Nilpotent Approximations of SubRiemannian Distances for Fast Perceptual Grouping of Blood Vessels in 2D and 3D
Abstract
We propose an efficient approach for the grouping of local orientations (points on vessels) via nilpotent approximations of subRiemannian distances in the 2D and 3D rototranslation 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 subRiemannian distances, and we discuss their relations to the fundamental solution of the subLaplacian 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) subRiemannian geometry is essential in achieving top performance and 2) that grouping via the fast analytic approximations performs almost equally, or better, than dataadaptive fast marching approaches on and .
Keywords:
SubRiemanian geometry Rototranslation group SE(2) SE(3) Nilpotent approximation Geodesic vessel tracking Perceptual grouping∎
1 Introduction
In this paper we derive analytic formulas for approximations of subRiemannian 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 subRiemannian distances on , see Fig. 1 for an illustration.
1.1 Nilpotent Approximation
The subRiemannian 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 subRiemannian geometry in nilpotent Lie groups, in which such homogenous norms provide exact fundamental solutions to subLaplacians.
The vectors obtained by the logarithmic map, expressed in a leftinvariant basis, are the socalled 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 BakerCampbellHausdorf (BCH) formula) a global isomorphism to the group. In our setting we have to truncate the commutator series in the BCH formula due to nonvanishing (higherorder) commutators, yielding a corresponding Heisenberg type approximation which we denote with . The obtained Taylor development of the group product and associated leftinvariant vector fields gives rise to a local approximation of the (subRiemannian) flows on in the sense of Rothschild and Stein rothschildstein ().
We then define a norm on based on the FollandKaplanKorányi gauge, which is known for its relation to the fundamental solution of the subLaplacian on the Heisenberg group folland1973fundamental (); kaplan1980fundamental (); koranyi1982kelvin (). We reason that the FollandKaplanKorá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 subLaplacian on , as an alternative to the works Duits2010 (); portegies2016arxiv (); citti2006cortical ().
The distance associated with the FollandKaplanKorányi type norm on is locally equivalent to the subRiemannian distance on , as was formally proven in full generality in the seminal work by Nagel, Stein and Wainger nagelsteinwainger (). In this paper, we show by qualitative and quantitative comparison that the norm on indeed provides a sharp local approximation of the subRiemannian 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 (); AbbasiSureshjani2017 (); favali2016analysis (). 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.
SubRiemannian 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 subRiemannian 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 apriori . There, a (sub)Finsler metric on position orientation space was used, similar to the subRiemannian 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 subRiemannian 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 dataadaptive subRiemannian distances computed via fast marching.
1.3 Paper Outline
In Sec. 2 and Sec. 3 we derive approximations for subRiemannian distances in respectively and . There, for each Lie group we first provide the preliminaries, then define the subRiemannian 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 SubRiemannian Distance and its Approximation in
2.1 The Lie Group
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 semidirect product of the group of planar translations and rotations , and its group product and inverse are respectively defined via:
(1) 
with group elements . The group acts on the (coupled) space of positions and orientations via
Since , we can uniquely identify the rototranslation 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
(2) 
which define corresponding leftinvariant vector fields
(3) 
via the pushforward of leftmultiplication, 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 leftinvariant 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
(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 nonzero commutators of are
(5) 
2.2 SubRiemannian Geometry in SE(2)
We consider a subRiemannian 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 subBundle of the full tangent bundle . Lengths of horizontal curves with are measured by the subRiemannian metric tensor^{1}^{1}1Due 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 subRiemannian tensor when (ThesisChenDa2016, , App. A) and (duits2016optimal, , Thm. 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 subRiemannian distances between two points is then given by
(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 leftinvariant, however, when substituting by in (7) we do have leftinvariance (i.e., then ).
2.3 A Nilpotent Approximation of
2.3.1 A Local Approximation via the BakerCampbellHausdorff Formula
Consider the exponential map from Lie algebra to the group
(8) 
with the basis vectors of given in (2), and with the canonical coordinates of the first kind given by
(9) 
For two leftinvariant vector fields and the BakerCampbellHausdorff (BCH) formula (see e.g. rossmann2002lie ()) gives:
(10) 
where denotes higher order nested brackets. Since the Lie algebra is not nilpotent it has nonvanishing 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 as^{2}^{2}2Note that such approximations of the BCH formula were already introduced in (nagelsteinwainger, , Thm. 2.22) in the general setting by Nagel, Stein and Wainger in 1985.
(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
(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 local^{3}^{3}3With 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 step2 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
(13) 
2.3.2 Homogeneous Norms on and the Fundamental Solution of the subLaplacian
In our approximation of the subRiemannian distance of Eq. (7) we make use of the following homogenous norm on :
(14) 
with constant a relative penalty for the nonhorizontal part. For this norm coincides with the wellknown FollandKaplanKorányi gauge, which is a widely studied norm on Carnot groups due to its relation to fundamental solutions of subLaplacians bonfiglioli2007stratified (): Folland found that , with homogeneous dimensions , is (a constant multiple of) the fundamental solution of the canonical subLaplacian on the Heisenberg group folland1973fundamental (); Kaplan showed that this relation more generally holds for Htype (Carnot) groups kaplan1980fundamental (); Korányi derived many more of its properties in relation to harmonic analysis and potential theory koranyi1982kelvin ().
In relation to subRiemannian geometry on and its subLaplacian , we find that the fundamental solution of can be approximated by the (explicit) fundamental solution of the canonical subLaplacian , 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 3dimensional 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 subRiemannian distance approximations we found that even sharper estimates could be obtained by relying on the explicit formula for the fundamental solution of the (Kohn) subLaplacian 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 subRiemannian distance
Finally we arrive at the subRiemannian distance approximations. By the BallBox theorem (see e.g. Bellaiche1996 ()) and equivalence of homogeneous norms, there exists a constant such that
with defined by Eq. (9). The logarithmic norm is locally equivalent to the subRiemannian distance, which was proved in full generality in (nagelsteinwainger, , Thm. 2 & 4). Via a scaling of the generators and we define the isotropic norm
(15) 
with and , and the given in (9). The norm closely approximates the subRiemannian distance for (no dataadaptivity) via
(16) 
with the coordinates of the first kind obtained via (9). In view of the FollandKaplanKorá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 subRiemannian distance and Appendix A for a quantitative comparison. The setting is used in all experiments on .
3 SubRiemannian Distance and its Approximation in
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 subRiemannian distance, which allows us to do perceptual grouping in 3D images as well.
3.1 The Lie Group
3.1.1
The Lie group is the semidirect product of the group of 3D translations and the group of 3D rotations . The group product and inverse for elements are defined by
(17) 
In the 3D case, we define the space of coupled positions and orientations as a Lie group quotient of :
The group action of onto is defined by
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:
A basis for the corresponding matrix Liealgebra is given by
(18)  
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 leftinvariant vector fields are again obtained using the pushforward of the leftmultiplication , but they depend on the choice of coordinates. In this paper we mostly rely on Euler angles in the parameterization of , i.e.,
(19) 
with a rotation with angle around . Then, the leftinvariant vector fields are given by
(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
(21) 
we obtain a rotation using the exponential map of matrices, i.e., . The relation between the spatial coefficients and is given by
(22) 
where and such that . Now
(23)  
using the relations above.
3.2 SubRiemannian Geometry in SE(3)
In the SE(3) case, a horizontal path is a curve with tangent vectors , where is now the subbundle of full tangent bundle spanned by . In this case we have one spatial control and two ‘angular’ controls and , so that the subRiemannian metric tensor becomes:
(24) 
The subRiemannian 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 of and the Approximated subRiemannian Distance
It is important to realize that the logarithmic map is only welldefined 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:
(25)  
where according to (23).
Also here, the FollandKaplanKorányitype norm can be used to approximate the fundamental solutions of the subLaplacian 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 subLaplacian on by , with homogeneous dimensions , as the exact solution of the subLaplacian 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
(26) 
with coordinates of the first kind given by the logarithmic map (22). This new group is a freenilpotent group of rank 3 and step 2.
We approximate the subRiemannian distance on via the norm (25). I.e.,
(27) 
and as such again obtain an approximation of the distance in the sense of Rothschild and Stein rothschildstein (). Based on the quantitative comparison to the subRiemannian distances 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 subRiemannian 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 apriori 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 nonbifurcating vessel segments. The graph is build up by inserting onebyone 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 apriori 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 subRiemannian 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 subgraph.

group crossing lines without prespecifying the number of groups.
In particular, it is the use of a subRiemannian 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 positionorientation 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 apriori 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 nonclosed contours. Here, we focus on the grouping of nonclosed 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
(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
(29) 
with the intrinsic gradient with inverse metric and the differential of , and the norm with respect to the metric tensor. In the standard (dataadaptive) 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
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 subRiemannian manifolds in and respectively.
4.3 Generating Key Points
The key point method is based on keeping track of a spatial arclength 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 arclength 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 arclength map is defined as
(30) 
with