Introduction to Basic Tracking

Local fiber tracking is an approach used to model white matter fibers by creating streamlines from local directional information. The idea is as follows: if the local directionality of a tract/pathway segment is known, one can integrate along those directions to build a complete representation of that structure. Local fiber tracking is widely used in the field of diffusion MRI because it is simple and robust.

In order to perform local fiber tracking, three things are needed: 1) A method for getting directions from a diffusion data set. 2) A method for identifying different tissue types within the data set. 3) A set of seeds from which to begin tracking. This example shows how to combine the 3 parts described above to create a tractography reconstruction from a diffusion data set.

To begin, let’s load an example HARDI data set from Stanford. If you have not already downloaded this data set, the first time you run this example you will need to be connected to the internet and this dataset will be downloaded to your computer.

from dipy.data import read_stanford_labels

hardi_img, gtab, labels_img = read_stanford_labels()
data = hardi_img.get_data()
labels = labels_img.get_data()
affine = hardi_img.affine

This dataset provides a label map in which all white matter tissues are labeled either 1 or 2. Lets create a white matter mask to restrict tracking to the white matter.

white_matter = (labels == 1) | (labels == 2)

1. The first thing we need to begin fiber tracking is a way of getting directions from this diffusion data set. In order to do that, we can fit the data to a Constant Solid Angle ODF Model. This model will estimate the Orientation Distribution Function (ODF) at each voxel. The ODF is the distribution of water diffusion as a function of direction. The peaks of an ODF are good estimates for the orientation of tract segments at a point in the image.

from dipy.reconst.shm import CsaOdfModel
from dipy.data import default_sphere
from dipy.direction import peaks_from_model

csa_model = CsaOdfModel(gtab, sh_order=6)
csa_peaks = peaks_from_model(csa_model, data, default_sphere,
                             relative_peak_threshold=.8,
                             min_separation_angle=45,
                             mask=white_matter)

2. Next we need some way of restricting the fiber tracking to areas with good directionality information. We’ve already created the white matter mask, but we can go a step further and restrict fiber tracking to those areas where the ODF shows significant restricted diffusion by thresholding on the general fractional anisotropy (GFA).

from dipy.tracking.local import ThresholdTissueClassifier

classifier = ThresholdTissueClassifier(csa_peaks.gfa, .25)

3. Before we can begin tracking is to specify where to “seed” (begin) the fiber tracking. Generally, the seeds chosen will depend on the pathways one is interested in modeling. In this example, we’ll use a \(2 \times 2 \times 2\) grid of seeds per voxel, in a sagittal slice of the corpus callosum. Tracking from this region will give us a model of the corpus callosum tract. This slice has label value 2 in the labels image.

from dipy.tracking import utils
import numpy as np

seed_mask = labels == 2
seeds = utils.seeds_from_mask(seed_mask, density=[2, 2, 2], affine=np.eye(4))

Finally, we can bring it all together using LocalTracking. We will then display the resulting streamlines using the dipy.viz package.

from dipy.tracking.local import LocalTracking
from dipy.viz import window, actor, colormap as cmap, have_fury
from dipy.tracking.streamline import Streamlines

# Enables/disables interactive visualization
interactive = False

# Initialization of LocalTracking. The computation happens in the next step.
streamlines_generator = LocalTracking(csa_peaks, classifier, seeds,
                                      affine=np.eye(4), step_size=.5)

# Generate streamlines object
streamlines = Streamlines(streamlines_generator)

# Prepare the display objects.
color = cmap.line_colors(streamlines)

if have_fury:
    streamlines_actor = actor.line(streamlines, cmap.line_colors(streamlines))

    # Create the 3D display.
    r = window.Renderer()
    r.add(streamlines_actor)

    # Save still images for this static example. Or for interactivity use
    window.record(r, n_frames=1, out_path='deterministic.png', size=(800, 800))
    if interactive:
        window.show(r)

Corpus Callosum Deterministic

We’ve created a deterministic set of streamlines, so called because if you repeat the fiber tracking (keeping all the inputs the same) you will get exactly the same set of streamlines. We can save the streamlines as a Trackvis file so it can be loaded into other software for visualization or further analysis.

from dipy.io.streamline import save_trk
save_trk("CSA_detr.trk", streamlines, affine,
         shape=labels.shape,
         vox_size=labels_img.header.get_zooms())

Next let’s try some probabilistic fiber tracking. For this, we’ll be using the Constrained Spherical Deconvolution (CSD) Model. This model represents each voxel in the data set as a collection of small white matter fibers with different orientations. The density of fibers along each orientation is known as the Fiber Orientation Distribution (FOD). In order to perform probabilistic fiber tracking, we pick a fiber from the FOD at random at each new location along the streamline. Note: one could use this model to perform deterministic fiber tracking by always tracking along the directions that have the most fibers.

Let’s begin probabilistic fiber tracking by fitting the data to the CSD model.

from dipy.reconst.csdeconv import (ConstrainedSphericalDeconvModel,
                                   auto_response)

response, ratio = auto_response(gtab, data, roi_radius=10, fa_thr=0.7)
csd_model = ConstrainedSphericalDeconvModel(gtab, response, sh_order=6)
csd_fit = csd_model.fit(data, mask=white_matter)

Next we’ll need to make a ProbabilisticDirectionGetter. Because the CSD model represents the FOD using the spherical harmonic basis, we can use the from_shcoeff method to create the direction getter. This direction getter will randomly sample directions from the FOD each time the tracking algorithm needs to take another step.

from dipy.direction import ProbabilisticDirectionGetter

prob_dg = ProbabilisticDirectionGetter.from_shcoeff(csd_fit.shm_coeff,
                                                    max_angle=30.,
                                                    sphere=default_sphere)

As with deterministic tracking, we’ll need to use a tissue classifier to restrict the tracking to the white matter of the brain. One might be tempted to use the GFA of the CSD FODs to build a tissue classifier, however the GFA values of these FODs don’t classify gray matter and white matter well. We will therefore use the GFA from the CSA model which we fit for the first section of this example. Alternatively, one could fit a TensorModel to the data and use the fractional anisotropy (FA) to build a tissue classifier.

classifier = ThresholdTissueClassifier(csa_peaks.gfa, .25)

Next we can pass this direction getter, along with the classifier and seeds, to LocalTracking to get a probabilistic model of the corpus callosum.

streamlines_generator = LocalTracking(prob_dg, classifier, seeds,
                                      affine=np.eye(4), step_size=.5,
                                      max_cross=1)

# Generate streamlines object.
streamlines = Streamlines(streamlines_generator)

if have_fury:
    streamlines_actor = actor.line(streamlines, cmap.line_colors(streamlines))

    # Create the 3D display.
    r = window.Renderer()
    r.add(streamlines_actor)

    # Save still images for this static example.
    window.record(r, n_frames=1, out_path='probabilistic.png', size=(800, 800))
    if interactive:
        window.show(r)

Corpus callosum probabilistic tracking.

save_trk("CSD_prob.trk", streamlines, affine,
         shape=labels.shape,
         vox_size=labels_img.header.get_zooms())