Association Fields via Cuspless Sub-Riemannian Geodesics in SE(2)

To model association fields that underly perceptional organization (gestalt) in psychophysics we consider the problem P curve of minimizing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\int _{0}^{\ell} \sqrt{\xi^{2} +\kappa^{2}(s)} {\rm d}s $\end{document} for a planar curve having fixed initial and final positions and directions. Here κ(s) is the curvature of the curve with free total length ℓ. This problem comes from a model of geometry of vision due to Petitot (in J. Physiol. Paris 97:265–309, 2003; Math. Inf. Sci. Humaines 145:5–101, 1999), and Citti & Sarti (in J. Math. Imaging Vis. 24(3):307–326, 2006). In previous work we proved that the range \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal{R} \subset\mathrm{SE}(2)$\end{document} of the exponential map of the underlying geometric problem formulated on SE(2) consists of precisely those end-conditions (x fin,y fin,θ fin) that can be connected by a globally minimizing geodesic starting at the origin (x in,y in,θ in)=(0,0,0). From the applied imaging point of view it is relevant to analyze the sub-Riemannian geodesics and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal{R}$\end{document} in detail. In this article we show that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal{R}$\end{document} is contained in half space x≥0 and (0,y fin)≠(0,0) is reached with angle π, show that the boundary \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\partial\mathcal{R}$\end{document} consists of endpoints of minimizers either starting or ending in a cusp, analyze and plot the cones of reachable angles θ fin per spatial endpoint (x fin,y fin), relate the endings of association fields to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\partial\mathcal {R}$\end{document} and compute the length towards a cusp, analyze the exponential map both with the common arc-length parametrization t in the sub-Riemannian manifold \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$(\mathrm{SE}(2),\mathrm{Ker}(-\sin\theta{\rm d}x +\cos\theta {\rm d}y), \mathcal{G}_{\xi}:=\xi^{2}(\cos\theta{\rm d}x+ \sin\theta {\rm d}y) \otimes(\cos\theta{\rm d}x+ \sin\theta{\rm d}y) + {\rm d}\theta \otimes{\rm d}\theta)$\end{document} and with spatial arc-length parametrization s in the plane \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathbb{R}^{2}$\end{document}. Surprisingly, s-parametrization simplifies the exponential map, the curvature formulas, the cusp-surface, and the boundary value problem, present a novel efficient algorithm solving the boundary value problem, show that sub-Riemannian geodesics solve Petitot’s circle bundle model (cf. Petitot in J. Physiol. Paris 97:265–309, [2003]), show a clear similarity with association field lines and sub-Riemannian geodesics.

• show that the boundary ∂R consists of endpoints of minimizers either starting or ending in a cusp, • analyze and plot the cones of reachable angles θ fin per spatial endpoint (x fin , y fin ), • relate the endings of association fields to ∂R and compute the length towards a cusp, • analyze the exponential map both with the common arclength parametrization t in the sub-Riemannian manifold (SE (2), Ker(− sin θ dx + cos θ dy), G ξ := ξ 2 (cos θ dx + sin θ dy)⊗ (cos θ dx + sin θ dy)+ dθ ⊗ dθ) and with spatial arc-length parametrization s in the plane R 2 . Surprisingly, s-parametrization simplifies the exponential map, the curvature formulas, the cusp-surface, and the boundary value problem, • present a novel efficient algorithm solving the boundary value problem, Keywords Sub-Riemannian geometric control · Association fields · Pontryagin's maximum principle · Boundary value problem · Geodesics in roto-translation space

Introduction
Curve optimization plays a major role both in imaging and visual perception. In imaging there exist many works on snakes and active contour modeling, whereas in visual perception illusionary contours arise in various optical illusions [48,52]. Mostly, these optimal curve models rely on Euler's Fig. 1 Stationary curves of the elastica problem ( 0 κ 2 (s) + ξ 2 ds → min) do not need to be global minimizers, cf. [54,66]. E.g. the nondashed elastica is a global minimum (for ξ = 1), whereas in dashed lines we have depicted a local minimum connecting the same boundary conditions elastica curves [33] (minimizing (κ 2 + ξ 2 )ds) to obtain extensions where typically external forces to the data are included, cf. [5,18,21,60,61]. The elastica problem suffers from the well-known fact that not every stationary curve is a global minimizer, e.g. many local minimizers exist, cf. Fig. 1. Stationarity of a curve can be reasonably checked by the visual system using local perturbations, whereas checking for (global) optimality [54,66] is much more difficult. Some visual illusions (e.g. the Kanisza triangle) involve corners requiring abrupt resetting of initial and ending conditions, which are difficult to explain in the elastica model. Another problem with elastica is that it is very hard to solve the boundary value problem analytically [4,6] (due to a highly non-linear ODE for curvature [48]) and this requires efficient numerical 3D shooting schemes.
On top of that elastica curves relate to modes of the direction process (for contour-completion [24]) where the direction of an oriented random walker is deterministic and its orientation is random. Such deterministic propagation only makes sense when the initial orientation is sharply defined. Instead Brownian motion with random behavior both in spatial propagation direction and in orientation direction [1,22,25], relates to hypo-elliptic diffusion on the planar roto-translation group. Such a Brownian motion models contour enhancement [25] rather than contour completion [24], see [28] for a short overview. The corresponding Brownian bridge measures [27,67] (relating to so-called completion fields in imaging [4,24,63,64]) tend to concentrate towards optimal sub-Riemannian geodesics [12,15,22,26,47,56]. So both elastica curves and sub-Riemannian geodesics relate to two different fundamental left-invariant stochastic processes [28] on sub-Riemannian manifolds on the 2D-Euclidean motion group SE (2), (respectively to the direction process [24,48] and to hypo-elliptic Brownian motion [1,22,25]).
In short, advantages of the sub-Riemannian geodesic model over the elastica model are: Fig. 2 An example of a smooth sub-Riemannian geodesic γ = (x(·), y(·), θ (·)) (in purple) in auxiliary problem P MEC , Eq. (12), whose spatial projection (in black) shows a cusp (red point). A cusp point is a point (x, y, θ) on γ such that the velocity (black arrow)ẋ of the projected curve x(·) = (x(·), y(·)) switches sign at (x, y). At such a point in SE(2) ≡ R 2 S 1 the tangent vector points (blue arrow) in θ -direction (Color figure online) • Every cuspless sub-Riemannian geodesic (stationary curve) is a global minimizer [15,16]. • The Euler-Lagrange ODE for normalized curvature z = κ/ κ 2 + ξ 2 can be reduced to a linear one. • The boundary value problem can be tackled via effective analytic techniques. • The locations where global optimality is lost can be derived explicitly. • Sub-Riemannian geodesics are parametrization independent in the roto-translation group SE (2), which is encoded via a pinwheel structure of cortical columns in the primary visual cortex [50,51]. However, the practical drawback of sub-Riemannian geodesics compared to elastica is that their spatial projections may exhibit cusps and it is hard to analyze when such a cusp occurs. See Fig. 2. Therefore, in this article we provide a complete analysis of such sub-Riemannian geodesics, their parametrization, solving the boundary value problem, and we show precisely when a cusp occurs. See Fig. 3.
This variational problem was studied as a possible model of the mechanism used by the visual cortex V1 to reconstruct curves which are partially hidden or corrupted. This model was initially due to Petitot (see [50,51] and references therein). Subsequently, the sub-Riemannian structure was introduced in the problem by Petitot [52] for the contact geometry of the fiber bundle of the 1-jets of curves in the plane (the polarized Heisenberg group), whereas Citti and Sarti [22] introduced the sub-Riemannian structure in SE (2) in problem P. The group of planar rotations and translations SE(2) is the true symmetry group underlying problem P. Therefore, we build on the SE(2) sub-Riemannian viewpoint first proposed by Citti and Sarti [22], and we solve their cortical model for all appropriate end-conditions. The stationary curves of problem P were derived by the authors of this paper in [12,26]. The problem was also studied by Hladky and Pauls in [40], and by Ben-Yosef and Ben-Shahar in [11].
In this article we will show that the model coincides 1 with the circle bundle model by Petitot [52] and that its min-Remark 1.1 Problem P is well-posed if and only if, 2 where R θ in denotes the counterclockwise rotation over θ in in the spatial plane and where R is a particular subset R 2 × S 1 (equal to the range of the underlying exponential map of P curve which we will define and derive later in this article), cf. [15,16].
We will see in the following that this set R is the set of all endpoints in R 2 × S 1 that can be connected with a cuspless stationary curve of problem P, starting from (0, 0, 0).

Remark 1.2 The physical dimension of parameter ξ is
[Length] −1 . From a physical point of view it is crucial to make the energy integrand dimensionally consistent. However, the problem with (x(0), θ (0)) = (0, 0, 0) and ξ > 0 is equivalent up to a scaling to the problem with ξ = 1: The minimizer x of P with ξ > 0 and boundary conditions (0, 0) and (x 1 , θ 1 ) relates to the minimizer x of P with ξ = 1 and boundary conditions (0, 0) and (ξ x 1 , θ 1 ), by spatial rescaling: x(s) = ξ −1 x(s). Therefore, in the remainder of this article we just consider the case ξ = 1 for simplicity.
It is not straightforward to derive the exact Euler-Lagrange equations together with a necessary geometric study of the set of all possible solution curves. The exact solutions to the problem can be derived using 3 types of techniques: 1. Direct derivation of the Euler-Lagrange equation. E.g.
the approach by Mumford [48], yielding a direct approach to the ODE for the curvature, see Appendix A. 2. The Pontryagin Maximum principle: A geometrical control theory approach based on Hamiltonians, cf. [3,12,47,53] and Appendix D. 3. The Bryant and Griffith's approach (based on the works by Marsen-Weinstein on reduction in theoretical mechanics [44]) using a symplectic differential geometrical approach based on Lagrangians [26, App. A], cf. [19].
In this article we will apply all three techniques as they are complementary. Furthermore, we aim to provide a complete overview on the surprisingly tedious problem (many inaccurate and/or incomplete results on the stationary curves have appeared in the mathematical imaging literature). Finally, we want to connect remarkably different approaches in previous works [11,14,22,26,47,58] on the topic. The first approach very efficiently produces only the Euler-Lagrange equation for the curvature of stationary curves, but lacks integration of a single curve and lacks a geometric study of the continuum of all stationary curves that arise by varying the possible boundary conditions.
The second approach includes profound geometrical understanding from a Hamiltonian point of view and deals with local optimality [3] of stationary curves.
The third approach 3 takes a Lagrangian point of view and provides additional differential geometrical tools from theoretical mechanics that help integrating and structuring the canonical equations. These additional techniques will be of use in deriving semi-analytic solutions to the boundary value problem and in the modeling of association fields.
All three approaches provide, among other results, the following linear hyperbolic ODË for normalized curvature where s denotes spatial arc-length and κ(s) denotes curvature of the spatial part r → x(r) of a geodesic γ = (x, θ) : [0, ] → R 2 S 1 , with θ(s) = arg(ẋ(s) + iẏ(s)). Such geodesics are globally minimizing, cf. [15,16] and Theorem 1 below). Furthermore, denotes sub-Riemannian arclength t as a function of s along a sub-Riemannian geodesic. Recall that spatial arclength s and sub-Riemannian arclength t are respectively determined by As a particular case of Eq. (6), the total sub-Riemannian arc-length T of the lifted curve s → γ = (x(s), θ (s)) with θ(s) = arg(ẋ(x) + iẏ(s)), relates to the total length of the spatial curve s → x(s) via T = t ( ). Firstly, application of Mumford's approach for deriving the ODE for curvature of elastica, to problem P is relatively straightforward, see Appendix A, but does not explicitly involve geometrical control and the Frenet formula still needs to be integrated.
Secondly, in our previous work [16] we considered an extended mechanical problem P MEC related to P. This problem P MEC will soon be explained in detail in Sect. 1.1, and is completely solved by Sachkov et al. in [47,55,56]. Application of the Pontryagin maximum principle to this related problem P MEC (after squaring the Lagrangian and constraining the total time to a fixed 4 T ) yields for ξ = 1 the maximized Hamiltonian 5 with momentum p = p 1 dθ + p 2 dx + p 3 dy and the induced canonical equations which via re-parametrization of cylinder H (p) = 1 2 sin(ν/2) = p 2 cos(θ ) + p 3 sin(θ ), produces the mathematical pendulum ODË For details on the involved computation see [16,47]. Thirdly, application of the Bryant and Griffith's (Lagrangian) approach to problem P will yield a canonical Pfaffian system on an extended manifold whose elements involve both position, orientation, control (curvature and length), spatial momentum and angular momentum. We will show that the essential part of this Pfaffian system is equivalent to ∇γ p = 0 where ∇ denotes a Cartan connection and p denotes momentum as a co-vector within T * (R 2 S 1 ). This fundamental identity allows us to analytically solve the boundary value problem.
1.1 Lift problem P to the roto-translation group Problem P relates to two different geometric control problems (P curve and P MEC ): • P curve : Fix ξ > 0 and boundary conditions (x in , y in , θ in ), In the space of integrable (possibly nonsmooth) controls v(·) : [0, ] → R, we aim to solve: → min (here T ≥ 0 is free) Problem P MEC has a solution by Chow's and Fillipov's theorems [3] regardless the choice of end-condition and has been completely solved in a series of papers by one of the authors (see [47,55,56]). It gives rise to a sub-Riemannian distance on the sub-Riemannian manifold within SE(2) as we will explain next.
The space R 2 × S 1 can be equipped with a natural group product where R θ denotes a counter-clockwise rotation over angle θ ∈ (−π, π] and with x = (x, y) T and x = (x , y ) T so that it becomes isomorphic to the 2D (special) Euclidean motion group consisting of rotations and translations in the plane, also known as roto-translation group, and commonly denoted by SE (2). As SE(2) acts transitive and free on the set of positions and orientations R 2 × S 1 we can identify point on orbits (x, y, θ) starting from the unity (0, 0, 0) with the corresponding group elements (x, y, R θ ). Therefore we write R 2 S 1 ≡ SE (2) to stress that the set R 2 × S 1 is equipped with a (semi-direct) group product (13). Now both problems P curve and P MEC are invariant with respect to rotations and translations so we may as well set (x in , y in , θ in ) = (0, 0, 0). Indeed, given a problem with general boundary conditions (x in , y in , θ in ) and (x f in , y f in , θ f in ), its minimizer γ opt (when it exists) is (x in , y in , θ in ) ·γ opt , whereγ opt is the minimizer from (0, 0, 0) to Throughout this article we use the following notation for the moving frame {A 1 , A 2 , A 3 } of left-invariant vector fields where on the right we consider vector fields as differential operators, for details on such identification see e.g. [3,7]. The corresponding co-frame of left-invariant dual basis vectors will be denoted bŷ X 1 = (0, 0, 1) ↔ ω 1 := dθ, X 2 = (cos θ, sin θ, 0) ↔ ω 2 := cos θ dx + sin θ dy, where frame and dual frame relate viâ where in the righthand side we have the Kronecker symbols δ i j = 1 if i = j and 0 else. Problem P MEC can now be reformulated as the computation of where d denotes the sub-Riemannian distance 6 on the sub-Riemannian manifold with sub-Riemannian metric tensor Problem P MEC is to be considered as an auxiliary mechanical problem (of optimal path planning of a moving car carrying a steering wheel and the ability to drive both forwardly and backwardly) associated to P curve . To this end we stress that P MEC cannot be interpreted as a problem of reconstruction of planar curves, [14]. The problem is that the minimizing curve γ = (x, θ) : [0, T ] → SE(2) may have a vertical tangent vector (i.e. in θ -direction) in between the ending conditions, which causes a cusp in the corresponding projected curve t → x(t) in the plane, see Fig. 2. Such a cusp corresponds to a point on an optimal path where the car is suddenly set in reverse gear.
Problem P MEC is invariant under monotonic re-parameterizations and at a cusp spatial arc-length parametrization breaks down. If (x f in , y f in , θ f in ) ∈ R no such cusps arise and P MEC and P curve are equivalent [15,16] and we can use arclength parametrization also in P MEC (in which case the first control-variable is set to 1, since ω 2 | γ (s) ,γ (s) = 1). In [16] we have proven the following Theorem.

Definition 1
Let R ⊂ SE(2) denote the set of end-points in SE(2) that can be reached from e with a stationary curve of problem P curve .

Theorem 1
In P curve we set initial condition (x in , y in , θ in ) = e = (0, 0, 0) and consider (x f in , y f in , θ f in ) ∈ R 2 S 1 . Then • (x f in , y f in , θ f in ) ∈ R if and only if P curve has a unique minimizing geodesic which exactly coincides with the unique minimizer of P MEC .
∈ R if and only if problem P curve is illdefined (i.e. P curve does not have a minimizer). 7 6 Usually the minimization in Eq. (16) is made in the space of Lipschitz functions, to guarantee the existence of minimizers via PMP. However, a posteriori one verifies that these minimizers are indeed C ∞ . 7 For end-condition (x f in , y f in , θ f in ) / ∈ R problem P MEC has a minimizer with internal cusp (and thereby violating the natural settings of P curve ). Such a minimizer of P MEC can be approximated by smooth curves satisfying the constraints of problem P curve . In these cases P curve does not allow local or global minimizers, nor does it allow a stationary curve [16]. Cuspless sub-Riemannian geodesics (projected on the plane) for admissible boundary conditions modeling the association field as in Fig. 8. According to Theorem 1 they are global minimizers. Remarkably the tangent vector to these geodesics (e.g. the red geodesic) is nearly vertical at the end condition and the large curvature at the end condition at the association field, indicate the association field lines end at close vicinity of cusps (Color figure online) As a result, for the case g in = (0, 0, 0), we say g f in ∈ SE (2) is an admissible end-condition for P curve if g f in ∈ R, as only for such end-conditions we have existence of a (smooth) global minimizer, see also [12]. See Fig. 4.

Structure of the Article
Firstly, in Sect. 3 we consider the origin of the problem of finding cuspless sub-Riemannian geodesics in (SE (2), , G β ), which includes cortical modeling of the primary visual cortex and association fields.
In Sect. 4 we provide a short road map on how to connect two natural parameterizations. The cuspless sub-Riemanian geodesics in the sub-Riemannian manifold (SE (2), , G β ) can be properly parameterized by the sub-Riemannian arclength parametrization (via t) or by spatial arclength parametrization (via s). Parametrization via t yields the central part of the mathematical pendulum phase portrait (recall Eq. (10)), whereas parametrization via s yields a central part of a hyperbolic phase portrait (recall Eq. (4)). The hyperbolic phase portrait does not coincide with a local linearization approximation (as in Hartman-Grobman's theorem [38]). In fact, it is globally equivalent to the relevant part of the pendulum phase portrait (i.e. the part associated to cuspless sub-Riemannian geodesics). The involved coordinate transforms are global diffeomorphisms.
In Sect. 5 we define the exponential map [2,47] for P curve and P MEC . Then we show that the set R ⊂ SE(2) (consisting of admissible end-conditions) equals the range of the exponential map of P curve . We will provide novel explicit formulas for the exponential map for P curve using spatial arc length parametrization s and moreover, for completeness and comparison, in Appendix B we will also provide explicit formulas for the exponential map of P MEC that were previously derived in previous work [47] by one of the authors.
We show that the exponential map of P curve follows by restriction of P MEC to the strip (ν, c) ∈ [0, 2π] × R, see Fig. 9. A quick comparison in Appendix B learns us that spatial arclength parametrization (also suggested in [22]) simplifies the formulas of the (globally minimizing, cuspless) geodesics of P curve considerably.
As the set of admissible end-conditions R equals the range of the exponential map of P curve , we analyze this important set R carefully in Sect. 6. More precisely, we 1. show that R is contained in half space x ≥ 0 and (0, y fin ) = (0, 0) is reached with angle π , 2. show in Theorem 6 that the boundary ∂R consists of the union of endpoints of minimizers either starting or ending in a cusp and a vertical line l above (0, 0, 0), and we compute the total spatial arc-length towards a cusp, 3. analyze and plot the cones of reachable angles θ fin per spatial endpoint (x fin , y fin ), 4. prove homeomorphic and diffeomorphic properties of the exponential map in Theorem 6, 5. show in Lemma 8 that geodesics that end with a cusp at θ f in = π 2 are precisely those with stationary curvature (κ(0) = 0) at the origin.
In Sect. 7 we solve the boundary value problem, where we derive a (semi)-analytic description of the inverse of the exponential map and present a novel efficient algorithm to solve the boundary value problem. This algorithm requires numerical shooting only in a small sub-interval of [−1, 1], rather than a numerical shooting algorithm in R 2 × S 1 .
In Sect. 8 we show a clear similarity of cuspless sub-Riemannian geodesics and the association field lines from psychophysics [34] and neuro-physiology [52]. This is not surprising as we will show that sub-Riemannian geodesics allowing x-parametrization, exactly solve the circle bundle model for association fields by Petitot, cf. [52]. It is remarkable that the endings of association fields are close to the cusp-surface ∂R, which we underpin with Lemma 8 and Remark 8.1.
For a concise overview of previous mathematical models for association fields and their direct relation to the cuspless sub-Riemannian geodesic model proposed in this article we refer to the final subsection in Appendix G.

Origin of Problem P: Cortical Modeling
In a simplified model (see [51, p. 79]), neurons of V1 are grouped into orientation columns, each of them being sensitive to visual stimuli at a given point of the retina and for a given direction on it. The retina is modeled by the real plane.
Orientation columns are connected between them in two different ways. The first kind is given by vertical connections, which connect orientation columns belonging to the The human visual system not only performs a score of local orientations (organized by a pinwheel structure in V1). It also checks (a priori) for alignment of local orientations in the enhancement and detection of elongated structures. In modeling both procedures it is crucial that one does not consider R 2 × S 1 as a flat Cartesian space. See Fig. 7.
The Euclidean motion group acts transitively and free on the space of positions and orientations, allowing us to identify the coupled space of positions and orientations R 2 S 1 with the roto-translation group SE(2) = R 2 SO(2). This imposes a natural Cartan connection [26,52] on the tangent Positions and orientations are coupled. The spatial and angular distance between (x 1 , θ 1 ) and (x 0 , θ 0 ) is the same as the spatial and angular distance of (x 2 , θ 1 ) between (x 0 , θ 0 ). However, (x 1 , θ 1 ) is much more aligned with (x 0 , θ 0 ) than (x 2 , θ 1 ) is. The left-invariant sub-Riemannian structure on the space R 2 S 1 takes this alignment into account. The connecting curves are spatial projections of sub-Riemannian geodesics in SE (2) Besides the non-commutative group structure on R 2 S 1 ≡ SE(2), contact geometry plays a major role in the functional architecture of the primary visual cortex (V1) [41], and more precisely its pinwheel structure, cf. [52]. In his paper [52] Petitot shows that the horizontal cortico-cortical connections of V1 implement the contact structure of a continuous fibration π : R × P 1 → P 1 with base space the space of the retina and P 1 the projective line of orientations in the plane. He applies his model to the Field's, Hayes' and Hess' physical concept of an association field, to several models of visual hallucinations [32] and to a variational model of curved modal illusory contours [42,48,65]. Such association fields reflects the propagation of local orientations in the primary visual cortex. For further remarks on the concept of an association field and its mathematical models see Appendix G. Intuitively, the tangents to the field lines of the association field provide expected local orientations, given that a local orientation is observed at the center of the field in Fig. 8). These association fields have been confirmed by Jean Lorenceau et al. [43] via the method of apparent speed of fast sequences where the apparent velocity is overestimated when the successive elements are aligned in the direction of the motion path and underestimated when the motion is orthogonal to the orientation of the elements. They have also been confirmed by electrophysiological methods measuring the velocity of propagation of horizontal activation [37]. There exist several other interesting low-level vision models and psychophysical measurements that have produced similar fields of association and perceptual grouping [39,49,68], for an overview see [52,Chaps. 5.5,5.6]. Remarkably, psychological physics experiments based on multiple Gabor patch-stimuli indicate a thresholding effect in contour recognition, if the slope variation in two subsequent elements (Gabor patches) is too large no alignment is perceived and if the orientations are no longer tangent but transverse to the curve no alignment is perceived, cf. [52].
In this article we will show that sub-Riemannian geodesics closely model the association fields from psychophysics and Fig. 8 Modeling the association field with sub-Riemannian geodesics and exponential curves, (a) the association field [34,52]. Compare the upper-right part of the association field to the following lines: in (b) we impose the end condition (blue arrows) for the SR-geodesic model in black and the end condition (red arrows) for the horizontal exponential curve model [57], Eq. (72), in grey; (c) comparison of sub-Riemannian geodesics with exponential curves with the same (co-circularity) ending conditions; (d) as in (b) including other ending conditions (Color figure online) that the location of cusps seems to provide a reasonable grouping criterium to connect two local orientations (consistent with endings of the association field), see Fig. 8. Next we will show that it does not matter whether one lifts problem P (given by Eqs. (1) and (2)) to the projective line bundle or to the group of rotations and translations in the plane.

No Need for Projective Line Bundles in P curve
The can as well be formulated on the projective line bundle P 1 [14,52] where antipodal points on the sphere S 1 are identified. See also [13].
In the setting of P curve , we then can study the problem with initial condition in the set and similarly for the final condition. Nevertheless, the structure of solutions does not change with respect to the solutions of the standard problem P curve . Indeed such flips are either not allowed or they do not produce new curves: • Flipping only one of the boundary conditions is not possible as in this article we shall show that if ( is not admissible. • If we both flip (i.e. θ → θ + π ) and switch both the initial and ending condition we get the same curve (in opposite direction).
So when insisting on cuspless solution curves in our central problem P, lifting problem P to the projective bundle R 2 P 1 is equivalent to lifting P to SE(2) ≡ R 2 S 1 .
In fact, identification of antipodal points does not make any difference when considering cuspless sub-Riemannian geodesics in (SE (2), , G ξ ). Therefore, in this article we will not identify antipodal points and we focus on problem P curve and its corresponding admissible boundary conditions (i.e. an explicit description of the set R ⊂ SE (2)).

Parametrization of Curves in P curve
The natural parametrization for sub-Riemannian geodesics in P MEC is the sub-Riemannian arclength parametrization. However, when considering only those sub-Riemannian geodesics in (SE (2), , G ξ ) without cusps (as in P curve ), i.e. the cuspless sub-Riemannian geodesics, the problem is actually a planar curve problem (as in P) and there it is more natural 8 to use spatial arclength parametrization.

Cusps and the Exponential Map Associated to P curve and P MEC
In order to express the exponential map associated to P curve (for ξ = 1) in spatial arclength parametrization we apply Bryant & Griffith's approach [20], which was previously successfully applied to the elastica problem [19]. Here we will also include an additional viewpoint on this technical approach via the Cartan connection. In case the reader is not so much interested in the geometrical details and underpinnings, it is also possible to skip the following derivations and to continue reading starting from the formulas for the sub-Riemannian geodesics γ (s) in Theorem 3.
To avoid large and cumbersome computations we first need some preliminaries on moving frames of references and Cartan connections. Recall to this end our notations for left-invariant frame {A i } 3 i=1 given by Eq. (14), and leftinvariant co-frame {ω i } 3 i=1 given by Eq. (15). The leftinvariant vector fields generate a Lie algebra where the non-zero structure constants are c 3 12 = −c 3 21 = −c 2 13 = c 2 31 = 1. This Lie-algebra serves as the moving frame of reference in R 2 S 1 ≡ SE(2). The Cartan connection ∇ on T (SE (2)) is given by where we used the following definitionṡ As a result (for details see Eq. (93) and Theorem 12 in Appendix C) covariant differentiation of a momentum covector field along a curve γ : [0, ] → SE (2) yields withλ k (s) = dλ k ,γ (s) . Finally we mention the Cartan's structural formula Step 1: Extend the manifold SE(2) with geometric control variables Consider the extended manifold Q = SE(2) × R + × R × R with coordinates (x, y, e iθ , σ, κ, r), where σ = x (r) so that ds = σ dr, where r → x(r) is some parametrization of the spatial part of the lifted curve r → γ (r) = (x(r), θ (r)) in SE (2). In order to extend the sub-Riemannian manifold (SE (2), Ker(ω 3 ), G ξ =1 ) such that the concept of horizontal curves is preserved we impose These equations determine the horizontal part Step 2: Include momentum Include the Lagrange multipliers as local momentum vectors in our target space. Therefore we extend Q to a larger space Z. We define Z as the affine sub-bundle Step 3: Minimization on extended space Z Consider a one parameter family {N r } of horizontal vector fields on SE (2) and compute the variation of the integrated Lagrangian-form ψ along such a N r : where we used the Stokes Theorem N r d( ∂ ∂r ψ) = ∂ Nr ∂ ∂r ψ = 0 and the formula for Lie derivatives of volume forms along vector fields L X A = X dA + d(X A) and where X A := A(X, ·) denotes the insert operator. Consequently, we must solve the canonical ODE system Γ (r) dψ| Γ (r) = 0 for all r > 0. (26) where Γ (r) ≡ (γ (r), κ(r), σ (r), r, p(r)). This boils down to Now by means of the Cartan structural formula (22), and Eq. (27) we obtain the Pfaffian system The first three equations represent the horizontality restriction. The two equations in the middle represent the Euler-Lagrange optimization of the energy and show that {λ 1 , λ 2 , λ 3 } are components of momentum with respect to the dual frame (under identification (24)). It is readily deduced that with γ a cuspless sub-Riemannian geodesic can be rewritten as where ∇ denotes the Cartan connection on the co-tangent bundle T * (SE (2)).
Proof The last 3 equations in (28) provide the momentum covector. They can be written as which by Eq. (21) can be rewritten as To this end we note that Finally, with respect to the second part of Eq. (30): from which the result follows.
The second part allows us to interpretate p = 3 i=1 λ i ω i as a momentum covector.

Remark 5.3
In contrast to Levi-Civita connections on Riemannian manifolds, the Cartan connection ∇ has torsion and thereby auto-parallel curves do not coincide with geodesics. In fact, Theorem 12 in Appendix C shows that auto-parallel curves are (horizontal) exponential curves.
Step 4: Integrate the Pfaffian system To integrate ∇γ p = 0 we resort to matrix-representation m : SE(2) → R 3×3 given by and express dual-vectors (covectors) as row vectors. Analogously to Bryant's work on elastica [19] we express equation (32) in explicit coordinates where we use short-notation for the row-vector from which we deduce that Before we will derive γ from Eq. (37) we will need the following lemma based on Noether's theorem. Formally, one can avoid this general abstract lemma (as in [19]) by observing Lemma 2 Cuspless sub-Riemannian geodesics are contained within the co-adjoint orbits for all s ∈ [0, s max ], with s max given by Eq. (41).
Proof According to Noether's theorem (i.e. conservation law on momentum) the moment map m : (2)), for all Ξ ∈ T (SE (2)) is constant along the characteristic curves Ξ =γ . The co-adjoint representation of SE (2) acting on the dual of its Lie-algebra (T (SE (2))) * is given by We have m(η g (q, p)) = (Ad g −1 ) * m(q, p), where the group action g → η g is given by As a result the co-adjoint orbits of SE(2) coincide with the cylinders in Eq. (38).
The minimizers of P curve are cuspless geodesics and their total length (towards a cusp) equals The curvature of orbits with c < 1 and z 0 > 0 is strictly positive. The curvature of orbits with c < 1 and z 0 < 0 is strictly negative. The curvature of orbits with c > 1 switches sign once at Proof Follows directly from the hyperbolic phase portrait induced byz = z and Theorem 2, and solving for respectively |z(s)| = 1 and z(s) = 0.
After these results on sub-Riemannian geodesics, we continue with solving for ∇p = 0, Eq. (37). Problem P curve is left-invariant and in the next lemma we select a suitable point on each co-adoint orbit to simplify the computations considerably.
Applying the above Lemma and Eq. (29) provides the next theorem, Theorem 3, where we provide explicit analytical formulae for the geodesics by integration of the Pfaffian system. To this end we first need a formal definition of the operator that integrates the Pfaffian system Eq. (28) and produces the corresponding geodesic of P curve in SE (2).
This operator needs initial momentum p 0 and total spatial length > 0 as input and produces the corresponding geodesic of P curve as output. By Eqs. (29) and (32) initial momentum equals with initial normalized curvature z 0 = κ 0 / κ 2 0 + 1. As a result, we have The Hamiltonian at the unity element, evaluated at initial momentum is given by Now let us use arclength parameterization (so set r = s and σ = 1) in the canonical ODE system (26) on Z. Via identification Eq. (24) this gives rise to an equivalent ODE system on Q × T * (SE (2)) with unity element e = (0, 0, 0) ∈ SE(2).
Remark 5.4 For sober notation we omit index e and write Exp = Exp e and H (p) = H (e, p) for exponential map and Hamiltonian. Furthermore, we include a tilde in this exponential map associated to the geometrical control problem of P curve to avoid confusion with the exponential map Exp : T e (SE(2)) → SE(2) from Lie-algebra to Lie group.  (2) propagate only in vertical direction, not allowing spatial arc-length parameterization. See also [16,Remark 31].
Note that the cuspless geodesic γ follows from cuspless geodesicγ = h 0 γ via the rigid body motion

Corollary 2
The end-point g f in of a cuspless sub-Riemannian geodesic is given by Proof From the previous Theorem 3 we deduce from which the result follows.

Corollary 3
The (x, y)-coordinates of the Exponential map involve one elliptic integral and the tangent vectors along geodesics do not involve any special functions. Furthermore, from −ẏ(s) ≥ 0 it follows that the spatial part of the geodesics is monotonically increasing along the Geodesics with c = 1 admit simple formulas:
For a plot of the critical surface see Fig. 10. In Theorem 3 we have derived the exponential map of P curve in terms of spatial arc-length parametrization s, whereas in previous work [15] the exponential map of P MEC is expressed in sub-Riemannian arc-length t. For comparison see Appendix B.
On the one hand one observes that the exponential map of P curve is much simpler when expressed in s and it is easier to integrate in current active shape models in imaging where the same kind of parametrization is used. On the other hand for P MEC it is more natural to choose t-parametrization as this parametrization does not beak down at cusps. The following theorem relates the exponential mappings for P curve and P MEC .

Theorem 4 Let
Exp denote 9 the exponential map of P curve . Let EXP denote the exponential map of P MEC . Then these exponential maps satisfy the following relation for all p 0 ∈ C ⊂ T * e (SE (2)), and all 0 < ≤ s max , (so that (p 0 , ) ∈ D, recall Eq. (46)), where t ( , p 0 ) is given by Eq. (6).
Proof We note that ≤ s max implies that the orbits do not hit the cusp lines in the pase portraits (i.e. |z| = 1 and ν = 0, 2π ) so that (ν(t), c(t)) stays within the central strip (i.e. ν(t) ∈ [0, 2π]) indicated in Fig. 9. The rest follows by Lemma 1. 9 For the sake of simplicity we do not index Exp the exponential map with the initial condition g in , as throughout this article we set g in = e = (0, 0, 0).

The Set R and the Cusp-Surface ∂R
According to Theorem 1 the set of points in SE(2) that can be reached with a global minimizer from unity element g in = e = (0, 0, 0) is equal to R given in Definition 1. Therefore, we first need to investigate this set in order to apply cuspless sub-Riemannian geodesics in vision applications. First of all we have the following characterization.
Theorem 5 Let s max (p 0 ) be given by Eq. (41). Let C be given by Eq. (46). The range of the exponential map given by coincides with the set R, consisting of points in SE (2) that can be reached with (globally minimizing) geodesics of P curve departing from e.
Proof Apply Theorems 1 and 3, where the analytic stationary solution curves of P curve break down iff = s max (p 0 ) in which case tangents to geodesics are vertical due to |z( )| = dθ dt (T ) = 1.
The exponential map of P curve coincides with the exponential map of P MEC [2,47] restricted to the strip ν(t (s)) ∈ [0, 2π] (in between the blue lines in Fig. 9), where we exclude the points (ν, c) = (0, 0) and (ν, c) = (2π, 0) from the strip (recall Remark 5.5) so that in the range we exclude the vertical line The exponential map of P MEC restricted to this strip is a homeomorphism (as follows by the results in [56]) thereby the exponential map of P curve is a homeomorphism as well. As a result (for formal proof see Appendix F) we have Theorem 6 Let D, R denote respectively the domain and range of the exponential map of P curve (recall Eqs. (46), (54)). Then with the subspace topology. 10 • Exp :D →R is a diffeomorphism. 10 As D and R are not open sets within the standard topologies on the embedding spaces T e (SE(2)) × R + and R 2 × S 1 . These subspace topologies do not coincide with the induced topology imposed by the embedding via the identity map, as such identity map is not continuous. However, with respect to the subspace topologies the set D, respectively R are open sets and the homeomorphism Exp : D → R is well-defined.
Finally, the boundary ∂R is given by These results can be observed in Fig. 11, which shows a well-posed, smooth, bijective relation between smooth regions in the phase portrait (i.e. D) and smooth regions in R ⊂ SE(2) and where the union of the blue and red surfaces form the cusp-surface adjacent to the line l. Subsequently, we provide some theorems on R and ∂R to get a better grip on the existence set of P curve , recall Eqs. (3) and (11).

The Elliptic Integral in the Exponential Map
In this section we will first express the single elliptic integral arising in the exponential map in Theorem 3 in a standard elliptic integral and then we provide bounds for this integral from which one can deduce bounds on the set R. For explicit bounds for the elliptic integral for the cases c < 1, where the sub-Riemannian geodesics are U-shaped, see Appendix H.

Observations and Theorems on R
In Theorem 3 we have derived the exponential map of P curve in explicit form. Before we derive some results on the range R of the exponential map we refer to Fig. 11 where we have depicted the set R using Theorem 3. In Fig. 11 we observe: 1. The range R of the exponential map is a connected, noncompact set and its piecewise smooth boundary coincides with the cusp-surface, Eq. (55). 2. The range of the exponential map produces a reasonable criterium (namely condition (3)) to connect two local orientations. Consider the set of reachable cones depicted in Fig. 14. 3. The range of the exponential map of P curve is contained in the half-space x f in ≥ 0 and |θ f in | = π can only be attained at x = 0 and y = 0 where geodesics arrive at a cusp. 4. The cone of reachable angles θ f in per position (x f in , where E(z, m) is given by Eq. (56). 6. The critical surface splits the range of the exponential map into four disjoint parts, cf. Fig. 11. These parts C 1 1 , C 0 1 , C + 2 and C − 2 directly relate to the splitting of the phase space, into the four parts C 1 1 , C 0 1 , C + 2 and C − 2 .
Let's underpin these observations with theorems.

Theorem 7
The range R of the Exponential map of P curve is contained within the half space x ≥ 0. In particular, its boundary ∂R (i.e. the cusp-surface) is contained within x ≥ 0.
Proof From Theorem 3 we deduce that One has (see Fig. 9) In the other cases in the phase portrait where the result is obvious. Via symmetry considerations one only needs to consider the case where z(s max ) = 1. Then we apply Lemma 5 (with a = −ż 0 and b = z 0 ) from which we deduce that In the remainder of this proof we will show that For analysis of R and ∂R and for (semi-)analytically solving of the boundary value problem the following identities (due to Theorem 3) come at hand.

Lemma 6
We have the following relation between the momentum at s = 0 p 0 = z 0 ω 1 + 1 − |z 0 | 2 ω 2 +ż 0 ω 3 and the end-condition g f in = (x f in , y f in , θ f in ): This yields a quadratic polynomial equation inż 0 : and whose solutions are expressed in z 0 viȧ  Fig. 13 for an illustration of such geodesics.

The Cones of Reachable Angles
We will provide a formal theorem that underpins our observations of the cone of reachable angles θ f in per end-position (x f in , y f in ), recall (57). Recall that θ endcusp (x f in , y f in ) denotes the final angle of the geodesic ending in (x f in , y f in , ·) with a cusp and where θ begincusp (x f in , y f in ) denotes the final angle of a geodesic ending in (x f in , y f in , ·) starting with a cusp. In case there exist two geodesics ending with a cusp at (x f in , y f in ) we order their end-angles by writing

Theorem 9
Let (x f in , y f in , θ f in ) ∈ R. If and then we have For a direct graphical validation of Theorem 9 see Fig. 11 (in particular the top view along θ ), where we note that the bound in (65) relates to the spatial projection of the curve that arises by taking the intersection of the blue and red surface on ∂R at θ = 0 (the thick black line in Fig. 11 at θ = 0). For more details on the proof see Appendix E.
As already mentioned in Sect. 3.1, it does not matter if one considers problem P curve on the projective line bundle R 2 P 1 or on R 2 S 1 ≡ SE(2). This is due to the following theorem.

Theorem 10
If then (x f in , y f in , θ f in + π) / ∈ R.

Solving the Boundary Value Problem
In order to explicitly solve the boundary-value problem for P curve for admissible boundary conditions (Eq. (3)) we can apply left-invariance (i.e. rotation and translation invariance) of the problem and consider the case g in = e = (0, 0, 0) and g f in ∈ R.
Recall from Eq. (20) that initial momentum p 0 is determined by z 0 andż 0 : Now solving the boundary value problem boils down to expressing (p 0 , ) directly into since when we achieve to do so we have and the globally minimizing curve of P curve is given by In fact, this means we must find the inverse of the exponential map Exp. The inverse of this exponential map exists due to Theorems 6, 4 and 1.
We invert the boundary value problem for a very large part analytically, yielding a novel very fast and highly accurate algorithm to solve the boundary value problem. In comparison to previous work on this topic [45], we have less parameters to solve (and moreover, our proposed optimization algorithm involves less parameters).
First of all we directly deduce from Theorem 3, Lemma 6 and Eq. (40) that where v, w, c are given by Now we have already expressed two of the three unknowns in the end condition = (z 0 ,ż 0 , g f in ) given by Eq. (67), z 0 =ż 0 (z 0 , g f in ) given by Eq. (64).
The remaining unknown variable z 0 ∈ [−1, 1] can be found via a simple numerical algorithm to find the unique root of a function F : I → R + , where I ⊂ [−1, 1] is a known and determined by g f in .
However, before we can formulate this formally there is a technical issue to be solved first, which is the choice of sign in Eq. (64).  (2) splits R into to parts and intersects the cusp surface ∂R at θ = π 2 and θ = − π 2 . The upper part requires positive sign in the formula forż 0 whereas the lower part requires a negative sign. Right: cross sections for x = 1 and x = 2 where V is contained within C 1 1 and C 0 1 and splits C + 2 and C − 2 , cf. Fig. 11 given by

Lemma 7 Let surface V ⊂ SE(2) be given by
(71) Proof The Exp is a (global) homeomorphism and its orbits s → Exp(p 0 , s) are analytic for each p 0 ∈ T * e (SE (2)). Thereby the sign cannot switch along orbits (unless D = 0, which only occurs at θ f in = ±π at ∂R). Furthermore, since Exp is a homeomorphism sign switches (in Eq. (64)) between neighboring orbits are not possible unless it happens across an orbit s → (z(s),ż(s)) withż 0 = 0. Now from the phase portrait it is clear that orbits in phase space s → (z(s),ż(s)) withż(s) > 0 and c > 1, i.e. orbits in C + 2 need a plus sign, whereas orbits in C − 2 need a minus sign in Eq. (64). The lineż 0 = 0 splits the phase portrait in two parts, and by the results in Theorem 6 this means that the surface V splits the set R into two parts. Now Exp maps C + 2 onto C + 2 and it maps C − 2 onto C − 2 , and C − 2 lies beneath V and C + 2 lies above V , from which the result follows.

Remark 7.1
The surface V is depicted in Fig. 15. Lemma 7 is depicted in Fig. 16, where we used Theorem 3 to compute for each point in 2] in phase space the sign of 2aż 0 + b at respectively s = 0, 1 2 s max (z 0 , 0), 3 4 s max (z 0 , 0) and s = s max (z 0 , 0). We see that the black points (where the sign is positive) lies above the orbits family of orbits with z 0 ∈ [−1, 1] andż 0 = 0.

Remark 7.2
The explicit parametrization for plane V is given by the union of the x-axis and the surface parame- The next theorem reduces the boundary value problem to finding the unique root of a single positive real-valued function.
Theorem 11 Let g f in ∈ R. The inverse of the exponential map in Definition 2 is given by given in Lemma 7 and with discriminant D(z 0 , g f in ) given by Eq. (63) and where z 0 denotes the unique zero F (z 0 ) = 0 of function F : I → R + defined on given by where · denotes the Euclidean norm on R 2 × S 1 .
Proof By Theorem 1 there is a unique stationary curve connecting e and g f in ∈ R. The exponential map of P curve is a homeomorphism by Theorem 6 and thereby the continuous function F has a unique zero, since andż 0 are already determined by z 0 and g f in via Theorem 3 and Lemma 7. Remark 7.3 Theorem 11 allows fast and accurate computations of sub-Riemannian geodesics, see Fig. 12 where the computed geodesics are instantly computed with an accuracy of relative L 2 -errors in the order of 10 −8 . Finally, we note that Theorem 6 implies that (our approach to) solving the boundary-value problem is well-posed (i.e. the solutions are both unique and stable).

Modeling Association Fields with Solutions of P curve
Contact geometry plays a major role in the functional architecture of the primary visual cortex (V1) and more precisely in its pinwheel structure, cf. [52]. In his paper [52] Petitot shows that the horizontal cortico-cortical connections of V1 implement the contact structure of a continuous fibration π : R × P → P with base space the space of the retina and P the projective line of orientations in the plane. This model was refined by Citti and Sarti [22], who formulated the model as a contact structure within SE(2) producing problem P curve given by Eq. (11). Petitot applied his model to the Field's, Hayes' and Hess' physical concept of an association field, to several models of visual hallucinations [32] and to a variational model of curved modal illusory contours [42,48,65].
In their paper, Field, Hayes and Hess [34] present physiological speculations concerning the implementation of the association field via horizontal connections. They have been confirmed by Jean Lorenceau et al. [43] via the method of apparent speed of fast sequences where the apparent velocity is overestimated when the successive elements are aligned in the direction of the motion path and underestimated when the motion is orthogonal to the orientation of the elements. They have also been confirmed by electrophysiological methods measuring the velocity of propagation of horizontal activation [37].
There exist several other interesting low-level vision models and psychophysical measurements that have produced similar fields of association and perceptual grouping [39,49,68], for an overview see [52, Chaps. 5.5, 5.6].

Three Models and Their Relation
Subsequently, we discuss three models of the association fields: horizontal exponential curves, Legendrian geodesics, and cuspless sub-Riemannian geodesics (which for many boundary conditions coincide with Petitot's circle bundle model, as we will explain below).
To this end Petitot [52] also proposed the "circle bundle model" which has the advantage that it is coordinate independent. Its energy integral can be expressed as 0 √ 1 + κ 2 ds, where s ∈ [0, ] denotes spatial arclength-parametrization. As long as the curve can be well-parameterized by x → (x, y(x), θ (x)) this model coincides 13 with sub-Riemannian geodesics.
For the explicit connections between each of the 3 mathematical models we refer to Appendix G.
On the one hand, a serious drawback arising in the co-circularity model for association fields is that the only the spatial part (x f in , y f in ) of the end-condition can be prescribed (the angular part is imposed by co-circularity), whereas with geodesics one can prescribe (x f in , y f in , θ f in ) (as long as the ending condition is contained within R). This drawback is clearly visible in Fig. 8, where the association 11 The dual basis in (SE(2)) 0 is equal to (dθ, dx, −θ dx + dy) and thereby the sub-Riemannian metric on (SE(2)) 0 does not include the |y (x)| 2 term. 12 The corresponding minimization problem (and induced sub-Riemannian distance) is left invariant in (SE(2)) 0 and not left-invariant in SE(2). 13 The preservation law and curvature in [52, Eq. (87)] does not fully match our results in [12,26,47,55,56], Appendix A and Theorem 2. On the other hand, the sub-Riemannian geodesic model has more difficulty describing the association field by Field and co-workers in the almost circular connections to the side (where the co-circularity model is reasonable). To this end we note that circles are not sub-Riemannian geodesics as the ODEz = ξz does not allow z to be constant.
This difficulty, however, can be tackled by variation of ξ in Problem P curve . Our algorithm explained in Sect. 5, combined with the scaling homothety described in Remark 1.2, is well-capable of reconstructing the almost circular field line cases as well. This can be observed in Fig. 17.

Variation of ξ and Association Field Modeling
See Fig. 17 to see the effect of ξ > 0 on the modeling of association fields. The larger ξ the shorter the spatial part of the paths, and the more bending we see in the vicinity of the end-points. The smaller the ξ the more circular the shape becomes at the sides of the association field model. Here we note that for these smaller values of ξ , the endpoints of the more straight association field lines become problematic. In Fig. 17 one can see that when choosing ξ too small the end-point of the most straight field line even lies outside the range R of the exponential map. This effect is due to the fact that the boundary ∂R of the range of the exponential map, depicted in Figs. 11 and 14, scales with ξ > 0 in spatial direction.
Varying of ξ 2 > 0 also takes into account a wellknown parameter in completion; namely the area of the completed figures (see e.g. [52]). This area equals A = (x f in − x in )(y f in − y in ). By Remark 1.1 we can as well set x in = y in = θ in = 0 and then as explained in Remark 1.2 solving P curve with ξ > 0 amounts to solving P curve with ξ = 1 with scaled end-conditions (x f in ξ, y f in ξ). In fact, such rescaling of end-conditions rescales the area as follows A → Aξ 2 .

A Conjecture and Its Motivation
The shape of the association field lines is well captured by the sub-Riemannian geodesics with ξ = 1, in comparison to e.g. the exponential curves as can be observed in part b) of Fig. 8. See also Fig. 17. On top of that, the field curves of the association field end with vertical tangent vectors, and these end-points are very close to cusp points in the sub-Riemannian geodesics modeling these field lines. This can be observed both in Fig. 4 and in Fig. 17, where the sub-Riemannian geodesics ending at the end-points of the association field is nearly vertical. We will underpin this observation also mathematically in Lemma 8 and Remark 8.1.
Apparently, both the shape of the association field lines and their ending is well-expressed by the sub-Riemannian geodesics model P curve , which was proposed by Citti and Sarti [22]. Therefore, following the general idea of Petitot's work [50] (in particular, his circle bundle model) and the results in this article on the existence set R this puts the following conjecture:

Conjecture 1
The criterium in our visual system to connect two local orientations, say g 0 = (x 0 , y 0 , θ 0 ) = (0, 0, 0) and g f in = (x f in , y f in , θ f in ) ∈ SE(2), could be modeled by checking whether g f in is within the range R of the exponential map.
Here we recall that from the results in [16] (summarized in Theorem 1) it follows that the set R consists precisely of those points in SE(2) that are connected to the origin by a unique global minimizer of P curve . This conjecture needs further investigation by psycho-physical and neurophysiological experiments. In any case, within the model P curve (relating to Petitot's circle bundle model [52] and the sub-Riemannian model by Citti and Sarti [22]) a curve is optimal if and only if it is stationary. Furthermore, the sub-Riemannian geodesics strongly deviate from horizontal exponential curves even if the end condition is chosen such that the co-circularity condition is satisfied (this can be observed in item c) of Fig. 8). This discrepancy between horizontal exponential curves and cusp-less sub-Riemannian geodesics in (SE(2), , G ξ ) is also intruiging from the differential geometrical viewpoint: see Theorem 12 in Appendix C.
In the remainder of this section we will mathematically underpin our observation that end-points of association fields are close to cusps.
Then forż 0 < 0 small we have Furthermore, under the conditions in Eq. (74), two of the following statements imply the remaining third one.
Proof If θ f in = θ( ) = π 2 then by Eq. (61) we have thaṫ z 0 = − 1 − |z( )| 2 , so that The rest follows by the fact that the second statement is equivalent to |z( )| = 1 and the formula for θ f in in Eq. (52).

Remark 8.1
The curves in the association field have θ f in = π 2 and relatively small initial curvature so that |ż 0 | 1 and therefore they end very close to cusps, i.e. ≈ s max .

Conclusion and Future Work
Under conditions (3) on the boundary conditions cuspless sub-Riemannian geodesics in SE(2), span{cos θ∂ x + sin θ∂ y , ∂ θ }, G ξ coincide with the lifts of global minimizers of P curve (i.e. curves optimizing 0 κ 2 + ξ 2 ds with free length and given boundary conditions).
As the derivation of these cuspless geodesics is much less trivial than it seems (many conflicting results have appeared in the imaging literature on this topic), we derived them via 3 different mathematical approaches producing the same results from different perspectives. There are two ways to reasonably parameterize such curves, via spatial arclength and sub-Riemannian arclength and in this article we explicitly relate these parameterizations. The phase portrait in momentum space induced by sub-Riemannian arclength parametrization corresponds to (a strip within) the phase portrait of the mathematical pendulum, whereas the phase portrait in momentum space induced by spatial arclength parametrization is a hyperbolic phase portrait associated to a linear ODE for normalized curvature z = κ/ κ 2 + ξ 2 . Using the latter approach we have analyzed and computed the existence set R for P curve (where every stationary curve is globally minimizing!). We have also solved the boundary value problem, where the numerics is reduced to finding the unique root of a continuous explicit real-valued function on a small subset of [−1, 1].
As such cuspless sub-Riemannian geodesics provide a suitable alternative to (involved and not necessarily optimal) elastica curves in computer vision. Moreover, they seem to provide a very adequate model for association fields and they are the solutions to Petitot's circle bundle model. They also relate to previous models for association fields based on horizontal exponential curves (i.e. "co-circularity") via the Cartan connection: Along horizontal exponential curves tangent vectors are parallel transported, whereas along sub-Riemannian geodesics momentum is parallel transported.
Our solutions, analysis and geometric control for the sub-Riemannain geodesics presented in this article form the venture point for data-dependent active contour models in SE(2) (in combination with contour-enhancement [1,14,22,26,29,35,36] and contour completion PDE's [4,8,30,48]) we are currently developing and applying in various applied imaging problems. Applications include extraction of the vascular tree in 2D-retinal imaging [10] and fiber-tracking in diffusion weighted magnetic resonance imaging [23,62] (where we use sub-Riemannian geodesics in SE(3) solving the 3D-version of P curve ). In these applications one replaces the constant measure on SE(2) in P curve by a data-dependent measureC : SE(2) → [1, ∞) in P curve , producing external force terms in the Euler-Lagrange equations that pull the geodesics towards the data.
Finally, future work will include comparison of numerical algorithms for P MEC and P curve .
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.
Remark Lemma 4 expresses the integral forỹ in a singe elliptic function.

B.2 The Geodesics for P MEC (and P curve ) Parameterized by Sub-Riemannian Arclength
Here we distinguish between different cases in the phase portrait of the mathematical pendulum, recall Fig. 9.

(82)
Remark 11.1 The geodesics of P MEC are defined for t ∈ R and every pair of points in SE(2) can be connected by a smooth geodesic.

Appendix C: The Cartan Connection and Its Relation to Geodesics and Exponential Curves
In this section we show that last three equations of the Pfaffian system Eq. (28) can be summarized in a single simple formula ∇γ p = 0, where ∇ denotes the Cartan connection on the co-tangent bundle T * (SE (2)) of SE (2).

])
Definition 3 The Cartan connection ∇ on the tangent bundle (SE (2), T (SE(2))) is given by the covariant derivatives withȧ k (t) =γ i (t)(A i | γ (t) a k ), for all tangent vectors X| γ (t) =γ i (t)A i | γ (t) along a curve t → γ (t) ∈ SE(2) and all sections μ(γ (t)) = 3 k=1 a k (γ (t))A k (γ (t)). The Christoffel symbols in (88) coincide with the structure constants of the Lie-algebra Theorem 12 Exponential curves are auto-parallel with respect to the Cartan connection ∇. Horizontal exponential curves are auto-parallel with respect to Cartan connection D. In fact, where Exp : T e (SE(2)) → SE(2) denotes the exponential map from Lie algebra to Lie group. Along an exponential  (2)) by means of the infinitesimal generator of the right-regular representation acting on C 1 (SE(2), R) via Via the identity e tdR(c i A i ) = R e 3 i=1 t(c i A i ) and application of the method of characteristics to linear, left-invariant convection systems on C 1 (SE (2) Thereby, for any Now, as the Christoffel symbols are anti-symmetric (see Eq. (89)) we have which implies Eq. (90). With respect to Eq. (91) we note that by duality Now Eq. (94) and Eq. (31) imply Eq. (91). z(l) = −z(0) = ±1 and c > 1) or at |θ f in | = π (cases z(l) = z(0) = ±1 and c > 1). This directly follows by the formula for θ f in in Eq. (52). Next we show that if a geodesic starts with a cusp it will have x f in < 2. This follows by the formula for x f in in Eq. (52) where the second term vanished if |z 0 | = 1 and where the first term is less than 2.
we again have More precisely, a geodesic starting with a cusp will satisfy Eq. (65), since the maximum value for x f in under the condition |z 0 | = 1 is obtained at |z( )| = 1 = |z 0 | and in these cases we find |ż 0 | =  (2) and consider the corresponding trajectory in phase space and assume that the trajectory has c < 1 andż 0 < −z 0 , then sign(z(s)) = sign(κ(s)) = sign(θ(s)) = −1 is for all s ∈ [0, ] and minimum value for |θ(s)| is obtained s = . We can extend this geodesic to [s min , s max ] with |z(s min )| = 1 and s min ≤ 0. Then the extended geodesic starts and ends in a cusp and we have θ f in ∈ [θ endcusp (x f in , y f in ), θ begincusp (x f in , y f in )] ⊂ R − due to the monotonic decrease of s → θ(s). This can be observed in Fig. 11 where all extended orbits in C 0 1 go down in theta direction converging towards the θ -minimum at the blue surface at s ↑ s max .
For geodesics with c > 1 the curvature κ(s) =θ(s) can switch sign only once at s B where z(s B ) = 0. By the assumption y f in ≤ 0 we can restrict ourselves to the nontrivial case z 0 < 0,ż 0 > 0 and c > 1, i.e. ((z 0 ,ż 0 ) ∈ C + 2 ) where we assume z( ) > 0 in order to consider the nonstraightforward case where the curvature switches sign. Again we extend the geodesic to [s min , s max ] (starting from a cusp ending at a cusp). The extended geodesic's angle initially decreases initially on [0, s B ], but then by reflection symmetry it will gain more than the initial decrease during s ∈ [s B , s max ]. Thereby, if a geodesics is within C + 2 its extension will stay in C + 2 , see Fig. 11 (as it cannot cross the critical surface) and converge (upwards in θ ) towards the blue surface (where geodesics end with a cusp) and we have θ f in ∈ [θ endcusp (x f in , y f in ), θ begincusp (x f in , y f in )]. Conclusion, if the bound (65) holds (x f in , y f in ) ∈ R 2 can be connected both with a geodesic starting with a cusp (at (0, 0) and a geodesic ending in a cusp at (x f in , y f in ), then the end-angles of these curves determine the cone of reachable orientations. In case Eq. (65) does not hold, there are no geodesics starting from a cusp that reach (x f in , y f in ) and in these cases the cone of reachable orientations is just bounded by the end-angles of two geodesics ending in a cusp (the two blue surfaces in Fig. 11).

Appendix F: Proof of Theorem 6
Within this proof we will rely on previous results by Sachkov in [47,55,56]. As these works rely on sub-Riemannian arclength parametrization and the pendulum phase portrait (Fig. 9), we will do the same in this proof. To this end (via Lemma 1) we re-express domain and range of the exponential map Exp of P curve in this parametrization where t cusp (p 0 ) = t (s max (z 0 ,ż 0 )) denotes sub-Riemannian arclength until the first cusp-time formally defined by t cusp (p 0 ) = inf t > 0 |ẋ(t) =ẏ(t) = 0 .
The function p 0 → t cusp (p 0 ) is continuous on C and it is a uniform lower bound [55] for the continuous function p 0 → t cut (p 0 ) which assigns to each initial momentum the corresponding time where global optimality of P MEC is lost along the corresponding stationary curve. Now that the preliminaries are done, let us start with the proof.
We will first show the mapping Exp :D →R is a diffeomorphism. Consider to this end the setŇ = {(p 0 , t) ∈ C × R + | t < t cut (p 0 )}. Since the function t cut : C → (0, ∞] is continuous, the setŇ is open. Let EXP denote the exponential map of P MEC . It follows from Th. 3.1 [55] that the mapping EXP|Ň is injective. Moreover, it was shown in Th. 2.5 [56] that this mapping is a local diffeomorphism. Thus the restriction EXP|Ň is a global diffeomorphism. Now since t cusp (p 0 ) ≤ t cut (p 0 ), we haveD ⊂Ň . So Exp|D ≡ EXP|D is a global diffeomorphism as well.
Regarding the final remark: Consider a sequence (g n ) n∈N in R that converges to some g ∈ SE(2) \R. We must show g = (x, y, θ) ∈ l or g is the end-point of a geodesic starting with a cups or g is the end-point of a geodesic ending at a cusp.
By the homeomorphic and diffeomorphic properties of Exp this sequence relates to a sequence in the phase portrait converging to a point (ν, c, t) ∈ ∂D. Let us consider all possible cases: • If (ν, c) = (0, 0) or (ν, c) = (2π, 0) we have g ∈ l.
• If (ν, c) = (0, c) with c > 0 or if (ν, c) = (2π, c) with c < 0 (i.e. the cases corresponding to the red lines in the phase portrait in Fig. 11b), this means g is the end-point of a geodesic starting with a cusp. • If (ν, c) = (0, c) with c > 0 or if (ν, c) = (2π, c) with c < 0 (i.e. the cases corresponding to the blue lines in the phase portrait in Fig. 11b), this means g is the end-point of a geodesic ending with a cusp.

Appendix G: Definition of Cusps, Geodesics and Association Field (Models)
We use the following definition of a geodesic.
Definition 4 A geodesic of P curve (respectively P MEC ) is a stationary curve γ of the corresponding geometric control problem formulated in Sect. 1.1. Geodesics of P curve are called cuspless sub-Riemannian geodesics.
Remark It can be shown that such a geodesic γ has the property that for every sufficiently small interval I = (t 1 , t 2 ) in the domain of such a curve the restriction γ | I is a minimizer between γ (t 1 ) and γ (t 2 ).
Smooth sub-Riemannian geodesics in P MEC may exhibit cusps when projected to the spatial plane, as we recall from Fig. 2. Roughly speaking, cusps are singular points in which spatial velocity changes its sign.
Remark Formally, in the minimizers of P curve , cusps do not occur, as the solutions break down at cusps. However, when |z(0)| = 1, where z(s) denotes normalized curvature Eq. (5) we say that a geodesic in P curve departs from a cusp. If |z( )| = 1, i.e. = s max (z 0 ,ż 0 ) given by Eq.(78), we say the geodesic of P curve ends in a cusp.

G.1 Association Field
The term association field comes from modeling contour integration in the human visual cortex by psycho-physical experiments [34,52]. The general idea of an association field is to provide an a priori link between relative positions and orientations within the sensorium of cortical columns in the primary visual cortex, Fig. 5. Intuitively, the tangents to the field lines of the association field provide expected local orientations, given that an local orientation is observed at the center of the association field.
Field and his co-workers psychophysical measurements have resulted in the association field depicted in item a) of Fig. 8. They relied on the hypothesis that the visual system can solve the continuity problem separately at each scale. 14 G.2 Models of the Association Field Within this article, we consider the cuspless sub-Riemannian geodesic model P curve , cf. [11,22,40], which is a natural extension of Petitot's circle bundle model [52, Chap. 6.6.5].
These other models relate to the cuspless sub-Riemannian geodesics as follows: • The Legendrian geodesics follow from the cuspless sub-Riemannian geodesic model by contracting (e.g. [25]) the sub-Riemannian manifold on (SE (2), , G ξ ), Eq. (17), towards its nilpotent approximation, cf. Eq. (73). • The horizontal exponential curves keeps the control variable in P curve constant and they can be considered as a rough local approximation of sub-Riemannian geodesics, see item (c) in Fig. 8. Finally, we recall Theorem 12.

Appendix H
We provide some estimates for the single elliptic integral appearing in the spatial arclength parametrization of cuspless sub-Riemannian geodesics, recall Theorem 3 and Lemma 4.

Lemma 9
We have the following lower and upper bound for the elliptic integral at s = s max , for the case z 0ż0 < 1 and c < 1: 1 − |z(s)| 2 ds ≤ s max , 14 Within the association field model P curve scaling of the endconditions amounts to scaling of ξ . are contained in the epigraph. Their slope equals −z 0ż0 √ 1−|z 0 | 2 . The rest follows by Fig. 18.