Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2720421

Formats

Article sections

Authors

Related links

Neuroimage. Author manuscript; available in PMC 2010 August 1.

Published in final edited form as:

Published online 2009 February 11. doi: 10.1016/j.neuroimage.2009.01.065

PMCID: PMC2720421

NIHMSID: NIHMS96332

Ilya Eckstein,^{a,}^{*} David W. Shattuck,^{a} Jason L. Stein,^{a} Katie L. McMahon,^{b} Greig de Zubicaray,^{b} Margaret J. Wright,^{c} Paul M. Thompson,^{a} and Arthur W. Toga^{a}

The publisher's final edited version of this article is available at Neuroimage

See other articles in PMC that cite the published article.

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.

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 [2–5]. 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.

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.

Let us first consider a simpler case of a single-fiber template, modeled as a curve : [0,1] → ^{3}. 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 Φ() = , where is a curve *maximizing* the alignment functional (a.k.a. anisotropic Dirichlet energy):

$${E}_{D}(\mathcal{Y})={\int}_{0}^{1}{\mathcal{Y}}^{\prime}{(t)}^{T}D(\mathcal{Y}(t)){\mathcal{Y}}^{\prime}(t)dt\phantom{\rule{0.16667em}{0ex}}\text{subject}\phantom{\rule{0.16667em}{0ex}}\text{to}\phantom{\rule{0.16667em}{0ex}}\mathcal{Y}(0)={y}_{0},\mathcal{Y}(1)={y}_{1},$$

(1)

where ′(*t*) denotes the tangent vector to
(*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 *E _{D}* is unbounded, we must constrain the solution by adding a smoothness term as follows:

$${E}_{\mathit{match}}(\mathcal{X},D)={E}_{D}(\mathcal{X})-\alpha \underset{{E}_{B}(\mathcal{X})}{\underbrace{{\int}_{0}^{1}\mid \mathrm{\Delta}\mathcal{X}(t){\mid}^{2}dt}},$$

(2)

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 *E _{B}*(), corresponds to the

Note that *E _{D}* is conceptually similar to the

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** = {**x*** _{i}*}

$${E}_{D}(\mathbf{X})=\sum _{i}\frac{{{\mathbf{e}}_{i}}^{T}{\overline{D}}_{i}{\mathbf{e}}_{i}}{\mid {\mathbf{e}}_{i}{\mid}^{2}}{l}_{i},\phantom{\rule{0.38889em}{0ex}}\text{where}\phantom{\rule{0.16667em}{0ex}}{\mathbf{e}}_{i}={\mathbf{x}}_{i+1}-{\mathbf{x}}_{i},$$

(3)

*D̄ _{i}* is the average value of

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 *x _{i}*

$$\mathbf{K}({x}_{i})=\frac{1}{{m}_{i}}\left(\frac{{e}_{i}}{\mid {e}_{i}\mid}-\frac{{e}_{i-1}}{\mid {e}_{i-1}\mid}\right),$$

(4)

where *m _{i}* = (

$${E}_{B}(\mathbf{X})={(\mathbf{KX})}^{T}\mathbf{M}\phantom{\rule{0.16667em}{0ex}}(\mathbf{KX})={\mathbf{X}}^{T}({\mathbf{K}}^{T}\mathbf{MK})\mathbf{X},$$

(5)

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

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

$${\mathbf{X}}_{t+dt}={\mathbf{X}}_{t}+dt\phantom{\rule{0.16667em}{0ex}}{\mathbf{M}}^{-1}\frac{\partial {E}_{\mathit{match}}}{\partial \mathbf{X}}({\mathbf{X}}_{t}),$$

(6)

with the derivation of ${\scriptstyle \frac{\partial {E}_{\mathit{match}}}{\partial \mathbf{X}}}$ found in Appendix B.

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 *X _{t}*

$$\begin{array}{l}{E}_{IB}(t+dt)\approx \parallel {\mathbf{K}}_{t}({\mathbf{X}}_{t+dt}-{\mathbf{X}}_{t})\mid {\mid}_{{L}^{2}}^{2}\\ =\parallel {\mathbf{K}}_{t}{\mathbf{U}}_{t}\mid {\mid}_{{L}^{2}}^{2}={{\mathbf{U}}_{t}}^{T}({\mathbf{K}}^{T}\mathbf{MK})\phantom{\rule{0.16667em}{0ex}}{\mathbf{U}}_{t}.\end{array}$$

(7)

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 **U*** _{t}*, as opposed to a shape. Furthermore, it represents a norm in the space of vector fields in

$${\mathbf{X}}_{t+dt}={\mathbf{X}}_{t}+dt\phantom{\rule{0.16667em}{0ex}}{\mathcal{L}}^{-1}\frac{\partial {E}_{\mathit{match}}}{\partial \mathbf{X}}({\mathbf{X}}_{t}),$$

(8)

where = **M**+λ**K**^{T}**MK**. 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.

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 *T _{i}* is computed (using CGAL

$$\mathcal{L}=\mathbf{M}+\lambda {\mathbf{K}}^{T}\mathbf{MK}-\mu \stackrel{\sim}{\mathbf{L}}.$$

(9)

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.

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.

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
* _{i}* is initially chosen simply as the curve-skeleton [27] of the corresponding probability map

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
* _{i}* and

Furthermore, the evolution procedure can benefit from having a (deformed) probability map *P*(
* _{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

$${E}_{\mathit{prior}}(\mathcal{X})={\int}_{0}^{1}\overline{p}(\mathcal{X}(t))dt\approx \sum _{i}{\int}_{\omega ({\mathbf{x}}_{i})}\overline{p}(x)dx,$$

where (*x*) is the probability of no fibers passing through *x* and *ω*(**x*** _{i}*) is the Voronoi (dual) region of the vertex

$$\frac{\partial {E}_{\mathit{prior}}}{\partial {\mathbf{x}}_{i}}={\int}_{\omega ({\mathbf{x}}_{i})}\nabla \overline{p}(x)dx,$$

where (*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 *E _{match}* is in fact sensitive to diffusion tensor magnitude, not just orientation.

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 (*b*_{0}) images with no diffusion sensitization (i.e., T2-weighted images) and 94 diffusion-weighted (DW) images (b-value 1159 *s*/*mm*^{2}) 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].

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.

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.

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.

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.

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 http://nihroadmap.nih.gov/bioinformatics. 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).

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

$$\mathbf{L}({\mathbf{x}}_{i})=\left(\frac{{\mathbf{e}}_{i}}{\mid {\mathbf{e}}_{i}\mid}-\frac{{\mathbf{e}}_{i-1}}{\mid {\mathbf{e}}_{i-1}\mid}\right)=(\frac{1}{\mid {\mathbf{e}}_{i}\mid}+\frac{1}{\mid {\mathbf{e}}_{i-1}\mid}){x}_{i}-\frac{{\mathbf{x}}_{i-1}}{\mid {\mathbf{e}}_{i-1}\mid}-\frac{{\mathbf{x}}_{i}}{\mid {\mathbf{e}}_{i}\mid}.$$

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:

$${E}_{B}(\mathbf{X})={\mathbf{X}}^{T}({\mathbf{K}}^{T}\mathbf{MK})\mathbf{X}={\mathbf{X}}^{T}({\mathbf{L}}^{T}{\mathbf{M}}^{-1}\mathbf{L})\mathbf{X},$$

where **L**^{T}**M**^{−1}**L** is, in fact, the discrete biLaplacian operator. The same simplification applies to Eq. 8 and 9.

From Eq. 2 we know: *E _{match}* =

$$\begin{array}{l}\frac{\partial {E}_{\mathit{match}}}{\partial {\mathbf{x}}_{i}}={\theta}_{i-1}({\mathbf{x}}_{i}-{\mathbf{x}}_{i-1})-{\theta}_{i}({\mathbf{x}}_{i+1}-{\mathbf{x}}_{i}),\\ \text{where}\phantom{\rule{0.16667em}{0ex}}{\theta}_{i}=2\frac{\mid {\mathbf{e}}_{i}{\mid}^{2}{\overline{D}}_{i}-({{\mathbf{e}}_{i}}^{T}{\overline{D}}_{i}{\mathbf{e}}_{i}){\mathbf{I}}_{3}}{\mid {\mathbf{e}}_{i}{\mid}^{4}}.\end{array}$$

The discrete biharmonic energy is simply *E _{B}* is defined in Appendix A above, and the corresponding partial derivatives can be written simply as
${\scriptstyle \frac{\partial {E}_{B}}{\partial \mathbf{X}}}=({\mathbf{L}}^{T}{\mathbf{M}}^{-1}\mathbf{L})\mathbf{X}$.

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

^{2}CGAL, Computational Geometry Algorithms Library, http://www.cgal.org

^{3}This 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]

32. Fillard P, Toussaint N, Pennec X. MEDINRIA: DT-MRI PROCESSING AND VISUALIZATION SOFTWARE. 2006.

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]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |