Abstract
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.
Similar content being viewed by others
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:{\mathbb {R}}^ 3\rightarrow {\mathbb {R}}\) with a set of oriented filters to construct a 3D orientation score \(U:{\mathbb {R}}^ 3\times S ^ 2 \rightarrow {\mathbb {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 coherence-enhancing 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.
2D orientation score for an exemplary image. In the orientation score crossing structures are disentangled because the different structures have a different orientation
A schematic view of image processing via invertible orientation scores
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:{\mathbb {R}}^ 2 \times S^1 \rightarrow {\mathbb {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.
Overview of applications of processing via orientation scores. Top: We reduce noise in real medical image data (3D rotational Xray) of the abdomen containing renal arteries by applying diffusions via 3D orientation 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 basisfollowed 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 coherence-enhancing 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.
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.
Construction of a 3D orientation score. Top: The data f are correlated with an oriented filter \(\psi _{{\mathbf {e}}_x}\) to detect structures aligned with the filter orientation \({\mathbf {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
2 Invertible Orientation Scores
2.1 Continuous Orientation Score Transform
Throughout this article we use the following definition of the Fourier transform on \({\mathbb {R}}^3\):
An invertible orientation score \({\mathcal {W}}_{\psi }[f]:{\mathbb {R}}^3 \times S^2 \rightarrow {\mathbb {C}}\) is constructed from a given ball-limited 3D dataset
with ball \(B_{\varrho } = \{\varvec{\omega }\in {\mathbb {R}}^3 \;|\; \Vert \varvec{\omega }\Vert < \varrho \}\) of radius \(\varrho >0\), by correlation \(\star \) with an anisotropic kernel
This is illustrated in Fig. 4. Here \(\psi \in {\mathbb {L}}_2({\mathbb {R}}^3) \cap {\mathbb {L}}_1({\mathbb {R}}^3)\) is a wavelet aligned with and rotationally symmetric around the z-axis, and \(\psi _{{\mathbf {n}}}\in {\mathbb {L}}_2({\mathbb {R}}^3)\) the rotated wavelet aligned with \({\mathbf {n}}\) given by
Here \({\mathbf {R}}_{{\mathbf {n}}} \in SO(3)\) is any 3D rotation which rotates the z-axis onto \({\mathbf {n}}\) where the specific choice of rotation does not matter because of the rotational symmetry of \(\psi \). The overline denotes a complex conjugate. The exact reconstruction formula [25] for this transformation is given by
with \({\check{\psi }}_{\mathbf {n}}({\mathbf {x}}) =\psi _{\mathbf {n}}(-{\mathbf {x}})\). The function \(M_\psi : {\mathbb {R}}^3 \rightarrow {\mathbb {R}}^+ \) is given by
and vanishes at \(\infty \), where the circumflex \((\,\hat{}\,)\) again denotes Fourier transformation. Due to our restriction to ball-limited data (2), this does not cause problems in reconstruction (5). The function \(M_\psi \) quantifies the stability of the inverse transformation [25], since \(M_\psi (\varvec{\omega })\) specifies how well frequency component \(\varvec{\omega }\) is preserved by the cascade of construction and reconstruction when \(M_\psi ^{-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_\psi (\varvec{\omega }) \approx 1,\) in view of the condition number of transformation \({\mathcal {W}}_\psi :{\mathbb {L}}_2^\varrho ({\mathbb {R}}^3)\rightarrow {\mathbb {L}}_2^\varrho ({\mathbb {R}}^3\times S^2)\) given by:
where M and \(\delta \) are assumed to be tight bounds in (7). In the codomain spatial frequencies are again limited to the ball:
Also, in the case we have \(M_\psi (\varvec{\omega }) = 1\) for \(\varvec{\omega }\in B_{\varrho }\) we have \({\mathbb {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 \in {\mathbb {L}}_{2}^{\varrho }({\mathbb {R}}^3)\) to an orientation integrated score
Its condition number is given by \({{\mathrm{cond}}}(A_\psi ) = \frac{\max \limits _{\varvec{\omega }\in B_{\varrho }}N_{\psi }(\varvec{\omega })}{\min \limits _{\varvec{\omega }\in B_{\varrho }}N_{\psi }(\varvec{\omega })} \).
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_\psi \) 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_\psi \) leads to restrictions on our wavelets which are captured in the following definition:
Definition 1
(Proper wavelet) Let us set a priori boundsFootnote 1 \(\delta ,M>0,\ 0 < \varepsilon \ll 1\). Furthermore, let \(\varrho \) be an a priori maximum frequency of our ball-limited image data. Then, a wavelet \(\psi \in {\mathbb {L}}_{2}({\mathbb {R}}^{3}) \cap {\mathbb {L}}_{1}({\mathbb {R}}^{3})\) is called a proper wavelet if
where \({\mathbf {R}}_{{\mathbf {e}}_z, \alpha } \in SO(3)\) is the 3D rotation around axis \({\mathbf {e}}_z\) over angle \(\alpha \).
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.
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
with Gaussian window in the Fourier domain given by
After splitting, \(\psi _0\) contains the average and low-frequency components and \(\psi _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 \(\psi _{i,{\mathbf {n}}}({\mathbf {x}})=\psi _{i}({\mathbf {R}}_{{\mathbf {n}}}^T {\mathbf {x}})\), as in Eq. (4). The vector transformation is then defined as
For this transformation we have the exact reconstruction formula
with
Again, \(M_{{\underline{\psi }}}\) 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 \(\psi \in {\mathbb {L}}_2({\mathbb {R}}^3) \cap {\mathbb {L}}_1({\mathbb {R}}^3)\) such that Eq. (7) holds, \(\delta =\min _{\varvec{\omega }\in B_\varrho } M_\psi (\varvec{\omega })\) and \(M=\max _{\varvec{\omega }\in B_\varrho } M_\psi (\varvec{\omega })\). Then the condition number of \({\mathcal {W}}_\psi : {\mathbb {L}}_2^\varrho ({\mathbb {R}}^3) \rightarrow {\mathbb {L}}_2^\varrho ({\mathbb {R}}^3 \times S^2)\) is given by
The condition number of \({\underline{{\mathcal {W}}}}_{{\underline{\psi }}}: {\mathbb {L}}_2^\varrho ({\mathbb {R}}^3) \rightarrow {\underline{{\mathbb {L}}}}_2^\varrho ({\mathbb {R}}^3 \times S^2)\) obtained from \({\mathcal {W}}_\psi \) 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 \({\mathcal {W}}_\psi \):
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_{\psi }\). Similarly, we get \(\max _{\varvec{\omega }\in B_\varrho } M_{\psi }(\varvec{\omega })\) for the second factor in Eq. (27). Then we obtain
Similarly the condition number of \({\underline{{\mathcal {W}}}}_{{\underline{\psi }}}\) is given by
Next we express \(M_{{\underline{\psi }}}\) in \(M_{\psi }\) of the original wavelet as
So it remains to quantify \(I(\varvec{\omega })\). For a wavelet splitting according to (18) we have
Hence
And since \(\frac{1}{2} \le 1-2x(1-x)) \le 1\) for \(0 \le x \le 1\) we have for \(M_\psi \) satisfying (7) the following bounds on \(M_{{\underline{\psi }}}\):
thereby guaranteeing stability after the splitting (18). As we cannot guarantee that \(\delta /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). \(\square \)
For this vector transformation we can also use the approximate reconstruction by integration (for \(N_\psi \approx 1\)) over orientations. Thus we have
As said we are only interested in processing of \({\mathcal {W}}_{\psi _1}[f]\) and not in processing of \({\mathcal {W}}_{\psi _0}[f]\), and so we directly calculate \(L_{\psi _0}[f]\) via
For a design with \(N_\psi (\varvec{\omega })=1\) for all \(\varvec{\omega }\in B_\varrho \), we have \({\hat{\phi }}_0={\hat{G}}_{s_{\rho }}\) and so
Then Eq. (35) becomes
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 simpleFootnote 2 electrostatic repulsion model [18].
Assume we have a number \(N_o\) of orientations \({\mathcal {V}}=\{{\mathbf {n}}_1,{\mathbf {n}}_2,\ldots ,{\mathbf {n}}_{N_o}\}\subset S^2\), and define the discrete invertible orientation score \({\mathcal {W}}_\psi ^d[f]:{\mathbb {R}}^3\times {\mathcal {V}}\rightarrow {\mathbb {C}}\) by
The exact reconstruction formula is in the discrete setting given by
with \({\varDelta }_i\) the discrete spherical area measure (\(\sum \nolimits _{i=1}^{N_o} {\varDelta }_i =4\pi \)) which for reasonably uniform spherical sampling can be approximated by \({\varDelta }_i\approx \frac{4 \pi }{N_o}\), and
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 \({\varDelta }_i\) more accurately. See for example [29, Eq. 83]).
Again, an exact reconstruction is possible iff \(0<\delta \le M_\psi ^d (\varvec{\omega })\le M<\infty \) and we have norm preservation when \(M_\psi ^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 \in {\mathbb {L}}_{2}^{\varrho }({\mathbb {R}}^3)\) to an orientation integrated score
This transformation has condition number \( {{\mathrm{cond}}}(A_\psi ^d) = \frac{\max \limits _{\varvec{\omega }\in B_{\varrho }}N_{\psi }^d(\varvec{\omega })}{\min \limits _{\varvec{\omega }\in B_{\varrho }}N_{\psi }^d(\varvec{\omega })} \).
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 \(\delta ,M>0,\ 0 < \varepsilon \ll 1\). Let \(\varrho \) be an a priori maximum frequency of our ball-limited image data. Then, a set of wavelets \(\{ \psi _{{\mathbf {n}}_i} \in {\mathbb {L}}_{2}^{\varrho }({\mathbb {R}}^{3}) \cap {\mathbb {L}}_{1}({\mathbb {R}}^{3}) \,|\, i=1,\dots ,N_o \}\), with a reasonable uniform spherical sampling (\({\varDelta }_i \approx \frac{4\pi }{N_o}\)), constructed as rotated versions of \(\psi \) is called a proper wavelet set if
where \({\mathbf {R}}_{{\mathbf {e}}_z, \alpha } \in SO(3)\) is a 3D rotation around axis \({\mathbf {e}}_z\) over angle \(\alpha \).
If, moreover, one has
then we speak of a proper wavelet with fast reconstruction property, cf. (43).
2.2.1 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
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 (\(N_o\)) of orientations.
-
2.
The transformation should be invertible, and reconstruction should be achieved by summation. Therefore, we aim for \(N_\psi ^d \approx 1\). Additionally, to guarantee all frequencies are transferred equally to the orientation score domain we aim for \(M_\psi ^d \approx 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 \(({\mathcal {F}}\psi ) (\varvec{\omega }) =g(\rho ) h (\vartheta ,\varphi )\), with
$$\begin{aligned} \begin{aligned} \varvec{\omega }&= (\omega _x,\omega _y,\omega _z) \\&= (\rho \sin \vartheta \cos \varphi ,\rho \sin \vartheta \sin \varphi ,\rho \cos \vartheta ). \end{aligned} \end{aligned}$$(49)Because by definition the wavelet \(\psi \) has rotational symmetry around the z-axis we have
. -
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.
.
Radial part g of \({\hat{\psi }}\), see Eq. (50) and radial parts \(g_0\) and \(g_1\) of \({\hat{\psi }}_0\) and \({\hat{\psi }}_1\) after splitting according to Sect. 2.1.1. The parameter \(\gamma \) controls the inflection point of the error function, here \(\gamma =0.8\). The steepness of the decay when approaching \(\rho _{\mathcal {N}}\) is controlled by the parameter \(\sigma _{erf}\) with default value \(\sigma _{erf} = \frac{1}{3} (\rho _{\mathcal {N}}-\varrho )\). At what frequency the splitting of \({\hat{\psi }}\) in \({\hat{\psi }}_0\) and \({\hat{\psi }}_1\) is done is controlled by parameter \(\sigma _\rho = \sqrt{2 s_\rho }\), see Eq. (18)
When directly setting orientation distribution A of Eq. (35) as angular part of the wavelet h we construct plate detectors. From left to right: Orientation distribution A, wavelet in the Fourier domain, the plate detector (real part) and the edge detector (imaginary part). Orange: Positive iso-contour. Blue: Negative iso-contour. Parameters used: \(s_o=\frac{1}{2}(0.25)^2,\sigma _{{{\mathrm{erf}}}}=3,\gamma =0.85\) and evaluated on a grid of \(51\times 51\times 51\) pixels (Color figure online)
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 (18) in Sect. 2.1.1 takes place. Following requirement 4 we only consider polar separable wavelets in the Fourier domain, so that 
. To satisfy requirement 2 we should choose radial function \(g (\rho ) = 1\) for \(\rho \in [0,\varrho ]\). In practice, this function should go to 0 when \(\rho \) tends to the Nyquist frequency \(\rho _{\mathcal {N}}\) to avoid long spatial oscillations. For the radial function \(g (\rho )\) we use,
with \({{\mathrm{erf}}}(z)= \frac{2}{\sqrt{\pi }} \int _{0}^{z} e^{-x^2} \;\mathrm{d}x\), 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 \ll \gamma < 1\). The steepness of the decay when approaching \(\rho _{\mathcal {N}}\) is controlled by the parameter \(\sigma _{erf}\) which we by default set to \(\sigma _{erf} = \frac{1}{3} (\rho _{\mathcal {N}}-\varrho )\). 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 \(\varrho \approx \rho _{\mathcal {N}}\)), and because radial function g causes \(M_\psi ^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_\psi ^d\) by \(\max (M_\psi ^d,\epsilon )\) in Eq. (5), with \(\epsilon \) small. Both make the reconstruction stable at the cost of not completely reconstructing the highest frequencies which causes a small amount of blurring.
We now need to find an appropriate angular part 
for the cake wavelets. First, we specify an orientation distribution \(A:S^2\rightarrow {\mathbb {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 \({\mathbf {n}}(\vartheta ,\varphi ) = (\sin \vartheta \cos \varphi ,\sin \vartheta \sin \varphi ,\cos \vartheta )\). 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_\psi ^d\) at the cost of less directionality.
First consider setting \(h=A\) so that \(\psi \) 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,Footnote 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) is the same for all values of \(\varphi \). As to the construction of \(h_ {Re}\), there is the general observation that 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
where integration is performed over \(S_p ({\mathbf {n}})\) denoting the great circle perpendicular to \({\mathbf {n}}\). This transformation preserves the symmetry of A, so we have 
. Thus, we finally set
For an overview of the transformations, see 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: \(s_o=\frac{1}{2}(0.25)^2,\gamma =0.85\) and evaluated on a grid of \(81\times 81\times 81\) pixels
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(\rho )A({\mathbf {n}}(\vartheta ,\varphi ))\) for the different filters, before applying the Funk transform to the orientation distribution A
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.
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 \(\{Y_l^m\}\) up to order L:
The spherical harmonics are given by
where \(P_l^m\) is the associated Legendre function, \(\epsilon = (-1)^m\) for \(m<0\) and \(\epsilon = 1\) for \(m>0\) and with integer \(l\ge 0\) and integer m satisfying \(-l \le m \le 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
3.2.1 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_l^m \rightarrow P_l(0)\,a_l^m\).
3.2.2 Anti-symmetrization
We have \(Y_l^m(\pi -\vartheta ,\varphi +\pi )=(-1)^l Y_l^m(\vartheta ,\varphi )\). We therefore anti-symmetrize the orientation distribution, see Eq. (53), via \(a_l^m\rightarrow \frac{(1-(-1)^l)}{2} a_l ^ m\).
3.2.3 Making Rotated Wavelets
To make the rotated versions \(\psi _{\mathbf {n}}\) of wavelet \(\psi \), we have to find \(h_{\mathbf {n}}\) in \({\hat{\psi }}_{\mathbf {n}}(\varvec{\omega }) =g (\rho )\,h_{\mathbf {n}}(\vartheta ,\varphi )\). 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_{m,m '} ^ l (\gamma ,\beta ,\alpha )\) (Wigner-D functions [41]):
Here \(\alpha , \beta \) and \(\gamma \) denote the Euler angles with counterclockwise rotations, where we rely on the convention \({\mathbf {R}}_{\gamma ,\beta ,\alpha } ={\mathbf {R}}_ {{\mathbf {e}}_z,\gamma } {\mathbf {R}}_ {{\mathbf {e}}_y,\beta } {\mathbf {R}}_ {{\mathbf {e}}_z,\alpha }\). This gives
where \(c_l^m\) are the coefficients of h given by
Because both anti-symmetrization and Funk transform preserve the rotational symmetry of A, we have \(h(\vartheta ,\varphi ) =\sum _{l = 0} ^L c_l^0 Y_l ^ 0 (\vartheta ,\varphi )\), and Eq. (62) reduces to
The filters from this section are summarized in the following result:
Result 1
Let \(A:S^2 \rightarrow {\mathbb {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 \({\hat{\psi }}\) in the Fourier domain via
with \(\varvec{\omega }=\rho \, {\mathbf {n}}_{\varvec{\omega }} = \rho \, {\mathbf {n}}(\vartheta ,\varphi )\). The real part of \(\psi \) is a tube detector given by
The imaginary part of \(\psi _{{\mathbf {n}}}\) is an edge detector given by
When expanding the angular part in spherical harmonics up to order L and choosing \(A = G_{s_o}^{S^2}\):
we have the following wavelet in the Fourier domain
and the coefficients of A and \({\hat{\psi }}\) relate via
We obtain rotated versions of our filter via
with \({\mathbf {n}}={\mathbf {n}}(\beta ,\gamma )\).
As we do not have analytical expressions for the spatial wavelets \(\psi _{\mathbf {n}}\), we sample the filter in the Fourier domain using Eq. (71) and apply a DFT afterward. The wavelet \(\psi \) 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_L^0/a_0^0<10^{-3}\) (e.g., \(L = 21\) for \(s_o = \frac{1}{2} (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_\psi ^d \approx 1\). We now focus on finding bounds for \(N_\psi ^d\) such that we can choose our parameters in a deliberate way.
Proposition 1
Let \(\{\psi _{{\mathbf {n}}(\beta _i,\gamma _i)} \,|\, i=1 \dots N_o\}\) be a set of wavelets constructed via the procedure in Result 1. Then we have bounds on \(N_{\psi }^d\) given by
with \({\mathbf {d}}_l=(d_{l}^{m})_{m=-l}^l\), \(d_{l}^{m}= \sum \limits _{i=1}^{N_o} c_{l}^{0} \cdot {\varDelta }_i \cdot D^{l}_{0,m}(0,\beta _i,\gamma _i)\) and here the norm is the \(\ell _2\)-norm on \({\mathbb {C}}^{2l+1}\).
Proof
First we expand function \(N_\psi ^d\) in spherical harmonics:
We have \(g (\rho ) = 1\) for \(\rho =\Vert \varvec{\omega }\Vert \le \rho _0\), but we still need to quantify the angular part. We define \({\mathbf {Y}}_{l}^N= (Y_{l}^{-l},Y_{l}^{-l+1},\dots ,Y_{l}^{l-1},Y_{l}^{l})\), so that
This varying component should remain small. We use the Cauchy-Schwarz inequality for each order l:
from which (73) follows. \(\square \)
See Fig. 9 for visual inspection of bounds of \(M_\psi ^d\) and \(N_\psi ^d\), and numerical results for the bounds of \(N_\psi ^d\).
Corollary 1
Given our analytical bounds (73) from Proposition 1 and \(N_o =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 \(\epsilon = 0.05\) when choosing parameter \(s_0 \gtrapprox 0.04\). In practice we have a proper wavelet set with fast reconstruction property already for smaller values of \(s_o\) (see Fig. 9).
Inspection of the stability of the transformation for different values of \(s_o\) given an orientation distribution \(A = G_{s_o}^{S^2}\) and for \(N_o=42\). Left: Spherical plot of A and the angular part of polar separable function \(N_\psi ^d\) and \(M_\psi ^d\). Orientation coverage is more uniform as the plots for\(N_\psi ^d\) and \(M_\psi ^d\) look more like a ball. Right: The upper and lower bounds of \(N_\psi ^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 (\(\approx 500\) orientations). Furthermore, we show \(1+\epsilon \) and \(1-\epsilon \) (orange dashed lines) for \(\epsilon =0.05\) (Color figure online)
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.
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=\Vert {\mathbf {x}}\Vert ^2-{\varDelta }\), 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 (\({\mathcal {F}}\circ H= H \circ {\mathcal {F}}\)) and eigenfunctions of H are also eigenfunctions of \({\mathcal {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_l^m\) the spherical harmonics, \((r,\theta ,\phi )\) and \((\rho ,\vartheta ,\varphi )\) spherical coordinates for \({\mathbf {x}}\) and \(\varvec{\omega }\), respectively, i.e.,
and \(g_n^l\) given by
where \(L_{n}^{(\nu )}(\rho )\) is the generalized Laguerre polynomial. We then choose the case with least radial oscillations \(\alpha _l^n=\alpha ^l \delta _n^0\). If we then choose
we have that \(M_\psi (\varvec{\omega })\) approximates 1 for all \(\varvec{\omega }\in {\mathbb {R}}^3\) in the Fourier domain as \(L \rightarrow \infty \) and we get the following waveletFootnote 4 [27]:
For this wavelet we have an analytical description in both spatial and Fourier domain.
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)
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 \({\mathbb {L}}_2({\mathbb {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 \in {\mathbb {L}}_{2}^{\varrho }({\mathbb {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].
4.2.1 The 3D Generalized Zernike Basis
The generalized Zernike functions are given by
with spherical coordinates
integer \(n,l\ge 0\) such that \(n=l+2p\), integer \(p\ge 0\) and \(m=-l,-l+1,\ldots ,l\) and \(\alpha >0\). The angular part is the spherical harmonic function and the radial part is given by
where \(P_{p}^{(\alpha ,l+\frac{1}{2})}\) denotes the Jacobi polynomial. The generalized Zernike functions are orthogonal on the unit ball
with \(\delta \) the Kronecker delta and with normalization factor
in which \((x)_\alpha =\frac{{\varGamma }(x+\alpha )}{{\varGamma }(x)}\) is the (generalized) Pochhammer symbol.
Fourier Transform
The inverse Fourier transform of the generalized Zernike function
is given by
with \({\mathbf {x}}= r \,{\mathbf {n}}(\theta ,\phi )\) and
Here \(J_a\) and \(j_a\) are the Bessel functions and spherical Bessel functions [55]. For integer \(\alpha \), 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., \(({\mathcal {F}}\psi )(\varvec{\omega })= F(\varvec{\omega }) = A (\vartheta ,\varphi ) B(\rho )\). When expanding such a function in the generalized Zernike basis,
we can split the coefficients in radial coefficients and angular coefficients
where
The coefficients \(c_{n,l}^{m,\alpha }\) 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,\dots \)
For each l, the radial functions \(R_n^{l,\alpha }\) with n varying are a basis for functions defined on the interval [0, 1]. For separable functions, we expand the same radial function \(B(\rho )\) for each l, and it can be shown that there is a recursion formula for the radial coefficients [43].
4.2.2 Wavelets
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({\mathbf {n}}(\vartheta ,\varphi )) = G_{s_o}^{S^2}({\mathbf {n}}(\vartheta ,\varphi ))\) 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 \(\psi _0\) and \(\psi _1\). In this section we design a radial function for \(\psi _1\) which is relevant for further processing. Furthermore, we already have an analytical description for \(\phi _0\), which we set to \(\phi _0=G_{s_\rho }\) (see Eq. (37)).
The radial function of \(\psi _0\) should therefore approximate \(B(\rho )=1-G_{s_{\rho }}(\rho )\) on the interval \([0,\varrho ]\) and should smoothly go to zero when approaching the edges of the interval. For the moment, we set \(\varrho =1\) and we include the scaling later. To start, we define the function
see Fig. 11a for the case \(\alpha =6,\beta =2\). For this function we have the following coefficients
To obtain a flatter function we multiply the function \(B_{\alpha ,\beta }\) with a second-order Taylor expansion of the reciprocal function \(\rho \mapsto (B_{\alpha ,\beta }(\rho ))^{-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 \(\beta \), so we can find the coefficients \(b_{n=l+2p}^{l,\alpha }\) for the flattened function as well. For the specific case \(\beta =2\) we get the following flattened function
with \(B_{\text {max}} = B_{\alpha ,2}(\rho _{\text {max}})\). For this flattened function the coefficients are given by
with \(c_i\) the coefficients of \(\rho ^0\), \(\rho ^2\) and \(\rho ^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:
Result 2
(Analytic 3D-wavelets in Zernike basis) Let \(\alpha >0\) and let \(A:S^2 \rightarrow {\mathbb {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 \({\hat{\psi }}\) in the Fourier domain via Eq. (65). The real part of \(\psi \) is a tube detector and the imaginary part of \(\psi \) is an edge detector, see Fig. 7. We choose radial function \(g(\rho )=B^{\text {flat}}_{6,2}(\frac{\rho }{\rho _{\mathcal {N}}})\) in Eq. (18) for \(\psi _1\) and angular function \(A({\mathbf {n}}(\vartheta ,\varphi )) = G_{s_o}^{S^2}({\mathbf {n}}(\vartheta ,\varphi ))\) and expand in the generalized Zernike basis:
The coefficients \(c_{n,l}^0\) follow by expanding A in spherical harmonics and \(B^{\text {flat}}_{6,2}\) in the radial Zernike polynomials, recall Eqs. (95) and (96). This yields \(a_l^0\) (Eq. (58)) and \(b_{n,l}^{\alpha }\) (Eq. (101)) and coefficients \(c_{n,l}^0\):
The spatial wavelet is given by
with \(S_{n,l}^{\alpha }\) given by Eq. (89) and \(Y_{l}^m\) the spherical harmonics of Eq. (57). Then we obtain rotated filters via
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).
Left: Function \(B_{6,2}(\rho )=(1-\rho ^2)^{6 } \rho ^{2}\). Right: Flattened function which are obtained from \(B_{\alpha ,\beta }(\rho )\) by multiplying with the second-order Taylor approximation of its reciprocal around the function maximum: \(B_{6,2}^{\text {flat}}(\rho )=\frac{1}{B_\text {max}}(1+\frac{7^3}{12}(\rho ^2-\frac{1}{7})^2)B_{6,2}(\rho )\)
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.
Comparison of the filters obtained by sampling in the Fourier domain and performing an inverse DFT (Result 1 in Sect. 3) and the filters expressed in the generalized Zernike basis (Result 2 in Sect. 4). Left: Iso-contour plot of the filter aligned with the x-axis showing one positive iso-contour (orange) and one negative iso-contour (blue). Middle: Cross section of the filter for \(z=0\). Right: The low-pass filter. Top Filters according to Result 1 with parameters \(s_{\rho }=\frac{1}{2} (1.9)^2\) and \(\gamma =0.85\). Bottom: The filters according to Result 2 with \(\alpha =3\) and \(\beta =2\). Both have \(s_{o}=\frac{1}{2} (0.4)^2\) and are evaluated on a grid of \(31 \times 31 \times 31\) voxels (Color figure online)
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. We see no qualitative differences between the two reconstructions.
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
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.
Overview of the datasets used in our experiments
6.1 Diffusion Via Invertible Orientation Scores
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 coherence-enhancing diffusion via 2D orientation scores have been developed for 2D data [37, 64]. Here, we propose an extension 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 \({\mathbb {R}}^ 3\times 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 \({\mathbf {R}}_{\mathbf {n}}\) any rotation for which \({\mathbf {R}}_{\mathbf {n}}\cdot {\mathbf {e}}_z = {\mathbf {n}}\). The operators \({\varPhi } \) which we consider are scale spaces on SE(3) (diffusions), and are given by \({\varPhi }={\varPhi }_t\) with
Here \({\tilde{W}}\) 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 condition \({\tilde{W}}(g, 0)\) based on exponential curve fit (see [32]), and with \({\mathcal {A}}_i|_{g}=(L_g)_* {\mathcal {A}}_i|_e\) the left-invariant vector fields on SE(3). The diffusion is better understood in locally adaptive frame \(\{{\mathcal {B}}_i\}_{i=1}^6\). Here \({\mathcal {B}}_3\) follows from an exponential curve fits and points along our structure. \({\mathcal {B}}_1\) and \({\mathcal {B}}_2\) span the plane spatially perpendicular to our structure, and \({\mathcal {B}}_4, {\mathcal {B}}_5\) and \({\mathcal {B}}_6\) correspond to angular diffusion (embedding in SE(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}\).
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
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
with \(c_1\) a constant automatically set to the \(50\%\) quantile of \(s({\tilde{U}})\) and \(s({\tilde{U}})(g)\) a measure for orientation confidence given by minus the Laplacian in the space orthogonal to the structure orientation \({\mathcal {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 \(|{\mathcal {B}}_3|_{g}({\tilde{U}})|\).
We then obtain Euclidean invariant image processing via
which includes inherent projection \({\mathbb {P}}_\psi \) of orientation scores, even if \({\varPhi }={\varPhi }_t\) maps outside of the space of orientation scores. We write \({\mathcal {W}}_\psi ^ {-1, \text {ext}}\) because we extend the inverse to \({\mathbb {L}}_2({\mathbb {R}}^3 \times S^2)\).
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:
Given such data, the noise is quantified via the contrast-to-noise ratio (CNR):
where \(\sigma (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 \({\varOmega }_S\) and background regions given by \({\varOmega }_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 \({\varOmega }_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.
6.1.4 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=\frac{1}{2}(0.45)^2\) and \(s_\rho = \frac{1}{2}(16)^2\). For CED we used the following settings: \(\alpha =0.2\) and c the \(50\%\)-quantile of \(\kappa \) (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 \(\alpha \) 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).
Results of coherence-enhancing diffusion via orientation scores on dataset B.1, see Eq. (110). Top: CEDOS for different amounts of diffusion time. Middle: Result for CED. Bottom: Result for isotropic Gaussian regularization. For all datasets, we show one iso-contour at \(\mu _\text {BG}+ 0.7 (\mu _\text {FG}-\mu _\text {BG})\), where \(\mu _\text {BG}\) and \(\mu _\text {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
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 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.
Measurement of vessel radius in vessel cross sections after different amounts of diffusion in dataset B.1. Top left: 3D Visualization of the data with the selected slices. Top middle: Radii measurements for increasing diffusion time for CEDOS. Top right: Radii measurements for increasing diffusion time for Gaussian regularization. The detected 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)
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.
Schematic visualization of the edges used in the tubularity measure \(V({\mathbf {x}}, {\mathbf {n}}, r)\). Left: A 3D iso-contour visualization of a vessel with orientation \({\mathbf {n}}\). The coordinates \((\theta ,r)\) are polar coordinates for the plane perpendicular to orientation \({\mathbf {n}}\) spanned by \({\mathbf {e}}_1\) and \({\mathbf {e}}_2\). Right: In this plane opposite edges are multiplied in \(E_{\text {prod}}({\mathbf {x}}, {\mathbf {n}}, r, \theta )\). The edge is expected to have outward orientation \({\mathbf {n}}^\perp (\theta )\)
Tubularity on artificial datasets and comparison of measured radius against ground truth radius for datasets C.1, C.2 and C.3 (top to bottom). Left: The data. Middle left: The centerline with ground truth radius in color. Middle right: The tubularity confidence \(s^t({\mathbf {x}})\) (max of tubularity over radius and orientation) with radius of max response \(r^*({\mathbf {x}})\) in color. Right: Measured radius \(r^*({\mathbf {x}})\) against ground truth radius \(r_\text {GT}({\mathbf {x}})\) on the ground truth centerline. The opacity of the plotted points is linearly scaled with the tubularity confidence
6.2.2 Tubularity Via Orientation Scores
For the tubularity measure we detect edges in the plane perpendicular to orientation \({\mathbf {n}}\). Within this plane, the product of two opposite edges at radius r and in-plane orientation \(\theta \) at position \({\mathbf {x}}\in {\mathbb {R}}^3\) is given by
where \({{\mathrm{Im^+}}}(z) = \max \{0,{{\mathrm{Im}}}(z)\}\) and \({\mathbf {n}}^\bot (\theta ) = \cos \theta \,{\mathbf {e}}_1 + \sin \theta \,{\mathbf {e}}_2\), with \(\{{\mathbf {e}}_1,{\mathbf {e}}_2\}\) an orthogonal basis for the orthogonal complement of \(\langle {\mathbf {n}}\rangle ={\text {span}}\{{\mathbf {n}}\}\). The product of the two \({{\mathrm{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 \(\theta \). This is done as follows
For the angular regularization kernel we use the simple \(2\pi \) periodic 1D-diffusion kernel
\(K^\text {or}_{{\sigma _o}^V}(\theta ) = \sum \limits _{k \in {\mathbb {Z}}} G_{S^1}^{\sigma _{o}^{V}}(\theta + 2k \pi )= \sum \limits _{k \in {\mathbb {Z}}} \frac{1}{\sqrt{2\pi } \sigma _o^v} e^{- \frac{|\theta +2 k \pi |^2}{2 |\sigma _o^V|^2}}\) over \(\theta \), instead of the true diffusion kernel (72), where for practical values we use \(\sigma _o^V=\pi /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({\mathbf {x}})\) is the tubularity confidence which is a measure for how certain we are at least one tubular structure is present at position \({\mathbf {x}}\). The features \({\mathbf {n}}^*({\mathbf {x}})\) and \(r^*({\mathbf {x}})\) are the orientation and radius of optimal tubularity response at position \({\mathbf {x}}\).
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: \(\sigma _{o}^{V} = \pi /8\) and we discretized the \(\theta \)-integral in Eq. (121) using 8 orientations.
Measured radius \(r^*({\mathbf {x}})\) against ground truth radius on the ground truth centerline for all 18 datasets. The opacity of the plotted points is linearly scaled with the tubularity confidence
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({\mathbf {x}},\cup B_{{\mathbf {c}}_i,r^*({\mathbf {c}}_i)})=\min _{{\mathbf {c}}_i} \left\{ \Vert {\mathbf {x}}-{\mathbf {c}}_i \Vert -r^*({\mathbf {c}}_i) \right\} \), where \({\mathbf {c}}_i\) are the positions given by the 1% quantile of highest responses
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 Gaussians, and we define the ground truth radius \(r_{GT}({\mathbf {x}})=\sigma ({\mathbf {x}})\) at the standard deviation \(\sigma ({\mathbf {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: \(\sigma _{o}^{V} = \pi /8\) and we discretized the \(\theta \)-integral in Eq. (121) using 12 orientations. For the orientation score transformation we used \(s_o=\frac{1}{2}(0.4)^2, s_\rho = \frac{1}{2}(5)^2\) and evaluated on a grid of \(21 \times 21 \times 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 \({\mathbf {n}}^*({\mathbf {x}})\) and a simple segmentation given by the 0-iso-contour of the distance map \(d({\mathbf {x}},\cup B_{{\mathbf {c}}_i,r^*({\mathbf {c}}_i)})=\min _{{\mathbf {c}}_i} \left\{ \Vert {\mathbf {x}}-{\mathbf {c}}_i \Vert -r^*({\mathbf {c}}_i) \right\} \), where \({\mathbf {c}}_i\) 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.
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 diffusion-weighted MRI and to compare with existing approaches [58, 68].
Notes
In practice we choose the default values \(\delta =\frac{1}{8}\) and \(M=1.1\) and \(\varepsilon =0.01\) and note it is actually the ratio \(\frac{M}{\delta }\) that determines the condition number. It is just a convenient choice to set the upper bound close to 1.
This can be observed in Fig. 6, where we recognize a planar derivative filter in the direction of \({\mathbf {n}}\) because it is anti-symmetric and has a mainly positive and mainly negative side.
The series in (81) converges point-wise but not in \({\mathbb {L}}_2\)-sense. So for taking the limit \(L \rightarrow \infty \) one must rely on a distributional wavelet transform, since \(\psi _H^{L \rightarrow \infty } \notin {\mathbb {L}}_2({\mathbb {R}}^3)\) but \(\psi _H^{L \rightarrow \infty } \in H_{-4}({\mathbb {R}}^3)\), 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.
References
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)
Ali, S.T.: A general theorem on square-integrability: vector coherent states. J. Math. Phys. 39(8), 3954 (1998)
Ali, S.T., Antoine, J.P., Gazeau, J.P.: Coherent States, Wavelets, and Their Generalizations. Graduate Texts in Contemporary Physics. Springer, New York (2000)
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)
Barbieri, D., Citti, G.: Reproducing kernel Hilbert spaces of CR functions for the Euclidean motion group. Anal. Appl. 13(03), 331–346 (2014)
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)
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)
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)
Bekkers, E.J.: Retinal Image Analysis Using Sub-Riemannian Geometry in SE(2). Ph.D. Thesis, University of Technology, Eindhoven (2017)
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)
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)
Born, M., Wolf, E.: Principles of Optics, 7th edn. Cambridge University Press, Cambridge (1999)
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)
Bouix, S., Siddiqi, K., Tannenbaum, A.: Flux driven automatic centerline extraction. Med. Image Anal. 9(3), 209–221 (2005)
Brunner, T.A.: Impact of lens aberrations on optical lithography. IBM J. Res. Dev. 41(1.2), 57–67 (1997)
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)
Burgeth, B., Pizarro, L., Didas, S., Weickert, J.: 3D-coherence-enhancing 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)
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)
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)
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)
Citti, G., Franceschiello, B., Sanguinetti, G., Sarti, A.: Sub-Riemannian mean curvature flow for image processing. SIAM J. Imaging Sci. 9(1), 212–237 (2016)
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)
Descoteaux, M., Angelino, E., Fitzgibbons, S., Deriche, R.: Regularized, fast, and robust analytical Q-ball imaging. Magn. Reson. Med. 58(3), 497–510 (2007)
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)
Duits, R.: Perceptual Organization in Image Analysis. Ph.D. Thesis, Technische Universiteit Eindhoven (2005)
Duits, R., Burgeth, B.: Scale spaces on lie groups. In: Scale Space and Variational Methods in Computer Vision, pp. 300–312. Springer, Berlin (2007)
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)
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)
Duits, R., Franken, E.: Left-invariant diffusions on the space of positions and orientations and their application to crossing-preserving smoothing of HARDI images. Int. J. Comput. Vis. 92(3), 231–264 (2010)
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)
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)
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)
Dunkl, C.F.: Orthogonal Polynomials of Several Variables. Cambridge University Press, Cambridge (2014)
Faraut, J., Harzallah, K.: Deux Cours d’Analyse Harmonique. Birkhauser, Basel (1987)
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)
Franken, E.: Enhancement of Crossing Elongated Structures in Images. Ph.D. Thesis, Technische Universiteit Eindhoven (2008)
Franken, E., Duits, R.: Crossing-preserving coherence-enhancing diffusion on invertible orientation scores. Int. J. Comput. Vis. 85(3), 253–278 (2009)
Freeman, W.T., Adelson, E.H.: The design and use of steerable filters. IEEE Trans. Pattern Anal. Mach. Intell. 13(9), 891–906 (1991)
Führ, H.: Abstract Harmonic Analysis of Continuous Wavelet Transforms. No. 1863 in Lecture Notes in Mathematics. Springer, Berlin (2005)
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)
Griffiths, D.J.: Introduction to Quantum Mechanics, 2nd edn. Cambridge University Press, Cambridge (2016)
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)
Janssen, A.J.: Generalized 3d Zernike Functions for Analytic Construction of Band-Limited Line-Detecting Wavelets. ArXiv preprint arXiv:1510.04837 (2015)
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)
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)
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)
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)
Lee, T.S.: Image representation using 2D Gabor wavelets. IEEE Trans. Pattern Anal. Mach. Intell. 18(10), 959–971 (1996)
Louis, A.: Orthogonal function series expansions and the null space of the Radon transform. SIAM J. Math. Anal. 15(3), 621–633 (1984)
Louis, A., Maass, D., Rieder, A.: Wavelets: Theory and Application. Wiley, New York (1997)
Mallat, S.: A Wavelet Tour of Signal Processing. Academic Press, London (1999)
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)
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)
Muhlich, M., Aach, T.: Analysis of multiple orientations. IEEE Trans. Image Process. 18(7), 1424–1437 (2009)
Olver, F.W.J.: NIST Handbook of Mathematical Functions. Cambridge University Press, Cambridge (2010)
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)
Petitot, J.: The neurogeometry of pinwheels as a sub-Riemannian contact structure. J. Physiol. 97(2–3), 265–309 (2003)
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)
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)
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)
Reisert, M.: Group Integration Techniques in Pattern Analysis—A Kernel View. Ph.D. Thesis, Albert-Ludwigs- University (2008)
Ruijters, D., Vilanova, A.: Optimizing GPU volume rendering. J. WSCG 14(1–3), 9–16 (2006)
Scharr, H.: Diffusion-like reconstruction schemes from linear data models. In: Pattern Recognition, pp. 51–60. Springer, Berlin (2006)
Sharma, U., Duits, R.: Left-invariant evolutions of wavelet transforms on the similitude group. Appl. Comput. Harmon. Anal. 39(1), 110–137 (2015)
Sifre, L., Mallat, P.S.: Rigid-Motion Scattering For Image Classification. Ph.D. Thesis, Ecole Polytechnique, CMAP, Paris (2014)
Steidl, G., Teuber, T.: Anisotropic smoothing using double orientations. In: Scale Space and Variational Methods in Computer Vision, pp. 477–489 (2009)
van Almsick, M.: Context Models of Lines and Contours. Ph.D. Thesis, Technische Universiteit Eindhoven (2007)
Vogt, T., Lellmann, J.: Measure-Valued Variational Models with Applications to Diffusion-Weighted Imaging. ArXiv preprint arXiv:1710.00798 (2017)
Weickert, J.: Coherence-enhancing diffusion filtering. Int. J. Comput. Vis. 31(2–3), 111–127 (1999)
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)
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)
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).
Author information
Authors and Affiliations
Corresponding author
Appendices
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 \(\{g_{n}^{l}(r)\}_{n=0}^{\infty }\) for a radial orthonormal basis for \({\mathbb {L}}_{2}({\mathbb {R}}^+, r^2 \mathrm{d}r)\) associated to each \(2l+1\)-dimensional SO(3)-irreducible subspace indexed by l. Then by the celebrated Bochner–Hecke theorem [34], this induces a corresponding orthonormal radial basis \(\{{\tilde{g}}_n^l(\rho )\}_{n=0}^{\infty }\) in the Fourier domain which can be obtained by a Hankel-type of transform [25, Ch. 7]. We expand our wavelets in spherical harmonics \(Y_{l}^{m}\) and ball coordinates (cf. (78)) accordingly:
Remark 4
In Sect. 4, we consider the modified Zernike basis in which case \(g^{l}_n\) and \({\tilde{g}}^{l}_n\) are given by, respectively, (85) and (82), whereas for the harmonic oscillator basis one has \(g^{l}_{n}= i^l (-1)^{n+l} {\tilde{g}}^{l}_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}(\gamma ,\beta ,0) \in {\mathbb {R}}\) (cf. Sect. 3.2.3) and one obtains the following steerable implementations of orientation scores:
where \({\mathbf {n}}=(\cos \gamma \sin \beta , \sin \gamma \sin \beta , \cos \beta )^T\), \(\star \) denotes correlation, the overline denotes complex conjugation and with function product \(({\tilde{g}}_{n}^l \otimes Y_{l}^{m})({\mathbf {x}}) = {\tilde{g}}_{n}^l(r)\; Y_{l}^{m}(\theta ,\phi )\).
B Table of Notations
Symbol | Explanation | Reference |
|---|---|---|
B.1 Spaces and Input Data | ||
\({\mathbb {R}}^{3} \times S^{2}\) | Space of positions and orientations | Page 1 |
\({\mathbb {L}}_2^{\varrho } ({\mathbb {R}}^3)=\{f \in {\mathbb {L}}_2 ({\mathbb {R}}^3)| \text {supp}({\mathcal {F}}f) \subset B_{\varrho }\}\), with \(\varrho >0\) | Space of frequency ball-limited images | (2) |
\(B_{\rho }=\{{\mathbf {x}}\in {\mathbb {R}}^3 \big | \Vert {\mathbf {x}}\Vert < \rho \}\), with \(\rho >0\) | The ball around the origin with radius \(\rho \) | (2) |
B.2 Orientation Score Transformation | ||
\({\mathcal {W}}_\psi \) | Orientation score transformation | (3) |
\({\mathcal {W}}_\psi ^{-1}\) | Data reconstruction via the inverse orientation score transformation | (5) |
\(\psi \) | Wavelet used for the orientation score transformation | (3) |
\({\hat{\psi }}\) | Fourier transform of the wavelet used for the orientation score transformation | (6) |
\(M_{\psi }\) | Factor used to quantify the stability of the transformation | (6) |
\(N_{\psi }\) | Factor used to quantify the stability of the transformation when using the simplified reconstruction by integration | (12) |
\({\mathcal {W}}_\psi ^d, ({\mathcal {W}}_\psi ^d)^{-1}, M_{\psi }^d, N_{\psi }^d\) | Discretized versions of \({\mathcal {W}}_\psi , ({\mathcal {W}}_\psi )^{-1}, M_{\psi }, N_{{}\psi }\) | |
\(\mathrm {d}\sigma \) and \({\varDelta }_i\) | Spherical area measure and discretized spherical area measure | |
B.3 Coordinates | ||
\({\mathbf {x}}=(x,y,z)\) | Cartesian coordinates real space | - |
\(\varvec{\omega }=(\omega _x,\omega _y,\omega _z)\) | Cartesian coordinates Fourier domain | - |
\((r,\theta ,\phi ), \quad {\mathbf {x}}= (r \sin \theta \cos \phi ,r \sin \theta \sin \phi ,r \cos \theta )\) | Spherical coordinates real space | (49) |
\((\rho ,\vartheta ,\varphi ),\quad \varvec{\omega }= (\rho \sin \vartheta \cos \varphi ,\rho \sin \vartheta \sin \varphi ,\rho \cos \vartheta )\) | Spherical coordinates Fourier domain | (49) |
B.4 Wavelets | ||
g | Radial function of the cake filters | (50) |
A | Orientation distribution used in wavelet construction | (52) |
\(Y_l^m\) | Spherical harmonics | (57) |
\(Z_{n,l}^{m,\alpha }\) | 3D generalized Zernike functions | (82) |
\(R_n^{l,\alpha }\) | Radial part of the 3D generalized Zernike functions | (84) |
\(S_{n,l}^{\alpha }\) | Radial part of the inverse Fourier transform 3D generalized Zernike functions | (89) |
B.5 Applications | ||
\(V({\mathbf {x}}, {\mathbf {n}}, r)\) | Tubularity measure | (121) |
\(s^t({\mathbf {x}})\) | Tubularity confidence | (122) |
\({\mathbf {n}}^*({\mathbf {x}})\) | Orientation of maximum tubularity response | (122) |
\(r^*({\mathbf {x}})\) | Radius of maximum tubularity response | (122) |
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.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.
About this article
Cite this article
Janssen, M.H.J., Janssen, A.J.E.M., Bekkers, E.J. et al. Design and Processing of Invertible Orientation Scores of 3D Images. J Math Imaging Vis 60, 1427–1458 (2018). https://doi.org/10.1007/s10851-018-0806-0
Received:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s10851-018-0806-0