New exact and numerical solutions of the (convection–)diffusion kernels on SE(3)

Direction process We consider hypo-elliptic diﬀusion and convection–diﬀusion on R 3 (cid:2) S 2 , the quotient of the Lie group of rigid body motions SE(3) in which group elements are equivalent if they are equal up to a rotation around the reference axis. We show that we can derive expressions for the convolution kernels in terms of eigenfunctions of the PDE, by extending the approach for the SE(2) case. This goes via application of the Fourier transform of the PDE in the spatial variables, yielding a second order diﬀerential operator. We show that the eigenfunctions of this operator can be expressed as (generalized) spheroidal wave functions. The same exact formulas are derived via the Fourier transform on SE(3). We solve both the evolution itself, as well as the time-integrated process that corresponds to the resolvent operator. Furthermore, we have extended a standard numerical procedure from SE(2) to SE(3) for the computation of the solution kernels, that is directly related to the exact solutions.© 2017 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

In [10], a Fourier-based method is presented to compute the kernel of the hypo-elliptic heat equation on unimodular Lie groups and examples are provided for various Lie groups.Three approaches to derive the exact solutions for the heat kernel for SE (2) have been proposed in SE (2) in [18,6], of which the first is equivalent to the approach in [10] and provides a solution in terms of a Fourier series.This approach will be extended in this paper to SE (3) and the connection to the Fourier transform on SE(3) will be made.The second approach of [18,6] provides a solution in terms of rapidly decaying functions and in the third approach this series is explicitly computed in terms of Mathieu functions.
In image analysis of 2D images, both the convection-diffusion and the diffusion PDEs can be used for crossing-preserving enhancement of elongated structures, after the image is lifted from R 2 to the group SE(2) via an invertible orientation score [19], other types of lifting operators [3,20], or using a semi-discrete variant [21].Enhancement of such structures is important in for example retinal imaging to improve the detection of blood vessels, which is challenging due to the low contrast and noise in the images.
Recently, both the diffusion and the convection-diffusion PDE have become of interest for enhancement of diffusion-weighted Magnetic Resonance Imaging (dMRI).The diffusion orientation distribution function (ODF) or the fiber orientation distribution (FOD) function can be considered to be a function on a domain in R 3 × S 2 , indicating at each position the estimated diffusion profile or distribution of fibers on the sphere S 2 .Such functions are usually position-wise estimated from the dMRI data.The advantage of processing with these PDEs is that they induce a better alignment between local orientations and their surrounding.In [5] the 3D extension of Mumford's direction process was used enhance dMRI data with the aid of stochastic completion fields.In [22] it is shown that the convection-diffusion kernel can be used to obtain asymmetric, regularized FODs.The practical advantages of the diffusion process for regularization of dMRI data, in particular better fiber tractography results and improved connectivity measures, are given in [23,24,25,26].
Although finite-difference implementations [27] and kernel approximations [28] of the PDEs on R 3 S 2 (for which the formal definition follows later) exist, so far no exact expressions are known.The derivation of these exact solutions will be one of the main results of this paper.Other contributions of this work are summarized at the end of this introduction.We first provide more details on the mathematical setting, established in previous work [28].

Left-Invariant Convection-Diffusion Operators on SE(3)
Let Ũ : SE(3) → R + be a square integrable function defined on the Lie group of rigid body motions SE(3) = R 3 SO(3), which is the semi-direct product of the translation group R 3 and the rotation group SO (3).We consider Ũ to be a given input, that requires regularization or enhancement of certain features.We use a tilde throughout the paper to indicate functions and operators on the group SE(3) (in contrast to functions/operators on a group quotient later in the paper, that we denote without a tilde).A possible approach is to use particular evolution equations that are special cases of the following general evolution process: , for all g ∈ SE(3), t ≥ 0, W (g, 0) = Ũ (g), for all g ∈ SE(3). ( Here Q is the generator of the evolution on the group, where we restrict ourselves to generators such that the evolution becomes a linear, second order convection-diffusion process.Moreover, Q should be a left-invariant operator and is composed of the left-invariant differential operators A = (A 1 , . . ., A 6 ).These vector fields are obtained via a pushforward A i | g = (L g ) * A i , where {A i } 6 i=1 is a choice of basis in the tangent space T e (SE(3)) at the unity element e = (0, I), where A 1 , A 2 , A 3 are spatial tangent vectors and A 4 , A 5 , A 6 are rotational tangent vectors.Here L g is the left-multiplication L g q = gq.If Q is a generator of a linear convection-diffusion PDE, then its general form can be written as: Here D = [D ij ] ∈ R 6×6 is a symmetric positive semi-definite 6 × 6-matrix, and a = (a i ) 6 i=1 ∈ R 6 .The explicit formulas for the A i require two charts [28], but these are not of crucial importance in this work, because of the particular configuration in which the A i 's occur in the two generators (two special cases of (2)) that we will consider.We use Υt ( Ũ ) = W (•, t) to denote the solution of the above evolution equation, so Υt : L 2 (SE(3)) → L 2 (SE(3)).The evolution can then formally be solved by convolution with the integrable solution kernel or impulse response Kt : SE(3) → R + : Υt ( Ũ )(g) = (e t QD,a Ũ )(g) = ( Kt * SE(3) Ũ )(g) = where we use the Haar measure on SE (3), which is the product of the Lebesgue measure on R 3 and the Haar measure on SO (3).From the dMRI application, we often encounter functions defined on R 3 × S2 instead of the functions on SE(3) as above.These functions are called Fiber Orientation Distributions (FODs) and they estimate at each position in the space R 3 the probability of having a fiber in a certain orientation.As shown in [25], to appropriately regularize these functions, a coupling between positions and orientations is required.This can be obtained by embedding the space R 3 × S 2 in SE (3).To this end, we define the space of coupled positions and orientations as a Lie group quotient in SE(3): R 3 S 2 := SE(3)/({0} × SO(2)). ( For two elements g = (x, R), g = (x , R ) ∈ SE(3), the group product is defined as gg = (x + Rx , RR ).Because this group product influences the group action in R 3 S 2 , we use the semi-product notation (even though this is usually reserved for the semi-direct product of groups).The group action of g Within the quotient structure R 3 S 2 two elements g, g (as above) are equivalent if with e z = (0, 0, 1) T the reference axis1 .Using the group action, equivalence relation ( 6) amounts to: Thereby, an arbitrary element in R 3 S 2 can be considered as the equivalence class of all rigid body motions that map reference position and orientation (0, e z ) onto (x, n) (with arbitrary x ∈ R 3 and n ∈ S 2 ).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).
Functions U : R 3 S 2 → R + are identified with functions Ũ : SE(3) → R + by the relation where R n is any rotation matrix such that R n e z = n (not to be confused with the notation R n,φ , that specifies the rotation by its axis n and angle of rotation φ).This means that the functions Ũ that we consider have the invariance Ũ (y, R) = Ũ (y, RR ez,α ), for all α ∈ (0, 2π].

Legal Diffusion and Convection-Diffusion Operators
For functions defined on the quotient R 3 S 2 the same machinery as in Section 1.2 can be used, although some restrictions apply.In [29] it is shown that operators on Ũ are legal (i.e., they correspond to well-defined operators on U that commute with rotations/translations) if and only if the following holds: with left-regular action (L g Ũ )(g ) = Ũ (g −1 g ) and right-regular action (R g Ũ )(g ) = Ũ (g g).For the diffusion case, an equivalent statement is that the diffusion matrix D is invariant under conjugation with Z α = R ez,α ⊕ R ez,α ∈ SO(6), i.e.Z −1 α DZ α = D.Among the few legal, left-invariant diffusion and convectiondiffusion generators on R 3 S 2 , recall (2), are the pure hypo-elliptic diffusion case and the hypo-elliptic convection-diffusion case, for details see [28,29].The first case corresponds to the forward Kolmogorov equation for hypo-elliptic Brownian motion, the second case corresponds to the forward Kolmogorov equation for the direction process (in 3D), as illustrated in Fig. 1.We denote the evolution generators of these two cases on the group and on the quotient with Qi and Q i , i = 1, 2, respectively.They are defined as: with D 33 > 0, D 44 = D 55 > 0. The generator Q 1 acts on sufficiently smooth functions Ũ and can be identified with a generator Q 1 acting on functions U : Here ∇ R 3 denotes the gradient on R 3 , and ∆ S 2 is the Laplace-Beltrami operator on the sphere S 2 .Also the following (hypo-elliptic) convection-diffusion generator can be identified with a legal generator on R 3 S 2 : again with D 44 = D 55 > 0. The corresponding generator Q 2 acting on sufficiently smooth functions U is defined as: Remark 1.1.The way Qi and Q i act on functions W and W , respectively, is related by the following identities: As a result, the following relation holds: regardless of the choice of rotation for R n (mapping e z onto n).
By the previous remark, we are no longer concerned with the left-invariant vector fields A i , and we focus on the following two PDE systems using (11) and (13): In particular, we want to find an exact solution and approximations for the impulse response or convolution kernel of these PDEs, i.e., the solution with initial condition U (y, n) = δ (0,ez) (y, n).We denote the linear and bounded evolution operator that maps U (y, n) on W (y, n, t) by

Relation with Stochastic Processes, Time-Integrated Processes and Resolvent kernels
Any evolution equation corresponding to a generator Q D,a as above, is in fact a Fokker-Planck equation for the evolution of a probability density function.The (convection)-diffusion processes generated by Q 1 and Q 2 were called contour enhancement and contour completion in [28], respectively, because of the stochastic processes they relate to, see Fig. 1.The random walks in the figure are simulated according to the formal definition of the stochastic processes in Appendix E. The convolution kernels of the PDEs can also be obtained with a Monte Carlo simulation by accumulating infinitely many random walks starting at (0, e z ) that obey the underlying stochastic process and have traveling time t, as illustrated in Fig. 1.The figure also displays the SE(2) equivalents of the processes, because they are easier to visualize and interpret.The SE(2) process for contour completion is better known as Mumford's direction process [12] and finds its application in computer vision.
Remark 1.2 (Glyph field visualization).In Fig. 1 and in later figures in the paper, we make use of so-called glyph field visualizations of probability distributions U on R 3 S 2 .A glyph at a grid point y ∈ c Z 3 , c > 0, is given by the surface {y + νU (y, n)n | n ∈ S 2 }, for a suitable choice of ν ∈ R. Usually we choose ν > 0 depending on the maximum of U , such that the glyphs do not intersect.
When we are interested in the position of a random walker regardless of the traveling time, we can condition on the traveling time and integrate.If we assume to have exponentially distributed traveling time t ∼ Exp(α), with mean 1/α, α > 0, then the probability density function of a random particle given its initial distribution U can be written as: This shows that the time-integrated process relates to the resolvent operator of Q D,a , in particular for Q 1 and Q 2 in Eq. ( 11) and ( 13), respectively.Therefore, apart from the time evolutions, we are also concerned with the corresponding resolvent equation: with (y, n) ∈ R 3 S 2 , α > 0, i = 1, 2. Spectral decomposition of generator Q i (with a purely discrete spectrum) therefore also yields the spectral decomposition of both the (compact) operator Υ t and the (compact) resolvent operator (Q i − αI) −1 .It follows that Monte Carlo simulations of random walks with exponentially distributed traveling times can be used to approximate resolvent kernels.Furthermore, it can be shown that the resolvent operator occurs in Tikhonov regularization of the input [28].For the convection-diffusion case, the time-integrated process is practically even more useful than the time-dependent case.

Contributions of the Paper
In this paper, we provide exact expressions for the kernels corresponding to the generators Q 1 and Q 2 for diffusion and convection-diffusion, in terms of eigenfunctions of the operator, after applying a Fourier transform in the spatial coordinates.As such we solve the 3D version of Mumford's direction process [12], and the 3D version of the (2D) hypo-elliptic Brownian motion kernel, considered for orientation processing in image analysis in [3,10,6].The main results can be found in Theorems 2.3, 3.3 and 3.4 and in Corr.2.5.
A similar approach was used in [17] for the SE(2)-case, but for the SE(3)-case exact solutions are not known, to the best of our knowledge.For SE(2), solutions were expressed in [17] in terms of Mathieu functions.In SE(3), we encounter (generalized) spheroidal wave functions, that can be written as a series of associated Legendre functions.We use the eigenfunctions to give expressions for both the time-dependent and the time-integrated process, associated with a resolvent equation.Using the results for the exact solutions for the time-integrated case, we derive in Thm 3.6 the minimal amount of repeated convolutions of the resolvent kernels required to remove the singularities in the origin.The proof of this result is given in Appendix C. We also provide two numerical methods to approximate the kernels.First, we extend the numerical algorithm by August [1] from SE(2) to SE(3).This algorithm has the advantage that it is closely related to the exact solutions and can therefore be shown to yield convergence to the exact solutions.As a second approach we provide a novel approximation of the kernels, that provide correct symmetries in contrast to existing approximations provided in [28].This rough and local approximation is convenient as it allows for straightforward and fast implementations, but no convergence results can be obtained.A numerical approximation based on finite differences was presented in [27], and is therefore not considered here.A full comparison between exact solutions, stochastic process limits and kernel approximations for these processes on SE(2) has been done in [16].In the present paper, the focus will be on the derivation of the exact solutions and its connection to a numerical approximation method that generalizes the results [1,17] from SE(2) to SE (3).
Finally, in order to make the connection to earlier work on exact solutions of heat kernels in the SE(2)-case [10,17] via harmonic analysis and the Fourier transform on SE(2), we also derive equivalent representations of the kernels using the SE(3) Fourier transform.The details of this equivalent representation (and alternative roadmap to the solutions) is presented in Appendix D, where we follow an algebraic Fourier theoretic approach.We stress that the body of this article relies on classical, geometrical and functional analysis.However, the approach in Appendix D provides the reader further insight in specific choices in the geometric analysis in the main body of the article.

Outline of the Paper
The paper is structured as follows.In Sections 2 and 3 we derive exact expressions for the convolution kernels of the differential equations in (16), with the main results in Theorem 2.3 and Theorem 3.4.We emphasize that the roadmap and computations in these sections are very similar, only in the convergence proofs of the encountered series of eigenfunctions we need different theory for each case.In Section 4 we present a matrix representation of the evolution in a Fourier basis and provide an algorithm to numerically compute a truncation of the exact series solution.Section 5 shows the Gaussian kernel approximations.We summarize our findings and conclude in Section 6.The derivation of equivalent solutions for the main PDEs using the SE(3) Fourier transform can be found in Appendix D.

Derivation of the Exact Solutions for Hypo-Elliptic Diffusion on R 3 S 2
In this section we derive the exact solutions for the hypo-elliptic diffusion case.We first set out the formal procedure for finding these solutions, before we present the details for this particular case and the specific eigenfunctions that we encounter.The evolution process for hypo-elliptic diffusion on R 3 S 2 , i.e., with generator Q 1 as in (11), is written as follows: As in (11), (13), we use ∇ R 3 to indicate the gradient with respect to the spatial variables, and ∆ S 2 to denote the Laplace-Beltrami operator on the sphere.Parameters D 33 , D 44 > 0 influence the amount of spatial and angular regularization, respectively.We use both a subscript and superscript 1 throughout this section for operators that arise from this evolution, to distinguish these from the operators corresponding to the convection-diffusion that we encounter in Section 3. The domain of the generator Q 1 equals where we use H 2 to denote a Sobolev space, although in D(Q 1 ) both H 2 (R 3 ) and H 2 (S 2 ) are equipped with the usual L 2 -norm.For details on Sobolev spaces see e.g.[30,31].By linearity and left-invariance of the differential equation, the solution can be found by R 3 S 2 -convolution with the corresponding integrable kernel K 1 t : R 3 S 2 → R + : The specific choice for the rotation matrix R n does not matter, since the left-invariance of the PDE implies that , see [28].This α-invariance is important because of the equivalence relation described in (6).Our approach for finding the exact kernel K 1 t is inspired by the approach for the SE(2) case in [17].We first apply a Fourier transform with respect to the spatial variables: The hat is used to indicate that a function has been Fourier transformed.The PDE (20) in terms of Ŵ then becomes: We fix ω ∈ R 3 and we define the operator ) as follows: We use the subscript ω to explicitly indicate that the operator depends on the frequency vector ω, as will the eigenfunctions of this operator.When we write then the correspondence between the operator B and the generator Q 1 can be written as: The heat kernel K1 ) of the PDE in the Fourier domain should then satisfy: with δ ez the δ-distribution on S 2 at e z .We show later that the L 2 (S 2 )-normalized eigenfunctions of the operator B 1 ω form an orthonormal basis for L 2 (S 2 ) and that, similar to the enumeration of spherical harmonics, these functions are indexed with integers l and m, |m| ≤ l.For the eigenfunctions, that we denote with Φ ω l,m , we have: The kernel K1 t can then be written in terms of these eigenfunctions as The solution of the differential equation ( 24) in the Fourier domain is given by: where we rely on the following inner product convention on L 2 (S 2 ): Thereby, the solution of Eq. ( 20) is given by: This expression should coincide with the convolution expression in Eq. ( 22).The following lemma gives us a useful identity for the eigenfunctions Φ ω l,m , that allows us to connect the series expression for the solution of ( 20) with the convolution form in (22).
Proof.We can write from which the result follows.
From the lemma we conclude that the eigenvalues λ l,m ω only depend on the norm r = ||ω|| of the frequency ω, so from now on we write λ l,m r .Moreover, we have the following relation between eigenfunctions: We can combine the previous considerations and Lemma 2.1 to give two equivalent series expressions for the kernel K1 satisfying Eq. ( 28): Finally, the solution of (20) can then be written as Now that we have formal expressions for the solution of the hypo-elliptic diffusion PDE, we can focus on finding the eigenfunctions Φ ω l,m .

Frequency-Dependent Choice of Variables
So far we have not specified a choice of spherical variables for the functions to which the operator B 1 ω is applied.This is needed in order to derive expressions for the eigenfunctions Φ ω l,m .As in the case for spherical harmonics, we want to use separation of variables for each fixed spatial frequency ω.To be able to do this, we choose to parameterize the orientation of n using angles dependent on ω.This choice should be such that the variables can be separated in both the Laplace-Beltrami operator ∆ S 2 and in the multiplication operator (ω • n) 2 .
The Laplace-Beltrami operator on a Riemannian manifold with metric tensor G is defined as: The G ij are matrix elements of the inverse G −1 of G.The Laplace-Beltrami operator on the sphere ∆ S 2 is, when we take standard spherical coordinates, given by: with x = sin φ cos θ, y = sin φ sin θ and z = cos φ.However, in these coordinates the term (ω • n) 2 is not separable.
To this end, we choose spherical coordinates with respect to the (normalized) frequency r −1 ω, with r = ||ω||, and a second axis perpendicular to ω.The specific choice for the latter axis is not important, since in our PDE only the angle between ω and n plays a role.For convenience we take as second axis ω×ez ||ω×ez|| , and we let β, γ denote the angles of rotation about axes ω×ez ||ω×ez|| and r −1 ω, respectively.For r −1 ω = e z , β and γ are just the standard spherical coordinates.Every orientation n ∈ S 2 can now be written in the form see Fig. 2.
Lemma 2.2.The Laplace-Beltrami operator on S 2 with the choice of variables as in Eq. ( 40) is given by: The full operator B 1 ω as defined in (25) with this choice of coordinates takes the separable form Proof.The result follows from direct computation of the metric tensor G w.r.t. the coordinates in in (40).

Separation of Variables
Thanks to the choice of coordinates (40) and Lemma 2.2, we can apply the method of separation of variables to solve the diffusion equation (24).We first look for solutions of the PDE without initial conditions, and we take the solutions of the separable form T (t)Φ(β, γ): It follows that we get 1 T (t) where λ r is the separation constant.We get that T (t) = Ke λrt , K constant and λ r an eigenvalue of B 1 ω .For finding eigenfunctions Φ(β, γ) of B 1 ω , we assume that these functions can be written as We then find: The resulting equation for C(γ) can be solved straightforwardly.Taking into account the 2π-periodicity of C, we find that C(γ) is a multiple of e imγ for some m ∈ Z. Normalization then gives:

Spheroidal Wave Equation
The equation for B(β), with separation constant m as above, now becomes: With the substitution x = cos β, y(x) = B(β) (which is commonly done for equations of this type), we get: We introduce two more parameters to bring this equation in a standard form: We then find: This equation is known as the spheroidal wave equation (SWE) [32, 33, Eq. 30.2.1], for which the eigenvalues and eigenfunctions, commonly referred to as spheroidal eigenvalues and spheroidal wave functions, are known and can be found up to arbitrary accuracy.In the remainder of this article, we denote them with λl,m ρ and S l,m ρ , respectively, for l ∈ N 0 , m ∈ Z, |m| ≤ l.For explicit analytic representations, see (A.11) in Appendix A.
The (normalized) spheroidal wave functions and the spheroidal eigenvalues as a result of our computations are displayed in Figs. 3 and 4 for a selection of values of l and m.It can be seen that all eigenvalues are real and all eigenfunctions are real-valued and vary continuously with parameter ρ.In the next section, we show that the spheroidal wave functions form a complete orthonormal basis for L 2 ([−1, 1]), from which it follows that the eigenfunctions Φ ω l,m form a complete orthonormal basis for L 2 (S 2 ).

Sturm-Liouville Form
The spheroidal wave equation can be written in Sturm-Liouville form: For this we choose p , and we have weight function w(x) = 1.In this formulation, p(x) vanishes at the boundary of the interval, which makes our problem a singular Sturm-Liouville problem (on a finite interval).It is sufficient to require boundedness of the solution and its derivative at the boundary points to have nonnegative, distinct, simple eigenvalues and existence of a countable, complete orthonormal basis of eigenfunctions {y k } ∞ k=1 [34] for the spheroidal wave equation.Since the weight function w(x) = 1, orthogonality is understood in the following sense: For our choice of p(x), q(x) and w(x), with m fixed, we have y k (cos β) = S l,m ρ (cos β) (with k = l 2 + l + 1 + m) and corresponding eigenvalues λl,m ρ .

Main Theorem
From the above considerations, we can come to the main result of this section.
Theorem 2.3.The normalized eigenfunctions Φ ω l,m of the operator B 1 ω in (42) are given by:  51).Function C m is as in Eq. (47).The solution of the hypo-elliptic diffusion on R 3 S 2 Eq. ( 20) is given by: Proof.From the separation of variables approach above, and Eq. ( 53), it follows that the Φ ω l,m are indeed normalized eigenfunctions.From the fact that the resolvent of operator B 1 ω is compact, which follows from Sturm-Liouville theory, we obtain, using the spectral decomposition of compact self-adjoint operators, a complete orthonormal basis of eigenfunctions of L 2 (R 3 S 2 ).As dim(S 2 ) = 2, we can number the orthonormal basis with two indices l, m, similar to the spherical harmonics, that appear in the special case ω = 0.This allows us to use the expression in (30) and the result follows.
In the particular case of ω = 0, the operator B 1 ω reduces to B 1 0 = D 44 ∆ S 2 , for which the spherical harmonic functions are the eigenfunctions.For the spherical harmonics, we use the following convention: with the associated Legendre polynomials P m l as defined in (A.1), and In the following corollary we write the eigenfunctions Φ ω l,m in (54) directly in terms of the spherical harmonics.

Time-Integrated Processes and Resolvent Kernels
Under the assumption of exponentially distributed traveling times, the probability of finding a particle at a certain position and orientation can be expressed in terms of the resolvent (Q 1 − αI) −1 of the generator Q 1 , recall Eq. ( 11) and Section 1.5.It can be seen that for eigenfunctions of B 1 ω , we have: It follows that the resolvent kernel is given by Thereby the probability density of finding a random walker at a certain position y with orientation n, regardless of the traveling time, is given by: 2.7.From the Hypo-Elliptic Diffusion Kernel to the Elliptic Diffusion Kernel So far in this section we have restricted our diffusion by choosing only D 33 and D 44 = D 55 as nonzero entries in D, recall (2), motivated from the use of this process in applications.However, even for elliptic diffusion it is still possible to obtain exact solutions, with just a simple transformation from the hypo-elliptic case.We are still required to use only legal generators, in the sense discussed in Section 1.4.Furthermore, for all differentiable functions Ũ : Hereby, also the case of elliptic diffusion can be considered, i.e., D = diag(D 11 , D 11 , D 33 , D 44 , D 44 , 0) > 0, such that the generator QE of the evolution on SE(3) and the generator Q E on R 3 S 2 become: Recall Remark 1.1 for the relation between the group and quotient generators.The operator B ω on H 2 (S 2 ), obtained as before from applying a Fourier transform in the spatial variables, changes accordingly: It is now fairly straightforward to obtain from Theorem 2.3 the following corollary: Corollary 2.5.Let D 33 > D 11 > 0 and t > 0. Then the elliptic heat kernel K E t = e tQ E δ (0,ez) is given by with Φ E,ω l,m (n ω (β, γ)) = S Finally, the limit D 11 ↓ 0 can be interchanged with the sum in the series, since by application of the Weiertrass criterion (and the existence of a uniform bound on all eigenfunctions Φ E,ω l,m (n)) the series is uniformly converging for all t > 0.

Derivation of the Exact Solutions for Convection-Diffusion on R 3 S 2
In this section we consider the second central equation of the paper, Eq. ( 16) with i = 2, for the 3D direction process or convection-diffusion process.We focus here on the time-integrated process, as this has proven to be more useful in applications.We show how exact solutions for the resolvent kernel can be found.Similar to the approach in Section 2, we derive eigenfunctions for the corresponding evolution operator in the Fourier domain.However, this operator can no longer be transformed into the standard Sturm-Liouville form.We therefore use the framework of perturbations of self-adjoint operators [35,36] to prove important properties of the eigenvalues and to prove completeness of the eigenfunctions.
The convection-diffusion system that we consider is the following: In this section, we focus on the time-integrated, resolvent process (although the solution strategy is the same): Again we fix ω ∈ R 3 and the operator B 2 ω (superscript 2) corresponding to Q 2 now becomes: When we express n in spherical coordinates β, γ with respect to ω as was done in (40) and Fig. 2, the differential operator in β, γ becomes: Since the Laplace-Beltrami operator is symmetric and the multiplication operator has a purely imaginary symbol, the operator B 2 ω in this case is not symmetric and not self-adjoint, but does satisfy: Here we note that the domain of the closed unbounded operator B 2 ω is the Sobolev space H 2 (S 2 ), equipped with the L 2 (S 2 )-norm.It will turn out to be useful to regard B 2 ω as the sum of a self-adjoint operator and a bounded operator M : L 2 (S 2 ) → L 2 (S 2 ), that just applies a multiplication, (Mf )(n ω (β, γ)) = cos β • f (n ω (β, γ)): So in particular, as before, the operator B 2 0 = D 44 ∆ S 2 has the spherical harmonics as eigenfunctions.We denote the eigenfunctions of B 2 ω with Ψ ω l,m and the corresponding eigenvalues with λ l,m r , even though the eigenvalues are not the same as in Section 2. Again we assume that the eigenfunctions can be written as the product-form Ψ ω l,m (n ω (β, γ)) = B(β)C(γ), leading to two ordinary differential equations.The equation with variable γ is the same as before, with solutions e imγ .The equation for variable β can be written as i.e., as now with We define the differential operator B m ρ as with slight abuse of notation, since now M : L 2 ([0, π]) → L 2 ([0, π]).Then Eq. ( 72) can be rewritten as: We explicitly denote the dependence on ρ and the separation constant m, as we need it later in our spectral analysis of the operator.

The Generalized Spheroidal Wave Equation
After applying the transformation x = cos β, y(x) = B(β) to Eq. ( 72), we get: This equation now has the form of a specific case of the generalized spheroidal wave equation (GSWE) [37,33,Sec. 30.12].
Remark 3.1.In literature [37,38], the generalized spheroidal wave equation also appears in the following form: In this equation C i , ω, η are constants and the equation has singularities at x = 0 and x = x 0 .With an appropriate choice of constants, the GSWE can be brought to the form of Equation ( 76).However, this does require taking the limit ω → 0, such that 2ηω stays bounded and nonzero.According to [38], this type of limit is considered by both Whittaker and Ince.We are only interested in the solution on the interval [0, x 0 ].In [37,38], solutions are provided in various series expansions, but not all of them converge on [0, x 0 ] and not all of them allow for taking the limit ω → 0 as above.Analogous to SWE case, we derive solutions that are series of associated Legendre functions.This is different from the solutions in [37,38], but our series is also suited for the case ω = 0.
We refer to Appendix B for the derivation of the eigenvalues λl,m ρ of the GSWE and the corresponding eigenfunctions that we denote with GS l,m ρ .Here we just state that the eigenfunctions of B 2 ω are given by in which we used the same substitution x = cos β as before.The functions GS l,m ρ for certain ρ, l and m are shown in Fig. 5. Recall (73) for the relation λ l,m r = −D 44 λl,m ρ between the eigenvalues corresponding to Ψ ω l,m and GS l,m ρ , respectively.From property (69) the following can be derived: This implies that As a result, we see that if {Ψ ω l,m } is complete and it admits a reciprocal basis {Ψ l,m ω } such that (Ψ l,m ω , Ψ ω l ,m ) = δ l l δ m m , then for the reciprocal basis functions we have Ψ l,m ω = Ψ ω l,m .In the next section, this completeness is further discussed.

The Time-Integrated Process
The kernel for the time-integrated process in the spatial Fourier domain corresponds to: However, there are conditions on the convergence of this series expression.Since operator B 2 ω , in contrast to B 1 ω , is no longer self-adjoint, the standard Sturm-Liouville theory, that ensures completeness of the eigenfunctions with negative, real eigenvalues, cannot be applied.In the following we formulate a lemma, on the eigenvalues of B m ρ , and a theorem, on the eigenfunctions, that combined imply that the convergence holds almost everywhere.Only for particular radii ||ω|| = D 44 ρ m n , for some ρ m n in the frequency domain there is no convergence, but it can be shown that for any m this happens only on a countable set of {ρ m n } n∈N0 that has no accumulation point.As a result, the series in (81) converges almost everywhere in the Fourier domain, and thereby the inverse Fourier transform, similar to Eq. (55) in Theorem 2.3, is still well-defined.
In the next lemma we prove properties of the eigenvalues λ l,m r of the operator B 2 ω that are necessary to have convergence of the series in (81).For this we also need to consider the operator B m ρ for fixed m, recall the definition in (74).Lemma 3.2 (Eigenvalues of the operator B 2 ω ).We have the following properties for the eigenvalues of B 2 ω : 1.Let m ∈ Z, then there exists a ρ m * > 0 such that B m ρ has real eigenvalues.Moreover, there is at most a countable set {ρ m n } ∞ n=1 where two eigenvalues collide and branch into a complex conjugate pair of eigenvalues.
2. For all r ≥ 0 the real part of the eigenvalues of B 2 ω is negative.Proof.We prove the two points subsequently: 1.For ρ = 0, the operator B m 0 , recall (74), is self-adjoint, negative semi-definite and therefore all eigenvalues λl,m 0 = −l(l + 1), l ≥ |m| are real, negative and simple.From the spectral inclusion theorem [39] it follows for the spectrum σ(B m ρ ) that: The operator norm is ||iρM|| = ρ and for B m 0 the minimal distance between two eigenvalues is the distance between the two smallest eigenvalues, when l = |m| and l = |m|  76) that when y(x) is an eigenfunction for λ, that y(−x) is an eigenfunction for λ.It cannot happen that branching of eigenvalues occurs without two eigenvalues colliding, since the multiplicity of λ(ρ) depends continuously on ρ.Since we have shown that no eigenvalues can collide for ρ < ρ m * , we are guaranteed to have real eigenvalues in this case.Now according to [40], there exists a nonzero analytical function F (λ, ρ), such that the equation F (λ, ρ) = 0 defines the eigenvalues λ m i , m fixed, as functions of ρ.We define ρ m n to be those values for ρ, in increasing order, for which λ m i (ρ) = λ m j (ρ) for some i = j.Due to the analyticity of F , the set {ρ m n } ∞ n=0 is countable and cannot have an accumulation point [40].We specify the analytic function whose zeros provides for given m ∈ Z the values (ρ m n ) ∞ n=0 later (in (Eq.84)).2. To show that all eigenvalues have a negative real part, it is sufficient to show that the symmetric part of the operator B 2 ω is negative definite.Indeed for all f ∈ H 2 (S 2 ) (dense in L 2 (S 2 )), we have: The dependency of the eigenvalues on ρ is displayed for two different values of m in Fig. 6.In the figure, the points ρ m n for which B m ρ has two colliding eigenvalues are indicated with red dots.The points ρ m n are in fact zeros of the analytic function where the right hand side only depends on ρ = D 44 ω .Moreover, for the behavior of the eigenvalues we have In particular, for m fixed and ρ ≤ ρ m 0 , all eigenvalues are real, and hence (MΨ ω l,m , Ψ ω l,m ) = 0. We use Lemma 3.2 in the next theorem, which proves the solution of the time-integrated differential equation is unique.Proof.
1. (Existence of the resolvent) Let ω ∈ R 3 and α > 0 given, r = ||ω||.Injectivity of (B 2 ω − αI) follows from the fact that the resolvent operator is bounded from below.For f ∈ H 2 (S 2 ), we have: Note that all eigenvalues are real for sufficiently small ρ.When ρ increases, each time two eigenvalues collide and branch into two complex conjugate eigenvalue pairs.Comparing the left and right figure, we note that the higher m, the higher the values for ρ where this branching occurs.Moreover, we have that Im(λ ρD 44 ) ∼ O(ρD 44 ).

(Completeness)
We first consider the operator B m ρ .By direct computation it can be shown that the multiplication operator ρM is bounded, with ρ as a bound for the operator norm.For ρ = 0, B m 0 has simple eigenvalues.Thereby, according to Kato [35, Ch.V, Sect.5, Th. 4.15a] there exists a complete basis of generalized eigenfunctions.Then B m ρ is closed with compact resolvent (B m ρ − αI) −1 .In the case of ρ = ρ m n , we still need to show that the basis of generalized eigenfunctions correspond to the actual eigenfunctions Ψ ω l,m as computed above.For ω < D 44 ρ * this is clear.For ρ ≥ ρ * , it follows by analytic extension in ρ, which is exactly what happens when we write down the eigenfunctions as a series of Legendre functions as in (B.5), as this boils down to a Taylor series in ρ.Furthermore, the functions GS l,m ρ are uniformly bounded on [−1, 1] and thereby form a Riesz basis [36], which makes the reciprocal basis unique.The reciprocal basis has the property that In fact we have Ψ l,m ω = S −1 Ψ ω l,m , where S : L 2 (S 2 ) → L 2 (S 2 ) denotes the frame operator, given by From the properties (69) and ( 80) it follows that the reciprocal basis Ψ l,m ω of Ψ ω l,m is linearly proportional to the conjugate basis, which implies that the reciprocal basis {Ψ l,m ω } is complete.Therefore, for any f ∈ L 2 (S 2 ), there is the convergent series representation Hence for ρ = ρ m n , (B 2 ω − αI) −1 is diagonalizable with eigenfunctions Ψ ω l,m and eigenvalues 1/(λ l,m r − α) with strictly negative real part.

Main Theorem
The following theorem summarizes the result regarding eigenfunctions of the operator B 2 ω corresponding to the generator Q 2 : Theorem 3.4.The eigenfunctions Ψ ω l,m of the operator B 2 ω in (67) are given by: Here GS l,m ρ is the eigenfunction, given by (B.5), of the generalized spheroidal wave equation.For almost every ω ∈ R 3 , these eigenfunctions form a complete bi-orthogonal system.Therefore the solution of the convection-diffusion equation on R 3 S 2 is given by: where the series converges in L 2 (R 3 S 2 )-sense.With λ l,m r = −D 44 l(l + 1) + O(r) we denote the countable eigenvalues of B 2 ω , with ω = r = ρD 44 .The eigenvalues are disjoint for ρ = 0, and ρ = ρ m n .Remark 3.5.
• For all ω ∈ R 3 , including the cases ||ω|| = D 44 ρ m n , we can decompose the resolvent into a complete basis of generalized eigenfunctions.In fact the projection onto the generalized eigenspace E λl,m ρ is given by with J a Jordan curve enclosing only λl,m ρ in positive direction, where we recall definition (74).
• Now if ρ = ρ n m , we have one-dimensional eigenspaces, spanned by GS l,m ρ : In the case that ρ = ρ m n , we have instead: At ρ = ρ m n , the algebraic multiplicity of λl,m ρ is 2, whereas the geometric multiplicity is 1.

Time integration with Gamma-distributed traveling times
The kernel for the time-integrated process with exponentially distributed traveling times has a singularity at the origin (y, n) = (0, e z ).It is possible to derive a relation between the resolvent kernel and the process with Gamma-distributed traveling times, that does not suffer from this singularity.Thanks to the exact solution representation of the kernels for both processes, we obtain the following refinement and generalization of [29,Thm. 12].
Theorem 3.6.The kernel of the time-integrated diffusion (i = 1) and convection-diffusion (i = 2) process with the assumption of Γ-distributed evolution time T , i.e., T ∼ Γ(k, α), k ∈ N, with E(T ) = k α , is related to repeated convolutions of the resolvent kernel as follows: where Γ(t; k, α) is the pdf of T .For the case i = 1, this kernel does not have a singularity in the origin when k ≥ 2. For the case i = 2, this holds when k ≥ 4.
Proof.We refer to Appendix C for the proof.
When comparing kernels for varying k, while keeping the expected value of the traveling time T fixed, using a k higher than the bounds given in the theorem above gives better shaped kernels in practice, with more outward mass and dampened singularities.We use this idea in the visualization of a time-integrated contour completion kernel in Section 4, Fig. 8.

Matrix Representation of the Evolution and Resolvent in a Fourier Basis
In [17] a connection between the exact solutions and a numerical algorithm proposed by August and Zucker in [1,2] was established for the SE(2) = R 2 S 1 case.In this algorithm the Fourier transform on L 2 (R 2 ) and L 2 (S 1 ) are applied subsequently.Here we again establish such a connection between the exact solutions and such a numerical algorithm in the 3D case.
We have seen in the previous sections that the eigenfunctions Φ ω l,m and Ψ ω l,m are closely related to the spherical harmonic functions.Using the Fourier transform on L 2 (S 2 ), we naturally obtain an ordinary differential equation in terms of spherical harmonic coefficients.In Sections 4.1 and 4.2 we derive this ODE for diffusion and convection-diffusion, respectively.At the end of Section 4.1 we include two remarks: one that shows that in deriving the ODEs we encounter a matrix representation of the evolution operator B i ω and its resolvent, and one that makes the connection with the Fourier transform on SE(3) as presented in Appendix D. Finally, we show in Section 4.3 that with the procedure presented in this section it is straightforward to compute a numerical solution to the PDEs, since it only requires truncation of the order of the spherical harmonics.

Hypo-Elliptic Diffusion
Recall that after a Fourier transform in the spatial coordinates, we have the following system: Since the spherical harmonics form a basis for L 2 (S 2 ), we can expand Ŵ (ω, n, t) for fixed ω and t in the spherical harmonic basis.Instead of doing this with standard spherical harmonics, we use a specific type of reoriented spherical harmonics, that we define as follows: and m ∈ Z, l ∈ N 0 , |m| ≤ l.Recall Fig. 2 and Eq. ( 40).The rotation R r −1 ω in (99) is defined through its matrix This rotation maps e z onto r −1 ω and e y onto ω × e z /||ω × e z ||, such that for every β, γ we have that Now for fixed ω and t we develop Ŵ (ω, where we define Ŵ l,m (ω, t) (and similarly Û l,m (ω)): Our goal is the recursion in Eq. ( 112) for Ŵ l,m (ω, t), that we can use to obtain a solution for these coefficients.To this end we start by substituting (103) into our differential equation (98): We aim at rewriting the term cos 2 β Y l,m (β, γ), see (110) below.In [33] the following identity for Legendre functions is given: Using this identity twice, we get: ) However, this identity only holds for the Legendre functions where the normalization factor N l,m is not included in the definition.The identity needs to be adapted, which is done as follows: This can directly be applied to the term cos 2 βY l,m (β, γ).We define the tridiagonal matrix M m 1 as follows: with Then we can write the relation for spherical harmonics in the form: Because matrix M m 1 has only three nonzero diagonals, at most three terms appear in the sum on the right.Since we consider these matrices for fixed m, it is more natural to change the order of summation in Eq. (104).With (110) we can rewrite (104) as: Equating the coefficients in this equation, the functions Ŵ l,m can be expressed recursively according to: ) In other words, for m fixed, we need to solve the following system: with Λ m = diag(|m|(|m| + 1), (|m| + 1)(|m| + 2), . . .).The solution of this system, being a matrix-vector representation of (98), is given by Remark 4.1.There is a direct connection between the exact solution presented in Theorem 2.3/Corollary 2.4 and the solution found via Eq.( 114).In fact, in case of the exact solutions, the generator and thereby also the evolution operator in (114), are diagonalized by the solutions of Eq. (A.10) in Appendix A. To clarify this observation, We note that Eq. (A.10) is the eigenvector problem corresponding to the eigenfunction problem in Eq. ( 29), while restricting operator B ω 1 to the span of spherical harmonics Y l,m , l ≥ |m|, for m ∈ Z fixed.
Remark 4.2.Eq. ( 112) follows from Eq. ( 20) by application of the operator: that can be roughly formulated as There is a close connection between this operator and the Fourier transform and irreducible representations on SE(3), see Appendix D.

Resolvent of the Convection-Diffusion
The same idea of substituting directly a series of spherical harmonics into the equation (as done in the previous section for the pure diffusion case) can be used for the resolvent case of the convection-diffusion equation.After the Fourier transform in the spatial coordinates, this equation reads Again, we express Ŵ in terms of spherical harmonics as Substituting this into Eq.( 116), we obtain: We use the relation in Eq. ( 105) once, and we get where matrix M m 2 only has non-zero elements on the upper and lower diagonal: with ξ l,m and ν l,m as in (105) and N m as in (109).Using this matrix and proceeding analogously to the diffusion case, the solution to Eq. ( 118) is found by solving:

Numerical Implementation of the Discrete Spherical Transform
When using the above procedure to compute the kernel, there are two places where the numerics differ from the exact solution: the series of spherical harmonics is truncated and the Fourier transform is carried out discretely.We introduce a parameter l max to indicate the maximal order of spherical harmonics that is taken into account.Furthermore, we take discrete values for ω on an equidistant cubic grid, say ω ijk , such that for each ω ijk ∈ R 3 the component ŵm of the solution requires, for the pure-diffusion case, solving the ODE: where M m 1,lmax , Λ m lmax ∈ R (lmax−|m|+1)×(lmax−|m|+1) .In the convection-diffusion, resolvent case, it comes down to solving: for all ω = ω ijk on an equidistant grid: where N denotes the number of samples.We use a discrete centered inverse Fourier transform to go back to the R 3 S 2 domain.
Remark 4.3.In our implementation of the algorithm for Eq. ( 122), we work with the basis Y l,m 0 = Y l,m .We use Wigner D-matrices to bring the coefficients with respect this basis to coefficients with respect to the rotated spherical harmonics Y l,m ω and vice versa.This is done for all m simultaneously to transform the coefficients of the initial condition.Then Eq. ( 114) is solved for each m, and the solution is transformed back, again for all m at once.

Analytical Approximations Using Logarithm on SE(3)
For performing the enhancement of an image U : R 3 S 2 via the PDE framework as before, it is beneficial to have a good approximation of the convolution kernel.We have derived the exact solutions and related them to Fourier based algorithms, but they are not feasible for fast computation of the enhancement of the image.Numerical approximations of the kernels exist [28], for which improvements have been proposed in our conference paper [42].A fast implementation for neuroimaging purposes was made by [26].In the following lemma we show three properties of the solution operator and the kernel of the pure diffusion case on the group SE(3), from which we derive later a number of properties for the kernel on the group quotient )) be a linear and legal kernel operator, and assume it maps L 2 (SE(3)) into L ∞ (SE(3)).Then we have:  k(g, q) Ũ (q) dq = ( Ǩ * SE(3) Ũ )(g), with k(g, q) = k(e, (g −1 q)) =: Ǩ(q −1 g) and with K(g) := k(g, e) = Ǩ(g −1 ).

preservation of correspondences:
The proof of this lemma was given in [42].Now we consider the special case where K represents the solution kernel K 1 t of the diffusion process with generator Q 1 .Corollary 5.2.The exact kernels K1 t : SE(3) → R + and K 1 t : R 3 S 2 → R + satisfy the following symmetries Now the Gaussian kernel approximation, based on the logarithm on SE(3) and the theory of coercive weighted operators on Lie groups [8], presented in earlier work [29,Ch:8,Thm.11]does not satisfy this symmetry.Next we will improve it by a new practical analytic approximation which does satisfy the property.

Comparison with Earlier Kernel Approximations for the Diffusion Kernel
We first present two existing approximations for K 1 t and K1 t , that do not automatically carry over the properties we have shown in the previous section.A possible approximation kernel for K 1 t is based on a direct product of two SE(2)-kernels, see [43, Eq.( 10), Eq.( 11)], to which we refer as K 1,App 1 t , with App being short for approximation.This approximation is easy to use since it is defined in terms of the Euler angles of the corresponding orientations.However, the symmetries described before are not preserved by K 1,App 1 t and errors tend to be larger when D 44 and t increase.Therefore we move to an approximation for the SE(3)-kernel K1 t , for which we show that it can be adapted such that it has the important symmetries.We need the logarithm on SE(3) for this approximation, so first consider an exponential curve in SE(3), given by γ(t) = g 0 exp t and we can relate to this the vector of coefficients c(g) = (c (1) (g), c (2) (g)) = (c 1 (g), • • • , c 6 (g)) T .We define matrix Ω as follows: We can write R in terms of Euler angles, R = R ez,γ R ey,β R ez,α .Let the matrix Ω γ,β,α be such that exp x,γ,β,α are given by the following equation: where q γ,β,α is the (Euclidean) norm of c (2) .Then an approximation for the kernel K1 t is, up to a rescaling of time t → ct with c ≥ 1 due to an equivalence of two norms 2 [8, Prop.6.1],given by [28]: with the smoothed variant of the weighted modulus, [29, Eq.78,79], given by where we included a scaling with ξ > 0 in the commutator term.In view of computations in [6, ch.5.2] we set ξ = 16.First numerical investigations indicate this is needed to get Gaussian kernels close to the exact kernels (compare the result in Figure 9 to the exact result in Figure 7).Now the difficulty lies in the fact that the function U is defined on the quotient for elements (y, R n ), where the choice for α in the rotation matrix R n (mapping e z onto n ∈ S 2 ) is not of importance.However, the logarithm is only well-defined on the group SE(3), not on the quotient R 3 S 2 , and explicitly depends on the choice of α.It is therefore not straightforward to use this approximation kernel such that the invariance properties in Corollary 5.2 are preserved.In view of Corollary 5.2, we need both left-invariance and rightinvariance for K1 t (g) under the action of the subgroup H.As right-invariance is naturally included via K1 t (y, R) = K 1 t (y, Re z ), left-invariance is equivalent to inversion invariance.In previous work the choice of α = 0 is taken, giving rise to the approximation However, this section is not invariant under inversion.In contrast, we propose to take the section α = −γ, which is invariant under inversion (since (R ez,γ R ey,β R ez,−γ ) −1 = R ez,γ R ey,−β R ez,−γ ).Moreover, this choice for estimating the kernels in the group is natural, as it provides the weakest upper bound kernel since by direct computation one has α = −γ ⇒ c 6 = 0. Finally, this choice indeed provides us the correct symmetry for the Gaussian approximation of K 1 t (y, n) as stated in the following theorem, that is proven in [42].
Theorem 5.4.When the approximate kernel K 1,new t on the quotient is related to the approximate kernel on the group K1,log i.e. we make the choice α = −γ, we have the desired α-left-invariant property and the symmetry property

Comparison Logarithmic Approximation to Exact Solutions
Qualitatively, the approximation kernel can be used effectively in applications.Here D 44 > 0 tunes the angular diffusion, t > 0 regulates the size and D 33 tunes the spatial diffusion.However, from the general theory [8] we can only expect locally good approximations.In comparison to the exact kernels, these approximation kernels do not accurately solve the PDEs for these parameters.This is mainly due to the fact that the non-commutative structure of SE( 3) is only encoded in the logarithm (128), which is not sufficient.On top of this, we recall the required rescaling of time, which requires a smaller time in the Gaussian approximations than in the exact solutions.A visualization of the approximation kernel can be found in Fig. 9 which is qualitatively comparable to the exact kernel depicted in Fig. 7.

Conclusion
We have provided the explicit solutions for the hypo-elliptic diffusion process for (restricted) Brownian motion on SE(3) and for the convection-diffusion process on SE (3), that is an extension of Mumford's 3Ddirection process [12].The solutions were derived by applying a Fourier transform in the spatial variables and a particular (frequency dependent) choice of coordinates, yielding for both processes a separable second order differential operator.Using the spectral decomposition of this operator, we have obtained a series expression for the solution kernel of the evolution equation and the kernel for the time-integrated process, related to the resolvent of this operator.
In the diffusion case, the eigenfunctions encountered in the spectral decomposition are similar to spherical harmonics, but require spheroidal wave functions instead of Legendre functions.Convergence of the series expressions in terms of the eigenfunctions was shown using Sturm-Liouville theory.The final expression for the exact solution is given in Theorem 2.3.
For the convection-diffusion case, generalized spheroidal wave functions are needed, that we derive applying a non-standard expansion in Legendre functions.In this case, in contrast to the diffusion case, the considered operator is no longer self-adjoint and as a result, standard Sturm-Liouville cannot be used for proving convergence of the series.Instead we prove completeness of the eigenfunctions and specific properties of the spectrum using perturbation theory of linear (self-adjoint) operators.The exact solution for the resolvent case is given in Theorem 3.4.
We have also established a numerical algorithm on SE(3) that generalizes the algorithm in [1, 2] for Mumford's 2D direction process to the 3D case.Using this method, approximate solutions can be found by rewriting the equation into an ODE in terms of expansion coefficients w.r.t. a rotated spherical harmonic basis.We connect this numerical algorithm to the exact solutions, showing that it in essence just diagonalizes (the matrix of) the resolvent operator.Truncation of the exact series representation of the kernels yield the same result.Results of the truncated kernels are shown in Figs.7 and 8.
Building on previous work on approximations of the kernels, we have included Gaussian estimates for the hypo-elliptic diffusion kernel and we have compared them qualitatively to the exact solution.The new approximation has the same symmetries as the exact kernel and by making use of a scaling factor, the approximation is good for reasonable parameter settings in practice.Although the parameters have to be chosen differently, the result is qualitatively comparable to the kernel in Fig. 7.
Finally, we have algebraically derived the same exact solutions via harmonic analysis, using the Fourier transform on SE(3), see Theorem Appendix D.3.
Future work will include extensive analysis of the Monte-Carlo simulations, briefly described in Appendix E. Furthermore, we will put connections of hypo-elliptic diffusion (and Brownian bridges) on SE(3) and recently derived sub-Riemannian geodesics in SE(3) [44].
where F G maps L 2 (G) unitarily onto the direct integral space ⊕ Ĝ B(H σ ) dν(σ) with ν the Plancherel measure on Ĝ.For more details see [45,46,47].One has the inversion formula: (D.2) In our Lie group case of SE(3) we identify all irreducible representations σ p,s having non-zero dual measure with the pair (p, s) ∈ R + × Z.This identification is commonly applied, see e.g.[9].All unitary irreducible representations (UIR's) of G, up to equivalence, with non-zero Plancherel measure are given by [48,46]: For an explicit formula for the modified spherical harmonics, see [9], where they are denoted with h l m,s .Moreover, we have inversion formula ([9, Eq.10.46]):Now in our position and orientation analysis we must constrain ourselves to a special type of functions f = Ũ ∈ L 2 (G), namely the ones that have the property that f (x, R) = Ũ (x, R) = U (x, Re z ).Then it follows by [9, eq.10.35] and by (D.7) that the only non-zero matrix elements (D.6) have the property that both m = 0 and m = 0.This can also be observed in the analytical examples in [9, ch:10.10].This reduces the five-fold sum in (D.5) to a three-fold sum, and the modified spherical harmonics become standard spherical harmonics.
Next we introduce some notation for generators of group representations of SE(3).

1. 3 .
The Space of Coupled Positions and Orientations R 3 S 2 Embedded in SE(3)

Figure 1 :
Figure 1: Various visualizations of the diffusion process (for contour enhancement, left, generated by Q 1 ) and the convectiondiffusion process (for contour completion, right, generated by Q 2 ).For the SE(2) case, both the projections of random paths on R 2 are shown, as well as an isocontour in SE(2) of the limiting distribution.For the SE(3) case, we show spatial R 3 -projections of random paths in SE(3), as given in Appendix E, and visualize the limiting distribution as a glyph field, see Remark 1.2.

Figure 6 :
Figure 6: Plot of the real and imaginary parts of the first 6 eigenvalues of B m ρ for m = 0 (left) and m = 3 (right).Note that all eigenvalues are real for sufficiently small ρ.When ρ increases, each time two eigenvalues collide and branch into two complex conjugate eigenvalue pairs.Comparing the left and right figure, we note that the higher m, the higher the values for ρ where this branching occurs.Moreover, we have that Im(λ ρD 44 ) ∼ O(ρD 44 ).

Figure 7 :
Figure 7: Glyph field visualization (as explained in Section 1.5) of the kernel K 1 t=2 (ω, n), with a higher resolution on the right.For this kernel, D 33 = 1, D 44 = 0.1.

Figure 8 :
Figure 8: Glyph field visualization of the time dependent kernel K 2 t=1 for the convection-diffusion case, with D 44 = 0.5 (left) and the time-integrated kernel, where we used a Gamma-distribution with k = 4 and α = .25.

)Remark 5 . 3 .
Proof.The symmetry (126) is due to Lemma 5.1.The second symmetry (127) is due to reflectional invariance (A 3 , A 4 , A 5 ) → (−A 3 , −A 4 , −A 5 ) in the diffusion(20) and reflection on the Lie algebra corresponds to inversion on the group.Note that (127) in terms of k would be: k(y, n, 0, e z ) = k(0, e z , y, n), by the relation k

Figure 9 :
Figure 9: Glyph field visualization of the logarithmic approximation kernel K 1,new t , with D 33 = 1., D 44 = 0.24, t = 0.7, ξ = 16.Although the parameters have to be chosen differently, the result is qualitatively comparable to the kernel in Fig.7.