Noname manuscript No.
(will be inserted by the editor)

Nilpotent Approximations of Sub-Riemannian Distances for
Fast Perceptual Grouping of Blood Vessels in 2D and 3D

arXiv:1707.02811v2 [] 8 Nov 2017

Erik J. Bekkers · Da Chen · Jorg M. Portegies

Received: date / Accepted: date

Abstract We propose an efficient approach for the
grouping of local orientations (points on vessels) via
nilpotent approximations of sub-Riemannian distances
in the 2D and 3D roto-translation groups SE(2) and
SE(3). In our distance approximations we consider homogeneous norms on nilpotent groups that locally approximate SE(n), and which are obtained via the exponential and logarithmic map on SE(n). In a qualitative validation we show that the norms provide accurate approximations of the true sub-Riemannian distances, and we discuss their relations to the fundamental solution of the sub-Laplacian on SE(n). The quantitative experiments further confirm the accuracy of
the approximations. Quantitative results are obtained
by evaluating perceptual grouping performance of retinal blood vessels in 2D images and curves in challenging 3D synthetic volumes. The results show that 1)
sub-Riemannian geometry is essential in achieving top
performance and 2) that grouping via the fast analytic approximations performs almost equally, or better, than data-adaptive fast marching approaches on
Rn and SE(n).
Keywords Sub-Riemanian geometry · Rototranslation group · SE(2) · SE(3) · Nilpotent approximation · Geodesic vessel tracking · Perceptual
grouping
E.J. Bekkers · J.M. Portegies
Centre for Analysis, Scientific computing and Applications
(CASA), Eindhoven University of Technology, the Netherlands
E-mail: e.j.bekkers@tue.nl, j.m.portegies@tue.nl
D. Chen
University Paris Dauphine, PSL Research University
CNRS, UMR 7534, CEREMADE, 75016 Paris, France
E-mail: chenda@ceremade.dauphine.fr

1 Introduction
In this paper we derive analytic formulas for approximations of sub-Riemannian distances on the 2D and 3D
rotation translation groups, denoted respectively with
SE(2) and SE(3). Additionally, we extend the perceptual grouping algorithm [13] for clustering of local orientations (points on blood vessels). Here clustering is
based on alignment of local orientations, which is quantified using sub-Riemannian distances on SE(n), see
Fig. 1 for an illustration.

1.1 Nilpotent Approximation
The sub-Riemannian distances on SE(n) are approximated via norms on the vectors obtained from the
logarithmic map (from group elements to the Lie algebra). This approach is motivated by problems from subRiemannian geometry in nilpotent Lie groups, in which
such homogenous norms provide exact fundamental solutions to sub-Laplacians.
The vectors obtained by the logarithmic map, expressed in a left-invariant basis, are the so-called exponential coordinates of the first kind. For a nilpotent
group of step two, like the Heisenberg group, these coordinates define (together with a group product defined
via the Baker-Campbell-Hausdorf (BCH) formula) a
global isomorphism to the group. In our SE(n) setting
we have to truncate the commutator series in the BCH
formula due to non-vanishing (higher-order) commutators, yielding a corresponding Heisenberg type approximation which we denote with (SE(n))0 . The obtained
Taylor development of the group product and associated left-invariant vector fields gives rise to a local ap-

2

Fig. 1 The red and green arrow have equal spatial and angular distance to the origin (black arrow). In a flat geometry on
R2 × S 1 the distance between the red and green arrow and
the source would be equal, and the geodesics straight lines
(see dashed lines). In sub-Riemannian geometry on SE(2)
the green arrow has a shorter distance to the source. The left
image shows 2D projections of the sub-Riemannian geodesics
in solid black, the right image shows their paths in SE(2).

proximation of the (sub-Riemannian) flows on SE(2)
in the sense of Rothschild and Stein [50].
We then define a norm on (SE(n))0 based on the
Folland-Kaplan-Korányi gauge, which is known for its
relation to the fundamental solution of the sub-Laplacian
on the Heisenberg group [29, 33, 35]. We reason that the
Folland-Kaplan-Korányi provides an accurate approximation to the fundamental solution on SE(n) as well,
as it provides the exact fundamental solution on the
Heisenberg type approximation (SE(n))0 . As such, we
provide an approach to approximating the heat kernel and fundamental solution of the sub-Laplacian on
SE(n), as an alternative to the works [22, 47, 12].
The distance associated with the Folland-KaplanKorányi type norm on (SE(n))0 is locally equivalent to
the sub-Riemannian distance on SE(2), as was formally
proven in full generality in the seminal work by Nagel,
Stein and Wainger [43]. In this paper, we show by qualitative and quantitative comparison that the norm on
(SE(n))0 indeed provides a sharp local approximation
of the sub-Riemannian distances on SE(n).

1.2 Perceptual Grouping
The motivation for perceptual grouping of local orientations comes from problems in medical image analysis in which the topologically correct reconstruction of
vessel (and pulmonary) trees is of great importance in
biomarker research and surgery planning. Knowing the
correct connectivity in tree structures not only allows
for local biomarker analysis (e.g., studies on bifurcation
and crossing properties [37]), but also allows for higher
level biomarker research via statistics on tree structures
[28]. Topological knowledge of vessel trees is also essential in determining artery/vein classification problems [26, 25, 15]. Finally, in many medical applications

Erik J. Bekkers et al.

involving vessel analysis, including topological tree reconstruction, distances between local orientations play
a crucial role [16, 57, 39, 55, 1, 27]. The approximate subRiemannian distance in this paper is analytic, fast, and
easy to implement, and as such may be a useful tool for
algorithms that rely on local orientation analysis.
Sub-Riemannian models are shown to be effective in
both image processing and in neuropsychological models for line perception in the primary visual cortex [45,
12, 53, 51, 40, 20, 7, 3, 48, 27]. In this paper we indeed observe by quantitative validation of automatic connectivity analysis that sub-Riemannian distances are preferred over their (full) Riemannian counter parts.
The approach taken in this paper for doing connectivity analysis is based on the perceptual grouping algorithm proposed by Cohen [13]. This algorithm turns
a set of key points into a graph by iteratively adding
edges between nodes based on their geodesic distances
while putting constraints on the number of connections
per node. The input set of key points may be obtained
via key point tracking algorithms [5, 34, 10], as is done
also in this paper, see Fig. 2.
In [13] an isotropic metric was used to define the
geodesic distances. Later, the perceptual grouping algorithm was adapted for use with anisotropic Riemannian metrics by Bougleux et al. [8]. In recent work [11],
it was further extended for the grouping of n closed contours for an a-priori n. There, a (sub-)Finsler metric on
position orientation space was used, similar to the subRiemannian metric used in this paper. As in [8] and
[11], we use the main algorithm of [13] as a backbone,
but we change the metric used for perceptual grouping
and we impose an additional constraint to avoid closed
loops (which are physically not realistic in the vessel
networks of interest).
With quantitative experiments we show that perceptual grouping with sub-Riemannian distances on SE(n)
is preferred over the use of (full) Riemannian distances
on SE(n), which is in turn preferred over grouping with
distances on Rn . Furthermore, the analytic approximations allow for fast perceptual grouping with competitive performance compared to data-adaptive subRiemannian distances computed via fast marching.

1.3 Paper Outline
In Sec. 2 and Sec. 3 we derive approximations for subRiemannian distances in respectively SE(2) and SE(3).
There, for each Lie group we first provide the preliminaries, then define the sub-Riemannian distance, and
then describe the proposed approximations. In Sec. 4
the algorithms (perceptual grouping, fast marching and

Title Suppressed Due to Excessive Length

3

key point tracking) are described, including an overview
of the different distances used in this paper. In Sec. 5 we
then compare the performance of the perceptual grouping algorithm using different distances, first on R2 and
SE(2) in Subsec. 5.1, then on R3 and SE(3) in Subsec. 5.2. General conclusions are provided in Sec. 6.

at g = (0, 0, 0) to an element in the group SE(2) by
following an integral curve along the left-invariant vector field (Lg )∗ X. The logarithmic map Log : SE(2) →
se(2) defines the mapping from group element to tangent vector at g = (0, 0, 0).
The Lie bracket for vector fields is defined as follows

2 Sub-Riemannian Distance and its
Approximation in SE(2)

[X, Y ] : = lim

2.1 The Lie Group SE(2)

I.e., it describes the infinitesimal displacement by following a path moving forth and back in X and Y directions. The Lie bracket of two vectors defines a new
vector (the commutator) and the Lie bracket of two
vector fields defines a new vector field. The non-zero
commutators of se(2) are

γ(t) − e
,
with
t2
(4)
γ(t) = Exp(−tY ) Exp(−tX) Exp(tY ) Exp(tX).
t→0

2.1.1 SE(2)
In order to measure distances between local orientations
we will consider the Lie group SE(2) as our base manifold. The group SE(2) = R2 o SO(2) is the semi-direct
product of the group of planar translations R2 and rotations SO(2), and its group product and inverse are
respectively defined via:
g · g 0 = (x, Rθ ) · (x0 , Rθ0 ) = (Rθ x0 + x, Rθ+θ0 ),
−1
g −1 = (−R−1
θ x, Rθ ),

(1)

with group elements g, g 0 ∈ SE(2). The group acts on
the (coupled) space of positions and orientations R2 o
S 1 via
g · (x0 , θ0 ) = (Rθ x0 + x, θ + θ0 ).
Since (x, Rθ ) · (0, 0) = (x, θ), we can uniquely identify the roto-translation group SE(2) with the space of
positions and orientations R2 o S 1 .
2.1.2 The Lie Algebra, Exponential Map and
Commutators

[A1 , A2 ] = −[A2 , A1 ] =

A3 ,

[A1 , A3 ] = −[A3 , A1 ] = −A2 .

(5)

2.2 Sub-Riemannian Geometry in SE(2)
We consider a sub-Riemannian geometry on SE(2) by
measuring distances between two points in SE(2) via
the lengths of shortest horizontal paths. A horizontal
path is a curve γ : [t0 , t1 ] ⊂ R → SE(2) with tangent
vectors γ̇(t) ∈ ∆|γ(t) := span{ A1 |γ(t) , A2 |γ(t) }, where
∆ denotes the sub-Bundle of the full tangent bundle
T (SE(2)) := span{A1 , A2 , A3 }. Lengths of horizontal
curves with γ̇(t) = u1 (t) A1 |γ(t) +u2 (t) A2 |γ(t) are measured by the sub-Riemannian metric tensor1
G ξ,C

γ(t)

(γ̇(t), γ̇(t)) := C(γ(t))2 (|u1 (t)|2 + ξ|u2 (t)|2 ),
(6)

The Lie algebra associated with SE(2) is the real vector
space se(2) = span{A1 , A2 , A3 } together with a bilinear operator [·, ·] : se(2) × se(2) → se(2) called the Lie
bracket (which we define below in Eq. (4)). The generators of the Lie algebra are given by the differential
frame {∂θ , ∂x , ∂y }|(0,0,0) at the origin
A1 = ∂θ |(0,0,0) , A2 = ∂x |(0,0,0) , A3 = ∂y |(0,0,0) , (2)
which define corresponding left-invariant vector fields
A1 |g = (Lg )∗ A1 =

∂θ |g ,

A2 |g = (Lg )∗ A2 =

cos θ ∂x |g + sin θ ∂y |g ,

in which C : SE(2) → R+ is an external cost which
penalizes the curves to move through certain regions in
SE(2), ξ is a parameter which balances the penalty of
motion in the angular and spatial directions and has
dimensions [1/length], and u1 and u2 are the control
parameters of the curve γ.
The sub-Riemannian distances between two points
g1 , g2 ∈ SE(2) is then given by
Z 1 q

ξ,C
d0 (g1 , g2 ) := inf
G |γ(t) (γ̇(t), γ̇(t))dt ,
(7)
0

(3)

A3 |g = (Lg )∗ A3 = − sin θ ∂x |g + cos θ ∂y |g
via the push-forward of left-multiplication, denoted by
(Lg )∗ , and with g = (x, y, θ) ∈ SE(2).
The exponential map Exp : se(2) → SE(2) defines a
mapping from a vector X ∈ se(2) in the tangent space

1
Due to the fact the the metric tensor is degenerate in the
A3 direction (tangent vectors are always contained within ∆)
it is not possible to represent the metric tensor in a standard
form as an invertible symmetric 3 × 3 matrix. This is however
possible when including an additional term −2 ξ 2 |u3 |2 (t) after which the tensor becomes (anisotropic) Riemannian [12,
52]. This Riemannian approximation converges to the subRiemannian tensor when  → 0 [9, App. A] and [24, Thm. 2].

4

Erik J. Bekkers et al.

Fig. 2 The pipeline for grouping vessel segments consists of 2 steps. First, key points are generated (from a single source point)
using minimal path tracking with key points [5]. Second, the automatically generated key points, with estimated orientations,
are grouped based on an adaption of the perceptual grouping algorithm [13] with the use of sub-Riemannian distances on
SE(2). The result on the right is obtained with the nilpotent approximations of the sub-Riemannian distances in SE(2).

where the infimum is taken over Lipschitz continuous
curves γ ∈ Lip([0, T ], SE(2)) with γ(0) = g1 , γ(1) =
g2 , and γ̇(t) = u1 (t) A1 |γ(t) + u2 (t) A2 |γ(t) . Note that
due to the inclusion of an external cost function C the
distance d is not strictly left-invariant, however, when
substituting C by Cg := C(g −1 h) in (7) we do have
left-invariance (i.e., then d(g · g1 , g · g2 ) = d(g1 , g2 )).
2.3 A Nilpotent Approximation (SE(2))0 of SE(2)
2.3.1 A Local Approximation via the
Baker-Campbell-Hausdorff Formula
Consider the exponential map from Lie algebra se(2) to
the group SE(2)
(c1 , c2 , c3 ) 7→ (x, y, θ) = Exp(c1 A1 + c2 A2 + c3 A3 ), (8)
with {Ai }3i=1 the basis vectors of se(2) given in (2),
and with (c1 , c2 , c3 ) the canonical coordinates of the
first kind given by
1
θ(y + x cot θ2 )
if θ =
6 0
1
2
c = θ,
c = 2
,
if θ = 0
x1
(9)
θ(−x + y cot θ2 ) if θ 6= 0
c3 = 2
.
y
if θ = 0
P3
For two left-invariant vector fields X = i=1 xi Ai and
P3
Y = i=1 y i Ai the Baker-Campbell-Hausdorff (BCH)
formula (see e.g. [49]) gives:
Log(Exp(X) Exp(Y )) = X + Y + 12 [X, Y ]
1
+ 12
([X, [X, Y ]] + [Y, [Y, X]])
+O([·, [·, [·, ·]]]),
(10)
where O([·, [·, [·, ·]]]) denotes higher order nested brackets. Since the Lie algebra se(2) is not nilpotent it has
non-vanishing Lie brackets of order ≥ 2 (cf. the commutator relations in (5)) the BCH formula gives an infinite
series of nested Lie brackets.

Here, we approximate the BCH formula SE(2) as2
1
Log(Exp(X) Exp(Y )) ≈ X + Y + [X, Y ],
2

(11)

by omitting the Lie brackets of order 2 (once nested
brackets) and higher, as if our Lie algebra se(2) is nilpotent of step 2. Then, together with the commutator relations [Ai , Ai ] = 0, A3 = [A1 , A2 ], and again omitting Lie brackets of order 2 (i.e., setting [A1 , A3 ] =
[A1 , [A1 , A2 ]] = 0), the BCH formula defines a group
product on the vector space R3 of the canonical coordinates of the first kind via
(x1 , x2 , x3 ) · (y 1 , y 2 , y 3 )


1 1 2
2 1
1
1
2
2
3
3
= x + y , x + y , x + y + (x y − x y ) .
2
(12)
The new group product (12), where the elements are
expressed in coordinates of the first kind (cf. Eq. (8)),
gives rise to a nilpotent Heisenberg group. It is a local3 approximation
of the
P3
P3 true group product g1 · g2 =
Exp( i=1 xi Ai ) · Exp( i=1 y i Ai ) given by (1). We denote this group by (SE(2))0 = H(3), with H(3) the 3
dimensional (nilpotent) Heisenberg group. Note that if
(x1 , x2 , x3 ) and (y 1 , y 2 , y 3 ) were coordinates of the first
kind for a group with a step-2 nilpotent algebra, then
this new group would be globally isomorphic to that
group. The new group (SE(2))0 defines a homogeneous
Carnot group with respect to the dilations
δs (c) = (s c1 , s c2 , s2 c3 ).

(13)

2
Note that such approximations of the BCH formula were
already introduced in [43, Thm. 2.22] in the general setting
by Nagel, Stein and Wainger in 1985.
3
With g1 , g2 ∈ SE(2) chosen close enough such that higher
order terms in (10) can be neglected.

Title Suppressed Due to Excessive Length

5

Fig. 3 Distances on SE(2) for ξ = 1, C = 1. Top row: Level sets of the distance volumes on SE(2). Bottom row: Minimum
intensity projections of the distances to the plane R2 with level set contours. From left to right: The sub-Riemannian distance
d0 (e, ·), see Eq. (7); Homogenous norms k·kξ,ζ , see Eq. (15), of the nilpotent approximation (SE(2))0 for respectively ζ = 44,
ζ = 16 (Folland-Kaplan-Korányi gauge) and c = ζ; The (ξ-isotropic) Riemannian distance d1 (e, ·) on SE(2), see Table 1 for
an overview of the different distances.

2.3.2 Homogeneous Norms on (SE(2))0 and the
Fundamental Solution of the sub-Laplacian
In our approximation of the sub-Riemannian distance
d0 of Eq. (7) we make use of the following homogenous
norm on (SE(2))0 :
p
kckζ := 4 (|c1 |2 + |c2 |2 )2 + ζ |c3 |2 ,
(14)
with constant ζ > 0 a relative penalty for the nonhorizontal part. For ζ = 16 this norm coincides with
the well-known Folland-Kaplan-Korányi gauge, which
is a widely studied norm on Carnot groups due to its
relation to fundamental solutions of sub-Laplacians [6]:
Folland found that kck2−Q
16 , with homogeneous dimensions Q, is (a constant multiple of) the fundamental solution of the canonical sub-Laplacian on the Heisenberg
group [29]; Kaplan showed that this relation more generally holds for H-type (Carnot) groups [33]; Korányi
derived many more of its properties in relation to harmonic analysis and potential theory [35].
In relation to sub-Riemannian geometry on SE(2)
and its sub-Laplacian L := A21 + A22 , we find that
the fundamental solution Γ of L can be approximated
by the (explicit) fundamental solution of the canonical sub-Laplacian L0 := X12 + X22 , with Jacobian basis
2
1
X1 = ∂c1 + c2 ∂c3 , X2 = ∂c2 − c2 ∂c3 on (SE(2))0 . This

solution in fact coincides with one of the approximations of Γ found by Duits & Franken [22]. There, the
fundamental solution of L was first approximated by
relying on a contraction of SE(2) to a 3-dimensional
Heisenberg group (via dilations on the group SE(2)),
and then derived the Gaussian estimates based on the
homogeneous norm k·k1 , i.e., ζ = 1, with exponential
coordinates derived from the contraction.
In our study on the sub-Riemannian distance approximations we found that even sharper estimates could
be obtained by relying on the explicit formula for the
fundamental solution of the (Kohn) sub-Laplacian on
H(3) (which is up to a constant given by kck−2
16 ). In
this context we thus obtain an estimate of the fundamental solution of L by estimating it with kck−2
16 , which
is proportional to the exact fundamental solution of L0
on our approximated group (SE(2))0 .
2.3.3 Approximation of the sub-Riemannian distance
Finally we arrive at the sub-Riemannian distance approximations. By the Ball-Box theorem (see e.g. [4])
and equivalence of homogeneous norms, there exists a
constant c such that
c−1 kLog(g)kζ ≤ d0 (e, g) ≤ ckLog(g)kζ ,

6

Erik J. Bekkers et al.

with Log(g) defined by Eq. (9). The logarithmic norm
is locally equivalent to the sub-Riemannian distance,
which was proved in full generality in [43, Thm. 2 &
4]. Via a scaling of the generators Ã2 = ξ −1 A2 and
Ã3 = ξ −1 A3 we define the ξ-isotropic norm
p
4
kckξ,ζ : = p
(|c1 |2 + |c̃2 |2 )2 + ζ|c̃3 |2
(15)
4
= (|c1 |2 + ξ 2 |c2 |2 )2 + ζ ξ 2 |c3 |2 ,
with c̃2 = ξc2 and c̃3 = ξc3 , and the ci given in (9). The
norm k·kξ,ζ closely approximates the sub-Riemannian
distance d0 (e, ·) for C = 1 (no data-adaptivity) via
d0 (g, h) ≈ | Log(g −1 h)|ξ,ζ ,

| Log(g)|ξ,ζ := kckξ,ζ (16)

with c the coordinates of the first kind obtained via
(9). In view of the Folland-Kaplan-Korányi gauge setting ζ = 16 in k·kξ,ζ would be a sensible choice. We
do observe however that ζ = 44 gives an even sharper
approximation, see Fig. 3 for a visual comparison to
the sub-Riemannian distance d0 and Appendix A for a
quantitative comparison. The setting ζ = 44 is used in
all experiments on SE(2).

In this section we extend the concepts of the previous
section to the group SE(3) of 3D translations and rotations. In the end we again obtain an approximation
for the sub-Riemannian distance, which allows us to do
perceptual grouping in 3D images as well.

3.1 The Lie Group SE(3)
3.1.1 SE(3)
The Lie group SE(3) = R3 o SO(3) is the semi-direct
product of the group of 3D translations R3 and the
group of 3D rotations SO(3). The group product and
inverse for elements g = (x, R), g 0 = (x0 , R0 ) ∈ SE(3)
are defined by

(17)

In the 3D case, we define the space of coupled positions
and orientations as a Lie group quotient of SE(3):
R3 o S 2 := SE(3)/(0 × SO(2)).
The group action of g ∈ SE(3) onto (y, n) ∈ R3 × S 2
is defined by
g · (y, n) = (x, R) · (y, n) = (x + Ry, Rn).

3.1.2 The Lie Algebra, Exponential Map and
Commutators
Analogously as in the SE(2) case, we associate with the
group SE(3) the Lie algebra se(3) using the exponential
and logarithmic maps. This is most easily done using
an isomorphism with the corresponding matrix group:


(x, Rγ,β,α ) ↔

Rγ,β,α xT
0
1


.

A basis for the corresponding matrix Lie-algebra is given
by

3 Sub-Riemannian Distance and its
Approximation in SE(3)

g · g 0 = (x, R) · (x0 , R0 ) = (x + Rx0 , RR0 ),
g −1 = (−R−1 x, R−1 ).

We can identify the element (x, n) ∈ R3 ×S 2 with group
elements (x, Rn ) ∈ SE(3)/(0 × SO(2)), where Rn is
any rotation matrix such that Rn ez = n.


0001
0 0 0 0

X1 = 
0 0 0 0 ,
0000


0000
0 0 0 0

X3 = 
0 0 0 1 ,
0000


0 010
 0 0 0 0

X5 = 
−1 0 0 0 ,
0 000




0000
0 0 0 1

X2 = 
0 0 0 0 ,
0000


00 0 0
0 0 −1 0

X4 = 
0 1 0 0 ,
00 0 0


0 −1 0 0
1 0 0 0

X6 = 
0 0 0 0 ,
0 0 00

(18)

and their equivalents Ai in the tangent space of SE(3)
span the Lie algebra se(3). Since it will be clear from
the context if we are in the SE(2) or SE(3) case, we
use the same notation for the generators of the Lie algebra as previously. Now the left-invariant vector fields
are again obtained using the push-forward of the leftmultiplication (Lg )∗ , but they depend on the choice of
coordinates. In this paper we mostly rely on ZY Z-Euler
angles in the parameterization of SO(3), i.e.,

Rγ,β,α = Rez ,γ Rey ,β Rez ,α ,

(19)

Title Suppressed Due to Excessive Length

7

with Rn,α a rotation with angle α around n. Then, the
left-invariant vector fields are given by

In the SE(3) case, a horizontal path is a curve γ :
R → SE(3) with tangent vectors γ̇(t) ∈ ∆|γ(t) :=
span{ A3 |γ(t) , A4 |γ(t) , A5 |γ(t) }, where ∆ is now the subbundle of full tangent bundle spanned by {Ai }6i=1 . In
this case we have one spatial control u3 and two ‘angular’ controls u4 and u5 , so that the sub-Riemannian
metric tensor becomes:

A1 |g = (cos α cos β cos γ − sin α sin γ)∂x
+(sin α cos γ + cos α cos β sin γ)∂y − cos α sin β∂z
A2 |g = (− sin α cos β cos γ − cos α sin γ)∂x
+(cos α cos γ − sin α cos β sin γ)∂y + sin α sin β∂z ,
A3 |g = sin β cos γ∂x + sin β sin γ∂y + cos β∂z
cos α
∂γ ,
A4 |g = cos α cot β∂α + sin α∂β −
sin β
sin α
A5 |g = − sin α cot β∂α + cos α∂β +
∂γ
sin β
A6 |g = ∂α ,

G ξ,C

(20)
for β 6= 0, π.
Remark 1 A second coordinate chart is needed to cover
the entire SO(3), for which for example ZY X-angles
can be used, as is done in e.g. [23], where also the expressions for the vector fields in this alternative coordinate chart are given. In fact, the basis elements Ai of
the Lie algebra correspond to partial derivatives with
respect to the ZY X-angles, similar to the SE(2)-case.
We can express each element se(3) in terms of the basis with coefficients c = (c1 , . . . , c6 )T . Furthermore, we
define c(1) := (c1 , c2 , c3 )T and c(2) := (c4 , c5 , c6 )T , the
spatial and rotational coefficients, respectively. We can
make the exponential map ExpSE(3) : se(3) → SE(3)
and logarithmic map LogSE(3) : SE(3) → se(3) explicit
using these coefficients. For a 3×3 matrix Ω of the form

0 −c6 c5
Ω :=  c6 0 −c4  ,
−c5 c4 0

3.2 Sub-Riemannian Geometry in SE(3)

γ(t)

(γ̇(t), γ̇(t)) := C(γ(t))2 ξ|u3 (t)|2 +

|u4 (t)|2 + |u5 (t)|2 ,

(24)

The sub-Riemannian distance between two elements
g1 , g2 ∈ SE(3) is still defined as in (7), but now the
infimum is taken over Lipschitz continuous curves γ ∈
Lip([0, T ], SE(3)) with γ(0) = g1 , γ(1) = g2 , and γ̇(t) =
u3 (t) A3 |γ(t) + u4 (t) A4 |γ(t) + u5 (t) A5 |γ(t) .

3.3 A Nilpotent Approximation (SE(3))0 of SE(3)
and the Approximated sub-Riemannian Distance
It is important to realize that the logarithmic map is
only well-defined on the group SE(3) and not on the
quotient R3 oS 2 , i.e., different choices for α in the rotational part result in different values for the coefficients
ci . Here, we choose the approach of [46] and set α = −γ
such that expected symmetries are preserved. With that
choice, the logarithm (23) gives for each (x, n) ∈ RoS 2
a unique vector c, on which we can put a norm:



(21)

we obtain a rotation using the exponential map of matrices, i.e., R = exp(Ω). The relation between the spatial coefficients c(1) and (x, R) is given by
c(1) =



1
q
I − Ω + q −2 1 − cot
2
2
2

 q 





(Ω)2 x, (22)

where q = ||c(2) || and Ω such that R = exp(Ω). Now

LogSE(3) (g) =

6
X

ci (g)Ai ,

i=1

ExpSE(3)

6
X

(23)

!
ci (g)Ai

and

= g,

i=1

using the relations above.

| LogSE(3) (g)|ξ,ζ := ||c||ξ,ζ :=
p
4
(ξ 2 |c3 |2 + |c4 |2 + |c5 |2 )2 + ζ (ξ 2 (|c1 |2 + |c2 |2 ) + |c6 |2 ),
(25)
where c = c(g) according to (23).
Also here, the Folland-Kaplan-Korányi-type norm
can be used to approximate the fundamental solutions
of the sub-Laplacian on SE(3). The norm ||c||ξ,ζ with
ζ = 1 was for example used in [23] approximations of
the heat kernel and the fundamental solution on SE(3),
of which only recently exact solutions were found in
[47]. In the context of this paper, we can approximate
the exact solutions of the sub-Laplacian on SE(3) by
kck2−Q
1,16 , with homogeneous dimensions Q = 9, as the
exact solution of the sub-Laplacian on the approximation group (SE(3))0 . The group (SE(3))0 that locally
approximates SE(3), is again obtained via a nilpotent

8

Erik J. Bekkers et al.

step 2 approximation of the BCH formula, and is defined by the group product
(x1 , x2 , x3 , x4 , x5 , x6 ) · (y 1 , y 2 , y 3 , y 4 , y 5 , y 6 )
 1
T
x + y 1 + 12 (x5 y 3 − x3 y 5 )
 x2 + y 2 + 1 (x3 y 4 − x4 y 3 ) 
2




x3 + y 3

 , (26)
=
4
4

x
+
y


5
5


x +y
1
6
6
4 5
5 4
x + y + 2 (x y − x y )
with xi , y i coordinates of the first kind given by the
logarithmic map (22). This new group is a free-nilpotent
group of rank 3 and step 2.
We approximate the sub-Riemannian distance d0 on
SE(3) via the norm (25). I.e.,
d0 (g, h) ≈ | LogSE(3) (g −1 h)|ξ,ζ ,

(27)

and as such again obtain an approximation of the distance in the sense of Rothschild and Stein [50]. Based on
the quantitative comparison to the sub-Riemannian distances d0 in Appendix A and the visualizations in Fig. 4
of the level sets we conclude that the approximated subRiemannian distance of (27) quite accurately approximates the true sub-Riemannian distance on SE(3). In
our analysis we found that the logarithmic norm with
ζ = 100 gave the best approximation, and as such we
used this norm in the perceptual grouping experiments
of Sec. 5.2.
Remark 2 The glyph at each grid point y in Fig. 4 is
given by the surface {y + νU (y, n)n|n ∈ S 2 }, for a
specific choice ν > 0, and with density U : R3 × S 2 →
R+ . The color of each orientation n = (n1 , n2 , n3 ) ∈ S 2
on the glyph is defined by the RGB color (n1 , n2 , n3 ).

4 Perceptual Grouping, Fast Marching and
Key Point Tracking
In this section the algorithms used in this paper are
explained. Our main application of interest is that of
grouping/clustering of points on blood vessels via the
perceptual grouping algorithm, which is explained in
Subsec. 4.1. The perceptual grouping algorithm takes
as input a set of key points that are obtained via the
minimal path tracking with key points algorithm [5], explained in Subsec. 4.3, which is an adaptation of the fast
marching algorithm, explained in Subsec. 4.2. Finally
since different metrics are used throughout the experiments (both for generating key points and for perceptual grouping) we end this section with an overview of
the used metrics in this paper in Subsec. 4.4.

4.1 The Perceptual Grouping Algorithm
The perceptual grouping algorithm presented in this
paper is a modification of the original algorithm proposed by Cohen [13], and which was later adapted for
perceptual grouping based on anisotropic distances [8].
In recent work [11], the perceptual grouping algorithm
was extended for the grouping of n closed contours for
an a-priori specified n. Like in [8] and [11], we use the
main algorithm of [13] as a backbone, but we change
the metric used for perceptual grouping and we impose
an additional constraint to avoid closed loops (which
are physically not realistic in the vessel networks of interest). Our adapted perceptual grouping algorithm is
given in pseudo code in Algorithm 1.
The goal of the perceptual grouping algorithm is
to construct a graph out of a set S of points of interest in which the edges DS are true connections (represented by geodesics) between points. Following the
terminology of [17, 5, 10] we will refer to the points of
interest as key points. Each key point is only linked to
at most 2 other key points (i.e., node degree δi is 2 at
most). The final graph thus only contains sets of nonbifurcating vessel segments. The graph is build up by
inserting one-by-one the edges which have the shortest
geodesic distance (if the node degree allows). As such,
only the strongest connections (shortest geodesics) appear in the final graph network. Since the original algorithm in [13] (and also [8]) does not include a mechanism to avoid closed loops we include an additional
check in the main algorithm to prevent this. Finally, in
order to avoid connecting key points which are too far
apart from each other we only consider edges of which
the spatial arc length of the connecting geodesic does
not exceed a certain a-priori threshold smax .
In summary, our changes relative to the works [13,
8, 11], is that we
– keep the choice for distance d(xi , xj ) open. In our
experiments the distances d will be mainly based on
sub-Riemannian geometry in SE(n).
– explicitly avoid making long distance connections by
filtering out such possible connections in an initialization step.
– avoid closed loops by not making connections between nodes that belong to the same sub-graph.
– group crossing lines without pre-specifying the number of groups.
In particular, it is the use of a sub-Riemannian metric
on SE(n) that allows for the grouping of crossing lines.
A first (successful) feasibility study on the possibility of
perceptual grouping of crossing lines was performed by
Chen et al. [11] using a (sub-)Finsler metric (based on
the Euler elastica model) on position-orientation space.

Title Suppressed Due to Excessive Length

9

Fig. 4 Distances on SE(3) for ξ = .1, C = 1, with the origin placed at e = (0, ex ). Top row: Level sets of the spatial
projections (minimum intensity projections over S 2 ) of the distance volumes on SE(3). Rows two to four: Glyph visualizations
2
in which each distance volume d is visualized with a ”Gaussian” density U (g) = e−d(e,g) . For an interpretation of the glyphs
see Remark 2. Row two: Glyph visualizations of sub-volume. Row three: Glyph visualization of slice at z = 0. Row four:
Zoomed in glyph visualization of the slice a z = 0. From left to right: The sub-Riemannian distance d0 (e, ·) on SE(3), see
Eq. (7) and (24); Homogenous norms k·kξ,ζ , see Eq. (25), of the nilpotent approximation (SE(3))0 for respectively ζ = 100,
ζ = 16 (Folland-Kaplan-Korányi gauge) and ζ = 1; The (ξ-isotropic) Riemannian distance d1 (e, ·) on SE(3), see Table 1 for
an overview of the different distances.

There it was successfully demonstrated on phantom images that their algorithm is able to deal with crossing
closed contours, however, it required specification of
the number of contours (which is not always a-priori
known). Furthermore, their metric relies on a notion of
directionality (instead of just orientations) which is useful in grouping closed contours, but may be disadvantages for grouping non-closed contours. Here, we focus
on the grouping of non-closed crossing contours without specifying the number of contours. Furthermore,
we quantify the performance of perceptual grouping of
crossing lines on a large set of both retinal images in
2D, and phantom images in 3D.

4.2 Fast Marching
Most of the distances (except for the fast analytic approximations) and the geodesics used in this paper are
computed via the fast marching algorithm, which is an
efficient numerical solver of the eikonal equation and
which can be used to obtain globally optimal geodesics
[14]. Let g0 be an arbitrary source point in a domain M
of interest, let G|g : Tg (M) × Tg (M) → R+ be a metric
tensor defined on the tangent space Tg (M) at g ∈ M,
and let
Z 1q
G|γ(t) (γ̇(t), γ̇(t))dt (28)
U (g):=d(g0 , g)= inf
γ∈S(g0 ,g)

0

its associated distance map, where the infimum is taken
over the set S(g0 , g) of Lipschitz continuous curves with

10

Erik J. Bekkers et al.

fast marching, and to [52] and [24] for fast marching in
sub-Riemannian manifolds in SE(2) and SE(3) respectively.

: S: a set of key points;
d(gi , gj ): distances between gi , gj ∈ S;
smax : max spatial length of geodesics;
variables: D̃S : set of possible edges;
δi : node degree of xi ;
output : DS : final set of edges;

input

4.3 Generating Key Points

Initialization:
Compute the distances d(gi , gj ) (and corresponding
geodesics) between all key points gi , gj ∈ S.
Initialize D̃S with the set of all edges between each
gi , gj ∈ S whose connecting geodesic has spatial arc
length smaller then smax , and set DS = Ø.
Main algorithm:
while D̃S 6= Ø do
1. Select edge and remove it from D̃S :
(gi , gj ) ← argmin d(g, h);
(g,h)∈D̃S

D̃S = D̃S − (xi , xj );
2. Check topology and update network :
if δi < 2 and δj < 2 and gi , gi are not
already in the same sub-graph in DS
then DS = DS + (gi , gj );
δi = δi + 1;
end

Algorithm 1: Perceptual grouping.

γ(0) = g0 , γ(1) = g, and with γ̇(t) ∈ Tγ(t) (M). Then
the distance map U is the unique viscosity solution of
the eikonal equation
p

k∇G U (g)kG = 1,
G (∇G U (g), ∇G U (g)) = 1,
↔
U (g0 ) = 0,
U (g0 ) = 0,
(29)

The key point method is based on keeping track of a
spatial arc-length map Ul (in which the spatial lengths
of the minimizing geodesics γ defining U are stored),
and stop as soon as a certain distance threshold is
passed [17]. The rationale behind this algorithm is that
among all points with equal geodesic distance values
U , the points reached by geodesics γ that best follow
the data (paths along which C is low) have maximum
spatial distance l(γ). Such a point maximizing spatial
distance in a given level set in U is called a key point.
The fast marching algorithm is highly suited for keeping track of a spatial arc-length map Ul , in addition
to U , due to the local updating approach (wavefront
propagation). Moreover, the algorithm can stop early if
one is only interested in finding the first key point with
length larger than lmax [17].
In summary, a key point is detected as follows. The
spatial arc-length map is defined as
Ul (g) := l(γg0 ,g ),
with γg0 ,g = argmin
γ∈S(g0 ,g)

For details on the fast marching algorithm on isotropic
manifolds we refer to [56, 54], to [41, 32] for anisotropic

R1q
0

G|γ(t) (γ̇(t), γ̇(t))dt the min-

imizing geodesic in (28), and with
Z

1

kẋ(t)kdt

l(γ) =

(31)

0

−1

with ∇G := G dU the intrinsic gradient with inverse
metric G −1 and dU the differential of U , and k·kG the
norm with respect to the metric tensor. In the standard (data-adaptive) Euclidean case with M = R2 ,
g0 = 0, g = x, γ̇(t) = u1 (t)∂x +u2 (t)∂y ∈ Tγ(t) (R2 ), and
with G|γ(t) (γ̇(t), γ̇(t)) = C(γ(t))2 (|u1 (t)|2 +|u2 (t)|2 ) the
eikonal equation is given by k∇U (x)k = C(x).
The fast marching algorithm efficiently solves the
eikonal equation in a one pass algorithm. It computes
the values of U in increasing order (starting with U (g0 ) =
0) based on the Bellman principle of optimality, in a
manner very similar to the Dijkstra algorithm for shortest paths on graphs [18]. The minimal geodesic connecting g0 with g is then obtained via a gradient descent on
U from g back to the origin g0 , i.e., solving the ODE

γ̇(t) ∝ −G −1 dU (γ(t)),
γ(0) = g0 .

(30)

the spatial arc-length of γ, with ẋ(t) = PRn γ̇(t) ∈ Rn
the spatial components of the tangents γ̇(t)4 . The fast
marching algorithm stops as soon as there is a g for
which Ul (g) ≥ lmax , and g will be called a key point.
With the above criteria one can iteratively detect
new key points based on geodesic distances to previously found key points, a method known as minimal
path tracking with key point detection [5]. One can make
several choice on when to stop the key point tracking algorithm [5, 34, 10]. In this work, we rely on the approach
by Chen et al. [10], where we only add key points on
locations which lie in a masked region (we use a binary
vessel centerline mask m : M → {0, 1}), i.e., we only
add a key point when both Ul (g) ≥ lmax and m(g) = 1.
The algorithm is stopped as soon Ul (g) ≥ 3 lmax .
4

In the lifted problem SE(2) the spatial components are
for example given by ẋ(t) = u2 (t)A2 |γ(t) + u3 (t)A3 |γ(t) ,
and in the SE(3) case ẋ(t) = u1 (t)A1 |γ(t) + u2 (t)A2 |γ(t) +
u3 (t)A3 |γ(t) .

Title Suppressed Due to Excessive Length

4.4 Overview of Distances Used in This Paper
Table 1 gives an overview of the different distances
discussed in this paper and used in the experiments.
The isotropic Euclidean metrics are used the generate
key points in R2 and R3 using the algorithm of Subsec. 4.3. The isotropic Euclidean distances are also used
in comparison to the other distances on SE(n) in the
perceptual grouping experiments. The sub-Riemannian
distances on SE(2) and SE(3) are explained respectively in Subsec. 2.2 and Subsec. 3.2. In the Riemannian
distances the full tangent bundle on SE(n) is considered. This means that now also non-horizontal curves in
SE(n) are considered, i.e., points on the curves γ are allowed to move sideways by the non-horizontal controls
u3 (t) in the SE(n) case, and u1 (t), u2 (t) in the SE(3)
case. Recall that in this case the blue and red oriented
particles in Fig. 1 do have the same distance to the
source (black arrow). Finally, the sub-Riemannian distance approximations, denoted with | LogSE(n) (g −1 h)|ξ,ζ
are discussed and defined in respectively Subsec. 2.3
and Eq. (16) for SE(2) and Subsec. 3.3 and Eq. (25)
for SE(3).
4.4.1 The Cost C
The cost functions C are constructed from functions
Uf : Rn × S n−1 → R on the orientation-lifted space.
These functions Uf are obtained via an orientation score
transform [21, 31] of image f : Rn → R by correlating
the image with a set of anisotropic wavelets ψ : Rn →
R:
Uf (g) = (Ug ψ, f )L2 (Rn ) ,
(32)
R
with (f, g)L2 (Rn ) = Rn f (x)g(x)dx the standard inner
product on L2 (Rn ), with the overline denoting complex
conjugation, and where Ug denotes the left-regular representation of the Lie group on images f . For the group
SE(2) acting on images f ∈ L2 (R2 ) it is defined as
(Ug f )(y) := f (R−1
θ (y − x))
with g = (x, θ) ∈ SE(2) (recall the group definitions in
Subsec. 2.1.1). For the group SE(3) acting on images
f ∈ L2 (R3 ) it is defined as
(Ug f )(y) := f (R−1
n (y − x))
with g = (x, Rn ) ∈ SE(3) (recall the group definition
in Subsec. 3.1.1).
The wavelets used in the orientation score transform
[21,31] are designed in such a way that all rotated version together cover the full Fourier spectrum. With this
design no data is lost in the transformation and a stable invertible transform (from orientation score) back to

11

image exists. For details on this wavelet-type transform
for lifting 2D images to functions on SE(2) we refer to
[21], and for lifting 3D images to 3D orientation scores
we refer to [31]. In all experiments we define the cost in
the following form
C(g) =

1
,
1 + λV(g)p

(33)

with V a vessel (or center line) enhancement obtained
by processing of the orientation score Uf , and which
is normalized between 0 and 1. Parameters λ and p
then control respectively the influence of the cost (dataadaptivity) and p the contrast.
Good choices for V for tracking of vessels in 2D position orientation space may be via the vessel enhancements of [58] or [30], similar to the SE(2) tracking experiments in [3]. For tracking in 3D orientation scores
V may be obtained via the crossing preserving vessel
enhancements of [19]. In related tracking problems in
lifted spaces the lifts are obtained via tubularity measures [11, 36, 38], or by correlating the image with a set
of rotated templates [44].
4.4.2 Projective Line Bundle
Finally, we remark that when dealing with geodesic distances in SE(n) we have to take into account that these
are defined for positions and orientations on the full
sphere S n−1 . The distances discussed in this paper thus
make a distinction between forward and backward arrival directions, i.e., d(e, (x, θ)) 6= d(e, (x, θ + π)).
In practice, and in particular in our perceptual grouping problem, we often do not know the direction of the
vessel, but we only have orientations. As such, we would
actually want to compute distances on the projective
line bundle Rn × P n−1 , with P n−1 := S n−1 / ∼ with
identification of antipodal points n1 ∼ n2 ↔ n1 = ±n2 .
We define the distances d˜ on the projective line bundle
by distances d on SE(n) via
˜ (x, n)) = min {d(g, (x, n)), d(g, (x, −n))} ,
d(g,

(34)

with n ∈ S n−1 , and g, (x, ±n) ∈ SE(n). Note that in
the SE(2) case we have with n(θ) = (cos θ, sin θ) ↔ θ
and −n(θ) = n(θ + π). For a more detailed analysis
on data-adaptive sub-Riemannian geodesics on the 2D
projective line bundle we refer [2].

5 Experiments
In the experiments we aim to quantify the performance
of perceptual grouping with different distances. For a
fair comparison we therefore generate automatically the

12

Erik J. Bekkers et al.

Table 1 Overview of the metrics used in this paper.
Distance
notation

Manifold
M

Tangent b.
T (M)

Metric tensor G

Tangent vectors γ̇
Isotropic Euclidean∗

kg − hk

R

2

2

1

T (R )

γ̇(t) = u (t)∂x + u2 (t)∂y

kg − hk

R3

T (R3 )

γ̇(t) =

G|γ(t) = C(γ(t))2 (|u1 (t)|2 + |u2 (t)|2 )

u1 (t)∂x + u2 (t)∂y
+u3 (t)∂z

G|γ(t) = C(γ(t))2



|u1 (t)|2 + |u2 (t)|2
+|u3 (t)|2



3
P



(Full) Riemannian SE(n)
d1 (g, h)

SE(2)

T (SE(2))

γ̇(t) =

3
P

ui (t)Ai |γ(t)

G|γ(t) = C(γ(t))2

i=1

i=1


ξi2 |ui (t)|2 ,

with ξ2 = ξ3 = ξ and ξ1 = 1
d1 (g, h)

SE(3)

T (SE(3))

γ̇(t) =

5
P

ui (t)Ai |γ(t)

G|γ(t) = C(γ(t))2



5
P

i=1

i=1


ξi2 |ui (t)|2 ,

with ξ1 = ξ2 = ξ3 = ξ and ξ4 = ξ5 = ξ6 = 1
Sub-Riemannian SE(n)
1

d0 (g, h)

SE(2)

∆

γ̇(t) = u (t)A1 |γ(t) + u2 (t)A2 |γ(t)

d0 (g, h)

SE(3)

∆

γ̇(t) =

G|γ(t) = C(γ(t))2 (|u1 (t)|2 + ξ 2 |u2 (t)|2 )

u3 (t)A3 |γ(t) + u4 (t)A4 |γ(t)
+u5 (t)A5 |γ(t)

G|γ(t) = C(γ(t))2



ξ 2 |u3 (t)|2 + |u4 (t)|2
+|u5 (t)|2



Sub-Riemannian approximation
−1

| LogSE(2) (g h)|ξ,ζ
| LogSE(3) (g −1 h)|ξ,ζ

Approximation of the sub-Riemannian distance on SE(2), cf. Eq. (16)
Approximation of the sub-Riemannian distance on SE(3), cf. Eq. (25)

∗

The isotropic Euclidean distances are used in key point generation and perceptual grouping. The other distances
are only used in the perceptual grouping algorithm.

most reasonable key points by using a vessel center line
mask m : Rn → [0, 1] (see Subsec. 4.3) based on the
ground truth data. Moreover, this guarantees that the
key points are always located on the ground truth center lines, which allows us to quantify performance using the ground truth data. In both the 2D and 3D case
the key points are then generated using the isotropic
Euclidean metric tensor, and with V(x) = m(x) (see
Subsec. 4.4.1). In all experiments we set p = 1, λ = 100
to compute the cost (cf. Eq. (33)).
In the perceptual grouping experiments the cost functions are constructed from orientation score transforms
Uf of the mask m on Rn . The costs on SE(n) are then
constructed via the modulus of the score:
V(g) = VSE(n) (g) := |Uf (g)|.

(35)
n

For equal comparison the costs on R are then constructed via V (x) = max
VSE(n) (x, n), i.e., via a maxn−1
n∈S

imum intensity projection over orientations n.
5.1 Perceptual Grouping in SE(2)
5.1.1 Experimental Setup
The data for the 2D retinal vessel grouping experiments
consists of 52 retinal image patches in which the ves-

sels have complicated topologies (each patch contains
at least 1 crossing, and at least 1 bifurcation). For each
retina patch the center lines were semi-automatically
traced, after which the connectivity (bifurcation relations) between the vessel segments were manually determined. The set of images contained in total 313 separate vessel segments. A connection between two nodes
was determined to be a true positive if both nodes lie
on the same vessel tree.

The minimum distance between key points in the
retina experiments (with patch sizes of ≈ 400 × 400
pixels) was set to lmax = 30 pixels. The maximum
geodesic arc length distance in the perceptual grouping
algorithm was set to smax = 80 pixels. The orientations
θ at each key point x was estimated by the orientation
that gave maximum response in the orientation score,
i.e., θ = argmax VSE(2) (x, θ). The circle S 1 was samθ∈S 1

pled with Nθ = 32. All distances were computing via
the fast marching algorithm of [41, 42] except for the
sub-Riemannian approximations, which were computed
directly using (9) and (15). The position-orientation
balancing parameter was set to ξ = 0.01.

Title Suppressed Due to Excessive Length

13

Fig. 5 Example 1 of the retinal vessel grouping experiments. Each connected component has its own color (note that the colors
might not match between experiments as the number of recovered components may differ), false connections are indicated
in red. Top row: experiments with data-adaptive distances (C 6= 1), and the ground truth vessel components including the
automatically generated key points. Bottom row: experiments without data-adaptive distance (C = 1).

Fig. 6 Example 2 of the retinal vessel grouping experiments. Each connected component has its own color (note that the colors
might not match between experiments as the number of recovered components may differ), false connections are indicated
in red. Top row: experiments with data-adaptive distances (C 6= 1), and the ground truth vessel components including the
automatically generated key points. Bottom row: experiments without data-adaptive distance (C = 1).

14

Erik J. Bekkers et al.

Table 2 Perceptual grouping performance for the 2D retinal
image experiments in terms of percentage of correct key point
connections (# of false connections in parenthesis).
Distance
kx − yk
d1 (g, h)
d0 (g, h)

C=1

89.99%
97.51%
(SE(2)|Sub-Riem.)
99.75%
| Log(g −1 h)|ξ,ζ (SE(2)| ≈Sub-Riem.) 99.72%
(R2

|Eucl.)

(SE(2)|Riem.)

(362)
(90)
(9)
(10)

Table 3 Perceptual grouping performance for the 3D synthetic volume experiments in terms of percentage of correct
key point connections (# of false connections in parenthesis).

C 6= 1

Distance

95.96% (146)
99.64% (13)
99.83% (6)
-

kx − yk
d1 (g, h)
d0 (g, h)

5.1.2 Results

C 6= 1

C=1

89.99%
93.02%
(SE(3)|Sub-Riem.)
96.79%
| Log(g −1 h)|ξ,ζ (SE(3)| ≈Sub-Riem.) 97.17%
(R3

|Eucl.)

(SE(3)|Riem.)

(78)
(54)
(25)
(22)

97.97% (16)
98.32% (13)
98.32% (13)
-

5.2 Perceptual Grouping in SE(3)
5.2.1 Experimental setup

Table 2 gives a quantitative overview of the results,
Figs. 5 and 6 shows the results on two of the 52 retina
patches. From Table 2 we make the following observations and conclusions:
1. Perceptual grouping is preferred in the lifted domain SE(2) instead of in the based domain R2 . This
suggest that taking orientation into account in the
grouping is essential.
2. A sub-Riemannian geometry on SE(2) is preferred
over a (ξ-isotropic) Riemannian geometry. This suggests that a sub-Riemannian geometry is necessary
to deal with the complex geometry at crossings and
parallel tracks (cf. Figs. 5 and 6).
3. The results obtained with the sub-Riemannian distances on SE(2) for C = 1 are almost equal. This
suggests that the approximations are quite accurate,
and that for C = 1 the analytic approximations may
be preferred due to speed and algorithm complexity
considerations.
4. Overall, results for C 6= 1 are better than for C = 1.
Note however, that the sub-Riemannian distances
on SE(2) for C = 1 are still better then the Euclidean distance on R2 and Riemannian distance on
SE(2) for C 6= 1, and only slightly under performs
relative to the sub-Riemannian C 6= 1 case. This
again shows that sub-Riemannian geometry is preferred, whether data is included in the metric tensors or not.
We conclude that in perceptual grouping of 2D vessels a sub-Riemannian geometry in SE(2) is preferred
over a Euclidean geometry in R2 , or a Riemannian geometry in SE(2). When accurate vesselness maps are
available, it is preferable to use these in the distances.
Furthermore, if one aims to design a easy to implement
and efficient perceptual grouping pipeline, approximate
sub-Riemannian distances should be used. With only a
2D key point tracking algorithm, a method for estimating orientations, and the analytic approximate distances (15) one obtains very accurate grouping results.

To quantify and study the influence of different distances in perceptual grouping algorithms for 3D vessels
we make use of synthetic 3D images. For these experiments 10 volumes were generated, each with 6 random
paths. Each path was generated with a Monte-Carlo
simulation of a random walk in SE(3) (see e.g. [59,
Ch. 3.5]). Due to the random construction it might occur that 2 paths cross each other. This is physiologically
unrealistic (vessels in 3D might bifurcate or touch, but
never grow through each other), but it does make the
experiments more challenging.
For each volume a binary center line mask was constructed using the generated ground truth paths. The
volumes were of size 51×51×51 voxels. The distance between key points was set to lmax = 5 voxels. The maximum geodesic arc length distance in the perceptual
grouping algorithm was set to smax = 15 voxels. The
orientation at each key point was again estimated as
the orientation that gave maximum response in VSE(n)
(Eq. (35)). The sphere S 2 was sampled with 200 orientations using Euler angles with n(β, γ) = Rγ,β,α .ez ,
π
π
π
with β ∈ { 2N
, 2 2N
, ..., π − 2N
}, γ ∈ {0, Nπβ , ..., 2π −
β
β
β
π
Nβ }, with Nβ = 10, and with Rγ,β,α given by (19). In
the lifted metric tensor we set ξ = 1.
5.2.2 Results
Table 3 gives a quantitative overview of the results,
Figs. 7 shows the results on one of the ten synthetic volumes. From Table 3 we can draw the same conclusions
as for the SE(2) case (using a sub-Riemannian geometry and including data adaptivity improves results).
Here, however, we make two additional observations
1. Data-adaptive fast marching seems less sensitive to
the choice of metric, but tracking in the lifted domain SE(3) still improves results. This can be explained by the fact that the volume is relatively
sparse, and by the fact that the cost function C is
constructed from ground truth data (the best possible cost). If the cost function dominates the metric,

Title Suppressed Due to Excessive Length

15

then the intrinsic energy/geometry has a smaller influence.
2. Out of all C = 1 distances (no-data adaptivity)
the grouping via the nilpotent distance approximations on SE(3) give best performance, even better then for the true sub-Riemannian distance. This
can be explained by the fact that for long distances
from the origin, the approximation gradually loose
their sub-Riemannian nature and allows more nonhorizontal behavior, as in the Riemannian case. It
could be that, due to the discrete sampling of the
sphere, not all orientations are accurately estimated.
The grouping based on the sub-Riemannian distance
approximations seems less sensitive to such errors.

using sub-Riemannian distances in such problems is
demonstrated by the experiments of this paper.

6 Conclusion

A Optimization of the Folland-Kaplan-Korányi
gauge parameter ζ

In this paper we have proposed an efficient approach
for perceptual grouping of local orientations via nilpotent approximations of sub-Riemannian distances in the
roto-translation group SE(n). The quantitative experiments on grouping of retinal blood vessels in 2D images, and perceptual grouping in challenging 3D synthetic volumes, showed that 1) sub-Riemannian geometry is essential in achieving top performance and 2)
that the grouping approach via the fast analytic approximations performs almost equally, or better, than
the data-adaptive fast marching approaches.
The sub-Riemannian distances on SE(2) and SE(3)
were approximated via norms on exponential coordinates of the first kind (obtained via the logarithmic
map). In both quantitative and visual comparison it
was found that the approximations accurately follow
the true sub-Riemannian distances, a conclusion which
was further supported by the equal performance in quantitative perceptual grouping experiments. We also numerically showed that the weighted logarithmic norms
used in this paper provide a more accurate approach
for approximating the heat kernel and fundamental solution of the sub-Laplacian on SE(n), compared to previous approaches [22, 47, 12].
Since the sub-Riemannian distance approximations
are analytic, they are easy to implement and fast to
compute. An interesting line of further research would
be to embed the sub-Riemannian distances in other
algorithms that rely on the quantification of the distance between local orientations. The results of this
paper could be further improved by augmenting the
sub-Riemannian distances with additional features (like
cross-sectional profile descriptors) and use a global graph
optimization approach as in [57, 26]. The potential of

Acknowledgements The following are gratefully acknowledge for their influence on the manuscript: Remco Duits at
Eindhoven University of Technology for suggestions on logarithmic approximations of the heat kernel and the fundamental solution on SE(2), SE(3) and the Heisenberg group;
Laurent Cohen at University Paris Dauphine for fruitful discussions on geodesic methods and perceptual grouping; JeanMarie Mirebeau at Laboratoire de mathématiques dOrsay,
Université Paris-Saclay for providing efficient and generic fast
marching code. The reviewers are gratefully acknowledged for
their valuable feedback on the manuscript. The research leading to the results of this article has received funding from the
European Research Council under the European Communitys 7th Framework Programme (FP7/20072014)/ERC grant
agreement No. 335555 (Lie Analysis).

In Figs. 3 and Fig. 4 we visually compared the nilpotent approximations of the sub-Riemannian distance on SE(2) and
SE(3) respectively. In this appendix we support by means of
a quantitative comparison our choices for ζ = 44 and ζ = 100
which appear in the logarithmic approximations (FollandKaplan-Korányi gauge) of | Log(g −1 h)|ξ,ζ as defined in (16)
and (27) on respectively SE(2) and SE(3).

A.1 Optimization of ζ for the SE(2) approximations
In the quantitative comparison on SE(2) we computed the L2
error between d0 (g, h) and the approximation | Log(g −1 h)|ξ,ζ
with ξ = 1 on a grid with a varying spatial domain size (from
x, y ∈ [−0.5, 0.5] to x, y ∈ [−4, 4]), and with varying choices
of ζ. The results are shown in Fig. 8 and are computed as
follows.
The reference sub-Riemannian distance d0 on SE(2) was
computed once via an anisotropic fast marching algorithm
[41, 42] on a grid which sampled x, y ∈ [−4, 4] at a sub-pixel
resolution of 0.01 with 128 orientations. The numerically computed sub-Riemannian distance volume was thus of dimensions 801 × 801 × 128.
In each experiment (with fixed spatial range and ζ) the
squared error between d0 and its approximation was sampled
on a regular grid that covered the specified domain with 41 ×
41 × 128 points. The averaged errors are plotted in Fig. 8.
Here we see that the approximation becomes more accurate towards the origin (x, y ∈ [−0.5, 0.5]) and that parameter
ζ has to be chosen larger in order to keep the anisotropy for
longer distances from the origin. The choice ζ = 44 generally
gave the best approximations and we rely on this setting in
the experiments on SE(2).

A.2 Optimization of ζ for the SE(3) approximations
In the quantitative comparison on SE(3) we computed the L2
error between d0 (g, h) and the approximation | Log(g −1 h)|ξ,ζ
with ξ = 1 on a grid with a varying spatial domain size (from

16

Erik J. Bekkers et al.

Fig. 7 Example 1 of the 3D synthetic vessel grouping experiments. Each connected component has its own color (note that
the colors might not match between experiments as the number of recovered components may differ), false connections are
indicated in red. Top row: experiments with data-adaptive distances (C 6= 1), and the ground truth vessel components including
the automatically generated key points. Bottom row: experiments without data-adaptive distance (C = 1).

Mean squared error

0.14
0.12
0.10

range = ±0.25

0.08

range = ±0.5

0.06

range = ±1.

0.04

range = ±2.

0.02

range = ±4.

0.00

0

20

40

60

80

100

ζ

Fig. 8 Mean squared errors between the sub-Riemannian
distance on SE(2) (see Eq. (7)) and its approximation (see
Eq. (16)). The errors is computed for varying choices of ζ and
on a varying grid size (from x, y ∈ [−0.5, 0.5] to x, y ∈ [−4, 4]).

x, y, z ∈ [−0.5, 0.5] to x, y, z ∈ [−4, 4]), and with varying
choices of ζ. The results are shown in Fig. 9 and are computed
as follows.
The reference sub-Riemannian distance d0 on SE(3) was
also computed once via an anisotropic fast marching algorithm [41, 42] on a grid which sampled x, y, z ∈ [−4, 4] at
a sub-pixel resolution of 0.1 with 31 × 62 Euler angles (cf.
Sec. 5.2). The numerically computed sub-Riemannian distance volume was thus of dimensions 201×201×201×31×62.
In each experiment (with fixed spatial range and ζ) the
squared error between d0 and its approximation was sampled
on a regular grid that covered the specified domain with 21 ×

21 × 21 × 31 × 62 points. The averaged errors are plotted in
Fig. 9.
Here we see that the approximation becomes more accurate towards the origin (x, y, z ∈ [−0.5, 0.5]). However, compared to the SE(2) experiments we do see a less stable localization of the optimal parameter ζ with varying spatial
resolutions. This behaviour can be explaind by 1) the subRiemannian distances are numerically computed via a fast
marching algorithm using Euler angles (which do not uniformly sample the sphere), and 2) the spatial resolution of
the computed reference sub-Riemannian distance volume was
only 0.1 (due to computer memory constraints). Although
very accurate from an application point of view, the subRiemannian distances on SE(3) are not exact, and the numerical errors induced by the algorithm may explain the variation in optimal ζ (in particular for the region close to the
origin). Overall, the choice ζ = 100 seems to be reasonable
in all ranges, and this was confirmed by visual comparison of
the distance maps in Fig. 4.

References
1. Abbasi-Sureshjani, S., Zhang, J., Duits, R., ter
Haar Romeny, B.: Retrieving challenging vessel connections in retinal images by line co-occurrence statistics.
Biological Cybernetics (2017)

Title Suppressed Due to Excessive Length

17

Mean squared error

0.20

0.15

range = ±0.5
0.10

range = ±1.
range = ±2.
range = ±4.

0.05

0.00
60

80

100

120

140

ζ

Fig. 9 Mean squared errors between the sub-Riemannian
distance on SE(3) (see Eq. (7) and (24)) and its approximation (see Eq. (27)). The errors is computed for varying choices
of ζ and on a varying grid size (from x, y, z ∈ [−0.5, 0.5] to
x, y, z ∈ [−4, 4]).

2. Bekkers, E., Duits, R., Mashtakov, A., Sachkov, Y.: Vessel tracking via sub-riemannian geodesics on R2 × P 1 .
arXiv preprint arXiv:1704.04192 (2017)
3. Bekkers, E.J., Duits, R., Mashtakov, A., Sanguinetti,
G.R.: A PDE approach to data-driven sub-Riemannian
geodesics in SE(2). SIAM Journal on Imaging Sciences
8(4), 2740–2770 (2015). DOI 10.1137/15M1018460
4. Bellaı̈che, A.: The tangent space in sub-Riemannian geometry, pp. 1–78. Birkhäuser Basel, Basel (1996)
5. Benmansour, F., Cohen, L.D.: Fast object segmentation
by growing minimal paths from a single point on 2d or
3d images. Journal of Mathematical Imaging and Vision
33(2), 209–221 (2009)
6. Bonfiglioli, A., Lanconelli, E., Uguzzoni, F.: Stratified
Lie groups and potential theory for their sub-Laplacians.
Springer Science & Business Media (2007)
7. Boscain, U., Duits, R., Rossi, F., Sachkov, Y.: Curve
cuspless reconstruction via sub-Riemannian geometry.
ESAIM: Control, Optimisation and Calculus of Variations 20(3), 748–770 (2014)
8. Bougleux, S., Peyré, G., Cohen, L.: Anisotropic geodesics
for perceptual grouping and domain meshing. In: European Conference on Computer Vision, pp. 129–142.
Springer (2008)
9. Chen, D.: New minimal path models for tubular structure
extraction and image segmentation. Ph.D. thesis, Universite Paris Dauphine, PSL Research University (2016)
10. Chen, D., Mirebeau, J.M., Cohen, L.D.: Vessel tree extraction using radius-lifted keypoints searching scheme
and anisotropic fast marching method. Journal of Algorithms & Computational Technology 10(4), 224–234
(2016)
11. Chen, D., Mirebeau, J.M., Cohen, L.D.: Global minimum for a finsler elastica minimal path approach. International Journal of Computer Vision 122(3), 458–483
(2017). DOI 10.1007/s11263-016-0975-5. URL http:
//dx.doi.org/10.1007/s11263-016-0975-5
12. Citti, G., Sarti, A.: A cortical based model of perceptual completion in the roto-translation space. Journal of
Mathematical Imaging and Vision 24(3), 307–326 (2006)
13. Cohen, L.D.: Multiple contour finding and perceptual
grouping using minimal paths. Journal of Mathematical Imaging and Vision 3(14), 225–236 (2001)
14. Cohen, L.D., Kimmel, R.: Global minimum for active
contour models: A minimal path approach. International
journal of computer vision 24(1), 57–78 (1997)

15. Dashtbozorg, B., Mendonça, A.M., Campilho, A.: An automatic graph-based approach for artery/vein classification in retinal images. IEEE TIP 23(3), 1073–83 (2014)
16. De, J., Cheng, L., Zhang, X., Lin, F., Li, H., Ong, K.H.,
Yu, W., Yu, Y., Ahmed, S.: A graph-theoretical approach
for tracing filamentary structures in neuronal and retinal
images. IEEE transactions on medical imaging 35(1),
257–272 (2016)
17. Deschamps, T., Cohen, L.D.: Fast extraction of minimal
paths in 3d images and applications to virtual endoscopy.
Medical image analysis 5(4), 281–299 (2001)
18. Dijkstra, E.W.: A note on two problems in connexion
with graphs. Numerische mathematik 1(1), 269–271
(1959)
19. Duits & Janssen, Hannink, J., Sanguinetti, G.: Locally
adaptive frames in the roto-translation group and their
applications in medical imaging. Journal of Mathematical
Imaging and Vision pp. 1–36 (2016)
20. Duits, R., Boscain, U., Rossi, F., Sachkov, Y.: Association
fields via cuspless sub-Riemannian geodesics in SE(2).
Journal of Mathematical Imaging and Vision 49(2), 384–
417 (2013). DOI 10.1007/s10851-013-0475-y. URL http:
//dx.doi.org/10.1007/s10851-013-0475-y
21. Duits, R., Duits, M., Almsick, M., Haar Romeny, B.: Invertible orientation scores as an application of generalized wavelet theory. PRIA 17(1), 42–75 (2007). DOI
10.1134/S1054661807010063
22. Duits, R., Franken, E.: Left-invariant parabolic evolutions on SE(2) and contour enhancement via invertible
orientation scores part I: Linear left-invariant diffusion
equations on SE(2). Quarterly of Applied Mathematics
68(2), 255–292 (2010)
23. Duits, R., Franken, E.: Left-invariant diffusions on the
space of positions and orientations and their application
to crossing-preserving smoothing of hardi images. International Journal of Computer Vision 92(3), 231–264
(2011)
24. Duits, R., Meesters, S.P., Mirebeau, J.M., Portegies,
J.M.: Optimal paths for variants of the 2d and 3d reedsshepp car with applications in image analysis. arXiv
preprint arXiv:1612.06137 (2016)
25. Eppenhof, K., Bekkers, E., Berendschot, T.T., Pluim,
J.P., ter Haar Romeny, B.M.: Retinal artery/vein classifcation via graph cut optimization. In: E. Trucco,
X. Chen, M.K. Garvin, J.J. Liu, X.Y. Frank (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. 121–128. Iowa Research Online (2015)
26. Estrada, R., Allingham, M.J., Mettu, P.S., Cousins,
S.W., Tomasi, C., Farsiu, S.: Retinal artery-vein classification via topology estimation. IEEE transactions on
medical imaging 34(12), 2518–2534 (2015)
27. Favali, M., Abbasi-Sureshjani, S., Romeny, B.t.H., Sarti,
A.: Analysis of vessel connectivities in retinal images by
cortically inspired spectral clustering. Journal of Mathematical Imaging and Vision 56(1), 158–172 (2016)
28. Feragen, A., Petersen, J., Owen, M., Lo, P., Thomsen, L.H., Wille, M.M.W., Dirksen, A., de Bruijne, M.:
Geodesic atlas-based labeling of anatomical trees: Application and evaluation on airways extracted from ct. IEEE
transactions on medical imaging 34(6), 1212–1226 (2015)
29. Folland, G.: A fundamental solution for a subelliptic operator. Bulletin of the American Mathematical Society
79(2), 373–376 (1973)
30. Hannink, J., Duits, R., Bekkers, E.: Crossing-preserving
multi-scale vesselness. In: P. Golland, N. Hata, C. Barillot, J. Hornegger, R. Howe (eds.) MICCAI 2014, LNCS,

18

31.

32.

33.

34.

35.

36.

37.

38.

39.

40.

41.

42.

43.

44.

45.

46.

Erik J. Bekkers et al.
vol. 8674, pp. 603–610. Springer International Publishing (2014). DOI 10.1007/978-3-319-10470-6\ 75. URL
http://dx.doi.org/10.1007/978-3-319-10470-6_75
Janssen, M., Duits, R., Breeuwer, M.: Invertible orientation scores of 3D images. In: J.F. Aujol, M. Nikolova,
N. Papadakis (eds.) SSVM, LNCS, pp. 563–575. Springer
International Publishing (2015, Corrected version on
arXiv:1505.07690)
Jbabdi, S., Bellec, P., Toro, R., Daunizeau, J., PélégriniIssac, M., Benali, H.: Accurate anisotropic fast marching for diffusion-based geodesic tractography. Journal of
Biomedical Imaging 2008, 2 (2008)
Kaplan, A.: Fundamental solutions for a class of hypoelliptic pde generated by composition of quadratic
forms. Transactions of the American Mathematical Society 258(1), 147–153 (1980)
Kaul, V., Yezzi, A., Tsai, Y.: Detecting curves with unknown endpoints and arbitrary topology using minimal
paths. IEEE Transactions on Pattern Analysis and Machine Intelligence 34(10), 1952–1965 (2012)
Korányi, A.: Kelvin transforms and harmonic polynomials on the heisenberg group. Journal of Functional Analysis 49(2), 177–185 (1982)
Law, M.W., Chung, A.C.: Three dimensional curvilinear structure detection using optimally oriented flux. In:
European conference on computer vision, pp. 368–382.
Springer (2008)
Leontidis, G., Al-Diri, B., Hunter, A.: Exploiting the retinal vascular geometry in identifying the progression to
diabetic retinopathy using penalized logistic regression
and random forests (2016)
Li, H., Yezzi, A.: Vessels as 4-d curves: Global minimal 4D paths to extract 3-D tubular surfaces and centerlines.
IEEE Transactions on Medical Imaging 26(9), 1213–1223
(2007)
Lo, P., Sporring, J., Ashraf, H., Pedersen, J.J., de Bruijne, M.: Vessel-guided airway tree segmentation: A voxel
classification approach. Medical image analysis 14(4),
527–538 (2010)
Mashtakov, A.P., Duits, R.: A cortical based model for
contour completion on the retinal sphere. Program Systems: Theory and Applications 7(4), 231–247 (2016)
Mirebeau, J.M.: Anisotropic fast-marching on cartesian
grids using lattice basis reduction. SIAM Journal on
Numerical Analysis 52(4), 1573–1599 (2014). DOI
10.1137/120861667
Mirebeau, J.M.: Fast Marching methods for Curvature
Penalized Shortest Paths (2017). URL https://hal.
archives-ouvertes.fr/hal-01538482. Working paper
or preprint
Nagel, A., Stein, E.M., Wainger, S.: Balls and metrics
defined by vector fields i: Basic properties. Acta Mathematica 155(1), 103–147 (1985)
Péchaud, M., Keriven, R., Peyré, G.: Extraction of tubular structures over an orientation domain. In: Computer
Vision and Pattern Recognition, 2009. CVPR 2009. IEEE
Conference on, pp. 336–342. IEEE (2009)
Petitot, J.: The neurogeometry of pinwheels as a subRiemannian contact structure. Journal of PhysiologyParis 97(23), 265 – 309 (2003). Neurogeometry and visual
perception
Portegies, J., Sanguinetti, G., Meesters, S., Duits, R.:
New Approximation of a Scale Space Kernel on SE(3)
and Applications in Neuroimaging.
In: J.F. Aujol,
M. Nikolova, N. Papadakis (eds.) SSVM, no. 9087 in
LNCS, pp. 40–52. Springer International Publishing
(2015)

47. Portegies, J.M., Duits, R.: New Exact and Numerical Solutions of the (Convection-)Diffusion Kernels on SE(3).
arXiv:1604.03843 [math] (2016). ArXiv: 1604.03843
48. Prandi, D., Remizov, A., Chertovskih, R., Boscain,
U., Gauthier, J.P.: Highly corrupted image inpainting through hypoelliptic diffusion.
arXiv preprint
arXiv:1502.07331 (2015)
49. Rossmann, W.: Lie groups: an introduction through linear groups, vol. 5. Oxford University Press on Demand
(2002)
50. Rothschild, L.P., Stein, E.M.: Hypoelliptic differential
operators and nilpotent groups. Acta Mathematica
137(1), 247–320 (1976)
51. Sachkov, Y.L.: Cut locus and optimal synthesis in the
sub-Riemannian problem on the group of motions of a
plane. ESAIM: Control, Optimisation and Calculus of
Variations 17(2), 293–321 (2011)
52. Sanguinetti, G., Bekkers, E., Duits, R., Janssen,
M., Mashtakov, A., Mirebeau, J.M.: Sub-Riemannian
fast marching in SE(2).
In: A. Pardo, J. Kittler (eds.) Progress in Pattern Recognition, Image
Analysis, Computer Vision, and Applications, Lecture
Notes in Computer Science, vol. 9423, pp. 366–374.
Springer International Publishing (2015). DOI 10.1007/
978-3-319-25751-8\ 44
53. Sarti, A., Citti, G.: The constitution of visual perceptual units in the functional architecture of v1. Journal of
Computational Neuroscience 38(2), 285–300 (2015)
54. Sethian, J.: Level Set Methods and Fast Marching
Methods: Evolving Interfaces in Computational Geometry, Fluid Mechanics, Computer Vision, and Materials Science. Cambridge Monographs on Applied
and Computational Mathematics. Cambridge University
Press (1999). URL https://books.google.nl/books?
id=ErpOoynE4dIC
55. Shang, Y., Deklerck, R., Nyssen, E., Markova, A., de Mey,
J., Yang, X., Sun, K.: Vascular active contour for vessel
tree segmentation. IEEE Transactions on Biomedical Engineering 58(4), 1023–1032 (2011)
56. Tsitsiklis, J.N.: Efficient algorithms for globally optimal
trajectories. Automatic Control, IEEE Transactions on
40(9), 1528–1538 (1995)
57. Türetken, E., Benmansour, F., Andres, B., Glowacki, P.,
Pfister, H., Fua, P.: Reconstructing curvilinear networks
using path classifiers and integer programming. IEEE
transactions on pattern analysis and machine intelligence
38(12), 2515–2530 (2016)
58. Zhang, J., Dashtbozorg, B., Bekkers, E., Pluim, J., Duits,
R., ter Haar Romeny, B.: Robust retinal vessel segmentation via locally adaptive derivative frames in orientation
scores. IEEE Transactions on Medical Imaging PP(99),
2631–2644 (2016). DOI 10.1109/TMI.2016.2587062
59. Zhang, J., Duits, R., Sanguinetti, G., ter Haar Romeny,
B.M.: Numerical approaches for linear left-invariant diffusions on SE(2), their comparison to exact solutions, and
their applications in retinal imaging. Numerical Mathematics: Theory, Methods and Applications 9(01), 1–50
(2016)

