-
Notifications
You must be signed in to change notification settings - Fork 429
/
tracking_bootstrap_peaks.py
139 lines (109 loc) · 4.61 KB
/
tracking_bootstrap_peaks.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
"""
====================================================
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 :ref:`example_tracking_introduction_eudx`
example. Let's start by loading the necessary modules for executing this
tutorial.
"""
# Enables/disables interactive visualization
interactive = False
from dipy.data import read_stanford_labels, small_sphere
from dipy.direction import BootDirectionGetter, ClosestPeakDirectionGetter
from dipy.io.stateful_tractogram import Space, StatefulTractogram
from dipy.io.streamline import save_trk
from dipy.reconst.csdeconv import (ConstrainedSphericalDeconvModel,
auto_response)
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
hardi_img, gtab, labels_img = read_stanford_labels()
data = hardi_img.get_data()
labels = labels_img.get_data()
affine = hardi_img.affine
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(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)
"""
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:
r = window.Renderer()
r.add(actor.line(streamlines, colormap.line_colors(streamlines)))
window.record(r, out_path='tractogram_bootstrap_dg.png', size=(800, 800))
if interactive:
window.show(r)
"""
.. figure:: tractogram_bootstrap_dg.png
:align: center
**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:
r = window.Renderer()
r.add(actor.line(streamlines, colormap.line_colors(streamlines)))
window.record(r, out_path='tractogram_closest_peak_dg.png',
size=(800, 800))
if interactive:
window.show(r)
"""
.. figure:: tractogram_closest_peak_dg.png
:align: center
**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
.. include:: ../links_names.inc
"""