Tracking of Lines in Spherical Images via Sub-Riemannian Geodesics in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\text {SO(3)}}$$\end{document}SO(3)

In order to detect salient lines in spherical images, we consider the problem of minimizing the functional \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\int \limits _0^l \mathfrak {C}(\gamma (s)) \sqrt{\xi ^2 + k_g^2(s)} \, \mathrm{d}s$$\end{document}∫0lC(γ(s))ξ2+kg2(s)ds for a curve \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document}γ on a sphere with fixed boundary points and directions. The total length l is free, s denotes the spherical arclength, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k_g$$\end{document}kg denotes the geodesic curvature of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document}γ. Here the smooth external cost \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathfrak {C}\ge \delta >0$$\end{document}C≥δ>0 is obtained from spherical data. We lift this problem to the sub-Riemannian (SR) problem in Lie group \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\text {SO(3)}}$$\end{document}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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi > 0, \mathfrak {C} \ne 1$$\end{document}ξ>0,C≠1. For \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathfrak {C}=1$$\end{document}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 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathfrak {C} \ne 1$$\end{document}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.


Introduction
In computer vision, it is common to extract salient curves in flat images via data-driven minimal paths or geodesics [1][2][3][4][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]. 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 S 2 .
Here we study the problem P curve (S 2 ) of finding a smooth curve n(·) on a unit sphere S 2 that satisfies given boundary conditions (both positions and velocities) n(0) = n 0 , n(l) = n 1 , n (0) = n 0 , n (l) = n 1 and minimizes the functional l 0 C(n(s)) ξ 2 + k 2 g (s) ds, where k g (·) 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 : S 2 → R + adding for adaptation to a given spherical image data. We state this problem P curve (S 2 ) more explicitly in Sect. 2.2. The problem P curve (S 2 ) is a spherical analog of a wellknown problem P curve (R 2 ) [30,31] (cf. Fig. 2) of finding a smooth curve x(·) on a plane R 2 that satisfies given boundary conditions and minimizes the functional l 0 c(x(s)) ξ 2 + k 2 (s) ds, Fig. 2 Left: problem P curve (S 2 ): 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 P curve (R 2 ) [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 : R 2 → 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,ȳ) 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 P curve (S 2 ) can be considered as a spherical extension of a (flat) cortical-based model P curve (R 2 ) 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][20][21][22][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 P curve on the set S 2 to a SR problem on the group SO(3), which we will call P mec , 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 P curve (S 2 ) for the uniform cost case C = 1. • Analysis and parametrization of cusp points arising in spherical projections of geodesics of P mec (SO (3)). • A new HJB-PDE theory for the numeric computation of globally optimal geodesics of P curve (S 2 ) and P mec (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.

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). 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 a max = 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 a max = 6.5 mm 10.5 mm · R = 13 21 in dimensionless coordinates.
In order to compute the maximum absolute angleȳ max , let us express the angle |ȳ| with respect to center of the eyeball (see right Fig. 3). Expressing the squared length of segment N Q yields (a + R cos |ȳ|) 2 cos −2 ψ = (a + R cos |ȳ|) 2 + R 2 sin 2 |ȳ|.
Solving this equation with respect to cos |ȳ|, we obtain unique nonnegative solution: A standard fundus camera used for producing of the retinal images has the angular range ψ ∈ [− π 8 , π 8 ]. Thus, substitution R = 1, a = 13 21 , and ψ = π 8 in Eq. (1) gives the maximum anglē We rely on the following parameterization of the image sphere S 2 and the retinal sphereS 2 (see Fig. 4): and y,ȳ ∈ [−π, π]. Next we present the explicit relation between the object coordinates (x,ȳ)-spherical coordinates on a unit sphereS 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 S 2 is obtained by reflection of the physical retinal sphereS 2 through the point (−a, 0, 0) T ∈ R 3 . In this case, we have 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 wherē In these formulas, we need to substitute a = a max = 13 21 and c = 4 5 < R = 1, depicted in Fig. 4, where we work in dimensionless coordinates.

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 In mathematical analysis, we can set η = 1; however, in experiments η is to be taken into consideration. Note that for η = 1 we have 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)| |y| , and it has a maximum when y = y max . We have 0 ≤ GD(y) ≤ GD(y max ) ≈ 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 P curve (S 2 ) or P curve (R 2 ) in the retinal imaging application. In the sequel, we will write P curve instead of P curve (S 2 ) as we will always be concerned with the case where the base manifold equals S 2 .

Problem P curve on S and P mec 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 P curve on S 2 , which we formulate next.

Mathematical Foundation and Notations
The Lie group SO (3) is the group of all rotations about the origin in R 3 . We shall denote a counterclockwise rotation around axis a ∈ S 2 with angle ϕ via R a,ϕ , in particular for rotations around standard axes where we denote cx = cos x, cy = cos y, cθ = cos θ, sx = sin x, sy = sin y, s θ = sin θ , and where The Lie group SO(3) defines an associated Lie algebra where T Id (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 A i follows by 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) We express L in matrix form as Remark 2 Note that at the unity one has The formulas for X i in (10) follow by (cf. Remark 1) We also use the isomorphism between so(3) and R 3 where A i ∈ so(3), R ∈ SO(3), e i ∈ R 3 , i = 1, 2, 3. Note that (6) is a product of matrix exponentials: We choose to rely on this parameterization to keep the analogy with previous SE(2) models [28,31,40] mentioned in introduction.

Statement of the Problem P curve
Let S 2 = {n ∈ R 3 n = 1} be a sphere of unit radius. We consider the problem P curve (see Fig. 2), which is for given boundary points n 0 , n 1 ∈ S 2 and directions n 0 ∈ T n 0 (S 2 ), n 1 ∈ T n 1 (S 2 ), n 0 = n 1 = 1 to find a smooth curve n(·) : [0, l] → S 2 that satisfies the boundary conditions n(0) = n 0 , n(l) = n 1 , n (0) = n 0 , n (l) = n 1 , and for given ξ > 0 minimizes the functional L(n(·)) := l 0 C(n(s)) ξ 2 + k 2 g (s) ds, where k g (s) denotes the geodesic curvature on S 2 of n(·) evaluated in time s, and C : S 2 → [δ, +∞), δ > 0, is a smooth function that we call "external cost." Here the total length l is free and s = s 0 1 dσ = s 0 n (σ ) dσ denotes the spherical arclength. Thus, we have n (s) = 1 and the Gauss-Bonnet formula 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 since the problem is left-invariant w.r.t. the natural left action of SO(3) onto S 2 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 P curve .

Remark 5
Following the previous works [30,31] in SE(2) group, we do not expect existence of minimizers for the problem P curve (S 2 ) in general.

Statement of P mec Problem in SO(3)
We call P mec the following SR problem in SO (3): 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 S 2 to the group SO(3), i.e., C(R) = C(R e 1 ).
We study the problem P mec for C = 1 (case of uniform external cost) in Sect. 3, but let us first consider some preliminaries.

Remark 6
In P mec , we only have two velocity controls u 1 and u 2 in a 3D tangent space T R(t) (SO(3)) (cf. Fig. 5).

Remark 7
The existence of minimizers for the problem P mec is guaranteed by Chow-Rashevskii and Filippov theorems on sub-Riemannian manifolds [37]. Moreover, the geodesics of P mec are smooth, since there are no abnormal extremals (see Example 1.3.13 in [49]).
where the controls u 1 and u 2 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) = e 1 and n (0) = e 3 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 dt → min, (22) with T fixed (see, e.g., [33]).

Remark 10
In analogy with the SR problem P mec in SE(2), cf. [30,31], we sometimes call X 1 = −R A 2 the "spatial generator" and X 2 = R A 1 the "angular generator", despite the fact that X 1 and X 2 are both angular generators on S 2 . The problem P mec (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 u 1 controls the motion along X 1 , and the input u 2 controls the motion along X 2 , see Fig. 5.

Relation Between the Problems P curve and P mec
We call a spherical projection the following projection map from SO(3) onto S 2 (see Fig. 6): In coordinates (x, y, θ) defined by (6), we have So we see that (x, y) are spherical coordinates on S 2 . The controls u 1 and u 2 along the "spatial generator" X 1 and the "angular generator" X 2 (cf. Remark 10)

Fig. 6
Left: illustration of the parameterizations used in P mec and P curve . 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 S 2 In problem P curve , one is interested in a curve n(s) = R(t (s)) e 1 , which satisfies (13) and minimizes (14). Here Next we show that the spherical projection (23) of certain minimizers of P mec provides solution of problem P curve . More precisely, this only holds for the minimizers whose spherical projection does not have a cusp.

Definition 1 The spherical projection of a minimizer of
That is, if the control in front of the "spatial" generator X 1 switches sign locally. We are interested in the first cusp time t cusp = min n∈N {t n cusp > 0}, and we call s max the corresponding value of spherical arclength, s.t. t cusp = t (s max ) via (25).
Notice that if u 1 ≡ 0 then the trajectory of P mec is projected at a single point on S 2 which does not provide a solution to P curve . This allows us to define s max as so that t cusp = t (s max ).
The following theorem states that minimizers n(s) of P curve for s ∈ [0, s max ) are given by spherical projection of the minimizers R(t) of P mec for t ∈ [0, t cusp ). Here s denotes the spherical arclength parameter, defined by d ds n(s) ≡ 1, and t denotes the SR arclength, defined by C(R(t)) ξ 2 u 2 be a minimizer of P mec parameterized by SR arclength, and let the corresponding Then for such boundary conditions, P curve has a minimizer n(s), along which we have for 0 ≤ s ≤ l < s max , and T = t (l).
Proof See Appendix 1.

The Sub-Riemannian Problem P 1 mec in SO(3)
In this section, we study the problem P mec in the special case of uniform external cost C = 1, which we call P 1 mec . 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][21][22][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 P mec , 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 P 1 mec 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.

P 1 mec as Optimal Control Problem on S 2 × S 1
In this section, we introduce the chart in SO(3) and formulate problem P 1 mec given by (17)- (20), see Sect. 2.3, as an optimal control problem on Z} are the spherical coordinates, 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 × S 1 θ ) 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 Therefore, the control system (17) can be rewritten aṡ Vector fields near the controls u i are two of three of the basis left-invariant vector fields in SO(3) expressed in our chart. We denote them by X 1 , X 2 , and the third basis left-invariant vector field is given by their commutator Thus we have Then problem P 1 mec given by (17)- (20) for C = 1, taking into account Remark 9, is equivalent to the following optimal control problem: Remark 11 A projective version of problem P 1 mec 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 = (x 1 , y 1 , θ 1 ) andh = (x 1 , y 1 , θ 1 + π). The associated SR distance is given by The model in [16] is connected to problem P 1 mec via right The existence of optimal trajectories whose spherical projection contains a cusp is mentioned in [16].

Application of the Pontryagin Maximum Principle
In this subsection, we apply the Pontryagin maximum principle (PMP) to problem P 1 mec 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 Here ν = (x, y, θ), X i are given by (32), ·, · denotes the action of a covector on a vector, and (d ν 1 The maximization condition of PMP reads as where u(t) denotes an extremal control and (λ(t), γ (t)) is an extremal.
The maximization condition gives the expression for the extremal controls Then the maximized Hamiltonian, which we call "the Hamiltonian", reads as The vertical part (for momentum components h i ) of the Hamiltonian system of PMP is given bẏ where {·, ·} denotes the Poisson bracket (see [37]). Using the standard relation between Poisson and Lie brackets Thus the Hamiltonian system with the Hamiltonian H follows by the last three identities and (38): with the boundary conditions 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].

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 h i (see illustration in Fig. 7) Different values of the Hamiltonian H correspond to different speed along extremal trajectories γ (·). By fixing the value  (41), we use SR arclength parameter can be seen as the curve formed by intersection of the cylinder H = 1 2 and the sphere M = const (see the yellow line in Fig. 7).
Elimination of h 1 from (41) and (42) yields Since h 2 2 ≤ 1, there exist the following types of solution of the vertical part of (39): 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.

Geodesics in SR Arclength Parametrization
In this subsection, we provide explicit formulas for the SR geodesics, where we reexamine results from [19,22].
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
where (β(t), c(t)) is a trajectory of (44), and Proof See Appendix 2. Note that (β(t), c(t)), andỹ(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.

Sub-Riemannian Wavefronts and Spheres
A useful tool for studying geodesics in left-invariant optimal control problems is the exponential map 2 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)]: Now we explain the notion of wavefront. By definition, the wavefront consists of endpoints of all the geodesics of the same length T : The outer surface of a wavefront forms a sub-Riemannian sphere, which is a set of endpoints equidistant from e: (3)), where d(q 1 , q 2 ) denotes the sub-Riemannian distance between q 1 and q 2 and t cut denotes the cut time, i.e., an instance of time, when a geodesic loses its optimality.

What Is Known About 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].

Remark 12
There are two reasons for a geodesic to lose optimality: 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 P 1 mec , 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.

Sub-Riemannian Geodesics in SO(3) with Cuspless Spherical Projection
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 X 1 (i.e., u 1 (t) > 0) before the first cusp time t < t cusp . Thus, by virtue of (36) the cuspless constraint is given by We shall often rely on short notation h i (s) := h i (t (s)), where we stress our notations in order to avoid confusion with the chain law for differentiation. Recall (26) and note that s max for a given geodesics is determined by its initial momentum. We write s max (h(0)) when we want to stress this dependence. We derive an explicit formula for s max (h(0)) in Sect. 4.1.

Theorem 3
For any s ∈ [0, s max (h(0))), h 1 (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 Here we recall that u 1 = ds dt , see (28), and h 1 = ξ 2 u 1 , see (36). Thus we obtained the vertical part of (49). Doing the same for the horizontal part in (39), we obtain the full system 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): (49) for all ξ = 1 is given by

Theorem 4 The general solution of the vertical part in
and for the case ξ = 1, we find straight lines parallel to h 2 axis in the (h 2 , h 3 )-phase portrait: Proof The solution to the system of ODEs is unique and it is readily checked (52) is a solution of the vertical part of (49) for all ξ = 1. For ξ < 1, we take main values for the complex square root and complex sinh and cosh. For both ξ = 1 and ξ ↓ 1 and ξ ↑ 1, we have h 2 (s) = h 0 2 + h 0 3 s and h 3 (s) = h 0 3 .

Computation of the First Cusp Time
To analyze the cusp points in problem P curve , we need to determine s max (h(0)). It is given, recall (26) and (36), by the minimal positive root of equation h 1 (s) = 0: (h(0)), recall (25), where
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].

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, s max (h(0))], where s max (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
where , and Proof It follows from Theoremes 2, 3, and 4. Momentum component h 1 is expressed in h 2 via the Hamiltonian (37) which equals to 1 2 along SR geodesics. The sign of h 1 (0) is equal to the sign of u 1 (0) = ds dt (0) and therefore nonnegative. Regarding the integration constraints, we note that M 2 = 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) where F denotes an elliptic integral of the first kind and Π denotes an elliptic integral of the third kind.

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 X i as differential operators of the first order, and we write X i (W) for the derivative of a function W : SO(3) → R along X i . We aim to solve the following optimal control problem: Here the terminal time T is free; and C : SO(3) → [δ, +∞), δ > 0 is the external cost. By rescaling of time τ = t T ∈ [0, 1] simultaneously with controlsū i (τ ) = T u i (T τ ), we write down the explicit solutions as Here the sub-Riemannian metric tensor is defined only on the distribution Δ, recall Remark 8.

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 fields 3 X i in SO(3) are given by the following differential operators: and corresponding basis left-invariant one forms ω i , satisfying ω i , X j = δ j i , 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: with Riemannian metric tensor G ε given by

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 γ
is not a conjugate point for all τ ∈ [0, 1] and all ε ≥ 0 sufficiently small.
with (ū 1 (τ ),ū 2 (τ )) = W(g) and with W(g) denoting the viscosity solution of the following boundary value problem: W(e) = 0. (62) For an outline of the proof, see Appendix 4, where we rely on a similar approximation approach as in [48, ch. 5 Here the diagonal matrix in the middle encodes the anisotropy between the X i 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.

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).

Verification of the Fast Marching Method in the Case of Uniform External Cost
In We have performed a series of such experiments and always obtained similar results when the geodesics γ i (s) were optimal for s ∈ [0, s i end ]. It was also remarkable that SR-FM resulted into different curves (length minimizers) when the geodesic γ i (s) was not optimal for s ∈ [0, s i end ]. 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 SR-FM 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( 15 32 π) and 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 .

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. (3) and SE(2) for C = 1

Comparison of SR Wavefronts in SO
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 S 2 x,y × S 1 θ was chosen in the way to obtain the analogy with (X, Y, Θ) ∈ R 2 × S 1 . 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(h 0 , ·), recall (47), i.e., by the projection in SO(3) of the solution of the Hamiltonian system (39), with initial con- = (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: 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: 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 WF SE (2) are the SR wavefronts of SE(2) and the green surfaces WF SO(3) are the SR wavefronts of SO(3). We see a very similar shape of WF SO(3) and WF SE(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 points 4 of WF SO(3) are located near singular points of WF SE (2) . This leads us to a conjecture that the location of conjugate points (open problem) can be estimated by the location of conjugate points 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 WF SE (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].

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 O XY , 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.

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 where we use standard multiscale vesselness [50] with λ i (X, Y ) eigenvalues of the Gaussian Hessian of image F : R 2 → R (maximized over scales) ordered by |λ 1 (X, Y )| < |λ 2 (X, Y )|, with β = c = 0.3, and with unit step function U (λ) = 1, for λ ≥ 0, 0, otherwise. Here we note that the Gaussian Hessian is given by H(G s * F) and computed by Gaussian derivatives [51] at multiscales s = 1 2 σ 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 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 is implied by Eq. (28). The second equality follows from application of PMP to problem P mec , which gives u 1 = h 1 ξ 2 C 2 and u 2 = h 2 C 2 . 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 h 1 = X 1 (W), h 2 = X 2 (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 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(·), θ (·)). in object coordinates on S 2 rather than planar curvature κ SO(3) in photo coordinates on projection on R 2 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 κ SO(3) g is not very significant. Thus we see the importance of including correct spherical geometry in vessel tracking algorithm in retinal images.

Vessel Analysis via SO(3) SR Geometry and SO(3) Riemannian Geometry
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.

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) ≡ R 2 S 1 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 S 2 , cf. Fig. 6, instead of a frame bundle above R 2 . 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][10][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.
But also the dynamics coincides, in the following sense. If n(s), s ∈ [0,l], is a smooth curve on the sphere S 2 with the initial conditionsñ(0) = e 1 ,ñ (0) = e 3 and a geodesic curvaturek g (s), then it can be lifted to a curveR(s) in SO (3) that is a trajectory of control system (17) with the controls u 1 (s) ≡ 1,ũ 2 (s) =k g (s) and the initial conditionR(0) = Id. The curveR(s)e 1 on S 2 has the same initial conditions and geodesic curvature asñ(s), thusR(s)e 1 ≡ñ(s). Summing up, since R(t) is a minimizer of P mec , then its projection n(s) to the sphere S 2 is a minimizer of P curve .

Appendix 3: Proof of Theorem 5
Recall that we consider the case (see Remark 4) where the curve starts from e 1 ∈ S 2 and goes in the direction of upper half of the sphere with tangent vector e 3 . Therefore, we have h 1 (0) = u 1 (0) ≥ 0, h 1 (0) > 0. Denote κ = (h 0 3 ) 2 + 1 − (h 0 2 ) 2 χ 2 . Notice that when (h 0 3 ) 2 a 2 + (h 0 2 ) 2 < 1 ⇔ κ < 0 we have s max (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.

Remark 18
The idea for the convergence (77)  where the powerset P(T g (SO(3))) of each tangent space T g (SO(3)) ≡ R 3 is equipped with the metric topology induced by the Hausdorff distance. The proof of this continuity is straightforward and is therefore omitted here.
The idea behind this is that for cuspless SR geodesics one has u 1 positive which holds by if the corresponding momentum component h 1 = X 1 (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.