Tracking of Lines in Spherical Images via Sub-Riemannian Geodesics in \({\text {SO(3)}}\)

Journal of Mathematical Imaging and Vision, Feb 2017

In order to detect salient lines in spherical images, we consider the problem of minimizing the functional \(\int \limits _0^l \mathfrak {C}(\gamma (s)) \sqrt{\xi ^2 + k_g^2(s)} \, \mathrm{d}s\) for a curve \(\gamma \) on a sphere with fixed boundary points and directions. The total length l is free, s denotes the spherical arclength, and \(k_g\) denotes the geodesic curvature of \(\gamma \). Here the smooth external cost \(\mathfrak {C}\ge \delta >0\) is obtained from spherical data. We lift this problem to the sub-Riemannian (SR) problem in Lie group \({\text {SO(3)}}\) and show that the spherical projection of certain SR geodesics provides a solution to our curve optimization problem. In fact, this holds only for the geodesics whose spherical projection does not exhibit a cusp. The problem is a spherical extension of a well-known contour perception model, where we extend the model by Boscain and Rossi to the general case \(\xi > 0, \mathfrak {C} \ne 1\). For \(\mathfrak {C}=1\), we derive SR geodesics and evaluate the first cusp time. We show that these curves have a simpler expression when they are parameterized by spherical arclength rather than by sub-Riemannian arclength. For case \(\mathfrak {C} \ne 1\) (data-driven SR geodesics), we solve via a SR Fast Marching method. Finally, we show an experiment of vessel tracking in a spherical image of the retina and study the effect of including the spherical geometry in analysis of vessels curvature.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

https://link.springer.com/content/pdf/10.1007%2Fs10851-017-0705-9.pdf

Tracking of Lines in Spherical Images via Sub-Riemannian Geodesics in \({\text {SO(3)}}\)

Tracking of Lines in Spherical Images via Sub-Riemannian Geodesics in SO(3) A. Mashtakov 0 1 2 3 4 R. Duits 0 1 2 3 4 Yu. Sachkov 0 1 2 3 4 E. J. Bekkers 0 1 2 3 4 I. Beschastnyi 0 1 2 3 4 R. Duits 0 1 2 3 4 0 MAMA, International School for Advanced Studies , Trieste , Italy 1 BMIA & CASA, Eindhoven University of Technology , Eindhoven , The Netherlands 2 RUDN University , Moscow , Russia 3 CASA & BMIA, Eindhoven University of Technology , Eindhoven , The Netherlands 4 CPRC, Program Systems Institute of RAS , Pereslavl-Zalessky , Russia In order to detect salient lines in spherical images, we consider the problem of minimizing the functional l ξ > 0, C = 1. For C = 1, we derive SR geodesics and evaluate the first cusp time. We show that these curves have a simpler expression when they are parameterized by spherical arclength rather than by sub-Riemannian arclength. For case C = 1 (data-driven SR geodesics), we solve via a SR Fast Marching method. Finally, we show an experiment of vessel tracking in a spherical image of the retina and study the effect of including the spherical geometry in analysis of vessels curvature. Sub-Riemannian geodesics; Geometric control; Spherical image; Lie group SO(3); Vessel tracking - 0 boundary points and directions. The total length l is free, s denotes the spherical arclength, and kg denotes the geodesic curvature of γ . Here the smooth external cost C ≥ δ > 0 is obtained from spherical data. We lift this problem to the sub-Riemannian (SR) problem in Lie group SO(3) and show that the spherical projection of certain SR geodesics provides a solution to our curve optimization problem. In fact, this holds only for the geodesics whose spherical projection does not exhibit a cusp. The problem is a spherical extension of a well-known contour perception model, where we extend the model by Boscain and Rossi to the general case B A. Mashtakov 1 Introduction In computer vision, it is common to extract salient curves in flat images via data-driven minimal paths or geodesics [1–5]. The minimizing geodesic is defined as the curve that minimizes the length functional, which is typically weighted by a cost function with high values at image locations with high curve saliency. Another set of geodesic methods, partially inspired by the psychology of vision, was developed in [13,14]. In these articles, sub-Riemannian (SR) geodesics in respectively the Heisenberg H(3) and the Euclidean motion group SE(2) are proposed as a model for contour perception and contour integration. The combination of such contour perception models with data-driven geodesic methods has been presented in [40]. There, a computational framework for tracking of lines via globally optimal data-driven sub-Riemannian geodesics on the Euclidean motion group SE(2) has been presented with comparisons to exact solutions [28]. Fig. 1 Photography of the retina. A part of the retina is projected onto the image plane. The camera coordinates are denoted by (X, Y ), and object coordinates are denoted by (x¯, y¯) In this work, we extend this framework for tracking of lines in spherical images (e.g., images of the retina, cf. Fig. 1). This requires a sub-Riemannian manifold structure in a different Lie group, namely the group SO(3) (consisting of 3D rotations) acting transitively on the 2-sphere S2. Here we study the problem Pcurve(S2) of finding a smooth curve n(·) on a unit sphere S2 that satisfies given boundary conditions (both positions and velocities) n(0) = n0, n(l) = n1, n (0) = n0, n (l) = n1 and minimizes the functional and minimizes the functional where kg (·) denotes the geodesic curvature of n(·), s denotes the spherical arclength, and total length l is free, see Fig. 2. Furthermore, we include a curve-stiffness parameter ξ > 0. In the optimization functional, we also include an external cost factor C : S2 → R+ adding for adaptation to a given spherical image data. We state this problem Pcurve(S2) more explicitly in Sect. 2.2. The problem Pcurve(S2) is a spherical analog of a wellknown problem Pcurve(R2) [30, 31] (cf. Fig. 2) of finding a smooth curve x(·) on a plane R2 that satisfies given boundary conditions x(0) = x0 = ( X0, Y0), x (0) = x0 = (cos Θ0, sin Θ0), x(l) = x1 = ( X1, Y1), x (l) = x1 = (cos Θ1, sin Θ1), Fig. 2 Left: problem Pcurve(S2): for given boundary conditions on a 2D sphere (both positions and velocities), we aim to find a curve minimizing the functional compromising length and geodesic curvature. In the optimization functional, we also include an external cost induced by spherical image data. Right: problem Pcurve(R2) [30,31]: for given boundary conditions on a 2D plane, to find a curve minimizing the compromise between length and curvature. The external cost factor is added for adaptation to flat image data (see [40]). where k(s) denotes the curvature, l denotes the total length, and curve stiffness is regulated by ξ > 0. The external cost factor c : R2 → R+ is added for adaptation to a given flat image data (see [40]). There are three motivations for our study. The first motivation comes from a medical image analysis application, where automatic extraction of the vascular retinal tree on images is helpful for early detection of many diseases, such as diabetic retinopathy, glaucoma, atherosclerosis, and others (see, e.g., [6, 18]). Optical retinal images are mostly acquired by flat cameras, and as a result, distortion appears. Such distortion could lead to questionable (distorted) geometrical features (vessel curvature, thickness, etc.) that are used as biomarkers [6, 10, 11] for different diseases. We will show that the distortion that appears near the boundary of a flat image can play a significant role in the quantitative analysis of the vascular structure and its curvature. As the retina itself is not flat, we should include the spherical geometry both in image processing and in image representation in order to avoid distortion. Such distortion comes from the central projection of the physical retinal surface to the image plane. It is illustrated in Fig. 1, where ( X, Y ) denote Cartesian camera coordinates and (x¯ , y¯) denote spherical object coordinates. This motivates us to study data-driven SR geodesics on the rigid body rotation group SO(3) and their spherical projections; likewise, it was done (see [39, 40]) for flat images and lifting to the roto-translation group SE(2). The second motivation comes from models of the visual system of mammals. As mentioned by Boscain and Rossi [15], the problem Pcurve(S2) can be considered as a spherical extension of a (flat) cortical-based model Pcurve(R2) for perceptual completion, proposed by Citti and Sarti [14], and Petitot [13]. Such a spherical extension is again motivated by the fact that the retina is not flat. By the same argument, cuspless SR geodesics in SO(3) could provide a model of association fields in the psychology of vision (see [31]). Note that nonuniform distribution of photoreceptors on the retina can be included in our model by putting an external cost C on a hemisphere. The third motivation for this study is that in geometric control theory optimal synthesis for the SR problem in SO(3) has not been achieved in the general case (not even for the case of uniform cost C = 1), despite many strong efforts in this direction [15,19–23,44,45]. In this paper, we will not provide optimal synthesis analytically, but instead we do provide a HJB theory for computing globally optimal (data-driven) geodesics. In previous works [39,40], we achieved this for SR geodesics in SE(2), used for tracking of blood vessels in flat 2D images. In view of these three motivations, we lift the problem Pcurve on the set S2 to a SR problem on the group SO(3), which we will call Pmec, and we provide explicit formulas for SR geodesics. This allows us to describe the set of endpoints in SO(3) reachable by geodesics whose spherical projections do not contain cusps. Furthermore, we present a Hamilton– Jacobi–Bellman PDE theory that allows us to numerically compute the SR distance map, from which a steepest descent backtracking (via the Pontryagin maximum principle) provides only the globally optimal geodesics for general external cost and general ξ > 0. We verify our numerical solution, by comparison with exact geodesics in the case C = 1. Finally, we use these results in a vessel tracking algorithm in spherical images of the retina, without central projection distortion. The main contributions of this article are: • New formulas for the geodesics of Pcurve(S2) for the uniform cost case C = 1. • Analysis and parametrization of cusp points arising in spherical projections of geodesics of Pmec(SO(3)). • A new HJB-PDE theory for the numeric computation of globally optimal geodesics of Pcurve(S2) and Pmec(SO(3)) for the general external cost case, with verification for C = 1. • Tracking of lines in spherical image data (e.g., retinal images) without central projection distortion, with comparison to tracking of lines in flat image data. 1.1 Vessel Tracking in Spherical Images of the Retina In this subsection, we first describe the mapping from object coordinates on the retina to camera coordinates, and then we discuss the relevance of considering spherical images of the retina rather than flat images, which are commonly used in the medical application (see Fig. 1). We show that distortion appears inevitably on flat images, with a significant relative error (up to 7%) in length measures. Even larger relative errors (over 20%) appear in the application of differential operators (used for vessel detection). Fig. 3 Top: schematic eye [42] and central projection of images onto the retina. Here R ≈ 10.5mm. Bottom: schematic eye, enlarged to support Eq. (1) We base our computations on the reduced schematic eye model (see Fig. 3), which is commonly used in clinical ophthalmology (see, e.g., [42]). Let R be the radius of an eyeball, a be the distance from the nodal point N to the center C , and ψ be the angle between visual axis and a light ray passing through N . Here we consider a simplified model, where the optical axis (the best approximation of a line passing through the optical center of the cornea, lens, and fovea) coincides with the visual axis (the line connecting fixation point and the fovea).1 The average radius of a human eye is R ≈ 10.5 mm, and the maximum distance between nodal point N and the central point C is amax = 17 mm − 10.5 mm = 6.5 mm. Now we switch to mathematical object coordinates of the retina where we use the eyeball radius to express lengths, i.e., we have R = 1 and amax = 10.5 mm · R = 2113 in dimensionless 6.5 mm coordinates. In order to compute the maximum absolute angle y¯max, let us express the angle |y¯| with respect to center of the eyeball (see right Fig. 3). Expressing the squared length of segment N Q yields (a + R cos |y¯|)2 cos−2 ψ = (a + R cos |y¯|)2 + R2 sin2 |y¯|. Solving this equation with respect to cos |y¯|, we obtain unique nonnegative solution: 1 There is small difference between these two axes (c.f. Fig. 33 [42]) which we neglect in our basic model. y¯max ≈ 0.63rad ≈ 36◦. We rely on the following parameterization of the image sphere S2 and the retinal sphere S¯2 (see Fig. 4): n(x , y) = (cos x cos y, cos x sin y, sin x )T , (−2a, 0, 0)T − (cos x¯ cos y¯, cos x¯ sin y¯, sin x¯)T , π , π2 ] and y, y¯ ∈ [−π, π ]. with x , x¯ ∈ [− 2 Next we present the explicit relation between the object coordinates (x¯ , y¯)—spherical coordinates on a unit sphere S¯2 representing the surface of the eyeball; flat photo coordinates (X, Y )—Cartesian coordinates on the image plane and the spherical coordinates (x , y) on the image sphere, see Figs. 1 and 4. To take into account the distance from the eyeball to the camera in our model, we introduce a parameter η > 0. In Fig. 4, by setting η = 1 we fix the distance equal to (a + c) radiuses of the eyeball. This corresponds to the case when the image sphere S2 is obtained by reflection of the physical retinal sphere S¯2 through the point (−a, 0, 0)T ∈ R3. In this case, we have (x , y) = (x¯, y¯), and we will always rely on this identification in the sequel. The general case η > 0 can be taken into account by congruency and scaling X → η R X, Y → η RY . The central projection Π (cf. Fig. 1) from (x , y) to (X, Y ) including the scaling factor η > 0 (with physical dimension length in units of R) is given by The inverse mapping Π −1 from (X, Y ) to (x , y) for η = 1 is given by x = arcsin(X p¯(X, Y )), y = arg( p1(X, Y ) + i Y p¯(X, Y )), a(a + c) + Ξa,c(X, Y ) p¯(X, Y ) = (X 2 + Y 2) + (a + c)2 Fig. 4 Spherical object coordinates (x¯, y¯) on a retina, Cartesian camera coordinates (X, Y ) on a flat image of the retina, and spherical camera coordinates (x, y) on a spherical image of the retina p1(X, Y ) = with Ξa,c(X, Y ) = (X 2 + Y 2)(1 − a2) + (a + c)2. 13 In these formulas, we need to substitute a = amax = 21 and c = 45 < R = 1, depicted in Fig. 4, where we work in dimensionless coordinates. 1.1.1 Quantification of Local and Global Deformation The local deformation from spherical photo coordinates (x , y) to planar photo coordinates (X, Y ) is now given by the Jacobian J (x , y) = det ∂∂Xx (x , y) ∂∂Xy (x , y) ∂∂Yx (x , y) ∂∂Yy (x , y) (a + c)2 cos x (1 + a cos x cos y) 2 = (a + cos x cos y)3 η . In mathematical analysis, we can set η = 1; however, in experiments η is to be taken into consideration. Note that for η = 1 we have 0.77 ≈ J (0, 0) ≤ J (x , y) ≤ J (ymax, ymax) ≈ 1.1, showing that local deformation plays a considerable role and varies in (x , y). Next we consider the global distortion along the line x = 0. It is defined as GD(y) = |y−Y (0,y)| , and it has a maximum |y| when y = ymax. We have 0 ≤ GD(y) ≤ GD(ymax) ≈ 0.07, and we see the distortion up to 7% along the line x = 0. The same holds along the line y = 0. We conclude that it makes a considerable difference to study Pcurve(S2) or Pcurve(R2) in the retinal imaging application. In the sequel, we will write Pcurve instead of Pcurve(S2) as we will always be concerned with the case where the base manifold equals S2. 2 Problem Pcurve on S2 and Pmec in SO(3) In this section, we first provide preliminaries on group theoretical tools and notation, and then we show that the spherical projection of certain SR geodesics γ (·) on the Lie group SO(3) provides solution curves to the problem Pcurve on S2, which we formulate next. 2.1 Mathematical Foundation and Notations The Lie group SO(3) is the group of all rotations about the origin in R3. We shall denote a counterclockwise rotation around axis a ∈ S2 with angle ϕ via Ra,ϕ , in particular for rotations around standard axes e1 = (1, 0, 0)T , e2 = (0, 1, 0)T , e3 = (0, 0, 1)T . We use representation of SO(3) by 3 × 3 matrices ⎛ cx cy −s x cy sθ − sy cθ sy sθ − s x cy cθ ⎞ = ⎝ cx sy cy cθ − s x s y sθ −cy sθ − s x s y cθ ⎠ ,(6) s x cx sθ cx cθ where we denote cx = cos x , cy = cos y, cθ = cos θ , s x = sin x , sy = sin y, s θ = sin θ , and where π π (x , y, θ ) ∈ − 2 , 2 The Lie group SO(3) defines an associated Lie algebra so(3) = TId(SO(3)) = span( A1, A2, A3), ⎛ 0 0 0 ⎞ ⎛ 0 0 1 ⎞ ⎛ 0 −1 0 ⎞ A1 = ⎝ 0 0 −1 ⎠ , A2 = ⎝ 0 0 0 ⎠ , A3 = ⎝ 1 0 0 ⎠ , 0 1 0 −1 0 0 0 0 0 where TId(SO(3)) denotes the tangent space at the unity element e ∈ SO(3), represented by identity matrix Id, which corresponds to the origin in the parameter space e ∼ Id ∼ {(x , y, θ ) = (0, 0, 0)}. Remark 1 The formula for Ai follows by Id, A2 = − ∂ R(∂xx,y,θ) |Id, A3 = ∂ R(∂xy,y,θ) |Id. The nonzero Lie brackets are given by There is a natural isomorphism between so(3) and the Lie algebra L of left-invariant vector fields in SO(3), where commutators of the vector fields in L correspond to the matrix commutators in so(3) [R A, R B] = R[ A, B], A, B ∈ so(3), R ∈ SO(3) . We express L in matrix form as ⎧ X1(x, y, θ ) = −R(x, y, θ ) A2, ⎪ L = span(X1, X2, X3), ⎨ X2(x, y, θ ) = R(x, y, θ ) A1, Remark 2 Note that at the unity one has X1 Id = − A2, X2 Id = A1, X3 Id = A3. The formulas for Xi in (10) follow by (cf. Remark 1) X1 R = R (∂x R)|Id, X2 R = R (∂θ R)|Id, X3 R = R (∂y R)|Id. We also use the isomorphism between so(3) and R3 Ai ∼ ei , R Ai R−1 ∼ Rei , where Ai ∈ so(3), R ∈ SO(3), ei ∈ R3, i = 1, 2, 3. Note that (6) is a product of matrix exponentials: R(x , y, θ ) = ey A3 e−x A2 eθ A1 . We choose to rely on this parameterization to keep the analogy with previous SE(2) models [28,31,40] mentioned in introduction. 2.2 Statement of the Problem Pcurve Let S2 = {n ∈ R3 n = 1} be a sphere of unit radius. We consider the problem Pcurve (see Fig. 2), which is for given boundary points n0, n1 ∈ S2 and directions n0 ∈ Tn0 (S2), n1 ∈ Tn1 (S2), n0 = n1 = 1 to find a smooth curve n(·) : [0, l] → S2 that satisfies the boundary conditions n(0) = n0, n(l) = n1, n (0) = n0, n (l) = n1, where kg(s) denotes the geodesic curvature on S2 of n(·) evaluated in time s, and C : S2 → [δ, +∞), δ > 0, is a smooth function that we call “external cost.” Here the total length l is free and s = 0s 1 dσ = 0s n (σ ) dσ denotes the spherical arclength. Thus, we have n (s) = 1 and the Gauss–Bonnet formula kg(s) = n (s) · (n(s) × n (s)). Remark 3 In introduction, we provided two motivations for this problem, where the external cost C plays a different role. Firstly, there is a motivation coming from retinal images, where C is adapted to spherical image data. Secondly, there is a motivation from the modeling of human vision, where the nonuniform distribution of photoreceptors can be modeled by the external cost. Remark 4 Without loss of generality, we can fix the boundary conditions as n(0) = e1, n(l) = n1, n (0) = e3, n (l) = n1, since the problem is left-invariant w.r.t. the natural left action of SO(3) onto S2 for C = 1. But also for C = 1 we can fix these boundary conditions by shifting the external cost. This boundary value convention is used throughout this article for Remark 5 Following the previous works [30,31] in SE(2) group, we do not expect existence of minimizers for the problem Pcurve(S2) in general. 2.3 Statement of Pmec Problem in SO(3) We call Pmec the following SR problem in SO(3): R˙ = −u1 R A2 + u2 R A1, R(T ) = R f in, L(R(·)) := R ∈ SO(3), (u1, u2) ∈ R2, ξ > 0, with T > 0 free. The external cost C : SO(3) → [δ, +∞), δ > 0, is a smooth function that is typically obtained by lifting the external cost C from the sphere S2 to the group SO(3), i.e., C (R) = C(R e1). We study the problem Pmec for C = 1 (case of uniform external cost) in Sect. 3, but let us first consider some preliminaries. Remark 6 In Pmec, we only have two velocity controls u1 and u2 in a 3D tangent space TR(t)(SO(3)) (cf. Fig. 5). Remark 7 The existence of minimizers for the problem Pmec is guaranteed by Chow-Rashevskii and Filippov theorems on sub-Riemannian manifolds [37]. Moreover, the geodesics of Pmec are smooth, since there are no abnormal extremals (see Example 1.3.13 in [49]). Remark 8 Sub-Riemannian manifolds are commonly defined by a triplet (M, Δ, G), with manifold M, distribution Δ ⊂ T (M) and metric tensor G. In our case M = SO(3), Δ = span{R A1, R A2}, G(R˙ , R˙ ) = C 2(·) ξ u1 + u22 , 2 2 where the controls u1 and u2 are components of velocity vector w.r.t. moving frame of reference, see (17). The choice of the sub-Riemannian structure (21) is determined by the initial conditions n(0) = e1 and n (0) = e3 in (16). Remark 9 By virtue of Cauchy–Schwarz inequality, the problem of minimizing the sub-Riemannian length L is equivalent to the problem of minimizing the action functional E (R(·)) = with T fixed (see, e.g., [33]). Remark 10 In analogy with the SR problem Pmec in SE(2), cf. [30,31], we sometimes call X1 = −R A2 the “spatial generator” and X2 = R A1 the “angular generator”, despite the fact that X1 and X2 are both angular generators on S2. The problem Pmec(SO(3)) can be seen as a model of the Reeds-Shepp car on a sphere. The Reeds-Shepp car can move forward and backward and rotate on a place. The input u1 controls the motion along X1, and the input u2 controls the motion along X2, see Fig. 5. 2.4 Relation Between the Problems Pcurve and Pmec We call a spherical projection the following projection map from SO(3) onto S2 (see Fig. 6): R → R e1 ∈ S2. = n(x , y) ∈ S2. So we see that (x , y) are spherical coordinates on S2. Fig. 5 The controls u1 and u2 along the “spatial generator” X1 and the “angular generator” X2 (cf. Remark 10) Fig. 6 Left: illustration of the parameterizations used in Pmec and Pcurve. The rotations are parameterized by (12), i.e., by angles x , y, and θ of rotation about basis axes, see (6). Right: illustration of the Hopf fibration [45]—circle bundle with the base S2 In problem Pcurve, one is interested in a curve n(s) = R(t (s)) e1, which satisfies (13) and minimizes (14). Here R(t ) = R(x (t ), y(t ), θ (t )), and t (s) = Next we show that the spherical projection (23) of certain minimizers of Pmec provides solution of problem Pcurve. More precisely, this only holds for the minimizers whose spherical projection does not have a cusp. Notice that if u1 ≡ 0 then the trajectory of Pmec is projected at a single point on S2 which does not provide a solution to Pcurve. This allows us to define smax as smax = inf {s > 0 u1(t (s)) = 0}, so that tcusp = t (smax). The following theorem states that minimizers n(s) of Pcurve for s ∈ [0, smax) are given by spherical projection of the minimizers R(t ) of Pmec for t ∈ [0, tcusp). Here s denotes the spherical arclength parameter, defined by dds n(s) ≡ 1, and t denotes the SR arclength, defined by C ( R(t )) ξ 2u21(t ) + u22(t ) ≡ 1. Theorem 1 Let R(t ), t ∈ [0, T ] be a minimizer of Pmec parameterized by SR arclength, and let the corresponding optimal control (u1(t ), u2(t )) satisfy the inequality u1(t ) > 0 for all t ∈ [0, T ]. Set Then for such boundary conditions, Pcurve has a minimizer n(s), along which we have n(s) = R(t (s)) e1, t (s) = u1(t ) = ddts (t ), u2(t ) = kg (s(t )) ddts (t ), for 0 ≤ s ≤ l < smax, and T = t (l). Proof See Appendix 1. 3 The Sub-Riemannian Problem Pm1ec in SO(3) In this section, we study the problem Pmec in the special case of uniform external cost C = 1, which we call Pm1ec. This problem can be seen as a left-invariant SR problem in Lie group SO(3), and it was tackled by many authors (see [15, 20–23, 44, 45]) in different contexts. It is remarkable that the statement of the problem is very simple, but due to the high complexity of the formulas describing geodesics, the minimizers are still unknown for general ξ > 0. The case ξ = 1 corresponds to a symmetric SR structure in SO(3), and it was completely solved in [15, 20, 30, 44, 45]. In this section, we consider the general case ξ > 0 and derive analytic formulas for the geodesics using the parametrization (6) for SO(3). Further, in Sect. 5 these formulas allow us to verify our numerical approach to compute the global minimizers of Pmec, presented for the first time by knowledge of the authors. In our analytical study, first we introduce a coordinate chart in SO(3) and write the optimal control problem Pm1ec in this chart. Then we apply Pontryagin maximum principle (PMP), which is the first-order necessary condition for optimality [37], and we derive the explicit formulas for the sub-Riemannian geodesics. Afterward we explain the notion of sub-Riemannian wavefront and sphere. Finally, we discuss what is known about optimality of the geodesics. 3.1 Pm1ec as Optimal Control Problem on S2 × S1 In this section, we introduce the chart in SO(3) and formulate problem Pm1ec given by (17)–(20), see Sect. 2.3, as an optimal control problem on M = Sx2,y × Sθ1, where π , π2 ] × R/{2π Z} are the spherical coordinates, (x , y) ∈ [− 2 cf. (24), and θ ∈ R/{2π Z}. We use this specific chart to stress the analogy with the closely related problem in SE(2) group (see [30,31]). Consider a smooth curve γ (·) = (x (·), y(·), θ (·)) ∈ C ∞(R, Sx2,y × S1) departing from the origin, i.e., γ (0) = θ e = (0, 0, 0). In Sect. 2.1, we parameterized the group SO(3) by the angles x , y, θ , recall (6). A smooth curve γ (·) corresponds to a smooth curve R(·) in SO(3) defined by the one-parameter family of matrices depending smoothly on the parameter t ∈ R and satisfying the initial condition R(0) = Id. A tangent vector of R(t ) = R(x (t ), y(t ), θ (t )) is expressed as ∂ R ∂ R ∂ R R˙ = ∂ x x˙ + ∂ y y˙ + ∂θ θ˙. Therefore, the control system (17) can be rewritten as ⎛ x˙ ⎞ γ˙ = ⎝ y˙ ⎠ θ ˙ ⎛ cos θ ⎞ ⎛ 0 ⎞ = ⎝ − sec x sin θ ⎠ u1 + ⎝ 0 ⎠ tan x sin θ 1 Vector fields near the controls ui are two of three of the basis left-invariant vector fields in SO(3) expressed in our chart. We denote them by X1, X2, and the third basis left-invariant vector field is given by their commutator X3 = [X1, X2]. Thus we have ⎛ cos θ ⎞ ⎛ 0 ⎞ ⎛ X1 = ⎝ − sec x sin θ ⎠ , X2 = ⎝ 01 ⎠ , X3 = ⎝ tan x sin θ Then problem Pm1ec given by (17)–(20) for C = 1, taking into account Remark 9, is equivalent to the following optimal control problem: γ˙ = u1 X1 + u2 X2, (u1, u2) ∈ R2 γ (0) = e, γ (T ) = ν f in = (x1, y1, θ1) ∈ S2 × S1, E = 21 0 T ξ 2u21 + u22 d t → min, ξ > 0. Remark 11 A projective version of problem Pm1ec with ξ = 1 is studied in [16]. It is obtained by identification of antipodal directions θ and θ + π ; thus, its solution is given by minimum between two geodesics with terminal configurations h = (x1, y1, θ1) and h¯ = (x1, y1, θ1 + π ). The associated SR distance is given by The model in [16] is connected to problem Pm1ec via right lsihniefatrgtra→nsfgor·meaπtiAo3ne−aπ2=A2 .π2T−hexl,obca=lchya,rθt i=s oξb+taiπne.dTbhye existence of optimal trajectories whose spherical projection contains a cusp is mentioned in [16]. 3.2 Application of the Pontryagin Maximum Principle In this subsection, we apply the Pontryagin maximum principle (PMP) to problem Pm1ec given by (33)–(35) and derive the Hamiltonian system of PMP. Here we follow standard geometric control theory (see, e.g., [37]). Due to the absence of abnormal extremals (see Remark 7), we consider only the normal case, where the controldependent Hamiltonian of PMP reads as k=1 Here ν = (x , y, θ ), Xi are given by (32), ·, · denotes the action of a covector on a vector, and (d ν1, d ν2, d ν3) = (d x , d y, d θ ) are basis one forms. The maximization condition of PMP reads as where u(t ) denotes an extremal control and (λ(t ), γ (t )) is an extremal. Introduce left-invariant Hamiltonians hi = λ, Xi , i = 1, 2, 3. The maximization condition gives the expression for the extremal controls u2 = h2. Then the maximized Hamiltonian, which we call “the Hamiltonian”, reads as h˙i = {H, hi } = The vertical part (for momentum components hi ) of the Hamiltonian system of PMP is given by j=1 where {·, ·} denotes the Poisson bracket (see [37]). Using the standard relation between Poisson and Lie brackets {hi , h j } = λ, [Xi , X j ] , we get Thus the Hamiltonian system with the Hamiltonian H follows by the last three identities and (38): ⎧⎪ h˙1 = −h2h3, ⎧⎪ x˙ = hξ21 cos θ, ⎪⎨ h˙2 = ξ12 h1h3, ⎪⎨ y˙ = − hξ21 sec x sin θ, ⎪⎪⎩ h˙3 = 1 − ξ12 h1h2 ⎪⎩⎪ θ˙ = hξ21 sin θ tan x + h2 — vertical part, — horizontal part, with the boundary conditions h1(0) = h01, h2(0) = h02, h3(0) = h30, x (0) = 0, y(0) = 0, θ (0) = 0. In (39), the horizontal part follows from (36) and (31). It is known that the Hamiltonian system (39) is Liouville integrable [25,27], and it was explicitly integrated in [19,22]. In the next subsections, we classify the possible solutions by values of the parameter ξ , and we adapt the explicit solution to our coordinate chart (x , y, θ ) ∈ M, where we follow the analogy to the closely related problem in SE(2). This allows us to obtain simpler formulas than in [19,22]. 3.3 Classification of Different Types of Solutions Here we describe the domains of parameter ξ > 0, which correspond to different dynamics of the Hamiltonian system (39). It is well known (see, e.g., [22]) that this system has two integrals of motion that depend only on momentum components hi (see illustration in Fig. 7) M 2 = h12 + h22 + h32 - Casimir function. t = Fig. 7 Top: integral manifolds of the vertical part of (39). The red sphere is a coadjoint orbit M2 = h21 + h22 + h32, and the green cylinder is the Hamiltonian 2H = hξ212 + h22 = 1. Bottom: vector field plot on the (h2, h3) plane (Color figure online) H = 21 in (41), we use SR arclength parameter along extremal trajectories. The momentum trajectory h(t ) can be seen as the curve formed by intersection of the cylinder H = 21 and the sphere M = const (see the yellow line in Fig. 7). Elimination of h1 from (41) and (42) yields Since h2 2 ≤ 1, there exist the following types of solution of the vertical part of (39): 1) Elliptic for ξ < 1, M 2 > ξ 2; 2) Linear for ξ = 1, M 2 > ξ 2; 3) Hyperbolic for ξ > 1, M 2 = ξ 2, M > 1; 4) A point (0,0) for ξ < 1, M 2 = ξ 2; 5) A segment of a straight line for ξ = 1, M = 1; 6) ξTw2;o crossing segments of straight lines for ξ > 1, M 2 = 7) No solution otherwise. According to this classification, we will refer to the case ξ < 1 as the elliptic case, ξ = 1 as the linear case, and ξ > 1 as the hyperbolic case. 3.4 Geodesics in SR Arclength Parametrization Different values of the Hamiltonian H correspond to different speed along extremal trajectories γ (·). By fixing the value In this subsection, we provide explicit formulas for the SR geodesics, where we reexamine results from [19,22]. In the coordinates (β, c) ∈ R/{4π Z} × R, defined by ⎪⎩ c˙ = −r sin β, c(0) = 2ξh03 , 1 where r = ξ2 − 1 ∈ (−1, +∞). Thus r ∈ (−1, 0) corresponds to the elliptic, r = 0 to the linear, and r ∈ (0, +∞) to the hyperbolic case. The system (44) has a symmetry (β˜(t ) = β(t ) + π, r˜ = −r ), which allows us to study the elliptic case and obtain the hyperbolic case by symmetry observations. Note also that the linear case should be treated separately, but the equations in this case are much simpler. Therefore, in the remainder of this manuscript we restrict ourselves to the case r ≥ 0 ⇔ 0 < ξ ≤ 1. In the following theorem, we provide explicit formulas for SR geodesics in SR arclength parametrization. Theorem 2 A solution of the system (39) reads as h1(t ) = ξ cos β2(t) , h2(t ) = sin β2(t) , h3(t ) = ξc2(t) , x (t ) = arg( R121(t ) + R221(t ) + iR31(t )), y(t ) = arg(R11(t ) + iR21(t )), θ (t ) = arg(R33(t ) + iR32(t )), ⎛ R11(t) R12(t) R13(t) ⎞ R(t) = ⎝ RR2311((tt)) RR2322((tt)) RR2333((tt)) ⎠ = D0T ey˜(t)A3 e−x˜(t)A2 eθ˜(t)A1 , M 2 − (h02)2, x˜(t ) = arg (h01)2 + (h02)2 + (h03)2, t − Proof See Appendix 2. Note that (β(t ), c(t )), and y˜(t ) admit explicit expression in Jacobi elliptic functions and elliptic integrals of the first, second, and third kind [19]. We present plots of spherical projections (23) of extremal trajectories in Fig. 8. Here it is remarkable that the spherical projections of SR geodesics in SO(3) can represent cusps, which can be undesirable for image analysis applications. This motivates us to study SR geodesics before the first cusp in their spherical projections. For short, we call such curves “cuspless geodesics.” In Sect. 4, we describe the possible end conditions, which can be connected by a cuspless geodesic, and in Sect. 5, we provide PDE-based approach for solving the boundary value problem (BVP) for data-driven SR geodesics, where we reformulate the problem as a solution to the SR-eikonal equation. 3.5 Sub-Riemannian Wavefronts and Spheres A useful tool for studying geodesics in left-invariant optimal control problems is the exponential map2 that maps an initial momentum h(0) and a time t to the endpoint of corresponding geodesic γ (·) [i.e., the exponential map integrates the Hamiltonian system (39)]: Te∗(SO(3)) × R → SO(3), (h(0), t ) → γ (t ). Now we explain the notion of wavefront. By definition, the wavefront consists of endpoints of all the geodesics of the same length T : WF(T ) = Exp(h(0), T ) h(0) ∈ Te∗(S1O(3)), H (h(0)) = 2 2 Not to be confused with the exponential map from Lie algebra to Lie The outer surface of a wavefront forms a sub-Riemannian sphere, which is a set of endpoints equidistant from e: 4 Sub-Riemannian Geodesics in SO(3) with Cuspless Spherical Projection S(T ) = Exp(h(0), T ) | h(0) ∈ Te∗(SO(3)), = {g ∈ SO(3) | d(e, g) = T } , where d(q1, q2) denotes the sub-Riemannian distance between q1 and q2 and tcut denotes the cut time, i.e., an instance of time, when a geodesic loses its optimality. In the previous subsection, we computed SR geodesics or in other words extremal trajectories. It is known (see, e.g., [33,37]) that sufficiently short arcs of SR geodesics are SR length minimizers (optimal trajectories). It is also known that in general a geodesic loses its optimality after a cut point. The corresponding instance of time is called cut time, and the set of all cut points in configuration space forms the so-called cut locus. In general, it is complicated to derive cut loci [37]. 1) Local optimality is lost at a conjugate point (critical point of exponential map that integrates the Hamiltonian system); 2) Global optimality is lost at a first Maxwell point (when two geodesics meet with the same length for the first time). In problem Pm1ec, both the Maxwell set and the conjugate locus are not known. The estimation of the first Maxwell time was given in [19], but still neither the cut locus, nor the SR length minimizers are known yet for general ξ > 0. The cut locus in a special case of the bi-invariant metric (for ξ = 1) was obtained in [15]. In [22], conjugate locus in the Riemannian case was constructed, and by considering a SR metric as a limiting case of the Riemannian metric, the corresponding formulas for the conjugate locus in SR case could be obtained. Here we do not provide such a limiting procedure. Instead, motivated by applications, in Sect. 5 we propose a numerical solution to compute the SR length minimizers and spheres. In Sect. 4, motivated by the study [31], we also give some statements on optimality of cuspless geodesics. We show that in contrast to SE(2) case [30,31] there exist nonoptimal cuspless SR geodesics. In this section, we study cuspless SR geodesics. Such geodesics allow parametrization by spherical arclength, which leads to simpler formulas than using SR arclength parametrization (see Sect. 3). Spherical arclength parameterization breaks down, when the spherical projection of a geodesic exhibits a cusp for the first time. So a natural question arises to characterize the set of end conditions R ⊂ SO(3) reachable by cuspless geodesics, similar to the SE(2) case studied in [31]. The cuspless SR geodesics are projections in SO(3) of the trajectories of the Hamiltonian system (39) that are going in positive direction of X1 (i.e., u1(t ) > 0) before the first cusp time t < tcusp. Thus, by virtue of (36) the cuspless constraint is given by h1(t ) > 0, for all t ∈ [0, T ]. We shall often rely on short notation hi (s) := hi (t (s)), where we stress our notations in order to avoid confusion with the chain law for differentiation. Recall (26) and note that smax for a given geodesics is determined by its initial momentum. We write smax(h(0)) when we want to stress this dependence. We derive an explicit formula for smax(h(0)) in Sect. 4.1. Theorem 3 For any s ∈ [0, smax(h(0))), h1(0) > 0 the system (39) is equivalent to the following system (see corresponding vector flow in Fig. 7): with the boundary conditions Proof Express the dynamics (39) in the spherical arclength parameter s. We have h2(s) = dhd2s(s) = hh˙21((tt((ss)))) ξ 2 = h3(s), h2(0) = h20, h3(s) = dhd3s(s) = hh˙31((tt((ss)))) ξ 2 = ξ 2 − 1 h2(s), h3(0) = h30. Here we recall that u1 = ddst , see (28), and h1 = ξ 2u1, see (36). Thus we obtained the vertical part of (49). Doing the same for the horizontal part in (39), we obtain the full system (49), where the Hamiltonian H (s) = ξ−2h21(s2)+h22(s) = 21 together with the cuspless constraint h1(s) > 0 allows us to write h1(s) = ξ 1 − h22(s). Finally, the boundary conditions (50) follow from (40), since t = 0 ⇔ s = 0. Note that in contrast to the mathematical pendulum ODEs (44) expressed in SR arclength t , where the Jacobi elliptic functions appear, the vertical part in s parameterization now is a simple linear system of ODEs integrated in elementary functions. To obtain simpler formulas, we define the parameter χ ∈ C (where one should take principal square root): Theorem 4 The general solution of the vertical part in (49) for all ξ = 1 is given by 0 ⎧⎪ h2(s) = h20 cosh sχ + hχ3 sinh sχ , ⎪ ⎨ h3(s) = h30 cosh sχ + χ h02 sinh sχ , and for the case ξ = 1, we find straight lines parallel to h2 axis in the (h2, h3)-phase portrait: ⎪⎪⎩ h1(s) = 1 − h22(s). 4.1 Computation of the First Cusp Time To analyze the cusp points in problem Pcurve, we need to determine smax(h(0)). It is given, recall (26) and (36), by the minimal positive root of equation h1(s) = 0: smax(h(0)) = min{s > 0 | h1(s) = 0}. Theorem 5 When moving along a SR geodesic t → γ (t ), the first cusp time is computed as tcusp(h(0)) = t (smax(h(0)), recall (25), where with s1 = sgn Re h20χ + h30 and κ = (h03)2 + 1 − (h02)2 χ 2 ∈ R. As a result, in the SR manifold (SO(3), Δ, G), recall (21), there do exist nonoptimal cuspless geodesics. Proof See Appendix 3. The presence of nonoptimal cuspless geodesics is remarkable, as in the SE(2) group every cuspless geodesic is globally optimal [31]. 4.2 Geodesics in Spherical Arclength Parametrization In this subsection, we derive the formulas for cuspless SR geodesics in s parameterization. Theorem 6 The unique solution of (49) is defined for s ∈ [0, smax(h(0))], where smax(h(0)) is given by (54). The solution to the vertical part is given by Theorem 4, and the solution to the horizontal part is given by ⎛ R11(s) R12(s) R13(s) ⎞ R(s) = ⎝ RR2311((ss)) RR2322((ss)) RR2333((ss)) ⎠ = D0T ey˜(s)A3 e−x˜(s)A2 eθ˜(s)A1 , x˜(s) = arg M 2 − h22(s) + i h2(s) , θ˜(s) = arg h3(s) − i ξ 1 − h22(s) , M2 − (h02)2, M = Proof It follows from Theoremes 2, 3, and 4. Momentum component h1 is expressed in h2 via the Hamiltonian (37) which equals to 21 along SR geodesics. The sign of h1(0) is equal to the sign of u1(0) = ddst (0) and therefore nonnegative. Fig. 9 Spherical projection of cuspless SR geodesics in SO(3) in elliptic, linear, and hyperbolic cases Regarding the integration constraints, we note that M 2 = h1 + h2 + h32 and 0 ≤ h1 = ξ 1 − h22 and 2 2 M 2 − h22 from which the result follows. Remark 13 In contrast to general SR geodesics in SO(3) given in t parameterization, where Jacobi elliptic functions appear, the cuspless SR geodesics in SO(3) given in s parameterization involve only a single elliptic integral s 1−h22(σ ) 0 M2−h22(σ ) dσ . This integral can be expressed in terms of standard elliptic integrals. For example, in the elliptic case ξ < 1, M 2 > ξ 2 we have M1−2−ξξ22 , Ψ (s) = s 1 − ξ 2 + Ψ0, with Ψ0 = = √ where F denotes an elliptic integral of the first kind and Π denotes an elliptic integral of the third kind. See plots of projected cuspless geodesics n(s) = R(s) e1 for s ∈ (0, smax) in Fig. 9. 5 PDE Approach for Data-Driven SR Geodesics in SO(3) In this section, we adapt the PDE approach for data-driven SR geodesics in SE(2) [39,40] to the SO(3) group. Here we consider the basis left-invariant vector fields Xi as differential operators of the first order, and we write Xi (W) for the derivative of a function W : SO(3) → R along Xi . We aim to solve the following optimal control problem: i=1 γ (t ) ∈ SO(3), γ (0) = e, γ (T ) = g1, (u1, u2) ∈ R2, T Here the terminal time T is free; and C : SO(3) → [δ, +∞), δ > 0 is the external cost. t By rescaling of time τ = T ∈ [0, 1] simultaneously with controls u¯i (τ ) = T ui (T τ ), we write down the explicit solutions as Here the sub-Riemannian metric tensor G = C 2(·) ξ 2ω1 ⊗ ω1 + ω2 ⊗ ω2 5.1 Sub-Riemannian Fast Marching in SO(3) Here we propose a SR-FM method (sub-Riemannian fast marching) for the computation of data-driven SR length minimizers (not necessarily cuspless) in SO(3) group, as a solution to the SR-eikonal system (62). This method was successfully used in [41] for the computation of data-driven SR length minimizers in SE(2) group. The method is based on a Riemannian approximation of sub-Riemannian manifold, and computing Riemannian geodesics in highly anisotropic space, which becomes the SR manifold in the limiting case as anisotropy tends to infinity. Here we follow the explanation in [41], where we work now in new settings of the SO(3) group and use the coordinate chart (x , y, θ ) defined in Sect. 3.1. Recall that the basis leftinvariant vector fields3 Xi in SO(3) are given by the following 3 In previous works [39,40], Xi was denoted by Ai differential operators: X1 = cos θ ∂x − sec x sin θ ∂y + tan x sin θ ∂θ , X3 = sin θ ∂x + sec x cos θ ∂y − tan x cos θ ∂θ , and corresponding basis left-invariant one forms ωi , satisfying ωi , X j = δij , are expressed as We approximate the sub-Riemannian manifold by a Riemannian manifold by fixing a small ε > 0. Moreover, the SR-eikonal equation (62) is well defined and it can be derived as a limiting case of the eikonal equation on a Riemannian manifold via the inverse metric tensor, see Appendix 4. Let us denote the Riemannian distance as follows: Gε := C2(·) ξ 2ω1 ⊗ ω1 + ω2 ⊗ ω2 + ξε22 ω3 ⊗ ω3 . (60) Remark 14 In our approach, we shall rely on standard notion of viscosity solution [7,8,12] of the Riemannian eikonal equation, which admits a generalization to the sub-Riemannian case. See details in Appendix 4. The following theorem summarizes our approach for the computation of data-driven sub-Riemannian length minimizers in SO(3). Theorem 7 Let g = e ∈ SO(3) be chosen such that there exists a unique minimizer γε∗ : [0, 1] → SO(3) of dSεO(3)(g, e) such that γε∗(τ ) is not a conjugate point for all τ ∈ [0, 1] and all ε ≥ 0 sufficiently small. 0 Then τ → dSO(3)(e, γ0∗(τ )) is smooth and the minimizer γ0∗(τ ) is given by γ0∗(τ ) = γb∗(1 − τ ), with γ˙b(τ ) = −u¯1(τ ) X1|γb(τ ) − u¯2(τ ) X2|γb(τ ) , and with W(g) denoting the viscosity solution of the following boundary value problem: ⎩ W(e) = 0. 2 = C (g), for g = e, For an outline of the proof, see Appendix 4, where we rely on a similar approximation approach as in [48, ch. 5, app. A] for an elastica functional and in [36, Thm. 2, Cor. 2, app. A] for a SR problem on S E (2). Remark 15 The approach in Theorem 7 can be adapted for producing only the cuspless minimizers. See Appendix 4. Remark 16 The SR spheres of radius t centered at e are given by St = {g ∈ SO(3) | W(g) = t }. Remark 17 Recall that conjugate points are points where local optimality is lost. For a formal definition and characterization, see [35, def. 8.4.3, cor. 8.4.5]. We would conjecture that the assumption on the conjugate points in the above theorem is not really needed, as even the conjugate points that are limits of first Maxwell points do not seem to cause problems in the backtracking procedure, akin to the S E (2) case [40, App. D], but nevertheless our proof in Appendix 4 does rely on this assumption. From the derivations above Theorem 7, we see that the fast marching approach for computing SR geodesics in SE(2), cf. [41], is easily generalized to the SO(3) case. To this end, we replace the matrix representation for Gε expressed in the fixed (x , y, θ ) Cartesian coordinate frame. In the SO(3) case, it equals ⎛ cos θ 0 sin θ ⎞ R = ⎝ − cos x sin θ sin x cos x cos θ ⎠ . 0 1 0 Here the diagonal matrix in the middle encodes the anisotropy between the Xi directions, while the matrix R is the basis transformation from the moving coframe {ω1, ω2, ω3} to the fixed coframe {dx , dy, dθ }, recall (58), in which the fast marching implementation via special anisotropic stencils [5] is used. In Sect. 6.1, we will show that the thereby obtained fast marching approach already presents reasonable precision for ε = 0.1. Experiments in Sect. 6.3 show the application of the method (with data-adaptive nonuniform cost) to tracking of blood vessels in retinal images. 6 Experiments In Sect. 6.1, we verify the SR-FM method by comparison of SR length minimizers/spheres obtained via SR-FM with the exact SR geodesics/wavefronts (cf. Sect. 3) for the case of uniform external cost (i.e., C = 1). In Sect. 6.2, we compare SR geodesics and wavefronts in the groups SE(2) and SO(3) (Sect. 6.2.1) for the uniform external cost case. In Sect. 6.3, we provide experiments of vessel tracking by SR geodesics in SO(3) when the external cost C is induced by spherical data, and compare them to the result of vessel tracking on the corresponding flat image by SR geodesics in SE(2). 6.1 Verification of the Fast Marching Method in the Case of Uniform External Cost In this subsection, we perform the experiments to validate the SR Fast Marching (SR-FM) method proposed in Sect. 5.1. The goal of the experiments is to check that the method produces an accurate approximate solution to the SR problem in SO(3) group in the case of uniform external cost C = 1. In all the experiments, we fixed the anisotropy parameter of the Riemannian approximation as ε = 0.1. In the first experiment, we compare the geodesics γ F M (·) obtained via SR-FM with the exact cuspless geodesics γ (·) computed via analytic formulas in Theorem 6. We perform the comparison as follows: 1) Fix ξ > 0 and the initial momenta h(0). 2) Compute the first cusp time smax (h(0)) corresponding to h(0), and set send = min{smax (h(0)), π2 }. 3) Compute the geodesic γ (s), s ∈ [0, send] via Thm. 6. 4) Compute the distance function W(g) in the domain π , π2 ] × [−π, π ] × [−π, π ] via SRg = (x , y, θ ) ∈ [− 2 FM. Here we compute the distance function in the grid of 201 × 401 × 401 points and then interpolate it using third-order Hermite interpolation. 5) Compute the geodesic γ F M (t ), t ∈ [0, W(gend)] via backtracking (61) from the endpoint g1 = γ (send). 6) Plot the spherical projections of γ (s) and γ F M (t ) and compare them. A typical result of the comparison is shown in Fig. 10, where we set ξ = 1.5, h3(0) = 0 for all the curves and varied initial momentum h2(0) ∈ {−0.99, −0.81, −0.63, −0.45, −0.27, −0.09, 0.09, 0.27, 0.45, 0.63, 0.81, 0.99}. As a result, we see that the geodesics computed numerically via SR-FM accurately follow the exact geodesics. We have performed a series of such experiments and always obtained similar results when the geodesics γ i (s) were optimal for s ∈ [0, seind]. It was also remarkable that SRFM resulted into different curves (length minimizers) when Fig. 10 Top: the spherical projection of cuspless SR geodesics in SO(3) computed by exact formulas in Theorem 6 (green dashed lines) and our numerical fast marching approximations (red lines). Bottom: example of nonoptimal cuspless geodesics (Color figure online) the geodesic γ i (s) was not optimal for s ∈ [0, seind]. Such an example is illustrated in Fig. 10 (bottom). The question of optimality of SR geodesics in SO(3) in the general case ξ > 0 is an open important problem, recall Sect. 3.5. Here we provide a numerical SR-FM method for computing only the optimal geodesics. In analogy with how it was done in [40], it is possible to compute Maxwell sets numerically. We have a strong conjecture that optimality is lost at the first Maxwell point induced by reflectional symmetries along the geodesic. We leave the computation of Maxwell points and analysis of optimality of the geodesics as a direction for future work. To verify the SR-FM method we also perform experiments with comparison of the geodesics obtained via SR-FM and the general geodesics (not necessarily cuspless) given by Theorem 2. A typical result is presented in Fig. 11 (top). Here we again observe an accurate result of SR-FM, but now for the geodesics whose spherical projections have cusps. To conclude this section, we provide one more experiment, where we compare the exact sub-Riemannian wavefronts and SR spheres computed via SR-FM. We show that the SRFM method provides a distance function W(g) that closely approximates the SR distance d(e, g) ≈ W(g). In Fig. 11, we show the comparison of the exact wavefront WF( 3125 π ) and Fig. 11 Top: comparison of exact SR geodesics in SO(3) (green dashed lines) obtained via analytic formulas (Theorem 2) and SR geodesics computed numerically via SR-FM (red lines). Here ξ = 4.5, tend = 32π , and initial momenta (hi2(0), hi3(0)) = {(−5, −0.99), (−5, 0.99), (5, −0.99), (5, 0.99)}. Bottom: Comparison for ξ = 1 of wavefront in SO(3) (green transparent surface) obtained via the analytic formulas of Theorem 2 for T = 3125 π and the SR sphere (orange solid surface) for the same T , obtained via SR-FM (Color figure online) isosurface S F M ( 3125 π ) = {g ∈ SO(3) |W (g) = 3125 π }. Here we see that the isosurface computed via SR-FM accurately follows the outer surface of the exact wavefront (i.e., exact SR sphere). A similar picture was obtained for different radii T . 6.2 Comparison of SR Geodesics and Wavefronts in SE(2) and SO(3) for C = 1 In this subsection, we compare SR geodesics and wavefronts in SE(2) and SO(3) and analyze their applicability in image processing. We also include a short discussion about optimality of geodesics, which is closely related to the analysis of self-intersections of wavefronts and the study of cut loci. To this end, we provide accurate plot of the wavefronts near their singularities. 6.2.1 Comparison of SR Wavefronts in SO(3) and SE(2) for C = 1 In this subsection, we show a comparison of the SR wavefronts in SO(3) and SE(2) in the case of uniform external cost C = 1. Here we employ the fact that the coordinate chart (x , y, θ ) in Sx2,y × Sθ1 was chosen in the way to obtain the analogy with ( X, Y, Θ ) ∈ R2 × S1. This allows us to plot SR wavefronts of SO(3) and SE(2) in the same 3D plot. In this comparison, we use the SR arclength parameterization t , where the geodesic γ SO(3)(·) is given by exponential map Exp(h0, ·), recall (47), i.e., by the projection in SO(3) of the solution of the Hamiltonian system (39), with initial condition h(0) = h0, γ SO(3)(0) = e = (0, 0, 0). To establish the comparison of wavefronts we switch to polar coordinates (β, c), recall (43), where the Hamiltonian system in SO(3) reads as: — vertical part, A solution to this system is given by Theorem 2. The Hamiltonian system for SR geodesics in SE(2) (see, e.g., [29] ) reads as: — vertical part, A solution to this system is given in [29]. Indeed, we observe a clear similarity between the SO(3) and SE(2) case, when using parametrization (6). In Fig. 12, we plot the SR wavefronts in SO(3) and SE(2) for several values of end time T . In these plots, we identify (x , y, θ ) and ( X, Y, Θ ), so that the red surfaces WFSE(2) are the SR wavefronts of SE(2) and the green surfaces WFSO(3) are the SR wavefronts of SO(3). We see a very similar shape of WFSO(3) and WFSE(2) for small radii T , but the difference increases when T increases. Thus, we conclude that the SR geodesics in SO(3) can be locally approximated by the SR geodesics in SE(2), but globally they are considerably different. One can also observe that the singular points4 of WFSO(3) are located near singular points of WFSE(2). This leads us to a conjecture that the location of conjugate points (open problem) can be estimated by the location of conjugate points 4 By singular points, we mean either conjugate or Maxwell points (recall Remark 12). Fig. 12 Comparison of SR wavefronts in SE(2) (red) and SO(3) (green) for T = 0.5, ξ = 1, η = 1. A zoomed in picture for the blue square is provided in Fig. 13 (Color figure online) of WFSE(2). Note also that in the general case ξ > 0 the set of singular points has a complicated shape. In Fig. 13, we present more detailed plot of the singularities, with the depicted Maxwell set and special cases of conjugate points, which are limit points of the Maxwell set on SR sphere (outer part of wavefront). In Fig. 13, we observe that the wavefront in SO(3) has a very special symmetry when ξ = 1. This is not present for ξ = 1, and this never happens in the SE(2) group. The case ξ = 1 for SR manifold in SO(3) was completely examined in [15], where it was shown that locally the conjugate locus is an interval, and globally it is a circle without a point. Changing ξ destroys the symmetry, conjugate and Maxwell points are getting separated, and the conjugate locus has an astroidal shape [34]. 6.2.2 Comparison of SR Geodesics in SO(3) and SE(2) for C = 1 In this subsection, we again consider the case C = 1 and compare SR geodesics γ SO(3)(·) = (x (·), y(·), θ (·)) and γ SE(2)(·) = (X (·), Y (·), Θ(·)) in the image plane. The SR-FM method is used for computation of the geodesics parameterized by SR arclength. Here we prepare background for comparison of the geodesics in retinal images via the schematic eye model, recall Sect. 1.1, where as a departure point we use an image (white for C = 1) on a plane OXY , recall Fig. 4. See Fig. 14, where we compare SE(2) and SO(3) SR geodesics projected on the plane and on the sphere (via mappings Π and Π −1). For details, see Appendix 5. 6.3 Vessel Analysis via SO(3) SR Geometry and SE(2) SR Geometry As explained in introduction in Sect. 1.1, we need to include the spherical geometry of the retina rather than the flat geometry of the flat image. This spherical geometry is encoded in our spherical image model, see Fig. 4. Next we will analyze the effect of including this geometry in the SR-FM vessel tracking method along data-driven SR geodesics in SO(3). More precisely, we propose vessel tracking in object coordinates (or spherical image coordinates) via SR geometry in SO(3) as an extension of vessel tracking in flat images [39,40] along SR geodesics in SE(2). Therefore, we want to investigate whether including the correct spherical geometry makes a difference in the vessel tracking in practice. Although a complete detailed comparison on large data sets is left for future work, we present preliminary experiments which indeed indicate considerable differences in both tractography and curvature measurements. These experiments are shown in Figs. 15 and 16 and next we explain them. We apply the same scheme as in Sect. 6.2.2, but now we compute data-driven geodesics, where the external cost is induced by image data. Next we explain the construction of the external cost that is applied in all experiments. For the sake of simple comparison, we restrict ourselves to a cost depending on the spherical coordinates only, and we set G(x , y) = V F (Π (x , y)) −1 where we use standard multiscale vesselness [50] with λi (X, Y ) eigenvalues of the Gaussian Hessian of image F : R2 → R (maximized over scales) ordered by |λ1(X, Y )| < |λ2(X, Y )|, with β = c = 0.3, and with unit 1, for λ ≥ 0, step function U (λ) = 0, otherwise. Here we note that the Gaussian Hessian is given by H(Gs ∗ F ) and computed by Gaussian derivatives [51] at multiscales s = 21 σ 2 ∈ {2, 3, 4, 5} in term of pixel sizes. In the experiment in Figure 15, we show that there is a considerable difference between SE(2) SR geodesics and SO(3) SR geodesics. We see that when internal geometry is dominant over the external cost (λ small) the SR geodesics in SO(3) are more stiff than SR geodesics in SE(2), and therefore in the boundary value problem they are less eager to take shortcuts and better follow the vessel structure. In case λ is large (external cost is dominant than the internal geometry), we see only small differences in the overall locations of the SE(2) curves and SO(3) curves. The results are stable w.r.t. choice of 1 ≤ η ≤ 2 (which controls the distance from the camera to the eye ball, relative to eye ball radius, recall Fig. 4). In the next experiment, we measure the curvature of the curves obtained by the vessel tracking method via SE(2) geometry and via SO(3) geometry. For this experiment, we Fig. 13 Comparison of self-intersections of SR wavefronts in SE(2) (red) and SO(3) (green) for ξ = 1 (linear case) and ξ = 0.5 (elliptic case). Here T = 0.5 and η = 1. The viewpoint is taken from the inside of the SR sphere. In this figure, we plot the part of wavefront depicted by the blue square in Fig. 12. If ξ = 1 both the wavefronts in SE(2) and in SO(3) do not intersect at a single point. The first Maxwell sets are depicted by dashed violet lines. The first conjugate points on the SR sphere are depicted by yellow dots (Color figure online) Fig. 14 Comparison of an SE(2) SR geodesic (blue) and an SO(3) SR geodesic (red), from left to rightin (x, y) spherical image coordinates, (X, Y ) flat image coordinates, and plotted on the sphere S2 (Color figure online) used the values ξ = 3, λ = 50 and η = 2. Although in this case the result of tractography is very similar for the SE(2) and SO(3) curves, we show that there is a considerable difference in their curvature. Corollary 1 (from Theorem 1) The geodesic curvature of a spherical projection of data-driven geodesic γ SO(3)(·) satisfies Proof The first equality in the chain 2 X2|γ (·)(W ) X1|γ (·) (W ) is implied by Eq. (28). The second equality follows from h1 application of PMP to problem Pmec, which gives u1 = ξ 2C2 and u2 = Ch22 . The third equality follows from the fact that in the points where the Bellman function W is differentiable (almost everywhere in our case) its derivatives are given by components of momentum covector h1 = X1(W ), h2 = X2(W ), see [37]. By the same argument, it can be checked (see [31] and [40]) that the planar curvature of spatial projections of SR geodesics in SE(2) satisfies κ SE(2)(·) = ξ 2 A2|γ SE(2)(·) W with A1 = cos Θ ∂X +sin Θ ∂Y , A2 = ∂Θ basis left-invariant vector fields SE(2). Thus the curvature analysis can be simply done based on vessel tracking, and this shows the benefit of our algorithm. In Fig. 16, we show an experiment of vessel curvature measurement based on tracking via SO(3) geometry and via SE(2) geometry. For completeness, we added also a comparison with a planar curvature κ SO(3)(·) of a planar projection Γ SO(3)(·) := Π (x (·), y(·)) of γ SO(3)(·) = (x (·), y(·), θ (·)). Fig. 15 Comparison of vessel tracking via γ SO(3) SR geodesics in SO(3) (green solid lines) in object coordinates and via γ SE(2) SR geodesics in SE(2) (red dashed lines) in the planar camera coordinates. Here the planar projection Γ SO(3) of γ SO(3) and spatial projection Γ SE(2) of γ SE(2) are depicted in the same flat image (Color figure online) Fig. 16 Left: two curves from the experiment in Fig. 15 (right-bottom figure) are depicted with slight shift. The upper curve is a spatial projection Γ SE(2) of the data-driven SR geodesic γ SE(2) with depicted (in color) planar curvature κSE(2) on top of the curve. The lower curve Γ SO(3) is the planar projection of the SR geodesic γ SO(3) with depicted geodesic curvature κgSO(3) on top of it. Right: three graphs are shown in the same plot: planar curvature κSE(2) of Γ SE(2); geodesic curvaΓ SO(3). The effect of considering geodesic curvature κgSO(3) in object coordinates on S2 rather than planar curvature κSO(3) in photo coordinates on projection on R2 is visible (compare the green solid and dashed graphs). A bigger difference comes from using SO(3) SR geometry than SE(2) SR geometry (compare red and green graphs) (Color figure online) We can see a considerable difference in curvature measurement via SO(3) geometry and SE(2) geometry. It is also seen that the difference between κ SO(3) and κgSO(3) is not very significant. Thus we see the importance of including correct spherical geometry in vessel tracking algorithm in retinal images. Fig. 17 Comparison of vessel tracking in a spherical image of the retina via SR geodesic (green solid line) and Riemannian geodesic (red dashed line) in SO(3) (Color figure online) Next we address the general benefit of using SR geodesics rather than (isotropic) Riemannian geodesics in the tracking of salient lines (blood vessels) in spherical images. To this end, we show a typical example where the tracking induced by a sub-Riemannian geodesic gives better result than the tracking induced by a Riemannian geodesic. The experiment in Fig. 17 is performed similar to Sect. 6.3, but now we compare the result of vessel tracking via SR geodesics (obtained by SR-FM) and Riemannian geodesics (isotropic metric with ε = 1). From this experiment, we see that similarly to the SE(2) case [40] the tracking via Riemannian geometry suffers from incorrect jumps toward nearly parallel neighboring vessels and yields nonsmooth curves, the tracking via SR geometry gives the desirable result. 7 Conclusion Data-driven sub-Riemannian geodesics in 3D Lie groups are a suitable tool for tractography of blood vessels in retinal imaging. In previous works on the SE(2) case [40,41], practical advantages have been shown in comparison with the (isotropic) Riemannian case, and geodesic methods in the image domain. However, these models included a SR geometry on SE(2) ≡ R2 S1 based on lifts of flat images, which does not match the actual object geometry: The retina is spherical rather than planar, cf. Fig. 1, Fig. 3, and Fig. 4. A similar observation holds for models in the psychology of contour perception in human vision [30]. Therefore, for geometric tracking we propose a frame bundle above S2, cf. Fig. 6, instead of a frame bundle above R2. Geometric tracking of geodesics is done along globally optimal data-driven SR geodesics in SO(3) (and their spherical projections) by our new numerical wavefront propagation method. The method was validated for the uniform cost case by comparisons with exact geodesics which we derived in Sect. 3, cf. Figs. 10 and 11. Here, in contrast to the SE(2) case [31], we do not have a scaling homothety. As a result, the parameter ξ has a considerable effect on the geometry, and for ξ = 1 we identify the linear case, for ξ > 1 we identify the hyperbolic case, and for ξ < 1 we identify the elliptic case, cf. Figs. 7 and 9. For all of these cases, we have computed the first cusp time in Theorem 5 (the first time where the spherical projection of a SR geodesic exhibits a cusp). Also, we have presented new formulas for such “cuspless” SR geodesics in Theorem 6. These formulas only involve a single elliptic integral thanks to spherical arclength parametrization, and for ξ = 1 our formulas are simpler than the general formulas for SR geodesics in SO(3). Furthermore, we used a specific parametrization of Lie group SO(3) that allowed us to compare between SR geodesics / wavefronts in SO(3) to SR geodesics / wavefronts in SE(2), cf. Figs. 12 and 13. In our comparison we took into account a standard optical model for the mapping between object coordinates on the retina and camera coordinates in the acquired planar retinal image. In our experiments, the differences between the SO(3) case and the SE(2) case are considerable, both for the case of uniform cost, cf. Fig. 14, and for the data-driven case in the retinal image analysis application, cf. Fig. 15. In general, we see that for realistic parameter settings (in optics) the SO(3) geodesics have a slower variation in curvature and are less eager to take shortcuts, see, e.g., Fig. 15 and 16. Furthermore, there are visible differences between geodesic curvature of data-driven SR geodesics on the sphere and the curvature of their planar projections. As in retinal imaging applications, curvature is considered as a relevant biomarker [9–11] for detection of diabetic retinopathy and other systemic diseases, the data-driven SR geodesic model in SO(3) is a relevant extension of our data-driven geodesic model in SE(2). Here we restricted ourselves to feasibility studies. More extensive comparisons between SO(3) geodesics and SE(2) geodesics on large retinal imaging benchmark sequences are beyond the scope of this article and are left for future work. Finally, we note that the computation time for data-driven SR geodesics in SO(3) is exactly the same as for the SE(2) case. Our specific choice of coordinates of SO(3) allowed us to modify the very efficient fast marching approach [41], with a simple replacement of the metric tensor matrix. 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, a.n. 335555. The authors gratefully acknowledge EUMarie Curie project MANET ag. no. 607643, and RFBR research Project No. 16-31-60083 mol_a_dk for financial support. The publication was supported by the Ministry of Education and Science of the Russian Federation (the Agreement number 02.a03.21.0008). The authors gratefully acknowledge dr. G.R. Sanguinetti and Dr. J.-M. Mirebeau on their efforts on anisotropic fast marching implementation of SR geodesics in SE(2) [41] that allowed for our generalization to the SO(3) case. Furthermore, the authors thank G.R. Sanguinetti on helpful comments and suggestions on earlier versions of this article. We thank Dr. Tos Berendschot for fruitful discussions on retinal imaging and for his encouragement to take into account the spherical nature of the retina in our geometric retinal image processing. The authors thank anonymous reviewers whose valuable comments improved the exposition of this work. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Appendix 1: Proof of Theorem 1 The function s(t ) := increasing for t ∈ [0, T ]; then there exists an inverse function t (s), s ∈ [0, l], l = s(T ), differentiable and increasing. Notice that ddst (t ) = u1(t ), which gives the second equality Define R¯ (s) := R(t (s)), n(s) := R¯ (s)e1, n(t ) := R(t )e1 = n(s(t )). From the dynamics of Pmec, we obtain identities Then the Gauss–Bonnet formula (15) gives which implies the third equality in (28). The boundary conditions (27) for the curve n(s) follow from its definition and the boundary conditions for R(t ). Recall that by definition C ( R(t )) = C(n(s(t ))). Since the minimizer R(t ) is parameterized by SR arclength, we have C 2( R(t )) ξ 2u21(t ) + u22(t ) ≡ 1, whence ddst (s) = u1(t1(s)) = implies (29) after integration w.r.t. s. Thus, we can see that the optimization functionals of Pcurve and Pmec coincide: 0T C (R(t )) ξ 2u21(t ) + u22(t )d t = 0l C(n(s)) ξ 2 + kg2(s)d s. But also the dynamics coincides, in the following sense. If n˜ (s), s ∈ [0, l˜], is a smooth curve on the sphere S2 with the initial conditions n˜ (0) = e1, n˜ (0) = e3 and a geodesic curvature k˜g (s), then it can be lifted to a curve R˜ (s) in SO(3) that is a trajectory of control system (17) with the controls u˜ 1(s) ≡ 1, u˜ 2(s) = k˜g (s) and the initial condition R˜ (0) = Id. The curve R˜ (s)e1 on S2 has the same initial conditions and geodesic curvature as n˜ (s), thus R˜ (s)e1 ≡ n˜ (s). Summing up, since R(t ) is a minimizer of Pmec, then its projection n(s) to the sphere S2 is a minimizer of Pcurve. Appendix 2: Integration of the Hamiltonian System Pm1ec To integrate the Hamiltonian system (39), we follow the idea of V. Jurdjevic [24], where one can employ left invariance and introduce initial rotation D0 of momentum space, such that the initial momentum transforms to h(0) = (0, 0, M ), solve the problem in this simple case (i.e., to find the trajectory R˜ (t ) that corresponds to h(0) = (0, 0, M )), and obtain general solution R(t ) by a backward transformation Proof of Theorem 4.1 Substitution of (36) to (31) gives Here P (t ) = −h1(t )ω2 + h2(t )ω1 + h3(t )ω3 is a momentum covector expressed in the basis ωi , dual to Ai , i.e., ωi , A j = δi j . The vertical part (64) has P (t ) = R−1(t ) P (0) R(t ) as a solution. Thus, the Hamiltonian system (64)—(65) preserves the norm of momentum covector h(t ) = (h2(t ), −h1(t ), h3(t )), i.e., the value M 2 = h21(t ) + h22(t ) + h23(t ), whose isosurface is a coadjoint orbit of SO(3). The Killing form allows to identify so(3) and so(3)∗ (see [19]), and by isomorphism (11), we write P(t ) = R−1(t ) P(0)R(t ) ∼ h(t ) = h(0) R(t ). Next we apply a left action of SO(3) on solution curve R(·) R˜ (·) = D0 R(·), where D0 is chosen such that the initial momentum matrix P(0) and initial covector h(0) = (−h2(0), h1(0), h3(0)) are reduced to a simple form P(0) = M A3 and h(0) = (0, 0, M ), as the nicest possible representant within the coadjoint orbit of SO(3). Next we represent the rotation matrix R˜ in form (6), i.e., R˜ = ey˜ A3 e−x˜ A2 eθ˜ A1 , where we parameterize the rotations π π via 3 angles x˜ ∈ [− 2 , 2 ], y˜ ∈ R/{2π Z}, θ˜ ∈ R/{2π Z}. Substitution of (68) in (66) gives the momentum matrix P(t) = M e−θ˜(t)A1 ex˜(t)A2 e#−y˜(t)A3 $A%3 ey˜(t)A&3 e−x˜(t)A2 eθ˜(t)A1 , and the equivalent relation for the momentum covector (h2, −h1, h3) = (0, 0, M ) ⎝ ⎛ cos x˜ 0 − sin x˜ ⎞ ⎛ 1 0 0 ⎞ 0 1 0 ⎠ ⎝ 0 cos θ˜ − sin θ˜ ⎠ . sin x˜ 0 cos x˜ 0 sin θ˜ cos θ˜ By multiplying the matrices in the right-hand side, we obtain From (69) we immediately have θ˜(t ) = arg (h3(t ) − ih1(t )) , π π and since x˜(t ) ∈ [− 2 , 2 ] ⇒ cos x˜(t ) ≥ 0, we can also express x˜(t ) = arg To obtain y˜(t ) first notice that due to left invariance, the matrix R˜ satisfies the same equation as R, i.e., R˜˙ = R˜ Ω. Thus x˜, y˜, θ˜ satisfy the same equations as x , y, θ , i.e., the horizontal part of (39). Thus we have y˙˜(t ) = − By virtue of (69), we can express (70) as Further, since M 2 = h12(t ) + h22(t ) + h32(t ) = const, we have t − Thus we obtained solution for x˜(t ), y˜(t ), θ˜(t ). To finish the proof, it only remains to obtain solution for x (t ), y(t ), θ (t ). Recall that the expression (46) for R(t ) follows from orthogonality of matrix D0, i.e., D0−1 = D0T , and parameterization of R(t ) by x (t ), y(t ), θ (t ) as R(t ) = ey(t)A3 e−x(t)A2 eθ(t)A1 . To find x (t ), y(t ), θ (t ) notice that action of R(t ) on e1 and R−1(t ) = RT (t ) on e3 gives (R11(t ), R21(t ), R31(t ))T (R31(t ), R32(t ), R33(t ))T = (cos x (t ) cos y(t ), cos x (t ) sin y(t ), sin x (t ))T , = (sin x (t ), cos x (t ) sin θ (t ), cos θ (t ) cos x (t ))T . π , π2 ], first two equations in (45) follow Since x (t ) ∈ [− 2 from (71) and the third equation in (45) follows from (72). Appendix 3: Proof of Theorem 5 Recall that we consider the case (see Remark 4) where the curve starts from e1 ∈ S2 and goes in the direction of upper half of the sphere with tangent vector e3. Therefore, we have h1(0) = u1(0) ≥ 0, h1(0) > 0. = 21 , we Via the Hamiltonian H (s) = h1(s) = 0 ⇔ h2(s) = ±1. In the linear case χ = 0, Eq. (73) reads as h02 + h3s = ±1, 0 that has the minimal positive root In the elliptic and the hyperbolic cases, we need to find a minimal positive root of In the elliptic case, we have χ = ia, where a = (0, 1), yielding solution = − ai log ⎝ (h03)2−a2 1−(h20)2 +ia ⎞ Denote κ = (h03)2 + 1 − (h02)2 χ 2. Notice that when (ha032)2 + (h02)2 < 1 ⇔ κ < 0 we have smax(h(0)) = ∞ and therefore we have cuspless trajectory of infinite SR length, recall (25). This fact is also clear from phase portrait of dynamic of vertical part (52), see Fig. 7 (bottom, left), where intersection of the integral curve with the red straight line indicates the moment, when cusp appears. We see that some geodesics have no cusps in their spherical projection up to infinity. Notice that since SO(3) is compact, it has bounded diameter (i.e., there exists D > 0, such that the SR distance between any two elements of SO(3) does not exceed D). Thus in contrast to SE(2) group, where every cuspless geodesic is optimal (see [31]), we observe that there exist nonoptimal cuspless geodesics in SO(3) group. In the hyperbolic case χ = a = ξ 2 − 1, ξ > 1, we have κ ≥ 0 and the minimal positive root of (74) reads as smHax (a, h02, h03) = a1 log Finally, using the parameter χ = formula (54) for all cases. a2 1−(h20)2 +(h03)2+a Appendix 4: Outline of the Proof of Theorem 7 By the result in [35, Thm 11.15], points g ∈ SO(3) where 0 the SR distance dSO(3)(e, ·) given by (56) is nonsmooth are either conjugate points, Maxwell points, or abnormal points. So by assumption (and absence of abnormal geodesics due to a 2-bracket generating distribution Δ, cf. [37, ch:20.5.1]), 0 we see that dSO(3)(e, ·) is smooth at {γ0∗(τ ) | 0 < τ ≤ 1}. As a result, the mapping is smooth as it is a smooth composition of maps. Similarly for the Riemannian case ε > 0, recall (59), we have by the assumptions that is smooth for all ε > 0. Let G−1 denote the inverse metric tensor associated with G given by (57), and let Gε−1 denote the inverse metric tensor associated with Gε given by (60). Now we rely on standard results [7,12] on backtracking of optimal Riemannian geodesics in Riemannian manifolds. This means that we find the smooth optimal geodesics via an intrinsic gradient descent on the Riemannian distance map Wε(g) := dSεO(3)(g, e). This yields the identity where Wε is the unique viscosity solution [7,12] of the corresponding Riemannian eikonal equation given by Gε|g ( Gε|g)−1dWε(g), ( Gε|g)−1dWε(g) = 1, for g = e, and Wε(e) = 0. Next we transfer these standard results toward the subRiemannian case, by a limiting procedure. Firstly, we note that following a fully tangential approach to [36, App.A,Thm.2] yields (see Remark 18 below) the following respectively pointwise and uniform convergence: lεi↓m0 Wε(g) = W0(g), and γε∗ → γ0∗ , as ε ↓ 0. Secondly, one has for any smooth function f : SO(3) → R that the intrinsic gradients converge: = ξ −2 X1|g f X1|g + X2|g f X2|g , for all g ∈ SO(3) . (78) Thirdly, as γε∗ and γ0∗ are solutions to the Hamiltonian system of the Pontryagin Maximum Principle, the trajectories are continuously depending on ε > 0 and so are their derivatives. As a result, we have (75) γ˙0∗(τ ) = lεi↓m0 γ˙ε∗(τ ) = lεi↓m0 Wε(g) (Gε−1dWε)(γε∗(τ )) Furthermore, from (76) it follows that Gε|g Gε|−g1 dWε(g), Gε|−g1 dWε(g) G|g G|−g1 dW(g), G|−g1 dW(g) where we note that the third term under the square root, after substitution of (78), vanishes. Finally, the viscosity solution property [7,12] is also naturally carried over due to continuous ε dependence of the Hamiltonian and the convergences (77), (79). Remark 18 The idea for the convergence (77) which is tangential to the SE(2) = R2 S1 case considered in [36, App.A,Thm.2] is based on a similar convergence result by Jean-Marie Mirebeau & Chen Da of a related elastica model [48, app.A]. It relies on closeness of controllable paths [48, Cor.A.5] and Arzela-Ascoli’s theorem. Another ingredient is the continuity of the mapping of the pair (ε, g) onto the corresponding indicatrix Bε(g) := {v ∈ Tg(SO(3)) | G|g (v, v) ≤ 1}, i.e., the continuity of where the powerset P(Tg(SO(3))) of each tangent space Tg(SO(3)) ≡ R3 is equipped with the metric topology induced by the Hausdorff distance. The proof of this continuity is straightforward and is therefore omitted here. The SR-eikonal and Backtracking Equations for Cuspless SR Geodesics Now in this article, we are primarily interested in the SR geodesics whose spherical projection does not have cusps. This can be taken into account by modifying the standard SR-eikonal equation (62) as 2 = C (g), for g = e, (80) ⎩ W(e) = 0, and backtracking equation (61) as where u1(t ) = max(ξ02,CX21(|γγbb((tt))()W)) , u2(t ) = XC2|2γ(bγ(tb)((tW)) ) , and where the system is integrated for t ∈ [0, W(g1)]. The idea behind this is that for cuspless SR geodesics one has u1 positive which holds by if the corresponding momentum component h1 = X1(W) is positive. Note that in the eikonal equation one substitutes momentum h = dW into the Hamiltonian to achieve equidistant wavefront propagation [12]. For further details and analysis on the positive control restriction on a similar problem on SE(2), see [36]. Similar analysis benefits may be expected on the SO(3) case, but this is beyond the scope of this article. Appendix 5: Comparison of the SE(2) and SO(3) Geodesics in the Image Plane We organize the comparison in the following way: 1) Choose an initial condition V0 = (X0, Y0, Θ0) and a terminal condition V1 = (X1, Y1, Θ1) in D = [Xmin, Xmax ] × [Ymin, Ymax ] × [−π, π] ⊂ R2 × S1 ∼= SE(2), where Xmax = Ymax = −Xmin = −Ymin = (c + a) tan ψmax = 110459 tan π8 . These values are obtained from schematic eye model for η = 1, recall Sect. 1.1, where ψmax = π8 is a maximum scanning angle via a standard fundus cameras. 2) Compute the SR distance WSE(2)(V ) = dSE(2)(V0, V ) in the volume of interest D V via SR-FM-SE(2) (subRiemannian fast marching in SE(2), see [41]). 3) Find via backtracking in SE(2) (see [41]) a geodesic γ SE(2)(·) satisfying γ SE(2)(0) = V0 and γ SE(2)(T ) = V1, with T = WSE(2)(V1). 4) Find the initial condition ν0 = (x0, y0, θ0) and the termi nal condition ν1 = (x1, y1, θ1) in SO(3), obtained from V0 and V1 via planar projection Π −1, recall Eq. (5), as (xi , yi ) = Π −1(Xi , Yi ) =: (Φ1(X, Y ), Φ2(X, Y )) and ∂Φ1 ∂Φ1 x˙(Xi , Yi , Θi ):= ∂ X |(Xi ,Yi ) cos Θi + ∂Y |(Xi ,Yi ) sin Θi , ∂Φ2 ∂Φ2 y˙(Xi , Yi , Θi ):= ∂ X |(Xi ,Yi ) cos Θi + ∂Y |(Xi ,Yi ) sin Θi . 5) Compute the SR distance WSO(3)(ν) = dSO(3)(ν0, ν) in the domain ν ∈ D = [xmin , xmax ] × [ymin, ymax ] × [−π, π ] via SR-FM in SO(3). Here we set xmax = ymax = −xmin = −ymin = 0.63, recall Eq. (1). 6) Find via backtracking (61) a geodesic γ SO(3)(·) satisfying γ SO(3)(0) = ν0 and γ SO(3)(T ) = ν1. 7) Plot in the image plane both the spatial projection Γ SE(2)(·) = (X (·), Y (·)) of the geodesic γ SE(2)(·) and the planar projection Γ SO(3)(·) = Π (x (·), y(·)) of the geodesic γ SO(3)(·). R. Duits received his M.Sc. degree (cum laude) in Mathematics in 2001 at the TU/e, The Netherlands. He received his Ph.D. degree (cum laude) at the Department of Biomedical Engineering at the TU/e. Currently he is an associate professor at the Department of Applied Mathematics & Computer Science and affiliated at the Biomedical Image Analysis group at the Biomedical Engineering Department. He has received an ERCStG 2013 personal grant (Lie analysis, no. 335555). His research interests include group theory, functional analysis, differential geometry, PDEs, harmonic analysis, Yu. Sachkov received his Ph.D. degree in Mathematics from Lomonosov Moscow State University, Russia, in 1992 and Doctor of Science degree from Steklov Mathematical Institute, Moscow, Russia, in 2008. He is Chief of Control Processes Research Center at Program Systems Institute, Pereslavl-Zale ssky (Russian Academy of Sciences), and Professor at University of Pereslavl, Russia. His research interests include geometric control theory, subRiemannian geometry on Lie groups, and their applications to mechanics, robotics, and vision. E. J. Bekkers is as a Ph.D. student at the Medical Image Analysis Group at Eindhoven University of Technology (TU/e). He has a former education in Biomedical Engineering (BSc+MSc from TU/e) and Mechanical Engineering (P from TU/e). His PhD project is entitled “Sub-Riemannian Geometry in Retinal Image Analysis” and is part of the academic– industrial research grant NWO Hé-Programme of innovation (co-funded by i-Optics). Within the application of retinal image analysis, the main focus is on the modeling and detection of key anatomical landmarks and the extraction of biomarkers for diabetes research. From a methodological point of view, focus is on the development of image analysis methods based on subRiemannian geometry and Lie group analysis (where in particular the Euclidean motion group is considered). His research interests include medical image analysis, differential geometry, convex optimization, and biomarker research. I. Beschastnyi received his M.S. at Kaliningrad State Technical University in 2011. He was trained in graduate program in Program Systems Institute, Pereslavl-Zalessky (Russian Academy of Sciences) from 2011 till 2014. Currently he is a Ph.D. student at International School of Advanced Studies (SISSA). His research interests include sub-Riemannian geometry on Lie groups and secondorder optimality conditions. 1. Peyré , G. , Péchaud , M. , Keriven , R., Cohen , L.D.: Geodesic methods in computer vision and graphics . Found. Trends. Comp. Comput. Graph. Vis . 5 ( 34 ), 197 - 397 ( 2010 ) 2. Sethian , J.A.: Level Set Methods and Fast Marching Methods . Cambridge University Press , Cambridge ( 1999 ) 3. Caselles , V. , Kimmel , R. , Sapiro , G. : Geodesic active contours . Int. J. Comp. Vis . 22 ( 1 ), 61 - 79 ( 1997 ) 4. Osher , S. , Fedkiw , R.P.: Level Set Methods and Dynamic Implicit Surfaces . Applied Mathematical Science, Springer, New York ( 2003 ) 5. Mirebeau , J.-M.: Anisotropic fast-marching on cartesian grids using Lattice Basis Reduction . SIAM J. Num. Anal . 52 ( 4 ), 1573 ( 2014 ) 6. Ikram , M.K. , Ong , Y.T. , Cheung , C.Y. , Wong, T.Y. : Retinal vascular caliber measurements: clinical significance, current knowledge and future perspectives . Ophthalmologica 229 ( 3 ), 125 - 136 ( 2013 ) 7. Lions , P.L. : Generalized Solutions of Hamilton-Jacobi Equations . Pitman , p. 317 ( 1982 ) 8. Crandall , M.G. , Lions , P.L. : Viscosity solutions of HamiltonJacobi equations . Trans. Am. Math. Soc . 277 ( 1 ), 1 - 42 ( 1983 ) 9. Sasongko , M.B. , Wong, T.Y. , Nguyen, T.T. , Cheung , C.Y. , Shaw , J.E. , Wang , J.J. : Retinal vascular tortuosity in persons with diabetes and diabetic retinopathy . Diabetologia 54 , 2409 - 2416 ( 2011 ) 10. Kalitzeos , A.A. , Lip , G.Y. , Heitmar , R.: Retinal vessel tortuosity measures and their applications . Exp. Eye Res . 106 , 40 - 46 ( 2013 ) 11. Bekkers , E. , Zhang , J. , Duits , R., ter Haar Romeny , B. : Curvature based biomarkers for diabetic retinopathy via exponential curve fits in SE(2) . In: Trucco, E. , Chen , X. , Garvin , M.K. , Liu , J.J. , Frank , X.Y . (eds.) Proceedings of the Ophthalmic Medical Image Analysis Second International Workshop , OMIA 2015 , Held in Conjunction with MICCAI 2015 , Munchen, Germany, October 9 , 2015 , pp. 113 - 120 . Iowa Research Online ( 2015 ) 12. Evans , L.C. : Partial Differential Equations. Graduate Studies in Mathematics vol. 19 , AMS, Providence ( 1998 ) 13. Petitot , J. : The neurogeometry of pinwheels as a sub-Riemannian contact structure . J. Physiol . Paris 97, 265 - 309 ( 2003 ) 14. Citti , G. , Sarti , A. : A cortical based model of perceptual completion in the roto-translation space . J. Math. Imaging Vis . 24 , 307 - 326 ( 2006 ) 15. Boscain, U. , Rossi , F. : Invariant Carnot-Caratheodory metrics on S3 , S O ( 3 ), S L ( 2 ) and lens spaces , SIAM. J. Control Optim . 47 , 1851 - 1878 ( 2008 ) 16. Boscain, U. , Rossi , F. : Projective Reed-Shepp car on S2 with quadratic cost . ESAIM: COCV 16 , 275 - 297 ( 2010 ) 17. Boscain, U. , Gauthier, J.-P., Chertovskih , R. , Remizov , A. : Hypoelliptic diffusion and human vision: a semidiscrete new twist . SIAM J. Imaging Sci . 7 ( 2 ), 669 - 695 ( 2014 ) 18. Bekkers , E. , Duits , R. , Berendschot , T. , Romeny , B.H. : A multiorientation analysis approach to retinal vessel tracking . J. Math. Imaging Vis . ( 2014 ) 19. Beschastnyi , I., Sachkov , Y. : Sub-Riemannian and almostRiemannian geodesics on SO(3) and S2 , ArXiv: 1409 . 1559 , ( 2014 ) 20. Berestovskii , V. : Geodesics of a left-invariant nonholonomic Riemannian metric on the group of motions of the Euclidean plane . Sib. Math. J . 35 ( 6 ), 1083 - 1088 ( 1994 ) 21. Berestovskii , V. , Zubareva , I.: Sub-Riemannian distance in the Lie groups SU (2) and S O ( 3 ). Sib. Adv. Math. 26 ( 2 ), 77 - 89 ( 2016 ) 22. Bonnard , B. , Cots , O. , Pomet , J.-B., Shcherbakova , N. : Riemannian metrics on 2D-manifolds related to the Euler-Poinsot rigid body motion . ESAIM: Control Optim. Calc. Var . 20 , 864 - 893 ( 2014 ) 23. Bonnard , B. , Chyba M. : Two applications of geometric optimal control to the dynamics of spin particle , 2014 , Preprint HAL Id: hal- 00956828 24. Jurdjevic , V. : Geometric Control Theory, Cambridge Studies in Advanced Mathematics (No. 52) . Cambridge University Press, Cambridge ( 1996 ) 25. Kruglikov , B. , Vollmer , A. , Lukes-Gerakopoulos , G. : On integrability of certain rank 2 sub-Riemannian structures . ArXiv: 1507. 03082 ( 2015 ) 26. Mashtakov , A.P. , Ardentov , A.A. , Sachkov , Y.L. : Parallel algorithm and software for image inpainting via sub-Riemannian minimizers on the group of rototranslations . Numer. Math. Theory Methods Appl . 6 ( 1 ), 95 - 115 ( 2013 ) 27. Mashtakov , A.P. , Sachkov, YuL: Superintegrability of left-invariant sub-Riemannian structures on unimodular three-dimensional lie groups . Differ. Equ . 51 (11), 1476 - 1483 ( 2015 ) 28. Sachkov , Yu.L.: Cut locus and optimal synthesis in the subRiemannian problem on the group of motions of a plane . ESAIM: COCV , 17 , 293 - 321 ( 2011 ) 29. Moiseev , I., Sachkov , Yu.L.: Maxwell strata in sub-Riemannian problem on the group of motions of a plane . ESAIM: COCV 16 , 380 - 399 ( 2010 ) 30. Boscain, U. , Duits , R. , Rossi , F. , Sachkov , Y.L. : Curve cuspless reconstruction via sub-Riemannian geometry . ESAIM: Control Optim. Calc. Var . 20 , 748 - 770 ( 2014 ). doi:10.1051/cocv/2013082 31. Duits , R. , Boscain, U. , Rossi , F. , Sachkov , Y.L. : Association Fields via Cuspless Sub-Riemannian Geodesics in SE(2) . JMIV 49 ( 2 ), 384 - 417 ( 2014 ). doi:10.1007/s10851- 013 - 0475 -y, http:// bmia.bmt.tue.nl/people/RDuits/cusp 32. Duits , R. , Ghosh , A. , Dela Haije , T.C.J., Mashtakov , A. : On subRiemannian geodesics in S E ( 3 ) whose spatial projections do not have cusps . J. Dyn. Control Syst . 22 ( 4 ), 771 - 805 ( 2016 ) 33. Montgomery , R.: A Tour of Subriemannian Geometries, Their Geodesics and Applications. Mathematical Surveys and Monographs ( 2002 ) 34. Agrachev , A.A.: Exponential mappings for contact subRiemannian structures . J. Dyn. Control Syst . 2 ( 3 ), 321 - 358 ( 1996 ) 35. Agrachev , A.A. , Barilari , D. , Boscain, U.: Introduction to Riemannian and Sub-Riemannian Geometry from the Hamiltonian Viewpoint . Preprint SISSA 09/ 2012 /M ( 2016 ) 36. Duits , R. , Meesters , S. , Mirebeau , J.M. , Portegies , J. : Optimal paths for variants of the 2D and 3D Reeds-Shepp car with applications in image analysis . ArXiv:1612 . 06137 , ( 2016 ) 37. Agrachev , A.A. , Sachkov, YuL: Control Theory from the Geometric Viewpoint . Springer, New York ( 2004 ) 38. Duits , R. , Ghosh , A. , Dela Haije , T.C.J., Sachkov , Y.L. : Cuspless sub-Riemannian geodesics within the Euclidean motion group S E (d) . In: Neuromathematics of Vision, Springer Series Lecture Notes in Morphogenesis , vol 1 , pp. 173 - 240 ( 2014 ) 39. Bekkers , E.J. , Duits , R, Mashtakov, A. , Sanguinetti , G.R. : (joint main authors). Data-driven sub-Riemannian geodesics in SE(2) . In: Aujol, Nikolova, Papadakis (eds.) Proceedings of SSVM 2015 , LNCS 9087 , pp. 613 - 625 ( 2015 ) 40. Bekkers E.J. , Duits R. , Mashtakov A. , Sanguinetti , G.R. : (joint main authors). A PDE approach to data-driven sub-Riemannian geodesics in SE(2) . SIAM J. Imaging Sci . 8 ( 4 ), 2740 - 2770 ( 2015 ) 41. Sanguinetti , G. , Duits , R. , Bekkers , E. , Janssen , M.H.J. , Mashtakov , A. , Mirebeau , J.M. : Sub-Riemannian fast marching in SE(2) . In: Progress in Pattern Recognition, Image Analysis , Computer Vision, and Applications. Lecture Notes in Computer Science , vol. 9423 , pp. 366 - 374 ( 2015 ) 42. American Academy of Ophthalmology, Clinical Optics , 2013 - 2014 . Basic and Clinical Science Course. American Academy of Ophthalmology , San Francisco ( 2013 ) 43. Rouy , E. , Tourin , A. : A viscosity solutions approach to shape-fromshading . SIAM J. Num. Anal . 29 , 867 - 884 ( 1992 ) 44. Calin , O. , Chang , D.-C. , Markina , I.: Sub-Riemannian geometry on the sphere S3 . Can . J. Math. 61 ( 4 ), 721 - 739 ( 2009 ) 45. Chang , D.-C. , Markina , I. , Vasilev , A. : Sub-Riemannian geometry on the 3-D sphere: complex analysis and operator theory , 3 , pp. 361 - 377 ( 2009 ) 46. Field , D.J. , Hayes , A. , Hess , R.: Contour integration by the human visual system: evidence for a local “association field” . Vis. Res . 33 ( 2 ), 173 - 193 ( 1993 ) 47. Chen , D. , Mirebeau , J.-M. , Cohen , L.D.: Global minimum for curvature penalized minimal path method . In: Xie, X. , Jones , M.W. , Tam , G.K.L. (eds.) Proceedings of the British Machine Vision Conference (BMVC) , pp. 86 . 1 - 86 .12. BMVA Press, Guildford ( 2015 ) 48. Chen , D. : New Minimal Path Models for Tubular Structure Extraction and Image Segmentation . PHD thesis , École Doctorale de Dauphine Paris, ( 2016 ) 49. Rifford , L.: Sub-Riemannian geometry and optimal transport . Springer Briefs in Mathematics. Springer International Publishing, New York ( 2014 ) 50. Frangi , A. , Niessen , W. , Vincken , K. , Viergever , M. : Multiscale vessel enhancement filtering . In Proc. of the MICCAI , Lecture Notes in Computer Science, Cambridge, pp. 130 - 137 , ( 1998 ) 51. ter Haar Romeny , B.M.: Front-End Vision and Multi-Scale Image Analysis . Springer, 484 pp, ( 2004 )


This is a preview of a remote PDF: https://link.springer.com/content/pdf/10.1007%2Fs10851-017-0705-9.pdf

A. Mashtakov, R. Duits, Yu. Sachkov, E. J. Bekkers, I. Beschastnyi. Tracking of Lines in Spherical Images via Sub-Riemannian Geodesics in \({\text {SO(3)}}\), Journal of Mathematical Imaging and Vision, 2017, 239-264, DOI: 10.1007/s10851-017-0705-9