-
Notifications
You must be signed in to change notification settings - Fork 429
/
bundle_extraction.py
215 lines (154 loc) · 5.63 KB
/
bundle_extraction.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
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
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
"""
==================================================
Automatic Fiber Bundle Extraction with RecoBundles
==================================================
This example explains how we can use RecoBundles [Garyfallidis17]_ to
extract bundles from tractograms.
First import the necessary modules.
"""
import numpy as np
from dipy.segment.bundles import RecoBundles
from dipy.align.streamlinear import whole_brain_slr
from dipy.viz import window, actor
from dipy.io.streamline import load_trk, save_trk
"""
Download and read data for this tutorial
"""
from dipy.data.fetcher import (fetch_target_tractogram_hcp,
fetch_bundle_atlas_hcp842,
get_bundle_atlas_hcp842,
get_target_tractogram_hcp)
target_file, target_folder = fetch_target_tractogram_hcp()
atlas_file, atlas_folder = fetch_bundle_atlas_hcp842()
atlas_file, all_bundles_files = get_bundle_atlas_hcp842()
target_file = get_target_tractogram_hcp()
atlas, atlas_header = load_trk(atlas_file)
target, target_header = load_trk(target_file)
"""
let's visualize atlas tractogram and target tractogram before registration
"""
interactive = False
ren = window.Renderer()
ren.SetBackground(1, 1, 1)
ren.add(actor.line(atlas, colors=(1,0,1)))
ren.add(actor.line(target, colors=(1,1,0)))
window.record(ren, out_path='tractograms_initial.png', size=(600, 600))
if interactive:
window.show(ren)
"""
.. figure:: tractograms_initial.png
:align: center
Atlas and target before registration.
"""
"""
We will register target tractogram to model atlas' space using streamlinear
registeration (SLR) [Garyfallidis15]_
"""
moved, transform, qb_centroids1, qb_centroids2 = whole_brain_slr(
atlas, target, x0='affine', verbose=True, progressive=True)
"""
We save the transform generated in this registration, so that we can use
it in the bundle profiles example
"""
np.save("slr_transform.npy", transform)
"""
let's visualize atlas tractogram and target tractogram after registration
"""
interactive = False
ren = window.Renderer()
ren.SetBackground(1, 1, 1)
ren.add(actor.line(atlas, colors=(1,0,1)))
ren.add(actor.line(moved, colors=(1,1,0)))
window.record(ren, out_path='tractograms_after_registration.png',
size=(600, 600))
if interactive:
window.show(ren)
"""
.. figure:: tractograms_after_registration.png
:align: center
Atlas and target after registration.
"""
"""
Read AF left and CST left bundles from already fetched atlas data to use them
as model bundles
"""
from dipy.data.fetcher import get_two_hcp842_bundles
model_af_l_file, model_cst_l_file = get_two_hcp842_bundles()
"""
Extracting bundles using recobundles [Garyfallidis17]_
"""
model_af_l, hdr = load_trk(model_af_l_file)
rb = RecoBundles(moved, verbose=True)
recognized_af_l, af_l_labels = rb.recognize(model_bundle=model_af_l,
model_clust_thr=5.,
reduction_thr=10,
reduction_distance='mam',
slr=True,
slr_metric='asymmetric',
pruning_distance='mam')
"""
let's visualize extracted Arcuate Fasciculus Left bundle and model bundle
together
"""
interactive = False
ren = window.Renderer()
ren.SetBackground(1, 1, 1)
ren.add(actor.line(model_af_l, colors=(.1,.7,.26)))
ren.add(actor.line(recognized_af_l, colors=(.1,.1,6)))
ren.set_camera(focal_point=(320.21296692, 21.28884506, 17.2174015),
position=(2.11, 200.46, 250.44) , view_up=(0.1, -1.028, 0.18))
window.record(ren, out_path='AF_L_recognized_bundle.png',
size=(600, 600))
if interactive:
window.show(ren)
"""
.. figure:: AF_L_recognized_bundle.png
:align: center
Extracted Arcuate Fasciculus Left bundle and model bundle
"""
"""
Save the bundle as a trk file. Rather than saving the recognized streamlines
in the space of the atlas, we save the streamlines that are in the original
space of the subject anatomy.
"""
save_trk( "AF_L.trk", target[af_l_labels], hdr['voxel_to_rasmm'])
model_cst_l, hdr = load_trk(model_cst_l_file)
recognized_cst_l, cst_l_labels = rb.recognize(model_bundle=model_cst_l,
model_clust_thr=5.,
reduction_thr=10,
reduction_distance='mam',
slr=True,
slr_metric='asymmetric',
pruning_distance='mam')
"""
let's visualize extracted Corticospinal Tract (CST) Left bundle and model
bundle together
"""
interactive = False
ren = window.Renderer()
ren.SetBackground(1, 1, 1)
ren.add(actor.line(model_cst_l, colors=(.1,.7,.26)))
ren.add(actor.line(recognized_cst_l, colors=(.1,.1,6)))
ren.set_camera(focal_point=(-18.17281532, -19.55606842, 6.92485857),
position=(-360.11, -340.46, -40.44),
view_up=(-0.03, 0.028, 0.89))
window.record(ren, out_path='CST_L_recognized_bundle.png',
size=(600, 600))
if interactive:
window.show(ren)
"""
.. figure:: CST_L_recognized_bundle.png
:align: center
Extracted Corticospinal Tract (CST) Left bundle and model bundle
"""
"""
Save the bundle as a trk file:
"""
save_trk("CST_L.trk", target[cst_l_labels], hdr['voxel_to_rasmm'])
"""
References
----------
.. [Garyfallidis17] Garyfallidis et al. Recognition of white matter
bundles using local and global streamline-based registration
and clustering, Neuroimage, 2017.
"""