Morphological and Linear Scale Spaces for Fiber Enhancement in DW-MRI

We consider morphological and linear scale spaces on the space ℝ3⋊S 2 of 3D positions and orientations naturally embedded in the group SE(3) of 3D rigid body movements. The general motivation for these (convection-)diffusions and erosions is to obtain crossing-preserving fiber enhancement on probability densities defined on the space of positions and orientations. The strength of these enhancements is that they are expressed in a moving frame of reference attached to fiber fragments, allowing us to diffuse along the fibers and to erode orthogonal to them. The linear scale spaces are described by forward Kolmogorov equations of Brownian motions on ℝ3⋊S 2 and can be solved by convolution with the corresponding Green’s functions. The morphological scale spaces are Bellman equations of cost processes on ℝ3⋊S 2 and we show that their viscosity solutions are given by a morphological convolution with the corresponding morphological Green’s function. For theoretical underpinning of our scale spaces on ℝ3⋊S 2 we introduce Lagrangians and Hamiltonians on ℝ3⋊S 2 indexed by a parameter η∈[1,∞). The Hamiltonian induces a Hamilton-Jacobi-Bellman system that coincides with our morphological scale spaces on ℝ3⋊S 2. By means of the logarithm on SE(3) we provide tangible estimates for both the linear- and the morphological Green’s functions. We also discuss numerical finite difference upwind schemes for morphological scale spaces (erosions) of Diffusion-Weighted Magnetic Resonance Imaging (DW-MRI), which allow extensions to data-adaptive erosions of DW-MRI. We apply our theory to the enhancement of (crossing) fibres in DW-MRI for imaging water diffusion processes in brain white matter.


Introduction
Diffusion-Weighted Magnetic Resonance Imaging (DW-MRI) involves magnetic resonance techniques for noninvasively measuring local water diffusion in tissue. Local water diffusion profiles reflect underlying biological fiber structure. For instance in the brain, diffusion is less constrained parallel to nerve fibers than perpendicular to them.
The diffusion of water molecules in tissue over time t is described by a transition density function p t , cf. [6]. Diffusion Tensor Imaging (DTI), introduced by Basser et al. [11], assumes that p t can be described for each position y ∈ R 3 by an anisotropic Gaussian. If {Y t } denotes the stochastic process describing the movement of water-molecules in R 3 , then one has p Y t = y | Y 0 = y = e − (y −y) T (D(y)) −1 (y −y) 4t (4πt) 3 2 √ det(D(y)) , (1) where D is a tensor field of positive definite symmetric tensors on R 3 estimated from the MRI data. In a DTIvisualization one usually plots the surfaces (the so-called DTI-glyphs) where μ > 0 is fixed and y ∈ Ω with Ω some compact subset of R 3 . The corresponding probability density U : R 3 × S 2 → R + on positions and orientations is given by with n ∈ S 2 = {n ∈ R 3 | n = 1} and y ∈ R 3 . Here ∞ 0 p(Y t = y + ρn | Y 0 = y)ρ 2 dρ denotes the Orientation Density Function (ODF), cf. [1,26], and we used the spatial probability density proportional to the volume of the DTI-glyphs. Due to the modeling assumption in Eq. (1) DTI is not capable of representing crossing fibers [6].
High Angular Resolution Diffusion Imaging (HARDI) is another recent DW-MRI technique for imaging water diffusion processes in fibrous tissues, cf. [29,61]. HARDI provides for each position in R 3 and for each orientation in S 2 an MRI signal attenuation profile, which can be related to the local diffusivity of water molecules in the corresponding direction. As a result, HARDI images are distributions (y, n) → U(y, n) over positions and orientations. HARDI is not restricted to functions on S 2 induced by a quadratic form and is thus capable of reflecting crossings.
We visualize these distributions U : R 3 × S 2 → R + as follows. Definition 1 A glyph visualization of the distribution U : R 3 × S 2 → R + is a visualization of a field y → S μ (U )(y) of glyphs, where each glyph is given by the surface S μ (U )(y) = y + μU (y, n)n | n ∈ S 2 ⊂ R 3 , for some y ∈ R 3 , and some suitably chosen μ > 0.
See Fig. 1, where a HARDI data set is depicted using a glyph visualization. In HARDI modeling the Fourier transform of the estimated transition densities is typically considered at a fixed characteristic radius known as the b-value, cf. [29].
In the remainder of this article we will visualize all DW-MRI images (including DTI) via glyph visualizations as defined in Definition 1.
To reduce noise and to infer information about fiber crossings, contextual information can be used to enhance the data. This enhancement is useful both for visualization purposes and as a preprocessing step for other algorithms, such as fiber tracking algorithms, which may have difficulty in noisy and/or incoherent regions. Recent studies indicate the increasing relevance for enhancement techniques in clinical applications [75,76,87].
Promising research has been done on constructing diffusion/regularization processes on the 2-sphere defined at each spatial locus separately [28,29,41] as an essential preprocessing step for robust fiber tracking. In these approaches position space R 3 and orientation space S 2 are decoupled, and diffusion is only performed over the angular part, disregarding spatial context. Consequently, these methods cannot propagate fiber fragments through complex fiber structures such as crossings, bifurcations and interruptions. These are precisely the interesting locations where the initial data (at least in case of DTI) fails to represent the underlying fiber structure.
Therefore, in contrast to previous work on enhancement of DW-MRI [18,20,28,29,41,45,46,80], we consider both the spatial and the orientational part to be included in the domain. That is, we consider a HARDI dataset as a function U : R 3 × S 2 → R + . Furthermore, we explicitly employ the proper underlying group structure, that arises by embedding the coupled space of positions and orientations R 3 S 2 := SE(3)/ {0} × SO (2) as the partition of left cosets into the group SE(3) = R 3 SO(3) of 3D-rigid motions. This group structure is necessary for alignment of fiber fragments. The general advantage of our approach on SE(3) is that we can enhance the original HARDI/DTI data using orientational and spatial neighborhood information simultaneously.
This allows us to extrapolate crossings in distributions on R 3 S 2 created from the original DTI data to [65] via Eq. (3). This may be used to reduce the number of scanning directions in areas where the (hypo-elliptic) diffusion [36,Chap. 4.2] on R 3 S 2 yield reasonable fiber extrapolations, cf. [36,65,67,68]. This can also be observed in Fig. 5 where the black-boxes indicate areas where fibers of the corpus callosum and the corona radiata cross.
HARDI already produces more detailed information about complex fiber-structures. Application of the same diffusion on HARDI [68] then removes spurious crossings that are not aligned with the surrounding glyphs, see

Contextual Processing of DW-MRI
Our diffusions relate to Brownian motion of oriented random walkers. If we constrain each random walker say at (y, n) to either proceed spatially in direction n or to change its orientation n we get Brownian motions with hypo-elliptic generators, with a constrained coupling between spatial and angular movements. Figure 2 shows that such a hypo-elliptic diffusion preserves the fiber-structure better than spatial/angular diffusion. The key idea is that our hypo-elliptic diffusions naturally include the context of fiber-fragments in the data. These diffusions are solved by shift-twist convolutions with the corresponding Green's functions, cf. [36]. This generalizes the well-known, efficient framework of tensor voting for contextual processing of tensors [57,58] in a PDE (scale space) setting on R 3 S 2 . Tensor voting applies a shift-twist convolution with a (singular) kernel representing alignment with an a priori voting field [43,Chap. 4], [7,84], [31,Chap. 5.4], [66]. Typi- (c) Alternation of linear angular and linear spatial diffusion visually destroys the anisotropic fiber structure. In contrast the hypo-elliptic diffusion in (d) Preserves the fiber-structure much better cally, such a voting field is based on the co-circularity principle [57,58]. This principle is naturally included in the leftinvariant PDE's on R 3 S 2 , as propagation locally takes place along horizontal exponential curves whose spatial projections have constant curvature, cf. [36]. To get a quick intuition of our hypo-elliptic diffusions on R 3 S 2 , one can regard such a Green's function as a "voting field" of glyphs (see Fig. 3). The shift-twist convolution aligns this voting field of glyphs with each position and orientation in the DW-MRI data. A shift-twist convolution is the example of a linear, left-invariant operator on L 2 (R 3 S 2 ) (the space of DW-MRI data), cf. [36, Lemma 1, Corollary 1]. A left-invariant operator on L 2 (R 3 S 2 ) is an operator that commutes with rotations and translations.

Why Morphological Scale Spaces on R 3 S 2 ?
Typically, if linear diffusions are directly applied to DTIdata, glyphs are propagated in too many directions. There- Fig. 3 The Green's functions for line propagation (completion and enhancement) in R 2 S 1 ≡ SE (2), left, (in dashed lines their Heisenberg approximations [32,35]) and in R 3 S 2 := SE(3)/({0}×SO(2)) depicted on the right using glyph-visualizations according to Definition 1. They provide a probabilistic "voting field" of glyphs fore we combined these diffusions with monotonic transformations in the codomain R + , such as squaring input and output cf. [36,65,68]. Visually, this produces anatomically plausible results (Figs. 5 and 6), but this fails if there are large global variations present in the data. This is often problematic around ventricle areas in the brain, where the glyphs are usually larger than those along the fibers, see Fig. 7(toprow). In order to achieve a better way of sharpening the data where global maxima do not dominate the sharpening of the data, cf. Fig. 4, we propose morphological scale spaces (erosions) on R 3 S 2 where transport takes place orthogonal to the fibers. The result of such an erosion after application of a linear diffusion is depicted down left in Fig. 7, where the diffusion has created crossings and where the erosion has visually sharpened the fibers.

Objectives
We provide a list of 9 issues that arise from previous work [36,65,68], organized in 4 categories: • Morphological scale space theory on R 3 S 2 .
1. Can we generalize the left-invariant convolutions on R 3 S 2 to left-invariant morphological convolutions?
2. Can we replace the grey-scale transformations [65,36,68] by Hamilton-Jacobi-Bellman (HJB) equations (erosions) on R 3 S 2 . Can we derive a morphological scale space theory on R 3 S 2 that allows us to visually sharpen the fibers in the DW-MRI data in a geometrical way? 3. Can we find the unique viscosity solutions of these HJB-equations? 4. Can we find analytic approximations for the viscosity solutions of these left-invariant HJB-equations on R 3 S 2 , akin to the analytic approximations of linear left-invariant diffusions, cf. [36, Chap. 6.2]?
5. Can we construct basic numerical finite difference schemes and analytic convolution schemes for morphological scale spaces on R 3 S 2 ?
• Combining morphology and diffusion.
6. Can we combine left-invariant diffusions and leftinvariant dilations in a pseudo-linear scale space [40] on R 3 S 2 ?
7. Can we combine diffusions and erosions for fiberenhancement?
• Probabilistic models of scale spaces on R 3 S 2 .
8. Can we provide a probabilistic interpretation of both linear and morphological scale spaces on R 3 S 2 ? 9. The contour completion kernel on R 3 S 2 suffers from a severe singularity at the origin. Can we get contour completion to work via iteration of resolvent processes? Does this relate to a new probabilistic model for contour completion?
The main objectives on which this article concentrates are issues 2, 3, 5 and 9. This article addresses all nine questions with affirmative answers and formal theorems. The approach dealing with the first issue fits within general group morphology [69], where we tackle some additional issues related to the group quotient structure of R 3 S 2 . The approach dealing with the sixth issue is placed in Appendix C as it is less effective than the approach dealing with the seventh issue that is applied in our final experiments. To achieve our objectives, we introduce linear, morphological and pseudo-linear scale spaces defined on where the input DW-MRI image serves as initial condition W (y, n, 0) = U(y, n). For a formal definitions of these scale spaces, see Appendix D.
In nearly all experiments we considered DTI-data as input. In this case the initial condition is either given by Eq. (3) or by for all y ∈ R 3 , n ∈ S 2 as proposed in the previous works [36,43,65,68]. In practice option (4) tends to make glyphs more isotropic. Furthermore, we introduce a theoretical Lagrangian and Hamiltonian framework on R 3 S 2 . Although not considered here, this framework could also be used for extension of geometric control and fiber-tracking [53] to geometrical control problems on sub-Riemannian manifolds [4] within the joint domain of positions and orientations, relating to Finsler geometry [3]. Moreover, our framework may also be employed in the modeling from the physically scanned data and the DTI and HARDI modeling via inverse problems on R 3 S 2 . Figure 7 shows a preview of how our scale space evolutions perform on the same neural DTI dataset considered in [65], where we used for the sake of comparison, as this was also applied in the previous works [65,68]. This min-max normalization induces a field of glyphs with similar sizes. In the final DTIexperiments of this article we will not apply this operator: in Figs. 21, 23 we directly visualize the concatenation of diffusion, erosion and minimum subtraction. This article primarily focusses on a novel mathematical theory for enhancement of probability distributions on R 3 × S 2 . Our experiments only serve as first feasibility studies supporting its potential value in DW-MRI applications. For validations that our enhancements improve subsequent fiber-tracking in a relevant neuro-imaging application (detection of the optic radiation for planning temporal lobe resection in the treatment of epilepsy) we refer to related works [27,77].

Structure of the Article
In Sect. 2 we motivate that the domain of DW-MRI is not to be considered as a flat Cartesian product of R 3 and S 2 . Instead we consider the coupled space R 3 S 2 of positions and orientations as a group quotient in the 3D-Euclidean motion group SE(3).
In Sect. 3 we show that operators on DW-MRI images that commute with rotations and translations must be leftinvariant. Furthermore, in Theorem 1 we classify all leftinvariant operators on L 2 (SE (3)) that allow well-defined restrictions to L 2 (R 3 S 2 ) (the space of DW-MRI data). Such operators will be called "legal" operators.
In Sect. 4 we consider both linear and morphological convolutions on DW-MRI. According to Theorem 2 in DTI and HARDI data containing fibers of the corpus callosum and the corona radiata in a human brain, with b-value 1000 s/mm 2 on voxels of (2mm) 3 , cf. [65]. We visualize a 10 × 16-slice of interest (162 samples on S 2 using icosahedron tessellations) from 104 × 104 × 10 × (162 × 3) datasets. Top row: region of interest with fractional anisotropy intensities with colorcoded DTI-principal direc-tions. Middle row, DTI data U visualized according to Eq. (2) resp. Definition 1. Bottom row: HARDI data (Q-ball with l ≤ 4, [29]) of the same region, hypo-elliptically diffused DTI data (y, n) → W (y, n, t), Eq. (53), using Eq. (4) as input. We applied min-max normalization of W (y, ·, t) for all positions y Sect. 4.1 all (reasonable) linear left-invariant operators are R 3 S 2 -convolutions. We show in Theorem 2 that R 3 S 2convolutions can be implemented without a goniometric parametrization of S 2 . Subsequently, we extend linear convolutions on R 3 S 2 to morphological convolutions and in Theorem 3 we show that their implementation is similar.
In this article we aim for both linear and morphological scale space PDE's on R 3 S 2 . This requires an introduction to the Lie algebra of left-invariant vector fields which serves as a local frame of reference attached to fiber-fragments in the DW-MRI data. This is done in Sect. 5, where in Theorem 4 it follows that the choice of legal scale space generators is quite limited.
Subsequently, in Sect. 6 we introduce two different sub-Riemannian manifolds within SE(3). One of them underlies hypo-elliptic diffusions and one of them underlies the erosions. The practical motivation for this, is that we want diffusion along fiber fragments, whereas we want to erode orthogonal to them. Figure 10 provides the crucial geometric intuition. Theorems 5 and 7 relate the left-invariant vector fields to the Frenet frame of the spatial part of the curves in the two sub-Riemannian manifolds. We introduce metrics and geodesics on both sub-Riemannian manifolds and we relate them to a basic curve optimization problem in R 3 , where we pay attention to the local existence of minimizing geodesics. Furthermore, in Theorem 8 we relate the Lagrangians and Hamiltonians on the sub-Riemannian manifold(s) via a Fenchel transform.
In Sect. 7 we consider both linear and morphological scale spaces where we employ the geometric tools of the  53), is applied to the HARDI dataset previous section. In Sect. 7.1 we classify all legal linear convection-diffusion. Then in Sect. 7.2 we introduce morphological scale spaces via a Hamilton-Jacobi-Bellman equation induced by the Hamiltonian of the previous section. Subsequently, in Theorem 9 we show that iso-contours of morphological scale spaces are geodesically equidistant and in Theorem 10 we show that their viscosity solutions are given by convolution with the corresponding morphological Green's functions. Additionally, in Appendix C we consider pseudo-linear scale spaces on R 3 S 2 which form a continuous transition between the linear and the morphological scale spaces on R 3 S 2 .
In Sect. 8 we apply the Ball-Box theorem and we obtain tangible estimates for both the linear and morphological Green's functions.
In Sect. 9 we discuss two different implementations of morphological scale spaces: 1 via morphological convolution and 2 via finite differences. The first approach directly relates to the unique viscosity solutions of the Hamilton-Jacobi-Bellmann equations (Theorem 10) and allows fast parallel algorithms via look-up tables, 1 whereas the second approach allows data-adaptation (Sect. 9.5) via fast upwind schemes using small finite difference stencils.
In Sect. 10 we provide a concise overview of the underlying probability theory of both linear and morphological scale spaces on R 3 S 2 . The linear scale spaces on 1 Akin to accelerated implementations [68] of linear shift-twist convolutions on R 3 S 2 R 3 S 2 are Fokker-Plank equations of stochastic processes for fiber completion/enhancement, recall Fig. 3. In Theorem 12 we present a practical and natural extension of the 3D contour-completion process. The morphological scale spaces are Hamilton-Jacobi-Bellmann equations describe cost-processes on R 3 S 2 .
Finally, in Sect. 11 we provide experiments where diffusions and erosions are combined.

The Space R 3 S and Its Embedding in the Rigid Body Motion Group SE(3)
For the enhancement of data, it is crucial to determine which fiber fragments are well aligned with each other. To find out which fiber fragments are well aligned, you must consider the spatial and orientational information together, as they interact with each other. Figure 8 shows how (x 0 , n 0 ) and (x 1 , n 1 ) are better aligned (by forming a more natural connecting curve) than (x 0 , n 0 ) and (x 2 , n 1 ) are, even though both pairs have the same spatial distance from each other, and the same difference in angles. This shows us that the space of positions and orientations is coupled and it is therefor conceptually wrong to consider the space of positions and orientations as a flat Euclidean space. Suppose we have a fiber fragment at position y ∈ R 3 with orientation n ∈ S 2 . This fiber fragment can be rotated and translated via the action of SE(3) (rigid body motion group) on R 3 × S 2 :  [25,26]. Evolutions are implemented by finite difference schemes [26,36], with step size t (Color figure online) for all x, y ∈ R 3 , R ∈ SO(3) and all n ∈ S 2 . Subsequent alignment reveals the non-commutative product of SE (3): where we see that SE(3) = R 3 SO(3) is the semi-direct product of SO(3) and R 3 as it involves a homomorphism This homomorphism is responsible for the non-commutative structure of the group SE(3). Considered as a set we have SE(3) = R 3 × SO(3), but considered as a group we must write R 3 SO(3) := R 3 τ SO(3) as its non-commutative group product, Eq. (7), is not equal to the commutative group product ( The group SE(3) acts transitively on R 3 × S 2 as every position and orientation can be reached from the unit position orientation (0, e z ), with e z = (0, 0, 1) T : where R n denotes any rotation such that R n e z = n.
To keep track of rigid body motions acting on R 3 S 2 it will be useful to embed the space of positions and orientations in SE (3). However one should be careful that identifying n ↔ R n via (8) is not unique. Indeed if we multiply R n with a rotation around the z axis via an angle α from the right, the identification Eq. (8) is still valid. To get rid of this problem The spatial and angular distance between (x 1 , n 1 ) and (x 0 , n 0 ) is the same as the spatial and angular distance of (x 2 , n 1 ) between (x 0 , n 0 ). However, (x 1 , n 1 ) is much more aligned with (x 0 , n 0 ) than (x 2 , n 1 ) is. The space R 3 S 2 takes this alignment into account (in contrast to R 3 × S 2 ). The connecting curves are the spatial projections of the geodesics in R 3 S 2 for β = 1 2 in Eq. (43) (with x 0 = (0, 0, 0), x 1 = (0, 0, 5) and x 2 = (3, 0, 4), n 0 = (0, 0, 1) and n 1 = (sin 2π 5 , 0, cos 2π 5 )) we consider the coupled space of positions and orientations as the group quotient 2 where we identify SO(2) with all rotations around the z-axis. Elements in R 3 S 2 are equivalence classes of rigid body motions under the equivalence relation where from now on R a,ψ denotes the counter-clockwise rotation around axis a ∈ S 2 by angle ψ ∈ [0, 2π). So each equivalence class [(y, To keep a sober notation we denote these equivalence classes with (y, n) ∈ R 3 S 2 where we use Eq. (8) to uniquely identify (y, n) ↔ (y, R n ) . 2 This involves a slight abuse of notation as S 2 = SO(3)/SO(2) is not a group (due to the hairy-ball/Brouwer's Theorem) and hence R 3 S 2 is not formally a semi-direct product of groups.
Here we would like to keep track of functions on the equivalence classes U : R 3 S 2 → R given by (y, n) → U(y, n) and the corresponding functions on the group Note that U(y, n) =Ũ(y, R n ) for all y ∈ R 3 and n ∈ S 2 . Instead of applying a rigid body motion to a single fiber fragment we can also apply a rigid body motion on a whole distribution on positions and orientations (in particular a DW-MRI image). This is done via with L g 1 g 2 = L g 1 L g 2 for any pair g 1 , g 2 ∈ SE(3).

Legal Operators
In the previous section we introduced a group structure on the domain of DW-MRI, where the rigid body motion group SE(3) acts transitively on the space of positions and orientations. In this section we will classify the operators U →Φ(Ũ) that uniquely relate to operators U → Φ(U ) defined on the group quotient R 3 S 2 , such that Φ commutes with rotations and translations. Later on Φ, will be an enhancement operator on DW-MRI, in which case the commutation requirement is natural. Mathematically, this is achieved by imposing left-invariance, i.e.
on an operator Φ : Occasionally, we shall resort in our design to the space of square integrable functions on the group. Here we need to consider the left and right regular action of SE(3) onto L 2 (SE(3)): for all V ∈ L 2 (SE (3)). If V corresponds to U : R 3 S 2 → R + , i.e. V =Ũ (Eq. (11)) then we have R h V = V for all h = (0, R e z ,α ) and all α ∈ [0, 2π). Therefore we need the following definition.

Definition 2
We define the spaceL 2 (SE(3)) of functions on the group that uniquely relate to functions on the quotient viã The next theorem specifies necessary and sufficient condition for operators onŨ : SE(3) → R + to be legal in the sense that they correspond to well-defined operators on U : R 3 S 2 → R + that commute with rotations and translations.
Theorem 1 Consider an operatorΦ on square integrable functions defined on the groupΦ :L 2 (SE(3)) → L 2 (SE (3)). Then such an operator uniquely relates to a left-invariant operator Φ : if and only if with tilde-operator U →Ũ given by Eq. (11).
Proof The natural identification between Φ andΦ in Eq. (14) is unique if and only if the choice of R n (which by definition is any rotation that maps e z to n) in the second equality in Eq. (14) does not matter and that is the case if and only ifΦ = R h •Φ for all h ∈ {0} × SO(2). With respect to the left-invariance we note that left action of SE(3) on the left cosets coincides with left action on the group: (2) and consequently we haveΦ

Operator
is legal and coincides with the identity operator. 2. Operator L g given by Eq. (91) is illegal as L is a group representation and SE(3) is not commutative. 3. One can construct legal operators from illegal operators by addition and concatenation. This will occur frequently in this article. The common way to parameterize SO(3) is via the Euler angles R(γ , β, α) = R e z ,γ R e y ,β R e z ,α , with α ∈ [0, 2π), β ∈ [0, π) and γ ∈ [0, 2π). This yields the following parametrization of S 2 n(β, γ ) = R(γ , β, α)e z = (cos γ sin β, sin γ sin β, cos β) T Like all parameterizations of S 2 = SO(3)/SO(2), the Euler angle parametrization suffers from the problem that there does not exist a global diffeomorphism from a sphere to a plane. We need a second chart to cover SO(3); which again parameterizes S 2 via different ball-coordinates β ∈ [−π, π),γ ∈ (− π 2 , π 2 ), n(β,γ ) = R e x ,γ R e y ,β R e z ,α e z = (sinβ, − cosβ sinγ , cosβ cosγ ) T , (17) but which has ambiguities at the intersection of the equator with the x-axis, cf. [36] and Fig. 9. We will use the second chart in our analysis as it does not have a singularity at the unity element. Operators on functions on the group SE(3), expressed in these Euler angle parameterizations on the group SE(3), are legal if they are independent of the final angle (respectively α andα), see Theorem 1. This holds in particular for scale space operators on SE(3) that we will consider later in Sect. 7.

Convolutions on R 3 S 2
In the previous section we classified the legal operatorsΦ acting on functions on the group. Such operators are leftinvariant. If in addition such operator is linear it is a convolution on the group [21,43]. The quotient structure requires some additional symmetry constraints on the convolution kernel. To avoid these technicalities we directly consider the corresponding operator Φ : L 2 (R 3 S 2 ) → L 2 (R 3 S 2 ) in the cases where Φ is a linear, respectively morphological, convolution on R 3 S 2 .

Linear Convolutions
According to the theorem below, all reasonable linear, left-invariant operators on DW-MRI images are R 3 S 2convolutions.
. Then there exists an integrable kernel k y, n; y , n 2 dy dσ n , is finite and we have (KU)(y, n) = R 3 S 2 k y, n; y , n U y , n dy dσ n , for almost every (y, n) ∈ R 3 S 2 and all U ∈ L 2 (R 3 S 2 ).
Then to each positive left-invariant kernel with R 3 S 2 k(0, e z ; y, n)dσ (n)dy = 1 we associate a unique probability density p : R 3 S 2 → R + with the invariance property by means of p(y, n) = k(y, n; 0, e z ). The convolution now reads where σ is the surface measure on S 2 and where R n ∈ SO(3) s.t. n = R n e z .
For a proof we refer to [36,Chap. 3]. For parallel implementation of R 3 S 2 -convolution (where we pre-compute all rotated and translated reflected kernels) we rely on the following Lemma.
Then for all y ∈ R 3 and all n ∈ S 2 we have where (·, ·) L 2 (R 3 S 2 ) denotes the L 2 -inner product and with the reflected kernelp : R 3 S 2 → R + given by for all y ∈ R 3 and all n ∈ S 2 .
For more details on efficient implementation of such a convolution see [36,68]. There are two relevant issues regarding the implementation we did not explicitly address in our previous works and which are also relevant in implementing the morphological convolutions that we shall introduce in the next subsection.
Lemma 2 Implementation of a R 3 S 2 -convolution via the pre-computed kernels in Lemma 1 can be done without goniometric formulas, that is R T n v , with v = y − y, in Lemma 1 (and Theorem 3) can be computed via: Proof Because of Eq. (19) we are free in the choice of R n ∈ SO(3) such that Eq. (9) holds. We take the counterclockwise rotation in the plane spanned by e z and n such that e z is mapped on to n. To this end we use the Rodrigues' formula where we take a = e z × n and cos ψ = n · e z = n 3 and sin ψ = 1 − |n 3 | 2 from which the result follows after taking the transpose of the corresponding matrix (expressed in the standard basis).
Remark For the pre-computation of the analytical approximations of the Green's functions [36,68] we do have to invert Eq. (17), which yields: with sign(x) = 1 if x > 0 and zero else, and with sign(x) = 1 if x ≥ 0 and zero else.

Morphological Convolutions
Dilations on the joint space of positions and orientations R 3 S 2 are obtained by replacing the (+, ·)-algebra by the where k : R 3 S 2 → R − denotes a morphological kernel. If this morphological kernel is induced by a semigroup (or evolution) then we write k t for the kernel at time t. Erosions on R 3 S 2 are given by: Dilation kernels are negative and erosion kernels are positive 3 and therefore we write for the corresponding erosion kernels. Implementation of morphological convolutions (i.e. dilations/erosions) is entirely similar to linear convolution implementation. Indeed in analogy with Lemma 1 we have the following result.
Then for all y ∈ R 3 and all n ∈ S 2 we have with reflection operator k →ǩ as in Eq. (21).
A tangential result can be obtained for dilations. A particular example of a kernel k is the morphological delta distribution which has the following reproducing properties for all bounded U : R 3 S 2 → R + .

A Moving Frame of Reference and Metric Tensor for Scale Spaces on R 3 S 2
In this article we aim for particular cases for operators Φ andΦ in Sect. 3. Namely, we aim for evolutions Φ = Φ t with fixed time t > 0 satisfying the semi-group property where we want to generalize both morphological scale spaces [17,40,82,83,86] and linear scale spaces [8, 30, 47, 48, 50-52, 54, 55, 85] on R 3 to R 3 S 2 . Generally speaking, this mission is achieved by replacing the bi-invariant generators ∂ x , ∂ y , ∂ z on the commutative group R 3 by leftinvariant generators on the Euclidean motion group SE(3). In this section we will derive these left-invariant generators and consider their important geometric interpretation. Evolutions on DW-MRI images must commute with rotations and translations. Therefore our evolutions on DW-MRI images (and the underlying metric-tensor) are expressed in a local frame of reference attached to fiber fragments. This frame of reference {A 1 , . . . , A 6 } consists of 6 left-invariant vector fields on SE(3) given by where {A 1 , . . . , A 6 } is the basis for the Lie-algebra at the unity element and (3) is the exponential map in SE (3). (3), using the pushforward (L g ) * of the leftmultiplication which is defined by, i.e. .
where e = (0, I ) denotes the unity element in SE (3), which is the case iff with A i given by Eq. (27) and constant c 1 , . . . , c 6 .
. Therefore, we obtain a local basis of left-invariant differential operators via the infinitesimal generators associated to R, that is via Equation (27) implies imply that condition (28) is necessary and sufficient. Let These left-invariant vector fields can be expressed in our second chart, Eq. (16), as forβ = π 2 andβ = − π 2 . However, these technical formulas are only needed for analytic approximation of Green's functions, see [36,Chap. 6]. In practice one uses finite difference approximations [36,Chap. 7], where spherical interpolation in between higher order tessellation of the icosahedron can be done by means of the discrete spherical harmonic transform [36,Chap. 7.1] or preferably by triangular interpolation (for full details see [26]).
The left-invariant vector fields span the Lie-algebra L(SE (3)) of the group and we have the following the fol- with structure constants So far we considered operators on SE(3). However, our MRI-datasets are defined on the quotient R 3 S 2 given by Eq. (10). So let us figure out what operators give rise to welldefined operators on the quotient, i.e. what operators actually make sense on DW-MRI data. Recall the definition of legal operators on the group, Definition 3. Recall also that legal operators are precisely those operators on the group that allow extension to the quotient. Intuitively, one has to ensure that the operators are isotropic in the grey-planes in Fig. 10.

Theorem 4
We have the following list of legal operators: 1. The only left-invariant vector field that is legal are A 3 and A 6 . 2. The only left-invariant 2nd order linear differential operators that are legal are given bỹ

The only left-invariant metric tensor
Proof Linearity and left-invariance immediately implies that a and D should be constant. By Theorem 1 we have one more condition remaining, for all g ∈ SE(3), with Zα = R e z ,α ⊕ R e z ,α ∈ SO(6). So application of Schur's Lemma [73] applied to the standard irreducible matrix representation of SO(3) into the space of linear automorphisms on R 3 , implies that if a matrix commutes with all rotations it is a multiple of the identity. This yields the only option Eq. (34). W.r.t. the metric tensor we have two sufficient and necessary constraints for all g, q, h ∈ SE(3) and all vector fieldsX,Ỹ ∈ T (SE (3)).
The first requirement (left-invariance) requires the metric tensor components to be constant, the second requires the matrix [g ij ] to commute with Zα from which the result follows.
Remark We will restrict ourselves to the spaceL 2 (SE(3)), Eq. (13). On this space A 6 vanishes so from now on we exclude all operators involving A 6 or dA 6 . Furthermore, restriction to the spaceL 2 , where S 2 denotes the Laplace-Beltrami operator on the sphere.
The well-defined metric tensor on the quotient is thereby parameterized by the values of D 11 , D 33 and D 44 , and we write the metric as a tensor product of left-invariant covectors: The metric tensor on the quotient R 3 S 2 is thereby given by where vector fields are described by the differential opera- for j = 1, 2, 3, where R e j ,h denotes the counter-clockwise rotation around axis e j by angle h, with e 1 = (1, 0, 0) T , e 2 = (0, 1, 0) T , e 3 = (0, 0, 1) T .

Sub-Riemannian Manifolds within SE(3) and Sub-Riemannian Metrics on R 3 S 2
In the previous section we have classified all legal metric tensors on SE(3) and we have parameterized all legal Riemannian spaces (SE (3), G) by the triplet (D 11 , D 33 , D 44 ) ∈ (R + ) 3 . Depending on the choices for D ii , the metric tensor punishes certain tangent directions more than others. Often one would like to impose that flows in our evolutions along specific directions within the tangent space are prohibited (i.e. have infinite cost). See for example Fig. 10 To ensure that certain tangent vector directions are prohibited within the flow of the evolutions U → Φ t (U ) we need the concept of sub-Riemannian manifolds. This section studies the differential geometrical consequences of such a prohibition. Together with the preceding group theoretical chapters, it forms the mathematical foundation in the design and analysis of our linear and morphological scale space operators U → Φ t (U ) which we will introduce in Sect. 7, analyze in Sect. 8, implement in Sect. 9, and apply in Sect. 11. When accepting the formula's for  (27). To be well-defined on R 3 S 2 we must impose isotropy in the tangent planes span{A 1 , A 2 } and span{A 4 , A 5 }. Diffusion/convection primarily takes place along A 3 in space and (outward) in the plane span{A 4 , A 5 } tangent to S 2 . Erosion takes place both inward in the tangent plane span{A 1 , A 2 } in space and inward in the plane span{A 4 , A 5 } the Hamiltonian and Lagrangian derived in this section the reader may choose to skip the mathematical details in this section and continue with the next section. In this case, Fig. 10 provides a quick summary of this chapter.
with the extra constraint that certain subspaces of the tangent space are prohibited, i.e. for all curves γ in (M, θ 1 , . . . , θ n ) we have The above definition and notation is a suitable simplification of a more accurate 4 notation and definition of a sub-Riemannian manifold. Curves in M satisfying (38) are called horizontal curves.
Next we consider two specific sub-Riemannian manifolds of (SE (3),G), that we need when defining respectively linear scale spaces and morphological scale spaces. In both examples we ensure that the available tangent vectors together with their commutators fill the full tangent space T (SE (3)), so that every pair of points in SE(3) can be connected by a horizontal curve.
for all s ∈ [0, L]. Curves satisfying Eq. (39) are called horizontal curves in SE (3) and their tangent vectors are given bẏ where we use spatial arc length parametrization and short notatioṅ 4 In geometric control theory is the curvature of the spatial part of the curve. The normal N and binormal B of the spatial part s → x(s) of a horizontal curve are given by Proof Recall from Eq. (30) that A 3 is aligned withñ(β,γ ), Eq. (17). Furthermore for spatial arc length parametrization one has ẋ(s) =γ 3 (s) = 1, since in the sub-Riemannian manifold we haveγ 1 =γ 2 = 0. W.r.t. the remainder we note that 6 where c k ij denote the structure constant of the Lie-algebra. Differentiation ofẋ(s) = A 3 (where one applies the product rule and Eq. (40)) produces the required results.
In previous work [36] we have shown that linear leftinvariant diffusion in the group takes place along exponential curves. 7 When restricting to horizontal diffusion 5 We take a different convention here than in our technical report [37, [36,Chap. 5]. It corresponds to the usual differentiation of moving frames in rigid body mechanics. 7 The exponential curves are auto-parallels of the Cartan connection. This connection has torsion and as a result the minimizing curves (geodesics) of (42) do not coincide with the auto-parallels.
The only legal metric tensors on this sub-Riemannian manifold are given by The corresponding metric distance on this sub-Riemannian manifold is given by with s the spatial arc length parameter and with

Theorem 6
The problem (42) is well-defined for g −1 2 g 1 close enough 8 to e. The equivalent problem on R 3 S 2 (or rather R 3 ) is: where s, L > 0, and κ(s) are respectively spatial arc length, total length, and curvature of the spatial part of the curve.
Proof The metric tensor is left-invariant so that d( is within the range of the exponential map of the corresponding geometric control problem with particles in SE(3) that move forwardly (i.e.ẋ(s) = +A 3 ) along horizontal curves one does not cross cusps [14], similar to the (SE (2), − sin θ dx + cos θ dy)case [15]. In such a case it follows that every stationary If it is cheaper to set the car in reverse to reach the final destination a smooth minimizer does not exist [14]. Bottom right: In black the planar boundary conditions for which the sub-Riemannian geodesics do not suffer from cusps (for β = 1) Eq. (43) has a unique smooth minimizer), for details see [15] curve (satisfying the Euler-Lagrange equations on the sub-Riemannian manifold) is a global optimizer. These Euler-Lagrange equations can be found by the Pontryagin maximum principle [4] or by Marsden-Weinstein reduction [16] as we did in [37,Appendix G]. For a visualization and characterization of this exponential map in respectively SE (2) and SE(3) and further details see [15,44]. The equivalence now directly follows from the identities in Theorem 5 (where the distance is independent of the choice of R n satisfying Eq. (9)).
These sub-Riemannian geodesics have as a natural counter-part the elastica curves, where in stead of κ(s) 2 + β 2 ds one minimizes the elastica functional κ(s) 2 + β 2 ds, cf. [38,59,72]. The advantage of elastica is that they do not involve cusps, so one does not have to constrain the end-points. Their disadvantage is the vast involvement of special functions and their loss of both global and local optimality, which is fully analyzed for the 2D-case in [72]. See Fig. 11 for an example of geodesic and an elastica and an illustration of the cusp-problem as reported by Boscain and Rossi [14]. For more details see [15,37,44].
If initial point (x 0 , n 0 ) = (x(0), n(0)) and endpoint (x 1 , n 1 ) = (x(L), n(L)) are chosen such that (43) is equivalent to the sub-Riemannian distance on (SE(2), − sinθ dx + cosθdỹ). Right: full field of reachable angles. Bottom left, the set of end-points in SE (2) where a global minimizer exists is an unbounded orientable volume with the cusp-surface as boundary the sub-Riemannian distance on R 3 S 2 given by Eq. (43) is equivalent to the sub-Riemannian distance on (SE (2), − sinθ dx + cosθ dỹ) that has been studied in full detail by Sachkov and Moiseev [70,71]. See Fig. 12(top left) for the convention of the (x,ỹ)-coordinate system.
In Fig. 8 we have depicted such minimizing curves and in Fig. 12 we have plotted the cones of reachable angles (which is precisely the set of end-points in SE(2) ≡ R 2 S 1 where the sub-Riemannian problem has a (global) minimizer). In orientation space this reachable area is a non-compact volume bounded by the cusp-surface (where geodesics arise with infinite curvature, which correspond to the boundaries of the hyperbolic phase portrait in [33, Appendix A]). For more details and formal proofs see [14,15,44].
for all s ∈ [0, L]. These horizontal curves satisfẏ x(s)⊥n(s) ∈ S 2 where the spherical part n(s) is not to be mistaken with the normal N(s) to the spatial part of the curve. More precisely, Theorem 7 Along the spatial part of a horizontal curveγ in (SE(3), dA 3 , dA 6 ) we have the following Frenet frame with curvature magnitude where we differentiate w.r.t. to spatial arc length.
Remark In our sub-Riemannian metric definitions (Eq. (46) and Eq. (42)) we always impose dA 6 | γ ,γ = 0 to avoid 9 artificial curvature and torsion in the spatial part of the geodesics, for details see cf.

Hamiltonians on (SE(3), dA 3 , dA 6 )
Recall, that in the sub-Riemannian manifold (SE(3), dA 1 , dA 2 , dA 6 ) we considered elastica besides the geodesics. We obtained the elastica functional by squaring the integrand I(s) while parameterizing by spatial arc length. In (SE(3), dA 3 , dA 6 ) we do a similar thing, that is we consider powers of the integrands while parameterizing by spatial arc length s ∈ [0, L]: Note that the metric-integrand relates to the Lagrangian via The Lagrangian can also be parameterized by the arc length parameter, say p, in the sub-Riemannian manifold, since by Eq. (48) we have with L the spatial length and t the sub-Riemannian length of the curve where we note that Now that we have introduced a Lagrangian on (SE(3), dA 1 , dA 2 , dA 6 ) we construct the Hamiltonian by means of the Fenchel transform on the Lie-algebra of left-invariant vector fields L(SE (3)). This construction is a generalization of the relation between Lagrangian and Hamiltonian on R n , see [39,Chap. 3.2.2]. In the subsequent section, the Hamiltonian will serve in the generator of morphological scale spaces, where we will generalize morphological scale spaces on R n to morphological scale spaces on R 3 S 2 .
Definition 6 Let X be a normed space with dual space X * , and let L : X → R ∪ {∞} be a real-valued function on X, then the Legendre-Fenchel transform L → FL on X is given by where x, y = x(y) and x ∈ X * .

Theorem 8 The Hamiltonian associated to the Lagrangian (47) is given by
Proof The first identity follows by application of the Pontryagin Maximum Principle [4,64] to the sub-Riemannian geodesic problem given in Eq. (46) on (SE (3), dA 1 , dA 2 , dA 6 ). The second identity follows from Eq. (50).

Left-Invariant Scale Spaces on R 3 S 2
In this section we consider two special cases of evolutions (enhancement) on DW-MRI. In both cases (y, n, t) → (Φ t (U ))(y, n) is the solution of a PDE system (evolution) defined on (R 3 S 2 ) × R + with initial condition Φ 0 (U ) = U . For smoothing and fiber propagation in the DW-MRI data U : R 3 S 2 → R + we consider linear scale spaces on R 3 S 2 , whereas for well-posed sharpening we consider morphological scale spaces (erosions) on R 3 S 2 .
For the sake of sober notation we write W (y, n, t) := (Φ t (U ))(y, n) for all t ≥ 0.

Linear Scale Spaces on R 3 S 2
The spherical and the spatial Laplacian can be expressed in terms of the left-invariant vector fields: where we restrict the rotation generators to the spacẽ L 2 (SE(3)) ⊂ L 2 (SE(3)), Eq. (13). These Laplacians generate diffusion over S 2 and R 3 separately and are thereby likely to destroy the fiber structure in DW-MRI, recall Fig. 7. Recall our result on legal operators, Theorem 4 and recall Definition 2. From Theorem 4, Eq. (52) and A 6 |L 2 (SE(3)) = 0, we deduce the following result:

Corollary 2
The only legal linear left-invariant convectiondiffusion equations on R 3 S 2 are the following scale spaces: 1. Angular diffusion generated by S 2 , 2. Spatial diffusion generated by R 3 , 3. Contour completion generated by ±A 3 + D 44 S 2 , 4. Hypo-elliptic contour enhancement generated by In the final option we set 0 < D 11 D 33 , as otherwise fibers propagate too much to the side. To this end we note that for D 11 ≥ D 33 we have The scale space evolutions for hypo-elliptic linear contour enhancement are given by: with D 11 = D 22 , D 44 = D 55 denotes the generator of the linear scale spaces on R 3 S 2 , with differential operators A i on R 3 × S 2 given by Eq. (37). Next we address the probabilistic interpretation of the left-invariant (legal) linear scale spaces. The density (y, n) → W (y, n, t) represents the probability density of finding a random orientated particle (Y t , N t ) at position y ∈ R 3 with orientation n at time t > 0 given that the initial distribution of random particles is given by (y, n) → U(y, n). Now in a Markov-process traveling time is memoryless and the only continuous memoryless distribution is the negative exponential distribution P(T = t) = λe −λt . When we condition out the traveling time we obtain the resolvent processes This resolvent process coincides with Tikhonov-regularization [36, Theorem 2]. It also appears in implicit finite difference implementations of the linear evolutions [26], when setting λ −1 = t. For contour completion the resolvent process is crucial since it integrates the propagation fronts to a smooth integrable density with a singularity at the origin [36]. However, in contrast to the time-dependent process, it does not satisfy the semigroup property (a crucial scale space axiom [30,85]). We have p t * R 3 S 2 p s = p s+t (57) which has the following implications on the underlying independent random variables: for all t, s > 0. However, for resolvents one has Thereby, when alternating the resolvent direction processes (e.g. in an implicit finite difference scheme of the time evolutions), the following relevant questions rise: • What happens from a probabilistic point of view when alternating these resolvent direction processes? • What happens with the effective probability density and the singularity at the origin? • Can we still relate the iterated resolvent process to the original process?
We will answer these questions in the probability section, Sect. 10.2.
In our previous work [36] one can find suitable approximations for all the Green's functions. See Fig. 3 for an overview with comparison to 2D-equivalents of the Green's functions. Reasonable parameter settings for 2Dcontour completion are (λ, D 11 1.25, 1, 0.04). The analogy in parameter-settings comes from the fact that 3Dcompletion and enhancement kernels can be approximated by direct products of their 2D-counterparts [36]. Increase t enlarges the time-dependent kernel. Increase of λ −1 enlarges the resolvent kernel. By a simple re-scaling argument we can set D 33 = 1 in which case the parameter D 44 is equal to β in the underlying metric-tensor Eq. (43) and thereby tunes the bending of the fiber-propagation kernel. In [5,17] the authors put an interesting connection between morphological and linear scale spaces on R n via the Cramer transform. The Cramer transform does not generalize easily to SE(3), but nevertheless the analogy and probabilistic interpretations of both type of scale spaces remain. = ±H η dW (y, R n , t) where the left-invariant gradient of the corresponding func-tionW (·, t) on SE(3) (recall Eq. (11)) is given by Remark For η → ∞ we naturally arrive at the timeindependent HJB-equation (Eikonal Eq.) where the Hamiltonian H ∞ relates to the homogeneous Lagrangian One solution is W (y, n) = δ((0, I ), (y, R n )) where δ denotes the sub-Riemannian distance on (SE(3), dA 3 , dA 6 ), cf. Eq. (46).
The HJB-equations on R 3 S 2 have an important geometric property (similar to HJB-equations on R n , cf. [60,74,81]) as we will show next.

Definition 7 A family {S t } t∈R + of surfaces
is geodesically equidistant if for all geodesics passing these surfaces transversally with The equivalence in Eq. (62) denotes equality up to a monotonic re-parametrization of time. The morphological scale spaces (59) do not have unique solutions. This is well-known for the special case where one has only angular erosion or only spatial erosion [24,39]. Based on these works we have the following definition of viscosity solutions on R 3 S 2 .

Definition 8 Suppose that the Hamiltonian is convex and
Then a viscosity solution is a bounded and continuous (not necessarily differentiable) weak solution W : (R 3 S 2 ) × R + → R of (59) such that
To get some intuition on the concept of viscosity solutions on R n we note that the name viscosity solutions refers to the propagation of viscous fluids, where propagation fronts can not propagate through each other, cf. Fig. 13 According to the following theorem, we have that the unique viscosity solutions of our morphological scale spaces are given by erosions/dilations with the corresponding Green's functions.

are given by (− case) left-invariant erosion
and with the corresponding morphological Green's function.
The morphological Green's functions are given by −k with spatial arc length s given by s(τ ) := τ 0 ẋ(τ ) dτ , and with R 3 S 2 -"erosion arc length" given by The exact erosion/dilation kernels satisfy the semigroup property for all t, s > 0, where we recall Eq. (25).
Proof See Appendix A.
In Eq. (70) we can use any parametrization of a smooth curve τ → γ (τ ) in the sub-Riemannian manifold (SE(3), dA 3 , dA 6 ). For example if τ = s (spatial arc-length) we find dp ds (s) = |γ 1 (s)| 2 + |γ 2 (s)| 2 where L denotes the total length of x(·) and where the Lagrangian L η is evaluated along the curve where the components of the tangent vector w.r.t. our moving left-invariant frame of reference are given bẏ Proof This follows by the homogeneity of the Lagrangian, the chain law, Eqs. (72) and (49).
Spatial arc length parametrization of the curves in sub-Riemannian manifolds within SE(d), d = 2, 3 is often preferable over sub-Riemannian arc length parametrization, as it involves less special functions, cf. [33, Appendix A], [70], and the parametrization breaks down precisely at the cusps (which are exactly the points where both global and local optimality is lost [15]) recall Fig. 11. However, for generalizing the results on viscosity solutions of HJB-equations on R 3 , [5,17,24,39] to R 3 S 2 the p-parametrization seems to be more suitable.

Tangible Approximations of Linear and Morphological Green's Functions
The derivation of the exact formula for the morphological kernel in Eq. (69) is very complicated: One first would have to derive the curve minimizers via reduction of a complicated Pfaffian system (following a Lagrangian approach similar to the derivation of the sub-Riemannian geodesics in [37, Appendix G], [16], [44]) or the Pontryagin maximum principle (following a Hamiltonian approach [4]). A direct classical variational Lagrangian approach (similar to [59]), without relying on symplectic geometry can also be done. It yields the same curves, but only after highly cumbersome bookkeeping and computations, [44]. On top of that, one has the problem that once the stationary curves are derived by integration of the Lagrangian or Hamiltonian system, one must solve a boundary value problem that involves at least 1 special function (similar to the SE(2)-case [15]) and in order to compute the morphological kernel in Eq. (69) one must resort to a quite involved numerical algorithm [44] for every (y, n) ∈ R 3 S 2 .
A more reasonable, but less accurate, alternative would be to apply an upwind finite difference scheme of (59) that is explained in Sect. 9.3. However, this would require interpolations in the rotated, translated and reflected kernels in Theorem 3. Therefore, if it comes to fast practical kernel implementations (both for the morphological and linear case) and analysis we prefer to use the approximations presented in the next theorem.
Proof We apply the theory of weighted sub-coercive operators on Lie groups here [78,Proposition] and/or the general Ball-box theorem [12,13] that provides a local uniform estimate for sub-Riemannian metrics. For hypo-elliptic diffusion we must consider the sub-Riemannian manifold: (SE(3), dA 1 , dA 2 , dA 6 ) (recall Fig. 10). This induces the following filtration: Now there exist Gaussian approximations (77) for 0 < a 1 ≤ a 2 and b 1 ≥ b 2 > 0 for hypo-elliptic Green's functionsP D 33 ,D 44 ,D 55 t with respect to the sub-Riemannian metric [78]. Now by the ball-box theorem the sub-Riemannian metric d on (SE(3), dA 1 , dA 2 , dA 6 ) can be estimated by the weighted modulus with A i = A i | e , c > 1, for all g close to e. For the explicit expression of the log-coefficients c i in terms of our Euler angles (17) see [36,Eq. (76)]. Here we stress that the exponential map exp : T e (SE(3)) → SE (3) is surjective. In case of the nilpotent Heisenberg approximation of the Euclidean motion one can get a grip on the constant c > 1 (which is not that far from 1 see [32,36]). The precise value of c > 1 is not crucial for this theorem as it can be taken into account by a rescaling of time t → c −1 t. Now the sub-Riemannian balls in (SE(3), dA 1 , dA 2 , dA 6 ) are similar to their ball-box estimates non-smooth, whereas the heat-kernels are smooth due to the Hormander Theorem [49]. Therefore we use the locally equivalent smooth modulus Now the result (75) follows by Eq. (78), (77), (79) and rescaling of the Lie-algebra generators. Note that a subtlety arises here: the logarithm is welldefined on SE(3) but it is not well-defined on R 3 S 2 , it is therefore important that in our approximations one takes a consistent section in the principal fiber bundle (SE(3), (2)) with total space SE(3), base manifold R 3 S 2 and structure group {0} SO(2) acting via the right action R h g = gh and with projection g → [g] = g({0} SO (2)). Since we identified (0, e z ) with the unity e = (0, I ) ∈ SE(3) we must setα = 0 in all Lie-algebra coefficients c i , since in our approximation we should stick to the consistent cross section, cf. [36,Chap. 7].
For the result (74) on approximations of morphological Green's functions, we first consider the case η → ∞, D 11 = D 44 = 1. In this case the Lagrangian is homogeneous and we consider the sub-Riemannian distance δ given by Eq. (46) defined on the sub-Riemannian manifold (SE(3), dA 3 , dA 6 ): Note that for erosions we need a different sub-Riemannian manifold than for diffusion, recall Fig. 10. Now again by the Ball-box theorem and the equivalent smooth modulus trick δ is locally equivalent to the following weighted modulus on the Lie-algebra 4 c 1 2 + c 2 2 + c 4 2 + c 5 2 2 + c 3 2 + c 6 2 .
The result for 1 2 < η ≤ 1 follows by the Fenchel-transform on L (SE(3)) and Theorem 8. The result (74) for general g ii = (D ii ) −1 > 0 follows by re-scaling of the Lie-algebra basis. Finally, the positive Green's function for erosion is minus the negative Green's function for dilation.

Remarks
• Erosion ball in SE(3) around the unity, i.e. the ball g ∈ SE(3) | δ(g, e) < r , with r > 0 small, are equivalent to balls g ∈ SE(3) | 4 |c 1 | 2 + |c 2 | 2 + |c 4 | 2 + +|c 5 | 2 2 + |c 3 | 2 + |c 6 | 2 < r induced by weighted modulus on the Lie-algebra of SE(3) (with c i = c i (g) the logarithm coefficients). This underlies our asymptotical formula (74). This gives us a simple analytic approximation formula for the balls in the sub-Riemannian manifold (SE(3), dA 3 , dA 6 ), without determining the minimizing curves (i.e. geodesics) in Eq. (46). • Next we provide some background in the weighting in the weighted modulus. In general it is not possible to connect any two points, say e and g, via a horizontal exponential curve. This means that not every element (y, R e z ,γ R e y ,β ) ∈ SE (3) which explains the different weighting of the missing Liealgebra directions. • From a practical point of view one typically wants to erode in the sub-Riemannian manifold (SE(3), dA 3 , dA 6 ) and to dilate in sub-Riemannian manifold (SE(3), dA 1 , dA 2 , dA 6 ). In this case the erosion process sharpens the fibers, whereas the latter process would dilations of fibers (rather than shortening them). One can combine these processes in a single erosion-dilation process on SE(3). • The underlying idea in the logarithmic approximation within a (sub-)Riemannian manifold in SE (3) is that one keeps the non-zero coefficientsγ i constant. This is equivalent to considering auto-parallels of the Cartan connection on the underlying principal fiber bundle [33,36]. The auto-parallels of the Cartan connection are the (horizontal) exponential curves, which due to the torsion of the Cartan connection do not coincide with the (sub)-Riemannian geodesics. However, due to the commutator relations and the CBH-formula this torsion kicks in as a second order effect and locally we can approximate with an aligned exponential curve. 10 9 Implementation of (Morphological) Scale Spaces on R 3 S 2 We distinguish between two implementations for morphological scale spaces: such that the spatial erosions are given by Advantages kernel implementations:

Geometric Erosion Parameters
The erosions given by Eq. (59) involve four fundamental parameters, t, D 11 , D 44 and η. However, the erosion equations can be re-scaled such that D 11 = 1, leaving only 3 degree of freedom. The fraction D 11 /D 44 determines the relative strength of angular erosion compared to spatial erosion and this parameter has physical dimension [Length] 2 . After rescaling the underlying sub-Riemannian metric tensor becomes A large value of D 11 /D 44 yields a domination of spatial erosion, a small value of D 11 /D 44 yields a domination of angular erosion. The stopping time parameter t > 0 indicates the size of the erosion-kernels. Although not considered here, statistical approaches could be used to estimate a suitable stopping time from a given DW-MRI data-set. To get a visual impression of the typical effect of (D 11 /D 44 , t) on the erosion output, we refer to Fig. 14. The parameter η ∈ ( 1 2 , 1] determines the power of the Hamiltonian. For η ↓ 1 2 we get a homogeneous Hamiltonian and a flat morphological approximation kernel (take η ↓ 1 2 in Theorem 11): The typical effect of including η is depicted in Fig. 15 where we applied angular erosion only (i.e. D 11 = ∞), see Eq. (82). Due to the fact that the approximation of the erosion kernel in Theorem 11 is a quadratic form on the logarithmic coefficients that depend on goniometric functions of the anglesβ andγ along the cross-sectionα = 0 the erosion kernel is not very suited for a glyph-field visualization. To get an impression of the shape of the balls in (SE(3), dA 3 , dA 6 ) it is more intuitive to plot the corresponding Gaussian estimation of the hypo-elliptic diffusion

Kernel Implementations
Kernel implementation for dilations/erosions on R 3 S 2 are entirely analogous to kernel implementations of linear diffusion on R 3 S 2 , except for a change in the algebra. For dilations one uses the (max, +)-algebra whereas with diffusion one uses the (+, ·)-algebra. For erosions one uses the (min, +)-algebra. This is directly observed by comparing Lemmas 1 and 5. We use the fast tangible morphological kernel approximations in Theorem 11.

Finite Difference Implementations
Similar to our finite difference schemes for left-invariant diffusions [26,36] we use a Euler-forward first order time integration scheme. However, for morphological scale spaces we use a so-called upwind-scheme, where the sign of the central difference determines whether a forward or backward difference is taken.
For further details on implementation of the left-invariant finite differences on R 3 S 2 and the involved interpolation we refer to [26].

Adaptive Morphological Scale Spaces
It is very interesting to consider adaptive dilations on R 3 S 2 , in a similar way as with adaptive diffusions [42]. This would generalize the efficient dilation [19,63] of DTI-tensors to (diffused) high angular resolution DW-MRIimages defined on R 3 S 2 . This could be done by multiplying the left-invariant gradient dW , Eq. (60), with a data adaptive positive definite help-matrix in SO(3) ⊕ SO(3), before substituting it into the Hamiltonian in Eq. (59). Here the idea of gauging the local reference frame to include curvature adaptation and deviation for horizontality based on locally optimal exponential-curves in SE(2) [33, Chap. 3.4], could be generalized to SE(3). However, adaptive dilations and erosions on the full space R 3 S 2 are beyond the scope of this article.
Here we just provide a simple mechanism for dataadaptive angular erosion and dilation as a first step towards adaptive erosions on R 3 S 2 . In the erosion evolution (59) one can include adaptivity by making D 44 dependent on the local Laplace-Beltramioperator with φ a monotonic, odd function and The intuitive idea here is to create a soft transition between dilation and erosion on a glyph x + μ{U(x, x)n | n ∈ S 2 }. In the vicinity of a spherical maximum the S 2 is positive, whereas in the vicinity of a spherical minimum the S 2 is negative. The constant c determines where the transition should take place. See Fig. 17 where we set φ(x) = D 44 x

The Underlying Probability Theory
In this section we address the underlying probability of scale spaces on R 3 S 2 . We first consider probability theory underlying linear scale spaces on R 3 S 2 , and then a probabilistic interpretation of morphological scale spaces on R 3 S 2 .

Brownian Motions on R 3 S 2
Next we formulate a left-invariant discrete Brownian motion on SE(3). The left-invariant vector fields {A 1 , . . . , A 6 } form a moving frame of reference to the group. There are two ways of considering vector fields. Either one considers them as differential operators on smooth locally defined functions, or one considers them as tangent vectors to equivalent classes of curves. These two viewpoints are equivalent, for formal proof see [9,Prop. 2.4]. So far we have used the first way of considering vector fields, but in this section we prefer to use the second way. We will write {e 1 (g), . . . , e 6 (g)} for the left-invariant vector fields (as tangent vectors to equivalence classes of curves) rather than the differential operators {.A 1 | g , . . . , A 6 | g }. We obtain the tangent vector e i from A i by replacing ∂ x ↔ (1, 0, 0, 0, 0, 0), 0, 1, 0, 0, 0, 0), 0, 0, 1, 0, 0, 0), 0, 0, 0, α cos β cos γ, α cos β sin γ, −α sin β), 0, 0, 0, α cos γ, α sin γ, 0), 0, 0, 0, cos γ sin β, sin γ sin β, cos β), (86) where we identified SO(3) with a ball with radius 2π whose outer-sphere is identified with the origin, using Euler an-gles R e z ,γ R e y ,β R e z ,α ↔ αn(β, γ ) ∈ B 0,2π . Next we formulate left-invariant discrete random walks on SE(3) expressed in the moving frame of reference {e i } 6 i=1 given by (86): , ε i,n+1 ∼ N (0, 1) and where the stepsize equals s = s N and where a := 5 i=1 a i e i denotes an apriori spatial velocity vector having constant coefficients a i with respect to the moving frame of reference {e i } 5 i=1 (just like in (34)). Now if we apply recursion and let N → ∞ we get the following continuous Brownian motion processes on SE (3): σ ji e j | (Y (τ ),N (τ )) dτ, with ε i ∼ N (0, 1) and (X(0), N (0)) ∼ U and where σ = √ 2D ∈ R 6×6 , σ > 0. Now if we set U = δ 0,e z (i.e. at t = 0, we have only one orientated particle located at (0, e z )) then suitable averaging of infinitely many random walks of this process yields the transition probability density with Fokker-Planck equation (53), for D 11 = 0. Instead of considering the probability density W (y, n, t) of finding a random walker at time t at position y and orientation n given that it started from the initial distribution U , one can also consider the probability density P λ (y, n) of finding a random walker at position y with orientation n given that it started from the initial distribution U regardless its traveling time T ∼ NE(λ). Then by Eq. (7.1) we must apply the resolvent operator to the initial condition: The (resolvent) Green's functions and the corresponding stochastic processes for contour-completion and contour enhancement are illustrated in Fig. 3.

The Time Integrated Direction Process
In the previous work [36] the resolvent processes of contour enhancement are studied and according to [36, Theorem 2] they are related to Tikhonov regularization on R 3 S 2 . The corresponding resolvent kernels do not satisfy the semi-group property (cf. Eq. (58)) and the question rises what happens when alternating Tikhonov regularization steps. The same question is particularly important with contour completion, where the resolvent processes integrate the moving propagation fronts of the convection-diffusions. The solution P λ (y, n) = (R D,a λ * R 3 S 2 U)(y, n) of (90) represents the probability density of finding a random walker at position y with orientation n given that it started from the initial distribution U regardless its traveling time. In a Markov process, traveling time is memoryless and thereby negatively exponentially distributed T ∼ NE(λ).
There is however, a practical drawback induced by this memoryless property: Both the time-integrated resolvent kernel of the direction process and the time-integrated resolvent kernel of the enhancement process suffer from a serious singularity at the unity element. For example, in [37, Appendix D] we derived an asymptotical formula for the resolvent contour enhancement kernel, from which we deduce R D 33 ,D 44 λ (y, n) ∼ 1 |g| 6 for |g| D 33 ,D 44 1, and where |g| D 33 ,D 44 = |(y, R n )| D 33 ,D 44 is the weighted modulus on SE(3) given by with c i the logarithmic coefficients given by (76), recall Theorem 11. Consequently, the kernels are in practice too much concentrated around the singularity for visually appealing results. Therefore, we introduce a k-fold iterative approach and address the three related questions raised at the end of Sect. 7.1.

The Iteration of Resolvent Operators
The sum T = k i=1 T i of k independent negatively exponentially distributed random variables T i ∼ NE(λ) is Gamma distributed: with temporal convolutions given by Now for the sake of illustration we set k = 2 for the moment and we compute the probability density of finding a random walker G T ∈ R 3 S 2 with traveling time T = T 1 + T 2 at position y with orientation n given that it at t = 0 started at (0) with orientation e z . Basically this means that the path of the random walker has two stages, first one with time T 1 ∼ NE(λ) and subsequently one with traveling time T 2 ∼ NE(λ) and we have R D,a λ,k=2 (y, n) where * = * R 3 S 2 denotes convolution on R 3 S 2 as in Theorem 2. By induction this can easily be generalized to the general k-step approach: An alternative derivation of the result arises by which only holds in the distributional sense, where with respect to the crucial second equality we recall Eq. (91). We have proved the following result.

Theorem 12
The probability density function (Green's function) of this stochastic process is given by Figure 20 shows some experiments of contour completion on an artificial data set containing fibers with a gap we would like to complete, for various parameter settings of (λ, k), where T i ∼ NE(λ) and k denotes the iteration index of the time integrated contour-completion process. In principle this scheme boils down to R 3 S 2 -convolution (applying Lemma 1 and Theorem 2) with the Green's functions depicted in Fig. 19. For a fair comparison we kept the expected value constant, that is E(T ) = E(T 1 + · · · + T k ) = k λ = 4 in Fig. 19 and E(T ) = E(T 1 + · · · + T k ) = k λ = 2 in Fig. 20, which roughly coincides with half the size of the gap. Note that the results hardly change after two iterations, as the graphs of the scaled Gamma distributions Γ (·; k 2 , k)/Γ ((2(k − 1)/k); k 2 , k) are similar for k = 2, 3, 4. We conclude that the k-step approach with k > 1 provides much more effective completion of gaps within fibers in R 3 S 2 than the case k = 1, as can be observed in Figs. 19 and 20. This mainly due to the fact that the order of the singularity at the origin decreases (it takes 5 steps to get rid of the singularity [37, Chap. 10.1.1, Appendix D]).

Bellman Processes on R 3 S 2
In this subsection we present a short overview of cost processes on SE(3). The mapping I → lim ε→0 log I defines a morphism of the (+, ·)-algebra to the (min, +)-algebra on R + (the co-domain of our DW-MRI images) and it is indeed readily verified that lim ε→0 log a · b = a + b and Using this map, various notions of probability calculus can be mapped to their counterparts in optimization problems. Next we mention the following definitions as given by Akian, Quadrat and Viot in [5] adapted to our case.

Definition 9
The decision space is the triplet (SE (3) One can formulate related concepts such as independent decision variables, conditional cost, mean and characteristic function of a decision variable, in the same way as in probability theory via the morphism of the (+, ·)-algebra to the (min, +)-algebra. The Laplace or Fourier transform in (+, ·)-algebra corresponds to the Frenchel transform in the (min, +)-algebra. Now we present the decision counterpart of the Markov processes, namely Bellman processes.
where c is called the transition cost which is a map from for all (t, g) ∈ R + × SE(3), and where c 0 is some cost density on SE(3).
We set the following transition cost density: Then the marginal cost for a Bellman process G t on (SE (3), O(SE(3)), C) to be in a state g at time t given an initial cost c 0 ,m(g, t) := C(G t = g) satisfies the Bellman equation where F L(SE(3)) denotes the Fenchel transform (recall Definition 6) on the Lie-Algebra L (SE(3)). Furthermore, we used This Bellman equation for the cost process coincides with the Hamilton-Jacobi-Bellman equation, Eq. (59) on R 3 S 2 , whose viscosity solution is given by morphological convolution with the corresponding morphological Green's function as shown in Appendix A.

Experiments
First we provide a summary of the experiments we considered so far. Then we show two more experiments on diffusion-weighted MRI-images of the brain of a healthy volunteer. Finally, we draw some general conclusions from the experiments.
Throughout this article, DW-MRI images are scanned in N o ≈ 50 gradient directions while using a b-value of 1000 s/mm 2 . We restricted ourselves to these values, since our evolutions perform similarly for other reasonable bvalues and other number of gradients N o used in the MRIscanning, as shown by Prckovska et al. [65]. In general DTIimages do not require N o ≈ 50, due to the Gaussian assumption Eq. (1), but for HARDI-imaging with b = 1000 s/mm 2 , N o ≈ 50 is reasonable.
In Figs. 2, 4 and 7 one can see that high-frequency components are removed by left-invariant diffusion on DW-MRI while preserving crossings, allowing a well-posed subsequent erosion. Such a subsequent erosion is depicted in Fig. 7 (lower-left figure). The advantage of a subsequent erosion is also depicted in Fig. 21 where we show multiple viewpoints in 3D.
Left-invariant convection-diffusion for contour-completion can be used to close gaps in DW-MRI data as can be observed in Fig. 20 (using the iterated resolvent kernels depicted in Fig. 19 according to Theorem 12).
Left-invariant diffusions for contour enhancement, Eq. (53) (in combination with grey-value transformations) can be used to extrapolate fiber crossings from DTI (see Fig. 5), whereas they remove spurious crossings in HARDI as can be observed in Fig. 6. However, these grey-value transformations suffer from the drawback of explosive enlargement of large glyphs illustrated in Fig. 4 and typically require masking of these large glyphs (typically present at ventricles). To avoid the drawbacks of greyvalue transformations we include a final experiment of concatenating leftinvariant erosion and diffusion on DW-MRI without the involvement of greyvalue transformations. See Fig. 23 where we compare the concatenation of squaring, diffusions and squaring to the concatenation of diffusion and erosion. For the erosions we applied the upwind-scheme of Sect. 9.3, for the diffusions we used the finite-difference schemes explained in our numerical work [26]. Summarizing, we conclude: • R 3 S 2 -erosions are both useful for visual sharpening of the data (Fig. 21) and for pre-processing before applying R 3 S 2 -diffusion (Fig. 23). • Concatenation of linear and morphological scale spaces on R 3 S 2 provide a more effective sharpening of the fibers than a pseudo-linear scale space approach, e.g. compare Fig. 23

Conclusion
We have developed crossing preserving, rotation-and translation covariant scale spaces on DW-MRI images. The underlying evolutions are convection-diffusion equations and Hamilton-Jacobi-Bellman equations of respectively stochas-tic and cost processes, on the space of positions and orientations R 3 S 2 . These scale spaces are expressed in a moving frame of reference allowing (hypo-elliptic) diffusion along fibers and erosion orthogonal to fibers. They extrapolate complex fiber-structures (crossings) from DTI, while reducing non-aligned crossings in HARDI. They can be implemented by finite difference methods [26], or by convolutions with analytic kernels. We have addressed the nine objectives arising from previous work which are listed in the introduction. The main results of this article are Theorems 10, 11 and 12. This is where Theorem 11 comes into play. The logarithmic modulus approximations provide lower and upper-bounds for the exact morphological Green's functions. By Theorem 9 the exact Green's functions arise by geodesically equidistant propagation of the morphological delta, δ C given by Eq. (25). Moreover, the logarithmic modulus approximations arise by freezing the coefficients p →γ i (p). These approximations are accurate for small h. In fact, we have the following generalization of Eq. (95) to R 3 S 2 : Lemma 5 Let c i ∈ R, i = 1, 2, 4, 5 and η > 1 2 and define L η :  6 ). This map is a diffeomorphism until the first cusp-time arises (likewise in Fig. 11). Application of first order Taylor-approximation around s = 0 yields and η > 1 2 implies that this is of the order where ψ −1 (0) = 0 and in the interior of a minimizing curve we have finite curvature so ψ −1 (h) = O(h) and the first term vanishes as h ↓ 0 and the result follows from the definition of the morphological kernel, Eq. (63).

Remark
In the subsequent proofs we need a sober notation and we just writẽ for all g ∈ SE(3), (y, n) ∈ R 3 S 2 and all t ≥ 0, for the corresponding morphological Green's function on the group SE(3), given in Eq (96).
A second issue is the semi-group property. The PDE for morphological scale spaces is a non-linear evolution and the question whether the semi-group property (57) that holds for linear left-invariant evolutions has a morphological counterpart, i.e. whether Eq. (71) holds. In the R n case, this directly follows from the isomorphic properties of the Cramer transform [5,17], but in R 3 S 2 this is difficult. for all g ∈ SE(3) and t ≥ 0, then the semi-group property follows from for all τ ∈ [0, t] and the associativity of morphological convolutions. So it remains to prove Eq. (97). To this end we follow a similar approach as in [56, Theorem 1] (avoiding specific scaling properties of the Lagrangian). By Theorem 9 it follows that the Green's functions are induced by geodesically equidistant propagation of the morphological delta, δ C given by Eq. (25). More precisely, they are given by Eq. (69).
In general one has k t v −1 g ≤k τ v −1 q +k t−τ q −1 g for all g, h, v ∈ SE(3) and all τ ∈ [0, t], since the concatenation of two optimal curves yields an admissible curve over which the left hand side is optimized. Here equality is obtained if v −1 q is on the minimizing curve between e and v −1 g, i.e. by left-invariance g is on the minimizing curve between e and v −1 g. Now due to continuity and convexity of (c 1 , c 2 , c 4 , c 5 ) → L η (c 1 , c 2 , c 4 , c 5 ) the infimum in (97) is actually a minimum and therefore we can choose v ∈ SE(3) such that W (q, τ ) =k τ v −1 q +Ũ(v).
Then by Eq. (98) one has W (g, t) ≤k t v −1 g +Ũ(v) ≤k τ v −1 q +k t−τ q −1 g +Ũ(v) =W (q, τ ) +k t−τ q −1 g for all q ∈ SE(3) and thereby (by taking the infimum over all q) we obtain In order to prove (57) it remains to be shown that Let w ∈ SE(3) such that Now take q as a point along the minimizing curve between g and w. Then we havẽ W (g, t) =k τ w −1 q +k t−τ q −1 g +Ũ(w) ≥k t−τ q −1 g +W (q, τ ) ≥ min r∈SE (3) k t−τ r −1 g +W (r, τ ) from which the result follows.
To complete the proof of Theorem 10 we must prove the following theorem.
Proof The proof consists of two parts, first we must show that they are indeed solutions and then we show that they are viscosity solutions. We will only consider the erosion case since the dilation case can be treated analogously.
Part II Next we verify that erosion (102) with the Green's function k t indeed satisfies (64) at a location (y 0 , n 0 , t 0 ) ∈ (R 3 S 2 ) × R + where W − V attains a local maximum. The dilation case follows analogously. First of all it follows by left-invariance (and the fact that SE(3) acts transitively on R 3 S 2 , cf. Eq. (6)) that without loss of generality we can restrict ourselves to the case (y 0 , n 0 ) = (0, e z ). Secondly, by the semigroup property in Lemma 6, we only need to consider the case t 0 = 0, while dropping the constraint t > 0.

Appendix B: Propagation of geodesically equidistant surfaces
In this section we provide a proof of Theorem 9. First of all we note that Eq. (62) in Definition 7 can be written out in more explicit form via the chain-rule for differentiatioñ from which we deduce that D XLη (·, X) = dW (·, t) where D X denotes the derivative w.r.t. a vector field X tangential to a transversal minimizing geodesic γ given by

Appendix C: Pseudo-linear Scale Spaces
So far we have considered anisotropic diffusions aligned with fibers and erosions orthogonal to the fibers. As these two types of left-invariant evolutions can be concatenated, cf. Fig. 4, the natural question arises whether there exist a single evolution process that combines erosion/dilation and diffusion. Moreover, a simple alternative to visually sharpen distributions on positions and orientations is to apply monotonic greyvalue transformations, for instance via the power operatorχ q χ q (U ) (y, n) = U(y, n) q , where q ≥ 1, where we recall its drawback illustrated in Fig. 4. Conjugation of the diffusion operator with a monotonically increasing grey-value transformation χ : is related to simultaneous erosion/dilation and diffusion, cf. [40]: with C ∈ R, t ≥ 0, are solved by where χ C : R + → R + denotes the monotonic greyvalue transformation in the co-domain given by whose inverse is given by Proof This is a straightforward generalization of the derivations in [40] where one replaces the shift-invariant generators ∂ x i on R n by the left-invariant generators {A i } 6 i=1 on SE (3). The derivations in [40] rely on the chain-law for monotonic transformations in the codomain and the product rule for differentiation and these rules also apply to {A i } 6 i=1 . Finally, the solution of the pseudo-linear scale space is legal (Definition 3) iff the diffusion operator e tQ D,a (A) is legal.
Remark The drawback of this correspondence, is that our diffusions primarily take place along the fibers, whereas our erosions take place orthogonal to the fibers, recall Fig. 10. The pseudo-linear scale spaces always erode/dilate (depending on sign(C)) and diffuse in the same direction.