Quantitative evaluation of 10 tractography algorithms on a realistic diffusion MR phantom.

TL;DR: A common dataset with known ground truth and a reproducible methodology to quantitatively evaluate the performance of various diffusion models and tractography algorithms is used and evidence that diffusion models such as (fiber) orientation distribution functions correctly model the underlying fiber distribution is provided.

About: This article is published in NeuroImage.The article was published on 2011-05-01 and is currently open access. It has received 410 citations till now. The article focuses on the topics: Diffusion MRI & Tractography.

Summary

1 Introduction

• The unique ability of Diffusion Weighted MRI (DW-MRI) based fiber tractography to map, in vivo, the architecture of white matter pathways has ignited strong interest in clinical and neuroscience research.
• Notably, one objective of these methods is to model and track in the presence of complex fiber configurations such as crossings or kissings, but objective comparison of the performance of each has not been done yet.
• The increasing number of diffusion models and tractography algorithms is both a blessing and a curse: Unfortunately, the fiber configuration was sparse as compared with the brain.
• Results were analyzed and ranked based on several metrics.

2 Material and Methods

• The authors will first briefly review the construction of a realistic diffusion MR ground truth dataset.
• The authors will then detail the rules of the tractography competition and the methodology developed to quantitatively evaluate and compare contributions.

2.1.1 Design of a Realistic Diffusion MR Phantom

• The construction of diffusion phantoms is a challenging task involving the following two steps [Poupon et al., 2008]: Design of a realistic and practically feasible fiber configuration.
• An adequate recipe to fill the container with a MRI compatible solution is also required.
• The configuration should be as realistic possible - containing crossing and kissing fibers as well as bundles of different curvatures.
• A polyurethane negative and positive prints of the target bundles were manufactured and used to strongly tighten the fibers together.
• First, a layer of bundles was placed everywhere in the phantom.

2.1.2 Diffusion-Weighted MRI Acquisitions

• Diffusion-weighted data of the phantom were acquired on the 3T Tim Trio MRI systems of the NeuroSpin centre, equipped with a whole body gradient coil (40 mT/m, 200 T/m/s), and using a 12-channel receive only head coil, in combination with the whole body transmit coil of the MRI system.
• A single-shot diffusion-weighted twice refocused spin echo echoplanar pulse sequence was used to perform the acquisitions, while compensating for the first order Eddy currents.
• Two datasets were acquired at two different spatial resolutions: 3 mm isotropic and 6 mm isotropic.
• Parameters for the 3 mm isotropic acquisition were as follows: field of view FOV=19.2cm, matrix 64x64, slice thickness TH=3mm, read bandwidth RBW=1775 Hz/pixel, partial Fourier factor 6/8, parallel reduction factor GRAPPA=2, repetition time TR=5s, 2 repetitions.
• The diffusion sensitization was applied along a set of 64 orientations uniformly distributed over the sphere.

2.1.3 Estimation of a Ground Truth Dataset

• To facilitate the evaluation of the different results submitted during the contest, the authors chose to restrict the analysis to a set of 16 fibers traversing 16 manually identified voxels, or seeds.
• The authors chose to 2We define SNR as the ratio between the signal magnitude and the noise power (i.e., standard deviation) [Kaufman et al., 1989].
• Competitors were asked to return a single representative fiber of the bundle traversing each seed voxel.
• Only S9 is ambiguous since two solutions are possible, but this ambiguity was detected soon after the contest started and could not be removed.
• Then, lines were smoothed using approximating cubic b-spline to remove any sampling noise, and fibers were resampled using 1000 uniformly distributed points, which formed the ground truth (Fig. 4 c).

2.2 Contest Rules

• Participants were free to use any combination of tools and algorithms that lead to the best result.
• Participants were free to use any, or several, of these datasets.
• The 64 diffusion gradients used during acquisition were provided as a text file.
• Each fiber has to come in a separate file named with the seed label it originates.
• This design favors deterministic tracking algorithms and this issue will be discussed later in the conclusions.

2.2.1 Common Fiber File Format

• Due to the increasing availability of tractography softwares (DtiStudio, Brainvisa, TrackVis, MedINRIA, Slicer to quote just a few), and by extension to the existence of numerous fiber file formats, one could not reasonably rely on one of them, mainly because they can be quite complex to produce, especially for those who are not familiar with them.
• Instead, the authors chose to rely on the simplest existing format: the text file.
• Participants were asked to return a single text file per fiber, where the fiber coordinates are listed in sequential order (i.e., x y z coordinates of the first point, x y z coordinates of the second point, etc.), one point per line.
• Thus, the number of lines corresponds exactly to the number of points of a fiber.

2.2.3 Pre-processing

• Another important issue to take care of when evaluating tractography results from different participants is the fiber sampling.
• The sampling is very likely to differ from one submission to another: some methods produce highly sampled fibers with several hundreds or thousands of points, while others only provide a dozen of points.
• To normalize this, fibers were parametrized by interpolating cubic b-splines.
• Interpolation was chosen in order not to alter the fiber coordinates as returned by the 12 participants.
• In the next section, the authors present the evaluation methodology used to compare tractography results with the ground truth.

2.3 Evaluation Methodology

• Evaluation was performed on a per-fiber basis.
• The authors recall that each participant had to return a dataset composed of 16 candidate fibers matching 16 ground truth fibers.
• Thus, the candidate fiber passing through seed N can be compared to the ground truth fiber going through the same seed.
• The evaluation methodology narrows down to the evaluation of differences between pairs of curves.
• In the following, the authors describe the evaluation measures of curve matching they used for this contest.

2.3.1 Generic Score of Fiber Match

• The optimal result is realized when the candidate fiber perfectly matches the ground truth, i.e., when both fibers are superimposed.
• The choice of c is obviously not unique, and without any prior knowledge there is no best choice for it.
• Then, following [Fillard et al., 2007], dynamic programming is used to determine the path of minimal cost within this distance matrix, which gives us the final correspondences between the arc length of both fibers.
• On the contrary, by taking the angular difference between tangents, the sRMSE will be low (resp. high) when fibers are parallel (resp. orthogonal).
• In the following, the authors express three metrics that were used for the contest: the spatial metric, the tangent metric and the curve metric.

2.3.2 The Contest Metrics

• The spatial metric is simply the L2 norm between two corresponding fiber positions.
• (3) The sRMSE endowed with the spatial metric is expressed in mm and ranges from 0 (overlapping points) to infinity.
• V1 and v2 are normalized tangent vectors to fiber points.
• The curve metric is expressed as the absolute difference of the curvature between two fiber points: dist(κ1, κ2) = |κ2 − κ1| (5) The sRMSE endowed with the curve metric is expressed in mm−1 and ranges from 0 to infinity.
• In the next section, the authors present the results of the qualitative and quantitative evaluation of the 10 contributions received during the Fiber Cup.

3 Results

• A total of 9 individual submissions were received, including one with 2 results, making a total of 10 tractography results.
• Results were analyzed following the methodology described in the previous section.
• Computation of quantitative metrics was performed on a regular PC (Intel Core 2 Duo, 2Gb of memory).
• For the sake of completeness, the authors also included the result of the probabilistic tractography algorithm implemented in FSL [Behrens et al., 2007a], as this is one of most widely used algorithm within the neuroscience community.
• The authors first summarize the 10 contributions in terms of diffusion model and tractography algorithm chosen.

3.1 Summary of Contributions

• An overview of the 10 tractography methods evaluated during the contest is given in Table 1.
• More precisely, the choice of the diffusion model appeared as more important than the tractography algorithm itself, often reduced to a streamline approach, although some variety can be noted.
• The choice of the higher resolution / reduced SNR dataset is interesting since a common problem in real acquisitions is to know how much of the SNR should be sacrificed in favor of the spatial resolution.
• Streamline with propagation direction following tensor PDD 3 × 3 × 3, b = 1500 Tracking is performed from all phantom voxels and only fibers going through seeds were kept.
• A scoring function determines the most likely fiber to return a single fiber / seed.

3.2 Qualitative Evaluation

• The authors present on Figure 5 an overlap of the 10 contributions for each ground truth fiber (one image corresponding to one seed location), and on Figure 6 the individual results for each contribution (one image corresponds to the 16 candidate fibers of one method).
• Note that the image number does correspond to the method Id of Table 1.
• From Figure 5, the authors can conclude that, except for S13 and S14 that are located on the isolated U-shape structure (Region 7, Fig. 1 right), at least one contribution per seed fails at reconstructing the correct pathway.
• Very often the algorithm chose the wrong direction when going through crossing regions).
• Besides crossings, most fibers appear as nicely reconstructed.

3.3 Quantitative Evaluation

• One sRMSE comprises the evaluation of a metric on two times 1000 points, the authors end up with a total of 960 000 point-to-point metrics being tested.
• The tangent metric evaluates whether the fiber trajectory correctly follows the ground truth.
• Note that method 2 obtains bad scores for S3 and S14.
• Indeed, from Fig. 6 (3), one can notice that S3 has a nonsmooth trajectory around branching 2 (Fig. 1 right) that appears like an inflection point.
• Other methods, in particular method 10, are penalized by this metric due to the high frequency noise of their fibers.

3.4 Ranking

• Tractography results were ranked according to the following rule.
• For each fiber and each metric, the method realizing the best score (i.e., the lowest metric value) was attributed 3 points.
• Obviously, improvements are possible since one may not desire to give the same importance to all metrics.
• To further illustrate the performance of the tested algorithms in real situations, the authors performed tractography on a brain dataset with the top two methods (methods 7 and 2) and compared them to a single-DT streamline tractography algorithm (method 4).
• Results are presented in Supplementary Section 3.

4.1 Comments on the Methods

• As expected, single tensor-based methods (Fig. 5 (3), (4) and (8)) seem to perform worse than others in crossing regions for the obvious reason that a single tensor is unable to correctly characterize the two-fiber compartment specific of those regions.
• In particular, in the lower crossing area (region 1, Fig. 1 right) methods 3 [Basser et al., 2000] and 4 [Lazar et al., 2003] chose to avoid it by contouring it, while method 8 [Fillard et al., 2003] stopped the tracking, very probably because the crossing yield a fiber curvature greater than the maximum angular deviation 22 authorized.
• Multi-tensor based approaches (Fig. 5 (1) and (5)) are clearly a big improvement compared to single-tensor methods.
• Method 10 estimates the ODF using a probability density constraint and a spatial regularity prior.
• Such a tractography technique is noise-sensitive and at the same time, highly dependent on accurate ODF estimation.

4.2 Recommendations

• It is still possible to make a few recommendations about methods which should be used and those which should be avoided in tractography.
• The recommendations that follow are based on the tested implementations of each method.
• Moreover, there is no guarantee that the results obtained on the phantom dataset can be directly transposed to real situations.
• Notably, the DT model with only few degrees of freedom is by essence less sensitive to noise than more complex models, which often makes it the unique alternative in clinical applications.
• Method 10 explicitly imposes a spatial regularity when estimating the ODF, which eventually leads to good fiber pathways even using a streamline tractography algorithm, which give some evidence that the fiber directions were correctly modeled by the ODF.

Procedure:

• FOD are estimated using the odfestimate matlab function as follows: odfestimate -input data -SHorder 4 -reconstructiontess 4 -output ODFSH.mat.
• Fibers were tracked using the fibertrack matlab function: fibertrack -seedposition [x,y,z] -eliminatemaxangle 45 -incominglocalanglethreshold 20 -possiblecurveanglethreshold 60 -possiblecurvemagthreshold 0.8 -nomaxmagthreshold 0.7 -output possiblefibertracts.mat.
• The angle threshold was set to 60 degrees.
• Finally, fibers were further refined using the calculatepathcrossedmaxtimes matlab function: calculatepathcrossedmaxtimes -input possiblefibertracts.mat -output finaltracts.mat 34.

Brain Dataset

• In order to illustrate the performance difference of the winning methods, the authors performed tractography on a real brain dataset with method 7 and 2 and compared them to a single-DT method (method 4).
• To present the results, two ROIs were manually defined and only fibers passing through the ROIs were retained.
• The first ROI delineates the corpus callosum (Figure 9 left): the bundle should not only contain the U-shaped cortico-cortical fibers, but also the longer projection fibers that pass several regions of crossing fibers in the brain (including the corona radiata).
• Fibers passing through the pons should go up to the motor cortex via the corticospinal tract and fan along the right and left hemisphere.
• Only Method 7 (and partly method 2) was able to detect the fanning structure of the cortico-spinal tract.

Figures

Figure 3: ADC and FA images of the 6x6x6 phantom dataset. Top: The single slice of the ADC image is shown for the three b-values used for acquisition. The values reported are the mean (standard deviation) of the ADC values of the phantom. Bottom: The single slice of the FA image is shown for the three b-values used for acquisition. The values reported are the mean (standard deviation) of the FA values of the phantom.

Figure 10: Examples of tractography on a real dataset of three methods. For each method, the entire white matter was reconstructed and only fibers passing trough the ROIs shown on Figure 9 are shown. Left: Fibers passing through the corpus callosum (Fig. 9 left). Right: Fibers passing through the pons (Fig. 9 right). Results are shown for: top: single-DT and streamline tractography (method 4), middle: FOD with constrained spherical deconvolution and streamline tractography (2nd best method), and bottom: for global tractography (winning method). The color codes the main fiber orientation: a left-right oriented fiber is colored in red, an infero-superior fiber is blue and an antero-posterior fiber is green.

Figure 8: Probabilistic tractography of the phantom dataset using FSL. The connectivity maps are shown on top of the b0 image. We refer the reader to Figure 4 for the seed locations. The color range goes from yellow (low probability of connexion) to red (high probability). The obtained CMs indicate that probabilistic tractography remains feasible on such synthetic dataset, which could be further used to evaluate this class of algorithm.

Figure 7: Quantitative comparison of the 10 contributions to the ground truth using the spatial metric of Eq. 3 (left), the tangent metric of Eq. 4 (middle), and the curve metric of Eq. 5 (right). Top: The metric scores for each method and each fiber. Units are in mm for the spatial metric, in degree for the tangent metric and in mm −1 for the curve metric. Bottom: Points attributed to each method. Color-coding is the following: dark blue: 0 point, light blue: 1 point, yellow: 2 points, and brown: 3 points. The total number of points for each method and each metric is indicated below the x-axis.

Table 3: Seed-dependent parameters for method 8.

Figure 2: ADC and FA images of the 3x3x3 phantom dataset. Top: The middle slice of the ADC image is shown for the three b-values used for acquisition. The values reported are the mean (standard deviation) of the ADC values of the phantom. Bottom: The middle slice of the FA image is shown for the three b-values used for acquisition. The values reported are the mean (standard deviation) of the FA values of the phantom.

Table 2: Final score ranking of the 10 contributions. See Table 1 for more information about each method.

Figure 1: Sketch of the phantom that mimics a coronal section of the human brain. Left: Fiber pathways are highlighted in colors. Arrows indicate the directions of the synthetic fiber bundles. Right: The various crossing, splitting and kissing fiber configurations are numbered for a fast text referencing. Note that angles between crossing fibers were carefully determined, although not used during evaluation.

Figure 5: Panel of images showing the reconstructed fiber of all contributions passing through each seed selected for the context. The name below the image indicates the seed fibers originate (a color-method correspondence table is given below). Such overview allows to understand the variability of the results: some seeds (such as S3 and S4) were more successful than others (such as S7) to be reconstructed.

Figure 4: The 16 seed voxels chosen for the contest. Top: Seeds defined on the 3× 3× 3 mm (a) and 6× 6× 6 mm (b) datasets. Bottom: Ground truth fibers for both datasets. For each figure, the b0 image is shown.

Figure 9: ROI used for tractography on a real dataset. Left: Corpus Callosum. Right: Pons. Top: Axial view. Middle: Coronal view. Bottom: Sagittal view. ROI were manually drawn on the FA map.

Figure 6: Each of the 10 contributions is shown individually. Image (11) represents the ground truth for visual comparison. Fibers are colored by the seed they originate (see Fig. 4 for more information on the seeds location). Visually, methods 2 and 7 seem to perform the best. We refer the reader to Table 1 for the method - Id correspondence.

Frequently asked questions

Q1. What are the contributions mentioned in the paper "Quantitative evaluation of 10 tractography algorithms on a realistic diffusion mr phantom" ?

In this work, the authors use a common dataset with known ground truth and a reproducible methodology to quantitatively evaluate the performance of various diffusion models and tractography algorithms. The results provide evidence that: 1. For high SNR datasets, diffusion models such as ( fiber ) orientation distribution functions correctly model the underlying fiber distribution and can be used in conjunction with streamline tractography, and 2. 

Q2. What are the future works mentioned in the paper "Quantitative evaluation of 10 tractography algorithms on a realistic diffusion mr phantom" ?

In the future, new evaluation criteria will be proposed. Another possibility is to evaluate whether the boundaries of a bundle are correctly reconstructed by measuring the spatial distance in-between two tracts delimiting the bundle. The authors believe that such a common dataset along with the methodology proposed here can serve as an evaluation basis for existing and new algorithms. New results can be submitted for evaluation by emailing them to fibercup09 @ gmail. 

Q3. What was the role of the diffusion model in motivating the project?

The potential of tractography to help map anatomical connections played a significant role in motivating an ambitious project to map the human ”connectome” 1. 

Q4. Why did the authors choose to use 16 spatial positions?

The choice of those 16 spatial positions was made to ensure that a single fiber bundle passes through each of them to avoid ambiguity on the result and to facilitate the evaluation. 

Q5. What are the popular tractography algorithms?

Among deterministic tractography algorithms, streamline algorithms were developed first [Mori et al., 1999b,Conturo et al., 1999,Basser et al., 2000], followed by more elaborated tensor deflection algorithms [Weinstein et al., 1999,Lazar et al., 2003] or more global approaches [Poupon et al., 2001, Mangin et al., 2002]. 

Q6. How many fibers were used to fill the phantom?

the positive and negative prints were squeezed while keeping fibers strongly tightened until the openings (i.e, where the fibers enter/leave the phantom) are exactly 1cm thick. 

Q7. How many fibers are captured in the phantom?

Compression is carefully controlled to make sure that fibers are captured in 1mm2 crosssection everywhere throughout the phantom. 

Q8. Why did the authors exclude probabilistic tractography from the evaluation panel?

the nature of the ground truth itself prevents the inclusion of probabilistic tractography algorithms into the evaluation panel, since those output gener-9ally connectivity maps (CM) and not fiber pathways. 

Q9. What are the popular probabilistic tractography methods?

Probabilistic tractography methods include DT-based algorithms [Parker et al., 2003,Behrens et al., 2003,Lazar and Alexander, 2005,Friman et al., 2006,RamirezManzanares and Rivera, 2006,Savadjiev et al., 2008,Koch et al., 2002,Zhang et al., 2009], calculation of geodesics in a DT-warped space [Lenglet, 2006, Jbabdi et al., 2004], and numerous HARDI-based methods [Parker and Alexander, 2005, Perrin et al., 2005a,Seunarine et al., 2006,Behrens et al., 2007b,Jbabdi et al., 2007,Savadjiev et al., 2008, Chao et al., 2007a, Seunarine et al., 2007, Haroon and Parker, 2007,Kaden et al., 2007,Jeurissen et al., 2010]. 

Q10. How many fibers are found in the phantom?

Since the authors know the number of fibers and the space they are captured in, the authors can deduce the density of fibers, which was close to 1900 fibers/mm2 everywhere, including in the crossings. 

Q11. What are the objectives of this study?

The objectives of this study are to provide a qualitative and quantitative comparison of several tractography methods on the same realistic dataset with known ground truth and to freely distribute this dataset along with the evaluation methodology so that new methods can be easily evaluated and compared to existing ones. 

Q12. What is the procedure for comparing two consecutive fibers?

This procedure ensures that the function c is monotonically increasing, i.e., if s1 >= s2, c(s1) >= c(s2), which13ensures that two consecutive points of a fiber are associated to two other consecutive points.