Design and Processing of Invertible Orientation Scores of 3D Images

The enhancement and detection of elongated structures in noisy image data are relevant for many biomedical imaging applications. To handle complex crossing structures in 2D images, 2D orientation scores \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$U: {\mathbb {R}} ^ 2\times S ^ 1 \rightarrow {\mathbb {C}}$$\end{document}U:R2×S1→C were introduced, which already showed their use in a variety of applications. Here we extend this work to 3D orientation scores \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$U: {\mathbb {R}} ^ 3 \times S ^ 2\rightarrow {\mathbb {C}}$$\end{document}U:R3×S2→C. First, we construct the orientation score from a given dataset, which is achieved by an invertible coherent state type of transform. For this transformation we introduce 3D versions of the 2D cake wavelets, which are complex wavelets that can simultaneously detect oriented structures and oriented edges. Here we introduce two types of cake wavelets: the first uses a discrete Fourier transform, and the second is designed in the 3D generalized Zernike basis, allowing us to calculate analytical expressions for the spatial filters. Second, we propose a nonlinear diffusion flow on the 3D roto-translation group: crossing-preserving coherence-enhancing diffusion via orientation scores (CEDOS). Finally, we show two applications of the orientation score transformation. In the first application we apply our CEDOS algorithm to real medical image data. In the second one we develop a new tubularity measure using 3D orientation scores and apply the tubularity measure to both artificial and real medical data.


Introduction
The enhancement and detection of elongated structures are important in many biomedical image analysis applications. These tasks become problematic when multiple elongated structures cross or touch each other in the data. In these cases it is useful to work with multi-orientation representations of image data. Such multi-orientation representations can be made using various techniques, such as invertible orientation scores (which is obtained via a coherent state transform) [3,5,10,30,36,42], continuous wavelet transforms [10,28,30,64], orientation lifts [13,71], or orientation channel representations [35]. Here we focus on constructing an invertible orientation score. In order to separate the crossing or touching structures (Fig. 1), we extend the domain of the data to include orientation. This is achieved by correlating our 3D data f : R 3 → R with a set of oriented filters to construct a 3D orientation score U : R 3 × S 2 → C. As the transformation between image and orientation score is stable, due to our design of anisotropic wavelets, we can robustly relate operators on the score to operators on images. To take advantage of the multi-orientation decomposition, we consider operators on orientation scores and process our data via orientation scores (Fig. 2).
Regarding the invertibility of the transform from image to orientation score, we note that in comparison to continuous wavelet transforms (see, e.g., [4,48,50,51] and many others) on the group of 3D rotations, translations and scalings, we use all scales simultaneously and exclude the scaling group from the wavelet transform and its adjoint, yielding a coherent state type of transform [2]. This makes it harder to design appropriate wavelets, but has the computational advantage of only needing one all-scale transformation.
The 2D orientation scores have already showed their use in a variety of applications. In [37,64] the orientation scores were used to perform crossing-preserving coherenceenhancing diffusions. These diffusions greatly reduce the noise in the data, while preserving the elongated crossing structures. Next to these generic enhancement techniques, the orientation scores also showed their use in retinal vessel tracking [8,10,19], in vessel segmentation [70] and biomarker analysis [11,60], where they were used to better handle crossing vessels. Here we aim to extend such techniques to 3D data.
To perform detection and enhancement operators on the orientation score, we first need to transform a given grayscale image or 3D dataset to an orientation score in an invertible way. In previous works various wavelets were introduced to perform a 2D orientation score transform. Some of these wavelets did not allow for an invertible transformation (e.g., Gabor wavelets [48]). A wavelet that allows an invertible transformation was proposed by Kalitzin et al. [46]. A generalization of these wavelets was found by Duits [25] who derived formal unitarity results and expressed the wavelets in a basis of eigenfunctions of the harmonic oscillator. This type of wavelet was also extended to 3D. This wavelet, however, has some unwanted properties such as poor spatial localization (oscillations) and the fact that the maximum of the wavelet does not lie at its center. In [25] a class of 2D cake wavelets was introduced that have a cake-piece-shaped form in the Fourier domain. The cake wavelets simultaneously detect oriented structures and oriented edges by constructing a complex orientation score U : R 2 × S 1 → C. Because the family of rotated cake wavelets cover the full Fourier spectrum, invertibility is guaranteed.
In this article we propose a 3D version of the cake wavelets. A preliminary attempt to generalize these filters was done in [44], where the plate detectors in [25] were extended to complex-valued cake wavelets with a line detector in the real part. Compared to these previous works, the filters in this work are now exact until sampling in the Fourier domain. For these filters we have no analytical description in the spatial domain as filters are obtained via a discrete inverse Fourier transform (DFT). Therefore, we additionally consider expressing filters of this type in the 3D generalized Zernike basis. For this basis we have analytical expressions for the inverse Fourier transform, allowing us to find analytical expressions for the spatial filters. This has the additional advantage that they allow for an implementation with steerable filters, see App. A. These analytical expressions are then used to validate the filters obtained using the DFT method. We also show applications of these filters and the orientation score transformation in 3D vessel analysis. That is, we present crossing-preserving diffusions for denoising 3D rotational Xray of blood vessels in the abdomen and we present a tubularity measure via orientation scores and features based on this tubularity measure, which we apply to cone beam CT data of the brain. An overview of the applications is presented in Fig. 3. Regarding our nonlinear diffusions of 3D rotational Xray images via invertible orientation scores, we observe that complex geometric structures in the vasculature (involving multiple orientations) are better preserved than with nonlinear diffusion filtering directly in the image domain. This is in line with previous findings for nonlinear diffusion filtering of 2D images [37] and related works [54,63,66] that rely on other more specific orientation decompositions.
For the sake of general readability, we avoid Lie group theoretical notations, until Sect. 6.1 where it is strictly needed. Let us nevertheless mention that our work fits in a larger Lie group theoretical framework, see, for example, [3,7,25,26,39] that has many applications in image processing. Besides the special cases of the Heisenberg group [6,31,57], the 2D Euclidean motion group [5,9,21,22,30,53,59], the similitude group [4,56,64] and the 3D rotation group [52], we now consider invertible orientation scores on the 3D Euclidean motion group (in which the coupled space of positions and orientations is embedded). The invertible orientation scores relate to coherent states from physics [42] for n = 3, with a specific (semi-)analytic design for our image processing purposes. Top: We reduce noise in real medical image data (3D rotational Xray) of the abdomen containing renal arteries by applying diffusions via 3D ori-entation scores. Bottom: Geometrical features can be directly extracted from our tubularity measure via 3D orientation scores. We apply this method to cone beam CT data of the brain

Contributions of the Article
The main contributions per section of the article are: -In Sect. 2 we give an overview of the discrete and continuous 3D orientation score transformation. Additionally, we present a transformation which is split in low and high frequencies and quantify the stability of the transformation in Lemma 1. item In Sect. 3 we present the cake wavelets obtained using the DFT method and give an efficient implementation using spherical harmonics which is summarized in Result 1. Furthermore, we analyze the stability of the transformation for these filters (Proposition 1). -In Sect. 4 we present the analytical versions of the cake wavelets obtained by expansion in the Zernike basis followed by a continuous Fourier transform which is summarized in Result 2. -In Sect. 5 we compare the two types of filters and show that the DFT filters approximate their analytical counterparts well. -In Sect. 6 we show two applications of the orientation score transformation: 1. We propose an extension of coherence-enhancing diffusion via our invertible orientation scores of 3D images. Compared to the original idea of coherenceenhancing diffusion acting directly on image data [16,17,69], there is the advantage of preserving crossings. Here we applied our method to real medical image data (3D rotational Xray) of the abdomen containing renal arteries. We show quantitatively that our method effectively reduces noise (quantified using contrast-to-noise ratios (CNR)) while preserving the complex vessel geometry and the vessel widths. Furthermore, qualitative assessment indicates that our denoising method is very useful as preprocessing for 3D visualization (volume rendering). 2. We develop a new tubularity measure in 3D orientation score data. This extends previous work on tubularity measures using 2D orientation scores [19] [9, Ch. 12] to 3D data. We show qualitatively that our measure gives sharp responses at vessel centerlines and show its use for radius extraction and complex vessel segmentation in cone beam CT data of the brain.

Outline of the Article
To summarize, we give a global outline of the article: First, we discuss the theory of invertible orientation score transforms in Sect. 2. Then we construct 3D cake wavelets and give a new efficient implementation using spherical harmonics in Sect. 3, followed by their analytical counterpart expressed in the generalized Zernike basis in Sect. 4. Then we compare the two types of filters and validate the invertibility of the orientation score transformation in Sect. 5. Finally, we address two application areas for 3D orientation scores in Sect. 6 and show results and practical benefits for both of them. In the first application (Sect. 6.1), we present a natural extension of the crossing-preserving coherence-enhancing diffusion on invertible orientation scores (CEDOS) [37] to the 3D setting.
In the second application (Sect. 6.2) we use the orientation score to define a tubularity measure and show experiments applying the tubularity measure to both synthetic data and brain CT data. Both application sections start with a treatment of related methods.

Continuous Orientation Score Transform
Throughout this article we use the following definition of the Fourier transform on R 3 : An invertible orientation score W ψ [ f ] : R 3 × S 2 → C is constructed from a given ball-limited 3D dataset with ball B = {ω ∈ R 3 | ω < } of radius > 0, by correlation with an anisotropic kernel This is illustrated in Fig. 4.
is a wavelet aligned with and rotationally symmetric around the z-axis, and ψ n ∈ L 2 (R 3 ) the rotated wavelet aligned with n given by Here R n ∈ SO(3) is any 3D rotation which rotates the zaxis onto n where the specific choice of rotation does not matter because of the rotational symmetry of ψ. The overline denotes a complex conjugate. The exact reconstruction formula [25] for this transformation is given by withψ n (x) = ψ n (−x). The function M ψ : R 3 → R + is given by and vanishes at ∞, where the circumflex (ˆ) again denotes Fourier transformation. Due to our restriction to ball-limited data (2), this does not cause problems in reconstruction (5). The function M ψ quantifies the stability of the inverse transformation [25], since M ψ (ω) specifies how well frequency component ω is preserved by the cascade of construction and reconstruction when M −1 ψ would not be included in Eq. (5). An exact reconstruction is possible as long as In practice it is best to aim for M ψ (ω) ≈ 1, in view of the condition number of transformation W ψ : where M and δ are assumed to be tight bounds in (7). In the codomain spatial frequencies are again limited to the ball: Construction of a 3D orientation score. Top: The data f are correlated with an oriented filter ψ ex to detect structures aligned with the filter orientation e x . Bottom left: This is repeated for a discrete set of filters with different orientations. Bottom right: The collection of 3D datasets constructed by correlation with the different filters is an orientation score and is visualized by placing a 3D dataset on a number of orientations Also, in the case we have M ψ (ω) = 1 for ω ∈ B we have L 2 -norm preservation and reconstruction Eq. (5) simplifies to We can further simplify the reconstruction for wavelets for which the following additional property holds: In that case the reconstruction formula is approximately an integration over orientations only: For the reconstruction by integration over angles only we can analyze the stability via the condition number of the mapping that maps an image f ∈ L 2 (R 3 ) to an orientation integrated score Its condition number is given by cond( Nψ (ω) . In practice, we always use this last reconstruction because practical experiments show that performing an additional convolution with the wavelet as done in reconstruction (5) after processing the score can lead to artifacts. It is, however, important to also consider the reconstruction (5) and M ψ because it is used to quantify the stability and norm preservation of the transformation from image to orientation score.
The fact that we use reconstruction by integration while still taking into account norm preservation by controlling M ψ leads to restrictions on our wavelets which are captured in the following definition: Definition 1 (Proper wavelet) Let us set a priori bounds 1 δ, M > 0, 0 < ε 1. Furthermore, let be an a priori maximum frequency of our ball-limited image data. Then, a wavelet ψ ∈ L 2 (R 3 where R e z ,α ∈ SO(3) is the 3D rotation around axis e z over angle α.
If, moreover, one has then we speak of a proper wavelet with fast reconstruction property, cf. (13).

Remark 1
The 1st condition (symmetry around the z-axis) allows for an appropriate definition of an orientation score rather than a rotation score. The 2nd condition ensures invertibility and stability of the (inverse) orientation score transform. The 3rd condition allows us to use the approximate reconstruction by integration over angles only.

Remark 2
Because of finite sampling in practice, the constraint to ball-limited functions is reasonable. The constraint is not a necessary one when one relies on distributional transforms [10, App. B], but we avoid such technicalities here. 1 In practice we choose the default values δ = 1 8 and M = 1.1 and ε = 0.01 and note it is actually the ratio M δ that determines the condition number. It is just a convenient choice to set the upper bound close to 1.

Low-Frequency Components
In practice we are not interested in the zero and lowest frequency components since they represent average value and global variations which appear at scales much larger than the structures of interest. We need, however, to store these data for reconstruction. Therefore, we perform an additional splitting of our wavelets into two parts with Gaussian window in the Fourier domain given bŷ After splitting, ψ 0 contains the average and low-frequency components and ψ 1 the higher frequencies relevant for further processing. In continuous wavelet theory it is also common to separately store very low-frequency components separately, see, e.g., [51,65]. In this case we construct two scores: one for the high-frequency components and one for the low-frequency components Here we again have ψ i,n (x) = ψ i (R T n x), as in Eq. (4). The vector transformation is then defined as For this transformation we have the exact reconstruction formula with Again, M ψ quantifies the stability of the transformation. The next lemma shows us that the stability of the transformation is maintained after performing the additional splitting.
such that Eq. (7) holds, δ = min ω∈B M ψ (ω) and M = max ω∈B M ψ (ω). Then the condition number of W ψ : The condition number of W ψ : obtained from W ψ by performing an additional splitting in low-and high-frequency components is given by thereby guaranteeing that stability is maintained after performing the splitting.
Proof First, we find the condition number of W ψ : For the first factor in Eq. (27) we find where in the last step the supremum is attained by a sequence of images whose Fourier transform concentrates around the maximum of M ψ . Similarly, we get max ω∈B M ψ (ω) for the second factor in Eq. (27). Then we obtain Similarly the condition number of W ψ is given by Next we express M ψ in M ψ of the original wavelet as So it remains to quantify I (ω). For a wavelet splitting according to (18) we have Hence And since 1 2 ≤ 1 − 2x(1 − x)) ≤ 1 for 0 ≤ x ≤ 1 we have for M ψ satisfying (7) the following bounds on M ψ : thereby guaranteeing stability after the splitting (18). As we cannot guarantee that δ/2 and M are tight bounds in Eq. (34), combining it with (29) will only give us an upper bound for the condition number in Eq. (26).
For this vector transformation we can also use the approximate reconstruction by integration (for N ψ ≈ 1) over orientations. Thus we have As said we are only interested in processing of W ψ 1 [ f ] and not in processing of W ψ 0 [ f ], and so we directly calculate For a design with N ψ (ω) = 1 for all ω ∈ B , we havê φ 0 =Ĝ s ρ and so Then Eq. (35) becomes

Discrete Orientation Score Transform
In the previous section, we considered a continuous orientation score transformation. In practice, we have only a finite number of orientations. To determine this discrete set of orientations we uniformly sample the sphere using a simple 2 electrostatic repulsion model [18].
The exact reconstruction formula is in the discrete setting given by with Δ i the discrete spherical area measure ( N o i=1 Δ i = 4π ) which for reasonably uniform spherical sampling can be approximated by For spherical samplings that are constructed by triangularization (such as tessellations of the icosahedron), one can rely on surface areas of spherical triangles to compute Δ i more accurately. See for example [29,Eq. 83]). Again, an exact reconstruction is possible iff 0 < δ ≤ M d ψ (ω) ≤ M < ∞ and we have norm preservation when M d ψ =1. Again, for the wavelets for which the image reconstruction can be simplified by a summation over orientations: For this reconstruction by summation we can analyze the stability via the condition number of the mapping that maps an image f ∈ L 2 (R 3 ) to an orientation integrated score This transformation has condition number cond(A d . Similar to Definition 1 for the continuous case, the reconstruction properties of a set of filters is captured in the following definition: Definition 2 (Proper wavelet set) Let us again set a priori bounds δ, M > 0, 0 < ε 1. Let be an a priori maximum frequency of our ball-limited image data. Then, a set of wavelets N o ), constructed as rotated versions of ψ is called a proper wavelet set if where R e z ,α ∈ SO(3) is a 3D rotation around axis e z over angle α.
If, moreover, one has then we speak of a proper wavelet with fast reconstruction property, cf. (43).

Low-Frequency Components
For the discrete transformation we will also perform a splitting in low-and high-frequency components as explained in Sect. 2.1.1. The reconstruction formula by summation in Eq. (43) is now given by

Steerable Orientation Score Transform
Throughout this article we shall rely on a spherical harmonic decomposition of the angular part of proper wavelets in the spatial and the Fourier domain. This has the benefit that one can obtain steerable [36,38,61] implementations of orientation scores, where rotations of the wavelets are obtained via linear combination of the basis functions. As such, computations are exact and no interpolation (because of rotations) takes place. Details are provided in Appendix A.

Wavelet Design Using a DFT
A class of 2D cake wavelets, see [10,27,37], was used for the 2D orientation score transformation. We now present 3D versions of these cake wavelets.
Because by definition the wavelet ψ has rotational symmetry around the z-axis we have h(ϑ, ϕ) = h(ϑ). 5. The kernel should be localized in the spatial domain, since we want to pick up local oriented structures. 6. The real part of the kernel should detect oriented structures, and the imaginary part should detect oriented edges. The constructed orientation score is therefore a complex-valued orientation score, as the wavelets are complex-valued. For an intuitive preview, see the boxes in Fig. 7.

Construction of Line and Edge Detectors
We now discuss the procedure used to make 3D cake wavelets before splitting in low and high frequencies according to To satisfy requirement 2 we should choose radial function g(ρ) = 1 for ρ ∈ [0, ]. In practice, this function should go to 0 when ρ tends to the Nyquist frequency ρ N to avoid long spatial oscillations. For the radial function g(ρ) we use, with erf(z) = 2 √ π z 0 e −x 2 dx, which is approximately equal to one for largest part of the domain and then smoothly goes to 0 when approaching the Nyquist frequency. We fix the inflection point of this function g and set the fundamental parameter for ball limitedness to with 0 γ < 1. The steepness of the decay when approaching ρ N is controlled by the parameter σ er f which we by default set to σ er f = 1 3 (ρ N − ). The additional splitting in low and high frequencies according to Sect. 2.1.1 effectively causes a splitting of the radial function, see Fig. 5.
In practice the frequencies in our data are limited by the Nyquist frequency (we have ≈ ρ N ), and because radial function g causes M d ψ to become really small close to the Nyquist frequency, reconstruction Eq. (40) becomes unstable. We solve this by using approximate reconstruction Eq. (13). Alternatively, one could replace M d ψ by max(M d ψ , ) in Eq. (5), with small. Both make the reconstruction stable at the cost of not completely reconstructing We now need to find an appropriate angular part h for the cake wavelets. First, we specify an orientation distribution A : S 2 → R + , which determines what orientations the wavelet should measure. To satisfy requirement 3 this function should be a localized spherical window, for which we propose the spherical diffusion kernel [20]: with s o > 0 and n(ϑ, ϕ) = (sin ϑ cos ϕ, sin ϑ sin ϕ, cos ϑ).
The parameter s o determines the trade-off between requirements 2 and 3 listed in the beginning of Sect. 3, where higher values give a more uniform M d ψ at the cost of less directionality.
First consider setting h = A so that ψ has compact support within a convex cone in the Fourier domain. The real part of the corresponding wavelet would, however, be a plate detector and not a line detector (Fig. 6). The imaginary part is already an oriented edge detector, 3

and so we set
where the real part of the earlier found wavelet vanishes by anti-symmetrization of the orientation distribution A while the imaginary part is unaffected. Note that the right-hand side of (52) and (53) where integration is performed over S p (n) denoting the great circle perpendicular to n. This transformation preserves the symmetry of A, so we have h Re (ϑ, ϕ) = h Re (ϑ). Thus, we finally set For an overview of the transformations, see Fig. 7. As discussed before, the additional splitting in low and high frequencies as described in Sect. 2.1.1 effectively causes a splitting in the radial function. How this affects the coverage of the Fourier domain is visualized in Fig. 8.

Efficient Implementation Via Spherical Harmonics
In Sect. 3.1 we defined the real part and the imaginary part of the wavelets in terms of a given orientation distribution.
In order to efficiently implement the various transformations (e.g., Funk transform) and to create the various rotated versions of the wavelet, we express our orientation distribution A in a spherical harmonic basis {Y m l } up to order L:  The "sharp" parts when the cones reach the center, however, cause the filter to be non-localized, which was solved in earlier works by applying a spatial window after filter construction. Right: By splitting the filter in lower and higher frequencies we solve this problem. In the figure we show g(ρ)A(n(ϑ, ϕ)) for the different filters, before applying the Funk transform to the orientation distribution A The spherical harmonics are given by where P m l is the associated Legendre function, = (−1) m for m < 0 and = 1 for m > 0 and with integer l ≥ 0 and integer m satisfying −l ≤ m ≤ l. For the diffusion kernel, which has symmetry around the z-axis we only need the spherical harmonics with m = 0, and we have the coefficients [20]: and Eq. (56) reduces to

Funk Transform
According to [23], the Funk transform of a spherical harmonic equals with P l (0) the Legendre polynomial of degree l evaluated at 0. We can therefore apply the Funk transform to a function expressed in a spherical harmonic basis by a simple transformation of the coefficients a m l → P l (0) a m l .

Making Rotated Wavelets
To make the rotated versions ψ n of wavelet ψ, we have to find h n inψ n (ω) = g(ρ) h n (ϑ, ϕ). To achieve this we use the steerability of the spherical harmonic basis. Spherical harmonics rotate according to the irreducible representations of the SO(3) group D l m,m (γ, β, α) (Wigner-D functions [41]): Here α, β and γ denote the Euler angles with counterclockwise rotations, where we rely on the convention R γ,β,α = R e z ,γ R e y ,β R e z ,α . This gives where c m l are the coefficients of h given by The filters from this section are summarized in the following result: Result 1 Let A : S 2 → R + be a function supported mainly in a sharp convex cone around the z-axis and symmetrically around the z-axis and g as radial function of Eq. (50). Then A provides our waveletψ in the Fourier domain viâ with ω = ρ n ω = ρ n(ϑ, ϕ). The real part of ψ is a tube detector given by The imaginary part of ψ n is an edge detector given by When expanding the angular part in spherical harmonics up to order L and choosing A = G S 2 s o : we have the following wavelet in the Fourier domain and the coefficients of A andψ relate via We obtain rotated versions of our filter viâ with n = n(β, γ ).
As we do not have analytical expressions for the spatial wavelets ψ n , we sample the filter in the Fourier domain using Eq. (71) and apply a DFT afterward. The wavelet ψ is a proper wavelet with fast reconstruction property (Definition 1).

Remark 3
The heat kernel on S 2 is given by where we recall Eq. (68). Because of the exponential decay with respect to l, we can describe the diffusion kernel well with the first few coefficients. In all experiments we truncate at smallest L such that a 0 L /a 0 0 < 10 −3 (e.g., L = 21 for

Stability of the Discrete Transformation with Fast Reconstruction for Filters of Result 1
To make a fast reconstruction by summation possible (requirement 2), we need a proper wavelet set with the fast reconstruction property (Definition 2) with N d ψ ≈ 1. We now focus on finding bounds for N d ψ such that we can choose our parameters in a deliberate way.
here the norm is the 2 -norm on C 2l+1 .
Proof First we expand function N d ψ in spherical harmonics: We have g(ρ) = 1 for ρ = ω ≤ ρ 0 , but we still need to quantify the angular part. We define This varying component should remain small. We use the Cauchy-Schwarz inequality for each order l:  Fig. 9).

Wavelet Design with Continuous Fourier Transform and Analytical Description in the Spatial Domain
In the previous section we described wavelets which were analytical in the Fourier domain and were sampled and inverse discrete Fourier transformed to find the wavelets in the spatial domain. To get more control on the wavelet properties in both the spatial and Fourier domain, it would be convenient to have an analytical description of the wavelets in both domains. This could be achieved by expressing the wavelets in a basis for which we have analytical expressions for the Fourier transform. We will now discuss 2 such options for the basis and describe filters expressed in them.

A Review on Expansion in the Harmonic Oscillator Basis
The first basis in which we could expand our wavelets are the eigenfunctions of the harmonic oscillator H = x 2 − Δ, which was also studied in [25,67]. We will quickly review this work and show the problems which were encountered when using this basis, before moving onto an alternative basis in the next section which aims to solve these problems. When using the eigenfunctions of the harmonic oscillator as a basis, the idea is that operator H and the Fourier transform commute (F • H = H • F) and eigenfunctions of H are also eigenfunctions of F. We then expand our wavelets in these eigenfunctions restricting ourselves to eigenfunctions which are symmetric around the z-axis: the spherical harmonics with m = 0. The wavelet is then given by with Y m l the spherical harmonics, (r , θ, φ) and (ρ, ϑ, ϕ) spherical coordinates for x and ω, respectively, i.e., x = (r sin θ cos φ, r sin θ sin φ, r cos θ), ω = (ρ sin ϑ cos ϕ, ρ sin ϑ sin ϕ, ρ cos ϑ), and g l n given by where L (ν) n (ρ) is the generalized Laguerre polynomial. We then choose the case with least radial oscillations α n l = α l δ 0 n . If we then choose we have that M ψ (ω) approximates 1 for all ω ∈ R 3 in the Fourier domain as L → ∞ and we get the following wavelet 4 [27]: 4 The series in (81) converges point-wise but not in L 2 -sense. So for taking the limit L → ∞ one must rely on a distributional wavelet transform, since ψ L→∞ i.e., it is contained in the dual of the 4 th -order Sobolev space. Such technicalities do not occur in the Zernike basis as we will see later in Sect. 4.
For this wavelet we have an analytical description in both spatial and Fourier domain.

Expansion in the Zernike Basis
The wavelets from the previous subsection have some unwanted properties such as poor spatial localization (long oscillations) and the fact that the wavelets maximum does not lie at the wavelets center, see Fig. 10. A possible explanation is that the basis used is orthogonal on the full L 2 (R 3 ) space and not limited to the ball in the Fourier domain, and truncation of this basis at the Nyquist frequency could lead to oscillations. An alternative basis for the unit ball is the Zernike basis which we can scale to be a basis in the Fourier domain for ball-limited images f ∈ L 2 (R 3 ), recall Eq. (2). The 2D Zernike basis is often used in applications as optics, lithography and acoustics [1,12,15], since efficient recursions can be used for calculating the basis functions and analytic formulas exist for many transformations among which the Fourier transform. The basis is therefore highly suitable for problems such as aberration retrieval where a wave function in the Fourier domain should be estimated from measurements in the image domain [24].
Orthogonal polynomials of several variables on the unit ball were considered by Louis [49] in the context of inversion of the Radon transform for tomographic applications. For a modern treatment of orthogonal polynomials on the unit ball, see [33,Ch. 4]. Here we will use the generalized Zernike functions [43], which vanish to a prescribed degree at the boundary of the unit ball, with explicit expansion results for particular functions supported by the unit ball. Since this basis is orthogonal on the unit ball and has explicit results for the Fourier transform it is highly suitable for the application in mind. Derivations for the results used in the following section can be found in [43].

The 3D Generalized Zernike Basis
The generalized Zernike functions are given by with spherical coordinates integer n, l ≥ 0 such that n = l + 2 p, integer p ≥ 0 and m = −l, −l + 1, . . . , l and α > 0. The angular part is the spherical harmonic function and the radial part is given by where P (α,l+ 1 2 ) p denotes the Jacobi polynomial. The generalized Zernike functions are orthogonal on the unit ball with δ the Kronecker delta and with normalization factor in which (x) α = Γ (x+α) Γ (x) is the (generalized) Pochhammer symbol.

Fourier Transform The inverse Fourier transform of the generalized Zernike function
is given by with x = r n(θ, φ) and Here J a and j a are the Bessel functions and spherical Bessel functions [55]. For integer α, the expression in Eq. (89) for q > 0 reduces to Expansion of Separable Functions An additional constraint for the wavelets is that they should be separable in the Fourier domain, i.e., (Fψ)(ω) = F(ω) = A(ϑ, ϕ)B(ρ). When expanding such a function in the generalized Zernike basis, we can split the coefficients in radial coefficients and angular coefficients The coefficients c m,α n,l in (92) reflect the separation of F as a product of an angular and radial factor as well as a corresponding separation of the generalized Zernike basis functions in (82). In the latter, the index l appears both in the angular and radial factor. Thus we have while for all l = 0, 1, . . .
For each l, the radial functions R l,α n with n varying are a basis for functions defined on the interval [0, 1]. For separable functions, we expand the same radial function B(ρ) for each l, and it can be shown that there is a recursion formula for the radial coefficients [43].

Wavelets
We now choose appropriate radial and angular functions for our wavelets expressed in the generalized Zernike basis. A(n(ϑ, ϕ) n(ϑ, ϕ)) for which the spherical harmonic coefficients are given by (58). After this we apply the same transformations (Funk transform and anti-symmetrization) to obtain the angular part of the wavelet.

Flat Radial Profile for All-Scale Transform
Recall the procedure of splitting of the lowest frequencies as described in Sect. 2.1.1 resulting in filters ψ 0 and ψ 1 . In this section we design a radial function for ψ 1 which is relevant for further processing. Furthermore, we already have an analytical description for φ 0 , which we set to φ 0 = G s ρ (see Eq. (37)).
The radial function of ψ 0 should therefore approximate B(ρ) = 1 − G s ρ (ρ) on the interval [0, ] and should smoothly go to zero when approaching the edges of the interval. For the moment, we set = 1 and we include the scaling later. To start, we define the function To obtain a flatter function we multiply the function B α,β with a second-order Taylor expansion of the reciprocal function ρ → (B α,β (ρ)) −1 around the function's maximum obtained at see Fig. 11b. The resulting function is again a sum of functions of type (97) with different values for β, so we can find the coefficients b l,α n=l+2 p for the flattened function as well. For the specific case β = 2 we get the following flattened function with B max = B α,2 (ρ max ). For this flattened function the coefficients are given by with c i the coefficients of ρ 0 , ρ 2 and ρ 4 in the second-order Taylor series of the reciprocal. These coefficients follow from (100) and are given by The filters from this section are summarized in the following result:  ψ 1 and angular function A(n(ϑ, ϕ)) = G S 2 s o (n(ϑ, ϕ)) and expand in the generalized Zernike basis: The coefficients c 0 n,l follow by expanding A in spherical harmonics and B The spatial wavelet is given by

Experiments with the Filters and the Transform Between a 3D Image and Its Orientation Score
Before considering applications of the filters in the next section, we first compare filters obtained by DFT (Sect. 3) to filters expressed in the generalized Zernike basis (Sect. 4) and inspect the quality of the reconstruction.

Comparison of Wavelets Obtained via DFT and Analytical Expressions Using the Zernike Basis
First we compare the filters obtained by sampling in the Fourier domain followed by a DFT (Sect. 3) to the filters obtained by expansion in the Zernike basis (Sect. 4). Settings were chosen such that the radial functions of both wavelets matched best and the same settings for the angular function were used. In Fig. 12 we show that the filters are very similar in shape. We see no major artifacts caused by sampling followed by an inverse DFT.

Quality of the Reconstruction
A visual inspection of the reconstruction after the transformation and reconstruction procedure can be found in Fig. 13. As expected, a small amount of regularization is observed. We see no qualitative differences between the two reconstructions.

Applications
Next we present two applications of 3D image processing via invertible orientation score transforms. Figure 14 summarizes the different datasets on which we validate this kind of processing. In the invertible orientation score transform used in all experiments, we choose to use the wavelets from Result 1 and the default values of Table 1, unless stated otherwise. The reason for only using the wavelets from Result 1 is that the wavelets from Result 2 are over 10 times slower in our current implementations, partly due to the fact that the series in (106) is converging, but unfortunately not rapidly converging.

Background and Related Methods
Many methods exist for enhancing elongated structures based on nonlinear diffusion equations. Coherence-enhancing diffusion (CED) filtering [69] uses the structure tensor to steer the diffusion process to mainly apply diffusion along the elongated structures, therefore preserving the edges. One downside of this method is that at situations where multiple oriented structures occur at the same position, one of the structures gets destroyed. This renders this method not suitable for crossing structures and in 3D data bifurcating vessels. Interesting extensions dealing with crossings by analyzing the environment using higher-order derivatives have been proposed [63].
Methods that deal with crossings by applying coherenceenhancing diffusion via 2D orientation scores have been developed for 2D data [37,64]. Here, we propose an extension Fig. 13 Comparison of construction and reconstruction of data A.1 using the different types of filters with the same settings as in Fig. 12. In each row, from left to right, an iso-contour of the data and 3 slices through the center of the data along the three principal axis. Top: The original data. Middle: the data after construction and reconstruction using the filters from Result 1. Bottom: the data after construction and reconstruction using the filters from Result 2 of coherence-enhancing diffusion via 3D orientation scores to enhance elongated structures, while preserving crossings and bifurcating vessels. Preliminary results on artificial data have been shown in [32]. Here we show first results on real data, quantify the results and furthermore add additional adaptivity to the diffusion equation.

CEDOS
We now use the invertible orientation score transformation to perform data enhancement according to Fig. 2. Because R 3 × S 2 is not a Lie group, it is common practice to embed the space of positions and orientations in the Lie group of positions and rotations SE(3) by setting with R n any rotation for which R n · e z = n. The operators Φ which we consider are scale spaces on SE(3) (diffusions), and are given by HereW is the solution of a nonlinear diffusion equation: where in coherence-enhancing diffusion on orientation scores (CEDOS) D i j is adapted locally to initial conditionW (g, 0) based on exponential curve fit (see [32]), and with A i | g = (L g ) * A i | e the left-invariant vector fields on SE(3). The diffusion is better understood in locally adaptive frame {B i } 6 i=1 . Here B 3 follows from an exponential curve fits and points along our structure. B 1 and B 2 span the plane spatially perpendicular to our structure, and B 4 , B 5 and B 6 correspond to angular diffusion (embedding in S E(3) leads to a third angular dimension). Our diffusion then takes the diagonal form where we limit ourselves to diffusion of type D 11 = D 22 , and D 44 = D 55 = D 66 to preserve the data symmetry of Eq. (107). The diffusion system in Eq. (110) has been set up in previous work [45] with constant D 11 and D 33 .
In this work we include further adaptivity by making the constants D 11 and D 33 data adaptive. We aim to enhance oriented structures and reduce noise as much as possible. Therefore, non-oriented regions are smoothed isotropically by setting Discrete wavelet size 11 × 11 × 11 with c 1 a constant automatically set to the 50% quantile of s(Ũ ) and s(Ũ )(g) a measure for orientation confidence given by minus the Laplacian in the space orthogonal to the structure orientation B 3 : Since we want to stop diffusion when reaching the end of a structure, we set with c 2 a constant automatically set to the 50% quantile of We then obtain Euclidean invariant image processing via which includes inherent projection P ψ of orientation scores, even if Φ = Φ t maps outside of the space of orientation scores. We write W −1,ext ψ because we extend the inverse to L 2 (R 3 × S 2 ).

Quantification Via Contrast-to-Noise Ratio
In the next section, we quantify the denoising capabilities of the diffusion in Eq. (110) with the contrast-to-noise ratio. Next we will explain how we calculate this measure for our real medical image data.
For signal f with noise N the noisy data is given by: Given such data, the noise is quantified via the contrast-tonoise ratio (CNR): where σ ( f − f N ) denotes the standard deviation of the difference signal f − f N , and where the numerator denotes the contrast in our data. For real data we do not have a ground truth f but only noise signal f N and we will use the following estimation for the standard deviation of the noise and the contrast of the signal. First we estimate the contrast by determining the average value over manually segmented parts of the vessel given by region Ω S and background regions given by Ω B : For estimating the noise of the signal, we select regions for which the signal f can be expected to be constant (see Fig. 15). Given such a region Ω B we estimate the noise standard deviation by The contrast-to-noise ratio (CNR) is then given by where the numerator denotes the contrast in our data.

Results on Cone Beam CT Data of the Abdomen
We tested our method on real Cone Beam CT data of the abdomen (Fig. 14b). The data were acquired using a Philips Allura Xper FC20 system, using a Cone Beam CT backprojection algorithm (XperCT) to generate the final volumetric image.
To quantify our method, we segmented the vessels and selected background regions, see Fig. 15. We then applied CEDOS with different end times and computed the CNR for these different end times ranging from 1 to 6, see Fig. 17. Diffusion constants D 11 and D 33 were determined using Eqs. (111), (113), and we set D 11 = 0.001. For the orientation score transformation, we used s o = 1 2 (0.45) 2 and s ρ = 1 2 (16) 2 . For CED we used the following settings: α = 0.2 and c the 50%-quantile of κ (see [69]).
As one can expect, in all cases we recognize a peak since initially noise is reduced whereas later also contrast is reduced. Compared to Gaussian diffusion and CED, we reach a higher CNR (Fig. 16). Furthermore, in the CEDOS and CED case, the CNR does not decrease much when applying more diffusion making it more robust with respect to choice in diffusion time. The fact that CED does not achieve high CNR ratios and performs relatively bad in this test is that the diffusion matrix is not designed to reduce background noise (which is used here to quantify noise) but mainly to enhance orientated patterns; for this reason we also set α relatively high to still achieve noise reduction in the background regions.
For 3D visualization of the diffusion results for optimal diffusion time (determined from the CNR), see Fig. 18.
Here we see that compared to Gaussian diffusion our anisotropic diffusion reduced more noise while still maintaining the important structures. A similar thing is achieved by anisotropic diffusion in CED but bifurcating vessels are destroyed by this method, as in this method the diffusion is mainly performed along the orientation of one of the vessels at the bifurcation (see the black circles in Fig. 18).

Influence of CEDOS on Vessel Edge Location
In order to test the influence of our regularization method on vessel features, we implemented a simple edge detection algorithm. We manually selected positions in the data and detect the edges in the vessel cross sections by extracting radial profiles from the centerline outward and looking for  Fig. 17 the minimum in first-order gaussian derivative. Compared to Gaussian regularization, the vessel edge position is not influenced by our regularization method (Fig. 19), which is highly important for applications which rely on accurate vessel lumen measurements, e.g., stent positioning and navigation of endovascular devices. The key explanation for this benefit is that at the vessel we get a very low D 11 (Eq. 111) and therefore we smooth only along the vessel.

Background and Related Methods
In this section we propose a tubularity measure based on the edge information in our orientation scores. This tubularity measure is then used for vessel width measurements and could be used for vessel segmentation. Tubularity measures are designed to have a high response on the centerline of tubular structures. One such tubularity measure is the image gradient flux filter [14] which is used for vessel segmentation. An extension of the gradient flux filter is the optimally oriented flux filter [47] which introduces the notion of oriented flux making the filter orientation sensitive.
A 2D tubularity measure based on orientation scores was proposed in [9]. The advantage of this tubularity measure is that it included nonlinearity and that the implementation via orientation scores still has a response at crossing vessels. Here we propose an extension of this tubularity measure to 3D making use of 3D orientation scores.

Tubularity Via Orientation Scores
For the tubularity measure we detect edges in the plane perpendicular to orientation n. Within this plane, the product of two opposite edges at radius r and in-plane orientation θ at For all datasets, we show one iso-contour at μ BG + 0.7(μ FG − μ BG ), where μ BG and μ FG are the mean of the background and foreground in the (processed) data determined using the selected regions used for determining the CNR. For a better impression of the full volume, see Fig. 18. We plot results for three different times according to Fig. 16, where case (a) corresponds to optimal diffusion time for Gaussian regularization and case (c) corresponds to optimal CEDOS which is also approximately equal to optimal CED time. We see that CEDOS preserves the complex vascular geometry position x ∈ R 3 is given by where Im + (z) = max{0, Im(z)} and n ⊥ (θ ) = cos θ e 1 + sin θ e 2 , with {e 1 , e 2 } an orthogonal basis for the orthogonal complement of n = span{n}. The product of the two Im + edge responses in Eq. (120) yields a better performance than taking the sum, as is done in [9,Fig. 12.2] and [19]. The idea behind taking the product instead of the sum is that we need a high edge response in both directions. See Fig. 20 for a schematic visualization. Since for a real tube all edges should be present, and we do no want any response for, e.g., plate structures, we take a minimum over the perpendicular orientations parametrized by θ . This is done as follows For the angular regularization kernel we use the simple 2π periodic 1D-diffusion kernel over θ , instead of the true diffusion kernel (72), where for practical values we use σ V o = π/8 the series of rapidly converging kernels can be reasonably truncated already at the first term k = 0.
From the tubularity measure (121), we extract the following features: Here s t (x) is the tubularity confidence which is a measure for how certain we are at least one tubular structure is present at position x. The features n * (x) and r * (x) are the orientation and radius of optimal tubularity response at position x.  x, n, r , θ). The edge is expected to have outward orientation n ⊥ (θ)

Results on Artificial Data
For the validation of our tubularity measure, we constructed 18 artificial datasets with a random tubular structure with randomly varying radius (Fig. 14c). For the tubularity measure we used the following settings: σ V o = π/8 and we discretized the θ -integral in Eq. (121) using 8 orientations.
As validation we compared the optimal radius to the ground truth radius and inspect the tubularity confidence. The transversal profiles in the tubes are modeled by Gaus- The tubularity confidence s t (x) (max of tubularity over radius and orientation) with radius of max response r * (x) in color. Right: Measured radius r * (x) against ground truth radius r GT (x) on the ground truth centerline. The opacity of the plotted points is linearly scaled with the tubularity confidence sians, and we define the ground truth radius r GT (x) = σ (x) at the standard deviation σ (x) of the Gaussian. Thereby, the tube boundaries are at the inflection point of the transversal Gaussian profiles.
The tubularity confidence is selective on the vessel centerline, and the found optimal radius is a reasonable estimation of the ground truth radius, see Figs. 21 and 22, and follows the correct trend although outliers typically tend to produce an overestimate.

Results on 3D Rotational Angiography Data of the Brain in Patients with Arteriovenous Malformation (AVM)
We applied the tubularity measure to 3D rotational angiography of the brain in patients with arteriovenous malformation (Fig. 14a). The data were acquired using a Philips Allura Xper FC20 system, using a 3D rotational angiography backprojection algorithm (3DRA) to generate the final volumetric image. For the tubularity measure we used the following settings: σ V o = π/8 and we discretized the θ -integral in Eq. (121) using 12 orientations. For the orientation score transformation we used s o = 1 2 (0.4) 2 , s ρ = 1 2 (5) 2 and evaluated on a grid of 21 × 21 × 21 pixels and default settings for the other parameters.
In Fig. 23 we show our tubularity measure for this medical data. The tubularity measure gives sharp responses for vessel centerlines. We also show optimal orientation n * (x) and a simple segmentation given by the 0-iso-contour of the where c i are the positions given by the 1% quantile of highest responses.

Conclusion
We presented theory and filters for the 3D orientation score transformation which is valuable in handling complex oriented structures in 3D image data. Then we showed applications of this transformation. First, we proposed filters for a 3D orientation score transform. We presented two types of filters: The first uses a discrete Fourier transform and the second is designed in the 3D generalized Zernike basis which allowed us to find analytical expressions for the spatial filters. Both filters allowed for an invertible transformation. The filters and the quality of the reconstruction are assessed in Sect. 5, where we showed that the discrete filters approximate their analytical counterparts well. We also verified the invertibility of our transformation by showing data reconstructions of real medical data.
The orientation score transform was then used in two different applications. In the first we presented an extension of coherence-enhancing diffusion via 3D orientation scores which we applied to real 3D medical data and showed our method effectively reduced noise while maintaining the complex vessel geometry. In the second application we propose a new nonlinear tubularity measure via 3D orientation scores. The tubularity measure has sharp responses for vessel centerlines, and we showed its use in radius extraction and complex vessel segmentation.
In this work basic applications of the tubularity measure are shown. Future work would include using the tubularity measure in more advanced vessel segmentation procedures. Furthermore, many other applications exist for the 3D orientation scores and many techniques developed for 2D orientation scores can now be extended to 3D. First steps are presented in this paper, where the extension of 2D CEDOS and a 2D tubularity measure via 2D orientation scores are given. It is also interesting to explore the nonlinear diffusion procedure (Eq. (110)) for contextual processing of diffusionweighted MRI and to compare with existing approaches [58,68].
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

(125)
Remark 4 In Sect. 4, we consider the modified Zernike basis in which case g l n andg l n are given by, respectively, (85) and (82), whereas for the harmonic oscillator basis one has g l n = i l (−1) n+lgl n given, respectively, by (77) and (79). We obtain steerability via finite series truncation at n = N and l = L. Then we rotate the steerable kernels via the Wigner-D functions D l 0,m (γ , β, 0) ∈ R (cf. Sect. 3.2.3) and one obtains the following steerable implementations of orientation scores: where n = (cos γ sin β, sin γ sin β, cos β) T , denotes correlation, the overline denotes complex conjugation and with function product (g l n ⊗ Y m l )(x) =g l n (r ) Y m l (θ, φ).