Skip to main content
NCBI home page
NIHPA Author Manuscripts logoLink to NIHPA Author Manuscripts
. Author manuscript; available in PMC: 2010 Mar 1.
Published in final edited form as: Neuroimage. 2008 Nov 13;45(1 Suppl):S111–S122. doi: 10.1016/j.neuroimage.2008.10.054

Mathematical Methods for Diffusion MRI Processing

C Lenglet a,b,*,1, JSW Campbell c,d, M Descoteaux e, G Haro f, P Savadjiev d, D Wassermann g,h, A Anwander i, R Deriche g, GB Pike c, G Sapiro b, K Siddiqi d, P Thompson j
PMCID: PMC2678879  NIHMSID: NIHMS100395  PMID: 19063977

Abstract

In this article, we review recent mathematical models and computational methods for the processing of diffusion Magnetic Resonance Images, including state-of-the-art reconstruction of diffusion models, cerebral white matter connectivity analysis, and segmentation techniques. We focus on Diffusion Tensor Images (DTI) and Q-Ball Images (QBI).

Keywords: diffusion MRI, Diffusion Tensor Imaging (DTI), High Angular Resolution Diffusion Imaging (HARDI), Q-Ball Imaging (QBI), Orientation Distribution Function (ODF), Spherical Harmonics (SH), Funk-Radon Transform (FRT), sub-voxel fiber configurations, tractography, manifold learning, N-cuts, stratification, segmentation, non-uniform complexity

1. Introduction

Diffusion MRI is an exciting extension of MRI introduced in the mid-1980s (Le Bihan and Breton, 1985; Merboldt et al., 1985; Taylor and Bushell, 1985). It provides a unique and sensitive probe for the architecture of biological tissues. The key idea behind diffusion MRI is that water diffusion in structured tissues, such as the brain’s white matter, reflects its architecture at a microscopic scale. This is because molecular motion is favored in directions aligned with fiber bundles and hindered in the orthogonal orientations. Measuring, at each voxel, the effect of water diffusion on the MR signal in a number of directions provides exquisite insight into the local orientation of fibers. Shortly after the first acquisitions of images characterizing the anisotropic diffusion of water molecules in vivo (Moseley et al., 1990; Osment et al., 1990), the Diffusion Tensor (DT) formalism was introduced (Basser et al., 1994). Diffusion Tensor Imaging (DTI) is an analytic means to precisely describe the three-dimensional nature of anisotropy in tissues. The diffusion tensor model encapsulates the average diffusion properties of water molecules (inside a typical 1–3mm sized voxel), as the 3×3 covariance matrix of a Gaussian distribution. DTI has been extensively used to study a wide range of neurological disorders such as cerebro-vascular diseases, multiple sclerosis, Alzheimer’s and Parkinson’s disease, schizophrenia and brain tumors. It has also been very useful for studying brain development, effects of aging, and the structure of anatomical and functional networks in brain substructures, such as the thalamus, striatum and various fiber bundles. Tractography and connectivity mapping techniques are at the core of many of these studies.

However, between one-third to two-thirds of imaging voxels in the human brain’s white matter are thought to contain multiple fiber bundle crossings (Behrens et al., 2007), in which case the Diffusion Tensor model breaks down. High Angular Resolution Diffusion Imaging (HARDI) techniques (Tuch, 2002) such as Diffusion Spectrum Imaging (DSI) (Callaghan et al., 1988; Wiegell et al., 2000) or Q-Ball Imaging (QBI) (Tuch, 2004), have therefore been proposed to overcome the limitations of the diffusion tensor model and recover fiber crossing information. With QBI, model-free mathematical approaches can be developed to reconstruct the angular profile of the diffusion displacement probability density function (PDF) of water molecules, known as the diffusion orientation distribution function (ODF). The underlying fiber distribution (fODF) can also be estimated, which is fundamental for tractography.

The estimation of diffusion tensors or diffusion/fiber ODFs is challenging given the complexity of the diffusion MRI data, the mathematical tools used to describe them and the computational tools used to process them. Here we address several different theoretical and computational issues that arise in processing diffusion MRI. Ultimately, the goal is to recover the underlying fiber distribution so the white matter architecture can be inferred with tractography and segmentation methods. Therefore, we guide the reader from local diffusion model reconstruction in Section 2, to tractography algorithms in Section 3, and finally, to fiber bundle segmentation methods, from DTI or QBI, in Section 4.

2. Local Diffusion Models

To measure the molecular motion in the direction of a given diffusion gradient g, the Stejskal-Tanner imaging sequence (Stejskal and Tanner, 1965) is commonly used and relates the MR signal attenuation S(q, τ) to the statistical properties of the net displacement vector R by

S(q,τ)=S0∫ℝ3p(R∣τ)e−2πiqTRdR=F[p(R∣τ)], (1)

where S0 is a reference signal acquired with no diffusion gradient, τ is the molecular diffusion time, q = γδg/2π is the displacement reciprocal vector (with γ the gyromagnetic ratio of water protons and δ the duration of the diffusion gradients), and p(R|τ) is the ensemble average propagator (EAP). S(q, τ) is thus expressed as the 3D Fourier transform ℱ of the EAP. This function is ultimately the function we are looking to reconstruct in diffusion MRI. Intuitively, one has to sample the diffusion PDF along many q vectors to be able to reconstruct the diffusion PDF. The space of all possible 3D q vectors is called q-space. This is the idea behind q-space imaging (Callaghan, 1991). If the diffusion process is assumed to be Gaussian, the Stjeskal-Tanner equation (1) boils down to

S(g,b)=S0e−bgTDg, (2)

where g is the unit vector q/|q|, b is the so-called b-value given by τ|q|2 and D is the 3 dimensional diffusion tensor. In this case, p takes the simple form

p(R∣τ)=1(4πτ)3∣D∣e−RTD−1R4τ.

In contrast, HARDI acquisition schemes are model-free. They do not make any assumption about the form of the EAP but rather sample q-space along as many directions and q-magnitudes as possible, to reconstruct p as accurately as possible. Typically, there are two strategies used in HARDI: 1) sampling of the whole q-space 3D Cartesian grid and estimation of the EAP by inverse Fourier transformation (DSI), 2) single shell spherical sampling and estimation of fiber distributions from the diffusion/fiber ODF (QBI), Persistent Angular Structure (Jansons and Alexander, 2003) or Diffusion Orientation Transform (Ozarslan et al., 2006). Recent work (Khachaturian et al., 2007) has proposed to improve the accuracy of QBI by fusing information from multiple q-shells. In this section, we focus on the estimation of diffusion tensors and diffusion/fiber ODFs from diffusion MRI.

2.1. Estimation of Diffusion Tensor Images

Several authors have already studied the properties of the non-linear space of diffusion tensors (i.e., symmetric positive-definite matrices) in the context of diffusion tensor processing (Basser and Pajevic, 2003; Fletcher and Joshi, 2004; Batchelor et al., 2005; Pennec et al., 2006; Lenglet et al., 2006b; Fillard et al., 2007; Schwartzman, 2006). Here we briefly describe the theoretical tools necessary for understanding the tensor estimation procedure introduced in this subsection.

2.1.1. Computational Framework

We consider the family of 3D Gaussian distributions with 0-mean as the 6-dimensional parameter space of variances and covariances. We identify it with S+(3), the set of 3 × 3 symmetric positive-definite matrices. A Riemannian metric can be introduced for S+(3) in terms of the Fisher information matrix (Burbea and Rao, 1982), and it has the form gij=g(Ei,Ej)=〈Ei,Ej〉∑=12tr(∑−1Ei∑−1Ej),∀∑∈S+(3), where {Ei} i, j = 1, …, 6 denotes the basis of the tangent space TΣS+(3) = SΣ(3) at Σ ∈ S+(3). The geodesic distance Inline graphic g induced by this metric was investigated by several authors and the original theorem can be found in Atkinson and Mitchell (1981). We recall that the distance between two elements Σ1 and Σ2 of S+(3) is the length of the minimizing geodesic curve between these two points. Calvo and Oller derived an explicit solution of the geodesic equations for general multivariate Gaussian distributions (Calvo and Oller, 1991). The geodesic starting at Σ1 in the direction V ∈ TΣ1S+(3), the tangent space of S+(3) at Σ1, is

∑(t)=∑11/2et∑1−1/2V∑1−1/2∑11/2,∀t∈[0,1], (3)

and the geodesic distance between the two endpoints Σ1 and Σ2 is

Dg(∑1,∑2)=12trace(log2(∑1−1/2∑2∑1−1/2)).

Using these concepts, we can now describe a diffusion tensor estimation algorithm that naturally enforces their positive-definiteness (Lenglet et al., 2006b).

2.1.2. Tensor Estimation

The estimation of a diffusion tensor image relies on the Stejskal-Tanner equation (2). At least N =30 (and not just the theoretical minimum of N =6) diffusion gradients g are typically necessary to robustly estimate the apparent diffusion coefficient, fractional anisotropy and tensor orientations (Jones, 2004). The classical technique for tensor estimation relies on a least-squares procedure where equation (2) is rewritten as a linear system, which can be solved efficiently. However, it has the disadvantage of potentially producing tensors with negative eigenvalues, which are physically impossible. We therefore seek to minimize the following objective function at each voxel of the acquired volume:

E(S0,S(g1,b),…,S(gN,b))=∑i=1Nψ(1bln(S(gi,b)S0)+giTDgi), (4)

where ψ: ℝ ↦ ℝ is a real-valued function which reduces the effect of outliers by replacing the classical least-squares residual by a function such as the Cauchy, Fair, Huber or Tukey M-estimators. We can minimize this energy by an intrinsic gradient descent procedure that naturally evolves on S+(3). The gradient of ℰ is (Lenglet et al., 2006b)

∇E=∑i=1Nψ′(1bln(S(gi,b)S0)+giTDgi)Dgi(DgiT),

where we recall that gi is known and given by the diffusion gradient direction, S is a diffusion weighted image and D is the unknown diffusion tensor. ∇ℰ can be used as the velocity V in the geodesics equation (3) to minimize ℰ while remaining on the manifold of interest S+(3). In Table 1, we compare the performance of a least-squares estimation procedure with the gradient descent technique outlined above and further detailed in Lenglet et al. (2006b). 2500 tensors were generated to span a wide range of configurations and an artificial set of 12 diffusion-weighted images was created from these tensors by using the Stejskal-Tanner equation. Gaussian noise was added and tensors were re-estimated with the two methods. The gradient descent technique clearly outperforms the least-squares approach as it is able to avoid non-positive definite tensors and produces results much closer to the ground-truth data. In practice, although this method is more time-consuming, it has proved to be very useful to avoid degenerate tensors in areas such as the genu or splenium of the corpus callosum. As recently proposed, the Log-Euclidean framework (Fillard et al., 2007) may also be used to estimate diffusion tensors robustly and efficiently. We now move to the estimation of higher-order models of diffusion properties.

Table 1.

Least-squares (LS) and gradient descent estimation: Number of non-positive definite (NPD) tensors (col. 1), and the difference between estimated and ground-truth tensor fields (col. 2–5)

NPD tensors Mean Variance Minimum Maximum
Least Squares 225 12.9 40.4 0.00199 227
Gradient Descent 0 0.658 5.78 0.00409 143
Open in a new tab

2.2. Analytical Reconstruction of the Diffusion ODF

In contrast with DTI, QBI (Tuch, 2004) is a model-independent method to estimate the diffusion ODF, which contains the full angular information of the diffusion PDF and is defined using spherical coordinates as Ψ(θ,φ)=∫0∞p(r,θ,φ)dr, where (θ, φ) obey the physics convention (θ ∈ [0, π], φ ∈ [0, 2π]). A smoothed version of the diffusion ODF can be directly reconstructed from a single-shell HARDI acquisition using the Funk-Radon transform (FRT) (Tuch, 2004). Intuitively, the FRT value at a given spherical point is the great circle integral of the signal on the sphere defined by the plane through the origin perpendicular to the point of evaluation. The original QBI has a numerical solution (Tuch, 2004). More recent methods (Anderson, 2005; Hess et al., 2006; Descoteaux et al., 2007a) have introduced an analytical spherical harmonic reconstruction solution that is faster, more robust to noise and does not require as many gradient-encoding directions. To develop the analytical solution to QBI, the HARDI signal first needs to be represented using spherical harmonics (SH) and then, the FRT can be solved analytically using the SH basis. Letting Yℓm denote the SH of order ℓ and degree m (m = −ℓ, …, ℓ), a modified SH basis that is real and symmetric is defined. For even order ℓ, a single index j in terms of ℓ and m is used such that j(ℓ, m) = (ℓ2 + ℓ + 2)/2 + m. The modified basis is given by

Yj={2Re(Yℓ∣m∣),ifm<0,Yℓm,ifm=0,2(−1)m+1Im(Yℓm),ifm>0, (5)

where Re(Yℓm) and Im(Yℓm) represent the real and imaginary parts of Yℓm respectively. The basis is designed tobe symmetric, real and orthonormal. It is then possible to obtain an analytical diffusion ODF estimate, Ψ, with

Ψ(θ,φ)=∑j=1L2πPℓ(j)(0)cj︸cj′Yj(θ,φ), (6)

where L = (ℓ + 1)(ℓ + 2)/2 is the number of elements in the spherical harmonic basis, cj are the SH coefficients describing the input HARDI signal, Pℓ(j) is the Legendre polynomial of order ℓ(j)2 and cj′ are the coefficients describing the ODF Ψ. Here, the cj coefficients are estimated with the solution presented in Descoteaux et al. (2006) with a Laplace-Beltrami regularization of the SH coefficients cj to obtain a more robust ODF estimation. The detailed implementation of the Laplace-Beltrami regularization and a comparison with other state-of-the art methods (Anderson, 2005; Hess et al., 2006) are presented in Descoteaux et al. (2006, 2007a).

2.3. Analytical Reconstruction of the Fiber ODF

The relation between the measured diffusion ODF and the underlying fiber distribution, the fiber ODF, is still an important open question in the field (Tuch, 2002; Perrin et al., 2005), the answer to which depends on the physics of diffusion, cell membrane permeability, the free diffusion coefficients, axonal packing, the distribution of axonal diameters, the degree of myelination in the underlying fiber bundles, and other parameters. The diffusion ODF is thus a blurred version of the “true” fiber ODF. Because of this blurring effect, the extracted maxima of the diffusion ODF are often used for fiber tractography. An alternative is to use spherical deconvolution (SD) methods that provide an estimate of the fiber ODF (Tournier et al., 2004; Jian and Vemuri, 2007), also called the fiber orientation density (FOD). These techniques have better angular resolution than QBI and produce sharper fiber ODF profiles than the q-ball diffusion ODF. Smaller fiber compartments with smaller volume fractions may be visible in the fiber ODF but not in the diffusion ODF. SD and fiber ODF estimation are currently topics of active research. Here, we use a simple linear transformation of the analytical QBI solution. A schematic view of the spherical deconvolution method is shown in Fig. 1.

Figure 1.

Figure 1

Open in a new tab

Left: Convolution between the diffusion ODF kernel, R, and true fiber ODF produces a smooth diffusion ODF estimate, Ψ. Right: The Funk-Radon transform of the HARDI signal, S, produces a smooth diffusion ODF, Ψ, which is transformed into a sharper fiber ODF estimate, F, by the deconvolution.

The fiber ODF is reconstructed in three steps:

  1. The regularized diffusion ODF coeffcients cj′ are reconstructed using equation (6) of the previous section, cj′=2πPℓ(j)(0)cj/S0.

  2. The single fiber diffusion ODF, R, used as a deconvolution kernel, is estimated from the real data. We assume an axially symmetric diffusion tensor model with eigenvalues (e2, e2, e1) and e1 ≫ e2 for the underlying single fiber diffusion model (Tournier et al., 2004). The values of e1 and e2 are estimated from 300 voxels with highest FA value in our real dataset, as these voxels can each be assumed to contain a single fiber population. The single fiber diffusion ODF kernel has an analytical expression (Descoteaux et al., 2007b) and is given by R(t)=(1−αt2)−1/28πbe1e2, where α = (1 −e2/e1), b is the b-value of the real dataset and t ∈ [−1, 1] is the variable that represents the dot product between the direction of the fiber and the point of evaluation (θ, φ) on the sphere.

  3. The SH coefficients of the fiber ODF, fj, are then obtained by a simple linear transformation, fj=cj′/rℓ(j), with rℓ(j)=2π∫−11R(t)Pℓ(j)(t)dt, which can be solved analytically by taking the power expansion of Pℓ(j)(t) and integrating rℓ(j) term by term. As for the analytical diffusion ODF solution, the spherical deconvolution is obtained with the Funk-Hecke theorem (Descoteaux et al., 2007a). Therefore, the fiber ODF, expressed in terms of the HARDI signal, is

    fj=8πbe1e2Pℓ(j)(0)S0Aℓ(α)cj,withAℓ(α)=∫−11(1−αt2)−1/2Pℓ(t)dt (7)

The final fiber ODF can be reconstructed for any (θ, φ) as F(θ,φ)=∑j=1RfjYj(θ,φ). F provides a valid choice for the fiber ODF (Descoteaux et al., 2007b), in close agreement with the SD method Tournier et al. (2004).

Diffusion tensors and ODFs are at the heart of the white matter connectivity and complexity analysis methods of the next sections.

3. White Matter Connectivity Analysis

3.1. Identification of Sub-voxel Fiber Bundle Configurations

Despite the many advantages of HARDI reconstructions over the diffusion tensor model, they can still be ambiguous and difficult to interpret in the presence of complex sub-voxel fiber tract configurations (Le Bihan et al., 2006; Parker and Alexander, 2005) and thus confound fiber tracking algorithms. Different fiber geometries can yield similar diffusion/fiber ODFs, but require different decisions in tractography. To illustrate this, consider the two types of sub-voxel fiber structures depicted in Fig. 2. Both the single curving fiber tract (left) and the fanning fiber tract (right) are likely to result in an almost indistinguishable ODF with a single broad peak oriented in the vertical direction (middle panel in left and right subfigures). This is due to the presence of a relatively wide array of fiber tangent directions within one voxel, which results in a broad fiber ODF profile. Although they yield a similar ODF, because of a similiar sub-voxel fiber tangent distribution, each of these structures requires a different decision from a fiber tracking algorithm. For a curving fiber bundle, only one path should be recovered, with a local tangent given by the medial direction of the broad maximum (left). But for a fanning configuration, multiple paths should be followed when propagating in one direction (polarity vector represented by a green arrow in the right subfigure), and only one direction should be followed when propagating in the other (blue arrow). This illustrates the importance of recovering the polarity of the fanning in addition to its extent. Hence, in order to take the appropriate action, it is crucial for tracking methods to be able to differentiate between these types of sub-voxel configurations. The distinction between a fanning fiber tract and a single fiber tract with high curvature is known to be particularly challenging (Le Bihan et al., 2006; Parker and Alexander, 2005). It can be addressed by relating fiber ODF data to the underlying white matter fiber tracts, modeled as 3D curves (Savadjiev et al., 2008). This approach is based on the 3D curve inference method (Savadjiev et al., 2006) described below, which infers the curves that best describe the underlying white matter fibers in each voxel. By integrating information over a local neighborhood, this method can resolve ambiguities that cannot be clarified by only considering individual ODFs.

Figure 2.

Figure 2

Open in a new tab

Sub-voxel fiber configurations that can cause ambiguous fiber ODF shapes. The red curves denote fibers in a specific voxel, and the yellow glyph shows the reconstructed ODF. The arrows show the directions of a typical fiber path traversing this voxel; configurations differ despite having the same ODF.

3.1.1. 3D Curve Inference

3D curve inference is a differential geometric method for inferring of helical arcs as osculating approximations to arbitrary 3D curves, based on a support measure defined over a local neighborhood. In the context of diffusion MRI orientation data, this enables the local curvature and torsion of white matter fibers to be estimated. Distinct sub-voxel fiber bundle configurations that share the same tangents (orientations) at a particular voxel can be distinguished from one another. As input, the algorithm requires a discretized regular (typically rectangular) 3D lattice, with a fiber ODF defined at each lattice location (voxel). Each of these ODFs is then sampled along several orientations. A notion of co-helicity can be formally defined (Savadjiev et al., 2006) and relates individual orientations at distinct voxels through a geometrical constraint. In particular, the conditions under which three orientations defined at three distinct locations in space can be tangent to a helix are determined. Based on these conditions, a measure of the support that a given orientation (at a given voxel) receives from co-helical configurations of neighboring orientations (at neighboring voxels) is calculated. This measure is weighted by the ODF value along these orientations at the voxels of interest. We then discretize the parameter space describing helical curves. A best-fit helix is determined for each orientation, based on the support obtained from the neighborhood. 3D curve inference can be used to perform ODF regularization (Savadjiev et al., 2006). ODF values in orientations that are not aligned with the inferred curves are discarded, whereas those that are aligned with the inferred curves are supported. The inferred best-fit helices can be used to disambiguate curving vs. fanning subvoxel fiber configurations (Savadjiev et al., 2008), as described next.

3.1.2. Labeling of Ambiguous Sub-voxel Fiber Tract Configurations

To motivate the approach, consider Fig. 3, which shows a schematic of the inferred curves in the case of a fanning fiber tract and of a single curving fiber tract. For simplicity, only a 2D case is illustrated, but the technique is applicable to any 3D ODF dataset. In the general 3D case, the inferred curves will be helical, i.e., they will have both curvature and torsion. Local helix approximations to fibers are constructed by searching for co-helical triplets of fiber ODF maximum orientations in a local neighborhood. A co-helical triplet is interpolated to a helix which is used as a local approximation to the (arbitrary) 3D curve that represents an underlying white matter fiber tract. Thus, a given fiber ODF orientation presents evidence for an underlying curve (helix) if it is the central element of a co-helical triplet in a spherical neighborhood centered on that voxel and if its parameterization inferred through 3D curve inference agrees with the parameters (curvature, torsion) of the cohelical triplet. As an example, the three sampling orientations corresponding to the three blue maxima in Fig. 3 (right) form a co-helical triplet of orientations. Similarly, in Fig. 3 (left), the groups of red, green and blue maxima all form co-helical triplets with the black maximum, which is common to all three groups. One or more such helices can pass through a given voxel, as one helix is assigned to each ODF orientation in that voxel. For example, three such helices pass through the central ODF in Fig. 3 (left), associated with the red, green and blue orientations, respectively. The number and the configuration of these helices are used to label the voxel as belonging to a fanning, crossing, or single fiber tract configuration, as described in Table 2. This labeling uses two types of information: 1) ODF shape information at each voxel. (The number of local maxima can be used to distinguish crossings from the other two cases. The crossing case is included for completeness.) 2) A geometric model inferred from a neighborhood of voxels. The helices inferred by the 3D curve inference algorithm are used to distinguish fanning configurations from single, possibly curving, fiber tracts. We emphasize that the inferred helical curves are local approximations to more global 3D curves, obtained independently at each voxel. Thus the local helical approximations will typically differ from one location (voxel) to the next. In summary, the approach uses evidence from inferred helix curves and their local configurations to disambiguate the three cases outlined in Table 2. Since helices are parametric curves, and since they are represented by co-helical triplets of tangents, it is straightforward to check the number and local configuration of helices that pass through any given voxel. These ideas are developed into an algorithm, described in pseudocode in Savadjiev et al. (2008), where implementation details are also discussed. As we show in the next section, this labeling information turns out to be very important when performing fiber tractography.

Figure 3.

Figure 3

Open in a new tab

Co-helical triplets obtained by the 3D curve inference algorithm. In the fanning case (left), multiple ODFs on the fanning side can have their maxima (red, green, blue) cohelical with different orientations in the central (ambiguous) ODF and with the same maximum (black) of an ODF on the merging side. The inferred curves are shown with dashed lines. In the single curving tract case (right), only one set of co-helical triplets of ODF maxima (blue) exists. Only one curve is inferred.

Table 2.

Labels obtained from local helix curves configuration and fiber ODF maxima.

Label ODF Maxima Curves Configuration
Single fiber tract 1 One helix.
Fanning 1 Two or more distinct and diverging helices.
Crossing ≥ 2 Number of helices = number of ODF maxima with sufficient angular separation.
Open in a new tab

3.2. Diffusion MRI Tractography

3.2.1. Overview

The previous sections have discussed how to robustly extract information about fiber orientations, at the voxel scale, using diffusion MRI. One application of such techniques is to infer global connectivity in the central nervous system. Fiber tractography is used to integrate the voxel-scale fiber orientations in order to create maps of connectivity between distant areas. The delineation of these pathways is useful in determining whether specific areas of the brain are connected, the course of these connections, and how these connections change in diseases. The extracted pathways can be used as regions of interest in which to investigate other scalar parameters, such as fractional anisotropy or measures extracted from other types of MRI contrast. The pathways may also be used for parcellation of given brain regions based on differences in connectivity to and from them (Behrens et al., 2003; Frey et al., 2006). Numerous fiber tractography algorithms exist, but different integration and interpolation schemes, and varying step sizes (Lori et al., 2000), and seeding protocols (Campbell et al., 2005), can greatly influence the streamline propagation results. Integration approaches include Fiber Assignment by Continuous Tracking (FACT) (Mori et al., 1999), Euler-Lagrange (Conturo et al., 1999), and Runge-Kutta (Basser et al., 2000). Each technique has benefits, such as robustness to high curvature (Lazar and Alexander, 2003), as well as drawbacks. Front evolution approaches, including level-set techniques, have also been investigated for tractography (Parker et al., 2002; Campbell et al., 2005; Tournier et al., 2003; Jackowski et al., 2005; Prados et al., 2006). There are clear theoretical benefits of incorporating information about complex sub-voxel fiber geometry in tractography. We have investigated the improvements in vivo and in phantoms. Fig. 4 shows an example tractography comparison in the human brain. Fig. 4a is the result of using the diffusion tensor model and propagating along the principal eigenvector. Fig. 4b is the result of using the q-ball model-free reconstruction method (Tuch, 2004), propagating along the fiber ODF maxima. Using DTI, generally only the most medial projections of the corpus callosum are seen. Q-Ball Imaging is more capable of picking up the projections that cross through the cortico-spinal tract and superior longitudinal fasciculus. The tractography shown in Fig. 4c is performed by exploiting the labeling of fanning and curving fibers described in Section 3.1. Using the labeling information generally results in subtle improvements in the sensitivity of the tractography, especially for fiber bundles fanning toward the cortex. This can be a great benefit, e.g., in assessing connectivity between cortical areas, such as co-activated areas in functional studies.

Figure 4.

Figure 4

Open in a new tab

Effect of incorporating information on complex sub-voxel fiber geometry in tractography. Fiber tracking result in human brain using (a) DTI, (b) QBI, and (c) the fiber geometry information described in Section 3.1. Tracking results are displayed as isosurfaces encompassing all voxels connected to the seed ROI (green)

3.2.2. Details and Comparison of Tractography Methods

For the tracking results shown in Fig. 4b/c, we proceeded as follows. In Fig. 4b, the deterministic method described in Savadjiev et al. (2008) was used. No interpolation was performed. Streamlines were propagated using FACT integration along the ODF maxima. The tracking was stopped if the fractional anisotropy (FA) was less than 0.1, the mean diffusivity was greater than 10−6 mm2/ms, or the angular difference in the orientation of the tract from one voxel to the next was greater than 80° (these parameters also apply to Fig. 4a). An alternate solution for deterministic tracking using fiber ODFs was proposed to take into account multiple maxima at each step (Descoteaux et al., 2007b). This is illustrated in Fig. 5a. Similar initialization and parameters to those of Savadjiev et al. (2008) were used and trilinear interpolation was performed to obtain diffusion ODF, fiber ODF and DT at subvoxel precision.

Figure 5.

Figure 5

Open in a new tab

Deterministic and probabilistic tracking of the projections of the corpus callosum to Broca’s area. (a) Deterministic Split-Streamline tractography, (b) Probabilistic tractogram, (c) Probabilistic fibers colored by their end point projections. In subfigures (b) & (c), colors indicate probability of connection.

In Fig. 4c, except where fanning was indicated, the direction of propagation was also given by the ODF maximum. However, when fanning was indicated (i.e., the dot product of the incoming direction with the fanning polarity vector was positive), the tracking algorithm followed all fiber directions from the fiber ODF, thus exploiting the rich information provided by the labeling of ambiguous sub-voxel fiber tact configurations of Section 3.1. Tracking was initiated in all voxels in a small region of interest in the corpus callosum. At voxels labeled as single curves, the propagation direction was given by the fiber ODF maximum closest to the incoming direction. At voxels labeled fanning fibers, the direction of propagation depended on whether the incoming direction was in the direction of the fanning, or in the direction of the merge. For the former, all directions in the extent of the fanning area, given by the fiber ODF distribution, were followed; for the latter, the fiber ODF maximum only was followed. For voxels with multiple fiber ODF maxima, the maximum closest to the incoming direction was followed. The fanning was accomplished by running the entire tracking process iteratively and randomly selecting a direction at each iteration. 10000 iterations were used. For all starting ROI voxels, the tracking was initiated on a 3×3×3 grid of start points in order to facilitate branching. Streamlines were propagated using FACT integration. Tracts that erroneously turned down the cortical-spinal tract were excluded. A connectivity profile of all voxels reached by the tracking was saved and an example result is presented in Fig.4c.

Finally, the rich information of fiber ODFs can be exploited in a probabilistic way (Descoteaux et al., 2007b). A random walk method (Koch et al., 2002) was extended as follows: A large number of particles is typically started from a seed point. The particles randomly propagate according to our local fiber ODF estimate, F, and the number of times a voxel is reached by the path of a particle is counted. This yields higher transitional probabilities along the main fiber directions. For each elementary transition of a particle, the probability of a movement from the seed point x to the target point y in direction uvxy is computed as the product of the local fiber ODFs in direction uvxy, i.e. P(x → y) = F(uvxy)x · F (uvxy)y, where P (x → y) is the probability for a transition from point x to point y and F(uvxy)x is the fiber ODF at point x in direction xy (by symmetry, direction xy and yx are the same). This method is illustrated in Fig. 5b/c. A tractogram (i.e., the 3D distribution of connected voxels to the seed voxel) of a voxel of the corpus callosum and sample fiber tracts included in the probability map are presented.

3.2.3. Sensitivity and Error Analysis

Fiber tractography is rarely fully automatic. It is prone to false positive and false negative results, and the precise tracking protocol has a large effect on the end results. As described above, we can initiate tracking on a grid within each start voxel. This approach facilitates branching and improves the sensitivity of the technique. We have found that the density of seeding in each start voxel impacts the tracking results: the point at which the result converges depends on the system under investigation. For any given pathway, the regions of interest selected to delineate the pathway will impact the result. The user can also choose between initiating tracking only in these regions of interest, or initiating tracking everywhere and retaining those tracts that pass through the ROIs (the “brute force” approach (Conturo et al., 1999)). Fig. 6a shows the result using brute force seeding for the same tract-delineating ROI used in Fig. 4. Brute force tractography generally reduces false negatives and seed point dependence of the results, and enables fanning to some extent, as can be seen by comparison to Fig. 4c. It also guarantees reversible tractography results: the connectivity between point x and point y should be the same as the connectivity between point y and point x. Tractography is usually supervised, in that exclusion masks are used to remove false positives. At the resolution of typical diffusion imaging protocols, it is probable that inferred fiber trajectories will jump from one pathway to another, so some prior knowledge of the anatomy is necessary. One common artifact, when reconstructing the corpus callosum for instance, is to jump onto the cortico-spinal tract (CST), as shown in Fig. 6b. Masks were used to remove the CST from the other tracking results in Fig. 4.

Figure 6.

Figure 6

Open in a new tab

Error and uncertainty in tractography. (a) Brute force seeding strategy, using identical parameters to Fig. 4b. (b) False positive tracking results in the absence of exclusion masks. (c) Maximum intensity projection of probabilistic connectivity values. All experiments used the same tract-delineating ROI in the corpus callosum.

As with any measurement, uncertainty in the tractography result should be quantified. Uncertainty in tractography arises from uncertainty in the directions of propagation in all of the voxels that constitute a tract. The uncertainty should reflect our confidence that there is a fiber in a given direction, and the confidence in the direction itself. Uncertainty in the direction arises from imaging noise and from limitations of the chosen reconstruction technique. In the case of DTI, it is relatively easy to acquire multiple datasets in order to estimate the standard deviation σθ of the distribution for the angle of propagation. For high angular resolution, in which the minimal data acquisition is much larger, statistical techniques such as bootstrap methods (Haroon and Parker, 2007; Berman et al., 2008) and Markov Chain Monte Carlo in a Bayesian framework (Fonteijn et al., 2007; Behrens et al., 2007) are being explored. Here, we illustrate this concept of uncertainty in probabilistic q-ball tractography in Fig. 6c, using the finite angular resolution of the acquisition to determine σθ (Frey et al., 2006). The confidence in the direction of propagation is given by a truncated Gaussian profile with maximum at the maximum of the diffusion ODF and standard deviation σθ. The connectivity map is shown as a maximum intensity projection. The connectivity value of a voxel to a given ROI is given by the lowest confidence value of all tract segments along the tract between the voxel the ROI (Parker et al., 2002; Campbell et al., 2005). An advantage of this “weakest link” approach over counting the number of times a voxel is reached by a random tracking process (Behrens et al., 2007) is that if a voxel is reached many times but by different routes, it still has low probability of connection to the ROI. It is also more amenable to the brute force tracking approach. Here, tracking was performed using the brute force approach, with 1000 starts per voxel in order to randomly sample the cones of uncertainty around each maximum.

In the human brain, there is no gold standard tractography result with which to evaluate tractography algorithms. Anatomy varies from brain to brain, and our understanding of human neuroanatomy is still incomplete. There is therefore a need for phantoms for evaluation of tractography. There has been considerable work on synthetic phantoms, e.g., (Lazar and Alexander, 2003; Tournier et al., 2002; Perrin et al., 2005; Close et al., 2008). Physical phantoms are useful for evaluating MRI sequences and evaluating post-processing techniques in the presence of imaging artifacts, such as eddy-current induced distortions, and noise. These include biological phantoms (Boujraf et al., 2001) and phantoms made from textiles (Watanabe et al., 2006). We have created a physical phantom made from excised rat spinal cords (Campbell et al., 2005) (see Fig. 7). While simple, this phantom provides a gold standard tractography result for quantitative comparisons between algorithms, seeding strategies, and other tracking parameters. A physical phantom should have structures that restrict diffusion in a timeframe compatible with an MRI experiment, and diffusion properties and complex fiber geometries comparable to those in vivo. It is also desirable to match the MR relaxation properties of human tissue. Fig. 7 illustrates the use of the phantom to compare tracking approaches. Both diffusion tensor and q-ball reconstruction were performed on the same dataset. Fig. 7a and 7b show tracking results using the diffusion tensor model and q-ball imaging, respectively, using the same tract-delineating ROI on the curved tract and brute force seeding. The tensor model fails to capture the crossing fibers necessary to reconstruct the whole tract. Fig. 7c illustrates ODF tracking without brute force seeding: the differences between using brute force seeding and not are often evident near curves. This example also illustrates that connections to regions more distant from the ROI are sparser.

Figure 7.

Figure 7

Open in a new tab

Fiber tracking result in phantom using (a) diffusion tensor and brute force seeding, (b) q-ball and brute force seeding, (c) q-ball with seeding only in the tract-delineating ROI, located in the bottom curved part of the tract. A transparent surface indicating the gold standard cord segmentation is shown for reference.

4. White Matter Segmentation and Complexity Analysis

Clustering methods for diffusion MRI have been recently introduced and provide a complementary point of view to the analysis of the white matter architecture. They typically rely on some metric between diffusion tensors or ODFs and allow us to identify various fiber bundles or regions of the white matter with different diffusion pro-files. While many techniques have been proposed to classify the gray matter, white matter and cerebrospinal fluid from T1-weighted MR images, the literature addressing the clustering of white matter and sub-cortical structures from diffusion MRI is fairly recent. In this context, two main features identify each element of a diffusion image: the position on the image and the diffusion characteristics. To perform effective clustering, the contribution of these two features must be carefully exploited. Quantifying the similarity between the diffusion features (tensors, ODFs) is still a subject of current investigations. In the following, we denote the position in the diffusion image by xi and Di stands for the diffusion characteristic (either the diffusion tensor or some representation of the ODF).

4.1. Methods Based on DTI

The first approach that used DTI to elucidate structure in the brain by means of clustering was designed to identify the different nuclei of the thalamus (Wiegell et al., 2003). It uses a k-means algorithm. The spatial metric is the Mahalanobis distance with respect to each cluster and the feature metric is the Frobenius norm of the difference between tensors. The choice of this last metric is crucial and discussed in the following, where we focus on fiber bundle segmentation.

One of the very first approaches to fiber bundle segmentation (Zhukov et al., 2003) was able to cluster white matter structures by only using the fractional anisotropy as the diffusion characteristic, in a surface evolution framework (which is well-suited for controlling the shape and smoothness of the resulting clusters). A 3D surface S is represented by the zero level set of a 4-dimensional function φ, S = {x ∈ ℝ3: φ (x, t) = 0}, and φ is evolved according to the differential equation ∂ φ(x, t)/∂t = −F (x)||∇ φ (x, t)||, where F is a scalar-valued function which drives the evolution of φ, and implicitely deforms the surface S along its normals. F is usually made of two terms F = Fc +βFs. Fc quantifies characteristics of the regions to segment and Fs drives the smoothness of the surface; β is a user-selected weight. Fs and Fc can be respectively chosen as the mean curvature of the surface S and an edge detector function ℝ ↦ [0, 1] applied to a smoothed FA map (Zhukov et al., 2003). To generalize this approach and take advantage of the full tensor information, the function Fc was adapted to work on a generalized structure tensor for diffusion tensor fields (Feddern et al., 2003). This approach allowed some improvements over the previous work and was among the first to focus on the definition of an adequate metric between diffusion tensors. Later, a statistical level-set segmentation method was introduced (Wang and Vemuri, 2004). In this method, Fc is based on a regional description of the inner and outer compartments:

Fc=−Df(D,D¯in)+Df(D,D¯out)=−∥D−D¯in∣∣F2+∥D−D¯out∣∣F2, (8)

where D̄in and D̄out are the Fréchet means of side and outside the surface S. The Fréchet mean D̄ of a set of N tensors {Di}i=1,…, N is analytically or iteratively computed as the minimizer of the tensors’ variance, depending on the choice of metric. Such a regional approach allows the tensors in the inner and outer regions to vary in a piecewise constant manner, contrary to approaches which only search for sharp variations of the FA, or other anisotropy maps (Zhukov et al., 2003). Thus, the algorithm presented in Wang and Vemuri (2004) is capable of detecting fiber bundles where the tensor’s shape changes smoothly. However, because of the use of the Euclidean distance Inline graphic f between tensors (Frobenius norm), the Fréchet mean D̄in/out is not guaranteed to be positive-definite, thus generating artifacts and incorrect segmentations in regions where tensors’ variation is large. To overcome this problem, several authors have studied the influence of the metric. Considering the diffusion tensor as the covariance matrix of a zero-mean Gaussian distribution, the symmetrized Kullback-Leibler divergence or J-Divergence was introduced (Wang and Vemuri, 2005):

Dj(D1,D2)=12trace(D1−1D2+D1D2−1)−6,

and used to extend previous work (Wang and Vemuri, 2004). It is a natural metric between probability distributions, which turns out to have a closed form expression in the Gaussian case. It also has a closed form expression for the mean tensor. Next, two different distances between tensors were introduced in a similar level set formulation (Jonasson et al., 2005). This was later extended to prevent overlapping when propagating multiple surfaces for the segmentation, for instance, of the thalamic nuclei (Jonasson et al., 2007b). The first distance, called Integrated Similarity, compares the diffusion properties from two different voxels. It is expressed as Dis(D1,D2)=14π∫S2min(d1(r)d2(r),d2(r)d1(r))dr, where d1(r) is the diffusion coefficient in direction r for the tensor D1. This metric compares diffusion coefficients over all possible directions and is very sensitive to small differences between the shapes of the tensors. It has, however, a high computational cost. Another metric is used to calculate the empirical mean of a set of tensors. It measures the overlap between two tensors and is defined as Do(D1,D2)=trace(D1D2). It has also been proposed to use a Riemannian metric derived from the Fisher information matrix in an extended statistical framework (Lenglet et al., 2006a). This metric yields a geodesic distance on the manifold of zero-mean Gaussian distributions S+(3), as presented in Section 2.1.1. This distance can also be expressed as Dg(D1,D2)=12∑i=13log2(λi) where the scalars λi are the eigenvalues of the matrix D1−1D2D1−1. The authors also demonstrated how to approximate a Gaussian distribution on S+(3) and to exploit this information in the segmentation procedure. The practical differences of using the three different distances, Inline graphic f/j/g, in the extended statistical surface evolution frame-work of Lenglet et al. (2006a) are illustrated in Fig. 8. Recently, a distance with similar properties to those of Inline graphic g was introduced (Arsigny et al., 2006). This Log-Euclidean distance

Figure 8.

Figure 8

Open in a new tab

DTI statistical segmentation using different distances: Inline graphicf,Inline graphicj, Inline graphicg Differences are observed in the splenium of the corpus callosum.

Dle(D1,D2)=trace((log(D1)−log(D2))2)

has the advantage to be simple to implement and fast to compute. However, it has not been extensively studied for segmentation tasks yet, aside from its use in two recent papers (Weldeselassie and Hamarneh, 2007; Malcolm et al., 2007). In both works, the image is segmented into two parts by minimizing an energy functional similar to the one of Wang and Vemuri (2005). A non-parametric approach relying on the Log-Euclidean distance and a Markov random field framework was also recently described (Awate et al., 2007). Finally, a graph-theoretical approach, known as N-Cuts (Shi and Malik, 2000), was used (Ziyan et al., 2006). This graph partitioning technique is based on the link between the second smallest eigenvector of the Laplacian matrix of a graph and optimal partitions. The nodes of the graph are the voxel xi of the image and the weights of the edges between those nodes are obtained from similarities between neighboring tensors. The similarity between tensors can be any of the previously described distances or restricted to the directional information of the principal eigenvectors. The outline of the procedure proposed in Ziyan et al. (2006) is as follows: First, a matrix Ws is built to encode local similarity between tensors. It is only non-zero for neighboring voxels:

{Ws}ij={exp(−D2(Di,Dj)σ2),if∥xi−xj∥≤10,otherwise. (9)

σ is a chosen scale parameter. Next, local similarities are propagated to a full affinity matrix W by converting Ws into a one-step transition probability matrix whose rows and columns sum to one. Markovian relaxation (Tishby and Slonim, 2000) is used to generate the n-step transition probability matrix. Finally, this matrix is recursively partitioned in two clusters by using eigenvalue decomposition and thresholding the eigenvector with the second smallest eigenvalue. This produces a hierarchical clustering. The number of recursions is ultimately chosen to obtain the desired number of clusters. The main limitation of this algorithm is the need for a uniform sampling of the (xi, Di) at the end of the Markovian relaxation. An extension of this approach to ODFs was recently proposed (Wassermann et al., 2008), as we will describe in the next section. As we are now going to discuss, level set, Markov random fields and graph-theoretic segmentation frameworks have also recently been extended to HARDI datasets.

4.2. Methods based on HARDI

As described in Section 2, the diffusion tensor model cannot describe complex white matter fiber configurations, and HARDI techniques like QBI were introduced to overcome this issue. It is thus natural to exploit this information to improve white matter segmentation results. The 5D space defined by the location of the ODFs on the acquisition grid and their orientational information can be used (Hagmann et al., 2006; Jonasson et al., 2007a). These segmentation procedures are respectively implemented with a hidden Markov random field or level set framework. Mixtures of von Mises-Fisher distributions were proposed to model the ODFs (McGraw et al., 2006) and segmentation was also performed with a hidden Markov random field scheme in this work. It is possible to use a spherical harmonics (SH) decomposition of the ODF at each voxel xi. The diffusion characteristic Di is then replaced by a vector of SH coefficients. Although the most appropriate metric between ODFs is an open area of research, the L2 norm can be used and efficiently computed. Two other clustering techniques have been proposed that take advantage of the SH representation. First, in Descoteaux and Deriche (2008), we generalized the level set algorithm presented in Wang and Vemuri (2004); Lenglet et al. (2006a) to the HARDI case. Equation (8) is modified such that the distance is the L2 norm of the difference of SH coefficients. As can be seen in the two images on the left of Fig. 9, the projections of the corticospinal tract to the cortex can only be successfully segmented by using HARDI data. Second, in Wassermann et al. (2008), we extended the work presented in Ziyan et al. (2006) in three ways: 1) The matrix presented in equation (9) was rewritten to use the L2 norm between the SH coefficients of the ODFs. 2) The matrix obtained after Markovian relaxation is normalized in order to relax the hypothesis on the sampling of the (xi, Di). 3) This matrix is finally used as an input for the diffusion maps method (Lafon and Lee, 2006). This method produces a mapping, where each element (xi, Di) is represented as a point in Euclidean space, and an estimation of the number of clusters existing in the set. A k-means algorithm is used to automatically find the clusters. The two images on the right of Fig. 9 show that segmentation results are coherent with white matter anatomical knowledge.

Figure 9.

Figure 9

Open in a new tab

Left: Comparison between DTI and HARDI level set segmentation methods. Right: Diffusion Maps clustering of HARDI data

In the following section, we show how HARDI clustering techniques can be extended to quantify the geometrical complexity of segmented fiber bundles. Beyond the assignment of labels, it is possible to provide a scalar measure that correlates with the expected variations of the configurations of the white matter fiber tracts.

4.3. Quantification of the Non-uniform Complexity of the White Matter

In this section, we describe how the white matter microstructure complexity (or dimensionality) may be studied. Here, the complexity is understood as meaning the minimum number of parameters needed to represent the diffusion MRI data (tensors, ODFs) in an underlying sub-manifold of ℝm (m ≥ 6 depends on the order of the SH approximation of ODFs). Regions with or without fiber crossings clearly belong to manifolds with different complexity, and we use stratifications (Haro et al., 2008b) (i.e., the union of manifolds with different dimensions and densities) to quantify the local complexity of DTI and HARDI datasets and relate it to known features of neuroanatomy (Haro et al., 2008a).

A geometric and probabilistic method was recently proposed to estimate the local dimension and density of point clouds in ℝm (Levina and Bickel, 2005). It was then extended by modeling high-dimensional sample points as a process of translated Poisson mixture, with regularizing restrictions (Haro et al., 2008b). Noise is naturally handled and it is possible to identify the underlying manifolds with different dimensions and densities. The outline of the method is as follows: If we sample an m-dimensional manifold with T points, the proportion of points falling into a ball around xt is kT≈ρ(xt)V(m)Rk(xt)m (Levina and Bickel, 2005). The point cloud of interest, embedded in high dimension D, is X = {xt ∈ ℝD; t = 1,…, T}, k is the number of points inside the ball, ρ(xt) is the local sampling density at point xt, V(m) is the volume of the unit sphere in ℝm, and Rk(xt) is the Euclidean distance from xt to its k-th nearest neighbor (kNN). The inhomogeneous process N(R, xt), which counts the number of points falling into a D-dimensional sphere B(R, xt) of radius R centered at xt, can be approximated by a Poisson process with rate λ(R, xt) = ρ(xt)V (m)mRm−1. The local intrinsic dimension and density estimators at each point xt are obtained from the Maximum Likelihood (ML) estimator based on a Poisson distribution with this rate.

Noise usually contaminates the point cloud, so the observed point process (i.e., the ODFs) is not a sampling of a low-dimensional manifold but rather a perturbation of this sample process. This is modeled with a translated Poisson process (Snyder and Miller, 1991), which translates an underlying (unobservable) point process into an output (observable) point process according to a conditional probability density f called the transition density. If each point is translated independently and no deletion or insertion occur, any translated Poisson process with an integrable intensity function λ on the input space X is also a Poisson process with intensity μ(z) = ∫Xf(z|x)λ(x)dx, ∀ z ∈Z, on the output space. The intensity of our observable process λ(r, xt) is parameterized by the Euclidean distance between points so the density f(s|r) transforms a distance r in the input space to a distance s in the observable space. The intensity of the Poisson process in the output space is μ(s,xt)=∫0R′f(s∣r)ρ(xt)V(m)mrm−1dr, where R′ > R since points originally at a distance greater than R from xt can be placed within a distance less than R after the translation process.

Maximizing the likelihood of the new translated Poisson process, we obtain a nonlinear recursive expression for the local dimension m(xt) at point xt, which is difficult to solve. We approximate it by an easier to compute closed expression, with explicit bounds on the approximation (see Haro et al. (2008b) for details),

m(xt)≈[1k−1∑i=1k−1∫0R′f(Ri∣r)logRkrdr∫0R′f(Ri∣r)dr]−1. (10)

The associated density estimator is ρ(xt)=(k−1)/(V(m(xt))Rkm(xt)). However, these estimators are local. We propose to compute an ML on the whole point cloud simultaneously (not just for each point independently) by using a mixture of translated Poisson distributions which accommodates noise and different classes (characterized by their own dimension and density) (Haro et al., 2008b). This Translated Poisson Mixture Model (TPMM) is solved with an Expectation Maximization algorithm, which leads to explicit estimations of each cluster dimensionality and density, as well as a soft clustering according to these parameters; see Haro et al. (2008b) for details.

We recently applied this technique to HARDI datasets. ODFs were estimated with the techniques described in Section 2. We examined the complexity of the raw Diffusion Weighted Images (points in ℝ30) as well as that of 4th and 6th order ODFs (ℝ15 and ℝ28 respectively), and corresponding (sharpened) fiber ODFs (Fig. 10, top). Clusterings from 4th and 6th order ODFs are almost identical, as 30 gradients may be insufficient to fit a detailed 6th order model. Clusterings obtained from the ODFs are clearly better than those from the raw HARDI data and we can readily distinguish (Fig. 10, top) the gray matter in green, complex white matter in purple (e.g., forceps minor/major, anterior/posterior corona radiata or superior longitudinal fasciculus), anisotropic white matter in light blue (e.g., genu/splenium of the corpus callosum or internal capsule), and highly anisotropic white matter in blue (e.g., genu of the corpus callosum, cortico-spinal tract). Using 4th order fiber ODFs had little effect. 6th order fiber ODFs decreased the clustering accuracy, perhaps by enhancing high-frequency noise in the higher-order model. We also compared our complexity/dimension estimates to the known complexity of white matter configurations, in the genu of the corpus callosum and forceps minor. Callosal fibers are tightly packed at the interhemispheric plane, but diverge and mingle with other fiber bundles as they progress toward the frontal lobes. Our method identifies and quantifies this increase in complexity. The dimension and density of the four submanifolds both increase (Fig. 10, bottom) as fibers leave the highly anisotropic genu region.

Figure 10.

Figure 10

Open in a new tab

Top: Influence of the input model on the labeling of 2 axial slices. Bottom: Increasing complexity in the forceps minor

5. Conclusion

Diffusion MRI and variations such as DTI or QBI have opened up a landscape of extremely exciting discoveries for medicine and neuroscience. The development of powerful analysis tools for these modalities has occupied the medical image analysis community for about a decade now and has already resulted in fundamental advances in research on various neurological disorders such as stroke, cancer or neurodegenerative diseases. Here we have briefly reviewed state-of-the-art mathematical models and computational techniques for processing diffusion MRI data. We showed how to efficiently estimate local diffusion models such as the diffusion tensor or diffusion/fiber orientation distribution functions. We then described a framework to identify of sub-voxel fiber bundle configurations (crossing, fanning, etc.) from QBI data. Along with an overview of current approaches for white matter tractography, we showed how this framework can be used to extend QBI deterministic tractography. We finally introduced a set of techniques based on surface evolution or graph-theoretic methods for clustering white matter structures in DTI or QBI. We also demonstrated how machine learning techniques can be applied to such datasets to quantify the non-uniform complexity of the cerebral white matter.

Acknowledgments

This work was partially supported by ONR, NGA, NSF, NIH (grants P41 RR008079, P30 NS057091, R01 EB007813, R01 HD050735, and P41 RR013642), DARPA, the PAI Procope, the Diffusion MRI INRIA ARC, NSERC Canada, FQRNT Québec and the Spanish Ministerio de Ciencia e Innovación, under program Juan de la Cierva and project TEC2007-66858/TCM.

Footnotes

2

ℓ (j) is the order associated with the jth element of the SH basis, i.e., for j = 1, 2, 3, 4, 5, 6, 7, … ℓ(j) = 0, 2, 2, 2, 2, 2, 4, …

Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

References

  1. Anderson A. Measurements of fiber orientation distributions using High Angular Resolution Diffusion Imaging. Magn Reson Med. 2005;54:1194–1206. doi: 10.1002/mrm.20667. [DOI] [PubMed] [Google Scholar]
  2. Arsigny V, Fillard P, Pennec X, Ayache N. Log-Euclidean metrics for fast and simple calculus on diffusion tensors. Magnetic Resonance in Medicine. 2006;56:411–421. doi: 10.1002/mrm.20965. [DOI] [PubMed] [Google Scholar]
  3. Atkinson C, Mitchell A. Rao’s distance measure. Sankhya: The Indian Journal of Stats. 1981;43(A):345–365. [Google Scholar]
  4. Awate S, Hui Z, Gee J. A fuzzy, nonparametric segmentation framework for DTI and MRI analysis: With applications to DTI-tract extraction. IEEE Transactions on Medical Imaging. 2007;26 (11):1525–1536. doi: 10.1109/TMI.2007.907301. [DOI] [PubMed] [Google Scholar]
  5. Basser P, Mattiello J, Bihan DL. Estimation of the effective self-diffusion tensor from the NMR spin echo. Journal of Magnetic Resonance B. 1994;(103):247–254. doi: 10.1006/jmrb.1994.1037. [DOI] [PubMed] [Google Scholar]
  6. Basser P, Pajevic S. A normal distribution for tensor-valued random variables: Applications to diffusion tensor MRI. IEEE Transactions on Medical Imaging. 2003;22 (7):785–794. doi: 10.1109/TMI.2003.815059. [DOI] [PubMed] [Google Scholar]
  7. Basser P, Pajevic S, Pierpaoli C, Duda J, Aldroubi A. In vivo fiber tractography using DT-MRI data. Magn Reson Med. 2000;44:625–632. doi: 10.1002/1522-2594(200010)44:4<625::aid-mrm17>3.0.co;2-o. [DOI] [PubMed] [Google Scholar]
  8. Batchelor P, Moakher M, Atkinson D, Calamante F, Connelly A. A rigorous framework for diffusion tensor calculus. Magnetic Resonance in Medicine. 2005;53 (1):221–225. doi: 10.1002/mrm.20334. [DOI] [PubMed] [Google Scholar]
  9. Behrens T, Johansen-Berg H, Woolrich M, Smith S, Wheeler-Kingshott C, Boulby P, Barker G, Sillery E, Sheehan K, Ciccarelli O, Thompson A, Bardy J, Matthews P. Non-invasive mapping of connections between human thalamus and cortex using diffusion imaging. Nature Neuroscience. 2003;6 (7):750–757. doi: 10.1038/nn1075. [DOI] [PubMed] [Google Scholar]
  10. Behrens TEJ, Johansen-Berg H, Jbabdi S, Rushworth MFS, Woolrich MW. Probabilistic diffusion tractography with multiple fibre orientations. what can we gain? NeuroImage. 2007;34 (1):144–155. doi: 10.1016/j.neuroimage.2006.09.018. [DOI] [PMC free article] [PubMed] [Google Scholar]
  11. Berman J, Chung S, Mukherjee P, Hess C, Han E, Henry R. Probabilistic streamline q-ball tractography using the residual bootstrap. NeuroImage. 2008;39 (1):215–222. doi: 10.1016/j.neuroimage.2007.08.021. [DOI] [PubMed] [Google Scholar]
  12. Boujraf S, Luypaert R, Eisendrath H, Osteaux M. Echo planar magnetic resonance imaging of anisotropic diffusion in asparagus stems. MAGMA. 2001;13 (2):82–90. doi: 10.1007/BF02668156. [DOI] [PubMed] [Google Scholar]
  13. Burbea J, Rao C. Entropy differential metric, distance and divergence measures in probability spaces: A unified approach. Journal of Multivariate Analysis. 1982;12:575–596. [Google Scholar]
  14. Callaghan PT. Principles of nuclear magnetic resonance microscopy. Oxford University Press; Oxford: 1991. [Google Scholar]
  15. Callaghan PT, Eccles CD, Xia Y. Rapid Communication: NMR microscopy of dynamic displacements: k-space and q-space imaging. Journal of Physics E Scientific Instruments. 1988;21:820–822. [Google Scholar]
  16. Calvo M, Oller J. An explicit solution of information geodesic equations for the multivariate normal model. Statistics and Decisions. 1991;9:119–138. [Google Scholar]
  17. Campbell J, Siddiqi K, Rymar V, Sadikot A, Pike G. Flow-based fiber tracking with diffusion tensor and q-ball data: Validation and comparison to principal diffusion direction techniques. NeuroImage. 2005;27 (4):725–736. doi: 10.1016/j.neuroimage.2005.05.014. [DOI] [PubMed] [Google Scholar]
  18. Close T, Tournier J, Calamante F, Johnson L, Mareels I, Con-nelly A. Software tool to generate complex structures for validation of fibre tracking. Proc. of the International Society for Magnetic Resonance in Medicine; 2008. p. 431. [Google Scholar]
  19. Conturo T, Lori N, Cull T, Akbudak E, Snyder A, Shimony J, McKinstry R, Burton H, Raichle M. Tracking neuronal fiber pathways in the living human brain. Proc of the National Academy of Sciences. 1999;96:10422–10427. doi: 10.1073/pnas.96.18.10422. [DOI] [PMC free article] [PubMed] [Google Scholar]
  20. Descoteaux M, Angelino E, Fitzgibbons S, Deriche R. Apparent diffusion coefficients from High Angular Resolution Diffusion Imaging: Estimation and applications. Magnetic Resonance in Medicine. 2006;56:395–410. doi: 10.1002/mrm.20948. [DOI] [PubMed] [Google Scholar]
  21. Descoteaux M, Angelino E, Fitzgibbons S, Deriche R. Regularized, fast, and robust analytical q-ball imaging. Magnetic Resonance in Medicine. 2007a;58 (3):497–510. doi: 10.1002/mrm.21277. [DOI] [PubMed] [Google Scholar]
  22. Descoteaux M, Deriche R. High angular resolution diffusion MRI segmentation using region-based statistical surface evolution. Journal of Mathetical Imaging in Vision 2008 [Google Scholar]
  23. Descoteaux M, Deriche R, Anwander A. Deterministic and probabilistic q-ball tractography: from diffusion to sharp fiber distributions. Tech Rep. 2007b;6273 INRIA Sophia Antipolis. [Google Scholar]
  24. Feddern C, Weickert J, Burgeth B. Level-set methods for tensor-valued images. Proc. of the Second IEEE Workshop on Geometric and Level Set Methods in Computer Vision; 2003. pp. 65–72. [Google Scholar]
  25. Fillard P, Arsigny V, Pennec X, Ayache N. Clinical DT-MRI estimation, smoothing and fiber tracking with Log-Euclidean metrics. IEEE Transactions on Medical Imaging. 2007;26 (11):1472–1482. doi: 10.1109/TMI.2007.899173. [DOI] [PubMed] [Google Scholar]
  26. Fletcher P, Joshi S. Principal geodesic analysis on symmetric spaces: Statistics of diffusion tensors. Proc. Computer Vision Approaches to Medical Image Analysis; Prague. 2004. pp. 87–98. [Google Scholar]
  27. Fonteijn HJ, Verstraten FJ, Norris D. Probabilistic inference on q-ball imaging data. IEEE Trans in Med Imaging. 2007;26 (11):1515–1524. doi: 10.1109/TMI.2007.907297. [DOI] [PubMed] [Google Scholar]
  28. Frey S, Campbell J, Pike G, Petrides M. Dissociation of areas 44 and 45 (broca’s area) using diffusion tensor fiber tractog-raphy. Society for Neuroscience Meeting; 2006. p. 159. [Google Scholar]
  29. Hagmann P, Jonasson L, Deffieux T, Meuli R, Thiran JP, Wedeen VJ. Fibertract segmentation in position orientation space from High Angular Resolution Diffusion MRI. Neu-roImage. 2006;32:665–675. doi: 10.1016/j.neuroimage.2006.02.043. [DOI] [PubMed] [Google Scholar]
  30. Haro G, Lenglet C, Sapiro G, Thompson P. On the nonuniform complexity of brain connectivity. Proc. IEEE Intl. Symposium on Biomedical Imaging: From Nano to Macro; Paris. 2008a. pp. 887–890. [Google Scholar]
  31. Haro G, Randall G, Sapiro G. Translated poisson mixture model for stratification learning. International J of Computer Vision 2008b [Google Scholar]
  32. Haroon HA, Parker GJ. Using variants of the wild bootstrap to quantify uncertainty in fibre orientations from q-ball analysis. Proc. of the Organization for Human Brain Mapping annual meeting; 2007. p. 273. [Google Scholar]
  33. Hess C, Mukherjee P, Han E, Xu D, Vigneron D. Q-ball reconstruction of multimodal fiber orientations using the spherical harmonic basis. Magnetic Resonance in Medicine. 2006;56:104–117. doi: 10.1002/mrm.20931. [DOI] [PubMed] [Google Scholar]
  34. Jackowski M, Kao C, Qiu M, Constable R, Staib L. White matter tractography by anisotropic wavefront evolution and diffusion tensor imaging. Medical Image Analysis. 2005;9 (5):427–440. doi: 10.1016/j.media.2005.05.008. [DOI] [PMC free article] [PubMed] [Google Scholar]
  35. Jansons K, Alexander D. Persistent angular structure: new insights from diffusion magnetic resonance imaging data. Inverse Problems. 2003;19:1031–1046. [Google Scholar]
  36. Jian B, Vemuri B. A unified computational framework for deconvolution to reconstruct multiple fibers from diffusion weighted MRI. IEEE Transactions on Medical Imaging. 2007;26 (11):1464–1471. doi: 10.1109/TMI.2007.907552. [DOI] [PMC free article] [PubMed] [Google Scholar]
  37. Jonasson L, Bresson X, Hagmann P, Cuisenaire O, Meuli R, Thiran JP. White matter fiber tract segmentation in DT-MRI using geometric flows. Medical Image Analysis. 2005;9:223–236. doi: 10.1016/j.media.2004.07.004. [DOI] [PubMed] [Google Scholar]
  38. Jonasson L, Bresson X, Thiran JP, Wedeen VJ, Hagmann P. Representing diffusion MRI in 5-D simplifies regularization and segmentation of white matter tracts. IEEE Trans on Med Imaging. 2007a;26 (11):1547–1554. doi: 10.1109/TMI.2007.899168. [DOI] [PubMed] [Google Scholar]
  39. Jonasson L, Hagmann P, Pollo C, Bresson X, Wilson C, Meuli R, Thiran JP. A level set method for segmentation of the thalamus and its nuclei in DT-MRI. Signal Processing. 2007b;87:309–321. [Google Scholar]
  40. Jones D. The effect of gradient sampling schemes on measures derived from diffusion tensor MRI: A Monte Carlo study. Magnetic Resonance in Medicine. 2004;51:807–815. doi: 10.1002/mrm.20033. [DOI] [PubMed] [Google Scholar]
  41. Khachaturian M, Wisco J, Tuch D. Boosting the sampling efficiency of q-ball imaging using multiple wavevector fusion. Magnetic Resonance in Medicine. 2007;57:289–296. doi: 10.1002/mrm.21090. [DOI] [PubMed] [Google Scholar]
  42. Koch M, Norris D, Hund-Georgiadis M. An investigation of functional and anatomical connectivity using magnetic resonance imaging. NeuroImage. 2002;16:241–250. doi: 10.1006/nimg.2001.1052. [DOI] [PubMed] [Google Scholar]
  43. Lafon S, Lee A. Diffusion maps and coarse-graining: A uni-fied framework for dimensionality reduction, graph partitioning, and data set parameterization. IEEE Pattern Analysis and Machine Intelligence. 2006;28 (9):1393–1403. doi: 10.1109/TPAMI.2006.184. [DOI] [PubMed] [Google Scholar]
  44. Lazar M, Alexander A. An error analysis of white matter tractography methods: synthetic diffusion tensor field simulations. NeuroImage. 2003;20 (2):1140–1153. doi: 10.1016/S1053-8119(03)00277-5. [DOI] [PubMed] [Google Scholar]
  45. Le Bihan D, Breton E. Imagerie de diffusion in vivo par résonance magnétique nucléaire. CR de l’Académie des Sciences. 1985;301:1109–1112. [Google Scholar]
  46. Le Bihan D, Poupon C, Amadon A, Lethimonnier F. Artifacts and pitfalls in diffusion MRI. J Magn Reson Imaging. 2006;24(3):478–88. doi: 10.1002/jmri.20683. [DOI] [PubMed] [Google Scholar]
  47. Lenglet C, Rousson M, Deriche R. DTI segmentation by statistical surface evolution. IEEE Transactions on Medical Imaging. 2006a June;25 (6):685–700. doi: 10.1109/tmi.2006.873299. [DOI] [PubMed] [Google Scholar]
  48. Lenglet C, Rousson M, Deriche R, Faugeras O. Statistics on the manifold of multivariate normal distributions: Theory and application to diffusion tensor MRI processing. J Mathematical Imaging and Vision. 2006b;25 (3):423–444. [Google Scholar]
  49. Levina E, Bickel P. NIPS. Vol. 17. Vancouver: 2005. Maximum likelihood estimation of intrinsic dimension. [Google Scholar]
  50. Lori N, Akbudak E, Snyder A, Shimony J, Conturo T. Diffusion tensor tracking of human neuronal fiber bundles: Simulation of effects of noise, voxel size and data interpolation. Proc. International Society for Magnetic Resonance in Medicine; 2000. p. 775. [Google Scholar]
  51. Malcolm J, Rathi Y, Tannenbaum A. A graph cut approach to image segmentation in tensor space. IEEE Conference on Computer Vision and Pattern Recognition; 2007. pp. 1–8. [DOI] [PMC free article] [PubMed] [Google Scholar]
  52. McGraw T, Vemuri B, Yezierski R, Mareci T. Segmentation of High Angular Resolution Diffusion MRI modeled as a field of von mises-fisher mixtures. Proc. European Conf. on Computer Vision; 2006. pp. 463–475. [Google Scholar]
  53. Merboldt K, Hanicke W, Frahm J. Self-diffusion NMR imaging using stimulated echoes. J Magn Reson. 1985;64:479–486. [Google Scholar]
  54. Mori S, Crain B, Chacko V, van Zijl P. Three dimensional tracking of axonal projections in the brain by magnetic resonance imaging. Annals of Neurology. 1999;45:265–269. doi: 10.1002/1531-8249(199902)45:2<265::aid-ana21>3.0.co;2-3. [DOI] [PubMed] [Google Scholar]
  55. Moseley M, Cohen Y, Mintorovitch J, Kucharczyk J, Tsuruda J, Weinstein P, Norman D. Diffusion-weighted MR imaging of anisotropic water diffusion in cat central nervous system. Radiology. 1990;176 (2):439–445. doi: 10.1148/radiology.176.2.2367658. [DOI] [PubMed] [Google Scholar]
  56. Osment P, Packer K, Taylor M, Attard JJ, Carpenter TA, Hall LD, Doran SJ, Herrod NJ. NMR imaging of fluids in porous solids. Phil Trans Roy Soc. 1990;333:441–452. [Google Scholar]
  57. Ozarslan E, Shepherd T, Vemuri B, Blackband S, Mareci T. Resolution of complex tissue microarchitecture using the diffusion orientation transform (DOT) NeuroImage. 2006;31 (3):1086–1103. doi: 10.1016/j.neuroimage.2006.01.024. [DOI] [PubMed] [Google Scholar]
  58. Parker G, Alexander D. Probabilistic anatomical connectivity derived from the microscopic persistent angular structure of cerebral tissue. Phil Trans R Soc B. 2005;360:893–902. doi: 10.1098/rstb.2005.1639. [DOI] [PMC free article] [PubMed] [Google Scholar]
  59. Parker GJM, Wheeler-Kingshott CAM, Barker GJ. Estimating distributed anatomical connectivity using fast marching methods and diffusion tensor imaging. IEEE Transactions Medical Imaging. 2002;21:505–512. doi: 10.1109/TMI.2002.1009386. [DOI] [PubMed] [Google Scholar]
  60. Pennec X, Fillard P, Ayache N. A Riemannian framework for tensor computing. International Journal of Computer Vision. 2006;66 (1):41–66. [Google Scholar]
  61. Perrin M, Poupon C, Cointepas Y, Rieul B, Golestani N, Pallier C, Riviere D, Constantinesco A, Bihan DL, Mangin J-F. Fiber tracking in q-ball fields using regularized particle trajectories. Information Processing in Medical Imaging; 2005. pp. 52–63. [DOI] [PubMed] [Google Scholar]
  62. Prados E, Lenglet C, Pons J, Wotawa N, Deriche R, Faugeras O, Soatto S. Control theory and fast marching methods for brain connectivity mapping. Proc. IEEE Conf. CVPR; 2006. pp. 1076–1083. [Google Scholar]
  63. Savadjiev P, Campbell JSW, Descoteaux M, Deriche R, Pike GB, Siddiqi K. Labeling of ambiguous sub-voxel fibre bundle configurations in high angular resolution diffusion MRI. NeuroImage. 2008;41 (1):58–68. doi: 10.1016/j.neuroimage.2008.01.028. [DOI] [PubMed] [Google Scholar]
  64. Savadjiev P, Campbell JSW, Pike GB, Siddiqi K. 3D curve inference for diffusion MRI regularization and fibre tractog-raphy. Medical Image Analysis. 2006;10 (5):799–813. doi: 10.1016/j.media.2006.06.009. [DOI] [PubMed] [Google Scholar]
  65. Schwartzman A. PhD thesis. Stanford University; 2006. Random ellipsoids and false discovery rates; statistics for diffusion tensor imaging data. [Google Scholar]
  66. Shi J, Malik J. Normalized cuts and image segmentation. IEEE Pattern Analysis and Machine Intelligence. 2000;22 [Google Scholar]
  67. Snyder DL, Miller MI. Random Point Processes in Time and Space. Springer-Verlag; 1991. [Google Scholar]
  68. Stejskal E, Tanner J. Spin diffusion measurements: spin echoes in the presence of a time-dependent field gradient. Journal of Chemical Physics. 1965;42:288–292. [Google Scholar]
  69. Taylor D, Bushell M. The spatial mapping of translational diffusion coefficients by the NMR imaging technique. Physics in Medicine and Biology. 1985;30 (4):345–349. doi: 10.1088/0031-9155/30/4/009. [DOI] [PubMed] [Google Scholar]
  70. Tishby N, Slonim N. NIPS. Vol. 14 2000. Data clustering by markovian relaxation and the information bottleneck method. [Google Scholar]
  71. Tournier J, Calamante F, Gadian D, Connelly A. Diffusion-weighted magnetic resonance imaging fibre tracking using a front evolution algorithm. NeuroImage. 2003;20 (1):276–288. doi: 10.1016/s1053-8119(03)00236-2. [DOI] [PubMed] [Google Scholar]
  72. Tournier J, Calamante F, Gadian D, Connelly A. Direct estimation of the fiber orientation density function from diffusion-weighted MRI data using spherical deconvolution. NeuroImage. 2004;23:1176–1185. doi: 10.1016/j.neuroimage.2004.07.037. [DOI] [PubMed] [Google Scholar]
  73. Tournier J, Calamante F, King M, Gadian D, Connelly A. Limitations and requirements of diffusion tensor fiber tracking: an assessment using simulations. Magnetic Resonance in Medecine. 2002;47 (4):701–708. doi: 10.1002/mrm.10116. [DOI] [PubMed] [Google Scholar]
  74. Tuch D. PhD thesis. Harvard University and Massachusetts Institute of Technology; 2002. Diffusion mri of complex tissue structure. [Google Scholar]
  75. Tuch D. Q-ball imaging. Magn Reson Med. 2004;52(6):1358–1372. doi: 10.1002/mrm.20279. [DOI] [PubMed] [Google Scholar]
  76. Wang Z, Vemuri B. Proc ECCV. 2004. Tensor field segmentation using region based active contour model; pp. 304–315. [Google Scholar]
  77. Wang Z, Vemuri B. DTI segmentation using an information theoretic tensor dissimilarity measure. IEEE Trans Med Imag. 2005;24(10):1267–1277. doi: 10.1109/TMI.2005.854516. [DOI] [PubMed] [Google Scholar]
  78. Wassermann D, Descoteaux M, Deriche R. Diffusion maps clustering for magnetic resonance Q-Ball imaging segmentation. International Journal of Biomedical Imaging. 2008;2008(526906) doi: 10.1155/2008/526906. [DOI] [PMC free article] [PubMed] [Google Scholar]
  79. Watanabe M, Aoki S, Masutani Y, Abe O, Hayashi N, Ma-sumoto T, Mori H, Kabasawa H, Ohtomo K. Flexible ex vivo phantoms for validation of diffusion tensor tractography on a clinical scanner. Radiat Med. 2006;24(9):605–609. doi: 10.1007/s11604-006-0076-4. [DOI] [PubMed] [Google Scholar]
  80. Weldeselassie Y, Hamarneh G. DT-MRI segmentation using graph cuts. SPIE Medical Imaging. 2007;6512:65121. [Google Scholar]
  81. Wiegell M, Larsson H, Wedeen VJ. Fiber crossing in human brain depicted with diffusion tensor MR imaging. Radiology. 2000;217:897–903. doi: 10.1148/radiology.217.3.r00nv43897. [DOI] [PubMed] [Google Scholar]
  82. Wiegell M, Tuch D, Larsson H, Wedeen VJ. Automatic segmentation of thalamic nuclei from diffusion tensor magnetic resonance imaging. NeuroImage. 2003;19 (2):391–401. doi: 10.1016/s1053-8119(03)00044-2. [DOI] [PubMed] [Google Scholar]
  83. Zhukov L, Museth K, Breen D, Whitaker R, Barr A. Level set modeling and segmentation of DT-MRI brain data. Journal of Electronic Imaging. 2003;12 (1):125–133. [Google Scholar]
  84. Ziyan U, Tuch D, Westin C-F. Segmentation of thalamic nuclei from DTI using spectral clustering. Proc. MICCAI; 2006. pp. 807–814. [DOI] [PubMed] [Google Scholar]

RESOURCES