Design and Processing of Invertible Orientation Scores of 3D Images.

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

Introduction

The enhancement and detection of elongated structures are important in many biomedical image analysis applications. These tasks become problematic when multiple elongated structures cross or touch each other in the data. In these cases it is useful to work with multi-orientation representations of image data. Such multi-orientation representations can be made using various techniques, such as invertible orientation scores (which is obtained via a coherent state transform) [3, 5, 10, 30, 36, 42], continuous wavelet transforms [10, 28, 30, 64], orientation lifts [13, 71], or orientation channel representations [35]. Here we focus on constructing an invertible orientation score. In order to separate the crossing or touching structures (Fig. 1), we extend the domain of the data to include orientation. This is achieved by correlating our 3D data 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 multi-orientation decomposition, we consider operators on orientation scores and process our data via orientation scores (Fig. 2).

Fig. 1

2D orientation score for an exemplary image. In the orientation score crossing structures are disentangled because the different structures have a different orientation

Fig. 2

A schematic view of image processing via invertible orientation scores

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.

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 . 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.

Fig. 3

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

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 , with a specific (semi-)analytic design for our image processing purposes.

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:

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).

We develop a new tubularity measure in 3D orientation score data. This extends previous work on tubularity measures using 2D orientation scores [19] [9, Ch. 12] to 3D data. We show qualitatively that our measure gives sharp responses at vessel centerlines and show its use for radius extraction and complex vessel segmentation in cone beam CT data of the brain.

Outline of the Article

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

Invertible Orientation Scores

Continuous Orientation Score Transform

Throughout this article we use the following definition of the Fourier transform on :An invertible orientation score is constructed from a given ball-limited 3D datasetwith ball of radius , by correlation with an anisotropic kernelThis is illustrated in Fig. 4. Here is a wavelet aligned with and rotationally symmetric around the z-axis, and the rotated wavelet aligned with given byHere is any 3D rotation which rotates the z-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 [25] for this transformation is given bywith . The function is given byand vanishes at , where the circumflex again denotes Fourier transformation. Due to our restriction to ball-limited data (2), this does not cause problems in reconstruction (5). The function quantifies the stability of the inverse transformation [25], 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 asIn practice it is best to aim for in view of the condition number of transformation given by:where M and 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 for we have -norm preservationand reconstruction Eq. (5) simplifies toWe 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 to an orientation integrated scoreIts condition number is given by .

Fig. 4

Construction of a 3D orientation score. Top: The data f are correlated with an oriented filter to detect structures aligned with the filter orientation . 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

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 norm preservation by controlling leads to restrictions on our wavelets which are captured in the following definition:

Definition 1

(Proper wavelet) Let us set a priori bounds1 . Furthermore, let be an a priori maximum frequency of our ball-limited image data. Then, a wavelet is called a proper wavelet if where is the 3D rotation around axis over angle .

If, moreover, one hasthen 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.

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 partswith Gaussian window in the Fourier domain given byAfter 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., [51, 65]. In this case we construct two scores: one for the high-frequency componentsand one for the low-frequency componentsHere we again have , as in Eq. (4). The vector transformation is then defined asFor this transformation we have the exact reconstruction formulawithAgain, 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 byThe condition number of obtained from by performing an additional splitting in low- and high-frequency components is given bythereby guaranteeing that stability is maintained after performing the splitting.

Proof

First, we find the condition number of :For the first factor in Eq. (27) we findwhere in the last step the supremum is attained by a sequence of images whose Fourier transform concentrates around the maximum of . Similarly, we get for the second factor in Eq. (27). Then we obtainSimilarly the condition number of is given byNext we express in of the original wavelet asSo it remains to quantify . For a wavelet splitting according to (18) we haveHenceAnd since for we have for satisfying (7) the following bounds on :thereby guaranteeing stability after the splitting (18). As we cannot guarantee that and M are tight bounds in Eq. (34), combining it with (29) will only give us an upper bound for the condition number in Eq. (26).

For this vector transformation we can also use the approximate reconstruction by integration (for ) over orientations. Thus we haveAs said we are only interested in processing of and not in processing of , and so we directly calculate viaFor a design with for all , we have and soThen Eq. (35) becomes

Discrete Orientation Score Transform

In the previous section, we considered a continuous orientation score transformation. In practice, we have only a finite number of orientations. To determine this discrete set of orientations we uniformly sample the sphere using a simple2 electrostatic repulsion model [18].

Assume we have a number of orientations , and define the discrete invertible orientation score byThe exact reconstruction formula is in the discrete setting given bywith the discrete spherical area measure () which for reasonably uniform spherical sampling can be approximated by , andFor spherical samplings that are constructed by triangularization (such as tessellations of the icosahedron), one can rely on surface areas of spherical triangles to compute more accurately. See for example [29, Eq. 83]).

Again, an exact reconstruction is possible iff and we have norm preservation when =1. Again, for the wavelets for whichthe 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 to an orientation integrated scoreThis 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 ball-limited 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 where is a 3D rotation around axis over angle .

If, moreover, one hasthen we speak of a proper wavelet with fast reconstruction property, cf. (43).

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

Steerable Orientation Score Transform

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

Wavelet Design Using a DFT

A class of 2D cake wavelets, see [10, 27, 37], was used for the 2D orientation score transformation. We now present 3D versions of these cake wavelets. 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]:

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 (Definition 2).

The kernel should be strongly directional.

The kernel should be separable in spherical coordinates in the Fourier domain, more explicitly , with Because by definition the wavelet has rotational symmetry around the z-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-valued orientation score, as the wavelets are complex-valued. For an intuitive preview, see the boxes in Fig. 7.

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: and evaluated on a grid of pixels

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 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,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 g and set the fundamental parameter for ball limitedness towith . 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 Sect. 2.1.1 effectively causes a splitting of the radial function, see Fig. 5.

Fig. 5

Radial part g of , see Eq. (50) and radial parts and of and after splitting according to Sect. 2.1.1. The parameter controls the inflection point of the error function, here . The steepness of the decay when approaching is controlled by the parameter with default value . At what frequency the splitting of in and is done is controlled by parameter , see Eq. (18)

In practice the frequencies in our data are limited by the Nyquist frequency (we have ), and because radial function g 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 cake wavelets. 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 [20]:with and . The parameter determines the trade-off between requirements 2 and 3 listed in the beginning of Sect. 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,3 and so we setwhere the real part of the earlier found wavelet vanishes by anti-symmetrization of the orientation distribution A while the imaginary part is unaffected. Note that the right-hand side of (52) and (53) is the same for all values of . As to the construction of , 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 definewhere integration is performed over denoting the great circle perpendicular to . This transformation preserves the symmetry of A, so we have . Thus, we finally setFor an overview of the transformations, see Fig. 7.

Fig. 6

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: and evaluated on a grid of pixels (Color figure online)

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.

Fig. 8

Coverage of the Fourier domain before and after splitting according to Sect. 2.1.1. Left: The different wavelets cover the Fourier domain. The “sharp” parts when the cones reach the center, however, cause the filter to be non-localized, which was solved in earlier works by applying a spatial window after filter construction. Right: By splitting the filter in lower and higher frequencies we solve this problem. In the figure we show for the different filters, before applying the Funk transform to the orientation distribution A

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 up to order L:The spherical harmonics are given bywhere is the associated Legendre function, for and for and with integer and integer m satisfying . For the diffusion kernel, which has symmetry around the z-axis we only need the spherical harmonics with , and we have the coefficients [20]:and Eq. (56) reduces to

Funk Transform

According to [23], the Funk transform of a spherical harmonic equalswith 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 .

Anti-symmetrization

We have . We therefore anti-symmetrize the orientation distribution, see Eq. (53), via .

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 (Wigner-D functions [41]):Here and denote the Euler angles with counterclockwise rotations, where we rely on the convention . This gives where are the coefficients of h given byBecause both anti-symmetrization and Funk transform preserve the rotational symmetry of A, we have , and Eq. (62) reduces toThe 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 z-axis and symmetrically around the z-axis and g as radial function of Eq. (50). Then A provides our wavelet in the Fourier domain viawith . The real part of is a tube detector given byThe imaginary part of is an edge detector given byWhen expanding the angular part in spherical harmonics up to order L and choosing :we have the following wavelet in the Fourier domainand the coefficients of A and relate viaWe obtain rotated versions of our filter viawith .

As we do not have analytical expressions for the spatial wavelets , we sample the filter in the Fourier domain using Eq. (71) and apply a DFT afterward. The wavelet is a proper wavelet with fast reconstruction property (Definition 1).

Remark 3

The heat kernel on is given bywhere 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 (e.g., for ).

Stability of the Discrete Transformation with Fast Reconstruction for Filters of Result 1

To make a fast reconstruction by summation possible (requirement 2), we need a proper wavelet set with the fast reconstruction property (Definition 2) with . We now focus on finding bounds for such that we can choose our parameters in a deliberate way.

Proposition 1

Let be a set of wavelets constructed via the procedure in Result 1. Then we have bounds on given bywith , and here the norm is the -norm on .

Proof

First we expand function in spherical harmonics:We have for , but we still need to quantify the angular part. We define , so thatThis varying component should remain small. We use the Cauchy-Schwarz inequality for each order l:from which (73) follows.

See Fig. 9 for visual inspection of bounds of and , and numerical results for the bounds of .

Fig. 9

Inspection of the stability of the transformation for different values of given an orientation distribution and for . Left: Spherical plot of A and the angular part of polar separable function and . Orientation coverage is more uniform as the plots for and look more like a ball. Right: The upper and lower bounds of . 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 ( orientations). Furthermore, we show and (orange dashed lines) for (Color figure online)

Corollary 1

Given our analytical bounds (73) from Proposition 1 and , 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 when choosing parameter . In practice we have a proper wavelet set with fast reconstruction property already for smaller values of (see Fig. 9).

Wavelet Design with Continuous Fourier Transform and Analytical Description in the Spatial Domain

In the previous section we described wavelets which were analytical in the Fourier domain and were sampled and inverse discrete Fourier transformed to find the wavelets in the spatial domain.

To get more control on the wavelet properties in both the spatial and Fourier domain, it would be convenient to have an analytical description of the wavelets in both domains. This could be achieved by expressing the wavelets in a basis for which we have analytical expressions for the Fourier transform. We will now discuss 2 such options for the basis and describe filters expressed in them.

A Review on Expansion in the Harmonic Oscillator Basis

The first basis in which we could expand our wavelets are the eigenfunctions of the harmonic oscillator , 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 () and eigenfunctions of H are also eigenfunctions of . We then expand our wavelets in these eigenfunctions restricting ourselves to eigenfunctions which are symmetric around the z-axis: the spherical harmonics with . The wavelet is then given bywith the spherical harmonics, and spherical coordinates for and , respectively, i.e.,and given bywhere is the generalized Laguerre polynomial. We then choose the case with least radial oscillations . If we then choosewe have that approximates 1 for all in the Fourier domain as and we get the following wavelet4 [27]:For this wavelet we have an analytical description in both spatial and Fourier domain.

Expansion in the Zernike Basis

The wavelets from the previous subsection have some unwanted properties such as poor spatial localization (long oscillations) and the fact that the wavelets maximum does not lie at the wavelets center, see Fig. 10. A possible explanation is that the basis used is orthogonal on the full 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 , recall Eq. (2).

Fig. 10

Wavelet expanded in the harmonic oscillator basis according to Eq. (81) for . Left: 3D visualization showing one negative (blue) and one positive (orange) iso-contour. Right: Cross section of the wavelet at (Color figure online)

The 2D Zernike basis is often used in applications as optics, lithography and acoustics [1, 12, 15], since efficient recursions can be used for calculating the basis functions and analytic formulas exist for many transformations among which the Fourier transform. The basis is therefore highly suitable for problems such as aberration retrieval where a wave function in the Fourier domain should be estimated from measurements in the image domain [24].

Orthogonal polynomials of several variables on the unit ball were considered by Louis [49] in the context of inversion of the Radon transform for tomographic applications. For a modern treatment of orthogonal polynomials on the unit ball, see [33, Ch. 4]. Here we will use the generalized Zernike functions [43], which vanish to a prescribed degree at the boundary of the unit ball, with explicit expansion results for particular functions supported by the unit ball. Since this basis is orthogonal on the unit ball and has explicit results for the Fourier transform it is highly suitable for the application in mind. Derivations for the results used in the following section can be found in [43].

The 3D Generalized Zernike Basis

The generalized Zernike functions are given bywith spherical coordinatesinteger such that , integer and and . The angular part is the spherical harmonic function and the radial part is given bywhere denotes the Jacobi polynomial. The generalized Zernike functions are orthogonal on the unit ballwith the Kronecker delta and with normalization factorin which is the (generalized) Pochhammer symbol.

Fourier Transform

The inverse Fourier transform of the generalized Zernike functionis given bywith andHere and are the Bessel functions and spherical Bessel functions [55]. For integer , the expression in Eq. (89) for reduces toExpansion of Separable Functions

An additional constraint for the wavelets is that they should be separable in the Fourier domain, i.e., . When expanding such a function in the generalized Zernike basis,we can split the coefficients in radial coefficients and angular coefficientswhere The coefficients 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 havewhile for all For each l, the radial functions with n varying are a basis for functions defined on the interval [0, 1]. For separable functions, we expand the same radial function for each l, and it can be shown that there is a recursion formula for the radial coefficients [43].

Wavelets

We now choose appropriate radial and angular functions for our wavelets expressed in the generalized Zernike basis.

Angular Function for the Zernike Wavelets

For the angular functions we again choose orientation distribution 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 and . In this section we design a radial function for which is relevant for further processing. Furthermore, we already have an analytical description for , which we set to (see Eq. (37)).

The radial function of should therefore approximate on the interval and should smoothly go to zero when approaching the edges of the interval. For the moment, we set and we include the scaling later. To start, we define the functionsee Fig. 11a for the case . For this function we have the following coefficientsTo obtain a flatter function we multiply the function with a second-order Taylor expansion of the reciprocal function around the function’s maximum obtained atsee Fig. 11b. The resulting function is again a sum of functions of type (97) with different values for , so we can find the coefficients for the flattened function as well. For the specific case we get the following flattened functionwith . For this flattened function the coefficients are given bywith the coefficients of , and in the second-order Taylor series of the reciprocal. These coefficients follow from (100) and are given byThe filters from this section are summarized in the following result:

Fig. 11

Left: Function . Right: Flattened function which are obtained from by multiplying with the second-order Taylor approximation of its reciprocal around the function maximum:

Result 2

(Analytic 3D-wavelets in Zernike basis) Let and let be a function supported mainly in a sharp convex cone around the z-axis and symmetrically around the z-axis. Then A provides our wavelet in the Fourier domain via Eq. (65). The real part of is a tube detector and the imaginary part of is an edge detector, see Fig. 7. We choose radial function in Eq. (18) for and angular function and expand in the generalized Zernike basis:The coefficients follow by expanding A in spherical harmonics and in the radial Zernike polynomials, recall Eqs. (95) and (96). This yields (Eq. (58)) and (Eq. (101)) and coefficients :The spatial wavelet is given bywith given by Eq. (89) and the spherical harmonics of Eq. (57). Then we obtain rotated filters viaSince 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).

Experiments with the Filters and the Transform Between a 3D Image and Its Orientation Score

Before considering applications of the filters in the next section, we first compare filters obtained by DFT (Sect. 3) to filters expressed in the generalized Zernike basis (Sect. 4) and inspect the quality of the reconstruction.

Comparison of Wavelets Obtained via DFT and Analytical Expressions Using the Zernike Basis

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

Fig. 12

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 . Right: The low-pass filter. Top Filters according to Result 1 with parameters and . Bottom: The filters according to Result 2 with and . Both have and are evaluated on a grid of voxels (Color figure online)

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.

Fig. 13

Comparison of construction and reconstruction of data A.1 using the different types of filters with the same settings as in Fig. 12. In each row, from left to right, an iso-contour of the data and 3 slices through the center of the data along the three principal axis. Top: The original data. Middle: the data after construction and reconstruction using the filters from Result 1. Bottom: the data after construction and reconstruction using the filters from Result 2

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.

Fig. 14

Overview of the datasets used in our experiments

Table 1

Default values for the parameters of the orientation score transform used in the application section

Parameter	Default value	Defining Eq.
N_o	42	(39)
γ	0.85	(51)
σ_erf	13(ρ_N-ϱ)	(50)
s_ρ	12(16)^2	(18)
s_o	12(0.45)^2	(52)
Discrete wavelet size	11×11×11	–

Diffusion Via Invertible Orientation Scores

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.

CEDOS

We now use the invertible orientation score transformation to perform data enhancement according to Fig. 2. Because 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 settingwith any rotation for which . The operators which we consider are scale spaces on SE(3) (diffusions), and are given by withHere is the solution of a nonlinear diffusion equation:where in coherence-enhancing diffusion on orientation scores (CEDOS) is adapted locally to initial condition based on exponential curve fit (see [32]), and with the left-invariant vector fields on SE(3). The diffusion is better understood in locally adaptive frame . Here follows from an exponential curve fits and points along our structure. and span the plane spatially perpendicular to our structure, and and correspond to angular diffusion (embedding in SE(3) leads to a third angular dimension). Our diffusion then takes the diagonal formwhere we limit ourselves to diffusion of type , and to preserve the data symmetry of Eq. (107). The diffusion system in Eq. (110) has been set up in previous work [45] with constant and .

In this work we include further adaptivity by making the constants and data adaptive. We aim to enhance oriented structures and reduce noise as much as possible. Therefore, non-oriented regions are smoothed isotropically by settingwith a constant automatically set to the quantile of and a measure for orientation confidence given by minus the Laplacian in the space orthogonal to the structure orientation :Since we want to stop diffusion when reaching the end of a structure, we setwith a constant automatically set to the quantile of .

We then obtain Euclidean invariant image processing viawhich includes inherent projection of orientation scores, even if maps outside of the space of orientation scores. We write because we extend the inverse to .

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 denotes the standard deviation of the difference signal , 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 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 and background regions given by :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 we estimate the noise standard deviation byThe contrast-to-noise ratio (CNR) is then given bywhere the numerator denotes the contrast in our data.

Fig. 15

Selected regions for determining the contrast-to-noise ratios for dataset B.1. Left: Grid containing slices of the data (top row), the same slices with segmented vessel parts (second row) and the slices with three selected background regions (third row), the slices after applying CEDOS. Right: 3D visualization of the data

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 and were determined using Eqs. (111), (113), and we set . For the orientation score transformation, we used and . For CED we used the following settings: and c the -quantile of (see [69]).

Fig. 17

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 , where and 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

As one can expect, in all cases we recognize a peak since initially noise is reduced whereas later also contrast is reduced. Compared to Gaussian diffusion and CED, we reach a higher CNR (Fig. 16). Furthermore, in the CEDOS and CED case, the CNR does not decrease much when applying more diffusion making it more robust with respect to choice in diffusion time. The fact that CED does not achieve high CNR ratios and performs relatively bad in this test is that the diffusion matrix is not designed to reduce background noise (which is used here to quantify noise) but mainly to enhance orientated patterns; for this reason we also set relatively high to still achieve noise reduction in the background regions.

Fig. 16

Contrast-to-noise ratio (CNR) against diffusion time for CEDOS compared to Gaussian regularization and CED of data B.1 depicted in Fig. 15. The times denoted by a, b and c correspond to the diffusion times shown in Fig. 17

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).

Fig. 18

Volume rendering of the diffusion results for datasets B.1, B.2 and B.3 (from top to bottom) visualized in the Philips viewer [62] using default settings in all cases. For all cases we used optimized diffusion time according to Fig. 16

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 (Eq. 111) and therefore we smooth only along the vessel.

Fig. 19

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)

Tubularity Measure

In this section we propose a tubularity measure based on the edge information in our orientation scores. This tubularity measure is then used for vessel width measurements and could be used for vessel segmentation.

Tubularity measures are designed to have a high response on the centerline of tubular structures. One such tubularity measure is the image gradient flux filter [14] which is used for vessel segmentation. An extension of the gradient flux filter is the optimally oriented flux filter [47] which introduces the notion of oriented flux making the filter orientation sensitive.

A 2D tubularity measure based on orientation scores was proposed in [9]. The advantage of this tubularity measure is that it included nonlinearity and that the implementation via orientation scores still has a response at crossing vessels. Here we propose an extension of this tubularity measure to 3D making use of 3D orientation scores.

Tubularity Via Orientation Scores

For the tubularity measure we detect edges in the plane perpendicular to orientation . Within this plane, the product of two opposite edges at radius r and in-plane orientation at position is given bywhere and , with an orthogonal basis for the orthogonal complement of . The product of the two edge responses in Eq. (120) yields a better performance than taking the sum, as is done in [9, Fig. 12.2] and [19]. The idea behind taking the product instead of the sum is that we need a high edge response in both directions. See Fig. 20 for a schematic visualization. Since for a real tube all edges should be present, and we do no want any response for, e.g., plate structures, we take a minimum over the perpendicular orientations parametrized by . This is done as followsFor the angular regularization kernel we use the simple periodic 1D-diffusion kernel

Fig. 20

Schematic visualization of the edges used in the tubularity measure . Left: A 3D iso-contour visualization of a vessel with orientation . The coordinates are polar coordinates for the plane perpendicular to orientation spanned by and . Right: In this plane opposite edges are multiplied in . The edge is expected to have outward orientation

over , instead of the true diffusion kernel (72), where for practical values we use the series of rapidly converging kernels can be reasonably truncated already at the first term .

From the tubularity measure (121), we extract the following features: Here is the tubularity confidence which is a measure for how certain we are at least one tubular structure is present at position . The features and are the orientation and radius of optimal tubularity response at position .

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: and we discretized the -integral in Eq. (121) using 8 orientations.

As validation we compared the optimal radius to the ground truth radius and inspect the tubularity confidence. The transversal profiles in the tubes are modeled by Gaussians, and we define the ground truth radius at the standard deviation 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.

Fig. 21

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 (max of tubularity over radius and orientation) with radius of max response in color. Right: Measured radius against ground truth radius on the ground truth centerline. The opacity of the plotted points is linearly scaled with the tubularity confidence

Fig. 22

Measured radius 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

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: and we discretized the -integral in Eq. (121) using 12 orientations. For the orientation score transformation we used and evaluated on a grid of 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 and a simple segmentation given by the 0-iso-contour of the distance map , where are the positions given by the quantile of highest responses.

Fig. 23

Tubularity on real data. Top: Dataset A.1. Bottom: Dataset A.2. Left: The data. Middle left: The projected tubularity (max over radius and orientation) colored by radius. Middle right: Orientation of max response for 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 , where are the positions given by the 1% quantile of highest responses

Conclusion

We presented theory and filters for the 3D orientation score transformation which is valuable in handling complex oriented structures in 3D image data. Then we showed applications of this transformation.

First, we proposed filters for a 3D orientation score transform. We presented two types of filters: The first uses a discrete Fourier transform and the second is designed in the 3D generalized Zernike basis which allowed us to find analytical expressions for the spatial filters. Both filters allowed for an invertible transformation. The filters and the quality of the reconstruction are assessed in Sect. 5, where we showed that the discrete filters approximate their analytical counterparts well. We also verified the invertibility of our transformation by showing data reconstructions of real medical data.

The orientation score transform was then used in two different applications. In the first we presented an extension of coherence-enhancing diffusion via 3D orientation scores which we applied to real 3D medical data and showed our method effectively reduced noise while maintaining the complex vessel geometry. In the second application we propose a new nonlinear tubularity measure via 3D orientation scores. The tubularity measure has sharp responses for vessel centerlines, and we showed its use in radius extraction and complex vessel segmentation.

In this work basic applications of the tubularity measure are shown. Future work would include using the tubularity measure in more advanced vessel segmentation procedures. Furthermore, many other applications exist for the 3D orientation scores and many techniques developed for 2D orientation scores can now be extended to 3D. First steps are presented in this paper, where the extension of 2D CEDOS and a 2D tubularity measure via 2D orientation scores are given. It is also interesting to explore the nonlinear diffusion procedure (Eq. (110)) for contextual processing of diffusion-weighted MRI and to compare with existing approaches [58, 68].

Symbol	Explanation	Reference
B.1 Spaces and Input Data
R^3×S^2	Space of positions and orientations	Page 1
L2ϱ(R^3)={f∈L_2(R^3)|supp(Ff)⊂B_ϱ}, with ϱ>0	Space of frequency ball-limited images	(2)
B_ρ={x∈R^3|‖x‖<ρ}, with ρ>0	The ball around the origin with radius ρ	(2)
B.2 Orientation Score Transformation
W_ψ	Orientation score transformation	(3)
Wψ-1	Data reconstruction via the inverse orientation score transformation	(5)
ψ	Wavelet used for the orientation score transformation	(3)
ψ^	Fourier transform of the wavelet used for the orientation score transformation	(6)
M_ψ	Factor used to quantify the stability of the transformation	(6)
N_ψ	Factor used to quantify the stability of the transformation when using the simplified reconstruction by integration	(12)
Wψd,(Wψd)^-1,Mψd,Nψd	Discretized versions of W_ψ,(W_ψ)^-1,M_ψ,N_ψ	(39), (40), (41), (42)
dσ and Δ_i	Spherical area measure and discretized spherical area measure	(5), (40)
B.3 Coordinates
x=(x,y,z)	Cartesian coordinates real space	-
ω=(ω_x,ω_y,ω_z)	Cartesian coordinates Fourier domain	-
(r,θ,ϕ),x=(rsinθcosϕ,rsinθsinϕ,rcosθ)	Spherical coordinates real space	(49)
(ρ,ϑ,φ),ω=(ρsinϑcosφ,ρsinϑsinφ,ρcosϑ)	Spherical coordinates Fourier domain	(49)
B.4 Wavelets
g	Radial function of the cake filters	(50)
A	Orientation distribution used in wavelet construction	(52)
Ylm	Spherical harmonics	(57)
Zn,lm,α	3D generalized Zernike functions	(82)
Rnl,α	Radial part of the 3D generalized Zernike functions	(84)
Sn,lα	Radial part of the inverse Fourier transform 3D generalized Zernike functions	(89)
B.5 Applications
V(x,n,r)	Tubularity measure	(121)
s^t(x)	Tubularity confidence	(122)
n^∗(x)	Orientation of maximum tubularity response	(122)
r^∗(x)	Radius of maximum tubularity response	(122)