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.

*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.

*R*

^{2}×

*S*

^{1}, 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.

- Filter each image of the database with a set of oriented edge detection kernels.
- 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. - 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*^{2}× S^{2}. - Project the data of the 4D histogram to a 3D histogram where the third coordinate is the relative 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).

*θ*. 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.

*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*.

*x*, Δ

*y*,

*θ*

_{ i },

*θ*

_{ j }) of the 4D histogram on the rotated point:

*x*, Δ

*y*) of an angle Δ

*θ*.

*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).

*H*(

*η*,

*ξ*, Δ

*θ*).

*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*

^{0}and

*G*

^{ π/2}are the first

*x*-derivative and

*y*-derivative of a Gaussian function, respectively:

*G*

^{0}and

*G*

^{ π/2}. Therefore, only two 2D convolutions are needed to generate the whole set of filtered images. Finally, we observe that

*G*

^{0}and

*G*

^{ π/2}are both separable kernels:

*σ*= 1 and choose a support of 7 pixels.

*π*/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.

*x*,

*y*,

*θ*) ∈

*R*

^{2}×

*S*

^{1}, 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.

*R*

^{2}×

*S*

^{1}

*γ*

_{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*

^{1}.

*S*

^{1}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*

^{2}×

*S*

^{1}, and the 2D curve

*γ*

_{2D}(

*t*) is lifted to a new curve

*γ*(

*t*) in the 3D space:

*R*

^{2}×

*S*

^{1}if it is the lifting of an edge (identified with a planar curve).

*γ*(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

*γ*′(

*t*) can be written as a linear combination of these two vectors:

*X*_{1}+

*k*

*X*_{2}. We have chosen to set the coefficient of

*X*_{1}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*_{1},

*X*_{2}. 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*

^{2}×

*S*

^{1}will be identified with the rototranslation group, the group of rotations and translations.

*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.

*R*

^{2}×

*S*

^{1}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*

_{0}= (

*x*

_{0},

*y*

_{0},

*θ*

_{0}) and

*p*

_{1}= (

*x*

_{1},

*y*

_{1},

*θ*

_{1}):

*θ*depends on (

*x*,

*y*), as shown in Equation 6.

*N*(0,

*σ*

^{2}) is a normally distributed variable with zero mean and variance equal to

*σ*

^{2}. 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.

*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

*σ*

^{2}. When

*σ*

^{2}= 0, the motion is completely deterministic and particles never deviate from straight paths. When

*σ*

^{2}→ ∞, 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):

*X*_{1}= (cos(

*θ*), sin(

*θ*), 0), the direction tangent to the path and a diffusion term on the orientation variable

*X*_{2}= (0, 0, 1).

*σ*, the standard deviation in curvature.

*R*

^{2}×

*S*

^{1}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}= cos(

*θ*) ∂

_{ x }+ sin(

*θ*) ∂

_{ y }and

*X*

_{2}= ∂

_{ θ }since are the partial derivatives in the direction of the vectors

*X*_{1}and

*X*_{2}, respectively.

*σ*in edges belonging to natural images.

*H*(

*x*,

*y*,

*θ*).

*H*, i.e., given a pair of contour points it is sufficient to switch the reference between them to verify the central symmetry.

- The boundary conditions are Neumann in the spatial coordinates and periodic in the directional coordinate.
- 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.
- The corresponding linear system solver uses a direct method based on the library UMFPACK.

*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*

^{2}×

*S*

^{1}. The best fitting between the model and the estimated probability allows to measure the parameter

*σ*of the model.

*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*

_{0},

*y*

_{0},

*θ*

_{0}) = (

*x*,

*y*,

*θ*):

*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*

^{2}×

*S*

^{1}for different values of

*k*are shown. The projection on the space

*R*

^{2}is visualized in Figure 6.

*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

*θ*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.

*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).

*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.

*θ*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).

*σ*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

^{−1}. 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).