Documentation 1.7.0 >
Examples >
Fiber Tracking >
Bootstrap and Closest Peak Direction Getters Example
Note
Click here 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 example_tracking_introduction_eudx 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=6)
csd_fit = csd_model.fit(data, mask=white_matter)
0%| | 0/58788 [00:00<?, ?it/s]
1%|1 | 665/58788 [00:00<00:08, 6642.17it/s]
3%|2 | 1531/58788 [00:00<00:07, 7828.00it/s]
4%|4 | 2411/58788 [00:00<00:06, 8271.13it/s]
6%|5 | 3299/58788 [00:00<00:06, 8509.40it/s]
7%|7 | 4191/58788 [00:00<00:06, 8654.36it/s]
9%|8 | 5101/58788 [00:00<00:06, 8805.25it/s]
10%|# | 6000/58788 [00:00<00:05, 8863.37it/s]
12%|#1 | 6915/58788 [00:00<00:05, 8952.84it/s]
13%|#3 | 7820/58788 [00:00<00:05, 8982.87it/s]
15%|#4 | 8719/58788 [00:01<00:05, 8970.62it/s]
16%|#6 | 9642/58788 [00:01<00:05, 9048.76it/s]
18%|#7 | 10558/58788 [00:01<00:05, 9081.19it/s]
20%|#9 | 11467/58788 [00:01<00:05, 9072.76it/s]
21%|##1 | 12385/58788 [00:01<00:05, 9101.37it/s]
23%|##2 | 13320/58788 [00:01<00:04, 9174.78it/s]
24%|##4 | 14238/58788 [00:01<00:04, 9108.06it/s]
26%|##5 | 15159/58788 [00:01<00:04, 9135.84it/s]
27%|##7 | 16096/58788 [00:01<00:04, 9204.41it/s]
29%|##8 | 17017/58788 [00:01<00:04, 9188.48it/s]
31%|### | 17959/58788 [00:02<00:04, 9257.33it/s]
32%|###2 | 18885/58788 [00:02<00:04, 9180.69it/s]
34%|###3 | 19804/58788 [00:02<00:04, 9133.31it/s]
35%|###5 | 20727/58788 [00:02<00:04, 9161.76it/s]
37%|###6 | 21644/58788 [00:02<00:04, 9127.46it/s]
38%|###8 | 22581/58788 [00:02<00:03, 9187.66it/s]
40%|###9 | 23500/58788 [00:02<00:03, 9166.85it/s]
42%|####1 | 24417/58788 [00:02<00:03, 9164.59it/s]
43%|####3 | 25349/58788 [00:02<00:03, 9210.42it/s]
45%|####4 | 26271/58788 [00:02<00:03, 9127.94it/s]
46%|####6 | 27184/58788 [00:03<00:03, 9050.06it/s]
48%|####7 | 28090/58788 [00:03<00:03, 8994.39it/s]
49%|####9 | 28990/58788 [00:03<00:03, 8908.74it/s]
51%|##### | 29882/58788 [00:03<00:03, 8844.98it/s]
52%|#####2 | 30767/58788 [00:03<00:03, 8790.85it/s]
54%|#####3 | 31647/58788 [00:03<00:03, 8746.60it/s]
55%|#####5 | 32540/58788 [00:03<00:02, 8799.29it/s]
57%|#####6 | 33449/58788 [00:03<00:02, 8882.81it/s]
58%|#####8 | 34362/58788 [00:03<00:02, 8955.76it/s]
60%|#####9 | 35258/58788 [00:03<00:02, 8891.18it/s]
62%|######1 | 36181/58788 [00:04<00:02, 8989.83it/s]
63%|######3 | 37093/58788 [00:04<00:02, 9028.20it/s]
65%|######4 | 38023/58788 [00:04<00:02, 9108.31it/s]
66%|######6 | 38954/58788 [00:04<00:02, 9168.38it/s]
68%|######7 | 39877/58788 [00:04<00:02, 9185.36it/s]
69%|######9 | 40807/58788 [00:04<00:01, 9217.98it/s]
71%|####### | 41737/58788 [00:04<00:01, 9241.37it/s]
73%|#######2 | 42662/58788 [00:04<00:01, 9213.28it/s]
74%|#######4 | 43600/58788 [00:04<00:01, 9261.21it/s]
76%|#######5 | 44527/58788 [00:04<00:01, 9172.68it/s]
77%|#######7 | 45452/58788 [00:05<00:01, 9193.13it/s]
79%|#######8 | 46372/58788 [00:05<00:01, 9174.07it/s]
80%|######## | 47290/58788 [00:05<00:01, 9141.84it/s]
82%|########2 | 48209/58788 [00:05<00:01, 9155.59it/s]
84%|########3 | 49148/58788 [00:05<00:01, 9225.22it/s]
85%|########5 | 50071/58788 [00:05<00:00, 9182.48it/s]
87%|########6 | 50990/58788 [00:05<00:00, 9165.69it/s]
88%|########8 | 51907/58788 [00:05<00:00, 9155.53it/s]
90%|########9 | 52823/58788 [00:05<00:00, 9134.20it/s]
91%|#########1| 53737/58788 [00:05<00:00, 9095.91it/s]
93%|#########2| 54647/58788 [00:06<00:00, 9095.53it/s]
95%|#########4| 55562/58788 [00:06<00:00, 9109.65it/s]
96%|#########6| 56473/58788 [00:06<00:00, 9003.52it/s]
98%|#########7| 57374/58788 [00:06<00:00, 8955.18it/s]
99%|#########9| 58270/58788 [00:06<00:00, 8782.56it/s]
100%|##########| 58788/58788 [00:06<00:00, 8991.24it/s]
we use the CSA fit to calculate GFA, which will serve as our stopping criterion.
csa_model = CsaOdfModel(gtab, sh_order=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)
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)
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 44.033 seconds)