A Cortical-Inspired Contour Completion Model Based on Contour Orientation and Thickness

An extended four-dimensional version of the traditional Petitot–Citti–Sarti model on contour completion in the visual cortex is examined. The neural configuration space is considered as the group of similarity transformations, denoted as M=SIM(2). The left-invariant subbundle of the tangent bundle models possible directions for establishing neural communication. The sub-Riemannian distance is proportional to the energy expended in interneuron activation between two excited border neurons. According to the model, the damaged image contours are restored via sub-Riemannian geodesics in the space M of positions, orientations and thicknesses (scales). We study the geodesic problem in M using geometric control theory techniques. We prove the existence of a minimal geodesic between arbitrary specified boundary conditions. We apply the Pontryagin maximum principle and derive the geodesic equations. In the special cases, we find explicit solutions. In the general case, we provide a qualitative analysis. Finally, we support our model with a simulation of the association field.


Introduction
A mathematical description of the functioning of the human body is a pressing problem in the modern world.Particularly, the specification of cerebration and neuron operation is of particular interest.In this paper, we model the visual information processing by the visual cortex.The complete mechanism of the visual signal processing is not fully studied; however, there is a profound understanding [1] of how image processing is carried out via the information accumulated in light-sensitive receptors, bipolar and ganglion cells of the eye retina.Such information includes the spatial coordinates of the image.
After the retina, the visual signal passes through LGN cells of the thalamus and arrives in the visual cortex.The visual cortex has a multi-layered organization and consists of billions of neural cells.Neurons are connected in a complex network, which is extremely difficult to analyze due to the huge number of elements and even more connections between them.The direct simulation approach to modeling such systems faces inevitable obstacles.However, there are some fundamental principles that are used in network configuration, e.g., the principle of minimum energy spent on establishing communication between two excited neurons of the network.A promising direction for studying such complex systems is to understand such principles and propose simple mathematical models based on these principles.Further mathematical analysis of such models can deepen the understanding of the original systems.
In [2], a mathematical model based on scale space theory was proposed to describe the primary processing mechanism.The model is based on the properties of Gaussian kernels and their derivatives as regularized differential operators, as well as solutions to the linear diffusion equation.The model is supported by experimental study [3], where the shape of the receptive fields (RFs) of bipolar and ganglion cells was established.The authors concluded that RFs are mathematically well approximated by filter profiles based on the Gaussian kernel and Gaussian derivatives.Then, their functioning can be represented as the action of a filter on the input signal (convolution of two functions).
In [4], Marr put forward the idea that retinal ganglion cells perform a convolution capable of extracting qualitative breaks encoded in the signal.Marr also posed that higher levels of visual processing are rooted in this first level of morphological organization of the retinal image, which he called the 2D-primal sketch.He discovered that the convolution of a signal with a receptive profile is a wavelet analysis of the signal, i.e., a spatially localized multiscale Fourier analysis capable of detecting discontinuities.An alternative transformation tool can be an internal compression of the image in accordance with its geometric structure.There is Marr's hypothesis that the image can be reconstructed from its different-scale edges.This contour reconstruction is extremely accurate as only finer details such as textures are smoothed out.In this case, a problem arises with the origin of the breaks.Discontinuities that are robust to large-scale changes can be mistaken for the edges of external objects.Therefore, compression of visual information, i.e., information limitation, is identified with morphological analysis, i.e., geometric constraint.
Further research required a more comprehensive approach.Hubel and Wiesel [5] understood the principles of the primary visual cortex V1 processing.They showed that some parts of the brain react not only to the spatial position of the visible image, but also to its orientation in space.They found that the receptive fields of V1 neurons were elongated rather than rounded.This indicated the ability of V1 cells to detect contour segments with different orientations throughout the image.Mathematically, the operation of V1 cells can be modeled as lifting a two-dimensional input image into an expanded space of positions and orientations.
Even though V1 is physically a two-dimensional neural layer, it implements more than two degrees of freedom.Hubel [1] called the difference between physical and abstract dimensions the grafting of variables.In his pioneering work [6], Hoffmann introduced differential fiber bundles to describe the visual cortex.Here, the base of the fiber bundle represents the retinal plane and the fibers represent the engrafted variables.Further development of this model was performed by Petitot and Tondut [7] (see also the honorable work [8] by Petitot).They described V1 cells as a fiber bundle equipped with a contact structure, and the neural long-range connections were identified as integral curves in the Heisenberg group.In their model, contour completion is carried out along integral curves by minimizing a suitable length functional.
Petitot [8] described the biological functioning of the visual cortex V1 as a sub-Riemannian (SR) structure on the Heisenberg group.Citti and Sarti [9] took into account the nature of the orientation angle and proposed the SR structure on the Lie group SE(2) = R 2 × SO(2).The roto-translation group SE(2) models the configuration space of V1 neurons as the space of R 2 positions and SO(2) orientations.In this model, the reconstruction of a hidden contour occurs by minimizing the excitation energy of neurons that perceive visual information.Such a process is interpreted as the action of the hypoelliptic diffusion operator studied in [10][11][12].Then, the reconstructed parts of the contour are the SR length minimizers on SE (2).The exact expression of the minimizers was found in [13].Such curves are used to render images [14] and to explain some visual illusions [15,16].
Modern computer vision is actively developing based on the principles of biological systems.They are used for image analysis tasks such as enhancement, segmentation, shading, and feature detection.In [17,18], the authors provided a mathematical background for image processing in an extended space of positions and orientations.Salient lines are tracked by SR length minimizers (see [19,20]) or by optimal trajectories in modified models [21][22][23][24], where the spatial propagation along the minimizers is restricted to avoid cusp points and to normalize the curvature of the detected salient lines.
The classic Petitot-Citti-Sarti model has been widely developed.It is the cornerstone of an entire scientific direction-neuromathematics of vision [25,26].This model was modified in subsequent works by many authors (e.g., [27][28][29][30][31][32][33][34][35]), taking into account some aspects of the physiology of vision and based on the needs in the field of image processing.A recent work [36] showed an overview of different models.
The work [27] bridged the classic Petitot-Citti-Sarti model with the edge co-occurrence statistics in natural images.The model's applicability for the association field construction was studied in [28], where the authors showed that the boundary conditions connected by the SR geodesic with cuspless planar projection match the criteria of good continuation.In [29], the authors modified the model to account for the spherical nature of the retina.Later, in [30], this spherical model was extended by considering the non-uniform distribution of photoreceptors on the retina.A link to the well-known Bresloff-Cowan spherical model [31] of V1 hypercolumns was studied in [32].A link to the widely used Wilson-Cowan equations of neural dynamics was studied in [37].In [33], a four-dimensional model accounting for the contour curvature was studied.In [38], the authors proposed a five-dimensional model considering the duration and velocity of visual stimuli.Another five-dimensional extension was proposed in [39], where orientation, frequency, and phase selective behavior of the V1 simple cells are analyzed.Based on this model, a contour completion method was developed in [40].A semidiscrete modification and its application to image processing was studied in [34].
In [35], the authors expanded their classic model by adding a scale parameter and introducing a symplectic structure to describe the structure of neural connectivity.In the current work, we take this model as a basis and develop it by explicitly introducing the length functional and considering the problem of finding length minimizers between arbitrarily specified boundary conditions.
Neurophysiological studies show that spatial hypercolumns also accumulate secondary information about the visible image, such as ocular dominance [41], contour curvature [42,43], contour thickness (scale) [44], and other features.It is known that neurons (simple cells) have different sizes in different areas of the visual cortex.Hence, from V1 to V2, we find simple cells sensitive to objects of different scales (see, e.g., Figure 37 in [25]).In this paper, inspired by [35], we consider a four-dimensional model, which is an extension of the classic Petitot-Citti-Sarti model by adding the contour thickness parameter.The configuration space of neurons is interpreted as a group of similarity transformations M = SIM(2).The left-invariant distribution [45] of the tangent subspaces models the possible directions of establishing a neural connection.The sub-Riemannian distance is proportional to the energy expended in interneuron activation between two excited border neurons.According to the model, the contours of the damaged image are restored using sub-Riemannian geodesics in the space M of positions, orientations and thicknesses (scales).This extension is also intended for image processing tasks to find salient lines (see Figure 1), and to restore damaged image contours (see Figure 2).Image enhancement via left-invariant evolution in the SIM(2) group is studied in [46].In this paper, we consider the problem of SR geodesics in SIM (2).In Section 2, we formulate the model and state the sub-Riemannian problem.In Section 3, we prove the 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 study the Hamiltonian PMP system.We obtain an explicit expression for the abnormal geodesics and provide a qualitative analysis of the Hamiltonian system for normal geodesics.In Section 5, we discuss the boundary value problem.In Section 6, we provide a simulation to construct the association field using the sub-Riemannian geodesics in our problem.
The main contributions of our research are the following: • An optimal control formulation of the contour completion problem; see ( 15) and ( 16).

Problem Formulation
The classic works of Petitot, Citti, and Sarti [8,9] present a model of V1 as a threedimensional Lie group SE(2) of positions and orientations.In their model, horizontal long-range connections between cells of V1 are represented by smooth curves adhering to a nonholonomic constraint: the curves must be tangent to the distribution ∆ = Ker ω, where ω ∈ Λ 1 SE(2) is a given left-invariant differential one-form.Through this approach, a horizontal connection between two neurons is established based on the principle of minimum energy spent on its creation.This leads to the natural modeling of the space V1 by the sub-Riemannian manifold (SE(2), ∆, g), where the metric g specifies the distance encoding the expended energy.According to the model, the visual system performs contour completion (restoration of a corrupted or partially hidden from observation contour) by finding a sub-Riemannian length minimizer between two configurations on the boundary of the damaged area (see Figure 2).
In the consequent work [44], the same authors introduced a new variable σ and defined a symplectic structure in the extended space SIM(2) of positions, orientations and scales.The symplectic structure generates nonholonomic constraints for establishing a long-range neural connection.Note that the explicit form of a metric encoding the energy spent to create the connection was not considered in [44].In our work, we explicitly present this sub-Riemannian metric and formulate an optimal control problem for finding length minimizers.In [44], the authors presented special types of integral curves with a fixed parameter of scale or orientation, which correspond to constant controls in our model.We are also motivated by image analysis applications, where the thickness of contours varies during tracking (see Figure 1).
According to the model, the contour completion mechanism by V1 is invariant under parallel translations, rotations, and scaling of the image on the retina.Such transformations constitute the group of orientation-preserving similarity transformations on the plane The retinal plane is a homogeneous space of the Lie group SIM(2), which acts transitively on it.Thus, SIM(2) models the configuration space of simple cells V1.Now, we explain the lifting of an observable image from the retinal plane to the group SIM (2), which represents V1.The set of receptive profiles of V1 simple cells over a retinal point is formed from the "mother" Gabor function (see [44]), by rotations on an angle θ, and dilations on e σ : where The lifted image O(x, y, θ, σ) : SIM(2) → R + is obtained by probing the observable image I(x, y) : R 2 → R + on the retinal plane with a family of two parametric Gabor filters, A selection from all different cells in a fiber is performed by maximum response selection The values ( θ, σ) are actual values of engrafted variables θ and σ associated with a retinal point (x, y).Such values are used as boundary conditions for contour completion problem that we formulate at the end of this paragraph as an optimal control problem.
Beforehand, we explain a nonholonomic constraint on a long-range (horizontal) neural connection.A horizontal connection is modeled by a smooth curve tangent to the distribution ∆ = Ker ω, where the one-form ω ∈ Λ 1 SIM(2) is given by (see [44]) ω = e −σ (− sin θ dx + cos θ dy).(7) Notice that the distribution ∆ is given by where X i are basis left-invariant vector fields on SIM(2) where Id is a unit element, and L q h = qh is the left translation on SIM(2) (see Appendix A).
By a horizontal curve, we call a Lipschitz curve tangent to ∆ at almost every point where We construct the sub-Riemannian metric by requiring that X 1 , X 3 , and X 4 be orthogonal.Thus, the length of a horizontal curve is given by where the parameters α > 0, β > 0 are coefficients of the sub-Riemannian metric that encode the balance between penalties for motion in a plane along the contour and changing its orientation and thickness.Further, for simplicity, we consider the model case α = β = 1.Any horizontal curve γ(t) of positive length can be reparameterized by arc length ∥ γ(t)∥ = 1 (see Lemma 3.15 in [47]).Thus, the problem of length minimization l(γ) → min is equivalent to time-optimal problem T → min.
Finally, we formulate a contour completion problem as the optimal control problem.Consider the following control system: For given boundary conditions q 0 , q 1 ∈ SIM(2), we aim to find the controls u 1 (t), u 3 (t), u 4 (t) ∈ L ∞ ([0, T], R), such that the corresponding trajectory q : [0, T] → SIM(2) transfers the system from the initial configuration q 0 to the final configuration q 1 by minimum time: Remark 1.The problem is invariant under the left action of SIM(2) since the vector fields X 1 , X 3 , and X 4 are left-invariant.Due to this property, without loss of generality, we set q(0) = Id.

Existence of Solutions
When studying Problems ( 15) and ( 16), the natural question arises about the existence of an admissible trajectory connecting boundary conditions (16).The control system is called completely controllable if an admissible trajectory exists for any q 0 , q 1 ∈ SIM(2).We study the complete controllability of System (15) using the technique of geometric control theory [48].
We have the following nonzero Lie brackets of the controlled vector fields: System ( 15) is symmetric with respect to the controls and it satisfies Hormander condition (i.e., the Lie algebra of the controlled vector fields spans at every point q ∈ SIM(2) the entire tangent space): By Chow-Rashevsky theorem [48], these two conditions plus connectedness of SIM(2) guarantee complete controllability.
Existence of an optimal admissible trajectory that satisfies conditions ( 16) is ensured by Filippov's theorem [48,49].In such a way, we proved the following theorem.
Theorem 1. Solutions to the optimal control problems ( 15) and ( 16) exist for any boundary condition.

Pontryagin Maximum Principle
The necessary condition for optimality is given by the Pontryagin maximum principle (PMP).In this section, we apply PMP to our problem (15) and (16).
The Hamiltonian system The maximum condition The function H(p, q) being maximized is called the Hamiltonian.This is the first integral of the Hamiltonian system.The case H = 0 is called abnormal, and the case H > 0 is called normal.The normal case is reduced to H = 1 by time reparameterization.
The Hamiltonian system (20) in our problems (15) and ( 16) is given by Natural coordinates for left-invariant systems [45] are given by h i (p, q) = ⟨p, X i (q)⟩: The inverse transformation from h to p is given by The Pontryagin function (19) is expressed in the coordinates (23) as follows: The Hamiltonian system (22) takes the following form in the coordinates (23): with the boundary conditions The subsystem for the configuration variables x, y, θ, σ is called the horizontal part, and the subsystem for adjoint variables h i is called the vertical part of the Hamiltonian system.Remark 2. The vertical part of the Hamiltonian system (26) can be alternatively derived by computing the Poisson brackets: ḣi = {H u , h i } (see [47]).
In PMP formulation for Problems (15) and ( 16) of searching for a non-trivial (of non-zero length) optimal curve without loss of generality, we choose the arc length parameterization Next, we consider two cases: H = 0 (the abnormal case), and H = 1 (the normal case).

Abnormal Case H = 0
The Pontryagin function is given by Since the set of admissible controls U : Then, the vertical part of ( 26) is reduced to The horizontal part of ( 26) is reduced to ẋ = ẏ = θ = 0, σ = u 4 .Taking into account the boundary condition x(0) = y(0) = θ(0) = σ(0) = 0 leads to the following theorem.Theorem 2. Abnormal extremal trajectories have the following form: where u 4 (•) is an integrable function with values ±1.
Abnormal optimal trajectories parameterized by arc length have the following form: The expression (31) of optimal trajectories immediately follows from (30) since, in a time-optimal problem, a motion of a system in opposite directions is not optimal.Therefore, the sign of u 4 is not changed.

Normal Case H = 1
The Pontryagin function H u = u 1 h 1 + u 3 h 3 + u 4 h 4 can be considered as a scalar product H u = ⟨(u 1 , u 3 , u 4 ), (h 1 , h 3 , h 4 )⟩.It reaches maximum on the control set U : u 2 1 + u 2 3 + u 2 4 = 1 when the vector (u 1 , u 3 , u 4 ) is collinear to (h 1 , h 3 , h 4 ) and has unit length.Note that, due to the choice of H = 1, this implies the following relation for extremal controls: The Hamiltonian system takes the form (33 This system has a first integral: the Hamiltonian H.We can find more first integrals by considering right-invariant vector fields Y i = R q * A i , R q h = h • q, and associated linear on the fibers of the cotangent bundle Hamiltonians g i = ⟨p, Y i ⟩, p ∈ T * q SIM(2).Among these four right-invariant Hamiltonians, only two (g 1 , g 2 ) are in involution (i.e., the Poisson bracket {g 1 , H} = {g 2 , H} = {g 1 , g 2 } = 0) and functionally independent.Thus, we found the following set of first integrals: To prove Liouville integrability, it is necessary to find four functionally independent first integrals in involution.Three of them we found above.We could not find the remaining first integral.Thus, the question of Liouville integrability remains open.Now, we focus on the vertical part and describe the coadjoint orbits [50].Consider the Poisson bivector, which is given by a matrix P = P ij with the components P ij = {h i , h j }.The structure of Poisson brackets coincides with the structure of Lie brackets (17).Thus, we have only the following nonzero Poisson brackets: We have det P = (h 2 1 + h 2 2 ) 2 ; thus rank P = 0 if h 2 1 + h 2 2 = 0, and rank P = 4 otherwise.In the case h 2 1 + h 2 2 = 0, the coadjoint orbit is zero dimensional.The explicit expression for the extremals are given by the following theorem.

Theorem 3. For the initial covector values
normal extremal trajectories have the following form: They are optimal on a time interval t ∈ [0, π h 30 ], when h 30 ̸ = 0; and up to infinity, when h 30 = 0.
Proof.The expression (37) is obtained by integration of the Hamiltonian system (33) with the boundary conditions ( 27) that, due to the condition (36), is reduced to Optimality for h 30 = 0 holds, since the corresponding extremal trajectory is a straight line passing at maximum speed.In the time-optimal problem, this trajectory is optimal, since any other trajectory of the control system requires more time to reach a point on this line.
Optimality for h 30 ̸ = 0 for t ∈ [0, π h 30 ] holds for the same reason.The trajectory for h 30 ̸ = 0 is not optimal for t > π h 30 , since it has a Maxwell point [45] at t = π h 30 .A Maxwell point is a point where two distinct geodesics meet with the same time.After a Maxwell point, an extremal loses its optimality.A reason for the Maxwell point in our case is periodicity of the angle θ ∈ S 1 .Indeed, consider two trajectories with the initial values h 1 (0) = h 2 (0) = 0, h 3 (0) = ±h 30 , h 4 (0) = h 40 .They reach the same configuration Remark 3. The abnormal optimal trajectories (31) coincide with normal optimal trajectories (37) when h 30 = 0. Thus, they are not strictly abnormal.
In the general case h 2 1 + h 2 2 > 0, the coadjoint orbit is four dimensional.We performed a qualitative analysis of the Hamitonian system leading to the following theorem.
Note that the condition h 40 < 0 is technical, and we use it in the proof.Based on the numerical experiments, we formulate the following conjecture.
Hypothesis 1.Any solution to the vertical part corresponding to the initial covector h 2 10 + h 2 20 > 0 has the following asymptotic behavior:

Approach to the Boundary Value Problem
A geodesic on a sub-Riemannian manifold is a horizontal curve whose sufficiently short arcs are length minimizers.In the optimal control formulation, SR geodesics are the Pontryagin extremal trajectories.Note that PMP is only a necessary, but not sufficient, condition for optimality.It is an infinite-dimensional analogy of the zero derivative condition when minimizing a smooth function in R n .One needs higher-order conditions to find the minimum among all the critical points.The Pontryagin extremals are first-order candidates for being optimal among all admissible trajectories of a control system.Sufficiently short arcs of SR geodesics are optimal, since they satisfy the Legendre condition [48].An extremal trajectory loses optimality at a so-called cut point [47].
There are two types of Pontryagin extremals, called abnormal and normal.In the previous section, we found an explicit expression for the abnormal extremals and derived the Hamiltonian system (33) for the normal extremals in SIM (2).Note that the set of reachable end conditions for arbitrary time (the attainable set) by the abnormal extremals is a one-dimensional subspace in the four-dimensional space SIM(2).In contrast, the attainable set by the normal extremals is the entire group SIM(2), as we proved in Theorem 1.
By varying the initial value of covector h(0) over the set C = {h ∈ R 4 | H = 1}, we obtain a three-parameter family of the normal extremals.Consider a so-called exponential map Exp(h(0), t) : C → SIM(2), which maps the initial covector and the instance of time t > 0 to the point γ(t) of the corresponding geodesic.It holds for general sub-Riemannian manifolds that the exponential map is not injective.For example, consider the initial covector h(0) = (0, 0, 1, 0), then the corresponding extremal trajectory is given by θ(t) = Mod(t, 2π) (see (33)), which is periodic with a period 2π.
Cut points are singularities of the exponential map.There are typically two reasons for a cut point [47]: a conjugate point (a point where the exponential map is degenerate), and a Maxwell point (a point where two distinct geodesics meet at the same time).Finding cut points is a hard mathematical problem, and its solution relies on explicit formulas for the geodesics and analysis of symmetries of the exponential map.Thus, to solve the boundary value problems ( 15) and ( 16) of finding a length minimizer between two given configurations, one needs to restrict the preimage of the exponential map to the domain corresponding to only optimal geodesics.Further, the shooting method can be applied.

Modeling of Association Field
The problem of contour completion (integration) by human visual system was investigated by psychophysicists.Gestalt laws have been proposed for several phenomena of visual perception.Among them, the law of good continuation plays the central role for perceptual completion.The principle of good continuation is found in the experiments of Field, Hayes and Hess [51].Those experiments have resulted in the notion of association field, which describes the set of possible subjective contours starting from a given initial configuration.The role of the scale in contour integration process was stressed in [52].
Inspired by Figure 16 in [51], we provide a simulation of association field by sub-Riemannian geodesics in SIM(2) (see Figure 3).A remarkable property of this model is that the further spatial propagation of the present geodesics does not appear with growing time, which corresponds to Hypothesis 1.This gives a natural bound for the spacial distance between given boundary configurations, which corresponds to the Field model.In Figure 4, we provide a simulation showing that the sub-Riemannian distance in SIM(2) can be used as a criterion for perceptual grouping of the patterns with different positions, orientations, and sizes.In this experiment, we show the points in SIM(2) that are equidistantly (with the distance d = 0.1) distributed along the given three segments of geodesics.They are plotted over the background consisting of points in SIM(2) on a regular spatial (x, y) grid and with randomly chosen orientation θ ∈ S 1 and the scale σ ∈ [0, 1.8].The grid is constructed in a way to guarantee that the distance between the background elements is greater than 0.18.One can see that the equidistantly distributed points are grouped in contrast to the points in the background.
We performed the above simulations in Wolfram Mathematica by numerical integration of the normal Hamiltonian system (33).In the experiments, we relied on the local properties of the exponential map.A further detailed study of the feasibility of SIM(2) model for contour completion on real images requires software to solve a boundary value problem, as discussed in Section 5.This will be a topic of our future research.

Conclusions
In this paper, we considered the sub-Riemannian problem in the Lie group SIM(2) of orientation-preserving similarity transformations of the plane.This problem arises when modeling the mechanism of completing contours by the visual cortex.We considered an extended Petitot-Citti-Sarti model, where the thickness of the contours is taken into account.Based on the principle of minimum energy in the contour completion process, we proposed a sub-Riemannian metric encoding the energy.We stated the contour completion problem as the problem of finding a length-minimizer with given boundary conditions.We reformulated the problem as a time-optimal problem.We proved a solution's existence and applied a necessary optimality condition: PMP.The Hamiltonian PMP system for the geodesics was derived.We found explicit parameterization of abnormal trajectories and provided a qualitative analysis of the normal Hamiltonian system.Finally, we presented simulations on constructing the association field by sub-Riemannian geodesics in SIM (2).

R
Set of real numbers.

Z
Set of integer numbers.S 1  One dimensional sphere (a circle): Lie group of rotations of the plane R 2 .This group is parameterized by angle θ ∈ S 1 .

SE(2)
Lie group of proper motions of the plane R 2 .This group is topologically equivalent to the manifold R 2 × S 1 and parameterized by a vector of parallel translation (x, y) ∈ R 2 and an angle θ ∈ S .

I
Image on the retinal plane.

[•, •]
Commutator (also known as Lie bracket).Commutator of two matrices A i and A j is defined by [A i , A j ] = A i A j − A j A i .Commutator of two vector fields X i , X j ∈ TM at a point q ∈ M is defined by [X i , X j ] q = ∂ ∂s∂t t=s=0 e −tX i • e sX j • e tX i (q), where e tX denotes a flow generated by the vector field X, see [47].sim (2) Lie algebra of the Lie group SIM(2).This is a vector space sim(2) = T Id SIM(2) spanned by the basis vectors A 1 , . . ., A 4 , see (A4) for their matrix representation.Lie Lie algebra generated by the given vector fields and all their commutators, see (18)  Denote by Id the unit element, which is given by identity matrix and corresponds to the origin, and denote by L q the left translation on element q ∈ SIM(2), which is given by matrix multiplication L q q ′ = qq ′ .(A3) The canonical basis for the Lie algebra sim(2) = T Id SIM(2) is given by Every left-invariant vector field is obtained by push-forward of a Lie algebra element under the left translation.The basis left-invariant vector fields are defined as X i (q) = L q * A i (see Figure A1).Explicitly, this construction is the following: Let γ(t) : R → SIM(2) be a smooth curve such that γ(0) = Id, d dt t=0 γ(t) = A i .Then, X i is given by X i (q) = d dt t=0 L q (γ(t)), or X i (q) = qA i in matrix representation (A5) Thus, we have X 1 (q) = e σ cos θ ∂ ∂x + sin θ ∂ ∂y , X 3 (q) = ∂ ∂θ , X 2 (q) = e σ − sin θ ∂ ∂x + cos θ ∂ ∂x , X 4 (q) = ∂ ∂σ . (A6) . Construction of left-invariant vector fields.Here, γ is a smooth curve passing through the unit element Id, L q is the left translation on element q ∈ SIM(2), and A is the tangent vector of γ.

Figure 1 .
Figure 1.Finding blood vessels (salient lines) in the fundus photography of the human retina.Specifications: x, y are spatial coordinates, θ is the orientation, and κ = e σ is the thickness of lines.

Figure 2 .
Figure 2. Restoration of image contours.From left to right: original image; corrupted image (the damaged area is a white disc); recovering contours via the classical model (sub-Riemannian geodesics in SE(2)); restoration via geodesics in SIM(2), taking into account the thickness of contours.

Figure 4 .
Figure 4. Perceptual grouping of the elements in SIM(2) having small SR distance between them.They are plotted over the background consisting of the elements located from each other on a bigger distance.

X 1 (− e σ sin θ − e σ cos θ 0 e σ cos θ − e σ sin θσ cos θ − e σ sin θ 0 e σ sin θ e σ cos θ
The parameters of SIM(2): (x, y) ∈ R 2 is a vector of parallel translation.In our model, (x, y) are coordinates in image plane.θTheparameter of SIM(2): θ ∈ S 1 is an angle of rotation.In our model, θ is the orientation angle between the abscissa and the tangent vector to the contour.
σThe parameter of SIM(2): σ ∈ R is a scaling parameter.In our model, e σ is the thickness of the contour.IdUnit element of the group SIM(2).It is given unit matrix and corresponds to the parameters values x . p Adjoint covector in the Darboux coordinates.pi i-th component of the covector p. h Adjoint covector in the left-invariant coordinates.hi Left-invariant Hamiltonian corresponding to the basis vector field X i , see (19).h i is the i-th component of the covector h.h i0 Initial (for t = 0) value of h i .hi1Terminal (for t = T) value of h i .
Scalar product of a covector and a vector, ⟨∑ p i dx i , ∑ v i∂ ∂x i ⟩ = ∑ p i v i .