Geometrical Optical Illusion via Sub-Riemannian Geodesics in the Roto-Translation Group
Benedetta Franceschiello, Alexey Mashtakov, Giovanna Citti, Alessandro Sarti
GGeometrical Optical Illusion via Sub-RiemannianGeodesics in the Roto-Translation Group
B. Franceschiello , A. Mashtakov , G. Citti , and A. Sarti Fondation Asile des Aveugles and Laboratory for Investigative Neurophisiology,Department of Radiology, CHUV - UNIL, Lausanne Program Systems Institute of RAS, Russia, CPRC, Department of Mathematics, University of Bologna, Italy CAMS, Center of Mathematics, CNRS - EHESS,Paris, France [email protected], [email protected],[email protected], [email protected]
Abstract.
We present a neuro-mathematical model for geometrical op-tical illusions (GOIs), a class of illusory phenomena that consists in amismatch of geometrical properties of the visual stimulus and its associ-ated percept. They take place in the visual areas V1/V2 whose functionalarchitecture have been modelled in previous works by Citti and Sarti asa Lie group equipped with a sub-Riemannian (SR) metric. Here we ex-tend their model proposing that the metric responsible for the corticalconnectivity is modulated by the modelled neuro-physiological responseof simple cells to the visual stimulus, hence providing a more biologicallyplausible model that takes into account a presence of visual stimulus. Il-lusory contours in our model are described as geodesics in the new metric.The model is confirmed by numerical simulations, where we compute thegeodesics via SR-Fast Marching.
Geometrical-optical illusions (GOIs) have been discovered in the XIX century byGerman psychologists (Oppel 1854 [50], Hering, 1878,[33]) and have been definedas situations in which there is an awareness of a mismatch of geometrical prop-erties between an item in the object space and its associated percept [68]. Theseillusions induce a misjudgement of the geometrical properties of the stimulus, dueto the perceptual difference between the features of the presented stimulus andits associated perceptual representation. An historical survey of the discovery ofgeometrical-optical illusions is included in Appendix I of [68] and a classificationof these phenomena, can be found in Coren and Girgus, 1978, [15]; Robinson,1998, [54]; Wade, 1982, [65].The aim of this paper is to propose a mathematical model for GOIs basedon the functional architecture of low level visual cortex (V1/V2). This neuro-mathematical model will allow to interpret at a neural level the origin of GOIsand to reproduce the arised percept for this class of phenomena. The mainidea is to adapt the model for the functional geometry of V1 provided in [14] a r X i v : . [ q - b i o . N C ] F e b or perceptual completion. Here we extend it introducing a new metric for theconnectivity of the visual cortex, wich takes into account the output of simplecells in V1/V2, as a coefficient modulating the sub-Riemannian metric. We alsopostulate that geometrical optical illusory curves arise as geodesics in this newconnectivity metric between two given sets. Then we will adapt to this definitionthe SR Fast-Marching (SR-FM) algorithm introduced in [21,57] as tool for thecomputation of geodesics with fixed two-point boundary conditions (extremapoints). As a result we will be able to explain the perceptual phenomena by thegeometry of V1 and SR differential geometry instruments.The paper is organised as follows. In Section 2, we introduce the perceptualproblem of geometrical optical illusion and we review the state of the art concern-ing the previously proposed mathematical models. In Section 3, we briefly recallthe functional architecture of the visual cortex and the cortical based modelintroduced by Citti and Sarti in [14]. In Section 4, we introduce the neuro-mathematical model proposed for GOIs, taking into account the modulationof the functional architecture induced by the stimuli. In Section 5, we discusssub-Riemannian geodesics and the sub-Riemannian Fast-Marching. Finally, wedescribe the implementation in Section 6 and discuss the results in Section 7. Fig. 1.
From left to right. Hering illusion: two straight vertical lines in front of a radialbackground appear as if they were bowed outwards. Ehm–Wackermann illusion: thecontext of concentric circles bends inwards the edges of the square. Poggendorff illusion:the presence of a central surface induces a misalignment of the crossing transversals.
The phenomena we consider here consist in misperception effects induced byelements of the image. The Hering illusion, introduced by Hering in 1861 [33], ispresented in Figure 1, left. In this illusion two vertical straight lines are presentedn front of a radial background, so that the lines appear as if they were bowedoutwards. A similar effect is observable in the Ehm–Wackermann illusion [26],i.e. a square on a background composed by concentric circles, Figure 1, center.One more famous GOI is the Poggendorff illusion, which consists in an apparentmisalignment of two collinear, oblique, transversals separated by a rectangularsurface (Figure 1, right). For the latter, psychological elements contributing tothis misperception have been presented in [19,54] and [62].The interest in GOI comes from the chance to provide a better explanationof these phenomena, helping to understand the unrevealed mechanisms of vi-sion ([25]). Many studies, which rely on neuro-physiological and imaging data,show the evidence that neurons in at least two visual areas, V1 and V2, carrysignals related to illusory contours, and those signals in V2 are more robust thanin V1 ([64], [48], reviews [25], [47]). A more recent study measured the activatedconnectivity in and between areas of early visual cortices ([61]). To integrate themathematical model with the recent findings, we propose a neural-based modelto interpret GOI.
The first models of GOI are purely phenomenological and provide quantitative analysis of the perceived geometrical distortion, such as the angle deformation,which is the attitude of perceiving acute angles larger and obtuse ones as smaller.Models of this type have been proposed in 1971 by Hoffmann in terms of orbitsof a Lie group acting on the plane [35], and by Smith [60] in terms of differen-tial equations. More recently Ehm and Wackermann in [26] and [27] provided avariational approach expressed by a functional dependent on length of the curveand the deflection from orthogonality along the curve. These approaches do nottake into account the underlying neurophysiological mechanisms.On the other hand an entire branch for modelling neural activity, the Bayesianframework, had its basis in Helmholtz theory [32]: our percepts are our best guessas to what is in the world, given both sensory data and prior experience.
Thedescribed idea of unconscious inference is at the basis of the Bayesian statisticaldecision theory, a principled method for determining optimal performance in agiven perceptual task ([31]). An application of this theory to motion illusions hasbeen provided by Weiss et al in [67], by Geisler and Kersten in [31], by Fermüllerand Malm in [28]. In our model, we aim to combine psycho-physical evidenceand neurophysiological findings, in order to provide a neuro-mathematical de-scription of GOIs. It is inspired by the celebrated models of Hoffman [36] andPetitot [52,51], who have founded a discipline now called neuro-geometry, aimedto describe the functional architecture of the visual cortex with geometrical in-struments in order to explain phenomenological evidence. More recent contribu-tions are due to August and Zucker [3], Sarti and Citti [14,58], Duits et al. [23,24].A recent work trying to integrate the neuro-physiology of V1/V2 for explainingsuch phenomena has been presented in [29].
The classical neuromathematical model of V1/V2
The retina, identified as M ⊂ R , is the first part of the visual path initiatingthe signal transmission, which passes through the Lateral Geniculate Nucleusand arrives in the visual cortex, where it is processed.Let us consider a visual stimulus, i.e. an image I : M ⊂ R → R + . (1)The receptive field (RF) of a cortical neuron is the portion of the retina which theneuron reacts to, and the receptive profile (RP) ψ ( χ ) is the function that modelsits activation when a point χ = ( χ , χ ) ∈ R of the retinal plane is elicited bya stimulus at that point. To be specific, ( χ , χ ) are the local coordinates ofthe neighbourhood centered at x = ( x, y ) ∈ R , to which the neuron reactsto, while ( x, y ) refers to the global coordinates system of the retina R . Simplecells of V1 are sensitive to position and orientation of the contrast gradient of animage ∇ I . Their properties have been discovered by Hubel and Wiesel in[37] andexperimentally described by De Angelis in [20]. Considering a basic geometricmodel, the set of simple cells RPs can be obtained via translation on the vector x and rotation on the angle θ ∈ S (cid:39) SO of a unique mother profile ψ ( χ ) . Receptive fields have been modelled as oriented filters in the middle of 80’sand since then the orientation extraction in image analysis has been subject ofseveral works. The first models have been presented by Daugman [18] (1985),Jones and Palmer [39] (1987) in terms of Gabor filters. In the same years Youngin [69] (1987) and Koenderink in [40] (1990) proposed to model RPs as Gaussianderivatives (DoG). We also refer to [55] (2008) and [51] (2008) for further expla-nations and details. Recently a new class of multi-orientation filters have beenintroduced by Duits et al. in [22] (2007): cake-wavelets. A comparison betweencake-wavelets and Gabor filters efficiency has been presented in [6]. Having thescope of modelling the functionality of the visual cortex, we chose Gabor filters,proved to be a good model of receptive profiles and their spiking responses [53].We will consider odd and even part of Gabor filters in order to measure θ cor-rectly for both contours and lines. Definition 1.
A mother Gabor filter is given by ψ ( χ ) = α πσ e − ( χ α χ σ e iχ λ , χ = ( χ , χ ) ∈ R , (2) where λ > is the spatial wavelength, α > is the spatial aspect ratio and σ > is the standard deviation of the Gaussian envelope. s discovered by Lee [41], the set of simple cells RPs can be obtained viatranslation on the vector ( x, y ) ∈ R and rotation on the angle θ ∈ S of aunique mother profile ψ ( χ ) . Since the set of parameters ( x, y, θ ) describes theLie group SE of rotations and translations, we identify the set of receptiveprofiles (RPs) with this group. Let η = ( x, y, θ ) ∈ SE . Then an action of theLie group SE on the homogeneous space R is given by η (cid:12) χ = (cid:18) xy (cid:19) + (cid:18) cos θ − sin θ sin θ cos θ (cid:19) (cid:18) χ χ (cid:19) . (3)Denote η − = ( − x cos θ − y sin θ, x sin θ − y cos θ, − θ ) ∈ SE the inverse elementto η . A general RP can be expressed as ψ η ( χ ) = ψ ( η − (cid:12) χ ) . (4) The retina is the light-sensitive layer of shell tissue of the eye, where the visualstimulus is first imprinted. In our model, we follow [14] and identify the retinawith the planar area M ⊂ R .A visual stimulus is modelled by intensity function (1): I : M ⊂ R → R + : ( x, y ) (cid:55)→ I ( x, y ) . It activates the retinal layer of photoreceptors. Then the cortical neurons whosereceptive fields intersect the activated layer spike. We model their spiking re-sponses (or spiking activity) O I ( η ) = O I ( x, y, θ ) by a Gabor transform. Definition 2.
Let ψ η ∈ L ( R ) be a Gabor filter, given by (4) , modelling thereceptive profile of simple cells of the primary visual cortex. The continuousGabor transform of a signal I ∈ L ( R ) is defined as: O I ( η ) = (cid:90) R I ( χ ) ψ η ( χ ) dχ. (5)In the previous definition the output O I ( η ) of receptive profiles of simple cellsin V1 in response to a visual stimulus I ( x, y ) is mathematically described asa convolution. Let us note that Gabor filters are complex valued: the real andimaginary parts have a different role and detect different features. The real partis even and spikes maximally along lines, while the imaginary part is odd anddetects the presence of surfaces, i.e. contours. The term functional architecture refers to the organisation of cells of the pri-mary visual cortex in structures. The hypercolumnar structure, discovered bythe neuro-physiologists Hubel and Wiesel in the 60s ([38]), organizes the cellsf V1/V2 in columns (called hypercolumns) covering a small part of the visualfield R and corresponding to parameters such as orientation, scale, direction ofmovement, color, for a fixed retinal position ( x, y ) , Figure 2 (top). Fig. 2.
Top: representation of hypercolumnar structure, for the orientation parameter,where L and R represent the ocular dominance columns (Petitot [51]). Bottom: the setof all possible orientations for each position of the retina ( x, y ) . In our framework, over each retinal point we consider a whole hypercolumnof cells, sensitive to all possible orientations, see Figure 2 (bottom). Hence foreach position ( x, y ) of the retina M ⊂ R we associate a whole set of filters RP ( x,y ) = { ψ ( x,y,θ ) : θ ∈ S } . (6)This expression defines a fiber { θ ∈ S } over each point ( x, y ) ∈ R .In this framework the hypercolumnar structure is described in terms of differ-ential geometry, but further explanations are requested to model the orientationselectivity process performed by the cortical areas in the space of feature S ([14]). From the physiological point of view the orientation selectivity is the action ofshort range connections between simple cells belonging to the same hypercolumnto select the most probable instance from the spiking activation of receptiveprofiles in response to a stimulus.Mathematically, this process is modelled by assignment to every point x =( x, y ) ∈ R the angle ¯ θ ∈ S — the orientation of a line passing through thepoint x . It is found as the element of fiber that gives the maximal response of (5): ¯ θ ( x, y ) = argmax θ ∈ S | O I ( x, y, θ ) | . (7)his process is called lifting and it associates to each retinal point ( x, y ) thecorresponding maximal output ¯ θ ( x, y ) , denoting the selected orientation (tangentdirection) to the visual stimulus at point x .The other connectivity, responsible for the formation of contours in the cortexby given a retinal stimulus, is called horizontal connectivity. Horizontal connec-tions are long ranged and connect cells of approximately the same orientation ,belonging to different hypercolumns. Modelling this behaviour requires to endowV1 with a differential structure, see [14], where horizontal curves are the liftingof retinal curves to the extended space of positions and orientations — the Liegroup SE ∼ = R (cid:111) SO .The horizontal connectivity is therefore modelled as adiffusion along the integral curves of the left invariant vector fields on the group.The basis of left-invariant vector fields on SE is given by X = cos θ ∂∂x + sin θ ∂∂y , X = ∂ θ , X = − sin θ ∂∂x + cos θ ∂∂y . In order to model the propagation of the horizontal connectivity in R (cid:111) SO ,in [14] Citti and Sarti proposed to endow R (cid:111) SO with a sub-Riemannian metric. Definition 3.
A SR manifold is given by a triple ( M, ∆, G ) , where M is a con-nected, simply connected smooth manifold, ∆ is a smooth subbundle of the tan-gent bundle to M , and G is a metric defined on ∆ . In particular, in [14], the horizontal connectivity in the cortex is modelled bymeans of distance function, defined on the SR manifold ( M, ∆, G ) , where M = SE , ∆ = span( X , X ) , G = ω ⊗ ω + ω ⊗ ω . (8)Here ⊗ denotes the Kronecker product, and ω i ∈ T ∗ M denotes the basis oneform dual to X i , i.e. (cid:10) ω i , X j (cid:11) = δ ij , where δ ij is the Kronecker delta. In the previous section we provided neuro-geometrical tools that, starting fromthe neural counterpart of the visual cortex, explain the behaviour of V1/V2 inpresence of visual stimuli. The original contribution of this paper is to extendthe previous model [14] in the following way: starting from the sub-Riemannianmetric G in (8), we weight the long range connectivity taking into account theintra-columnar response of simple cells in V1/V2.We will obtain a new polarized metric G in (11). Illusory countours in ourmodel arise as local minimizers (i.e. geodesics) of this polarized metric. Here we introduce the idea that the isotropic cortical metric G defined on thehorizontal subbundle ∆ , see (8), can be modulated by the output of simple cells ig. 3. Top Left: Hering illusion, cf. Fig. 1.
Top Right:
Level line of the externalcost C ( x, y, θ ) for θ = 2 . rad. Bottom Left:
Ehm–Wackermann illusion, cf. Fig. 1.
Bottom Right:
Level line of the external cost C ( x, y, θ ) for θ = 2 . rad. O I ( η ) , induced by the visual stimulus I . This phenomena is a weak type of learn-ing, or pre-activation, where the activated cells are more sensitive to the corticalpropagation. The proposed modulation P ( O I ) of the metric induced by the vi-sual stimulus I is maximal in correspondence of the edges, and is expressed as P ( O I ) = P ( η ) = Re( O I ( η )) + Im( O I ( η )) , (9)where O I ( η ) is the output of simple cells as defined in (5).This formulation allows to detect both the presence of lines (first term) andthe presence and polarity of contours (second term). Once computed O I ( η ) fromthe initial images, we add to P ( η ) a positive values as the values of the output O I range from negative to positive value. The modification will make P ( η ) a positivevalue, which can be considered a polarization term for the metric. Finally, wenormalize it, obtaining the following external cost: C ( O I ) = C ( η ) = c + P ( η ) (cid:113) c + P ( η ) , (10)where c is a suitable positive constant.We define the polarized metric on the distribution ∆ = span( X , X ) as G = (cid:32) C ξ ( O I ) C ( O I ) (cid:33) , G − = (cid:18) C ξ ( O I ) 00 C ( O I ) (cid:19) , (11)where C ξ ( O I ) = ξ C ( O I ) . Here ξ > is a real parameter, which will be fixedand kept constant. It allows to weight differently the translational, i.e. X , X , ig. 4. Left: Poggendorff illusion, cf. Fig. 1.
Center: the initial stimulus with a secondtransversal corresponding to the perceptual completion.
Right:
A level set of C ( x, y, θ ) for θ = 2 . rad. The saturation is slightly lowered to show detected contours and lines. and rotational, i.e. X components of the metric. In other words it allows tomodulate the anisotropy between the retinical (on R ) and the hypercolumnar( S ) components. The choice of the constant ξ is discussed across experiments inSection 6. Let us notice we denoted the metric as G , where coefficient indicatesthat the metric has no contribution in the direction X . A curve γ : [0 , T ] (cid:55)→ SE is an integral curve of a vector field X starting fromthe point a iff ˙ γ ( t ) = X | γ ( t ) and γ (0) = a . In this last case γ will be also denoted exp( tX )( a ) = γ ( t ) . Denote M = (SE , ∆, G ) the SR manifold with polarized metric (11). Definition 4.
A Lipschitzian curve γ : [0 , T ] → SE on M is called horizontal,if ˙ γ ( t ) ∈ ∆ | γ ( t ) for a.e. t ∈ [0 , T ] . In other words, a horizontal curve γ on M is an integral curve of ˙ γ ( t ) = u ( t ) X | γ ( t ) + u ( t ) X | γ ( t ) , where u , u are real-valued functions from L ∞ ([0 , T ]) . Definition 5.
The SR length of a horizontal curve γ on M is defined as l ( γ ) := T (cid:82) (cid:112) G ( ˙ γ ( t ) , ˙ γ ( t )) d t. (12)or a general introduction to sub-Riemannian geometry see [45]. Note, that thevector fields X , X satisfy the Hörmander condition [45], i.e. their generated Liealgebra coincides with the tangent space at every point. Due to this property,the Rashevski (1938)–Chow (1939) theorem guarantees that any two elements η , η ∈ SE can be connected by a horizontal curve.Hence in a connected manifold M for any couple of points η and η thefollowing set is not empty: Γ η ,η = { γ horizontal curve , γ (0) = η , γ ( T ) = η } . As a consequence it is possible to define a distance on a connected manifold M . Definition 6.
The Carnot-Carathéodory distance on the sub-Riemannian man-ifold M between two points η and η is defined as d ( η , η ) = inf γ ∈ Γ η ,η l ( γ ) . (13)Note that Filippov’s theorem [1] implies existence of length-minimizers andinfimum in (13) can be replaced by minimum. Computation of sub-Riemannian (Carnot–Carathéodory) distance in general isa very difficult problem. For example, even in the case of left-invariant SR struc-tures on Lie groups the length-minimizers are known only in several simplestcases: the Heisenberg group [11,63], the groups SO , SU , and SL with theKilling metric [10], SE [56], SH [12], the Engel group [2], and 2-step corank2 nilpotent SR problems [4]. Our case M = (SE , ∆, G ) is much more diffi-cult than the case of a left-invariant SR metric, since the metric G dependson the functional parameter – the external cost C . Thus, to obtain an analyticexpression for SR distance (13) does not seem possible.Instead, based on idea of Riemannian approximation [57], we build a numer-ical method to compute the SR-distance as a limiting case of the Riemanniandistances, when the penalization of movement in the direction X (forbidden inSR case) tends to infinity. In other words, the Riemannian approximation relaxesthe horizontality constraint ˙ γ ∈ ∆ and extends the SR metric G to the highlyanisotropic Riemannian metric G defined in whole tangent bundle T SE .There are several possible definitions for Riemannian distance functions whichapproximate a Carnot-Carathéodory distance in the Gromov-Hausdorff sense.We use the following Riemannian approximation for the SR metric G . Definition 7 (Riemannian approximation of the SR metric).
A Rieman-nian approximation G (cid:15) of the sub-Riemannian metric G is defined over the wholetangent bundle T SE and has the following expression in the frame ( X , X , X ) : G (cid:15) = diag (cid:16) C ξ , C , (cid:15) C ξ (cid:17) , where ξ > , (cid:15) > are parameters of the metric anisotropy, C = C ( O I ) isexternal cost (10) and C ξ = ξ C .emark 1. Note that the metric blows up as (cid:15) tends to . On the other handthe inverse of the metric, computed as G − (cid:15) = diag (cid:16) C ξ , C , (cid:15) C ξ (cid:17) formally tends to G − .For every (cid:15) > we denote as d (cid:15) ( · , · ) the Riemannian distance associated to themetric G (cid:15) . The following result, proved in [30, Theorem 1.1] in general settings,provide a relationship between Riemannian and sub-Riemannian distance: Theorem 1.
Let M = SE . The sequence ( M, d (cid:15) ) converges to ( M, d ) as (cid:15) → in the Gromov-Hausdorff sense. See also [46, Theorem 1.2.1] for another related Riemannian approximationscheme and [21] for an extension to more general Finsler metrics by Duits etal. In order to better understand this assertion we provide explicit estimates ofthe approximated distance. To this end, let us recall that the exponential mapis a local diffeomorphism around each point (see [49]).
Proposition 1 ([49]).
For every fixed point η ∈ SE the function Φ η ( ζ ) = exp (cid:32) (cid:88) i =1 ζ i X i (cid:33) ( η ) , ζ = ( ζ , ζ , ζ ) ∈ R , is a local diffeomorphism around the point η and its inverse: Φ − η defines localcoordinates in a neighborhood of η . To better describe the dependence of the distance d (cid:15) on the parameter (cid:15) , letus define the regularized gauge function: N ( ζ ) = (cid:113) ζ + ζ + | ζ | , N (cid:15) ( ζ ) = (cid:113) ζ + ζ + min (cid:0) | ζ | , (cid:15) − ζ (cid:1) , for (cid:15) > and the associated pseudo-distance function d (cid:15) ( η , η ) = N (cid:15) ( Φ − η ( η )) which provides an estimate of the distance d (cid:15) for both (cid:15) = 0 and (cid:15) > (see [13]). Lemma 1 ([13]).
There exists
A > , independent of (cid:15) , such that for all η , ηA − d (cid:15) ( η, η ) ≤ d (cid:15) ( η, η ) ≤ A d (cid:15) ( η, η ) . (14)The dependence from (cid:15) becomes now quite clear. Indeed for every fixed ζ : (cid:15) → ⇒ N (cid:15) ( ζ ) → N ( ζ ) , which is independent of (cid:15) , and provides an estimate for the SR distance. GOIs as Sub-Riemannian Geodesics
A second original aspect of our model is the way we recover subjective contours.Many results using sub-Riemannian (SR) geodesics as model for subjective con-tours in SE are already present in literature. SR geodesics and their applicationto image analysis were also studied in [9,34,43], e.g. to retinal vessel tracking byBekkers et al. [8,42,7]. Explicit formulas for SR geodesics and optimal synthesisin SE are obtained by Sachkov [56]. Let us also underline that in [14] geodesicsarise as foliation of subjective surfaces.The problem we face here is more general, since we do not know the exactposition of the geodesic extrema. Let us consider an example, the Hering illusion:it presents a misperception of two vertical straight lines, perceived as bowedoutwards. The perceived curves are modelled as geodesics of the polarized metricin SE , but we only know the spatial R -components of its extrema ( ( x , y ) and ( x , y ) ), while the angles θ and θ are unknown. As a result the reconstructedperceptual curve is described as the minimizing geodesic between two a prioriknown sets, obtained by fixing the spatial component ( x, y ) ∈ R and varyingthe orientation θ ∈ S .We can generalize this problem by means of the following definition. Definition 8.
Let K , K ⊂ M = SE be compact and non empty sets.For fixed (cid:15) ≥ , consider the (Riemannian, if (cid:15) > , or SR, if (cid:15) = 0 ) metric G (cid:15) .We call (cid:15) -minimizing geodesic with extrema in the sets K and K the curve γ (cid:15) (horizontal, if (cid:15) = 0 ), such that l (cid:15) ( γ (cid:15) ) = min ˜ γ { l (cid:15) (˜ γ ) | ˜ γ : [0 , T ] → M, ˙˜ γ ( t ) ∈ ∆ (cid:15) , ˜ γ (0) ∈ K , ˜ γ ( T ) ∈ K } . (15) Here l (cid:15) denotes the length in G (cid:15) , and ∆ (cid:15) = span( X , X , (cid:15)X ) ⊆ T M , (cid:15) ≥ . Note that minimum in (15) exists due to compactness of the sets K and K and existence of a minimizing geodesic connecting any two given points η ∈ K and η ∈ K , as we will formally show in Proposition 2.In other words, an (cid:15) -minimizing geodesic realizes the distance d (cid:15) between twosets K and K , where the distance is defined in the following. Definition 9.
Let K , K ⊂ SE be compact non empty sets.The distance function from K (Riemannian if (cid:15) > or SR if (cid:15) = 0 ) is defined as d (cid:15), K ( η ) = inf η ∈ K d (cid:15) ( η, η ) . Hence, the distance function between two sets K , K is defined as d (cid:15) ( K , K ) = inf η ∈ K d (cid:15), K ( η ) . urthermore, if γ is (cid:15) -minimizing geodesic with extrema in K and K , then γ (0) ∈ K , γ ( T ) ∈ K , l (cid:15) ( γ ) = d (cid:15) ( K , K ) . Remark 2.
The special case in which K = { η } and K = { η } contain only onepoint, the previous problem reduces to find the length minimizing curve between η and η . The existence of a minimum is well known and it is called minimizinggeodesic. We refer to [45] for general properties of minimizing geodesics.From the existence of a geodesic with fixed extrema, we can deduce the existenceof (cid:15) -minimizing geodesic in the sense of Definition 8. Proposition 2.
In the assumptions of Definition 8 the minimum in (15) exists.Proof.
Indeed we can find a sequence η ,n in K and a sequence η ,n in K such that d (cid:15) ( η ,n , η ,n ) tends to d (cid:15) ( K , K ) . Since η ,n and η ,n are boundedthey have a limit, respectively η and η . The geodesic curve between these twopoints exists by Remark 2 and gives the minimum in (15).From the convergence of the distance d (cid:15) to the distance d as (cid:15) → , we candeduce the following proposition: Proposition 3.
Let K , K ⊂ SE be compact non empty sets.The following convergence result holds: d (cid:15) ( K , K ) → d ( K , K ); d (cid:15),K ( η ) → d ,K ( η ) as (cid:15) → Proof.
By definition of distances d (cid:15) and d we immediately see that ∀ (cid:15) > : ≤ d (cid:15) ( K , K ) ≤ d ( K , K ) . For any sequence (cid:15) j → as j → ∞ consider the sequence d (cid:15) j ( K , K ) . Since d (cid:15) ( K , K ) is bounded, d (cid:15) j ( K , K ) is converging. Clearly the limit l satisfies l ≤ d ( K , K ) . Let us assume by contradiction that d (cid:15) j ( K , K ) converges to alimit l < d ( K , K ) . For every j there exist η ,(cid:15) j ∈ K and η ,(cid:15) j ∈ K such that d (cid:15) j ( η ,(cid:15) j , η ,(cid:15) j ) = d (cid:15) j ( K , K ) . Since ( η ,(cid:15) j ) j ∈ N and ( η ,(cid:15) j ) j ∈ N are bounded, then they have a convergent subse-quences: η ,j → η and η ,j → η . Now we note that | d (cid:15) j ( η ,(cid:15) j , η ,(cid:15) j ) − d (cid:15) j ( η ,(cid:15) j , η ) | + | d (cid:15) j ( η ,(cid:15) j , η ) − d (cid:15) j ( η , η ) | → , since d (cid:15) j is Lipshitz continuous, with Lipshitz constant . Hence d (cid:15) j ( η ,(cid:15) j , η ,(cid:15) j ) ≤ | d (cid:15) j ( η ,(cid:15) j , η ,(cid:15) j ) − d (cid:15) j ( η ,(cid:15) j , η ) | ++ | d (cid:15) j ( η ,(cid:15) j , η ) − d (cid:15) j ( η , η ) | + d (cid:15) j ( η , η ) → d ( η , η ) and as a result d ( K , K ) ≤ d ( η , η ) = lim (cid:15) j → d (cid:15) j ( η ,(cid:15) j , η ,(cid:15) j ) , and this is a contradiction, so that d (cid:15) j ( K , K ) → d ( K , K ) . Since any se-quence has the same limit, then d (cid:15) ( K , K ) → d ( K , K ) as (cid:15) → . .2 Riemannian and sub-Riemannian eikonal equation Now we show that the distance function from a set satisfies a first order PDEcalled eikonal equation. We first recall the notion of a (sub-)Riemannian gradient.
Definition 10.
The Riemannian gradient of a function f in the metric G (cid:15) isthe vector ∇ (cid:15) f = (cid:88) j =1 G ij(cid:15) ( X j f ) X i . For (cid:15) = 0 we obtain the sub-Riemannian gradient ∇ f = (cid:88) j =1 G ij ( X j f ) X i . In the Riemannian setting it is well know that the distance function from aset is a viscosity solution (in the sense by [17,16]) of the eikonal equation:
Proposition 4 ([17,16]).
Let K be a compact non empty set with C ∞ boundaryand (cid:15) > be a constant. Then in the points of differentiability outside the set K the Riemannian distance function d (cid:15),K ( η ) satisfies the equation (cid:107)∇ (cid:15) d (cid:15),K ( η ) (cid:107) G (cid:15) = 1 , and d (cid:15),K ( η ) vanishes in K . In Carnot groups with SR metric the same assertion has been proved in [46].This result can be extended to the present setting:
Proposition 5.
The distance function d ,K given by Definition 9 satisfies thefollowing eikonal equation: (cid:26) (cid:107)∇ d ,K ( η ) (cid:107) G = 1 , for η / ∈ K,d ,K ( η ) = 0 , for η ∈ ∂K. (16)We omit the proof which is similar to the one contained in [46]. Viceversa thefollowing result holds: Proposition 6.
The problem (cid:26) (cid:107)∇ u (cid:107) G = 1 , for η / ∈ K,u = 0 , for η ∈ ∂K. (17) has a unique viscosity solution, u which coincides with the distance function d ,K ( η ) from the set K in the sub-Riemannian setting (for (cid:15) = 0 ) or Riemanniansetting (for (cid:15) > ). The proof can be obtained working as in [5] .3 Sub-Riemannian Fast Marching
One of the most efficient method to compute geodesics in the Euclidean settingis Fast-Marching, introduced by Sethian in [59]. In case of geodesics between twopoints it has been extended by Mirebeau (in the case of Riemannian metric [44]),and by Sanguinetti et al in [57] (in the SE setting with a general SR metric).The Fast-marching method works as follows: – First the distance map, from the initial point is computed as viscosity solu-tion of the eikonal equation. – Then a backtracking procedure is applied. The latter is based on the re-lationship between the gradient of the distance function and the directionof the geodesics far from cut points, i.e. points where the geodesics losesminimality (see [1]).In particular, in [8] it was shown that if η (cid:54) = η ∈ SE and the unique minimizinggeodesic γ (cid:15) : [0 , T ] → SE joining η and η does not contain cut points, then ˙ γ ( t ) = ∇ (cid:15) d (cid:15) ( γ ( t ) , η ) . As a consequence, they show that the geodesic can berecovered with the following backtracking procedure: Proposition 7.
By given two distinct points η (cid:54) = η ∈ SE consider a geodesic γ (cid:15) : [0 , T ] → SE joining η and η . If γ does not contain cut points then γ (cid:15) ( t ) = γ (cid:15)b ( T − t ) , where γ b,(cid:15) is a solution of the problem (cid:26) ˙ γ (cid:15)b ( t ) = −∇ (cid:15) d (cid:15) ( γ (cid:15)b ( t ) , η ) , t ∈ [0 , T ] ,γ (cid:15)b (0) = η. (18) Here we extend backtraking procedure (18) in Proposition 7 to the geodesicswith extrema in a set, introduced in Definition 8. Before that we make a remark.
Remark 3. If K and K are compact sets with smooth boundary and d (cid:15),K attains a minimum at the point η on the set K , then the minimizing geodesicwith extrema in K and K coincides with the minimizing geodesic with firstextremum in K and second extremum η .As a consequence, the minimizing geodesic can be found via backtracking ofthe distance from the starting set: Proposition 8.
Let (cid:15) > . If the following assumptions are satisfied:1. K , K ⊂ SE are compact non empty sets with smooth boundary;2. d (cid:15),K attains a minimum at the point η on the set K ;3. minimizing geodesic γ (cid:15) joining K and η does not contain cut points;then γ (cid:15) ( t ) = γ (cid:15)b ( T − t ) , where γ (cid:15)b is a solution of the problem: (cid:26) ˙ γ (cid:15)b ( t ) = −∇ (cid:15) d (cid:15),K ( γ (cid:15)b ( t )) , t ∈ [0 , T ] γ (cid:15)b (0) = η . (19) roof. We can assume that γ (cid:15) is parametrized by arclength. Then we have t = d (cid:15),K ( γ (cid:15) ( t )) and differentiating with respect to t we get (cid:104)∇ (cid:15) d (cid:15),K ( γ (cid:15) ( t )) , ˙ γ (cid:15) ( t ) (cid:105) (cid:15) ≤ (cid:107)∇ (cid:15) d (cid:15),K ( γ (cid:15) ( t )) (cid:107)(cid:107) ˙ γ (cid:15) ( t ) (cid:107) ≤ . Since equality holds, then ˙ γ (cid:15) ( t ) is parallel to ∇ (cid:15) d (cid:15),K ( γ (cid:15) ( t )) . But they have thesame norm, hence ˙ γ (cid:15) ( t ) = ∇ (cid:15) d (cid:15),K ( γ (cid:15) ( t )) . Since γ b,(cid:15) ( t ) = γ (cid:15) ( T − t ) , then ˙ γ (cid:15) ( t ) = −∇ (cid:15) d (cid:15)K ( γ (cid:15) ( t )) . In the limit for (cid:15) → we recover a minimizing geodesic for the SR problem. Corollary 1. If K and K are compact non empty sets with smooth boundary,and for every (cid:15) > the (cid:15) -minimizing geodesic γ (cid:15) : [0 , T ] → SE joining K and K does not contain cut points, then there exists lim (cid:15) → γ (cid:15) = γ and γ is the minimizing SR geodesic joining K and K . Clearly the minimizing geodesic between two sets can be as well computed viaminimization on the complete set of geodesics connecting points of K and K . The imaginary part of the output corresponds to the response of odd Gaborfilters (contours polarity), while the real part performs the orientation detectionof lines (line detection). We assume that the orientation domain θ takes values in [ − π, π ) . As Re ( O I ( x, y, θ )) has π -periodicity, the energy volume is duplicated toinsure a well definition of the activation responses over the whole angular domain [ − π, π ) . Let us notice that while combining odd and even receptive profile, sinceodd receptive profiles over a straight line do not produce any response, we tookdifferent scales for even and odd filters: even Gabor profiles need to be sharpto detect line orientation, while odd Gabor profiles need to be wider to detectwhether they are aligned or not along a surface contour. The first step performed consists into the convolution of the initial image with abank of even and odd Gabor filters. A response O I is produced and opportunelycombined to obtain P , as described in (9). P , corresponding to the polarization ofour SR metric, is shifted to positive values and normalized to obtain C ( x, y, θ ) ,finally used as weight for the connectivity (figure 3, right column). The SRgeodesic that solves (13) is obtained in two steps:. Computation of the distance map solving (16) via Sub-Riemannian Fast-Marching, see figure 5, left;2. Computation of the geodesic by gradient descent (19).The constructed metric in R (cid:111) SO is a Riemannian approximation of the SRmetric, weighted by the external cost C ( x, y, θ ) . When switching from imagecoordinates to mathematical coordinates one should take care of correctly eval-uating ξ , which represents the anisotropy between the two horizontal direction, ξ∆x = ∆θ , where ∆x, ∆θ are the discretization steps along x and θ . In theexperiments for Hering and Ehm–Wackermann illusions, we set (cid:15) = 0 . , ξ = 7 ,while it varied proportionally to the geometrical elements of the image (entrytransversal and width of the surface) in the experiments for Poggedorff illusion.As was shown in [57], (cid:15) is taken sufficiently small to give an accurate approxi-mation of the SR-case. We processed the initial stimuli of the illusions through the method presentedin Section 4.1 and implemented in Section 6.2.
Fig. 5. From left to right : (1) minimum of distance map d (cid:15),K ( η ) from the boundaryvalue condition (initial seed) of equation (16), along the direction θ , computed throughSR-Fast-Marching. (2): 2D projection of the computed geodesics. The perceptual curveis cyan. (3): 3D plot of the computed geodesics. ξ = 4 . The Hering illusion, introduced by Hering in 1861 [33] is presented in Figure 1,left. In this illusion two vertical straight lines are presented in front of radialbackground, so that the lines appear as if they were bowed outwards. As de-scribed in the previous sections, we first convolve the distal stimulus with thentire bank of Gabor filters to compute the polarization of the metric P ( x, y, θ ) :we take 32 orientations selected in [0 , π ) , σ = 7 . pixels, α = 0 . pixels. The re-sulting computed perceptual curves are shown in figure 6. In order to determinethe perceptual angle, we varied θ between (0 , π/ as starting set, using − θ asending set. Fig. 6.
Representation of the computed perceptual curves.
This illusion, introduced by Ehm and Wackermann in [26], consists in presentinga square over a background of concentric circles, Figure 1, center. This context,the same we find in Ehrenstein illusion, bends the edges of the square toward thecenter of the image. Here we take the same number of orientations, 32, selectedin [0 , π ) , σ = 10 pixels, α = 0 . pixels. In order to determine the perceptualangle, we varied θ between (0 , π/ as starting set, using − θ as ending set. Theresulting computed perceptual curves are shown in Figure 7. Manipulating the elements of Poggendorff to understand how to magnify theillusory phenomenon has been done in many works [19,66]. In [66], the authorsperformed psychophysical experiments to obtain quantitative measures of themagnitude of the illusion: the illusory effects increased with increasing separationbetween the parallels as well as increasing the width of the obtuse angle formedby the transversal. We were not able to estimate computationally the effectinduced by obtuse angles because of the interaction of the parameter ξ in thesub-Riemannian Fast-Marching.Here we consider odd Gabor filters with the following values: α = 1 , θ ∈ ( − π, π ) ( values for Even Gabor profiles, 62 for Odd Gabor filters), σ = 3 . pixels (Even) and σ = 7 . pixels (Odd). The dimensions of images are × ig. 7. Representation of the computed perceptual curves. pixels. The scale parameter σ is chosen in relationship with image resolution andis taken smaller for even Gabor filters, to construct filter sharp enough to detectlines. On the other hand, σ for odd Gabor filters is taken bigger, to detect thepresence of surface and obtaining vanishing integral along lines. For the entrytransversal, we chose θ = π/ , π/ , π/ , π/ and width = , , pixels. Wecomputed the distance between the entry trasversal and the set containing theending points on the right side of the surface. The shortest curves computedthrough this model are in accord with the perceptual expectation. The anglevariation of the transversal, create an increased obtuse angle effect ( θ = π/ )and a non illusory effect ( θ = π/ ). In Figure 8, all the 2D projections of thecomputed geodesics for varying transversal entry angle and surface width ispresented. When θ = π/ , no illusion is shown and the geodesic is a straight line. Discussion
In this paragraph we show a table reporting the collected data con-cerning the SR lengths of the computed curve. It refers to the change of lengthvarying the widths and angles, underlining the observed phenomena.Type of curve Width = 7 pixels Width = 15 pixels Width = 25 pixelsPercep. curve θ = π/ θ = π/ θ = π/ θ = π/ θ = π/ θ = π/ Now we consider a variant of the Poggendorff illusion, called Round Poggendorff,see Fig. 10, left. The presence of the central surface induces a misperception: the ig. 8.
Poggendorff stimuli processed with their corresponding computed geodesicsoverlapped. In red we show an undersampling of geodesics computed from the leftentry transversal to the right side of the central surface. Varying the orientation fromLeft to Right: θ = π/ , π/ , π/ , π/ and with fixed width, Top: pixels, Central: pixels, Bottom: pixels arches do not seem cocircular and the left arc seems to be projected to somepoint with a certain orientation on the left bar. As in the previous example, weprovide a terminal set to the SR-Fast Marching: the seed is fixed at the crossingpoint between the right arc and the right bar, ξ = 2 . and possible terminalorientations belong to [ − π/ , , where θ = 0 is the angle corresponding to theorthogonal projection over the left bar and θ = − π/ is the boundary conditionof the circle at crossing point with the left bar. The SR length of the geodesic is1.32668 and the corresponding computed end point is { . , . , − . } . In this paper a neuro-mathematical model for the geometrical optical illusions ispresented, based on the functional architecture of V1. Perceptual curves arise asgeodesics of a polarized metric in SE , directly induced by the visual stimulus.The geodesics are computed through SR-FM and the perceptual curves result to ig. 9. Reconstruction of the perceptual bars, from the y coordinate corresponding tothe cyan curve of 8. Varying the orientation from Left to Right : θ = π/ , π/ , π/ for fixed width Top: pixels, Central: pixels, Right: pixels ig. 10. Left: Round Poggendorff illusion, courtesy of Talasli et al. see [62].
Center:
A family of geodesics starting from ( x , y , θ ) with multiple endpoints. The aim is todetermine ( y, θ ) minimizing the length of the perceptual curve. Right:
A minimizerhas end point ( y, θ ) = (0 . , − . be shorter (w.r.t. SR-metric) than the corresponding geometrical continuation.The model has been compared with psychophysical evidences which explain howthe effect varies depending on the width of the central surface and the angle of thetransversal. Improving the understanding of these phenomena is very importantbecause it can lead to insights about the behaviour of the visual cortex [25],allowing new applications to be developed. Acknowledgment and Contributions
The research leading to these results has received funding from the People Pro-gramme (Marie Curie Actions) of the European Union’s Seventh FrameworkProgramme FP7/2007-2013/ under REA grant agreement n607643.The work of Mashtakov is supported by the Russian Science Foundation un-der grant 17-11-01387 and performed in Ailamazyan Program Systems Instituteof Russian Academy of Sciences. The work on the revised manuscript was par-tially supported by the European Community’s Seventh Framework Programme(FP7/2007-2013) / ERC grant Lie Analysis, agr. n. 335555, and by the Euro-pean Community’s project Marie SkÅĆodowska-Curie grant agreement GHAIAn. 777822.We thank the anonymous reviewer for the valuable comments, which helpedto considerably improve the quality of the manuscript. We thank A. Ardentovfor the fruitful discussions on the revised manuscript.