diff --git a/data/sub-02_dwi.bval b/data/sub-02_dwi.bval deleted file mode 100644 index 1bab7ac..0000000 --- a/data/sub-02_dwi.bval +++ /dev/null @@ -1 +0,0 @@ -5 5 610 1240 615 1225 625 1230 620 1230 610 1240 615 1225 620 1220 620 5 1225 615 1240 625 1230 620 1240 615 1225 615 1230 625 1230 610 1230 5 615 1235 615 1220 625 1245 620 1240 620 1225 615 1235 620 1225 605 5 1225 625 1235 615 1230 620 1230 610 1235 620 1230 620 1220 625 1235 5 610 1230 620 1240 615 1230 625 1230 610 1235 610 1240 620 1235 620 5 1220 615 1230 615 1230 620 1240 610 1215 620 1240 620 1225 620 1225 615 1230 5 890 1800 895 1785 910 1790 900 1795 890 1800 895 1785 905 1780 900 5 1785 895 1805 905 1795 900 1800 895 1780 895 1790 905 1790 890 1790 5 895 1795 895 1775 905 1805 900 1805 900 1785 895 1795 900 1785 885 5 1785 905 1795 895 1795 900 1790 890 1795 900 1790 900 1780 910 1795 5 890 1790 900 1805 895 1790 905 1790 890 1795 890 1800 905 1800 900 5 1780 895 1790 895 1790 900 1805 890 1775 900 1805 900 1785 900 1785 895 1790 diff --git a/data/sub-02_dwi.bvec b/data/sub-02_dwi.bvec deleted file mode 100644 index c5968be..0000000 --- a/data/sub-02_dwi.bvec +++ /dev/null @@ -1,3 +0,0 @@ -0.531695 0.531695 -0.545977 0.483628 0.691191 -0.105459 -0.294575 -0.867939 -0.957356 0.14789 0.296586 -0.263343 -0.167652 -0.81816 0.568106 0.825166 -0.707535 0.531695 0.491263 -0.697671 0.247618 -0.297736 -0.778689 -0.942033 0.368183 -0.0645985 0.450837 -0.227713 0.988353 0.633661 0.611956 0.869951 -0.394799 0.531695 0.211295 -0.79625 -0.71349 0.0577723 -0.110309 -0.12057 -0.296779 -0.638925 -0.381349 -0.362386 -0.90571 -0.164736 -0.764349 0.881573 0.582101 0.531695 -0.511588 0.235808 -0.693941 0.173357 0.958996 0.766879 0.0119803 -0.409544 0.425315 0.982639 -0.73073 0.408617 -0.317282 0.0222084 0.561343 0.531695 -0.122333 0.658063 0.839668 -0.129015 -0.51021 -0.696646 0.536149 -0.163584 0.56363 0.902852 -0.724038 0.0750007 0.476689 0.68334 0.136501 0.531695 -0.186741 -0.903589 0.524559 -0.0124602 -0.977455 -0.802906 -0.315849 0.0979143 0.363383 -0.876238 0.453462 0.0771761 0.69704 -0.546704 0.946959 0.457953 -0.113123 0.531695 -0.544936 0.483667 0.690919 -0.105761 -0.295221 -0.868069 -0.957646 0.147654 0.296574 -0.263495 -0.167325 -0.818163 0.568269 0.824569 -0.707712 0.531695 0.491245 -0.696686 0.24725 -0.297772 -0.778943 -0.941885 0.368626 -0.0641064 0.450583 -0.228784 0.988407 0.633947 0.611876 0.86899 -0.394911 0.531695 0.211367 -0.796546 -0.713488 0.0576796 -0.110199 -0.120379 -0.296724 -0.639359 -0.381307 -0.362002 -0.904865 -0.164665 -0.765339 0.881319 0.581293 0.531695 -0.511303 0.235777 -0.694147 0.173465 0.959197 0.767218 0.0121282 -0.409605 0.425303 0.982425 -0.730638 0.408316 -0.317067 0.0218063 0.561006 0.531695 -0.123679 0.65794 0.840407 -0.129021 -0.509844 -0.696571 0.536687 -0.164075 0.563762 0.903298 -0.723724 0.0749389 0.476648 0.683608 0.135636 0.531695 -0.187123 -0.903731 0.524375 -0.0130766 -0.977257 -0.803017 -0.316274 0.0977579 0.363638 -0.876749 0.45367 0.0759775 0.697124 -0.547355 0.946716 0.457396 -0.113138 -0.58283 0.58283 0.121865 0.854567 -0.512577 0.288723 -0.941474 0.445659 -0.288893 -0.594316 0.26337 0.886353 -0.748207 -0.384613 -0.770115 0.0628862 -0.472155 0.58283 0.0442891 -0.447981 -0.942525 0.934804 0.482426 0.301846 0.799217 0.702844 0.462096 0.340428 0.150971 0.762401 -0.45612 -0.0786779 0.547192 0.58283 0.149505 -0.564176 0.433314 0.183068 0.955914 -0.989788 0.734683 0.769266 -0.706198 -0.489619 -0.0967702 -0.823982 0.63408 -0.0583448 -0.0405124 0.58283 0.079565 0.931015 -0.693473 -0.34494 -0.241181 0.622141 -0.723782 -0.423991 -0.75564 -0.111409 -0.402825 -0.652493 0.251743 -0.998248 -0.746255 0.58283 0.024262 0.290716 0.479636 -0.960859 0.424062 0.274791 -0.84412 -0.515091 0.350898 0.420225 -0.0727105 0.917063 0.736002 -0.653662 0.870117 0.58283 -0.221191 -0.243925 0.702476 -0.548266 0.126745 0.550065 0.946883 0.495776 0.149623 0.39115 0.878861 -0.89406 -0.164084 -0.788454 0.094938 -0.400783 0.532024 0.58283 0.121912 0.854112 -0.511209 0.287513 -0.940903 0.444569 -0.287904 -0.593459 0.262493 0.885784 -0.74656 -0.383118 -0.769135 0.0622394 -0.471158 0.58283 0.0438219 -0.446851 -0.942305 0.934524 0.481515 0.301062 0.798308 0.701373 0.461051 0.339449 0.150393 0.761973 -0.455249 -0.0768719 0.545989 0.58283 0.148737 -0.563409 0.431571 0.182186 0.955546 -0.989721 0.733506 0.768908 -0.704729 -0.488417 -0.0955991 -0.823485 0.632655 -0.057818 -0.040301 0.58283 0.0782924 0.930443 -0.69274 -0.344176 -0.239751 0.621009 -0.72284 -0.423305 -0.755001 -0.109278 -0.401929 -0.651377 0.250936 -0.998154 -0.745804 0.58283 0.0241904 0.28941 0.477565 -0.960595 0.421882 0.274129 -0.843759 -0.513902 0.349883 0.419104 -0.0723075 0.916624 0.734631 -0.652827 0.869069 0.58283 -0.220578 -0.240917 0.70134 -0.547476 0.125888 0.549227 0.946691 0.493615 0.148325 0.389783 0.878583 -0.893047 -0.162313 -0.7874 0.0941528 -0.399917 0.531251 -0.6145 0.6145 -0.82889 -0.189258 -0.509432 0.951587 0.163865 -0.219247 0.00330035 0.790517 -0.917972 0.380827 -0.641933 -0.427418 0.290135 -0.561379 0.525798 0.6145 0.869885 -0.559078 0.224348 0.193634 0.401135 -0.146504 0.475073 -0.708405 -0.763684 0.91228 -0.019154 0.131215 0.646115 -0.486822 0.738048 0.6145 0.965921 0.218384 -0.550611 -0.981401 -0.272139 -0.0760449 0.610052 0.00228013 0.596538 -0.793064 -0.412704 0.54214 0.117104 0.468428 -0.812106 0.6145 0.855539 0.278576 -0.193757 0.922477 -0.148859 -0.1576 -0.689925 -0.807778 0.498112 -0.148356 0.551149 0.638189 -0.914307 -0.0548447 -0.357767 0.6145 0.992193 0.694577 0.254768 0.245163 0.748236 0.662703 -0.00251417 0.841381 -0.747791 0.0909314 -0.685917 0.391625 0.480696 0.325226 -0.473566 0.6145 -0.957184 0.352176 -0.481005 -0.836211 -0.168868 -0.229719 0.060429 -0.862913 -0.919546 0.28144 0.148242 -0.441249 0.698005 0.281879 -0.30701 0.793506 0.839139 0.6145 -0.829567 -0.191203 -0.511171 0.951919 0.165973 -0.220939 0.00493084 0.791204 -0.918227 0.382044 -0.643934 -0.428754 0.292408 -0.562327 0.526454 0.6145 0.869918 -0.561206 0.225676 0.194926 0.401736 -0.149046 0.476255 -0.709906 -0.764465 0.912377 -0.0208293 0.13232 0.646804 -0.488822 0.738879 0.6145 0.966024 0.219285 -0.55198 -0.981571 -0.273475 -0.0772054 0.611493 0.000303432 0.598298 -0.79398 -0.414825 0.542916 0.11834 0.468971 -0.812696 0.6145 0.855827 0.28051 -0.195628 0.922742 -0.149867 -0.160387 -0.69091 -0.808107 0.499089 -0.151323 0.551925 0.63952 -0.914604 -0.0566807 -0.359235 0.6145 0.992027 0.695239 0.256218 0.246193 0.749717 0.663055 -0.00617345 0.842012 -0.748167 0.09167 -0.686291 0.392664 0.48283 0.326339 -0.475734 0.6145 -0.957251 0.353878 -0.482859 -0.83672 -0.170647 -0.231331 0.0612066 -0.864169 -0.919656 0.281746 0.149252 -0.443503 0.698335 0.283557 -0.308001 0.794264 0.839626 diff --git a/data/sub-02_dwi.nii.gz b/data/sub-02_dwi.nii.gz deleted file mode 100644 index 50a2489..0000000 Binary files a/data/sub-02_dwi.nii.gz and /dev/null differ diff --git a/docs/_config.yml b/docs/_config.yml index 8fcbea1..80bb366 100644 --- a/docs/_config.yml +++ b/docs/_config.yml @@ -7,8 +7,6 @@ copyright: "2021" sphinx: config: html_extra_path: ["assets"] - extra_extensions: - - sphinx_exercise # Execution settings execute: diff --git a/docs/extra/nifti.md b/docs/extra/nifti.md index 3168571..40f9c40 100644 --- a/docs/extra/nifti.md +++ b/docs/extra/nifti.md @@ -15,7 +15,7 @@ kernelspec: ## Investigating NIfTI images with *NiBabel* [NiBabel](https://nipy.org/nibabel/) is a Python package for reading and writing neuroimaging data. -To learn more about how NiBabel handles NIfTIs, check out the [Working with NIfTI images](https://nipy.org/nibabel/nifti_images.html) page of the NiBabel documentation. +To learn more about how NiBabel handles NIfTIs, check out the [Working with NIfTI images](https://nipy.org/nibabel/nifti_images.html) page of the NiBabel documentation, from which this lesson is adapted from. ```{code-cell} python import nibabel as nib @@ -105,178 +105,25 @@ dwi_affine = dwi_img.affine dwi_affine ``` -To explain this concept, recall that we referred to coordinates in our data as *(i,j,k)* coordinates such that: +To explain this concept, recall that we referred to coordinates in our data as *voxel coordinates (i,j,k)* coordinates such that: * i is the first dimension of `dwi_data` * j is the second dimension of `dwi_data` * k is the third dimension of `dwi_data` Although this tells us how to access our data in terms of voxels in a 3D volume, it doesn't tell us much about the actual dimensions in our data (centimetres, right or left, up or down, back or front). -The affine matrix allows us to translate between *voxel coordinates* in (i,j,k) and *world space coordinates* in (left/right,bottom/top,back/front). +The affine matrix allows us to translate between *voxel coordinates* and *world space coordinates* in (left/right,bottom/top,back/front). + An important thing to note is that in reality in which order you have: * left/right * bottom/top * back/front -Depends on how you've constructed the affine matrix, but for the data we're dealing with it always refers to: +depends on how the affine matrix is constructed. For the data we're dealing with, it always refers to: * Right * Anterior * Superior -Applying the affine matrix is done through using a *linear map* (matrix multiplication) on voxel coordinates (defined in `dwi_data`). - -## Diffusion gradient schemes - -In addition to the acquired diffusion images, two files are collected as part of the diffusion dataset. -These files correspond to the gradient amplitude (b-values) and directions (b-vectors) of the diffusion measurement and are named with the extensions `.bval` and `.bvec` respectively. - -```{code-cell} python -bvec = "../../data/sub-01_dwi.bvec" -bval = "../../data/sub-01_dwi.bval" -``` - -The b-value is the diffusion-sensitizing factor, and reflects both the timing & strength of the gradients (measured in s/mm^2) used to acquire the diffusion-weighted images. - -```{code-cell} python -!cat ../../data/sub-01_dwi.bval -``` - -The b-vector corresponds to the direction of the diffusion sensitivity. Each row corresponds to a value in the i, j, or k axis. The numbers are combined column-wise to get an [i j k] coordinate per DWI volume. - -```{code-cell} python -!cat ../../data/sub-01_dwi.bvec -``` - -Together these two files define the dMRI measurement as a set of gradient directions and corresponding amplitudes. - -In our example data, we see that 2 b-values were chosen for this scanning sequence. -The first few images were acquired with a b-value of 0 and are typically referred to as b=0 images. -In these images, no diffusion gradient is applied. -These images don't hold any diffusion information and are used as a reference (for head motion correction or brain masking) since they aren't subject to the same types of scanner artifacts that affect diffusion-weighted images. - -All of the remaining images have a b-value of 1000 and have a diffusion gradient associated with them. -Diffusion that exhibits directionality in the same direction as the gradient results in a loss of signal. -With further processing, the acquired images can provide measurements which are related to the microscopic changes and estimate white matter trajectories. - - - -We'll use some functions from [Dipy](https://dipy.org), a Python package for pre-processing and analyzing diffusion data. -After reading the `.bval` and `.bvec` files with the `read_bvals_bvecs()` function, we get both in a numpy array. Notice that the `.bvec` file has been transposed so that the i, j, and k-components are in column format. - -```{code-cell} python -from dipy.io import read_bvals_bvecs - -gt_bvals, gt_bvecs = read_bvals_bvecs(bval, bvec) -gt_bvecs -``` - -Below is a plot of all of the diffusion directions that we've acquired. - -```{code-cell} python -import matplotlib.pyplot as plt - -fig = plt.figure() -ax = fig.add_subplot(111, projection="3d") -ax.scatter(gt_bvecs.T[0], gt_bvecs.T[1], gt_bvecs.T[2]) -plt.show() -``` - -It is important to note that in this format, the diffusion gradients are provided with respect to the image axes, not in real or scanner coordinates. Simply reformatting the image from sagittal to axial will effectively rotate the b-vectors, since this operation changes the image axes. Thus, a particular bvals/bvecs pair is only valid for the particular image that it corresponds to. - -The diffusion gradient is critical for later analyzing the data - -```{code-cell} python -:tags: [output_scroll] - -import numpy as np - -if dwi_affine.shape == (4, 4): - dwi_affine = dwi_affine[:3, :3] - -rotated_bvecs = dwi_affine[np.newaxis, ...].dot(gt_bvecs.T)[0].T -rotated_bvecs -``` - -```{code-cell} python - -fig = plt.figure() -ax = fig.add_subplot(111, projection="3d") -ax.scatter(rotated_bvecs.T[0], rotated_bvecs.T[1], rotated_bvecs.T[2]) -plt.show() -``` - -Inspired by MRtrix3 and proposed in the [BIDS spec](https://github.com/bids-standard/bids-specification/issues/349), dMRIPrep also creates an optional `.tsv` file where the diffusion gradients are reported in scanner coordinates as opposed to image coordinates. -The [i j k] values reported earlier are recalculated in [R A S]. - -```{code-cell} python -:tags: [output_scroll] - -rasb = np.c_[rotated_bvecs, gt_bvals] - -rasb -``` - -We can write out this `.tsv` to a file. - -```{code-cell} python -np.savetxt(fname="../../data/sub-01_rasb.tsv", delimiter="\t", X=rasb) -``` - -```{code-cell} python -from dipy.core.gradients import gradient_table - -gtab = gradient_table(gt_bvals, rotated_bvecs) -``` - -## Brain Masking - -One of the first things we do before image registration is brain extraction, separating any non-brain material from brain tissue. -This is done so that our algorithms aren't biased or distracted by whatever is in non-brain material and we don't spend extra time analyzing things we don't care about - -As mentioned before, the b0s are a good reference scan for doing brain masking. lets index them. - -```{code-cell} python -gtab.b0s_mask -``` - -```{code-cell} python -bzero = dwi_data[:, :, :, gtab.b0s_mask] -``` - -skullstrip -5 volumes in b0 - -```{code-cell} python -bzero.shape -``` - -take median image - -```{code-cell} python - -median_bzero = np.median(bzero, axis=-1) -``` - -```{code-cell} python -from dipy.segment.mask import median_otsu - -b0_mask, mask = median_otsu(median_bzero, median_radius=2, numpass=1) -``` - -```{code-cell} python -from dipy.core.histeq import histeq - -sli = median_bzero.shape[2] // 2 -plt.subplot(1, 3, 1).set_axis_off() -plt.imshow(histeq(median_bzero[:, :, sli].astype("float")).T, - cmap="gray", origin="lower") - -plt.subplot(1, 3, 2).set_axis_off() -plt.imshow(mask[:, :, sli].astype("float").T, cmap="gray", origin="lower") - -plt.subplot(1, 3, 3).set_axis_off() -plt.imshow(histeq(b0_mask[:, :, sli].astype("float")).T, - cmap="gray", origin="lower") -``` +Applying the affine matrix is done through using a *linear map* (matrix multiplication) on the voxel coordinates (defined in `dwi_data`). diff --git a/docs/nipreps/dmriprep.md b/docs/nipreps/dmriprep.md index f142209..d98336f 100644 --- a/docs/nipreps/dmriprep.md +++ b/docs/nipreps/dmriprep.md @@ -6,7 +6,7 @@ dMRIPrep is an analysis-agnostic tool that addresses the challenge of robust and ## Development -In his 2019 ISMRM talk [^veraart2019], Jelle Veraart polled the developers of some of the major dMRI analysis software packages. +In his 2019 ISMRM talk [^veraart2019], Dr. J. Veraart polled the developers of some of the major dMRI analysis software packages. The goal was to understand how much consensus there was in the field on whether to proceed with certain pre-processing steps. The poll showed that consensus is within reach. However, different pre-processing steps have varying levels of support. diff --git a/docs/nipreps/nipreps.md b/docs/nipreps/nipreps.md index b8f0204..3273dc3 100644 --- a/docs/nipreps/nipreps.md +++ b/docs/nipreps/nipreps.md @@ -83,7 +83,7 @@ This eases the burden of maintaining these tools but also helps focus on standar *NiPreps* only support BIDS-Derivatives as output and so are agnostic to subsequent analysis. *NiPreps* also aim to be robust in their codebase. -The pipelines are modular and rely on widely-used tools such as AFNI, ANTs, FreeSurfer, FSL, Nilearn, or Dipy and are extensible via plug-ins. +The pipelines are modular and rely on widely-used tools such as AFNI, ANTs, FreeSurfer, FSL, Nilearn, or DIPY and are extensible via plug-ins. This modularity in the code allows each step to be thoroughly tested. Some examples of tests performed on different parts of the pipeline are shown below: ```{tabbed} unittest diff --git a/docs/references.bib b/docs/references.bib deleted file mode 100644 index a845151..0000000 --- a/docs/references.bib +++ /dev/null @@ -1,2 +0,0 @@ ---- ---- diff --git a/docs/tutorial/data.md b/docs/tutorial/data.md index ed06174..a93adf3 100644 --- a/docs/tutorial/data.md +++ b/docs/tutorial/data.md @@ -16,43 +16,49 @@ kernelspec: :tags: [hide-cell] import warnings -from eddymotion.viz import plot_dwi warnings.filterwarnings("ignore") - -def _data_repr(value): - if value is None: - return "None" - return f"<{'x'.join(str(v) for v in value.shape)} ({value.dtype})>" - ``` -Diffusion imaging probes the random, microscopic motion of water protons by using MRI sequences that are sensitive to the geometry and environmental organization surrounding these protons. +Diffusion imaging probes the random, microscopic movement of water molecules by using MRI sequences that are sensitive to the geometry and environmental organization surrounding these protons. This is a popular technique for studying the white matter of the brain. The diffusion within biological structures, such as the brain, are often restricted due to barriers (eg. cell membranes), resulting in a preferred direction of diffusion (anisotropy). A typical dMRI scan will acquire multiple volumes (or ***angular samples***), each sensitive to a particular ***diffusion direction***. + + + +*Sourced from Dr. A. Rokem, DIPY Workshop 2021* + These *diffusion directions* (or ***orientations***) are a fundamental piece of metadata to interpret dMRI data, as models need to know the exact orientation of each angular sample. ```{admonition} Main elements of a dMRI dataset - A 4D data array, where the last dimension encodes the reconstructed **diffusion direction *maps***. -- Tabular data or a 2D array, listing the **diffusion directions** and the encoding gradient strength. +- Tabular data or a 2D array, listing the **diffusion directions** (`.bvec`) and the encoding **gradient strength** (`.bval`). ``` In summary, dMRI involves ***complex data types*** that, as programmers, we want to access, query and manipulate with ease. ## Python and object oriented programming -Python is an [object oriented programming](https://en.wikipedia.org/wiki/Object-oriented_programming) language, which represent and encapsulate data types and corresponding behaviors into programming structures called *objects*. +Python is an [object oriented programming](https://en.wikipedia.org/wiki/Object-oriented_programming) language. +It allows us to represent and encapsulate data types and corresponding behaviors into programming structures called *objects*. Therefore, let's leverage Python to create *objects* that contain dMRI data. -In Python, *objects* can be specified by defining a class with name `DWI`. -To simplify class creation, we'll use the magic of a Python library called [`attrs`](https://www.attrs.org/en/stable/). +In Python, *objects* can be specified by defining a class. +In the example code below, we've created a class with the name `DWI`. +To simplify class creation, we've also used the magic of a Python library called [`attrs`](https://www.attrs.org/en/stable/). ```{code-cell} python + """Representing data in hard-disk and memory.""" import attr +def _data_repr(value): + if value is None: + return "None" + return f"<{'x'.join(str(v) for v in value.shape)} ({value.dtype})>" + @attr.s(slots=True) class DWI: @@ -78,147 +84,229 @@ class DWI: ``` -This first code implements several *attributes* and the first *behavior* - the `__len__` *method*. +This code implements several *attributes* as well as a *behavior* - the `__len__` *method*. The `__len__` method is special in Python, as it will be executed when we call the built-in function `len()` on our object. -Let's test this memory structure with some *simulated* data: + +Let's test out the `DWI` data structure with some *simulated* data: ```{code-cell} python -# NumPy is a fundamental Python library + +# NumPy is a fundamental Python library for working with arrays import numpy as np -# Let's create a new DWI object, with only gradient information that is random +# create a new DWI object, with only gradient information that is random dmri_dataset = DWI(gradients=np.random.normal(size=(4, 64))) -# Let's call Python's built-in len() function +# call Python's built-in len() function print(len(dmri_dataset)) + ``` -The output of this `print()` statement is telling us that this (simulated) dataset has 64 diffusion weighted samples. +The output of this `print()` statement is telling us that this (simulated) dataset has 64 diffusion-weighted samples. ## Using the new data representation object -For simplicity, we will be using the full implementation of the `DWI` class from our [`eddymotion` package](https://github.com/nipreps/EddyMotionCorrection/blob/57c518929146b23cc9534ab0b2d024aa136e25f8/emc/dmri.py) +The code shown above was just a snippet of the `DWI` class. For simplicity, we will be using the full implementation of this class from our [`eddymotion` package](https://github.com/nipreps/EddyMotionCorrection/blob/main/eddymotion/dmri.py) Under the `data/` folder of this book's distribution, we have stored a sample DWI dataset with filename `dwi.h5`. -Please note that the file has been minimized by zeroing all but two diffusion weighted orientation maps. +Please note that the file has been minimized by zeroing all but two diffusion-weighted orientation maps. + Let's get some insights from it: ```{code-cell} python -# Import the class from the library + +# import the class from the library from eddymotion.dmri import DWI -# Load the sample file +# load the sample file dmri_dataset = DWI.from_filename("../../data/dwi.h5") print(len(dmri_dataset)) + ``` -In this case, the dataset is reporting to have 102 diffusion weighted samples. +In this case, the dataset is reporting to have 102 diffusion-weighted samples. -Python will automatically generate a summary of this object if we just issue the name of our new object. -This pretty-printing of the object informs us about the data and metadata that together compose this particular DWI dataset: +Python will automatically generate a summary of this object if we just type the name of our new object. +This pretty-printing of the object informs us about the data and metadata that, together, compose this particular DWI dataset: ```{code-cell} python + dmri_dataset + ``` -### Visualizing the data +We'll go over some of the components of `dmri_dataset` through this lesson. + +## Visualizing the data +Let's start out by seeing what the data looks like. The fully-fledged `DWI` object has a convenience function to plot the dataset: ```{code-cell} python -dmri_dataset.plot_mosaic(); + +dmri_dataset.plot_mosaic() + ``` -When calling `plot_mosaic()` without arguments, the *b=0* reference is plotted. +When calling `plot_mosaic()` without any arguments, the *b=0* reference is plotted. This *b=0* reference is a map of the signal measured ***without gradient sensitization***, or in other words, when we are not measuring diffusion in any direction. The *b=0* map can be used by diffusion modeling as the reference to quantify the signal drop at every voxel and given a particular orientation gradient. -We can also get an insight of how a diffusion weighted orientation looks like by selecting them with the argument `index`: +We can also get some insight into how a particular diffusion-weighted orientation looks like by selecting them with the argument `index`. + +```{admonition} Exercise +Try calling `plot_mosaic` with an index of 10 or 100. +``` + +**Solution** ```{code-cell} python -dmri_dataset.plot_mosaic(index=10); +:tags: [hide-cell] + +dmri_dataset.plot_mosaic(index=10, vmax=5000) + ``` or: ```{code-cell} python -dmri_dataset.plot_mosaic(index=100); +:tags: [hide-cell] + +dmri_dataset.plot_mosaic(index=100, vmax=5000) + ``` -As we can see, ***diffusion weighted*** images consistently drop almost all signal in voxels filled with cerebrospinal fluid because there, water-diffusion is free regardless of the direction you are measuring in. +Diffusion that exhibits directionality in the same direction as the gradient results in a loss of signal. +As we can see, ***diffusion-weighted*** images consistently drop almost all signal in voxels filled with cerebrospinal fluid because there, water diffusion is free (isotropic) regardless of the direction that is being measured. + +We can also see that the images at `index=10` and `index=100` have different gradient strengths. +The higher the magnitude of the gradient, the more diffusion that is allowed to occur, indicated by the overall decrease in signal intensity. +There is also a lot more noise. +## Visualizing the gradient information -### Visualizing the gradient information +Our `DWI` object stores the gradient information in the `gradients` attribute. -Our `DWI` object stores the gradient information in the `gradients` attribute: +```{admonition} Exercise +Let's see the shape of the gradient information +``` + +**Solution** ```{code-cell} python +:tags: [hide-cell] + dmri_dataset.gradients.shape + +``` + +We get a $4\times102$. + +```{admonition} Exercise +Try printing the gradient information to see what it contains. +Remember to transpose (`.T`) the array +``` + +**Solution** + +```{code-cell} python +:tags: [hide-cell] + +print(dmri_dataset.gradients.T) + +``` + +Later, we'll refer to this array as the gradient table. + +It consists of one row per diffusion-weighted image, with each row consisting of 4 values corresponding to [ R A S+ b ]. + +[ R A S+ ] are the components of the **gradient direction**. +Note that the directions have been re-oriented with respect to *world space coordinates*. +For more information on this, refer to {doc}`../extra/nifti`. + +The last column, b, reflects the **timing and strength of the gradients** in units of s/mm². + +To get a better sense of which gradient directions were sampled, let's plot them! + +```{code-cell} python + +dmri_dataset.plot_gradients(cmap="viridis") + ``` -<<<<<<< HEAD -### The *LOGO* (leave-one-gradient-out) splitter +We've projected all of the gradient directions onto the surface of a sphere, with each unique gradient strength colour-coded. + +## The *LOGO* (leave-one-gradient-out) splitter + One final behavior that will make our endeavor easier in the long run is a convenience method for data splitting. In particular, we are implementing some sort of cross-validation scheme where we will iterate over different data splits. In this case, the splitting strategy is a simple leave-one-out. -Because one "*datapoint*" in our DW dataset corresponds to one gradient, we will refer to this partitioning of the dataset as *leave-one-gradient-out (LOGO)*: +Because one "*datapoint*" in our DWI dataset corresponds to one gradient, we will refer to this partitioning of the dataset as *leave-one-gradient-out (LOGO)*: ```{code-cell} python - def logo_split(self, index, with_b0=False): - """ - Produce one fold of LOGO (leave-one-gradient-out). - - Parameters - ---------- - index : :obj:`int` - Index of the DWI orientation to be left out in this fold. - with_b0 : :obj:`bool` - Insert the *b=0* reference at the beginning of the training dataset. - - Return - ------ - (train_data, train_gradients) : :obj:`tuple` - Training DWI and corresponding gradients. - Training data/gradients come **from the updated dataset**. - (test_data, test_gradients) :obj:`tuple` - Test 3D map (one DWI orientation) and corresponding b-vector/value. - The test data/gradient come **from the original dataset**. - - """ - dwframe = self.dataobj[..., index] - bframe = self.gradients[..., index] - - # if the size of the mask does not match data, cache is stale - mask = np.zeros(len(self), dtype=bool) - mask[index] = True - - train_data = self.dataobj[..., ~mask] - train_gradients = self.gradients[..., ~mask] - - if with_b0: - train_data = np.concatenate( - (np.asanyarray(self.bzero)[..., np.newaxis], train_data), - axis=-1, - ) - b0vec = np.zeros((4, 1)) - b0vec[0, 0] = 1 - train_gradients = np.concatenate( - (b0vec, train_gradients), - axis=-1, - ) - - return ( - (train_data, train_gradients), - (dwframe, bframe), + +def logo_split(self, index, with_b0=False): + """ + Produce one fold of LOGO (leave-one-gradient-out). + + Parameters + ---------- + index : :obj:`int` + Index of the DWI orientation to be left out in this fold. + with_b0 : :obj:`bool` + Insert the *b=0* reference at the beginning of the training dataset. + + Return + ------ + (train_data, train_gradients) : :obj:`tuple` + Training DWI and corresponding gradients. + Training data/gradients come **from the updated dataset**. + (test_data, test_gradients) :obj:`tuple` + Test 3D map (one DWI orientation) and corresponding b-vector/value. + The test data/gradient come **from the original dataset**. + + """ + dwframe = self.dataobj[..., index] + bframe = self.gradients[..., index] + + # if the size of the mask does not match data, cache is stale + mask = np.zeros(len(self), dtype=bool) + mask[index] = True + + train_data = self.dataobj[..., ~mask] + train_gradients = self.gradients[..., ~mask] + + if with_b0: + train_data = np.concatenate( + (np.asanyarray(self.bzero)[..., np.newaxis], train_data), + axis=-1, + ) + b0vec = np.zeros((4, 1)) + b0vec[0, 0] = 1 + train_gradients = np.concatenate( + (b0vec, train_gradients), + axis=-1, ) + + return ( + (train_data, train_gradients), + (dwframe, bframe), + ) + ``` -This function will allow us to easily partition the dataset as follows: +This function is contained in the `DWI` class shown earlier and will allow us to easily partition the dataset as follows: ```{code-cell} python + +from eddymotion.viz import plot_dwi + data_train, data_test = dmri_dataset.logo_split(10) -plot_dwi(data_test[0], dmri_dataset.affine, gradient=data_test[1]); +plot_dwi(data_test[0], dmri_dataset.affine, gradient=data_test[1]) + ``` -where `data_train` is a tuple containing all DW volumes and corresponding gradient table excluding the left-out, which is store in `data_test`. -Consequently, `data_test[0]` contains the held-out DW map and `data_test[1]` the corresponding gradient vector (RAS+B format). -======= ->>>>>>> fix spacing +`data_train` is a tuple containing all diffusion-weighted volumes and the corresponding gradient table, excluding the left-out, which is stored in `data_test`. +`data_test[0]` contains the held-out diffusion-weighted volume and `data_test[1]`, the corresponding gradient table. + +## Next steps: diffusion modeling + +By modeling the diffusion signal, the acquired images can provide measurements which are related to the microscopic changes and estimate white matter trajectories. diff --git a/docs/tutorial/intro.md b/docs/tutorial/intro.md index c3ad68c..2a44326 100644 --- a/docs/tutorial/intro.md +++ b/docs/tutorial/intro.md @@ -20,6 +20,6 @@ While we can address the misalignment, it is really problematic to overcome the ## Objective: Implement a head-motion estimation code This tutorial focuses on the misalignment problem. -We will build from existing software (Dipy, for diffusion modeling) and ANTs (for image registration), as well as commonplace Python libraries (NumPy) a software framework for head-motion estimation in diffusion MRI data. +We will build from existing software (Dipy for diffusion modeling and ANTs for image registration), as well as commonplace Python libraries (NumPy), a software framework for head-motion estimation in diffusion MRI data. The algorithmic and theoretical foundations of the method are described by Dr. M. Cieslak in his [OHBM 2019 poster](https://github.com/mattcieslak/ohbm_shoreline/blob/master/cieslakOHBM2019.pdf). diff --git a/docs/welcome.md b/docs/welcome.md index 4b1763e..a7daeeb 100644 --- a/docs/welcome.md +++ b/docs/welcome.md @@ -1,6 +1,6 @@ # Welcome! -## *Implementing a head-motion correction algorithm for diffusion MRI in Python, using Dipy and NiTransforms* +## *Implementing a head-motion correction algorithm for diffusion MRI in Python, using DIPY and NiTransforms* **Summary**. This tutorial walks attendees through the development of one fundamental step in the pre-processing of diffusion MRI data using a community-driven approach and relying on existing tools. diff --git a/requirements.txt b/requirements.txt index e81392a..116737c 100644 --- a/requirements.txt +++ b/requirements.txt @@ -10,5 +10,4 @@ nilearn nitransforms niworkflows numpy -requests -sphinx-exercise \ No newline at end of file +requests \ No newline at end of file