A Cortical-Inspired Sub-Riemannian Model for Poggendorff-Type Visual Illusions

We consider Wilson-Cowan-type models for the mathematical description of orientation-dependent Poggendorff-like illusions. Our modelling improves two previously proposed cortical-inspired approaches, embedding the sub-Riemannian heat kernel into the neuronal interaction term, in agreement with the intrinsically anisotropic functional architecture of V1 based on both local and lateral connections. For the numerical realisation of both models, we consider standard gradient descent algorithms combined with Fourier-based approaches for the efficient computation of the sub-Laplacian evolution. Our numerical results show that the use of the sub-Riemannian kernel allows us to reproduce numerically visual misperceptions and inpainting-type biases in a stronger way in comparison with the previous approaches.

1. Introduction.The question of how we perceive the world around us has been an intriguing topic since the ancient times.For example, we can consider the philosophical debate around the concept of entelechy, which started with the early studies of the Aristotelian school to answer this question, while, on the side of phenomenology and on its relation to natural sciences, we can think of the theory started by Husserl.A well-known and accepted theory of perception is the one formulated within Gestalt psychology [1,2].
The Gestalt psychology is a theory for understanding the principles underlying the configuration of local forms giving rise to a meaningful global perception.The main idea of the Gestalt psychology is that the mind constructs the whole by grouping similar fragments rather than simply summing the fragments as if they were all different.In terms of visual perception, such similar fragments correspond to point stimuli with same (or very close) valued features of the same type.As an enlightening example from vision science, we tend to group the same colored objects in an image and to perceive them as an ensemble rather than the objects with different colors.There have been many psychophysical studies which attempted to provide quantitative parameters describing the tendencies of the mind in visual perception based on Gestalt psychology.A particularly important one is the pioneering work of Field et.al. [3], where the authors proposed a representation, called association field, modelling specific Gestalt principles.Furthermore, they also showed that it is more likely that the brain perceives together fragments which are similarly oriented and aligned along a curvilinear path than the ones rapidly changing orientations.
The presented model for neuronal activity is a geometrical abstraction of the orientation-sensitive V1 hypercolumnar architecture observed by Hubel and Wiesel [4][5][6].This abstraction generates a good phenomenological approximation of the V1 neuronal connections existing in the hypercolumnar architecture as reported by Bosking et.al. [7].In this framework, the corresponding projections of the neuronal connections in V1 onto a 2D image plane are considered to be the association fields described above and the neuronal connections are modeled as the horizontal integral curves generated by the model geometry.The projections of such horizontal integral curves were shown to produce a close approximation of the association fields, see Figure 1.For this reason, the approach considered by Citti, Petitot and Sarti and used in this work is referred to be cortically-inspired.
We remark that the presented model for neural activity is a phenomenological model which provides a mathematical understanding of early perceptual mechanisms at the cortical level by starting from very structure of receptive profiles.Nevertheless, it has been very useful for many image-processing applications, see for example, [8,9].Fig. 1: Projections of horizontal integral curves approximate the association fields from the experiment of Field, Hayes and Hess [3].They are generated by the sub-Riemannian model geometry proposed by Citti and Sarti [10].Figures are adapted from [3] and [10].
In this work, we follow this approach for a better understanding of the visual perception biases due to visual distortions often referred to as visual illusions.Visual illusions are described as the mismatches between the reality and its visual perception.They result either from a neural conditioning introduced by external agents such as drugs, microorganisms, tumors [11,12] or from self-inducing mechanisms evoking visual distortions via proper neural functionality applied on specific stimulus [13,14].The latter type of illusions are due to the effects of neurological and biological constraints on the visual system, [15].
In this work, we focus on illusions induced by contrast induction and orientation misalignments, with a particular focus on the well-known Poggendorff illusion and its variations, see Figure 2.This is a geometrical optical illusion [16,17] in which a misaligned oblique perception is induced by the presence of a central bar [18].Fig. 2: The original Poggendorff illusion: The red colored line is aligned with the black line although the blue one is falsely perceived as its continuation.Source: Wikipedia.
1.1.The functional architecture of the primary visual cortex .It has been known since the celebrated experiments of Hubel and Wiesel [4][5][6] that neurons (simple cells) in the primary visual cortex (V1) perform boundary (hence orientation) detection and propagate their activations through the cortical connectivity, in accordance with the psychophysical results of Fields and Hayes [3].Hubel and Wiesel showed that simple cells have a spatial arrangement based on the so-called hypercolumns in V1.In this arrangement, simple cells which are sensitive to different orientations at the same retinal location are found in the same vertical column constructed on the cortical surface.Adjacent columns contain the simple cells which are sensitive to close positions.
Several models have been proposed to describe the functional architecture of V1 and the neural connectivity within it.Koenderink et.al. [19,20] focused on differential geometric approaches to study the visual space where they modelled the invariance of simple cells with respect to suitable symmetries in terms of a family of Gaussian functions.Hoffman [21,22] provided the basic framework of vision models by interpreting the hypercolumn architecture of V1 as a fiber bundle.Following a similar reasoning, Petitot and Tondut [23] further developed this modelling, providing a new model, coherent both with the structure of orientation sensitive simple cells and the long range neural connectivity between V1 simple cells.In their model, they first observed that the simple cell orientation selectivity induces a contact geometry (associated to the first Heisenberg group) rendered by the fibers of orientations.Moreover, they showed that a specific family of curves found via a constrained minimization approach in the contact geometry fits the aforementioned association fields reported by Field et.al. [3].In [10,24], Citti and Sarti further developed the model of Petitot and Tondut, by introducing a group based approach, then refined by Boscain, Gauthier et al. [25,26].See also the monograph [27].The so-called Citti-Petitot-Sarti (CPS) model exploits the natural sub-Riemannian (sR) structure of the group of rotations and translations SE(2) as the V1 model geometry.
In this framework, simple cells are modelled as points of the three-dimensional group M = R 2 ×P 1 .Here, P 1 is the projective line, obtained by identifying antipodal points in S 1 .The response of simple cells to the two-dimensional visual stimuli is identified by lifting them to M via a Gabor wavelet transform.Neural connectivity is then modelled in terms of horizontal integral curves given by the natural sub-Riemannian structure of M. Activity propagation along neural connections can further be modelled in terms of diffusion and transport processes along the horizontal integral curves.
In recent years, the CPS model have been exploited as a framework for several cortical-inspired image processing problems by various researchers.We mention the large corpus of literature by Duits et al., see e.g., [28][29][30] and the state-of-the-art image inpainting and image recognition algorithms developed by Boscain, Gauthier, et al. [8,31].Some extensions of the CPS model geometry and its applications to other image processing problems can be found in [32][33][34][35][36][37][38][39].
1.2.Mean-field neural dynamics & visual illusions.Understanding neural behaviors is in general a very challenging task.Reliable responses to stimuli are typically measured at the level of population assemblies comprised by a large number of coupled cells.This motivates to reduce, whenever possible, the dynamics of a neuronal population to a neuronal mean-field model which describes largescale dynamics of the population as the number of neurons goes to infinity.These mean-field models, inspired by the pioneering work of Wilson and Cowan [40,41], and Amari [42], are low dimensional in comparison to their corresponding ones based on large-scale population networks.Yet, they capture the same dynamics underlying the population behaviors.
In the framework of the CPS model for V1 discussed above, several mathematical models were proposed to describe the neural activity propagation favoring the creation of visual illusions, including Poggendorff type illusions.In [37], for instance, illusions are identified with suitable strain tensors, responsible of the perceived displacement from the grey levels of the original image.In [43], illusory patterns are identified by a suitable modulation of the geometry of SE(2) = R 2 × S 1 and computed as the associated geodesics via the fast-marching algorithm.
In [44][45][46], a variant of the Wilson-Cowan (WC) model based on a variational principle and adapted to the M geometry of V1 was employed to model the neuronal activity and generate illusory patterns for different illusion types.The modelling considered in these works is strongly inspired by the integrodifferential model firstly studied in [47] for perception-inspired Local Histogram Equalization (LHE) techniques and later applied in a series of work, see, e.g., [48,49] for the study of contrast and assimila-tion phenomena.By further incorporating a cortical-inspired modelling, the authors showed in [44][45][46] that cortical LHE models are able to replicate visual misperceptions induced not only by local contrast changes, but also by orientation-induced biases similar as the ones in Figure 2. Interestingly, the cortical LHE model [44][45][46] was further shown to outperform both standard and cortical-inspired WC models and rigorously shown to correspond to the minimization of a variational energy, which suggests better efficient representation properties [50,51].One major limitation in the modelling considered in these works is the use of neuronal interaction kernels (essentially, isotropic 3D Gaussian) which are not compatible with the natural sub-Riemannian structure of V1 proposed in the CPS model.

Main contributions.
In this work, we encode the sub-Riemannian structure of V1 into both WC and LHE models by using a sub-Laplacian procedure associated with the geometry of the space M described in Section 1.1.Similarly as in [44][45][46], the lifting procedure associating to a given two dimensional image the corresponding neuronal response in M is performed by means of all-scale cake wavelets, introduced in [52,53].A suitable gradient-descent algorithm is applied to compute the stationary states of the neural models.
Within this framework, we study the family of Poggendorf visual illusions induced by local contrast and orientation alignment of the objects in the input image.In particular, we aim to reproduce such illusions by the proposed models in a way which is qualitatively consistent with the psychophysical experience.
Our findings show that it is possible to reproduce Poggendorff-type illusions by both the sR cortical-inspired WC and LHE models.This, compared with the results in [44,45] where the cortical WC model endowed with a Riemannian (isotropic) 3D kernel was shown to fail to reproduce Poggendofftype illusions, shows that adding the natural sub-Laplacian procedure in the computation of the flows improves the capability of those cortical-inspired models in terms of reproducing orientation-dependent visual illusions.

Cortical-inspired modelling.
In this section we recall the fundamental features of CPS models.The theoretical criterion underpinning the model relies on the so-called neurogeometrical approach introduced in [10,23,54].According to this model, the functional architecture of V1 is based on the geometrical structure inspired by the neural connectivity in V1.

Receptive profiles.
A simple cell is characterized by its receptive field, which is defined as the domain of the retina to which the simple cell is sensitive.Once a receptive field is stimulated, the corresponding retinal cells generates spikes which are transmitted to V1 simple cells via retinogeniculo-cortical paths.
The response function of each simple cell to a spike is called receptive profile (RP), and denoted by ψ (ζ,θ) : Q → C. It is basically the impulse response function of a V1 simple cell.Conceptually it is the measurement of the response of the corresponding V1 simple cell to a stimulus at a point1 ζ = (x, y) ∈ Q.
In this study we assume the response of simple cells to be linear.That is, for a given visual stimulus f : Q → R we assume the response of the simple cell at V1 coordinates (ζ, θ) to be This procedure defines the cortical stimulus a 0 : M → C associated with the image f .We note that receptive field models consisting of cascades of linear filters and static non-linearities, although not perfect, may be more adequate to account for responses to stimuli [20,55,56].Several mechanisms such as, e.g., response normalization, gain controls, cross-orientation suppression, or intra-cortical modulation, might intervene to change radically the shape of the profile.Therefore, the above static and linear model for the receptive profiles should be considered as a first approximation of the complex behavior of a real dynamic receptive profile, which cannot be perfectly described by static wavelet frames.
Regarding the form of the RP, in [10], a simplified basis of Gabor functions were proposed as good candidates for modelling the position-orientation sensitive receptive profiles for neuro-physiological reasons [57,58].This basis has then been extended to take into account additional features such as scale [54], velocity [33], and frequency-phase [39].On the other hand, Duits et al. [53] proposed so-called cake kernels as a good alternative to Gabor functions, and showed that cake kernels were adequate for obtaining simple cell output responses which were used to perform certain image processing tasks such as image enhancement and completion based on sR diffusion processes .
In this study, we employ cake kernels as the models of position-orientation RPs obtaining the initial simple cell output responses to an input image, and we use the V1 model geometry M to represent the output responses.We model the activity propagation along the neural connectivity by using the combination of a diffusion process based on the natural sub-Laplacian and a Wilson-Cowan type integro-differential system.

Horizontal connectivity and sub-Riemannian diffusion.
Neurons in V1 present two type of connections: local and lateral.Local connections connect neurons belonging to the same hypercolumn.On the other hand, lateral connections account for the connectivity between neurons belonging to different hypercolums, but along a specific direction.In the CPS model these are represented2 by the vector fields (2.2) The above observation yield to the modelling of the dynamic of the neuronal excitation {Z t } t≥0 starting from a neuron (ζ, θ) via the following stochastic differential equation where u t and v t are two one-dimensional independent Wiener processes.As a consequence, in [25] the cortical stimulus a 0 induced by a visual stimulus f 0 is assumed to evolve according to the Fokker-Planck equation Here, β > 0 is a constant encoding the unit coherency between the spatial and orientation dimensions.The operator L is the sub-Laplacian associated to the sub-Riemannian structure on M with orthonormal frame {X 1 , X 2 }, as presented in [10,25].It is worth mentioning that this operator is not elliptic, since {X 1 , X 2 } is not a basis of T M.However, span{X 1 , X 2 , [X 1 , X 2 ]} = T M. Hence, {X 1 , X 2 } satisfies the Hörmander condition and L is an hypoelliptic operator [59] which models the activity propagation between neurons in V1 as the diffusion concentrated to a neighborhood along the (horizontal) integral curves of X 1 and X 2 .
A direct consequence of hypoellipticity is the existence of a smooth kernel for (2.4).That is, there exists a function (t, ξ, ν An analytic expression for k t can be derived in terms of Mathieu functions [9,29].This expression is however cumbersome to manipulate, and it is usually more efficient to resort to different schemes for the numerical implementation of (2.4), see, e.g., Section 4.

Reconstruction on the retinal plane.
Activity propagation evolves the lifted visual stimulus in time.In order to obtain a meaningful result, which is represented on a 2-dim image plane, we have to transform back the evolved lifted image to the 2-dim image plane.We achieve this by using the projection given by where f : R 2 × (0, T ] → R and 0 < T < ∞ denote the processed image and the final time of the evolution, respectively.One easily checks that this formula yields f (•, 0) = f 0 under the assumption (2.7) 3. Describing neuronal activity via Wilson-Cowan-type models.In neurophysiological experiments, reliable neural responses to visual stimuli are generally observed at the neuronal population level: the information processing and the response produced are obtained by integrating the individual dynamics of the neurons interacting within the population.Modelling neuronal populations can be done via coupled differential systems (networks) consisting of a large number of equations, and the average behavior can in principle be used to represent the population behavior.This requires high computational power and the use of challenging analytical approaches due to the high dimension of the network.A different mesoscopic approach consists in considering the average network behavior as the number of neurons in the network is let to infinity.The asymptotic limit of the network can thus be written in terms of the probability distribution (density) of the state variables.This asymptotic limit is the so-called mean-field limit.It has been successfully used as a reference framework in several papers, see, e.g., [60][61][62] and will also be the approach considered in this work.
Here µ : Q → R is a smoothed version of the simple cell output response a 0 via a Gaussian filtering, while parameters λ > and M > 0 are fixed positive constants.Following the standard formulation of WC models studied, e.g., in [60,63] we have that the role of the time-independent external stimulus h : Q × [0, π) → R is played here by h(ξ) := λa 0 (ξ) + µ(ξ) while model parameters can be set as β := 1 + λ and ν := 1/2M .The function σ : R → [−1, 1] stands for a nonlinear saturation function, which we choose as the sigmoid: The connectivity kernel ω ξ models the interaction between neurons in M. Its definition should thus take into account the different type of interactions happening between connected neurons in V1, e.g. it should model at the same time both local and lateral connections via the sub-Riemannian diffusion described in Section 2.2.
In [44,45] the authors showed that (3.1) does not arise from a variational principle.That is, it there exists no energy function E : L 2 (M) → R such that (3.1) can be recast as the problem Under this formulation, stationary states a * of (3.1) are (local) minima of E.
The interest of considering an evolution model following a variational principle in the sense (3.3) is given by its connection with the optimization-based approaches considered in [64] to describe the efficient coding problem as an energy minimization problem which involves natural image statistics and biological constraints which force the final solution to show the least possible redundancy.Under this interpretation, the non-variational model (3.1) is suboptimal in reducing redundant information in visual stimuli, see [44, Section 2.1] for more details.

Local
Histogram Equalization (LHE) model.In order to build a model which complies with the efficient neural coding described above, in [44,45], the authors showed that (3.1) can be transformed into a variational problem by replacing the term σ(a(η, t)) with σ(a(ξ, t) − a(η, t)) for a suitable choice of the nonlinear sigmoid function σ, thus enforcing non-linear activations on local contrast rather than on local activity.The corresponding model reads: where σ(r) := −σ(r + 1/2), and σ as in (3.2).This model has been first introduced in [47] as a variational reformulation of the Local Histogram Equalization (LHE) procedure for RGB images.The corresponding energy E : L 2 (M) → R for which (3.3) holds is: where Σ : R → R is any (even) primitive function for σ.As it is clear from (3.5), the local histogram equalization properties of the model are due here to the activation averaging which is localized by the kernel ω ξ , which should thus be adapted to the natural geometry of M (see Section 3.3 for a more detailed discussion).

3.3.
A sub-Riemannian choice of the interaction kernel ω ξ .In (3.1) and (3.4), the geometric structure of the underlying space M is captured by the connectivity kernel ω ξ , which characterizes the activity propagation along neural connections in V1.In [44,45], simple 3-dimensional Gaussiantype kernels were considered.This choice was shown to be good enough in these works to reproduce a large number of contrast-and orientation-dependent Poggendorff-like illusions via the LHE model in (3.4), but not by the WC one (3.1).
Here, motivated by he discussion in Section 2.2, we study the effect of a more natural choice for the interaction kernel ω ξ , which we set as ω ξ (η) = k τ (ξ, η), where k τ : M × M → R is the sub-Riemannian heat kernel evaluated at time τ > 0. Indeed, 3-dimensional isotropic Gaussian kernels are obtained via the Euclidean heat equation are not coherent with the intrinsically anisotropic neuronal connectivity structure of V1.Recalling (2.5), this choice of ω ξ allows to rewrite the WC equations (3.1) as Using this formulation, the evaluation of the interaction term at point (ξ, t) ∈ M × (0, T ] can be done by solving the sub-Riemannian heat equation and let it evolve for a certain inner-time τ > 0. This avoids to deal directly with the explicit expression of k τ whose numerical implementation is very delicate, as explained, e.g., in [9].A similar simplification is not readily available for the LHE equation (3.4), due to the dependence on ξ of the integrand function.In this setting, we follow the discussion in [47] and replace the non-linearity σ by a polynomial approximation of sufficiently large order n.Namely, we look for a polynomial approximation of σ of the form σ(r) = c 0 + . . .+ c n r n , which allows us to write σ a(ξ, t) − a(η, t) This allows to approximate the interaction term in (3.4) as Finally, the resulting (approximated) sub-Riemannian LHE equation reads: 4. Discrete modelling and numerical realisation.In this Section, we report a detailed description of how models (sR-WC) and (sR-LHE) can be formulated in a complete discrete setting, providing, in particular, some insights on how the sub-Riemannian evolution can be realised.We further add a self-contained section regarding the gradient-descent algorithm used to perform the numerical experiments reported in Section 5, for more details see [44,46].4.1.Discrete modelling and lifting procedure via cake wavelets.The sub-Riemannian diffusion e τ L is discretised by a final time τ = m∆τ where m and ∆τ denote the number of iterations and the time-step, respectively.For N ∈ N + and ∆x, ∆y ∈ R + denoting the spatial sampling size, we then discretise the given gray-scale image function f 0 associated to the retinal stimulus on a uniform square spatial grid Q := {(x i , y j ) = (i∆x, j∆y) : i, j = 1, 2, . . ., N } ⊂ R 2 and denote, for each i, j = 1, 2, . . ., N , the brightness value at point ζ i,j := (x i , y j ) ∈ Q by (4.1) As far as the orientation sampling is concerned, we use a uniform orientation grid with points Θ := {θ k := k∆θ, k = 1, . . ., K}, K ∈ N + and ∆θ = π/K .We can then define the discrete version of the simple cell response a 0 (x i , y j , θ k ) to the visual stimulus located at ζ i,j ∈ Q with local orientation θ k ∈ Θ at time t = 0 of the evolution as where L : Q → Q × Θ is the lifting operation to be defined.
To do so, we consider in the following the image lifting procedure based on cake kernels introduced in [53] and used, e.g., in [32,44,46].We write the cake kernel centered at ζ i,j and rotated by θ k as where , m ∈ {1, 2, . . ., N }.We can then write the lifting operation applied to the initial image f 0 for all ζ i,j ∈ Q and θ k ∈ Θ as: Finally, for P ∈ N + we consider a time-discretisation of the interval (0, T ] at time nodes T := {t p := p∆t, p = 1, . . .P }, P ∈ N + with ∆t := T /P . The resulting fully-discretised neuronal activation at ζ i,j = (x i , y j ) ∈ Q, θ k ∈ Θ and t p ∈ T will be thus denoted by: (4.5) A p [i, j, k] = a(ζ i,j , θ k , t p ).
4.2.Sub-Riemannian heat diffusion.Let g : M → R be a given cortical stimulus, and denote and set G[i, j, k] = g(ξ i,j , θ k ).In this section we describe how to compute The main difficulty here is due the degeneracy arising from the anisotropy of the sub-Laplacian.Indeed, developing the computations in (2.4), we have In particular, it is straightforward to deduce that the eigenvalues of are (0, β 2 , 1).
The discretisation of such anisotropic operators can be done in several ways, see for example [29,30,39,65].In our implementation, we follow the method presented in [26] which is tailored around the group structure of SE(2), the universal cover of M, and based on the non-commutative Fourier transform, see also [8].
It is convenient to assume for the following discussion ∆x = ∆y = √ N and ∆θ = π/K.The "semi-discretised" sub-Laplacian L K can be defined by where by Λ K we denote the central difference operator discretising the derivatives along the θ direction, i.e. the operator (4.9) The operator D is the diagonal operator defined by (4.10) The full discretisation is then achieved by discretising the spatial derivatives as Under the discretisation L K of L defined in (4.8), we now resort to Fourier methods to compute efficiently the solution of the sub-Riemannian heat equation (4.13) In particular, let Ĝ[r, s, k] be the discrete Fourier transform (DFT) of G w.r.t. the variables i, j, i.e.
A straightforward computation shows that Hence, (4.13) is mapped by the DFT to the following completely decoupled system of N 2 ordinary linear differential equations on C K : (4.16) which can be solved efficiently through a variety of standard numerical schemes.We chose the semiimplicit Crank-Nicolson method [66] for its good stability properties.Let us remark that the operator at the r.h.s. of the above equations are periodic tridiagonal matrices, i.e., tridiagonal matrices with additional non-zero values at positions (1, K) and (K, 1).Thus, the linear system appearing at each step of the Crank-Nicolson method can be solved in linear time w.r.t.K via a variation of the Thomas algorithm.
The desired solution exp τ G can be then simply recovered by applying the inverse DFT to the solution of (4.16) at time τ .

Discretisation via gradient descent.
We follow [44,45,47] and discretise both models (sR-WC) and (sR-LHE) via a simple explicit gradient descent scheme.Denoting the discretised version of the local mean average µ(ξ) appearing in the models by U [i, j, k] = µ(i∆x, y∆j, k∆θ), we have that the the time stepping reads for all p ≥ 1 (4.17 where SA p−1 is defined depending on the model by: (4.18) with C ,p−1 being the discretised version of the coefficient C in (3.6) at time t p−1 .A sufficient condition on the time-step ∆t guaranteeing the convergence of the numerical scheme (4.17) is ∆t ≤ 1/(1 + λ) (see [47]).4.4.Pseudocode.Our algorithmic procedure consists of three main numerical sub-steps.The first one is the lifting of the two dimensional input image f 0 to the space M via (4.4).The second one is the Fourier-based procedure described in Section 4.2 to compute the sub-Riemannian diffusion (4.6) which can be used as kernel to describe the neuronal interactions along the horizontal connection.This step is intrinsically linked to the last iterative procedure, based on computing the gradient descent update (4.17)-(4.18)describing the evolution of neuronal activity in the cortical framework both for the (sR-WC) and (sR-LHE).
We report the simplified pseudo-code below.
5. Numerical Experiments.In this section we present the results obtained by applying models (sR-WC), (sR-LHE) via Algorithm 1 to two Poggendorf-type illusions reported in Figure 3.Our results are compared to the ones obtained by applying the corresponding WC and LHE 3-dimensional models with a 3D-Gaussian kernel as described in [44,45].The objective of the following experiments is to understand whether the output produced by applying (sR-WC) and (sR-LHE) to the images in Figure 3 agrees with the illusory effects perceived.Since the quantitative assessment of the strength of these effects is a challenging problem, the outputs of Algorithm 1 have too be evaluated by visual inspection.Namely, for each output, we consider whether the continuation of a fixed black stripe on one side of a central bar connects with a segment on the other side.Differently from inpainting-type problems, we stress that for these problems the objective is to replicate the perceived wrong alignments due to contrast and orientation effects rather than its collinear prosecution and/or to investigate when both type of completions can be reproduced.
Testing data: Poggendorff-type illusions.We test the (sR-WC) and (sR-LHE) models on a grayscale version of the Poggendorff illusion in Figure 2 and on its modification reported in Figure 3b where the background is constituted by a grating pattern: in this case the perceived bias depends also on the contrast between the central surface and the background lines.Parameters.Images in Figure 3 have size N × N pixels, with N = 200.The lifting procedure to the space of positions and orientations is obtained by discretising [0, π) into K = 16 orientations (this is in agreement with the standard range of 12-18 orientations typically considered to be relevant in literature [67,68]).The relevant cake wavelets are then computed following [32], setting the frequency band bw = 5 for all experiments.The scaling parameter β appearing in (2.4) 2), and the parameter M appearing in (sR-WC), (sR-LHE) is set to M = 1.
Parameters varying from test to test are: the slope α > 0 of the sigmoid functions σ in (3.2) and σ, the fidelity weight λ > 0, the variance of the 2D Gaussian filtering σ µ use to compute the local mean average µ in (sR-WC) and (sR-LHE), the gradient descent time-step ∆t, the time step ∆τ and the final time τ used to compute the sub-Riemannian heat diffusion e τ L .

Poggendorff gratings.
In Figure 4 we report the results obtained by applying (sR-WC) to the Poggendorff grating image in Figure 3b.We compare them with the ones obtained by the cortical-inspired WC model considered [44,45] where the sR heat-kernel is an isotropic 3D Gaussian which are reported in Figure 4a.In Figure 4b, we observe that the sR diffusion encoded in (sR-WC) favours the propagation of the grating throughout the central gray bar so that the resultant image agrees with our perception of misalignment.We stress that such illusion could not be reproduced via the cortical-inspired isotropic WC model proposed in [44,45].The use of the appropriate sub-Laplacian diffusion is thus crucial in this example to replicate the illusion.We further report in Figure 5 the result obtained by applying (sR-LHE) on the same image.We observe that in this case both (sR-LHE) model and the LHE cortical model introduced in [44,45] reproduce the illusion.
Note that both (sR-WC) and (sR-LHE) further preserve fidelity w.r.t. the given image outside the target region, which is not the case in the LHE cortical model presented in [44,45].
5.2.Dependence on parameters: inpainting vs. perceptual completion.The capability of the (sR-LHE) model to reproduce visual misperceptions depends on the chosen parameters.This fact was already observed in [45] for the cortical-inspired LHE model proposed therein endowed by a standard Gaussian filtering.There, LHE was shown to reproduce illusory phenomena only in the case where the chosen standard deviation of the Gaussian filter was set to be large enough (w.r.t. the overall size of the image).On the contrary, the LHE model was shown to perform geometrical completion (inpainting) for small values of the standard deviation.Roughly speaking, this corresponds to the fact that perceptual phenomena -such as geometrical optical illusions -can be modelled only when the interaction kernel is wide enough for the information to cross the central gray line.This is in agreement with psycho-physical experiences in [17] where the width of the central missing part of the Poggendorff illusion is shown to be directly correlated with the intensity of the illusion.
In the case under consideration here, the parameter encoding the width of the interaction kernel is the final time τ of the sub-Riemannian diffusion used to model the activity propagation along neural connections.To support this observation, in Figure 6 we show that the completion obtained via (sR-LHE) shifts from a geometrical one (inpainting), where τ is small, to a perceptual one, for τ sufficiently big.
As far as (sR-WC) model is concerned, we observed that despite the improved capability of replicating the Poggendorf gratings, the transition from perceptual completion to inpainting could not be reproduced.In agreement with the efficient representation principle, this supports the idea that visual perceptual phenomena are better encoded by variational models as (sR-LHE) than by non-variational ones as (sR-WC).

Poggendorff illusion.
In Figure 7 we report the results obtained by applying LHE methods to the standard Poggendorff illusion in Figure Figure 3a.In particular, in Figure 7a we show the result obtained via the LHE method of [44,45], while in Figure 7b we show the result obtained via (sR-LHE), with two close-ups in Figures 7c and 7d showing a normalized detail of the central region onto the set of values [0, 1].As shown by these preliminary examples, the prosecutions computed by both (LHE) models agree with our perception as the reconstructed connection in the target region links the two misaligned segments, while somehow 'stopping' the connection of the collinear one.
This phenomenon as well as a more detailed study on how the choice of the parameters used to generate Figure 3a (such as the incidence angle, the width of the central gray bar, the distance between lines), in a similar spirit to [69] where psysho-physicis experiments were performed on analogous images, is an interesting topic for future research.(resp.Figure 7a).

Conclusion.
In this work we presented a sub-Rimannian version of the Local Histogram Equalization mean-field model previously studied in [44,45] and here denoted by (sR-LHE).The model considered is a natural extension of existing ones where the kernel used to model neural interactions was simply chosen to be 3D Gaussian kernel, while in (sR-LHE) this is chosen as the sub-Riemannian kernel of the heat equation formulated in the space of positions and orientations given by the primary visual cortex (V1).A numerical procedure based on Fourier expansions is described to compute such evolution efficiently and in a stable way and a gradient-descent algorithm is used for the numerical discretisation of the model We tested the (sR-LHE) model on orientation-dependent Poggendorff-type illusions and showed that (i) in presence of a sufficiently wide interaction kernel, model (sR-LHE) is capable to reproduce the perceptual misalignments perceived, in agreement with previous work (see Figures 5 and 7), (ii) when the interaction kernel is too narrow, (sr-LHE) favors a geometric-type completion (inpainting) of the illusion (see Figure 6) due to the limited amount of diffusion considered.
We also considered (sR-WC), a similar model obtained by using the sub-Riemannian interaction kernel in the standard Wilson-Cowan equations.We showed that the introduction of such corticalbased kernel improves the capability of WC-type models of reproducing Poggendorff-type illusions, in comparison to the analogous results reported [44,45] where the cortical version of WC with a standard 3D Gaussian kernel was shown to fail to replicate the illusion.
Finally, we stress that, in agreement with the standard range of 12-18 orientations typically considered to be relevant in literature [67,68], all the aforementioned results have been obtained by considering K = 16 orientations.The LHE and WC models previously proposed were unable to obtain meaningful results with less that K = 30 orientations.
(a) Association fields (b) Projected horizontal integral curves

Fig. 3 :
Fig. 3: Gray-scale Poggendorff-type illusions.Figure (3a) is the standard 200×200 Poggendorff illusion with a 30 pixel-wide central and an incidence angle of π/3 drawn by the black lines with the central bar .Figure (3b) is a variation of the classical Poggendorff illusion where a further background grating is present.