/
plot_003_viz.py
474 lines (392 loc) · 19.1 KB
/
plot_003_viz.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
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
"""
============================
Visualizing AFQ derivatives
============================
Visualizing the results of a pyAFQ analysis is useful because it allows us to
inspect the results of the analysis and to communicate the results to others.
The pyAFQ pipeline produces a number of different kinds of outputs, including
visualizations that can be used for quality control and for quick examination
of the results of the analysis.
However, when communicating the results of pyAFQ analysis, it is often useful
to have more specific control over the visualization that is produced. In
addition, it is often useful to have visualizations that are visually appealing
and striking. In this tutorial, we will use the `fury <https://fury.gl/>`_
library [1]_ to visualize outputs of pyAFQ as publication-ready figures.
"""
import os
import os.path as op
import nibabel as nib
import numpy as np
from dipy.io.streamline import load_trk
from dipy.tracking.streamline import transform_streamlines
from fury import actor, window
from fury.colormap import create_colormap
import AFQ.data.fetch as afd
from AFQ.viz.utils import PanelFigure
##############################################################################
# Get some data from HBN POD2
# ----------------------------
# The Healthy Brain Network Preprocessed Open Diffusion Derivatives (HBN POD2)
# is a collection of resources based on the Healthy Brain Network dataset
# [2, 3]_. HBN POD2 includes data derivatives - including pyAFQ derivatives -
# from more than 2,000 subjects. The data and the derivatives can be browsed at
# https://fcp-indi.s3.amazonaws.com/index.html#data/Projects/HBN/BIDS_curated/
#
# Here, we will visualize the results from one subject, together with their
# anatomy and using several variations. We start by downloading their
# pyAFQ-processed data using fetcher functions that download both the
# preprocessed data, as well as the pyAFQ-processed data (Note that this
# will take up about 1.75 GB of disk space):
afd.fetch_hbn_preproc(["NDARAA948VFH"])
study_path = afd.fetch_hbn_afq(["NDARAA948VFH"])[1]
deriv_path = op.join(
study_path, "derivatives")
afq_path = op.join(
deriv_path,
'afq',
'sub-NDARAA948VFH',
'ses-HBNsiteRU')
bundle_path = op.join(afq_path,
'bundles')
#############################################################################
# Read data into memory
# ----------------------
# The bundle coordinates from pyAFQ are always saved in the reference frame of
# the diffusion data from which they are generated, so we need an image file
# with the dMRI coordinates as a reference for loading the data (we could also
# use `"same"` here).
fa_img = nib.load(op.join(afq_path,
'sub-NDARAA948VFH_ses-HBNsiteRU_acq-64dir_space-T1w_desc-preproc_dwi_model-DKI_FA.nii.gz'))
fa = fa_img.get_fdata()
sft_arc = load_trk(op.join(bundle_path,
'sub-NDARAA948VFH_ses-HBNsiteRU_acq-64dir_space-T1w_desc-preproc_dwi_space-RASMM_model-CSD_desc-prob-afq-ARC_L_tractography.trk'), fa_img)
sft_cst = load_trk(op.join(bundle_path,
'sub-NDARAA948VFH_ses-HBNsiteRU_acq-64dir_space-T1w_desc-preproc_dwi_space-RASMM_model-CSD_desc-prob-afq-CST_L_tractography.trk'), fa_img)
#############################################################################
# Transform into the T1w reference frame
# --------------------------------------
# Our first goal is to visualize the bundles with a background of the
# T1-weighted image, which provides anatomical context. We read in this data and
# transform the bundle coordinates, first into the RASMM common coordinate frame
# and then subsequently into the coordinate frame of the T1-weighted data (if
# you find this confusing, you can brush up on this topic in the
# `nibabel documentation <https://nipy.org/nibabel/coordinate_systems.html>`_).
t1w_img = nib.load(op.join(deriv_path,
'qsiprep/sub-NDARAA948VFH/anat/sub-NDARAA948VFH_desc-preproc_T1w.nii.gz'))
t1w = t1w_img.get_fdata()
sft_arc.to_rasmm()
sft_cst.to_rasmm()
arc_t1w = transform_streamlines(sft_arc.streamlines,
np.linalg.inv(t1w_img.affine))
cst_t1w = transform_streamlines(sft_cst.streamlines,
np.linalg.inv(t1w_img.affine))
#############################################################################
#
# .. note::
# A virtual frame buffer is needed if you are running this example on
# a machine that is not connected to a display ("headless"). If this is
# the case, you can either set an environment variable called `XVFB` to `1`
# or you can deindent the following code (and comment out the `if` statement)
# to initialize the virtual frame buffer.
if os.environ.get("XVFB", False):
print("Initializing XVFB")
import xvfbwrapper
from xvfbwrapper import Xvfb
vdisplay = Xvfb()
vdisplay.start()
#############################################################################
# Visualizing bundles with principal direction coloring
# -----------------------------------------------------
# The first visualization we will create will have the streamlines colored
# according to their direction. The color of each streamline will be RGB encoded
# according to the RAS of the average orientation of its segments. This is the
# default behavior in fury.
#
# Fury uses "actors" to render different kind of graphics. For the bundles, we
# will use the `line` actor. These objects are wrappers around the `vtkActor
# class <https://vtk.org/doc/nightly/html/classvtkActor.html>`_, so methods of
# that class (like `GetProperty()`) are available to use. We like to set the
# aesthetics of the streamlines, so that they are rendered as tubes and with
# slightly thicker line-width than the default. We create a function that sets
# these properties of the line actor via the `GetProperty` method. We will
# reuse this function later on, also setting the key-word arguments to the call
# to `actor.line`, but for now we use the default setting, which colors each
# streamline based on the RAS orientation, and we set the line width to 8.
def lines_as_tubes(sl, line_width, **kwargs):
line_actor = actor.line(sl, **kwargs)
line_actor.GetProperty().SetRenderLinesAsTubes(1)
line_actor.GetProperty().SetLineWidth(line_width)
return line_actor
arc_actor = lines_as_tubes(arc_t1w, 8)
cst_actor = lines_as_tubes(cst_t1w, 8)
#############################################################################
# Slicer actors
# -------------
# The anatomical image is rendered using `slicer` actors. These are actors that
# visualize one slice of a three dimensional volume. Again, we create a helper
# function that will slice a volume along the x, y, and z dimensions. This
# function returns a list of the slicers we want to include in our
# visualization. This can be one, two, or three slicers, depending on how many
# of {x,y,z} are set. If you are curious to understand what is going on in this
# function, take a look at the documentation for the
# :met:`actor.slicer.display_extent` method (hint: for every dimension you
# select on, you want the full extent of the image on the two *other* two
# dimensions). We call the function on the T1-weighted data, selecting the # x
# slice that is half-way through the x dimension of the image (`shape[0]`) and
# the z slice that is a third of a way through that x dimension of the image
# (`shape[-1]`).
def slice_volume(data, x=None, y=None, z=None):
slicer_actors = []
slicer_actor_z = actor.slicer(data)
if z is not None:
slicer_actor_z.display_extent(
0, data.shape[0] - 1,
0, data.shape[1] - 1,
z, z)
slicer_actors.append(slicer_actor_z)
if y is not None:
slicer_actor_y = slicer_actor_z.copy()
slicer_actor_y.display_extent(
0, data.shape[0] - 1,
y, y,
0, data.shape[2] - 1)
slicer_actors.append(slicer_actor_y)
if x is not None:
slicer_actor_x = slicer_actor_z.copy()
slicer_actor_x.display_extent(
x, x,
0, data.shape[1] - 1,
0, data.shape[2] - 1)
slicer_actors.append(slicer_actor_x)
return slicer_actors
slicers = slice_volume(t1w, x=t1w.shape[0] // 2, z=t1w.shape[-1] // 3)
#############################################################################
# Making a `scene`
# -----------------
# The next kind of fury object we will be working with is a `window.Scene`
# object. This is the (3D!) canvas on which we are drawing the actors. We
# initialize this object and call the `scene.add` method to add the actors.
scene = window.Scene()
scene.add(arc_actor)
scene.add(cst_actor)
for slicer in slicers:
scene.add(slicer)
#############################################################################
# Showing the visualization
# -------------------------
# If you are working in an interactive session, you can call::
#
# window.show(scene, size=(1200, 1200), reset_camera=False)
#
# to see what the visualization looks like. This would pop up a window that will
# show you the visualization as it is now. You can interact with this
# visualization using your mouse to rotate the image in 3D, and mouse+ctrl or
# mouse+shift to pan and rotate it in plane, respectively. Use the scroll up and
# scroll down in your mouse to zoom in and out. Once you have found a view of
# the data that you like, you can close the window (as long as its open, it is
# blocking execution of any further commands in the Python interpreter!) and
# then you can query your scene for the "camera settings" by calling::
#
# scene.camera_info()
#
# This will print out to the screen something like this::
#
# # Active Camera
# Position (238.04, 174.48, 143.04)
# Focal Point (96.32, 110.34, 84.48)
# View Up (-0.33, -0.12, 0.94)
#
# We can use the information we have gleaned to set the camera on subsequent
# visualization that use this scene object.
scene.set_camera(position=(238.04, 174.48, 143.04),
focal_point=(96.32, 110.34, 84.48),
view_up=(-0.33, -0.12, 0.94))
#############################################################################
# Record the visualization
# -------------------------
# If you are not working in an interactive session, or you have already figured
# out the camera settings that work for you, you can go ahead and record the
# image into a png file. We use a pretty high resolution here (2400 by 2400) so
# that we get a nice crisp image. That also means the file is pretty large.
out_folder = op.join(afd.afq_home, "VizExample")
os.makedirs(out_folder, exist_ok=True)
window.record(
scene,
out_path=op.join(out_folder, 'arc_cst1.png'),
size=(2400, 2400))
############################################################################
# Setting bundle colors
# ---------------------
# The default behavior produces a nice image! But we often want to be able to
# differentiate streamlines based on the bundle to which they belong. For this,
# we will set the color of each streamline, based on the bundle. We often use
# the Tableau20 color palette to set the colors for the different bundles. The
# `actor.line` initializer takes `colors` as a keyword argument and one of the
# options to pass here is an RGB triplet, which will determine the color of all
# of the streamlines in that actor. We can get the Tableau20 RGB triplets from
# Matplotlib's colormap module. We initialize our bundle actors again and clear
# the scene before adding these new actors and adding back the slice actors. We
# then call `record` to create this new figure.
from matplotlib.cm import tab20
color_arc = tab20.colors[18]
color_cst = tab20.colors[2]
arc_actor = lines_as_tubes(arc_t1w, 8, colors=color_arc)
cst_actor = lines_as_tubes(cst_t1w, 8, colors=color_cst)
scene.clear()
scene.add(arc_actor)
scene.add(cst_actor)
for slicer in slicers:
scene.add(slicer)
window.record(
scene,
out_path=op.join(out_folder, 'arc_cst2.png'),
size=(2400, 2400))
#############################################################################
# Adding core bundles with tract profiles
# ---------------------------------------
# Next, we'd like to add a representation of the core of each bundle and plot
# the profile of some property along the length of the bundle. Another option
# for setting streamline colors is for each point along a streamline to have a a
# different color. Here, we will do this with only one streamline, which
# represents the core of the fiber -- it's median trajectory, and with a profile
# of the FA. This requires a few distinct steps. The first is adding another
# actor for each of the core fibers. We determine the core bundle as the median
# of the coordinates of the streamlines after each streamline is resampled to
# 100 equi-distant points. In the next step, we extract the FA profiles for each
# of the bundles, using the `afq_profile` function. Before we do this, we make
# sure that the streamlines are back in the voxel coordinate frame of the FA
# volume. The core bundle and FA profile are put together into a new line actor,
# which we initialize as a thick tube (line width of 40) and with a call to
# create a colormap from the FA profiles (we can choose which color-map to use
# here. In this case, we chose `'viridis'`). Note that because we are passing
# only one streamline into the line actor, we have to pass it inside of square
# brackets (`[]`). This is because the `actor.line` initializer expects a
# sequence of streamlines as input. Before plotting things again, we initialize
# our bundle line actors again. This time, we set the opacity of the lines to
# 0.1, so that they do not occlude the view of the core bundle visualizations.
#
# .. note::
# In this case, the profile is FA, but you can use any sequence of 100
# numbers in place of the FA profiles: group differences, p-values, etc.
#
from dipy.tracking.streamline import set_number_of_points
core_arc = np.median(np.asarray(set_number_of_points(arc_t1w, 100)), axis=0)
core_cst = np.median(np.asarray(set_number_of_points(cst_t1w, 100)), axis=0)
from dipy.stats.analysis import afq_profile
sft_arc.to_vox()
sft_cst.to_vox()
arc_profile = afq_profile(fa, sft_arc.streamlines, affine=np.eye(4))
cst_profile = afq_profile(fa, sft_cst.streamlines, affine=np.eye(4))
core_arc_actor = lines_as_tubes(
[core_arc],
40,
colors=create_colormap(arc_profile, 'viridis')
)
core_cst_actor = lines_as_tubes(
[core_cst],
40,
colors=create_colormap(cst_profile, 'viridis')
)
scene.clear()
arc_actor = lines_as_tubes(arc_t1w, 8, colors=color_arc, opacity=0.1)
cst_actor = lines_as_tubes(cst_t1w, 8, colors=color_cst, opacity=0.1)
scene.add(arc_actor)
scene.add(cst_actor)
for slicer in slicers:
scene.add(slicer)
scene.add(core_arc_actor)
scene.add(core_cst_actor)
window.record(
scene,
out_path=op.join(out_folder, 'arc_cst3.png'),
size=(2400, 2400))
#############################################################################
# Adding ROIs
# -----------
# Another element that we can add into the visualization are volume renderings
# of regions of interest. For example, the waypoint ROIs that were used to
# define the bundles. Here, we will add the waypoints that were used to define
# the arcuate fasciculus in this subject. We start by reading this data in.
# Because it is represented in the space of the diffusion data, it too needs to
# be resampled into the T1-weighted space, before being visualized. After
# resampling, we might have some values between 0 and 1 because of the
# interpolation from the low resolution of the diffusion into the high
# resolution of the T1-weighted. We will include in the volume rendering any
# values larger than 0. The main change from the previous visualizations is the
# adition of a `contour_from_roi` actor for each of the ROIs. We select another
# color from the Tableau 20 palette to represent this, and use an opacity of
# 0.5.
#
# .. note::
# In this case, the surface rendered is of the waypoint ROIs, but very
# similar code can be used to render other surfaces, provided a volume that
# contains that surface. For example, the gray-white matter boundary could
# be visualized provided an array with a binary representation of the volume # enclosed in this boundary
from dipy.align import resample
waypoint1 = nib.load(
op.join(
afq_path,
"ROIs", "sub-NDARAA948VFH_ses-HBNsiteRU_acq-64dir_space-T1w_desc-preproc_dwi_desc-ROI-ARC_L-1-include.nii.gz"))
waypoint2 = nib.load(
op.join(
afq_path,
"ROIs", "sub-NDARAA948VFH_ses-HBNsiteRU_acq-64dir_space-T1w_desc-preproc_dwi_desc-ROI-ARC_L-2-include.nii.gz"))
waypoint1_xform = resample(waypoint1, t1w_img)
waypoint2_xform = resample(waypoint2, t1w_img)
waypoint1_data = waypoint1_xform.get_fdata() > 0
waypoint2_data = waypoint2_xform.get_fdata() > 0
scene.clear()
arc_actor = lines_as_tubes(arc_t1w, 8, colors=color_arc)
cst_actor = lines_as_tubes(cst_t1w, 8, colors=color_cst)
scene.add(arc_actor)
scene.add(cst_actor)
for slicer in slicers:
scene.add(slicer)
surface_color = tab20.colors[0]
waypoint1_actor = actor.contour_from_roi(waypoint1_data,
color=surface_color,
opacity=0.5)
waypoint2_actor = actor.contour_from_roi(waypoint2_data,
color=surface_color,
opacity=0.5)
scene.add(waypoint1_actor)
scene.add(waypoint2_actor)
window.record(
scene,
out_path=op.join(out_folder, 'arc_cst4.png'),
size=(2400, 2400))
#############################################################################
# Making a Figure out of many fury panels
# ---------------------------------------
# We can also make a figure that contains multiple panels, each of which
# contains a different visualization. This is useful for communicating the
# results of an analysis. Here, we will make a figure with four panels, using
# some of the visualizations we have already created. We will use some
# convenient methods from pyAFQ.
pf = PanelFigure(3, 2, 6, 9)
pf.add_img(op.join(out_folder, 'arc_cst1.png'), 0, 0)
pf.add_img(op.join(out_folder, 'arc_cst2.png'), 1, 0)
pf.add_img(op.join(out_folder, 'arc_cst3.png'), 0, 1)
pf.add_img(op.join(out_folder, 'arc_cst4.png'), 1, 1)
pf.format_and_save_figure(f"arc_cst_fig.png")
#############################################################################
#
# .. note::
# If a virtual buffer was started before, it's a good idea to stop it.
if os.environ.get("XVFB", False):
print("Stopping XVFB")
vdisplay.stop()
#############################################################################
# References
# ----------
# .. [1] Garyfallidis et al., (2021). FURY: advanced scientific visualization.
# Journal of Open Source Software, 6(64), 3384,
# https://doi.org/10.21105/joss.03384
#
# .. [2] Alexander LM, Escalera J, Ai L, et al. An open resource for
# transdiagnostic research in pediatric mental health and learning
# disorders. Sci Data. 2017;4:170181.
#
# .. [3] Richie-Halford A, Cieslak M, Ai L, et al. An analysis-ready and
# quality controlled resource for pediatric brain white-matter research.
# Scientific Data. 2022;9(1):1-27.