A model of natural image edge co-occurrence in the rototranslation group

Abstract

In this paper, we propose to model the edge information contained in natural scenes as points in the 3D space of positions and orientations. This space is equipped with a strong geometrical structure and it is identified as the rototranslation group. In this space, we compute a histogram of co-occurrence of edges from a database of natural images and show that it can be interpreted as a probability density function, expressed by the fundamental solution of a suitable Fokker–Planck equation defined in the 3D structured space. Both estimated statistics and model predictions are reconsidered and compared with the partial gestalt association fields proposed by D. J. Field, A. Hayes, and R. F. Hess (1993). Finally, parametric identification allows to estimate the variance of the co-occurrence random process in natural images.

Introduction

Early visual cortex provides low-level feature extraction in terms of local frequency and orientation analysis. It is evident that cell responses with Gabor-like profile in the visual cortex do exist (Jones & Palmer, 1987). They provide an estimation of the position and local orientation of the image gradient and a subsequent process link this information to form integral contours. 

From the psychophysical point of view, the integration problem has been faced for example by Yen and Finkel (1998), who proposed a facilitation field for contour integration, and by Field, Hayes, and Hess (1993), who proposed an association field to model the link between position-orientation information. 

From a different point of view, many experiments have been conducted to establish the relationship between position and orientation of gradients in natural images. Experiments aimed to estimate the statistical distribution of co-occurrence of position-orientation couples from large data base of natural images. In Elder and Goldberg (2002), it has been shown that statistics of co-occurrence is invariant with respects to scale and that a co-circularity constraint underlies the distribution. Many subsequent works have been done in this direction (August & Zucker, 2000; Elder & Goldberg, 2002; Geisler, Perry, Super, & Gallogly, 2001; Gilbert, Sigman, Cecchi, & Magnasco, 2001; Orabona, Metta, & Sandini, 2006). Two good review papers are by Simoncelli and Olshausen (2001) and by Geisler (2008). 

Mathematical representations of edge distributions and link have been provided by many authors (August & Zucker, 2000; Mumford, 1994; Petitot & Tondut, 1999; Williams & Jacobs, 1997, among others). In Citti and Sarti (2006), a couple of position orientation is naturally described as a point in the rototranslation group (SE(2)), and the links between this points are the integral curves in the group, providing a simple model for association fields. Stochastic models for boundary completion have been proposed in Franken, Duits, and ter Haar Romeny (2007), Mumford (1994), Williams and Jacobs (1997), and August and Zucker (2000), substituting random walks to the deterministic integral curves. 

In this paper, we first define an experimental setup for the estimation of edge statistics in natural images and compute histograms of co-occurrence of edges from a database of natural images. Then, we propose to model the experimental data as points in the rototranslation group. This allows us to describe the link between edges in a very precise way with sophisticated mathematical instruments. From a deterministic point of view, different edges are connected by the integral curves that are solutions of the natural advection (transport) equation in the group. In a stochastic setting, connectivity of points in the group becomes a probability of co-occurrence, regulated by the stochastic advection equation in the group, i.e., a Fokker–Planck equation. The model equations, while discretized and simulated, are able to predict qualitative and quantitative behavior of the probability density of edges co-occurrence in natural images. 

Parametric identification allows to estimate the variance of the co-occurrence random process in natural images. 

In the Natural image edge co-occurrence statistics section, we introduce a novel methodology to compute the edge co-occurrences probability in natural images. In The mathematical model section, we describe the deterministic and stochastic mathematical models for the edge co-occurrences. The results and comparisons between experimental data and numerical simulation are presented in the Results section. Parameter identification of the variance of the stochastic process is also provided. 

Natural image edge co-occurrence statistics

In this section, we describe a new technique for computing the co-occurrence probability of edge orientations in a database of natural scenes. We have improved a standard approach used by other research groups over the last few years (August & Zucker, 2000; Elder & Goldberg, 2002; Geisler et al., 2001; Gilbert et al., 2001; Krüger, 1998; Orabona et al., 2006). In all these works, authors use pairs of quadrature filters (Morrone & Burr, 1988) to obtain an energy function for each image. The image energy is the sum of the square of the response of an odd and an even filter, resulting in a value independent of the polarity. 

In the sequel, we will first describe the mechanism for estimating the probability of co-occurrences, and then we will provide the implementation details. 

Our approach

The probability of co-occurrences is obtained computing a 3D histogram where the two first dimensions correspond to the relative position of two edges and the final one to their relative orientation. Then, the histogram lies in the space R ² × S ¹, the same space we will later use to model contours in The mathematical model section. In contrast with previous papers on co-occurrence statistics, we focus only on the image contours and their orientations and use oriented derivative filter so that the polarity of the contrast of contours is considered. 

Our method can be summarized in four steps:

1.  
   Filter each image of the database with a set of oriented edge detection kernels.
2.  
   For each image, perform non-maximal suppression with a threshold in order to obtain a list of points (x, y, θ) corresponding to edges with their respective orientations.
3.  
   Count how many times two detected edges with relative position (Δx, Δy) have orientations (θ _c, θ _p) and store the data in a 4D histogram in R ² × S².
4.  
   Project the data of the 4D histogram to a 3D histogram where the third coordinate is the relative orientation.

Step 1: We focus our study only on the image contours and their orientations; therefore, the edge detection filter will be an odd kernel acting as an oriented derivative filter. This choice allows us to consider the polarity of the contrast of contours, distinguishing if the image gradient goes from a darker zone to a lighter one or in the opposite sense. Therefore, the considered set of filters is indexed by their orientation θ ∈ [0, 2π). Among the options present in literature, one of the most common is the imaginary component of Gabor filters. Instead, we used directional derivatives of a Gaussian filter (DoG) with orientation θ since using the steerable filter architecture we achieve very efficient computational implementation (Freeman & Adelson, 1991). 

Step 2: For every image the filtering process produces a stack of filtered images, each one obtained by convolution with the kernel of orientation θ. We perform non-maximal suppression consisting in selecting for each pixel the maximum output of the filters as θ varies. Then we construct a list of triplets containing the pixels (x, y) where the maximum exceeded a fixed threshold, and the corresponding orientation θ where the maximum is achieved. A careful analysis ensures that the results are almost independent from the threshold. The same happens with the size and variance of the kernel used. 

This edge detection mechanism is inspired from the architecture of the primary visual cortex since we try to reproduce the hyper-columnar architecture using the bank of DoG filters. Each filter imitates the impulse response of the simple cells of V1. Also the non-maximal suppression is present in the hyper-columns (Hubel, 1995) (Figure 1). 

Figure 1

 

Geometric configuration of two edges in an image. The black bars represent the first edge at position (x _i , y _i ) with orientation θ _i and the second one at position (x _j , y _j ) with orientation θ _j . The co-occurrence final histogram is defined respect to the (η, ξ, Δθ).

Step 3: For each image, we have computed a list of points (x _i , y _i , θ _i ) that represents the set of positions of contours and their orientations, and we can now estimate statistics of co-occurrences. Given any pair of points (x _i , y _i , θ _i ) and (x _j , y _j , θ _j ) in the list, we say that we have a co-occurrence of two points with relative positions Δx = x _j − x _i , Δy = y _j − y _i and with orientations θ _i and θ _j . Hence, every co-occurrence is represented by a quadruplet (Δx, Δy, θ _i , θ _j ), which we store in a 4D histogram. In other words, we count how many times the quadruplets occur. In this procedure, we only take into account couples of oriented points satisfying ∣Δx∣, ∣Δy∣ < d. 

The same procedure is repeated for all considered images, accumulating in a 4D histogram. 

Step 4: We consider relative orientations and not absolute ones, projecting the points (Δx, Δy, θ _i , θ _j ) of the 4D histogram on the rotated point: 

$(\eta ,\xi ,\mathrm{\Delta }\theta )=({\mathrm{R}}_{\mathrm{\Delta }\theta }(\mathrm{\Delta }x,\mathrm{\Delta }y),{\theta }_j-{\theta }_j),$

(1)

where is the matrix of rotation of the vector (Δx, Δy) of an angle Δθ. 

The assumption that the position differences and the orientation differences suffice for compressing what is a 6D distribution (two (x, y) positions and two orientations) to a 3D histogram amounts to a spatial-orientation homogeneity assumption formulated and tested in August and Zucker (2000). 

This step has advantages for the computational complexity of the algorithm since it reduces the dimension of the data to be stored. As we will see in the next section, it also simplifies the geometrical tools to be used in the model. 

Finally, we normalize the full histogram over the total number of occurrences to get the probability of observing an edge element at every possible relative position and orientation difference from a given (reference) edge element. This normalized histogram estimating the probability density function of co-occurrences of a particular geometrical configuration of contours will be called from now on H(η, ξ, Δθ). 

Numerical implementation

The image database was obtained from the Web site: http://hlab.phys.rug.nl/imlib/index.html (van Hateren, 1997; van Hateren & Snippe, 2001) and it has been many times used in literature to compute natural image statistics (Gilbert et al., 2001; Kalkan, Wörgötter, & Krüger, 2007). It consists on 4000 high-quality gray scale digital images, 1536 × 1024 pixel and 12 bits in depth. Every convolution kernel G ^θ is an oriented derivative of Gaussian (DoG) filter with direction θ. They can be written as a linear combination of two derivative of Gaussian filters (Freeman & Adelson, 1991): 

$G^{\theta }=\cos \left(\theta \right)G^0+\sin \left(\theta \right)G^{\pi /2},$

(2)

where G ⁰ and G ^π/2 are the first x-derivative and y-derivative of a Gaussian function, respectively: 

$G^0=-2x{e}^{-(x^2+y^2)},G^{\pi /2}=-2y{e}^{-(x^2+y^2)}.$

(3)

 

This means we can obtain a filtered image corresponding to any oriented kernel by taking linear combinations of the images filtered with G ⁰ and G ^π/2. Therefore, only two 2D convolutions are needed to generate the whole set of filtered images. Finally, we observe that G ⁰ and G ^π/2 are both separable kernels: 

$G^0=-2\left(x{e}^{-x^2}\right)e^{-y^2},G^{\pi /2}=-2e^{-x^2}\left(y{e}^{-y^2}\right),$

(4)

reducing the 2D convolutions to the concatenation of two 1D filters. Making a quick analysis, it turns out that the resulting histogram is practically independent from the size and variance of the kernel used. We have set σ = 1 and choose a support of 7 pixels. 

The orientations were discretized at 32 different values (−15π/16, −14π/16, …, −π/16, 0, π/16, … 15π/16, π) and considered contours far up to 32 pixels so the final dimension of the histogram is 65 × 65 × 32. 

The total computation time for each image depends on the quantity of contours found on it. The mean value registered for the 4000 images in the database was 10 seconds so the total computational time was around 10 hours. This computation was performed on a Core 2 Duo 2.1-GHz laptop computer and the programs were implemented in C language. 

Visualization of the histogram is provided in the Results section. 

The mathematical model

The goal of this section is to present a new model for the probability of co-occurrence of edges in natural images. We first represent edges as curves in 3D space of positions and orientations (x, y, θ) ∈ R ² × S ¹, and recall that it can be identified with the rototranslation group with a suitable geometrical structure. Then we refer to classical probabilistic models of edge description in this space and finally we present our model. 

Representing edges in R ² × S ¹

We will introduce a natural scenes edge representation using instruments of differential geometry. A contour is represented in the 2D plane as a regular curve γ _2D(t) = (x(t), y(t)), and almost everywhere we can assume that its tangent vector is non vanishing, so that it is identified by an orientation θ(t): D ⊂ R → S ¹. S ¹ represents [0, 2π), isomorphic to the unit circle. The function θ takes values in the whole circle in order to represent polarity of the contours: two contours with the same orientation but with opposite contrast are represented through opposite angles on the unit circle. In this way, the points of the curve will be described in the space (x, y, θ) ∈ R ² × S ¹, and the 2D curve γ _2D(t) is lifted to a new curve γ(t) in the 3D space: 

$(x\left(t\right),y\left(t\right))\to (x\left(t\right),y\left(t\right),\theta \left(t\right)).$

(5)

We will call admissible curve a curve in R ² × S ¹ if it is the lifting of an edge (identified with a planar curve). 

In Figure 2, we illustrate the lifting process. By definition, the tangent vector to the curve γ _2D (the blue curve in Figure 2) can be expressed as (cos(θ), sin(θ)), so that 

$\theta =\arctan (y^{\prime }/x^{\prime }).$

(6)

 

Figure 2

 

A contour represented by the blue curve γ _2D(t) is lifted into the rototranslation group obtaining the red curve γ(t). The tangent space to the rototranslation group is spanned by the vectors X ₁ and X ₂.

Let us now consider its lifting γ (the red curve in Figure 2). By definition, the projection on the x, y plane of its tangent vector is tangent to the curve γ _2D and it has direction (cos(θ), sin(θ)). γ′(t) has a non vanishing component in the direction 

$X_1=(\cos \left(\theta \right),\sin \left(\theta \right),0)$

(7)

and a second component in the direction of the lifting: 

$X_2=(0,0,1).$

(8)

Then γ′(t) can be written as a linear combination of these two vectors: X ₁ + k X ₂. We have chosen to set the coefficient of X ₁ to one since it cannot vanish. In particular, admissible curves are integral curves of two vector fields in a 3D space and cannot have components in the orthogonal direction (Figure 3). 

$X_3=(-\sin \left(\theta \right),\cos \left(\theta \right),0).$

(9)

We will call horizontal tangent space the space generated by the vectors X ₁, X ₂. In differential geometry (Do-Carmo, 1976), these vectors can be identified by the generators of the Lie algebra of rotations and translations. Hence, the space R ² × S ¹ will be identified with the rototranslation group, the group of rotations and translations. 

Figure 3

 

The horizontal tangent planes in each point of the rototranslation group determine the differential structure of the space.

An admissible curve in this group is an integral curve of the vector field X ₁, X ₂ and is defined as the solution of the following ordinary differential equation: 

${\gamma }^{\prime }\left(t\right)=X_1\left(t\right)+k\left(t\right)X_2\left(t\right).$

(10)

 

Writing the same equation in coordinates we get: 

$\begin{array}{c}{x}^{\prime }=\cos \left(\theta \right)\\ y^{\prime }=\sin \left(\theta \right)\\ {\theta }^{\prime }=k\left(t\right)\end{array}.$

(11)

 

It is well known that k(t) is the curvature of γ _2D(t). Writing the curve in this way, it becomes obvious that the shape of the curve is completely defined by the function k: its curvature or angular velocity. 

If we equip the tangent planes with an Euclidean metric, then the length of an admissible curve can be computed as usual: 

$\lambda \left(\gamma \right)\left(t\right)=\underset{0}{\overset{t}{\int }}\parallel {\gamma }^{\prime }\left(s\right)\parallel ds=\underset{0}{\overset{t}{\int }}\parallel X_1+k{X}_2\parallel ds=\underset{0}{\overset{t}{\int }}\sqrt{1+k^2}ds.$

(12)

In order to define a distance in term of the length, we need to answer the following question: Is it possible to connect every pairs of points of R ² × S ¹ using an admissible curve? This is not a simple question since admissible curves are integral curves of a linear combination of two vectors fields in a 3D space. However, the answer is yes. See Citti and Sarti (2006) for a detailed justification. Consequently, it is possible to define a notion of distance between two points p ₀ = (x ₀, y ₀, θ ₀) and p ₁ = (x ₁, y ₁, θ ₁): 

$d(p_0,p_1)=\inf \{\lambda \left(\gamma \right):\gamma \mathrm{i}\mathrm{s}\mathrm{a}\mathrm{n}\mathrm{a}\mathrm{d}\mathrm{m}\mathrm{i}\mathrm{s}\mathrm{s}\mathrm{i}\mathrm{b}\mathrm{l}\mathrm{e}\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{v}\mathrm{e}\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{n}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{n}\mathrm{g}{p}_0\mathrm{a}\mathrm{n}\mathrm{d}{p}_1\}.$

(13)

In the Euclidean case, this infimum is realized by a geodesic that is a segment. Also here the infimum is realized on a curve called geodesics, which however can be locally curvilinear. 

The metric induced by Equation 13 is clearly non-Euclidean. Moreover, it is not even Riemannian and it is called Sub-Riemannian. We refer to Citti and Sarti (2006) for a detailed technical explanation. 

Metrics of this type often appear when one of the dimensions is a state variable depending on the others, and in this case the state variable θ depends on (x, y), as shown in Equation 6. 

The functional Equation 12 has the same asymptotic behavior as the modified elastica functional introduced by Nitzberg, Mumford, and Shiota (1993). It is quadratic when the curvature tends to zero and linear when it tends to infinity. The same functional is used in perceptual completion problems to complete occluded parts, to perform image inpainting, or to retrieve subjective contours (Ballester, Bertalímo, Caselles, Sapiro, & Verdera, 2001; Citti & Sarti, 2006; Masnou & Morel, 1998; Nitzberg & Mumford, 1990; Sarti, Citti, & Petitot, 2008). 

A classical stochastic model for edge description

David Mumford in Field et al. (1993) proposed to model natural scene edges as random walks in the plane where the curvature is taken as a random variable. This process is modeled by the following stochastic differential equation (SDE) (also used in August and Zucker, 2003, and Williams and Jacobs, 1995): 

$\begin{array}{c}{x}^{\prime }=\cos \left(\theta \right)\\ y^{\prime }=\sin \left(\theta \right)\\ {\theta }^{\prime }=N(0,{\sigma }^2).\end{array}$

(14)

where N(0, σ ²) is a normally distributed variable with zero mean and variance equal to σ ². Note that this is the probabilistic counter part of the deterministic Equation 11, naturally defined in the group structure. Indeed, both systems are represented in terms of left invariant operators of the Lie group (Citti & Sarti, 2006), the first one with deterministic curvature and the second with normal random variable curvature. 

These equations describe the motion of a particle moving with constant speed in a direction randomly changing accordingly with the stochastic process N. The effect is that particles tend to travel in straight lines, but over time, they drift to the left or right by an amount proportional to σ ². When σ ² = 0, the motion is completely deterministic and particles never deviate from straight paths. When σ ² → ∞, the motion is completely random, and the stochastic process becomes a 2D isotropic random walk. For this reason, the stochastic process defined above is called diffusion and the deterministic process advection. We denote p the probability density to find a particle at the point (x, y) moving with direction θ at the instant of time t conditioned by the fact that it started from a given location with some known velocity. This probability density satisfies a deterministic equation known in literature as the Kolmogorov forward equation or Fokker–Planck equation (FP) (Oksendal, 2005): 

${\partial }_{t}p=-\cos \left(\theta \right){\partial }_{x}p-\sin \left(\theta \right){\partial }_{y}p+({\sigma }^2/2){\partial }_{\theta }^{2}p.$

(15)

 

On this formulation, the FP equation consists on an advection term in the direction X ₁ = (cos(θ), sin(θ), 0), the direction tangent to the path and a diffusion term on the orientation variable X ₂ = (0, 0, 1). 

This equation has been largely used in computer vision and applied to perceptual completion related problems. It was first used by Williams and Jacobs (1997) to compute stochastic completion field, by August (2001) and August and Zucker (2003), to define the curve indicator random field, and more recently by Duits and Franken (2007) and Franken et al. (2007), applying it to perform contour completion, denoising, and contour enhancement. 

The only unknown parameter of the model is the variance σ, the standard deviation in curvature. 

A time independent Fokker–Planck equation for the co-occurrence of oriented edges

We propose here to use the FP equation for modeling co-occurrence probability of edge orientation. The main novelty of our model is that the proposed FP equation independent of the time variable, together with the fact that we consider both forward and backward propagations. 

The fundamental solution of the proposed FP equation in R ² × S ¹ express the probability to find in natural images an oriented contour at a certain position conditioned by the fact that there was an oriented contour in a reference point. On the other hand, each point of a contour can be described by a particle that can go either backward or forward with identical probability in the direction of the contour. The model for the probability density propagation is then the sum of two FP Green functions, one corresponding to a particle moving forward, 

$-X_{1}p(x,y,\theta )+({\sigma }^2/2)X_{22}p(x,y,\theta )=(1/2)\delta (x,y,\theta ),$

(16)

and the other corresponding to a particle moving backward: 

$X_{1}p(x,y,\theta )+({\sigma }^2/2)X_{22}p(x,y,\theta )=(1/2)\delta (x,y,\theta ).$

(17)

 

In the two equation above, we have written the Fokker–Planck operator using the notation of Citti and Sarti (2006), calling the differential operators X ₁ = cos(θ) ∂ _x + sin(θ) ∂ _y and X ₂ = ∂ _θ since are the partial derivatives in the direction of the vectors X ₁ and X ₂, respectively. 

One of the main goals of the present study is to estimate the value of the standard deviation in curvature σ in edges belonging to natural images. 

Numerical computation of the fundamental solution of the FP equation

We have numerically computed the fundamental solution of the FP Equations 16 and 17, using COMSOL Multiphisics v3.5, a Finite Element Method solver. Then we have compared the sum of the 2 fundamental solutions Equations 16 and 17 with the histogram H(x, y, θ). 

Let's notice that the sum of fundamental solutions are symmetrical respects to the origin by construction. This property is also verified by the histogram H, i.e., given a pair of contour points it is sufficient to switch the reference between them to verify the central symmetry. 

In the following, we list some details of the numerical simulation:

1.  
   The boundary conditions are Neumann in the spatial coordinates and periodic in the directional coordinate.
2.  
   The distribution δ(x, y, θ) is numerically approximated by a Gaussian normalized function: 

   ${\delta }_{\epsilon }(x_i,y_i,{\theta }_i)={\left(\frac{1}{\epsilon \sqrt{\pi }}\right)}^3{e}^{-\frac{1}{{\epsilon }^2}(x_i^2+y_i^2+{\theta }_i^2)}$

   (18)

   where (x _i , y _i , θ _i ) are mesh points.
3.  
   The solution is computed on an adaptive mesh. Later we interpolate the mesh in a rectangular grid (65 × 65 × 32 large) in order to compare it with the histogram.
4.  
   The corresponding linear system solver uses a direct method based on the library UMFPACK.

Results

The histogram constructed in the Natural image edge co-occurrence statistics section and the theoretical model presented in The mathematical model section have been introduced in a completely independent way. The main scope of this section is to compare the two approaches from both a qualitative and a quantitative point of view. We will show that the geometry formally introduced in The mathematical model section is naturally encoded in the histogram, first deducing from the histogram the integral curves of the vector fields X _i . After that, as proposed in The mathematical model section, we show that the fundamental solution of the FP operator models the histogram in the rototranslation group, isomorphic to R ² × S ¹. The best fitting between the model and the estimated probability allows to measure the parameter σ of the model. 

The co-circularity pattern

One of the most important and studied pattern in natural image statistics of the possible geometrical contour configuration is the co-circular pattern, introduced in this field for the first time by Parent and Zucker (1989). This property has an interesting psychophysic counterpart, i.e., the association fields, introduced by Field et al. (1993) as a model explaining the good continuation law of Gestalt psychology. They used oriented Gabor-like patches distributed over a grid to prove that alignment of the elements along a path plays a large role in the ability to detect the path. Small variations in the alignment with respect to the path significantly reduced the ability of the observer to detect it. From these experiments, they construct the field shown in Figure 4. 

Figure 4

 

The field lines of the association field (Field et al., 1993).

On the other hand, from the neurophysiological point of view, it is largely accepted that the neural correlate of these association fields are the long range intra-columnar connections (Ben-Shahar & Zucker, 2004; Bosking, Zhang, Schofield, & Fitzpatrick, 1997; Petitot & Tondut, 1999). These connections exist between simple cells of V1 who share almost the same orientation and are responsible for contour integration and completion of subjective contours. Many perceptual completion models inspired by this cortical architecture explain the association fields as co-circular connections. These associations fields have been modeled as integral curves of the rototranslation group in Citti and Sarti (2006). 

Integral curves in the rototranslation group

The expression of the integral curves of the structure is provided by Equation 11. Choosing the curvature k constant over time, one can integrate the system of equations to obtain an explicit formula for the curve γ(t) starting from the origin (x ₀, y ₀, θ ₀) = (x, y, θ): 

$\begin{array}{c}x=(1/k)\sin (kt)\\ y=(1/k)(1-\cos (kt))\\ \theta =kt.\end{array}$

(19)

 

For each value of the parameter k, it is obtained a different curve with helicoidal path since the projection over the plane is a circle with curvature k and it moves with constant speed in direction θ. On Figure 5, some integral curves in the 3D space of positions and orientations R ² × S ¹ for different values of k are shown. The projection on the space R ² is visualized in Figure 6. 

Figure 5

 

The fan of integral curves given by Equation 19 by varying the curvature k, visualized in the R ² × S ¹ space. The set of curves defines a surface with a singular point in the origin.

Figure 6

 

The fan of integral curves projected on the xy plane. On the left are visualized the projections of the curves defined by Equation 19 showing clearly the co-circularity constraint. On the right is visualized the vector field of unitary vectors oriented with the maximal probability (red) with superimposed its integral curves (blue). Again the co-circularity constraint is evident.

Circular curves deduced by the histogram

The co-circular property of the edge statistics was observed by several groups (August & Zucker, 2000; Elder & Goldberg, 2002; Geisler et al., 2001; Gilbert et al., 2001; Orabona et al., 2006). In our work, we strongly focus on curves, and the co-circularity pattern is shown by means of a family of curves reproducing the association fields. Indeed, we will deduce a 2D vector field from the histogram, compute its integral curves, and compare them with the previously detected family of curves, again through integral curves. 

The histogram is a discrete approximation of the probability density function for the co-occurrence of edges. In each point (x, y, θ), it expresses the probability to have a contour conditioned to the fact that there exists a horizontal contour at the reference point. As in Gilbert et al. (2001), we select for each point (x, y) the most probable orientation θ _m , such that 

$H(x,y,{\theta }_m(x,y))=\max \{H(x,y,\theta ),\theta \in {\mathrm{S}}^1\}.$

(20)

 

This orientation set the direction of the unitary vector field: 

$\mathbf{V}(x,y=(\cos \left({\theta }_m(x,y)\right),\sin \left({\theta }_m(x,y)\right)),$

(21)

depicted in Figure 6. Its integral curves can be computed solving the differential equation: 

$\begin{array}{c}{\gamma }^{\prime }\left(t\right)=\mathbf{V}(x,y)\\ \gamma \left(0\right)=(x_0,y_0),\end{array}$

(22)

and they are visualized in Figure 6. 

As it is clear from the result, the integral curves of the vector field V(x, y) (Figure 6 right) optimally approximate the 2D projection of the integral curves of the vector fields X ₁, X ₂ (Figure 6 left), and both show a co-circular pattern modeling the association fields of Figure 4. 

The probabilistic interpretation

Qualitative comparison between the histogram and the fundamental solution of the FP equation

Figure 8 shows a visualization of the resulting histogram (blue) and the computed fundamental solution of the FP (red). The isosurface, i.e., the surface defined by a constant intensity level of the 3D function, corresponding to the 2% of the maximum has been selected and a surface rendering visualization has been adopted. The qualitative resembling of the two data sets is evident. Both the distributions are thick versions of the surface generated by the set of the integral curves shown in Figure 5. In Figure 9 and in supplementary material, the torsion of the isosurfaces is visible in both data sets. 

The comparison between the two distributions is visualized with a different technique in Figures 10 and 11, where isocontour plots at constant θ section are shown. Every θ constant slice corresponds to a configuration of two contours differing by θ in their orientations. The drop-off in correlation strength with distance present in Figures 10 and 11 is strictly related to both the spread of the integral curves of the deterministic model, and the diffusion term in defined by the variance σ. In Figure 12, it is shown a diagram similar to the one used in Gilbert et al. (2001) to interpret the distribution of edges for each possible configuration. Each black bar corresponds to an edge supposing there was a horizontal edge in the center. If the hypothesis of co-circularity is verified, edges with the same orientation are distributed around a straight line whose angle direction is half the edge orientation angle. In contrast with the earlier work, there are only two branches in the cut planes of the histogram instead of four since we are considering polarity of the contours. 

In Figure 7, it is plotted in blue the amount of co-occurrences of the histogram H for each θ constant slice and in red the integral of the FP distribution at a constant θ slice after performing a suitable rescaling explained below. Both distributions are mainly superimposed for all the values of the θ angle strictly less than ∣π/2∣ and different than 0 (collinear configuration). 

Figure 7

 

Comparison between the number of co-occurrences in every θ constant slice in H by varying the angle (blue) and the corresponding value in the FP distribution.

Over 10% of the co-occurrences of H corresponds to collinear edges (θ = 0). This is explained by the fact that in the natural image database, there are many sharp straight edges corresponding mostly to buildings. The same argument explains the two peaks present at θ = ±π/2, orthogonal contours. 

As outlined by Figure 6 in the paper of Simoncelli and Olshausen (2001) where they modeled the same type of distribution, the distributions had sharp peak at 0 and much longer tails than a Gaussian density. Field (1987) argued that the representation corresponding to these densities, in which most neurons had small amplitude responses, had an important neural coding property, which he termed sparseness. Our invariant model still underestimates the peak in 0, while it perfectly fits the tails of low intensities. The underestimation in 0 can be probably be explained with a violation of the rotation invariance of the image we postulated here. Indeed, the probability density is stronger in the collinear direction. The range of considered values of θ excludes automatically the presence of parallelism effects, which is statistical significant for larger value of theta, as shown in Krüger (1998) as well in Geisler et al. (2001). The modelization of parallelism in terms of Lie groups theory has been faced at a purely theoretical level in Sarti et al. (2008) and in terms of image statistics will be studied in Sanguinetti, Citti, and Sarti (in preparation). 

Best fitting and parameter identification

The best fitting between the experimental and simulated distributions has been accomplished by minimizing the mean square error by varying the parameter σ of Equations 16 and 17 (the variance of the stochastic process). With a discretization step of 0.01 into a range of 0 < σ < 10, the minimum error value results E _m < 2% showing at a quantitative level that the model represent accurately the experimental distribution. This prove the Mumford hypothesis on contour reappearance (Mumford, 1994). Moreover, the minimizer corresponding to E _m results σ = 1.73 pixels⁻¹. In the estimate of this parameter σ, we filtered out the artifact in the statistics, which is the value attained in 0, because it violates the invariance assumption we made, as discussed above. This value can be considered as the natural constant of the stochastic process underlying edge distribution in natural images (Figures 8, Figures 9, Figures 10, Figures 11, and Figures 12). 

Figure 8

 

Isosurface visualization of the histogram H (blue) and the FP fundamental solution (red) on orthogonal projection on the xθ plane (left) and yθ plane (right).

Figure 9

 

Isosurface visualization of the histogram H (blue) and the FP fundamental solution (red) with orthographic projection.

Figure 10

 

Isocontour visualization of the histogram H at different angles θ.

Figure 11

 

Isocontour visualization of the FP fundamental solution at different angles θ.

Figure 12

 

The co-circularity pattern as the most probable configuration. In contrast with the results presented in Gilbert et al. (2001), only two branches have been in the cut planes of the histogram in Figure 10 instead of four since we are considering polarity of the contours.

Conclusions

In this paper, we developed the connection between the rototranslation model (Citti & Sarti, 2006) for the interaction between oriented neurons in visual cortex and edge co-occurrence statistics. Inspired by Mumford (1994) and Williams and Jacobs (1995), a Fokker–Planck equation is introduced and its fundamental solution is used to bridge between the geometric model and the edge statistics. This approach can be formally extended to different Lie Groups corresponding to other low-level visual features like for example: curvature detection, ladder configuration to describe parallelism, or scale. Future research will check the relationship between these geometric structures and statistical properties of natural images. Some simplifications introduced in the present model to preserve local invariance can be relaxed in future studies to take into account the un-isotropy of different image features. These studies aim to achieve a deep understanding of the mutual relationship between statistics of images and geometrical models of the primary visual cortex in the perspective to identify underlying process of Gestalt principles. 

Supplementary Materials

Supplementary Movie - Supplementary Movie 

Supplementary Movie. Isosurface visualization of the histogram H (blue) and the FP fundamental solution (red). 

Supplementary Movie - Supplementary Movie 

Acknowledgments

This work was partially supported by GALA project Geometric Analysis in Lie groups and Applications and by MONESIA project MObility Network Europe-Southamerica: an Institutional Approach, both supported by the European Commission. We thank Steven Zucker for fruitful discussions. 

Commercial relationships: none. 

Corresponding author: Gonzalo Sanguinetti. 

Email: [e-mail]. 

Address: Viale Risorgimento 2, DEIS, CSITE Bologna, Bologna, Italy.