Estimation of fiber Orientation Probability Density ... - Semantic Scholar

NeuroImage 47 (2009) 638–650

Contents lists available at ScienceDirect

NeuroImage j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / y n i m g

Estimation of fiber Orientation Probability Density Functions in High Angular Resolution Diffusion Imaging Antonio Tristán-Vega a,⁎, Carl-Fredrik Westin b, Santiago Aja-Fernández c a b c

Laboratory of Image Processing, University of Valladolid, Spain Laboratory for Mathematics in Imaging, Harvard Medical School, USA Laboratory of Image Processing, University of Valladolid, Spain

a r t i c l e

i n f o

Article history: Received 19 December 2008 Revised 20 March 2009 Accepted 9 April 2009 Available online 22 April 2009 Keywords: HARDI Funk–Radon Transform Orientation distributions Angular PDF Spherical harmonics

a b s t r a c t An estimator of the Orientation Probability Density Function (OPDF) of fiber tracts in the white matter of the brain from High Angular Resolution Diffusion data is presented. Unlike Q-Balls, which use the Funk–Radon transform to estimate the radial projection of the 3D Probability Density Function, the Jacobian of the spherical coordinates is included in the Funk–Radon approximation to the radial integral. Thus, true angular marginalizations are computed, which allows a strict probabilistic interpretation. Extensive experiments with both synthetic and real data show the better capability of our method to characterize complex microarchitectures compared to other related approaches (Q-Balls and Diffusion Orientation Transform), especially for low values of the diffusion weighting parameter. © 2009 Elsevier Inc. All rights reserved.

Introduction Diffusion MRI is a recent MR modality that has made possible to characterize the white matter architecture in the brain in vivo (Basser et al., 1996). A popular model of the diffusion profile is based on the Gaussian assumption, which allows the diffusion to be modeled with a single covariance matrix, namely the diffusion tensor. With the advent of parallel imaging and better scanners, it is today feasible to acquire a larger number of diffusion weighted images in clinical time, the so-called High Angular Resolution Diffusion Imaging (HARDI). An advantage of HARDI is that it can better model more complex fiber architectures others than one single fiber bundle at each image voxel, such as crossing, bending, or kissing fibers. A number of HARDI techniques have recently appeared, and among them Q-Ball imaging, which is based on the Funk–Radon Transform (FRT, (Tuch, 2004)), has gained popularity mainly because it can be robustly computed with closed-form expressions (Descoteaux et al., 2007) and does not require any assumptions about the behavior of the diffusion signal outside of the sampled sphere (Campbell et al., 2005; Perrin et al., 2008; Poupon et al., 2008). Q-Balls use the FRT to approximate the Orientation Distribution Function (ODF) as the radial projection (integral) of the 3D Probability ⁎ Corresponding author. E-mail addresses: [email protected] (A. Tristán-Vega), [email protected] (C.-F. Westin), [email protected] (S. Aja-Fernández). 1053-8119/$ – see front matter © 2009 Elsevier Inc. All rights reserved. doi:10.1016/j.neuroimage.2009.04.049

Density Function (PDF). An alternative approach to represent the orientation distribution is to marginalize the radial part of the PDF. Although these two approaches at first seem equivalent, a major difference for the latter is that the Jacobian of the spherical coordinates is included in the radial integral defining the marginalization, computing an angular function which is a true marginal PDF with a valid probabilistic interpretation. An estimator of such a function (the Orientation Probability Density Function, OPDF) is presented here. Like Q-Balls, our estimator is based on the FRT, so only very weak assumptions have to be made about the behavior of the diffusion signal outside the sampled sphere. As an additional contribution, we propose a matrix implementation that allows to compute it in a very fast and robust way. Our preliminary results indicate that this probabilistic approach better resolves fiber crossings compared to other related approaches (Q-Balls and Diffusion Orientation Transform), especially for low values of the diffusion weighting parameter. Results on brain HARDI data also show promising results indicating the ability to visualize regions of crossing fibers. Background The generalized gradient diffusion equation An assumption inherent to diffusion imaging is that the probability of the displacement of a water molecule in a given angular direction is related to the presence of a fiber bundle in this same direction. Under

A. Tristán-Vega et al. / NeuroImage 47 (2009) 638–650

the assumption of narrow pulses, i.e., the duration of the sensitizing gradients δ is much smaller than the time between pulses Δ1, the PDF of the displacement of a water molecule to a position R ∈ R3 is given by the Fourier transform of the attenuation signal E(q) (Callahan, 1991): Z Z Z EðqÞexpð−2πiq  RÞdq ð1Þ P ðRÞ = FfEðqÞgðRÞ = ℝ3

where q = γδ = 2πG = qG = OGO with γ the gyromagnetic ratio, and g = G = OGO is the unitary direction of the magnetic field gradient. In Diffusion Spectral Imaging (DSI), the whole q-space or at least a representative subspace of ℝ3 is sampled, so the Fourier transform may be numerically computed to calculate P(R) and infer the underlying fiber population (Tuch et al., 2003; Wedeen et al., 2005). This technique requires the sampling of all directions g and all magnitudes q, which clearly reduces its clinical usefulness. Sphere sampling of the diffusion signal To avoid the need to sample the whole q-space, the main assumption of Diffusion Tensor Imaging (DTI) is that the PDF of the displacement may be accurately described by a Gaussian process, so Eq. (1) reduces to the well-known Stejskal–Tanner equation (Stejskal and Tanner, 1965): T

1 −R D P ðRÞ = qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi exp 4τ ð4πτÞ3 jDj

1

R

!

  T fEðqÞ = exp −bg Dg

ð2Þ

where τ = Δ − δ = 3 is the effective diffusion time and b = 4πτq2 is the diffusion weighting parameter. The symmetric, positive-definite covariance matrix D is the Diffusion Tensor (DT). The only unknowns are the six free components (due to the symmetry) of the DT, so it is enough to know the value of E(q) for six independent gradient directions g and one single q, since D does not depend on q. More gradient directions may be acquired to palliate the effect of noise, having a regular sampling of the sphere S ≡ {q = qg|4πτq2 = b0}. For a large number of sensitizing gradients, the term HARDI is often used. If only one single fiber bundle is present in each voxel, the Gaussian model is accurate enough; for F fiber bundles without water exchange, the linearity of the Fourier transform yields the multi-tensor model: EðqÞ =

F X

  T pf exp −bg Df g

ð3Þ

f =1

literature differ in the way they manage this incomplete information. A restricted model for E(q) may be assumed: multi-tensor models (Kreher et al., 2005; Peled et al., 2006) use Eq. (3) to characterize E (q); the number of fiber bundles F must be known beforehand or at least be estimated. Continuous Mixtures Of Gaussians (CMOG; Jian et al., 2007) overcome this problem assuming a continuous Wishart distribution for the parameters of the tensors. Higher order tensors (Descoteaux et al., 2006) are used to model the complex ADC for a given q; this yields the implicit assumption that the ADC does not depend on q, which is only true for the Gaussian model. The same assumption is used in the Diffusion Orientation Transform (DOT; Özarslan et al., 2006). In other approaches the estimation of P(R) is treated as an inverse problem: Persistent Angular Structures (PAS; Jansons and Alexander, 2003) assume that the diffusion of water molecules is restricted to a given radius R0(P(R) = Θ(θ, ϕ)δ(R − R0) where now δ is Dirac's delta function) and infer orientation information with minimum complexity; spherical deconvolution (Anderson, 2005; Descoteaux et al., 2009; Tournier et al., 2007, 2008) uses Eq. (3) for a given basic tensor (convolution kernel) in each spatial direction and infer a Point Spread Function p(θ, ϕ) for these directions. Finally, as shown in the Funk–Radon Transform section, Q-Balls obviate any assumption on the behavior of E(q) outside S at the expense of a blurring in the computation of orientation information (Tuch, 2004; Tuch et al., 2003). Inference of orientation information Although E(q) is in general not completely characterized in HARDI, it is not necessary either to completely characterize P(R) (the probability density of each displacement), but only its underlying orientation information (displacements in the same directions are associated to the same fiber bundle); for each orientation r given by (θ, ϕ), the expression for the marginal probability is: pðrÞupðθ; /Þ =

Z ∞

Z Z Z

ℝ3

Z

π

P ðRÞdR = Z

ð4Þ

where the positive function D is the Apparent Diffusion Coefficient (ADC) and (θ, ϕ) are the angular coordinates in the spherical system (we use physics convention, so g = ½sinθcos/; sinθsin/; cosθT . For the Gaussian model the positive function D(q, g) reduces to the positivedefinite quadratic form gTDg, so only in this case the ADC does not depend on q and the sampling of S is enough to completely characterize E(q). State of the art For non-Gaussian diffusion profiles E(q) is not completely characterized by HARDI data sets, so the previous approaches in the 1

2

P ðRrÞR sinθdR

Δ and δ will be used later on to denote the Laplacian operator and Dirac's delta function. We keep here this notation since it is the most widely accepted.

ð5Þ

0

with R = ‖R‖ so R = Rr. The term R2 sin θ is the Jacobian of the transformation to spherical coordinates and it is therefore required to compute actual probabilities so that the integral of P(Rr) equals 1:

0

=

where pf is the partial volume fraction of fiber f and Df is the corresponding diffusion tensor. For more complex architectures, more general expressions for E(q) have to be used. In general, we can always write:   2 EðqÞ = exp −4πτq Dðq; θ; /Þ b1

639

π

Z

π

Z

0

Z =

0

Z

2π 0 2π

Z ∞ 0

Z ∞ 0

0 2π

0

2

P ðRrÞR sinðθÞdRd/dθ  2 P ðRrÞR dR sinðθÞd/dθ

Φðθ; /Þ sinððθÞÞd/dθ = 1;

ð6Þ

where Φ(θ, ϕ) ≡ Φ(r) is the OPDF defined as: Φðθ; /ÞuΦðrÞ =

Z ∞

2

P ðRrÞR dR:

ð7Þ

0

In the last integral of Eq. (6) the term sin (θ) is once again due to the spherical coordinates system, i.e., it is required to compute the integral over the surface of the sphere. The orientation information for each direction r is characterized by the OPDF, which is a true probability density since its integral over all possible directions is 1. This definition has been already proposed before in Wedeen et al., (2005) as a “weighted radial summation”, although this approach, based on DSI, highly differs from the work here presented. Nevertheless, it turns out that it is usually easier to compute the ODF (Tuch et al., 2003) as the radial projection of P(R): Wðθ; /ÞuWðrÞ =

Z ∞ 0

P ðRrÞdR:

ð8Þ

640

A. Tristán-Vega et al. / NeuroImage 47 (2009) 638–650

Fig. 1. Examples of theoretically computed orientation functions for a simple tensor model: (a) Ψ(r), computed as C1(gTD− 1g)− 1/2; (b) Φ(r), computed as C2(gTD− 1g)− 3/2; and (c) ϒ(r), computed as C3 exp(−RTD− 1R / 4τ) with typical values R0 = 12 μm and τ = 20 ms.

The Jacobian of the spherical coordinates is dropped from the definition of Ψ(r), so this function does not represent a true Probability Density Function, as has been noted before (Tuch, 2004; Tuch et al., 2003). In practice, this has the effect of blurring the orientation information. Other orientation functions, such as the probability profile for a given radius R0, ϒ(r) ≡ P(R0r), have been used in the literature (Özarslan et al., 2006). In the case of the tensor model in Eq. (2) all these functions may be explicitly computed; in Fig. 1 we show rendered surfaces of Ψ(r), Φ(r), and ϒ(r) for a typical tensor configuration with eigenvalues [1.7, 0.3, 0.3] U 10− 3 mm2/s, where the value of the orientation function at each direction is represented as the distance of the surface to the origin. As may be seen, both Φ(r) and ϒ (r) have sharper profiles than Ψ(r). The blurring in the orientation information provided by the ODF (Ψ(r)) is a source of uncertainty in the estimation of fiber directions: two fiber bundles crossing in close directions would be more difficult to distinguish with the ODF than with the OPDF, since the two local maxima of the ODF could hinder each other. This statement will be thoroughly validated in the Results and discussion section.

The advantage of using the ODF relies on the fact that it may be accurately approximated by the Funk–Radon transform G of the attenuation signal, defined for a given unitary direction r as (Tuch, 2004; Tuch et al., 2003): Z q = q0

= 2πq0 gZ

Z ∞ −∞

δðr  qÞEðqÞdq Z ∞ Z ∞Z −∞

0

Spherical Harmonics Expansion Since both the attenuation signal E and the OPDF Φ are functions defined on a sphere, it is very convenient to represent them in the basis of Spherical Harmonics (SH), see for example Frank (2002). SH are defined as the eigenfunctions of the Laplace– Beltrami operator: Δb Y ðθ; /Þ =

The Funk–Radon Transform

GfEðqÞgðrÞJ

J0: the larger the intensity of the sensitizing gradients b0 = 4πτq20 , the narrower the Bessel kernel and the better the approximation, see Fig. 2. Z is a normalization constant. In this case only the values of E along the equators of the sphere S are required, so no assumption has to be made about the behavior of E outside S. This is the principle of Q-Ball imaging as introduced in Tuch, (2004); Tuch et al. (2003). This technique was designed to estimate the ODF but not the OPDF, so the aim here is to use the FRT to build an estimator of the OPDF. The blurring due to the Bessel kernel will be discussed in the Methods section.

2π 0

P ðRrÞdR = 2ZWðrÞ

and may be put in the form: m Yl ðθ; /Þ

mm l

P ðRrÞJ0 ð2πq0 ρÞρdudρdR ð9Þ

where ρ and φ are the coordinates of a cylindrical system with the z axis aligned with r: the integral of E(q) along the equator perpendicular to r (with radius q0) is proportional to the integral of P(R) inside a tube along r given by the lobes of the Bessel function

  2 1 A AY ðθ; /Þ 1 A Y ðθ; /Þ = mY ðθ; /Þ sin θ + 2 sin θ Aθ Aθ sin θ A/2 ð10Þ

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2l + 1 ðl − mÞ! m im/ = P ðcos θÞe 4π ðl + mÞ! l

= − lðl + 1Þ;

ð11Þ

m = − l N l; l = 0 N ∞

where Plm are the associated Legendre polynomials; since the functions we are representing have all radial symmetry, it is enough to consider even orders, l = 0, 2, 4…, see Descoteaux et al. (2006) and Appendix B for more details. It has been proven by Descoteaux et al. (2007) that SH are eigenfunctions as well for the FRT. As a final remark, it has been proven that the SH decomposition of the ADC is equivalent to its representation with higher order tensors (Descoteaux et al., 2006).

Fig. 2. Illustration of the FRT. The integral of E (which has been sampled in the sphere) along the equator q0 is equivalent to the integral of P inside the tubes along u; each tube represent a local extreme of the Bessel function J0 (right). The larger q0, the narrower the tube.

A. Tristán-Vega et al. / NeuroImage 47 (2009) 638–650

tensors or DOT, may be used to analytically compute the first derivative of E:

Methods Estimation of the OPDF from the attenuation signal Our goal is to estimate the OPDF (Φ(r)) as defined in Eq. (7). From Eq. (9), we know that the integral of a function along the direction r may be approximated by the integral on an equator perpendicular to r (i.e., the FRT evaluated at r) of its inverse Fourier transform. It is clear from Eq. (7) that the function to integrate is R2P(R), so: Z ΦðrÞ =

∞

2

R P ðRrÞdR =

0

1 2

Z

∞

−∞

2

R P ðRrÞdRg

1 n˜ o G EðqÞ ðrÞ 2Z

! A2 A2 A2 EðqÞ + + Aq21 Aq22 Aq23

n o −1 −1 2 R P ðRÞ ðqÞ = E˜ðqÞ = F 4π 2

ð13Þ

T

where now q = [q1,q2,q3] ; the right hand side of Eq. (13) is clearly the Laplacian of E. Using this result in Eq. (12), our estimator may be written:   1 −1 ΔE ð q Þ ðrÞ = C  GfΔEðqÞgðrÞ G 2Z 4π2

ð14Þ

for a negative constant C. Once again, due to the spherical symmetry of the problem, it makes sense to use the spherical coordinates representation of the Laplacian operator: ΔE = =

    1 A 1 A AE 1 A2 E 2 AE q + sin θ + 2 2 Aq 2 2 Aq Aθ Aθ q q sin θ q sin θ A/2 1 1 Δq E + 2 Δb E q2 q

ð15Þ

The problem reduces to estimate the Laplacian of the attenuation signal E(q) in the sampled sphere at b0 to compute its FRT, which in turn may be decomposed in a radial part and a part which is proportional to the Laplace–Beltrami operator. Expressions for the tensor model may be found in Appendix A. Estimation of Δb As said before, SH are eigenfunctions of the Laplace–Beltrami operator, so once the attenuation signal has been expressed in the basis of SH the second term in Eq. (15) is trivial to compute:

Eðq0 gÞ =

∞ X

l X

m

l=0 m= −l l even

=

∞ 1 X 2 q l=0

m

Cl Yl ðgÞZ

l X

1 Δb Eðq0 gÞ q2 m

m

lðl + 1ÞCl Yl ðgÞ

  AEðq; gÞ A 2 = exp −4πτq Dðq; gÞ Aq Aq     A = − 4πτq 2Dðq; gÞ + q Dðq; gÞ exp −4πτq2 Dðq; gÞ Aq   2 g − 8πτqDðq; gÞexp −4πτq Dðq; gÞ ; ð17Þ ð17Þ

ð12Þ

where we have used the approximation given by Eq. (9) to integrate R2 P(R). The unknown function E˜ is easy to determine, since it has to be the inverse Fourier transform of R2 P(R); from basic Fourier analysis theory:

ΦðrÞg

641

ð16Þ

with the assumption of a constant ADC, ∂D/∂q = 0 since it is assumed that D(q,g) = D(q0,g) = D(g) for all q; on the contrary, we only need to assume that ∂D / ∂q ≪ 2D / q for q0 so that we may neglect the derivative ∂D / ∂q and therefore estimate ∂E / ∂q from one single q0. Qualitatively speaking, this is equivalent to assume that the attenuation signal for a given spatial direction g changes much faster than the shape of the ADC at this particular direction does. This assumption is checked in the final part of the numerical validation. Finally, we may estimate Δb as:   Δq Eðq0 gÞg − 8πτq20 Dðq0 ; gÞ 3 − 8πτq20 Dðq0 ; gÞ Eðq0 gÞ = 2logEðq0 gÞð3 + 2logEðq0 gÞÞEðq0 gÞ

Issues in the computation of the FRT The blurring due to the Bessel kernel in the FRT (see the Funk– Radon Transform section) has an additional side effect in our case due to J0 is not positive (see Fig. 2, right). Consider the expression of the Laplacian for the tensor model given in Appendix A by Eq. (A.3) (for more complex models, the linearity of the problem may be used to infer identical conclusions); due to the positive-definite character of D, it may be possible to find that, for the directions of maximum diffusion (largest ‖Dq‖):   2 −8πτ −8πτODqO + traceðDÞ EðqÞN 0

ð19Þ

and so certain directions of q, once the factor −1 / (4π2) is applied, contribute to the FRT integral with negative values. This means that the integrals corresponding to certain directions r of the OPDF orthogonal to the directions q of maximum diffusion may be negative. Therefore, the modulus of the OPDT estimation has to be computed in all cases. The Orientation Probability Density Transform An overview of the whole estimation scheme may be seen in Fig. 3. The radial part of the Laplacian is expanded in the SH basis so that we may use the linearity of this expansion to compute the SH coefficients of the whole estimator. The synthesis formula of the SH, see Appendix B, may be used to compute the estimated probability density for any arbitrary direction r, not necessarily corresponding to one of the original sampled directions. As an analogy

m= −l

l even

where Clm are the coefficients of the SH expansion, computed as in Descoteaux et al. (2007), and in practice the sum is truncated to a few terms l ≤ L. Estimation of Δq The radial term in Eq. (15) is more difficult to compute since we have information only in the sphere of radius q0. Going back to Eq. (4), the assumption of a nearly constant ADC, underlying in higher order

ð18Þ

Fig. 3. Summary of the proposed OPDF estimator, the so-called OPDT.

642

A. Tristán-Vega et al. / NeuroImage 47 (2009) 638–650

Table 1 The coefficients fj−1 for different ratios of e2/e1 (Descoteaux et al., 2009 Eq. (13)); in the last row, fj−1 are computed based on the Laplace–Beltrami operator.

e2/e1 = 0.5 e2/e1 = 0.2 e2/e1 = 0.1 Δb

lj = 0

lj = 2

lj = 4

lj = 6

lj = 8

0.159 0.159 0.159 0

3.53 1.61 1.20 6

36.16 7.44 4.07 20

301.96 27.95 11.24 42

2292.1 95.41 28.23 72

The lj denotes the order of the SH (see Appendix B for details on notation).

with the DOT, we define the Orientation Probability Density Transform (OPDT) as the estimator described in Fig. 3. OPDT is based on the Q-Balls implementation by Descoteaux et al. (2007), so it is fast and robust, and may be represented as a set of matrix computations as described in Appendix B. On the other hand, OPDT may be seen as a kind of contrast-enhanced Q-Ball. OPDT is computed as the FRT of the diffusion signal corrected with the Laplacian operator, which has the effect of sharpening it. Given this interpretation, our work is related to the technique introduced in Descoteaux et al., (2009), where the attenuation signal is preconditioned to obtain a sharpened ODF, namely the fiber Orientation Distribution Function (fODF). This estimator may be written, with the notation used in Appendix B, as: Wsharp ðθ; /ÞuWsharp ðrÞ =

HX −1 j=0

fF gjj

Cj Y ðrÞ fj j

ð20Þ

where Cj are the coefficients of the SH expansion of E(q), Yj are each of the H SH basis functions and {F}jj are the coefficients resulting from the computation of the FRT of Yj; the sharpening effect is due to fj, which come from the deconvolution of the ODF with a Gaussian (tensor) kernel, i.e., the shape of the ODF in case one single fiber direction is present. These terms depend on the ratio e2/e1 between the eigenvalues of the tensor used as deconvolution kernel, which has

to be estimated from data. Now, neglecting the radial part Δq of the Laplacian (as described later on, this is a good approximation is some cases), OPDT fits Eq. (20), with coefficients fj−1 which are the eigenvalues of the SH basis for the Laplace–Beltrami operator and therefore do not need to be computed. To give a deeper insight in the performance of OPDT as a contrast enhancement of Q-Balls, in Table 1 we show the values of fj−1 for several configurations. The SH coefficients of higher orders are multiplied by a factor greater than 1 to enhance fast variations on the diffusion profile, hence improving the angular contrast. Nevertheless, our approach strongly differs from the one in Descoteaux et al. (2009); note first that OPDT may be casted into Eq. (20) only if the radial part of the Laplacian is ignored, which is not always valid (this topic is discussed in the experimental section). Otherwise, OPDT performs a highly non-linear correction of the attenuation signal. Second, note that the orientation function estimated by the OPDT, the OPDF, has a probabilistic interpretation, while the estimator in Eq. (20) does not. Third, note that we do not need to estimate any deconvolution kernel, which makes the work by Descoteaux et al. (2009) closer to the spherical deconvolution method in Anderson (2005) (based as well in SH expansions) than it is to Q-Balls, DOT or OPDT. Results and discussion We study the capability to resolve complex micro-architectures in different scenarios. The parameters of interest are the value of b0, the number of sampling gradients N in the sphere S and the noise power σ 2n. We compare OPDT with Q-Balls to show the different representation capabilities of ODF and OPDF. We test the DOT as well; although it is not an estimator neither for the ODF nor for the OPDF (DOT estimates ϒ(r) ≡ P(R0r) for a given R0), this estimator is closely related to OPDT and Q-Balls (see the discussion in the previous section). SH expansions up to order L = 6 (H = 28 SH coefficients) have been taken. For the regularization parameter (see Appendix B) we use λ = 0.006 (Descoteaux et al., 2007).

Fig. 4. Mean error in the angle of the detected fibers as a function of the original angle between fibers. The curves stop when the estimators are not able to detect the fiber crossings (there are not two maxima in the orientation function).

A. Tristán-Vega et al. / NeuroImage 47 (2009) 638–650

For DOT we use R0 = 12 μm. Although a higher value (R0 = 25 μm) may yield more accurate estimates for noise-free scenarios (Özarslan et al., 2006); (Prčkovska et al., 2008), in the presence of noise its performance is highly degraded. According to our experience, R0 = 12 μm is the best trade-off between accuracy and robustness to noise in a real situation, so we keep this value in all our experiments to achieve a fair comparison. For the noisy scenarios, we corrupt the simulated signal E(q) with Rician noise (Drumheller, 1993). Although the attenuation signal is computed as the quotient between the DWI signals and the baseline image, i.e., as the quotient between two Rician-distributed signals (Gudbjartsson and Patz, 1995), DWI are always much more noisy than the baseline; moreover, several baseline images are often available in HARDI data sets, so we consider the baseline is almost noise-free and therefore E(q) is closely enough to Rician distributed. Capability of resolving two crossing fibers The capability of correctly resolving two crossing fibers and the accuracy in the determination of their directions is a standard benchmark for the assessment of micro-architecture resolution capabilities (Descoteaux et al., 2007; Özarslan et al., 2006; Prčkovska et al., 2008; Tournier et al., 2007, 2008; Tuch, 2004). To generate the attenuation signal, we use two synthetic tensors with eigenvalues [1.7, 0.3, 0.3] U 10− 3 mm2/s and principal eigenvectors in the XY plane: for the first tensor, the principal eigenvector forms an angle of (90 − α) / 2° with the X axis, and for the second one an angle of (90 − α) / 2° with the Y axis, so that the directions of both fiber bundles form an angle of α degrees. Eq. (3) is used to generate the attenuation signal for N antipodal symmetric directions, with partial volume fractions pf = 1/2 each. Five different scenarios are considered: S-1, with b0 = 1200 s/mm2 and N = 50; S–2, with b0 = 1200 s/mm2 and N = 100; S–3, with b0 = 3000 s/mm2 and N = 50; S–4, with b0 = 3000 s/mm2 and

643

N = 100; S–5, with b0 = 3000 s/mm2 and N = 200. Results may be found in Fig. 4. We represent the average angular error in the detected maxima only if the two maxima are found. OPDT is able to recover closer crossing fibers in all cases, and the difference is more important with lower b0 values. OPDT is more accurate than Q-Balls in all cases. DOT is more accurate than OPDT only for higher values of b0 and fibers crossing in large angles, and yet its advantage is quite subtle. Comparing Q-Balls and DOT, their performance is quite similar for low values of b0, but DOT performs better for higher b0 (similar results have been reported in Prčkovska et al. (2008)). Note that all estimators are able to better resolve crossings in closer directions for higher b0, since a greater b0 augments the contrast in the diffusion profile. On the contrary, using more gradient directions improves the results only marginally: the minimum detectable angle is practically the same, and the improvement on the angular error is negligible. As a final comment, note the anomalous behavior of OPDT near 45°, when the angular error decreases as the fibers bundles get closer; this phenomenon is thoroughly discussed later on in this Section. As an additional experiment, we show in Fig. 5 the result of the three estimators for scenario S-4. OPDT yields more accurate estimates and is able to detect the two fibers when the other estimators show only one maximum (Q-Balls for 45° and 40°) or a maximum which does not correspond to any of the true fiber populations (DOT for 45° and 40°). From this section we may conclude that OPDT is preferable: first to resolve fibers crossing in low angles, and second for low values of b0; the third conclusion is that the OPDF is able to yield more accurate results and a higher angular contrast in most of cases. Behavior in the presence of noise Although conventional MRI shows very good Signal to Noise Ratio (SNR), this is not the case in DTI/HARDI data sets; moreover,

Fig. 5. Estimates of the fiber population for S-4 with Q-Balls (top), DOT (middle) and OPDT (bottom), for crossing angles, from left to right: 90°, 60°, 45° and 40°. Green lines correspond to the ground-truth, while red lines correspond to the estimated directions.

644

A. Tristán-Vega et al. / NeuroImage 47 (2009) 638–650

Fig. 6. Mean error in the angle of the detected fibers as a function of the original angle between fibers in the presence of noise. Results are the average of 100 Montecarlo trials. In this case we consider a failure when the estimator is not able to detect the two fibers in the 50% of the trials.

higher values of b0 yield stronger attenuation which for the same noise power σ 2n produces poor SNR. In this section we study the effect of varying the Peak Signal to Noise Ratio (PSNR, defined here as the maximum value of the baseline divided by σn) of the

Rician-distributed signal E(q) in the accuracy of the estimators. Analogous results to that of Fig. 4 may be found in Fig. 6 for S-1 (PSNR = 13,33), S-2 (PSNR = 13,33), S-4 (PSNR = 5) and S-5 (PSNR = 5).

Fig. 7. Estimates of the fiber population for S-4 with Q-Balls (top), DOT (middle) and OPDT (bottom), for crossing angles, from left to right: 90°, 75°, 65° and 60°. Rician noise with PSNR = 5 has been added to E(q).

A. Tristán-Vega et al. / NeuroImage 47 (2009) 638–650

645

Fig. 8. Estimates of the fiber population with three crossing fibers for S-2 with PSNR = 13.3 (top) and S-5 with PSNR = 5 (bottom).

Obviously, noise worsens the accuracy of all the estimators. A number of techniques to palliate this effect have been successfully used in the literature (Clarke et al., 2008; Tristán-Vega and AjaFernández, 2008). In this case using more gradient directions improves both the minimum angle detectable and the accuracy. Note that increasing b0 improves the result even when the SNR dramatically decreases. The advantage of using OPDT is evident for lower b0; for higher b0, both OPDT and DOT show the same capability to

resolve fiber crossings, but OPDT is slightly more accurate. Contrary to the behavior shown in Fig. 4, Q-Ball is now more accurate for higher b-values, especially for N = 200. Qualitatively speaking, Q-Ball works by averaging the attenuation signal in a whole equator perpendicular to the direction of interest; therefore this estimator is indeed more robust to noise than OPDT, which is based on derivatives computed on the attenuation signal, accentuating the effect of noise. Besides, it has been previously described that there is a side effect in OPDT for noisy environments: the

Fig. 9. Pseudo-color representation of (2D/q)/(∂D/∂q) (error in the estimation of Δq, left) and Δb/Δq (relative importance of Δb, right) for 2 (top) and 3 (bottom) crossing fibers. We represent half the (θ, ϕ) space since the attenuation signal is symmetric. Fiber directions are represented as black spots, and values greater than 10 have been clipped to this value in all cases.

646

A. Tristán-Vega et al. / NeuroImage 47 (2009) 638–650

than the effect of noise, and both OPDT and DOT perform better than Q-Ball for all crossing angles. First, note that a b0 value of 1200 s/mm2 is more realistic in practical applications than 3000 s/mm2, so the advantage of OPDT remains clear; second, for b0 = 3000 s/mm2 OPDT has an accuracy close to that of Q-Balls with the same capability as DOT to resolve crossings. An analogous representation to that of Fig. 5 is shown in Fig. 7 for scenario S-4 with PSNR = 5. Consistently with the results in Fig. 6, Q-Balls give closer estimates to the true fibers, but its representation of the fiber populations is less intuitive and more blurred, failing to find the two populations before DOT or OPDT do, so the third important point is the better representation capability of OPDT and DOT. Resolution of more complex architectures Fig. 10. Relative error (averaged for N = 100 gradient directions) in the estimation of the Laplacian as a function of the angle between 2 crossing fibers.

OPDF estimated may yield negative values which have to be corrected by taking the modulus. This is a clear limitation of OPDT. Nevertheless, the signal average inherent to Q-Balls produce an angular blurring of the orientation information, which justifies the fact that OPDT is able to recover lower crossing angles even in noisy scenarios. For low b-values, this blurring is more important

For illustrative purposes we show in Fig. 8 a situation with three crossing fiber bundles; directions are: (θ1, ϕ1) = (π/2, 0), (θ2, ϕ2) = (π/2, π/2) and (θ3, ϕ3) = (π3/20, π/2); eigenvalues of Di are: [2, 0.2, 0.3]·10− 3, [1.8, 0.4, 0.3]·10− 3 and [2, 0.1, 0.1]·10− 3 mm2/s. We test scenarios S-2 with PSNR = 13.3 and S-5 with PSNR = 5. Q-Balls completely blurs the orientation information and are not able to distinguish the three fibers. DOT fails to recover the three fibers for S-2, but for S-5 it is able to detect all of them with a mean error of 18°. OPDT is able to yield acceptable results in both situations, with respective mean errors of 10° and 16°.

Fig. 11. Axial slice in the upper brain (the FA map of the whole slice is represented for illustrative purposes). Several important tracts may be identified: the superior longitudinal fasciculus (slf), the superior corona radiata (scr), the cingulus (cg) and the corpus callosum (cc). Glyph surfaces are colored here with the common convention of color per orientation.

A. Tristán-Vega et al. / NeuroImage 47 (2009) 638–650

647

Fig. 12. Detail of the regions highlighted in Fig. 11. OPDT is able to detect fiber crossings even for b0 = 700 s/mm2; note that region 2 is especially difficult due to the small angles of crossing.

We may conclude that OPDT is more accurate and robust as well for complex micro-architectures. As mentioned before, its advantage is clear for lower b0 values. Accuracy in the approximation of the Laplacian We show in Fig. 9 a pseudo-color representation of (2D/q)/ (∂D/∂q) (see Eq. (17)) and Δb/Δq (see Eq. (15)), computed as described in Appendix A for 2 and 3 crossing fibers. We show averaged values of these quotients for all bo between 1000 and 3000 s/mm2. In most of the (θ, ϕ) space (2D/q)/(∂D/∂q) is greater than 10, and besides Δb represents more than the 90% of the total value of Δ, so the error is negligible. The greatest errors (small (2D/q)/(∂D/∂q)) are committed in the directions of the fiber bundles. However, with 2 fibers these directions correspond to

zones where Δb is clearly greater, so the overall error is small. With 3 fibers there are directions where (2D/q)/(∂D/∂q) is small and Δb/Δq is less than 10, but minima of (2D/q)/(∂D/∂q) never overlap with minima of Δb/Δq. Therefore, we may conclude first that more complicated micro-architectures yield greater errors, and second that the relative error due to Δq remains moderated compared to

the total value of Δ = q20 Δq + Δb . To support these conclusions, we show in Fig. 10 the relative error in the estimation of the Laplacian as a function of the angle between 2 crossing fibers. The error decreases in all cases as the fibers get closer, since the structure of E(q) is simpler. Compared to Fig. 4, note that with OPDT the angular error anomaly decreases for fibers in angles below 50° or 60°, which may be explained by the fact that from 50° and below the relative error in Δq dramatically decreases. When the fibers get too close, however, the estimator is

Fig. 13. Axial slice of the midbrain, where the middle cerebellar peduncle (mcp), the pontine crossing tract (pct), the corticopontine tract (cpt), the corticospinal tract (cst), the medial lemnicus (ml) and the inferior cerebellar peduncle (icp) may be identified.

648

A. Tristán-Vega et al. / NeuroImage 47 (2009) 638–650

not able to resolve them and the error increases not depending on the error in Δq. Note that the error in Δq is obviously more important with lower b0: our assumption implies that the exponential decay of E(q) is much faster than the change in the shape of the ADC, and decreasing b0 slows the decay of E(q). Nevertheless we have previously shown that OPDT is more advantageous for lower b0, so we may conclude that the error in the estimation of Δq is negligible; moreover, note that in the noisy case (see Fig. 6) the aforementioned artifact near 50°–60° is not present, so we may infer that the effect of noise is more important than the error in the approximation of Δq itself. In vivo experiments To test the behavior of the proposed estimator on real data, we use a HARDI volume of a volunteer, scanned in a 1.5 T GE Echospeed system (scanning sequence: Max. gradient amplitudes: 40 mT/M. Rectangular FOV 220 × 165 mm. 128 × 96 scan matrix (256 × 192 image matrix). 4 mm slice thickness, 1 mm interslice distance. Receiver bandwidth 6 kHz. TE 70 ms; TR 80 ms (effective TR 2500 ms). Scan time 60 s/slice). It comprises 8 non-weighted baseline images and 51 independent gradient directions with diffusion weighting parameter b0 = 700 s/mm2. We deliberately choose such a small b value to demonstrate the representation capabilities of OPDT even for very extreme conditions. The data set has been denoised with the filter proposed in Tristán-Vega and Aja-Fernández (2008). In Fig. 11 a representative slice is shown; as can be seen, OPDT is able to yield anatomically correct results (Mori et al., 2005): the longitudinal fasciculus (green, left hand) goes longitudinally from top to bottom, as well as the cingulus (green, right hand), while the corpus callosum (red, right hand) crosses the cingulum from left to right and the corona radiata (blue, center) is perpendicular to the plane of the slice. More interestingly, OPDT is able to resolve the fiber crossings inside the highlighted zones (the horizontal red parts of the longitudinal fasciculus cross transversely with its longitudinal part in 1 and 2, and the corpus callosum crosses the cingulus in 3. These regions are shown in detail in Fig. 12.

To finish, we show in Fig. 13 an axial slice of the midbrain, where the fiber configuration is especially complicated. The main orientations given by OPDT are still anatomically correct (Mori et al., 2005), and yet the especially complicated crossings between the middle cerebellar peduncle and the pontine crossing tract in the highlighted zone are correctly detected. A detailed view of this zone is shown in Fig. 14. Conclusion We have introduced a simple, robust, yet very accurate estimator, the OPDT, for the probability density of fiber populations in a given direction. Its complexity is comparable to that of Q-Balls or DOT, but it can better resolve fiber crossings and represent fiber populations in many situations. Besides, it has a correct probabilistic interpretation: Q-Balls only estimate radial projections of true probability densities (the ODF), and DOT estimates only the value of the probability density at a given radius R0. On the contrary, OPDT is an estimator of true marginal probability densities (OPDF), which turns out to be an important property in the representation of fiber populations. Compared to other related methods, the main assumption of our approach, i.e., the slow variation of the ADC, is rather conservative and has been proven not to be a limitation (yet its effect is hindered by the effect of noise), which gives OPDT the advantage of Q-Balls of not needing to assume any kind of underlying model, therefore showing general validity whenever Eq. (1) holds. The proposed estimator is more advantageous for low b0 values, so its computation is not very restrictive on the machinery used to acquire the diffusion images. In fact, the main utility of OPDT is in scenarios of low b0. On the contrary, for higher values of b0 the accuracy in the estimation is not so clearly improved by OPDT compared to Q-Balls or DOT. In this sense, OPDT has two main drawbacks: • The computation of the FRT produces negative values of the orientation information which have to be corrected by taking the modulus. This is an issue especially for very noisy scenarios, for example in images acquired with a high b0 where the SNR is highly degraded. In these cases Q-Balls may be preferable, since their accuracy is similar to that of OPDT for high b0 and they are more robust to noise due to the FRT averaging. • Although only in a local sense, we still need to make an assumption on the behavior of the ADC outside the sampled sphere. This assumption is very weak but compared to Q-Balls, for which absolutely no assumption is required, it is still a drawback. As a final remark, the value of R0 for DOT used in our experiments is not optimal for noise-free scenarios, where taking R0 ⋍ 25 μm may highly improve the accuracy of the estimation for high b0. This improvement may even outperform OPDT in some situations. The problem is that using a high b0 (which is necessary for DOT to work properly) results in an exponential worsening of the SNR. As a consequence, these two conditions (high b0 and high SNR) are not compatible in general, unless an adequate denoising scheme (whose necessity in HARDI has been widely reported previously; Clarke et al., 2008; Descotaux et al., 2008) is implemented. Acknowledgments

Fig. 14. Detail of the highlighted region in Fig. 13.

The authors would like to acknowledge grant number TEC200767073/TCM from the Comisión Interministerial de Ciencia y Tecnología (Spain) and National Institute of Health grant NIH R01 MH074794 (CFW). Antonio Tristán-Vega would like to especially acknowledge Dr. Yogesh Rathi for his valuable help in the preparation of the experiments with real data.

A. Tristán-Vega et al. / NeuroImage 47 (2009) 638–650

A. Expressions of the Laplacian for the tensor model A.1. Laplacian in Cartesian coordinates With the notation in Eq. (13), the Stejskal–Tanner equation may be written in the form:   EðqÞuEðq1 ; q2 ; q3 Þ = exp −4πτqT Dq   = exp −4πτ½q1 ; q2 ; q3 D½q1 ; q2 ; q3 T

ðA:1Þ

and so, due to the symmetry of the diffusion tensor:   A T T EðqÞ = − 4πτ ½1; 0; 0Dq + q D½1; 0; 0 EðqÞ Aq1 = − 8πτ ½1; 0; 0DqEðqÞ   A 2 T EðqÞ = − 8πτ −8πτ ð½1; 0; 0DqÞ + ½1; 0; 0D½1; 0; 0 EðqÞ 2 Aq1 2

ðA:2Þ Computing the corresponding terms for q2 and q3, the expression for the Laplacian is:   2 ΔE = − 8πτ −8πτODqO + traceðDÞ EðqÞ

ðA:3Þ

A.2. Radial projection of the Laplacian From Eq. (2) and taking into account that b = 4πτq2, Δq, may be written as:   A 2 T q2 EðqÞ = q −8πτqg Dg EðqÞ Aq    2 2 T 2 T Δq EðqÞ = q −8πτqg Dg − 24πτq g Dg EðqÞ

649

P H = ðL = 2 + 1ÞðL + 1Þ. The set of N sampling direcH= l = 0L 2l tions of E(q) is denoted as Θ = fðθ1 ; /1 Þ; ðθ2 ; /2 Þ; N ; ðθN ; /N Þgu  g1 ; g2 ; N ; gN ; the set of N' sampling directions for n    which  the ODF  θ1V ; /e1 ; θ2V ; /2V ; N ; θNV V ; /NV V gu or the OPDF is computed is Θ V = fr1 ; r2 ; N ; rN V g. Therefore, we may define the N × 1 vectors E (attenuation signals) and D (log-signal); the R × 1 vectors C (SH coefficients of E or ΔE) and C' (the SH coefficients of the ODF or the OPDF); the N' × 1 vector P (the ODF or OPDF at each direction in Θ') as:



fEgi = E q0 gi ; fDgi = logE q0 gi ; ðN × 1Þ n o V = CjV ; ðH × 1Þ fCgj = Cj ; C j   V fPgi V = Wðri V Þ or ΦððÞÞ; N × 1

We define the matrix of eigenvalues of the SH basis L and the FRT matrix F as H × H diagonal matrices of the form (Descoteaux et al., 2007):   fLgjj = − lj lj + 1 ; fF gjj 8 1; > > <   = 1  3  5  N  lj − 1 > lj = 2 > : 2πð −1Þ ; 2  4  6  N  lj

  = 8q2 πτgT Dg 8πτqT Dq − 3 EðqÞ

j=0 ðB:4Þ jN0

Finally, we define the matrix B (resp. B') as the N × H (resp. N' × H) matrix of the evaluations of the H SH basis functions in the set Θ (resp. Θ'): n o   V V V fBgij = Yj ðθi ; /i Þ; B V = Yj θi V ; /i V ; i j

ðA:4Þ

ðB:3Þ

ðB:5Þ

B.2. Q-Balls computation The SH analysis equation is posed by (Descoteaux et al., 2007) as a regularized Least Squares (LS) problem, so:

Note that this last equation is formally identical to the approximation in Eq. (18). Moreover, for the tensor model q2 Dðq; gÞ = q2 gT DgZ Dðq; gÞ = DðgÞ = gT Dg (D does not depend on q), and Eq. (18) is the exact expression of Δq. B. Matrix computation of the OPDT The analysis and synthesis equations of the SH expansion, as well as the expression of the FRT, may be written as matrix operations (Descoteaux et al., 2007). In this Appendix we extend this work to the OPDT; we derive the matrix computation of the Laplacian and then exploit the properties of the resulting matrices to give a closed-form representation.

  T 2 −1 T B E C = B B + λL

ðB:6Þ

where λ is the regularization parameter; the product with L2 is used to penalize higher values of the Laplacian, i.e. higher energies of the second derivatives. Since SH are eigenfunctions for the FRT, the following expression holds:   V T 2 −1 T B E C = F C = F B B + λL

ðB:7Þ

It only remains to evaluate the SH expansion at the directions in Θ' to obtain:   V V V T 2 −1 T B E P = B C = B F B B + λL

ðB:8Þ

B.1. Notation Like in Descoteaux et al. (2007), we change the indexation of the SH basis to represent them with one single index j. For each even order l, it is easy to note that we have Hl = 2l + 1 basis functions Ylm , so the following indexation may be used:   n o l l −l l Yj ; j = ðl − 1Þ N ðl + 3Þ = Yl N Yl ; l even ðB:1Þ 2 2 Besides, we use the notation: lj = l :

l l ðl − 1ÞVjV ðl + 3Þ; l even 2 2

ðB:2Þ

for the order of the j-th SH. It is easy to note that the total number of SH basis functions for an expansion up to order L = 2L̂ is

Note that the whole Eq. (B.8) may be precomputed for all voxels except for the product with E, so the estimation of Ψ(r) requires only a product of an N' × N matrix by an N × 1 vector. B.3. OPDT computation From the scheme in Fig. 3 and the LS representation in Eq. (B.6), together with the definition of D in Eq. (B.4) it follows: " #   T 2 −1 T 2 B D  ð 3 + 2D Þ  E C = B B + λL q20   1 T 2 −1 T B E + 2 L B B + λL q0

ðB:9Þ

650

A. Tristán-Vega et al. / NeuroImage 47 (2009) 638–650

where the dot products · must be understood as element-wise operations. The first term in Eq. (B.9) is the SH decomposition of the approximation of Δb as in Eq. (18). In the second term, we first compute the SH expansion of E and then use the property of the SH of being eigenfunctions of the Laplace–Beltrami operator to compute Δb. Taking into account that both L and the LS matrix (BTB + λL2)− 1 are symmetric matrices and therefore they commute, we may write: C=

   1 T 2 −1 T T B B + λL 2B ½D  ð3 + 2DÞ  E + LB E q20

ðB:10Þ

And finally the estimator of Φ(r) is written: P=

    1 0 T 2 −1 T T B F B B + λL 2B ½D  ð3 + 2DÞ  E + LB E ðB:11Þ 2 2 4π q0

Note that we do not have to know the exact value of q0 if we normalize Φ(r), as opposed to the case of the DOT, see (Özarslan et al., 2006, Appendix A, Eq. (28)). The N' × H matrix B'F (BTB + λL2)− 1 and the N × H matrix LBT may be precomputed, so for each voxel we have to compute N logarithms, 2N products (the radial part of the Laplacian), two products of an N × H matrix by an N × 1 vector, and a product of an N' × 1 matrix by an N × 1 vector. References Anderson, A.W., 2005. Measurement of fiber orientation distributions using high angular resolution diffusion imaging. Magn. Reson. Med. 54 (5), 1194–1206. Basser, P.J., Pierpaoli, C., 1996. Microstructural and physiological features of tissues elucidated by quantitative-diffusion-tensor MRI. Magn. Reson. Med. 111 (3), 209–219. Callahan, P.T., 1991. Principles of Nuclear Magnetic Resonance Microscopy. Clarendon Press, Oxford. Campbell, J.S., Siddiqi, K., Rymar, V.V., Sadikot, A.F., Pike, G.B., 2005. Flow-based fiber tracking with diffusion tensor and Q-ball data: validation and comparison to principal diffusion direction techniques. NeuroImage 27 (4), 725–736. Clarke, R.A., Scifo, P., Rizzo, G., Dell'Acqua, F., Scotti, G., Fazio, F., Sep. 2008.. Noise correction on Rician distributed data for fibre orientation estimators. IEEE Trans. Med. Imag. 27 (9), 1242–1251. Descoteaux, M., Angelino, E., Fitzgibbons, S., Deriche, R., 2006. Apparent diffusion profile estimation from high angular resolution diffusion images: estimation and applications. Magn. Reson. Med. 56 (2), 395–410. Descoteaux, M., Angelino, E., Fitzgibbons, S., Deriche, R., 2007. Regularized, fast, and robust analytical Q-Ball imaging. Magn. Reson. Med. 58, 497–510. Descotaux, M., Wiest-Daesslé, N., Prima, S., Barillot, C., Deriche, R., 2008. Impact of Rician adapted non-local means filtering on HARDI. Medical Image Computing and Computer-Assisted Intervention. Vol. 5242 of Lecture Notes in Computer Science. Springer, Verlag, pp. 122–130.

Descoteaux, M., Deriche, R., Knösche, T.R., Anwander, A., Feb. 2009. Deterministic and probabilistic tractography based on complex fibre orientation distributions. IEEE Trans. Med. Imag. 28 (2), 269–286. Drumheller, D., Apr. 1993. General expressions for Rician density and distribution functions. IEEE Trans. Aerosp. Electron. Syst. 29 (2), 580–588. Frank, L.R., Jan. 2002.. Characterization of anisotropy in high angular resolution diffusion-weighted MRI. Magn. Reson. Med. 47 (6), 1083–1099. Gudbjartsson, H., Patz, S., 1995. The Rician distribution of noisy MRI data. Magn. Reson. Med. 34, 910–914. Jansons, K.M., Alexander, D.C., 2003. Persistent angular structures: new insights from diffusion magnetic resonance imaging data. Inverse Problems 19, 1031–1046. Jian, B., Vemuri, B., Özarslan, E., Carney, P.R., Marecy, T.H., Aug. 2007. A novel tensor distribution model for the diffusion-weighted MR signal. NeuroImage 37 (1), 164–176. Kreher, B., Schneider, J., Mader, I., Martin, E., Henning, J., Il'yasov, K., Sep. 2005. Multitensor approach for analysis and tracking of complex fiber configurations. Magn. Reson. Med. 54 (5), 1216–1225. Mori, S., Wakana, S., Nagae-Poetscher, L.M., van Zijl, P.C., 2005. MRI Atlas of Human White Matter. Elsevier, Amsterdam, The Netherlands. Özarslan, E., Sepherd, T.M., Vemuri, B.C., Blackband, S.J., Mareci, T.H., 2006. Resolution of complex tissue microarchitecture using the diffusion orientation transform (DOT). NeuroImage 31, 1086–1103. Peled, S., Friman, O., Jolesz, F., Westin, C.-F., 2006. Geometrically constrained two-tensor model for crossing tracts in DWI. Magn. Reson. Med. 24 (9), 1263–1270. Perrin, M., Cointepas, Y., Cachia, A., Poupon, C., Thirion, B., Riviére, D., Cathier, P., Kouby, V.E., Constantinesco, A., Bihan, D.L., Mangin, J.-F., 2008. Connectivity-based parcellation of the cortical mantle using Q-ball diffusion imaging. J. Biomed. Imaging 8 (3), 1–10. Poupon, C., Roche, A., Dubois, J., Mangin, J.-F., Poupon, F., 2008. Real-time MR diffusion tensor and Q-ball imaging using Kalman filtering. Med. Image Anal. 12 (5), 527–534 special issue on the 10th international conference on medical imaging and computer assisted intervention — MICCAI 2007. Prčkovska, V., Roebroeck, A., Pullens, W., Vilanova, A., ter Haar Romeny, B., Sep. 2008. Optimal acquisition schemes in high angular resolution diffusion weighted imaging. Medical Image Computing and Computer-Assisted Intervention. Vol. 5242 of Lecture Notes in Computer Science. Springer, Verlag, pp. 9–17. Stejskal, E.-O., Tanner, J.-E., 1965. Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient. J. Chem. Phys. 42, 288–292. Tournier, J.-D., Calamante, F., Connelly, A., 2007. Robust determination of the fibre orientation distribution in diffusion MRI: non-negativity constrained superresolved spherical deconvolution. NeuroImage 35, 1459–1472. Tournier, J.-D., Yeh, C.-H., Calamante, F., Cho, K.-H., Connelly, A., Lin, C.-P., 2008. Resolving crossing fibres using constrained spherical deconvolution: validation using diffusion-weighted imaging phantom data. NeuroImage 42, 617–625. Tristán-Vega, A., Aja-Fernández, S., Sep. 2008. Joint LMMSE estimation of DWI data for DTI processing. Medical Image Computing and Computer-Assisted Intervention. Vol. 5242 of Lecture Notes in Computer Science. Springer, Verlag, pp. 27–34. Tuch, D.S., 2004. Q-Ball imaging. Magn. Reson. Med. 52, 1358–1372. Tuch, D.S., Reese, T.G., Wiegell, M.R., Wedeen, V.J., 2003. Diffusion MRI of complex neural architecture. Neuron 40, 885–895. Wedeen, V.J., Hagmann, P., Tseng, W.-Y.I., Reese, T.G., Weisskoff, R.M., 2005. Mapping complex tissue architecture with diffusion spectrum imaging. Magn. Reson. Med. 54, 1377–1386.