Geodesic Tracking via New Data-driven Connections of Cartan Type for Vascular Tree Tracking

We introduce a data-driven version of the plus Cartan connection on the homogeneous space $\mathbb{M}_2$ of 2D positions and orientations. We formulate a theorem that describes all shortest and straight curves (parallel velocity and parallel momentum, respectively) with respect to this new data-driven connection and corresponding Riemannian manifold. Then we use these shortest curves for geodesic tracking of complex vasculature in multi-orientation image representations defined on $\mathbb{M}_{2}$. The data-driven Cartan connection characterizes the Hamiltonian flow of all geodesics. It also allows for improved adaptation to curvature and misalignment of the (lifted) vessel structure that we track via globally optimal geodesics. We compute these geodesics numerically via steepest descent on distance maps on $\mathbb{M}_2$ that we compute by a new modified anisotropic fast-marching method. Our experiments range from tracking single blood vessels with fixed endpoints to tracking complete vascular trees in retinal images. Single vessel tracking is performed in a single run in the multi-orientation image representation, where we project the resulting geodesics back onto the underlying image. The complete vascular tree tracking requires only two runs and avoids prior segmentation, placement of extra anchor points, and dynamic switching between geodesic models. Altogether we provide a geodesic tracking method using a single, flexible, transparent, data-driven geodesic model providing globally optimal curves which correctly follow highly complex vascular structures in retinal images. All experiments in this article can be reproduced via documented Mathematica notebooks available at GitHub (https://github.com/NickyvdBerg/DataDrivenTracking).


Introduction
Retinal images are often used to examine the vascular system with optical scanning devices that image the vasculature in the retina noninvasively.The vasculature in the eye is known to be typically representative of the vasculature throughout the body.This allows doctors to monitor the circulatory system and aids in the diagnosis of different kinds of diseases like diabetes, hypertension [2][3][4] and Alzheimer's disease [5].Typically, high levels of tortuosity in the vasculature are a biomarker for such diseases [6][7][8][9].Successful automatic vessel tracking detects complex vasculature and aids the effective diagnosis of such diseases.Here, geometric models come into play via geodesic tracking methods where geodesics are the shortest paths that follow the biological blood vessels.They help in tracking and subsequent analysis of the vascular tree in the retina originating from the optic nerve [10][11][12][13].
Geodesic tracking has been extensively studied where many prevalent approaches perform the tracking in the standard 2D image domain [14][15][16][17][18][19].For many methods in this category, calculating the geodesics in R 2 leads to certain difficulties in accurately following the blood vessel.For example, one common difficulty is the inaccurate tracking of crossing structures and bifurcations.This has motivated methods that aim to lift the image function to higher dimensional spaces.For example, the space of positions and orientations [20,21] or radius-lifted spaces [22] where the lifting yields the benefit of disentangling seemingly complex crossing structures in the retinal images.In this article, we focus on the methods [13,[22][23][24] that perform the geodesic tracking in the 3D-space of positions and orientations M 2 .This is based on lifted images, or socalled orientation scores.The well-known benefit of this lifted approach is that lines involved in crossings are manifestly disentangled in M 2 .As visualized in Figure 1, the crossing circles in the image become disjoint spirals (cf.Fig. 1b) in the homogeneous space of positions and orientations.
However, practical considerations of working in (the domain M 2 of) orientation scores, like memory reduction and enabling low computation times, result in some undesirable effects.For example, using a limited number of orientations leads to imperfections in the computation of the orientation scores.Hence, some vessels can be assigned a near angular coordinate that may not reflect their true orientation, and therefore does not align with the vessel data correctly.We denote this problem as 'misalignment' (also referred to as deviation from horizontality [25]).Moreover, considering a limited number of orientations results in a sampling bias on orientations, and thereby the possibility of missing high curvature regions yielding poor curvature adaptation (cf.Fig. 2).
In this article, we provide a novel, data-driven tracking model that improves upon existing geodesic tracking methods.Our model demonstrates an improved curvature adaptation, reduces misalignment, and exhibits a high degree of geometric interpretability.
We will aim for a single geometric Finslerian model to deal with complex vasculature without requiring heavy preprocessing (e.g.placement of anchor points, pre-skeletonization) and associated extra parameters, and without suffering from the 'cusp problems' reported in [13,26,27].Fig. 2 Orientations sampling bias in geodesic tracking.Sampling bias can lead to wrong tracking results, and our new model will overcome this as we will show later in more detail (Fig. 9).
The cusp problem is tackled by creating an asymmetric Finslerian model1 (M 2 , F  ) extension of the data-driven Riemannian manifold similar to the much less data-driven techniques in [29].For a quick impression of such a 'cusp' in a spatially projected sub-Riemannian geodesic, see Fig. 19a in Appendix D. Clearly, cusps are undesirable for vascular tracking, and an asymmetric Finslerian version of the Riemannian manifold tackles this problem.Intuitively, cusps in spatial projections of sub-Riemannian geodesics arise sometimes as optimal paths of a 'Reed-Shepp' car (imagine a car driving along the geodesic track) [29,30] where the car was required to use its reverse gear to follow the optimal path.In the asymmetric Finslerian model we turn off the reverse gear of the car, while allowing for 'in-place rotations' see Fig. 19b.
Pre-processing techniques for geodesic tracking such as pre-skeletonization and iterative placement of anchor/key points are typically used in conjunction with Bézier curves [31] or splines on Lie-groups [31,32], but often require additional parameters and fine-tuning.Specifically, extensive use of anchor points implies that anchors get relatively close to each other, and then the choice of geometric model in between becomes increasingly less relevant (even non-data-driven sub-Riemannian distance approximations suffice as shown in the work by Bekkers et al. [32]).As a result, this reduces the geometric interpretability of the overall model.In this work, we therefore aim for a single geometrical model.Therefore, we will not use pre-processing, pre-skeletonization [22], multiple anchor points [33], and connectivity by perceptional grouping [16,32,34], even though these techniques are theoretically interesting and applicable.
In tracking an entire vascular tree, we limit the number of anchor points to at most one (which is computed without explicit manual supervision) and only use the boundary conditions for each vessel (see Fig. 3).Thanks to our new modified version of the anisotropic fast marching algorithm [35], we can now better adapt to curvature and spatial misalignment efficiently (see Fig. 4).We also address common pitfalls at Fig. 3 This tracking result (bottom) of a vascular tree in an optical image of the eye (top) is calculated with only two runs of the anisotropic fast marching algorithm.In the images, seeds, bifurcations, and tips are indicated by green, purple, and red points respectively.The white and cyan lines denote the tracking results obtained in the first and second run respectively.Details follow in the experimental section 6.
Fig. 4 Tracking result of the previous left-invariant model [29] (red), and the new data-driven left-invariant model (green).The tracking performed in the lifted space of positions and orientations is projected back onto the input 2D image.Our proposed model (in green) demonstrates a significantly improved accuracy in adapting to curvature of blood vessels in the optical image.complex overlapping structures, where one must impose additional constraints to avoid taking wrong exits in the tracking.The implementation of such constraints is easily accounted for in our model, as we will see.
Similarly to the plus-control variant of previous geometric curve optimization models [29] we replace symmetric, anisotropic, (sub-)Riemannian geodesic tracking by asymmetric, anisotropic Finslerian models that avoid cusps and allow for automatic placement of in-place rotations (recall Fig. 19 and see [29,36]).As a result, our model automatically accounts for bifurcations, thereby reducing the number of anchor points.
In summary, the main contributions of this article are: 1. We introduce a new geodesic tracking model that uses a crossing-preserving approach for tracking complex vasculatures in M 2 .Our model uses a new anisotropic fast marching algorithm to compute cusp-free data-driven geodesics.The induced geometric vessel tracking better adapts for vessel curvature and orientation sampling biases, compared to the previous model in [29].2. We mathematically analyze these solutions (the family of all geodesics) via our data-driven version of the plus Cartan connection (Section 3) that underlies the Hamiltonian flow as we will show in Theorem 1. 3. Finally, we demonstrate our method on highly challenging examples of retinal images with complex vasculature where adequate tracking results are obtained with only two runs of the proposed anisotropic fast marching algorithm.

Structure of the Article
In Section 2, we provide background of the geometrical tools underlying our method.We explain the space of positions and orientations M 2 , and why it is beneficial to apply tracking in this 3D space rather than in 2D position space.
In Sections 3 and 4, we describe our model.We begin in Section 3 by introducing a new data-driven Cartan connection ∇  associated with a data-driven left-invariant metric tensor field G  .These geometric tools allow for curvature adaptation and correction of misalignments in existing geodesic tracking algorithms in M 2 .
In Section 4 we use the data-driven Cartan connections and data-driven left-invariant metric tensor fields.We present our main theoretical result in Theorem1 where we characterize 'straight curves' and 'shortest curves' in data-driven left-invariant Riemannian manifolds on a finite-dimensional Lie group , analyze the Hamiltonian flow of all geodesics together, provide the geodesic backtracking formula of the new geodesic tracking model, address the symmetries of the geodesics and connections of the new model.
Then in Section 5 we employ the geometrical models and tools and present a numerical algorithm to compute the distance map for the special case where the Lie group equals the roto-translation group  =  (2).Additionally, we explain how to compute the backtracking of geodesics from end to source point.We present a new version of the anisotropic fast marching algorithm [37] that applies to our new data-driven model.
In Section 6, we report an extensive experimental evaluation of geodesic tracking in retinal images from the annotated STAR dataset [38,39], and show that our new model allows for adequate geometric tracking of highly complex vasculatures.
Finally, in Section 7, we end with a brief discussion of future work and conclude.

Lifted space of positions and orientations M 2
We begin by introducing the lifted space of positions and orientations M 2 .As motivated earlier, working in this space allows for convenient ways to separate difficult crossing structures in R 2 .In this article, we specifically focus on the challenging problem of vessel tracking in retinal images having complex vasculature.

Definition 1
The space of two-dimensional positions and orientations M 2 is defined as a smooth manifold where  1 ≡ R/(2Z) using the identification n = (cos , sin ) ↔ . ( Elements in the homogeneous space are denoted by ordered pairs (x, ) ∈ R 2 ×  1 but to stress the semidirect product structure of the roto-translation group  (2 The space M 2 is a homogeneous space under the transitive action of roto-translations given by the following mapping:

Definition 2 (Lie group action on domain)
For each roto-translation g = (y,   ) ∈  (2) the mapping of  g : M 2 → M 2 is given by where   ∈  (2) is the matrix associated with a counterclockwise rotation with rotation angle  ∈  1 .
Clearly, the concatenation of two rigid body motions is again a rigid body motion and indeed one has for all g 1 , g 2 ∈  (2), i.e.  is a group representation.After setting a reference point p 0 = (0, 0) ∈ M 2 we have a 1-to-1 relation between the roto-translation (x,   ) mapping p 0 to p = (x, ) and the point in the homogeneous space p = (x, ) ∈ M 2 .In particular, p 0 is then identified with the unity element e := (0, ) 2. In short, we identify Under this identification, the product of two elements, say (x 1 ,  1 ), (x 2 ,  2 ) ∈ M 2 , in the space of positions (x 1 , x 2 ∈ R 2 ) and orientations ( 1 ,  2 ∈  1 ) is given by where henceforth   2 ∈  (2) denotes the counter-clockwise rotation matrix with rotation angle  2 ∈ R/(2Z).

Definition 3
Let  be a Banach space, then () denotes the space of bounded linear operators on .
As mentioned we lift the image from R 2 to  (2) to separate crossing structures in the corresponding orientation score (Fig. 1).Next we explain how this is precisely done.

Definition 4 (Orientation Scores)
The orientation score transform   : L 2 (R 2 ) → L 2 ( (2)) using anisotropic wavelet  maps an image  to an orientation score  =    .The orientation score is given by where the rotated and translated mother wavelets are obtained by the Lie group action U :  (2) → (L 2 (R 2 )) given by for all g = (x,   ) ∈  (2) and all y ∈ R 2 .Definition 4 inputs an image  ∈ R 2 and yields a function  ∈  (2).This is achieved by taking a convolution with a rotated wavelet filter, where the canonical/mother wavelet function  is rotated counter-clockwise with the angle , as we can see in (3).By varying  over all orientations R/(2Z), the image is lifted from R 2 to M 2 .We use the real part of the cake wavelets [24,40] depicted in Fig. 1c as then the space of orientation scores is naturally embedded in L 2 ( (2)) [40], and gives practically informative orientation scores [2].More information on the orientation score transform, its range, its invertibility, and the choice of wavelets  can be found in previous works [2,[40][41][42].
We can rotate and translate images via  ↦ → U g  .This corresponds to a left action on the orientation score: for all g ∈  (2).Left actions are defined as follows:
As shown in [40,Thm. 21], by construction of (3), orientation score processing must be left-invariant (i.e., equivariant with respect to left actions L g ) and not right-invariant.The key reason for this is the fundamental relation ( 4).This leftinvariance will also be crucial in our geometric tracking which needs to be left-invariant (and not right-invariant).Indeed rotating and translating the image should yield to an equally rotated and translated (lifted) tracking output curve.

Metric Tensor Fields and Finsler Functions
To calculate shortest paths in orientation scores, we need to establish local costs on tangents (velocities).We do so by assigning a metric tensor G p (•, •) to every point p = (x, ) in the lifted space of positions and orientations.It is beneficial to design this metric tensor depending on the specific application.Typically, this choice of the metric tensor field establishes the geometric model.Additionally, diagonalizing this tensor at every point p provides a local frame of reference.
First, we introduce the static frame denoted by {  ,   ,   }, induced by the coordinates , ,  for all points in M 2 .Its dual frame, for the cotangent bundle  * (M 2 ), is denoted by {d, d, d}, and can be used to express the metric tensor field G.
It is advantageous to use left-invariant vector fields for our application, since it guarantees that tracking results are equivariant to the group of roto-translations.More specifically, tracking is independent of the roto-translation of the image, meaning that tracking on a roto-translated image is identical to tracking on the original and roto-translating the result.

Definition 6 (Frame of Left-Invariant Vector Fields)
The frame of left-invariant vector fields (left-invariant frame) is obtained by a pushforward of the static frame at the origin p 0 .We define the pushforward ( g ) * :  h () →  gh () by After computations, we obtain the left-invariant vector fields The corresponding dual frame is given by {  } 3 =1 where   A  =    .A brief computation gives  1 = cos  d + sin  d,  2 = − sin  d + cos  d and  3 = d.
In addition to that, we define the left-invariant metric tensor field: for all p ∈ M 2 , all p ∈  p (M 2 ) and all g ∈  (2).

Remark 1 (Left-invariant Metric Tensor Field)
Let G denote a left-invariant metric tensor field on .Then there exists a unique constant matrix where ⊗ denotes the usual tensor product.
In the standard left-invariant model we restrict ourselves to the case    =      , and then G is diagonal with respect to the co-frame Often we do not want to work with symmetric Riemannian metric tensor fields (for instance to avoid cusps, and to ensure that fronts only move forward, see Fig. 5), and then we resort to the general Finsler geometry as done in [29] and [35].Essentially, this means that we replace the symmetric norm √︃ G|  ( ) ( (), ()) in the Riemannian distance/metric:
Next, we define some relevant geometrical tools associated to (7) and (8).

Remark 2 (The Space of Curves over which we optimize)
To adhere to standard conventions in Riemannian geometry we optimize over the space of piecewise continuously differentiable curves in M 2 (indexed by  > 0): In (7) we set  = 1, as there the choice of  is irrelevant by parameterization independence of the functional.

Remark 3 (Control Sets)
The control set in the tangent bundle  (M 2 ) is defined as with p ∈ M 2 and G the underlying Riemannian metric tensor field.The corresponding asymmetric control set in the tangent bundle  (M 2 ) is defined as with p = (x, n) ∈ M 2 and F the underlying Finslerian metric tensor field.In the limiting case where backward motions become prohibited as  ↓ 0 (i.e. the sub-Finslerian setting) we only get half of the Riemannian control sets

Cartan Connections
The theory of Cartan connections was developed by Élie Cartan.His viewpoint on differential geometry relies on moving frames of reference (repère mobile).The idea is to connect tangent spaces by group actions on homogeneous spaces.This geometric tool allows us to understand the geodesic flow associated to (7) and its data-driven extensions.
The homogeneous space that we will use for crossingpreserving 2D image processing is the homogeneous space of positions and orientations M 2 .Here, the pushforward ( g ) * of the left-multiplication connects  e (M 2 ) to  g (M 2 ) as it maps  e (M 2 ) (isometrically3) onto  g (M 2 ) and ( −1 g ) * , known as the Cartan-Ehresmann form, maps  g () back to  e ().
First, we introduce the general definition of Cartan connections, after which we also introduce the Cartan plus connection [44].In this article, we will introduce a datadriven version of the Cartan plus connection, leading to a generalization of the existing theory on shortest and straight curves in M 2 .
We use the following special case of a Cartan connection to define shortest and straight curves.Note that this Cartan connection is easily expressed in the left-invariant frame [22,[44][45][46][47][48].Then the Cartan plus connection is given by Remark 4 Note that the • symbol denotes the composition of functions such that for example For explicit coordinate expressions, see [44].Next, we explain how to read and compute (12).The covariant derivative ∇   of a vector field  with respect to a vector field  is again a vector field.Indeed the above formula gives so that it becomes clear where  and  typically enter in the open slots of the expression (12).Note that vector field A  in (13) is a differential operator applied to the smooth function  ∋ g ↦ →   g ( g ) ∈ R.

Remark 5
The connection ∇ [+] is called 'Cartan plus connection' as we add the two sums between the two large round brackets.In differential geometry one also has Cartan connections with a realvalued scalar factor in front of the second term, but this does not serve our applications [49].
Now that we explained Cartan connections, let us return to the core purpose of designing a geometric model such that projected geodesics follow the blood vessels.

Geodesic Tracking in the Space of Positions and Orientations M 2
Geodesic tracking methods aim to find the shortest paths following the underlying biological blood vessels in the retinal image.Such shortest paths are obtained by finding minimizing geodesics, which are defined as curves with the shortest length functionals.Typically, such length functionals are driven by a cost function that is small at locations of the blood vessels and high at all other places.Many different approaches to determine the minimizing geodesics have been proposed over the years, ranging from classical geodesic tracking in the image domain [15] to tracking in higher-dimensional homogeneous spaces [2,13,23,29].
As already explained, we lift the input image  from R 2 to M 2 using orientation scores in homogeneous spaces [41,42] (see Fig. 1).The tracking process involves computing a geodesic distance map in M 2 and then using steepest descent to find the shortest curve in the lifted space.Finally, we project the curve back onto the input image in R 2 to get the final tracking result, see the examples in Figure 4.
Over time, different models have been introduced that describe how the geodesics should behave.Imagine a car moving along such geodesics.Then the Reeds-Shepp car model [30] which describes the problem of shortest paths for cars between an initial and final point, and the Reeds-Shepp forward model [29] turns off the reverse gear of the car.In both cases the spatially projected geodesics (optimal paths) tend to follow blood vessels in medical images well.

Symmetric Reeds-Shepp Car Model
The left-invariant metric tensor field of the symmetric Reeds-Shepp car model, G, is given by the symmetric tensor field The anisotropy parameter  penalizes vectors with large sideways components.Note that the classical sub-Riemannian model corresponds to the limit  ↓ 0. For formal convergence results of the Riemannian model to the sub-Riemannian model see [29,Thm.2].In practice, choosing  = 0.1 usually provides a good enough approximation of the sub-Riemannian model for our purposes.The last parameter , a weighting parameter, influences the flexibility of the tracking.It either stimulates or discourages angular movement over spatial movement [29].
The cost function  : M 2 → [, 1], with  > 0, discourages movement at specific locations, e.g.outside vessel structures.In this article, the smooth costfunction (, , ) ↦ →  (, , ) is typically a version of the multi-scale crossing preserving vesselness map [50] explained in Appendix D. For an impression of what such a map  looks like see the 3D visualization in Figure 18.As we consider rather complex vasculatures it is often more intuitive to display their minimum projections over , see for example Figure 13.

Asymmetric Reeds-Shepp Car Model
Besides the symmetric version of the left-invariant metric tensor field of the Reeds-Shepp car model, an asymmetric version has been introduced in [29].The forward gear leftinvariant metric tensor field of this model is given by the asymmetric Finsler norm/function for all p = (x, n), p = ( x, n), with  − := min{0, }.Eq. ( 15) coincides with (8) with The parameters  and  and the cost function  have the same meaning as in the symmetric model.However, we consider an extra variable  ∈ (0, 1] in the asymmetric Reeds-Shepp car model.This parameter determines how strongly the model needs to adhere to the forward gear.Note that when  = 1, we find the symmetric Reeds-Shepp car model, and when  → 0, backward movement becomes prohibited.In that case, we move from cusps to change orientation to inplace rotations visualized (cf.Fig. 19 in Appendix D).These asymmetric Finslerian models are also highly beneficial in image segmentation as shown by Chen and Cohen [51].

Anisotropic Fast Marching
We provide here a brief overview of the partial differential equation (PDE) framework associated with geodesic distance maps, and of their numerical computation, see Section 5 for further details.We already mentioned that it is common to calculate minimizing geodesics in two steps; first calculating the geodesic distance map, then calculating the shortest curve using steepest descent.To get a first impression of how this looks like in practice, see Fig. 5.The geodesic distance map is characterized, in the PDE framework, as the viscosity solution of a static first-order Hamilton-Jacobi-Bellman equation, known as the Eikonal Equation.For numerically solving the Eikonal equation, it is discretized using for instance finite differences [13], leading to a coupled non-linear system of equations, which is typically solved using a front propagation method such as the fast marching algorithm (FMM).Classical references on the FMM include [18,52], anisotropic variants are presented in [35,37,53], and the details related to our new model will follow in Section 5.1.
The fast marching algorithm is a numerical method for solving the coupled system of equations discretizing the eikonal PDE.The algorithm proceeds in only one pass over the domain hence providing significant efficiency gains, but also requiring that the numerical scheme satisfies two conditions (monotonicity and causality), see [35,37,Def. 2.1].The proposed variant of this method uses Selling's algorithm [54] to calculate in a preliminary step a decomposition of the quadratic forms defining the dual metric.This dual metric suitably only involves positive weights and vectors with integer coordinates, see Proposition 2. These ingredients are used to devise an adaptive finite differences scheme, discretizing the anisotropic Eikonal PDE and obeying the required conditions, see Section 5.2.The eikonal PDE is solved via an given by ( 17)-does not, and moreover adapts for curvature in M 2 , cf.Fig. 7).Explicit formulas for ∇  and ∇ [+] will follow later in Table 2.
anisotropic Fast Marching algorithm, and its solution provides the desired distance map.Finally, the minimizing geodesic is calculated by solving an ordinary differential equation defined in terms of the distance map [35,37].
In previous studies of the Reeds-Shepp model and variants [29,35], the geodesic metric tensor matrix featured a block diagonal structure, which was exploited in the discretization.However, while working with data-driven metric tensor fields, this block format does not apply!Therefore we adapt the anisotropic fast marching algorithm to cope with the general setting.In this article, we will briefly discuss the changes that were necessary to solve data-driven metric tensor fields.Such data-driven geometric models (that we explain in the next section) give better tracking results than the previous model, as one can see in Fig. 5.In addition to that, using the anisotropic fast marching algorithm to calculate the geodesics, only a limited number of runs (one for a single vessel, Fig. 5, and only two for a full vasculature, Fig. 3) are needed to correctly track the vascular structures.

Flowchart and Overview of the Methodology
Before we dive into the details of our method, we provide a sequential flowchart of our methodology, and more information on where more details of each part will be addressed.

Data-Driven Metric & Data-Driven Cartan Connection
In multi-orientation image processing, it is beneficial (for vessel segmentation [38]) to rely on locally adaptive frames [55,56].However, the locally adaptive frames in M 2 typically require a stable selection of the principal eigenvector (eigenvector corresponding to the largest eigenvalue) of a symmetrized Hessian of the function  :  (2) → R. Recall that a Hessian is defined by a (dual) connection.Even if one uses the left Cartan connection, selecting the principal eigenvector can be locally unstable [56] and the largest eigenvalue may not be unique.For instance, if line structures are not locally present at all.Therefore, in this article, we take a slightly different approach by creating an unconditionally stable data-driven left-invariant metric tensor field.

Definition 10 (Data-Driven Left-Invariant Metric)
Let  be a Lie group.Then the metric tensor field G  is datadriven left invariant when it satisfies for all (g, g) ∈  () and all q ∈ : Recall that in our case of interest where  =  (2) and where  = W   is an orientation score of the image  , the equivariance relation (4) holds, so roto-translation of an image  ↦ → U g  is equivalent to roto-translation  ↦ → L g  of the score.
Consequently (as will follow in Theorem 1) if a metric tensor field is data-driven left invariant then a roto-translation U g  of the input image  yields a new geodesic   that is rotated and translated accordingly: Thus, Definition 10 is a valid constraint in our application as we want the vessel tracking along geodesics to be equivariant with respect to roto-translations.
By creating such a data-driven metric tensor field G  on our Lie group of interest  =  (2) ≡ M 2 , data-driven corrections are made for spatial and angular misalignment in existing models relying on the standard left-invariant frame [25,56,57].We will see that a better fitted metric tensor field G  has a significant impact on the tracking results for very tortuous vessels, as shown in Figure 5.For our case of interest M 2 , a reasonable choice that satisfies the constraint, and that we use in our experiments, is given by: , where G and F are given in ( 14) and ( 15), ( 17) Here the Hessian field  is defined in Lemma 4 in Appendix C, and ∥ • ∥ * the dual norm corresponding to the primal norm given by √︁ G p ( p, p) with  = 1 in ( 14).Parameter  > 0 regulates inclusion of data-driven 2nd order line-adaptation to the orientation score data , cf.Fig. 1.
Finally, the data-driven left-invariant metric tensor field relies on the usual Reeds-Shepp car models G respectively F with external smooth cost  (p) satisfying: computed from the orientation score , as explained in Appendix D. There we combine ideas on crossing-preserving vesselness maps from [29,38,50].
Remark 6 Within G and F in (17) we set  2 = 0.01 =  11 / 22 as relative costs for sideward motion, recall (14).Ideally we want this to be high, but as we will prove in the numerics section (Subsection 5.2), a spatial anisotropy of  2 = 0.01 still guarantees numerical accuracy.We follow [13,29] and set bending stiffness  2 =  11 = 0.01 and  33 = 1.
Proof See Lemma 6 in Appendix C.
Next, we list a few remarks that underpin and explain our specific choice of metric tensor field.

Remark 7
In geometric image analysis [58], eigenvectors of the Hessian typically provide a local coordinate frame along lines.In orientation scores, this is not different [25].In M 2 , Hessians  = ∇ [+], * d are not symmetric and we rely on a singular value decomposition via the dual norm in (17) which only relies on the symmetric product, see Remark 9.
Remark 8 Formally speaking, the (old) metric tensor fields G and asymmetric version F are also data-driven if it comes to scalar-adaptation via cost function , but as they do not adapt for any kind of directional data-adaptation (as illustrated in Fig. 6) we do not refer to them as 'the data-driven model'.  in (17) then boils down to a straightforward Euclidean norm: where For details of Hessians of functions on manifolds with a connection, see Appendix C. For now let us focus on the notion of data-driven left-invariant frames, where we improve upon the 'Locally Adaptive Derivatives (LADs)' in [38,56].

Definition 11 (Data-Driven Left-Invariant Frame)
Any data-driven metric tensor field G  can be diagonalized: and this defines the positively oriented data-driven left-

Remark 10 (Advantages of our data-driven metric and frame)
The local frame of reference {A   } depends on the image data, cf.Fig. 6.In fact, Eq. ( 20) is used to define the dual of the data-driven left-invariant frame via diagonalisation.This deviates from LADs in previous work [38,44].We now have the advantage of coercivity recall (18), independent of the orientation score data , which makes the tracking algorithms unconditionally stable.Furthermore, we now have another advantageous property over LADs, namely that A   = A  for  constant.∇  Table 1 Comparison of (notation of) current and previous work.Diagonalization is w.r.t.dual frame associated to the frames depicted in Fig. 6.
In order to calculate distances using the new data-driven metric tensor field, we need to introduce the data-driven Riemannian distance.

Definition 12 (Data-Driven Riemannian Distance)
The data-driven Riemannian distance  G  from a point p ∈ M 2 to a point q ∈ M 2 is given by where Γ 1 was defined in (9), and () := d d ().
Remark 12 Note that this distance can always be transformed to a quasi-distance when we are working in with the forward gear version of the model: Using the new data-driven metric frame, recall Def.11, we introduce the data-driven Cartan plus connection, which will be used to express 'short' and 'straight' curves in Section 4.

Definition 13 (Data-Driven Cartan Plus Connection)
The data-driven Cartan plus connection is given by Explicit coordinate expressions will follow in Lemma 1.
In Table 1, an overview of the notation used for the new concepts introduced in this work and concepts introduced in earlier work is given.
In Figure 7, the exponential curves and the control sets for both discussed Cartan connections, ∇ [+] and ∇  are visualized.In addition to that, the tracking results relying on different models are plotted.One sees that the data-driven Cartan connection better adapts for curvature leading to more accurate tracking results.
Cartan Plus Connection ∇ [+]  Data-driven Cartan Plus Connection ∇  Fig. 7 The advantage of using data-driven Cartan connections ∇  instead of the non-data-driven Cartan plus connection ∇ [+] .In grey, a shortest curve (geodesic) between two points in M 2 is visualized, along with its spatial projection.The left geodesic has parallel momentum w.r.t.∇  (cf.Thm 1) and the right w.r.t.∇ [+] [44, Thm.1].The new geodesic better adapts for curvature (and spatial alignment).This is also visible in the corresponding control sets (11) depicted by the white closed surfaces above at several green points on the geodesics.The red arrow indicates the principal direction of the local metric tensor (left: G, right: data-driven G  ).The control sets belonging to ∇ [+] are only aligned to the underlying structure in the spatial domain, whereas the control sets belonging to ∇  align with the appropriate curvature in the tangent space as well.In the bottom row, we depict exponential curves through the green points with a tangent in the principal direction (left of G, right of G  ).They are straight-curves of ∇ [+] (left) and ∇  (right), and 'steer' the geodesic tracking as we will show in Thm 1.

The New Geometric Tracking Model: Asymmetric Finsler Functions steered by Locally Adaptive Frames
We discuss a new data-driven version of the Cartan connection.This result applies to all Lie-groups  of finite dimension dim() = .Note that  (2) ≡ M 2 , but the result does not apply to all homogeneous spaces (like M  for  > 2).The notation used in this section is summarized in Table 2.
We consider a locally adaptive frame A    =1 with dual frame     =1 .This can be any well-defined frame that depends on the underlying data.The (data-driven) metric tensor field that is considered, is given by (20).The datadriven terms can adapt for curvature and deviation from horizontality where the direction of the left-invariant frame deviates from the underlying line structure.Relation between bases  () and  * () ⟨d  , LIV metric tensor field For specific choices of coefficients   , see (14).
Data-driven LIV metric tensor field that must satisfy (16).We choose to use (17) and write   :=    .For its asymmetric Finslerian extension F  see also (17). Hamiltonian Fundamental Symplectic form Structure Functions Data Driven Cartan Connection see also (31) Dual Data Driven Cartan Connection In previous works, the Cartan plus connection, which relies on the left-invariant frame, has been used to describe straight and shortest curves in Lie groups [49].However, this frame is not always adequate in multi-orientation image processing as it does not always align perfectly with the underlying line structures in the orientation scores (see Figure 6).To improve the tracking results, we, therefore, switch to using a data-driven Cartan connection associated with the data-driven metric tensor field G  given by (17).Let us first define what we mean by a 'data-driven Cartan connection'.

Definition 14
The data-driven Cartan connection and its corresponding dual are given by where ∇  *   := ∇  * (, ) and ∇    := ∇  (,  ).Remark 13 The relation between ∇ and its dual ∇ * is for all vector fields ,  and all covector fields  on , which may be interpreted as a product rule for the pairing between the vectors and co-vectors.In particular for  = A   ,  = A   and  =    we get c = − c   .
In the next lemma, we will express the data-driven Cartan connection and its corresponding dual explicitly in coordinates, which will provide us an expression on which we will build in the proof of our main theorem, Theorem 1.
Lemma 1 When expressing Eq. ( 28) and (29) more explicitly in data-driven left-invariant frame components (gauge frame components for short), one finds and for the dual connection where  =

Main Theorems
Our goal is to analyze and structure the Hamiltonian flow belonging to the new data-driven geometric model determined by a data-driven metric tensor field G  .For convenience, we restrict ourselves in our main theorem to the case where the homogeneous space equals a full finite-dimensional Lie group  as the base manifold.We are mainly interested in the case  =  (2) ≡ M 2 with  = 3 and with data-driven metric tensor field G  given by (17).
Before stating the main theoretical result, that generalizes [44, Thm.1&2] to data-driven metric tensor fields G  , we introduce the concepts of parallel momentum and parallel velocity.They are now determined by the data-driven Cartan connection ∇  and its dual ∇  * .Definition 15 (Parallel momentum) Let  : [0, 1] →  be a smooth curve.Then, the curve (•) has ∇  parallel momentum (•) Definition 16 (Parallel velocity) Let  : [0, 1] →  be a smooth curve.Then, the curve (•) has parallel velocity (•) Remark 14 Using the antisymmetry of the structure functions ( 26) and (31) in Lemma 1 we can rewrite Eq. ( 34) to Next, we formulate the main theoretical results where we generalize the main results [44, Thm.1&2] from the leftinvariant setting (, G) with connection ∇ [+] , to the new data-driven geometric models (, G  ) with connection ∇  .In more detail, the next theorem shows 1. how to compute globally optimal distance minimizers in a geometry that relies on data-driven left-invariant frames: These geodesics have parallel momentum w.r.t.connection ∇  (Def.15). 2. that the locally optimal straight curves are the straight curves w.r.t.connection ∇  : These curves have parallel velocity (i.e. are auto-parallel) w.r.t.∇  (Def.16).
Note that the equation for geodesics of the new data-driven model (M 2 , G  ) gives a wild expression in the left-invariant frame.In the fixed frame it is even worse.However, our new tool of the connection ∇  expressed in the locally adaptive frame    allows us to describe these geodesic equations (and the Hamiltonian flow) concisely and intuitively by the next theorem, using the tools listed in Table 2.
The shortest curve  : [0, 1] →  with (0) = g and (1) = g 0 may be computed by steepest descent backtracking on the distance map  (g) =  G  (g, g 0 ) where Exp integrates the following vector field on : and where  is the viscosity solution of the eikonal PDE system 5 For the definition of pullback of a dual connection, see Remark 26 in Appendix B.
Remark 15 (Assumptions on point g in backtracking ( 37)) For the geodesic backtracking formulated above, we need a differentiable distance map along the path and a well-posed Hamiltonian flow (i.e. the mapping from (0) to () arising from (35) must have a non-vanishing Jacobian) along the path.If g would be a first Maxwell-point (i.e. two distinct geodesics meet for the first time at g) the distance map is not differentiable at g.If g would be a conjugate point (often limits of first Maxwell points [13]) then the Hamiltonian flow breaks down.In the latter case, local optimality is lost.
In the first case, global optimality is lost.Fortunately the viscosity property of the viscosity solution [67] of (38) kills non-optimal fronts [13] and one may resort to multi-valued backtracking via sub-gradient backtracking.
The next 3 remarks show how incredibly simple the Hamiltonian flow, the eikonal PDE, and the backtracking of geodesics become when expressed in the data-driven leftinvariant frame.

Remark 16
Eq. ( 33) is in gauge frame components simply: Remark 17 Eq.( 38) is in gauge frame components simply: Remark 18 Eq.( 37) is in gauge frame components simply: This explains the definition of () below (37).A more explicit integration formula for (37) can be obtained in a similar way as in [27,64] (where exact solutions are derived for  =  = 1) via momentum preservation laws.

Asymmetric Finsler Extension to Automatically Deal with Bifurcations
One can always decide to include a positive control variant, to avoid cusps in the spatial projection of geodesics in  =  (2).This is done by considering the geodesics in the asymmetric Finslerian manifold (M 2 , F  ), recall (17), rather than the geodesics in the Riemannian manifold (M 2 , G  ).
It is not too hard in practice: a slight adaptation of the eikonal PDE (taking the positive part of one momentum component) will guarantee that all optimal geodesic wavefronts propagate in the preferred forward direction around the source point, as can be observed in Figure 5 where the fronts initially move 'down-left' (and not 'up-right').
The asymmetric Finslerian model (M 2 , F  ) is still wellposed (controllable and piecewise regular geodesics) even if  ↓ 0. In fact, such asymmetric Finslerian geodesics are by construction piecewise concatenations of Riemannian geodesics and in-place rotations.These observations follow by a straightforward generalization of [29,Thm.1,2,4] to the data-driven setting, where the control set formulation of the geodesic distances, still applies: where Γ  was defined in (9).Moreover, the field of control sets given by p ↦ → B F  (p) recall (10) remains continuous when using F  or G  (instead of F or G), which directly follows by [29,Lemma 6] in conjunction with the important coercivity property (21).The nice thing in practice is that in-place rotations are automatically placed at optimal locations by the eikonal PDE system (solved by the anisotropic fast-marching algorithm that we discuss in the next section).It is not surprising that, when using a reasonable cost function  (see Figs. 12 & 13), these in-place rotations are automatically placed at bifurcations in complex vascular trees as can be observed in the upcoming Figure 16.

Numerical Solutions to the Eikonal PDE System: Extension of the Anisotropic Fast-Marching Algorithm that allows for all Left-Invariant Data-Driven Metrics
In this section we describe the computation of globally minimizing geodesics for the new models F  considered in this paper, whose fundamental ingredient is the numerical solution to an anisotropic eikonal PDE.The Reeds-Shepp forward optimal control model F , of which a variant F  is considered in this paper, has already been addressed numerically using a Semi-Lagrangian [29] and Eulerian [35] discretization of the corresponding eikonal PDE.Both works however take advantage of the fact that the original geodesic model F regards the tangent spaces to the physical R 2 and the angular  1 domains as orthogonal to each other.However, in our case of interest (with model F  given by ( 17)), we cannot expect such a block-matrix structure in the fixed coordinate system (, , ).
To overcome this problem, we describe below the extension of [35] to the adaptive frame setting considered here, where this orthogonality relation is lost; in contrast, [29] could not be generalized in an efficient manner.

Remark 19 (Convenience notations for the numerical section)
Throughout this section, we label the dependence w.r.t. the relaxation parameter  ∈ (0, 1], so as to analyze it more easily, and to investigate the limit  → 0. In contrast, we often drop the dependence w.r.t. the data , which is regarded as fixed.
We also take advantage of the fact that the manifold As a result, by identifying co-vectors and vectors by their components in R 3 , the functional brackets ⟨•, •⟩ below boil down to the ordinary dot product on R 3 .Similarly, the tensor product  ⊗  boils down to  ⊤ .

Asymmetric quadratic eikonal PDE
The Reeds-Shepp forward model, is defined through a sub-Finslerian quasi-metric6, relaxed by a small parameter  > 0, recall F was given by Eq. ( 8) and its data-driven version F  was given by (17).Throughout this section, and in our vessel tracking experiments, we are concerned with the case where sideward motions and backward motions become very expensive: we set spatial anisotropy parameter  =  with 0 <  ≪ 1 in the Finsler norm F  given by (17).
The generic form of the data-driven Finsler metric function considered in this paper (17) reads: for any point p ∈ Ω, within a given bounded connected domain Ω ⊂ M 2 , and any tangent vector p ∈  p (M 2 ) ≡ R 3 , and where the two small parameters ,  relate via In the following analysis, we only use the property that the tensor field  0 is pointwise positive definite, that the differential forms  1 and  2 are pointwise linearly independent, and that  0 : Ω →  ++ 3 , and  1 ,  2 : Ω → R 3 (following the conventions of Remark 19) have Lipschitz regularity.Here  ++ 3 denotes the space of 3 × 3 symmetric positive definite real-valued matrices.

Remark 20
In the normal left-invariant setting F =1 = F the asymmetric metric expressed in the fixed frame ( , , ), of the tangent space at any coordinates (, , ), has a block diagonal structure.In contrast the data-driven metrics G  , F  , in general, does not have a block-matrix structure, as the independent elements {   } 3 =1 may point anywhere due to their data-driven nature, as can be seen in Figure 6, keeping in mind the duality (23).Therefore, we must introduce a new modification of the anisotropic fast-marching algorithm.
The purpose of the second term in ( 42) is to increase the cost of sideways motions, whereas the final term prevents motions in reverse gear; both are excluded in the genuine Reeds-Shepp forward car model obtained in the limit (akin to [29, Thm.2]) as  → 0, which is non-holonomic.
The distance map  : Ω → R from a given point p 0 ∈ Ω and w.r.t. the Finsler function F  , is the unique viscosity solution to the following anisotropic eikonal system where the dual Finsler function equals by definition ).The structure of the metric (44), referred to as asymmetric quadratic, allows to compute a closed form expression of the dual metric F *  , and thereby the eikonal PDE ( 44), as we will see below.Then F is a quasi-norm (i.e. an asymmetric norm), whose dual quasi-norm reads for all p ∈ R 3 with  := ( +   ⊤ ) −1 and  :=  −1 / Proof See [29, Lemma 4].
For concreteness, we apply Lemma 2 to our Finsler function F  of interest (44), defined pointwise by the parameters This then results in the following dual Finsler functions:

Lemma 3 (Dual Finsler Functions)
With our choice (42) of Finsler function F  used in (44), the dual Finsler function F *  is given for all p ∈  * p (M 2 ) ≡ R 3 by where we used the shorthand notation  −1 := ( 0 ) −1 , the cross product  :=  1  ×  2  , and the orthogonalization coefficient Proof Follows by Lemma 2 and Taylor expansion, for details, see Appendix F.
Lemma 3 shows that one can define an ideal sub-Finsler function F * 0 that arises in the limiting case  ↓ 0, and that Our goal, achieved in Sections 5.1 and 5.2 below, is to design a numerical scheme that is consistent with the sub-Finslerian eikonal PDE F * 0 (p, d (p)) = 1, and which satisfies the conditions that make the fast marching algorithm applicable.

Discretization and consistency
We discretize the eikonal PDE (44), which has an asymmetric quadratic structure (45), using adaptive finite differences on the Cartesian grid Ω ℎ := Ω ∩ ℎZ 3 of the domain Ω, where ℎ > 0 denotes the grid scale.Note that 2/ℎ must be an integer, for consistency with the periodic boundary conditions in the angular coordinate.The numerical scheme construction relies on Selling's decomposition of positive definite matrices [54] and on a corollary related to the approximation of the squared positive part of a linear form.The versions of these results presented in [35,Corollary 4.12,Corollary 4.13] are gathered in the following proposition.

Proposition 2 (Selling matrix decomposition, see [35]) For any
For any In addition The above holds with the constants  := 6,  := 4 √

3.
Remark 21 A key aspect of Proposition 2 is that the vectors ( e  ) and ( f  ) have integer coordinates, hence we avoid (offgrid) interpolations in our discretization.
We next analyze the approximation error towards the ideal model as  → 0 and ℎ → 0 suitably.To this end we denote by  ℎ  the finite differences scheme on Ω ℎ associated with the parameters   and   of our application (47).Note that the stencil radius is (46).Now combining ( 50) with ( 52), we obtain the overall consistency error The optimal order of the consistency error O (ℎ 3 ) is achieved with the scaling  = ℎ 1 3 .In our practical experiments however, we rather use the small fixed value  =  = 0.1 which by (53) indeed yields a sufficiently accurate scheme [35]!

Anisotropic Fast-Marching for Computing Distance Maps of Data-driven Left-invariant Finsler Models
In fast-marching methods (FMM), one divides the grid points into 3 categories: Far, Trial, and Accepted.In each step of the FMM, one selects the trial point p whose function value  (p) is minimal.The point p is moved into the accepted set, and  (p) is frozen.In contrast, all the trial or far points whose stencil contains p see their function value updated, and they are moved into the trial set.This procedure generalizes and abstracts the classical FMM [52], for details see [22].When all points have moved to the accepted category, the FMM terminates, and a geodesic can be easily calculated by steepest descent as described in Section 5.4.
The update of a single function value  (p) is defined as follows: isolate  (p) in the numerical scheme expression (51), and express it by the values of its neighbors so as to obey F (p) = 1.The latter equation admits by construction a unique solution, which is obtained as the largest root of a quadratic equation.
There are two crucial properties of the discretization : -Discrete Degenerate ellipticity: (p) is a non-decreasing function of the finite differences in (51).-Causality:  (p) only depends on the positive part of all finite differences in (51).
These are the two key7 assumptions of [53,Theorem 2.3], implying that the discretized PDE (53) admits a unique solution   ℎ , which is computable in a single pass over the discretization domain Ω ℎ , using anisotropic fast-marching.
Following the steps of the proof [35] associated with the standard Reeds-Shepp forward model, and with straightforward adaptations, we obtain that   ℎ converges uniformly as  → 0 and ℎ/ → 0 to the solution  of the sub-Finslerian Eikonal PDE F * 0 (p, d (p)) = 1.

Steepest Descent for the Finslerian Geodesics
In previous work [29,Thm.4],standard Riemannian steepest descent tracking on the distance map , recall (37) in Theorem 1, is generalized from the symmetric Riemannian setting to the (possibly asymmetric) Finsler model setting.
That idea also transfers to the new data-driven left-invariant model as the Finslerian back-tracking [29, App.B] still applies.Steepest descent tracking (37) from p to source point p 0 again becomes (using the embedding of M 2 ⊂ R 3 ) with 0 <  ≪ 1 and where the derivative is taken with respect to the second entry of the dual Finsler function, and where  is the viscosity solution of the eikonal PDE system (44).
In the practical implementations we use a second order Runge-Kutta method for time integration approximations and at time  = 1 we arrive at the source-point p 0 .
This geodesic backtracking in (M 2 , F  ) again boils down to piecewise concatenations of Remark 22 Taking the negative part of the final term in (42) means that we switch between two Riemannian manifolds (one with the final term active and with the final term non-active).
At locations where  3  ≈  3 , for instance at locations where A  1 ≈ A 1 this means that one switches between anisotropic geodesics and spherical geodesics (in-place rotations).In the vessel tracking applications we want such in-place rotations to appear above bifurcations in the vasculature.
A closely related situation is discussed in [29,Thm.4],but now it is applied to the new data-driven model F  (17) with dual (F  ) * = lim  ↓0 F *  .By Theorem 1 the Riemannian geodesics have parallel momentum and their transitionstowards spherical geodesics is like  continuously differentiable.The surface S where Finslerian geodesics of F  in M 2 switch from one Riemannian manifold to the other is given by Interestingly, in the mixed-model F  that we will explain later in Section 6, the condition in Remark 22 above is satisfied at bifurcations.Then in-place rotations are indeed automatically placed at the bifurcations in the tracking of blood vessels, as can be seen in Figure 16.

Experiments
We choose the data-driven left-invariant metric tensor field with forward gear F  as given in Eq. ( 17).An elaboration on the calculation of the cost function can be found in Appendix D. We will discuss the new model's ability to adapt for curvature.Additionally, we show and discuss some full vascular tree tracking results.

Curvature Adaptation
The data-driven left-invariant metric tensor field G  and its asymmetric variant F  both consists of a "standard" leftinvariant metric tensor field to which a data-driven term is added, as introduced in (17).The idea of the second datadriven term in this equation is that it pushes the main direction of the model into the direction of the underlying vessel, as illustrated in Fig. 8, where no data-dependent cost function  = 1 was used to generate the tracking result.We see that the data-driven left-invariant metric tensor field takes the image data into account and steers the tracking in the direction of the underlying vascular structure, even when the cost function does not contain information about vessel locations and curvature.
The data-driven term also leads to better tracking results of very tortuous vascular structures as we see in Fig. 9.In Fig. 9a, the tracking results relying on (solely) the left-invariant metric tensor field F fail to follow the underlying vessel correctly, contrary to the new data-driven left-invariant model F  (17) which follows the vascular structure correctly.In addition, one sees that when using the left-invariant model, the wave fronts (indicated in orange) suffer from the discretization.In the data-driven left-invariant model, these discretization artifacts are gone, and the wavefronts follow the underlying vascular structure correctly.When only considering 8 orientations, as in Fig. 9b, the data-driven left-invariant frame is still able to follow the vascular structure correctly.Although the discretization is clearly visible in the tracking results, the data-driven left-invariant metric tensor field is still able to follow the vessel correctly.It is important to note that both models use the same cost function.The differences in the tracking results are due to the data-driven left-invariant metric tensor fields' ability to better adapt for: 1. gradual curvature change of blood vessels.The same applies to other applications such as the detecting of cracks, see Figure 10,  2. orientation biases by orientation score data-alignment as depicted in Figure 6.

Complete Vascular Tree Tracking
In the previous section, we discussed the curvature adaption feature of the new (asymmetric) data-driven left-invariant metric tensor field F  .This model also can automatically place 'in-place' rotations in globally optimal geodesics which are typically placed at bifurcations.).We see that the iso-contours of the data-driven metric tensor field follow the vessel structure better, and the tracking is correct (even with 8 orientations).55), is preferable as it only adapts for curvature (like in Fig. 9) in between those complex structures, and indeed provides correct tracts everywhere, as can be seen in Fig. 11c.The geodesics of both models are computed using  = 100.
However, these valuable abilities of the model may also lead to extremely complex structures to some undesirable behavior.In full vascular tree tracking, we see that the datadriven term may steer the tracking away from the actual vessel at crossings in extreme cases, as can be seen in Fig. 11.
To overcome this problem (see item c in Fig. 11), and to still take advantage of improved data-adaptation (like in Figs.5, 8 and 9) we introduce a (new) mixed model that prevents this undesirable behavior at (extreme) complex structures, where it locally relies on the old model.Then in between (extreme) complex structures we still benefit from the directional data adaptation in the orientation score.
The mixed metric tensor field G  (and its asymmetric version F  ) is given by the data-driven left-invariant metric tensor field away from the crossing structures, and by the The results are typically not sensitive to the choice of  and  in our application as long as  > 2. Therefore we always set  = 5 and  = 1 pixel-size in our experiments.This construction of the metric tensor field ensures that the metric tensor field is not tempted to move in the wrong direction in extreme cases where vessels cross each other.The tracking result relying on the mixed metric tensor field is visualized in Fig. 11c, and does not show the earlier mentioned undesirable behavior, as shown in Fig. 11b.Therefore, this new model will be used in all full vascular tree tracking results.All results presented in this section are calculated using parameters  11 = 0.01,  22 =  33 = 1.For the curve optimisation this is the same as setting  =  = 0.1 in (15) used in (17).Even for such extreme anisotropy settings, our numerical algorithm is appropriate as motivated in Section 5. We always observed that tracking is stable with respect to small variations in these parameters, so there was no point in fine-tuning them.

Asymmetric Double Step
The tracking results were computed in two steps; first tips are connected to the nearest bifurcation/seed (cost function visualized in Fig. 12), after which those points are connected to the nearest seed (cost function visualized in Fig. 13).The used cost functions (from tip to nearest bifurcation and from bifurcation to seed) support movement along the thin and thick vessels respectively.The tracking results that correspond to these cost functions are depicted in Fig. 14.The calculated geodesics are all correct, except for 2 difficulties when: 1. crossings and bifurcations are very close to each other, 2. vascular structures are kissing.
Next, we compare the results of the new mixed model F  and the left-invariant model F in Fig. 15.There are some visible differences between both tracking methods, marked in pink and blue.First, we see that the tracking results relying on the mixed method ensure that the centerline is better followed, and multiple geodesics are at approximately the same place (in blue).Second, we see the ability of the mixed method to adapt to the direction of the vascular structure (in pink).Instead of moving towards a bifurcation point away from the seed in the left-invariant metric tensor field, the curvature adaptation ensures that the tracking results immediately move towards the seed it is connected to.(55).The first step involves connecting the tips (marked in red) to the nearest bifurcation (marked in purple) using the cost function depicted in Fig. 12.Second, these bifurcations are now tracked to the seeds (marked in green) using the cost function depicted in Fig. 13.The first step involves connecting the tips (marked in red) to the nearest bifurcation (marked in purple).Second, these bifurcations are now tracked back to the seeds (marked in green).The cost functions used in the first (white) and second (blue) step are given in [29], with  = 800 and  = 4, and Fig. 13 respectively.The main differences between both models indicate that the mixed model F  follows the vascular structure better and is better able to follow the centerline of a given vascular tree (marked by pink and blue circles respectively).

Asymmetric Single Step with Prior Classification of Seeds and Tips
Common practical setups in vascular tracking of retinal images include the prior knowledge of the locations of tips and seeds of vessel structures.We implemented our datadriven model using a prior classification of the connected tips and seeds.More specifically, in every run of the fast marching algorithm, one of the seeds is considered together with its corresponding tips, and the connecting vessel structures are tracked.Figure 16 shows our result in this setup and demonstrates that our approach can determine the geodesics that accurately follow the vascular structure in the retinal image.

Accuracy of the Model
We now present some quantitative evaluations to measure the accuracy of our data-driven metrics for geodesic tracking.We measure the mistake ratio  for the images in the STAR dataset.For these images, the ground truth of the vessels is known, which allows us to calculate the percentage of the vessel that is not on the correct vascular structure, where  = # pixels not on correct vessel # pixels of all geodesics .
We have calculated this accuracy for images of the STAR dataset where we connect the tips to their nearest bifurcation, since one should use the new model away from crossing structures.The results are presented in Fig. 17.We see that for most tracks, the performance improves when switching to the new data-driven model, and in the cases where there is no improvement, the results do not get significantly worse.On average, we find an improvement of 10.7% of the calculated tracks for the considered images.

Conclusion and Future Work
In this article, we introduced the concept of a data-driven leftinvariant metric tensor field G  and its asymmetric variant F  .The metric tensor field is defined by the underlying image, where movement along line structures is encouraged by its design in (17).In addition, a data-driven version ∇  of the plus Cartan connection, relying on G  , was introduced.
We used these geometrical tools to formulate a challenging data-driven version of [44, Thm.1], which was stated and proved in Theorem 1.In this theorem, 'straight' and 'short' curves are described with respect to the data-driven Cartan connection.In particular, it describes the entire Hamiltonian flow of the new Riemannian manifold model (M 2 , G  ) in terms of the new data-driven Cartan connection ∇  , and   (55).Prior grouping of tips (in red) and seeds (in green) results in perfect tracking of the vessel tree, using only one efficient anisotropic fast-marching run via the numerical method in Section 5.Both results are calculated with the cost function visualized in Fig. 13.At the purple points, we have bifurcations and our tracking is solely based on the mixed model produced (spatially projected) geodesics  with (automatic) in-place rotations at such bifurcations.explains the backtracking procedure for backtracking datadriven left-invariant geodesics in (M 2 , G  ).As subsequently explained this can be extended to the asymmetric Finsler model (M 2 , F  ) that often yields the same geodesics, but also automatically places in-place rotations.The latter is beneficial at bifurcations in complex vasculature when using crossing-preserving vesselness costs for the costfunction .
The diagonalization of the new data-driven left-invariant models G  and F  provides locally adaptive frames that are beneficial over previous approaches to locally adaptive frames in M 2 [38,56,68] in the sense that: 1. they coincide with the usual left-invariant frame if the data is locally constant, 2. they are more stable as they are constructed by coercive metric tensor fields, recall (21).
To calculate the minimizing geodesics efficiently, an adaptation to the efficient anisotropic fast marching algorithm was required and presented in Section 5.The metric tensor component matrix was no longer of block form in the fixed coordinate system, and the necessary changes to overcome this have been discussed and analyzed in Section 5. We also provide an asymptotic error analysis of our numerical scheme.
To show the performance of the data-driven metric tensor field and the mixed metric tensor field, we have tested them on 2D images of the retina.All experiments confirm that the new model is better able to adapt for curvature.In addition to that, for the tracking of a single vessel, a low number of orientations is sufficient to find the correct minimizing geodesic, as can be seen in Figure 9. Full vascular tree tracking needs to be handled with care at difficult crossings structures, which is done in the mixed model F  introduced in (55).
In general, the tracking results perform very well in the discussed two-step approach (see e.g. Figure 15), where tips are first connected to the nearest bifurcation, after which the geodesics connecting these bifurcations to the corresponding seeds are calculated.After prior classification of seeds and tips belonging to the same vascular tree, the tracking results follow the vessels perfectly, recall Figure 16.
Despite some very appealing theoretical and practical advantages of our model, we still require considerable computation and runtime (tripling the overall processing time) to make the data-driven metric-tensor field and distance maps.Therefore, the exact usage of the proposed data-driven metric depends on the specific context of the tracking requirements.
For future work, it would be interesting to look into the possibilities to train the cost function  using PDE-G-CNNs [69], which is now geometrically computed as explained in Appendix D. In the past, this method had promising results for vessel segmentation.Besides using PDE-G-CNNs to construct the cost function, it would be worth looking into the possibilities to use neural networks to calculate the distance function as was done in [70].Secondly, we address the characterisation of shortest curves.
where  denotes the Hamiltonian flow, and where ∇  and ∇  denote the gradient with respect to  and  respectively.
This follows from the computation of the Hamiltonian via the Fenchel transform.The vertical part is given by λ In the above derivation, we have used that Additionally, it is important to note that Hence, we also have found that λ = γ and λ = γ , which we aimed to show.Note that reformulation in a coordinate-free matter yields ∀  ∈ {1,...,} γ = λ ⇔ |  ( ) = G −1   ( ) (), for all  ∈ [0,  (g)].
Remark 25 In Theorem 1 we give a backtracking formulation (where geodesics go 'down-hill' to the origin via steepest descent) where we rescaled time back to the interval  ∈ [0, 1] (as this is more practical).Similar to [29,Thm.4,Eq.28]this is done as follows: in the ODE backtracking for geodesics (37) we included an extra negative scaling factor − (g) in comparison to all the canonical ODEs above.
Thirdly, we address the symmetry (36) of the data-driven Cartan connection.
By construction of ( 17) and ( 28), we have the correct symmetry (36).Indeed from (17), it follows that A L q   gp = ( g ) * A   p and   L q  gp = ( g ) *    p and via (28) we get ( g ) * ∇ L g  * = (∇  ) * , where we use the fact that G  is diagonal w.r.t.basis {A   } and where we respectively applied the pushforward of a vector field, the pullback of a co-vector field, and the pullback of a dual connection.
Finally, we address the integration of geodesics and their symmetry.

Appendix C The Used Metric Tensor Field is indeed a Data-Driven Left Invariant Metric Tensor Field
We first rely on a convenient standard formula of the Hessian of smooth function on a manifold relative to a connection on that manifold in Lemma 4. Then we provide an alternative formulation of such a Hessian in Lemma 5 (via the notion of parallel transport).Finally, we prove that G  , that heavily relies (17) on a Hessian  of a sufficiently smooth orientation score  :  (2) → R, is indeed a data-driven left-invariant vector field in Lemma 6. i.e.X is the unique autoparallel curve through p with tangent vector  p .For all ,  ∈ [−, ] let  X , :  X ()  →  X ( )  be the parallel transport operator along the curve X, which is uniquely defined by the following properties 1.  X  , =  ∀ ∈ [−, ], 2.  X  2 , 3 •  X  1 , 2 =  X  1 , 3 , ∀ 1 ,  2 ,  3 ∈ [−, ], 3. smooth with respect to X,  and .
Then  ↦ →  X 0,  p ∈   ( )  gives a smooth vector field along the curve X that is unique parallel transport of  p along that curve, i.e. with the property ∇ X ( )  X 0,•  p = 0.In order to prove that G  is a data-driven left invariant metric tensor field, we need the following identities: ∇ (g) *   g *  =  g * ∇   , where (63) is the definition of the pushforward, and ( 64) is the equivariance of the Cartan connection ∇ = ∇ =1 .In addition, it is important that so ∥ ∥ = ∥ Ỹ ∥, where  ∈  p (M 2 ) and Ỹ ∈  gp (M 2 ).
We check the data-driven left invariance for the separate terms of the metric tensor field, starting with G p ( p, p): ∥| p ( q, •) ∥ 2 * .

( a )
Image  .(b) Orientation Score  belonging to image  .(c) Visualization of the real part of the cake-wavelet, with orientations  = 0 (left) and  = 3 /16 (right), used to lift image  to the orientation score.

Fig. 1
Fig. 1 (a) Visualization of a grayscale image  : R 2 → R and (b) its corresponding orientation score  : M 2 → R in the space of positions and orientations M 2 given by (3), using a standard cakewavelet as depicted in 1c.We use a volume-rendering where the orange spirals indicate data-points p = ( , ,  ) ∈ M 2 with high amplitudes | (p) |.

Fig. 5
Fig. 5 Top row: All points in the orange surface have the same distance to the seed.The isocontours are projected back onto the image, as depicted in the ground plane.Bottom row: Several isocontours are projected onto the image and a projection of the curve is visualized.Left: The Riemannian geodesic with parallel velocity to the Cartan plus connection ∇ [+] (red) takes a wrong shortcut.Right: The Riemannian geodesic (green) with parallel velocity to the Cartan plus connection ∇  in the Riemannian manifold (M 2 , G  ) -or more precisely the Finslerian manifold (M 2 , F  )given by (17)-does not, and moreover adapts for curvature in M 2 , cf.Fig.7).Explicit formulas for ∇  and ∇[+] will follow later in Table2.

Fig. 6
Fig.6 Visualizations of the left-invariant frame and the data-driven leftinvariant frame in M 2 .Locally along one of the spirals in the orientation score depicted in Fig.1.In Fig.6a, the main direction of A 1 is not properly aligned with the underlying 3D structure, whereas in Fig.6bA  1 is.

Fig. 8
Fig. 8 Influence of data-driven metric tensor fields: (top) Tracking with the vanilla left-invariant metric tensor field from (14).(bottom)Tracking with the proposed data-driven left-invariant metric tensor field from(17).To isolate its effect in the tracking process and record the effect of only directional adaptation of the underlying metric, we have set the cost function  = 1.We observe that the data-driven nature of our model allows for a better fidelity to the underlying vascular structure.The parameters are given by  11 = 0.01,  22 = 1,  33 = 1,  = 100.

( a )Fig. 9
Fig. 9 Comparison tracking results of left-invariant and data-driven left-invariant metric tensor field: Tracking results for left-invariant metric tensor field F (left) and data-driven left-invariant metric tensor field F  (right).The parameters for the (data-driven) left-invariant metric tensor field are given by  11 = 0.01,  22 = 0.16,  33 = 1,  = 100.The cost function is given by  = 1/(1 + 200|  | 2).We see that the iso-contours of the data-driven metric tensor field follow the vessel structure better, and the tracking is correct (even with 8 orientations).

( a )
Tracking results with the left-invariant metric tensor field.(b) Tracking results with the data-driven left-invariant metric tensor field with  = 100.

Fig. 10
Fig. 10 Application in crack detection: Tracking results for the leftinvariant and data-driven left-invariant metric tensor field on an image of cracks in a building.The presented results are calculated using 32 orientations, and parameter settings  11 = 0.01,  22 = 1,  33 = 1.In regions with high curvature, the data-driven model adapts more gradually for curvature to get more data evidence than the left-invariant model which tends to prefer in-place rotations.

( a )
Fig. 11 Motivation mixed model: Too much curvature adaptation at crossings is dangerous in extreme cases.The mixed model, introduced in Eq. (55), is preferable as it only adapts for curvature (like in Fig.9) in between those complex structures, and indeed provides correct tracts everywhere, as can be seen in Fig.11c.The geodesics of both models are computed using  = 100.

Fig. 12
Fig. 12 Projected cost functions for tracking in two steps -from bifurcation to tip: Cost used to connect tips to the nearest bifurcation.Black and white mean low and high costs respectively.This cost function supports movement along the thin vessels very well.The multiscale vesselness is computed as explained in App.D, and the considered spatial scales is   = 1.

Fig. 13
Fig. 13 Projected cost functions for tracking in two steps -from seed to bifurcation: Cost used to connect bifurcations to optic nerve.Black and white mean low and high costs respectively.This cost function supports movement along the thick vessels very well.The multiscale vesselness is computed as explained in App.D, and the considered spatial scales are   ∈ {1, 2}.

(a) Image 1 :
Tracking with the mixed model with  = 20.(b) Image 2: Tracking with the mixed model with  = 15.

Fig. 14
Fig. 14 Two step tracking of Vascular Tree structures: Tracking with mixed model (M 2 , F  ) proposed in(55).The first step involves connecting the tips (marked in red) to the nearest bifurcation (marked in purple) using the cost function depicted in Fig.12.Second, these bifurcations are now tracked to the seeds (marked in green) using the cost function depicted in Fig.13.

( a )
Tracking with the left-invariant model F. (b) Tracking with the mixed model F  proposed in (55) with  = 65.

Fig. 15
Fig. 15 Two-step tracking of vascular tree structures:The first step involves connecting the tips (marked in red) to the nearest bifurcation (marked in purple).Second, these bifurcations are now tracked back to the seeds (marked in green).The cost functions used in the first (white) and second (blue) step are given in[29], with  = 800 and  = 4, and Fig.13respectively.The main differences between both models indicate that the mixed model F  follows the vascular structure better and is better able to follow the centerline of a given vascular tree (marked by pink and blue circles respectively).

( a )
Image 1: Tracking with the mixed model with  = 50.(b) Image 2: Tracking with the mixed model with  = 15.

Fig. 16
Fig. 16 Tracking of Vascular Tree per Seed: Tracking with the mixed model (M 2 , F  ) proposed in(55).Prior grouping of tips (in red) and seeds (in green) results in perfect tracking of the vessel tree, using only one efficient anisotropic fast-marching run via the numerical method in Section 5.Both results are calculated with the cost function visualized in Fig.13.At the purple points, we have bifurcations and our tracking is solely based on the mixed model produced (spatially projected) geodesics  with (automatic) in-place rotations at such bifurcations.

Table 2
Table of Geometric Tools and Notations 4.1 Combine Optimally Straight and Short: A new Data-Driven Version ∇  of the Cartan Connection 1. symmetric Riemannian geodesics in M 2 with metric tensor field G p ( p, p) recall (17), and 2. symmetric Riemannian geodesics in M 2 with metric tensor field G  p ( p, p)+( −2 −1)| 1  ( p)| 2 that are in-place rotations, at locations where (1,2,8,9,13,15,16,24,38,48)rplot of the accuracy of the mixed model vs. the left-invariant model applied on images in the STAR dataset(1,2,8,9,13,15,16,24,38,48).The accuracy is calculated on the calculated tracks between the tips and the nearest bifurcation for all vessels in one single run.The red area marks where the former left-invariant model performs better than the new mixed model, incidences indicated by an 'x'.The green area marks where the new mixed model performs better than the former left-invariant model, incidences are indicated by an 'o'.Most measurements show the improved performance of the mixed models compared to the left-invariant model (LIF).