Tracking of Lines in Spherical Images via SubRiemannian Geodesics in \({\text {SO(3)}}\)
Tracking of Lines in Spherical Images via SubRiemannian 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 , PereslavlZalessky , 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 subRiemannian arclength. For case C = 1 (datadriven 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.
SubRiemannian 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
subRiemannian (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 wellknown 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 datadriven 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, subRiemannian (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
datadriven geodesic methods has been presented in [40].
There, a computational framework for tracking of lines via
globally optimal datadriven subRiemannian 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 subRiemannian manifold structure in a
different Lie group, namely the group SO(3) (consisting of 3D
rotations) acting transitively on the 2sphere 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 curvestiffness 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 datadriven 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 rototranslation 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) corticalbased 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 (datadriven)
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 HJBPDE 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 leftinvariant 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 leftinvariant 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 ChowRashevskii and Filippov theorems on
subRiemannian manifolds [37]. Moreover, the geodesics of
Pmec are smooth, since there are no abnormal extremals (see
Example 1.3.13 in [49]).
Remark 8 SubRiemannian 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 subRiemannian 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 subRiemannian 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
ReedsShepp car on a sphere. The ReedsShepp 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 SubRiemannian 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 leftinvariant 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 firstorder necessary condition for
optimality [37], and we derive the explicit formulas for the
subRiemannian geodesics. Afterward we explain the notion
of subRiemannian 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
oneparameter 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
leftinvariant vector fields in SO(3) expressed in our chart.
We denote them by X1, X2, and the third basis leftinvariant
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 leftinvariant 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
PDEbased approach for solving the boundary value problem
(BVP) for datadriven SR geodesics, where we reformulate
the problem as a solution to the SReikonal equation.
3.5 SubRiemannian Wavefronts and Spheres
A useful tool for studying geodesics in leftinvariant 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 subRiemannian
sphere, which is a set of endpoints equidistant from e:
4 SubRiemannian 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 subRiemannian 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 socalled
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 biinvariant 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 DataDriven SR Geodesics in
SO(3)
In this section, we adapt the PDE approach for datadriven
SR geodesics in SE(2) [39,40] to the SO(3) group. Here we
consider the basis leftinvariant 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 subRiemannian metric tensor
G = C 2(·) ξ 2ω1 ⊗ ω1 + ω2 ⊗ ω2
5.1 SubRiemannian Fast Marching in SO(3)
Here we propose a SRFM method (subRiemannian fast
marching) for the computation of datadriven SR length
minimizers (not necessarily cuspless) in SO(3) group, as a
solution to the SReikonal system (62). This method was
successfully used in [41] for the computation of datadriven SR
length minimizers in SE(2) group. The method is based on
a Riemannian approximation of subRiemannian 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 leftinvariant one forms ωi ,
satisfying ωi , X j = δij , are expressed as
We approximate the subRiemannian manifold by a
Riemannian manifold by fixing a small ε > 0. Moreover, the
SReikonal 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
subRiemannian case. See details in Appendix 4.
The following theorem summarizes our approach for the
computation of datadriven subRiemannian 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 dataadaptive nonuniform cost) to tracking
of blood vessels in retinal images.
6 Experiments
In Sect. 6.1, we verify the SRFM method by comparison of
SR length minimizers/spheres obtained via SRFM 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 (SRFM) 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 SRFM 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
thirdorder 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
SRFM 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 SRFM 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 SRFM method we also perform experiments
with comparison of the geodesics obtained via SRFM 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 SRFM, but now
for the geodesics whose spherical projections have cusps.
To conclude this section, we provide one more experiment,
where we compare the exact subRiemannian wavefronts and
SR spheres computed via SRFM. 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 SRFM (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 SRFM (Color figure
online)
isosurface S F M ( 3125 π ) = {g ∈ SO(3) W (g) = 3125 π }. Here
we see that the isosurface computed via SRFM 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 selfintersections 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
SRFM 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 SRFM vessel
tracking method along datadriven 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 datadriven 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 selfintersections 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 datadriven 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 leftinvariant
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 (rightbottom
figure) are depicted with slight shift. The upper curve is a spatial
projection Γ SE(2) of the datadriven 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 subRiemannian 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 SRFM) 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
Datadriven subRiemannian 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 datadriven 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 datadriven 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 datadriven 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 datadriven SR geodesic model in SO(3)
is a relevant extension of our datadriven 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 datadriven
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/20072013)/ERC grant
Lie Analysis, a.n. 335555. The authors gratefully acknowledge
EUMarie Curie project MANET ag. no. 607643, and RFBR research Project
No. 163160083 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 righthand 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 2bracket 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 X1g f X1g + X2g f X2g , 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)
Gg 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 JeanMarie Mirebeau & Chen Da of a related elastica
model [48, app.A]. It relies on closeness of controllable
paths [48, Cor.A.5] and ArzelaAscoli’s theorem. Another
ingredient is the continuity of the mapping of the pair
(ε, g) onto the corresponding indicatrix Bε(g) := {v ∈
Tg(SO(3))  Gg (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 SReikonal 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
SReikonal 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 ) = XC22γ(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 SRFMSE(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 SRFM 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, PereslavlZale
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 “SubRiemannian
Geometry in Retinal Image Analysis”
and is part of the academic–
industrial research grant NWO
HéProgramme of innovation
(cofunded by iOptics). 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, PereslavlZalessky
(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 subRiemannian
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 fastmarching 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 HamiltonJacobi 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 subRiemannian contact structure . J. Physiol . Paris 97, 265  309 ( 2003 )
14. Citti , G. , Sarti , A. : A cortical based model of perceptual completion in the rototranslation space . J. Math. Imaging Vis . 24 , 307  326 ( 2006 )
15. Boscain, U. , Rossi , F. : Invariant CarnotCaratheodory 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 ReedShepp 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. : SubRiemannian and almostRiemannian geodesics on SO(3) and S2 , ArXiv: 1409 . 1559 , ( 2014 )
20. Berestovskii , V. : Geodesics of a leftinvariant 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.: SubRiemannian 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 2Dmanifolds related to the EulerPoinsot 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. , LukesGerakopoulos , G. : On integrability of certain rank 2 subRiemannian structures . ArXiv: 1507. 03082 ( 2015 )
26. Mashtakov , A.P. , Ardentov , A.A. , Sachkov , Y.L. : Parallel algorithm and software for image inpainting via subRiemannian minimizers on the group of rototranslations . Numer. Math. Theory Methods Appl . 6 ( 1 ), 95  115 ( 2013 )
27. Mashtakov , A.P. , Sachkov, YuL: Superintegrability of leftinvariant subRiemannian structures on unimodular threedimensional 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 subRiemannian 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 subRiemannian 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 SubRiemannian 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 SubRiemannian 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 ReedsShepp 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 subRiemannian 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). Datadriven subRiemannian 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 datadriven subRiemannian 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. : SubRiemannian 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 shapefromshading . SIAM J. Num. Anal . 29 , 867  884 ( 1992 )
44. Calin , O. , Chang , D.C. , Markina , I.: SubRiemannian geometry on the sphere S3 . Can . J. Math. 61 ( 4 ), 721  739 ( 2009 )
45. Chang , D.C. , Markina , I. , Vasilev , A. : SubRiemannian geometry on the 3D 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.: SubRiemannian 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.: FrontEnd Vision and MultiScale Image Analysis . Springer, 484 pp, ( 2004 )