A PDE Approach to Data-driven Sub-Riemannian Geodesics in SE(2)

We present a new flexible wavefront propagation algorithm for the boundary value problem for sub-Riemannian (SR) geodesics in the roto-translation group $SE(2) = \mathbb{R}^2 \rtimes S^1$ with a metric tensor depending on a smooth external cost $\mathcal{C}:SE(2) \to [\delta,1]$, $\delta>0$, computed from image data. The method consists of a first step where a SR-distance map is computed as a viscosity solution of a Hamilton-Jacobi-Bellman (HJB) system derived via Pontryagin's Maximum Principle (PMP). Subsequent backward integration, again relying on PMP, gives the SR-geodesics. For $\mathcal{C}=1$ we show that our method produces the global minimizers. Comparison with exact solutions shows a remarkable accuracy of the SR-spheres and the SR-geodesics. We present numerical computations of Maxwell points and cusp points, which we again verify for the uniform cost case $\mathcal{C}=1$. Regarding image analysis applications, tracking of elongated structures in retinal and synthetic images show that our line tracking generically deals with crossings. We show the benefits of including the sub-Riemannian geometry.

1. Introduction. In computer vision, it is common to extract salient curves in images via minimal paths or geodesics minimizing a length functional. The minimizing geodesic is defined as the curve that minimizes the length functional, which is typically weighted by a cost function with high values on image locations with high curve saliency. To compute such datadriven geodesics many authors use a two step approach in which firstly a geodesic distance map to a source is computed and then steepest descent on the map gives the geodesics. In a PDE framework, the geodesic map is obtained via wavefront propagation as the viscosity solution of a Hamilton-Jacobi-Bellman equation (the Eikonal equation). For a review of this approach and applications see [37,27,31].
Another set of geodesic methods, partially inspired by the psychology of vision was developed in [15,30]. Here, the roto-translation group SE(2) = R 2 ⋊ S 1 endowed with a sub-Riemannian (SR) metric models the functional architecture of the primary visual cortex and geodesics are used for completion of occluded contours. A stable wavelet-like approach to lift 2D-images to functions on SE(2) was proposed in [16]. Within the SE(2) framework, images and curves are lifted to the 3D space R 2 ⋊ S 1 of coupled positions and orientations in which intersecting curves are disentangled. The SR-structure applies a restriction to so-called horizontal curves which are the curves naturally lifted from the plane (see Fig. 1A). For explicit formulas of SR-geodesics and optimal synthesis see [36]. SR-geodesics in SE (2) were also studied in [8,10,18,23,25]. Here, we propose a new wavefront propagation-based method for finding SR-geodesics within SE(2) with a metric tensor depending on a smooth external cost C : SE(2) → [δ, 1], δ > 0 fixed. Our solution is based on a Hamilton-Jacobi-Bellman 1 Figure 1. A: Every point in the planar curve γ2D(t) = (x(t), y(t)) is lifted to a point g = γ(t) = (x(t), y(t), θ(t)) ∈ SE(2) on an horizontal curve (solid line) by considering the direction of the tangent vectoṙ γ2D(t) of the planar curve as the third coordinate. Then, tangent vectorsγ(t) ∈ span{ A1| γ(t) , A2| γ(t) } = ∆| γ(t) , Eq. (2.1). B: In the lifted domain SE(2) crossing structures are disentangled. C: The SR-geodesic (green) better follows the curvilinear structure along the gap than the Riemannian geodesic (red).
(HJB) equation in SE(2) with a SR metric including the cost. It is of interest to interpret the viscosity solution of the corresponding HJB equation as a sub-Riemannian distance [38]. Using Pontryagin's Maximum Principle (PMP), we derive the HJB-system with an Eikonal equation providing the propagation of geodesically equidistant surfaces departing from the origin. We prove this in Thm. 3.1, and we show that SR-geodesics are computed by backtracking via PMP. In Thm. 3.2, we consider the uniform cost case (i.e. C = 1) and we show that the surfaces coincide with the SR-spheres, i.e. the surfaces from which every tracked curve is globally optimal. This uniform cost case has been deeply studied in [36] relying on explicit ODE-integration in PMP. In this article, we will rely on a PDE-approach and viscosity solutions of HJB-equations, allowing us to extend to the general case where C is a smooth cost uniformly bounded from below and above. We will often use the results in [36] as a golden standard to verify optimality properties of the viscosity solutions and accuracy of the involved numerics of our PDE-approach. We find a remarkable accuracy compared to exact solutions, 1st Maxwell sets (i.e. the location where for the first time two distinct geodesics with equal length meet), and to the cusp surface [18,10].
Potential towards applications of the method with non-uniform cost is demonstrated by performing vessel tracking in retinal images. Here the cost function is computed by lifting the images via oriented wavelets, as is explained in Section 4.1. Similar ideas of computing geodesics via wavefront propagation in the extended image domain of positions and orientations, and/or scales, have been proposed in [28,21,7]. In addition to these interesting works we propose to rely on a SR geometry. Let us illustrate some key features of our method. In Fig. 1B one can see how disentanglement of intersecting structures, due to their difference in orientations, allows to automatically deal with crossings (a similar result can be obtained with the algorithm in [28]). The additional benefit of using a SR geometry is shown in Fig. 1C where the SR-geodesic better follows the curvilinear structure along the gap.
1.1. Structure of the Article. The article is structured as follows. First, in Section 2, we give the mathematical formulation of the curve optimization problem that we aim to solve in this paper. We also provide embedding into previously studied geometrical control problems formulated on SE(2) and/or R 2 . In Section 3 we describe our PDE approach that provides the SR geodesic distance map as the viscosity solution of a boundary value problem (BVP) involving a sub-Riemannian Eikonal equation. Furthermore, in Theorem 3.1, we show that sub-Riemannian geodesics are obtained from this distance map by back-tracing (imposed by the PMP computations presented in Appendix A). In Theorem 3.2 we show that for the uniform cost case (i.e. C = 1) such back-tracing will never pass a Maxwell-points or conjugate points, and thereby for C = 1 our approach provides only the globally optimal solutions.
In Section 4 we describe an iterative procedure on how to solve the BVP by solving a sequence of initial value problems (IVP) for the corresponding HJB-equation. Before involvement of numerics, we express the exact solutions in concatenated morphological convolutions (erosions) and time-shifts in Appendix E. Here we rely on morphological scale space PDE's [12,2,17], and we show that solutions of the iterative procedure converges towards the sub-Riemannian distance map. Then in Section 4.1 we construct the external cost C, based on a lifting of the original image to an orientation score [16].
In Section 5, we describe a numerical PDE-implementation of our method by using leftinvariant finite differences [20] in combination with an upwind-scheme [33]. Finally, in Section 6, we present a verification of the proposed method by comparison with the exact solutions. We also provide simple numerical approaches (extendable to the non-uniform cost case) to compute 1st Maxwell points, conjugate points and cusp points [18], which we verify for the uniform cost case with results in [36]. Finally, in Subsection 6.2, we show application of the method to vessel tracking in optical images of the retina. We discuss the two main parameters that are involved: the balance between external and internal cost, and the balance between spatial and angular motion. First feasibility studies are presented on patches, and we discuss on how to proceed towards automated retinal vessel tree segmentation.
2. Problem Formulation. The roto-translation group SE(2) carries group product: where R θ is a counter-clockwise planar rotation over angle θ. This group can be naturally identified with the coupled space of positions and orientations R 2 ⋊ S 1 , by identifying R θ ↔ θ while imposing 2π-periodicity on θ. Then for each g ∈ SE(2) we have the left multiplication L g h = gh. Via the push-forward (L g ) * of the left-multiplication we get the left-invariant vector fields {A 1 , A 2 , A 3 } from the Lie-algebra basis {A 1 , A 2 , A 3 } = { ∂ x | e , ∂ θ | e , ∂ y | e } at the unity e = (0, 0, 0): Then all tangentsγ(t) ∈ T γ(t) (SE(2)) along smooth curves t → γ(t) = (x(t), y(t), θ(t)) ∈ SE(2) can be expanded asγ(t) = 3 k=1 u k (t) A k | γ(t) , where the contravariant components u k (t) of the tangents (velocities) can be considered as the control variables. Not all curves t → γ(t) in SE(2) are naturally lifted from the plane in the sense that θ(t) = arg(ẋ(t)+iẏ(t)). This holds for so-called horizontal curves which have u 3 = 0 and thusγ(t) The allowed (horizontal) directions in tangent bundle T (SE(2)) (see Fig. 1A) form a so-called distribution ∆ := span{A 1 , A 2 }.
We study the problem of finding SR minimizers, i.e. for given boundary conditions γ(0) = e, γ(T ) = g, we aim to find the horizontal curve γ(t) (havingγ ∈ ∆) that minimizes the total SR length If t is the SR arclength parameter, our default parameter, then G C | γ(t) (γ(t),γ(t)) = 1 and l = T . Then, SR minimizers γ are solutions to the optimal control problem (with free T > 0): In the naming of this geometric control problem we adhere to terminology in previous work [10,18]. Stationary curves of the problem (2.4) are found via PMP [1]. Existence of minimizers follows from Chow-Rashevsky and Filippov's theorem [1], and because of absence of abnormal trajectories (due to the 2-bracket generating distribution ∆) they are smooth.
Remark 2. The Cauchy-Schwarz inequality implies that the minimization problem for the SR length functional l is equivalent (see e.g. [26]) to the minimization problem for the action functional with fixed T : 2.1. Embedding into Geometric Control Theory. The problem (2.4) actually comes from a mechanical problem in geometric control, where a so-called Reeds-Shepp car [32] can proceed both forward and backward in the path-optimization. As pointed out in [9] such a problem, for certain end conditions, cannot be considered as a curve optimization problem on the plane. In fact in [17] the set R ⊂ SE(2) of allowable end-conditions g is explicitly determined for C = 1. The boundary of the set R is formed by end-points of geodesics whose spatial projections either depart or end in a cusp [17,Fig.12,Thm.9]. For g = (x, θ) ∈ R ⊂ SE(2) and C = 1 the following problem is well-posed where L denotes spatial length and κ curvature of smooth curve γ ∈ C ∞ ([0, L], R 2 ). For C = 1 the set of end-conditions (x, θ) for which the above problem is well-posed varies with C.
In this article we take P C mec (SE(2)) (and not (2.6)) as a venture point, as for the detection of elongated structures we must look both forward and backward. From the image analysis point of view, our primary interest is in the following refined problem of P C mec (SE(2)): (2.7) P C contour (SE(2)) : with curve γ : [0, T ] → SE(2), with controls: (u 1 (t), u 2 (t)) ∈ R 2 , and u 1 (t) does not change sign, where R C is the set of all g ∈ SE(2) such that the minimizing SR-geodesic(s) γ(·) = (x(·), θ(·)) do not exhibit a cusp in their spatial projections x(·). Remark 3.For existence results of the minimizers, a priori L 1 -conditions have to be imposed on the controls u 1 , u 2 . However, in retrospect, after application of a generalized version of PMP [39], [10, Thm 5.2 and App.B], the controls of the minimizers turn out to be smooth and bounded. Regarding extension to C = 1 (recall (2.2)) we note that δ ≤ C ≤ 1.
Remark 4.For C = 1 one has and with R ⊂ {(x, y, θ) ∈ SE(2) | x ≥ 0} described explicitly in [18], and partly depicted in item C of Figure 2. In Section 6.1 we shall provide a very simple numerical tool to compute the surface in SE(2) where cusps appear also for C = 1. This surface is a boundary of a volume in SE(2) that contains the set R C . Summarizing, in P C mec (SE(2)) the optimal control t → u 1 (t) can switch sign. • If g is chosen such that the optimal control u 1 ≥ 0 then the lift of problem P C curve (R 2 ) coincides with P C mec (SE(2)) and also with P C contour (SE(2)). Figure 2. 3 Sub-Riemannian geodesics for uniform cost C = 1 to reveal the differences of the 3 geometric control problems (2.4), (2.6) and (2.7), for β = 1. A: plots of spatial projections of sub-Riemannian geodesics in R 2 , B: SR-geodesics in SE(2), C: A part of R (the set of end-conditions for which Pcurve(R 2 ) is welldefined [18]) depicted as reachable cones around the origin. End condition g1 = (x1, y1, θ1) ∈ R yields the minimizing curve of Pcontour(SE(2)), Pmec(SE(2)) and Pcurve(R 2 ). End condition g2 = (−x1, y1, −θ1) yields the minimizing curve in Pcontour(SE(2)), Pmec(SE(2)), and it is invalid for Pcurve(R 2 ). End condition g3 is invalid for both Pcurve(R 2 ) and Pcontour(SE(2)), as it induces an internal cusp. For C = 1 the set of allowable end conditions for Pcontour(SE(2)) equals R C=1 = R ∪ Q. For C = 1 this set R C differs, and can be computed.

A B C
• If g is chosen such that the optimal control u 1 ≤ 0 then problem P C mec (SE(2)) and problem P C contour (SE(2)) coincide. • If g is chosen such that the optimal control u 1 (t) switches sign at some internal time t ∈ (0, T ), then g ∈ SE(2) \ R C and the spatial projection of the corresponding minimizing SR-geodesic(s) has an internal cusp, which we consider not desirable in our applications of interest. See Figure 2, where each of the 3 above cases is illustrated. Now let us return to our main objective, that is to derive solutions of (2.4) via a PDEwavefront propagation method for data-driven SR-geodesics. The word 'data-driven' refers to the fact that our PDE-approach is also applicable to C = 1, a property that does not apply to previous methodology following ODE-approaches [36,10,23,18,30] on SR-geodesics in vision/image analysis.
3. Solutions via Data-driven Wavefront Propagation. The following theorem summarizes our method for the computation of data-driven sub-Riemannian geodesics in SE (2). The idea is illustrated in Fig. 3.
Theorem 3.1. Let W (g) be a solution of the following boundary value problem (BVP) with Eikonal-equation Then the iso-contours are geodesically equidistant with speed dt dt = C(γ(t)) β 2 |u 1 (t)| 2 + |u 2 (t)| 2 = 1 and they provide a specific part of the SR-wavefronts departing from e = (0, 0, 0). A SR-geodesic departing from g ∈ SE(2) is found by backward integration Proof The next theorem provides our main theoretical result. Recall that Maxwell points are SE(2) points where two distinct geodesics with the same length meet. The 1st Maxwell set corresponds to the set of Maxwell-points where the distinct geodesics meet for the first time.
In the subsequent theorem we will consider a specific solution to (3.1), namely the viscosity solution as defined in Definition C.3 in Appendix C.
Theorem 3.2. Let C = 1. Let W (g) be the viscosity solution of the BVP (3.1). Then S t , Eq. (3.2), equals the SR-sphere of radius t. Backward integration via (3.3) provides globally optimal geodesics reaching e at t = W (g) = d(g, e) := = t} are non-smooth at the 1st Maxwell set M, cf. [36], contained in   [25] giving rise to interior folds (corresponding to multiple valued non-viscosity solutions of the HJB-equation). The Maxwell set M consists precisely of the dashed line on x cos θ 2 + y sin θ 2 = 0 and the red circles at |θ| = π. The dots are 2 (of the 4) conjugate points on St which are limits of 1st Maxwell points (but not Maxwell points themselves). In B we see the astroidal structure of the conjugate locus [35,13]. In A we see that the unique viscosity solutions stop at the 1st Maxwell set. Comparison of A and B shows the global optimality and accuracy of our method at A.
. From left to right: A: projection of γ1 and γ2 on the plane (x, y), B: 2D-slices (x = x * , y = y * ) of level sets of W (g) with distinguished value W (g) = t (again in orange). On top we plotted, the Maxwell point, the intersection of surface x cos θ 2 + y sin θ 2 = 0 (in purple, this set contains a part of the 1st Maxwell set) with the 2D-slices. C: The SR-sphere St in SE(2), D: section around g * revealing the upward kink due to the viscosity solution. From this kink we see that the tracking (3.3) does not cross a 1st Maxwell point as indicated in red, yielding global optimality in Thm. 3.2.
The Hamiltonian H f ree for the free time problem (2.4) minimizing l equals For details see Appendix A and B. Eq. (3.1) can be written as H f ree (g, p) = 0 with momentum 1  4. An Iterative IVP-procedure to Solve the SR-Eikonal BVP. To obtain an iterative implementation to obtain the viscosity solution of the SR-Eikonal BVP given by (3.1), we rely on viscosity solutions of initial value problems (IVP). In this approach we put a connection between morphological scale spaces [12,2], and morphological convolutions with morphological kernels, on the SR manifold (SE(2), ∆, G C ) and the SR Eikonal BVP.
In order to solve the sub-Riemannian Eikonal BVP (3.1) we resort to subsequent auxiliary IVP's on SE(2) for each r ∈ [r n , r n+1 ], with r n = nǫ at step n ∈ N ∪ {0}, ǫ > 0 fixed: n+1 (e, r n ) = 0 for n = 1, 2, . . ., and for n = 0, where δ M e is the morphological delta given by Let W ǫ n+1 denote the viscosity solution of (4.1) carrying the following support so in (4.1) at the n-th iteration (n ≥ 1) we use, for g = e, the end condition W ǫ n (g, r n ) of the n-th evolution as an initial condition W ǫ n+1 (g, r n ) of the (n + 1)-th evolution. Only for g = e we set initial condition W ǫ n+1 (e, r n ) = 0. Then we define the pointwise limit Finally, regarding the application of our optimality results, it is important that each IVPsolution W ǫ n+1 (g, r) is the unique viscosity solution of (4.1), as then via (4.3) the viscosity property for the viscosity solutions of the HJB-IVP problem naturally carries over to the viscosity property of the viscosity solutions of system (3.1). Thus we obtain W = W ∞ as the unique viscosity solution of the SR-Eikonal BVP.
Details on the limit (4.3), which takes place in the continuous setting before numeric discretization is applied, can be found in Appendix E. In Appendix E we provide solutions of (4.1) by a time shift in combination with a morphological convolution 2 with the corresponding morphological kernel, and show why the double limit is necessary. A quick intuitive explanation is given in Figure 6, where we see that for ǫ > 0 we obtain stair-casing (due to a discrete rounding of the distance/value function) and where in the limit ǫ ↓ 0 the solution W ∞ (g) = W (g) = d(g, e) is obtained.
As n grows the staircase grows, as ǫ → 0 the size of the steps in the staircase vanishes and we see W ∞ (g) = d(g, e) in the most right column.

This is important for representing viscosity solutions of left-invariant HJB-equations on SE(2)
by Lax-Oleinik [19] type of formulas (akin to the SE(3)-case [17]). Remark 8.The stair-casing limit depicted in Figure 6 is similar to the basic Eikonal BVP on R with solution d(x, 0) = |x|. On R the approach (4.1),(4.2) and (4.3) provides pointwise Remark 9.By general semigroup theory [2], one cannot impose both the initial condition and a boundary condition W ǫ (e, r) = 0 at the same time, which forced us to update the initial condition (at g = e) in our iteration scheme (4.1). The separate updating with value 0 for g = e in Eq.  C(x, y, θ) = 1 is a vessel detector where we use spatially isotropic Gaussian derivatives [20]. The vesseldetector gives responses only if there are convex variations orthogonal to the elongated- The lifting is done using real-valued anisotropic wavelets ψ: In this work we use the real part of so-called cake wavelets [16] to do the lifting. These wavelets have the property that they allow stable reconstruction and do not tamper data evidence before processing takes place in the SE (2) domain. Other type of 2D wavelets could be used as well. In related work by Péchaud et al. [28] the cost C was obtained via normalized cross correlation with a set of templates. In Eq. (4.5) two parameters, λ and p, are introduced. Parameter λ can be used to increase the contrast in the cost function. E.g., by choosing λ = 0 one creates a uniform cost function and by choosing λ > 0 data-adaptivity is included. Parameter p > 1 controls the steepness of the cost function, and in our experiments it is always set to p = 3.

Implementation.
To compute the SR geodesics with given boundary conditions we first construct the value function W ∞ in Eq. (3.1), implementing the iterations at Eq. (4.1), after which we obtain our geodesic γ via a gradient descent on W ∞ from g back to e, recall Thm. 3.1 (and Thm. 3.2). Throughout this section, we keep using the continuous notation g ∈ SE(2) although within all numerical procedures g is sampled on the following (2N + 1) × (2N + 1) × (2N θ ) equidistant grid: As a default we set N = 60, x max = 7, N θ = 32. The time-discretization grid is also chosen to be equidistant with time steps ∆r = ǫ.
On this grid we compute an iterative upwind scheme to obtain the viscosity solution W ǫ at iteration Eq. (4.1). Here we initialize W ǫ (·, 0) = δ M D e (·), with the discrete morphological delta, given by δ M D (g) = 0 if g = e and 1 if g = e, and iterate r)) for g = e W ǫ (e, r + ∆r) = 0, with free-time Hamiltonian (see Appendix A, Eq. (A.4)) given by until convergence. We set ∆r = ǫ in Eq. (4.1). In the numerical upwind scheme, the leftinvariant derivatives are calculated via where A + i and A − i denote respectively the forward and backward finite difference approximations of A i . Note that W ǫ in (5.2) is a first order finite difference approximation of W ǫ n+1 in (4.1) at time interval r ∈ [nǫ, (n + 1)ǫ] and we iterate until the subsequent L ∞ -norms differ less than 10 −6 . This upwind scheme is a straightforward extension of the scheme proposed in [33] for HJB-systems on R n . It produces sharp ridges at the 1st Maxwell set (cf. Fig. 4) as it is consistent at local maxima. For numerical accuracy and left-invariance we applied finite differences in the moving frame of left-invariant vector fields, using B-spline interpolation. This is favorable over finite differences in the fixed coordinate grid {x, y, θ}. For details on these kind of left-invariant finite differences, and comparisons to other finite difference implementations (e.g. in fixed coordinate grid) see [20].
In our implementation the origin e is treated separately as our initial condition is not differentiable. We apply the update W ǫ (e, r) = 0 for all r ≥ 0. We set step size ǫ = 0.1 min(s xy β, s θ ) with s xy and s θ step sizes in respectively the x-y-directions and θ-direction.
6. Experiments and Results. 6.1. Verification for the Uniform Cost Case. Throughout the paper we have illustrated the theory with figures obtained via our new wavefront propagation technique. In this subsection we go through the figures that support the accuracy of our method. As the problem (2.4) for C = 1 was solved [36,18], we use this as a golden standard for comparison. Figure 8. A: SR-geodesic example for the uniform cost case shows or PDE-discretizations (with 12 and 64 sampled orientations in red and green respectively) are accurate in comparison to analytic solutions in [36,18] (in black). B: The blue surface represents the cusp surface numerically computed via the proposed HJB-system (with C = 1) and subsequent calculation of the zero-crossings of A1W ∞ (x, y, θ). Indeed if a SR-geodesic (in green) passes this surface, it passes in θ-direction (with infinite curvature [10,18]), yielding a cusp on the spatial ground plane. The same blue surface is computed in [18, Fig. 11]. We even see the additional fold (top left passing the grey-plane) as some globally optimal SR-geodesics even exhibit 2 cusps.  Fig. 8A . Using the semi-analytic approach for solving the BVP in [18] an almost identical result is obtained. The curves computed with our method with angular step-sizes of 2π/12 and 2π/64 are shown in Fig. 8A in red and green respectively. Already at low resolution we observe accurate results. In Fig. 4 we compare one SR-sphere for T = 4 (Fig. 4A) found via our method with the SR-wavefront departing from e (Fig. 4B) computed by the method of characteristics [25]. We observe that our solution is non-smooth at the 1st Maxwell set M (3.5) and that the unique viscosity solution stops precisely there, confirming Theorem 3.2. Finally, the blue surface in Fig. 8B represents the cusp surface, i.e. the surface consisting of all cusp points. Cusps are points that can occur on geodesics when they are projected into the image plane (see Fig. 8B). This happens at points g c where the geodesic is tangent to ∂ θ | gc = A 2 | gc . Then, the cusp surface S cusp is easily computed as a zero-crossing:

B A
It is in agreement with the exact cusp surface S cusp analytically computed in [18, Fig. 11]. Remark 10.The simple key geometric idea behind (6.1) is that we have a cusp at time t if u 1 (t) = 1 C 2 (γ(t)) h 1 (t) = 1 C 2 (γ(t)) A 1 W ∞ (γ(t)) = 0 which directly follows from Appendix A.

Comparison and Computation of SR-spheres.
Numerical verification of the solution obtained by our PDE approach was also done by comparison with the exact SR-distance, that was computed by explicit formulas for SR-geodesics (given on p.386 in [25]) and the explicit formulas for cut time (that coincides with the first Maxwell time, given by (5.18)- (5.19) in [25]). The experiments were done in the following way: 1. Compute a set of points: of end points lying on the SR-sphere of fixed radius T by analytic formulas for the exponential map and 1st Maxwell time t M AX 1 [25]. The number of end points i max was chosen as i max = 72 T 2 . Here C is the cylinder in momentum space given by (2)) | H f ixed (e, p) = 1/2 , where we recall (3.6). The sampling points p i are taken by a uniform grid on the rectifying coordinates (ϕ, k) of the mathematical pendulum (the ODE that arises in the PMP procedure, cf. [25, ch:3.2]), both for the rotating pendulum case (p i ∈ C 2 , yielding S-curves) and the oscillating pendulum case (p i ∈ C 1 , yielding U -curves), where we note that C = C 1 ∪ C 2 . 2. Evaluate the distance function W ∞ (g i ) = W ∞ (x i , y i , θ i ) obtained by our numerical PDE-approach in section 5 for every point of the set EP (T ). We use 3rd order Hermite interpolation for W ∞ (x i , y i , θ i ) at g = g i ∈ EP (T ) in between the grid (5.1).

Compute the maximum absolute error max
|W (g i ) − T | and the maximum relative error max Remark 11.The exponential map Exp : C × R + → SE(2) provides the end-point (x(t), y(t), θ(t)) = γ(t) = Exp(p 0 , t) of the SR-geodesic γ, given SR-arclength t and initial momentum p 0 ∈ T * e (SE (2)). This exponential map integrates the PMP ODE's in Appendix A.
The absolute and relative errors of the SR-distance computations for each end points located on SR-spheres of fixed radii are presented in Figure 10. The red graph corresponds to a sampling of (N, N θ ) = (50, 64), recall Eq. (5.1), used in the SR-distance computation by our numerical PDE approach, and the blue graph corresponds to the finer sampling (N, N θ ) = (140, 128). We see that the maximum absolute error does not grow, and that the relative error decreases when increasing the radius of the SR-sphere. Increase of sampling rate improves the result. For the finer sampling case, neither the absolute errors nor the relative errors exceed 0.1.

Comparison and Computation of 1st Maxwell Set.
Considering forward and backward derivatives acting on the distance function W ∞ (x, y, θ), we can compute the 1st Maxwell set (recall eq. (3.5), see also Appendix D) as set of points where forward and backward left-invariant derivatives acting on distance function have different signs:  Here i = 1 corresponds to the local component of the 1st Maxwell set (i.e. the part around e), and i = 2 corresponds to the global component of the 1st Maxwell set. The local component consists of two connected components lying on the surface given by x cos θ 2 + y sin θ 2 = 0 (i.e. the purple surface in Figure 5), and the global component is a plane given by equation θ = π (for details, see [36]). In figure 11 we compare the local component of M num computed by our PDE approach with its exact counterpart M, presented in [36, thm 3.5]. It shows that M num is close to the exact M, except for the conjugate points at the boundary (where (6.2) is less accurate). Although note depicted here a similar picture was obtained for the global component, where M num indeed covers the plane θ = π. Summaring, this experiment verifies the correctness of the proposed method, but it also shows that the method allows to observe the behavior of the 1st Maxwell set. Eq. (6.2) allows us to numerically compute the Maxwell set for the data-driven cases C = 1 where exact solutions are out of reach.

Feasibility Study for Application in Retinal
Imaging. As a feasibility study for the application of our method in retinal images we tested the method on numerous image patches exhibiting both crossings, bifurcations, and low contrast, (Fig. 12, Fig. 13). For each seed point g 0 the value function g → W ∞ (g −1 0 g) was calculated according to the implementation Figure 11. Comparison of the 1st Maxwell set obtained by our numerical PDE-approach with the exact 1st Maxwell set [36]. Note that the local components of the 1st Maxwell set are part of the purple surface in For the construction of the cost function (see e.g. Fig. 7) we set p = 3, and the lifting was done using cake wavelets with angular resolution π/16. More precisely we used a cake-wavelet with standard parameters (N = 8, N θ = 32, s θ = π 8 , σ s = 20px, γ = 0.8), for details see cf. [4, ch:2]. The precise choice of anisotropic wavelet is not decisive for the algorithm (so other type of anisotropic wavelets and cost constructions could have been applied as well).
In all experiments we run with 4 settings for the two parameters (β, λ) determining the sub-Riemannian geodesics, we set β small = 0.05, β large = 0.1, λ small = 10, λ large = 100. The idea of these settings is to see the effect of the parameters, where we recall β controls global stiffness of the curves, and λ controls the influence of the external cost. We also include comparisons to a Riemannian wavefront propagation method on R 2 , and a Riemannian wavefront propagation method on SE(2). These comparisons clearly show the advantage of including the sub-Riemannian geometry in the problem. For results on two representative patches, see Figure 12. For results on 25 other patches see www.bmia.bmt.tue.nl/people/RDuits/Bekkersexp.zip.
Typically, the wavefront propagation tracking methods on (R 2 , G R 2 ) produces incorrect shortcuts at crossings and very non-smooth curves. The Riemannian wavefront propagation tracking method (with spatial isotropy) (SE (2), G C f ull ) often deals correctly with crossings, but typically suffers from incorrect jumps towards nearly parallel neighboring vessels. Also it yields non-smooth curves. This can be corrected for when including extreme anisotropy, see Remark 12 below. The Sub-Riemannian wavefront propagation method produces smooth curves that appropriately deals with crossings. For high contrast images with reliable cost C best results are obtained with low β and large λ. However, in low contrast images and/or patient data with severe abnormalities, low λ is preferable, as in these cases the cost function is less reliable see Figure 13. Remark 12. It is possible to construct a family of anisotropic Riemannian metric tensors, recall (3.7)): , which bridges the SR-metric G C of our method (obtained by ǫ → ∞) to the full Riemannian metric tensor G C f ull (obtained by ǫ → 1). For the values of β's considered here, Riemannian geodesics and smooth Riemannian spheres for highly anisotropic cases ǫ ≥ 10 approximate SR-geodesics and non-smooth SRspheres. In fact, with such extreme anisotropy in the Riemannian setting, the non-smooth ridges M in the SR spheres (see. e.g. the 1st Maxwell sets in Figure 4) are only little smoothed, and also the cusp-surface hardly changes. This observation allows to use the anisotropic fastmarching [24] as an alternative fast method for computing the solution of (3.1), instead of the iterative upwind finite difference approach in Section 4.
The experiments indicate β = 0.01 (small) in combination with λ = 100 (large) are preferable on our patches. This typically holds for good quality retinal images of healthy volunteers. In lesser quality retinal images of diabetic patients, however, the cost function is less reliable and here λ = 10 (small) can be preferable, see Figure 13.
However, it might not be optimal to set the β parameter globally, as we did in these experiments, as smaller vessels are often more tortuous and therefore require more flexibility, see e.g. [5, fig.7]. Furthermore, we do not include precise centerline extraction, which could e.g. be achieved by considering the vessel width as an extra feature (as in [7,28,21]).
In future work we will pursue on sub-Riemannian fast-marching methods for fully automated vascular tree extraction starting from an automatically detected optic nerve head via state-of-the art method [6] followed by piecewise sub-Riemannian geodesics (comparable to [14] in a Riemannian setting) in between boundary points detected by a SE(2)-morphological approach. First experiments show that such a SR fast-marching method leads to considerable decrease computation time, hardly reduces the accuracy of the method, and can be used to perform accurate and fast automatic full retinal tree segmentation. The advantage of such an approach over our previous work on automated vascular tree detection [4] is that each curve is a global minimizer of a formal geometric control curve optimization problem. However, the SR-fast marching and automatic detection of the complete vascular tree via piecewise sub-Riemannian geodesics is beyond the scope of this theoretically oriented paper.

Conclusion.
In this paper we propose a novel, flexible and accurate numerical method for computing solutions to the optimal control problem (2.4), i.e. finding SR-geodesics in SE(2) with non-uniform cost. The method consists of a wavefront propagation of geodesically equidistant surfaces computed via the viscosity solution of a HJB-system in (SE(2), ∆, G C ), Figure 12.
Data-adaptive sub-Riemannian geodesics obtained via our proposed tracking method (Thm. 3.2), with external cost (4.5), with p = 3, β equals β small = 0.01, β large = 0.1 and λ equals λ small = 10, λ large = 100. We applied tracking from 2 seed-points each with several end-points (to test the crossings/bifurcations). To distinguish between tracks from the two seed-points we plotted tracts in different lighting-intensity. We indicated the valid cases only if all trajectories are correctly dealt with. and subsequent backwards integration, which gives the optimal tracks. We used PMP to derive both the HJB-equation and the backtracking. We have shown global optimality for the uniform case (C = 1) and that our method generates geodesically equidistant surfaces. Compared to previous works regarding SR-geodesics in (SE(2), ∆, G 1 ) [36,18,23], we solve the boundary value problem without shooting techniques, using a computational method that Figure 13. Tractography (again for λ = λ small = 10, λ = λ large = 100 and β = 0.01, p = 3) in a patch of a challenging low-contrast retinal image of a diabetic patient. In case of low-contrast (and less reliable cost) it is better to keep λ small, in contrast to high contrast cases depicted in Figure 12. To distinguish between tracks from the two seed-points we plotted tracts in different lighting-intensity. We indicated the valid cases only if all trajectories are correctly dealt with.
always provides the optimal solution. Compared with wavefront propagation methods on the extended domain of positions and orientations in image analysis [28,29], we consider a SR metric instead of a Riemannian metric. Results in retinal vessel tracking are promising.
Fast, efficient implementation using ordered upwind schemes (such as the anisotropic Fast Marching method presented in [24]) is planned as future work as well as adaptation to other Lie groups such as SE(3), SO(3). Of particular interest in neuroimaging is application to high angular resolution diffusion imaging (HARDI) by considering the extension to SE(3) [17,29].
Acknowledgements. The research leading to these results has received funding from the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013) / ERC grant Lie Analysis, agr. nr. 335555. The authors gratefully acknowledge EU-Marie Curie project MANET ag. no.607643 for financial support.
with dual basis {ω i } for T * (SE(2)) defined by ω i , A j = δ i j . For a consistency check, we also apply the PMP-technique directly to Problem (2.4) with free terminal time T , where typically [1] the Hamiltonian vanishes. Then, using SR arclength parameter t, the control dependent Hamiltonian of PMP equals Optimization over all controls under SR arclength parametrization constraint C β 2 |u 1 | 2 + |u 2 | 2 = 1 produces via EL-optimization w.r.t. (u 1 , u 2 ) (via unit Lagrange multiplier) the (maximized) Hamiltonian: and by straightforward computations one can verify that both the horizontal part and the vertical part of PMP (but now applied to H f ree ) is exactly the same as (A.3) and (A.2).
Remark 13.The two approaches produce the same curves and equations, but their Hamiltonians are different. Nevertheless, we have H f ree = 0 ⇔ H f ixed = 1 2 .
Appendix B. Lemmas Applied in the Proof of Theorem 3.1.
In this section we consider preliminaries and lemmas needed for Thm. 3.1. Before we can make statements on SR-spheres we need to explain the notion of geodesically equidistant surfaces, and their connection to HJB-equations. In fact, propagation of geodesically equidistant surfaces in (SE(2), ∆, G C ) is described by a HJB-system on this SR-manifold.
Recall Remark 2. Also recall, that in Appendix A we have applied PMP to this problem yielding constant Hamiltonian Here P * h i ω i is a dual projection expressed in dual basis ω i given by Proof Substitute an arbitrary transversal minimizer γ(r) into V (·, r) and take the total derivative w.r.t. r: As a result we have ⇔ H(P * ∆ dV (γ(r), r)) = − ∂V ∂r (γ(r), r), Now every point g ∈ S r is part of a transversal minimizing curve γ(r) and the result follows. So the "⇒" is proven. Conversely, if the HJB-equation is satisfied it follows by the same computations (in reverse order) that L(γ(r),γ(r)) = d dr V (γ(r), r), which equals W ′ 0 (r). Remark 15.In PMP [1] (see also Appendix A) the controls are optimized to obtain the Hamiltonian H from the control dependent Hamiltonian H u . The first equivalence in (B.3) is due to the maximum condition of PMP. The second equivalence in (B.3), is by definition of the Hamiltonian, where by the convexity assumption of the Lagrangian the supremum is actually a maximum [19, ch:8].

Remark 16.The relation between the various Hamiltonians is
is the free-time Hamiltonian given by (3.8),

Remark 18.Combined
Cauchy-Dirichlet problems exist [38], but they are defined on (analytic) open and bounded domains. Thereby they cannot be applied to our set of interest SE(2) \ {e} as this would violate semigroup theory [2,19,34,12]. This is also clear in view of the Cramer transform [2], putting an isomorphism between HJB-and diffusion systems. symmetries (i.e. we include the constraint t ≤ min{t i max } where the minimum is taken over all positive Maxwell times along each trajectory), then we find M to be contained within the union of the following sets 3 : where [25, th:5.2] shows we must discard the first reflectional symmetry ǫ 1 as it does not produce Maxwell points. Now for generic geodesics (not passing the special conjugate points that are limit points of Maxwell points and not Maxwell points themselves) t cut = t 1 M AX , as proven in [35, th:3.3], where t 1 M AX > 0 denotes the first Maxwell time associated to the 8 discrete reflectional symmetries.
During the back-tracking the set M is never reached at internal times (only when starting at them, recall Remark 6), since they are "uphill" from all possible directions during dual steepest descent tracking (3.3), as we will shown next. As a result we have t ≤ t cut = t 1 M AX . Consider Fig. 5. At Maxwell points g * ∈ M due to the reflectional symmetries there exist two distinct directions in the 2D-horizontal part ∆ g * of the tangent space T g * (SE (2)) where the directional derivative is positive. If there would be a direction in the tangent space where the directional derivative is negative then there would be a direction in ∆ g * with zero directional derivative of W (·) at g * towards the interior of the sphere yielding contradiction. Here we note that due to the viscosity property of the HJB-solution, kinks at the Maxwell points are pointing upward (see Fig. 5 and Fig. 14) in the backward minimization tracking process [11, fig.30]. Furthermore, we note that SR spheres S t are continuous [36] and compact, as they are the preimage S t = d(·, e) ← ({t}) of compact set {t} under continuous mapping d(·, e). Continuity of d(·, e) implies the spheres are equal to the 2D-boundaries of the SR balls (w.r.t. the normal product topology on R 2 × S 1 ).
The algorithm also cannot pass conjugate points that are limits of 1st Maxwell points, but not Maxwell points themselves. See Fig. 4. Such points exist on the surface R 2 = 0 and are by definition within M \ M. Suppose the algorithm would pass such a point at a time t > 0 (e.g. there exist 4 such points on the sphere with radius 4, see Fig. 5) then due to the astroidal shape of the wavefront at such a point, cf. [35, Fig.11], there is a close neighboring tract that would pass a 1st Maxwell point which was already shown to be impossible (due to the upward kink-property of viscosity solutions).
Remark 20.The sub-Riemmannian spheres are non-smooth only at the 1st Maxwell set M. They are smooth at the conjugate points in M \ M (where the reflectional symmetry no longer produce two curves/fronts). In the other points on S t \ M the SR-spheres are locally smooth (by the Cauchy-Kovalevskaya theorem and the semigroup property of the HJB-equations). Figure 14. Overview of Maxwell points. Two Maxwell points on the purple surface x cos θ 2 + y sin θ 2 = 0 and two red Maxwell points on the surface |θ| = π (recall Figure 5). In all cases we see that local kinks in the viscosity solutions are upward, and the back-tracking algorithm can not pass these points. algebra, i.e. obtained by morphological convolution (erosions) with the morphological impulse response, cf. [12]. Now as the HJB-equations of such morphological scale spaces do not involve a global offset by 1 in the right-hand side of the PDE we need to combine such erosions with a time shift in order to take into account the global offset. It turns out that the combination of these techniques provide staircases with steps of size ǫ, so that we obtain the appropriate limit by taking the limit ǫ → 0 afterwards as is done in (4.3).
Morphological convolutions over the SE(2) group are obtained by replacing in linear leftinvariant convolutions (likewise the SE(3)-case studied in [17]), the usual (+, ·)-algebra by the (min, +)-algebra. Such erosions on SE(2) are given by Furthermore, to include the updating of the initial condition in (4.1) we define (E.2)W (g) := W (g) if g = e 0 if g = e.
where we use short notation for the sub-Riemannian derivative d SR V := 2 i=1 A i V ω i , recall (3.7) in Remark 5, and where HamiltonianH is given by (B.7). Now let us first consider the case n = 0. By the results in Appendix B the Hamiltonian system (4.2) provides geodesically equidistant wavefront propagation traveling with unit speed and departing directly from the unity element.As a result, we find V ǫ 1 (g, r new ) = k rnew (g) = 0 if d(g, e) ≤ r new ∞ else.
Corollary E.2. Let n ∈ N. Let ǫ > 0. The following identity holds (E.4) W ǫ n+1 (g, r n+1 ) = (k ǫ ⊖W ǫ n (·, r n ))(g) + ǫ Proof The first part follows by Lemma E.1 for r = r n+1 (i.e. r new = ǫ). The second part follows by induction. Recall from Lemma E.1 that W ǫ 1 (g, r) = k r (g) + r. Now application of (E.3) for n = 1 yields (E.5) W ǫ 2 (g, r 2 ) = (k ǫ ⊖W ǫ 1 (·, r 1 ))(g) + ǫ = ǫ + inf with B g,ǫ = {h ∈ SE(2) | d(g, h) ≤ ǫ}. This can intuitively be seen from the geometric meaning of an erosionW ǫ 1 → k ǫ ⊖W ǫ 1 where one drops cylinders from below on the graph of W ǫ 1 (·, r n ) and considering the new hull where cilinders get stuck. Eq. (E.5) can also be seen directly from the definition of k ǫ . Let us verify each case separately: • If d(g, e) > 2ǫ we have that the value must be infinite, since suppose it were finite then by the definition of the morphological kernel k ǫ we would need to have that d(g, e) ≤ d(g, h) + d(h, e) ≤ 2ǫ yielding contradiction. • If d(g, e) ≤ ǫ, then in the erosion-minimization we can take h = e and we obtain ǫ + 0.