A subRiemannian model of the visual cortex with frequency and phase
Abstract
In this paper we present a novel model of the primary visual cortex (V1) based on orientation, frequency and phase selective behavior of the V1 simple cells. We start from the first level mechanisms of visual perception: receptive profiles. The model interprets V1 as a fiber bundle over the 2dimensional retinal plane by introducing orientation, frequency and phase as intrinsic variables. Each receptive profile on the fiber is mathematically interpreted as a rotated, frequency modulated and phase shifted Gabor function. We start from the Gabor function and show that it induces in a natural way the model geometry and the associated horizontal connectivity modeling the neural connectivity patterns in V1. We provide an image enhancement algorithm employing the model framework. The algorithm is capable of exploiting not only orientation but also frequency and phase information existing intrinsically in a 2dimensional input image. We provide the experimental results corresponding to the enhancement algorithm.
Keywords: SubRiemannian geometry, neurogeometry differential geometry Gabor functions visual cortex image enhancement
1 Introduction
The question of how we perceive has been an intriguing topic for different disciplines. One of the first school which faced the problem is the Berlin school of experimental psychology, called Gestalt psychology school, [81], [52], [51] which formulates precise laws which can explain visual perception. The Gestalt psychology is a theory for understanding the principles underlying the emergence of perceptual units, as the result of a grouping process. The main idea is that perception is a global phenomenon, which considers the scene as a whole, and is much more than the pure sum of local perception. The first perceptual laws are of qualitative type, based on similarity, closure, good continuation, alignment. After that there have been many psychophysical studies which attempted to provide quantitative version of the grouping process. With the developments of neuroscience studies, researchers started to look for cortical implementation of Gestalt laws, with a particular attention to neural architectures of the visual cortex. A particularly important one for our study is the pioneering work of Field et. al. [27], which models Gestalt principles of good continuation and alignment. They experimentally proved that fragments aligned along a curvilinear path can be perceived as a unique perceptual unit much better than fragments with rapidly changing orientations. The results of their experiments were summarized in a representation, called association fields, which represent the complete set of paths with fixed initial position and orientation which can be perceived as perceptual units. The visual cortex is a part of the mammalian brain which is responsible for the first level processing tasks of perceptual organization of local visual features in a visual stimulus (two dimensional image). It is known from neurophysiological experiments that the visual cortex contains neurons (simple cells) which are locally sensitive to several visual features, namely, orientation [38], [39], [37], [41], spatial frequency [60], [42], [44], [43], [76], [77], [69], [70], phase [19], [67], [55], [62], scale [8] and ocular dominance [75], [54], [44]. The simple cells are organized in a hypercolumnar architecture, which was first discovered by Hubel and Wiesel [40]. In this architecture, a hypercolumn is assigned to each point of the retinal plane (if we disregard the isomorphic cortical mapping between retinal and cortical planes), and the hypercolumn contains all the simple cells sensitive to a particular value of the same feature type. Simple cells are able to locally detect features of the visual stimulus, and neural connectivity between the simple cells integrates them in a coherent global unity. Those two mechanisms, the feature detection and the neural connectivity, comprise the functional geometry of V1.
Several models were proposed for the functional geometry of V1 associated to the simple cells which were only orientation sensitive. Early models date back to ’80s. Koenderink and van Doorn [50], [49] revealed the similarity between Gaussian derivative functions and simple cell receptive profiles. They proposed visual models based on the functions of Gaussian derivatives as the mathematical representations of the receptive profiles. Their findings indeed encouraged many studies relying on the choice of a family of Gaussian derivative functions and Gaussian kernels, among which we would like to mention the works of Young [82] and Lindeberg [57], [59].
A different modeling approach from the above mentioned ones was to employ Gabor functions as the mathematical representations of the orientation sensitive simple cell receptive profiles. The motivation for this choice was relying on an uncertainty principle as was elaborated by Daugman [18] through a generalization of the hypothesis of Marĉelja [61] (see also [45] where Jones and Palmer compared statistically the results obtained via Gabor functions and the neurophysiological results collected from V1 of a cat). Furthermore Hoffman (see [34], [35]) proposed to model the hypercolumnar architecture of V1 as a fiber bundle. Following the second school (which uses the Gabor functions) and by further developing the model proposed by Petitot and Tondut [63] (see also [64] and [65] of Petitot), where hypercolumnar architecture was interpreted as a fiber bundle associated to a contact geometry, Citti and Sarti [13] introduced a group based approach. They proposed a new model of the functional geometry of V1, which considered the subRiemannian geometry of the rototranslation group () as the suitable model geometry. The main reason for employing geometry was due to that the corresponding Lie algebra to was providing a good model of the actual neural connectivity in V1. This model proposed in [13] has been extended to other visual features in addition to orientation, such as scale by Sarti et. al. [72], and to other cell types such as complex cells sensitive to velocitiy and movement direction by Barbieri et. al. [1] and Cocci et. al. [15]. Apart from those, a semidiscrete model was presented by Prandi et. al. in [68]. Furthermore, image processing applications employing Gabor transform in order to extract visual features from medical images were proposed in [24] by Duits and Sharma (see also [74]). Other applications in medical image analysis employing scale and orientation information can be found in [11] and [46], where Gabor transform is employed for the detection of local frequencies in tagging MRI (magnetic resonance imaging) images and thus for the computation of local frequency deformations in those images. Interested reader can also refer to [26] for different applications of geometric approach in general, in computer vision and robotics. Additionally to those studies, the models in terms of cortical orientation and orientationfrequency selectivity, which were provided by Bressloff and Cowan [10], [9], could be useful references for the reader. We refer to [14] for a review of several cortical models including many of the above mentioned ones.
The theoretical criterion underpinning the modeling we propose in this paper relies on the so called neurogeometrical approach described by Citti and Sarti [13], Petitot and Tondut [63], Sarti et. al. [72]. Following this approach, processing capabilities of sensorial cortices, and in particular of the visual cortex are modeled based on the geometrical structure of cortical neural connectivity. Global and local symmetries of the visual stimuli are captured by the cortical structure which is invariant under those symmetries (see Sanguinetti et. al. [71]). We will follow a similar framework and we will start from the first level perceptual tasks performed by the simple cells, from local feature extraction. This starting point will lead us to the model geometry of V1 associated to the simple cells sensitive to orientation, spatial frequency and phase information at each position in a given two dimensional image.
At the level of Gestalt organisation, the neurogeometrical architecture in [13] implements the psychophysical law of good continuation, the architecture in the affine group [72] implements good continuation and ladder, the architecture in the Galilean group [1], [15] implements common fate, the architecture we are considering here in a Gabor based subRiemannian geometry implements similarity between textures/patterns and contains all the previous models employing the neurogeometrical approach.
Once the light reflects from a visual stimulus and arrives to the retina, it evokes some spikes which are transmitted along the neural pathways to the simple cells in V1. Each simple cell gives a response called receptive profile to those spikes. In other words, receptive profile is the impulse response of a simple cell. The simple cells extract the information of local visual features by using their receptive profiles and it is possible to represent the extracted features mathematically in a higher dimensional space than the given two dimensional image plane. We will call this space the lifted space or the lifted geometry. We will use an extended Gabor function as the receptive profile of the simple cells. We will see that this choice naturally induces the corresponding Lie algebra of the subRiemannian structure, which is the corresponding lifted geometry to our model. The Lie algebra and its integral curves model neural connectivity between the simple cells. Moreover, since some pairs of the algebra are not commutative, it is possible to formulate an uncertainty principle and this principle is satisfied by the extended Gabor function. That is, the extended Gabor function minimizes uncertainties arising from simultaneous detection of frequencyphase and simultaneous detection of positionorientation (see also [25, Section 7.5], [2], [3], [4] and [74] for similar phenomena in different frameworks).
Concerning the question of which family of functions to use as receptive profiles, let us recall that receptive field models consisting of cascades of linear filters and static nonlinearities may be adequate to account for responses to simple stimuli such as gratings and random checkerboards, but their predictions of responses to complicated stimuli (such as natural scenes) are correct only approximately. A variety of mechanisms such as response normalization, gain controls, crossorientation suppression, intracortical modulation can intervene to change radically the shape of the profile. Then any static and linear model for the receptive profiles has to be considered just as a very first approximation of the complex behavior of a real dynamic receptive profile, which is not perfectly described by any of the static wavelet frames.
For example derivatives or difference of Gaussian functions are suitable approximations of the behavior of classical receptive profiles of the simple cells. In [58, 59], Lindeberg starts from certain symmetry properties of the surrounding world and derives axiomatically the functions of Gaussian derivatives obtained from the extension of the family of rotationally symmetric Gaussian kernels to the family of affine Gaussian kernels, and proposes to model the simple cell receptive fields in terms of those Gaussian derivatives (see also Koenderink [50], [49], Young [82], Landy and Movshon [53]). Indeed Gaussian functions are good models of the receptive profiles if we restrict ourselves to the visual features except for frequency and phase. They provide good results for orientation and scale detection as shown by the scalespace school (see, e.g., the works of Lindeberg [56], [57], [59], Florack [28], ter Haar Romeny [79], [78], Hannink et. al. [33]). However, we are interested here in two dimensional visual perception based on orientation, frequency and phase sensitive simple cells. Differently from the case with orientationscale sensitive simple cells, frequencyphase sensitive simple cells cannot be modeled in a straightforward way by Gaussian derivative functions. A different order Gaussian derivative must be used for the extraction of each frequency component of a given image. This requires the use of different functions of each one of them corresponds to a certain frequency, thus to a certain order derivative. In other words, frequency is not a parameter as in the case of scale but each frequency corresponds to a different function. It is not possible to derive a natural geometry starting from the derivatives of the Gaussian and it is rather required to assign an adequate geometric setting to the set of extracted feature values by the Gaussian derivatives in order to represent those values. At this point, a Gabor function seems to be a good candidate for the detection of different orientation, frequency and phase values in a two dimensional image, since orientation, frequency and phase are parameters of the Gabor function. In other words, instead of using different functions, we can use a single Gabor function corresponding to a set of parameter values in order to detect different feature values. In this way, we obtain a subRiemannian model geometry as the natural geometry induced directly by the Gabor function (i.e., by the receptive profile itself). Moreover, the Gabor function is able to model both asymmetric simple cells and even/odd symmetric simple cells thanks to its phase offset term appearing in its wave content while the functions of the Gaussian derivatives account only for the symmetric simple cells. Considering those points, we propose to use a Gabor function with frequency and phase parameters as the receptive profile model. The Gabor function allows to extend the model provided in [13] to the true distribution of the profiles in V1 (including the asymmetric receptive profiles with the phase shifts) in a straightforward way. Finally, we would like to refer to Duits and Franken [21], [23], [22], Franken and Duits [31], Sharma and Duits [74], Zhang et. al. [83], Bekkers et. al. [7] for information about applications which employ other wavelets corresponding to unitary transforms for feature extraction.
Here we consider the model framework provided in [13] as the departure point of our study. We extend this model from orientation selective framework to an orientation, frequency and phase selective framework. Furthermore we provide the neural connectivity among the simple cells not only orientation selective but also frequency selective with different phases. Thanks to the use of all frequency components of the Gabor functions, Gabor transform can be followed by an exact inverse Gabor transform, which was not the case in the model presented in [13] since a single frequency component of the Gabor function was used. The projection of our generalized model onto can be considered as equivalent to the model provided in [13]. The procedure that we use to obtain the extended framework can be employed for the extension to a model associated with orientationscale selective simple cells as well (see [5]).
We will see in Section 2 the model structure. We will show how the model geometry with the associated horizontal connectivity can be derived starting from the receptive profile model, i.e., from the Gabor function. Then in Section 3 we will provide the explicit expressions of the horizontal integral curves, which are considered as the models of the association fields in V1. Finally in Section 4, we will provide an image enhancement algorithm using the model framework together with the results obtained by applying a discrete vesion of the algorithm on some test images.
2 The model
The model is based on two mechanisms. The first one is the feature extraction linear mechanism. The second mechanism is the propagation along horizontal connectivity, which models the neural connectivity in V1. We describe the model by using those two mechanisms in terms of both a group structure and a subRiemannian structure.
2.1 Feature extraction and representation
2.1.1 Receptive profiles, symplectic structure and contact form
Being inspired by the receptive profile models proposed in [13] for the orientation selective behavior and in [20], [16], [1] for the spatiotemporal behavior of the simple cells , we propose to represent the receptive profile of a simple cell in our setting with the Gabor functions of the type
(1) 
with the spatial frequency^{1}^{1}1Spatial frequency refers to with a wavelength in our terminology. and , where we represent the coordinates associated to a 6dimensional space with .
In the case of V1 complex cells with spatiotemporal dynamics, the variable represents the velocity of a twodimensional plane wave propagation (see Barbieri et. al [1] for details). However, we are not interested in the complex cells and any temporal behavior, and we can choose . In our framework we interpret as the phase centered at ,
In this way we obtain a 5dimensional space
(2) 
where denotes the feature variables . Then we may write the associated Gabor function which is centered at and sensitive to feature values by using (1) as follows:
(3) 
The standard Liouville form reduces to
(4) 
Indeed is a contact form since
(5) 
is a volume form. In other words it is maximally nondegenerate and it does not vanish at any point on the manifold .
2.1.2 Set of receptive profiles
An important property of Gabor functions is that they are invariant under certain symmetries. Therefore any Gabor function can be obtained from a reference Gabor function (mother Gabor function), up to a certain transformation law.
Let us denote the origin for the layer of a frequency by . Then a suitable choice of the mother Gabor function with the frequency is
(6) 
We set
(7) 
which describes at each frequency the relation between a generic receptive profile centered at and the mother Gabor function through
(8) 
The set of all receptive profiles obtained from the mother Gabor function with all possible combinations of feature values at each possible frequency is called the set of receptive profiles.
2.1.3 Output of a simple cell
We obtain the output response of a simple cell (which is located at the point and sensitive to the feature values ) to a generic image as a convolution with Gabor filter banks:
(9) 
We apply the convolution for all feature values at every point in order to obtain the output responses of all receptive profiles in the set of receptive profiles. It is equivalent to applying a multifrequency Gabor transform on the given two dimensional image. Since we use all frequency components of the transform, we can employ the exact inverse Gabor transform in order to obtain the initial image:
(10) 
with denoting the complex conjugate. We will call the output response lifted image and the Gabor transform lifting.
We remark here that we consider the whole complex structure of the result of the convolution (9) as the output response of a simple cell. It is different from the cases of the previous visual cortex models which were choosing either real or imaginary part of the output responses obtained as the result of the convolution with corresponding Gabor filters (see for example [13], [72], [73]). In other words they were not taking into account the half of the information which they obtained from an image. Furthermore, inverse Gabor transform was not possible in the previous models of visual cortex given in [13], [72], [73] since in those models a single frequency Gabor transform was employed to obtain the output responses.
2.2 Horizontal vector fields and connectivity
Horizontal vector fields are defined as the elements of
(11) 
where denotes the tangent bundle of the 5dimensional manifold . They are induced naturally by the 1form given in (4). The horizontal vector fields are found explicitly as
(12) 
The corresponding horizontal distribution is therefore as follows:
(13) 
All nonzero commutators related to the horizontal vector fields given in (12) follow as
(14) 
Note that the horizontal vector fields are bracket generating since
(15) 
for all , where denotes the tangent space of at . Obviously (15) shows that the horizontal vector fields fulfill the Hörmander condition [36], and consequently they provide the connectivity of any two points on through the horizontal integral curves defined on due to the Chow’s theorem [12]. This connectivitiy property is particularly important since it guarantees that any two points in V1 can be connected via the horizontal integral curves, which are the models of the neural connectivity patterns in V1.
2.3 Functional architecture of the visual cortex
2.3.1 The architecture as a Lie group
Receptive profiles evoke a group structure at each frequency . We can describe the group structure underlying the set of receptive profiles by using the transformation law given in (7).
First we notice that the elements induce the group given by
(16) 
which is indeed a Lie group associated to a fixed frequency .
The differential of the lefttranslation
(19) 
is given by
(20) 
The vector fields , and are bracket generating due to that
(21) 
for every . Hence , and generate the Lie algebra corresponding to .
2.3.2 The architecture as a subRiemannian structure
The functional geometry is associated to a subRiemannian structure at each frequency . We denote by a submanifold of with points restricted to a fixed . In this case the horizontal distribution is found by
(22) 
Furthermore the induced metric is defined on and at every point will make orthonormal.
Finally the associated subRiemannian structure to the frequency is written as the following triple:
(23) 
3 Horizontal integral curves
The lifting mechanism leaves each lifted point isolated from each other since there is no connection between the lifted points. Horizontal vector fields endow the model with an integration mechanism which provides an integrated form of the local feature vectors obtained from the lifted image at each point on .
Once a simple cell is stimulated, its activation propagates between the simple cells along certain patterns which can be considered as the integrated forms of the local feature vectors. This propagation machinery is closely related to the association fields [27], which are the neural connectivity patterns between the simple cells residing in different hypercolumns (long range horizontal connections) within V1. The association fields coincide with the anisotropic layout of the long range horizontal connections at the psychophysical level. In the classical framework of [13], those association fields were modeled as the horizontal integral curves of . We follow a similar approach and propose to model the association fields in our model framework as the horizontal integral curves associated to the 5dimensional subRiemannian geometry of . We conjecture that those horizontal integral curves coincide with the long range horizontal connections between orientation, frequency and phase selective simple cells in V1.
We denote a time interval by with and then consider a horizontal integral curve associated to the horizontal vector fields given in (12) and starting from an initial point . Let us denote the velocity of by . At each time the velocity is a vector at . In order to compute the horizontal integral curves, we first consider the vector field which is given by
(24) 
with coefficients (which are not necessarily constants) where . Then we can write each component of as follows:
(25) 
In the case of that the coefficients are real constants and , we solve the system of the ordinary differential equations (given in (25)) of with the initial condition and find the solution as follows:
(26) 
If then the solution becomes
(27) 
Note that (26) and (27) describe the whole family of the horizontal integral curves described by the horizontal distribution
We are interested rather in two specific subfamilies corresponding to the horizontal vector fields which reside in either one of the two orthogonal subspaces which are defined at every point as
(28) 
satisfying
(29) 
Figure 1 gives an illustration of the orthogonal layout of and at points on an orientation fiber, i.e., on a horizontal integral curve along corresponding to some fixed and . See also Figure 2, where the integral curves along the vector fields and with varied and values, respectively, are presented.
We remark here that is the horizontal tangent space of at the point once frequency and phase are fixed. In other words at each point with and fixed on , one finds the submanifold which is the classical subRiemannian geometry corresponding to the model given in [13]. This property allows the simple cell activity to be propagated in each subspace corresponding to a frequencyphase pair separately and it will be important for image enchancement applications employing our model framework.
4 Enhancement
Image enhancement refers to smoothing a given input image, reducing the noise and at the same time preserving the geometric structures (edges, corners, textures and so on). We perform our image enhancement procedure on the output responses instead of on the input image. Since the output responses encode the local feature values of orientation, frequency and phase, this allows us to exploit the additional information obtained from those features. Our enhancement procedure is based on an iterative LaplaceBeltrami procedure on the simple cell output responses in the 5dimensional subRiemannian geometry and it results in a mean curvature flow in the geometry.
4.1 Laplace Beltrami procedure
Anisotropic metric on the space of simple cell output responses defines the subRiemannian Laplacian in the subRiemannian space generated by the simple cells:
(30) 
where coefficients are nonnegative constants representing the weights of the second order horizontal vector fields which are given in (12). The weights are used to adjust the operator to the subRiemannian homogeneity of . They are particularly important in the discrete case, where different dimensions of the space need not necessarily be sampled in the same way.
It has been proved by Franceschiello et. al. in [29] that the output induces a metric on the space of the model geometry proposed in [13] and the metric elicits certain visual illusions. In the article of Franceschiello et. al. [29] a simplified diagonal metric was used. On the other hand, following the approach of Kimmel et. al. [48], [47], we choose the metric induced by the output on and use a simplified version of this metric for the applications.
The metric induced by the output responses is defined as follows:
Definition 1
(31) 
with constants .
We denote the inverse metric by and its elements by .
Mean curvature flow provides an adapted enhancement to the surface underlying the image function since the flow is restricted to the evolving level sets of the image. LaplaceBeltrami operator is written as:
(32) 
where is the determinant of the induced metric. LaplaceBeltrami operator can be considered as the linearization of the motion by curvature explained in [6]. For practial reasons, we will use a LaplaceBeltrami process with the operator given in (32) associated to a reduced version of the metric provided in Definition 1.
The evolution equation for the enhancement via subRiemannian LaplaceBeltrami procedure is written as:
(33) 
for all and .
4.1.1 Reduced equation
It is possible to perform the LaplaceBeltrami procedure in each frequency and phase subspace separately in a reduced framework. In that case we choose and . In this way we apply the evolution equation on surfaces in each frequency and phase subspace, i.e., on each manifold, which is the submanifold with elements representing the points with fixed and . In this framework the metric boils down to
(34) 
We choose and suitably by regarding the fixed values.
The motivation for choosing is that we would like to avoid excessive diffusion in the direction of the vector field . We already have sufficient diffusion in this direction due to the commutator . Direct application of introduces additional diffusion in ortogonal directions to the object boundaries, which is not desired since it might destroy object boundaries and contour structures in the input image. Furthermore, the use of the reduced version lowers the computational load since now multiple LaplaceBeltrami procedures are applied in 3dimensional subRiemannian geometry at each frequency instead of in the 5dimensional subRiemannian geometry .
We remark that the vector field does not perform information flow only in the orthogonal direction
(35) 
to the boundaries but also in the direction of phase. However the elimination of from the LaplaceBeltrami procedure must be accordingly taken into account also in the metric given in Definition 1 in order to provide the coherency between the LaplaceBeltrami operator and the employed metric. This is the reason for that we fix in the reduced version of the metric given in (34).
We also choose . Indeed we assume that no information flow takes place along the vector field , i.e., in the phase direction. We notice that the Gabor transform produces a rotated version of the image by the angle for each phase (see [24] for more details). Hence the LaplaceBeltrami procedure is applied on the rotated versions of the same initial image and the result is the same but only rotated for each phase value .
Although in the present study we will not provide any results related to image inpainting task of the LaplaceBeltrami procedure, we would like to mention a few related points. The use of the vector field becomes important in texture image inpainting. In that case, on the contrary to the enhancement, we would like to have information flow in orthogonal directions to the object boundaries and reduce the flow along the boundaries. In that case, since also the spatial frequency of the texture patterns have a great importance, we would like to keep the track of the frequency as well as the phase of the evolving output responses, and we would need to fix and to zero instead of and in that case.
4.2 Implementation of the algorithm
4.2.1 The algorithm
We present the steps of our algorithm based on (33) by starting from the initial image function at .

Denote the discrete step in time by . At the iteration (i.e., ) compute the result of the discretized version (of the operator ) applied on the current value of at time instant as and update the solution and the value of by using (33) as follows:

Repeat step 2 until the final time is achieved.

Apply the inverse Gabor transform given by (10) on .
4.2.2 Discrete simple cell output responses
We discretize the image function on a uniform spatial grid as
(36) 
with ( is the number of samples in spatial dimensions) and denoting the pixel width (In general we use square images as input image and we fix in terms of pixel unit). Furthermore the discretized simple cell response of on uniform orientation, frequency and phase grids with points , and (, , (where we denote the number of samples in the orientation dimension by , in the frequency dimension by and in the phase dimension by , and the distances between adjacent samples in the orientation dimension by , in the frequency dimensions by and in the phase dimension by ) is denoted by
(37) 
where and .
In this case the discrete version of the Gabor function given by (8) is written as:
(38) 
where , , . Then we fix (i.e., ) in the reduced framework (which was explained in Section 4.1.1) and write the discrete cell response obtained from the image via the discrete Gabor transform as:
(39) 
The time correspondence in the discrete case is represented by the time index where the time interval is discretized by samples and represents the time instant with satisfying and . In this case the discretized Gabor coefficient is written as
(40) 
4.2.3 Explicit scheme with finite differences
Here we provide the discrete scheme related to the numerical approximation of the algorithm. We propose an explicit finite difference scheme in order to iterate the evolution equation given in (33). The reason for choosing explicit scheme is that implicit scheme requires large memory in our 4dimensional (reduced) anisotropic framework.
We obtain the explicit scheme first by writing (33) in terms of the horizontal vector fields , , and given in (12). Then following Unser [80] and Franken [32], we implement the horizontal vector fields by using central finite differences which are interpolated by Bsplines on a uniform spatial sample grid. Note that Bspline interpolation is required since not all horizontal vectors are aligned with the spatial sample grid.
The interpolation is achieved by determining the coefficients
(41) 
in such a way that the spline polynomial with the Bspline basis functions coincides with the horizontal derivatives of the output at the grid points. For example, in the case of the first horizontal derivative , the condition must hold if we consider a discrete output as explained in Section 4.2.2. We refer to the explanations of Unser [80] for details.
We fix and define
(42) 
See Figure 3 for an illustration of those vectors. We write the central finite differences of the first order horizontal derivatives as
(43) 
and of the second order horizontal derivatives which we use as
(44) 
Then the numerical iteration (discretized from step 2 of the algorithm provided in Section 4.2.1) with a time step is written as follows:
(45) 
where represents the discretized version of given in (32) (with coefficients ) in terms of the central finite differences.
4.2.4 Stability analysis
We must consider two points for the stability of our finite discrete scheme:

Suitable choice of the time step size ,

Preserving the space homogeneity during the LaplaceBeltrami evolution.
The stability analysis for the case is explained in [32] and [23] based on Gershgorin theory. We adapt this technique to our reduced framework and find the upper limit for the time step as:
(46) 
where is the sampling distance between adjacent orientation samples, denotes the number of the orientation samples and is the ratio between orientation and spatial samples. Parameter is either or in our experiments, yielding the condition for stable processes for both values. We refer to [32, Chapter 6] and [17] for details.
The second point is due to that we sample each dimension by using a different number of samples. In order to perform subRiemannian diffusion by regarding the sample unit coherency one must choose the parameters , of the operator in such a way that the space homogeneity of is preserved.
4.3 Experiments
4.3.1 Gabor transform
The delicate point related to the lifting and inversion process is that Gabor functions must be sampled (in orientation , frequency and phase dimensions) in such a way that they cover all the spectral domain (that is, they must fulfill the Plancherel’s formula [66]).
We present some results of the Gabor transforminverse transform procedure associated to our setting and the effects of number of samples in the orientation dimension in Figure 5. We use the Gabor filter banks obtained from (6) and (8) with scale value of 2 pixels (total filter size is 24 pixels) in order to lift the test images (see Figure 5 for some examples of those Gabor functions). On the top row, we see the results related to an artificial test image (left), and at the bottom we see the results related to a real test image (left) taken from Kimmel et. al. [47] We see in the middle and right columns those two images now transformed and then inverse transformed with different number of orientation samples. We sample the space at frequencies , orientations (middle), (right) and phases . We observe that the decrease in the number of orientation samples reduce the quality of the tansformation procedure noticeably in both test images.
4.3.2 Enhancement
The lifting procedure is performed by the Gabor filters of the type given by (6) and (8) with pixels (the filter size is pixels) and time step in the experiments.
In Figure 6, we see the results by of the enhancement procedure applied on an artificially produced gray scale test image with white noise. The lifting is achieved with frequency samples , phase samples and orientation samples . Note that , thus . In order to fulfill physical unit coherency we choose and . The experiments are done with 15 and 30 iterations.
We continue with Figure 7 where we apply our procedure on a real image taken from Kimmel et. al. [47]. In [47] they use a multiscale LaplaceBeltrami procedure with a fixed frequency. We use the same phase and orientation samples as in the case of Figure 6 while we employ the frequency samples for the lifting. Here the coefficients , are chosen as in the case of Figure 6. We perform the experiments with 30 and 50 iterations.
We show in Figure 8, the results related to our LaplaceBeltrami procedure applied on another real image, with dimensions , taken from Kimmel et. al. [47]. We use the same sampling parameters as in the previous case of Figure 7 for the lifting. We perform our LaplaceBeltrami procedure with 6 and 15 iterations. The results are presented together with the multiscale LaplaceBeltrami results obtained by Kimmel et. al [47] for a comparison. Our algorithm takes advantage of different frequencies present in images and therefore can preserve texture structures in such images as in Figure 8. Compare the elongated structures towards the right edge of the images corresponding to Kimmel et. al. [47] (middle right) and to our procedure (bottom right).
5 Conclusion
In this paper we have shown that the multifeature selective simple cells and the associated V1 functional geometry can be modeled starting from a suitably chosen receptive profile, which was the extended Gabor function in our framework. We have derived the whole model subRiemannian geometry and the corresponding horizontal connectivity directly from the receptive profile. In addition to this construction of the model, we have also provided an image processing application employing our model framework: image enhancement via a subRiemannian LaplaceBeltrami procedure. We have provided the algorithm and its discretization explicitly as well as some experimental results. We have also mentioned that in fact, the enhancement procedure could be switched to an image inpainting procedure via a modification of the reduced metric used for the enhancement.
Funding
G. Citti and A. Sarti are funded by the European Union’s Horizon 2020 research and innovation programme under the Marie SkłodowskaCurie grant agreement GHAIA, No 777822.
References
 [1] (2014) A corticalinspired geometry for contour perception and motion integration. Journal of Mathematical Imaging and Vision 49 (3), pp. 511–529. Cited by: §1, §1, §2.1.1, §2.1.1.
 [2] (2011) Coherent states of the Euclidean group and activation regions of primary visual cortex. arXiv preprint arXiv:1111.0669. Cited by: §1.
 [3] (2012) An uncertainty principle underlying the functional architecture of V1. Journal of PhysiologyParis 106 (5), pp. 183–193. Cited by: §1.
 [4] (2015) Reproducing kernel Hilbert spaces of CR functions for the Euclidean motion group. Analysis and Applications 13 (03), pp. 331–346. Cited by: §1.
 [5] (2018) A geometric model of multiscale orientation preference maps via Gabor functions. Journal of Mathematical Imaging and Vision 60 (6), pp. 900–912. Cited by: §1.
 [6] (2019) Uniqueness of viscosity mean curvature flow solution in two subriemannian structures. SIAM Journal on Mathematical Analysis 51 (3), pp. 2633–2659. Cited by: §4.1.
 [7] (2018) Rototranslation covariant convolutional networks for medical image analysis. In International Conference on Medical Image Computing and ComputerAssisted Intervention, pp. 440–448. Cited by: §1.
 [8] (1969) On the existence of neurones in the human visual system selectively sensitive to the orientation and size of retinal images. The Journal of Physiology 203 (1), pp. 237. Cited by: §1.
 [9] (2001) Geometric visual hallucinations, Euclidean symmetry and the functional architecture of striate cortex. Philosophical Transactions of the Royal Society B: Biological Sciences 356 (1407), pp. 299–330. Cited by: §1.
 [10] (2003) The functional geometry of local and horizontal connections in a model of V1. Journal of PhysiologyParis 97 (2), pp. 221–236. Cited by: §1.
 [11] (2013) Myocardial deformation from local frequency estimation in tagging MRI. In International Conference on Functional Imaging and Modeling of the Heart, pp. 284–291. Cited by: §1.
 [12] (2002) Über systeme von linearen partiellen differentialgleichungen erster ordnung. In The Collected Papers Of WeiLiang Chow, pp. 47–54. Cited by: §2.2.
 [13] (2006) A cortical based model of perceptual completion in the rototranslation space. Journal of Mathematical Imaging and Vision 24 (3), pp. 307–326. Cited by: §1, §1, §1, §1, §1, §2.1.1, §2.1.3, §3, §3, §4.1.
 [14] (2014) Neuromathematics of vision. Springer. Cited by: §1.
 [15] (2015) Cortical spatiotemporal dimensionality reduction for visual grouping. Neural Computation. Cited by: §1, §1.
 [16] (2012) Spatiotemporal receptive fields of cells in V1 are optimally shaped for stimulus velocity estimation. JOSA A 29 (1), pp. 130–138. Cited by: §2.1.1.
 [17] (2011) Numerical schemes for linear and nonlinear enhancement of DWMRI. In International Conference on Scale Space and Variational Methods in Computer Vision, pp. 14–25. Cited by: §4.2.4.
 [18] (1985) Uncertainty relation for resolution in space, spatial frequency, and orientation optimized by twodimensional visual cortical filters. JOSA A 2 (7), pp. 1160–1169. Cited by: §1.
 [19] (1983) Spatialfrequencyspecific inhibition in cat striate cortex cells.. The Journal of Physiology 336 (1), pp. 359–376. Cited by: §1.
 [20] (1993) Spatiotemporal organization of simplecell receptive fields in the cat’s striate cortex. i. general characteristics and postnatal development. Journal of neurophysiology 69 (4), pp. 1091–1117. Cited by: §2.1.1.
 [21] (2009) Line enhancement and completion via linear left invariant scale spaces on SE(2). In International Conference on Scale Space and Variational Methods in Computer Vision, pp. 795–807. Cited by: §1.
 [22] (2010) Leftinvariant parabolic evolutions on se (2) and contour enhancement via invertible orientation scores part ii: nonlinear leftinvariant diffusions on invertible orientation scores. Quarterly of Applied Mathematics, pp. 293–331. Cited by: §1.
 [23] (2010) Leftinvariant parabolic evolutions on SE(2) and contour enhancement via invertible orientation scores part i: linear leftinvariant diffusion equations on SE(2). Quarterly of Applied Mathematics, pp. 255–292. Cited by: §1, §4.2.4.
 [24] (2013) Evolution equations on Gabor transforms and their applications. Applied and Computational Harmonic Analysis 35 (3), pp. 483–526. Cited by: §1, §4.1.1.
 [25] (2005) Perceptual organization in image analysis: a mathematical approach based on scale, orientation and curvature. Vol. 68. Cited by: §1.
 [26] (1993) Threedimensional computer vision: a geometric viewpoint. MIT press. Cited by: §1.
 [27] (1993) Contour integration by the human visual system: evidence for a local \csq@thequote@oinit\csq@thequote@oopenassociation field\csq@thequote@oclose. Vision Research 33 (2), pp. 173–193. Cited by: §1, §3.
 [28] (1997) Image structure. Series in Mathematical Imaging and Vision. Springer Berlin. Cited by: §1.
 [29] (2018) A neuromathematical model for geometrical optical illusions. Journal of Mathematical Imaging and Vision 60 (1), pp. 94–108. Cited by: §4.1.
 [30] (2007) Nonlinear diffusion on the 2D Euclidean motion group. Scale Space and Variational Methods in Computer Vision, pp. 461–472. Cited by: Figure 3.
 [31] (2009) Crossingpreserving coherenceenhancing diffusion on invertible orientation scores. International Journal of Computer Vision 85 (3), pp. 253. Cited by: §1.
 [32] (2008) Enhancement of crossing elongated structures in images. Eindhoven University of Technology. Eindhoven, The Netherlands. Cited by: §4.2.3, §4.2.4.
 [33] (2014) Crossingpreserving multiscale vesselness. In International Conference on Medical Image Computing and ComputerAssisted Intervention, pp. 603–610. Cited by: §1.
 [34] (1970) Higher visual perception as prolongation of the basic Lie transformation group. Mathematical Biosciences 6, pp. 437–471. Cited by: §1.
 [35] (1989) The visual cortex is a contact bundle. Applied Mathematics and Computation 32 (23), pp. 137–167. Cited by: §1.
 [36] (1967) Hypoelliptic second order differential equations. Acta Mathematica 119 (1), pp. 147–171. Cited by: §2.2.
 [37] (1963) Shape and arrangement of columns in cat’s striate cortex. The Journal of Physiology 165 (3), pp. 559–568. Cited by: §1.
 [38] (1959) Receptive fields of single neurons in the cat’s striate cortex. The Journal of Physiology 148 (3), pp. 574–591. Cited by: §1.
 [39] (1962) Receptive fields, binocular interaction and functional architecture in the cat’s visual cortex. The Journal of Physiology 160 (1), pp. 106–154. Cited by: §1.
 [40] (1974) Uniformity of monkey striate cortex: a parallel relationship between field size, scatter, and magnification factor. Journal of Comparative Neurology 158 (3), pp. 295–305. Cited by: §1.
 [41] (1977) Ferrier lecture: functional architecture of macaque monkey visual cortex. Proceedings of the Royal Society of London B: Biological Sciences 198 (1130), pp. 1–59. Cited by: §1.
 [42] (1997) Spatial relationships among three columnar systems in cat area 17. The Journal of Neuroscience 17 (23), pp. 9270–9284. Cited by: §1.
 [43] (2008) Models and measurements of functional maps in V1. Journal of Neurophysiology 99 (6), pp. 2745–2754. Cited by: §1.
 [44] (2000) Spatial frequency maps in cat visual cortex. The Journal of Neuroscience 20 (22), pp. 8504–8514. Cited by: §1.
 [45] (1987) An evaluation of the twodimensional Gabor filter model of simple receptive fields in cat striate cortex. Journal of Neurophysiology 58 (6), pp. 1233–1258. Cited by: §1.
 [46] (2013) Direct myocardial strain assessment from frequency estimation in tagging MRI. In International Workshop on Statistical Atlases and Computational Models of the Heart, pp. 212–219. Cited by: §1.
 [47] (2000) Images as embedded maps and minimal surfaces: movies, color, texture, and volumetric medical images. International Journal of Computer Vision 39 (2), pp. 111–129. Cited by: Figure 7, §4.1, §4.3.1, §4.3.2, §4.3.2.
 [48] (2000) On the geometry of texture. Technical report TechnionIsrael Inst. of Tech. Haifa Dept. of Computer Science. Cited by: Figure 8, §4.1.
 [49] (1987) Representation of local geometry in the visual system. Biological cybernetics 55 (6), pp. 367–375. Cited by: §1, §1.
 [50] (1984) The structure of images. Biological cybernetics 50 (5), pp. 363–370. Cited by: §1, §1.
 [51] (2013) Principles of gestalt psychology. Routledge. Cited by: §1.
 [52] (1970) Gestalt psychology: an introduction to new concepts in modern psychology. WW Norton & Company. Cited by: §1.
 [53] (1991) Computational models of visual processing. MIT press. Cited by: §1.
 [54] (1978) Ocular dominance columns and their development in layer iv of the cat’s visual cortex: a quantitative study. Journal of Comparative Neurology 179 (1), pp. 223–244. Cited by: §1.
 [55] (1990) Spatiotemporal interactions and the spatial phase preferences of visual neurons. Experimental Brain Research 80 (2), pp. 441–445. Cited by: §1.
 [56] (1994) Scalespace theory: a basic tool for analyzing structures at different scales. Journal of Applied Statistics 21 (12), pp. 225–270. Cited by: §1.
 [57] (1998) Feature detection with automatic scale selection. International Journal of Computer Vision 30 (2), pp. 79–116. Cited by: §1, §1.
 [58] (2011) Generalized Gaussian scalespace axiomatics comprising linear scalespace, affine scalespace and spatiotemporal scalespace. Journal of Mathematical Imaging and Vision 40 (1), pp. 36–81. Cited by: §1.
 [59] (2013) A computational theory of visual receptive fields. Biological Cybernetics 107 (6), pp. 589–635. Cited by: §1, §1.
 [60] (1977) Spatial frequency rows in the striate visual cortex. Vision Research 17 (2), pp. 257–264. Cited by: §1.
 [61] (1980) Mathematical description of the responses of simple cortical cells. JOSA 70 (11), pp. 1297–1300. Cited by: §1.
 [62] (2002) Detection and discrimination of relative spatial phase by V1 neurons. Journal of Neuroscience 22 (14), pp. 6129–6157. Cited by: §1.
 [63] (1999) Vers une neurogéométrie. fibrations corticales, structures de contact et contours subjectifs modaux. Mathématiques informatique et sciences humaines (145), pp. 5–102. Cited by: §1, §1.
 [64] (2003) The neurogeometry of pinwheels as a subRiemannian contact structure. Journal of PhysiologyParis 97 (2), pp. 265–309. Cited by: §1.
 [65] (2008) Neurogéométrie de la vision. Modeles mathématiques et physiques des architectures fonctionelles. Paris: Éd. École Polytech. Cited by: §1.
 [66] (1910) Contribution à l’étude de la représentation d’une fonction arbitraire par des intégrales définies. Rendiconti del Circolo Matematico di Palermo (18841940) 30 (1), pp. 289–335. Cited by: §4.3.1.
 [67] (1988) Responses of simple and complex cells to compound sinewave gratings. Vision Research 28 (1), pp. 25–39. Cited by: §1.
 [68] (2015) Image processing in the semidiscrete group of rototranslations. In International Conference on Geometric Science of Information, pp. 627–634. Cited by: §1.
 [69] (2013) Organization and origin of spatial frequency maps in cat visual cortex. The Journal of Neuroscience 33 (33), pp. 13326–13343. Cited by: §1.
 [70] (2016) Pinwheeldipole configuration in cat early visual cortex. NeuroImage 128, pp. 63–73. Cited by: §1.
 [71] (2010) A model of natural image edge cooccurrence in the rototranslation group. Journal of Vision 10 (14), pp. 37–37. Cited by: §1.
 [72] (2008) The symplectic structure of the primary visual cortex. Biological Cybernetics 98 (1), pp. 33–48. Cited by: §1, §1, §1, §2.1.3.
 [73] (2009) Functional geometry of the horizontal connectivity in the primary visual cortex. Journal of PhysiologyParis 103 (1), pp. 37–45. Cited by: §2.1.3.
 [74] (2015) Leftinvariant evolutions of wavelet transforms on the similitude group. Applied and Computational Harmonic Analysis 39 (1), pp. 110–137. Cited by: §1, §1, §1.
 [75] (1978) Ocular dominance in layer iv of the cat’s visual cortex and the effects of monocular deprivation.. The Journal of Physiology 281 (1), pp. 267–283. Cited by: §1.
 [76] (2004) The organization of orientation and spatial frequency in primary visual cortex. Proceedings of the National Academy of Sciences 101 (48), pp. 16941–16946. Cited by: §1.
 [77] (2012) Parallel development of orientation maps and spatial frequency selectivity in cat visual cortex. European Journal of Neuroscience 35 (1), pp. 44–55. Cited by: §1.
 [78] (2010) Multiscale and multiorientation medical image analysis. In Biomedical Image Processing, pp. 177–196. Cited by: §1.
 [79] (2003) Frontend vision and multiscale image analysis multiscale computer vision theory and applications. Mathematica. Cited by: §1.
 [80] (1999) Splines: a perfect fit for signal and image processing. IEEE Signal Processing Magazine 16 (6), pp. 22–38. Cited by: §4.2.3, §4.2.3.
 [81] (1938) Laws of organization in perceptual forms.. Cited by: §1.
 [82] (1987) The gaussian derivative model for spatial vision. i retinal mechanisms. Spatial vision 2 (4), pp. 273–293. Cited by: §1, §1.
 [83] (2016) Robust retinal vessel segmentation via locally adaptive derivative frames in orientation scores. IEEE Transactions on Medical Imaging 35 (12), pp. 2631–2644. Cited by: §1.