VARIATIONAL DENOISING OF DIFFUSION WEIGHTED MRI

In this paper, we present a novel variational formulation for restoring high angular resolution diffusion imaging (HARDI) data. The restoration formulation involves smoothing signal measurements over the spherical domain and across the 3D image lattice. The regularization across the lattice is achieved using a total variation (TV) norm based scheme, while the finite element method (FEM) was employed to smooth the data on the sphere at each lattice point using first and second order smoothness constraints. Examples are presented to show the performance of the HARDI data restoration scheme and its effect on fiber direction computation on synthetic data, as well as on real data sets collected from excised rat brain and spinal cord.


1.
Introduction. Observing the directional dependence of water diffusion in the nervous system can allow us to infer structural information about the surrounding tissue. Axonal membranes and myelin sheath present a barrier to molecules diffusing in directions perpendicular to the white matter fiber bundles whereas in directions parallel to the fibers, the diffusion process is less restricted [10]. This results in anisotropic diffusion that can be observed using magnetic resonance (MR) measurements by the utilization of magnetic field gradients [47]. In general, the acquired MR signal depends on the strength and the direction of these diffusion sensitizing gradients. Repeated measurements of water diffusion in tissue with varying gradient directions provide a means to quantify the level of anisotropy as well as to determine the local fiber orientation within the tissue.
In a series of publications, Basser and colleagues [6,7,8] have formulated an imaging modality called "diffusion tensor MRI (DT-MRI or DTI)" that employs a second order, positive definite, symmetric diffusion tensor to represent the local tissue structure. They have proposed several rotationally invariant scalar indices that quantify different aspects of water diffusion observed in tissue, similar to different "stains" used in histological studies [4]. Under the hypothesis that the preferred orientations of water diffusion will coincide with the fiber directions, one can determine the directionality of neuronal fiber bundles. This fact has been exploited to generate fiber-tract maps that yield information on structural connections in human [8,34,38,22] as well as rat brains [42,63,55,40,39] and spinal cords [54]. Despite its apparent success, DT-MRI has significant shortcomings when the tissue of interest has a complicated geometry. This is due to the relatively simple tensor model that assumes a unidirectional -if not isotropic-local structure. In the case of orientational heterogeneity, DT-MRI technique is likely to yield incorrect fiber directions, and artificially low anisotropy values. This is due to the Gaussian model implicit in DTI that allows only one preferred direction for water diffusion. In order to overcome these difficulties several approaches have been taken. Qspace imaging, a technique commonly used to examine porous structures [13], has been suggested as a possible solution [59]. However this scheme requires strong gradient strengths and long acquisition times [5], or significant reduction in the resolution of the images. Q-space imaging requires many images to be acquired since the space of diffusion encoding gradients is sampled on a 3D lattice. As a more viable alternative Tuch et al. have proposed to do the acquisition such that the diffusion sensitizing gradients sample the surface of a sphere [53,52]. In this high angular resolution diffusion imaging (HARDI) method, one does not have to be restricted to the tensor model and instead, it is possible to calculate diffusion coefficients along many directions. This method does not require more powerful hardware systems than those required by DT-MRI. Several groups have already performed HARDI acquisitions in clinical settings and have reported 43 to 126 different diffusion weighted images acquired in 20 to 40 minutes of total scanning time [30,52,33] indicating the feasibility of the high angular resolution scheme as a clinical diagnostic tool.
In Figure 1, we present a matrix of simulated voxels showing renderings of DTIbased estimates of orientation and HARDI-based orientation estimates computed using the scheme we present in Section 2.1. The orientation heterogeneity is evident from the HARDI-based renderings at each voxel since HARDI measurements can resolve multiple dominant directions of molecular diffusion in a voxel, a lacking feature of DTI. Since the HARDI data acquisition is very nascent, not many techniques of processing the HARDI data have been reported in literature. In the following section we will review the recently reported techniques of HARDI data denoising, which may be done prior to further analysis or visualization.
1.1. Review. We will first briefly describe the physics of acquisition and then point to various recent restoration techniques followed by methods for computing anisotropy measures from HARDI. This will be followed by an overview of our method.
1.1.1. Physics of Diffusion MR and HARDI Acquisition. The random process of diffusion of water molecules is described by the diffusion displacement PDF p t (r). This is the probability that a given molecule has a diffusion displacement of r after time t. The relation between the measured MR signal, and the diffusion displacement PDF is given by [13] (1) where S(q) is the MR signal when a diffusion gradient pulse of strength G and duration δ is applied yielding the wave vector q = γδG where γ is the gyromagnetic ratio for protons. S 0 is the image acquired with no diffusion encoding gradient applied. The above formula indicates that water displacement probabilities are simply the Fourier transform of S(q)/S 0 . It is the orientational modes of p t (r) that are taken to be the underlying fiber directions. The HARDI processing proceeds by acquiring diffusion weighted images with many diffusion encoding gradient directions, effectively sampling a spherical shell of the q-space (the space of diffusion encoding gradients) as described by Tuch [51]. It is desired that this sampling minimize the average angle between gradient directions so that the diffusion signal may be accurately reconstructed. The gradient direction for each image has been chosen to correspond to the vertices of an icosahedron which has been repeatedly subdivided. Our data sets include diffusion-weighted images acquired with the application of diffusion gradients along 81 or 46 directions in addition to one image with no diffusion weighting. Since the process of diffusion is known to have antipodal symmetry [30], we need to sample only one of the hemispheres in q-space.
1.1.2. Restoration. Processing of HARDI data sets has received increased attention lately and a few researchers have reported their results in literature. The use of spherical harmonic expansions have been quite popular in this context since the HARDI data primarily consists of scalar signal measurements on a sphere located at each lattice point on a 3D image grid. Tuch et al. [53,52] developed the HARDI acquisition and processing and later Frank [30] showed that it is possible to use the spherical harmonics expansion of the HARDI data to characterize the local geometry of the diffusivity profiles. Although elimination of odd-ordered terms and the truncation of the Laplace series provide some level of smoothing, there is no discussion of smoothing the data across the lattice points. Chen et al. [20] find a regularized spherical harmonic expansion by solving a constrained minimization problem. However the expansion is a truncated spherical harmonic expansion of order four, restricting the level of complexity that can be modeled using this approach. In [33], Jansons and Alexander described a new statistic, persistent angular structure, which was computed from the samples of a 3D function. In this case, the function described displacement of water molecules in each direction. The goal in their work was to resolve voxels containing one or more fibers. However, there was no discussion on how to restore the noisy HARDI data prior to resolution of the fiber paths. More recently, Descoteaux et al., [23,24], proposed an analytical solution to the reconstruction of the diffusion orientation distribution function (ODF). They model the signal using a spherical harmonic function of order eight and fit this model to the noisy data using a regularization constraint involving the Laplace-Beltrami operator for smoothing the HARDI data over the sphere of directions at each voxel. Their analytic form for the ODF reconstruction requires a numerical solution to a linear system and they do not consider regularization across the 3D lattice which can be important in order to obtain a piecewise smooth representation of the given HARDI data. Wiest-Daessle et al. [62,61] described several variants of non-local mean denoising applied to diffusion MRI. The approach which is applicable to HARDI involved considering the dataset as a vector-valued image, however this approach does not respect the directional relationship among the images. Assemlal et al. also employ only spatial regularization approaches to robustly determine the diffusion ODF [1] and PDF [2] fields. Savadjiev et al. [46] formulate a novel spatial regularization in terms of the underlying 3D curves which represent neuronal fibers.
In contrast to HARDI denoising, DT-MRI denoising has been more popular and numerous techniques exist in literature. For sampling of the techniques used to denoise DT-MRI, we refer the reader to [50,17,60,57,58,19,27,9,28,3,32]. Most of these works use a linearized Stejskal-Tanner equation [47] describing the MR signal decay with the exception of Wang et al., [57,58]. Using the Stejskal-Tanner equation as is, is quite important in preserving the accuracy of the restored data and this was shown in the experiments in [58]. Another important constraint in the DTI restoration is the positive definiteness of the tensors, in this context, work in [18] introduced an elegant differential geometric framework to achieve the solution. The work in [57,58] and [60] chose alternative methods to impose the positive definiteness of the restored tensor fields namely, a linear algebraic and a PDE-based method respectively. Approaches to filtering based on the Riemannian geometry of the manifold of symmetric positive-definite matrices have been reported [31,14].

1.2.
Overview of our modeling scheme. In this section we present a novel and effective variational formulation that will directly estimate a smooth signal S(θ, φ) and the probability distribution of the water molecule displacement over all directions p(θ, φ), given the noisy measurement whereŜ is the signal measurement taken on a sphere of constant gradient magnitude over all (θ, φ), b is the diffusion weighting factor, D(θ, φ) is the apparent diffusivity as a function of the direction expressed by the polar and azimuthal angles on the sphere and η(θ, φ) is Rician noise. The noise is due to additive Gaussian noise corrupting the complex-valued k-space measurements. However, for high signal-tonoise ratios we may consider η to be Gaussian distributed. A variational formulation for denoising using a data constraint based on the Rician likelihood was given by Basu et al. [9]. However, this leads to a highly nonlinear evolution equation since it involves the ratio of two Bessel functions. A modification to the non-local means algorithm which can handle Rician noise was presented by Descoteaux et al. [25]. However, neither of these approaches address smoothing over the spherical domain. In contrast, Clarke et al. [21] propose a robust method for estimating fiber orientation distributions in the presence of Rician noise, but they do not consider smoothness constraints over the voxel lattice. In this work we will assume a high SNR so that the Gaussian additive noise is a good assumption. Since we are performing high-field ex-vivo experiments, we can acquire many images and use averaging to increase the SNR so that this assumption is valid.
The variational principle involves smoothing S values over the sphere and across the 3D image lattice. The key factor that complicates this problem is that the domain of the data at each voxel in the lattice is a sphere. One may use the level-set techniques developed by Tang et al., [48] to achieve this smoothing. However, when data sets are large, it becomes computationally impractical to apply the level-set technique at each voxel independently to restore these scalar-valued measurements on the sphere. Alternative approaches to solving variational problems over nonplanar domains have been described in recent literature. Cecil et al. [15] propose several numerical approaches to dealing with discontinuous derivatives due to periodic boundaries encountered when solving problems on S 1 and S 2 . Liu et al. [37] proceed by finding a conformal mapping from the surface to the plane, then solving the problem in the 2D parameter space. Bogdanova et al. [12] presented explicit formulations of differential operators on parametric surfaces in terms of the Riemannian metric. Since our input data are sparsely distributed over a triangulated sphere (gradient directions are computed by subdividing an icosahedron, we simply use the spherical triangles as our computational domain. We arrive at a computationally efficient solution to this problem by using the finite element method (FEM) on the sphere and choosing local basis functions for the data restoration. Unlike the reported work on spherical harmonic basis expansion of the diffusivity function on the sphere [29,43,20], the FEM basis functions have local support and are more stable to perturbations due to noise in the data. From the denoised data we will compute a probability, p t (θ, φ), of molecular diffusion over a sphere of directions.
The rest of this paper is organized as follows: Section (2) contains a variational formulation of the HARDI denoising problem including smoothing the scalar signal over a sphere of directions at each 3D lattice point and across lattice points, computation of probability of water molecular diffusion over the sphere of directions and several measures of anisotropy computed from the field of probability densities. In section (3), we present several experimental results depicting the performance of our algorithms on synthetically generated and real data sets. Finally, we conclude in section (4). Appendices A and B contain the details of the finite element basis used and the element as well as the global equations.
2. Formulation of the HARDI restoration. Normally, the diffusion weighted images are quite noisy especially when acquired using large field gradients. One can reduce some amount of noise by signal averaging for each gradient direction used. However, this by itself does not preserve the details in the data. We now present a variational formulation for effective denoising of the HARDI data.
2.1. Variational smoothing. We propose a membrane-spline deformation energy minimization for smoothing the measured imageŜ(x, θ, φ). The variational principle for estimating a smooth S(x, θ, φ) is given by where Ω is the domain of the image lattice and S 2 is the sphere on which the signal measurements are specified at each voxel. The first term of Equation (3) is a data fidelity term which makes the solution to be close to the given data. The degree of data fidelity can be controlled by the input parameter µ. The second term is a regularization constraint enforcing smoothness of the data over the spherical domain at each voxel. The minimizer of this energy term is a membrane spline over the sphere which is in Sobolev space H 1 (S 2 ) [49]. The third term is another regularization term which causes the solution to be piecewise smooth over the spatial domain (the 3D voxel lattice). The minimizer of this TV norm is in the space BV (R 3 ), functions of bounded variation [35]. g(x) inhibits smoothing across discontinuities in S over the lattice. More on this in section (2.3). The choice of membrane spline smoothness over S 2 is motivated by the partial volume effect in MRI. The signal at each voxel is the average over a volume much larger than a single axonal fiber. Within this volume there may be fibers of varying orientation and regions of isotropy. Though the diffusivity function may be nearly discontinuous over S 2 at a point near a fiber bundle, it is highly unlikely for the volume average to be so. For this reason, we do not use TV norm minimization over the spherical domain.

2.2.
Finite element method based smoothing of S(θ, φ). We will consider a deformation energy functional which is a weighted combination of the thin-plate spline energy and the membrane spline energy, which is commonly used in computer vision literature for smoothing scalar-valued data in ℜ 3 (see McInerney et al., [41], Lai et al., [36]). In our case, the data at each voxel is an image on the sphere, S(θ, φ), so the problem is inherently 2 dimensional.
The diffusion-encoding gradient directions are taken as the vertices of a subdivided icosahedron, to achieve a nearly uniform sampling of gradient directions over the sphere. We map this piecewise planar approximation of the sphere to the global FEM coordinate system (u, v) by setting (u = θ, v = φ) for each gradient direction. This domain is triangulated so each face of the subdivided icosahedron will have a corresponding triangle in the (u, v) domain. A periodic boundary condition is imposed so that S(2π, v) = S(0, v). The area element in the (u, v) domain is du dv = sin φdθdφ. A similar mapping was used by McInerney & Terzopoulos [41] and Vemuri and Guo [56] for finite elements over a spherical domain.
Note that, after mapping, the data can be seen as a height field over the (u, v) plane. The smoothness of the height function, z(u, v), will be enforced by the smoothing functional The weight on the membrane term is α and the weight on the thin-plate term is β.
Once we have computed a smooth z(u, v), the result will then be mapped back to the image on the sphere, S(θ, φ). The data energy due to virtual work of the data forces, f , and virtual displacement, z, is By the principal of virtual work, the spline system is in equilibrium when the total work done by all forces is zero for all virtual displacements. The restoration of S(θ, φ) at each voxel is formulated as the energy minimization defining the equilibrium condition of the system. We use polynomial shape functions, N i , as a basis for the unknown smooth approximation, z, of the data over the u, v plane. We may write z as where N is a (1 × n) row vector, and q is a column vector of nodal variables. The domain, Ω, is partitioned into triangular elements, Ω j , each with their own local shape functions. The shape functions, in terms of local (barycentric) coordinates are given in Appendix A. For each element j, we have, In the rest of this section we will derive linear equations for the element potential energies, E j p , and data energies, E j d , in terms of the coefficients q j . Finally, we will assemble a global linear system, and solve for q. This will allow us to evaluate z(u, v) using Equation 7.
The global potential energy is the sum of the energies of each finite element, where the local potential energy function for each element is given by The element strain vector (given by Dhatt and Touzot [26]) is which may be rewritten as Since q j is constant over each element we can derive the element stiffness matrix, K, in terms of D and B giving us the element strain energy as, We will model the data constraint as springs pulling z(u, v) toward the measured values z 0 (u, v), as illustrated in Figure(2). The force at each point will obey f = k(z − z 0 ), where k is the spring constant. For small displacements the spring constant, k = µ 2 where µ is the data constraint coefficient from Equation (3). The element deformation energy due to virtual displacement z(u, v) is given by We can split the deformation energy into two terms : We may now balance the deformation energy and data energy by solving the following linear system: The global linear system for smoothing the entire mesh may be obtained by appropriately summing the local element matrices, as detailed in Appendix B. The global system is symmetric, and has a sparse banded structure with 18 nonzero diagonal bands. Since the global matrix is positive-definite, an efficient solution to q is obtained via Cholesky factorization.
2.3. Spatial smoothing of S(x). We are now ready to describe the smoothing of the data across the 3D lattice. There are many existing methods that one can apply to this problem as discussed earlier.
Smoothing the raw vector-valued data, S(x), is posed as a variational principle involving a first order smoothness constraint on the solution to the smoothing problem. Note that the data at each voxel are m measurements of S over a sphere of directions and can be assembled into a vector after the smoothing on the spherical coordinate domain has been accomplished. We propose a weighted TV-norm minimization for smoothing this vector-valued image S. This smoothing scheme reduces the effect of inter-region blurring, a drawback Gaussian convolution and isotropic diffusion suffer. Our method is a modified version of the work in Blomgren et. al., [11]. The novelty here lies in the choice of the weighting i.e., the coupling term between the channels. The variational principle for estimating a smooth S(x) is given by where, Ω is the image domain, µ is a regularization factor and m is the number of images. The first term here is the regularization constraint on the solution to have a certain degree of smoothness. The second term in the variational principle makes the solution faithful to the data in the L 2 sense. We have used the coupling function g(x) = 1/(1 + ||∇GA(x)|| 2 ) for smoothing HARDI, where GA is the generalized anisotropy index defined inÖzarslan et al., [45] and is computed from the variance of normalized diffusivity. For a more detailed discussion on GA, we refer the reader to [45]. This selection criterion preserves edges in anisotropy while smoothing the rest of the data. This anisotropy measure is chosen since it can be computed without explicitly computing the ODF, and it is our goal to smooth the data prior to ODF computation. An image of the coupling term for a typical slice is shown in Figure  (3).
Here we have used a different TV-norm than the one used by Blomgren and Chan [11]. The TV n,m norm is an L 2 norm of the vector of TV n,1 norms ( Ω ∇S i (x) 2 dx) for each channel. We use the L 1 norm instead, which is known to have better discontinuity preservation properties. The gradient descent form of the above minimization is given by The use of a modified TV-norm in equation (20) results in a looser coupling between channels than when using the TV n,m norm. This reduces the numerical complexity of Equation (21) and makes solution for large data sets feasible.
The gradient descent of the vector-valued image smoothing using the T V n,mnorm TV n,m (S(x)) = m i=1 [TV n,1 (S i )] 2 presented in [11] is given by, Note that the TV n,m norm appears in the gradient descent solution of the vectorvalued minimization problem. Considering that our data sets consist of up to 82 images, corresponding to (magnetic field) gradient directions, calculating the TV n,m norm by numerically integrating over the 3-dimensional data set at each step of an iterative process would be prohibitively expensive. In contrast using our modified TV-norm described earlier leads to a more efficient solution. We are now ready to present the numerical solution to equation (21).

Fixed-Point Lagged-Diffusivity.
Since the m Equations(21) are coupled only through the function g, we can drop the subscript on S with no ambiguity (later the subscript will refer to spatial coordinates.) In this section we will discuss the numerical solution for each channel, S, of the vector-valued image S. Equation (21) is nonlinear due to the presence of ∇S i in the denominator of the first term. We linearize Equation (21) by using the method of "lagged-diffusivity" presented by Chan and Mulet [16]. By considering ∇S to be a constant for each iteration, and using the value from the previous iteration we can instead solve Here the superscript denotes iteration number. Equation (23) can be recast in the form We now discretize the above equation in the following.

Discretized Equations.
To write Equation (24) as a linear system (AS t+1 = f t ), we discretize the Laplacian and gradient terms. Using central differences for the Laplacian we have We define the standard central differences to be We can rewrite Equation (24) in discrete form using the definitions in Equation (26) −S x−1,y,z − S x,y−1,z − S x,y,z−1 This results in a sparse linear system. The matrix of coefficients of S t+1 has 7 nonzero bands, and is given by The matrix in Equation (28) is symmetric and diagonally dominant. We employ the conjugate gradient descent to solve this system. The solution of Equation (28) represents one fixed-point iteration. This iteration is continued until |S t −S t+1 | < c, where c is a small prespecified tolerance. (1) can now be evaluated by computing the quantity S(q)/S 0 and performing the FFT. If the signal, S, is assumed to decay mono-exponentially from the origin of q-space (where S(0) = S 0 ), one can interpolate the signal values for arbitrary q. It is then possible to extrapolate (using the monoexponential decay model) from the spherical coordinate locations to grid points in cartesian space and then perform the FFT on this extrapolated data. The result is a probability of water molecule displacement over a small time constant. Since the quantity of interest is primarily the direction of water displacement, one can integrate out the radial component of p t (r) to get p t (θ, φ). This is commonly referred to as the diffusion orientation distribution function or simply ODF. Computing the ODF with this method is computationally expensive since it requires a 3D FFT at each voxel, and then a numerical integration for each direction. For the sake of efficiency, we will compute a probability profile (not the ODF), which will make processing large datasets feasible. This probability profile, written as p t (r, θ, φ) quantifies the probability that a water molecule diffuses through a sphere of fixed radius, r. A more detailed treatment byÖzarslan et al.

Computing probabilities. The probability in Equation
can be found in [44]. This scheme provides a fast way to calculate the orientation profiles. In our implementation we have evaluated the series given in [44] up to l = 6 terms since the reconstructed surfaces have very simple shapes which can be accurately represented using a truncated spherical harmonics series, and r 0 was set to 17.5µm. An alternative approach is to use the Funk-Radon Transform proposed by Tuch [51], however this introduces smoothing due to a spherical convolution step which would make evaluation of our denoising algorithm more difficult.
To enhance the visual impact of the probability profiles we apply a sharpening transform to the distribution by subtracting a uniform distribution (sphere) from each profile, as shown in Figure (4). The radius of the sphere is the minimum of the probability over all directions. By performing this operation the direction of maximum probability becomes more apparent. 3. Experimental results. The denoising and rendering techniques described in the previous section were first applied to a synthetic HARDI dataset. This dataset was generated using the technique described byÖzarslan et al. in [45]. The dataset was designed to depict a region of curving fibers, a region of straight fibers, and a crossing between the two. A total of 81 acquisition directions are simulated with b = 1500 s/mm 2 .
A small sample of the probability surfaces p(θ, φ) computed from the synthetic data, taken from near the crossing region, is shown in Figure (5a). The real-valued synthetic data was corrupted with Gaussian noise of zero mean, and variance σ 2 = 0.005. p(θ, φ) surfaces computed from the noisy data (without any denoising) are shown in Figure (5b). The same voxels are shown -after smoothing over the spherical manifold at each voxel independently -in Figure ( figure 5e) depict better smoothing than those in either of 5c) or 5d), visually indicating that one needs to perform smoothing on the sphere and across the lattice and not just one or the other.
From Figure (5b), it can be seen that the noise has a large influence on the smoothness of the distribution. As expected from the variational formulation, the spikes of noise present in the raw data have been smoothed while preserving the overall shape of the S profile. This smoothness is evident in the computed probability profiles as well. Table 1. Error between ground-truth probabilitiesp and probabilities computed from restored synthetic data when SNR = 14.  Table 2. Error between ground-truth probabilitiesp and probabilities computed from restored synthetic data when SNR = 5. A quantitative evaluation can be obtained by comparing the distributions computed from the smoothed data with the ground-truth by using the square root of J-divergence (symmetrized KL-divergence) as a measure. This divergence is defined as (29) d(p, q) = J(p, q) where In Table (1) we compare the distances, d(p, p), between the densities computed from the original synthetic data, (p), and the unrestored data, the data restored only using the FEM method, the data restored using only the TV-norm minimization, and the data restored using both techniques. For each technique, the mean distance, µ(d(p, p)), between the densities in corresponding voxels and the standard deviation, σ(d(p, p)), is presented. As evident from Table (1), the TV restoration has superior performance over the FEM technique in terms of the mean error. The combination of techniques has a lower mean error and standard deviation of the error than either the L 2 -norm based or the TV-norm based restoration. Note also that the error achieved by applying smoothing over the sphere prior to smoothing over the voxel lattice is lower than when the order is reversed. Since TV-norm minimization can be seen as a nonlinear diffusion process, performing the denoising in this order propagates smoothed intensities within homogeneous regions. In subsequent experiments we perform the denoising in this order. The denoising algorithm was applied to a dataset consisting of one non-diffusion weighted image and 46 diffusion weighted images of a rat spinal cord. Our data were acquired using a 14.1 Tesla (600 MHz) Bruker Avance Imaging spectrometer system with a diffusion weighted spin echo pulse sequence. Imaging parameters were : TR = 1400 ms, TE = 25 ms, Delta = 17.5 ms, delta = 1.5 ms, b high = 1500s/mm 2 , b low = 0s/mm 2 , diffusion gradient strengths = 0 mT/m with 28 averages were measured for b = 0s/mm 2 and diffusion gradient strengths = 733s/mm 2 with 7 averages were measured for each of the 46 diffusion weighting-gradient directions. The 46 directions were derived from the tessellation of a hemisphere. The image field of view was 4.3×4.3×12mm 3 , acquisition matrix was 72×72×40. The approximate SNR for the S 0 and diffusion weighted images were 58 and 50 respectively. The parameter values used for the restoration were µ = 0.97, α = 0.02, β = 0.0, k = 1.0.
Axial slices before and after denoising are shown for the non-diffusion weighted image in Figure (6) and one diffusion weighted image in Figure (7). The ringing artifacts visible near the sample boundary in Figure (6) have been noticeably decreased. Note that the edges in the image have been well preserved.
Figures (8) and (9) show restored probability profiles from rat brain and spinal cord datasets. The brain data were acquired using a 17.6 Tesla (750 MHz) Bruker Avance Imaging spectrometer system with a diffusion weighted spin echo pulse sequence. Imaging parameters were : TR = 2000 ms, TE = 28 ms, Delta = 17.8 ms, delta = 2.2 ms, b high = 1500s/mm 2 , b low = 0s/mm 2 , 6 averages for each of the 81 diffusion weighting-gradient directions. The 81 directions were derived from the tessellation of a hemisphere. The image field of view was 150 × 150 × 300µm 3 , acquisition matrix was 100×100×60. The approximate SNR for the S 0 and diffusion weighted images were 206 and 177 respectively. Figure (8b) shows a detail from the rat hippocampus. The piecewise smoothing behavior of the algorithm is evident within the anisotropic hippocampus region. This region has been smoothed independently of the more isotropic surrounding regions. The spherical smoothing term has also suppressed some peaks of the distribution which were probably due to noise in the acquired data. Figure (8c) shows a detail from the rat corpus callosum. The data dependent coupling term in the restoration algorithm has permitted intraregion smoothing within the corpus callosum while preventing interregion smoothing. Note that the fiber directions within the corpus callosum have been well preserved.
Figure (9) shows details from the rat spinal cord dataset. Here the noise reduction can be seen to enhance the coherence of structures in the inner core of grey matter.
The data were processed by a MATLAB implementation of the algorithm running on a system with Intel Quad Core QX6700 2.66 GHz CPU and 4 GB RAM. The computation times for the finite element smoothing over the sphere depends on the number of diffusion-encoding gradient directions in the image acquisition. For the spinal cord data with 46 directions the time was 0.018 seconds per voxel, and for the brain dataset with 81 directions the time was 0.038 seconds per voxel. The computation time for the TV-norm minimization problem for each diffusion weighted image depends on the size of the acquisition matrix. For the spinal cord the resolution was 72 × 72 × 40 and the computation time per image was 28.3 seconds. For the brain dataset the resolution was 100 × 100 × 60 and computation required 82.4 seconds per image. 4. Conclusion. In this paper, we presented a new variational formulation for restoring HARDI data and an FEM technique for implementing the restoration. Our formulation of the HARDI restoration involves two types of smoothness constraints. The first is smoothness over the spherical domain of acquisition directions, and the second is smoothness between neighboring voxels in the Cartesian domain. The smoothing technique is capable of preserving discontinuities in the data. This was demonstrated on synthetic and real anatomical data. By using J-divergence as a measure of distance between distribution, we were able to show quantitatively that the combination of restoration techniques performs better than either technique alone.
Appendix A.
Local element coordinates. We now present the coordinate system for the local elements. For local elements, triangular elements are used with a barycentric coordinate system (γ, ξ, η). Each coordinate is in the range [0, 1] and γ + ξ + η = 1 for points on the triangle.
The global coordinates, (u, v), can be computed from the local coordinates by The Jacobian, J, of the transformation between coordinate systems is defined by Integrals over the (u, v) domain to be converted to integrals over the local (ξ, η) domain by Using the Gauss-Radau quadrature rules given in [26], we can approximate the integral in Equation (33) by the summation where η i,j = r i (1 − ξ j ), wj j = a j (1 − ξ j ), ξ j , and wi i are given in Table 3.
Derivatives over (u, v) can be written in terms of local coordinates by applying the chain rule: The partial derivatives of ξ and η with respect to u and v can be computed by inverting the Jacobian The inverse of J is given by We use the fifth order element shape functions given by Dhatt and Touzot [26]. This element guarantees C 1 (surface normal) continuity across triangles. The quintic basis functions are given by The quintic shape functions have nodal variables which can be written in terms of local or global coordinates as, The local and global nodal variables are related to each other by Appendix B.
Global matrices. We wish to construct global matrices so that the energy balance over the entire FEM mesh is given by the linear system where K is a (6n × 6n) matrix since we have 6 variables per node. We will consider the simple case of 2 elements. Expanding the element Equation (19) in terms of nodal variables for element 0 we get where each q j l is a (6 × 1) column vector of nodal variables. We expand each K j to be (6n × 6n) by inserting rows and columns of zeros corresponding to each node of the mesh. Also expand f j to (6n × 1). The global K and q are obtained by summing the expanded matrices from each element in the mesh. For our 2 element example we have