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

Emre Baspinar

► To cite this version:

Emre Baspinar. Multi-Frequency Image Completion via a Biologically-Inspired Sub-Riemannian Model with Frequency and Phase. Journal of Imaging, 2021, 7 (12), pp.271. 10.3390/jimag-ing7120271. hal-03509322

HAL Id: hal-03509322

https://hal.science/hal-03509322

Submitted on 22 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Article

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

Emre Baspinar

CNRS/NeuroPSI, Campus CEA, 91400 Saclay, France; emre.baspinar@cnrs.fr

**Abstract:** 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.

**Keywords:** image completion; visual cortex; sub-Riemannian geometry; neurogeometry; differential geometry; Gabor function

<figure><img src="image_3.png" /></figure>

check for
updates

**Citation:** Baspinar, E. Multi-Frequency Image Completion via a Biologically-Inspired Sub-Riemannian Model with Frequency and Phase. *J. Imaging* **2021**, *7*, 271. https://doi.org/10.3390/jimaging7120271

Academic Editor: Edoardo Provenzi

Received: 27 October 2021

Accepted: 1 December 2021

Published: 9 December 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_4.png" /></figure>

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

Visual perception has drawn the attention of experts from fields of philosophy, psychology and neuroscience, as well as the attention of mathematicians and physicists working on perceptual modeling. The question of how we perceive was studied by Edmund Husserl in their pioneering philosophical texts in phenomenology [1–3]. With regard to 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–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 grouping similar fragments rather than purely summing the fragments as if they were indifferent. In terms of visual perception, those similar fragments can be thought of as 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 aligned pieces rather 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 studied how the visual system captures aligned fragments constituting lines or cocircular curves on a background of randomly oriented fragments; see Figure 2. 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.
<figure><img src="image_5.png" /><figcaption>**Figure 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.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 2.** 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.</figcaption></figure>

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 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 fragments of the object evoke our
sensory receptors. Examples of those two categories of completion resulting in subjective contours are the Kanizsa triangle and Ehrenstein's illusion, which are illustrated in Figure 4. We focus 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.

<figure><img src="image_7.png" /><figcaption>**Figure 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.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 4.** **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).</figcaption></figure>

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–14], spatial frequency [15–22], phase [23–26], scale [27] and ocular dominance [17,28,29]. The simple cells are organized in a hypercolumnar architecture, which was first discovered by Hubel and Wiesel [13]. In this architecture, a hypercolumn is assigned to each point (x, y) of the retinal plane $M \approx \mathbb{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 extract the feature values of the visual stimulus, and the activity propagation along the neural connections between the simple cells integrates those values to a coherent global unity. Those two mechanisms, the feature detection and the neural connectivity, comprise the functional geometry of V1.

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 modeling 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 the 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–41] and orientation preference maps [42], as well as frameworks for image processing [43–50], for pattern recognition [51] and for medical applications [52,53], 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–12]. This description provides a good phenomenological approximation of the biologically implemented V1 neural connections, which were reported by Bosking et al. [54]. 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*.

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Left: real association fields. Right: projections of SE(2) horizontal integral curves. Figures are adapted from [7,33].</figcaption></figure>

The sub-Riemannian model geometry proposed in [33] was extended in a recent work to multi-frequency and multi-phase setting [55]. The extended model corresponds to a natural geometry which is *derived* [55,56] 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 a visual feature. The extended model geometry was applied to enhancement on images in which several spatial frequencies were equally present [55], such as texture images. In the present work, we employ the same extended multi-frequency sub-Riemannian setting presented in [55] and propose an image completion algorithm within, which is aimed for use with grayscale texture images containing multiple spatial frequencies.

In Section 2, we explain the model framework which was introduced in [55] 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 [48,50]. At the end, we give the main conclusions and some perspectives for related future research.

## 2. Model Framework

In this section, we explain the cortical model geometry in which our completion algorithm is defined. This model framework is composed of two mechanisms: feature value extraction and neural connectivity [55].

### 2.1. Feature Value Extraction

Each simple cell is sensitive to a specific part of the retina, which is called the *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, the 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 \in M$ on the image plane $M$ and selective to the set of feature values $z \in S^1 \times \mathbb{R}^+ \times S^1$ is denoted by $\Psi_{(q,z)}$, where $q = (q_1, q_2)$ and $z = (\theta, f, \phi)$ together denote a fixed point $(q, z)$ in $Q \simeq \mathbb{R}^2 \times S^1 \times \mathbb{R}^+ \times S^1$. 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,55,57]. In the orientation, frequency and phase selective model framework, the receptive profile of a simple cell is a Gabor function of the following type:

$$ \Psi_{(q,z)}(x, y, s) := \frac{1}{2\sigma^2} e^{-i(r \cdot (x-q_1, y-q_2) - (s-\phi))} e^{-\frac{|x-q_1|^2 + |y-q_2|^2}{2\sigma^2}}, \quad (1) $$

where $f > 0$ represents the spatial frequency (spatial frequency is found via $f = \frac{1}{\lambda}$, where $\lambda > 0$ denotes the wavelength), $r = (r_1, r_2) = (-f \sin \theta, f \cos \theta)$ and $\sigma > 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 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 $\theta$ is the orientation angle to which the simple cell receptive profile is sensitive. Parameter $\phi$ is the reference phase and it creates a phase shift in the waves of the Gabor function as it changes. Here, $s - \phi$ represents the phase centered at $\phi$ and parameter $\phi$ is the reference phase; it creates a phase shift in the waves of the Gabor function as it changes. In other words, the even and odd symmetricity of the receptive profile is modulated with phase shift value to take into account all possible receptive profile patterns which can occur at a given time instant. Note that $q = (q_1, q_2)$ and $(x, y)$ denote the spatial parameters and spatial variables, respectively. Therefore, the vector fields in the cortical space $Q$ are written in terms of derivatives with respect to the variables $x$ and $y$, as is explained further.

We always consider a two-dimensional grayscale static image as a visual input. Each pixel on the image is thought of as a point stimulus. Normally, once the subject is exposed to the image, each pixel is projected onto the retinal surface, from which it is mapped to the cortical surface via retino-cortical maps. For the sake of simplicity, 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) \in Q$ to a given two-dimensional grayscale image $I: M \to [0, 1]$ via the convolution with the Gabor filter:

$$O_I(q, z) = \int_M \Psi_{(q,z)}(x, y, 0) I(x, y) \, dx \, dy. \qquad (2)$$

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 sometimes call this set of output responses *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 representations of the feature values in the five-dimensional V1 model geometry $Q$.

<figure><img src="image_10.png" /><figcaption>**Figure 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.</figcaption></figure>

In general, static receptive profile models based on linear filter banks and static nonlinearities [33,55,58–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 normalization, 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:

$$I(q_1, q_2) = \frac{\sqrt{f}}{\|\Psi\|_{\mathbb{L}^2}} \int_Q O_I(x, y, z) \bar{\Psi}_{(x,y,z)}(q_1, q_2, 0) \, dx \, dy \, dz, \qquad (3)$$

with $\Psi$ denoting the complex conjugate of the corresponding Gabor function $\Psi$.

## 2.2. 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, 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
activity propagation, and therefore interactivity, between simple cells. The activity propagation provides an integrated form of the local feature vectors associated with 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 (1):

$$ \Theta_{(\theta,f)} = -f \sin(\theta) \, dx + f \cos(\theta) \, dy - ds. \qquad (4) $$

The one-form naturally induces the horizontal vector fields corresponding to the model geometry Q. The horizontal vector fields are formally defined as the elements of

$$ \ker(\Theta) := \{X \in TQ : \Theta(X) = 0\}, \qquad (5) $$

where TQ denotes the tangent bundle of Q. The horizontal vector fields corresponding to Q are found from (5) as

$$ \begin{aligned} X_1 &= \cos(\theta) \partial_x + \sin(\theta) \partial_y, & X_2 &= \partial_\theta, \\ X_3 &= -\sin(\theta) \partial_x + \cos(\theta) \partial_y + f \partial_s, & X_4 &= \partial_f. \end{aligned} \qquad (6) $$

Those horizontal vector fields endow Q with a sub-Riemannian structure, as explained in [55]. The sub-Riemannian structure is composed of the manifold Q, the horizontal tangent bundle HQ, which is a subbundle of the tangent bundle TQ, and a scalar product g describing a notion of distance on Q. It is expressed as (Q, HQ, g). The horizontal vector fields span the horizontal tangent space $H_{(q,z)}Q$ at each $(q, z) \in Q$. The horizontal tangent space can be thought of as the analog 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 (6). Despite the degenerate character of the horizontal tangent space $H_{(q,z)}Q$, they provide full connectivity in Q due to the fact that $X_1$ and $X_2$ do not commute, as we see below.

Nonzero commutators of the horizontal vector fields are found as

$$ \begin{aligned} [X_1, X_2] &= \sin(\theta) \partial_x - \cos(\theta) \partial_y, \\ [X_2, X_3] &= -\cos(\theta) \partial_x - \sin(\theta) \partial_y, \\ [X_3, X_4] &= -\partial_s. \end{aligned} \qquad (7) $$

The horizontal vector fields are bracket generating since

$$ T_{(q,z)}Q = \operatorname{span}(X_1, X_2, X_3, X_4, [X_1, X_2])|_{(q,z)}, \qquad (8) $$

for all $(q, z) \in Q$ where $T_{(q,z)}Q$ denotes the tangent space at $(q, z) \in Q$. Indeed, (8) shows that the horizontal vector fields fulfill the 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-65]. This connectivity property has particular importance since it ensures 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 are close approximations of the association fields at the psychophysical level.

We emphasize that the horizontal tangent space is a natural consequence of the choice of the receptive profile function, which is the Gabor function given in (1). The one-form given in (4) induces the horizontal vector fields, which characterize the horizontal tangent space at every point in the cortical space. The horizontal tangent space is endowed with
the aforementioned scalar product g, which is defined only in the horizontal tangent space. This scalar product provides a measure of the length ℓ of a horizontal curve γ connecting not only two spatial points, but also the corresponding orientation, frequency and phase values to those two spatial points. Therefore, the notion of length defined on the horizontal tangent space H(q,z) takes into account how close feature values two simple cells are sensitive to. The length of a horizontal curve is described via ℓ(γ) = ∫√g(γ'(t), γ'(t)) dt. The horizontal curve γ connecting two points in the cortical space Q, i.e., two simple cells, and minimizing the length ℓ is called geodesic. Geodesics can be used as models of the neural connectivity in which a principle of minimizing the energy of activated neurons can also be taken into account. However, the computation of geodesics and the geodesic flow in a high-dimensional sub-Riemannian geometry such as Q is cumbersome and too mathematically involved. Therefore, in the present paper, we restrict ourselves to horizontal curves in general, and consider establishing the neural connectivity based on geodesics in Q as the refinement of the presented framework here.

### 3. 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 [55], 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 [55] 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 $\mathcal{I} = [0, T]$ with $0 < T < \infty$ and consider a horizontal integral curve $(q_1, q_2, \theta, f, \phi) = \gamma : \mathcal{I} \to \mathcal{M}$ associated with the horizontal vector fields given in (6). We denote the initial point of $\gamma$ by $\hat{\alpha} = (\hat{q}_1, \hat{q}_2, \hat{\theta}, \hat{f}, \hat{\phi})$ and its velocity by $\gamma'$. At each time instant $t \in \mathcal{I}$, the velocity is written as a vector $\gamma'(t) \in \text{span}(X_1, X_2, X_3, X_4)(\gamma(t))$ at $\gamma(t) = (q_1(t), q_2(t), \theta(t), f(t), \phi(t)) \in Q$. One way to compute the horizontal integral curves starting from the initial point $\hat{\alpha}$ is to solve the following ODE system for all $t \in \mathcal{I}$:

$$ \gamma'(t) = (c_1 X_1 + c_2 X_2 + c_3 X_3 + c_4 X_4)|_{\gamma(t)}, \quad (9) $$

where $c_{2,3,4}$ denote the coefficients. In the activity propagation machinery which we propose here, the propagation is concentrated along a neighborhood of the horizontal integral curves with constant coefficients $c_{2,3,4}$ in the cortical space Q [55]; see Figures 7 and 8. However, the coefficients need not necessarily be constants in the generic framework of horizontal integral curves.

<figure><img src="image_11.png" /><figcaption>**Figure 7.** A horizontal integral curve along the vector field $X_1 + X_2$. It represents an orientation fiber once $f$ and $\phi$ 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.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 8. Horizontal integral curve fans corresponding to $X_1 + c_2X_2$ (left) and $X_3 + c_4X_4$ (right) where $c_2$ and $c_4$ are varied.</figcaption></figure>

## 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–45,50].

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 $\Sigma \subset Q$ the subspace in which all output responses $O_I$ are defined. These are the output responses obtained as the lifting of the two-dimensional input image $I$. We denote by $\Pi \subset \Sigma$ the subspace in which we only find the output responses corresponding to the part to be completed in the image $I$. The sub-Riemannian diffusion operator for all $(q, z) \in \Sigma$ is defined as

$$ \mathcal{L} := X_1^2 + \beta_2^2 X_2^2 + \beta_3^2 X_3^2 + \beta_4^2 X_4^2, $$

where $X_i^2$ denotes the ith second-order horizontal derivative, and $\beta_{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) \in \Sigma$ by

$$ \partial_t u(q, z, t) = \mathcal{L}u(q, z, t), \quad u(q, z, 0) = O_I(q, z), \quad t \in (0, T], \quad 0 < T < \infty, \quad (10) $$

with the corrupted region marked by the boundary conditions given by $u(\tilde{q}, \tilde{z}, t) = O_I(\tilde{q}, \tilde{z})$ for all $(\tilde{q}, \tilde{z}) \in \Sigma - \Pi$ and $t \in (0, T]$. Here, $T$ denotes a sufficiently large final time and $u: Q \times [0, T] \to \mathbb{C}$ stands for the output responses evolving in time. We denote the number of orientations, frequencies and phases of the $N \times N$ image $I$ by $K$, $L$ and $M$, respectively. We consider the cortical space to be periodic in orientation, frequency and phase dimensions to guarantee the diffusion to remain in axes boundaries corresponding to those dimensions. Then, the coefficients are found as:

$$ \beta_2 = \frac{K}{N\sqrt{2}}, \quad \beta_3 = \frac{L}{N\sqrt{2}}, \quad \beta_4 = \frac{M}{N\sqrt{2}}. $$

We note that $\mathrm{span}(X_1, X_2)$ and $\mathrm{span}(X_3, X_4)$ define two subspaces of the horizontal tangent space $H_{(q,z)}Q$ at each point $(q, z) \in 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 (10) 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,50] for each frequency and phase channel separately by using the SE(2) sub-Riemannian diffusion operator

$$ \tilde{\mathcal{L}} = X_1^2 + \beta_2^2 X_2^2, \qquad (11) $$

where $\beta_2$ provides the unit coherency between the diffusion components in spatial and orientation dimensions. The approximate sub-Riemannian procedure is described by

$$ \partial_t u(q,z,t) = \tilde{\mathcal{L}}u(q,z,t), \quad u(q,z,0) = O_I(q,z), \qquad (12) $$

with the boundary conditions given as $u(\tilde{q}, \tilde{z}, t) = O_I(\tilde{q}, \tilde{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.

## 5. 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). 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 (10) or by (12) if the approximate setting is used. In time, we integrate (10) for the exact framework and (12) for the approximate framework by applying it on the output responses via iterating it with the time step $\Delta t$ until the final time $T$, at which the steady state is reached, i.e., $\partial_t u = 0$. We employ an explicit method where we use B-splined interpolated finite differences to implement the horizontal vector fields given in (6). The final step is to transform back the evolved output responses to the two-dimensional image plane. This is achieved by the inverse Gabor transform given by (3).

### 5.1. Discretization of the Output Responses

We employ a uniform spatial grid to discretize the image plane such that

$$ I[i, j] = I(i\Delta x, j\Delta y), \qquad (13) $$

where $i, j \in \{1, 2, \dots, N\}$, with $N$ denoting the image size and $\Delta x, \Delta y \in \mathbb{R}^+$ denoting the pixel width. In our case, the input images are square images and $\Delta x = \Delta 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 $\Delta\theta$, in the frequency dimension with $\Delta f$ and in the phase dimension with $\Delta s$. Discretized output response $O_I(q_{1,i}, q_{2,j}, \theta_k, f_l, \phi_m)$ given to $I[i, j]$ on uniform orientation, frequency and phase grids with points $\theta_k = k\Delta\theta, f_l = l\Delta f$ and $\phi_m = m\Delta s$ ($k \in \{1, 2, \dots, K\}, l \in \{1, 2, \dots, L\}, m \in \{1, 2, \dots, M\}$) is denoted by

$$ O_I[i, j, k, l, m] = O_I(q_{1,i}, q_{2,j}, \theta_k, f_l, \phi_m), \qquad (14) $$

where $q_{1,i} = i\Delta x$ and $q_{2,j} = j\Delta y$.

The Gabor function given by (1) is written in the discrete setting as follows:

$$ \Psi_{[i,j,k,l,m]}[\tilde{i}, \tilde{j}, \tilde{n}] = \Psi_{(q_{1,i}, q_{2,j}, \theta_k, f_l, \phi_m)}(x_{\tilde{i}}, y_{\tilde{j}}, s_{\tilde{n}}), \qquad (15) $$

where $\tilde{i}, \tilde{j} \in \{1, 2, \dots, N\}, \tilde{k} \in \{1, 2, \dots, K\}, \tilde{n} \in \{1, 2, \dots, M\}$, for each orientation $\theta_k$, frequency $f_l$ and phase $\phi_m$. We fix $s_{\tilde{n}} = 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:

$$O_I[i, j, k, l, m] = \sum_{\tilde{i}, \tilde{j}} \Psi_{[i, j, k, l, m]}[\tilde{i}, \tilde{j}, 0] I[\tilde{i}, \tilde{j}]. \quad (16)$$

We discretize the time interval by $V \in \mathbb{N}^+$ samples and denote it by $h_v$. Here, $h_v$ is the time instant $h_v = v \Delta t$, with $\Delta t$ satisfying $T = V \Delta t$ and $v \in \{1, 2, \dots, V\}$. Discretized evolving output response is written as

$$U[i, j, k, l, m, v] = u(q_{1,i}, q_{2,j}, \theta_k, f_l, \phi_m, h_v) \quad (17)$$

Finally, the discrete inverse transform applied via a normalized kernel $\Psi$ on the evolved output responses until the final time $T$ gives the completed two-dimensional image $I_T$, which is found as follows:

$$I_T[i, j] = \sum_{\tilde{i}, \tilde{j}, \tilde{k}, \tilde{m}} \sqrt{f_{\tilde{i}}} \sum_{\tilde{l}} U_T[\tilde{i}, \tilde{j}, \tilde{k}, \tilde{l}, \tilde{m}] \Psi_{[\tilde{i}, \tilde{j}, \tilde{k}, \tilde{l}, \tilde{m}]}[i, j, 0]. \quad (18)$$

## 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 (10) and (12), 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 (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)$

$$\mathrm{sp}(x, y) = \sum_{i,j \in \mathbb{Z}} b(i, j) \rho(x - i, y - j). \quad (19)$$

The coefficients are determined such that the spline polynomial $\mathrm{sp}(x, y)$, together with the B-spline basis functions $\rho(x - i, y - j)$, coincides with the horizontal derivatives at the grid points. For example, the condition $\mathrm{sp}(i\Delta x, j\Delta 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

$$\begin{aligned} e_{\xi}^{k} &:= (\Delta x \cos(\theta_k), \Delta y \sin(\theta_k)), \\ e_{\eta}^{k} &:= (-\Delta x \sin(\theta_k), \Delta y \cos(\theta_k)), \end{aligned} \quad (20)$$

whose illustrations corresponding to $\Delta x = \Delta y = 1$ case are given in Figure 9. We abuse the notation to denote the evolving output responses:

$$U = U[i, j, k, l, m, v] = u(q_{1,i}, q_{2,j}, \theta_k, f_l, \phi_m, h_v), \quad (21)$$

and then, we write the central finite differences of the second-order horizontal derivatives as
$$X_1 X_1 U \approx \frac{1}{(\Delta x)^2} \left( u(q + e_{\xi}^k, \theta_k, f_l, \phi_m) - 2u(q, \theta_k, f_l, \phi_m) \right. \\ \left. \qquad + u(q - e_{\xi}^k, \theta_k, f_l, \phi_m) \right),$$

$$X_2 X_2 U \approx \frac{1}{(\Delta \theta)^2} \left( u(q, \theta_{k+1}, f_l, \phi_m) - 2u(q, \theta_k, f_l, \phi_m) \right. \\ \left. \qquad + u(q, \theta_{k-1}, f_l, \phi_m) \right),$$

$$X_3 X_3 U \approx \frac{1}{(\Delta x)^2} \left( u(q + e_{\eta}^k, \theta_k, f_l, \phi_m) - 2u(q, \theta_k, f_l, \phi_m) + u(q - e_{\eta}^k, \theta_k, f_l, \phi_m) \right) \\ + \frac{f_l^2}{(\Delta s)^2} \left( u(q, \theta_k, f_l, \phi_{m+1}) - 2u(q, \theta_k, f_l, \phi_m) \right. \\ \left. + u(q, \theta_k, f_l, \phi_{m-1}) \right) + \frac{f \cos(\theta_k)}{2 \Delta s \Delta_x} \left( u(q + e_{\eta}^k, \theta_k, f_l, \phi_{m+1}) \right. \\ \left. - u(q - e_{\eta}^k, \theta_k, f_l, \phi_{m+1}) - u(q + e_{\eta}^k, \theta_k, f_l, \phi_{m-1}) \right. \\ \left. + u(q - e_{\eta}^k, \theta_k, f_l, \phi_{m-1}) \right) - \frac{f \sin(\theta_k)}{2 \Delta s \Delta_x} \left( u(q + e_{\xi}^k, \theta_k, f_l, \phi_{m+1}) \right. \\ \left. - u(q - e_{\xi}^k, \theta_k, f_l, \phi_{m+1}) - u(q + e_{\xi}^k, \theta_k, f_l, \phi_{m-1}) \right. \\ \left. + u(q - e_{\xi}^k, \theta_k, f_l, \phi_{m-1}) \right) + \frac{f^2}{2 \Delta s} \left( u(q, \theta_k, f_l, \phi_{m+1}) \right. \\ \left. - 2u(q, \theta_k, f_l, \phi) + u(q, \theta_k, f_l, \phi_{m-1}) \right),$$

$$X_4 X_4 u[i, j, k, l, m] \approx \frac{1}{(\Delta f)^2} \left( u(q, \theta_k, f_{l+1}, \phi_m) - 2u(q, \theta_k, f_l, \phi_m) + u(q, \theta_k, f_{l-1}, \phi_m) \right).$$

Finally, we write the discretized numerical iteration for (10) and (12) as follows:

$$\begin{aligned} U[i, j, k, l, m, v] &= u(q_{i,1}, q_{j,2}, \theta_k, f_l, \phi_m, h_v) \\ &= u(q_{i,1}, q_{j,2}, \theta_k, f_l, \phi_m, h_{v-1}) + \Delta t \bar{\mathcal{L}} u(q_{i,1}, q_{j,2}, \theta_k, f_l, \phi_m, h_{v-1}), \end{aligned} \quad (23)$$

where $\bar{\mathcal{L}}$ represents the discretized version of either $\mathcal{L}$ or $\bar{\mathcal{L}}$, depending on which one between the exact and approximate frameworks is chosen. The discretization is achieved by replacing the second-order horizontal derivatives with their discrete versions given in (22).

<figure><img src="image_13.png" /><figcaption>**Figure 9.** Illustration of the vectors $e_{\xi}^k$ and $e_{\eta}^k$ at (0,0) with $\Delta x = \Delta y = 1$. The figure was modified and adapted from [44].</figcaption></figure>
### 5.3. Pseudocode of the Algorithm

We denote the processed two-dimensional image at the final time $T$ by $I_T$ and the evolving discrete output responses at the time instant $v \Delta t$ by $U_v$. Then, we provide a general scheme of the exact completion Algorithm 1 as follows:

<caption>**Algorithm 1:** Completion algorithm pseudocode.</caption>

**Data:** $N \times N$ input grayscale image $I$
**Parameters:** $\sigma, \beta_2, \beta_3, \beta_4, T, K, L, M, \Delta t, \Delta\theta, \Delta f, \Delta s, \text{tol}$
**Result:** Processed image $I_T$
Compute the lifting $O_I$ via (16);
Initialize the iteration index $v \leftarrow 0$;
repeat
  $v \leftarrow v + 1$;
  Compute the right hand side of (23);
  Update $U_v$ via (23);
until ($\|U_v - U_{v-1}\|_{\mathbb{L}^2} / \|U_v\|_{\mathbb{L}^2} < \text{tol}$;
Inverse transform the evolved output responses $U_T$ to obtain the completed image $I_T$ via (18).

The detailed Matlab® code can be found in the following link: https://drive.google.com/file/d/1YET47AxYo5_FQs3ObREAqvQp2KjuKh7n/view (accessible starting from 27 October 2021).

## 6. Numerical Experiments

We choose $\sigma = 2$ as the scale of the localizing Gaussian given in (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 $\theta \in \{0, \frac{\pi}{32}, \frac{2\pi}{32}, \frac{3\pi}{32}, \dots, \frac{31\pi}{32}\}$, $\phi \in \{0, \frac{\pi}{8}, \frac{2\pi}{8}, \frac{3\pi}{8}, \frac{4\pi}{8}\}$ and $\Delta 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, with a sinusoidal function whose frequency increases linearly as the radius increases. The arcs with sinusoidal patterns 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 \in \{2, 2.5455, 3.0909, \dots, 8\}$. We observe that both the approximate and exact frameworks show a similar performance and they provide proper completion.

<figure><img src="image_14.png" /><figcaption>**Figure 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.</figcaption></figure>

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 \in \{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 for $T = 15$ in Figure 12 on the same texture image but now occluded by vertical and
horizontal bars. In such a 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 in all cases, both the approximate and the exact algorithms are able to complete 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 perform 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 is used. The parameters are the same as in the bottom row results of Figure 11. We use *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 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, the Parseval formula associated with 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 \in \{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 using 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 display the single-frequency experiments by using the occluded image found in the bottom row of Figure 15. Similarly, in Figure 18, we show 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.

<figure><img src="image_15.png" /><figcaption>**Figure 11.** Completion of an occluded real texture image by two different arc patterns on the top and bottom rows. **Left:** Original image. **Middle left:** Image with occluding arcs. **Middle right:** Completed image via the approximate framework. **Right:** Completed image via the exact framework.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 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.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 13.** Single-frequency completion associated with the bottom row of Figure 11.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 14.** Single-frequency completion associated with the bottom row of Figure 12.</figcaption></figure>
Finally, in Figure 19, we compare the result of our algorithm to the results obtained by applying the algorithms explained in [48,50]. The method proposed in [50] 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 [48] uses a semidiscrete version of the classical model geometry SE(2), and this allows us to perform completion via the integration of a parallelizable 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 us 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 completion results comparable 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 the idea that our algorithm is better adapted to the texture images, such as the ones in Figures 11 and 12, than to the natural images, such as the one given in Figure 19. In our simulation in Figure 19, we use $f \in \{1.50, 2.09, 2.68, 3.27, 3.86, 4.45, 5.05, 5.64, 6.23, 6.82, 7.41, 8.00\}$ and $T = 50$.

<figure><img src="image_19.png" /><figcaption>**Figure 15.** Completion of an occluded real texture image by two different arc patterns on the top and bottom rows. **Left:** original image taken from [68]. **Middle left:** image with occluding arcs. **Middle right:** completed image via the approximate framework. **Right:** completed image via the exact framework.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure 16.** Completion of an occluded real texture image by two different line patterns on the top and bottom rows. **Left:** original image taken from [68]. **Middle left:** image with occluding arcs. **Middle right:** completed image via the approximate framework. **Right:** completed image via the exact framework.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 17.** Single-frequency completion associated with the bottom row of Figure 15.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 18.** Single-frequency completion associated with the bottom row of Figure 16.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**Figure 19.** Top left: test image from [69]. Top right: image to be completed. Bottom left: completion via our method. Bottom middle: completion via the method in [50]. Bottom right: completion via the method in [48].</figcaption></figure>

## 7. Conclusions

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 output responses are represented in the five-dimensional sub-Riemannian model geometry. The second mechanism, which is the sub-Riemannian diffusion, models the activity propagation between the simple cells in V1. It is concentrated in a neighborhood 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 [55]. 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 geometry. 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, are shown on the two-dimensional image plane. This final result is the completed image in which the initially occluded parts are revealed.

One of the novelties 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,50] and its variant [48]. 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,48,50]. 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. In the proposed completion algorithm, it is possible to embed 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 [55]. 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 similar way as was carried out in SE(2) [70,71], but this time, a study regarding 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 up new interesting questions, especially at theoretical level.

**Funding:** This research was funded by the Human Brain Project (European Union grant H2020-945539).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The code can be found in the following link: https://drive.google.com/file/d/1YET47AxYo5_FQs3ObREAqvQp2KjuKh7n/view (accessible starting from 27 October 2021).

**Acknowledgments:** I would like to thank Giovanna Citti and Alessandro Sarti for their inspiring discussions. I hereby acknowledge that I was supported by the Human Brain Project (European Union grant H2020-945539).

**Conflicts of Interest:** The author declares no conflict of interest.

## References

1. Husserl, E.; Holenstein, E. *Logische Untersuchungen*. Springer Netherlands: The Hague, The Netherlands, 1975.

2. Husserl, E.; Claesges, U. *U. Ding und Raum: Vorlesungen 1907*; Springer Netherlands: The Hague, The Netherlands, 1973.

3. Husserl, E.; Schumann K. *Ideen zu einer reinen Phänomenologie und phänomenologischen Philosophie*; Springer Netherlands: The Hague, The Netherlands, 1976.

4. Wertheimer, M. Laws of Organization in Perceptual Forms. In *A Source Book of Gestalt Psychology*; Ellis, W.D., Ed.; Trubner & Company: Kegan Paul, Trench, 1938; pp. 71–88.

5. Köhler, W. *Gestalt Psychology: An Introduction to New Concepts in Modern Psychology*; Liveright: New York, NY, USA, 1970.

6. Koffka, K. *Principles of Gestalt Psychology*; Routledge: Abingdon, UK, 2013; Volume 44.

7. Field, D.J.; Hayes, A.; Hess, R.F. Contour integration by the human visual system: Evidence for a local "association field". *Vis. Res.* **1993**, *33*, 173–193. [CrossRef]

8. Kanizsa, G. *Organization in Vision: Essays on Gestalt Perception*; Praeger Publishers: Westport, CT, USA, 1979.

9. Kanizsa, G. *Grammatica del Vedere: Saggi su Percezione e Gestalt*; Il Mulino: Bologna, Italy, 1980.

10. Hubel, D.H.; Wiesel, T.N. Receptive fields of single neurones in the cat's striate cortex. *J. Physiol.* **1959**, *148*, 574. [CrossRef] [PubMed]

11. Hubel, D.H.; Wiesel, T.N. Receptive fields, binocular interaction and functional architecture in the cat's visual cortex. *J. Physiol.* **1962**, *160*, 106. [CrossRef] [PubMed]

12. Hubel, D.H.; Wiesel, T.N. Shape and arrangement of columns in cat's striate cortex. *J. Physiol.* **1963**, *165*, 559. [CrossRef] [PubMed]

13. Hubel, D.H.; Wiesel, T.N. Uniformity of monkey striate cortex: A parallel relationship between field size, scatter, and magnification factor. *J. Comp. Neurol.* **1974**, *158*, 295–305. [CrossRef]

14. Hubel, D.H.; Wiesel, T.N. Ferrier lecture: Functional architecture of macaque monkey visual cortex. *Proc. R. Soc. Lond. Biol. Sci.* **1977**, *198*, 1–59.

15. Maffei, L.; Fiorentini, A. Spatial frequency rows in the striate visual cortex. *Vis. Res.* **1977**, *17*, 257–264. [CrossRef]

16. Hübener, M.; Shoham, D.; Grinvald, A.; Bonhoeffer, T. Spatial relationships among three columnar systems in cat area 17. *J. Neurosci.* **1997**, *17*, 9270–9284. [CrossRef]
17. Issa, N.P.; Trepel, C.; Stryker, M.P. Spatial frequency maps in cat visual cortex. *J. Neurosci.* **2000**, *20*, 8504–8514. [CrossRef]

18. Issa, N.P.; Rosenberg, A.; Husson, T.R. Models and measurements of functional maps in V1. *J. Neurophysiol.* **2008**, *99*, 2745–2754. [CrossRef]

19. Sirovich, L.; Uglesich, R. The organization of orientation and spatial frequency in primary visual cortex. *Proc. Natl. Acad. Sci. USA* **2004**, *101*, 16941–16946. [CrossRef] [PubMed]

20. Tani, T.; Ribot, J.; O'Hashi, K.; Tanaka, S. Parallel development of orientation maps and spatial frequency selectivity in cat visual cortex. *Eur. J. Neurosci.* **2012**, *35*, 44–55. [CrossRef] [PubMed]

21. Ribot, J.; Aushana, Y.; Bui-Quoc, E.; Milleret, C. Organization and origin of spatial frequency maps in cat visual cortex. *J. Neurosci.* **2013**, *33*, 13326–13343. [CrossRef]

22. Ribot, J.; Romagnoni, A.; Milleret, C.; Bennequin, D.; Touboul, J. Pinwheel-dipole configuration in cat early visual cortex. *NeuroImage* **2016**, *128*, 63–73. [CrossRef] [PubMed]

23. De Valois, K.; Tootell, R.B. Spatial-frequency-specific inhibition in cat striate cortex cells. *J. Physiol.* **1983**, *336*, 359–376. [CrossRef]

24. Pollen, D.A.; Gaska, J.P.; Jacobson, L.D. Responses of simple and complex cells to compound sine-wave gratings. *Vis. Res.* **1988**, *28*, 25–39. [CrossRef]

25. Levitt, J.B.; Sanchez, R.M.; Smith, E.L.; Movshon, J.A. Spatio-temporal interactions and the spatial phase preferences of visual neurons. *Exp. Brain Res.* **1990**, *80*, 441–445. [CrossRef]

26. Mechler, F.; Reich, D.S.; Victor, J.D. Detection and discrimination of relative spatial phase by V1 neurons. *J. Neurosci.* **2002**, *22*, 6129–6157. [CrossRef]

27. Blakemore, C.T.; Campbell, F.W. On the existence of neurones in the human visual system selectively sensitive to the orientation and size of retinal images. *J. Physiol.* **1969**, *203*, 237. [CrossRef]

28. Shatz, C.J.; Stryker, M.P. Ocular dominance in layer iv of the cat's visual cortex and the effects of monocular deprivation. *J. Physiol.* **1978**, *281*, 267–283. [CrossRef]

29. LeVay, S.; Stryker, M.P.; Shatz, C.J. Ocular dominance columns and their development in layer iv of the cat's visual cortex: A quantitative study. *J. Comp. Neurol.* **1978**, *179*, 223–244. [CrossRef]

30. Hoffman, W.C. Higher visual perception as prolongation of the basic lie transformation group. *Math. Biosci.* **1970**, *6*, 437–471. [CrossRef]

31. Hoffman, W.C. The visual cortex is a contact bundle. *Appl. Math. Comput.* **1989**, *32*, 137–167. [CrossRef]

32. Petitot, J.; Tondut, Y. Vers une neurogéométrie. fibrations corticales, structures de contact et contours subjectifs modaux. *Math. Sci. Hum.* **1999**, *145*, 5–101. [CrossRef]

33. Citti, G.; Sarti, A. A cortical based model of perceptual completion in the roto-translation space. *J. Math. Imaging Vis.* **2006**, *24*, 307–326. [CrossRef]

34. Sarti, A.; Citti, G.; Petitot, J. The symplectic structure of the primary visual cortex. *Biol. Cybern.* **2008**, *98*, 33–48. [CrossRef]

35. Barbieri, D.; Citti, G.; Cocci, G.; Sarti, A. A cortical-inspired geometry for contour perception and motion integration. *J. Math. Imaging Vis.* **2014**, *49*, 511–529. [CrossRef]

36. Cocci, G.; Barbieri, D.; Citti, G.; Sarti, A. Cortical spatiotemporal dimensionality reduction for visual grouping. *Neural Comput.* **2015**, *27*, 1252–1293. [CrossRef] [PubMed]

37. Franceschiello, B.; Mashtakov, A.; Citti, G.; Sarti, A. Geometrical optical illusion via sub-Riemannian geodesics in the roto-translation group. *Differ. Geom. Its Appl.* **2019**, *65*, 55–77. [CrossRef]

38. Franceschiello, B.; Sarti, A.; Citti, G. A neuromathematical model for geometrical optical illusions. *J. Math. Imaging Vis.* **2018**, *60*, 94–108. [CrossRef]

39. Bertalmío, M.; Calatroni, L.; Franceschi, V.; Franceschiello, B.; Villa, A.G.; Prandi, D. Visual illusions via neural dynamics: Wilson-Cowan-type models and the efficient representation principle. *J. Neurophysiol.* **2020**, *123*, 1606–1618. [CrossRef]

40. Bertalmio, M.; Calatroni, L.; Franceschi, V.; Franceschiello, B.; Prandi, D. Cortical-inspired Wilson-Cowan-type equations for orientation-dependent contrast perception modelling. *J. Math. Imaging Vis.* **2021**, *63*, 263–281. [CrossRef]

41. Baspinar, E.; Calatroni, L.; Franceschi, V.; Prandi, D. A cortical-inspired sub-Riemannian model for Poggendorff-type visual illusions. *J. Imaging* **2021**, *7*, 41. [CrossRef]

42. Baspinar, E.; Citti, G.; Sarti, A. A geometric model of multi-scale orientation preference maps via Gabor functions. *J. Math. Imaging Vis.* **2018**, *60*, 900–912. [CrossRef]

43. Duits, R.; Franken, E. Left-invariant parabolic evolutions on SE(2) and contour enhancement via invertible orientation scores part I: Linear left-invariant diffusion equations on SE(2). *Q. Appl. Math.* **2010**, *68*, 255–292. [CrossRef]

44. Duits, R.; Franken, E. Left-invariant parabolic evolutions on se(2) and contour enhancement via invertible orientation scores part ii: Nonlinear left-invariant diffusions on invertible orientation scores. *Q. Appl. Math.* **2010**, *68*, 293–331. [CrossRef]

45. Boscain, U.; Duplaix, J.; Gauthier, J.P.; Rossi, F. Anthropomorphic image reconstruction via hypoelliptic diffusion. *SIAM J. Control Optim.* **2012**, *50*, 1309–1336. [CrossRef]

46. Duits, R.; Führ, H.; Janssen, B.; Bruurmijn, M.; Florack, L.; van Assen, H. Evolution equations on Gabor transforms and their applications. *Appl. Comput. Harmon. Anal.* **2013**, *35*, 483–526. [CrossRef]

47. Boscain, U.; Gauthier, J.F.; Prandi, D.; Remizov, A. Image reconstruction via non-isotropic diffusion in Dubins/Reed-Shepp-like control systems. In *Proceedings of the 53rd IEEE Conference on Decision and Control, Los Angeles, CA, USA*, 15–17 December 2014; IEEE: Piscataway, NJ, USA, 2014; pp. 4278–4283.
48. Boscain, U.; Chertovskih, R.A.; Gauthier, J.F.; Remizov, A.O. Hypoelliptic diffusion and human vision: A semidiscrete new twist. *SIAM J. Imaging Sci.* **2014**, *7*, 669–695. [CrossRef]

49. Sharma, U.; Duits, R. Left-invariant evolutions of wavelet transforms on the similitude group. *Appl. Comput. Harmon. Anal.* **2015**, *39*, 110–137. [CrossRef]

50. Citti, G.; Franceschiello, B.; Sanguinetti, G.; Sarti, A. Sub-Riemannian mean curvature flow for image processing. *SIAM J. Imaging Sci.* **2016**, *9*, 212–237. [CrossRef]

51. Prandi, D.; Gauthier, J.-P. *A Semidiscrete Version of the Citti-Petitot-Sarti Model as a Plausible Model for Anthropomorphic Image Reconstruction and Pattern Recognition*; SpringerBriefs in Mathematics; Springer: Cham, Switzerland, 2018; [CrossRef]

52. ter Haar Romeny, B.M.; Bekkers, E.J.; Zhang, J.; Abbasi-Sureshjani, S.; Huang, F.; Duits, R.; Dashtbozorg, B.; Berendschot, T.T.J.M.; Smit-Ockeloen, I.; Eppenhof, K.A.J.; et al. Brain-inspired algorithms for retinal image analysis. *Mach. Vis. Appl.* **2016**, *27*, 1117–1135. [CrossRef]

53. Bekkers, E.; Duits, R.; Berendschot, T.; Romeny, B.H. A multi-orientation analysis approach to retinal vessel tracking. *J. Math. Imaging Vis.* **2014**, *49*, 583–610. [CrossRef]

54. Bosking, W.H.; Zhang, Y.; Schofield, B.; Fitzpatrick, D. Orientation selectivity and the arrangement of horizontal connections in tree shrew striate cortex. *J. Neurosci.* **1997**, *17*, 2112–2127. [CrossRef] [PubMed]

55. Baspinar, E.; Sarti, A.; Citti, G. A sub-Riemannian model of the visual cortex with frequency and phase. *J. Math. Neurosci.* **2020**, *10*, 1–31. [CrossRef] [PubMed]

56. Baspinar, E. Minimal Surfaces in sub-Riemannian Structures and Functional Geometry of the Visual Cortex. Dissertation Thesis, Alma Mater Studiorum Università di Bologna, Bologna, Italy, 2018.

57. Daugman, J.G. Uncertainty relation for resolution in space, spatial frequency, and orientation optimized by two-dimensional visual cortical filters. *JOSA A* **1985**, *2*, 1160–1169. [CrossRef]

58. Bekkers, E.J.; Lafarge, M.W.; Veta, M.; Eppenhof, K.A.J.; Pluim, J.P.W.; Duits, R. Roto-translation covariant convolutional networks for medical image analysis. In *International Conference on Medical Image Computing and Computer-Assisted Intervention*; Springer: Berlin/Heidelberg, Germany, 2018; pp. 440–448.

59. Duits, R.; Duits, M.; van Almsick, M.; Romeny, B.H. Invertible orientation scores as an application of generalized wavelet theory. *Pattern Recognit. Image Anal.* **2007**, *17*, 42–75. [CrossRef]

60. Koenderink, J.J.; van Doorn, A.J. Representation of local geometry in the visual system. *Biol. Cybern.* **1987**, *55*, 367–375. [CrossRef] [PubMed]

61. Lindeberg, T. A computational theory of visual receptive fields. *Biol. Cybern.* **2013**, *107*, 589–635. [CrossRef]

62. Hörmander, L. Hypoelliptic second order differential equations. *Acta Math.* **1967**, *119*, 147–171. [CrossRef]

63. Rashevskii, P.K. Any two points of a totally nonholonomic space may be connected by an admissible line. *Uch. Zap. Ped. Inst. Liebknechta, Ser. Phys. Math.* **1938**, *2*, 83–94.

64. Chow, W.-L. Über systeme von linearen partiellen differentialgleichungen erster ordnung. *Math. Ann.* **1940**, *117*, 98–105. [CrossRef]

65. Agrachev, A.; Barilari, D.; Boscain, U. *A Comprehensive Introduction to Sub-Riemannian Geometry*; Cambridge University Press: Cambridge, UK, 2019; Volume 181.

66. Unser, M. Splines: A perfect fit for signal and image processing. *IEEE Signal Process. Mag.* **1999**, *16*, 22–38. [CrossRef]

67. Franken, E.M. *Enhancement of Crossing Elongated Structures in Images*; Eindhoven University of Technology: Eindhoven, The Netherlands.

68. Kimmel, R.; Malladi, R.; Sochen, N. Images as embedded maps and minimal surfaces: Movies, color, texture, and volumetric medical images. *Int. J. Comput. Vis.* **2000**, *39*, 111–129. [CrossRef]

69. Boscain, U.V.; Chertovskih, R.; Gauthier, J.E.; Prandi, D.; Remizov, A. Highly corrupted image inpainting through hypoelliptic diffusion. *J. Math. Imaging Vis.* **2018**, *60*, 1231–1245. [CrossRef]

70. Duits, R.; Almsick, M.V. The explicit solutions of linear left-invariant second order stochastic evolution equations on the 2d euclidean motion group. *Q. Appl. Math.* **2008**, *66*, 27–67. [CrossRef]

71. Zhang, J.; Duits, R.; Sanguinetti, G.; ter Haar Romeny, B.M. Numerical approaches for linear left-invariant diffusions on se (2), their comparison to exact solutions, and their applications in retinal imaging. *Numer. Math. Theory Methods Appl.* **2016**, *9*, 1–50. [CrossRef]