Tractography#

Tractography follows the fibre directions of the FODs from seed points, step by step, into streamlines. Tracker has two algorithms: "det" moves along the FOD peak nearest the current direction, and "prob" draws each direction at random according to the FOD. By default each streamline is tracked both ways from its seed. Citing lists the papers, and Conventions the point order and the file formats.

The data on this page

The examples on this page run on a small synthetic FOD image, written to the current folder by the code below: fod.nii.gz, SH coefficients up to degree 8 in the tournier basis, and mask.nii.gz, the voxels of its two bundles. It has 20 x 20 x 4 voxels of 2 mm. One bundle runs along x and one along y, and where they cross the FOD has both lobes. To use your own data, use your files instead.

import numpy as np

from larmor.io import save_volume
from larmor.sd import fit_sh
from larmor.types import Volume

n = 600
z = 1 - (2 * np.arange(n) + 1) / n
azimuth = np.pi * (3 - np.sqrt(5)) * np.arange(n)
radius = np.sqrt(1 - z**2)
sphere = np.stack([radius * np.cos(azimuth), radius * np.sin(azimuth), z], axis=1)


def fod(*fibres):
    """SH coefficients of a narrow lobe along each fibre."""
    cosines = sphere @ np.array(fibres).T
    return fit_sh(np.exp(-8 * (1 - cosines**2)).sum(axis=1), sphere, lmax=8)


x, y = np.eye(3)[0], np.eye(3)[1]
along_x = np.zeros((20, 20, 4), bool)
along_x[:, 8:12] = True
along_y = np.zeros((20, 20, 4), bool)
along_y[8:12, :] = True
coeffs = np.zeros((20, 20, 4, 45), np.float32)
coeffs[along_x] = fod(x)
coeffs[along_y] = fod(y)
coeffs[along_x & along_y] = fod(x, y)

affine = np.diag([2.0, 2.0, 2.0, 1.0])
save_volume(Volume(coeffs, affine), "fod.nii.gz")
save_volume(Volume((along_x | along_y).astype(np.uint8), affine), "mask.nii.gz")

Tracking#

from larmor.io import load_volume
from larmor.tractography import MaskSeeds, Tracker

fod = load_volume("fod.nii.gz")
mask = load_volume("mask.nii.gz")
tracker = Tracker("det")
run = tracker.track(
    fod, "tracks.tck", MaskSeeds(mask), num_streamlines=200, mask=mask, seed=0
)
print(f"{run.num_streamlines} streamlines from {run.num_seeds} seeds")

track() writes the tractogram batch by batch and returns a TractographyRun: the parameters used, the number of seeds and streamlines, the seeds rejected in num_rejected (no_start where no direction above the cutoff was found, too_short), and, with None as the output, the streamlines themselves in run.tractogram (for small runs only).

lar track fod.nii.gz tracks.tck --seed-image mask.nii.gz --select 200 --mask mask.nii.gz --seed 0

The progress bar counts the seeds tracked and those rejected; --json run.json saves the record of the run: every parameter used, the counts and the versions.

Seeding#

Seeds

Python

Command line

Random positions in a mask’s voxels, until a number of streamlines are kept

MaskSeeds(mask) with num_streamlines

--seed-image mask.nii.gz --select 200

N random positions in every voxel of a mask

seeds_per_voxel(mask, n, rng=0)

--seed-random-per-voxel mask.nii.gz 2

Seeds saved from an earlier run

read_seeds("seeds.txt")

--seed-file seeds.txt

Fixed seeds (Seeds) are all tracked unless num_streamlines is given. The seeds of the kept streamlines are saved, in output order, with output_seeds= or --output-seeds; tracked again from that file with the same run seed, they give the same streamlines.

from larmor.tractography import seeds_per_voxel

seeds = seeds_per_voxel(mask, 2, rng=0)
run = Tracker("prob").track(fod, "prob.tck", seeds, mask=mask, seed=0)

lar track fod.nii.gz prob.tck --algorithm prob \
    --seed-random-per-voxel mask.nii.gz 2 --mask mask.nii.gz --seed 0

Stopping#

A streamline stops where it:

  • leaves the tracking mask (mask=, --mask), or the image;

  • reaches an FOD amplitude below the cutoff (cutoff, --cutoff): by default "auto", 0.2 times the mean largest peak of the FOD;

  • turns by more than the maximum angle between steps (max_angle, --angle): 60 degrees for det and 15 for prob by default;

  • reaches the maximum length (max_length, --max-length), 100 times the voxel size by default.

Streamlines shorter than min_length (--min-length, 5 times the voxel size by default) are rejected. The step size (step_size, --step-size) is 0.1 times the voxel size by default, and bidirectional=False (--unidirectional) tracks each seed one way only. The reason each end stopped is kept per streamline as end_status (EndStatus), which TRK files keep and TCK files cannot hold.

Reproducible runs#

The run seed (seed=, --seed, in [0, 2^64)) fixes the seeds drawn and every random number, so the same seed gives the same tractogram whatever the number of threads. Without one, a seed is drawn and recorded in the run record.

Output formats#

track writes TCK or TRK files. To work with other programs, convert the result, for example to VTK:

lar convert tracks.tck tracks.vtk

Loading and saving data covers the formats, and Isolating a bundle and Tract profiles the next steps.

Viewing the streamlines#

lar bv -t tracks.tck -f fod.nii.gz

See also#