Les cahiers du GREYC

Les cahiers du GREYC Année 2007 numéro 1 Haz-Edine Assemlal, David Tschumperlé, Luc Brun {Haz-Edine.Assemlal, David.Tschumperle, Luc.Brun}@greyc.ensicaen.fr A Variational Framework for the Robust Estimation of ODFs From High Angular Resolution Diffusion Images Groupe de Recherche en Informatique, Image, Instrumentation de Caen. CNRS - UMR 6072 Université de Caen - Campus II Bd du Maréchal Juin, 14032 Caen Cedex - FRANCE Tél : 02 31 56 73 31 - Fax : 02 31 56 73 30 E-mail : greyc@info.unicaen.fr Ecole Nationale Supérieure d’Ingénieur de Caen Bd du Maréchal Juin, 14050 Caen Cedex - FRANCE Tél : 02 31 45 25 04 - Fax : 02 31 45 26 98 E-mail : greyc@ensicaen.fr Abstract We address the problem of estimating complex diffusion models from high angular resolution diffusion MRI images (also known as HARDI datasets). Rather than considering a classical 2nd-order tensor to model the water molecule diffusion in tissues, we describe each voxel diffusion by a model-free orientation distribution function (ODF) expressed as a set of spherical harmonics coefficients. We propose to estimate the ODFs volume directly from the raw HARDI data by minimizing a nonlinear energy functional which considers the non-gaussianity of the MRI Rician noise as well as introduces a regularity constraint on the estimated field. The estimation is thus performed by a set of multi-valued partial differential equations composed of both robust estimation and discontinuity-preserving regularization terms. We show that fiber-tracking is more accurate when using this regularized estimation as opposed to non-regularized methods. We finally illustrate the importance of these constraints in the ODFs estimation process through both synthetic and real HARDI datasets. 1 Introduction Diffusion Magnetic Resonance Imaging (dMRI) [16] is a non-invasive method which allows to observe the Brownian motion of water molecules in the brain tissues in vivo. As this motion is constrained by nearby fibers in brain white matter, dMRI enables to collect local informations on the wiring structures inside the human brain. In this setting, Diffusion Tensor Imaging (DTI) is a well-known particular case of such a modality which maps each voxel signal to a 2nd-order tensor model [3, 19, 29]. Anyway, doing this implicitly assumes that the diffusion is locally Gaussian, which leads to limitations when estimating intra-voxel structures different from single fiber configurations : a Gaussian function has only a single directional maximum and therefore cannot retrieve possible multiple maxima of the diffusion function, including patterns like crossing or kissing fibers. As the brain white matter is known to have such complex arrangements, the inadequacy of the DTI classical tensor model leads to imprecise estimation of the underlying fiber map structure. In order to overtake this significant shortcoming, numbers of higher order models have been proposed in the literature. Actually, these models better consider the physics behind the MRI acquisition : the measured signal originates from protons of hydrogen nuclei, which are mostly found in water molecules. When presented to a specific magnetic field, the rotating magnetisation vectors of the spins induce an electromotive force which constitutes the MR signal. Let P (x, x0 ) be the conditional diffusion probability density function (PDF) which describes the probability for a spin to displace from position x0 to position x in the experimental diffusion time τ [5]. The observed signal is the average over all spins within a voxel. So, the resulting ensemble average propagator of a relative motion r can be expressed as Z P (r) = P (x − x0 ) = P (x, x0 )ρ(x0 )dx0 2 where ρ(x0 ) is the initial spin density [5]. Stejskal and Tanner [23] proposed a spin echo sequence which gives a relation between the signal attenuation S(q) and the diffusion probability density function P (r), Z S(q) T = P (r)e−2πiq rq dr = F[P (r)] (1) S0 where F denotes the 3D Fourier transform, q the diffusion wave-vector and S0 the baseline image taken without any gradient. The latter is defined as q = γδg/(2π) where γ is the gyromagnetic ratio for the proton nucleus, δ is the diffusion gradient duration, and g the diffusion gradient vector. S0 stands for the baseline image, which is a MR image acquired without any preferred diffusion gradient. Equation (1) naturally suggests to use the inverse Fourier transform to estimate the PDF. This technique proposed by Tuch in [25] is known as Q-Space Imaging (QSI), but has significant restrictions as it cannot retrieve the PDF from in vivo acquisition [5, 25]. First, phase diffusion signal is often corrupted by biological motion mostly due to cardiac pulsations. Tuch [25] proposed to solve this problem by the use of the modulus Fourier transform P (r) = F [|S(q)|] instead, since the diffusion signal is real and positive. Nevertheless, a sufficient sampling of the PDF needs numerous acquisitions in the diffusion signal space (known as Q-Space) and therefore requires a huge acquisition time meanwhile the patient should remain motionless. Finally, QSI requires higher q values than standard scanners values and consequently requires higher gradient g which creates eddy current distortions in the magnetic field. As a result of QSI limitations, High Angular Resolution Diffusion Imaging (HARDI) [25] comes as an interesting alternative. Instead of sampling the diffusion signal all over the space, the acquisition is made on the single sphere following ns discrete gradient directions. The radial component of the three-dimensional PDF is lost, but it provided only details on tissue microstructures, and did not give valuable informations about the diffusion orientation. Hence, diffusion orientation can be measured through the Orientation Distribution Function (ODF) defined as the radial projection of the spherical diffusion function. Other higher order models based also on the Gaussian assumption of the diffusion have been proposed in the literature: multi-fiber Gaussian tensors which model the signal as a finite number of Gaussian fibers [25] and spherical deconvolution techniques which use Gaussian kernels to estimate the diffusion signal [17, 24]. These approaches are model-based methods, implying a strong a priori knowledge about the local fiber configuration. On the contrary, Q-Ball Imaging (QBI) proposed by Tuch in [26] has the advantage of being model-independent. Q-Ball Imaging (QBI) seeks to reconstruct the ODF from HARDI data. Given a unit spatial direction u ∈ R3 , Ψ(u) is the radial projection of the PDF on the line directed by u. Thus, the exact ODF Ψ can be written without loss of generality with u taken as the z-axis, as Z Z ∞ P (αu)dα = P (r, θ, z)δ(θ, z)rdrdθdz (2) Ψ(u) = o Tuch [26] showed that the Funk-Radon transform (FRT) G from the raw HARDI data approxi3 mates the ODF on the Q-space single sphere: Z ′ Gq′ [S(q)](u) = 2πq P (r, θ, z)J0 (2πq ′ r)rdrdθdz (3) where J0 stands for the zeroth-order Bessel function. Consequently, the estimated ODF in a direction u is given by the integral over the diffusion signal on the plane orthogonal to u. This leads to an interesting model-free method for retrieving orientation diffusion informations. The FRT (3) is very close to the exact ODF (2), except for the Dirac function δ replaced by J0 . Then instead of projecting the diffusion function along a thin line, the projection is done along a Bessel beam and the larger q ′ is, the closer the J0 approximates δ. As a consequence, the FRT from the raw HARDI data is smoother than the exact ODF. Tuch proposed in [26] an algorithm to reconstruct the FRT, but it implies complex numerical computations such as diffusion signal interpolation on the sphere using a kernel fit. Yet, Descoteaux et al. proposed in [10] a very elegant analytical resolution of the ODF leading to a linear estimation technique. Once having estimated diffusion directions, an interesting application of diffusion MRI consists in retrieving neuronal fibers in brain white matter by the mean of a so called fiber-tracking algorithm. This is classically done by computing the integral curve of interpolated DTI dominant eigenvectors [7, 27]. However, these methods are very sensitive to noise since it always suppose that the dominant eigenvector is correct. Noise issue was tackled in [8, 9, 27] who proposed to apply regularization schemes on tensor or principal direction before applying the fiber-tracking step. One of the main limitation of the DTI model is that it is not able to retrieve several intravoxel fiber distributions, leading to wrong or biased estimation of dominant fiber directions. On the other hand, recent higher order models as ODF fields are promising for estimating correct neuronal fibers. In the following sections, we remind the linear estimation technique of the ODFs introduced by Descoteaux et al. [10] (section 2). In section 3, we present our contribution, i.e. a new variational framework for a more robust estimation of ODFs. It has the advantage of being nonlinear, allowing to estimate and regularize simultaneously a whole volume of ODFs. Then, we detailed our fiber-tracking algorithm in section 4. And we finally validate this model by illustrating results on synthetic and human brain HARDI data (section 5). 2 Linear Estimation Geometrically speaking, the FRT value at a given direction u is the great circle integral of the signal on the plane orthogonal to u. Complex numerical schemes [26] have been proposed to compute this integral and Descoteaux et al. [10] recently proposed an elegant analytical method based on Funk-Hecke theorem [2] to calculate this integral from a signal expressed in the spherical harmonics (SH) basis. The spherical harmonics Ylm of degree l > 0 and phase order m ∈ [−l, l] are a set of orthonormal functions and form a basis to describe complex functions defined on the unit sphere. 4 Imaginary and non-symmetric parts of the diffusion signal essentially capture noise, so it is interesting to force the fitted spherical function to be constraint to be symmetric and real. So we use a modified basis constrained to be describe orthonormal, symmetric and real spherical functions only [1, 10, 11]. Let Ylm be a spherical harmonic of degree l and order m in the standard basis, Yj (j = (l2 + l + 2)/2 + m denotes the order) be a spherical harmonic in the modified basis. √  √ 2 m if − l ≤ m < 0 ((−1)m Ylm + Yl−m ), 2 Re(Y ) =  l 2     Yl 0 , if m = 0 Yj = (4)     √  √ 2 Im(Ylm ) = 22i ((−1)m+1 Ylm + Yl−m ), if 0 < m ≤ l where Ylm is defined using associated Legendre polynomials Plm . Thus, any function χ defined on the unit sphere ∀θ, φ ∈ Ωχ = [0, π] × [0, 2π[ , χ : Ωχ → R can be described as: N X χ(θ, φ) = cj Yj (θ, φ) (5) j=1 where N corresponds to the highest degree of the decomposition into modified SH basis Yj . Let S : R3 → ΩS ∈ Rns be the vector field of diffusion signal in ns discrete directions on the sphere (typically S(p) (q) = S(q)/S0 ) and C : ΩC ∈ R3 → RN be the vector of coefficients of spherical harmonics at voxel p = (x, y, z). Descoteaux et al. [10] proposed to fit the signal with a continuous spherical function by a least square minimization minχ∈Ωχ ||S(p) (θi , φi ) − χ(θi , φi )||2 = (6) 2 minC∈ΩS ||S(p) (θi , φi ) − B̃C(p) (θi , φi )|| where θi , φi follow gradient discretization of the diffusion signal on the single sphere , and B̃ is a matrix of SH functions Yj . Best fitting coefficients C are then given by a modified Moore-Penrose pseudo-inverse scheme. C(p) = (B̃ T B̃ + λL̃)−1 B̃ T S(p) (7) where λ is the weight term on the regularization matrix L, which penalizes high degrees in the estimation under the assumption that high degrees SH are likely to capture noise. This is somehow similar to the low-pass filter introduced by Tournier et al. in [24]. At this point, we have a continuous spherical function fitting the diffusion signal. We want now to recover the ODF which gives the orientation of the diffusion. Descoteaux et al. [10] showed that the FRT approximating the ODF can be expressed using the SH basis, by: X Gq′ [S(p) (q)] = P̃ B̃C(p) = 2πPlj (0)cj(p) Yj(p) (8) j 5 where P̃ is a transition matrix from Q-space signal to diffusion probability space, Plj are associated Legendre polynomials at order lj (value of l knowing j). Both P̃ and B̃ are independent to the estimated voxel p. Descoteaux et al. actually pointed that the spherical harmonics are a powerful tool to recover an approximation of the ODF. However, it is important to notice that MRI noise distribution follows a Rice distribution not a Gaussian one. This arises from the Gaussian noise in the original complex frequency domain of the diffusion signal [12, 14, 15, 22]. Therefore, a least square fit is definitely not the best choice for such an estimation process. Furthermore, estimation is made voxel-by-voxel and does not reflect the spatial regularity of the diffusion function. Our contribution is to propose a variational framework which is adaptable to noise distribution and is able to use valuable informations given by the neighbour voxels. 3 Variational Framework 3.1 Theory The key idea is to estimate and regularize the whole volume of voxels at the same time. Indeed, it enables to take into account correlation between all parts of the processing pipeline instead of doing the different parts separately. From the decomposition of the observed dMRI volume S into a filtered dMRI volume Sr and noise ǫ, we seek to reconstruct Sr = B̃ P̃ C where C : ΩC ∈ R3 → RN stands for the volume of spherical harmonics coefficients. N ∈ R is the highest SH degree and P̃ the N -rank order matrix defined by Descoteaux et al. in [10]   .. .   2πPlj (0) (9) P̃ =   .. . where Plm are the Legendre polynomial of degree lj and phase order m, so that j = (l2 + l + j 2)/2 + m is the index in the modified SH basis [10]. Let ns ∈ R be the number of gradient directions and B̃ be the matrix of size (ns , N )   Y1 (θ1 , φ1 ) . . . YN (θ1 , φ1 )   .. .. .. B̃ =  (10)  . . . Y1 (θns , φns ) . . . YN (θns , φns ) with Yj the modified SH basis functions. We propose to robustly estimate and regularize the ODF field simultaneously by minimizing this nonlinear functional energy E defined as: # ( ) Z "X ns min E(C) = ψ(Srk ) + αϕ(||∇C||)dΩS (11) C∈ΩC ΩS k 6 where ψ(Srk ) stands for the likelihood term which measures the dissimilitudes at voxel p between the raw signal S and its ODF estimation Sr = B̃ P̃ −1 C at gradient direction k, ψ : R → R+ and ϕ : R → R+ are real and positive functions, α ∈ R is the regularization weight and ||∇C|| the gradient norm defined as X ||∇C|| = ||∇Cj || (12) j Note that if ψ(s) = s2 and α = 0 in (11), we minimize the least square criterion (7, corresponding to the Descoteaux’s method with λ = 0). As the minimization cannot be computed straightforwardly, the gradient descent coming from the Euler-Lagrange derivation of (11) leads to a set of multi-valued partial derivate equation (PDE) (13). In practice, we first set C(t=0) to U0 , an initial estimate of SH coefficients. In order to estimate a solution, SH coefficients velocity ∂C giving the direction from the current C to a solution is computed. The latter is done several ∂t < ε,). times until convergence (typically when ε ∈ R+ , ε → 0, ∂C ∂t   Ct=0 = U0  ∂Cj ∂t = Pj−1 (13) Pns k B̃k,j ψ ′ (Srk ) + α div(ϕ(||∇C||)) The initial estimate U0 is computed either by considering a random field or a more structured one. A good choice is to start from an initial set which is not so far from the global minimum; so the linear least square estimation (7) seems to be an adequate alternative. Indeed, least square minimization is the global minimum when ψ(s) = s2 and α = 0. One can expect the minimum to be close enough to the least square minimum through variations of ψ and ϕ; and should consequently bring down the number of iterations required to converge. It is worth to mention that similar methods have been proposed for the regularization of DTI [9, 18, 21] and apparent diffusion coefficient (ADC) [6, 28]. Yet none is able to take advantage of the informations provided by the ODFs, and we propose a new approach for ODFs estimation. 3.2 Likelihood function ψ The diffusion MR magnitude images are corrupted by noise ǫ. Yet, as MRI noise follows a Rician distribution, least square criterion is not the best choice. The ψ function should rather be defined to support a robust ODF estimation, increasing with the dissimilitudes between the raw signal S and its reconstruction Sr = B̃ P̃ −1 C and having a threshold r ∈ R+ to stage E for irrelevant data (cf Fig.1). Nonetheless, the best ψ function is the one specific to the MR scanners, that is to say the Rice distribution, which probability density function is     S S · Sr −(S2 + S2r ) p(S|Sr , σ) = 2 exp I0 (14) σ 2σ 2 σ2 7 Figure 1: Example of likelihood ψ functions for robust estimation S = 10 and σ = 5: dot-dashed is a generic function and solid line is the specific function to Rice noise. where σ is the standard deviation of the noise and I0 is the modified zeroth-order Bessel function of the first kind. We adapt the Rician bias correction filter introduced by Basu et al. in [4] from 2nd-order DTI to complex ODF. It is based on a maximum a posteriori approach so we construct the filtered volume Sr that maximizes the log-posterior probability log p(Sr |S) = log p(S|Sr ) + log p(Sr ) − log p(S) (15) where p(S|Sr ) is the likelihood term, p(Sr ) is the prior, or the regularization term, and p(S) is the normalizing constant. We are interested in the likelihood term, thus combining equation (14) and equation (15) the pointwise loglikelihood becomes   S (S2 + S2r ) S · Sr log p(S|Sr , σ) = log 2 − = ψ(Sr ) (16) + log I0 σ 2σ 2 σ2 Fig.1 illustrates variation of the opposite function with scalar values of Sr when S = 10 and σ = 5. The energy is low when S ≈ Sr and increases with their dissimilitudes. Note that σ has to be known a priori and can be either retrieved as a parameter specific to the MR scanner, or can be either computed from a uniform area as described by Sijbers et al. in [22]. Combining equation (13) and the derivative of equation (16) with respect to Cj gives the PDE adapted to Rician noise, P       Pj−1 nk s B̃k,j ∂Cj S · Sr S · Sr −Sr + S I1 /I0 + α div(ϕ(||∇C||)) = ∂t σ2 σ2 σ2 8 (17) 3.3 Regularization function ϕ Besides, regularization should go strong on homogeneous area (low ||∇C||), and preserve contours not only between isotropic and anisotropic regions but also among voxels with different number of fibers (large ||∇C||) (cf Fig.2). Figure 2: Example of possible regularization function for ϕ. Frequential regularization, as described by Descoteaux et al. in [10] was introduced as a frequential low-pass filter L̃ in (7). We drop this frequential regularization step since the spatial regularization puts sufficient constraints on the diffusion signal to be estimated. Our experiments have confirmed that using both spatial and frequential regularization is useless. 4 Fiber-tracking DTI-based fiber-tracking has been widely used [7, 8, 9, 27] but it has significant drawbacks when dealing with intra-voxel structures. Indeed, not only DTI cannot model crossing or kissing fibers but it also estimates wrong directions in the case of multiple fiber configurations. On the contrary, ODF does not fall into this restrictions. Nevertheless, although the issue of robust fiber-tracking has received numerous contribution with DTI model it is still an open problem when using ODFs. In order to illustrate the influence of a robust ODF estimation on fiber-tracking, we propose a model for retrieving neuronal fiber in brain white matter. A way to do fiber-tracking is to use estimated displacement due to diffusion which is given by the ODF in order to find dominant directions. Once directions are retrieved, only one is kept based on a a priori on the fibers distribution, resulting in a diffusion tensors field w. A line 9 integration scheme is needed to propagate a fiber along a curve C through the tensors volume (c.f . Fig.7). One may want to use Euler method Ca+h = Ca + hwa + O(h2 ) (18) where a is the current position in the curve C and h is the integration step. In practice, Euler’s method is not stable and precise and so Runge-Kutta [20] comes as an interesting alternative. This method can be seen as the result of reduction in precision of a curve C ′ more precise than C because of a smaller integration step Ca+h = Ca + k1 k2 k3 k4 + + + + O(h5 ) 6 3 3 6 (19) where ki are the slope estimated in a + i/4h. Actually, fourth-order Runge-Kutta is by far the most precise and is the one we used on our tests. Besides, we assume that there are no neuronal fibers in water regions of the brain, and consequently there is a need to identify this regions. Generalized Fractional Anisotropy (GFA) (c.f . bottom of Fig.7.a) as proposed by Tuch in [26] measures the variation within the diffusion as a spherical function. It can be expressed in the spherical harmonics basis which has the advantage to be much faster to compute. s c2 std(Ψ) GF A = = 1 − PN 0 2 (20) rms(Ψ) j=0 cj This gives a convenient way to measure apart isotropic from anisotropic area; therefore we used it to stop fiber line integration when arriving in water area, i.e. when GF A is below a threshold. 5 Numerical Experiments For all our experiments, we used the robust Rician estimation function ψ presented in section 3.2 and the discontinuity-preserving regularization function ϕ(s) = 1+s12 /κ2 , where κ1 and κ2 are two thresholds depending on the value range of the original HARDI dataset. We first present results of our variational framework on synthetic HARDI data. These data are created using a Gaussian multi tensor model [1] to simulate n fibers crossing: S(u) = n X TD u k pk e−bu (21) k=1 where u follows a discretization of the sphere (72 directions obtained from the subdivision of a regular icosahedron), pk corresponds to the proportion of the kth fiber in the voxel and Dk , a 3x3 10 matrix, stands for the eigen values of each individual fibers. Our synthetic data simulate horizontal and vertical fibers (respectively bottom left and upper right in Fig.3(a)) crossing (bottom left in Fig.3(a)) with surrounding water regions (upper left region in Fig.3(a)). For visualization purposes in Fig.5, ODF are slightly sharpened using Laplace-Beltrami operator ∆b as described by Descoteaux et al. in [10], resulting in C(p)sharp = C(p) − β∆b C(p) , where β is the sharpening intensity. As a side effect this sharpening tends to transform isotropic area into anisotropic ones. This is visible in the left upper part. It is worth noting that chosen gradient norm ||∇C|| is an adequate measure to set apart isotropic area from anisotropic area as the corresponding gradient has a strong value on signal discontinuities (cf Fig.3(b)). Subsequently, divergence div(ϕ||∇C||) performs well in regularizing homogeneous area without degrading the contours as shown in Fig.3(c). (a) Synthetic data (b) ||∇C|| (c) div(ϕ||∇C||) Figure 3: Gradient norm and its corresponding diffusion divergence for the synthetic example. The brighter the pixel is, the higher ||∇C|| is and the stronger the regularization is. In order to simulate dMRI acquisitions, we add Rician noise [12, 14, 15, 22] of variance σ to the signal. Let ζ be the PSNR between the noisy and original data. Generalized Fractional Anisotropy (GFA) (c.f . Fig.4) as proposed in equation 20 measures the variation within the diffusion as a spherical function. This comparison demonstrates the need for a regularization within the estimation process. Indeed, GFA is adequate to have hindsight on global coherence of the ODF estimation of the volume since every voxel is summarized by a scalar value. When it comes to noisy input data, regularization greatly improves the spatial coherence of the volume estimation as illustrated in Fig.4. Besides, this gives a convenient way to measure apart isotropic from anisotropic area; therefore we used it to stop fiber line integration when arriving in water area, i.e. when GF A is below a threshold. We computed statistics on performance of LS and PDE estimation submitted to input HARDI raw data of PSNR ζ ∈ [5, 30], and to several regularization weight α ∈ [0, 1]. For each ζ and α, nd = 5 trials have been made. Finally, we compute PSNR (P SN RLS and P SN RP DE ) between the exact data and its estimation. Fig.6 show influence of ζ and α on P SN RP DE − 11 (a) Anisotropy of the (b) Anisotropy of the (c) Anisotropy of the original data. Least Square estimation PDE estimation on on noisy data. noisy data. Figure 4: Comparison of generalized anisotropy. Noisy data with ζ = 20. (a) LS (b) PDE (c) Zoom. LS/PDE Figure 5: Crossing fibers surrounded by water. Comparison between Least Square (LS) and our variational framework estimations on noisy data, ζ = 20. Color indicates the diffusion direction. For visualisation purposes, ODF are sharpened. (a) Linear LS estimation on noisy data. (b) PDE estimation on noisy data. (c) Zoom on differences. Left is LS estimation, right is PDE estimation. P SN RLS and highlights the interest of spatial regularization, especially when having low PSNR input data. Out of the results, the larger the diffusion regularization weight α is, the better the results of PDE estimation over the linear LS estimation are. Besides, the lower ζ of input data is, the larger the difference between linear LS estimation and our variational approach using regularization. Though PSNR could not be considered as an optimal measure of the quality of estimation since noise describes a Rician distribution, subjective quality from estimated data (cf Fig.5) endorse it. Indeed, Fig.5(c) shows that even when a voxel is corrupted, the reconstruction using neighbourhood helps enhancing reconstruction. Fiber-tracking was tested on another synthetic data-set which simulate horizontal and vertical fibers (respectively right and top in Fig.7) merging into one horizontal fiber (left in Fig.7). From 12 PDE - LS PSNR(denoised, original) 14 12 10 8 6 4 2 00 0.8 0.6 5 10 15 20 25 PSNR(noised, original) 0.4 0.2 30 0 Regularization weight Figure 6: Influence of Diffusion weight and PSNR of input data on the difference of PSNR between LS estimation and original data and PSNR between PDE estimation and original data. the several fibers distributions estimated, we retrieved one using a simple a priori, i.e. to follow the direction which is the most vertical. As expected, DTI is not able to retrieve correctly the profile of any underlying fiber as shown in Fig.(7.b). Instead, it estimates a wrong direction, which is a mixture of the two main directions from each fiber distribution. Therefore the estimated fiber is a fictive one since a correct path in this dataset would be either horizontal or going vertical. Fibertracking on ODF does not have this problem, but it is sensitive to noise (c.f . Fig.7.c). However our variational method successfully estimates the ODFs field from noisy data (P SN R = 15dB), which leads to good fiber-tracking (c.f . Fig.7.d). Finally, we test our variational framework on a human brain HARDI dataset. The raw HARDI data were acquired using a 1.5T MRI scanner, with b = 500s/mm2 and 31 gradient directions resulting on a volume of size 256x256x16. Fig.8 shows primary results, verified against Gray’s atlas [13], and compared to DTI and LS estimation. The region shown is known to be the cross of several fibers. Though we used the same sharpening coefficient β using Laplace-Beltrami operator to LS (cf Fig.8(d)) and PDE (cf Fig.8(e)) vizualisation, PDE estimation is sharper in the presence of fibers crossing (cf Fig.8(f)). This originates in the three-dimensional spatial regularization, i.e. PDE estimation considers neighbour voxels in the same slice and from other slices which do not appear in Fig.8). 13 (a) GFA (b) DTI (c) LS (d) noisy PDE Figure 7: Crossing fibers distributions: estimation and fiber-tracking. 6 Discussion and Conclusion We have proposed a new and robust variational reconstruction framework to estimate ODFs from HARDI data in a robust way, and showed that this approach has numbers of advantages. These include model-independence and more robustness to MRI Rician noise distribution. Furthermore, the ability to reconstruct a voxel taking the whole neighbourhood informations into account clearly enhance the spatial coherence of the reconstruction. This last property greatly improve the ability to recover reliable and accurate intra-voxel fibers distributions within the human brain and opens promising perspectives for studying more precisely the neuronal fiber network. References [1] D. Alexander, G. Barker, and S. Arridge. Detection and modeling of non-gaussian apparent diffusion coefficient profiles in human brain data. Magn. Res. in Med., 48(2):331–340, 2002. 5, 10 14 (a) Anisotropy (d) LS (b) Brain atlas (e) PDE (c) DTI (f) Zoom. LS/PDE Figure 8: Estimation of human brain white matter fibers, in a region with numbers of fibers cross: bottom left of image is corpus callosum, upper right is cingulum gyrus and upper right is frontal gyrus. (a) Colored Fractional anisotropy from DTI. (b) Brain Atlas [13]. (c) DTI estimation. (d) ODF estimation using LS. (e) ODF estimation using PDE. (f) Zoom on differences. Left is LS estimation, right is PDE estimation. [2] G. Andrews, R. Askey, and R. Roy. Special Functions. Cambridge University Press, 1999. 4 [3] P. Basser, J. Mattiello, and D. LeBihan. Mr diffusion tensor spectroscopy and imaging. Biophysical Journal, 66(1):256–267, 1994. 2 [4] S. Basu, T. Fletcher, and R. Whitaker. Rician noise removal in diffusion tensor mri. MICCAI, pages 117–125, 2006. 8 [5] P. Callaghan. Principles of Nuclear Magnetic Resonance Microscopy. Oxford University Press, USA, 1991. 2, 3 [6] Y. Chen, W. Guo, Q. Zeng, X. Yan, F. Huang, H. Zhang, G. He, B. Vemuri, and Y. Liu. Estimation, smoothing, and characterization of apparent diffusion coefficient profiles from high angular resolution dwi. CVPR, 1:588–593, 2004. 7 15 [7] T. Conturo, N. Lori, T. Cull, E. Akbudak, A. Snyder, J. Shimony, R. McKinstry, H. Burton, and M. Raichle. Tracking neuronal fiber pathways in the living human brain. pages 10422–10427. NAS of the USA, 1999. 4, 9 [8] O. Coulon, D. Alexander, and S. Arridge. A regularization scheme for diffusion tensor magnetic resonance images. pages 92–105, Davis,USA, 2001. ICIP. 4, 9 [9] R. Deriche, D. Tschumperle, and C. Lenglet. Dt-mri regularization and fiber tractography. Arlington, VA, 2004. ISBI. 4, 7, 9 [10] M. Descoteaux, E. Angelino, S. Fitzgibbons, and R. Deriche. Apparent diffusion coefficients from high angular resolution diffusion images: Estimation and applications. Magn. Res. in Med., 56:395– 410, August 2006. 4, 5, 6, 9, 11 [11] L. Frank. Characterization of anisotropy in high angular resolution diffusion-weighted mri. Magn. Res. Med, 47:1083–1099, 2002. 5 [12] G. Granlund and H. Knutsson. Signal Processing for Computer Vision. Kluwer Academic Publishers, 1995. 6, 11 [13] H. Gray. Anatomy of the human body. Lea and Febiger, 1918. 13, 15 [14] H. Gudbjartsson and S. Patz. The rician distribution of noisy mri data. Magn. Reson. Med., 34(6):910–914, 1995. 6, 11 [15] R. Henkelman. Measurement of signal intensities in the presence of noise in mr images. Medical Physics, 12(2):232–233, 1985. 6, 11 [16] D. LeBihan, E. Breton, D. Lallemand, P. Grenier, E. Cabanis, and M. Laval-Jeantet. Mr imaging of intravoxel incoherent motions: Application to diffusion and perfusion in neurologic disorders. Radiology, pages 401–407, 1986. 2 [17] E. Ozarslan, B. Vemuri, and T. Mareci. Fiber orientation mapping using generalized diffusion tensor imaging. In Biomedical imaging: macro to nano, volume 1, pages 1036–1039. International Symposium, apr 2004. 3 [18] O. Pasternak, N. Sochen, and Y. Assaf. Visualization and Image Processing of Tensor Fields. J.Weickert, H.Hagen, Springer,2005. 7 [19] C. Poupon, C. Clark, V. Frouin, J. Régis, I. Bloch, D. Le Bihan, and J. Mangin. Regularization of diffusion-based direction maps for the tracking of brain white matter fascicles. NeuroImage, 12:184– 195, 2000. 2 [20] W. Press, B. Flannery, S. Teukolsky, and W. Vetterling. Numerical recipes in C: the art of scientific computing, chapter Runge-Kutta Method, pages 710–714. Cambridge University Press, 1992. 10 [21] A. Ramirez-Manzanares and M. Rivera. A method for estimating brain nerve bundles, by restoring and filtering intra-voxel information in diffusion tensor mri data. Nice, France, oct 2003. VLSM. 7 [22] J. Sijbers, A. den Dekker, J. van Audekerke, M. Verhoye, and D. van Dyck. Estimation of the noise in magnitude mr images. Magnetic Resonance Imaging, 16(1):87–90, 1998. 6, 8, 11 [23] E. Stejskal and J. Tanner. Spin diffusion measurements: spin echoes in the presence of a timedependent field gradient. Journal of Chemical Physics, 42(1):288–292, 1965. 3 16 [24] J. Tournier, F. Calamante, D. Gadian, and A. Connely. Direct estimation of the fiber orientation density function from diffusion-weighted mri data using spherical deconvolution. NeuroImage, 23:1179– 1185, 2004. 3, 5 [25] D. Tuch. Diffusion MRI of Complex Tissue Structure. PhD thesis, Massachusetts Institute Of Technology, January 2002. 3 [26] D. Tuch. Q-ball imaging. Magn. Res. in Med., 52(6):1358–1372, 2004. 3, 4, 10 [27] B. Vemuri, Y. Chen, M. Rao, T. McGraw, Z. Wang, and T. Mareci. Fiber tract mapping from diffusion tensor mri. Vancouver, Canada, 2001. VLSM. 4, 9 [28] Z. Wang, B. Vemuri, Y. Chen, and T. Mareci. A constrained variational principle for direct estimation and smoothing of the diffusion tensor field from complex DWI. Transactions on Medical Imaging, pages 930–939, 2004. 7 [29] C. Westin, S. Maier, H. Mamata, A. Nabavi, F. Jolesz, and R. Kikinis. Processing and visualization for diffusion tensor mri. Medical Image Analysis, 6:93–108, 2002. 2 17