Shape-Adaptive DCT for Denoising of 3D Scalar and Tensor Valued Images

Diffusion Tensor Imaging

Tensor-valued data occur in many branches of mathematical sciences; see, e.g., Weickert and Hagen²⁰. In this paper, the tensor-valued data comes from diffusion tensor MRI of the human brain. From a set of K direction sensitive magnetic resonance images a symmetric positive definite tensor is constructed in each voxel of the image domain. This matrix yields structural information of the tissue in each voxel.

The relationship between each direction weighted measurement and the diffusion tensor D is given by the Stejskal–Tanner equation ^21,22

2

where b is a positive scalar given by the measuring pulse sequence, and is one of the predefined selected directions for which measurement S_k is acquired. From K direction weighted measurements we obtain K equations that we use for estimating the six unknowns of the diffusion tensor D. This can be done, for example, by a linear least-squares method, or other more adaptive methods.⁵ We note that as the transformation (Eq. 2) is nonlinear, we do not know the distribution of the noise in each element of the tensor D. Hahn et al.²³ have studied how noise propagates through the estimation process.

In structured tissue such as in the heart muscle or in the white matter of the brain, the self-diffusion of water is highly anisotropic. In gray matter and in cerebrospinal fluid, the self-diffusion of water is almost isotropic. Based on knowledge of the diffusion tensor D, a model of the myelinated nerve fiber pathways in the white matter can be constructed via fiber tracking algorithms.^1,2,24,25

The quality of the estimated diffusion tensor depends on several parameters. One particular parameter is the number of acquisitions or excitations (NEX) performed in each of the diffusion sensitizing directions. A high number of acquisitions in each direction gives a tensor estimate of good quality, provided the patient does not move during the acquisitions, whereas a small number of acquisitions give a shorter examination time for each patient. Therefore, we have a compromise between efficiency of the acquisitions and quality of the resulting images. In this paper, we investigate the possibility of post-processing the data from a small number of acquisitions, and still being able to construct tensor estimates of high quality. Thus, a practical research goal is to decrease the scanner time required for each patient.

Shape-Adaptive Discrete Cosine Transform

The 2-dimensional discrete cosine transform (DCT) is extensively used in image science. In its original formulation, it transforms a quadratic region in the spatial domain into the frequency domain. Being a harmonic transform, the DCT has a compactification property, i.e., good approximations of the image can be constructed by employing only a few of the coefficients in the frequency domain.^26,27 However, when the image domain of the transform contains sharp edges and only a few of the coefficients in the frequency domain are employed for the reconstruction to the spatial domain, various artifacts such as smearing of edges and Gibbs phenomena occur. To avoid these artifact, the region should be as homogeneous as possible. This is achieved by replacing the static regions from the standard DCT by regions Ω_x, which adapts to the information in the image around a point x. We choose these regions in such a way that the data can be well approximated by a smooth, slowly varying function. Such a function is well approximated by few coefficients from the frequency domain. The regions Ω_x should ideally not contain any discontinuities.

In a series of papers, Katkovnik, Foi, Egiazarian, Astola, and others describe shape adaptive DCT (SA-DCT) for denoising of 2D grayscale and color images.^16–19 The algorithm can be divided into three different stages: (1) construction of an adaptive neighborhood for each point in the domain, (2) transformation and thresholding of each neighborhood, and (3) estimation of the noise-free image. The adaptive neighborhoods are constructed by local polynomial approximations (LPA) in combination with the intersection of confidence intervals (ICI) rule. The transformation of each neighborhood to the frequency domain is done by a DCT algorithm, and hard thresholding is applied on the coefficients in the frequency domain. The inverse DCT algorithm is then applied. This results in a denoised region. Since each pixel x has its own region Ω_x, and in general these regions may overlap, we get an overcomplete basis. This overcomplete basis is used to construct the final image by weighting the basis elements in a proper way.

The state-of-the-art results obtained by the SA-DCT methods in 2D as well as their efficiency makes them attractive for denoising of 3D scalar valued images and 3D matrix valued images. In this paper, we extend the framework of SA-DCT to both 3D scalar valued and 3D matrix valued images.

Shape-adaptive Discrete Cosine Transform

Let be a closed spatial domain of dimension N and denote the noisy data set, which is discretized on a uniform grid. In the rest of this, paper we restrict the attention to 2D and 3D data sets, i.e. . We refer to the resulting denoised data set as I⁻. Presently, we treat the standard deviation σ of the noise of I as an input parameter.

A main ingredient of the SA-DCT method is the adaptive neighborhood Ω_x surrounding each voxel x∈Ω. The idea is that this neighborhood should contain voxels that in some way are “similar”, or homogeneous. A neighboring point y∈Ω can either have an intensity I(y), which is close to the intensity I(x), or the intensities can differ substantially. In the case where I(y)≈I(x), we want to include the point y in the adaptive neighborhood of x, i.e., y∈Ω_x. To decide which voxels should belong to the adaptive neighborhood of a given point, we use local polynomial approximations (LPA) and the intersection of confidence intervals (ICI) rule.¹⁹

To construct the adaptive neighborhood, we consider a set of directions such that each component of θ_i is either −1, 0 or 1, but never all equal to zero. It follows that there must be 3^N −1 such directions in an N-dimensional data set. In 2D,¹⁸ there are eight such directions: the four cardinal and the four intermediate compass directions. In 3D, there are 26 unique directions following a similar pattern.

We span a star-shaped skeleton around each point x in the image domain by tracing the voxels along straight lines in the directions of θ_i. The length d_i corresponding to the straight line in the direction of θ_i in the skeleton is determined by the ICI algorithm. We close the skeleton such that it becomes a polygonal hull by joining neighboring endpoints of the vertices in the skeleton by line segments (in 2D) or triangles (in 3D). We denote the domain inside this closed polygonal hull by Ω_x. For each voxel in the image domain, such an adaptive neighborhood is constructed. The pseudocode for our SA-DCT denoising is given in Algorithm 1.

In the following section, we explain how we can use the LPA-ICI method to compute the length d_i of each branch in the star .

LPA-ICI

To span the region , we calculate the support of each branch in the star, i.e., number of voxels that should be included along each direction vector θ_i, i = 1,...,3^N − 1. The idea is that the voxels in should have intensity values that are close to the intensity value of the center voxel x. Variations in the included data should be caused by the noise level and small local variations, and not by edges in the image.

To achieve this, we filter each direction with LPA kernels of varying length (scale) h∈H={h₁,...,h_j}, where h₁<h₂<...<h_J. The generation of these filter kernels are described in Foi 2005.¹⁶ For each kernel g^(h) containing weights , where i = 1,..., h, we have the property that the center voxel x has the highest weight . In addition, the weights sum to 1 and decrease with the length of the filter.

We can consider this filtering as a convolution of the data with a filter kernel of varying length. When the kernel is applied to the voxels in direction θ_i, we get the filtered value

3

The standard deviation of the noise in μ^(h) is given through the relation

4

For each direction, we then get the confidence intervals

5

where Γ > 0 is a global parameter of the algorithm. A large Γ results in a large noise tolerance, and more voxels will be included in the regions, and vice versa.

The ICI rule states that along the direction vector θ_i we should choose the largest distance d_i∈H where we have intersection of all the confidence intervals; see Figure 1. More precisely,

6

Having determined the length of each branch in the star- shaped domain , we define the neighborhood Ω_x as all voxels inside the polygonal hull closing with branches d_iθ_i, where i = 1,...,3^N − 1. By construction, the intensities in this region should not contain large changes caused by edges in the image. The noise in this region can now easily be removed by thresholding small coefficients in the frequency domain. We use the discrete cosine transform for this purpose as described in the next subsection. In the upper left panel of Figure 2, we have visualized an 2D example of the output of the LPA-ICI algorithm. The arrows centered around voxel x indicate each of the eight directions θ_i employed by the algorithm and the length of each arrow indicates the distance d_i. The set of all voxels covered by an arrow in the figure is and all voxels inside the dashed polygonal hull spanned by the arrows is Ω_x.

Note that since we only perform a LPA-ICI estimation on the voxels that coincide with the skeletonized domain , we do not have direct control over the intensity values in the set . It has been shown that for scalar images, this approach is a good compromise between efficiency and accuracy.¹⁷

Before the region Ω_x is denoised in the frequency domain, the mean of the region is estimated and subtracted from each intensity in Ω_x. This process is referred to as DC separation, and it has been shown to improve certain pathological issues associated with the shape-adaptive DCT described in the next section.¹⁷ After the DCT denoising is complete and the coefficients have been transformed back to the spatial domain, the precomputed mean is added to shift the mean of the intensities back to the same level as before the DCT. Although this shift of intensities is not entirely justified from the approximation standpoint (since the precomputed mean also will contain noise), it has been shown to visually provide superior results with few adverse effects.¹⁷

The DCT Algorithm

The discrete cosine transform (DCT) is used extensively in signal and image processing. The one-dimensional DCT of a signal {Z₀,...,Z_M−1} of length M is defined as

7

for k = 0,..., M − 1, where and . Note that this transform can be expressed as a matrix–vector product

8

When c_k is defined as above, then the A matrix is orthogonal. This implies that the inverse DCT can be expressed as

9

Two- and three-dimensional DCT are usually achieved by successively applying the one-dimensional DCT along the coordinate axes (i.e., separability). However, note that Ω_x will in general not be rectangular. Sikora has developed an algorithm for discrete cosine transform on nonrectangular domains.²⁶ In this paper, we employ this algorithm. However, we still use the orthogonal transform (Eq. 7) and not the one originally presented in Sikora 1995.²⁶

In the following, we let denote the quadratic (in 2D) or cubic (in 3D) null-extension of Ω_x. Null-extension in this context means a padding of the region with a particular value, “null”. This value will only be used as a “place-holder”, and not in any direct calculations. Accordingly, it is never considered a coefficient in the discussion below. The purpose of the null-extension is to extend the size of the region into a more manageable form, as later described.

Note that when examining Ω_x along the coordinate axes, it may be noncontiguous. One approach to alleviate this problem could be by zero-padding the region instead. However, this causes problems when applying the traditional DCT algorithm as many components of the DCT domain will be needed to represent these high jumps in intensities introduced by zero-padding. Sikora’s approach avoids this problem by first shifting all nonnull values of along the first coordinate axis so that they become consecutive in . A one-dimensional DCT, where the length of the signal M is equal to the number of nonnull values, is then applied to the shifted data. The same procedure is then applied to each dimension in turn, by first shifting the data and then applying the 1D DCT. When applying the inverse DCT, we need to invert these shifts, so a record of the rearrangements must be maintained.

In Figure 2 we display a 2D example of how the null-padded voxels are shifted to produce consecutive values. The upper right panel shows the null-padded region, where black is used to indicate voxel intensities (and later DCT coefficients) and white is used to show null values. The lower left panel shows the intensities shifted in the x direction towards the origin in the lower left corner. A 1D DCT is then applied to the nonnull elements of each row containing nonnull values; first to row 2 containing M = 1 intensities, then to row 3 containing M = 4 consecutive intensities, and so forth. The lower right panel shows the results after the row-wise DCT, shifted along the y axis. Again a 1D DCT is applied, this time to each nonnull element of each column. When the data is 3D, this procedure has to be repeated again in the z-direction. The final configuration is a set of shifted DCT coefficients in the origin corner of the cube.

Thresholding in the DCT Domain

Let denote the domain transformed from Ω_x using the DCT algorithm described in the previous section and let denote a given coefficient in . In addition, let denote the number of coefficients in the neighborhood . The cutoff threshold f is given as¹⁹

10

and the hard thresholded coefficients are given as

11

for all .

The thresholded region is then transformed and shifted back into the spatial domain by the inverse DCT giving in the spatial domain.

Estimation from Overcomplete Basis

Notice that since we calculate a region Ω_x for every voxel in the image, we have extensive region overlap, i.e., we have an overcomplete basis. To reconstruct an image from this information, we weight the data together. We assign a weight to every region and use the information in overlapping regions to estimate the denoised image. It is a standard approach to use weights that are inversely proportional to the mean variance of the region. However, for adaptive regions this leads to oversmoothing.¹⁷ To compensate for this, we can divide the weights by the square of the size of the region.

The mean variance of the region is given by

12

where is the number of nonzero coefficients in . This gives the following weights for the regions

13

The regions can now be weighted together, giving the final recovered estimate using the relation

14

for all p∈Ω and where the sum is taken over all voxels x such that contains p.

Generalization to DTI

In the previous section, we extended the SA-DCT methods from 2D to 3D images. This code can be used directly for denoising of 3D MRI images. In this section, we will show how the code also can be applied for matrix-valued DTI images.

In a previous work, two of the authors have investigated total variation (TV) regularization of tensor-valued data,¹¹ where the estimated tensor is regularized. To ensure positive definiteness of the regularized tensor, it is represented implicitly by Cholesky factorization as D=LL^T, where L is a lower triangular matrix. In the present work, we have adopted this approach when regularizing tensor-valued data by SA-DCT methods. Thus, we apply the 3D algorithm to each of the elements of L.

We are aware that there are other ways to apply the 3D code on matrix-valued DTI images, e.g., we could denoise directly on the S_i images. However, the advantage of the presented method is that we are able to guarantee positive eigenvalues. In addition, we have conducted numerical tests that indicate that the performance is similar for both methods.

3D Scalar-valued Data

We have in this paper generalized the SA-DCT methods from 2D to 3D images. In the first example, we want to show the difference between the 2D SA-DCT algorithm applied slice by slice in a 3D data set, referred to as quasi-3D, and application of our genuine 3D algorithm. We use 3D data from the BrainWeb, a database of freely available semirealistic simulated MR images.²⁸ The true data set has a range from 0 to 1 and a mean of 0.26. The added noise is normal distributed with zero mean and a variance of 0.07.

In Figure 3 we show a coronal slice of the image comparing the performance of the quasi-3D SA-DCT algorithm and the genuine 3D version. Line artifacts between neighboring slices can be observed, since the 2D algorithm ignores information in the z-direction. An additional problem with the 2D approach, dealing with isotropic voxels, is determining which axis to slice across. We have arbitrarily chosen the z-direction, but a choice of x or y would in this example yield different suboptimal results. As we observe from Figure 3, the result from the full 3D algorithm proposed in this paper yields results in which the noise reduction is performed in a consistent manner across all three dimensions. In Figure 4 we show a zoom-in of a small portion of the result from Figure 3.

We define the error as the Euclidean 2-norm of the difference between the noise-free image and the denoised image. The error between the noisy image and the noise-free image was 140.13. The image denoised using the quasi-3D algorithm gave 47.60, and using the full 3D algorithm reduced the error to 41.92.

3D Tensor-valued Synthetic Data

The main motivation behind this paper is denoising of tensor-valued images, in particular diffusion tensor MR images. First, we denoise a synthetic DTI data set where the object is a simulated torus. The DTI torus has been generated using the software Teem, written by Gordon Kindlmann.²⁹ For visualization of the color-coded fractional anisotropy (FA) images, derived from the estimated tensor images, we have used DtiStudio developed by Susumu Mori and coworkers.³⁰

Our DTI phantom consists of a doughnut-shaped object with cigar-shaped diffusion in the direction of the main circumference. The baseline image S₀ is constant equal to 1 and the six direction sensitive measurements S₁,...,S₆ have the range 0 to 0.80. We added normal distributed noise with zero mean and a variance of 0.01 to S₀ and variance of 0.04 to S₁,...,S₆. The six gradient directions g₁,...,g₆ used are given by the columns of the matrix

The diffusion tensors computed from the clean DTI data is used as a reference and the denoised tensors are compared against this ground truth. We define the error as the sum of squared element-wise tensor differences

15

where D_ij(p) denotes the tensor element in position (i,j) at voxel p. The global error in the noisy data was found to be 95.86.

When first computing the diffusion tensors from the noisy DTI data and then apply the 3D SA-DCT algorithm to each of the six elements of the voxel-wise lower triangular L decomposition of the tensor, the global error was 18.03. The results of this denoising procedure are shown in Figure 5.

3D Tensor-valued Real Brain Data

Finally, we tested SA-DCT denoising on real diffusion tensor images from a healthy human brain. The human subject data were acquired using a 3.0-T scanner (Magnetom Trio, Siemens Medical Solutions, Erlangen, Germany) with a eight-element head coil array and a gradient subsystem with the maximum gradient strength of 40 mT m⁻¹ and maximum slew rate of 200 mT m⁻¹ ms⁻¹. The DTI data were based on spin-echo single-shot Echo Planar Imaging (EPI) acquired utilizing generalized auto calibrating partially parallel acquisitions (GRAPPA) technique with acceleration factor of 2, and 64 reference lines. The DTI acquisition consisted of one baseline EPI, S₀, and six diffusion weighted images S₁,...S₆ (b-factor of 1,000 s mm⁻²) along the same gradient directions as in the previous example. Each acquisition had the following parameters: TE/TR/averages was 91 ms/10,000 ms/2, field of view (FOV) was 256 × 256 mm, slice thickness/gap was 2 mm/0 mm, acquisition matrix was 192 × 192 pixels and partial Fourier encoding was 75%.

As we are working with real data, we do not have access to an exact solution of the denoising problem. Instead, we used a higher quality reference data set for comparison. This data set was obtained by registering and averaging 18 such acquisitions. The noisy input to the denoising algorithm was a data set with four averaged acquisitions consuming about 20% of the acquisition time, compared to the higher-quality one.

For better evaluation of our denoising algorithm, we have compared our 3D SA-DCT results with those obtained using the total variation PDE-based method reported in (cf. Fig. 6).¹¹ This PDE model is essentially a generalization of the well-known Rudin Osher Fatemi (ROF) model and the Blomgren Chan model.^31,32 The solution here is the minimizer u of an energy functional on the form

16

where R(u) is a regularization functional which measures the smoothness of u, and F(u,f) is a fidelity functional that measures the distance from the noisy data f and the solution u. The solution is a compromise between a completely smooth solution (λ = 0) and a solution that is close to the input data (λ≫0).

We calculated the error in the tensor in the same way as in the previous example. Using this measurement, we found the global error of the noisy image (four averages compared to 18 averages) to be 96.61. Denoising with the 3D SA-DCT algorithm was able to reduce the error to 77.07. The PDE denoising algorithm produced a solution with a global error of 76.19.