Design and Processing of Invertible Orientation Scores of 3D Images

Journal of Mathematical Imaging and Vision, Mar 2018

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 \(U: {\mathbb {R}} ^ 2\times S ^ 1 \rightarrow {\mathbb {C}}\) were introduced, which already showed their use in a variety of applications. Here we extend this work to 3D orientation scores \(U: {\mathbb {R}} ^ 3 \times S ^ 2\rightarrow {\mathbb {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.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

https://link.springer.com/content/pdf/10.1007%2Fs10851-018-0806-0.pdf

Design and Processing of Invertible Orientation Scores of 3D Images

Journal of Mathematical Imaging and Vision https://doi.org/10.1007/s10851 Design and Processing of Invertible Orientation Scores of 3D Images M. H. J. Janssen 0 1 2 3 A. J. E. M. Janssen 0 1 2 3 E. J. Bekkers 0 1 2 3 J. Oliván Bescós 0 1 2 3 R. Duits 0 1 2 3 M. H. J. Janssen 0 1 2 3 A. J. E. M. Janssen 0 1 2 3 0 Department of Mathematics and Computer Science, Eindhoven University of Technology , Eindhoven , The Netherlands 1 CASA, Eindhoven University of Technology , Eindhoven , The Netherlands 2 J. Oliván Bescós 3 Philips , Interventional X-ray, Eindhoven , The Netherlands 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 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 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. Orientation scores; 3D wavelet design; Zernike polynomials; Scale spaces on SE(3); Coherence-enhancing diffusion; Tubular structure detection; Steerable 3D wavelet 1 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 : R3 → R with a set of oriented filters to construct a 3D orientation score U : R3 × S2 → 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 : R2 × S1 → 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. 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 1.1 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. 1.2 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. 2 Invertible Orientation Scores 2.1 Continuous Orientation Score Transform Throughout this article we use the following definition of the Fourier transform on R3: fˆ(ω) = (F f )(ω) = R3 e−iω·x f (x)dx. with ball B correlation = {ω ∈ R3 | ω < } of radius with an anisotropic kernel > 0, by (Wψ [ f ])(x, n) = (ψn f )(x) = R3 ψn(x − x) f (x )dx . This is illustrated in Fig. 4. Here ψ ∈ L2(R3) ∩ L1(R3) is a wavelet aligned with and rotationally symmetric around the z-axis, and ψn ∈ L2(R3) the rotated wavelet aligned with n given by ψn(x) = ψ (RnT x). Here Rn ∈ S O( 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 ( 3 ) ( 4 ) ( 5 ) ( 6 ) ( 7 ) ( 8 ) f (x) = (Wψ−1[Wψ [ f ]])(x) ⎡ ⎡ = F−1 ⎢ Mψ−1F ⎢ x˜ → ⎣ ⎣ S2 ⎤ ⎤ (ψˇn Wψ [ f ](·, n))(x˜)dσ (n)⎥ ⎥ (x), ⎦ ⎦ with ψˇ n(x) = ψn(−x). The function Mψ : R3 → 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 ∃M>0,δ>0∀ω∈B : 0 < δ ≤ Mψ (ω) ≤ M < ∞. In practice it is best to aim for Mψ (ω) ≈ 1, in view of the condition number of transformation Wψ : L2 (R3) → L2 (R3 × S2) given by: An invertible orientation score Wψ [ f ] : R3 × S2 constructed from a given ball-limited 3D dataset → C is cond(Wψ ) = Wψ Wψ−1 M = δ , f ∈ L2 (R3) = { f ∈ L2(R3) | supp(F f ) ⊂ B }, where M and δ are assumed to be tight bounds in ( 7 ). In the codomain spatial frequencies are again limited to the ball: Fig. 4 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 ex . Bottom left: This is repeated for a discrete set of filters with different orientations. Bottom right: The collection of 3D L2 (R3 × S2) = {U ∈ L2(R3 × S2)|∀n∈S2 U (·, n) ∈ L2 (R3)}. ( 9 ) 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 We can further simplify the reconstruction for wavelets for which the following additional property holds: Also, in the case we have Mψ (ω) = 1 for ω ∈ B L2-norm preservation we have Nψ (ω) = ψˆ n(ω) dσ (n) ≈ 1. S2 ( 12 ) ( 13 ) f 2L2(R3) = Wψ f 2L2(R3×S2), for all f ∈ L2 (R3), ( 10 ) and reconstruction Eq. ( 5 ) simplifies to 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 ∈ L2 (R3) to an orientation integrated score Aψ ( f ) = S2 Wψ f (·, n) dσ (n). max Nψ (ω) ω∈B Its condition number is given by cond(Aψ ) = min Nψ (ω) . ω∈B 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 bounds1 δ, M > 0, 0 < ε 1. Furthermore, let be an a priori maximum frequency of our ball-limited image data. Then, a wavelet ψ ∈ L2(R3) ∩ L1(R3) is called a proper wavelet if where Rez,α ∈ S O( 3 ) is the 3D rotation around axis ez over angle α. If, moreover, one has 3.) ∃ 21 < 0< ∀ω∈B 0 : Nψ (ω) ∈ [1 − ε, 1 + ε], 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. ( 14 ) ( 15 ) ( 16 ) ( 17 ) 2.1.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 ψ = ψ0 + ψ1, with ψˆ 0 = Gˆ sρ ψˆ , ψˆ 1 = (1 − Gˆ sρ )ψˆ , ( 18 ) with Gaussian window in the Fourier domain given by Gˆ sρ (ω) = e−sρ ω 2 , 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 (Wψ1 [ f ])(x, n) = (ψ1,n f )(x), and one for the low-frequency components (Wψ0 [ f ])(x, n) = (ψ0,n f )(x). Here we again have ψi,n(x) = ψi (RnT x), as in Eq. ( 4 ). The vector transformation is then defined as Wψ [ f ] = (Wψ0 [ f ], Wψ1 [ f ]). For this transformation we have the exact reconstruction formula f (x) = (Wψ−1Wψ f )(x) = F −1 Mψ−1F x˜ → S2 (ψˇ 1,n Wψ1 [ f ](·, n))(x˜ ) + (ψˇ 0,n Wψ0 [ f ](·, n))(x˜ )dσ (n) (x) ( 23 ) ( 19 ) ( 20 ) ( 21 ) ( 22 ) with 1 In practice we choose the default values δ = 81 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. 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. Lemma 1 Let ψ ∈ L2(R3) ∩ L1(R3) such that Eq. ( 7 ) holds, δ = minω∈B Mψ (ω) and M = maxω∈B Mψ (ω). Then the condition number of Wψ : L2 (R3) → L2 (R3 × S2) is given by | cond(Wψ )|2 = Wψ 2 Wψ−1 2 M = δ . The condition number of Wψ : L2 (R3) → L2 (R3 × S2) obtained from Wψ by performing an additional splitting in low- and high-frequency components is given by | cond(Wψ )|2 = Wψ 2 Wψ−1 2 ≤ 2δM , thereby guaranteeing that stability is maintained after performing the splitting. Proof First, we find the condition number of Wψ : | cond(Wψ )|2 = f 2L2 sup sup f ∈L2(R3) Wψ f 2L2 · f ∈L2(R3) Wψf f2L22L2 . For the first factor in Eq. ( 27 ) we find sup f 2L2 sup F f 2L2 f ∈L2(R3) Wψ f 2L2 = f ∈L2(R3) F Wψ f 2L2 = = sup f ∈L2(R3) S2 R3 |ψˆ n(ω)|2| fˆ(ω)|2dωdσ (n) R3 | fˆ(ω)|2dω sup f ∈L2(R3) R3 Mψ (ω)| fˆ(ω)|2dω R3 | fˆ(ω)|2dω 1 = ωm∈aBx Mψ (ω) . 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 1 M cond(Wψ ) = ωm∈aBx Mψ (ω) · ωm∈aBx Mψ (ω) = δ . Similarly the condition number of Wψ is given by 1 cond(Wψ ) = ωm∈aBx Mψ (ω) · ωm∈aBx Mψ (ω). Next we express Mψ in Mψ of the original wavelet as 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 Lψ0 [ f ] via Lψ0 [ f ](x) = (φ0 f )(x), with φ0 = S2 ψ0,n dσ (n). ( 36 ) Lψ0 [ f ](x) ( 35 ) For a design with Nψ (ω) = 1 for all ω ∈ B , we have φˆ0 = Gˆ sρ and so N d ψ = ψˆ ni (ω)Δi ≈ 1, 1 x 2 φ0(x) = Gsρ (x) = (4π sρ )3/2 e− 4sρ . Then Eq. ( 35 ) becomes the image reconstruction can be simplified by a summation over orientations: S2 Wψ1 [ f ](x, n) dσ (n) + (Gsρ ∗ f )(x). 2.2 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 simple2 electrostatic repulsion model [18]. Assume we have a number No of orientations V = {n1, n2, . . . , nNo } ⊂ S2, and define the discrete invertible orientation score Wψd [ f ] : R3 × V → C by (Wψd [ f ])(x, ni ) = (ψni f )(x). The exact reconstruction formula is in the discrete setting given by f (x) = ((Wψd )−1[Wψd [ f ]])(x) = F −1 (Mψd )−1F ( 37 ) ( 38 ) ( 39 ) ( 40 ) No i=1 No i=1 No i=1 No i=1 x˜ → (ψˇ ni Wψd [ f ](·, ni ))(x˜ ) Δi (x), with Δi the discrete spherical area measure ( iN=o1 Δi = 4π ) which for reasonably uniform spherical sampling can be approximated by Δi ≈ 4Nπo , and Mψd (ω) = ψˆ ni (ω) 2 Δi . 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 2 For our applications this numerical approach in [18] was sufficient. There are alternative approaches [40] with profound analysis on optimized spherical samplings minimizing squared quadrature errors in Fourier transforms. For this reconstruction by summation we can analyze the stability via the condition number of the mapping that maps an image f ∈ L2 (R3) to an orientation integrated score No i=1 Adψ ( f ) = ψ ni f Δi . This transformation has condition number cond(Adψ ) = max Nψd (ω) ω∈B min Nψd (ω) . ω∈B 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 {ψni ∈ L2 (R3) ∩ L1(R3) | i = 1, . . . , No}, with a reasonable uniform spherical sampling (Δi ≈ 4Nπo ), constructed as rotated versions of ψ is called a proper wavelet set if ( 42 ) ( 43 ) ( 44 ) ( 45 ) ( 46 ) ( 47 ) ( 41 ) 3.) ∃ 21 < 0< ∀ω∈B 0 : Nψd (ω) ∈ [1 − ε, 1 + ε], where Rez,α ∈ S O( 3 ) is a 3D rotation around axis ez over angle α. If, moreover, one has 2.3 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. 3 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. Thanks to the splitting in Sect. 2.1.1 we no longer need the extra spatial window used there. Our 3D transformation using the 3D cake wavelets should fulfill a set of requirements, compare [37]: 1. The orientation score should be constructed for a finite number (No) of orientations. 2. The transformation should be invertible, and reconstruction should be achieved by summation. Therefore, we aim for N d ψ ≈ 1. Additionally, to guarantee all frequencies are transferred equally to the orientation score domain we aim for M d ψ ≈ 1. The set should be a proper wavelet set with fast reconstruction property (Definition 2). 3. The kernel should be strongly directional. 4. The kernel should be separable in spherical coordinates in the Fourier domain, more explicitly (F ψ )(ω) = g(ρ)h(ϑ, ϕ), with ω = (ωx , ωy , ωz ) = (ρ sin ϑ cos ϕ, ρ sin ϑ sin ϕ, ρ cos ϑ ). ( 49 ) 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. 3.1 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 Fig. 5 Radial part g of ψˆ , see Eq. ( 50 ) and radial parts g0 and g1 of ψˆ0 and ψˆ1 after splitting according to Sect. 2.1.1. The parameter γ controls the inflection point of the error function, here γ = 0.8. The steepness of the decay when approaching ρN is controlled by the parameter σer f with default value σer f = 31 (ρN − ). At what frequency the splitting of ψˆ in ψˆ0 and ψˆ1 is done is controlled by parameter σρ = 2sρ , see Eq. ( 18 ) ( 18 ) in Sect. 2.1.1 takes place. Following requirement 4 we only consider polar separable wavelets in the Fourier domain, so that (F ψ )(ω) = g(ρ)h(ϑ ). 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, 1 g(ρ) = 2 1 − erf ρ − σer f with erf(z) = √2π 0z e−x2 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 = γ ρN , 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 = 31 (ρ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 ( 50 ) ( 51 ) Orange: Positive iso-contour. Blue: Negative iso-contour. Parameters used: so = 21 (0.25)2, σerf = 3, γ = 0.85 and evaluated on a grid of 51 × 51 × 51 pixels (Color figure online) the highest frequencies which causes a small amount of blurring. We now need to find an appropriate angular part h for the cake wavelets. First, we specify an orientation distribution A : S2 → 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]: 2 A(n(ϑ, ϕ)) = GsSo (n(ϑ, ϕ)), with so > 0 and n(ϑ, ϕ) = (sin ϑ cos ϕ, sin ϑ sin ϕ, cos ϑ ). The parameter so 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 ( 52 ) 1 hI m (ϑ ) = 2 ( A(n(ϑ, ϕ)) − A(−n(ϑ, ϕ))) 1 = 2 2 2 GsSo (n(ϑ, ϕ)) − GsSo (−n(ϑ, ϕ) , 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 ) is the same for all values of ϕ. As to the construction of h Re, there is the general observation that 3 This can be observed in Fig. 6, where we recognize a planar derivative filter in the direction of n because it is anti-symmetric and has a mainly positive and mainly negative side. we detect a structure that is perpendicular to the shape in the Fourier domain, so for line detection we should aim for a plane detector in the Fourier domain. To achieve this we apply the Funk transform to A, and we define h Re(ϑ, ϕ) = F A(n(ϑ, ϕ)) 1 = 2π Sp(n(ϑ,ϕ)) A(n ) ds(n ), where integration is performed over Sp(n) denoting the great circle perpendicular to n. This transformation preserves the symmetry of A, so we have h Re(ϑ, ϕ) = hRe(ϑ ). Thus, we finally set h(ϑ ) = hRe(ϑ ) + hI m (ϑ ). 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. ( 54 ) ( 55 ) ( 53 ) 3.2 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 {Ylm } up to order L : A(n(ϑ, ϕ)) = Fig. 7 Cake Wavelets. Top: 2D cake wavelets. From left to right: Illustration of the Fourier domain coverage, the wavelet in the Fourier domain and the real and imaginary part of the wavelet in the spatial domain [10]. Bottom: 3D cake wavelets. Overview of the transformations used to construct the wavelets from a given orientation distribution. Upper part: The wavelet according to Eq. ( 53 ). Lower part: The wavelet according to Eq. ( 54 ). IFT: Inverse Fourier Transform. Parameters used: so = 21 (0.25)2, γ = 0.85 and evaluated on a grid of 81 × 81 × 81 pixels 2D Cake-wavelets 3D Cake-wavelets An symmetrize Funk Transform IFT Edge Detector [Im] Line Detector [Re] IFT IFT Edge Detector [Im] Tube Detector [Re] Fig. 8 Coverage of the Fourier domain before and after splitting according to Sect. 2.1.1. Left: The different wavelets cover the Fourier domain. 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 Ylm (ϑ, ϕ) = 2l + 1 4π (l − |m|)! eimϕ P |m|(cos ϑ ), (l + |m|)! l ( 57 ) where Plm 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]: alm = ⎩ ⎧⎨ 0 m = 0, 24l+π1 e−l(l+1)so m = 0, and Eq. ( 56 ) reduces to L l=0 with Pl (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 alm → Pl (0) alm . 3.2.2 Anti-symmetrization We have Ylm (π − ϑ, ϕ + π ) = (−1)l Ylm (ϑ, ϕ). We therefore anti-symmetrize the orientation distribution, see Eq. ( 53 ), via alm → (1−(2−1)l ) alm . 3.2.3 Making Rotated Wavelets To make the rotated versions ψn of wavelet ψ , we have to find hn in ψˆ n(ω) = g(ρ) hn(ϑ, ϕ). To achieve this we use the steerability of the spherical harmonic basis. Spherical harmonics rotate according to the irreducible representations of l the SO( 3 ) group Dm,m (γ , β, α) (Wigner-D functions [41]): ( 58 ) ( 60 ) where clm are the coefficients of h given by clm = Pl (0) alm + Because both anti-symmetrization and Funk transform preserve the rotational symmetry of A, we have h(ϑ, ϕ) = lL=0 cl0Yl0(ϑ, ϕ), and Eq. ( 62 ) reduces to L l ( 59 ) hn(ϑ, ϕ) = cl0 D0l,m (γ , β, 0) Ylm (ϑ, ϕ). l m =−l RRγ,β,α Ylm (ϑ, ϕ) = Dm,m (γ , β, α)Ylm (ϑ, ϕ). ( 61 ) l Here α, β and γ denote the Euler angles with counterclockwise rotations, where we rely on the convention Rγ,β,α = Rez,γ Rey,β Rez,α. This gives hn(ϑ, ϕ) = RRγ,β,α h (ϑ, ϕ) = clm Dml,m (γ , β, α)Ylm (ϑ, ϕ), ( 62 ) ( 63 ) ( 64 ) ( 65 ) ( 66 ) ( 67 ) ( 68 ) ( 69 ) ( 70 ) Result 1 Let A : S2 → 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 via ψˆ (ω) = g(ρ) (F A(nω) + A(nω) − A(−nω)) , with ω = ρ nω = ρ n(ϑ, ϕ). The real part of ψ is a tube detector given by Re(ψ ) = F −1 (ω → g(ρ)(F A)(nω)) . The imaginary part of ψn is an edge detector given by 1 Im(ψ ) = i F −1 (ω → g(ρ) ( A(nω) − A(−nω))) . When expanding the angular part in spherical harmonics up to order L and choosing A = GsSo2 : A(n(ϑ, ϕ)) = al0Yl0(ϑ, ϕ), al0 = 2l + 1 e−l(l+1)so , 4π we have the following wavelet in the Fourier domain L l=0 L l=0 ψˆ (ω) = g(ρ) cl0Yl0(ϑ, ϕ), and the coefficients of A and ψˆ relate via cl0 = We have g(ρ) = 1 for ρ = ω ≤ ρ0, but we still need to quantify the angular part. We define YlN = (74) (76) We obtain rotated versions of our filter via L 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 S2 is given by GsSo2 (n(ϑ, ϕ)) = al0Yl0(ϑ, ϕ) 24l+π1 e−so(l+1)l Yl0(ϑ, ϕ), (72) 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 aL0 /a00 < 10−3 (e.g., L = 21 for so = 21 (0.25)2). 3.3 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. Proposition 1 Let {ψn(βi ,γi ) | i = 1 . . . No} be a set of wavelets constructed via the procedure in Result 1. Then we have bounds on Nψd given by 1 − dl with dl = (dlm )lm=−l , dlm = No cl0 · Δi · D0l,m (0, βi , γi ) and i=1 here the norm is the 2-norm on C2l+1. Proof First we expand function Nψd in spherical harmonics: Nψd (ω) = F [ψni ](ω)Δi = g(ρ) hni (ϑ, ϕ)Δi No i=1 No i=1 l=1 from which (73) follows. See Fig. 9 for visual inspection of bounds of Mψd and N d , ψ and numerical results for the bounds of N d . ψ Corollary 1 Given our analytical bounds (73) from Proposition 1 and No = 42, we can guarantee that our set of wavelets from Result 1 is a proper wavelet set with fast reconstruction property according to Definition 2 with = 0.05 when choosing parameter s0 0.04. In practice we have a proper wavelet set with fast reconstruction property already for smaller values of so (see Fig. 9). 4 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. lower bounds of N d . Comparison of the bounds according to ψ Eq. (73) (filled blue line) and numerical results (dashed blue line) of the bounds by a very fine sampling of the sphere (≈ 500 orientations). Furthermore, we show 1 + and 1 − (orange dashed lines) for = 0.05 (Color figure online) with Ylm 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 ϑ ), Fig. 9 Inspection of the stability of the transformation for different values of so given an orientation distribution A = GsSo2 and for No = 42. Left: Spherical plot of A and the angular part of polar separable function Nψd and Md . Orientation coverage is more uniform as the plots for ψ Nψd and Mψd look more like a ball. Right: The upper and 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. 4.1 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 αln (−1)n+l i l gnl (ρ) Yl0(ϑ, ϕ), (77) n=0 l=0 ∞ L n=0 l=0 (78) (79) (80) and gnl given by gnl (ρ) = ρ1 Enl+ 21 (ρ), Enν (ρ) = 2(n!) Γ (n + ν + 1) where L (nν)(ρ) is the generalized Laguerre polynomial. We then choose the case with least radial oscillations αln = αl δn0. If we then choose αl = 3 Γ l + 2 , Γ (l + 1) we have that Mψ (ω) approximates 1 for all ω ∈ R3 in the Fourier domain as L → ∞ and we get the following wavelet4 [27]: 4 The series in (81) converges point-wise but not in L2-sense. So for taking the limit L → ∞ one must rely on a distributional wavelet transform, since ψHL→∞ ∈/ L2(R3) but ψHL→∞ ∈ H−4(R3), i.e., it is contained in the dual of the 4th -order Sobolev space. Such technicalities do not occur in the Zernike basis as we will see later in Sect. 4. ψH ψH (0, ·, ·) in mind. Derivations for the results used in the following section can be found in [43]. Fig. 10 Wavelet expanded in the harmonic oscillator basis according to Eq. (81) for L = 15. Left: 3D visualization showing one negative (blue) and one positive (orange) iso-contour. Right: Cross section of the wavelet at x = 0 (Color figure online) ψH (x) = L l=0 √1 rl e− r22 Yl0(θ , φ). l! For this wavelet we have an analytical description in both spatial and Fourier domain. 4.2 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 L2(R3) 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 ∈ L2 (R3), 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 (82) (83) (84) (85) (86) (87) (88) in which (x )α = ΓΓ(x(+x)α) is the (generalized) Pochhammer symbol. Fourier Transform The inverse Fourier transform of the generalized Zernike function (F −1 Znm,,lα)(2π x) = e2πi(ω·x) Rnl,α(ρ) Ylm (ϑ, ϕ)dω ω ≤1 is given by (F −1 Znm,,lα)(2π x) = 4π il Snα,l (2πr ) Ylm (θ , φ), 0 1 with x = r n(θ , φ) and Snα,l (q) = Rnl,α(ρ) jl (qρ) ρ2dρ 4.2.1 The 3D Generalized Zernike Basis The generalized Zernike functions are given by Znm,,lα(ω) = Rnl,α(ρ)Ylm (ϑ, ϕ), with spherical coordinates ω = ρ n(ϑ, ϕ), 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 (81) Rnl,α(ρ) = ρl (1 − ρ2)α Pp =α,nl+−2l21 (2ρ2 − 1), where Pp(α,l+ 12 ) denotes the Jacobi polynomial. The generalized Zernike functions are orthogonal on the unit ball ω ≤1 Nnα,l = Znm11,,l1α(ω)Znm22,,l2α(ω) dω Here Ja and ja are the Bessel functions and spherical Bessel functions [55]. For integer α, the expression in Eq. (89) for q > 0 reduces to 2α(−1) p( p + 1)α jn+qαα++11(q) . 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, F (ω) = cnm,,lα(F ) Znm,,lα(ω), n,l,m we can split the coefficients in radial coefficients and angular coefficients 1 cnm,,lα(F ) = Nnα,l ω ≤1 F (ω) Znm,,lα(ω) dω (1 − ρ2)α 1 = alm ( A) b˜nl,α(B), b˜nl,α(B) = Nnα,l bnl,α(B) where alm ( A) = bnl,α(B) = π 2π 0 0 1 0 A(n(ϑ, ϕ)) Ylm (ϑ, ϕ) sin ϑ dϑ dϕ, B(ρ) Rnl,α(ρ) ρ2dρ . The coefficients cnm,,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, . . . B(ρ) = b˜l,α Rnl,α(ρ). n n=l,l+2,... For each l, the radial functions Rnl,α 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 (96) each l, and it can be shown that there is a recursion formula for the radial coefficients [43]. We now choose appropriate radial and angular functions for our wavelets expressed in the generalized Zernike basis. Angular Function for the Zernike Wavelets For the angular functions we again choose orientation distribution A(n(ϑ, ϕ)) = GsSo2 (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 proce dure 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 = Gsρ (see Eq. ( 37 )). The radial function of ψ0 should therefore approximate B(ρ) = 1 − Gsρ (ρ) 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 Bα,β (ρ) = (1 − ρ2)αρβ , see Fig. 11a for the case α = 6, β = 2. For this function we have the following coefficients bl,α,β n=l+2 p = β−l 2 p (2α + β + l + 2 p + 3) 21 (β+l+1)+α+ p α+ p 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 ρmax = ! , 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 bl,α n=l+2 p for the flattened function as well. For the specific case β = 2 we get the following flattened function Bαfl,a2t(ρ) = Bα,2(ρ) Bαre,c2(ρ) (97) . (98) (99) = Bmax ρ2(1−ρ2)α 1+ (ρ2 −ρm2ax)2 , with Bmax = Bα,2(ρmax). For this flattened function the coefficients are given by bl,flat n=l+2 p = ci bnl,=α,l2++22pi , with ci 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 Bαre,c2(ρ) = ci ρ2i = c0 + c1ρ2 + c2ρ4, ⎛ c0⎞ with ⎝ c1⎠ c2 ⎛ 1 + (α+1)3 ρm4ax⎞ = ⎜ −2 (α+21α)3 ρmax ⎟⎠ . 2 ⎝ (α2+α1)3 2α The filters from this section are summarized in the following result: Result 2 (Analytic 3D-wavelets in Zernike basis) Let α > 0 and let A : S2 → R+ be a function supported mainly in a sharp convex cone around the z-axis and symmetrically around the z-axis. Then A provides our wavelet ψˆ in the Fourier domain via Eq. ( 65 ). The real part of ψ is a tube detector and the imaginary part of ψ is an edge detector, see Fig. 7. We choose radial function g(ρ) = B6fl,a2t( ρρ ) in Eq. ( 18 ) for ψ1 and angular function A(n(ϑ, ϕ)) N= GsSo2 (n(ϑ, ϕ)) and expand in the generalized Zernike basis: ψˆ 1(ω) = cn0,l Rnl,α Yl0(ϑ, ϕ). (103) ρ ρN n − l = 2p, l, n ≥ 0 The coefficients cn0,l follow by expanding A in spherical harmonics and B6fl,a2t in the radial Zernike polynomials, recall Eqs. (95) and (96). This yields al0 (Eq. ( 58 )) and bnα,l (Eq. (101)) and coefficients cn0,l : 0 cn,l = Fig. 11 Left: Function B6,2(ρ) = (1 − ρ2)6ρ2. Right: Flattened function which are obtained from Bα,β (ρ) by multiplying with the second-order Taylor approximation of its reciprocal around the function 3 maximum: B6fl,a2t(ρ) = Bm1ax (1 + 172 (ρ2 − 71 )2)B6,2(ρ) with Snα,l given by Eq. (89) and Ylm the spherical harmonics of Eq. ( 57 ). Then we obtain rotated filters via l ψ1,n(x) = n − l = 2p, m =−l l, n ≥ 0 with (cn)nm,l := (cn(β,γ ))nm,l = cn0,l D0l,m (γ , β, 0). (cn)nm,l 4π il Snα,l (2πr ρN ) Ylm (θ, φ), (106) Since now we do have analytical expressions for the spatial filter, in contrary to the filters from Sect. 3, we sample the filters in the spatial domain using Eq. (106). The filter is a proper wavelet with fast reconstruction property (recall Definition 1). 5 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. 5.1 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. 5.2 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. Cross section of the filter for z = 0. Right: The low-pass filter. Top Filters according to Result 1 with parameters sρ = 21 (1.9)2 and γ = 0.85. Bottom: The filters according to Result 2 with α = 3 and β = 2. Both have so = 21 (0.4)2 and are evaluated on a grid of 31 × 31 × 31 voxels (Color figure online) We see no qualitative differences between the two reconstructions. 6 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. 6.1.1 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 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. 6.1.2 CEDOS We now use the invertible orientation score transformation to perform data enhancement according to Fig. 2. Because R3 × S2 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 U˜ (x, R) = U (x, R · ez ), U (x, n) = U˜ (x, Rn), (107) with Rn any rotation for which Rn · ez = n. The operators Φ which we consider are scale spaces on SE( 3 ) (diffusions), and are given by Φ = Φt with 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 Φt (U )(y, n) = W˜ (y, Rn, t ). Here W˜ is the solution of a nonlinear diffusion equation: ⎧ ⎪⎪⎨ ⎪⎩⎪ ∂ W˜ ∂ t 6 i, j=1 (g, t ) = Ai |g Di j (U˜ )A j |g W˜ (g, t ), W˜ (g, 0) = Wψ1 [ f ](x, R · ez ), g = (x, R) (108) (109) where in coherence-enhancing diffusion on orientation scores (CEDOS) Di j is adapted locally to initial condition W˜ (g, 0) based on exponential curve fit (see [32]), and with Ai |g = (L g )∗Ai |e the left-invariant vector fields on SE( 3 ). The diffusion is better understood in locally adaptive frame {Bi }i6=1. Here B3 follows from an exponential curve fits and points along our structure. B1 and B2 span the plane spatially perpendicular to our structure, and B4, B5 and B6 correspond to angular diffusion (embedding in S E ( 3 ) leads to a third angular dimension). Our diffusion then takes the diagonal form Defining Eq. ( 39 ) ( 51 ) ( 50 ) ( 18 ) ( 52 ) – (112) where we limit ourselves to diffusion of type D11 = D22, and D44 = D55 = D66 to preserve the data symmetry of Eq. (107). The diffusion system in Eq. (110) has been set up in previous work [45] with constant D11 and D33. In this work we include further adaptivity by making the constants D11 and D33 data adaptive. We aim to enhance oriented structures and reduce noise as much as possible. Therefore, non-oriented regions are smoothed isotropically by setting ( D11(U˜ ))(g) = 1 − exp c1 s(U˜ )(g) 2" , (111) Since we want to stop diffusion when reaching the end of a structure, we set Table 1 Default values for the parameters of the orientation score transform used in the application section Parameter No γ σer f sρ so Discrete wavelet size Default value with c1 a constant automatically set to the 50% quantile of s(U˜ ) and s(U˜ )(g) a measure for orientation confidence given by minus the Laplacian in the space orthogonal to the structure orientation B3: s(U˜ )(g) = − Bi |2gU˜ (g). i∈{1,2,4,5,6} Fig. 15 Selected regions for determining the contrast-to-noise ratios for dataset B.1. Left: Grid containing slices of the data (top row), the same slices with segmented vessel parts (second row) and the slices with three selected background regions (third row), the slices after applying CEDOS. Right: 3D visualization of the data ( D33(U˜ ))(g) = 1 − exp ⎝ − c2 B3|gU˜ (g) "2⎞ ⎠ , with c2 a constant automatically set to the 50% quantile of |B3|g (U˜ )|. We then obtain Euclidean invariant image processing via Υ f = Wψ−1,ext ◦ Φ ◦ Wψ f = Wψ−1 ◦ Pψ Φ ◦ Wψ f , (114) 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 L2(R3 × S2). 6.1.3 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: (113) Given such data, the noise is quantified via the contrast-tonoise ratio (CNR): C N R( f N , f ) = max f (x) − min f (x) x x σ ( f − f N ) , (116) 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 : μS = f N |ΩS , μB = f N |ΩB . (117) 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 f N (x) = f (x) + N (x). (115) σN = σ ( f N |ΩB ). (118) The contrast-to-noise ratio (CNR) is then given by C N R( f N ) = μS − μB σN , (119) where the numerator denotes the contrast in our data. 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 D11 and D33 were determined using Eqs. (111), (113), and we set D11 = 0.001. For the orientation score transformation, we used so = 21 (0.45)2 and sρ = 21 ( 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). 6.1.5 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 8 6 4 2 1 2 3 4 5 6 Fig. 16 Contrast-to-noise ratio (CNR) against diffusion time for CEDOS compared to Gaussian regularization and CED of data B.1 depicted in Fig. 15. The times denoted by a, b and c correspond to the diffusion times shown in 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 D11 (Eq. 111) and therefore we smooth only along the vessel. 6.2 Tubularity Measure 6.2.1 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. 6.2.2 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 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 ∈ R3 is given by Eprod(x, n, r , θ ) = Im+[U (x + r n⊥(θ ), n⊥(θ ))] · Im+[U (x − r n⊥(θ ), −n⊥(θ ))], (120) where Im+(z) = max{0, Im(z)} and n⊥(θ ) = cos θ e1 + sin θ e2, with {e1, e2} 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 V (x, n, r ) = θ∈[0,2π) min 0 π KσorV (θ − θ ) Eprod(x, n, r , θ )dθ . o (121) For the angular regularization kernel we use the simple 2π periodic 1D-diffusion kernel σ V G So1 (θ + 2kπ ) = KσoorV (θ ) = k∈Z k∈Z √21πσov e over θ , instead of the true diffusion kernel (72), where for − |θ2+|σ2oVkπ|2|2 practical values we use σoV = π/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: st (x) = max n∈S2,r∈R+ V (x, n, r ), n∗(x) = arg max max V (x, n, r ), n∈S2 r∈R+ r ∗(x) = arg max max V (x, n, r ). r∈R+ n∈S2 (122) (123) (124) Here st (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. vessel width is not influenced by our regularization method while this does occur for Gaussian regularization. Bottom: Cross sections of one vessel for increasing diffusion time with detected vessel edge positions (green points) and search area for edge detection (red circle) (Color figure online) Fig. 20 Schematic visualization of the edges used in the tubularity measure V (x, n, r ). Left: A 3D iso-contour visualization of a vessel with orientation n. The coordinates (θ, r ) are polar coordinates for the plane perpendicular to orientation n spanned by e1 and e2. Right: In this plane opposite edges are multiplied in Eprod(x, n, r , θ ). The edge is expected to have outward orientation n⊥(θ ) 6.2.3 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: σoV = π/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 Gausof tubularity over radius and orientation) with radius of max response r ∗(x) in color. Right: Measured radius r ∗(x) against ground truth radius rGT(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 rGT (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. 6.2.4 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: σoV = π/8 and we discretized the θ -integral in Eq. (121) using 12 orientations. For the orientation score transformation we used so = 21 (0.4)2, sρ = 21 ( 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 distance map d(x, ∪ Bci ,r ∗(ci )) = minci { x − ci − r ∗(ci )}, where ci are the positions given by the 1% quantile of highest responses. 7 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. Fig. 23 Tubularity on real data. Top: Dataset A.1. Bottom: Dataset A.2. Left: The data. Middle left: The projected tubularity (max over radius and orientation) colored by radius. Middle right: Orientation of max response for 1% quantile of highest responses in the tubularity confidence. Right: Segmentation based on radius of max response for 1% quantile of highest responses in the tubularity confidence. The plotted surface is the 0-iso-contour of the distance map d(x, ∪Bci ,r∗(ci )) = minci { x − ci − r ∗(ci )}, where ci are the positions given by the 1% quantile of highest responses 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]. Acknowledgements We would like to thank professors Jacques Moret and Laurent Spelle from the NEURI center, Department of NRI, Bicêtre University Hospital Paris, France, for providing the AVM dataset used in this article. The research leading to the results of this article has received funding from the European Research Council (ERC) under the European Community 7th Framework Programme (FP7/20072014)/ERC Grant Agreement No. 335555 (Lie Analysis). A Steerable Orientation Score Transform In this article we often rely on a spherical harmonic decomposition of the angular part of proper wavelets in the spatial and the Fourier domain. As the choice of radial basis varies in Sect. 4, and since it does not affect steerable filter [36,38,61,67] properties, we simply write {gnl (r )}n∞=0 for a radial orthonormal basis for L2(R+, r 2dr ) associated to each 2l + 1-dimensional S O ( 3 )-irreducible subspace indexed by l. Then by the celebrated Bochner–Hecke theorem [34], this induces a corresponding orthonormal radial basis {g˜nl (ρ)}n∞=0 in the Fourier domain which can be obtained by a Hankeltype of transform [25, Ch. 7]. We expand our wavelets in spherical harmonics Ylm and ball coordinates (cf. (78)) accordingly: αln gnl (r ) Yl0(θ , φ), αln g˜nl (ρ) Yl0(ϑ, ϕ), n=0 with ω = (ρ cos ϕ sin ϑ, ρ cos ϕ sin ϑ, ρ cos ϑ ). Remark 4 In Sect. 4, we consider the modified Zernike basis in which case gnl and g˜nl are given by, respectively, (85) and (82), whereas for the harmonic oscillator basis one has gnl = i l (−1)n+l g˜nl given, respectively, by (77) and (79). We obtain steerability via finite series truncation at n = N and l = L . Then wel rotate the steerable kernels via the Wigner-D functions D0,m (γ , β, 0) ∈ R (cf. Sect. 3.2.3) and one obtains the following steerable implementations of orientation scores: (125) U (x, n) N L l n=N0 l=L0 m=l −l n=0 l=0 m=−l αln D0l,m (γ , β, 0) · ((gnl ⊗ Yl0) f )(x) αln D0l,m (γ , β, 0) · F −1 )g˜nl ⊗ Yl0 · fˆ* (x) (126) where n = (cos γ sin β, sin γ sin β, cos β)T , denotes correlation, the overline denotes complex conjugation and with function product (g˜nl ⊗ Ylm )(x) = g˜nl (r ) Ylm (θ , φ). B Table of Notations Symbol B.1 Spaces and Input Data R3 × S2 L2 (R3) = { f ∈ L2(R3)|supp(F f ) ⊂ B }, with > 0 Explanation Space of positions and orientations Space of frequency ball-limited images The ball around the origin with radius ρ Orientation score transformation Data reconstruction via the inverse orientation score transformation Wavelet used for the orientation score transformation Fourier transform of the wavelet used for the orientation score transformation Factor used to quantify the stability of the transformation Factor used to quantify the stability of the transformation when using the simplified reconstruction by integration Discretized versions of Wψ , (Wψ )−1, Mψ , Nψ Spherical area measure and discretized spherical area measure Cartesian coordinates real space Cartesian coordinates Fourier domain Spherical coordinates real space Spherical coordinates Fourier domain Radial function of the cake filters Orientation distribution used in wavelet construction Spherical harmonics 3D generalized Zernike functions Radial part of the 3D generalized Zernike functions Radial part of the inverse Fourier transform 3D generalized Zernike functions Tubularity measure Tubularity confidence Orientation of maximum tubularity response Radius of maximum tubularity response ( 39 ), ( 40 ), ( 41 ), ( 42 ) ( 5 ), ( 40 ) Page 1 ( 2 ) ( 2 ) ( 3 ) ( 5 ) ( 3 ) ( 6 ) ( 6 ) ( 12 ) ( 49 ) ( 49 ) ( 50 ) ( 52 ) ( 57 ) (82) (84) (89) (121) (122) (122) (122) E. J. Bekkers received his M.Sc. degree in Biomedical Engineering in 2012 at Eindhoven University of Technology (TU/e), the Netherlands. In January 2017 he received his Ph.D. degree (cum laude) for his thesis "Retinal Image Analysis using Sub-Rieman nian Geometry in SE( 2 )," which was conducted in a combined position at Biomedical Engineering at TU/e and at the medical device company i-Optics, the Hague, the Netherlands. He currently is a postdoc researcher at the Department of Mathematics and Computer Science at TU/e, working on the application of Lie group (and control) theory in medical image analysis. His research interests include differential geometry, PDEs, variational methods, machine learning, medical imaging and statistical biomarker analysis. J. Oliván Bescós (1976) received his Physics degree from University of Zaragoza, Spain, in 1999. In 2004 he finished his Ph.D. in Electrical Engineering at the University of Twente, The Netherlands, where he conducted research on medical visualization and image processing. Since 2004 he works at Philips Healthcare, where he participates in the design, development and validation of algorithms and medical software to support doctors in imaging-guided interventions. He has been part of several multi-disciplinary projects acting as a liaison between clinicians, image processing and visualization scientists and engineering teams. R. Duits received his M.Sc. degree (cum laude) in Mathematics in 2001 at the TU/e, The Netherlands. He received his Ph.D. degree (cum laude) at the Department of Biomedical Engineering at the TU/e. Now he is associate professor at the Department of Applied Mathematics & Computer Science and affiliated part time at the Biomedical Image Analysis group at the Biomedical Engineering Department. He has received an ERC-StG 2013 personal grant (Lie analysis, No. 335555). His research interests subtend group theory, functional analysis, differential geometry, PDE’s, harmonic analysis, geometric control and their applications to biomedical imaging and vision. 1. Aarts , R.M. , Janssen , A.J.E.M. : On-axis and far-field sound radiation from resilient flat and dome-shaped radiators . J. Acoust. Soc. Am . 125 ( 3 ), 1444 - 1455 ( 2009 ) 2. Ali , S.T.: A general theorem on square-integrability: vector coherent states . J. Math. Phys . 39 ( 8 ), 3954 ( 1998 ) 3. Ali , S.T. , Antoine , J.P. , Gazeau , J.P. : Coherent States, Wavelets, and Their Generalizations. Graduate Texts in Contemporary Physics . Springer, New York ( 2000 ) 4. Antoine , J.P. , Murenzi , R. , Vandergheynst , P. : Directional wavelets revisited: Cauchy wavelets and symmetry detection in patterns . Appl. Comput. Harmon. Anal. 6 , 314 - 345 ( 1999 ) 5. Barbieri , D. , Citti , G.: Reproducing kernel Hilbert spaces of CR functions for the Euclidean motion group . Anal. Appl . 13 ( 03 ), 331 - 346 ( 2014 ) 6. Barbieri , D. , Citti , G. , Cocci , G. , Sarti , A. : A cortical-inspired geometry for contour perception and motion integration . J. Math. Imaging Vis . 49 ( 3 ), 511 - 529 ( 2014 ) 7. Batard , T. , Sochen , N.: A class of generalized laplacians on vector bundles devoted to multi-channel image processing . J. Math. Imaging Vis . 48 ( 3 ), 517 - 543 ( 2014 ) 8. Bekkers , E. , Duits , R. , Mashtakov , A. , Sanguinetti , G.: A PDE approach to data-driven sub-Riemannian geodesics in SE(2) . SIAM J. Imaging Sci . 8 ( 4 ), 2740 - 2770 ( 2015 ) 9. Bekkers , E.J. : Retinal Image Analysis Using Sub-Riemannian Geometry in SE(2) . Ph.D. Thesis , University of Technology, Eindhoven ( 2017 ) 10. Bekkers , E.J. , Duits , R. , Berendschot , T., ter Haar Romeny, B.M.: A multi-orientation analysis approach to retinal vessel tracking . J. Math. Imaging Vis . 49 ( 3 ), 583 - 610 ( 2014 ) 11. Bekkers , E.J. , Zhang, J., Duits , R., ter Haar Romeny, B.M. : Curvature based biomarkers for diabetic retinopathy via exponential curve fits in SE(2) . In: Proceedings of the Ophthalmic Medical Image Analysis First International Workshop , pp. 113 - 120 ( 2015 ) 12. Born , M. , Wolf , E.: Principles of Optics, 7th edn . Cambridge University Press, Cambridge ( 1999 ) 13. Boscain , U. , Chertovskih , R.A. , Gauthier, J.P. , Remizov , A.O. : Hypoelliptic diffusion and human vision: a semidiscrete new twist . SIAM J. Imaging Sci . 7 ( 2 ), 669 - 695 ( 2014 ) 14. Bouix , S. , Siddiqi , K. , Tannenbaum , A. : Flux driven automatic centerline extraction . Med . Image Anal . 9 ( 3 ), 209 - 221 ( 2005 ) 15. Brunner , T.A. : Impact of lens aberrations on optical lithography . IBM J. Res. Dev . 41 ( 1 .2), 57 - 67 ( 1997 ) 16. Burgeth , B. , Didas , S. , Weickert , J.: A general structure tensor concept and coherence-enhancing diffusion filtering for matrix fields . In: Laidlaw, D. , Weickert , J . (eds.) Visualization and Processing of Tensor Fields, Mathematics and Visualization , pp. 305 - 323 . Springer, Berlin ( 2009 ) 17. Burgeth , B. , Pizarro , L. , Didas , S. , Weickert , J.: 3D-coherenceenhancing diffusion filtering for matrix fields . In: Florack, L. , Duits , R. , Jongbloed , G., van Lieshout, M.C. , Davies , L . (eds.) Mathematical Methods for Signal and Image Analysis and Representation , Computational Imaging and Vision , pp. 49 - 63 . Springer, London ( 2012 ) 18. Caruyer , E. , Lenglet , C. , Sapiro , G. , Deriche , R.: Design of multishell sampling schemes with uniform coverage in diffusion MRI . Magn. Reson. Med . 69 ( 6 ), 1534 - 1540 ( 2013 ) 19. Chen , D. , Cohen , L.D. : Automatic vessel tree structure extraction by growing minimal paths and a mask . In: 2014 IEEE 11th International Symposium on Biomedical Imaging (ISBI) , pp. 802 - 805 ( 2014 ) 20. Chung , M.K. : Heat kernel smoothing on unit sphere . In: 3rd IEEE International Symposium on Biomedical Imaging: Nano to Macro , 2006 , pp. 992 - 995 ( 2006 ) 21. Citti , G. , Franceschiello , B. , Sanguinetti , G. , Sarti , A. : SubRiemannian mean curvature flow for image processing . SIAM J. Imaging Sci . 9 ( 1 ), 212 - 237 ( 2016 ) 22. Citti , G. , Sarti , A. : A cortical based model of perceptual completion in the roto-translation space . J. Math. Imaging Vis . 24 ( 3 ), 307 - 326 ( 2006 ) 23. Descoteaux , M. , Angelino , E. , Fitzgibbons , S. , Deriche , R.: Regularized, fast, and robust analytical Q-ball imaging . Magn. Reson. Med . 58 ( 3 ), 497 - 510 ( 2007 ) 24. Dirksen , P. , Braat , J.J.M. , Janssen , A.J.E.M. , Juffermans , C.A.H. : Aberration retrieval using the extended Nijboer-Zernike approach . J. Microlithogr. Microfabr. Microsyst . 2 , 61 - 67 ( 2003 ) 25. Duits , R.: Perceptual Organization in Image Analysis . Ph.D. Thesis , Technische Universiteit Eindhoven ( 2005 ) 26. Duits , R. , Burgeth , B. : Scale spaces on lie groups . In: Scale Space and Variational Methods in Computer Vision , pp. 300 - 312 . Springer, Berlin ( 2007 ) 27. Duits , R. , Duits , M., van Almsick , M. , ter Haar Romeny, B.M.: Invertible orientation scores as an application of generalized wavelet theory . PRIA 17 ( 1 ), 42 - 75 ( 2007 ) 28. Duits , R. , Felsberg , M. , Granlund , G., ter Haar Romeny, B.M. : Image analysis and reconstruction using a wavelet transform constructed from a reducible representation of the euclidean motion group . Int. J. Comput. Vis . 72 ( 1 ), 79 - 102 ( 2007 ) 29. Duits , R. , Franken , E.: Left-invariant diffusions on the space of positions and orientations and their application to crossingpreserving smoothing of HARDI images . Int. J. Comput. Vis . 92 ( 3 ), 231 - 264 ( 2010 ) 30. Duits , R. , Franken , E.: Left-invariant parabolic evolutions on SE(2) and contour enhancement via invertible orientation scores Part I: linear left-invariant diffusion equations on SE(2). Q . Appl. Math. 68 ( 2 ), 255 - 292 ( 2010 ) 31. Duits , R. , Führ , H. , Janssen , B. , Bruurmijn , M. , Florack , L., van Assen , H.: Evolution equations on Gabor transforms and their applications . Appl. Comput. Harmon. Anal . 35 ( 3 ), 483 - 526 ( 2013 ) 32. Duits , R. , Janssen , M.H.J. , Hannink , J. , Sanguinetti , G.R. : Locally adaptive frames in the roto-translation group and their applications in medical imaging . J. Math. Imaging Vis . 56 ( 3 ), 367 - 402 ( 2016 ) 33. Dunkl , C.F. : Orthogonal Polynomials of Several Variables . Cambridge University Press, Cambridge ( 2014 ) 34. Faraut , J. , Harzallah , K. : Deux Cours d'Analyse Harmonique . Birkhauser, Basel ( 1987 ) 35. Felsberg , M. : Adaptive filtering using channel representations . In: Florack, L.M.J. , Duit , R. , Jongbloed , G., van Lieshout, M.C. , Davies , L . (eds.) Mathematical Methods for Signal and Image Analysis and Representation , vol. 41 , pp. 31 - 48 ( 2012 ) 36. Franken , E.: Enhancement of Crossing Elongated Structures in Images . Ph.D. Thesis , Technische Universiteit Eindhoven ( 2008 ) 37. Franken , E. , Duits , R.: Crossing-preserving coherence-enhancing diffusion on invertible orientation scores . Int. J. Comput. Vis . 85 ( 3 ), 253 - 278 ( 2009 ) 38. Freeman , W.T., Adelson , E.H. : The design and use of steerable filters . IEEE Trans. Pattern Anal. Mach. Intell . 13 ( 9 ), 891 - 906 ( 1991 ) 39. Führ , H.: Abstract Harmonic Analysis of Continuous Wavelet Transforms . No. 1863 in Lecture Notes in Mathematics. Springer , Berlin ( 2005 ) 40. Gräf , M. , Potts , D. : On the computation of spherical designs by a new optimization approach based on fast spherical fourier transforms . Numer . Math. 119 ( 4 ), 699 - 724 ( 2011 ) 41. Griffiths , D.J. : Introduction to Quantum Mechanics, 2nd edn . Cambridge University Press, Cambridge ( 2016 ) 42. Isham , C.J. , Klauder , J.R. : Coherent states for n-dimensional Euclidean groups E(n) and their application . J. Math. Phys . 32 ( 3 ), 607 - 620 ( 1991 ) 43. Janssen , A.J. : Generalized 3d Zernike Functions for Analytic Construction of Band-Limited Line-Detecting Wavelets . ArXiv preprint arXiv:1510.04837 ( 2015 ) 44. Janssen , M.H.J. , Duits , R. , Breeuwer , M. : Invertible orientation scores of 3D images . In: Scale Space and Variational Methods in Computer Vision, Lecture Notes in Computer Science , vol. 9087 , pp. 563 - 575 . Springer ( 2015 ) 45. Janssen , M.H.J. , Haije , T.C.J.D. , Martin , F.C. , Bekkers , E.J. , Duits , R. : The Hessian of axially symmetric functions on SE(3) and application in 3D image analysis . In: Scale Space and Variational Methods in Computer Vision , pp. 643 - 655 . Springer, Cham ( 2017 ) 46. Kalitzin , S.N. , ter Haar Romeny, B.M. , Viergever , M.A. : Invertible apertured orientation filters in image analysis . Int. J. Comput. Vis . 31 ( 2-3 ), 145 - 158 ( 1999 ) 47. Law , M.W.K. , Chung , A.C.S.: Three dimensional curvilinear structure detection using optimally oriented flux . In: Computer Vision-ECCV 2008 , pp. 368 - 382 . Springer, Berlin ( 2008 ) 48. Lee , T.S. : Image representation using 2D Gabor wavelets . IEEE Trans. Pattern Anal. Mach. Intell . 18 ( 10 ), 959 - 971 ( 1996 ) 49. Louis , A. : Orthogonal function series expansions and the null space of the Radon transform . SIAM J. Math. Anal . 15 ( 3 ), 621 - 633 ( 1984 ) 50. Louis , A. , Maass , D. , Rieder , A. : Wavelets: Theory and Application . Wiley, New York ( 1997 ) 51. Mallat , S.: A Wavelet Tour of Signal Processing . Academic Press, London ( 1999 ) 52. Mashtakov , A. , Duits , R. , Sachkov , Y. , Bekkers , E.J. , Beschastnyi , I. : Tracking of lines in spherical images via sub-Riemannian geodesics in SO(3) . J. Math. Imaging Vis . 58 ( 2 ), 239 - 264 ( 2017 ) 53. Mashtakov , A.P. , Ardentov , A.A. , Sachkov , Y.L. : Parallel algorithm and software for image inpainting via sub-Riemannian minimizers on the group of rototranslations . Numer. Math. Theory Methods Appl . 6 ( 01 ), 95 - 115 ( 2013 ) 54. Muhlich , M. , Aach , T. : Analysis of multiple orientations . IEEE Trans. Image Process . 18 ( 7 ), 1424 - 1437 ( 2009 ) 55. Olver , F.W.J. : NIST Handbook of Mathematical Functions . Cambridge University Press, Cambridge ( 2010 ) 56. Pechaud , M. , Keriven , R. , Peyre , G.: Extraction of tubular structures over an orientation domain . In: 2009 IEEE Conference on Computer Vision and Pattern Recognition , pp. 336 - 342 ( 2009 ) 57. Petitot , J.: The neurogeometry of pinwheels as a sub-Riemannian contact structure . J. Physiol . 97 ( 2-3 ), 265 - 309 ( 2003 ) 58. Portegies , J.M. , Fick , R.H.J. , Sanguinetti , G.R. , Meesters , S.P.L. , Girard , G. , Duits , R.: Improving fiber alignment in hardi by combining contextual pde flow with constrained spherical deconvolution . PLoS ONE 10 ( 10 ), 1 - 33 ( 2015 ) 59. Prandi , D. , Boscain , U. , Gauthier, J.P. : Image processing in the semidiscrete group of rototranslations . In: Geometric Science of Information , pp. 627 - 634 . Springer, Cham ( 2015 ) 60. Rajan , S. , Das , T. , Krishnakumar , R.: An analytical method for the detection of exudates in retinal images using invertible orientation scores . In: Proceedings of the World Congress on Engineering , vol. 1 ( 2016 ) 61. Reisert , M. : Group Integration Techniques in Pattern Analysis-A Kernel View . Ph.D. Thesis , Albert- Ludwigs- University ( 2008 ) 62. Ruijters , D. , Vilanova , A. : Optimizing GPU volume rendering . J. WSCG 14 ( 1-3 ), 9 - 16 ( 2006 ) 63. Scharr , H.: Diffusion-like reconstruction schemes from linear data models . In: Pattern Recognition , pp. 51 - 60 . Springer, Berlin ( 2006 ) 64. Sharma , U. , Duits , R.: Left-invariant evolutions of wavelet transforms on the similitude group . Appl. Comput. Harmon. Anal . 39 ( 1 ), 110 - 137 ( 2015 ) 65. Sifre , L. , Mallat , P.S. : Rigid-Motion Scattering For Image Classification . Ph.D. Thesis , Ecole Polytechnique, CMAP, Paris ( 2014 ) 66. Steidl , G. , Teuber , T. : Anisotropic smoothing using double orientations . In: Scale Space and Variational Methods in Computer Vision , pp. 477 - 489 ( 2009 ) 67. van Almsick, M. : Context Models of Lines and Contours . Ph.D. Thesis , Technische Universiteit Eindhoven ( 2007 ) 68. Vogt , T. , Lellmann , J.: Measure-Valued Variational Models with Applications to Diffusion-Weighted Imaging . ArXiv preprint arXiv:1710.00798 ( 2017 ) 69. Weickert , J.: Coherence-enhancing diffusion filtering . Int. J. Comput. Vis . 31 ( 2-3 ), 111 - 127 ( 1999 ) 70. Zhang , J., Dashtbozorg , B. , Bekkers , E. , Pluim , J.P.W. , Duits , R., ter Haar Romeny, B.M.: Robust retinal vessel segmentation via locally adaptive derivative frames in orientation scores . IEEE Trans. Med. Imaging 35 ( 12 ), 2631 - 2644 ( 2016 ) 71. Zweck , J. , Williams , L.R. : Euclidean group invariant computation of stochastic completion fields using shiftable-twistable functions . J. Math. Imaging Vis . 21 ( 2 ), 135 - 154 ( 2004 ) M. H. J. Janssen received his M.Sc. degree (cum laude) in Biomedical Engineering in 2014 at TU/e, The Netherlands . Afterward, he worked as a Ph.D. student in the Department of Applied Mathematics and Computer Science at TU/e, working on the topic of 3D multiorientation analysis for medical image processing. His research interests include geometrical methods for medical image analysis, multi-orientation analysis and differential equations on Lie groups . Currently, he is working for Vital- Health software, developing e-Health software applications . A. J. E. M. Janssen ( 1953 ) received his Engineering and Ph.D. degree in Mathematics from the Eindhoven University of Technology, The Netherlands , in 1976 and 1979, respectively. From 1979 to 1981 he was a Bateman Research Instructor at the Mathematics Department of California Institute of Technology. From 1981 until 2010 he worked with Philips Research Laboratories, Eind hoven, where he was a Research Fellow since 1999 and recipient of the Gilles Holst award in 2003. His principal responsibility at Philips Research was to provide highlevel mathematical service and consultancy in mathematical analysis. He continues his research and consultancy activities from his current affiliations TU/e, EURANDOM, with further contacts with TU-Delft (Optics) and ASML . In 2003 , he was elected Fellow of the IEEE for his contributions to mathematical time-frequency analysis and signal processing . Dr. Janssen has published over 200 papers in the fields of Fourier analysis, signal analysis, time-frequency analysis, electron microscopy, optical diffraction theory, acoustics and queueing theory .


This is a preview of a remote PDF: https://link.springer.com/content/pdf/10.1007%2Fs10851-018-0806-0.pdf

M. H. J. Janssen, A. J. E. M. Janssen, E. J. Bekkers, J. Oliván Bescós, R. Duits. Design and Processing of Invertible Orientation Scores of 3D Images, Journal of Mathematical Imaging and Vision, 2018, 1-32, DOI: 10.1007/s10851-018-0806-0