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
|
"""
=========================
Getting started with DIPY
=========================
In diffusion MRI (dMRI) usually we use three types of files, a Nifti file with
the diffusion weighted data, and two text files one with b-values and
one with the b-vectors.
In DIPY_ we provide tools to load and process these files and we also provide
access to publicly available datasets for those who haven't acquired yet
their own datasets.
Let's start with some necessary imports.
"""
from os.path import expanduser, join
import matplotlib.pyplot as plt
from dipy.core.gradients import gradient_table
from dipy.data import fetch_sherbrooke_3shell
from dipy.io import read_bvals_bvecs
from dipy.io.image import load_nifti, save_nifti
###############################################################################
# With the following commands we can download a dMRI dataset
fetch_sherbrooke_3shell()
###############################################################################
# By default these datasets will go in the ``.dipy`` folder inside your home
# directory. Here is how you can access them.
home = expanduser("~")
###############################################################################
# ``dname`` holds the directory name where the 3 files are in.
dname = join(home, ".dipy", "sherbrooke_3shell")
###############################################################################
# Here, we show the complete filenames of the 3 files
fdwi = join(dname, "HARDI193.nii.gz")
print(fdwi)
fbval = join(dname, "HARDI193.bval")
print(fbval)
fbvec = join(dname, "HARDI193.bvec")
print(fbvec)
###############################################################################
# Now, that we have their filenames we can start checking what these look like.
#
# Let's start first by loading the dMRI datasets. For this purpose, we
# use a python library called nibabel_ which enables us to read and write
# neuroimaging-specific file formats.
data, affine, img = load_nifti(fdwi, return_img=True)
###############################################################################
# ``data`` is a 4D array where the first 3 dimensions are the i, j, k voxel
# coordinates and the last dimension is the number of non-weighted (S0s) and
# diffusion-weighted volumes.
#
# We can very easily check the size of ``data`` in the following way:
print(data.shape)
###############################################################################
# We can also check the dimensions of each voxel in the following way:
print(img.header.get_zooms()[:3])
###############################################################################
# We can quickly visualize the results using matplotlib_. For example,
# let's show here the middle axial slices of volume 0 and volume 10.
axial_middle = data.shape[2] // 2
plt.figure("Showing the datasets")
plt.subplot(1, 2, 1).set_axis_off()
plt.imshow(data[:, :, axial_middle, 0].T, cmap="gray", origin="lower")
plt.subplot(1, 2, 2).set_axis_off()
plt.imshow(data[:, :, axial_middle, 10].T, cmap="gray", origin="lower")
plt.show()
plt.savefig("data.png", bbox_inches="tight")
###############################################################################
# .. rst-class:: centered small fst-italic fw-semibold
#
# Showing the middle axial slice without (left) and with (right) diffusion
# weighting.
#
#
#
# The next step is to load the b-values and b-vectors from the disk using
# the function ``read_bvals_bvecs``.
bvals, bvecs = read_bvals_bvecs(fbval, fbvec)
###############################################################################
# In DIPY, we use an object called ``GradientTable`` which holds all the
# acquisition specific parameters, e.g. b-values, b-vectors, timings and
# others. To create this object you can use the function ``gradient_table``.
gtab = gradient_table(bvals, bvecs=bvecs)
###############################################################################
# Finally, you can use ``gtab`` (the GradientTable object) to show some
# information about the acquisition parameters
print(gtab.info)
###############################################################################
# You can also see the b-values using:
print(gtab.bvals)
###############################################################################
# Or, for example the 10 first b-vectors using:
print(gtab.bvecs[:10, :])
###############################################################################
# You can get the number of gradients (including the number of b0 values)
# calling ``len`` on the ``GradientTable`` instance:
print(len(gtab))
###############################################################################
# ``gtab`` can be used to tell what part of the data is the S0 volumes
# (volumes which correspond to b-values of 0).
S0s = data[:, :, :, gtab.b0s_mask]
###############################################################################
# Here, we had only 1 S0 as we can verify by looking at the dimensions of S0s
print(S0s.shape)
###############################################################################
# Just, for fun let's save this in a new Nifti file.
save_nifti("HARDI193_S0.nii.gz", S0s, affine)
###############################################################################
# Now, that we learned how to load dMRI datasets we can start the analysis.
# See example :ref:`sphx_glr_examples_built_reconstruction_reconst_dti.py` to
# learn how to create FA maps.
|
