Locally Adaptive Frames in the Roto-Translation Group and their Applications in Medical Imaging

Locally adaptive differential frames (gauge frames) are a well-known effective tool in image analysis, used in differential invariants and PDE-flows. However, at complex structures such as crossings or junctions, these frames are not well-defined. Therefore, we generalize the notion of gauge frames on images to gauge frames on data representations $U:\mathbb{R}^{d} \rtimes S^{d-1} \to \mathbb{R}$ defined on the extended space of positions and orientations, which we relate to data on the roto-translation group $SE(d)$, $d=2,3$. This allows to define multiple frames per position, one per orientation. We compute these frames via exponential curve fits in the extended data representations in $SE(d)$. These curve fits minimize first or second order variational problems which are solved by spectral decomposition of, respectively, a structure tensor or Hessian of data on $SE(d)$. We include these gauge frames in differential invariants and crossing preserving PDE-flows acting on extended data representation $U$ and we show their advantage compared to the standard left-invariant frame on $SE(d)$. Applications include crossing-preserving filtering and improved segmentations of the vascular tree in retinal images, and new 3D extensions of coherence-enhancing diffusion via invertible orientation scores.


Introduction
Many existing image analysis techniques rely on differential frames that are locally adapted to image data. This includes methods based on differential invariants [54,39,32,47], partial differential equations [54,63,38], and non-linear and morphological scale spaces [13,12,64], used in various image processing tasks such as tracking and line detection [6], corner detection and edge focussing [39,8], segmentation [58], active contours [15,16], feature based clustering etc. These local coordinate frames provide differential frames directly adapted to the local image structure via a Hessian or a structure tensor. Their primary benefit is that they allow to include anisotropy and curvature in a Euclidean invariant way. They induce complete sets of differential invariants, [32], which are typically applied in greyscale imaging for segmentation [58,50], low-level object recognition, etc. See Fig.1, where we have depicted local adaptive frames (also known as 'gauge frames' according to the terminology used in [32,Section 3.3.3] and [10,39]) based on eigenvector decomposition of the Gaussian Hessian of the MR-image in the background. It is some- Hessian H s f (x) = (G s * Hf )(x) and its eigensystem provides us a basis { ∂ a | x , ∂ b | x } for T x (R 2 ) at each location x ∈ R 2 "gauged" with image f . Right: Such gauge frames can be used for adaptive anistropic diffusion and geometric reasoning. The problem however is that at isotropic, more complex structures such as blobstructures/crossings the gauge frames are ill-defined causing weird fluctuations.
times problematic that such locally adapted differential frames are directly placed in the image domain R d , as at the vicinity of complex anisotropic structures, e.g. crossings, textures, bifurcations, one typically requires multiple local spatial coordinate frames, cf.
embedded as the partition of all left-cosets in the rototranslation group SE(d) on R d . Such techniques rely on various kinds of lifting, such as coherent state transforms (orientation scores), continuous wavelet transforms [22,26,55,6], orientation lifts [65,11], or orientation channel representations [31]. The two key advantages of this image domain extension is that processing in the space R d S d−1 of positions and orientations allows 1. to disentangle oriented structures involved in crossings, and to include curvature, cf. Fig.2, Fig.4, Fig.5. 2. to include stochastic PDE models for the alignment of all local orientations in the image, see e.g. [49,65,4,25,19,26,11,5,53,55,62,17,56,46].
In this article we will not discuss in detail on how such a new image representation or lift is to be constructed from grey-scale image f : R d → R, and we assume it to be a sufficiently smooth given We aim for adaptive anisotropic diffusion of images (here d = 2) which takes into account curvature. At areas with low orientation confidence (in blue) isotropic diffusion is required, whereas at areas with high orientation confidence (in red) anisotropic diffusion with curvature adaptation is required. Application of locally adaptive frames directly in the image domain suffers from interference (3rd column), whereas application of locally adaptive frames in the domain R d S d−1 of distributions U (e.g. invertible orientation scores, U = W ψ f , cf. Fig. 3) allows for adaptation along all separate elongate structures (4th column). input. Here |U (x, n)| is to be considered a probability density of finding a local orientation (i.e. an elongated structure) at position x ∈ R d with orientation n ∈ S d−1 . In case U is obtained via a coherent state transform of a grey-scale image (also known as invertible orientation score), cf. [2,26,6,33], which we will take as default, one typically takes the real part or the modulus [33, ch.5.3] for locally optimal frames, see Fig. 3.
In case one has to deal with more complex diffusion weighted MRI-data, the function U can be obtained after a modeling procedure (e.g. via the common approaches [59,60,1,57]).

Goals
In this article, our quest is to find locally optimal differential frames associated to the smooth function U : R d S d−1 → R, relying on similar Hessian-and/or structure-tensor type of techniques for gauge frames on images, recall Fig.1. The overall pipeline of including locally adaptive frames in suitable operators Φ (e.g. diffusion operators Φ = Φ t with diffusion time t > 0) in orientation lifts U : R d S d−1 → R of images f : R d → R, and the resulting effective image processing Υ , is depicted, for d = 2, in Fig. 4. For d > 2 the same pipeline applies. With this overall pipeline in mind, the main goals of this article are: 1. To overcome the problem of single locally adaptive image frames per position on R d , and to generalize well-known approaches via structure tensor or Hessian to R d S d−1 , d ∈ {2, 3}. 2. To set up adaptive local frames for each pair (x 0 , n 0 ) in the joint space of positions and orientations, based on specific curves t → (x(t), n(t)) that locally fit the data (x, n) → U (x, n) best. 3. To formalize best local exponential curve fits and to include them in data-driven diffusions via efficient algorithms. 4. To extend adaptation for curvature and deviation from horizontality in 2D, cf.
[33], to 3D, and to investigate the usage of structure tensor or Hessian. 5. To employ our locally optimal frames in medical imaging in improved crossing-preserving diffusions U → Φ t (U ) and improved differential invariants.
In order to achieve these goals, one must rely on the semi-direct product group structure on the domain R d S d−1 , Eq. (1). Therefore, in the next subsection we provide some preliminaries on this group structure. Subsequently, in Subsection 1.3, we illustrate the notion of best exponential curve fits (see the blue line and blue formula in Fig. 4), where for the sake of illustration we consider the basic case d = 2 first.

Preliminaries on Domain R d S d−1
Now let us assume input U is given and let us first concentrate on its domain. This domain equals the joint space R d S d−1 of positions and orientations of dimension 2d − 1, and is embedded as a Lie group quotient (1) in the group of rotations and translations SE(d) = R d SO(d) of dimension n d = d(d + 1)/2. The group SE(d) is endowed with semi-direct group product for all g = (x, R) and g = (x , R ) ∈ SE(d). For effective alignment-modeling it is crucial not to consider the space of positions and orientations as a flat Cartesian product without group structure, see e.g. [24, fig.7] and [28]. Instead the non-commutative relation in (2) between rotations and translation is to be included in alignment modeling and visual perception. Therefore we model the joint space of positions and orientations as the Lie group quotient (1), where we identify with some a priori set reference axis a ∈ S d−1 , with S d−1 = {a ∈ R d | a = 1} to itself. Throughout this article we set this a priori reference axis as follows: Within the Lie group quotient structure (1), two rigid body motions g = (x, R), g = (x , R ) ∈ SE(d) are by definition equivalent if g ∼ g ⇔ h := (g ) −1 g ∈ {0} × SO(d − 1) ⇔ x − x = 0 and ∃ R h ∈SO(d−1) : (R ) −1 R = R h .
As a result it is readily verified that Thereby, a single element in R d S d−1 can be considered as the equivalence class of all rigid body motions that map reference position and orientation (0, a) onto (x, n). Similar to the common identification of S 2 ≡ SO(3)/SO(2), we denote elements of the Lie group quotient R 3 S 2 by (x, n). In the basic 2D case, d = 2, we have R 2 S 1 ≡ SE(2), and we do not have to worry about the quotient structure. In the 3D case d = 3 things become much more subtle: relevant subroutines will be expressed in the extended space SE (3), and this leaves the obligation to check if the overall routine in R 3 S 2 is legal [28], i.e. independent on choice of representant in each left coset. Fig. 4: The overall pipeline of image processing f → Υ f via left-invariant operators Φ, where we compute per element g = (x, y, θ) a best exponential curve fit γ c * g (t) (in blue, with spatial projection in red) with tangenṫ γ c * g (0) = c * (g) = (c 1 , c 2 , c 3 ) T , with c * (g) the optimal tangent at g. Based on this fit we construct for each g a local frame { B 1 | g , B 2 | g , B 3 | g }. This provides a field of local frames g → { B 1 | g , B 2 | g , B 3 | g } used in our operators Φ on the lift (here Φ is a non-linear adaptive diffusion operator with fixed diffusion time).

Preliminaries on Exponential Curve Fits
In case d = 2 we have the left-invariant vector fields As explained in previous works [26,22] they form a first natural frame of reference in the domain of U , that we aim to adapt later. The corresponding dual frame (2)) is given by ω i , A j = δ i j and thus ω 1 = cos θdx + sin θdy, ω 2 = − sin θdx + cos θdy, ω 3 = dθ.
Remark 1 Throughout this article we take the algebraic viewpoint on tangent vectorsγ(t) on SE(d), i.e. we consider them as differential operators on locally defined functions (as is commonly done in differential geometry, see e.g. [3]): and φ smooth and locally defined on an open set around , and with n d := dim(SE(d)) = d(d + 1)/2. Note that n 2 = 3 and n 3 = 6.
Let us first consider the notion of exponential curves in SE(2), which we will employ in our frame-adaptation.
Let (c 1 , c 2 , c 3 ) T ∈ R 3 be a given column vector. Then the unique exponential curve passing trough g = (x, y, θ) with initial velocity c(g) with A i = A i | e denoting the tangent vector at unity element e = (0, 0, 0) in respectively x, y, θ-direction: In fact such exponential curves satisfẏ and thereby have constant velocity in the moving frame of reference, i.e.γ i = c i in Eq. (6). Via Eq. (8) and Eq. (4), and the method of characteristics, the exponential curves can be expressed as: for the case c 1 = 0, and all t ≥ 0, and for the case c 1 = 0, where g 0 = (x 0 , y 0 , θ 0 ) ∈ SE(2). In both cases we must impose in order to ensure that t equals the arc length parameter in the Riemannian manifold (SE(2), G µ ) with with µ > 0 a stiffness parameter balancing between costs of spatial motion relative to cost of angular motion, and where γ denotes any smooth curve in SE(2).
In this article exponential curves (being 1-parameter subgroups) have our primary interest (rather than the geodesics) as left-invariant PDE-evolutions locally take place along these curves, even in the non-linear setting. Therefore, in this article, we consider fitting exponential  Fig. 3). Isosurfaces (in red) with relatively large iso-values in the orientation score concentrate around horizontal curves (i.e. lifted curves with θ(s) = arg{ẋ(s) + iẏ(s)}) lifted from the plane. Left: In case of crossing straight lines, the tangent c = c(x, y, θ) of local exponential curve fit at a centerline-point (x, y, θ) (where d H = 0) at each line, is aligned with A 1 | (x,y,θ) . Right: At centerlines of curved structures, the tangent c of the exponential curve fit is tilted in the horizontal tangent plane span{A 1 , A 3 } and includes curvature adaptation.
curves rather than geodesics. The precise details on how to compute such best exponential curve fits in SE(d) will follow in Section 3 (for d = 2) and in Section 4 (for d = 3).
In Fig. 5 (12) Their practical advantage in 2D imaging has been shown in previous work [33, ch:6.5], [26], and is quickly shown in Fig. 6. In this article we generalize to d ≥ 2 and we propose several optimization problems and formal variational methods solving them to find such best exponential curve fits. We distinguish between two methods: 1. First order variation of the data along the exponential curve should be minimal, solved by the eigenvector of a generalized structure tensor of the extended data representations on SE(d). 2. Second order variation of the data along the exponential curve should be minimal, solved by the eigenvector of a generalized regularized Hessian of the extended data representations on SE(d).  Fig. 3, is that orientation bias towards sampled orientations is reduced if only 4 orientations are used. Bottom: The advantage of including curvature κ in the crossingpreserving flows. From left to right: input image f , processing without (middle) and with (right) inclusion of (d H , κ) via locally adaptive frame (14) based on 2nd order local best exponential curve fit γ c g (·) to data |U |.

From Best Exponential Curve Fit to Locally Adaptive Frame
For now let us assume that for each g ∈ SE(d) we have computed the tangent vector c(g) of the best exponential curve fit to extended data representation U . Next we illustrate how we use c(g) to construct a normalized full differential frame with n d = dim(SE(d)) = d(d + 1)/2. We stress that despite the inclusion of local data adaptation, the formal left-invariance of our algorithms is maintained.
To get a first grasp on how to construct a locally adaptive frame from a single best exponential curve fit, we start with considering the basic case d = 2, where inclusion of second order best exponential curve fits have been studied in previous works [34,26], and where Fig. 7 illustrates the geometrical idea. Note that these techniques are recently improved/refined to the multiple scale setting in [55], which we will also consider later in the application section. However, here we will present the SE(2)-case in such a way that generalizations to d ≥ 2 are in reach. These generalizations are studied and employed later in this article (for a quick preview on the case d = 3 see e.g. Fig. 9, 21 and 24) Consider d = 2 and Fig. 7, where both the frame {A 1 , A 2 , A 3 } and the locally adaptive frame {B 1 , B 2 , B 3 } in the tangent bundle T (SE(2)) are depicted. The explicit relation between the normalized Gauge frame (2)) and the left-invariant vector field frame {A 1 , A 2 , A 3 } is given by with given by with deviation from horizontality d H given by Eq. (12) and with spherical angle ν = arg(µ + iκ) ∈ [−π/2, π/2]. The multiplication M −1 µ A ensures that each of the vector fields in the locally adaptive frame B = (B 1 , B 2 , B 3 ) T is normalized w.r.t. the G µ -metric, recall (10).
Remark 3 When imposing isotropy (w.r.t. metric G µ ) in the plane orthogonal B 1 (e.g. done in the crossingpreserving diffusions in [26,33]), the precise choice of R c ∈ SO(3) mapping (1, 0, 0) T onto (µc 1 , µc 2 , c 3 ) T is not relevant. However, in practice one may not insist Fig. 7: Locally adaptive frame { B 1 | g , B 2 | g , B 3 | g } (in blue) in T g (SE(2)) (with g placed at the origin) is obtained from frame { A 1 | g , A 2 | g , A 3 | g } (in red) and c(g), via normalization and two subsequent rotations R c = R 1 R 2 , see Eq. (14), revealing deviation from horizontality d H in R 1 , spherical angle ν, and curvature κ in Eq. (15). Vector field A 1 takes a spatial derivative in direction n, whereas B 1 takes a derivative along the tangent c of the local best exponential curve. on such isotropy, and then there exists only one natural and efficient choice: R c = R 1 R 2 which is the 3Drotation which maps (1, 0, 0) T onto (c 1 , c 2 , c 3 ) T in such a way that (0, 1, 0) T ↔ A 2 is mapped onto B 2 which is again within the spatial tangent plane. Note that The generalization to the d-dimensional case of the construction of locally adaptive frame and the tangent vector c (which we represent by a column vector) of a given best exponential curve fitγ c (·) to dataŨ (x, R) := U (x, Ra) is surprisingly simple as briefly explained in Appendix A. For explicit formulae of the left-invariant vector fields in the d-dimensional case see [30], and for d = 3 see [29]. For a quick illustration of locally adaptive frames in d = 3 see Figure 24 in Appendix A.

Structure of the Article
In the introduction (and Appendix A) we have explained how to derive a locally adaptive frame {B 1 , . . . B n d } in SE(d) from a single best exponential curve fit γ c g (·) to dataŨ : SE(d) → R. Therefore, the body of this article is devoted to theory and algorithms to construct these best exponential curve fits at each g ∈ SE(d).
As the construction of local best exponential curve fits to data on SE(d) is technical (in particular for d = 3), we first reformulate the well-known approach to Gauge frames in smooth images f : R d → R via the structure tensor from a group theoretical viewpoint in Section 2. This gives a roadmap towards SE(2)extensions explained in Section 3, where we deal with new exponential curve fits of the 1st order in Subsection 3.1 computed via a structure tensor, and exponential curves of 2nd order in Section 3.2. The 2nd order exponential curves are computed via a symmetric sum or a symmetric product of a non-symmetric Hessian.
The new technical extension to SE(3) is done in Section 4. In Section 4.2 we provide a detailed differential geometrical approach to first order best exponential curves in SE(3), including formal parallel transport, well-posed construction of optimal curves on the quotient R 3 S 2 = SE(3)/({0} × SO(2)), and global optimality w.r.t. the corresponding energy functional. We also prove that in order to obtain exponential curves which are both torsionfree and which allow for well-posed projection on R 3 S 2 , one must resign to a two-step optimization algorithm presented in Subsection 4.2.3. This provides excellent local exponential curve fits, even if the data includes elongated structures with torsion. In Section 4.3 we present second order best exponential curves in SE(3) minimizing an energy functional and computed via the symmetric sum of the Hessian. Again we must rely on a two-step optimization procedure for torsion-free exponential curve fits, similar to the algorithm presented in Subsection 4.2.3. We also indicate serious problems of extending the symmetric product approach [34,26] from SE(2) to SE (3).
Subsequently, in Section 5 we consider new experiments regarding medical imaging applications and feasibility studies that show the advantage of including locally adaptive frames over non-adaptive left-invariant frames in invertible, (multiple-scale) orientation scores. In this section we distinguish between the SE(2) case, with clear advantages in crossing-preserving multi-scale vesselness filters and subsequent vessel detection in retinal imaging as we will show in Subsection 6.1, and the SE(3)-case where we include first feasibility studies that show proof of concept of the advantage of crossingpreserving coherence enhancing diffusion steered by locally adaptive frames via invertible 3D orientation scores on artificial dataset defined on R 3 .
Finally, we provide conclusions and envisioned 3D medical imaging applications in future work.

Optimal Exponential Curves in R d
In this section we reformulate the classical construction of a locally adaptive frame to image f at location x ∈ R d , in a group-theoretical way. This reformulation seems technical at first sight, but helps in understanding the formulation of best projected exponential curves in the higher dimensional Lie group SE(d). We will take the structure tensor approach [9,45], which will be shown to yield first order best exponential curve fits.
The Gaussian gradient with Gaussian kernel G s (x) = (4πs) −d/2 e − x 2 4s , is used in the definition of structure matrix: with s = 1 2 σ 2 s , and ρ = 1 2 σ 2 ρ the scale of regularization typically yielding a non-degenerate and positive definite matrix. For all x ∈ R d we have optimal tangent vector The solution to optimization problem (17) is obtained by the Euler-Lagrange equation and the minimizer is found as the eigenvector c * (x) with the smallest eigenvalue. Now let us put Eq. (17) in group-theoretical form. This is helpful in our subsequent generalizations to SE(d).
On R d exponential curves are straight lines: and on T (R d ) we impose standard flat metric tensor (17) can be rewritten as where the curve is a parallel transported curve associated to t → γ c x (t).
Definition 1 Let c * (x) ∈ T x (R d ) be the minizer in (19). We say γ x (t) = x + e tc * (x) is the best exponential curve fit to image-data f : R d → R at location x.
3 Optimal Exponential Curves in SE (2) As mentioned in the introduction we distinguish between two approaches: a first order optimization approach based on a structure tensor on SE(2), and a second order optimization approach based on a Hessian on SE (2). A comparison of the two approaches (with the same amount of Gaussian blurring) on basic test images has been presented in previous work [33, ch:5.5.1, fig.5. 5-5.7]. These experiments show that for large curvatures the 2nd order approach is preferable, whereas for small curvatures the 1st order approach is preferable in terms of accuracy and stability. A formal theoretical exponential curve optimization formulation of the first order approach (as done for the second order approach in [26] based on a Hessian) via the structure tensor is missing both in [33] and in [37], and is presented here in the next subsection. It also serves as an introduction to its more technical SE(3)extension in Section 4.2.

Exponential Curve Fits in SE(2) of the 1st Order
According to 1 in [33, 5.3.2] the 3×3 structure matrix at g = (x, R θ ) ∈ SE(2) is, in the fixed coordinate system given by with s p , ρ p > 0 isotropic scales on R 2 and s o , ρ o > 0 scales on S 1 of separable Gaussian kernels G ρp,ρo (x, θ) = G ρp (x)G ρo (θ), and separable G sp,so on R 2 × S 1 . We assume that U , ρ 0 , ρ p , and g, are chosen such that S sp,so,ρp,ρo (g) is a non-degenerate matrix.

Remark 4
The scalings are naturally included in the gradient, since by definition the gradient ∇U is the Riesz representation vector of derivative dU and indeed where we recall Eq. (11). By Eq. (4), (79) we have Thus the structure tensor field on SE(2) equals 2 S sp,so,ρp,ρo (g) = Now finding the normalized eigenvector c * (g) with smallest eigenvalue of this structure matrix S sp,so,ρp,ρo (g) solves an optimization problem. Akin to (20) we shall rely on parallel transported exponential curve with c = (c 1 , c 2 , c 3 ) T ∈ R 3 , and for all g = (x, R) and all h = (x , R ), which starts at neighboring location h ∈ SE(2) (at t = 0) and which carries the same angular velocity and spatial velocity as the exponential curve t → γ g h (t) departing from g.

Theorem 1
The normalized eigenvector M µ c * (g) with smallest eigenvalue of the rescaled structure matrix M µ S sp,so,ρp,ρo (g)M µ provides the solution c * (g) of the following optimization problem: Proof By the construction of (24) and the fundamental property (8) of exponential curves we have that with γ c h,g (0) = h, and where because of normalization c µ = 1, t denotes the Riemannian arclength in Riemannian manifold (SE(2), G µ ). As a result the convex optimization functional can be written as whereas the boundary condition can be written as Euler-Lagrange gives ∇E(c * ) = λ 1 ∇ϕ(c * ), with λ 1 the smallest eigenvalue of M µ S sp,so,ρp,ρo (g)M µ and we have from which the result follows.

Exponential Curve Fits in SE(2) of the 2nd Order
We now discuss the second order optimization approach based on a Hessian matrix. At each g ∈ SE(2) we define a 3 × 3 non-symmetric Hessian matrix with V = G sp,s0 * U , and where i denotes the row index and where j denotes the column index, and where both G sp,so and G ρp,ρo are Gaussian kernels with isotropic spatial part so one can use 'flat' R 2 × S 1 -convolutions or 'shift-twist' SE(2)-convolutions.
Remark 5 In contrast to the structure tensor in the 1st order exponential curve fit, the mapping U → H sp,so,ρp,ρo in the 2nd order approach is linear, and thereby H ρp,ρo,sp,so = H ρp+sp , ρo+so ,0,0 , and thus we may exclude the (ρ p , ρ o ) and write H sp,so := H 0,0,sp,so .
Later, in a non-linear extension of the 2nd order approach the parameters will appear again.
Remark 6 As the Hessian matrix is non-symmetric, we will consider two optimizations next. The first one, in Theorem 2 is linear in U , and involves the symmetric sum. The second one, in Theorem 3, is non-linear in U and involves the symmetric product.
Theorem 2 Let g ∈ SE(2) be such that the eigenvalues of the symmetrized Hessian have the same sign. Then the normalized eigenvector M µ c * (g) with smallest eigenvalue of the symmetrized Hessian matrix provides the solution c * (g) of the following optimization problem: Proof Similar to the proof of Theorem 1 we have: In previous work of the first author and Franken [34] best exponential curve fits were based on Hessians of the modulus of an orientation score in order to have phase invariance across lines. In retrospect, this choice may not have been optimal as this imposes little variations across lines and if s p > 0 is chosen too small the eigenvector with minimal eigenvalue of the Hessians can point in the wrong direction. On top of this, extensive experiments with the Hessian based on the modulus U = |W ψ f | of orientation score W ψ f : SE(2) → C, did show us that the second order approach in Theorem 2 worked less effective than the following approach (based on the symmetric product).

Remark 7
In the next section we will generalize results to the case d = 3, i.e. R 3 S 2 → SE(3). It turns out that Theorem 1 and Theorem 2 have straightforward generalizations to d = 3, whereas we did not manage to generalize the approach in Theorem 3 to d = 3 in such a way that the group quotient structure is respected.

Optimal Exponential Curves in SE(3)
In this section we will generalize the best exponential curve fit theory from the preceding chapter on SE(2) to SE(3). Things become more technical both from the computational point of view and from the conceptual point of view (we have to deal with R 3 S 2 : (2)). Therefore, we will start in Subsection 4.1 with some notations on the basic tools of exponential curves, deviation from horizontality, (sub-)Riemmanian metrics and (sub-)Riemannian parametrization. These tools and preliminaries are used both in our best exponential curve fit theory of the first order on SE(3) (outlined in Subsection 4.2), and in our best exponential curve fit theory of the second order (outlined in Subsection 4.3). In Subsection 4.2 we set up the structure tensor on SE(3) in Subsection 4.2.1. This structure tensor is a key ingredient in Subsection 4.2.2 where we solve the variational formulation of best exponential curve fits of the first order to dataŨ : SE(3) → R and compute their minimizers by spectral decomposition of this structure tensor. In Subsection 4.2.3 we present a two-step approach for achieving torsion-free exponential curve fits respecting the quotient structure (1).
In Subsection 4.3 we first introduce a Hessian on SE(3) in Subsection 4.3.1. In contrast to the SE(2) case we only consider the variational formulation of best exponential curve fits of the second order to dataŨ : SE(3) → R that are solved by spectral decomposition of the symmetric sum of the Hessian. Again torsion-free exponential curve fits are accomplished via a two-step approach in Subsection 4.2.3. (3). (2) is identified with all rotations about reference axis a = e z = (0, 0, 1) T . Two rigid body motions

Preliminaries and Notations for Best Exponential Curve Fits in SE
where R a,φ ∈ SO(3) denotes a counterclockwise rotation about axis a with angle φ. Legal operatorsΦ : with right regular, and left regular representations: for all h, g ∈ SE(3). They relate one-to-one to operators Φ : Tildes indicate we are operating in the extended space SE(3) of rigid body motions: where R n is any rotation that maps e z onto n ∈ S 2 . The group SE(3) acts on R 3 S 2 from the left via and thereby we associate to each curveγ(t) = (x(t), Exponential curves in SE(3) are denoted by: where {A i } form a basis for the Lie algebra These exponential curves are auto-parallel w.r.t. (left) Cartan connection 3 [29, Thm 12, App.C], i.e.
Intuitively speaking, they are the "straight curves" in a curved geometry on SE (3), and they are parameter- ; (c 4 , c 5 , c 6 )) T ∈ R 6 , as their tangents have constant components w.r.t. the frame of left-invariant vector fields (obtained by pushforward (L g ) * of left-multiplication L g h = gh): It can be shown that [18,29] A and explicit formula's in several coordinate charts are given in [29,Eq. 24,25], [17].
In order to normalize the speed of which we move along our exponential curves we use the following (leftinvariant) metric tensor i A i |γ, and with stiffness parameter µ. Now, for exponential curves, one hasγ i = c i is constant, so we again impose the normalization constraint in order to ensure that our exponential curves are always parameterized by Riemannian arclength t. The precise formulae for exponential curves in SE (3) are well-known, see e.g. [29,Eq.54]. Their spatial part are circular spirals with constant curvature and torsion 3 Here we rely on a (left) Cartan connection carrying nonzero torsion [27, Thm 3.9 & App. B ]. This (left) Cartan connection is not to be mistaken with torsion-free biinvariant Cartan-Schouten connections on Lie groups, which also have various interesting applications in image analysis and statistics on Lie groups, cf. [51,52]. Within our orientation score framework, considered in the application section, right-invariance is undesirable, see [27,Lemma 3.3]. The (left) Cartan-connection also underlies "the shortest curves" i.e. (sub)-Riemannian geodesics within SE(d) [24,30].
with torsion τ (t) = |c (1) · c (2) |κ(t) and curvature Another feature is deviation from horizontality: They tell us how much our Riemannian manifold structure spatially deviates from sub-Riemannian manifold . [53,11,19,24] and along precisely those curvesγ = (x(·), R(·)) one has the natural correspondence between position and orientation parṫ where s denotes spatial arc length and t denoting sub-Riemannian arc length, with with curvature magnitude κ(s) = ẍ(s) . In particular, for horizontal exponential curves one has However, it is not always wise to enforce (37) as can be deduced from Fig. 8. In this Figure, we visualize a distribution U : R 3 S 2 → R + via 'glyph visualization'. If U : R 3 S 2 → C is not strictly positive we perform a glyph field visualization of the absolute value |U |.
for some y ∈ R 3 , and some suitably chosen ν > 0.
In Fig. 9 we show examples of how best exponential curve fits provide more suitable frames for line enhancement and tracking in distributions U : R 3 S 2 → R + .

Exponential Curve Fits in SE(3) of the 1st Order
Now let us generalize Definition 1 to the setting of R 3 S 2 . The first problem that arises is that many exponential curves are not well-defined on the quotient R 3 S 2 , Eq. (1). Also the left-invariant gradient, computed via M µ = diag{µ i } := diag{µ, µ, µ, 1, 1, 1}: is not well-defined on the quotient. Therefore, we must resign to the embedding space of SE(3) first.
Using Eq. (39), we rewrite the exponential curves in SE (3), Eq. (29), as The real part of the orientation score U : R 3 S 2 → R (using a cake-wavelet described in [43]) provides us a density from which exponential curve fits are computed via the algorithm in Section 4.2.3. Top: spatial parts of exponential curves aligned with spatial genera- These exponential curves are aligned with the highest angular response in the orientation score only and lack data-adaptation. In contrast to our best exponential curve fit t →γ (63), whose spatial parts are depicted as black lines.
constraining the gradient to unity element e = (0, I).
Firstly, both smoothing kernels are given bỹ with G R 3 sp the heat kernel on R 3 and G S 2 so the heat kernel on S 2 , yielding separability and SO(2)-invariance: The invariance is needed for well-posed convolution operators on R 3 S 2 [28, Cor.1], as we will employ soon. Secondly, regarding Eq. (41), we stress thatγ c h,g (t) is defined as follows.
is the unique exponential curve starting from h which is parallel transported from exponential curveγ c g (t), in the sense that it has both the same spatial velocity and the same angular velocity. This is similar to Eq. (20) where only spatial velocity is preserved. The next theorem tells us more explicitly what this means both on Lie group and Lie algebra level. Note thatγ c g,g =γ c g .
departing from g = (x, R), as defined above, then we have and we have the following fundamental relations Proof See Appendix B.
See Fig. 10 for an intuitive illustration of the theorem.

Remark 8
In order to achieve parallel transport of exponential curves in SE(3) one applies transformation in the Lie algebra. Indeed such a transformation preserves the left-invariant metric: for all h ∈ SE(3) and all t ∈ R. For further differential geometrical details see Appendix B.
Now that we have taken into account parallel transportation of curves in SE(3) let us include the appropriate structure tensor into problem (41) of finding best exponential curve fits of the first order.

Remark 9
Gaussian gradient (48) generalizes (16), and it respects the non-commutative group structure thanks to the isotropy of the spatial diffusion.
Remark 10 By construction (27) and (31) Remark 11 We assume that s p , s o > 0 and functionŨ are chosen in such a way that the null space of the structure matrix is precisely equal to N (and not larger).
Due the assumption in Remark 11 we need to impose the condition in our best exponential curve optimization to avoid nonuniqueness of solutions.
Theorem 5 Problem (41) can be rewritten as Proof By direct computations (where we recall that the gradient by definition is the Riesz representation vector of the exterior derivative d with respect to metric tensor G µ (39)) we have This yields providing us the final result.
So we have the following Euler-Lagrange equations.

Exponential Curve Fits in
In reducing the problem to R 3 S 2 we first note that with R ez,α ∈ SO(3) the counterclockwise rotation about e z where from now on h α := (0, R ez,α ). Eq. (54) follows from ∇Ṽ (gh α ) ≡ Z T α ∇Ṽ (g). We deduce that for all α ∈ [0, 2π), g ∈ SE (3). As a result we have and thereby we obtain: and we can define the following curve on R 3 S 2 : , as (57) is independent of the choice of R n (denoting any rotation which maps e z onto n).
Theorem 6 The structure tensor defined by (46) can be expressed as S sp,so,ρp,ρo (x, R n ) = 2π with R T n,n = (R T n R n ⊕ R T n R n ) ∈ SO(6). Its eigenvector c * (x, R n ) with smallest eigenvalue provides the solution of (41) and defines a curve in R 3 S 2 : Thereby integration over third Euler-angle α is no longer needed in the definition of the structure tensor (41) as it just produces a constant 2π factor. Using (42) we havẽ n) where we use isotropy of G R 3 sp on R 3 . Furthermore, apply Corollary 1. Finally, by Eq. (56) optimal curve (59) is independent of the choice of R n ∈ SO(3). .

Torsion-free Exponential Curve Fits of the 1st Order via a Two-step Approach
Theorem 5 provides us best exponential curve fits that possibly carry torsion. On the one hand, from Eq. (34) we deduce that along an exponential curve one has τ = (c 1 c 4 +c 2 c 5 +c 3 c 6 )κ. On the other hand, we need to exclude the null space N from our optimization domain and by Eq. (51) this is achieved by including constraint c 6 = 0. As a result we are insisting on zero torsion along horizontal exponential curves where c 1 = c 2 = 0. Along other exponential curves torsion appears if c 1 c 4 +c 2 c 5 = 0. On top of this, torsion is a higher order less-stable feature than curvature. Therefore we would like to exclude it altogether from our exponential curve fits presented in Theorem 5 and Theorem 6, by a different theory and algorithm. The results of the algorithm clearly show that even if structures do have torsion, the local best exponential curve fits do not need to to carry torsion in order to achieve good results in local frame adaptation, see e.g. Fig. 9.
The constraint of zero torsion forces us to split our best exponential curve fit algorithm into two parts: 1. Estimate at g ∈ SE(3) the spatial velocity part c (1) (g) from the spatial structure tensor. 2. Move to a different location g new ∈ SE(3) where a horizontal exponential curve fit makes sense and then estimate the angular velocity c (2) from the rotation part of the structure tensor over there.
This forced splitting is a consequence of the next lemma.
Lemma 1 Among the class of exponential curves with nonzero spatial velocity c (1) = 0 such that their spatial projections do not have torsion, we have that the constraint c 6 = 0 does not impose constraints on curvature if and only if the exponential curve is horizontal.
Proof For a horizontal curveγ c g (t) we have d H = 0 ⇔ c 1 = c 2 = 0 and indeed τ = c (1) · c (2) κ = (c 3 c 6 )κ = 0 and we see that constraints c 6 = 0 and τ = 0 reduce to only one constraint (as c (1) = 0 ⇔ c 3 = 0). The curvature magnitude stays constant along the exponential curve and the curvature vector at t = 0, recall Eq. (35), is in this case given by which can indeed be any vector orthogonal to spatial velocity c (1) = (0, 0, c 3 ) T . Now let us check whether the condition is necessary. Suppose t →γ c g (t) is not horizontal, and suppose it is torsion free with c 6 = 0. Then we have c 1 c 4 + c 2 c 5 = 0, as a result the initial curvature is both orthogonal to vector c (1) = (c 1 , c 2 , c 3 ) T and orthogonal to (−c 2 , c 1 , 0) T , and thereby constrained to a one dimensional subspace.
Corollary 2 In order to allow for all possible curvatures in our torsion-free exponential curve fits we must relocate the exponential curve optimization at g ∈ SE(3) in the scoreŨ : SE(3) → R to a position g new ∈ SE (3) where an horizontal exponential curve can be expected. Subsequently, we can use Theorem 4 to parallel transport the horizontal and torsion-free curve through g new to a torsion-free exponential curve through g.
Step 2: Find the optimal spatial velocity: which boils down to finding the eigenvector with minimal eigenvalue of the 3 × 3 spatial sub-matrix of the structure tensor (46).
Step 3: Given c (1) (g) we aim for an auxiliary set of coefficients, where we also take into account rotational velocity. To achieve this in a stable way we will transfer to a different location in the group: and apply parallel transport (of Theorem 4) afterwards. At g new , it is natural to enforce horizontality, and we consider auxiliary optimization problem where we note that zero deviation from horizontality (36) and zero torsion (34) is equivalent to d H = 0 and τ = 0 ⇔ c 1 = c 2 = c 6 = 0.
Step 4: The auxiliary coefficients c new (g new ) = (0, 0, c 3 (g new ), c 4 (g new ), c 5 (g new ), 0) T of a torsion-free, horizontal optimal exponential curveγ cnew gnew through g new . Now we apply parallel transport (via Theorem 4) of this horizontal optimal exponential curve towards the corresponding exponential curve through g: (62) This gives the final (not necessarily horizontal), torsionfree, best exponential curve fit t →γ

Remark 12
The parameters c * f inal (g) of the final best exponential curve fit follow by a two-step optimization and need not be equal to c * (g) in Theorem 5, where a single optimization over the full space SE(3) is done.

Exponential Curve Fits in SE(3) of the 2nd Order
In this section we will generalize Theorem 2 to the case d = 3, where again we include the restriction to torsionfree exponential curves.
Before continuing with this, we would like to stress that it is very hard to generalize the other approach presented in Theorem 3 to the case d = 3, as the corresponding non-symmetric Hessian on SE(3) contains a final zero row instead of a final zero column (that would correspond to N = span(0, 0, 0, 0, 0, 1) T ). This problem has been overlooked in previous work [29, ch:9].

The Hessian on SE(3)
We define the 6 × 6 non-symmetric Hessian matrix by withṼ =G sp,s0 * Ũ , and where i = 1, . . . , 6 denotes the row index and where j = 1, . . . , 6 denotes the column index. Note that the double scaling with the matrix M µ −2 is natural in view of (22).

Remark 13
To implement the Hessian efficiently we use c k ij A k , recall Remark 9. Furthermore, for the purpose of accuracy we apply Theorem 7 Let g ∈ SE(3) be such that the symmetrized Hessian matrix 1 2 M µ (H sp,so (g) + (H sp,so (g)) T )M µ has eigenvalues with the same sign. Then the normalized eigenvector M µ c * (g) with smallest absolute non-zero eigenvalue of the symmetrized Hessian matrix provides the solution c * (g) of the following problem: withṼ =G sp,so * Ũ .
Proof The proof is essentially the same as proof of Theorem 2 (now summations run from 1 to 5).

Remark 14
The restriction to g ∈ SE(3) such that the eigenvalues of the symmetrized Hessian carry the same sign is necessary for a unique solution of the optimization. Note that in case of our first order approach via the structure tensor, no such cumbersome constraints arise. In case g ∈ SE(3) is such that the eigenvalues of the symmetrized Hessian have different sign there are 3 natural options: 1. Move towards a neighboring point where the Hessian eigenvalues have the same sign and apply parallel transport (Theorem 4) of the best exponential curve fit at the neighboring point. 2. Construct within the set of multiple minimizers c * the one with minimal curvature via an additional Euler-Lagrange technique. 3. Take c * (g) still as the eigenvector with smallest absolute eigenvalue (representing minimal absolute principal curvature). This is no longer solves (65), but maximizes the energy in the orthogonal subspace.
Although it is not difficult to solve the equations in the second approach analytically, it does not seem to produce appropriate fits at curved structures. Therefore, we prefer the other options.

Now similar considerations as in Subsection 4.2.3 ap-
ply and again torsion-free exponential curves are obtained via the two-step approach. Here we follow the same algorithm as in Subsection 4.2.3, but now with the Hessian field H sp,so (64) instead of the structure tensor field. However, there is one additional complication to be taken into account compared to the first order approach via the structure tensor: the computation of mixed left-invariant derivatives for the horizontal 3 × 3-Hessian Here c (1) (g) is the spatial velocity vector that is found as the eigenvector with smallest eigenvalue of the leftinvariant Hessian computed with Gaussian spatial derivatives at scale s p . The solution of auxiliary problem (61) (with the Hessian H sp,so instead of structure tensor) now boils down to finding the eigenvector with smallest eigenvalue of the matrix in (66). The final solution is again found via Eq. (63).
So it just remains to compute (66). Expressed in standard Euler-angles R new = R ez,γ R ey,β R ex,α where we can set α = 0, we find that M hor µ H hor (g new )M hor with n = (sin β cos γ, sin β sin γ, cos β) T , which follows from the formulas of the left-invariant derivatives in standard Euler angles in [29,Eq. 24], that hold for β = 0, i.e. g new = (x, I), with x ∈ R 3 free. In the vicinity of g new = (x, I) one relies on different Euler angles, see [29, Eq. 5, Fig. 4].

Adaptation of Exponential Curve Fits of the 2nd Order via Factorization
Practical experiments (see e.g. Fig. 11) indicate it is helpful to use the basis symmetric Hessians where in every entry (both in the upper triangle and in the lower triangle) the angular derivative comes first. Such matrices have both a zero last column and zero last row as ∂ αṼ = 0. Again the null space is avoided by imposing c 6 = 0 in the minimization. However, the venture point is now different; instead of fitting exponential curves directly, we apply the following optimization to obtain c * (g).
Theorem 8 Let g ∈ SE(3) be such that Hessian matrix M µ (H(g))M µ has two eigenvalues with the same sign. Then the normalized eigenvector M µ c * (g) with smallest eigenvalue provides the solution c * (g) of the following optimization problem: Proof We have (for details see Appendix C) This can be decomposed in the two-step approach. Effectively, this means that in (66) the upper triangle of the Hessian is replaced by the lower triangle, whereas the lower triangle is maintained. See Fig. 11, where the results of the best exponential curve fits of second order (relying on Hessians) are similar to best exponential curve fits of first order (relying on the structure tensor) depicted in Fig. 9.

Embodiment of Applications
In Section 6 we present image analysis applications of the theory as developed according to the grand scheme depicted in Fig. 4. This requires specification of the orientation score W ψ , and, in particular, of the wavelet ψ, and the application dependent left-invariant operator Φ that together yield the image enhancing operator Υ acting on images f . To this end, we describe in this section the continuous-wavelet/coherent-state transform that we use, based on cake-wavelets appropriate for processing of disk-limited images f . Furthermore, we describe three types of left-invariant operators Φ, one for multiple-scale vessel detection, one for coherence enhancing diffusion via multiple-scale 2D orientation scores, and one for coherence enhancing diffusion via 3D orientation scores. Fig. 11: In black the spatially projected part of best local exponential curve fits t → γ c g (t) of the second kind (fitted to the real partŨ = Re{W ψ f } of 3D invertible orientation score W ψ f , cf. Fig. 12) via cake wavelet ψ, cf. In the image analysis applications discussed in this section our function U : R d S d−1 → R + is given by either the real part of the invertible orientation score, see Fig. 3 for d = 2, or the real part of a fixed scale layer of a continuous wavelet transform. An orientation score is given by where R n is any rotation mapping reference axis a onto n ∈ S d−1 , and where and for d > 2 we restrict ourselves to ψ satisfying and for all rotations R a,α ∈ Stab(a) (for d = 3 this means for all rotations about axis a, Eq. (3). As a result U is well-defined on the left cosets 1)) as the choice of R n ∈ SO(d) mapping a onto n is irrelevant. See Fig. 12 for an example of a 3D orientation score. We will rely on both the 2D 'cake-wavelets' used in [6,22] and their recent 3D-equivalents [42]. Detailed formulas and recipes to construct such wavelets efficiently can be found in [42], and are therefore omitted here. Nonetheless, in order to provide the global intuitive picture these wavelets are depicted in Fig. 13. In this article, in all experiments using 3D-invertible orientation scores, we used the following 3D cake-wavelet parameters: N 0 = 42, s φ = 0.7, k = 2, N = 20, γ = 0.85, L = 16 evaluated on a grid of 21x21x21 pixels, for details see [42]. Here we restrict ourselves to disk-limited images with ρ close to the Nyquist frequency. In this case admissible wavelets satisfy: uniformly for ω ∈ R d , for some a priori fixed δ, M > 0. Under this assumption, the condition number of the transform between image and orientation score equals [23,Thm.20] where the inequality is in fact an equality. Here the ideal case of condition number 1 is obtained if δ = M = 1, but if M/δ is not too large this is also fine.
For admissible wavelets exact reconstruction is performed via the adjoint:  Fig. 13: Visualization of cake-wavelets in 2D (top) and 3D (bottom). In 2D we fill up the 'pie' of frequencies with overlapping differentiable "pieces of cake", and application of an inverse DFT (for details see [6]) provides complex-valued wavelets whose real parts are line detectors, and whose imaginary parts are edge detectors. In 3D one must include anti-symmetrization and the Funk transform [21] on L 2 (S 2 ) to obtain the same, for details see [42]. The intuitive idea behind applying the Funk transform, is to redistribute spherical data from orientations towards great circles which lay in planes orthogonal to those orientations. After all we want the real part of our oriented wavelets to be line detectors (and not plate detectors) in the spatial domain.
This allows us to robustly relate operators Φ on orientation scores to operators on images Υ via where W * ,ext ψ is the natural extension of the adjoint W * ψ from the range R(W ψ ) to its embedding Hilbert space, for details see [22]. Recall Fig. 4. In general operator Φ does not map the space of orientation scores into itself, but that is not problematic as the orthogonal complement to the space of orientation scores is mapped to 0 by W * ,ext ψ . In order to obtain sensible Euclidean invariant operators Υ (i.e. operators Υ commuting with rotations and translation) operators Φ must be left-invariant and not right-invariant, see [23]. Operator Φ is left-invariant if Φ • L g = L g • Φ, for all g ∈ SE(d), where we recall L g is given by (26).
In the subsequent sections we consider two types of left-invariant operators: 1. for d = 2, 3, non-linear adaptive diffusions steered along the gauge frame {B 1 , . . . , B n d }, i.e.
where W (g, t), with t ≥ 0, is the solution of: where we recall {B 1 , . . . , B n d } denotes the gauge frame associated to locally best exponential curve fit to dataŨ at location g ∈ SE(d).
Occasionally, for d = 2 we extend the (all-scale) orientation scores to multiple-scale invertible orientation scores. Such multiple-scale orientation scores [55] coincide with wavelet transforms on the similitude group SIM (2) = R 2 SO(2)×R + , where one uses another Bspline [61,31] basis decomposition along the log-radial axis in the Fourier domain.
In our experiments for d = 2 we used N = 4, N = 12 or N = 20 orientation layers and a decomposition centered around M = 4 discrete scales a l given by l = 0, . . . , M − 1 where a max is inverse proportional to the Nyquist-frequency ρ n and a min close to the inner scale [32] induced by sampling. For details see [55], for the overall idea of splitting an invertible orientation score into multiple scales yielding a continuous wavelet transform on SIM (2), see Fig. 14. Such wavelet transform on W ψ f : SIM (2) → C is given by and we again set U := Re{W ψ f }. Then we define the total integrated multiple scale SIM (2)-vesselness by: where SE(2)-vesselness operator Φ is given by Eq. (74), and where µ ∞ and µ i,∞ denote maxima w.r.t. sup-norm · ∞ taken over the subsequent terms.

Image Analysis Applications
In case d = 2 and Φ is the SIM (2)-vesselness filter (77), we consider the application of enhancing and detecting the vascular tree structure in retinal images. Such image processing task is highly relevant as the retinal vasculature provides the only means for non-invasive observation of the human vascular system. A variety of eye-related and systematic diseases such as glaucoma, age-related macular degeneration, diabetes, hypertension, arteriosclerosis or Alzheimer's affect the vasculature and may cause functional or geometric changes [41]. Automated quantification of these defects promises massive screenings for systematic and eye-related vascular diseases on the basis of fast and inexpensive imaging modalities, i.e. retinal photography. To automatically assess the state of the retinal vascular tree, vessel segmentations and/or models have to be created and analyzed. Because retinal images usually suffer from low contrast on small scales, the vasculature needs to be enhanced prior to model creation/segmentation. One well-established approach is the Frangi vesselness filter [35]. It frequently occurs in rubust retinal vessel segmentation methods [14,48]. However, a known drawback of the Frangi filter is that it can not handle crossings or bifurcations that make up huge parts of the retinal vascular network. This is precisely where the orientation score framework and the presented locally adaptive frame theory comes into play. We will show SIM (2)-vesselness via multiple scale orientation scores has advantages over multiple-scale vesselness filtering acting directly in the image domain. Furthermore, we show that SIM (2)-vesselness (given by Eq. (77) and Eq. (74)) really benefits from including the left-invariant locally adaptive frame In case d = 2 and Φ = Φ t is a left-invariant diffusion steered by the gauge frame {B 1 , B 2 , B 3 } given by Eq. (72) and Eq. (73), we consider medical imaging examples where elongated structure need to be enhanced in noisy data. E.g. 2-photon microscopy 2Dimages containing collagen fibers which need crossing preserving enhancement, prior to coherence quantification, also investigated in 3D-imaging [36]. Here we show benefits of extending results [34,27] to a multiple scale setting, again employing locally adaptive frames in SE(2), adapted to each scale layer in the continuous wavelet transform (76).
In case d = 3 the envisioned applications include blood vessel detection in 3D MR-angiography, e.g. the detection of the Adamkiewicz vessel, relevant for surgery planning. Also extensions towards fiber-enhancement of diffusion-weighted MRI [29,28] the non-linear diffusions are of interest. Some preliminary practical results have been conducted on such 3D-datasets [42,22,20], but here we shall restrict ourselves to very basic artificial 3D-datasets to show proof of concept of our locally adaptive frames in the roto-translation space only, and leave these three applications for future work.

Improved Differential Invariants
As a basic example of a practical differential invariant [32,47,10,39], we consider our SIM (2)-vesselness (77). We will compare our method using the locally adaptive frame {B 1 , B 2 , B 3 }, to SIM (2)-vesselness relying on the left-invariant frame {A 1 , A 2 , A 3 }. Furthermore, we compare our method to the multiple scale vesselness filter acting directly in the image domain.
In order to perform a first, simple practical comparison, we test these 3 techniques on the publically available 4 High Resolution Fundus (HRF)-dataset [44], containing manually segmented vascular trees by medical experts. The HRF-dataset consists of wide-field fundus photographs for a healthy, diabetic retinopathy and a glaucoma group (15 images each).
These practical experiments on the HRF-dataset reveal two advantages, supporting practical relevance and proof of concept of our locally adaptive frame theory via invertible multiple scale orientation scores. In short, we show the following two advantages: For the comparison we rely on a basic standard morphological component segmentation that we will explain next. We could also have used more state-of-the-art vessel segmentation techniques, such as e.g. VETOS (Vessel Edge Tracking via Orientation Scores) [6] possibly followed by effective automatic optic disk detection [7], but in order to check for the effectivity of our SIM (2)extension of vesselness and the use of gauge frames, the simple morphological segmentation seems more suited here as this segmentation does not include any further contextual alignment modeling in SE (2). To show the advantages mentioned above, we devised a simple segmentation algorithm to turn a vesselness filtered image V(f ) into a segmentation (i.e. a binary map). First an adaptive thresholding is applied, yielding a binary image where Θ is the unit step function, G γ is a Gaussian of scale γ = 1 2 σ 2 1 and h is a threshold parameter. In a second step, the connected morphological components in f B are subject to size and elongation constraints. Components counting less than τ pixels or showing elongations below a threshold ν are removed. Parameters γ, τ and ν are fixed at 100 px, 500 px and 0.85 respectively. The vesselness map V(f ) : R 2 → R is either the multi-scale vesselness by Frangi et al. [35] or the SIM (2)-vesselness (77).
The segmentation algorithm described above is evaluated on the HRF dataset. Average sensitivity and accuracy over the whole dataset are shown in Fig. 18 as a function of the threshold value h.
Regarding the first advantage, our method via invertible scale-orientation scores performs considerably better than the one based on the multi-scale Frangi filter. The segmentation results obtained with SIM (2)vesselness (77) are more stable w.r.t variations in the threshold h and the performance on the small vasculature has improved as measured via the sensitivity. Average sensitivity, specificity and accuracy at a threshold of h = 0.05 resp. given by 0.786, 0.988, 0.969 (healthy), 0.811, 0.963, 0.953 (diabetic retinopathy) and 0.797, 0.976, 0.964 (glaucoma) compare well with other segmentation methods evaluated on the HRF dataset for the healthy cases (see [14,Tab. 5]). On the diabetic retinopathy and glaucoma group, our method even outperforms existing segmentation methods. Fig. 18 shows a full segmentation computed with the proposed method together with a more detailed patch. In Fig. 16 we see that our method improves sensitivity both at non-crossing structures (due to line propagation of the anisotropic wavelets ψ in the wavelet transform (76)) and at crossing/bifurcating structures. As expected, we see a larger improvement at crossings.
Finally, regarding the second advantage we refer to Fig. 15 and Fig. 19, where the SIM (2)-vesselnessfiltering via the locally adaptive frame produces a visually much more appealing soft-segmentation of the blood vessels than SIM (2)-vesselness filtering via the non-adaptive frame. It therefore also produces a more accurate hard-segmentation as can be deducted from the comparison in Fig. 18. For comparison, the multiscale Frangi vesselness filter is also computed via summation over single scale results and max-normalized. Generally, we conclude from the experiments that the locally adaptive frame approach better reduces background noise, showing much less false positives in the final segmentation results. See the typical segmentation results on relatively challenging patches in Fig. 17.

Improved Crossing-Preserving Flows on Continuous Wavelet Transforms
Here, we briefly show a multiple-scale extension of coherence enhancing diffusion [34] via multiple scale orientation scores. In fact, instead of constructing an invertible orientation score as done in [34] we apply the continuous wavelet transform (76) where we split the cake-wavelets in a B-spline basis along the log radial axis, cf.

Experiments in 3D Images
We now show first results of the extension of coherence enhancing diffusion via orientation scores (CEDOS [34]) to the 3D setting. Again, data is processed according to Fig. 4. First, we construct an orientation score according to (68), using the 3D cake wavelets (Fig. 13). For determining the gauge frame {B 1 , . . . , B 6 } we use the first order structure tensor method in combination with Eq. (79) in Appendix A. In CEDOS we have Φ = Φ t , as defined in (72) and (73), which is a diffusion along the gauge frame. In comparison to 2D CEDOS we did not yet include adaptivity (via structureness [34]) in the diffusion matrix.
These diffusion operators can enhance elongated structures in 3D data, see Fig. 21 and Fig. 23, while preserving the crossings as can be seen in Fig. 22.      38)) of the absolute value of the diffused orientation scores with and without the use of gauge frames in the invertible orientation score of the artificial dataset depicted in Fig. 21. Diffusion along the gauge frames include better adaptation for curvature. This is mainly due to the angular part in the B 3 -direction, cf. Fig. 24, which includes curvature, in contrast to A 3 -direction. The angular part in B 3 causes some additional angular blurring leading to more isotropic glyphs.

Conclusion
Locally adaptive frames ('gauge frames') on images based on the gradient, structure tensor or Hessian are illposed at the vicinity of complex structures such as bifurcations and crossings. Therefore we create locally adaptive frames within invertible orientation scores to allow for curvature adaptation and adaptation for deviation from horizontality. This gives rise to a whole family of local frames (indexed by S d−1 ) per position in the image domain, enabling us to deal with crossings and bifurcations. In order to generalize gauge frames in the image domain to gauge frames in the orientation scores, we rely on local best exponential curve fits at each position and orientation within the real part of the orientation score. In Appendix A we have shown how such a single best exponential curve fit gives rise to a suitable locally adaptive frame { B 1 | g , . . . B n d | g } in T g (SE(d)). We distinguished between best exponential curve fits of the first order and second order: 1. Along first order exponential curve fits, the first order variation of the data (on SE(d)) along the exponential curve is locally minimal. The Euler-Lagrange equations are solved by finding the eigenvector of the structure tensor with smallest eigenvalue. 2. Along second order exponential curve fits, a second order variation of the data (on SE(d)) along the exponential curve is locally minimal. The Euler-Lagrange equations are solved by finding the eigenvector of the corresponding Hessian with smallest eigenvalue.
In case d = 3, i.e. in SE (3), these two approaches are presented for the first time. Closer inspection revealed it is natural to include a restriction to torsion-free exponential curve fits in order to be both compatible with the null-space of the structure/Hessian tensors and the quotient structure of R 3 S 2 embedded in SE (3).
We have presented algorithms following an efficient two-step approach to compute such torsion-free exponential curve fits. Experiments on artificial datasets show that even if the elongated-structures have torsion, the frame { B 1 | g , . . . , B n d | g } is well-adapted to the local structure at g ∈ SE(d). This adequate adaptation of the frame field {B 1 , . . . , B n d } allows us to perform crossing-preserving coherence enhancing diffusion via invertible orientation scores (CEDOS). Regarding the CEDOS algorithm, we managed to improve it to a multiple scale setting for d = 2, dealing with (crossing) multiple scale structures, construct and implement CEDOS, for the first time, for 3D images (d = 3).
In the application section we performed experiments on biomedical images containing crossing collagen fibers to show advantages of our multiple scale extension of CE-DOS [34]. Finally, we considered the application of a differential invariant (SIM (2) As such, the generic applicability and advantage of locally adaptive frames in orientation scores in 2D medical imaging has been shown. In this work, the advantages of including locally adaptive left-invariant frames into adaptive diffusions (extending CEDOS [34] to 3D) in 3D-imaging has only been shown on artificial datasets. Therefore, in future work we will study the use of locally adaptive frames in real 3D medical imaging applications, e.g. in 3D MR angiography [43], or in improved enhancement and detection of Adamkiewicz vessel [22]. tive frames in SE (3)   be an exponential curve through g that fits dataŨ : SE(d) → R best at g ∈ SE(d) in Lie group SE(d) of dimension n d = d(d + 1)/2. In Section 3 (d = 2), and in Section 4 (d = 3), we provide theory and algorithms to derive such curves. In this section we assume γ c g (·) is given.
Let us write c = c (1) c (2) ∈ R n d , with spatial velocity c (1) ∈ R d and rotational velocity c (2) ∈ R r d , with r d = dim SO(d) = d(d − 1)/2 of the exponential curve γ c g (·). Akin to the case d = 2 discussed in the introduction we define the Gauge frame via Eq. (14), but now with B = (B 1 , . . . , B n d ) T , A = (A 1 , . . . , A n d ) T , M µ = µI d ⊕ I r d ∈ R n d ×n d , R c = R 2 R 1 ∈ SO(n d ).
(79) Again R 1 is the counter-clockwise rotation that rotates the spatial reference axis a 0 , recall our convention (3), onto µ c (1) a c (2) strictly within the 2D-plane spanned by these two vectors. Rotation R 2 is the counter-clockwise rotation that rotates µ c (1) a c (2) onto µc (1) c (2) strictly within the 2D-plane spanned by these two vectors, so that a 0 Concluding, the preferred spatial direction a 0 ·A is mapped onto a 0 · B = c · A. Furthermore, all other spatial generators stay in the spatial tangent bundle (and they remain unchanged iff the exponential curve γ c g is a horizontal curve). The next theorem shows us that our choice of assigning an entire gauge frame to a single best exponential curve fit is the right one. For explicit formulae of the left-invariant vector fields in the d-dimensional case we refer to [30], and for d = 3 see [29].
Theorem 9 Let c(g) be the local tangent of best exponential curve fit t →γ c g (t) at g ∈ SE(d) in the datã U (g) = U (g (0, a)), recall (28). Consider the mapping of the frame of left-invariant vector fields A| g to the locally adaptive frame: with R c = R 2 R 1 ∈ SO(n d ), with subsequent counter-clockwise planar rotations R 1 , R 2 given by (80). The mapping A g → B g has the following properties: -The main spatial tangent direction (L g ) * a 0 · A| e is mapped to best exponential curve fit direction c T (g)· A| g . -Spatial left-invariant vector fields that are G µ -orthogonal to this main spatial direction stay in the spatial part of the tangent space T g (SE(d)) under rotation R c and they are invariant up to normalization under the action (81) iff the best exponential curve fit is horizontal.
Proof Regarding the first property we note that (L g ) * a 0 · A| e = a 0 · A| g as left-invariant vector fields are obtained by push-forward of the left multiplication. Furthermore, by Eq. (81) and Eq. (80) we have Regarding the second property, we note that if b · a = 0 ⇒ andγ c g is horizontal iff c (1) c (1) = a in which case the planar rotation R 2 reduces to the identity and R 2 R 1 b 0 = b 0 T and only spatial normalization by µ −1 is applied.
Remark 15 For d = 2 and a = (1, 0) T the above theorem can be observed in Fig. 7, where main spatial direction A 1 = cos θ∂ x + sin θ∂ y is mapped onto B 1 = c · A and where A 2 is mapped onto B 2 = µ −1 (− sin d H A 1 + cos d H A 2 ).
Remark 16 For d = 3 and a = (0, 0, 1) T the above theorem can be observed in Fig. 24, where main spatial direction A 3 = n · ∇ R 3 is mapped onto B 3 = c · A, and where A 2 and A 3 are mapped to the strictly spatial generators B 2 and B 3 . For further details see [43, ch:6.4.1].
Remark 17 For d = 3 and a = (0, 0, 1) T the above theorem can be observed in Fig. 24, where main spatial direction A 3 = n · ∇ R 3 is mapped onto B 3 = c · A, and where A 2 and A 3 are mapped to the strictly spatial generators B 2 and B 3 . For further details see [43, ch:6.4.1].

Remark 18
It is important that the same µ is used both in the optimization procedures of Theorems 1, 2, 3, 5, 7, 8, and in the construction (79) of orthonormal gauge frame {B 1 , . . . , B n d }, n d in this appendix. This is consistent from the geometrical viewpoint. Furthermore all optimization procedures are of the type arg min

B Proof of Theorem 4
By definitionγ c h,g is the parallel transported curve fromγ c g s.t. it has the same spatial velocity and the same rotational (and angular) velocity for all times, i.e.: x g (t) =ẋ h (t), with Ω = c (2) · ∇| e . Consideringγ c h,g (0) = h = (x , R ) and γ c g (0) = g = (x, R) we deduce that x h (t) − x g (t) = x h (0) − x g (0) = x − x, (R g (t)) −1 R h (t) = (R g (0)) −1 R h (0) −1 = R −1 R from which (44) can be deduced. Furthermore, by separate application of standard parallel transport on both the spatial and the rotational part of the tangent bundle T (SE(3)) we have with unity element e = (0, I), and where we stress that we have equivalence (via standard parallel transport on T (R 3 ) and T (SO(3))). Therefore we have for initial velocity and initial rotational velocitẏ (Rc (2) ) i A i ≡Ṙ g (0), and thereby we find (2) for the exponential curveγ c h,g starting from h, and therebỹ γ c h,g (t) = h e t((R ) T Rc (1) ,(R ) T Rc (2) )·M µ 2 ∇| e . The final equality follows by direct computation (i.e. application of Eq. (2) and (x, R) −1 = (−R −1 x, R −1 )) which shows that (44) is equivalent to (45). The top black curve is the spatial projection ofγ c g (·), and the bottom black curve is the angular projection of the exponential curve. After application of R T 2 the exponential curve is horizontal w.r.t. the intermediate frame {Â i }, subsequently R T 1 leaves spatial generators orthogonal to a horizontal curve invariant, so thatÂ 1 = B 1 andÂ 2 = B 2 are strictly spatial, as is in accordance with Theorem 9. The angular part of B 3 is visually included via the curvature at the blue arrow. The spatial part of B 4 and B 5 is depicted in the center of the ball.