Research Article  Open Access
S. Jbabdi, P. Bellec, R. Toro, J. Daunizeau, M. PélégriniIssac, H. Benali, "Accurate Anisotropic Fast Marching for DiffusionBased Geodesic Tractography", International Journal of Biomedical Imaging, vol. 2008, Article ID 320195, 12 pages, 2008. https://doi.org/10.1155/2008/320195
Accurate Anisotropic Fast Marching for DiffusionBased Geodesic Tractography
Abstract
Using geodesics for inferring white matter fibre tracts from diffusionweighted MR data is an attractive method for at least two reasons: (i) the method optimises a global criterion, and hence is less sensitive to local perturbations such as noise or partial volume effects, and (ii) the method is fast, allowing to infer on a large number of connexions in a reasonable computational time. Here, we propose an improved fast marching algorithm to infer on geodesic paths. Specifically, this procedure is designed to achieve accurate front propagation in an anisotropic elliptic medium, such as DTI data. We evaluate the numerical performance of this approach on simulated datasets, as well as its robustness to local perturbation induced by fiber crossing. On real data, we demonstrate the feasibility of extracting geodesics to connect an extended set of brain regions.
1. Introduction
For decades, dissection, lesion studies, or axonal transport of tracers have been the only available techniques for studying the brain's anatomical connections. It is not surprising that due to the invasiveness of these methods, most of the data concerning the largescale, white matter tracts of the brain were collected on animals, for example, cats [1] or monkeys [2], while structural data for the human brain were largely missing [3]. Diffusion weighted MR imaging now offers a propitious and unique framework to explore noninvasively the organisation of white matter in the living human brain [4, 5]. Despite the poor spatial resolution of this technique, already diffusion data are beginning to inform us about human brain largescale connections [6–8] and how they relate to the functional role of cortical and subcortical networks [9, 10].
Inferring on white matter architecture from diffusion data relies on the properties of water diffusion in the tissues. Water molecules diffuse more easily along the fibre tracts than across them, and this anisotropy is captured by the diffusionweighted MR signal. Inferring on connexions given this local feature is challenging, since the observations (diffusion properties) are indirectly related to the actual structure (axonal orientations, size, and packing). The tractography algorithms use the information of directionality contained in diffusion data to infer connectivity between brain regions. Usually, information about the orientation of white matter fibres is estimated locally, via models (e.g., diffusion tensor imaging (DTI) [11], mixture models [12], or partial volume models [13, 14]) or in a modelfree manner (e.g., Qball imaging [15]). Fibre tracking consists then in inferring connexions between distant brain regions, given this local orientation. This can be done either in a deterministic way, by trusting the local orientation information and following these directions until reaching a target region (i.e., streamline tractography [16–19]), or in a probabilistic way, by building distributions of connexions, using local probabilistic models for fibre orientation distributions [13, 14, 20].
In both cases, when tracking a fibre between two regions of the brain, these algorithms start in one seed region, and try to find the tracts, or distribution of tracts, that will end up in the target region. In cases where the local orientation information present in the diffusion data is consistent with the presence of this pathway, then these tractography algorithms manage in general to recover the connexion between the seed and the target. However, it often happens that in some parts of the trajectory, the local diffusion information no longer supports the presence of the pathway. This can either be due to a high level of noise compared to the actual signal, or to the presence of a high number of crossing fibres heterogeneous in their orientations. This issue is crucial in streamlining algorithms, and is also met in probabilistic algorithms when a single orientation per voxel is modelled [21]. The problem with those algorithms is that when tracking from a seed, the algorithm has no information about the region it will end up in.
A possible solution to the problem of local perturbations in the diffusion data may be provided by global tractography, that is, optimising a global criterion while seeking for connexions. A global tractography algorithm can potentially overcome errors in estimating local structure, because its goal is to connect two given regions. In other words, if we tell the algorithm which connexion we are looking at, that is, which pair of regions is to be connected, it is better at finding it. Geodesic tractography (GT), first proposed by Parker et al. [22], falls into this category. GT is based on the hypothesis that brain fibers can be interpreted as minimal distance paths (geodesics) for a metric derived from the water diffusion profile. This distance criterion is global by definition.
The basic idea for constructing a geodesic in a metric space is to build a distance field from a seed region, the very same region one would use as a seed for streamline tractography. This is done by solving the socalled Eikonal equation, a partial differential equation (PDE) that describes the time of arrival at each point of the space, as a function of the local speed. In a constant speed field, this PDE can be easily integrated, and the geodesics are simply straight lines. When the speed varies across the space, the geodesics can curve, preferring high local speed locations to decrease the arrival time. Finally, if the speed depends on the direction of travel (e.g., along versus across a fibre tract), then the PDE is said to be anisotropic.
Solving the Eikonal equation in a heterogeneous and highly anisotropic medium, as is the human brain, is a technically challenging problem [23]. This is especially true if one uses singlepass algorithms, which is particularly important when dealing with data containing hundreds of thousands of voxels. There have been a few attempts at solving this problem in the context of diffusionbased tractography [22, 24–27].
We describe a method for constructing geodesics in an anisotropic medium, and apply it to the problem of DTIbased tractography. This method relies on works in optimal path planning [28] and, more recently, vessel extraction in 3D angiography images [29]. It has been shown to be very accurate in anisotropic media [29], and requires less computation than the exact method proposed in Sethian and Vladimirsky [30] in a general framework for anisotropic optimal path planning. The main contribution of this work is to show how this method applies to the case of an elliptic medium, where the algorithm performs extremely well both in terms of accuracy and efficiency, as shown in the simulations. We also show the feasibility of applying such method to the extraction of structural connectivity in an extended brain network using diffusion data from a healthy subject.
2. Methods
In this section, we will give some theoretical background on geodesics and the Eikonal equation, and describe a singlepass algorithm for building geodesics.
2.1. Geodesics and the Eikonal Equation
A geodesic is a pathway minimising an integral of the form where describes an infinitesimal distance along a pathway , relative to a metric tensor .
Now, let be the arrival time function starting from a location , that is, is equal to the minimum value of the integral along a geodesic connecting to . Then, the arrival time function and the geodesics satisfy these two fundamental equations: where is the spatial gradient of . Equation (2a) is the anisotropic version of the socalled Eikonal equation. In the isotropic case, this equation is usually written , where is the local speed. Hence, this equation tells us two things: (i) it is a generalisation of the speed equation, stating that the time of arrival is inversely proportional to the speed, and (ii) changing the local metric tensor can be seen as changing the local speed. Equation (2b) shows that the tangent of the geodesic lines is parallel to the gradient of the time of arrival function with respect to the inverse metric. This is very important because it gives us a convenient way to reconstruct geodesics from any point in space, given the solution to the Eikonal equation. Figure 1 shows example geodesics in an isotropic space composed of two subsets with different local speeds.
Proof. Recall that the function is the minimum value of along the geodesic from point to an arbitrary point : A general variation of (3) is given (see, e.g., [31]) as Since we have integrated along a geodesic, the second term on the righthand side of (4) equals zero (Euler condition). We obtain Equation (2b) directly follows. Finally, and using the symmetry of the metric tensor , we get the Eikonal equation:
Equations (2a) and (2b) summarise the two steps for building geodesics: (i) solve the Eikonal equation for , given a metric tensor and a starting point ; (ii) construct geodesics between any given point and the starting point by following the gradient of with respect to the inverse metric .
2.2. FastMarching Algorithm
A few algorithms have been proposed in the literature for computing the function on a discrete grid. The most popular are Tsitsiklis's method [28] and Sethian's method [32], which are based on the construction of the time of arrival function using front propagation. These methods are also referred to as fast marching methods because they construct the function in a singlepass through the grid nodes. Tsitsiklis's method relies on (1) while Sethian's method uses the Eikonal equation (2a). Both methods are suitable in the case of isotropic media, that is, where the metric is proportional to the identity matrix, but they fail in anisotropic media [23]. An exact scheme to deal with anisotropy has been proposed by Sethian and Vladimirsky [30], but while remaining a singlepass algorithm, it still requires a computational effort that is growing with the amount of anisotropy. A variant of the initial fastmarching algorithm of Tsitsiklis [28] has been proposed to deal with anisotropic media [29], which is more computationally efficient than the exact scheme of Sethian [30]. Yet, it relies on a generic optimisation procedure that was undocumented for the special case of the elliptical media we face with DTI tractography. We extended this method by deriving a solution to the optimisation procedure in this case.
The general idea of the fastmarching algorithm was borrowed from the graph theory. It is a direct extension of Dijkstra's algorithm for finding minimal paths in a graph [33]. The algorithm relies on a very simple observation: suppose that the time of arrival is known inside a close set of grid nodes (a set we will refer to as the known set). Then, the first nodes that will be encountered by the propagating front are the nodes on the edge of the known set (this narrow band of grid nodes will be called the trial set). Secondly, the first node that will be encountered by the propagating front is the closest one to known (in terms of geodesic distance), and crucially, there will be no other way to make this distance smaller after propagating the front further. This means that the arrival time at this voxel will not change, and can be frozen. In other words, the value of the time of arrival can be calculated, starting from , in a singlepass through the voxels, only by considering, at each iteration, the neighbouring voxels of the propagating front. The other voxels (the far set) are not examined. Figure 2(a) schematises this front propagation scheme. The fastmarching algorithm is summarised in the appendices.
(a)
(b)
(c)
The crucial step in this front propagation is the computation of the distance between the front and the neighbouring voxels in the trial set. In our case, this distance is anisotropic, and we cannot use the standard methods, because they rely on the assumption that the gradients of are parallel to its geodesic lines (see [23] for further details). To account for the anisotropy, we consider a set of simplexes (triangles) that cover the whole neighbourhood around a voxel of the narrow band [29], and minimise the distance function between the simplexes and that voxel (see Figures 2(b) and 2(c)). The introduction of these simplexes allows to describe the trajectories on a continuous rather than a discrete grid. The definition of a simplex neighbouring a point is simply a set of three points that are 26 neighbours of , defining a triangle that we denote . There are 48 such triangles around for the 26 connexities (Figure 2(c)). The procedure for computing the anisotropic distance between the propagating front and the voxels in the trial set is given in the appendices.
During the updating procedure, the time of arrival at a voxel of the trial set is calculated from its neighbours on a simplex using an approximation (strictly speaking, two approximations!). Normally, if the geodesic passing by comes from simplex , then the time of arrival is given by We use a parametric approximation to this formula, given by the minimisation of the following function: where is the quadratic norm with respect to the metric and . Equation (8) follows the approximations of Tsitsiklis [28]. Term (I) approximates the distance from the starting point to the simplex centre of mass as a weighted sum of the distances to the nodes of the simplex. Term (II) approximates the remaining distance by considering the local metric as being constant, equal to its value at .
Minimising in the simplex can be written as a constrained optimisation problem that can be solved explicitly, since and the simplex are convex. The analytical solution is detailed in the appendices.
2.3. How to Choose the Metric?
In the GT framework, we make the hypothesis that white matter fibres are geodesics with respect to a metric tensor. But so far, we have not specified which metric tensor we mean. In DTI, the inverse tensor () seems to be the natural choice. Intuitively, water molecules diffusion is faster along the tract than across them. When inverting the diffusion tensor, the highest eigenvalues become the lowest, and the shortest distance is parallel to the fibres. One can also notice that the inverse tensor defines a metric in a Riemannian space that induces a LaplaceBeltrami operator (generalisation of the Laplace operator) which is encountered in the diffusion equation [25, 34].
However, the inverse tensor is not suitable in all circumstances. Consider the situation described in Figure 3 were a circular tract of radius connects points A and B, with diffusion tensors tangent to the tract having the same shape. Suppose the rest of the space is isotropic, with the same mean diffusion as along the tract. If one considers the inverse tensor metric , the distance between A and B through the circular path is where is the largest eigenvalue of the tensors along the circular pathway. On the other hand, the straight line distance between A and B is equal to . Hence, a necessary condition for the circular tract to be a geodesic is that its length is smaller than a straight line, that is, which leads to , that is, a condition on the tensor shape to be peaky enough. Of course, one can imagine that even if this condition is satisfied, a geodesic path might certainly lie somewhere in between a straight line and the circular line, as shown in Figure 4. Which metric to choose is hence still debatable. Nonetheless, in our simulations and real data applications, we will use the inverse diffusion tensor as a metric for defining geodesics.
(a)
(b)
(c)
(d)
(e)
(f)
3. Applications
3.1. Simulations
We have evaluated the GT method on simulated data. The purpose of these simulations is twofold. First, they show how the anisotropic fastmarching algorithm performs on elliptic media, in both homogeneous field (where the analytical solution is available) and a heterogeneous field. Second, they allow to compare GT with streamlining in cases where the data present local perturbations (crossing fibres).
In a homogeneous medium, where the data support the same diffusion tensor in every voxel, the analytic solution to the Eikonal equation is given by It is easy to check that in this case, and . We generated a tensor where the two smaller eigenvalues are equal, and gradually increased the anisotropy. Figure 5 shows the level curves of the analytic versus the numerical solution to the Eikonal equation. The two solutions are very close even for a large anisotropy, corresponding to a ratio of 50 between the largest and the lowest tensor eigenvalues. Table 1 summarises the mean and standard deviations of the relative error for different values of the anisotropy, which is expressed both in terms of the ratio between the largest and the lowest tensor eigenvalue, or in terms of the more widely used fractional anisotropy (FA, see, e.g., [35]).

(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
In a heterogeneous medium, such an analytical solution does not exist. However, we can verify that the Eikonal equation is satisfied, that is, is equal to one. We used the same circular tensor field as shown in Figure 4. In Table 1, we show the mean and standard deviations of for different anisotropies. Notice that these are close to one, but with a higher deviation from one with increasing anisotropy.
Finally, we show results of GT in the case of local perturbations. We generated a tensor field simulating a crossing fibre situation. The zone where the two fibres cross has a diffusion tensor that is the average of the two crossing fibres' tensors. We increased the crossing fibre area and compared the behaviour of GT to streamlining tractography (Figure 6). As expected, because the streamlining simply follows the direction of highest diffusion given by the tensor, the fibre trajectory was deviated. In the case of GT, there was little, if any, deviation from the straight line.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
3.2. Real Data
Acquisition
Data from a single healthy subject were acquired at Service de Neuroradiologie (CHNO des
QuinzeVingts, Paris). Six gradient weighted and one T_{2}weighted images were acquired
on a 1.5 Tesla MR Scanner (GE Signa) using the following scan parameters:
image matrix, inplane pixel size; slice thickness;
; milliseconds; Number of averages = 8. Thirtysix contiguous
slices covering the whole brain were acquired. The total scanning time was
approximately 14 minutes.
Regions of Interest
Five hundred and sixtyseven () regions
covering the whole cortex were manually selected in the DTI space. Each region
was represented by a single voxel. The anatomical localization of these regions
is shown in Figure 7. We performed a front propagation from each region,
which provided the distance functions . Then back propagation allowed us to construct the geodesics
connecting the whole set of voxel pairs. We computed a heuristic connectivity
index consisting of the mean diffusivity along each geodesic, multiplied by the
mean FA along the pathways.
In order to better visualize this anatomical
connectivity index in a matrix form, the set of brain regions were grouped with
respect to their localization. The regions were divided into five groups,
including the frontal lobe (left: 99 voxels, right: 101 voxels), the limbic
cortex (left: 31, right: 30), the occipital lobe (left: 56; right: 54), the
parietal lobe (left: 64; right: 62), and the temporal lobe (left: 34; right:
36). This classification was based on an automatic labelling of the voxels
locations given by the Talairach Daemon (http://ric.uthscsa.edu/projects/tdc),
after registering the DTI data into the MNI standard space, and subsequent
correction from MNI to Talairach space (see, e.g., [36]). Figure 8 shows the distribution of the connectivity index, in
the matrix form, between any two regions, arranged by group and by hemisphere.
The matrix shown in Figure 8 reveals an organization of the connectivity index
that follows the anatomical organization of the brain regions regarding their
locations. Since the connectivity index encompasses the anisotropy factor, its
value highly depends on which regions we are connecting, which means which
global pathways the geodesics are close to.
First, the diagonal blocks of the matrix show clearly
a lower level of connectivity than the extradiagonal blocks. This seems to
indicate that the connectivity index penalises short fibers, and inversely
favors long fibers, especially interhemispheric fibers. Secondly, the blocks
that show the highest connectivity index are the blocks that connect the right
and left occipital lobes.
This result is not surprising since the fiber tracts
that connect right and left occipital lobes follow a trajectory through the
splenium of the corpus callosum (forceps major), which is a highly anisotropic
area.
Geodesics
We further
investigated which of the constructed geodesics may represent actual fiber
trajectories. To approach this question, we thresholded the connectivity matrix
in order to emphasize the geodesics with the largest connectivity indices.
Specifically, we considered the geodesics with the highest connectivity
indices for each interhemispheric block connecting symmetrical groups, taken
independently. Figure 9 represents each group of geodesics in different
colors. The most probable geodesics paths follow the principal long association
fasciculi. The frontal lobe is connected to the occipital lobe via the
frontooccipital fasciculus. The temporal lobe is connected to the occipital
via the inferior longitudinal fasciculus, and to the frontal lobe via the
uncinate fasciculus. All major long association tracts are represented by these
geodesics.
(a)
(b)
Geodesics versus Streamlining
Finally, in
order to compare the results of our method to a conventional fiber tracking
method, we performed a streamline tractography from the seed voxels,
with four tracts per voxel. As a stopping criterion, we chose a maximum step
angle of , and an anisotropy threshold of [19]. To
compare the results to GT, we selected the four geodesics, having the highest
probability index, for each voxel in the set of seed voxels. This way, we have
the same number of tracts using both methods ( tracts). Figure 10
shows the results of these two procedures. The
streamline method produces many incomplete tracts, especially association
tracts, while the proposed GT method succeeded in reconstructing the major
association and commissural tracts, including the uncinate, the inferior
frontotemporal, and the callosal fibers. Note that the frontooccipital tract
is not present at this level of threshold (we only considered four geodesics
per voxel).
(a)
(b)
4. Discussion
Global optimisation is a valuable strategy in the context of path planning. When one has the information of where to start and where to go, this information is used to overcome local poor optimality. In the context of white matter diffusionbased tractography, where we often have strong hypotheses about the localisation of the regions in the brain, global optimisation can overcome some serious weaknesses of the process. Mainly, uncertainty about local fibre orientation, reflecting partial volume effects caused by crossing fibres, or local low signal to noise, can be handled efficiently using GT.
We have presented here a method to perform such globalbased path planning in an anisotropic medium. The method is very robust to high anisotropy, and provides an extremely accurate numerical solution to the Eikonal equation.
On realdata experiments, the reconstructed geodesics that have a high connectivity index correspond to known fiber tract fasciculi connecting the cortex. These fasciculi can all be retrieved by other tractography methods that use DTI data, providing priors on their location using one or more regions of interest [37, 38], especially intermediate regions located in white matter. GT automatically depicted these fasciculi with no prior.
However, the Ushaped fibers, that is, the short association tracts, are not favored by our connectivity index. This can be easily seen by looking at the diagonal blocks of the matrix in Figure 8. The long association tracts, as well as the commissural fibers, are more present with a higher connectivity index.
GT also allows one to construct interhemispheric tracts between each pair of regions located in different hemispheres. These tracts include homotopic and heterotopic connexions, that is, tracts connecting, respectively, symmetrical and asymmetrical regions lying in different hemispheres. It is worth noting that standard tractography methods usually fail to recover most callosal connexions, apart from the medial ones. This is a good illustration of the problem of crossing fibres, as those connexions cross the superior longitudinal fasciculus. However, recent probabilistic tractography with more complex local models has successfully traced those types of connexions [14, 20, 21].
There is an intuitive relationship between geodesic, for the inverse tensor metric, and probabilistic tractographies. Probabilistic tractography consists of constructing a distribution of connexions, by sampling tracts using local orientation distributions. In the basic case where this local probability model for fibre orientations is defined using the tensor model (i.e., a Gaussian local model with a covariance matrix proportional to the diffusion tensor ), the probability of a tract following an orientation given by at a location writes then, for some pathway connecting to , and for some discretisation of this pathway, the probability of moving along is the product of the infinitesimal step probabilities: Maximising this probability could then be related to minimising the geodesic distance, relative to the inverse tensor metric. While the probabilistic method gives a distribution of connexions, GT gives the mode of this distribution, that is, the path with highest probability. Note also that the probabilistic model given by (12) can be improved to fit the data more accurately (e.g., multiple tensors, etc.), which can be seen as a change in the metric tensor in GT.
Using GT, it is possible to study the organisation of large brain networks in terms of their anatomical connexions. Such networks have been studied in terms of structural invariants in a graph theoretical framework by several authors [39–41]. These works have been conducted for studying the structural organisation of the cat or primate brain, as well as for the human functional brain organisation, but have never been applied to large human anatomical networks, because no method has been proposed to construct such networks. GT could provide this structural information, via a graph that has been thresholded or not, since the connectivity index in itself contains information about the connectional structure.
There are two major issues when using geodesics for the tractography. First, choosing a metric for which geodesics represent fibre pathway trajectories is not straightforward. The correct metric might show more anisotropy than the diffusion tensor, as discussed earlier. Also, the choice of the metric might depend on the white matter fibres under investigation. The second issue is that, for any pair of regions in the brain, there exists a geodesic between those regions. However, this is not true for white matter fibres. One then has to decide when a geodesic is a fibre trajectory, for example, by defining indices and performing statistical thresholding under some null hypothesis. This problem of thresholding tractography results is not specific to GT, but is met by any other tractography method. It is though a bigger problem in the case of GT because every pair of regions is potentially connected. Another problem with GT is that, in the presence of two separate connexions between two regions, we are only able to detect one of them (the shortest one in terms of geodesic distance).
One way to validate GT results would be by comparison with another measure of connectivity. For example, measures of functional connectivity using functional magnetic resonance imaging (fMRI) by means of correlations [42] or partial correlations [43] are thought to be closely linked to the anatomical structure sustaining the brain regions, seen as graph nodes. The GT technique provides a unique tool for performing a comparison between anatomical and functional connectivity, since it can apply to large networks, and provide a measure of anatomical connectivity between each pair of nodes of the brain network. It can readily be used to compare the architectures of brain networks that have been studied in humans from the functional perspective (e.g., Salvador et al. [44] used partial correlations of fMRI data on a set of 100 regions), or using voxelbased morphometry to correlate cortical thickness between different cortical areas (e.g., He et al. [45] used this technique to study 100 cortical areas in humans). Such investigations have considerable possible applications, both cognitive and clinical. On the one hand, this method could serve as a basis for comparing anatomical and functional connectivities, as said earlier, and could help to understand how the brain works as an evolving network. On the other hand, the structure of restricted networks has already helped to distinguish between healthy subjects and patients, for example, Alzheimer disease in the case of functional connectivity [46], and Schizophrenia in the case of white matter morphology [47]. The GT method could serve for the characterisation of the structural organisation of those brain networks in terms of their connectional fingerprints.
Appendices
A. Algorithms
Algorithm 1. Fast marching algorithm
Definition 1s. Let Known be the set of points whose
value has been computed and will not change. Let Trial be the set of voxels
that are being examined (26neighbourhood of Known), and let Far be the set of
voxels that have not been examined yet. Finally, if is a set of voxels, let denote the
number of voxels that belong to .
(i)Initialization:
(a)move to Known and
set ,(b)move to
Far every such that and set ,(c)update
in the
neighbourhood of using Algorithm
2,
(ii)While :
(a)search
for the voxel in Trial with
the smallest value of ,(b)move to Known,(c)update
in the
neighbourhood of using Algorithm 2.
Algorithm 2. Updating procedure for the distance function
at voxel
:
(i) for all in the
26neighborhood of and ,
B. Explicit Solution for the Updating Procedure
Here we provide an explicit solution for the minimisation problem formulated in (8). Recall that the problem was to find the minimum, inside a simplex, for the following expression: In order to simplify the notations, and considering that , we will use the following: The function depends simply on and : is differentiable and convex, it is then minimal when . When constraining the minimum to lie inside the simplex , the solution is either that for which the gradient is zero, if it lies inside , or it is on the edges of if the gradient is zero outside the simplex. In the latter case, the minimisation problem is 1D, and the solution simplifies greatly.
First, let us write the unconstrained solution: implies This equation means that the minimum of lies on the straight line defined by the equation . This simplifies the problem, as the problem is again 1D if we replace by the function which expression depends on the values of and .
(i)If , (ii)If , In the last case, the problem reduces to minimising a 1D function of the form , in which case the solution writes where and .
Finally, if the solution given by the above lies outside the simplex (i.e., ), then we minimise on the edges of the simplex, which is again a 1D problem. This is equivalent to setting one of the to zero, and keeping the results which minimises :
(i): (ii): (iii):Acknowledgments
The authors would like to acknowledge funding from the Association pour la Recherche contre le Cancer and the Dr. Hadwen Trust For Humane Research (SJ). They are extremely grateful to Dr. Thien Huong N'Guyen, Service de Neuroradiologie (CHNO des QuinzeVingts, Paris), for kindly providing the diffusion data.
References
 J. W. Scannell, G. A. P. C. Burns, C. C. Hilgetag, M. A. O'Neil, and M. P. Young, “The connectional organization of the corticothalamic system of the cat,” Cerebral Cortex, vol. 9, no. 3, pp. 277–299, 1999. View at: Publisher Site  Google Scholar
 D. J. Felleman and D. C. van Essen, “Distributed hierarchical processing in the primate cerebral cortex,” Cerebral Cortex, vol. 1, no. 1, pp. 1–47, 1991. View at: Publisher Site  Google Scholar
 F. Crick and E. Jones, “Backwardness of human neuroanatomy,” Nature, vol. 361, no. 6408, pp. 109–110, 1993. View at: Publisher Site  Google Scholar
 M. Catani, R. J. Howard, S. Pajevic, and D. K. Jones, “Virtual in vivo interactive dissection of white matter fasciculi in the human brain,” NeuroImage, vol. 17, no. 1, pp. 77–94, 2002. View at: Publisher Site  Google Scholar
 S. Wakana, H. Jiang, L. M. NagaePoetscher, P. C. M. van Zijl, and S. Mori, “Fiber tractbased atlas of human white matter anatomy,” Radiology, vol. 230, no. 1, pp. 77–87, 2004. View at: Publisher Site  Google Scholar
 T. E. J. Behrens, H. JohansenBerg, M. W. Woolrich et al., “Noninvasive mapping of connections between human thalamus and cortex using diffusion imaging,” Nature Neuroscience, vol. 6, no. 7, pp. 750–757, 2003. View at: Publisher Site  Google Scholar
 P. L. Croxson, H. JohansenBerg, T. E. J. Behrens et al., “Quantitative investigation of connections of the prefrontal cortex in the human and macaque using probabilistic diffusion tractography,” Journal of Neuroscience, vol. 25, no. 39, pp. 8854–8866, 2005. View at: Publisher Site  Google Scholar
 M. F. S. Rushworth, T. E. J. Behrens, and H. JohansenBerg, “Connection patterns distinguish 3 regions of human parietal cortex,” Cerebral Cortex, vol. 16, no. 10, pp. 1418–1430, 2006. View at: Publisher Site  Google Scholar
 H. JohansenBerg, T. E. J. Behrens, M. D. Robson et al., “Changes in connectivity profiles define functionally distinct regions in human medial frontal cortex,” Proceedings of the National Academy of Sciences of the United States of America, vol. 101, no. 36, pp. 13335–13340, 2004. View at: Publisher Site  Google Scholar
 H. W. Powell, G. J. Parker, D. C. Alexander et al., “Hemispheric asymmetries in languagerelated pathways: a combined functional MRI and tractography study,” NeuroImage, vol. 32, no. 1, pp. 388–399, 2006. View at: Publisher Site  Google Scholar
 P. J. Basser, J. Mattiello, and D. LeBihan, “Estimation of the effective selfdiffusion tensor from the NMR spin echo,” Journal of Magnetic Resonance, vol. 103, no. 3, pp. 247–254, 1994. View at: Publisher Site  Google Scholar
 K. M. Jansons and D. C. Alexander, “Persistent angular structure: new insights from diffusion magnetic resonance imaging data,” Inverse Problems, vol. 19, no. 5, pp. 1031–1046, 2003. View at: Publisher Site  Google Scholar
 T. E. J. Behrens, M. W. Woolrich, M. Jenkinson et al., “Characterization and propagation of uncertainty in diffusionweighted MR imaging,” Magnetic Resonance in Medicine, vol. 50, no. 5, pp. 1077–1088, 2003. View at: Publisher Site  Google Scholar
 T. Hosey, G. Williams, and R. Ansorge, “Inference of multiple fiber orientations in high angular resolution diffusion imaging,” Magnetic Resonance in Medicine, vol. 54, no. 6, pp. 1480–1489, 2005. View at: Publisher Site  Google Scholar
 D. S. Tuch, T. G. Reese, M. R. Wiegell, N. Makris, J. W. Belliveau, and V. J. Wedeen, “High angular resolution diffusion imaging reveals intravoxel white matter fiber heterogeneity,” Magnetic Resonance in Medicine, vol. 48, no. 4, pp. 577–582, 2002. View at: Publisher Site  Google Scholar
 M. Mori, B. J. Crain, V. P. Chacko, and P. C. M. van Zijl, “Threedimensional tracking of axonal projections in the brain by magnetic resonance imaging,” Annals of Neurology, vol. 45, no. 2, pp. 265–269, 1999. View at: Publisher Site  Google Scholar
 D. K. Jones, A. Simmons, S. C. R. Williams, and M. A. Horsfield, “Noninvasive assessment of axonal fiber connectivity in the human brain via diffusion tensor MRI,” Magnetic Resonance in Medicine, vol. 42, no. 1, pp. 37–41, 1999. View at: Publisher Site  Google Scholar
 T. E. Conturo, N. F. Lori, T. S. Cull et al., “Tracking neuronal fiber pathways in the living human brain,” Proceedings of the National Academy of Sciences of the United States of America, vol. 96, no. 18, pp. 10422–10427, 1999. View at: Publisher Site  Google Scholar
 P. J. Basser, S. Pajevic, C. Pierpaoli, J. Duda, and A. Aldroubi, “In vivo fiber tractography using DTMRI data,” Magnetic Resonance in Medicine, vol. 44, no. 4, pp. 625–632, 2000. View at: Publisher Site  Google Scholar
 G. J. M. Parker and D. C. Alexander, “Probabilistic anatomical connectivity derived from the microscopic persistent angular structure of cerebral tissue,” Philosophical Transactions of the Royal Society of London, Series B, Biological Sciences, vol. 360, no. 1457, pp. 893–902, 2005. View at: Publisher Site  Google Scholar
 T. E. J. Behrens, H. J. Berg, S. Jbabdi, M. F. S. Rushworth, and M. W. Woolrich, “Probabilistic diffusion tractography with multiple fibre orientations: what can we gain?” NeuroImage, vol. 34, no. 1, pp. 144–155, 2007. View at: Publisher Site  Google Scholar
 G. J. M. Parker, C. A. M. WheelerKingshott, and G. J. Barker, “Estimating distributed anatomical connectivity using fast marching methods and diffusion tensor imaging,” IEEE Transactions on Medical Imaging, vol. 21, no. 5, pp. 505–512, 2002. View at: Publisher Site  Google Scholar
 D. L. Chopp, “Replacing iterative algorithms with singlepass algorithms,” Proceedings of the National Academy of Sciences of the United States of America, vol. 98, no. 20, pp. 10992–10993, 2001. View at: Publisher Site  Google Scholar
 S. Jbabdi, P. Bellec, G. Marrelec, V. Perlbarg, and H. Benali, “A level set method for building anatomical connectivity paths between brain areas using DTI,” in Proceedings of the 2nd IEEE International Symposium on Biomedical Imaging: Macro to Nano (ISBI '04) , vol. 1, pp. 1024–1027, Arlington, Va, USA, April 2004. View at: Google Scholar
 C. Lenglet, R. Deriche, and O. Faugeras, “Inferring white matter geometry from diffusion tensor MRI: application to connectivity,” in Proceedings of the 8th European Conference on Computer Vision (ECCV '04), T. Pajdla and J. Matas, Eds., vol. 3024 of Lecture Notes in Computer Science, pp. 127–140, Prague, Czech Republic, April 2004. View at: Google Scholar
 M. Jackowski, C. Y. Kao, M. Qiu, R. T. Constable, and L. H. Staib, “White matter tractography by anisotropic wavefront evolution and diffusion tensor imaging,” Medical Image Analysis, vol. 9, no. 5, pp. 427–440, 2005. View at: Publisher Site  Google Scholar
 P. Staempfli, T. Jaermann, G. R. Crelier, S. Kollias, A. Valavanis, and P. Boesiger, “Resolving fiber crossing using advanced fast marching tractography based on diffusion tensor imaging,” NeuroImage, vol. 30, no. 1, pp. 110–120, 2006. View at: Publisher Site  Google Scholar
 J. N. Tsitsiklis, “Efficient algorithms for globally optimal trajectories,” IEEE Transactions on Automatic Control, vol. 40, no. 9, pp. 1528–1538, 1995. View at: Publisher Site  Google Scholar
 L. Qingfen, Enhancement, extraction, and visualization of 3D volume data, Ph.D. thesis, Linkóping University, Linkóping, Sweden, 2003.
 J. A. Sethian and A. Vladimirsky, “Ordered upwind methods for static HamiltonJacobi equations,” Proceedings of the National Academy of Sciences of the United States of America, vol. 98, no. 20, pp. 11069–11074, 2001. View at: Publisher Site  Google Scholar
 V. I. Smirnov, A Course on Higher Mathematics, vol. 4, Pergamon Press, New York, NY, USA, 1964.
 J. A. Sethian, Level Set Methods and Fast Marching Methods, Cambridge University Press, Cambridge, Mass, USA, 2002.
 E. W. Dijkstra, “A note on two problems in connexion with graphs,” Numerische Mathematik, vol. 1, no. 1, pp. 269–271, 1959. View at: Publisher Site  Google Scholar
 L. O'Donnell, S. Haker, and C.F. Westin, “New approaches to estimation of white matter connectivity in diffusion tensor MRI: elliptic PDEs and geodesics in a tensorwarped space,” in Proceedings of the 5th International Conference on Medical Image Computing and ComputerAssisted Intervention (MICCAI '02), T. Dohi and R. Kikinis, Eds., pp. 459–466, Tokyo, Japan, September 2002. View at: Google Scholar
 C.F. Westin, S. E. Maier, H. Mamata, A. Nabavi, F. A. Jolesz, and R. Kikinis, “Processing and visualization for diffusion tensor MRI,” Medical Image Analysis, vol. 6, no. 2, pp. 93–108, 2002. View at: Publisher Site  Google Scholar
 M. Brett, I. S. Johnsrude, and A. M. Owen, “The problem of functional localization in the human brain,” Nature Reviews Neuroscience, vol. 3, no. 3, pp. 243–249, 2002. View at: Publisher Site  Google Scholar
 S. Mori, K. Frederiksen, P. C. M. van Zijl et al., “Brain white matter anatomy of tumor patients evaluated with diffusion tensor imaging,” Annals of Neurology, vol. 51, no. 3, pp. 377–380, 2002. View at: Publisher Site  Google Scholar
 S. Mori, S. Wakana, L. M. NagaePoetscher, and P. C. M. van Zijl, MRI Atlas of Human White Matter, Elsevier Science, Amsterdam, The Netherlands, 2005.
 G. Tononi, A. R. McIntosh, D. P. Russell, and G. M. Edelman, “Functional clustering: identifying strongly interactive brain regions in neuroimaging data,” NeuroImage, vol. 7, no. 2, pp. 133–149, 1998. View at: Publisher Site  Google Scholar
 O. Sporns, G. Tononi, and G. M. Edelman, “Theoretical neuroanatomy and the connectivity of the cerebral cortex,” Behavioural Brain Research, vol. 135, no. 12, pp. 69–74, 2002. View at: Publisher Site  Google Scholar
 O. Sporns, D. R. Chialvo, M. Kaiser, and C. C. Hilgetag, “Organization, development and function of complex brain networks,” Trends in Cognitive Sciences, vol. 8, no. 9, pp. 418–425, 2004. View at: Publisher Site  Google Scholar
 M. McIntosh and F. GonzalezLima, “Structural equation modelling and its application to network analysis of functional brain imaging,” Human Brain Mapping, vol. 2, no. 12, pp. 2–22, 1994. View at: Google Scholar
 G. Marrelec, J. Daunizeau, M. PélégriniIssac, J. Doyon, and H. Benali, “Conditional correlation as a first step toward common framework for functional brain interactivity modeling in functional MRI and MEG/EEG,” IEEE Transactions on Signal Processing, vol. 53, no. 9, pp. 3503–3516, 2005. View at: Publisher Site  Google Scholar
 R. Salvador, J. Suckling, M. R. Coleman, J. D. Pickard, D. Menon, and E. Bullmore, “Neurophysiological architecture of functional magnetic resonance images of human brain,” Cerebral Cortex, vol. 15, no. 9, pp. 1332–2342, 2005. View at: Publisher Site  Google Scholar
 Y. He, Z. J. Chen, and A. C. Evans, “Smallworld anatomical networks in the human brain revealed by cortical thickness from MRI,” Cerebral Cortex, vol. 17, no. 10, pp. 2407–2419, 2007. View at: Publisher Site  Google Scholar
 M. D. Greicius, G. Srivastava, A. L. Reiss, and V. Menon, “Defaultmode network activity distinguishes Alzheimer's disease from healthy aging: evidence from functional MRI,” Proceedings of the National Academy of Sciences of the United States of America, vol. 101, no. 13, pp. 4637–4642, 2004. View at: Publisher Site  Google Scholar
 R. A. A. Kanaan, J.S. Kim, W. E. Kaufmann, G. D. Pearlson, G. J. Barker, and P. K. McGuire, “Diffusion tensor imaging in schizophrenia,” Biological Psychiatry, vol. 58, no. 12, pp. 921–929, 2005. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2008 S. Jbabdi et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.