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

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

Erik J. Bekkers E.J. Bekkers J.M. Portegies Centre for Analysis, Scientific computing and Applications (CASA), Eindhoven University of Technology, the Netherlands
33email: 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
Da Chen E.J. Bekkers J.M. Portegies Centre for Analysis, Scientific computing and Applications (CASA), Eindhoven University of Technology, the Netherlands
33email: 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
Jorg M. Portegies E.J. Bekkers J.M. Portegies Centre for Analysis, Scientific computing and Applications (CASA), Eindhoven University of Technology, the Netherlands
33email: 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
###### 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 and . In our distance approximations we consider homogeneous norms on nilpotent groups that locally approximate , and which are obtained via the exponential and logarithmic map on . 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 . 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 and .

###### Keywords:
Sub-Riemanian geometry Roto-translation group SE(2) SE(3) Nilpotent approximation Geodesic vessel tracking Perceptual grouping

## 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 and . Additionally, we extend the perceptual grouping algorithm Cohen2001multiple () 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 , see Fig. 1 for an illustration.

### 1.1 Nilpotent Approximation

The sub-Riemannian distances on 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 sub-Riemannian 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 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 . The obtained Taylor development of the group product and associated left-invariant vector fields gives rise to a local approximation of the (sub-Riemannian) flows on in the sense of Rothschild and Stein rothschild-stein ().

We then define a norm on based on the Folland-Kaplan-Korányi gauge, which is known for its relation to the fundamental solution of the sub-Laplacian on the Heisenberg group folland1973fundamental (); kaplan1980fundamental (); koranyi1982kelvin (). We reason that the Folland-Kaplan-Korányi provides an accurate approximation to the fundamental solution on as well, as it provides the exact fundamental solution on the Heisenberg type approximation . As such, we provide an approach to approximating the heat kernel and fundamental solution of the sub-Laplacian on , as an alternative to the works Duits2010 (); portegies2016arxiv (); citti2006cortical ().

The distance associated with the Folland-Kaplan-Korányi type norm on is locally equivalent to the sub-Riemannian distance on , as was formally proven in full generality in the seminal work by Nagel, Stein and Wainger nagel-stein-wainger (). In this paper, we show by qualitative and quantitative comparison that the norm on indeed provides a sharp local approximation of the sub-Riemannian distances on .

### 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 Leontidis2016 ()), but also allows for higher level biomarker research via statistics on tree structures feragen2015geodesic (). Topological knowledge of vessel trees is also essential in determining artery/vein classification problems estrada2015retinal (); Eppenhof2015 (); Dashtbozorg2014 (). Finally, in many medical applications involving vessel analysis, including topological tree reconstruction, distances between local orientations play a crucial role de2016graph (); turetken2016reconstructing (); lo2010vessel (); shang2011vascular (); Abbasi-Sureshjani2017 (); favali2016analysis (). The approximate sub-Riemannian 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 Petitot2003 (); citti2006cortical (); Sarti2015 (); yuriSE2FINAL (); mashtakov2016cortical (); Duits2013 (); boscain2014 (); Bekkers2015SIIMS (); prandi2015highly (); favali2016analysis (). 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 Cohen2001multiple (). 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 benmansour2009fast (); kaul2012detecting (); chen2016vessel (), as is done also in this paper, see Fig. 2.

In Cohen2001multiple () 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. Bougleux2008anisotropic (). In recent work Chen2017 (), it was further extended for the grouping of closed contours for an a-priori . There, a (sub-)Finsler metric on position orientation space was used, similar to the sub-Riemannian metric used in this paper. As in Bougleux2008anisotropic () and Chen2017 (), we use the main algorithm of Cohen2001multiple () 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 is preferred over the use of (full) Riemannian distances on , which is in turn preferred over grouping with distances on . Furthermore, the analytic approximations allow for fast perceptual grouping with competitive performance compared to data-adaptive sub-Riemannian distances computed via fast marching.

### 1.3 Paper Outline

In Sec. 2 and Sec. 3 we derive approximations for sub-Riemannian distances in respectively and . 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 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 and in Subsec. 5.1, then on and in Subsec. 5.2. General conclusions are provided in Sec. 6.

## 2 Sub-Riemannian Distance and its Approximation in Se(2)

### 2.1 The Lie Group Se(2)

#### 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 is the semi-direct product of the group of planar translations and rotations , and its group product and inverse are respectively defined via:

 g⋅g′=(x,Rθ)⋅(x′,Rθ′)=(Rθx′+x,Rθ+θ′),g−1=(−R−1θx,R−1θ), (1)

with group elements . The group acts on the (coupled) space of positions and orientations via

 g⋅(x′,θ′)=(Rθx′+x,θ+θ′).

Since , we can uniquely identify the roto-translation group with the space of positions and orientations .

#### 2.1.2 The Lie Algebra, Exponential Map and Commutators

The Lie algebra associated with is the real vector space together with a bilinear operator called the Lie bracket (which we define below in Eq. (4)). The generators of the Lie algebra are given by the differential frame 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,A3|g=(Lg)∗A3=−sinθ∂x|g+cosθ∂y∣∣g (3)

via the push-forward of left-multiplication, denoted by , and with .

The exponential map defines a mapping from a vector in the tangent space at to an element in the group by following an integral curve along the left-invariant vector field . The logarithmic map defines the mapping from group element to tangent vector at .

The Lie bracket for vector fields is defined as follows

 [X,Y]:=limt→0γ(t)−et2,withγ(t)=Exp(−tY)Exp(−tX)Exp(tY)Exp(tX). (4)

I.e., it describes the infinitesimal displacement by following a path moving forth and back in and 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 are

 [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 by measuring distances between two points in via the lengths of shortest horizontal paths. A horizontal path is a curve with tangent vectors , where denotes the sub-Bundle of the full tangent bundle . Lengths of horizontal curves with are measured by the sub-Riemannian metric tensor111Due to the fact the the metric tensor is degenerate in the direction (tangent vectors are always contained within ) it is not possible to represent the metric tensor in a standard form as an invertible symmetric matrix. This is however possible when including an additional term after which the tensor becomes (anisotropic) Riemannian citti2006cortical (); Sanguinetti2015CIARP (). This Riemannian approximation converges to the sub-Riemannian tensor when (ThesisChenDa2016, , App. A) and (duits2016optimal, , Thm. 2).

 Gξ,C∣∣γ(t)(˙γ(t),˙γ(t)):=C(γ(t))2(|u1(t)|2+ξ|u2(t)|2), (6)

in which is an external cost which penalizes the curves to move through certain regions in , is a parameter which balances the penalty of motion in the angular and spatial directions and has dimensions [1/length], and and are the control parameters of the curve .

The sub-Riemannian distances between two points is then given by

 d0(g1,g2):=inf{∫10√Gξ,C∣∣γ(t)(˙γ(t),˙γ(t))dt}, (7)

where the infimum is taken over Lipschitz continuous curves with , , and . Note that due to the inclusion of an external cost function the distance is not strictly left-invariant, however, when substituting by in (7) we do have left-invariance (i.e., then ).

### 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 to the group

 (c1,c2,c3)↦(x,y,θ)=Exp(c1A1+c2A2+c3A3), (8)

with the basis vectors of given in (2), and with the canonical coordinates of the first kind given by

 c1=θ,c2={12θ(y+xcotθ2)ifθ≠0xifθ=0,c3={12θ(−x+ycotθ2)ifθ≠0yifθ=0. (9)

For two left-invariant vector fields and the Baker-Campbell-Hausdorff (BCH) formula (see e.g. rossmann2002lie ()) gives:

 Log(Exp(X)Exp(Y))=X+Y+12[X,Y]+112([X,[X,Y]]+[Y,[Y,X]])+O([⋅,[⋅,[⋅,⋅]]]), (10)

where denotes higher order nested brackets. Since the Lie algebra is not nilpotent it has non-vanishing Lie brackets of order (cf. the commutator relations in (5)) the BCH formula gives an infinite series of nested Lie brackets.

Here, we approximate the BCH formula as222Note that such approximations of the BCH formula were already introduced in (nagel-stein-wainger, , Thm. 2.22) in the general setting by Nagel, Stein and Wainger in 1985.

 Log(Exp(X)Exp(Y))≈X+Y+12[X,Y], (11)

by omitting the Lie brackets of order 2 (once nested brackets) and higher, as if our Lie algebra is nilpotent of step 2. Then, together with the commutator relations , , and again omitting Lie brackets of order 2 (i.e., setting ), the BCH formula defines a group product on the vector space of the canonical coordinates of the first kind via

 (x1,x2,x3)⋅(y1,y2,y3)=(x1+y1,x2+y2,x3+y3+12(x1y2−x2y1)). (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 local333With chosen close enough such that higher order terms in (10) can be neglected. approximation of the true group product given by (1). We denote this group by , with the 3 dimensional (nilpotent) Heisenberg group. Note that if and 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 defines a homogeneous Carnot group with respect to the dilations

 δs(c)=(sc1,sc2,s2c3). (13)

#### 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 of Eq. (7) we make use of the following homogenous norm on :

 ∥c∥ζ:=4√(|c1|2+|c2|2)2+ζ|c3|2, (14)

with constant a relative penalty for the non-horizontal part. For this norm coincides with the well-known Folland-Kaplan-Korányi gauge, which is a widely studied norm on Carnot groups due to its relation to fundamental solutions of sub-Laplacians bonfiglioli2007stratified (): Folland found that , with homogeneous dimensions , is (a constant multiple of) the fundamental solution of the canonical sub-Laplacian on the Heisenberg group folland1973fundamental (); Kaplan showed that this relation more generally holds for H-type (Carnot) groups kaplan1980fundamental (); Korányi derived many more of its properties in relation to harmonic analysis and potential theory koranyi1982kelvin ().

In relation to sub-Riemannian geometry on and its sub-Laplacian , we find that the fundamental solution of can be approximated by the (explicit) fundamental solution of the canonical sub-Laplacian , with Jacobian basis , on . This solution in fact coincides with one of the approximations of found by Duits & Franken Duits2010 (). There, the fundamental solution of was first approximated by relying on a contraction of to a 3-dimensional Heisenberg group (via dilations on the group ), and then derived the Gaussian estimates based on the homogeneous norm , i.e., , 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 (which is up to a constant given by ). In this context we thus obtain an estimate of the fundamental solution of by estimating it with , which is proportional to the exact fundamental solution of on our approximated group .

#### 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. Bellaiche1996 ()) and equivalence of homogeneous norms, there exists a constant such that

 \gothicc−1∥Log(g)∥ζ≤d0(e,g)≤\gothicc∥Log(g)∥ζ,

with defined by Eq. (9). The logarithmic norm is locally equivalent to the sub-Riemannian distance, which was proved in full generality in (nagel-stein-wainger, , Thm. 2 & 4). Via a scaling of the generators and we define the -isotropic norm

 ∥c∥ξ,ζ:=4√(|c1|2+|~c2|2)2+ζ|~c3|2=4√(|c1|2+ξ2|c2|2)2+ζξ2|c3|2, (15)

with and , and the given in (9). The norm closely approximates the sub-Riemannian distance for (no data-adaptivity) via

 d0(g,h)≈|Log(g−1h)|ξ,ζ,|Log(g)|ξ,ζ:=∥c∥ξ,ζ (16)

with the coordinates of the first kind obtained via (9). In view of the Folland-Kaplan-Korányi gauge setting in would be a sensible choice. We do observe however that gives an even sharper approximation, see Fig. 3 for a visual comparison to the sub-Riemannian distance and Appendix A for a quantitative comparison. The setting is used in all experiments on .

## 3 Sub-Riemannian Distance and its Approximation in Se(3)

In this section we extend the concepts of the previous section to the group 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 is the semi-direct product of the group of 3D translations and the group of 3D rotations . The group product and inverse for elements are defined by

 g⋅g′=(x,R)⋅(x′,R′)=(x+Rx′,RR′),g−1=(−R−1x,R−1). (17)

In the 3D case, we define the space of coupled positions and orientations as a Lie group quotient of :

 R3⋊S2:=SE(3)/(0×SO(2)).

The group action of onto is defined by

 g⋅(y,n)=(x,R)⋅(y,n)=(x+Ry,Rn).

We can identify the element with group elements , where is any rotation matrix such that .

#### 3.1.2 The Lie Algebra, Exponential Map and Commutators

Analogously as in the case, we associate with the group the Lie algebra using the exponential and logarithmic maps. This is most easily done using an isomorphism with the corresponding matrix group:

 (x,Rγ,β,α)↔(Rγ,β,αxT01).

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

 X1 =⎛⎜ ⎜ ⎜⎝0001000000000000⎞⎟ ⎟ ⎟⎠,X2=⎛⎜ ⎜ ⎜⎝0000000100000000⎞⎟ ⎟ ⎟⎠, (18) X3 =⎛⎜ ⎜ ⎜⎝0000000000010000⎞⎟ ⎟ ⎟⎠,X4=⎛⎜ ⎜ ⎜⎝000000−1001000000⎞⎟ ⎟ ⎟⎠, X5 =⎛⎜ ⎜ ⎜⎝00100000−10000000⎞⎟ ⎟ ⎟⎠,X6=⎛⎜ ⎜ ⎜⎝0−100100000000000⎞⎟ ⎟ ⎟⎠,

and their equivalents in the tangent space of span the Lie algebra . Since it will be clear from the context if we are in the or 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 left-multiplication , but they depend on the choice of coordinates. In this paper we mostly rely on -Euler angles in the parameterization of , i.e.,

 Rγ,β,α=Rez,γRey,βRez,α, (19)

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

 A1|g=(cosαcosβcosγ−sinαsinγ)∂x+(sinαcosγ+cosαcosβsinγ)∂y−cosαsinβ∂zA2|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β∂zA4|g=cosαcotβ∂α+sinα∂β−cosαsinβ∂γ,A5|g=−sinαcotβ∂α+cosα∂β+sinαsinβ∂γA6|g=∂α, (20)

for .

###### Remark 1

A second coordinate chart is needed to cover the entire , for which for example -angles can be used, as is done in e.g. duits2011HARDI (), where also the expressions for the vector fields in this alternative coordinate chart are given. In fact, the basis elements of the Lie algebra correspond to partial derivatives with respect to the -angles, similar to the -case.

We can express each element in terms of the basis with coefficients . Furthermore, we define and , the spatial and rotational coefficients, respectively. We can make the exponential map and logarithmic map explicit using these coefficients. For a matrix of the form

 Ω:=⎛⎜⎝0−c6c5c60−c4−c5c40⎞⎟⎠, (21)

we obtain a rotation using the exponential map of matrices, i.e., . The relation between the spatial coefficients and is given by

 c(1)=(I−12Ω+q−2(1−q2cot(q2))(Ω)2)x, (22)

where and such that . Now

 LogSE(3)(g)=6∑i=1ci(g)Ai,and (23)

using the relations above.

### 3.2 Sub-Riemannian Geometry in SE(3)

In the SE(3) case, a horizontal path is a curve with tangent vectors , where is now the sub-bundle of full tangent bundle spanned by . In this case we have one spatial control and two ‘angular’ controls and , so that the sub-Riemannian metric tensor becomes:

 Gξ,C∣∣γ(t)(˙γ(t),˙γ(t)):=C(γ(t))2(ξ|u3(t)|2+|u4(t)|2+|u5(t)|2), (24)

The sub-Riemannian distance between two elements is still defined as in (7), but now the infimum is taken over Lipschitz continuous curves with , , and .

### 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 and not on the quotient , i.e., different choices for in the rotational part result in different values for the coefficients . Here, we choose the approach of portegies2015ssvm () and set such that expected symmetries are preserved. With that choice, the logarithm (23) gives for each a unique vector , on which we can put a norm:

 |LogSE(3)(g)|ξ,ζ:=||c||ξ,ζ:= (25) 4√(ξ2|c3|2+|c4|2+|c5|2)2+ζ(ξ2(|c1|2+|c2|2)+|c6|2),

where according to (23).

Also here, the Folland-Kaplan-Korányi-type norm can be used to approximate the fundamental solutions of the sub-Laplacian on SE(3). The norm with was for example used in duits2011HARDI () approximations of the heat kernel and the fundamental solution on , of which only recently exact solutions were found in portegies2016arxiv (). In the context of this paper, we can approximate the exact solutions of the sub-Laplacian on by , with homogeneous dimensions , as the exact solution of the sub-Laplacian on the approximation group . The group that locally approximates , is again obtained via a nilpotent step 2 approximation of the BCH formula, and is defined by the group product

 (x1,x2,x3,x4,x5,x6)⋅(y1,y2,y3,y4,y5,y6)=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝x1+y1+12(x5y3−x3y5)x2+y2+12(x3y4−x4y3)x3+y3x4+y4x5+y5x6+y6+12(x4y5−x5y4)⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠T, (26)

with 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 on via the norm (25). I.e.,

 d0(g,h)≈|LogSE(3)(g−1h)|ξ,ζ, (27)

and as such again obtain an approximation of the distance in the sense of Rothschild and Stein rothschild-stein (). Based on the quantitative comparison to the sub-Riemannian distances in Appendix A and the visualizations in Fig. 4 of the level sets we conclude that the approximated sub-Riemannian distance of (27) quite accurately approximates the true sub-Riemannian distance on . In our analysis we found that the logarithmic norm with 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 in Fig. 4 is given by the surface , for a specific choice , and with density . The color of each orientation on the glyph is defined by the RGB color .

## 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 benmansour2009fast (), 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 Cohen2001multiple (), and which was later adapted for perceptual grouping based on anisotropic distances Bougleux2008anisotropic (). In recent work Chen2017 (), the perceptual grouping algorithm was extended for the grouping of closed contours for an a-priori specified . Like in Bougleux2008anisotropic () and Chen2017 (), we use the main algorithm of Cohen2001multiple () 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 of points of interest in which the edges are true connections (represented by geodesics) between points. Following the terminology of deschamps2001fast (); benmansour2009fast (); chen2016vessel () 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 is 2 at most). The final graph thus only contains sets of non-bifurcating 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 Cohen2001multiple () (and also Bougleux2008anisotropic ()) 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 .

In summary, our changes relative to the works Cohen2001multiple (); Bougleux2008anisotropic (); Chen2017 (), is that we

• keep the choice for distance open. In our experiments the distances will be mainly based on sub-Riemannian geometry in .

• 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 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. Chen2017 () using a (sub-)Finsler metric (based on the Euler elastica model) on position-orientation space. 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 cohen1997global (). Let be an arbitrary source point in a domain of interest, let be a metric tensor defined on the tangent space at , and let

 U(g):=d(g0,g)=infγ∈S(g0,g)∫10√G|γ(t)(˙γ(t),˙γ(t))dt (28)

its associated distance map, where the infimum is taken over the set of Lipschitz continuous curves with , , and with . Then the distance map is the unique viscosity solution of the eikonal equation

 {√G(∇GU(g),∇GU(g))=1,U(g0)=0,↔{∥∇GU(g)∥G=1,U(g0)=0, (29)

with the intrinsic gradient with inverse metric and the differential of , and the norm with respect to the metric tensor. In the standard (data-adaptive) Euclidean case with , , , , and with the eikonal equation is given by .

The fast marching algorithm efficiently solves the eikonal equation in a one pass algorithm. It computes the values of in increasing order (starting with ) based on the Bellman principle of optimality, in a manner very similar to the Dijkstra algorithm for shortest paths on graphs dijkstra1959note (). The minimal geodesic connecting with is then obtained via a gradient descent on from back to the origin , i.e., solving the ODE

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

For details on the fast marching algorithm on isotropic manifolds we refer to Tsitsiklis1995 (); Sethian1999 (), to Mirebeau2014 (); Jbabdi2008 () for anisotropic fast marching, and to Sanguinetti2015CIARP () and duits2016optimal () for fast marching in sub-Riemannian manifolds in and respectively.

### 4.3 Generating Key Points

The key point method is based on keeping track of a spatial arc-length map (in which the spatial lengths of the minimizing geodesics defining are stored), and stop as soon as a certain distance threshold is passed deschamps2001fast (). The rationale behind this algorithm is that among all points with equal geodesic distance values , the points reached by geodesics that best follow the data (paths along which is low) have maximum spatial distance . Such a point maximizing spatial distance in a given level set in is called a key point. The fast marching algorithm is highly suited for keeping track of a spatial arc-length map , in addition to , 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 deschamps2001fast ().

In summary, a key point is detected as follows. The spatial arc-length map is defined as

 Ul(g):=l(γg0,g), (30)

with