{ "cells": [ { "cell_type": "markdown", "source": [ "## Diffusion preprocessing\n", "\n", "Diffusion preprocessing typically comprises of a series of steps, which may vary depending on how the data is acquired. Some consensus has been reached for certain preprocessing steps, while others are still up for debate. The lesson will primarily focus on the preprocessing steps where consensus has been reached. Preprocessing is performed using a few well-known software packages (e.g. FSL, ANTs). For the purposes of these lessons, preprocessing steps requiring these software packages has already been performed for the dataset `ds000221` and the commands required for each step will be provided. This dataset contains single shell diffusion data with 7 b=0 s/mm^2 volumes (non-diffusion weighted) and 60 b=1000 s/mm^2 volumes. In addition, field maps (found in the `fmap` directory are acquired with opposite phase-encoding directions).\n", "\n", "To illustrate what the preprocessing step may look like, here is an example preprocessing workflow from QSIPrep (Cieslak _et al_, 2020):\n", "![preprocess](../fig/preprocessing/preprocess_steps.jpg)\n", "\n", "dMRI has some similar challenges to fMRI preprocessing, as well as some unique [ones](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3366862/).\n", "\n", "Our preprocesssing of this data will consist of following steps and will make use of sub-010006:\n", "1. Brainmasking the diffusion data\n", "2. Applying FSL `topup` to correct for susceptibility induced distortions\n", "3. FSL Eddy current distortion correction\n", "4. Registration to T1w" ], "metadata": {} }, { "cell_type": "markdown", "source": [ "### Brainmasking\n", "\n", "The first step to the preprocessing workflow is to create an appropriate brainmask from the diffusion data! Start, by first importing the necessary modules. and reading the diffusion data! We will also grab the anatomical T1w image to use later on, as well as the second inversion from the anatomical acquisition for brainmasking purposes." ], "metadata": {} }, { "cell_type": "code", "execution_count": null, "source": [ "import bids\n", "from bids.layout import BIDSLayout\n", "\n", "bids.config.set_option('extension_initial_dot', True)\n", "\n", "layout = BIDSLayout(\"../data/ds000221\", validate=False)\n", "\n", "subj='010006'\n", "\n", "# Diffusion data\n", "dwi = layout.get(subject=subj, suffix='dwi', extension='nii.gz', return_type='file')[0]\n", "\n", "# Anatomical data\n", "t1w = layout.get(subject=subj, suffix='T1w', extension='nii.gz', return_type='file')[0]" ], "outputs": [], "metadata": {} }, { "cell_type": "code", "execution_count": null, "source": [ "import numpy as np\n", "import nibabel as nib\n", "\n", "dwi = nib.load(dwi)\n", "dwi_affine = dwi.affine\n", "dwi_data = dwi.get_fdata()" ], "outputs": [], "metadata": {} }, { "cell_type": "markdown", "source": [ "DIPY's `segment.mask` module will be used to create a brainmask from this. This module contains a function `median_otsu`, which can be used to segment the brain and provide a binary brainmask! Here, a brainmask will be created using the first non-diffusion volume of the data. We will save this brainmask to be used in our later future preprocessing steps. After creating the brainmask, we will start to correct for distortions in our images." ], "metadata": {} }, { "cell_type": "code", "execution_count": null, "source": [ "import os\n", "from dipy.segment.mask import median_otsu\n", "\n", "# vol_idx is a 1D-array containing the index of the first b0\n", "dwi_brain, dwi_mask = median_otsu(dwi_data, vol_idx=[0])\n", "\n", "# Create necessary folders to save mask\n", "out_dir = '../data/ds000221/derivatives/uncorrected/sub-%s/ses-01/dwi/' % subj\n", "\n", "# Check to see if directory exists, if not create one\n", "if not os.path.exists(out_dir):\n", " os.makedirs(out_dir)\n", "\n", "img = nib.Nifti1Image(dwi_mask.astype(np.float32), dwi_affine)\n", "nib.save(img, os.path.join(out_dir, \"sub-%s_ses-01_brainmask.nii.gz\" % subj))" ], "outputs": [], "metadata": {} }, { "cell_type": "markdown", "source": [ "### FSL [`topup`](https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/topup)\n", "\n", "Diffusion images, typically acquired using spin-echo echo planar imaging (EPI), are sensitive to non-zero off-resonance fields. One source of these fields is from the susceptibilitiy distribution of the subjects head, otherwise known as susceptibility-induced off-resonance field. This field is approximately constant for all acquired diffusion images. As such, for a set of diffusion volumes, the susceptibility-induced field will be consistent throughout. This is mainly a problem due to geometric mismatches with the anatomical images (e.g. T1w), which are typically unaffected by such distortions.\n", "\n", "`topup`, part of the FSL library, estimates and attempts to correct the susceptibility-induced off-resonance field by using 2 (or more) acquisitions, where the acquisition parameters differ such that the distortion differs. Typically, this is done using two acquisitions acquired with opposite phase-encoding directions, which results in the same field creating distortions in opposing directions.\n", "\n", "Here, we will make use of the two opposite phase-encoded acquisitions found in the `fmap` directory of each subjet. These are acquired with a diffusion weighting of b = 0 s/mm^2. Alternatively, if these are not available, one can also extract and make use of the non-diffusion weighted images (assuming the data is also acquired with opposite phase encoding directions).\n", "\n", "First, we will merge the two files so that all of the volumes are in 1 file." ], "metadata": {} }, { "cell_type": "code", "execution_count": null, "source": [ "%%bash\n", "\n", "mkdir -p ../data/ds000221/derivatives/uncorrected_topup/sub-010006/ses-01/dwi/work\n", "\n", "fslmerge -t \\\n", " ../data/ds000221/derivatives/uncorrected_topup/sub-010006/ses-01/dwi/work/sub-010006_ses-01_acq-SEfmapDWI_epi.nii.gz \\\n", " ../data/ds000221/sub-010006/ses-01/fmap/sub-010006_ses-01_acq-SEfmapDWI_dir-AP_epi.nii.gz \\\n", " ../data/ds000221/sub-010006/ses-01/fmap/sub-010006_ses-01_acq-SEfmapDWI_dir-PA_epi.nii.gz" ], "outputs": [], "metadata": { "jupyter": { "outputs_hidden": true } } }, { "cell_type": "markdown", "source": [ "Another file we will need to create is a text file containing the information about how the volumes were acquired. Each line in this file will pertain to a single volume in the merged file. The first 3 values of each line refers to the acquisition direction, typically along the y-axis (or anterior-posterior). The final value is the total readout time (from center of first echo to center of final echo), which can be determined from values contained within the json sidecar. Each line will look similar to `[x y z TotalReadoutTime]`. In this case, the file, which we created, is contained within the `pedir.txt` file in the derivative directory.\n", "\n", "```\n", "0 1 0 0.04914\n", "0 1 0 0.04914\n", "0 1 0 0.04914\n", "0 -1 0 0.04914\n", "0 -1 0 0.04914\n", "0 -1 0 0.04914\n", "```\n", "\n", "With these two inputs, the next step is to make the call to `topup` to estimate the susceptibility-induced field. Within the call, a few parameters are used. Briefly:\n", " * `--imain` specifies the previously merged volume\n", " * `--datain` specifies the text file containing the information regarding the acquisition.\n", " * `--config=b02b0.cnf` makes use of a predefined config file supplied with `topup`, which contains parameters useful to registering with good b=0 s/mm^2 images.\n", " * `--out` defines the output files containing the spline coefficients for the induced field, as well as subject movement parameters" ], "metadata": {} }, { "cell_type": "code", "execution_count": null, "source": [ "%%bash\n", "\n", "topup \\\n", " --imain=../data/ds000221/derivatives/uncorrected_topup/sub-010006/ses-01/dwi/work/sub-010006_ses-01_acq-SEfmapDWI_epi.nii.gz \\\n", " --datain=../data/ds000221/derivatives/uncorrected_topup/sub-010006/ses-01/dwi/work/pedir.txt \\\n", " --config=b02b0.cnf \\\n", " --out=../data/ds000221/derivatives/uncorrected_topup/sub-010006/ses-01/dwi/work/topup" ], "outputs": [], "metadata": {} }, { "cell_type": "markdown", "source": [ "Next, we can apply the correction to the entire diffusion weighted volume by using `applytopup` Similar to `topup`, a few parameters are used. Briefly:\n", " * `--imain` specifies the input diffusion weighted volume\n", " * `--datain` again specifies the text file containing information regarding the acquisition - same file previously used\n", " * `--inindex` specifies the index (comma separated list) of the input image to be corrected\n", " * `--topup` name of field/movements (from previous topup step\n", " * `--out` basename for the corrected output image\n", " * `--method` (optional) jacobian modulation (jac) or least-squares resampling (lsr)" ], "metadata": {} }, { "cell_type": "code", "execution_count": null, "source": [ "%%bash\n", "\n", "applytopup \\\n", " --imain=../data/ds000221/sub-010006/ses-01/dwi/sub-010006_ses-01_dwi.nii.gz \\\n", " --datain=../data/ds000221/derivatives/uncorrected_topup/sub-010006/ses-01/dwi/work/pedir.txt \\\n", " --inindex=1 \\\n", " --topup=../data/ds000221/derivatives/uncorrected_topup/sub-010006/ses-01/dwi/work/topup \\\n", " --out=../data/ds000221/derivatives/uncorrected_topup/sub-010006/ses-01/dwi/dwi \\\n", " --method=jac" ], "outputs": [], "metadata": {} }, { "cell_type": "markdown", "source": [ "### FSL [`Eddy`](https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/eddy)\n", "\n", "Another source of the non-zero off resonance fields is caused by the rapid switching of diffusion weighting gradients, otherwise known as eddy current-induced off-resonance fields. Additionally, the subject is likely to move during the diffusion protocol, which may be lengthy.\n", "\n", "`eddy`, also part of the FSL library, attempts to correct for both eddy current-induced fields and subject movement by reading the gradient table and estimating the distortion volume by volume. This tool is also able to optionally detect and replace outlier slices.\n", "\n", "Here, we will demonstrate the application of `eddy` following the `topup` correction step, by making use of both the uncorrected diffusion data, as well as distortion corrections from the previous step. Additionally, a text file, which maps each of the volumes to one of the corresponding acquisition directions from the `pedir.txt` file will have to be created. Finally, similar to `topup`, there are also a number of input parameters which have to be specified:\n", "\n", " * `--imain` specifies the undistorted diffusion weighted volume\n", " * `--mask` specifies the brainmask for the undistorted diffusion weighted volume\n", " * `--acqp` specifies the the text file containing information regarding the acquisition that was previously used in `topup`\n", " * `--index` is the text file which maps each diffusion volume to the corresponding acquisition direction\n", " * `--bvecs` specifies the bvec file to the undistorted dwi\n", " * `--bvals` similarily specifies the bval file to the undistorted dwi\n", " * `--topup` specifies the directory and distortion correction files previously estimated by `topup`\n", " * `--out` specifies the prefix of the output files following eddy correction\n", " * `--repol` is a flag, which specifies replacement of outliers" ], "metadata": {} }, { "cell_type": "code", "execution_count": null, "source": [ "%%bash\n", "\n", "mkdir -p ../data/ds000221/derivatives/uncorrected_topup_eddy/sub-010006/ses-01/dwi/work\n", "\n", "# Create an index file mapping the 67 volumes in 4D dwi volume to the pedir.txt file\n", "indx=\"\"\n", "for i in `seq 1 67`; do\n", " indx=\"$indx 1\"\n", "done\n", "echo $indx > ../data/ds000221/derivatives/uncorrected_topup_eddy/sub-010006/ses-01/dwi/work/index.txt\n", "\n", "eddy_openmp \\\n", " --imain=../data/ds000221/sub-010006/ses-01/dwi/sub-010006_ses-01_dwi.nii.gz \\\n", " --mask=../data/ds000221/derivatives/uncorrected/sub-010006/ses-01/dwi/sub-010006_ses-01_brainmask.nii.gz \\\n", " --acqp=../data/ds000221/derivatives/uncorrected_topup/sub-010006/ses-01/dwi/work/pedir.txt \\\n", " --index=../data/ds000221/derivatives/uncorrected_topup_eddy/sub-010006/ses-01/dwi/work/index.txt \\\n", " --bvecs=../data/ds000221/sub-010006/ses-01/dwi/sub-010006_ses-01_dwi.bvec \\\n", " --bvals=../data/ds000221/sub-010006/ses-01/dwi/sub-010006_ses-01_dwi.bval \\\n", " --topup=../data/ds000221/derivatives/uncorrected_topup/sub-010006/ses-01/dwi/work/topup \\\n", " --out=../data/ds000221/derivatives/uncorrected_topup_eddy/sub-010006/ses-01/dwi/dwi \\\n", " --repol" ], "outputs": [], "metadata": {} }, { "cell_type": "markdown", "source": [ "### Registration with T1w\n", "\n", "The final step to our diffusion processing is registration to an anatomical image (eg. T1-weighted). This is important because the diffusion data, typically acquired using echo planar imaging or EPI, enables faster acquisitions at the cost of lower resolution and introduction of distortions (as seen above). Registration with the anatomical image not only helps to correct for some distortions, it also provides us with a higher resolution, anatomical reference.\n", "\n", "First, we will create a brainmask of the anatomical image using the second inversion. To do this, we will use FSL `bet` twice. The first call to `bet` will create a general skullstripped brain. Upon inspection, we can note that there is still some residual areas of the image which were included in the first pass. Calling `bet` a second time, we get a better outline of the brain and brainmask, which we can use for further processing." ], "metadata": {} }, { "cell_type": "code", "execution_count": null, "source": [ "%%bash\n", "\n", "mkdir -p ../data/ds000221/derivatives/uncorrected/sub-010006/ses-01/anat\n", "\n", "bet \\\n", " ../data/ds000221/sub-010006/ses-01/anat/sub-010006_ses-01_inv-2_mp2rage.nii.gz \\\n", " ../data/ds000221/derivatives/uncorrected/sub-010006/ses-01/anat/sub-010006_ses-01_space-T1w_broadbrain \\\n", " -f 0.6\n", "\n", "bet \\\n", " ../data/ds000221/derivatives/uncorrected/sub-010006/ses-01/anat/sub-010006_ses-01_space-T1w_broadbrain \\\n", " ../data/ds000221/derivatives/uncorrected/sub-010006/ses-01/anat/sub-010006_ses-01_space-T1w_brain \\\n", " -f 0.4 -m\n", "\n", "mv \\\n", " ../data/ds000221/derivatives/uncorrected/sub-010006/ses-01/anat/sub-010006_ses-01_space-T1w_brain_mask.nii.gz \\\n", " ../data/ds000221/derivatives/uncorrected/sub-010006/ses-01/anat/sub-010006_ses-01_space-T1w_brainmask.nii.gz" ], "outputs": [], "metadata": { "jupyter": { "outputs_hidden": true } } }, { "cell_type": "markdown", "source": [ "Note, we use `bet` here, as well as the second inversion of the anatomical image, as it provides us with a better brainmask. The `bet` command above is called to output only the binary mask and the fractional intensity threshold is also increased slightly (to 0.6) provide a smaller outline of the brain initially, and then decreased (to 0.4) to provide a larger outline. The flag `-m` indicates to the tool to create a brainmask in addition to outputting the extracted brain volume. Both the mask and brain volume will be used in our registration step.\n", "\n", "Before we get to the registration, we will also update our DWI brainmask by performing a brain extraction using `dipy` on the eddy corrected image. Note that the output of `eddy` is not in BIDS format so we will include the path to the diffusion data manually. We will save both the brainmask and the extracted brain volume. Additionally, we will save a separate volume of only the first b0 to use for the registration." ], "metadata": {} }, { "cell_type": "code", "execution_count": null, "source": [ "from dipy.segment.mask import median_otsu\n", "\n", "# Path of FSL eddy-corrected dwi\n", "dwi = \"../data/ds000221/derivatives/uncorrected_topup_eddy/sub-010006/ses-01/dwi/dwi.nii.gz\"\n", "\n", "# Load eddy-corrected diffusion data\n", "dwi = nib.load(dwi)\n", "dwi_affine = dwi.affine\n", "dwi_data = dwi.get_fdata()\n", "\n", "dwi_brain, dwi_mask = median_otsu(dwi_data, vol_idx=[0])\n", "dwi_b0 = dwi_brain[:,:,:,0]\n", "\n", "# Output directory\n", "out_dir=\"../data/ds000221/derivatives/uncorrected_topup_eddy/sub-010006/ses-01/dwi\"\n", "\n", "# Save diffusion mask\n", "img = nib.Nifti1Image(dwi_mask.astype(np.float32), dwi_affine)\n", "nib.save(img, os.path.join(out_dir, \"sub-010006_ses-01_dwi_proc-eddy_brainmask.nii.gz\"))\n", "\n", "# Save 4D diffusion volume\n", "img = nib.Nifti1Image(dwi_brain, dwi_affine)\n", "nib.save(img, os.path.join(out_dir, \"sub-010006_ses-01_dwi_proc-eddy_brain.nii.gz\"))\n", "\n", "# Save b0 volume\n", "img = nib.Nifti1Image(dwi_b0, dwi_affine)\n", "nib.save(img, os.path.join(out_dir, \"sub-010006_ses-01_dwi_proc-eddy_b0.nii.gz\"))" ], "outputs": [], "metadata": {} }, { "cell_type": "markdown", "source": [ "To perform the registration between the diffusion volumes and T1w, we will make use of ANTs, specifically the `antsRegistrationSyNQuick.sh` script and `antsApplyTransform`. We will begin by registering the diffusion b=0 s/mm^2 volume to get the appropriate transforms to align the two images. We will then apply the inverse transformation to the T1w volume such that it is aligned to the diffusion volume\n", "\n", "Here, we will constrain `antsRegistrationSyNQuick.sh` to perform a rigid and affine transformation (we will explain why in the final step). There are a few parameters that must be set:\n", "\n", "* `-d` - Image dimension (2/3D)\n", "* `-t` - Transformation type (`a` performs only rigid + affine transformation)\n", "* `-f` - Fixed image (anatomical T1w)\n", "* `-m` - Moving image (DWI b=0 s/mm^2)\n", "* `-o` - Output prefix (prefix to be appended to output files)" ], "metadata": {} }, { "cell_type": "code", "execution_count": null, "source": [ "%%bash\n", "\n", "mkdir -p ../data/ds000221/derivatives/uncorrected_topup_eddy_regT1/sub-010006/ses-01/transforms\n", "\n", "# Perform registration between b0 and T1w\n", "antsRegistrationSyNQuick.sh \\\n", " -d 3 -t a \\\n", " -f ../data/ds000221/derivatives/uncorrected/sub-010006/ses-01/anat/sub-010006_ses-01_space-T1w_brain.nii.gz \\\n", " -m ../data/ds000221/derivatives/uncorrected_topup_eddy/sub-010006/ses-01/dwi/sub-010006_ses-01_dwi_proc-eddy_b0.nii.gz \\\n", " -o ../data/ds000221/derivatives/uncorrected_topup_eddy_regT1/sub-010006/ses-01/transform/dwi_to_t1_" ], "outputs": [], "metadata": {} }, { "cell_type": "markdown", "source": [ "The transformation file should be created which we will use to apply the inverse transform with `antsApplyTransform` to the T1w volume. Similar to the previous command, there are few parameters that will need to be set:\n", "\n", "* `-d` - Image dimension (2/3/4D)\n", "* `-i` - Input volume to be transformed (T1w)\n", "* `-r` - Reference volume (b0 of DWI volume)\n", "* `-t` - Transformation file (can be called more than once)\n", "* `-o` - Output volume in the transformed space.\n", "\n", "Note that if more than 1 transformation file is provided, the order in which the transforms are applied to the volume is in reverse order of how it is inputted (e.g. last transform gets applied first)." ], "metadata": {} }, { "cell_type": "code", "execution_count": null, "source": [ "%%bash\n", "\n", "# Apply transform to 4D DWI volume\n", "antsApplyTransforms \\\n", " -d 3 \\\n", " -i ../data/ds000221/derivatives/uncorrected/sub-010006/ses-01/anat/sub-010006_ses-01_space-T1w_brain.nii.gz \\\n", " -r ../data/ds000221/derivatives/uncorrected_topup_eddy/sub-010006/ses-01/dwi/sub-010006_ses-01_dwi_proc-eddy_b0.nii.gz \\\n", " -t [../data/ds000221/derivatives/uncorrected_topup_eddy_regT1/sub-010006/ses-01/transform/dwi_to_t1_0GenericAffine.mat,1] \\\n", " -o ../data/ds000221/derivatives/uncorrected_topup_eddy_regT1/sub-010006/ses-01/anat/sub-010006_ses-01_space-dwi_T1w_brain.nii.gz" ], "outputs": [], "metadata": {} }, { "cell_type": "markdown", "source": [ "![](../fig/preprocessing/transformed_volumes.png)" ], "metadata": {} }, { "cell_type": "markdown", "source": [ "Following the transformation of the T1w volume, we can see that anatomical and diffusion weighted volumes are now aligned. It should be highlighted that as part of the transformation step, the T1w volume is resampled based on the voxel size of the ference volume (e.g. DWI).\n", "\n", "### Preprocessing notes:\n", "\n", "1. In this lesson, the T1w volume is registered to the DWI volume. This method minimizes the manipulation of the diffusion data. It is also possible to register the DWI volume to the T1w volume and would require the associated diffusion gradient vectors (bvec) to also be similarly rotated. If this step is not performed, one would have incorrect diffusion gradient directions relative to the registered DWI volumes. This also highlights a reason behind not performing a non-linear transformation for registration, as each individual diffusion gradient direction would also have to be subsequently warped. Rotation of the diffusion gradient vectors can be done by applying the affine transformation to each row of the file. Luckily, there are existing scripts that can do this. One such Python script was created by Michael Paquette: [`rot_bvecs_ants.py`](https://gist.github.com/mpaquette/5d59ad195778f9d984c5def42f53de6e).\n", "\n", "2. We have only demonstrated the preprocessing steps where there is general consensus on how DWI data should be processed. There are also additional steps with certain caveats, which include denoising, unringing (to remove/minimize effects of Gibbs ringing artifacts), and gradient non-linearity correction (to unwarp distortions caused by gradient-field inhomogeneities using a vendor acquired gradient coefficient file.\n", "\n", "3. Depending on how the data is acquired, certain steps may not be possible. For example, if the data is not acquired in two directions, `topup` may not be possible (in this situation, distortion correction may be better handled by registering with a T1w anatomical image directly.\n", "\n", "4. There are also a number of tools available for preprocessing. In this lesson, we demonstrate some of the more commonly used tools alongside `dipy`." ], "metadata": {} }, { "cell_type": "markdown", "source": [ "### References\n", "\n", ".. [Cieslak2020] M. Cieslak, PA. Cook, X. He, F-C. Yeh, T. Dhollander, _et al_, \"QSIPrep: An integrative platform for preprocessing and reconstructing diffusion MRI\", https://doi.org/10.1101/2020.09.04.282269" ], "metadata": {} } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.4" } }, "nbformat": 4, "nbformat_minor": 4 }