1"""
2=========================
3Getting started with DIPY
4=========================
5
6In diffusion MRI (dMRI) usually we use three types of files, a Nifti file with
7the diffusion weighted data, and two text files one with b-values and
8one with the b-vectors.
9
10In DIPY_ we provide tools to load and process these files and we also provide
11access to publicly available datasets for those who haven't acquired yet
12their own datasets.
13
14With the following commands we can download a dMRI dataset
15"""
16
17from dipy.data import fetch_sherbrooke_3shell
18fetch_sherbrooke_3shell()
19
20"""
21By default these datasets will go in the ``.dipy`` folder inside your home
22directory. Here is how you can access them.
23"""
24
25from os.path import expanduser, join
26home = expanduser('~')
27
28"""
29``dname`` holds the directory name where the 3 files are in.
30"""
31
32dname = join(home, '.dipy', 'sherbrooke_3shell')
33
34"""
35Here, we show the complete filenames of the 3 files
36"""
37
38fdwi = join(dname, 'HARDI193.nii.gz')
39
40print(fdwi)
41
42fbval = join(dname, 'HARDI193.bval')
43
44print(fbval)
45
46fbvec = join(dname, 'HARDI193.bvec')
47
48print(fbvec)
49
50"""
51``/home/username/.dipy/sherbrooke_3shell/HARDI193.nii.gz``
52
53``/home/username/.dipy/sherbrooke_3shell/HARDI193.bval``
54
55``/home/username/.dipy/sherbrooke_3shell/HARDI193.bvec``
56
57Now, that we have their filenames we can start checking what these look like.
58
59Let's start first by loading the dMRI datasets. For this purpose, we
60use a python library called nibabel_ which enables us to read and write
61neuroimaging-specific file formats.
62"""
63
64from dipy.io.image import load_nifti
65data, affine, img = load_nifti(fdwi, return_img=True)
66
67"""
68``data`` is a 4D array where the first 3 dimensions are the i, j, k voxel
69coordinates and the last dimension is the number of non-weighted (S0s) and
70diffusion-weighted volumes.
71
72We can very easily check the size of ``data`` in the following way:
73"""
74
75print(data.shape)
76
77"""
78``(128, 128, 60, 193)``
79
80We can also check the dimensions of each voxel in the following way:
81"""
82
83print(img.header.get_zooms()[:3])
84
85"""
86``(2.0, 2.0, 2.0)``
87
88We can quickly visualize the results using matplotlib_. For example,
89let's show here the middle axial slices of volume 0 and volume 10.
90"""
91
92import matplotlib.pyplot as plt
93
94axial_middle = data.shape[2] // 2
95plt.figure('Showing the datasets')
96plt.subplot(1, 2, 1).set_axis_off()
97plt.imshow(data[:, :, axial_middle, 0].T, cmap='gray', origin='lower')
98plt.subplot(1, 2, 2).set_axis_off()
99plt.imshow(data[:, :, axial_middle, 10].T, cmap='gray', origin='lower')
100plt.show()
101plt.savefig('data.png', bbox_inches='tight')
102
103"""
104.. figure:: data.png
105   :align: center
106
107   Showing the middle axial slice without (left) and with (right) diffusion
108   weighting.
109
110The next step is to load the b-values and b-vectors from the disk using
111the function ``read_bvals_bvecs``.
112"""
113
114from dipy.io import read_bvals_bvecs
115bvals, bvecs = read_bvals_bvecs(fbval, fbvec)
116
117"""
118In DIPY, we use an object called ``GradientTable`` which holds all the
119acquisition specific parameters, e.g. b-values, b-vectors, timings and others.
120To create this object you can use the function ``gradient_table``.
121"""
122
123from dipy.core.gradients import gradient_table
124gtab = gradient_table(bvals, bvecs)
125
126"""
127Finally, you can use ``gtab`` (the GradientTable object) to show some
128information about the acquisition parameters
129"""
130
131print(gtab.info)
132
133"""
134B-values shape (193,)
135         min 0.000000
136         max 3500.000000
137B-vectors shape (193, 3)
138         min -0.964050
139         max 0.999992
140
141You can also see the b-values using:
142"""
143
144print(gtab.bvals)
145
146"""
147
148::
149
150    [    0.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.
151      1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.
152      1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.
153      1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.
154      1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.
155      1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.  1000.
156      1000.  1000.  1000.  1000.  1000.  2000.  2000.  2000.  2000.  2000.
157      2000.  2000.  2000.  2000.  2000.  2000.  2000.  2000.  2000.  2000.
158      2000.  2000.  2000.  2000.  2000.  2000.  2000.  2000.  2000.  2000.
159      2000.  2000.  2000.  2000.  2000.  2000.  2000.  2000.  2000.  2000.
160      2000.  2000.  2000.  2000.  2000.  2000.  2000.  2000.  2000.  2000.
161      2000.  2000.  2000.  2000.  2000.  2000.  2000.  2000.  2000.  2000.
162      2000.  2000.  2000.  2000.  2000.  2000.  2000.  2000.  2000.  3500.
163      3500.  3500.  3500.  3500.  3500.  3500.  3500.  3500.  3500.  3500.
164      3500.  3500.  3500.  3500.  3500.  3500.  3500.  3500.  3500.  3500.
165      3500.  3500.  3500.  3500.  3500.  3500.  3500.  3500.  3500.  3500.
166      3500.  3500.  3500.  3500.  3500.  3500.  3500.  3500.  3500.  3500.
167      3500.  3500.  3500.  3500.  3500.  3500.  3500.  3500.  3500.  3500.
168      3500.  3500.  3500.  3500.  3500.  3500.  3500.  3500.  3500.  3500.
169      3500.  3500.  3500.]
170
171Or, for example the 10 first b-vectors using:
172"""
173
174print(gtab.bvecs[:10, :])
175
176"""
177
178::
179
180    array([[ 0.        ,  0.        ,  0.        ],
181           [ 0.999979  , -0.00504001, -0.00402795],
182           [ 0.        ,  0.999992  , -0.00398794],
183           [-0.0257055 ,  0.653861  , -0.756178  ],
184           [ 0.589518  , -0.769236  , -0.246462  ],
185           [-0.235785  , -0.529095  , -0.815147  ],
186           [-0.893578  , -0.263559  , -0.363394  ],
187           [ 0.79784   ,  0.133726  , -0.587851  ],
188           [ 0.232937  ,  0.931884  , -0.278087  ],
189           [ 0.93672   ,  0.144139  , -0.31903   ]])
190
191``gtab``  can be used to tell what part of the data is the S0 volumes
192(volumes which correspond to b-values of 0).
193"""
194
195S0s = data[:, :, :, gtab.b0s_mask]
196
197"""
198Here, we had only 1 S0 as we can verify by looking at the dimensions of S0s
199"""
200
201print(S0s.shape)
202
203"""
204``(128, 128, 60, 1)``
205
206Just, for fun let's save this in a new Nifti file.
207"""
208
209from dipy.io.image import save_nifti
210save_nifti('HARDI193_S0.nii.gz', S0s, affine)
211
212"""
213Now, that we learned how to load dMRI datasets we can start the analysis.
214See example :ref:`example_reconst_dti` to learn how to create FA maps.
215
216.. include:: ../links_names.inc
217
218"""
219