Note
Go to the end to download the full example code
Bootstrap and Closest Peak Direction Getters Example#
This example shows how choices in direction-getter impact fiber tracking results by demonstrating the bootstrap direction getter (a type of probabilistic tracking, as described in Berman et al. (2008) [Berman2008] a nd the closest peak direction getter (a type of deterministic tracking). (Amirbekian, PhD thesis, 2016)
This example is an extension of the Introduction to Basic Tracking example. Let’s start by loading the necessary modules for executing this tutorial.
from dipy.core.gradients import gradient_table
from dipy.data import get_fnames, small_sphere
from dipy.direction import BootDirectionGetter, ClosestPeakDirectionGetter
from dipy.io.gradients import read_bvals_bvecs
from dipy.io.image import load_nifti, load_nifti_data
from dipy.io.stateful_tractogram import Space, StatefulTractogram
from dipy.io.streamline import save_trk
from dipy.reconst.csdeconv import (ConstrainedSphericalDeconvModel,
auto_response_ssst)
from dipy.reconst.shm import CsaOdfModel
from dipy.tracking import utils
from dipy.tracking.local_tracking import LocalTracking
from dipy.tracking.streamline import Streamlines
from dipy.tracking.stopping_criterion import ThresholdStoppingCriterion
from dipy.viz import window, actor, colormap, has_fury
# Enables/disables interactive visualization
interactive = False
hardi_fname, hardi_bval_fname, hardi_bvec_fname = get_fnames('stanford_hardi')
label_fname = get_fnames('stanford_labels')
data, affine, hardi_img = load_nifti(hardi_fname, return_img=True)
labels = load_nifti_data(label_fname)
bvals, bvecs = read_bvals_bvecs(hardi_bval_fname, hardi_bvec_fname)
gtab = gradient_table(bvals, bvecs)
seed_mask = (labels == 2)
white_matter = (labels == 1) | (labels == 2)
seeds = utils.seeds_from_mask(seed_mask, affine, density=1)
Next, we fit the CSD model.
response, ratio = auto_response_ssst(gtab, data, roi_radii=10, fa_thr=0.7)
csd_model = ConstrainedSphericalDeconvModel(gtab, response, sh_order_max=6)
csd_fit = csd_model.fit(data, mask=white_matter)
0%| | 0/58788 [00:00<?, ?it/s]
1%|█ | 461/58788 [00:00<00:13, 4291.62it/s]
2%|██▊ | 1251/58788 [00:00<00:09, 5968.83it/s]
3%|████▌ | 2019/58788 [00:00<00:08, 6448.60it/s]
5%|██████▏ | 2758/58788 [00:00<00:08, 6569.33it/s]
6%|████████▏ | 3642/58788 [00:00<00:07, 7114.51it/s]
8%|██████████ | 4422/58788 [00:00<00:07, 7120.94it/s]
9%|████████████ | 5306/58788 [00:00<00:07, 7434.41it/s]
11%|██████████████ | 6230/58788 [00:00<00:06, 7757.72it/s]
12%|███████████████▊ | 7006/58788 [00:00<00:06, 7550.48it/s]
13%|█████████████████▋ | 7795/58788 [00:01<00:06, 7446.56it/s]
15%|███████████████████▌ | 8630/58788 [00:01<00:06, 7503.35it/s]
16%|█████████████████████▏ | 9381/58788 [00:01<00:06, 7310.39it/s]
17%|██████████████████████▋ | 10113/58788 [00:01<00:07, 6555.28it/s]
18%|████████████████████████▎ | 10835/58788 [00:01<00:07, 6568.44it/s]
20%|██████████████████████████▎ | 11723/58788 [00:01<00:06, 7017.36it/s]
21%|████████████████████████████▏ | 12570/58788 [00:01<00:06, 7229.66it/s]
23%|█████████████████████████████▊ | 13304/58788 [00:01<00:06, 7076.27it/s]
24%|███████████████████████████████▌ | 14030/58788 [00:01<00:06, 6946.87it/s]
25%|█████████████████████████████████▍ | 14902/58788 [00:02<00:06, 7251.33it/s]
27%|███████████████████████████████████▍ | 15765/58788 [00:02<00:05, 7442.34it/s]
28%|█████████████████████████████████████▏ | 16571/58788 [00:02<00:05, 7420.85it/s]
29%|██████████████████████████████████████▉ | 17315/58788 [00:02<00:05, 7235.54it/s]
31%|████████████████████████████████████████▊ | 18158/58788 [00:02<00:05, 7376.24it/s]
32%|██████████████████████████████████████████▍ | 18897/58788 [00:02<00:05, 7191.10it/s]
34%|████████████████████████████████████████████▏ | 19700/58788 [00:02<00:05, 7236.31it/s]
35%|██████████████████████████████████████████████▎ | 20609/58788 [00:02<00:05, 7557.79it/s]
36%|████████████████████████████████████████████████▏ | 21438/58788 [00:02<00:04, 7564.27it/s]
38%|██████████████████████████████████████████████████▏ | 22335/58788 [00:03<00:04, 7755.24it/s]
39%|███████████████████████████████████████████████████▉ | 23111/58788 [00:03<00:04, 7247.16it/s]
41%|█████████████████████████████████████████████████████▌ | 23841/58788 [00:03<00:04, 7074.54it/s]
42%|███████████████████████████████████████████████████████▍ | 24689/58788 [00:03<00:04, 7278.40it/s]
43%|█████████████████████████████████████████████████████████▏ | 25466/58788 [00:03<00:04, 7228.15it/s]
45%|██████████████████████████████████████████████████████████▉ | 26275/58788 [00:03<00:04, 7278.76it/s]
46%|████████████████████████████████████████████████████████████▋ | 27042/58788 [00:03<00:04, 7199.37it/s]
47%|██████████████████████████████████████████████████████████████▍ | 27808/58788 [00:03<00:04, 7140.76it/s]
49%|████████████████████████████████████████████████████████████████ | 28526/58788 [00:03<00:04, 6968.04it/s]
50%|█████████████████████████████████████████████████████████████████▉ | 29384/58788 [00:04<00:04, 7230.88it/s]
51%|███████████████████████████████████████████████████████████████████▋ | 30128/58788 [00:04<00:04, 7102.52it/s]
52%|█████████████████████████████████████████████████████████████████████▏ | 30839/58788 [00:04<00:04, 6915.31it/s]
54%|██████████████████████████████████████████████████████████████████████▊ | 31531/58788 [00:04<00:04, 6745.36it/s]
55%|████████████████████████████████████████████████████████████████████████▌ | 32338/58788 [00:04<00:03, 6934.90it/s]
57%|██████████████████████████████████████████████████████████████████████████▋ | 33236/58788 [00:04<00:03, 7316.79it/s]
58%|████████████████████████████████████████████████████████████████████████████▋ | 34144/58788 [00:04<00:03, 7612.29it/s]
60%|██████████████████████████████████████████████████████████████████████████████▋ | 35041/58788 [00:04<00:03, 7788.92it/s]
61%|████████████████████████████████████████████████████████████████████████████████▋ | 35960/58788 [00:04<00:02, 7973.36it/s]
63%|██████████████████████████████████████████████████████████████████████████████████▊ | 36893/58788 [00:05<00:02, 8140.33it/s]
64%|█████████████████████████████████████████████████████████████████████████████████████ | 37864/58788 [00:05<00:02, 8361.53it/s]
66%|███████████████████████████████████████████████████████████████████████████████████████ | 38779/58788 [00:05<00:02, 8356.54it/s]
68%|█████████████████████████████████████████████████████████████████████████████████████████ | 39683/58788 [00:05<00:02, 8331.39it/s]
69%|███████████████████████████████████████████████████████████████████████████████████████████▏ | 40625/58788 [00:05<00:02, 8419.51it/s]
71%|█████████████████████████████████████████████████████████████████████████████████████████████ | 41473/58788 [00:05<00:02, 8220.26it/s]
72%|███████████████████████████████████████████████████████████████████████████████████████████████▏ | 42409/58788 [00:05<00:01, 8321.52it/s]
74%|█████████████████████████████████████████████████████████████████████████████████████████████████▎ | 43361/58788 [00:05<00:01, 8435.75it/s]
75%|███████████████████████████████████████████████████████████████████████████████████████████████████▍ | 44280/58788 [00:05<00:01, 8426.12it/s]
77%|█████████████████████████████████████████████████████████████████████████████████████████████████████▍ | 45190/58788 [00:06<00:01, 8394.26it/s]
78%|███████████████████████████████████████████████████████████████████████████████████████████████████████▎ | 46030/58788 [00:06<00:01, 8179.24it/s]
80%|█████████████████████████████████████████████████████████████████████████████████████████████████████████▎ | 46884/58788 [00:06<00:01, 8061.05it/s]
81%|███████████████████████████████████████████████████████████████████████████████████████████████████████████▏ | 47717/58788 [00:06<00:01, 7933.34it/s]
83%|████████████████████████████████████████████████████████████████████████████████████████████████████████████▉ | 48540/58788 [00:06<00:01, 7811.43it/s]
84%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████▏ | 49516/58788 [00:06<00:01, 8145.02it/s]
86%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████ | 50331/58788 [00:06<00:01, 7937.50it/s]
87%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████▊ | 51159/58788 [00:06<00:00, 7826.95it/s]
89%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████ | 52123/58788 [00:06<00:00, 8123.40it/s]
90%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████ | 53004/58788 [00:07<00:00, 8102.06it/s]
92%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▊ | 53828/58788 [00:07<00:00, 7932.16it/s]
93%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▋ | 54622/58788 [00:07<00:00, 7865.30it/s]
94%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▌ | 55487/58788 [00:07<00:00, 7744.00it/s]
96%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▍ | 56311/58788 [00:07<00:00, 7680.53it/s]
97%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▍ | 57190/58788 [00:07<00:00, 7787.62it/s]
99%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▏ | 57969/58788 [00:07<00:00, 7587.76it/s]
100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████▉| 58768/58788 [00:07<00:00, 7502.99it/s]
100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 58788/58788 [00:07<00:00, 7509.89it/s]
we use the CSA fit to calculate GFA, which will serve as our stopping criterion.
csa_model = CsaOdfModel(gtab, sh_order_max=6)
gfa = csa_model.fit(data, mask=white_matter).gfa
stopping_criterion = ThresholdStoppingCriterion(gfa, .25)
Next, we need to set up our two direction getters
Example #1: Bootstrap direction getter with CSD Model#
boot_dg_csd = BootDirectionGetter.from_data(data, csd_model, max_angle=30.,
sphere=small_sphere)
boot_streamline_generator = LocalTracking(boot_dg_csd, stopping_criterion,
seeds, affine, step_size=.5)
streamlines = Streamlines(boot_streamline_generator)
sft = StatefulTractogram(streamlines, hardi_img, Space.RASMM)
save_trk(sft, "tractogram_bootstrap_dg.trk")
if has_fury:
scene = window.Scene()
scene.add(actor.line(streamlines, colormap.line_colors(streamlines)))
window.record(scene, out_path='tractogram_bootstrap_dg.png',
size=(800, 800))
if interactive:
window.show(scene)
Corpus Callosum Bootstrap Probabilistic Direction Getter
We have created a bootstrapped probabilistic set of streamlines. If you repeat the fiber tracking (keeping all inputs the same) you will NOT get exactly the same set of streamlines.
Example #2: Closest peak direction getter with CSD Model#
pmf = csd_fit.odf(small_sphere).clip(min=0)
peak_dg = ClosestPeakDirectionGetter.from_pmf(pmf, max_angle=30.,
sphere=small_sphere)
peak_streamline_generator = LocalTracking(peak_dg, stopping_criterion, seeds,
affine, step_size=.5)
streamlines = Streamlines(peak_streamline_generator)
sft = StatefulTractogram(streamlines, hardi_img, Space.RASMM)
save_trk(sft, "closest_peak_dg_CSD.trk")
if has_fury:
scene = window.Scene()
scene.add(actor.line(streamlines, colormap.line_colors(streamlines)))
window.record(scene, out_path='tractogram_closest_peak_dg.png',
size=(800, 800))
if interactive:
window.show(scene)
Corpus Callosum Closest Peak Deterministic Direction Getter
We have created a set of streamlines using the closest peak direction getter, which is a type of deterministic tracking. If you repeat the fiber tracking (keeping all inputs the same) you will get exactly the same set of streamlines.
References#
[Berman2008]
Berman, J. et al., Probabilistic streamline q-ball
tractography using the residual bootstrap, NeuroImage, vol 39, no 1, 2008
Total running time of the script: (0 minutes 43.602 seconds)