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

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.

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,2,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,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 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.

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.

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,14], spatial frequency [15,16,17,18,19,20,21,22], phase [23,24,25,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\simeq {\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 exctract 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 ($\mathrm{SE}\left(2\right)$) 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,50], for pattern recognition [51] and for medical applications [52,53], were proposed in the sub-Riemannian geometry of $\mathrm{SE}\left(2\right)$.

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. [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 $\mathrm{SE}\left(2\right)$ 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.

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 $\mathrm{SE}\left(2\right)$ 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 ${\mathsf{\Psi }}_{(q,z)}$, where $q=(q_1,q_2)$ and $z=(\theta ,f,\varphi )$ together denote a fixed point $(q,z)$ in Q≃${\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:

$${\mathsf{\Psi }}_{(q,z)}(x,y,s):=\frac{1}{2{\sigma }^2}{e}^{-i\left(r·(x-q_1,y-q_2)-(s-\varphi )\right)}{e}^{-\frac{|x-q_1{|}^2+{\left|y-q_2\right|}^2}{2{\sigma }^2}},$$

(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)=(-fsin\theta ,fcos\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 $\varphi$ is the reference phase and it creates a phase shift in the waves of the Gabor function as it changes.Here, $s-\varphi$ represents the phase centered at $\varphi$ and parameter $\varphi$ 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:

$$\begin{array}{c}\hfill O_I(q,z)={\int }_M{\mathsf{\Psi }}_{(q,z)}(x,y,0)I(x,y)dxdy.\end{array}$$

(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.

In general, static receptive profile models based on linear filter banks and static nonlinearities [33,55,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 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}}{{\parallel \mathsf{\Psi }\parallel }_{{\mathbb{L}}^2}}{\int }_Q{O}_I(x,y,z){\overline{\mathsf{\Psi }}}_{(x,y,z)}(q_1,q_2,0)dxdydz,$$

(3)

with $\overline{\mathsf{\Psi }}$ denoting the complex conjugate of the corresponding Gabor function $\mathsf{\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):

$${\mathsf{\Theta }}_{(\theta ,f)}=-fsin\left(\theta \right)dx+fcos\left(\theta \right)dy-ds.$$

(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\left(\mathsf{\Theta }\right):=\{X\in TQ:\mathsf{\Theta }(X)=0\},$$

(5)

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

$$\begin{array}{cc}\hfill X_1=& cos\left(\theta \right){\partial }_x+sin\left(\theta \right){\partial }_y,X_2={\partial }_{\theta },\hfill \\ \hfill X_3=& -sin\left(\theta \right){\partial }_x+cos\left(\theta \right){\partial }_y+f{\partial }_s,X_4={\partial }_f.\hfill \end{array}$$

(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{array}{cc}\hfill [X_1,X_2]=& sin\left(\theta \right){\partial }_x-cos\left(\theta \right){\partial }_y,\hfill \\ \hfill [X_2,X_3]=& -cos\left(\theta \right){\partial }_x-sin\left(\theta \right){\partial }_y,\hfill \\ \hfill [X_3,X_4]=& -{\partial }_s.\hfill \end{array}$$

(7)

The horizontal vector fields are bracket generating since

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

(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,64,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 aformentioned 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 $\gamma$ 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 $\ell \left(\gamma \right)=\int \sqrt{g({\gamma }^{\prime }\left(t\right),{\gamma }^{\prime }\left(t\right))}dt$. The horizontal curve $\gamma$ 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 $\mathrm{SE}\left(2\right)$ 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 $\mathrm{SE}\left(2\right)$, 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,\varphi )=\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{\varphi })$ and its velocity by ${\gamma }^{\prime }$. At each time instant $t\in \mathcal{I}$, the velocity is written as a vector ${\gamma }^{\prime }\left(t\right)\in span(X_1,X_2,X_3,X_4)\left(\gamma \left(t\right)\right)$ at $\gamma \left(t\right)=(q_1\left(t\right),q_2\left(t\right),\theta \left(t\right),f\left(t\right),\varphi \left(t\right))\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}$:

$$\begin{array}{c}\hfill {\gamma }^{\prime }\left(t\right)=(c_1{X}_1+c_2{X}_2+c_3{X}_3+c_4{X}_4){|}_{\gamma \left(t\right)},\end{array}$$

(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 Figure 7 and Figure 8. However, the coefficients need not necessarily be constants in the generic framework of horizontal integral curves.

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 $\mathrm{SE}\left(2\right)$ 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,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 $\mathsf{\Sigma }\subset 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 denote by $\Pi \subset \mathsf{\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 \mathsf{\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 \mathsf{\Sigma }$ by

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

(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 \mathsf{\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}},{\beta }_3=\frac{L}{N\sqrt{2}},{\beta }_4=\frac{M}{N\sqrt{2}}.$$

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)\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 $\mathrm{SE}\left(2\right)$ [33,45,50] for each frequency and phase channel separately by using the $\mathrm{SE}\left(2\right)$ sub-Riemannian diffusion operator

$$\tilde{\mathcal{L}}=X_1^2+{\beta }_2^2{X}_2^2,$$

(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),u(q,z,0)=O_I(q,z),$$

(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 achived 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),$$

(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,{\varphi }_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 ${\varphi }_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,{\varphi }_m),$$

(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:

$${\mathsf{\Psi }}_{[i,j,k,l,m]}[\tilde{i},\tilde{j},\tilde{n}]={\mathsf{\Psi }}_{(q_{1,i},q_{2,j},{\theta }_k,f_l,{\varphi }_m)}(x_{\tilde{i}},y_{\tilde{j}},s_{\tilde{n}}),$$

(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 ${\varphi }_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}}{\mathsf{\Psi }}_{[i,j,k,l,m]}[\tilde{i},\tilde{j},0]I[\tilde{i},\tilde{j}].$$

(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,{\varphi }_m,h_v)$$

(17)

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

$$\begin{array}{c}\hfill I_T[i,j]=\sum _{\tilde{i},\tilde{j},\tilde{k},\tilde{m}}\sqrt{f_{\tilde{l}}}\sum _{\tilde{l}}{U}_T[\tilde{i},\tilde{j},\tilde{k},\tilde{l},\tilde{m}]{\overline{\mathsf{\Psi }}}_{[\tilde{i},\tilde{j},\tilde{k},\tilde{l},\tilde{m}]}[i,j,0].\end{array}$$

(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)$

$$sp(x,y)=\sum _{i,j\in Z}b(i,j)\rho (x-i,y-j).$$

(19)

The coefficients are determined such that the spline polynomial $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 $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{array}{c}\hfill \begin{array}{cc}\hfill e_{\xi }^k:=& (\Delta xcos\left({\theta }_k\right),\Delta ysin\left({\theta }_k\right)),\hfill \\ \hfill e_{\eta }^k:=& (-\Delta xsin\left({\theta }_k\right),\Delta ycos\left({\theta }_k\right)),\hfill \end{array}\end{array}$$

(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,{\varphi }_m,h_v),$$

(21)

and then, we write the central finite differences of the second-order horizontal derivatives as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill X_1{X}_{1}U\approx \frac{1}{{(\Delta x)}^2}& (u(q+e_{\xi }^k,{\theta }_k,f_l,{\varphi }_m)-2u(q,{\theta }_k,f_l,{\varphi }_m)\hfill \\ & +u(q-e_{\xi }^k,{\theta }_k,f_l,{\varphi }_m)),\hfill \\ \hfill X_2{X}_{2}U\approx \frac{1}{{(\Delta \theta )}^2}& (u(q,{\theta }_{k+1},f_l,{\varphi }_m)-2u(q,{\theta }_k,f_l,{\varphi }_m)\hfill \\ & +u(q,{\theta }_{k-1},f_l,{\varphi }_m)),\hfill \\ \hfill X_3{X}_{3}U\approx \frac{1}{{(\Delta x)}^2}& \left(u(q+e_{\eta }^k,{\theta }_k,f_l,{\varphi }_m)-2u(q,{\theta }_k,f_l,{\varphi }_m)+u(q-e_{\eta }^k,{\theta }_k,f_l,{\varphi }_m)\right)\hfill \\ \hfill & +\frac{f_l^2}{{(\Delta s)}^2}(u(q,{\theta }_k,f_l,{\varphi }_{m+1})-2u(q,{\theta }_k,f_l,{\varphi }_m)\hfill \\ \hfill & +u(q,{\theta }_k,f_l,{\varphi }_{m-1}))+\frac{fcos\left({\theta }_k\right)}{2\Delta s{\Delta }_x}(u(q+e_{\eta }^k,{\theta }_k,f_l,{\varphi }_{m+1})\hfill \\ \hfill & -u(q-e_{\eta }^k,{\theta }_k,f_l,{\varphi }_{m+1})-u(q+e_{\eta }^k,{\theta }_k,f_l,{\varphi }_{m-1})\hfill \\ \hfill & +u(q-e_{\eta }^k,{\theta }_k,f_l,{\varphi }_{m-1}))-\frac{fsin\left({\theta }_k\right)}{2\Delta s{\Delta }_x}(u(q+e_{\xi }^k,{\theta }_k,f_l,{\varphi }_{m+1})\hfill \\ \hfill & -u(q-e_{\xi }^k,{\theta }_k,f_l,{\varphi }_{m+1})-u(q+e_{\xi }^k,{\theta }_k,f_l,{\varphi }_{m-1})\hfill \\ \hfill & +u(q-e_{\xi }^k,{\theta }_k,f_l,{\varphi }_{m-1}))+\frac{f^2}{2\Delta s}(u(q,{\theta }_k,f_l,{\varphi }_{m+1})\hfill \\ \hfill & -2u(q,{\theta }_k,f_l,\varphi )+u(q,{\theta }_k,f_l,{\varphi }_{m-1})),\hfill \\ \hfill X_4{X}_{4}u[i,j,k,l,m]\approx \frac{1}{{(\Delta f)}^2}& \left(u(q,{\theta }_k,f_{l+1},{\varphi }_m)-2u(q,{\theta }_k,f_l,{\varphi }_m)+u(q,{\theta }_k,f_{l-1},{\varphi }_m)\right).\hfill \end{array}\end{array}$$

(22)

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

$$\begin{array}{cc}\hfill U[i,j,k,l,m,v]=& u(q_{i,1},q_{j,2},{\theta }_k,f_l,{\varphi }_m,h_v)\hfill \\ \hfill =& u(q_{i,1},q_{j,2},{\theta }_k,f_l,{\varphi }_m,h_{v-1})+\Delta t\overline{\mathcal{L}}u(q_{i,1},q_{j,2},{\theta }_k,f_l,{\varphi }_m,h_{v-1}),\hfill \end{array}$$

(23)

where $\overline{\mathcal{L}}$ represents the discretized version of either $\mathcal{L}$ or $\tilde{\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).

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:

Algorithm 1: Completion algorithm pseudocode.

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\times 128$ grayscale images in the experiments related to Figure 10, Figure 11, Figure 12, Figure 13, Figure 14, Figure 15, Figure 16, Figure 17 and Figure 18 and $256\times 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}\}$, $\varphi \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.

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 Figure 11 and Figure 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 Figure 13 and Figure 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 paramaters 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 Figure 13 and Figure 14.

In Figure 15 and Figure 16, we perform the same type of experiments as in Figure 11 and Figure 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.

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 $\mathrm{SE}\left(2\right)$ and it combines the $\mathrm{SE}\left(2\right)$ sub-Riemannian diffusion with a concentration mechanism resulting in a diffusion driven motion by curvature in $\mathrm{SE}\left(2\right)$. The algorithm explained in [48] uses a semidiscrete version of the classical model geometry $\mathrm{SE}\left(2\right)$, 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 Figure 11 and Figure 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$.

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 oputput 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 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, 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 $\mathrm{SE}\left(2\right)$ 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 $\mathrm{SE}\left(2\right)$ [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.