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
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
|
"""
====================
Tracking Quick Start
====================
This example shows how to perform fiber tracking using Dipy.
We will use Constrained Spherical Deconvolution (CSD) [Tournier07]_ for local
reconstructions and then generate deterministic streamlines using the fiber
directions (peaks) from CSD and fractional anisotropic (FA) as a
stopping criterion.
First, let's load the necessary modules.
"""
from dipy.reconst.dti import TensorModel, fractional_anisotropy
from dipy.reconst.csdeconv import (ConstrainedSphericalDeconvModel,
auto_response)
from dipy.reconst.peaks import peaks_from_model
from dipy.tracking.eudx import EuDX
from dipy.data import fetch_stanford_hardi, read_stanford_hardi, get_sphere
from dipy.segment.mask import median_otsu
from dipy.viz import fvtk
from dipy.viz.colormap import line_colors
"""
Load one of the available datasets with 150 gradients on the sphere and 10 b0s
"""
fetch_stanford_hardi()
img, gtab = read_stanford_hardi()
data = img.get_data()
"""
Create a brain mask. This dataset is a bit difficult to segment with the default
``median_otsu`` parameters (see :ref:`example_brain_extraction_dwi`) therefore we use
here a bit more advanced options.
"""
maskdata, mask = median_otsu(data, 3, 1, False,
vol_idx=range(10, 50), dilate=2)
"""
For the constrained spherical deconvolution we need to estimate the response
function (see :ref:`example_reconst_csd`) and create a model.
"""
response, ratio = auto_response(gtab, data, roi_radius=10, fa_thr=0.7)
csd_model = ConstrainedSphericalDeconvModel(gtab, response)
"""
Next, we use ``peaks_from_model`` to fit the data and calculated the fiber
directions in all voxels.
"""
sphere = get_sphere('symmetric724')
csd_peaks = peaks_from_model(model=csd_model,
data=data,
sphere=sphere,
mask=mask,
relative_peak_threshold=.5,
min_separation_angle=25,
parallel=True)
"""
For the tracking part, we will use the fiber directions from the ``csd_model``
but stop tracking in areas where fractional anisotropy (FA) is low (< 0.1).
To derive the FA, used here as a stopping criterion, we would need to fit a
tensor model first. Here, we fit the Tensor using weighted least squares (WLS).
"""
tensor_model = TensorModel(gtab, fit_method='WLS')
tensor_fit = tensor_model.fit(data, mask)
FA = fractional_anisotropy(tensor_fit.evals)
"""
In order for the stopping values to be used with our tracking algorithm we need
to have the same dimensions as the ``csd_peaks.peak_values``. For this reason,
we can assign the same FA value to every peak direction in the same voxel in
the following way.
"""
stopping_values = np.zeros(csd_peaks.peak_values.shape)
stopping_values[:] = FA[..., None]
"""
For quality assurance we can also visualize a slice from the direction field
which we will use as the basis to perform the tracking.
"""
ren = fvtk.ren()
slice_no = data.shape[2] / 2
fvtk.add(ren, fvtk.peaks(csd_peaks.peak_dirs[:, :, slice_no:slice_no + 1],
stopping_values[:, :, slice_no:slice_no + 1]))
print('Saving illustration as csd_direction_field.png')
fvtk.record(ren, out_path='csd_direction_field.png', size=(900, 900))
"""
.. figure:: csd_direction_field.png
:align: center
**Direction Field (peaks)**
``EuDX`` [Garyfallidis12]_ is a fast algorithm that we use here to generate
streamlines. If the parameter ``seeds`` is a positive integer it will generate
that number of randomly placed seeds everywhere in the volume. Alternatively,
you can specify the exact seed points using an array (N, 3) where N is the
number of seed points. For simplicity, here we will use the first option
(random seeds). ``a_low`` is the threshold of the fist parameter
(``stopping_values``) which means that there will that tracking will stop in
regions with FA < 0.1.
"""
streamline_generator = EuDX(stopping_values,
csd_peaks.peak_indices,
seeds=10**4,
odf_vertices=sphere.vertices,
a_low=0.1)
streamlines = [streamline for streamline in streamline_generator]
"""
We can visualize the streamlines using ``fvtk.line`` or ``fvtk.streamtube``.
"""
fvtk.clear(ren)
fvtk.add(ren, fvtk.line(streamlines, line_colors(streamlines)))
print('Saving illustration as csd_streamlines_eudx.png')
fvtk.record(ren, out_path='csd_streamlines_eudx.png', size=(900, 900))
"""
.. figure:: csd_streamlines_eudx.png
:align: center
**CSD-based streamlines using EuDX**
We used above ``fvtk.record`` because we want to create a figure for the tutorial
but you can visualize the same objects in 3D using ``fvtk.show(ren)``.
To learn more about this process you could start playing with the number of
seed points or even better specify seeds to be in specific regions of interest
in the brain.
``fvtk`` gives some minimal interactivity however you can save the resulting
streamlines in a Trackvis (.trk) format and load them for example with the
Fibernavigator_ or another tool for medical visualization.
Finally, let's save the streamlines as a (.trk) file and FA as a Nifti image.
"""
import nibabel as nib
hdr = nib.trackvis.empty_header()
hdr['voxel_size'] = img.get_header().get_zooms()[:3]
hdr['voxel_order'] = 'LAS'
hdr['dim'] = FA.shape[:3]
csd_streamlines_trk = ((sl, None, None) for sl in streamlines)
csd_sl_fname = 'csd_streamline.trk'
nib.trackvis.write(csd_sl_fname, csd_streamlines_trk, hdr, points_space='voxel')
nib.save(nib.Nifti1Image(FA, img.get_affine()), 'FA_map.nii.gz')
"""
.. [Garyfallidis12] Garyfallidis E., "Towards an accurate brain tractography", PhD thesis, University of Cambridge, 2012.
.. [Tournier07] J-D. Tournier, F. Calamante and A. Connelly, "Robust determination of the fibre orientation distribution in diffusion MRI: Non-negativity constrained super-resolved spherical deconvolution", Neuroimage, vol. 35, no. 4, pp. 1459-1472, 2007.
.. include:: ../links_names.inc
"""
|
