1. Introduction
The principles of biological visual systems are of great interest among researchers in many fields of science. An important problem is the development and study of realistic mathematical models describing a certain stage of visual signal processing. In this paper, we perform a study of a mathematical model that describes a mechanism of visual signal processing by the primary visual cortex of the human brain.
Human vision is a complex process that has not yet been fully understood. Perception of visual information starts in the eye. A ray of light hits the retina and causes the activation of light-sensitive receptors. The retina is a part of the brain hierarchically organized in several layers. The first layer consists of light-sensitive receptors. In the next layers of the retina, bipolar and ganglion cells are contained. Here, the primary processing of the visual signal takes place. Furthermore, the visual signal arrives through the optic nerve to the lateral geniculate nucleus (LGN) of the thalamus. Afterward, the signal is transmitted to the visual cortex of the brain.
The mechanism of visual signal processing by neurons of the retina and LGN have been extensively studied. The scale-space theory [
1] proposed a mathematical model for this stage. This theory was inspired by the properties of the Gaussian kernel and its derivatives as regularized differential operators as well as solutions to the linear diffusion equation. The receptive fields of bipolar and ganglion cells as well as LGN neurons are well approximated by filter profiles based on the Gaussian kernel and Gaussian derivatives. Their work is modeled by the action of the filter to the input signal. Mathematically, this operation is defined as the convolution of two functions.
After LGN, the signal is transmitted to the visual cortex of the brain. Physiological research indicates that the visual cortex is composed of multiple layers. The research by Hubel and Wiesel [
2] has made significant progress in understanding the principles of the primary visual cortex V1. In particular, they found that the receptive fields of V1 neurons are elongated rather than rounded. The cells of V1 are capable of detecting segments of contours with different orientations from the whole image. Mathematically, the work of V1 cells can be understood as a lift of a 2D input image to the extended space of positions and orientations
.
The mathematical model of V1 as a sub-Riemannian structure on the Heisenberg group was proposed by J. Petitot [
3]. Then, this model was refined by G. Citti and A. Sarti [
4] as a sub-Riemannian structure on the Lie group
. Here,
models the configuration space of neurons V1, which can be understood as the space of positions
and orientations
. According to this model, the process of completion of occluded contours occurs by minimizing the excitation energy of neurons that perceive visual information in the areas of the observed scene that are occluded. This process can be interpreted as the action of the hypoelliptic diffusion operator studied in [
5,
6,
7]. The resulting curves are sub-Riemannian length minimizers in
. Their exact parameterization was obtained in [
8]. Such curves are used for image inpainting [
9] and for the explanation of some visual illusions [
10].
The principles of biological visual systems are actively used in computer vision. Based on these principles, effective methods of image processing are created: enhancement, segmentation, inpainting, and feature detection. For example, in [
11], the authors describe an approach that is based on lifting an image into the extended space of positions and orientations. After such lifting, the sub-Riemannian length minimizers are used to detect the salient lines [
12,
13,
14]. This approach is actively used in medical image processing.
The experimental data [
15] (see also [
16]) suggest that not only detectors of orientation but also detectors of curvature exist in V1. In this paper, we consider the mathematical model of the visual cortex ([
17], p. 57) [
18] that is obtained by extension of the classical Petitot–Citti–Sarti model. The extension is performed by taking into account the curvature of contours of the observed image. This leads to a sub-Riemannian structure in the four-manifold
, where the fourth component means the curvature. In the four-dimensional model of the visual cortex, a partially occluded contour is completed via the planar projection of a sub-Riemannian length-minimizer in
M that satisfies the boundary conditions (position, orientation, and curvature) obtained from the boundary of the occluded region. See
Figure 1.
In [
18], the authors show that the 4D model allows for better perceptual grouping and completion of complex images than the classical 3D model. They consider the completion of images via minimal surfaces in the 4D model. In this paper, we study the related problem of finding length-minimizers of the associated sub-Riemannian manifold. In analogy with [
9], knowledge of length-minimizers provides a method for the completion of isophotes of corrupted images. Such a method is an alternative to the minimal surface method, proposed in [
18], where the linearized Euler–Lagrange equation is used to compute the surface. Such a linearization provides an approximation for the geodesic flow. We expect that the usage of the precise length-minimizers will provide a more accurate method for image completion since it does not involve linearization.
In this article, we consider the problem of sub-Riemannian geodesics in the space
M of positions, orientations, and curvatures. We formulate the problem as an optimal control problem in
Section 2. In
Section 3, we prove complete controllability of the system and the existence of optimal controls. Then, in
Section 4, we apply a necessary optimality condition—the Pontryagin maximum principle (PMP)—and examine the Hamiltonian system of PMP. We obtain explicit parametrization of abnormal extremals and provide a numerical investigation resulting in the Liouville integrability of the normal Hamiltonian system of PMP. In conclusion, we summarize the main results.
Note that the question of integrability of the geodesic flow in
M is important both for theory and for applications. Liouville integrability ensures the absence of undesirable chaotic behavior. It guarantees that the trajectories of the system remain close to each other under a small perturbation of the initial value. This property is highly important for the stability of numerical schemes of integration. Note that integrability is a rare phenomenon. In the common case, a system of ODE is not integrable. In [
19], the authors show that, even in the simplest case of Carnot groups, the sub-Riemannian geodesic flow is not integrable already in dimension 8 (and consequently in higher dimensions).
3. Complete Controllability of the System
The first question that arises when studying problems (
1)–(
3) is the existence of an admissible trajectory connecting boundary conditions (
3). If for any
,
the answer is positive, then the control system is called completely controllable. Let us investigate the controllability of system (
1) using the technique of geometric control theory [
20].
System (
1) has the following form:
where the vector fields near the controls are given by
We investigate the controllability using Chow–Rashevskii theorem. In our case, it suffices to check that the rank condition is satisfied. To do this, we calculate the following Lie brackets of the fields
:
For the matrix composed of the vector fields
we have
We conclude that the rank of the matrix is four, so the vector fields are linearly independent. Therefore, they define a basis of the tangent space at every point q. Thus, we see that all of the conditions of the Chow–Rashevskii theorem are satisfied and we obtain the following result.
Theorem 1. The control system (
1)
is completely controllable. Remark 2. Since the values are unbounded, condition (
4)
is equivalent to where the family of planes Δ
is called the distribution. Due to condition (
5),
the growth vector of the distribution Δ
equals . Such systems are called structures of Engel type. Furthermore, a question of the existence of optimal trajectories arises: does there always exist an admissible trajectory satisfying conditions (
3), on which the functional (
2) reaches its minimum value? For our problems (
1)–(
3), the answer is positive. The existence of optimal trajectories is guaranteed by the Filippov theorem [
20,
21].
4. Pontryagin Maximum Principle
Before proceeding to the examination of extremal trajectories, let us reduce the problem under consideration to a simpler one. By virtue of Cauchy–Schwarz inequality, the original problem is equivalent to the problem of minimizing the action functional
Apply to problems (
1), (
3), (
6) a necessary condition of optimality—Pontryagin maximum principle (PMP)—we introduce the Pontryagin function:
PMP states that, if , , is an optimal process, then the following conditions hold:
Hamiltonian system ;
Maximum condition ;
Nontriviality condition .
Denote
The Pontryagin function takes the form
In the formulation of PMP, without loss of generality, it suffices to consider two cases: , the abnormal case, and , the normal case. Next, we consider both cases in detail.
4.1. Abnormal Case
The Pontryagin function is This is a linear function, unbounded when . Thus, the maximum condition is satisfied if and only if . The maximized Hamiltonian in this case is
Let be the components of the covector p in canonical Darboux coordinates. They change by the law .
The Hamiltonian system of PMP has the form
The condition
implies
By virtue of system (
7), the second identity implies
, while differentiation of the first identity implies
. Thus, we have
Since every admissible curve of positive length is a Lipschitz reparameterization of an arclength parameterized admissible one (see ([
22], Lemma 3.16)), one can chose the natural parameterization
for every trajectory of system (
1) that is not a fixed point. Thus, identity (
9) is equivalent to
on the intervals of time where the trajectory is not a fixed point.
The identity
implies
. From the third equation of system (
7), it follows that
. On the other hand, the first identity of system (
8) takes the form
. Thus, we have
. Note that, if
, then
, which contradicts the nontriviality condition in PMP. Thus, the case
remains to be considered.
In this case, the covector p is constant and nonzero. We obtain that the abnormal extremals have the form , where , and is any real-valued function. It is easy to see that the trajectory is not optimal if changes its sign. This holds since, when changes its sign, the trajectory is followed in opposite directions. Choosing the natural parameterization on optimal trajectories, we obtain . This results in the optimal trajectory . Thus, we obtain the following result.
Theorem 2. The abnormal extremal trajectories in problems (
1)–(
3)
have the form , where is any real-valued function. Naturally parameterized abnormal optimal trajectories have the form . 4.2. Normal Case
The Pontryagin function
takes the form
The maximum condition gives the expression for the extremal controls:
The maximized Hamiltonian
takes the form
By definition
,
we have
The Hamiltonian system in canonical coordinates has the form
We rewrite this system in the coordinates
:
By choosing the natural parameterization
on the extremal trajectories, one fixes the level surface of the Hamiltonian
. A polar angle
is introduce into the plane
:
By rewriting Hamiltonian system (
10) we obtain the following result.
Theorem 3. The naturally parameterized normal extremal trajectories in problems (
1)–(
3)
are solutions to the system The question of the existence of an analytic expression for extremal trajectories arises: is the Hamiltonian system integrable? To prove the Liouville integrability of system (
10), it suffices to find four functionally independent first integrals in involution. One such integral is the Hamiltonian
H, and two more first integrals
a and
b follow directly from the representation of the Hamiltonian system in canonical coordinates. These three integrals are functionally independent and are in involution. The existence of the remaining first integral is under investigation.
Numerical experiments have been carried out, indicating the presence of the fourth integral. To this end, we consider the four-dimensional system, which is decoupled from the rest of the variables:
Numerical Simulations
The existence of the first integral of system (
12) was studied using the Poincare map. The method consists of three steps. In the first step, one needs to find a periodic trajectory of the system. The second step consists of the construction of a manifold transversal to the flow of the system in a neighborhood of a point of the periodic trajectory. In the third step, small perturbations of the initial point are considered and the corresponding trajectories are computed numerically until they intersect the transversal submanifold several times (
N times). Such points are called orbits of the Poincare map. If the dynamics exhibit chaotic behavior, the trajectories strongly diverge under a small perturbation of the initial point and the orbits form a set of points chaotically distributed in the transversal submanifold. In the case of an integrable system, an orbit forms a set of points lying in a submanifold of a smaller dimension in the transversal manifold. Next, we perform this method for system (
12).
There exists a one-parameter family of periodic trajectories
obtained by changing the initial value
. In this case, the period is
. With such initial values, the point
is fixed. For the next steps, we choose a periodic trajectory passing through the initial point
,
,
,
.
It can be checked that the hyperspace is transversal to the periodic trajectory at the initial point. Indeed, the tangent vector to the periodic trajectory at the initial point is , which is orthogonal to the hyperspace .
New initial points were chosen from a small neighborhood of the initial point of the periodic trajectory. The trajectories departing from these points were computed. Points of intersections of every such trajectory with the transversal hyperspace were computed (for the first time, the second time, …, the N-th time).
In
Figure 2, the orbits of the Poincare map in the space
are shown. The red dot corresponds to the orbit of a periodic trajectory with initial value
,
,
. For the rest of the trajectories, the following parameters were used. For the orange trajectory,
,
,
. For the green trajectory,
,
,
. For the black trajectory,
,
,
. For the blue trajectory,
,
,
. The number of iterations of the Poincare map for all trajectories was chosen as
.
The trajectories were computed by numerical integration (NDSolve in Wolfram Mathematica). The instances when the trajectory intersects the transversal hyperspace were determined using the numerical solution (FindRoot) of the equation along the given trajectory with initial approximation at each iteration .
It is remarkable that the points of the Poincare map fill continuous closed curves. This indicates that, for a small perturbation of the initial conditions, the trajectories remain close over a long time interval. This situation arises when the system is integrable (at least in some domain). Thus, the presented numerical experiments indicate the presence of a fourth independent first integral, which results in Liouville integrability of the Hamiltonian system (
11).
5. Conclusions
In this paper, we consider the problem of sub-Riemannian geodesics in the four-dimensional manifold . This problem arises when modelling the mechanism of occluded contours completion in the extended model of the visual cortex by J. Petitot, G. Citti, and A. Sarti. The extension of the classical 3D model is performed by taking the curvature of the contours into account.
The following main results are obtained in the article. Complete controllability of the system and existence of optimal controls are proven. The Hamiltonian system of PMP is derived. The explicit parametrization of abnormal trajectories is found. In the normal case, three functionally independent first integrals are found. A hypothesis is formulated, confirmed by numerical experiments: the normal Hamiltonian system of PMP is Liouville integrable.
Integrability is an important property of the model. This property indicates the adequateness of the model (the visual signal is processed in a deterministic way). It ensures the stability of numerical methods of integration, which appear in brain-inspired image processing algorithms (such as image enchantment, inpainting, and salient curve extraction).