Multi-Frequency Image Completion via a Biologically-Inspired Sub-Riemannian Model with Frequency and Phase

We present a novel cortically-inspired image completion algorithm. It uses five-dimensional sub-Riemannian cortical geometry, modeling the orientation, spatial frequency and phase-selective behavior of the cells in the visual cortex. The algorithm extracts the orientation, frequency and phase information existing in a given two-dimensional corrupted input image via a Gabor transform and represents those values in terms of cortical cell output responses in the model geometry. Then, it performs completion via a diffusion concentrated in a neighborhood along the neural connections within the model geometry. The diffusion models the activity propagation integrating orientation, frequency and phase features along the neural connections. Finally, the algorithm transforms the diffused and completed output responses back to the two-dimensional image plane.

1. Introduction.Visual perception has drawn attention of experts from philosophy, psychology, neuroscience, as well as the attention of mathematicians and physicist working on perceptual modeling.The question of how we perceive was studied by Edmund Husserl in his pioneering philosophical texts in phenomenology [1][2][3].In the side of psychology, we can think of the well known Berlin school of experimental psychology, Gestalt psychology school, which formulated what is known today as Gestalt psychology of perception [4][5][6].
Gestalt psychology is a theory which attempts to provide the principles giving rise to a meaningful global perception of a visual scene by starting from the local properties of the objects within the scene.The main idea of Gestalt psychology is that the mind constructs the global whole by rather grouping similar fragments than purely summing the fragments as if they were indifferent.In terms of visual perception, those similar fragments can be thought of as the point stimuli with the same (or closely) valued features of the same type.In the present paper, we consider orientation, spatial frequency and phase as features.We extract these features from a given two dimensional grayscale input image which is partially occluded and reveal those occluded parts via an integration of the extracted features.By doing so, we follow the Gestalt principle which is called the law of good continuity to provide a feature integration mechanism reconstructing the occluded parts in the image.
The law of good continuity states that we group rather aligned pieces than those with sharp abrupt directional changes when we perceive an object as a whole which is formed by fragments, see Figure 1.This law was studied by Field, Hayes and Hess [7] from a psychophysical point of view.They worked on how the visual system captures aligned fragments constituting lines or cocircular curves on a two dimensional background of randomly oriented fragments; see Figure 4.They noticed that the aligned patterns captured by the visual system were locally overlapping with what they called association fields, which are shown in Figure 3.Those fields provide a geometric characterization of the law of good continuity and they can be interpreted as the psychophysical representations of the neural connections which are biologically implemented in the visual cortex.
The visual system is capable of perceiving a lacunar curve as complete by capturing the whole curve pattern underlying the lacunar curve; see Figure 1.This is due to a general phenomena which is called perceptual completion.This phenomena provides the perception of contours and figures which are actually not present in the visual stimulus.We call such contours subjective contours.
Kanizsa [8,9] explains two categories of completion: modal completion and amodal completion.The first one refers to the completion following the modality of vision.There is no direct stimulus corresponding to the object, yet, the perceived object is indistinguishable from the real stimuli and we perceive it as a whole.In the second category, the completion does not make use of the modality of vision.In other words, the stimulus corresponding to the object is partial but we still perceive the object as complete.In this case, the object is recognized as a whole although only some specific Fig. 1: An example of the law of good continuity.We capture the curve on the bottom left as the curve underlying the aligned fragments on the top left, we do not capture any curve underlying the fragments on the top right due to the abruptly changing orientation angles of the fragments.fragments of the object evoke our sensory receptors.Examples of those two categories of completion resulting in subjective contours are Kanizsa triangle and Ehrenstein's illusion, which are illustrated in Figure 2. We focus rather on amodal completion in the present work, and provide an approach which reconstructs the occluded parts in a compatible way with the law of good continuity and association fields, where those two notions are considered in an extended way based on position, orientation, frequency and phase alignment.Fig. 2: Left: Kanizsa triangle.There is no direct stimulus, yet we perceive a white triangle on top of the rest (modal completion).We perceive another triangular whose border is marked by the black lines on the bottom layer (amodal completion).Right: Ehrenstein illusion.We perceive a white disk around the center despite the absence of a direct stimulus (modal completion).We recognize that each vertical, horizontal or diagonal line fragment comprises a whole line which is occluded by the white disk (amodal completion).
The primary visual cortex (V1) is the main area of the cerebral cortex which is responsible for the first step processing of visual input so that a proper visual perception is achieved at a higher perceptual level.V1 contains a particular family of neurons, simple cells.These neurons are locally sensitive to visual features, such as orientation [10][11][12][13], spatial frequency [14][15][16][17][18][19][20][21], phase [22][23][24][25], scale [26] and ocular dominance [16,27,28].The simple cells are organized in a hypercolumnar architecture, which was first discovered by Hubel and Wiesel [29].In this architecture, a hypercolumn is assigned to each point (x, y) of the retinal plane M R 2 , and the hypercolumn contains all the simple cells sensitive to a particular value of the same feature type.The simple cells are able to locally exctract the feature Fig. 3: Top: association fields aligned with a horizontal patch as shown by Field, Hayes and Hess [7].Bottom: solid curves represent the association fields between strongly associated fragments, and the dashed ones imply no fields between weakly associated fragments.A simple cell is identified by its receptive field, which is defined as the domain of the retina to which the cell is sensitive and connected through the retino-geniculo-cortical paths.Once a receptive field is stimulated, it evokes a spike transmitted to the corresponding simple cells.Each one of those simple cells produces a response to the spike.This response is what is known as receptive profile.
Hoffman proposed modelling the hypercolumnar architecture of V1 in terms of a contact bundle structure [30,31].This framework was followed by Petitot and Tondut, where they improved the contact bundle structure and proposed a boundary completion method within the corresponding contact geometry [32].Moreover, they geometrically interpreted the association fields as the integration along the vector fields generating the contact geometry.This setting was developed further by Citti and Sarti [33] to a framework in which they introduced a group based approach to study the geometric modeling of V1 hypercolumnar architecture and the functional connectivity.They used Gabor function as the receptive profile model and proposed the sub-Riemannian geometry of the group of rotations and translations (SE(2)) as the V1 model geometry.This framework was extended to higher dimensional geometries where scale [34], and velocity [35,36] were taken into account.Various biologically-inspired models for optical illusions [37][38][39][40][41] and orientation preference maps [42], as well as frameworks for image processing [43][44][45][46][47][48][49], for pattern recognition [50] and for medical applications [51,52], were proposed in the sub-Riemannian geometry of SE (2).
The model presented in [33] is an abstract geometric description of the orientation-sensitive V1 hypercolumnar architecture reported by Hubel and Wiesel [10][11][12] .This description provides a good phenomenological approximation of the biologically implemented V1 neural connections which were reported by Bosking et al. [53].In this model framework, the projections of a particular family of curves onto the two dimensional image plane provide good approximations of the association fields.In other words, these curves model the V1 neural connections; see Figure 5.They are called horizontal integral curves and they are obtained by integration along the vector fields generating the SE(2) sub-Riemannian geometry.For this reason, the approach considered by Citti, Petitot and Sarti and used in our present paper is referred to as biologically-inspired.[7,33].
The sub-Riemannian model geometry proposed in [33] was extended in a recent work to multifrequency and multi-phase setting [54].The extended model corresponds to a natural geometry which is derived [54,55] from one of the very first perceptual mechanisms of vision: receptive profile.It is different from the classical approach, in which a suitable geometry is assigned to the neural responses represented in terms of receptive profiles.This extended model takes advantage of orientations, spatial frequencies and phases in a given 2D input image to encode the visual information.This is not the case in the sub-Riemannian model proposed in [33], which uses the SE(2) geometry, and in which only orientation can be represented as visual feature.The extended model geometry was applied to enhancement on images in which several spatial frequencies were equally present [54], such as texture images.In the present work, we employ the same extended multi-frequency sub-Riemannian setting presented in [54], and propose an image completion algorithm within, which is aimed for grayscale texture images containing multiple spatial frequencies.
In Section 2, we explain the model framework which was introduced in [54] and its relation to our completion algorithm.In Section 3, we present a specific family of integral curves defined in the model geometry, which are the models of the neural connections in V1.Then in Section 5, we introduce our completion algorithm, its discrete scheme and the corresponding pseudocode.Finally, we provide our simulation results and their comparison to some results obtained previously in [56] and [49].At the end, we give the main conclusions and some perspectives for the related future research.
2. Model framework.In this section, we will explain the cortical model geometry in which our completion algorithm is defined.This model framework is composed by two mechanisms: feature value extraction and neural connectivity [54].

Feature value extraction.
Each simple cell is sensitive to a specific part of the retina, which is called receptive field.Once the receptive field is stimulated by visual stimulus, the retinal cells in the receptive field produce spikes which are transmitted through retino-geniculo-cortical pathways to the related simple cells in V1.Each simple cell generates a response to those spikes, which is the receptive profile corresponding to the simple cell.In other words, receptive profile is the impulse response function of the simple cell.The simple cell receptive profile which is sensitive to the stimulus located at q ∈ M on the image plane M and selective to the set of feature values z ∈ S1 × R + × S 1 is denoted by Ψ (q,z) , where q = (q 1 , q 2 ) and z = (θ, f, φ) together denote a fixed point Here Q represents the five dimensional sub-Riemannian V1 model geometry.
Simple cell receptive profiles can be modeled in terms of Gabor functions [33,41,54,57].In the orientation, frequency and phase selective model framework, receptive profile of a simple cell is a Gabor function of the following type: , where f > 0 represents the spatial frequency 1 , r = (r 1 , r 2 ) = (−f sin θ, f cos θ) and σ > 0 is the scale of the localizing Gaussian.The complex exponential is the wave content and it is the main component capturing the orientation, frequency and phase information of the objects in a given two dimensional image.The second exponential is a Gaussian which spatially localizes the receptive profile around the point (q 1 , q 2 ).Frequency f determines how many wave peaks are found within the localizing Gaussian window; see Figure 6.As the number of wave peaks increases, the Gabor function can detect higher frequencies.Orientation θ is the orientation angle to which the simple cell receptive profile is sensitive.Parameter φ is the reference phase and it creates a phase shift in the waves of the Gabor function as it changes.
We disregard the coordinate map between the image plane and the retinal surface, as well as the retino-cortical map between the retinal surface and the V1 surface.We assume that the image plane is identically mapped to the retinal and cortical surfaces.We assume the responses of the simple cells to be linear and we compute the output response of a simple cell located at (q, z) ∈ Q to a given two-dimensional grayscale image I : M → [0, 1] via the convolution with the Gabor filter: We apply the convolution for every feature value z and for every point q.Consequently, we obtain the output responses of all receptive profiles corresponding to the V1 simple cells.We will call this set of output responses sometimes, lifted image; and the Gabor transform, lifting.It is equivalent to the result of a multi-frequency Gabor transform applied to the given two dimensional input image.Those output responses are the representations of the feature values in the five dimensional V1 model geometry Q.
In general, static receptive profile models based on linear filter banks and static nonlinearities [33,54,[58][59][60][61] provide good responses to simple stimuli.However, their responses to complicated stimuli, such as natural images, are approximate up to a certain level.Several mechanisms such as response normalisation, gain controls, cross-orientation suppression and intra-cortical modulation, can result in radical changes in the receptive profile shape.Therefore, the aforementioned Gabor filter bank model for the receptive profiles should be considered as a first approximation of highly complex real dynamic receptive profile.We employ all frequency components of the Gabor transform during the lifting.Therefore, exact inverse Gabor transform is valid and we use it to obtain the corresponding two dimensional image to the output responses: with Ψ denoting the complex conjugate of the corresponding Gabor function Ψ.

Horizontal connectivity.
Lifting provides the output responses, which are complex valued functions in the five dimensional model geometry Q.Each output response encodes the feature values corresponding to the orientation, frequency and phase of a pixel defined on the two dimensional image plane.The output responses, hence the simple cells, are isolated from each other once lifting from the image plane M to the model geometry Q takes place.Therefore, the model geometry Q should be endowed with an integration mechanism which provides the activity propagation, therefore the interactiviy, between the simple cells.The activity propagation provides an integrated form of the local feature vectors associated to the lifted image.This propagation is concentrated along a specific family of integral curves, horizontal integral curves, corresponding to the model geometry Q.The horizontal integral curves can be thought of as the models of the long range lateral connections, which connect the simple cells residing in different hypercolumns but selective to the same (or close) feature values.
We may associate the following differential one-form to each receptive profile described by (2.1): The one-form induces naturally the horizontal vector fields corresponding to the model geometry Q.
The horizontal vector fields are formally defined as the elements of where T Q denotes the tangent bundle of Q.The horizontal vector fields corresponding to Q are found from (2.5) as Those horizontal vector fields endow Q with a sub-Riemannian structure as explained in [54].They span the horizontal tangent space H (q,z) Q at each (q, z) ∈ Q.The horizontal tangent space can be thought of as the analogue of the Euclidean tangent space.In other words, differential operators such as gradient and Laplacian are defined in terms of the horizontal vector fields in the sub-Riemannian geometry.We remark that the differential operators are degenerate since the horizontal tangent space is spanned by four vector fields although it corresponds to Q, which is a five dimensional geometry.The horizontal integral curves are defined as the integrated curves along the horizontal vector fields given in (2.6).Albeit the degenerate character of the horizontal tangent space H (q,z) Q, they provide the full connectivity in Q due to that X 1 and X 2 do not commute as we will see below.
Nonzero commutators of the horizontal vector fields are found as (2.7) The horizontal vector fields are bracket generating since (2.8) for all (q, z) ∈ Q where T (q,z) Q denotes the tangent space at (q, z) ∈ Q.Indeed, (2.8) shows that the horizontal vector fields fulfill Hörmander condition [62].Consequently, they provide the connectivity of any two points in Q through the horizontal integral curves due to the Chow-Rashevskii theorem [63][64][65].This connectivity property has particular importance since it assures that any two points in the V1 sub-Riemannian model geometry Q can be connected via the horizontal integral curves, which are the models of the neural connections implemented biologically in V1 and which are close approximations of the association fields at the psychophysical level.

Horizontal integral curves.
The association fields were proposed to be modeled by the horizontal integral curves of SE(2) in the classical orientation sensitive framework [33].A similar line of thought was followed in [54] and it was proposed to employ the horizontal integral curves corresponding to the the five dimensional sub-Riemannian geometry Q as the cortical counterparts of the association fields.The projection of those horizontal integral curves of Q are the same as the projections of the horizontal integral curves of SE (2), which are shown in Figure 5.It was conjectured in [54] that the horizontal integral curves of Q coincide with the long range lateral connections between orientation, frequency and phase selective simple cells in V1.
Let us denote a time interval by I = [0, T ] with 0 < T < ∞ and consider a horizontal integral curve (q 1 , q 2 , θ, f, φ) = γ : I → M associated to the horizontal vector fields given in (2.6).We denote the initial point of γ by α = (q 1 , q2 , θ, f , φ) and its velocity by γ .At each time instant t ∈ I, the velocity is written as a vector γ (t) ∈ span(X 1 , X 2 , X 3 , X 4 ) γ(t) at γ(t) = (q 1 (t), q 2 (t), θ(t), f (t), φ(t)) ∈ Q.One way to compute the horizontal integral curves starting from the initial point α is to solve the following ODE system for all t ∈ I: where c 2,3,4 denote the coefficients.In the activity propagation machinery which we propose here, the propagation is concentrated along a neighbourhood of the horizontal integral curves with constant coefficients c 2,3,4 in the cortical space Q [54]; see Figures 7 and 8.However, the coefficients need not necessarily be constants in the generic framework of horizontal integral curves.Fig. 7: A horizontal integral curve along the vector field X 1 + X 2 .It represents an orientation fiber once f and φ are fixed.The tangent planes spanned by X 1 , X 2 (left) and X 3 , X 4 (right) are shown at six points on the curve.4. Sub-Riemannian diffusion in the cortical space.The activity propagation was proposed to be modeled in terms of a sub-Riemannian diffusion process in the classical orientation sensitive SE(2) framework [33].This sub-Riemannian diffusion process can be interpreted as the model of the interacting neural dynamics defined in terms of the corresponding horizontal vector fields and it was applied to several biologically-inspired image processing algorithms [43][44][45]49].
We follow a similar approach as in [33] and define a sub-Riemannian diffusion procedure in the five dimensional model geometry Q.We denote by Σ ⊂ Q the subspace in which all output responses O I s are defined.These are the output responses obtained as the lifting of the two dimensional input image I.We will denote by Π ⊂ Σ the subspace in which we find only the output responses corresponding to the part to be completed in the image I.The sub-Riemannian diffusion operator for all (q, z) ∈ Σ is defined as where β 2,3,4 are coefficients assuring the unit coherency in spatial, orientation, frequency and phase dimensions.The sub-Riemannian diffusion is described for all (q, z) ∈ Σ by (4.1) ∂ t u(q, z, t) = Lu(q, z, t), u(q, z, 0) = O I (q, z), t ∈ (0, T ], 0 < T < ∞, with the corrupted region marked by the boundary conditions given by u(q, z, t) = O I (q, z) for all (q, z) ∈ Σ − Π and t ∈ (0, T ].Here T denotes a sufficiently large final time and u : Q × [0, T ] → C stands for the output responses evolving in time.We denote the number of orientations, frequencies and phases of the N × N image I by K, L and M , respectively.Then, the coefficients are found as: We note that span(X 1 , X 2 ) and span(X 3 , X 4 ) define two subspaces of the horizontal tangent space H (q,z) Q at each point (q, z) ∈ Q.This allows us to decompose the sub-Riemannian horizontal tangent space into two components of which each one is defined by the vector fields of those two subspaces of T (q,z) Q.Consequently, we may approximate the sub-Riemannian diffusion described by (4.1) as a diffusion applied in each frequency and phase channel separately.More precisely, we apply the classical sub-Riemannian diffusion procedure defined in SE(2) [33,45,49] for each frequency and phase channel separately by using the SE(2) sub-Riemannian diffusion operator where β 2 provides the unit coherency between the diffusion components in spatial and orientation dimensions.The approximate sub-Riemannian procedure is described by with the boundary conditions given as u(q, z, t) = O I (q, z).The advantage of such approximation is that we perform the diffusion procedure in L three dimensional spaces instead of in a five dimensional geometry by still taking advantage of the frequency information extracted from the input image.We implement both the exact and the approximate procedures by using a simple forward Euler scheme in which the derivatives are implemented in terms of B-spline interpolated central finite differences.

Algorithm.
Our algorithm is based on three steps.Once a two dimensional grayscale input image is given, the first step is lifting the image via the convolution with Gabor functions as given in (2.2).This provides the output responses which encode the orientation, frequency, phase values and which are represented in the five dimensional model geometry.The second step is the sub-Riemannian diffusion described by (4.1), or by (4.3) if the approximate setting is used.We integrate in time (4.1) for the exact framework, and (4.3) for the approximate framework, by applying it on the output responses via iterating it with the time step ∆t until the final time T , at which the steady-state is reached, i.e., ∂ t u = 0. We employ an explicit method where we use B-splined interpolated finite differences to implement the horizontal vector fields given in (2.6).The final step is to transform back the evolved output responses to the two dimensional image plane.This is achived by the inverse Gabor transform given by (2.3).

Discretization of the output responses.
We employ a uniform spatial grid to discretize the image plane such that (5.1) where i, j ∈ {1, 2, . . ., N } with N denoting the image size and ∆x, ∆y ∈ R + denoting the pixel width.
In our case, the input images are square images and ∆x = ∆y = 1 in terms of pixel unit.We denote the number of samples in the orientation dimension by K, in the frequency dimension by L and in the phase dimension by M .We express the distance between two adjacent samples in the orientation dimension with ∆θ, in the frequency dimension with ∆f and in the phase dimension with ∆s.Discretized output response O I (q 1,i , q 2,j , θ k , f l , φ m ) given to I[i, j] on uniform orientation, frequency and phase grids with points θ k = k ∆θ, f l = l ∆f and φ m = m ∆s (k ∈ {1, 2, . . ., K}, l ∈ {1, 2, . . ., L}, m ∈ {1, 2, . . ., M }) is denoted by where q 1,i = i∆x and q 2,j = j∆y.The Gabor function given by (2.1) is written in the discrete setting as follows: (5.3) where ĩ, j ∈ {1, 2, . . ., N }, k ∈ {1, 2, . . ., K}, ñ ∈ {1, 2, . . ., M }, for each orientation θ k , frequency f l and phase φ m .We fix s ñ = 0 and express the discretized cell response obtained from the input image I[i, j] via the lifting described by the discrete Gabor transform as follows: (5.4) We discretize the time interval by V ∈ N + samples, and denote it by h v .Here h v is the time instant h v = v ∆t with ∆t satisfying T = V ∆t and v ∈ {1, 2, . . ., V }.Discretized evolving output response is written as Finally, the discrete inverse transform applied via a normalized kernel Ψ on the evolved output responses until the final time T gives the completed two dimensional image I T , which is found as follows: (5.6) 5.2.Explicit scheme with finite differences.Here we provide the explicit numerical scheme which we employ to iterate the exact and approximate frameworks given in (4.1) and (4.3), respectively.Our motivation for choosing an explicit scheme rather than an implicit scheme is that the latter requires large memory and computational power in our multidimensional framework.
We follow [44,66] to implement the horizontal vector fields given in (2.6) via B-spline interpolated central finite differences.The interpolation takes place on a uniform grid.It is needed since the horizontal vectors are not always aligned with the spatial grid point samples.B-spline interpolation is based on the coefficients b(i, j) The coefficients are determined such that the spline polynomial sp(x, y), together with the B-spline basis functions ρ(x − i, y − j), coincides with the horizontal derivatives at the grid points.For example, the condition sp(i∆x, j∆y) = X 1 O I [i, j, k, l, m] must be satisfied once the correct coefficients b are determined; see [44,66,67] for more details.
We define e k ξ :=(∆x cos(θ k ), ∆y sin(θ k )), e k η :=(−∆x sin(θ k ), ∆y cos(θ k )), (5.8) whose illustrations corresponding to ∆x = ∆y = 1 case are given in Figure 9.We abuse the notation to denote the evolving output responses: (5.9) Fig. 9: Illustration of the vectors e k ξ and e k η at (0, 0) with ∆x = ∆ y = 1.The figure was modified and adapted from [44]. and then we write the central finite differences of the second order horizontal derivatives as (5.10) Finally we write the discretized numerical iteration for (4.1) and (4.3) as follows: (5.11)where L represents the discretized version of either L or L, depending on which one between exact and approximate frameworks is chosen.The discretization is achieved by replacing the second order horizontal derivatives with their discrete versions given in (5.10).
Inverse transform the evolved output responses U T to obtain the completed image I T via (5.6).
The detailed Matlab ® code can be found in the following link: https://drive.google.com/file/d/1YET47AxYo5 FQs3ObREAqvQp2KjuKh7n/view (accessible starting from 29 October 2021).
6. Numerical Experiments.We choose σ = 2 as the scale of the localizing Gaussian given in (2.1) for all experiments.We use 128 × 128 grayscale images in the experiments related to Figures 10-18 and 256 × 256 grayscale image in the experiment corresponding to Figure 19.We use θ ∈ {0, π 32 , 2π 32 , 3π 32 , . . ., 31π 32 }, φ ∈ {0, π 8 , 2π 8 , 3π 8 , 4π 8 } and ∆t = 0.1 in all the experiments.Our first results are obtained by using an artificial test image as shown in Figure 10.In the test image, we have arcs of circles centered at the top left corner and with different radii, where a sinusoidal function whose frequency increases linearly as the radius increases.The arcs with sinusoidal pattern are occluded by zero valued arcs belonging to the circles which are centered at the top right corner.We apply our completion procedure with T = 10 and f ∈ {2, 2.5455, 3.0909, . . ., 8}.We observe that both the approximate and exact frameworks show a similar performance and they provide proper completion.
In Figure 11, we present two cases of our completion algorithm applied to a real texture image which is occluded by arcs corresponding to circles with different radii.We perform our completion procedure with f ∈ {2.00, 2.55, 3.09, 3.64, 4.18, 4.73, 5.27, 5.82, 6.36, 6.91, 7.45, 8.00}.Here T = 10 for the upper row and T = 5 for the bottom row.We use the same parameters but now T = 15 in Figure 12 on the same texture image but now occluded by vertical and horizontal bars.In such set of test images, the challenge is that the bars cross each other.Our completion algorithms provide proper completion of such crossing areas.We observe that both the approximate and the exact algorithms are able to complete in all cases the occluded parts in Figures 11 and 12.
In order to see the impact of using multiple frequencies in the completion process via the sub-Riemannian diffusion, we performed the sub-Riemannian diffusion by using only one single frequency; see Figures 13 and 14.In other words, f is fixed to a constant both in the lifting procedure and in the sub-Riemannian diffusion.In Figure 13, the same corrupted image given in the bottom row of Figure 11 was used.The paramaters are the same as in the bottom row results of Figure 11.We used f = 2, 2.55, 3.64 from left to right.We observe that the performance of single frequency framework is visibly lower than the multi-frequency framework whose results are presented in Figure 11 due to the presence of a wide range of spatial frequencies in the image to be completed.We observe the same also in the single frequency results presented in Figure 14, which are obtained by using the bottom row corrupted image given in Figure 12.Finally, we note that the inverse Gabor transform in those single frequency experiments is not well defined since the information corresponding to the other frequency components in the input image is lost in the lifting procedure.In other words, Parseval formula associated to multi-frequency Gabor transform does not hold any longer [46,Equation 2].Therefore, we simply project the processed output responses onto the image plane via a sum over orientation and phase components in the single frequency results shown in Figures 13 and 14.
In Figures 15 and 16, we perform the same type of experiments as in Figures 11 and 12, respectively, but now with f ∈ {1.00, 1.27, 1.55, 1.82, 2.09, 2.36, 2.64, 2.91, 3.18, 3.45, 3.73, 4.00} and by employing a different type of texture image.Here T = 10 for both rows in Figure 15, and T = 15 for both rows in Figure 16.In Figure 17, we performed the single frequency experiments by using the occluded image found in the bottom row of Figure 15.Similarly, in Figure 18, we performed the same single frequency experiments but now using the occluded image given in the bottom row of Figure 16.Similarly to the case of the previous real test image, single frequency completion cannot perform a proper completion due to the loss of the information corresponding to the other frequencies present in the input image.
Finally, in Figure 19, we compare the result of our algorithm to the results obtained by applying the algorithms explained in [49] and [68].The method proposed in [49] is defined in the classical model framework SE(2) and it combines the SE(2) sub-Riemannian diffusion with a concentration    mechanism resulting in a diffusion driven motion by curvature in SE(2).The algorithm explained in [68] uses a semidiscrete version of the classical model geometry SE (2), and this allows to perform completion via the integration of parallelisable finite set of Mathieu-type diffusions combined with a dynamical restoring mechanism.The main difference between our method and those previously proposed algorithms is that our model uses a higher order sub-Riemannian geometry, this allows to take into account frequency and phase information as well.Moreover, we do not combine in our framework the diffusion procedure with a concentration or dynamical restoring mechanism.We see in Figure 19 that our algorithm produces comparable completion results to the other two methods.We observe that our algorithm is able to preserve the contextual information better than the other two, especially the high frequency structures, thanks to the use of multiple frequencies and the exact inverse Gabor transform.The trade off is that in the corresponding stationary state, especially in the low frequency parts such as the below eye region, the completion is weaker compared to the other two methods.This supports that our algorithm is better adapted to the texture images, such as the ones in Figures 11, 12, than to the natural images such as the one given in Figure 19.In our simulation in Figure 19      The second mechanism, which is the sub-Riemannian diffusion, models the activity propagation between the simple cells in V1.It is concentrated in a neighbourhood along the horizontal integral curves corresponding to the sub-Riemannian model geometry.Those horizontal integral curves are conjectured to be good approximations of the neural connections in V1 [54].Once they are projected onto the two dimensional image plane, their projections overlap closely with the association fields as was shown in Figure 5. Resulting from the sub-Riemannian diffusion, subjective contours of amodal completion are reconstructed in the five dimensional model goemetry.Finally, the inverse Gabor transform provides the representation of the evolved output responses, therefore the reconstructed subjective contours together with the rest of the lifted image, on the two dimensional image plane.This final result is the completed image in which the initially occluded parts are revealed.
One of the novelty of the algorithm is that it is not only an image completion algorithm but it takes into account neurophysiological and psychophysical orientation, frequency and phase constraints observed in the visual cortex.Therefore, it should not be considered as an highly specialized image processing algorithm such as those found in medical imaging, radar imaging, robotics and computer vision.It should be considered rather as an algorithm compatible with a natural geometry which was derived from one of the first step mechanisms of the mammalian visual perception, from receptive profile.In other words, it reflects the cortical architecture.Moreover, it uses multi-frequency and phase channels, which was not the case in the previously proposed completion algorithms using the classical orientation sensitive SE(2) sub-Riemannian framework [33,49] and its variant [56].This allows our algorithm to employ the inverse Gabor transform instead of projecting the processed output responses onto the two dimensional image plane to provide the completed final image; this was not possible in the aforementioned previous methods [33,49,56].This provides good preservation of the contextual information (contours, edges etc.) with different frequencies in the input image.
One of the interesting aspects for future work is to consider a concentration mechanism.It is possible to embed in the proposed completion algorithm a similar concentration mechanism to the one presented in [33], but now with a concentration in each frequency and phase channel.Moreover, the proposed completion algorithm uses the same model geometry as the enhancement algorithm which was presented in [54].Another interesting future work is to combine those two algorithms to perform both completion and enhancement at the same time by employing orientation, frequency and phase information existing in a two dimensional input image.Finally, the study of an analytical solution in a similary way as was done in SE(2) [71,72], but this time for the sub-Riemannian diffusion defined in the five dimensional model geometry, could provide new techniques to perform the completion task, as well as many other image processing applications in the model geometry.This would open new interesting questions especially at theoretical level.

Fig. 4 :
Fig. 4: Two experimental settings from Field, Hayes and Hess [7].A stimuli with aligned patches which we capture (left) and a stimuli plus the background with randomly oriented patches (right) are shown.Abrupt changes in the fragment orientations make it difficult to detect the aligned pattern in the bottom row.

Fig. 6 :
Fig. 6: Two Gabor functions with low (left column) and high (right column) spatial frequencies.Top row: even (real) component of the Gabor functions.Bottom row: odd (imaginary) components.

Fig. 10 :
Fig. 10: Completion of arcs with sinusoidal pattern.Left: original image.Middle left: corrupted image with occluding arcs.Middle right: completed image via the approximate method.Right: completion via the exact method.

Fig. 11 :
Fig. 11: Completion of an occluded real texture image by two different arc patterns on the top and bottom rows.Left: Original image taken from.Middle left: Image with occluding arcs.Middle right: Completed image via the approximate framework.Right: Completed image via the exact framework.

Fig. 12 :
Fig. 12: Completion of a real texture image occluded by two different line patterns on the top and bottom rows.Left: original image.Middle left: image with occluding vertical and horizontal lines.Middle right: completed image via the approximate framework.Right: completed image via the exact framework.

Fig. 13 :
Fig. 13: Single frequency completion associated to the bottom row of Figure 11

Fig. 14 :
Fig.14: Single frequency completion associated to the bottom row of Figure12

Fig. 15 :
Fig. 15: Completion of an occluded real texture image by two different arc patterns on the top and bottom rows.Left: original image taken from [69].Middle left: image with occluding arcs.Middle right: completed image via the approximate framework.Right: completed image via the exact framework.

Fig. 16 :
Fig. 16: Completion of an occluded real texture image by two different line patterns on the top and bottom rows.Left: original image taken from [69].Middle left: image with occluding arcs.Middle right: completed image via the approximate framework.Right: completed image via the exact framework.

Fig. 17 :
Fig.17: Single frequency completion associated to the bottom row of Figure15

Fig. 18 :
Fig.18: Single frequency completion associated to the bottom row of Figure16

Fig. 19 :
Fig.19: Top left: test image from[70].Top right: image to be completed.Bottom left: completion via our method.Bottom middle: completion via the method in[49].Bottom right: completion via the method in[56].

7 .
Conclusion.In this work, we presented a completion algorithm which uses multiple frequency and phase channels to take advantage of the spatial frequency and phase information of a given two dimensional grayscale image.The algorithm consists of three mechanisms: feature extraction, sub-Riemannian diffusion and inverse transform.The first one is a linear filtering of the given input image with the Gabor filter banks.The filtering encodes the visual feature values in the output responses, i.e., orientation, frequency and phase values, of each pixel in the two dimensional input image.The oputput responses are represented in the five dimensional sub-Riemannian model geometry.