Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Neuroimage. Author manuscript; available in PMC 2010 August 1.
Published in final edited form as:
PMCID: PMC2720421

Active Fibers: Matching Deformable Tract Templates to Diffusion Tensor Images


Reliable quantitative analysis of white matter connectivity in the brain is an open problem in neuroimaging, with common solutions requiring tools for fiber tracking, tractography segmentation and estimation of intersubject correspondence. This paper proposes a novel, template matching approach to the problem. In the proposed method, a deformable fiber-bundle model is aligned directly with the subject tensor field, skipping the fiber tracking step. Furthermore, the use of a common template eliminates the need for tractography segmentation and defines intersubject shape correspondence. The method is validated using phantom DTI data and applications are presented, including automatic fiber-bundle reconstruction and tract-based morphometry.

Keywords: diffusion tensor imaging, DTI tractography, template matching, deformable curve evolution

1 Introduction

Diffusion tensor magnetic resonance imaging (DT-MRI) is a relatively new technique that measures diffusion of water molecules in tissue [1]. This modality has become particularly important in neuroimaging, since it offers insight into the in vivo properties of the white matter (WM) tissue in the brain, and thus can help study connectivity of WM axons in human subjects. The relevance of DT imaging to such studies stems from the fact that, in axons, water molecules diffuse faster along the axonal directions than across their membranes. Thus, modeling water diffusion in each voxel (which is done in DT-MRI using a second-order tensor) provides a certain indication of the local orientation of the neural fibers in that voxel. The problem of tractography consists of inferring the complete shapes of those fibers based on local diffusion measurements [25]. In the simplest setting, these measurements are typically reduced to the principal eigenvectors of the diffusion tensors [2,3].

Still, reliable reconstruction of the neural fibers remains an open problem in neuroimaging. Importantly, accurate inference at the individual axon level is fundamentally infeasible due to the low resolution of today’s DT-MR images. In fact, the very notion of accuracy is ill-posed in this context, since precise validation is typically impossible due to the lack of ground truth data for human subjects. As a result, there is no gold standard in tractography algorithms to date. The most commonly used DTI streamline methods [2,3] suffer from several notorious shortcomings, such as the ambiguity of the diffusion tensor model with several fibers crossing a single voxel (partial volume effects) and early termination of fibers, due to low fractional anisotropy (FA) and curvature thresholds. Furthermore, fiber tractography alone, although instrumental for visualization, has only a limited anatomical and analytical applicability. In order to perform quantitative analysis, it is necessary to apply a post-processing step of tractography segmentation that involves filtering, labeling and grouping of the fibers into higher-level bundles that have a clear anatomical interpretation and are found consistently across subjects [6]. Finally, to conduct statistical studies, one needs a way of establishing shape correspondence between bundles in different subjects, which is an open problem as well.

These challenges motivate the use of algorithms that can leverage prior anatomical knowledge. To this end, in this paper we advocate a template matching, “one shot” approach to the reconstruction and subsequent connectivity analysis of fiber tracts of interest. Indeed, multiple WM tracts in the human brain can be identified consistently across subjects, thereby lending themselves to a template representation. Additionally, template matching methods offer several distinct advantages in geometric object reconstruction, and in medical imaging applications in particular (e.g., in brain mapping [7]). Typically, the template serves as a shape prior, modeling the expected geometric and topological (and, sometimes, semantic or anatomical) properties of the reconstructed object, whereas the variability is captured through the fitting process. Furthermore, template matching facilitates object reconstruction at a superior level of detail, as the template is typically represented at a higher resolution. Finally, the resulting mapping between the subject and the template shape provides a “built-in” way to establish a common intersubject coordinate system—a valuable property for statistical shape analysis and morphometric studies. Despite these advantages, template matching approaches in the context of WM tractography remain largely unexplored to date.

The contributions of this paper can be summarized as follows:

  • We propose a novel method for automatic matching of fiber tracts to diffusion tensor (DT) images. The method is hybrid in the sense that a geometric template is matched directly to the subject tensor field, skipping the fiber tracking step. Additionally, the use of a common template for all subjects facilitates the subsequent analysis, eliminating the problems of tractography segmentation and fiber correspondence.
  • The proposed method can benefit from any available prior knowledge. The use of a shape prior helps avoid some common problems of standard tractography algorithms, e.g., early termination of fibers due partial volume effects, low FA and/or high curvature. Furthermore, tract-specific anatomical constraints, as well as probabilistic priors, can be incorporated into the matching process to improve reliability of the result.
  • We present a deformable fiber evolution model cast in the framework of generalized geometric flows, based on our derivation of a discrete curve Laplacian operator. The basic model is then further extended to deformable fiber-bundles that can readily be used in practical scenarios.
  • We provide a validation of the method based on a digital DTI phantom. Additionally, we present several practical applications on human brain data, i.e., direct reconstruction of fiber-bundles and tract-based morphometry.

2 Related Work

To date, several knowledge-based approaches to tract-specific reconstruction and quantification have been proposed in the literature. Most commonly, interactive methods are employed, where one selects, from a whole-brain tractography result, a subset of fibers passing through user-defined regions of interest (ROIs). Specialized protocols, e.g., [8,9], have been suggested to dissect tracts of interest in this fashion, and a multi-subject WM atlas was constructed using this technique[10]. However, such interactive selection requires extensive expert knowledge of the 3D fiber tract anatomy. To facilitate the process, automated ROI placement has been suggested, based on affine registration of the structural MRI (T1) subject images to the above atlas [9]. Alternatively, automatic clustering methods [11,12] can be applied to label fibers without the use of ROIs, relying on a previously learned atlas instead. In either case, the problem of intersubject tract shape correspondence remains, although it can be somewhat facilitated with the latter methods by using smaller, very specific clusters. Additionally, the results are still affected by tractography artifacts. To remedy this problem, Hua et al. [8] proposed to skip the tractography step altogether and perform the segmentation directly using atlas-aligned tract probability maps. However, they acknowledge that their results are influenced by the atlas-subject registration quality. In fact, the reliance on affine registration of structural MRI for WM correspondence becomes a limitation in all the above methods, as such a registration, even if extended to a deformable solution, can never be guaranteed to correctly align WM tracts. This is due to the fact that fiber orientations can only be inferred from diffusion images. Ziyan et al. [13,14] tackle this problem by recovering the pairwise affine transformations between corresponding fiber-bundles and merging the transformation into a global poly-affine deformation. Alas, the affine transformation model for each bundle is still too restricted to capture intersubject variability and is still sensitive to tractography artifacts. We thus propose a fully deformable alignment approach, where geometric tract templates are registered directly to the subject DT field, as explained in the next section.

3 Algorithm

3.1 Problem Formulation

Let us first consider a simpler case of a single-fiber template, modeled as a curve An external file that holds a picture, illustration, etc.
Object name is nihms96332ig1.jpg: [0,1] → R3. Given a diffusion tensor field D(x), we wish to find the best deformable alignment of the template curve at each point to the tensor orientations, subject to certain smoothness and position constraints. More formally (and ignoring smoothness temporarily), we are looking for a deformation function Φ that satisfies Φ(An external file that holds a picture, illustration, etc.
Object name is nihms96332ig1.jpg) = An external file that holds a picture, illustration, etc.
Object name is nihms96332ig2.jpg, where An external file that holds a picture, illustration, etc.
Object name is nihms96332ig2.jpg is a curve maximizing the alignment functional (a.k.a. anisotropic Dirichlet energy):


where An external file that holds a picture, illustration, etc.
Object name is nihms96332ig2.jpg′(t) denotes the tangent vector to An external file that holds a picture, illustration, etc.
Object name is nihms96332ig2.jpg (t). Intuitively, the integrand measures the local alignment of the tangent at a given point to the tensor orientation at that point, using all the information present in the tensor without the need to recover the eigenvectors individually. With this measure, the total energy is mostly influenced by the anisotropic tensors, whereas the isotropic ones are “indifferent” to the local fiber orientation. As a consequence, this global energy is relatively insensitive to local artifacts such as partial volume effects.

Since the energy ED is unbounded, we must constrain the solution by adding a smoothness term as follows:


where Δ is the Laplace-Beltrami operator for curves (henceforth referred to simply as the Laplacian), and α is the desired smoothness factor. This new term, denoted by EB(An external file that holds a picture, illustration, etc.
Object name is nihms96332ig1.jpg), corresponds to the biharmonic energy of An external file that holds a picture, illustration, etc.
Object name is nihms96332ig1.jpg: while ED(An external file that holds a picture, illustration, etc.
Object name is nihms96332ig1.jpg) measures the alignment of An external file that holds a picture, illustration, etc.
Object name is nihms96332ig1.jpgt to D, the biharmonic energy serves as a smoothness measure of An external file that holds a picture, illustration, etc.
Object name is nihms96332ig1.jpgt.

Note that ED is conceptually similar to the path cost function An external file that holds a picture, illustration, etc.
Object name is nihms96332ig2.jpg′ (t)T D−1An external file that holds a picture, illustration, etc.
Object name is nihms96332ig2.jpg′ (t), where D−1 serves as a Riemannian metric tensor that can be used to measure intrinsic distances between two points in the white matter [4,5,15]. Indeed, we can also see the integral of this cost function along the curve as a misalignment energy (to be minimized), an option that reduces the dependence regularization. The drawback is potential noise amplification due to the tensor inversion and numerical instabilities in the case of very anisotropic tensors. In practice, the two formulations are almost identical in terms of implementation, and the optimal choice between the two may depend on the nature of the data. In either case, however, since diffusion tensors are never perfectly anisotropic in reality, the DT data is preprocessed with a normalized sharpening step, as suggested by Fletcher et al. [15]. Specifically, D^(x)=D(x)13(D(x)D(x)13)α, with α = 4 chosen empirically in our experiments. Throughout this paper, we refer to the sharpened tensor field simply as D.

3.2 Discrete Setup

In practice, the curves are given in a polygonal representation rather than a parametric one, and we will stick to the discrete setting throughout the rest of this paper. Updating the notation to Y = Φ(X) with X = {xi}i=0..N, the integral in (1) becomes a discrete weighted sum:


i is the average value of D over the edge element ei and li = |ei| is the integration weight associated with ei. This choice of using edges (as opposed to vertices) as basic discrete elements is not accidental: for a polyline edge e, we can directly define a discrete counterpart of the tangent vector at any point p on e simply as ee. Note that this definition converges to the continuous case as the sampling density approaches infinity.

To discretize the biharmonic energy, we first need to define a Laplacian operator for discrete curves. One way to do that is to consider first the discrete curvature at a vertex xi [set membership] X which can be defined straightforwardly as


where mi = (li + li−1)/2 is the “discrete mass” of vertex xi (here, simply the arc length distance between the midpoints of ei−1 and ei). In fact, K(X) can be seen as the 1D counterpart of the mean curvature of a surface—this analogy becomes obvious if one regards a curve as a ribbon surface with a fixed infinitesimal width. Thus, we can reuse here a fact known for discrete surfaces: the matrix K is nothing but the discrete Laplacian in its pointwise formulation [16,17] which, in turn, is the direct analog of the continuous Δ. (For a detailed discussion of discrete surface Laplacians see, e.g., [17,18].) As a result, we can write the discrete biharmonic energy of X as EB(X)=KXL22, where ||·||L2 denotes the discrete L2 norm of vector fields in R3, given by XL22=XTMX [19]. Here M denotes the diagonal (or lumped) “mass matrix” associated with X, defined as Mii = mi (recall Eq. 4). Essentially, the mass matrix is needed to account for irregular sampling of the geometry. We thus obtain:


with a further simplification of this equation available in Appendix A.

3.3 Optimization Procedure

The optimization algorithm is essentially an iterative evolution of the template curve X so as to maximize the combined energy Ematch, using a gradient descent approach. This process is known as discrete gradient flow and is given by:


with the derivation of EmatchX found in Appendix B.

Deformation Prior

So far, we only used the template fiber X (with new positions assigned to its endpoints) as the initialization for our optimizer. While such an initialization is helpful to avoid local minima, in reality we also want to ensure that the final curve stays, in some sense, as close as possible to the template. This will minimize the likelihood of fibers straying from their own trajectories due the influence of other nearby tracts. One way to implement this is by adding a fidelity term to the total energy, penalizing stretching and bending of the evolving curve with respect to the original. In practice, adding a bending term combined with endpoint constraints will bound the stretching as well. The task is, therefore, to design an energy that measures bending with respect to a curved reference shape X. Although this can be done by employing a discrete thin shell model [20], such a strategy would have several drawbacks. First, the resulting energy would be highly nonlinear, leading to slower convergence. Perhaps more importantly, we would like to keep the final mixture of different energies as simple as possible, as it can become increasingly delicate to find a proper balance between them. For these reasons, we opt for a different approach as follows. Consider the incremental bending energy of Xt+1 w.r.t. Xt that can be expressed as the total squared change of curvature: EIB(t+dt)=Kt+dtXt+dtKtXtL22. For a sufficiently small dt, we can assume an approximately isometric deformation, i.e., KtKt+dt, leading to


Note that this deformation energy is very similar to the biharmonic energy from Sec. 3.2, with the sole difference that it is defined on a deformation field Ut, as opposed to a shape. Furthermore, it represents a norm in the space of vector fields in R3, and we can now modify the gradient flow in Eq. 6 to penalize instantaneous bending according to this norm. This is done using a Sobolev-like deformation prior [19,21,22] of the form Id+λΔ2 which can be shown to correspond directly to EIB. Here λ is a parameter assigning a relative cost of the bending deformation component. The resulting discrete flow equation then becomes:


where L = MKTMK. This scheme has several key advantages. First, it helps avoid local minima and resist noise, without oversmoothing the final solution. Secondly, the deformation prior does not, in theory, prevent the evolution from reaching the original global optimum 1. The iteration is terminated when convergence becomes sufficiently slow, typically yielding a tradeoff between the objective function and the fidelity criterion. While this is an approximate solution, it performs reasonably well for most our cases, as reported in Sec. 4.1.3.

3.4 Generalization to a Fiber-Bundle Template

So far we have only considered a single-fiber template. In practice, however, we typically seek to reconstruct specific fiber-bundles rather than individual fibers. It is therefore more useful to perform template matching on the bundle level. Thus, we generalize the basic method to fiber-bundle templates in the following simple way. Since a bundle is essentially a set of curves, the naïve approach is to treat each curve independently. However, this may cause a bundle to disintegrate, since there is no “glue” to hold the fibers together as they evolve. We address this problem by extending the deformation prior from the previous section so as to enforce a smooth motion of the entire bundle, such that the neighboring fibers move in a coherent fashion. With the deformation prior approach, adding such a global smoothness criterion becomes straightforward when given a Laplacian operator for the underlying shape, i.e., the full bundle, as opposed to a set of “autonomous” curves. Indeed, the matching energy from Eq. 2 readily extends to the bundle case simply by summing up the individual fiber energies. Thus, all we need to do is generalize the deformation prior by defining a bundle Laplacian. However, this task involves a challenge. Since the Laplacian operator essentially encodes local spatial relationships, even its simplest discrete-case formulation (a.k.a. the graph Laplacian matrix, see e.g. [23]) typically relies on vertex connectivity information, i.e., a finite element mesh is required. A common way to obtain a good-quality mesh is by computing a 3D Delaunay triangulation (i.e., tetrahedralization) of the bundle vertex set. However, we are interested in a mesh whose edges preserve the given curve connectivities, which becomes a constrained Delaunay problem that is much harder to solve in 3D. The brute-force alternative is to establish vertex connectivity simply by finding k nearest neighbors for each for each point. However, this will often produce local meshes of poor numerical quality. Consequently, we establish inter-fiber connectivity using a hybrid of the two methods, in a way similar to the approach proposed in the context of point-based surfaces [24]. A local 3D Delaunay triangulation Ti is computed (using CGAL2) around each fiber vertex vi, based on a candidate set of nearest points from the neighboring fibers and its existing edge neighbors of vi —see Figure 2 for an illustration. The edges incident on vi in Ti are then used to establish the local connectivity which is integrated into a global graph Laplacian matrix — an approximation we retain due to its simplicity and robust numerical properties. Specifically, we define the bundle Laplacian L as follows: Li j = [sharp]ei j and Lii = −Σvj[set membership]N(vi) Li j, where [sharp]ei j counts the edges connecting vi and vj in both directions. Note that L is symmetric, which is consistent with the properties of continuous Laplace-Beltrami operators3. Finally, the updated deformation prior from Eq. 8 becomes:

Fig. 2
Establishing Inter-Fiber Connectivity Within a Bundle

Here, the second term of the right-hand side still acts on the bundle fibers separately, involving solely intra-fiber connectivity. The new term, however, relies on interfiber connectivity, with μ being the cost of diverging fiber motion. Thus, informally speaking, the fibers can still diverge to reach a global optimum, but only if there is no “cheaper” way. Otherwise, they will tend to move in concert, which also makes the evolution process more robust to noise. Note that the dimension of the system is now given by the total number of vertices in the entire bundle.

4 Applications

4.1 Fiber Tract Reconstruction

4.1.1 Phantom Data Experiments

We first evaluate the single-fiber template approach using a phantom DT image with a known ground truth (Figure 1). The phantom tensor field was simulated based on a fixed-width helical tube. In the shown experiment, the template was purposefully initialized to a helix with much larger radius, with only its endpoints constrained to approximately correct locations. Although this initial shape was a very different one, the algorithm was able to find a correct alignment, as demonstrated in figure 1.

Fig. 1
DTI Phantom Test

4.1.2 Human Brain Data

As stated in the introduction section, our motivation for using a template for fiber tract reconstruction is to be able to incorporate prior anatomical knowledge in the template matching process. Obviously, this approach hinges on the availability of an adequate template. One viable option is basing such a template on the postmortem WM fiber tract atlas [25] available with the SPM Anatomy Toolbox [26]. This atlas comprises probabilistic maps for ten WM tracts, based on histological staining in several postmortem human brains, and thus reflecting actual ground truth data.

Here we focus without loss of generality on two specific tracts, namely the cingulum and the fornix. For each of these tracts, the geometric fiber template An external file that holds a picture, illustration, etc.
Object name is nihms96332ig3.jpg i is initially chosen simply as the curve-skeleton [27] of the corresponding probability map P(An external file that holds a picture, illustration, etc.
Object name is nihms96332ig3.jpg i) (see Figure 3 for illustration), providing a compact representation of the tract shape. Now, given a pair of co-registered subject T1 and DT MR images, the task is to match the template to the subject data.

Fig. 3
Creation of a Fiber-Bundle Template

The matching procedure involves two steps as follows. The first step is a purely structural registration: since the atlas is aligned with the T1-weighted MNI reference brain [28], we register the reference brain with the subject T1 image using a non-linear method [29]. The associated deformation field is then applied to each An external file that holds a picture, illustration, etc.
Object name is nihms96332ig3.jpg i and P(An external file that holds a picture, illustration, etc.
Object name is nihms96332ig3.jpg i), generating their corresponding deformed instances Ti and P(An external file that holds a picture, illustration, etc.
Object name is nihms96332ig3.jpg i)′, respectively. Still, structural registration alone, however accurate, is not enough to ensure a correct tract reconstruction. Thus, it is followed by a second step in which Ti becomes the initial subject-specific template and is evolved to match the subject tensor data, using the deformable fiber evolution algorithm described in Section 3. Since the use of deformable bundle model, as opposed to individual fibers, improves robustness no noise, the fiber template is first converted to a bundle representation by sampling curves from a fixed-radius tube around each original fiber. This step also includes a moderate Laplacian smoothing applied to the fiber template (using the Laplacian operator from Section 3.2), to reduce artifacts and possible overfitting effects specific to Step 1.

Furthermore, the evolution procedure can benefit from having a (deformed) probability map P(An external file that holds a picture, illustration, etc.
Object name is nihms96332ig3.jpg i)′ to reduce the method’s sensitivity to the initial position, as well as its dependence on endpoint constraints. This is done by extending the notion of a template beyond merely a shape prior to also include a probabilistic prior, represented by a new energy term. This energy is based directly on P(An external file that holds a picture, illustration, etc.
Object name is nihms96332ig3.jpg i)′, simply measuring for each point the probability of not having a fiber passing through it, i.e.,


where [p with macron](x) is the probability of no fibers passing through x and ω(xi) is the Voronoi (dual) region of the vertex xi. One can easily verify that the minimizer of this energy is indeed the most likely curve given P(An external file that holds a picture, illustration, etc.
Object name is nihms96332ig3.jpg i)′, subject to endpoint constraints. The partial derivatives of Eprior are given by


where [nabla] [p with macron](x) can be approximated using finite differences.

Likewise, one can now incorporate a segmentation mask (either binary or probabilistic) of the cortico-spinal fluid (CSF) to prevent the fibers from drifting into the ventricles. This is especially useful since CSF-filled regions are characterized by high isotropic diffusion, and Ematch is in fact sensitive to diffusion tensor magnitude, not just orientation.

4.1.3 Results

Data Processing

HARDI data were acquired from 5 healthy adult subjects (age: 23.9 ± 1.4SD years) on a 4 Tesla Bruker Medspec MRI scanner using an optimized diffusion tensor sequence [31]. The study was approved by the institutional review boards (IRBs) at the Queensland Institute for Medical Research and at the University of California, Los Angeles; all subjects gave informed consent. 105 images were acquired: 11 baseline (b0) images with no diffusion sensitization (i.e., T2-weighted images) and 94 diffusion-weighted (DW) images (b-value 1159 s/mm2) in which gradient directions were evenly distributed on the hemisphere [2]. Imaging parameters were: TR/TE 92.3/8250 ms, 55 × 2mm contiguous slices, FOV = 23 cm. The reconstruction matrix was 128×128, yielding a 1.8×1.8 mm2 in-plane resolution. The total scan time was 14.5 minutes. The DT images were obtained from diffusion weighed data using the MedInria tool [32]; skull stripping in T1 MRI images was done in BrainSuite [33].

Reconstruction Results

Figure 4 shows several examples of the atlas-based tract reconstruction procedure described in Section 4.1.2. The superimposed T1 images suggest that the reconstructed bundles appear to follow anatomically valid trajectories, notwithstanding the proximity of e.g., the corpus callosum fibers which were not part of the template. Figure 5 provides an additional quantitative validation, depicting the effect of each of the algorithm’s steps, for one of the cases. Note that a relatively moderate deformation of the template shape may lead to a significantly better alignment to the diffusion tensor field, as illustrated by the corresponding plot of the misalignment energy. Furthermore, the mean energy value along a given pathway can be used as a likelihood measure, a claim supported by Ziyan et al. [13], who suggested a similar “Fiber-Tensor Fit” measure in a somewhat different context. Naturally, more conservative measures, such as FA, mean diffusivity (MD), etc., can be computed along the tracts as well.

Fig. 4
Cingulum and Fornix Bundles Reconstruction
Fig. 5
Deformable Alignment of the Fornix Template

In addition, Figure 6 provides a comparison of cingulum reconstruction with our method to streamline tachography performed in DTI Studio [30], using a specialized ROI placement protocol [8]. While both methods yielded approximately similar results for the upper part of the cingulum, the ROI-based approach failed to capture the entire C-shaped structure, which is likely to be due to lower FA values in those regions, and/or the high fiber curvature. In contrast, our algorithm does not involve any FA or curvature thresholds, and therefore was not subject to this shortcoming. While it could still be possible to capture the entire cingulum bundle using DTI tractography through a more complex combination of multiple ROIs, this would require even more extensive expert knowledge and manual labor, as opposed to the proposed automatic method. Note also the multiple “straying” fibers in the tractography result, requiring an additional filtering step and complicating intersubject comparisons, whereas the template matching approach avoids these issues by design. Finally, it should pointed out that the described procedure can be performed using any other probabilistic WM atlas, if available.

Fig. 6
Comparison to Fiber Tracking with Manual ROIs

4.2 Tract Based Morphometry

The above results are immediately suitable for tract-based morphometric analysis. While the actual analysis is out of the scope of this paper, several possibilities are available. Since each fiber-bundle result has a one-to-one correspondence to the template, the latter becomes a common intersubject “coordinate system”, to which each point of any reconstructed bundle. This property can be used to compute local tract-based maps of scalar measures such as FA, MD, etc., as well as for more complex comparisons of tract shapes and deformation-based morphometry. In fact, the deformable evolution approach makes it straightforward to quantify the deformation of the discrete bundle at each vertex w.r.t. the initial template shape. For instance, one can define the cumulative local deformation as the path integral of the incremental bending energy from Eq. 7 at a given vertex.

5 Discussion and Conclusions

5.1 Potential Clinical Applications

Besides neurological research, white matter connectivity analysis is becoming an important tool in neurosurgical planning. The key distinction in clinical applications is the requirement to correctly analyze pathological cases, e.g., detect WM connectivity in the presence of brain tumors. However, the underlying anisotropic diffusion assumptions no can no longer be relied on in the presence of tumor tissue (see, e.g. [34]). In fact, to our knowledge, there are no adequate methods to date to reliably reconstruct fiber tract connectivity in demyelinated tissue. Nevertheless, solutions based on basic streamline tractography approaches are increasingly employed in image-guided neurosurgery. While visually useful, such tractography results should be considered with great caution. For instance, if the detected fibers appear to be disrupted by a tumor, that is not yet a sufficient evidence that there are no unobserved fibers going though it, which poses a risk of visual misinterpretation [35]. By contrast, our proposed method does not rely solely on local diffusion measurements, and therefore will always produce a continuous fiber-bundle that may still pass through the tumor region. Then, complementary measures (e.g., FA, MD) could be used to locally assess the certainty of the result. While it is clearly premature to suggest the template-based technique for clinical practice at this point, we believe this approach could ultimately be well suited for this application. Furthermore, the flexibility to incorporate more sophisticated deformation priors, e.g., based on tumor growth models, suggests potential extensions to specific cases.

5.2 Conclusions

We have presented a novel approach to WM fiber-bundle reconstruction and connectivity analysis, addressing the problem from the template matching point of view. Given an adequate and sufficiently detailed template (or atlas) of the WM fiber tracts, the proposed method has the potential of simultaneously addressing the problems of fiber tracking, fiber-bundle segmentation and anatomical labeling. Reliance on a “good” template, while being a distinctive feature of this method, can also be a limitation. For instance, the reconstruction procedure of the cingulum and fornix bundles in Section 4.1.2 could be made more robust by explicitly incorporating e.g., a corpus callosum template into our prior. Thus, creation of a more complete atlas with built-in tract-specific knowledge becomes a major avenue of future work. Another potential extension of the method is template-based DTI clustering (or segmentation), as the proposed scheme does not attempt to accurately recover the bundle width. Furthermore, since the method is inherently suitable for tract-based morphometric analysis, we intend to extend it for such applications. Finally, incorporating the recent advances in DTI registration [36] into the template matching process is yet another interesting direction.

Obviously, this method has its limitations. Being a template-based approach, it can be useful in the reconstruction of known tracts, but not in the detection of pathways that are not expected a priori. Also, while the use of a template can help improve robustness to noise in the DW images, an overconstrained template may not fully capture subject-specific pathologies. Indeed, it is not always clear how the prior knowledge should be weighted in relation to the observed data. Thus, the confidence in the final match should always be assessed using the available measures, i.e., the misalignment energy, FA, MD, etc., with low confidence values indicating that further inspection may be required. A potential way to generalize such a strategy is to incorporate confidence measures in the algorithm directly as weighting criteria, using e.g., a Bayesian filtering framework. Another issue is the sensitivity to the choice of the initial template: while this problem is alleviated by the use of a probabilistic prior (when available), the use of multiple shape priors (e.g., as in [37]) is likely to help better capture subject variability. Finally, we hope that more DTI datasets with expert-delineated tracts will soon become available, which should also allow a more rigorous validation of the proposed method, by “learning” the best set of parameter values.

The “holy grail” of the template matching approach is a combined atlas of both white and grey matter structures, all integrated in a single geometric deformable template. With this goal in mind, we are currently experimenting with the integration of the proposed method with a deformable cortical surface matching framework, similarly based on geometric flows [38].


This work is funded by the National Institutes of Health through the NIH Roadmap for Medical Research, Grant U54 RR021813. Information on the National Centers for Biomedical Computing can be obtained from The authors are also thankful to Nathan Hageman for his DTI phantom data, as well as to Yonggang Shi and Boris Gutman for sharing their processing tools. Additional support was provided by the National Institute of Child Health and Human Development (RO1 HD050735), and the National Health and Medical Research Council, Australia (Grant 496682).

A Laplacian and biLaplacian Operators for Discrete Curves

Observe that the matrix K from Eq. 4 can be rewritten as K = M−1L, where M is the discrete mass matrix and


This relationship is consistent with what we know about discrete surfaces [16,17], with L being the corresponding definition of the discrete vertex-based Laplacian operator for curves. Consequently, we can simplify Eq. 5 as follows:


where LTM−1L is, in fact, the discrete biLaplacian operator. The same simplification applies to Eq. 8 and 9.

B Energy Gradient Derivation

From Eq. 2 we know: Ematch = EDαEB. In the discrete case, ED is given by Eq. 3 and its partial derivatives are as follows:


The discrete biharmonic energy is simply EB is defined in Appendix A above, and the corresponding partial derivatives can be written simply as EBX=(LTM1L)X.


1In practice, a very small time step may be required to reach the global optimum with explicit Euler integration.

2CGAL, Computational Geometry Algorithms Library,

3This symmetric property becomes important, e.g., in the definition of the Sobolev norm which makes use of the Laplace-Beltrami operator.

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.


1. Basser P, Pierpaoli C. Microstructural and Physiological Features of Tissues Elucidated by Quantitative-Diffusion-Tensor MRI. Journal of Magnetic Resonance, Series B. 1996;111(3):209–219. [PubMed]
2. Mori S, Crain B, Chacko V, van Zijl P. Three-dimensional tracking of axonal projections in the brain by magnetic resonance imaging. Ann Neurol. 1999;45(2):265–9. [PubMed]
3. Basser P, Pajevic S, Pierpaoli C, Duda J, Aldroubi A. In vivo fiber tractography using DT-MRI data. Magnetic Resonance in Medicine. 2000;44(4):625–632. [PubMed]
4. Lenglet C, Deriche R, Faugeras O. Inferring white matter geometry from diffusion tensor MRI: Application to connectivity mapping. ECCV. 2004:127–140.
5. O’Donnell L, Haker S, Westin C-F. New approaches to estimation of white matter connectivity in diffusion tensor MRI: Elliptic pdes and geodesics in a tensor-warped space. Fifth International Conference on Medical Image Computing and Computer-Assisted Intervention (MICCAI’02); Tokyo, Japan. 2002. pp. 459–466.
6. Heiervang E, Behrens T, Mackay C, Robson M, Johansen-Berg H. Between session reproducibility and between subject variability of diffusion MR and tractography measures. Neuroimage. 2006;33(3):867–877. [PubMed]
7. Toga AW, Mazziotta JC. Brain Mapping: The Methods. 2. Academic Press; 2002.
8. Hua K, Zhang J, Wakana S, Jiang H, Li X, Reich D, Calabresi P, Pekar J, van Zijl P, Mori S. Tract probability maps in stereotaxic spaces: Analyses of white matter anatomy and tract-specific quantification. Neuroimage. 2008;39(1):336–347. [PMC free article] [PubMed]
9. Zhang W, Olivi A, Hertig S, van Zijl P, Mori S. Automated fiber tracking of human brain white matter using diffusion tensor imaging. Neuroimage. 2008;42(2):771–777. [PMC free article] [PubMed]
10. Wakana S, Jiang H, Nagae-Poetscher L, van Zijl P, Mori S. Fiber Tract–based Atlas of Human White Matter Anatomy. Radiology. 2004;230:77. [PubMed]
11. Maddah M, Mewes A, Haker S, Grimson W, Warfield S. Automated atlas-based clustering of white matter fiber tracts from DTMRI. MICCAI. 2005:188–195. [PubMed]
12. O’Donnell L, Westin C. Automatic Tractography Segmentation Using a High-Dimensional White Matter Atlas, Medical Imaging. IEEE Transactions on. 2007;26(11):1562–1575. [PubMed]
13. Ziyan U, Sabuncu M, O’Donnell L, Westin C. Nonlinear registration of diffusion MR images based on fiber bundles. MICCAI. 2007:351–358. [PubMed]
14. Ziyan U, Sabuncu M, Eric L, Grimson W, Westin C-F. A Robust Algorithm for Fiber-Bundle Atlas Construction. ICCV. 2007:1–8.
15. Fletcher P, Tao R, Jeong W, Whitaker R. A Volumetric Approach to Quantifying Region-to-Region White Matter Connectivity in Diffusion Tensor MRI. Information Processing in Medical Imaging 2007 Conference Proceedings; 2007. p. 346. [PubMed]
16. Meyer M, Desbrun M, Schröder P. A. Barr, Discrete differential-geometry operators for triangulated 2-manifolds. Proc. of the Int. Workshop on Vis. and Math; 2002.
17. Wardetzky M, Bergou M, Harmon D, Zorin D, Grinspun E. Discrete quadratic curvature energies. Computer Aided Geometric Design. 2007;24(8–9):499–518.
18. Desbrun M, Kanso E, Tong Y. Discrete differential forms for computational modeling. International Conference on Computer Graphics and Interactive Techniques. 2006:39–54.
19. Eckstein I, Pons J-P, Tong Y, Kuo C-C, Desbrun M. Generalized surface flows for mesh processing. Symposium on Geometry Processing; 2007.
20. Grinspun E, Hirani AN, Desbrun M, Schröder P. Discrete shells. ACM SIGGRAPH Symposium on Computer Animation; 2003. pp. 62–67.
21. Charpiat G, Keriven R, Pons J-P, Faugeras O. Designing spatially coherent minimizing flows for variational problems based on active contours, in: ICCV. 2005;(2):1403–1408.
22. Sundaramoorthi G, Yezzi A, Mennucci A. Sobolev Active Contours. IJCV. 2007;73(3):345–366.
23. Zhou K, Huang J, Snyder J, Liu X, Bao H, Guo B, Shum H. Large mesh deformation using the volumetric graph Laplacian. ACM SIGGRAPH 2005. 2005;24(3):496–503.
24. Clarenz U, Rumpf M, Telea A. Finite elements on point based surfaces. Proc. Eurographics Symp. Point-Based Graphics; 2004. pp. 201–211.
25. Bürgel U, Amunts K, Hoemke L, Mohlberg H, Gilsbach J, Zilles K. White matter fiber tracts of the human brain: Three-dimensional mapping at microscopic resolution, topography and intersubject variability. Neuroimage. 2006;29(4):1092–1105. [PubMed]
26. Eickhoff S, Stephan K, Mohlberg H, Grefkes C, Fink G, Amunts K, Zilles K. A new SPM toolbox for combining probabilistic cytoarchitectonic maps and functional imaging data. Neuroimage. 2005;25(4):1325–1335. [PubMed]
27. Shi Y, Lai R, Krishna S, Sicotte N, Dinov I, Toga AW. Anisotropic Laplace-Beltrami Eigenmaps. Bridging Reeb Graphs and Skeletons; MMBIA. 2008. pp. 1–7. [PMC free article] [PubMed]
28. Evans A, Collins D, Mills S, Brown E, Kelly R, Peters T. 3D statistical neuroanatomical models from 305 MRI volumes. Nuclear Science Symposium and Medical Imaging Conference; 1993. pp. 1813–1817.
29. Leow A, Yanovsky I, Chiang M, Lee A, Klunder A, Lu A, Becker J, Davis S, Toga A, Thompson P. Statistical Properties of Jacobian Maps and the Realization of Unbiased Large-Deformation Nonlinear Image Registration. IEEE TMI. 26(06) [PubMed]
30. Jiang H, van Zijl P, Kim J, Pearlson G, Mori S. DtiStudio: Resource program for diffusion tensor computation and fiber bundle tracking. Computer Methods and Programs in Biomedicine. 2006;81(2):106–116. [PubMed]
31. Jones D, Horsfield M, Simmons A. Optimal Strategies for Measuring Diffusion in Anisotropic Systems by Magnetic Resonance Imaging. Magnetic Resonance in Medicine. 1999;42(3):515–525. [PubMed]
33. Shattuck D, Leahy R. BrainSuite: an automated cortical surface identification tool. Med Image Anal. 2002;6(2):129–142. [PubMed]
34. Yamasaki F, Kurisu K, Satoh K, Arita K, Sugiyama K, Ohtaki M, Takaba J, Tominaga A, Hanaya R, Yoshioka H, et al. Apparent Diffusion Coefficient of Human Brain Tumors at MR Imaging. Radiology. 2005:2353031338. [PubMed]
35. Kinoshita M, Yamada K, Hashimoto N, Kato A, Izumoto S, Baba T, Maruno M, Nishimura T, Yoshimine T. Fiber-tracking does not accurately estimate size of fiber bundle in pathological condition: initial neurosurgical experience using neuronavigation and subcortical white matter stimulation. Neuroimage. 2005;25(2):424–429. [PubMed]
36. Chiang M, Leow A, Klunder A, Dutton R, Barysheva M, Rose S, McMahon K, de Zubicaray G, Toga A, Thompson P. Fluid Registration of Diffusion Tensor Images Using Information Theory, Medical Imaging. IEEE Transactions on. 2008;27(4):442–456. [PMC free article] [PubMed]
37. Etyngier P, Segonne F, Keriven R. Shape Priors using Manifold Learning Techniques. ICCV. 2007:1–8.
38. Eckstein I, Joshi A, Leahy R, Kuo C, Desbrun M. Generalized Surface Flows for Deformable Registration and Cortical Matching; MICCAI. 2007. pp. 692–700. [PMC free article] [PubMed]