Design and Processing of Invertible Orientation Scores of 3D Images for Enhancement of Complex Vasculature
Abstract
The enhancement and detection of elongated structures in noisy image data is relevant for many biomedical imaging applications. To handle complex crossing structures in 2D images, 2D orientation scores were introduced, which already showed their use in a variety of applications. Here we extend this work to 3D orientation scores . 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 cakewavelets, which are complex wavelets that can simultaneously detect oriented structures and oriented edges. Here we introduce two types of cakewavelets, the first uses a discrete Fourier transform, the second is designed in the 3D generalized Zernike basis, allowing us to calculate analytical expressions for the spatial filters. Finally, we show two applications of the orientation score transformation. In the first application we propose an extension of crossingpreserving coherence enhancing diffusion via our invertible orientation scores of 3D images which we apply 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.
Keywords:
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 is 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 multiorientation representations of image data. Such multiorientation representations can be made using various techniques, such as invertible orientation scores (which is obtained via a coherent state transform) Ali2000 (); DuitsQAM2010Part1 (); BekkersJMathImagingVis2014 (); ThesisFranken (); IshamJMathPhys1991 (); BarbieriAnalAppl2014 (), continuous wavelet transforms DuitsIJCV2007 (); DuitsQAM2010Part1 (); SharmaACHA2015 (); BekkersJMathImagingVis2014 (), orientation lifts ZweckJMathImagingVis2004 (); BoscainSIAMJImagingSci2014 (), or orientation channel representations Felsberg2012 (). 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 with a set of oriented filters to construct a 3D orientation score . 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 multiorientation 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. Lee1996 (); Antoine1999 (); Mallat1999 (); Louis1997 () 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 AliJMathPhys1998 (). This makes it harder to design appropriate wavelets, but has the computational advantage of only needing one allscale transformation.
The 2D orientation scores have already showed their use in a variety of applications. In FrankenIJCV2009 (); SharmaACHA2015 () the orientation scores were used to perform crossingpreserving 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 BekkersJMathImagingVis2014 (); BekkersSIAMJImagingSci2015 (); Chen2014 (), in vessel segmentation Zhang2016 () and biomarker analysis Bekkers2015 (); Rajan2016 (), 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 greyscale 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 Lee1996 ()). A wavelet that allows an invertible transformation was proposed by Kalitzin KalitzinIJCV1999 (). A generalization of these wavelets was found by Duits ThesisDuits () 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 ThesisDuits () a class of 2D cakewavelets were introduced, that have a cakepiece shaped form in the Fourier domain. The cakewavelets simultaneously detect oriented structures and oriented edges by constructing a complex orientation score . Because the family of rotated cakewavelets cover the full Fourier spectrum, invertibility is guaranteed.
In this article we propose a 3D version of the cakewavelets. A preliminary attempt to generalize these filters was done in JanssenSSVM2015 (), where the platedetectors in ThesisDuits () were extended to complexvalued cakewavelets with a linedetector 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 a steerable implementation. 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 crossingpreserving 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 FrankenIJCV2009 () and related works Scharr2006 (); steidl2009anisotropic (); muhlich2009analysis () that rely on other more specific orientation decompositions.
For the sake of general readability we avoid Lie group theoretical notations, until Section 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 Ali2000 (); FuhrBook2005 (); DuitsSSVM2007 (); ThesisDuits (); BatardJMathImagingVis2014 () that has many applications in image processing. Besides the special cases of the Heisenberg group BarbieriJMIV2014 (); PetitotJPhysiolParis2003 (); DuitsApplComputHarmonAnal2013 (), the 2D Euclidean motion group CittiJMathImagingVis2006 (); BarbieriAnalAppl2014 (); Mashtakov2013 (); Prandi2015 (); CittiSIAMJImagingSci2016 (); DuitsQAM2010Part1 (); ThesisBekkers (), the similitude group Antoine1999 (); SharmaACHA2015 (); Pechaud2009 () and the 3D rotation group MashtakovJMathImagingVis2017 (), we now consider invertible orientation scores on the 3D Euclidean motion group (in which the coupled space of positions and orientations is embedded), which relate to coherent states from physics IshamJMathPhys1991 () for , with a specific (semi)analytic design for our image processing purposes.
1.1 Contributions of the Article
The main contributions per section of the article are:

In Section 5 we compare the two types of filters and show that the DFT filters approximate their analytical counterparts well.

In Section 6 we show two applications of the orientation score transformation:

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 imagedata WeickertIJCV1999 (); BurgethBook2009 (); Burgeth2012 (), 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).

We develop a new tubularity measure in 3D orientation score data. This extends previous work on tubularity measures using 2D orientation scores Chen2014 ()(ThesisBekkers, , 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
First, we discuss the theory of invertible orientation score transforms in Section 2. Then we construct 3D cakewavelets and give a new efficient implementation using spherical harmonics in Section 3, followed by their analytical counterpart expressed in the generalized Zernike basis in Section 4. Then we compare the two types of filters and validate the invertibility of the orientation score transformation in Section 5. Finally, we address two application areas for 3D orientation scores in Section 6 and show results and practical benefits for both of them. In the first application (Subsection 6.1), we present a natural extension of the crossing preserving coherence enhancing diffusion on invertible orientation scores (CEDOS) FrankenIJCV2009 () to the 3D setting. In the second application (Subsection 6.2) we use the orientation score to define a tubularity measure, and show experiments applying the tubularity measure to both synthetic data and brainCT data. Both application sections start with a treatment of related methods.
2 Invertible Transformation
2.1 Continuous Orientation Score Transform
Throughout this article we use the following definition of the Fourier transform on :
(1) 
An invertible orientation score is constructed from a given balllimited 3D dataset
(2) 
with ball of radius , by correlation with an anisotropic kernel
(3) 
Here is a wavelet aligned with and rotationally symmetric around the axis, and the rotated wavelet aligned with given by
(4) 
Here is any 3D rotation which rotates the axis onto 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 for this transformation is given by
(5) 
with . The function is given by
(6) 
and vanishes at , where the circumflex again denotes Fourier transformation. Due to our restriction to balllimited data (2) this does not cause problems in reconstruction (5). The function quantifies the stability of the inverse transformation ThesisDuits (), since specifies how well frequency component is preserved by the cascade of construction and reconstruction when would not be included in Eq. (5). An exact reconstruction is possible as long as
(7) 
In practice it is best to aim for in view of the condition number of transformation given by:
(8) 
where in the codomain spatial frequencies are again limited to the ball:
(9) 
Also, in the case we have for we have norm preservation
(10) 
and reconstruction Eq. (5) simplifies to
(11) 
We can further simplify the reconstruction for wavelets for which the following additional property holds:
(12) 
In that case the reconstruction formula is approximately an integration over orientations only:
(13) 
For the reconstruction by integration over angles only we can analyze the stability via the condition number of the mapping that maps an image to an orientation integrated score
(14) 
Its condition number is given by .
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 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 normpreservation by controlling leads to restrictions on our wavelets which are captured in the following definition:
Definition 1 (Proper Wavelet)
Let us set a priori bounds^{1}^{1}1In practice we choose the default values and and and note it is actually the ratio that determines the condition number. It is just a convenient choice to set the upper bound close to 1. . Furthermore, let be an a priori maximum frequency of our balllimited image data. Then, a wavelet is called a proper wavelet if
(15)  
(16) 
where is the 3D rotation around axis over angle .
If moreover, one has
(17) 
then we speak of a proper wavelet with fast reconstruction property, cf. (13).
Remark 1
The 1st condition (symmetry around the 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 balllimited functions is reasonable. The constraint is not a necessary one when one relies on distributional transforms (BekkersJMathImagingVis2014, , 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 this data for reconstruction. Therefore we perform an additional splitting of our wavelets into two parts
(18) 
with Gaussian window in the Fourier domain given by
(19) 
After splitting, contains the average and low frequency components and 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. Mallat1999 (); Sifre2014 (). In this case we construct two scores. One for the highfrequency components
(20) 
and one for the lowfrequency components
(21) 
Here we again have , as in Eq. (4). The vector transformation is then defined as
(22) 
For this transformation we have the exact reconstruction formula
(23) 
with
(24) 
Again, 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 such that Eq. (7) holds, and . Then the condition number of is given by
(25) 
The condition number of obtained from by performing an additional splitting in low and high frequency components is given by
(26) 
thereby guaranteeing that stability is maintained after performing the splitting.
Proof
First, we find the condition number of :
(27) 
For the first factor in Eq. (27) we find
(28) 
Similarly, we get for the second factor in Eq. (27). Then we obtain
(29) 
Similarly the condition number of is given by
(30) 
Next we express in of the original wavelet as
(31) 
So it remains to quantify . For a wavelet splitting according to (18) we have
(32) 
Hence
(33) 
And since for we have for satisfying (7) the following bounds on :
(34) 
thereby guaranteeing stability after the splitting (18).
For this vector transformation we can also use the approximate reconstruction by integration (for ) over orientations. Thus we have
(35)  
As said we are only interested in processing of and not in processing of , and so we directly calculate via
(36) 
For a design with for all , we have and so
(37) 
Then Eq. (35) becomes
(38) 
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 an electrostatic repulsion model CaruyerMagnResonMed2013 ().
Assume we have a number of orientations , and define the discrete invertible orientation score by
(39) 
The exact reconstruction formula is in the discrete setting given by
(40) 
with the discrete spherical area measure () which for reasonably uniform spherical sampling can be approximated by (otherwise one could use (DuitsIntJComputVis2010, , Eq. (83))), and
(41) 
Again, an exact reconstruction is possible iff and we have norm preservation when =1. Again, for the wavelets for which
(42) 
the image reconstruction can be simplified by a summation over orientations:
(43) 
For this reconstruction by summation we can analyze the stability via the condition number of the mapping that maps an image to an orientation integrated score
(44) 
This transformation has condition number .
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 . Let be an a priori maximum frequency of our balllimited image data. Then, a set of wavelets , with a reasonable uniform spherical sampling (), constructed as rotated versions of is called a proper wavelet set if
(45)  
(46) 
where is a 3D rotation around axis over angle .
If moreover, one has
(47) 
then we speak of a proper wavelet with fast reconstruction property, cf. (43).
2.2.1 Low Frequency Components
2.3 Steerable Orientation Score Transform
Throughout this article we shall rely on spherical harmonic decomposition of the angular part of proper wavelets in spatial and Fourier domain. This has the benefit that one can obtain steerable FreemanPAMI1991 (); ThesisFranken (); ThesisReisert () 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 cakewavelets, see BekkersJMathImagingVis2014 (); FrankenIJCV2009 (); DuitsPRIA2007 (), was used for the 2D orientation score transformation. We now present 3D versions of these cakewavelets. Thanks to the splitting in Section 2.1.1 we no longer need the extra spatial window used there. Our 3D transformation using the 3D cakewavelets should fulfill a set of requirements, compare FrankenIJCV2009 ():

The orientation score should be constructed for a finite number () of orientations.

The transformation should be invertible and reconstruction should be achieved by summation. Therefore we aim for . Additionally, to guarantee all frequencies are transferred equally to the orientation score domain we aim for . The set should be a proper wavelet set with fast reconstruction property (Def. 2)

The kernel should be strongly directional.

The kernel should be separable in spherical coordinates in the Fourier domain, i.e., , with
(49) Because by definition the wavelet has rotational symmetry around the axis we have .

The kernel should be localized in the spatial domain, since we want to pick up local oriented structures.

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 orientation score. 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 cakewavelets before splitting in low and high frequencies according to (18) in Section 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 for . In practice, this function should go to 0 when tends to the Nyquist frequency to avoid long spatial oscillations. For the radial function we use,
(50) 
with , 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 and set the fundamental parameter for balllimitedness to
(51) 
with . The steepness of the decay when approaching is controlled by the parameter which we by default set to . The additional splitting in low and high frequencies according to Section 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 ), and because radial function causes 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 by in Eq. (5), with 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 cakewavelets. First, we specify an orientation distribution , 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 Chung2006 ():
(52) 
with and . The parameter determines the tradeoff between requirements 2 and 3 listed in the beginning of Section 3, where higher values give a more uniform at the cost of less directionality.
First consider setting 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, and so we set
(53) 
where the real part of the earlier found wavelet vanishes by antisymmetrization of the orientation distribution while the imaginary part is unaffected. As to the construction of , 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 , and we define
(54) 
where integration is performed over denoting the great circle perpendicular to . This transformation preserves the symmetry of , so we have . Thus, we finally set
(55) 
For an overview of the transformations see Fig. 7.
3.2 Efficient Implementation via Spherical Harmonics
In Subsection 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 in a spherical harmonic basis up to order :
(56) 
The spherical harmonics are given by
(57) 
where is the associated Legendre function, for and for and with integer and integer satisfying . For the diffusion kernel, which has symmetry around the axis we only need the spherical harmonics with , and we have the coefficients Chung2006 ():
(58) 
and Eq. (56) reduces to
(59) 
3.2.1 Funk Transform
According to DescoteauxMagnResonMed2007 (), the Funk transform of a spherical harmonic equals
(60) 
with the Legendre polynomial of degree evaluated at . We can therefore apply the Funk transform to a function expressed in a spherical harmonic basis by a simple transformation of the coefficients .
3.2.2 AntiSymmetrization
We have . We therefore antisymmetrize the orientation distribution, see Eq. (53), via .
3.2.3 Making Rotated Wavelets
To make the rotated versions of wavelet we have to find in . 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 (WignerD functions Griffiths2016 ()):
(61) 
Here and denote the Euler angles with counterclockwise rotations, where we rely on the convention . This gives
(62) 
where are the coefficients of given by
(63) 
Because both antisymmetrization and Funk transform preserve the rotational symmetry of , we have , and Eq. (62) reduces to
(64) 
The filters from this section are summarized in the following result:
Result 1
Let be a function supported mainly in a sharp convex cone around the axis and symmetricaly around the axis and as radial function of Eq. (50). Then provides our wavelet in the Fourier domain via
(65) 
with . The real part of is a tube detector given by
(66) 
The imaginary part of is an edge detector given by
(67) 
When expanding the angular part in spherical harmonics up to order and choosing :
(68) 
we have the following wavelet in the Fourier domain