Basic Tracking with DIPY

Citation and Resources:¶

Tools included in this workflow¶

DIPY:

• Garyfallidis, E., Brett, M., Amirbekian, B., Rokem, A., van der Walt, S., Descoteaux, M., Nimmo-Smith, I., & Dipy Contributors (2014). Dipy, a library for the analysis of diffusion MRI data. Frontiers in neuroinformatics, 8, 8. Garyfallidis et al. (2014)

FURY - Free Unified Rendering in Python:

• Eleftherios Garyfallidis, Serge Koudoro, Javier Guaje, Marc-Alexandre Côté, Soham Biswas, David Reagan, Nasim Anousheh, Filipi Silva, Geoffrey Fox, and Fury Contributors. “FURY: advanced scientific visualization.” Journal of Open Source Software 6, no. 64 (2021): 3384. Garyfallidis et al. (2021)

• FURY Home

Dataset¶

• Stanford HARDI dataset, distributed with DIPY. Original acquisition described in: Rokem, A., Yeatman, J. D., Pestilli, F., Kay, K. N., Mezer, A., van der Walt, S., & Wandell, B. A. (2015). Evaluating the accuracy of diffusion MRI models in white matter. PLoS ONE, 10(4), e0123272. Rokem et al. (2015)

Educational resources¶

• Fiber Tracking Tutorials DIPY (v1.9.0)

Table of Contents¶

Introduction
Install and import Python libraries
DIPY command-line workflows
Data Preparation
Step 1: Diffusion model
Step 2: Stopping criterion
Step 3: Seeding
Step 4: Local Tracking: Probabilistic and Deterministic Approaches

Introduction¶

Local fiber tracking is a method for reconstructing white matter pathways from diffusion MRI data. It models the trajectories of fiber bundles by following local orientation information derived from the diffusion signal. The key idea is that if the dominant diffusion direction within each voxel is known, one can integrate along these directions to reconstruct streamlines representing possible neural tracts.

To perform local fiber tracking, four components are required:

1. A model to estimate local diffusion directions (e.g., DTI, CSA or CSD).

2. A stopping criterion that defines when a streamline should terminate (e.g., based on fractional anisotropy or tissue boundaries).

3. A set of seed points from which tracking begins.

4. Local tracking - the actual tractography execution.

In this notebook, we will combine these components to perform tractography using the DIPY library, demonstrating each step of the local tracking workflow.

Goal¶

Build practical understanding of local tractography workflows in DIPY, from diffusion modeling to streamline generation.

Learning Objectives¶

By the end of this notebook, you will be able to:

• Load and minimally preprocess diffusion MRI data in DIPY.

• Fit DTI, CSA and CSD models and visualize derived metrics (FA, FOD).

• Define stopping criteria and seed regions for tractography.

• Run both deterministic and probabilistic tractography.

• Compare outputs visually and quantitatively.

Install and import Python libraries¶

DIPY command-line workflows¶

This notebook installs DIPY 1.12.1. DIPY 1.12 introduced command-line entry points for several preprocessing and reconstruction tasks. See the DIPY Workflows Interfaces for the complete workflow catalog.

The Python examples in this notebook show how the individual steps work. The CLIs run the same workflows from shell scripts or batch jobs.

• dipy_brain_mask creates a brain mask with median Otsu, EVAC+, or SynthSeg. Median Otsu is the method that most closely matches this notebook’s DWI preprocessing.

• dipy_fit_fwdti fits the free-water elimination DTI model and writes its scalar maps. Use suitable multi-shell diffusion data for this model.

• dipy_fit_powermap creates an anisotropic powermap, which can serve as a pseudo-T1 image when an anatomical scan is unavailable.

Replace the placeholder paths in these command templates with files from your dataset:

dipy_brain_mask dwi.nii.gz --method median_otsu --bvalues_files dwi.bval --out_dir outputs/mask
dipy_fit_fwdti dwi.nii.gz dwi.bval dwi.bvec mask.nii.gz --out_dir outputs/fwdti
dipy_fit_powermap dwi.nii.gz dwi.bval dwi.bvec mask.nii.gz --out_dir outputs/powermap

The next cell runs --help for each command. This checks the installed entry points without starting a full reconstruction.

Data Preparation¶

Load Stanford HARDI dataset provided by DIPY¶

When we load the dataset, we get two key objects:

• img: the NIfTI image containing the 4D diffusion data (3D brain volume × diffusion-weighted directions).

• gtab (gradient table): encodes the acquisition scheme — specifically the b-values and gradient directions for each volume. B-values control the strength of diffusion weighting: b=0 means no diffusion weighting (a baseline image), while higher values (e.g., 2000 s/mm²) make the signal sensitive to water diffusion along the corresponding gradient direction.

We also define an affine matrix — a 4×4 transformation that maps between voxel indices and world (scanner) coordinates. Here we use an identity matrix so that tractography operates directly in voxel space.

This dataset provides a label map in which white matter tissues are labeled as either 1 (corpus callosum) or 2 (other white matter regions). Let’s create a white matter mask to restrict tracking to the white matter.

Minimal preprocessing¶

To keep things simple, we’ll compute a brain mask using median_otsu.

This function is inspired by MRtrix3’s bet which has default values median_radius=3, numpass=2. However, from tests on multiple 1.5T and 3T data from GE, Philips, and Siemens, the most robust choice is median_radius=4, numpass=4 (default).

Step 1: Diffusion model¶

To reconstruct white matter pathways, we first need to estimate the local fiber directions within each voxel. In diffusion MRI, water molecules diffuse more easily along axons than across them — this directional dependence can be modeled to infer underlying fiber orientations.

Several models can be used to represent this directional information. We will look into the following three models (for a complete overview of reconstruction models available in DIPY, see here):

A) DTI (Diffusion Tensor Imaging)

• Estimates a single principal diffusion direction per voxel.

• Computes scalar measures such as Fractional Anisotropy (FA) and Mean Diffusivity (MD).

• Simple and fast, but cannot resolve crossing fibers, which are present in a substantial portion of white matter voxels.

While FA is a scalar measure of anisotropy, you can also visualize a tensor ODF derived from the fitted tensor. This ODF is Gaussian-shaped and symmetric, representing the diffusion profile of a single tensor. It provides intuition about fiber orientation but is limited in regions with crossing or branching fibers.

B) CSA (Constant Solid Angle) ODFs:

• Estimates an Orientation Distribution Function (ODF) per voxel from the diffusion-weighted signal.

• Peaks of the ODF indicate likely fiber orientations.

• Can represent multiple orientations, but peaks may be less sharp than CSD, potentially leading to lower tracking precision.

The CSA model helps illustrate the transition from single-tensor representations to full ODFs, showing how we move toward richer, multi-directional voxel-wise information for tractography.

C) CSD (Constrained Spherical Deconvolution):

• Estimates the Fiber Orientation Distribution (FOD) by deconvolving the diffusion signal with a single-fiber response function. The FOD is represented using spherical harmonics (SH) — a mathematical basis for representing functions on a sphere, analogous to Fourier series on a circle.

• Resolves multiple distinct fiber populations per voxel (crossings, branchings).

• Produces sharper, more accurate orientation peaks, improving both deterministic and probabilistic tractography.

ODF vs FOD: An ODF describes the probability of diffusion in each direction, while a Fiber Orientation Distribution (FOD) specifically represents the underlying fiber orientations after deconvolving the single-fiber response from the signal. In practice, CSD produces FODs with sharper peaks than raw ODFs.

In short: while DTI and CSA provide useful but simplified diffusion representations, CSD offers a more precise reconstruction of fiber orientations, especially in regions with complex microstructure.

A) Reconstruction of the diffusion signal with DTI (single tensor) model¶

Let’s start with DTI, the simplest diffusion model, before exploring more complex representations with CSA and CSD.

In DTI, the diffusion profile within each voxel is modeled as an ellipsoid described by a 3×3 tensor. The tensor’s eigenvectors define the principal axes of diffusion, and the eigenvalues (evals) quantify the magnitude of diffusion along each axis. From these, we can derive scalar maps such as:

• Fractional Anisotropy (FA): ranges from 0 (isotropic — equal diffusion in all directions, e.g., CSF) to 1 (highly directional diffusion, e.g., along a coherent fiber bundle).

We will visualize DTI results in several ways:

• FA map — highlights anisotropic regions like white matter.

• Tensor ellipsoids — show the shape and orientation of diffusion tensors in small regions (e.g., corpus callosum).

• Tensor ODFs — visualize the diffusion orientation profile derived from the single tensor in each voxel.

The fitting will not be accurate in the background of the image as there is no signal and possibly we will find FA values with NaNs (not a number). We can easily remove these in the following way.

We can also compute the colored FA or RGB-map. First, FA needs to be scaled between 0 and 1, then the RGB map can be computed and visualized.

Let’s also visualize the tensor ellipsoids of a small rectangular area in an axial slice of the corpus callosum (CC):

Finally, we will visualize the tensor Orientation Distribution Functions for the same area as we did with the ellipsoids.

B) CSA (Constant Solid Angle) ODF¶

Here, we’ll use peaks_from_model to apply the CsaOdfModel to each white matter voxel and estimate fiber orientations that we can use for tracking. This function requires a reconstruction model, the diffusion data, and a sphere as input. The sphere defines a set of discrete directions on the unit sphere at which ODF values are evaluated — DIPY provides several pre-computed spheres (e.g., default_sphere with 362 points, repulsion724 with 724 points) where the points are distributed to provide uniform angular coverage.
The csa_peaks object contains several useful attributes: gfa (generalized fractional anisotropy — a multi-shell generalization of FA), peak_dirs (peak directions), peak_values (maxima values of the ODF), peak_indices (position on the discrete sphere), shm_coeff (spherical harmonic coefficients) and odf (the full orientation distribution function).

• DIPY’s example Reconstruction with Constant Solid Angle

For quality assurance, we can visualize an axial slice from the direction field estimated by the CSA model. Each glyph represents a local fiber orientation derived from the ODF peaks. The underlying ODFs are not shown here.

C) CSD (Constrained Spherical Deconvolution)¶

Constrained Spherical Deconvolution (CSD) reconstructs the fiber Orientation Distribution Function (fODF) by estimating how the diffusion signal would look for a single, coherently oriented fiber (the response function) and then deconvolving this response from the measured signal. This allows recovery of multiple fiber orientations within a voxel, making it particularly effective in regions with crossing fibers.

• DIPY’s CSD reconstruction example

First, let’s verify the b-values of the dataset by looking at the attribute gtab.bvals.

Given the multiple gradient directions in this dataset, we can proceed with the two steps of CSD:

1. Estimation of the fiber response function

2. Use the response function to reconstruct the fODF

Step 1. Estimation of the fiber response function from a local brain region¶

To estimate the single-fiber response function, we can use a brain region where fibers are known to be coherent and aligned, such as the corpus callosum. The auto_response_ssst function automates this process: it selects a small cuboid ROI (typically centered in the brain), computes FA values, and extracts the response function from voxels with FA > 0.7 — likely corresponding to single-fiber regions.

Internally, this function first creates a mask of suitable voxels (mask_for_response_ssst) and then computes the response function within that mask (response_from_mask_ssst). These two steps can also be performed manually; the resulting responses should be identical.

The response tuple contains two elements: an array with the eigenvalues of the response function and the average S0 for this response:

The response tensor should be prolate and anisotropic, with the axial diffusivity about five times greater than the radial diffusivity (i.e., the two smaller eigenvalues are equal and roughly one-fifth of the largest).

We can double-check that we have a good response function by visualizing the response function’s ODF. Here is how you would do that:

Step 2. fODF reconstruction¶

After estimating the response function, we can start the deconvolution process:

We can also find the peak directions (maxima) of the ODFs with peaks_from_model.

We now visualize both the ODFs and peaks in the same space.

Why DTI, CSA and CSD? DTI provides intuitive scalar maps (FA/MD) but only estimates a single fiber direction per voxel. CSA and CSD can both resolve multiple fiber orientations, with CSD producing sharper FODs that improve tractography accuracy in complex white matter regions.

Step 2: Stopping criterion¶

Before running tractography, we must define when to stop tracking. The stopping criterion determines whether tracking continues or terminates at each step along a streamline. Tracking typically stops when:

• The streamline reaches a region with unreliable diffusion information (e.g., low anisotropy)

• It exits the white matter into gray matter or CSF

• It reaches the image boundaries

In this notebook, we’ll explore two basic stopping criteria:

1. Threshold Stopping Criterion — stops when a scalar metric (FA or GFA) falls below a threshold

2. Binary Stopping Criterion — stops when the streamline exits a predefined white matter mask

Both approaches produce streamlines classified as either:

• Valid: ending at appropriate stopping points (e.g., at mask boundaries or image edges)

• Invalid: terminating prematurely in regions where tracking should continue

Note: More advanced stopping criteria, such as Anatomically-Constrained Tractography (ACT), will be covered in a separate notebook on advanced tracking methods.

1️⃣ Threshold Stopping Criterion

The threshold stopping criterion stops tracking when an interpolated scalar value falls below a specified threshold. We’ll demonstrate this using the FA map from DTI with a threshold of 0.15. At each tracking step, the FA value is estimated using trilinear interpolation at the streamline’s current position. If FA < 0.15, tracking stops.

Streamline stopping states:

• ENDPOINT: tracking stopped because the FA value fell below the threshold

• OUTSIDEIMAGE: tracking stopped because the position was outside the image domain

• TRACKPOINT: tracking stopped because no valid direction was available, even though FA ≥ threshold at that location

This criterion is data-driven and flexible, but streamlines may occasionally extend slightly beyond the anatomical white matter boundaries.

2️⃣ Binary Stopping Criterion

Instead of stopping based on a continuous scalar map (like FA), the binary criterion relies solely on a white matter mask. Tracking is allowed only inside the mask—as soon as the current position lies outside it, the streamline terminates.

This criterion uses nearest-neighbor interpolation to determine whether the current tracking location is inside (mask = 1) or outside (mask = 0) the white matter.

Streamline stopping states:

• ENDPOINT: tracking stopped because the mask value was 0 (outside white matter)

• OUTSIDEIMAGE: tracking stopped because the position was outside the image domain

• TRACKPOINT: tracking stopped due to no available direction at that location, even though it was inside the mask

This approach is simpler and ensures that streamlines remain strictly within the white matter mask.

Step 3: Seeding¶

Next, we define where tracking starts — the seed points. Each seed acts as a starting location for streamline propagation based on local fiber orientations.

You can choose seeds:

• In a specific region of interest (ROI), such as the corpus callosum.

• Or throughout the entire white matter to reconstruct whole-brain tractography.

For example, to seed in the corpus callosum:

If you want to do whole-brain tractography instead:

Step 4: Local Tracking: Probabilistic and Deterministic Approaches¶

With the diffusion model (Step 1), stopping criteria (Step 2), and seeds (Step 3) defined, we are ready to perform tractography.

DIPY supports two main approaches for local tracking:

• Deterministic: follows the principal fiber orientation at each step.

• Probabilistic: samples directions from the fiber orientation distribution, accounting for uncertainty.

Both approaches use the same seeds and stopping criteria but differ in how the next tracking direction is chosen.

Probabilistic tractography¶

Probabilistic fiber tracking is one approach for reconstructing white matter connections using diffusion MRI. Unlike deterministic tractography, which follows the peak fiber orientation at each step, probabilistic tractography samples the direction of each step from a distribution of possible orientations. In this approach, the Fiber Orientation Distribution (FOD) estimated by the CSD model encodes the relative strength of different fiber orientations within each voxel. The FOD can be discretized on a well-distributed sphere to obtain directional weights (akin to a probability mass function, PMF) for sampling tracking directions. Since a PMF cannot have negative values, any negative FOD amplitudes are clipped to zero.

This probabilistic sampling allows streamlines to explore multiple plausible pathways, capturing uncertainty in local fiber orientation. This is particularly valuable in complex white matter regions with crossing or branching fibers, where a single deterministic direction may oversimplify the underlying anatomy. As with deterministic tracking, probabilistic streamlines begin from seed points and are guided by stopping criteria such as white matter masks or FA thresholds.

Local Tracking is used for probabilistic tractography, which takes the direction getter along with the stopping criterion and seeds as input.

Since we seeded streamlines exclusively from the corpus callosum, all resulting fibers represent interhemispheric connections originating from this region. Let’s visualize them using probabilistic direction getter from SH (peaks_from_model):

Deterministic tractography¶

We now perform deterministic tractography using the fiber orientation distributions (FODs) estimated with the CSD model. Instead of sampling directions probabilistically, the Deterministic Maximum Direction Getter always follows the direction with the highest value in the distribution at each step, subject to tracking constraints (e.g., maximum turning angle of 30°).

This method uses the full orientation distribution rather than only peak directions, distinguishing it from approaches like EuDX, which rely on peak-following. The resulting streamlines represent the direction of maximum orientation likelihood starting from the defined seed mask.

Let’s visualize the Corpus Callosum using deterministic maximum direction getter:

The following alternative visualizations display the corpus callosum streamlines overlaid on T1-weighted anatomical slices (axial at z=35 and sagittal at x=40), with the seed ROI shown as a semi-transparent yellow contour.

The sagittal view illustrates these interhemispheric pathways connecting the two hemispheres.

Quantitative Comparison of Tracking Methods¶

As a final step, we compare the deterministic and probabilistic tractography results quantitatively. Since both methods used identical seeds, stopping criteria, and tracking parameters, we can directly compare their outputs in terms of streamline count and length distributions. Streamline length is a basic metric that provides insight into the spatial extent of reconstructed pathways. By examining the distribution of lengths across all streamlines, we can assess whether the two tracking approaches produce similar or systematically different patterns for this particular white matter structure.

Dependencies in Jupyter/Python¶

• Using the package watermark to document system environment and software versions used in this notebook, alongside the Neurodesktop version extracted from the JUPYTER_IMAGE or NEURODESKTOP_VERSION environment variables.