A cortical-inspired geometry for contour perception and motion integration
In this paper we develop a geometrical model of functional architecture for the processing of spatio-temporal visual stimuli. The model arises from the properties of the receptive field linear dynamics of orientation and speed-selective cells in the visual cortex, that can be embedded in the definition of a geometry where the connectivity between points is driven by the contact structure of a 5D manifold. Then, we compute the stochastic kernels that are the approximations of two Fokker Planck operators associated to the geometry, and implement them as facilitation patterns within a neural population activity model, in order to reproduce some psychophysiological findings about the perception of contours in motion and trajectories of points found in the literature.
keywords: Visual cortex Lie groups Contact geometry Galilean group Cognitive neuroscience Spatio-temporal models
The modular structure of the mammalian visual cortex has been discovered in the seventhies by the pioneeristic work of Hubel and Wiesel [31, 30]. Many families of cells contribute to sampling and coding the stimulus image. Each family is sensible to a specific feature of the image: position, orientation, scale, color, curvature, velocity, stereo.
The main behaviour of simple cells is that of detection of positions and local orientations via linear filtering of the stimulus, and the linear filter associated to a given cell is called its receptive profile. Daugman  proposed a criterion of minimal uncertainty for the shape receptive profile, which results in two dimensional spatial Gabor filters, that was later confirmed by Jones and Palmer  and Ringach . Simple cells also possess a temporal behaviour, studied by De Angelis et al. , and spatio-temporal receptive profiles could again be interpreted by Cocci et al.  in terms of minimal uncertainty, resulting in three dimensional Gabor filters. However, not all cell activity in response to stumuli can be justified in terms of linear filters, and nonlinearities can sometimes become relevant to model their behaviours .
Cells are spatially organized in such a way that for every point (x,y) of the retinal plane there is an entire set of cells, each one sensitive to a particular instance of the considered feature, giving rise to the so-called hypercolumnar organization. Hypercolumnar organization and neural connectivity between hypercolumns constitute the functional architecture of the visual cortex, that is the cortical structure underlying the processing of visual stimulus.
The mathematical modelling of the functional architecture of the visual cortex in terms of differential geometry was introduced with the seminal works of Koenderink  and Hoffmann [27, 28]. While the first author pointed out the differential action of perceptual mechanisms, in particular with respect to jet spaces arising from linear filters, the second author proposed to model the hypercolumnar organization in terms of a fiber bundle structure and pointed out the central role of symmetries in perception expressing them in terms of Lie groups and Lie algebras.
Problems of perceptions can also be afforded with a purely psychophysical approach. The study of Field Hayes and Hess  introduced the notion of association field, as path of information integration along images that can quantitatively satisfy the assumptions of the Gestalt principle of good continuation. The perceptual role of this mechanism was indeed that of contour integration, that typically occurs along field lines associated to locally coherent directions in images.
Almost simultaneously Mumford  proposed a variational approach to describe smooth edges, in terms of the elastica functional, that could be implemented with stochastic processes defining curves with random curvature at any point. Indeed, they produce probability distributions in the space of positions and orientations whose probability peaks follow elastica curves. Williams and Jacobs  could use such stochastic processes to implement a mechanism of stochastic completion, and interpreted the probability kernel they obtained as tensors representing geometric connections on the space of positions and orientations associated to the neural representation of images due to simple cells.
Many of such results dealing with differential geometry were given a unified framework under the new name of neurogeometry by Petitot and Tondud , who related the association fields of Field Hayess and Hess with the contact geometry introduced by Hoffmann and the elastica of Mumford.
The problem of edge organization in images was then addressed in terms of a stochastic process of the type of Mumford, introducing nonlinearities in order to take into account the role of curvature, by August and Zucker [5, 6], while the variational approach of Nitzberg, Mumford and Shiota was extended from edges to level lines of images by Ambriosio and Masnou .
Then, in , Citti and Sarti showed how the functional architecture could be described in terms of Lie groups structures. They interpreted the geometric action of receptive profiles as a lifting of level lines into the space of positions and orientations, and addressed the problem of occlusion with a nonlinear diffusion-concentration process in such a space of liftings. In particular, this approach allows to introduce orientation, instead of depth, as a third dimension for the disentanglement of crossing level lines. Moreover, by making use of the sub-Riemannian structure associated to the Lie symmetries of the Euclidean motion group they could prove that their geometrical nonlinear diffusion approximates the variational result of Ambrosio and Masnou. In their model, then, contour completion is justified as a propagation in the sub-Riemannian setting, and the integral curves of the vector fields that generate the Lie algebra can be considered as a mathematical representation of the association fields of Field, Hayes and Hess, hence proving the relation between neural mechanisms and image completion. This method was then concretely implemented in .
The problem of boundary completion was also addressed from a slightly different point of view by Zucker , who showed the role of Frenet frames.
Exact solutions to the Fokker-Planck equation associated to Mumford stochastic process were provided by Duits and van Almsick , and later Duits and Franken [15, 16] unified such stochastic approach with nonlinear mechanisms of the type of August and Zucker, keeping left invariance with respect to the Lie symmetry of the Euclidean motion group and yet allowing the invertibility of the whole process. Their result was applied to the problem of contour enhancement and contour completion, working on the whole Lie group by means of a representation via suitably defined linear filters. This approach was then extended to different geometric setting by Duits and Führ , again with applications to the processing of images.
We also mention the works of Hladky and Pauls  and Boscain et al , where the authors provided technical developments on the structure of the sub-Riemannian paths and minimal surfaces involved in the cited works.
In this paper we propose a mathematical model of cortical functional architecture for the processing of spatio-temporal visual information, that is compatible with both phenomenological experiments and neurophysiological findings. The adopted theoretical framework follows the outlined path of mathematical models of the activity of the visual cortex, and in particular it continues the geometric approach of Citti and Sarti  to a space of higher dimension, since it takes into account time and velocity of stimuli. We will extend the stochastic process of Mumford  to this setting, working in a space of liftings arising from the filtering with spatio-temporal receptive profiles, and make use of assumptions like the ones of Ermentrout and Cowan  in order to construct a population dynamics able to provide a new form of association fields adapted to the problem of motion integration and motion completion under occlusion. Moreover, the resulting kernels will be comparable to measured neural activities in the presence of stimuli characterized by their direction of motion.
Our starting point in Section 2 is a process of detection resulting from linear filtering with three dimensional Gabor functions with two spatial and one temporal dimensions, which have been proposed as a model of spatio-temporal receptive profiles of primary visual cortex simple cells (see , ). Spatio-temporal Gabor filters extend simple cells Gabor behavior as spatial filters [13, 51], that proved its usefulness for the classical task of edge detection, while its role for motion detection was already pointed out in . The Gabor transform takes in input a moving image , where are spatial coordinates and is time and provides in output a representation of the signal in the phase space
with and representing spatio-temporal position and representing, respectively, spatial and temporal frequency.
Since we are mainly interested in spatio-temporal dynamical aspects, we will assign to temporal frequency the meaning of velocity on the coordinates. We will also select the subset of all detected features corresponding to a fixed value of . This will end up to be a 5D manifold, with a contact structure, induced by a normalization of the Liuoville form.
In Section 3 we will outline that this constraint carries a notion of admissible curves  in a deterministic and a stochastic setting, allowing to compute the kernels connecting filters in the 5D manifold. The stochastic processes are completely determined by the described structures of the tangent space of the 5D manifold and they turn out to be described by the fundamental solution of a Fokker Planck equation . In fact the processes contain diffusion in the fiber variables and transport along the remaining generators of admissible tangent directions, in the spirit of [40, 69]. We will compute these kernels in two limit cases: the motion of a contour at a fixed time instant, and the motion of a point moving in time. The first one reduces to the geometry of contours, with a notion of instantaneous velocity, the second one corresponds to point trajectories in time, and can be related to the subset of the Galilei group  on the plane.
In Section 4 we will discuss the compatibility of the previously calculated kernels with psychophysiological findings reported in the recent literature. Then we will insert the connectivity kernels computed in Section 3 in a neural population activity model , regarding them as cortical facilitation patterns.
In section 5 we will use this population activity model equipped with the suitable connectivity kernels in two numerical simulations, comparing the results to recent phenomenological findings of the perception of contours in motion , and to fMRI measurements of cortical neural activity related to motion perception .
2 The geometry of spatio-temporal dynamics
In this section we extend an approach introduced in  that amounts to model each V1 simple cell in terms of its receptive profile, to interpret its action as a Gabor filtering, and to introduce a geometry of the space compatible with the properties of the output. In this paper we will consider each cell as sensitive to a local orientation and apparent velocity, that is the velocity orthogonal to a moving stimulus. The collected data lies on a five dimensional manifold of space, time, orientation and velocity, a single cell being represented as a point on . The geometric structure of this manifold will be described in terms of a contact structure, that provides a constraint on the tangent space and on admissible connectivity among cells.
2.1 Spatio-temporal receptive profiles
It is known that the visual cortex decomposes the visual stimulus by measuring its local features. Local orientation and direction of movement have been the first visual features of neurons in V1 that have been studied. Receptive profiles (RPs) are descriptors of the linear filtering behavior of a cell and they can be reconstructed by processing electrophysiological recordings . It has been shown that the spatial characteristics of these RPs can be modeled by 2-dimensional Gabor functions [13, 32].
However a large class of cells shows a very specific space-time behavior in which the spatial phase of the RP changes gradually as a function of time . Although many models have been proposed to reproduce these dynamics [14, 1], in  it has been shown that a 3-dimensional sum-of-inseparable Gabor model can fit very well experimental data of both separable and inseparable RPs (Fig. 1). Following this approach, we choose Gaussian Gabor filters centered at position on the image plane, activated around time , with spatial frequency , temporal frequency , spatial width (circular gaussians) and temporal width
where we have used the abbreviations and . Also note that the functions (1) correspond to the propagation of a two dimensional plane wave within the activating window with phase velocity
The variable is canonically associated  to the phase space
endowed with the symplectic structure compatible with the complex structure of of variable , that is
where we have denoted with the Liouville form
Since we are mainly interested in the geometry of level lines, our analysis can be restricted to the information captured by filters possessing a fixed central spatial frequency . This amounts to fixing a scale of oscillations, hence disregarding the harmonic content of the filtering by focusing on the features of orientations and velocity. We then obtain the reduced 1-form
where . This form is defined on the spatio-temporal phase space with fixed frequency, that is the 5-dimensional manifold
and it is associated to every Gabor filter as shown in Fig. 2. At this level, time is introduced as a base variable with the same role of . The dual variable of is , which has the same role of , both being the engrafted variables with respect to time and space. To every point of corresponds univocally a Gabor filter whose parameters are the coordinates of the point itself.
2.2 Admissible tangent space as constraint on the connectivity on
We will model the connectivity between points in the space in terms of admissible tangent directions of itself. From the geometric point of view, the presence of the 1-form (5) is equivalent to the choice of a vector field with the same coefficients as with respect to the canonical basis , that is its dual vector (see Fig. 2):
The kernel of , denoted (or ) is the space of vectors orthogonal to . A basis of this space is constituted of the so-called horizontal or admissible vectors
that defines the horizontal tangent space. In particular we nota that the Euclidean metric on the horizontal planes makes the vector fields orthogonal.
The only manifolds (curves or surfaces) admissible in this space are the ones whose tangent vectors are linear combinations of the horizontal vectors (6). The neural connectivity between the receptive fields at different points in the space will be defined in terms of these vectors in Section 3.
For reader convenience we compute here explicitly all the non-zero commutation relations between the vector fields 555To every vector field we can associate a directional derivative with the same coefficients. Then we will call commutator of and We say that and commute if . Note that partial derivatives always commute, while directional derivatives in general do not. It is important to note that, even though the commutator is expressed formally as a second derivative, it is indeed a first derivative, so that it has an associated vector. Hence we can as well define the commutator between the vectors and as the vector associated to . :
In Fig. 3a we visualize the structure restricted to the spatial and engrafted variables , spanned by and . The tilting of the planes is due to the non commutative relation (8). In Fig 3b we show the spatio-temporal structure restricted to the variables , where . Also in this figure the tilting of the planes spanned by and is due to the non vanishing commutator condition (9). It is worth noting that in this setting, the propagation is permitted only along the horizontal planes, while it is forbidden in their orthogonal direction . Let us note that is linearly independent of the horizontal tangent space at every point. Hence together with their commutators span the whole space at every point. This is the so called Hörmander condition, which will have a crucial role in studying the property of the space. The same condition can be expressed in terms of properties of the form , which is called contact form  since is never zero being the volume form of the space.
We also note that we can define on a smooth composition law
where is a counterclockwise rotation of an angle , that is such that the vector fields are left invariant with respect to this law. This implies that also the admissible curves of the structrue and its kernels will be invariant. The manifold together with this composition law can also be identified with a subset of the Galilei group (see Appendix).
2.3 The output of the receptive profiles
The output of the Gabor filters selects a set of points in corresponding to specific features of the image. Since the Gabor filters are always connected in terms of the vectors (6), also their output will inherit this structure, and will be concentrated around an admissible surface.
The energy output of a cell with receptive profile in presence of a spatio-temporal stimulus is given by
We note that taking the square modulus in (11) disregards the phase of the corresponding linear filtering, which for many applications (see e.g. ) is crucial. In this case however all the relevant information is encoded in this energy model. In particular the geometric quantities and are encoded in the points of maxima of the energy. To see why, we recall that the analogous of formula (11) with purely spatial Gabor filters was studied in , where the filtering output was a function of the variables and alone (see also ). In order to outline the goemetric meaning of the lifting, we consider its action on level lines which are smooth. In that case, the output to a stimulus takes its maximum around a value
where identifies the orientation orthogonal to the level lines of . As a consequence, in  it is proved that the level lines of are lifted to curves admissible in the sense that their tangent vector lies in the kernel of the form
which is a contact form that can be obtained by restricting (5) to the space .
Since here we take into account time, then the output is a function defined on . In perfect analogy with the lower dimensional case, and denoting with a level set of the stimulus , the output takes its maximum around the values and
such that the vector is orthogonal to the level set of . The vector is orthogonal to the spatial level line, and the scalar represents the apparent velocity in this direction.
is the 5D lifting of the level set of , and the orthogonality condition implies that it is admissible, in the sense that its tangent vectors lie in the horizontal space, kernel of the form (see Fig. 4).
3 Curves and kernels of connectivity
The lifted points of the spatio-temporal stimulus are connected by admissible integral curves. However we will see that not all admissible curves can be considered lifted ones, and we will describe which ones have this property. These curves will represent the association fields in space-time, analogously to the association fields of Field Hayes and Hess in the pure spatial case, and are reminiscent of the classical Gestalt concept of good continuation. We will discuss in Section 4 the compatibility of this structure with the cortical functional architecture. A probabilistic version of the connectivity field will be provided on the basis of Fokker Planck equations first introduced by D. Mumford in the spatial case , and interpreted as model of connectivity in Lie groups in several works such as [48, 73, 6, 15, 16, 57]. Here we will extend this stochastic approach to the space-time contact structure.
3.1 Generators of lifted curves
We have already seen that the output of RP filtering is concentrated around admissible submanifolds. Nevertheless not all admissible submanifolds can be lifting of images, since lifting are graphs of the functions and . For example the plane generated by the vectors and cannot be recovered by lifting. Hence we will study the possible linear combinations of vector fields , which can be tangent to lifted level lines. An -admissible curve is characterized by , i.e.
with not necessarily constant. Lifted curves depend on the vectors or , which are tangent to the base space , and their linear combination. To simplify the problem, we will consider separately two special cases of particlar interest for the model: the limit cases of contour motion detected at a fixed time , and motion of a point in time . The first one is described by integral curves of the vector and of the generators of the engrafted variables and :
These curves lie in the section of the contact structure depicted in Figure 3a, where represents Euclidean curvature and is the rate of change of local velocity along the curve. On the other hand, the motion of a point in time can be described as
that are suitable to describe spatio-temporal trajectories of points. The coefficient of the direction is the angular velocity, and the coefficient of is the tangential acceleration.
3.2 Curves and kernels for contours in motion
3.2.1 Moving contours as deterministic integral curves
We consider now the geometry of moving contours at a fixed time. These are generated by the integral curves (12) so to satisfy the system of equations
with variable coefficients , . In particular, due to condition (8) the set of points reacheable by piecewise constant integral curves of this type (for an explicit expression see e.g. ) starting from a fixed point is the space
where is fixed, so the space can be identified with the space of points .
As we did with (10) we can define on a smooth composition law
where is a counterclockwise rotation of an angle , that is such that the vector fields are left invariant with respect to this law. Since this is a group law, the manifold together with this (16) provides the Lie group , that is the direct product of the group and the group of the reals with addition.
3.2.2 Stochastic kernel
Let us consider a probabilistic counter part of trajectories (14). The vector field expresses a derivative in the direction of a variable coded on the retina, while the vectors , express derivatives in the direction of an eingrafted variable. Due to their different role, we will consider the following vector-valued stochastic process
where is a two dimensional Brownian motion. The distribution of these stochastic curves is mostly concentrated around the surface . The density of points reached by this stochastic kernel is then a candidate to implement the mechanism of association fields (see Fig. 6). This approach generalizes the approach of random paths introduced in [40, 69] for the problem.
If we call the density of points of reached at the value of the evolution parameter by the sample paths of the process (17), then can be obtained as the fundamental solution
where is the Fokker-Planck operator
containing a diffusion over the fiber variables and and a drift in the base variables, and with and we denote first and second order derivatives in direction . Since contains a set of vector fields that generates the whole tangent space of , by Hörmander theorem this operator is hypoelliptic . This means that is non-null in all for any , even if the operator contains only 3 linearly independent fields. Moreover, since it defines the Fokker-Planck equation for the stochastic process (17), is indeed a (conditional) probability density on , evolving with the parameter .
In order to characterize each point of in terms of the density of paths (17) that reach it, independently of the value of the evolution parameter, we need a notion corresponding to the fan , expressed in terms of . The density of points reached at any value of the evolution parameter by the stochastic dynamics (17) is given by
This derived quantity is actually the fundamental solution of the operator , so explicitly we have
We recall that in  several representations of the exact solution to a problem analogous to (20) were presented. In this work we will deal with the numerical implementation of the fundamental solution , developed with standard Markov Chain Monte Carlo methods (MCMC) . This is done by generating random paths obtained from numerical solutions of the system (17) and averaging their passages over discrete volume elements, and appears suitable to treat also the more involved case of subsequent equation (26). An example is shown in Fig. 6, where an isosurface plot of the kernel is plotted over the integral curves (14) with constant coefficients, already depicted in Fig. 5. By comparing such numerical approximations with , we can confirm the accuracy of this method. Moreover, from the figure it can be seen that the probability density is concentrated around the surface and decays rapidly away from it. This is reasonable since for this kind of hypoelliptic operators one can generally obtain estimates for the fundamental solution in terms of exponential decay with respect to a geodesic distance computed with respect to the minimal set of vector fields that, together with their commutators, span the whole tangent space . Such a distance is anisotropic, and its balls are squeezed in the directions of the commutators . Here the vector fields involved are , hence the concentration of the fundamental solution around the set defined by their integral curves follows by the sharper decay in the direction of the commutators. We also note that in this particular case, due to the availability of exact solutions and estimates from , this behaviour can also be checked directly.
3.3 Curves and kernels for point trajectories
3.3.1 Point trajectories as deterministic integral curves
The trajectory-type curves introduced in (13) are solutions to the system of ordinary differential equations
for given initial point . In general the coefficients need not to be constant, but we can have a local approximation of any curve in a neighborhood of the starting point if we consider the case of constant coefficients. In this model case when and are not zero, (22) is explicitly solved by
where and , and a reasonable choice is to set , in order to synchronize the evolution parameter with the time parameter .
The fan of such curves is depicted in Fig. 7. Each of them describes a motion on the plane that, for small times, consists approximately of arcs of circles with radius , while for sufficiently large times corresponds to enlarging spirals around a slightly moving center. Their instantaneous acceleration is given by , so that corresponds to the tangential acceleration666By direct computation, is the time derivative of the modulus of the velocity along direction :
while their curvature is
In particular, when we obtain straight lines along the direction, while for we obtain circular trajectories of radius .
3.3.2 Stochastic kernel
Similarly to what we have done in section 3.2.2, let us consider the vector-valued stochastic process
where is a two dimensional brownian motion.
The density of points of reached at the value of the evolution parameter by the sample paths of the process (23), is the fundamental solution of the equation
where is the Fokker-Planck operator
Equation (24) is not hypoelliptic, indeed its fundamental solution is concentrated on a submanifold of codimension 1 defined by the equation , from the system (23). However it is still a Fokker-Planck equation, hence is nonnegative, and integrating the density with respect to the evolution parameter we derive the fundamental solution of the operator . Explicitly we have
and we note in particular that this equation is now a hypoelliptic equation that is also the Fokker-Planck equation of a stochastic process defining a propagation in the physical time . Moreover, since the vector fields that consitute the operator (25) are left invariant with respect to the composition law (10), the solution (19) possesses the symmetry
where stands for the left inverse with respect to (10), see also Appendix.
Similarly to what we have done in the previous section, we can obtain at each point by solving numerically the system of stochastic differential equations (23) and applying standard Markov Chain Monte Carlo methods. An isosurface plot of the kernel is shown in Fig. 8, over the same integral curves of Fig. 7 and again we see the concentration around the fan .
4 Neural propagation of boundaries and trajectories
The curves (14) and (22) can be related to well-defined perceptual mechanisms. More precisely, we propose to consider them as association fields, in the sense of , devoted to the integration of contour in motion and trajectories. The perceptual tasks of subjective boundary completion and motion integration have been widely studied by both psychologists and physiologists, and their underlying physiological explanation continues to be an open issue of discussion. The perceptual bias towards collinear stimuli has classically been associated to the long-range horizontal connections linking cells in V1 sharing similar preferences in stimulus orientation. This specialized form of intra-striate connectivity pattern is found across many species, including cats, tree shrews and macaques [33, 9], the main difference being the specificity and the spatial extent of the connections. Furthermore, axons seem to follow the retinotopic cortical map anisotropically, with the axis of anisotropy being related to the orientation tuning of the originating cell . These connections have already been modeled in  by means of a sub-Riemannian diffusion process over the orientation space .
Similarly to what happens for the integration of merely spatial visual information, the brain is also capable to easily predict stimulus trajectories , and to group together boundary elements sharing similar motion or apparent motion paths [35, 50]. One possible explanation for these effects could be the existence of specialized facilitatory networks linking cells anisotropically and coherently with their preference in velocity and axis of motion direction. Supporting these assumptions, it has been found that the neural preferences in direction of movement are also structurally mapped in the cortical surface, with nearby neurons being tuned for similar motion direction , and it has been shown that excitatory horizontal connections in the V1 of the ferret are strictly iso-direction tuned . Furthermore, it is known that also extra-striate area MT/V5 is retinotopically organized, its horizontal connectivity pattern being highly structured, with connections reaching columns of cells tuned for similar orientation and direction preference anisotropically and asymmetrically [38, 2]. Moreover, striate and extra-striate cortical areas seem to cooperate, and surround modulation in V1 can be given by the connectivity patterns implemented in both areas by means of fast feedforward and feedback inter-areal projections .
Regarding the physiological correlates of motion integration, it was shown that the motion of an occluded object trajectory is significantly represented in the human brain by the same visual areas that process real motion . Moreover, a recent study based on electrophysiological recordings of the V1 of tree shrews showed some non-linear neural behaviors that are coherent with the phenomenological dynamics of motion integration . It is indeed a general assumption that some cortical area in the visual cortex is responsible for predicting future motion, a possible implementation being a specialized connectivity for spatio-temporal trajectory facilitation, that is different from the one in V1 responsible for contour integration [25, 62, 64, 68].
Here, we do not speculate on the exact physiological origin of these psychophysiological findings. However, we want to show that the connectivity kernels and arising from the geometry defined in Section 2 is capable to reproduce qualitatively some of the effects reported on the works that have been previously cited. We will embed this connectivity in the neural population activity model described in the next paragraph. Then, in Section 5 we will simulate the response of cortical visual neurons to artificial stimuli, comparing the results with some psychophysiological findings reported in the recent literature.
4.1 Modeling neural activity
The state of a population of cells can be characterized by a real-valued activity variable, which depends on the interaction of the feedforward input at different points, due to the cortical connectivity. The first population activity models are due to Wilson and Cowan , and Ermentrout and Cowan :
where is the neural activity of the population, is the feedforward input (11), is the facilitation kernel, is the facilitation strength and is the sigmoidal function
In the stationary case a first order approximation of the solution of (28) is
that is the activity formula that we will use in the experiments of Section 5. Let’s explicitly note that the term
represents the mean neural extra-cellular activity in response to a stimulus.
The geometry of the functional architectures is contained in the kernel , taking into account the deep structures of the connectivity space. Then the overall probability of activation can be obtained by convolution of the activity with the kernel , and the term
is the facilitation pattern resulting from the contribution of horizontal cortico-cortical or feedback inter-areal connectivities.
In the model of Ermentrout-Cowan only position and orientation were considered and symmetry properties were imposed to the facilitation patterns. In , for the same features it was proposed to choose as a kernel the fundamental solution of a Fokker Planck equation deduced from the Euclidean symmetries. In the next section, two numerical simulations will be performed using equation (30). In the first one, we will model the propagation of boundaries in motion using the connectivity kernel (19). In the second experiment we will model instead the propagation of trajectories, using the 5D kernel (26).
5 Numerical simulations
In this section some numerical simulations will be performed to test the reliability of the kernels computed in Section 3 when introduced in the activity model (30) and the results will be compared to some psychophysical and physiological experiments.
5.1 The feedforward and extracellular activity in response to a stimulus
For the subsequent numerical simulations, we measured in every discrete point of a stimulus of size , its local energy of orientation and speed by convolving the input image sequence with a pre-determined bank of Gabor filters of fixed spatial frequency centered at the points , thus discretizing Eq. (11) so to have
where and are the component frequencies and .
We chose the spatio-temporal frequencies and the maximum velocity value represented in the Gabor filter set depending on the stimulus, so that given the couple the filters have a maximum temporal frequency of . We normalized the filters so that the response could range from (in regions with no changes in luminance) to (square plane waves sharing the filter parameters and ), corresponding to a normalization of the stimulus contrast. The Gabor scale parameters will follow the relations
whose meaning is to approximately have spatial subregions under the Gabor Gaussian envelope, and a variable that allows the filter with maximum velocity to cover one wavelength within the Gabor active time interval (see [12, 14, 51] for the physiological justifications).
The thresholded feedforward input is then computed following Eq. (29), representing the mean neural extra-cellular activity in response to the stimulus. In the subsequent analyses, the values to which we set the parameters and in the above stated formulas will be explicitly specified. It is worth noting that these parameters are not due to physiological data, due to the difficulty in finding some reference in the literature, but are chosen as to have a computationally reasonable number of non-null measurements .
5.2 Experiment 1 - Contours in motion
For the first simulation we use a dashed circle in linear motion moving upwards within the visual space. The image size is pixels and the image sequence is composed of frames. The circle segments are approximately pixels wide and the circle is moving with a uniform speed of pixels per frame (see Fig. 9). To get a measurement of local spatio-temporal features, we convolve the stimulus with a set of 3-dimensional Gabor functions with (11). The spatial frequency parameter is set to , so that the Gabor moving subregions match the width of the segments. The maximum local velocity represented in the filter set is pixel per frame, while spatial and temporal scale parameters are calculated accordingly with (34). After having convolved the stimulus with the Gabors, we model the neural extra-cellular activity by using (29) with parameters , , obtaining .
We select the 4-dimensional subset corresponding to the nd frame , where , thus neglecting behaviors over time. This is because the model given by the stochastic model of connectivity distribution is stationary, even if it is relying on spatio-temporal information, and is compatible with the dimensionality of in (15).
Following the assumptions that the measurements can model the output of a cortical direction- and speed-selective cell (or of a neural population that are selective to the same visual features), the overall continuation probability can be obtained with a discretized version of the cortico-cortical facilitation pattern (32) with the kernel :
where is calculated for every fiber vector on the same discretized domain of the lifted stimulus, using the stochastic approach described in Section 3.3.2. The structure of (35) is that of a discrete version of a group convolution with respect to the composition law (16), due to the symmetry (21) of the kernel and the left invariance of the Lebesgue measure on .
It is worth noting that the parameters and governing the diffusion process are related to the maximum spatio-temporal curvature of an illusory contour. Due to the lack of reference physiological data, the values for the diffusion coefficients were chosen in such a way that at the final value taken by the evolution parameter of the stochastic paths, the mean square displacements of the fiber variables are , that means and .
Finally, the population activity is computed by
that is a discretized version of (30) where is a coefficient governing the total strength of the excitatory connections. In Fig. 10 are shown the results of the simulation. In the top row we have plotted a iso-level surface of the normalized and thresholded output measurements . The process correctly lifts the stimulus around its theoretical values in the domain . In the bottom row, the iso-level surface of the global cortical activity is shown. Note that continuation is performed within the whole manifold , propagating the local speed and direction cues at the end of the segments, and interpolating the data to define a 4-dimensional set that carries information even in the visual space between the segments, where no changes in contrast were originally present in the stimulus. These results are coherent with recent psychophysiological findings, where it is stated that a global shape in motion is better perceived if local velocity changes smoothly along its contour , and that the integration of the motion of a partially occluded object is facilitated when its visible contours define closed configurations .
5.3 Experiment 2 - Motion integration
The stimuli that we will use in this simulation will be several instances of an object in motion along a certain direction, that disappears at a given time position , and reappears at with a direction of motion changed by in a position that is coherent to a piecewise continuous trajectory (Fig. 11). It is a well documented fact that humans tend to perceive the two trajectories as pertaining to a single unit just for small values of and , an effect that is commonly referred to as motion integration. In particular, we know from many psychophysiological experiments that the chances of detecting a straight or curvilinear trajectory in noise increase with stimulus duration and with direction coherence [66, 63].
To do so, we created multiple instances from the stimulus paradigm described above, and shown in Fig. 11, assigning different values to the parameters and . It is worth noting that we are using an elongated object whose axis of motion is always orthogonal to its major axis of eccentricity. The choice of this particular stimulus is naturally explained by the geometry, as the points that are connected by its continuation property implicitly define a local orientation and a direction of motion that is orthogonal to that orientation. The Gabor filtering of a different kind of stimulus (for example a moving dot) would indeed measure high responses also for those velocity projections that are coherent with the real axis of motion, due to the motion streak effect , requiring a different neural model to detect unambiguously the direction of motion. Thus, even if the majority of the psychophysiological experiments use moving dots as stimuli, the validity of the following experiment is not influenced, as exploring the exact physiological implementation of motion integration is out of the scope of this paper.
The image sequences to process are pixels wide and are composed of frames. The value of eccentricity of the moving ellipsoidal object is set to , and its velocity is set to pixels per frame. Each stimulus instance is uniquely identified by the couple . We convolve each sequence with a set of Gabor filters with and set to match width of the object’s minor axis, and again, we threshold the output according to (29) using parameters , , obtaining the set of measurements . In the previous simulation we have taken a temporal slice of this output, and then we have propagated the activity using the connectivity kernel . Similarly, here we will propagate the activity present in , but without discarding time. The corresponding facilitation pattern is then obtained convolving it with the 5-dimensional kernel defined on the same discrete domain of the lifted stimulus using the same approach that led to (35)
The results obtained by processing two stimulus configurations are shown in Fig. 12, where only the central portion of the full temporal domain is plotted, to highlight the different effects that activity propagation has in the two cases. When the change in direction of motion of the object is under a certain threshold, the trajectory is completed smoothly, yielding a strong activation in the temporal interval when the object is not present in the scene. If the value of increases, however, the spatio-temporal curvature of the optimal subjective trajectory progressively becomes too high to be completed by the connectivity, even if a weak response linking the two parts of the stimulus could yet be present.
It is worth noting that due to the thresholding stage, this subjective interpolation effect is strongly non-linear, as it cannot be explained by the sum of the responses to the first () and second () part of the stimulus alone (see Fig. 13). This is coherent with the findings in the work of Wu et al, where electrophysiological experiments recorded a similar non-linear behavior of trajectory interpolation for small values of . In that work, the cause of this effect is left unexplained. Even if the stimulus paradigms are slightly different (they use a field of moving dots with a null value of ), the result of our simulations allow us to make some speculations.
A good qualitative description of the effects that the modeled excitatory connectivity has on stimulus response can be found in Fig. 14, where for every stimulus configuration we plot the difference
where is the full stimulus and and represent the first () and the second () part of the stimulus. The visualization of the output highlights the role that a trajectory-specialized cortical connectivity could have in performing tasks of motion integration. The spatio-temporal extension of over gets smaller for higher values of and . The decay of the facilitation effect is coherent with the observations made in  (where ) and  (where ), even if the experiments are methodologically different. Moreover, even if, coherently with our results, some psychophysiological experiments showed that a broken trajectory in noise is easily detectable [65, 59], as far as we know little has been done to explore the effect of changing the duration of the temporal occluding gap.
Regarding the implementation of our model, we would like to highlight functional meaning of the parameter couple , driving the diffusion along the fiber variables when calculating the kernels and . In particular, the parameter seems to be strongly related to the maximum perceived curvature of illusory contours and trajectories, while sets the maximum rate of change of local velocity along admissible subjective contours (trajectories). It would be interesting to try to fine-tune the parameters of the model in order to reproduce quantitatively as precisely as possible the psychophysiological findings found in literature. We aim to carry out this analysis in a future paper.
It is worth noting that the physiological correlates of the first simulation are not well documented and, as far as we know, no one has yet studied the neural activity in cortical regions responding to subjective contours in motion. We suppose that it would be a difficult issue to address, as multiple cortical layers, as V1, V4 and MT/V5, may be involved . Recently, some significative results have been obtained with static illusory contours, using promising electrophysiological techniques . Those kinds of experiments, targeted to the detection of illusory motion contours, could give additional clues about the neural computation that governs the influence of spatio-temporal features in the detection of moving shapes and boundaries.
As discussed in the beginning of Section 4, a possible physiological implementation of the geometry used in the second numerical simulation, could be a trajectory-specialized cortico-cortical connectivity between neurons in the higher visual areas. The advantages of having such an excitatory connectivity implemented in the visual cortex could be very important in performing complex cognitive tasks. For example, we assume that the brain could use the improvements in contour detection with respect to the only related functional architecture, due to an additional information, that of velocity, that indeed in practical situation is a feature that is coherent on objects. The influence that the connectivity arising from our geometrical model has in the processing of visual tasks such as spatio-temporal grouping or segmentation will be the subject of a future paper.
In this paper we proposed a model of cortical functional architecture for the processing of spatio-temporal visual information. Motion features are detected first by simple and complex cells with RP modelled by 2D+time Gabor filters. Linear filtering with such a profiles lifts the visual stimulus from the base space to the phase space comprising spatial and temporal frequencies, in which a Liouville form can be defined.
We defined a 5D phase space with fixed frequency by taking a reduction of the previous differential form. Then, exploiting the commutation properties of the horizontal basis, we regarded the tangent space of the contact manifold as a possible constraint acting on the connectivity between points, giving the definition of admissible integral, or horizontal, curves. Possible linear combinations of the horizontal basis that are compatible with the definition of admissible curve have been studied, in order to model the possible lifting of the visual stimuli as association fields, in the sense of . We considered the corresponding deterministic integral curves for two modeling limit cases: contours in motion and trajectories of a point in motion.
Then we considered horizontal stochastic paths, i.e. trajectories of points that always move along the tangent space of and are allowed to change the value associated to the fiber variables in a random, equidistributed way. We have seen that the resulting integrations over the evolution parameter of the paths coincide with the kernels of the Fokker Planck operators defined on the geometries. Both deterministic curves and stochastic kernels inherit the symmetries of the non associative reduced Galilean group as described in Appendix.
After a discussion about the compatibility of our theoretical framework with various phenomenological and psychophysiological findings on visual perception and cognition reported in the literature, we introduced the stochastic kernels as facilitation inducers in a neural population activity model.
Suitable numerical simulations have been carried out, by processing pre-determined artificial stimuli with the neural population activity model previously described, showing that the Fokker Planck stochastic kernels endow the model with the capability of completion and continuation of contours in motion and trajectories of points, coherently with the phenomenological experiments of Rainville  and Wu et al .
In conclusion, results have shown that the proposed functional geometry is compatible with existing psychophysical and physiological experiments, even if a complete knowledge about the effective neural implementation needs supplementary empirical data.
Appendix: Relations with the Galileian group
The geometry of can be related with the classical Galilei group of motions on the plane, when choosing an appropriate local frame of reference.
The two dimensional Galilei group is the semidirect product Lie group defined by the composition law 
where provides space-time translations, is a counterclockwise planar rotation of an angle and is the velocity vector responsible of the purely Galilean transformation, or boost, expressed in Cartesian coordinates.
The manifold can be identified with a subset of this group by identifying the direction of velocity with the angle of the rotation. Setting the velocity in polar coordinates
then we can identify with the set of points
that is, using the identification . Since the angle is interpreted as the local direction of a boundary, then this amounts to select the apparent velocity orthogonal to it.
This manifold is not invariant under the Galilei composition, but it is invariant with respect to the smooth composition law (10). This composition law is not associative, but possesses a trivial neutral element and each possesses a left inverse element, i.e. such that
and for these reasons we can consider it a reduced Galilei non associative group .
-  Adelson, E.H., Bergen, J.R.: Spatiotemporal energy models for the perception of motion. J. Opt. Soc. Am. A 2(2), 284-299 (1985)
-  Ahmed, B., Cordery, P.M., McLelland, D., Bair, W., Krug, K.: Long-range clustered connections within extrastriate visual area V5/MT of the rhesus macaque. Cereb. Cortex 22(1), 60-73 (2012)
-  Ambrosio, L., Masnou, S.: A direct variational approach to a problem arising in image recostruction. Interfaces and Free Boundaries 5 (1), 63-81 (2003)
-  Angelucci, A., Levitt, J.B., Walton, E.J.S., Hupe, J., Bullier, J., Lund, J.S.: Circuits for local and global signal integration in primary visual cortex. J. Neurosci. 22(19), 8633-46 (2002)
-  August, J., Zucker, S. W.: The curve indicator random field: curve organization via edge correlation. In: Boyer, K., Sarkar, S. (eds.) Perceptual organization for artificial vision systems, pp. 265-288. Kluwer (2000)
-  August, J., Zucker, S. W.: Sketches with curvature: the curve indicator random field and Markov processes. IEEE Trans. Pattern Analysis and Machine Intelligence 25(4), 387-400 (2003)
-  Boscain, U., Duits, R., Rossi, F., Sachkov, Y: Curve cuspless reconstruction via sub-Riemannian geometry
-  Bosking, W.H., Zhang, Y., Schofield, B., Fitzpatrick, D.: Orientation selectivity and the arrangement of horizontal connections in tree shrew striate cortex. J. Neurosci. 17(6), 2112-27 (1997)
-  Chisum, H.J., Mooser, F., Fitzpatrick, D.: Emergent properties of layer 2/3 neurons reflect the collinear arrangement of horizontal connections in tree shrew visual cortex. J. Neurosci. 23(7), 2947-60 (2003)
-  Citti, G., Manfredini, M., Sarti, A.: Neuronal oscillations in the visual cortex: convergence to the Riemannian Mumford-Shah functional. SIAM Journal of Mathematical Analysis 35 (6), 1394-1419 (2004).
-  Citti, G., Sarti, A.: A cortical based model of perceptual completion in the roto-translation space. J. Math. Imag. Vis. 24(3), 307-326 (2006)
-  Cocci, G., Barbieri, D., Sarti, A., Spatio-temporal receptive fields of cells in V1 are optimally shaped for stimulus velocity estimation. J. Opt. Soc. Am. A 29(1), 130-138 (2012)
-  Daugman, J., Uncertainty relation for resolution in space, spatial frequency, and orientation optimized by two-dimensional visual cortical filters. J. Opt. Soc. Am. A 2(7), 1160-1169 (1985)
-  DeAngelis, G.C., Ohzawa, I., Freeman, R.D.: Spatiotemporal organization of simple-cell receptive fields in the cat’s striate cortex. I. General characteristics and postnatal development. J. Neurophysiol. 69(4), 1091-1117 (1993)
-  Duits, R., Franken, E.M.: Left invariant parabolic evolution equations on SE(2) and contour enhancement via invertible orientation scores, part I: Linear left-invariant diffusion equations on SE(2). Q. Appl. Math. 68, 255-292 (2010)
-  Duits, R., Franken, E.M.: Left invariant parabolic evolution equations on SE(2) and contour enhancement via invertible orientation scores, part II: Nonlinear left-invariant diffusion equations on invertible orientation scores. Q. Appl. Math. 68, 293-331 (2010)
-  Duits, R., Führ, H., Janssen, B., Bruurmijn, M., Florack, L., Van Assen, H.: Evolution equations on Gabor transforms and their applications. arXiv:1110.6087 (2011)
-  Duits, R., van Almsick, M.: The explicit solutions of linear left-invariant second order stochastic evolution equations on the 2D Euclidean motion group. Quart. Appl. Math. 66, 27-67 (2008)
-  Ermentrout, G.B., Cowan, J.D.: Large scale spatially organized activity in neural nets. SIAM J. Appl. Math. 38(1), 1-21 (1980)
-  Field, D.J., Hayes, A., Hess, R.F.: Contour integration by the human visual system: evidence for a local association field. Vision Res. 33, 173-193 (1993)
-  Folland, G.B.: Harmonic analysis on phase space. Princeton University Press (1989)
-  Geiges, H.: An introduction to contact topology. Cambridge University Press (2008)
-  Geisler, W.S.: Motion streaks provide a spatial code for motion direction. Nature 400(6739), 65-69 (1999)
-  Graham, N.V.: Beyond multiple pattern analyzers modeled as linear filters (as classical V1 simple cells): Useful additions of the last 25 years. Vision Research 51, 1397-1430 (2011).
-  Grzywacz, N.M., Watamaniuk, S.N.J., McKee, S.P.: Temporal coherence theory for the detection and measurement of visual motion. Vision Res. 35(22), 3183-203 (1995)
-  Hladky, R.K., Pauls, S.D.: Minimal Surfaces in the Roto-Translation Group with Applications to a Neuro-Biological Image Completion Model. J. Math Imaging Vis 36, pp. 1-27 (2010)
-  Hoffman, W.C.: The Lie Algebra of Visual Perception. Journal of Mathematical Psychology 3, 65-98 (1966)
-  Hoffman, W.C.: The visual cortex is a contact bundle. Appl. Math. Comput. 32, 137-167 (1989)
-  Hörmander, L.: Hypoelliptic second order differential equations. Acta Math. 119(1), 147-171 (1967).
-  Hubel, D.H.: Eye, Brain and Vision. Scientific American Library, New York (1988)
-  Hubel, D.H., Wiesel, T.N.: Ferrier lecture: Functional architecture of macaque monkey visual cortex. Proceedings of the Royal Society of London Series B 198(1130), 1-59 (1977)
-  Jones, J.P., Palmer, L.A.: An evaluation of the two-dimensional Gabor filter model of simple receptive fields in cat striate cortex. J. Neurophysiol. 58(6), 1233-1258 (1987)
-  Kisvarday, Z.F., Toth, E., Rausch, M., Eysel, U.T.: Orientation-specific relationship between populations of excitatory and inhibitory lateral connections in the visual cortex of the cat. Cereb. Cortex 7(7), 605-18 (1997)
-  Koenderink, J.J., van Doorn, A.J: Representation of local geometry in the visual system. Biol. Cybernet. 55,367-375 (1987)
-  Ledgeway, T., Hess, R.F., Geisler, W.S.: Grouping local orientation and direction signals to extract spatial contours: empirical tests of ”association field” models of contour integration. Vision Res. 45(19), 2511-22 (2005)
-  Lee, T.S., Nguyen, M.: Dynamics of subjective contour formation in the early visual cortex. Proc. Natl. Acad. Sci. USA 98(4), 1907-11 (2001)
-  Lorenceau, J., Alais, D.: Form constraints in motion binding. Nat. Neurosci. 4(7), 745-51 (2001)
-  Malach, R., Schirman, T.D., Harel, M., Tootell, R.B., Malonek, D.: Organization of intrinsic connections in owl monkey area MT. Cereb. Cortex 7(4), 386-93 (1997)
-  Montgomery R., A tour of subriemannian geometries, their geodesics and applications. AMS 2002.
-  Mumford, D.: Elastica and computer vision. In: Bajaj, C.L. (eds.) Algebraic Geometry and its Applications, pp. 491-506. Springer-Verlag (1994)
-  Mumford, D., Shah, J.: Optimal approximation by piecewise smooth functions and associated variational problems. Comm. Pure Appl. Math. 17, 577-685 (1989)
-  Nagel, A., Stein, E.M., Wainger, S.: Balls and metrics defined by vector fields I: Basic properties. Acta Math. 155, 103-147 (1985)
-  Nitzberg, M., Mumford, D., Shiota, T.: Filtering, segmentation and depth. Lecture Notes in Computer Science 662, Springer 1993.
-  Øksendal, B.: Stochastic Differential Equations: an Introduction with Applications. Springer (2010)
-  Olson, I.R., Gatenby, J.C., Leung, H.C., Skudlarski, P., Gore, J.C.: Neuronal representation of occluded objects in the human brain. Neuropsychologia 42(1), 95-104 (2004)
-  Pan, Y., Chen, M., Yin, J., An, X., Zhang, X., Lu, Y., Gong, H., Li, W., Wang, W.: Equivalent representation of real and illusory contours in macaque V4. J. Neurosci. 32(20), 6760-70 (2012)
-  Pauls, S.D., Hladky, R.K.: Minimal surfaces in the roto-translation group with applications to a neuro-biological image completion model. J. Math. Imaging Vis. 36(1), 563-594 (2010)
-  Petitot, J., Tondut, Y.: Vers une neurogéométrie. Fibrations corticales, structures de contact et contours subjectifs modaux. Math. Sci. Hum. 145, 5-101 (1999)
-  Petkov, N., Subramanian, E.: Motion detection, noise reduction, texture suppression and contour enhancement by spatiotemporal Gabor filters with surround inhibition. Biol. Cybern. 97(5-6), 423-439 (2007)
-  Rainville, S.J.M., Wilson, H.R.: Global shape coding for motion-defined radial-frequency contours. Vision Res. 45(25-26), 3189-201 (2005)
-  Ringach, D.L.: Spatial structure and symmetry of simple-cell receptive fields in macaque primary visual cortex. J. Neurophysiol. 88(1), 455-463 (2002)
-  Roerig, B., Kao, J.P.: Organization of intracortical circuits in relation to direction preference maps in ferret visual cortex. J Neurosci. 19(24) (1999)
-  Rothschild, L.P., Stein, E.M.: Hypoelliptic differential operators and nilpotent Lie groups. Acta Math. 137, 247â320 (1977)
-  Robert, C.P., Casella, G.: Monte Carlo statistical methods. Springer, 2nd Edition (2004)
-  Sabinin, L.V.: Smooth quasigroups and loops. Kluwer (1999)
-  Sanguinetti, G., Citti, G., Sarti, A.: Image Completion Using a Diffusion Driven Mean Curvature Flowin A Sub-Riemannian Space. VISAPP (2)’08, 46-53 (2008).
-  Sanguinetti, G., Citti, G., Sarti, A.: A model of natural image edge co-occurrence in the rototranslation group. J. Vis. 10(14) (2010)
-  Sarti, A., Citti, G., Petitot, J.: The symplectic structure of the primary visual cortex. Biol. Cybern. 98, 33-48 (2008)
-  Scholl, B.J., Pylyshyn, Z.W.: Tracking multiple items through occlusion: clues to visual objecthood. Cogn. Psychol. 38(2), 259-290 (1999)
-  Sorba, P.: The Galilei group and its connected subgroups. J. Math. Phys 17(6), 941-953 (1976)