There are some important considerations for how your diffusion data will be preprocessed.
You must choose how to combine dwi scans within a session (Merging multiple scans from a session)
How to correct for susceptibility distortion.
How to perform motion correction
Building a pipeline¶
Merging multiple scans from a session¶
For q-space imaging sequences it is common to have multiple separate scans to acquire the entire sampling scheme. These scans get aligned and merged into a single DWI series before reconstruction. It is also common to collect a DWI scan (or scans) in the reverse phase encoding direction to use for susceptibility distortion correction (SDC).
This creates a number of possible scenarios for preprocessing your DWIs. These
scenarios can be controlled by the
combine_all_dwis argument. If your study
has multiple sessions, DWI scans will never be combined across sessions.
Merging only occurs within a session.
False (not present in the commandline call), each dwi
scan in the
dwi directories will be processed independently. You will have one
preprocessed output per each DWI file in your input.
combine_all_dwis is set to
True, two possibilities arise. If all DWIs in a session
are in the same PE direction, they will be merged into a single series. If there are
two PE directions detected in the DWI scans and
'fieldmaps' is not in
images are combined according to their PE direction, and their b0 reference images are used to
perform SDC. Further complicating this is the FSL workflow, which combines distortion correction
with eddy/motion correction and will merge scans with different PE directions.
If you have some scans you want to combine and others you want to preprocess separately, consider creating fake sessions in your BIDS directory.
Susceptibility correction methods¶
The are three kinds of SDC available in qsiprep:
Phase Encoding POLARity (PEPOLAR) techniques (also called blip-up/blip-down): This is the implementation from sdcflows, using 3dQwarp to correct a DWI series using a fieldmap in the fmaps directory [Jezzard1995]. The reverse phase encoding direction scan can come from the fieldmaps directory or the dwi directory.
1a. If using
TOPUPis used for this correction.
Phase-difference B0 estimation: Use a B0map sequence that includes at lease one magnitude image and two phase images or a phasediff image.
Fieldmap-less estimation (experimental): The SyN-based susceptibility distortion correction implemented in FMRIPREP. To use this method, include argument –use-syn-sdc when calling qsiprep. Briefly, this method estimates a SDC warp using ANTS SyN based on an average fieldmap in MNI space. For details on this method, see fmriprep’s documentation
qsiprep determines if a fieldmap should be used based on the
fields in the JSON sidecars in the
QSIPrep can be configured to produce a very similar pipeline to the HCP dMRI pipelines.
HCP and HCP-Lifespan scans acquire complete multi-shell sequences in opposing phase
encoding directions, making them a special case where Phase Encoding POLARity (PEPOLAR) techniques are used
and the corrected images from both PE directions are averaged at the end. To produce
qsiprep that is directly comparable to the HCP dMRI pipeline you
will want to include:
--distortion-group-merge average \ --combine-all-dwis \
If you want to disable the image pair averaging and get a result with twice as
many images, you can substitute
Outputs of qsiprep¶
qsiprep generates three broad classes of outcomes:
Visual QA (quality assessment) reports: one HTML per subject, depicting images that provide a sanity check for each step of the pipeline.
Pre-processed imaging data such as anatomical segmentations, realigned and resampled diffusion weighted images and the corresponding corrected gradient files in FSL and MRTrix format.
Additional data for subsequent analysis, for instance the transformations between different spaces or the estimated head motion and model fit quality calculated during model-based head motion correction.
Quantitative QA: A single-line csv file per subject summarizing subject motion, coregistration quality and image quality.
qsiprep outputs summary reports, written to
These reports provide a quick way to make visual inspection of the results easy. One useful
graphic is the animation of the q-space sampling scheme before and after the pipeline. Here is
a sampling scheme from a DSI scan:
Preprocessed data (qsiprep derivatives)¶
There are additional files, called “Derivatives”, written to
Derivatives related to T1w files are nearly identical to those produced by
can be found in the
*T1w_brainmask.nii.gzBrain mask derived using ANTs’
*T1w_dtissue.nii.gzTissue class map derived using FAST.
*T1w_preproc.nii.gzBias field corrected T1w file, using ANTS’ N4BiasFieldCorrection
_brainmaskabove, but in MNI space.
*T1w_space-MNI152NLin2009cAsym_class-WM_probtissue.nii.gzProbability tissue maps, transformed into MNI space
_dtissueabove, but in MNI space
_preprocabove, but in MNI space
*T1w_space-MNI152NLin2009cAsym_target-T1w_warp.h5Composite (warp and affine) transform to map from MNI to T1 space
*T1w_target-MNI152NLin2009cAsym_warp.h5Composite (warp and affine) transform to transform T1w into MNI space
Derivatives related to diffusion images are in the
*_confounds.tsvA tab-separated value file with one column per calculated confound and one row per timepoint/volume
Volumetric output spaces include
T1w (default) and
*dwiref.nii.gzThe b0 template
*desc-brain_mask.nii.gzThe generous brain mask that should be reduced probably
*desc-preproc_dwi.nii.gzResampled DWI series including all b0 images.
*desc-preproc_dwi.bvecFSL-style bvals and bvecs files. These will be incorrectly interpreted by MRTrix, but will work with DSI Studio. Use the
.bfile for MRTrix.
desc-preproc_dwi.bThe gradient table to import data into MRTrix. This and the
_dwi.nii.gzcan be converted directly to a
.miffile using the
mrconvert -grad _dwi.bcommand.
*b0series.nii.gzThe b0 images from the series in a 4d image. Useful to see how much the images are impacted by Eddy currents.
*bvecs.nii.gzEach voxel contains a gradient table that has been adjusted for local rotations introduced by spatial warping.
*cnr.nii.gzEach voxel contains a contrast-to-noise model defined as the variance of the signal model divided by the variance of the error of the signal model.
See implementation on
For each DWI processed by qsiprep, a
file will be generated.
These are TSV tables, which look like the example below:
framewise_displacement trans_x trans_y trans_z rot_x rot_y rot_z hmc_r2 hmc_xcorr original_file grad_x grad_y grad_z bval n/a -0.705 -0.002 0.133 0.119 0.350 0.711 0.941 0.943 sub-abcd_dwi.nii.gz 0.000 0.000 0.000 0.000 16.343 -0.711 -0.075 0.220 0.067 0.405 0.495 0.945 0.946 sub-abcd_dwi.nii.gz 0.000 0.000 0.000 0.000 35.173 -0.672 -0.415 0.725 0.004 0.468 1.055 0.756 0.766 sub-abcd_dwi.nii.gz -0.356 0.656 0.665 3000.000 45.131 0.021 -0.498 1.046 0.403 0.331 1.400 0.771 0.778 sub-abcd_dwi.nii.gz -0.935 0.272 0.229 3000.000 37.506 -0.184 0.117 0.723 0.305 0.138 0.964 0.895 0.896 sub-abcd_dwi.nii.gz -0.187 -0.957 -0.223 2000.000 16.388 -0.447 0.020 0.847 0.217 0.129 0.743 0.792 0.800 sub-abcd_dwi.nii.gz -0.111 -0.119 0.987 3000.000
The motion parameters come from the model-based head motion estimation workflow. The
hmc_xcorr are whole-brain r^2 values and cross correlation scores (using the ANTs definition)
between the model-generated target image and the motion-corrected empirical image. The final
columns are not really confounds, but book-keeping information that reminds us which 4d DWI series
the image originally came from and what gradient direction (
and gradient strength
bval the image came from. This can be useful for tracking down
mostly-corrupted scans and can indicate if the head motion model isn’t working on specific
gradient strengths or directions.
Quality Control data¶
A single-line csv file (
desc-ImageQC_dwi.csv) is created for each output image. This file is
particularly useful for comparing the relative quality across subjects before deciding who to
include in a group analysis. The columns in this file come from DSI Studio’s QC calculation and is
described in [Yeh2019]. Columns prefixed by
raw_ reflect QC measurements from the data before
preprocessing. Columns prefixed by
mni_ contain QC metrics calculated on the
preprocessed data. Motion parameter summaries are also provided, such as the mean and max of
framewise displacement (
max_fd). The max and mean absolute values for translation
and rotation are
max_rotation and the maxima of their derivatives are
max_rel_rotation. Finally, the difference in spatial overlap
between the anatomical mask and the anatomical brain mask and the DWI brain mask is calculated
using the Dice distance in
Confounds and “carpet”-plot on the visual reports¶
fMRI has been using a “carpet” visualization of the BOLD time-series (see [Power2016]), but this type of plot does not make sense for DWI data. Instead, we plot the cross-correlation value between each raw slice and the HMC model signal resampled into that slice. This plot is included for each run within the corresponding visual report. Examples of these plots follow:
Preprocessing pipeline details¶
qsiprep adapts its pipeline depending on what data and metadata are
available and are used as the input.
The anatomical sub-workflow begins by constructing an average image by conforming all found T1w images to LPS orientation and a common voxel size, and, in the case of multiple images, averages them into a single reference template (see Longitudinal processing).
Brain extraction, brain tissue segmentation and spatial normalization¶
Then, the T1w image/average is skull-stripped using ANTs’
which is an atlas-based brain extraction workflow.
Once the brain mask is computed, FSL
fast is utilized for brain tissue segmentation.
Finally, spatial normalization to MNI-space is performed using ANTs’
in a multiscale, mutual-information based, nonlinear registration scheme.
In particular, spatial normalization is done using the ICBM 2009c Nonlinear
Asymmetric template (1×1×1mm) [Fonov2011].
When processing images from patients with focal brain lesions (e.g. stroke, tumor
resection), it is possible to provide a lesion mask to be used during spatial
normalization to MNI-space [Brett2001].
ANTs will use this mask to minimize warping of healthy tissue into damaged
areas (or vice-versa).
Lesion masks should be binary NIfTI images (damaged areas = 1, everywhere else = 0)
in the same space and resolution as the T1 image, and follow the naming convention specified in
BIDS Extension Proposal 3: Common Derivatives
This file should be placed in the
sub-*/anat directory of the BIDS dataset
to be run through
In the case of multiple T1w images (across sessions and/or runs), T1w images are
merged into a single template image using FreeSurfer’s
This template may be unbiased, or equidistant from all source images, or
aligned to the first image (determined lexicographically by session label).
For two images, the additional cost of estimating an unbiased template is
trivial and is the default behavior, but, for greater than two images, the cost
can be a slowdown of an order of magnitude.
Therefore, in the case of three or more images,
templates aligned to the first image, unless passed the
flag, which forces the estimation of an unbiased template.
The preprocessed T1w image defines the
In the case of multiple T1w images, this space may not be precisely aligned
with any of the original images.
Reconstructed surfaces and functional datasets will be registered to the
T1w space, and not to the input images.
Preprocessing of DWI files is split into multiple sub-workflows described below.
Head-motion estimation (SHORELine)¶
A long-standing issue for q-space imaging techniques, particularly DSI, has
been the lack of motion correction methods. DTI and multi-shell HARDI have
eddy in FSL, but DSI has relied on aligning the
interleaved b0 images and applying the transforms to nearby non-b0 images.
qsiprep introduces a method for head motion correction that iteratively
creates target images based on
First, all b0 images are aligned to a midpoint b0 image (or the first b0 image
hmc_align_to="first") and each non-b0 image is transformed along with
its nearest b0 image.
Then, for each non-b0 image, a
model is fit to all the other images with that image left out. The model is then
used to generate a target signal image for the gradient direction and magnitude
(i.e. q-space coordinate) of the left-out image. The left-out image is registered
to the generated target
signal image and its vector is rotated accordingly. A new model is fit on the
transformed images and their rotated vectors. The leave-one-out procedure is
then repeated on this updated DWI and gradient set.
"none" is specified as the hmc_model, then only the b0 images are used
and the non-b0 images are transformed based on their nearest b0 image. This
is probably not a great idea.
Susceptibility distortion correction is run as part of this pipeline to be
consistent with the
Ultimately a list of 6 (or 12)-parameters per time-step is written and fed to the confounds workflow. These are used to estimate framewise displacement. Additionally, measures of model fits are saved for each slice for display in a carpet plot-like thing.
Head-motion estimation (TOPUP/eddy)¶
DTI and multi-shell HARDI can be passed to
head motion correction, susceptibility distortion correction and eddy current
qsiprep will use the BIDS-specified fieldmaps to configure and
TOPUP before passing the fieldmap to
eddy if fieldmaps are available.
DWI reference image estimation¶
This workflow estimates a reference image for a DWI series. This procedure is different from the DWI reference image workflow in the sense that true brain masking isn’t usually done until later in the pipeline for DWIs. It performs a generous automasking and uses Dipy’s histogram equalization on the b0 template generated during motion correction.
Susceptibility Distortion Correction (SDC)¶
The PEPOLAR and SyN-SDC workflows from FMRIPREP are copied here. They operate on the output of reference estimation, after head motion correction. For a complete list of possibilities here, see Merging multiple scans from a session.
Pre-processed DWIs in a different space¶
A DWI series is resampled to an output space. The
specified on the commandline call. All transformations, including head motion
correction, susceptibility distortion correction, coregistration and (optionally)
normalization to the template is performed in a single shot using a Lanczos kernel.
There are two ways that the gradient vectors can be saved. This workflow always
produces a FSL-style bval/bvec pair for the image and a MRTrix .b gradient table
with the rotations from the linear transforms applied. You can also write out
local_bvecs file that contains a 3d vector that has been rotated to account
for nonlinear transforms in each voxel. I’m not aware of any software that can
use these yet, but it’s an interesting idea.