fMRIPrep adapts its pipeline depending on what data and metadata are
available and are used as the input.
For example, slice timing correction will be
performed only if the SliceTiming metadata field is found for the input
dataset.
A (very) high-level view of the simplest pipeline (for a single-band dataset with only one task, single-run, with no slice-timing information nor fieldmap acquisitions) is presented below:
.. workflow::
:graph2use: orig
:simple_form: yes
from fmriprep.workflows.tests import mock_config
from fmriprep.workflows.base import init_single_subject_wf
with mock_config():
wf = init_single_subject_wf('01')
Note
Each node in this workflow is either a processing node or a sub-workflow. Several conventions appear in this workflow that will be apparent throughout fMRIPrep.
inputnodes are special nodes that provide the runtime-generated inputs to a workflow. These are like function "arguments". There are correspondingoutputnodes in most other workflows, which are like function return values.- Workflows end with
_wf, and are generated by a function of the forminit_{workflow}_wf. For example,anat_preproc_wfis a sub-workflow that is generated by the :func:`~smriprep.workflows.anatomical.init_anat_preproc_wf` (see below). Because each task and run of functional data is processed separately, :func:`~fmriprep.workflows.bold.base.init_func_preproc_wf` names the resulting workflows using input parameters, resulting infunc_preproc_task_{task}_run_{run}_wf. - Datasinks begin with
ds_, and save files to the output directory. This is in contrast to most nodes, which save their outputs to the working directory.ds_report_nodes indicate that the node is saving text and figures for generating reports, rather than processed data. - When a name appears in parentheses, such as
(reports)inabout (reports)it is the module where the interface is defined. In this case,aboutis an :class:`~fmriprep.interfaces.reports.AboutSummary`, found in the :mod:`fmriprep.interfaces.reports` module.
The anatomical sub-workflow begins by constructing an average image by conforming all found T1w images to RAS orientation and a common voxel size, and, in the case of multiple images, averages them into a single reference template (see Longitudinal processing).
.. workflow::
:graph2use: orig
:simple_form: yes
from niworkflows.utils.spaces import Reference, SpatialReferences
from smriprep.workflows.anatomical import init_anat_preproc_wf
wf = init_anat_preproc_wf(
bids_root='.',
freesurfer=True,
hires=True,
longitudinal=False,
omp_nthreads=1,
output_dir='.',
skull_strip_mode='force',
skull_strip_template=Reference('MNI152NLin2009cAsym'),
spaces=SpatialReferences([
('MNI152Lin', {}),
('fsaverage', {'density': '10k'}),
('T1w', {}),
('fsnative', {})
]),
skull_strip_fixed_seed=False,
t1w=['sub-01/anat/sub-01_T1w.nii.gz'],
t2w=[],
)
Important
Occasionally, openly shared datasets may contain preprocessed anatomical images
as if they are unprocessed.
In the case of brain-extracted (skull-stripped) T1w images, attempting to perform
brain extraction again will often have poor results and may cause fMRIPrep to crash.
fMRIPrep can attempt to detect these cases using a heuristic to check if the
T1w image is already masked.
This must be explicitly requested with ---skull-strip-t1w auto.
If this heuristic fails, and you know your images are skull-stripped, you can skip brain
extraction with --skull-strip-t1w skip.
Likewise, if you know your images are not skull-stripped and the heuristic incorrectly
determines that they are, you can force skull stripping with --skull-strip-t1w force,
which is the current default behavior.
See also sMRIPrep's :py:func:`~smriprep.workflows.anatomical.init_anat_preproc_wf`.
Then, the T1w reference is skull-stripped using a Nipype implementation of
the antsBrainExtraction.sh tool (ANTs), which is an atlas-based
brain extraction workflow:
.. workflow::
:graph2use: orig
:simple_form: yes
from niworkflows.anat.ants import init_brain_extraction_wf
wf = init_brain_extraction_wf()
Once the brain mask is computed, FSL fast is utilized for brain tissue segmentation.
fMRIPrep includes a single figure overlaying the brain mask (red), and tissue boundaries (blue = gray/white; magenta = tissue/CSF):
Brain extraction and segmentation report
Finally, spatial normalization to standard spaces is performed using ANTs' antsRegistration
in a multiscale, mutual-information based, nonlinear registration scheme.
See :ref:`output-spaces` for information about how standard and nonstandard spaces can
be set to resample the preprocessed data onto the final output spaces.
Animation showing spatial normalization of T1w onto the MNI152NLin2009cAsym template.
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 standard 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 T1w image, and use the _roi suffix,
for example, sub-001_label-lesion_roi.nii.gz.
This file should be placed in the sub-*/anat directory of the BIDS dataset
to be run through fMRIPrep.
Because lesion masks are not currently part of the BIDS specification, it is also necessary to
include a .bidsignore file in the root of your dataset directory. This will prevent
bids-validator from complaining
that your dataset is not valid BIDS, which prevents fMRIPrep from running.
Your .bidsignore file should include the following line:
*lesion_roi.nii.gz
Note
The lesion masking instructions in this section predate the release of BIDS Derivatives. As of BIDS 1.4.0, the recommended naming convention is:
manual_masks/
└─ sub-001/
└─ anat/
├─ sub-001_desc-tumor_mask.nii.gz
└─ sub-001_desc-tumor_mask.json
In an upcoming version of fMRIPrep, we will search for lesion masks as pre-computed
derivatives. Until this is supported, we will continue to look for the _roi suffix.
In the case of multiple T1w images (across sessions and/or runs), T1w images are
merged into a single template image using FreeSurfer's `mri_robust_template`_.
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,
but aligning three or more images is too expensive to justify being the default behavior.
For consistency, in the case of multiple images, fMRIPrep constructs
templates aligned to the first image, unless passed the --longitudinal
flag, which forces the estimation of an unbiased template.
Note
The preprocessed T1w image defines the T1w space.
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.
fMRIPrep uses FreeSurfer_ to reconstruct surfaces from T1w/T2w
structural images.
If enabled, several steps in the fMRIPrep pipeline are added or replaced.
All surface preprocessing may be disabled with the --fs-no-reconall flag.
Note
Surface processing will be skipped if the outputs already exist.
In order to bypass reconstruction in fMRIPrep, place existing reconstructed
subjects in <output dir>/sourcedata/freesurfer prior to the run, or specify
an external subjects directory with the --fs-subjects-dir flag.
fMRIPrep will perform any missing recon-all steps, but will not perform
any steps whose outputs already exist.
If FreeSurfer reconstruction is performed, the reconstructed subject is placed in
<output dir>/sourcedata/freesurfer/sub-<subject_label>/ (see :ref:`fsderivs`).
Surface reconstruction is performed in three phases.
The first phase initializes the subject with T1w and T2w (if available)
structural images and performs basic reconstruction (autorecon1) with the
exception of skull-stripping.
Skull-stripping is skipped since the brain mask :ref:`calculated previously
<t1preproc_steps>` is injected into the appropriate location for FreeSurfer.
For example, a subject with only one session with T1w and T2w images
would be processed by the following command:
$ recon-all -sd <output dir>/freesurfer -subjid sub-<subject_label> \
-i <bids-root>/sub-<subject_label>/anat/sub-<subject_label>_T1w.nii.gz \
-T2 <bids-root>/sub-<subject_label>/anat/sub-<subject_label>_T2w.nii.gz \
-autorecon1 \
-noskullstrip
The second phase imports the brainmask calculated in the Preprocessing of structural MRI sub-workflow. The final phase resumes reconstruction, using the T2w image to assist in finding the pial surface, if available. See :py:func:`~smriprep.workflows.surfaces.init_autorecon_resume_wf` for details.
Reconstructed white and pial surfaces are included in the report.
Surface reconstruction (FreeSurfer)
If T1w voxel sizes are less than 1mm in all dimensions (rounding to nearest
.1mm), `submillimeter reconstruction`_ is used, unless disabled with
--no-submm-recon.
If T2w or FLAIR images are available, and you do not want them included in
FreeSurfer reconstruction, use --ignore t2w or --ignore flair,
respectively.
lh.midthickness and rh.midthickness surfaces are created in the subject
surf/ directory, corresponding to the surface half-way between the gray/white
boundary and the pial surface.
The smoothwm, midthickness, pial and inflated surfaces are also
converted to GIFTI_ format and adjusted to be compatible with multiple software
packages, including FreeSurfer and the `Connectome Workbench`_.
Note
GIFTI surface outputs are aligned to the FreeSurfer T1.mgz image, which
may differ from the T1w space in some cases, to maintain compatibility
with the FreeSurfer directory.
Any measures sampled to the surface take into account any difference in
these images.
.. workflow::
:graph2use: orig
:simple_form: yes
from smriprep.workflows.surfaces import init_surface_recon_wf
wf = init_surface_recon_wf(omp_nthreads=1,
hires=True)
See also sMRIPrep's :py:func:`~smriprep.workflows.surfaces.init_surface_recon_wf`
Typically, the original brain mask calculated with antsBrainExtraction.sh
will contain some innaccuracies including small amounts of MR signal from
outside the brain.
Based on the tissue segmentation of FreeSurfer (located in mri/aseg.mgz)
and only when the :ref:`Surface Processing <workflows_surface>` step has been
executed, fMRIPrep replaces the brain mask with a refined one that derives
from the aseg.mgz file as described in
:py:class:`~niworkflows.interfaces.freesurfer.RefineBrainMask`.
:py:func:`~fmriprep.workflows.bold.base.init_func_preproc_wf`
.. workflow::
:graph2use: orig
:simple_form: yes
from fmriprep.workflows.tests import mock_config
from fmriprep import config
from fmriprep.workflows.bold.base import init_func_preproc_wf
with mock_config():
bold_file = config.execution.bids_dir / 'sub-01' / 'func' \
/ 'sub-01_task-mixedgamblestask_run-01_bold.nii.gz'
wf = init_func_preproc_wf(str(bold_file))
Preprocessing of :abbr:`BOLD (blood-oxygen level-dependent)` files is split into multiple sub-workflows described below.
:py:func:`~niworkflows.func.util.init_bold_reference_wf`
.. workflow::
:graph2use: orig
:simple_form: yes
from niworkflows.func.util import init_bold_reference_wf
wf = init_bold_reference_wf(omp_nthreads=1)
This workflow estimates a reference image for a :abbr:`BOLD (blood-oxygen level-dependent)` series. If a single-band reference ("sbref") image associated with the BOLD series is available, then it is used directly. If not, a reference image is estimated from the BOLD series as follows: When T1-saturation effects ("dummy scans" or non-steady state volumes) are detected, they are averaged and used as reference due to their superior tissue contrast. Otherwise, a median of motion corrected subset of volumes is used.
The reference image is then used to calculate a brain mask for the :abbr:`BOLD (blood-oxygen level-dependent)` signal using NiWorkflows' :py:func:`~niworkflows.func.util.init_enhance_and_skullstrip_bold_wf`. Further, the reference is fed to the :ref:`head-motion estimation workflow <bold_hmc>` and the :ref:`registration workflow to map BOLD series into the T1w image of the same subject <bold_reg>`.
The red contour shows the brain mask estimated for a BOLD reference volume. The blue and magenta contours show the tCompCor and aCompCor masks, respectively. (See :ref:`bold_confounds`, below.)
:py:func:`~fmriprep.workflows.bold.hmc.init_bold_hmc_wf`
.. workflow::
:graph2use: colored
:simple_form: yes
from fmriprep.workflows.bold import init_bold_hmc_wf
wf = init_bold_hmc_wf(
mem_gb=1,
omp_nthreads=1)
Using the previously :ref:`estimated reference scan <bold_ref>`,
FSL mcflirt is used to estimate head-motion.
As a result, one rigid-body transform with respect to
the reference image is written for each :abbr:`BOLD (blood-oxygen level-dependent)`
time-step.
Additionally, a list of 6-parameters (three rotations,
three translations) per time-step is written and fed to the
:ref:`confounds workflow <bold_confounds>`.
For a more accurate estimation of head-motion, we calculate its parameters
before any time-domain filtering (i.e., :ref:`slice-timing correction <bold_stc>`),
as recommended in [Power2017].
:py:func:`~fmriprep.workflows.bold.stc.init_bold_stc_wf`
.. workflow::
:graph2use: colored
:simple_form: yes
from fmriprep.workflows.bold import init_bold_stc_wf
wf = init_bold_stc_wf(
metadata={'RepetitionTime': 2.0,
'SliceTiming': [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]},
)
If the SliceTiming field is available within the input dataset metadata,
this workflow performs slice time correction prior to other signal resampling
processes.
Slice time correction is performed using AFNI 3dTShift.
All slices are realigned in time to the middle of each TR.
Slice time correction can be disabled with the --ignore slicetiming
command line argument.
If a :abbr:`BOLD (blood-oxygen level-dependent)` series has fewer than
5 usable (steady-state) volumes, slice time correction will be disabled
for that run.
One of the major problems that affects :abbr:`EPI (echo planar imaging)` data is the spatial distortion caused by the inhomogeneity of the field inside the scanner.
Applying susceptibility-derived distortion correction, based on fieldmap estimation.
Please note that all routines for susceptibility-derived distortion correction have been excised off of fMRIPrep for utilization on other projects (e.g., dMRIPrep). For more detailed documentation on :abbr:`SDC (susceptibility-derived distortion correction)` routines, check on the SDCFlows component.
Theory, methods and references are found within the SDCFlows documentation.
:py:func:`~fmriprep.workflows.bold.resampling.init_bold_preproc_trans_wf`
.. workflow::
:graph2use: orig
:simple_form: yes
from fmriprep.workflows.bold import init_bold_preproc_trans_wf
wf = init_bold_preproc_trans_wf(mem_gb=3, omp_nthreads=1)
A new preproc :abbr:`BOLD (blood-oxygen level-dependent)` series is generated from the slice-timing corrected or the original data (if :abbr:`STC (slice-timing correction)` was not applied) in the original space. All volumes in the :abbr:`BOLD (blood-oxygen level-dependent)` series are resampled in their native space by concatenating the mappings found in previous correction workflows (:abbr:`HMC (head-motion correction)` and :abbr:`SDC (susceptibility-derived distortion correction)` if executed) for a one-shot interpolation process. Interpolation uses a Lanczos kernel.
:py:func:`~fmriprep.workflows.bold.registration.init_bold_reg_wf`
.. workflow::
:graph2use: orig
:simple_form: yes
from fmriprep.workflows.bold import init_bold_reg_wf
wf = init_bold_reg_wf(
freesurfer=True,
mem_gb=1,
omp_nthreads=1,
use_bbr=True,
bold2t1w_dof=9,
bold2t1w_init='register')
The alignment between the reference :abbr:`EPI (echo-planar imaging)` image
of each run and the reconstructed subject using the gray/white matter boundary
(FreeSurfer's ?h.white surfaces) is calculated by the bbregister routine.
See :func:`fmriprep.workflows.bold.registration.init_bbreg_wf` for further details.
If FreeSurfer processing is disabled, FSL flirt is run with the
:abbr:`BBR (boundary-based registration)` cost function, using the
fast segmentation to establish the gray/white matter boundary.
See :func:`fmriprep.workflows.bold.registration.init_fsl_bbr_wf` for further details.
After either :abbr:`BBR (boundary-based registration)` workflow is run, the resulting affine transform will be compared to the initial transform found by FLIRT. Excessive deviation will result in rejecting the BBR refinement and accepting the original, affine registration.
:py:func:`~fmriprep.workflows.bold.resampling.init_bold_std_trans_wf`
.. workflow::
:graph2use: colored
:simple_form: yes
from niworkflows.utils.spaces import SpatialReferences
from fmriprep.workflows.bold import init_bold_std_trans_wf
wf = init_bold_std_trans_wf(
freesurfer=True,
mem_gb=3,
omp_nthreads=1,
spaces=SpatialReferences(
spaces=[('MNI152Lin', {}), ('MNIPediatricAsym', {'cohort': '6'})],
checkpoint=True),
multiecho=False,
)
This sub-workflow concatenates the transforms calculated upstream (see
Head-motion estimation, Susceptibility Distortion Correction (SDC) --if
fieldmaps are available--, EPI to T1w registration, and an anatomical-to-standard
transform from Preprocessing of structural MRI) to map the
:abbr:`EPI (echo-planar imaging)`
image to the standard spaces given by the --output-spaces argument
(see :ref:`output-spaces`).
It also maps the T1w-based mask to each of those standard spaces.
Transforms are concatenated and applied all at once, with one interpolation (Lanczos) step, so as little information is lost as possible.
The output space grid can be specified using modifiers to the --output-spaces
argument.
:py:func:`~fmriprep.workflows.bold.resampling.init_bold_surf_wf`
.. workflow::
:graph2use: colored
:simple_form: yes
from fmriprep.workflows.bold import init_bold_surf_wf
wf = init_bold_surf_wf(
mem_gb=1,
surface_spaces=['fsnative', 'fsaverage5'],
medial_surface_nan=False)
If FreeSurfer processing is enabled, the motion-corrected functional series (after single shot resampling to T1w space) is sampled to the surface by averaging across the cortical ribbon. Specifically, at each vertex, the segment normal to the white-matter surface, extending to the pial surface, is sampled at 6 intervals and averaged.
Surfaces are generated for the "subject native" surface, as well as transformed to the
fsaverage template space.
All surface outputs are in GIFTI format.
If CIFTI output is enabled, the motion-corrected functional timeseries (in T1w space) is first
sampled to the high resolution 164k vertex (per hemisphere) fsaverage. Following that,
the resampled timeseries is sampled to HCP Pipelines_'s fsLR mesh (with the left and
right hemisphere aligned) using `Connectome Workbench`_'s -metric-resample to generate a
surface timeseries for each hemisphere. These surfaces are then combined with corresponding
volumetric timeseries to create a CIFTI2 file.
:py:func:`~fmriprep.workflows.bold.confounds.init_bold_confs_wf`
.. workflow::
:graph2use: colored
:simple_form: yes
from fmriprep.workflows.bold.confounds import init_bold_confs_wf
wf = init_bold_confs_wf(
name="discover_wf",
mem_gb=1,
metadata={"RepetitionTime": 2.0,
"SliceTiming": [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]},
regressors_all_comps=False,
regressors_dvars_th=1.5,
regressors_fd_th=0.5,
)
Given a motion-corrected fMRI, a brain mask, mcflirt movement parameters and a
segmentation, the discover_wf sub-workflow calculates potential
confounds per volume.
Calculated confounds include the mean global signal, mean tissue class signal,
tCompCor, aCompCor, Frame-wise Displacement, 6 motion parameters, DVARS, spike regressors,
and, if the --use-aroma flag is enabled, the noise components identified by ICA-AROMA
(those to be removed by the "aggressive" denoising strategy).
Particular details about ICA-AROMA are given below.
:py:func:`~fmriprep.workflows.bold.confounds.init_ica_aroma_wf`
ICA-AROMA denoising is performed in MNI152NLin6Asym space, which is automatically
added to the list of --output-spaces if it was not already requested by the user.
The number of ICA-AROMA components depends on a dimensionality estimate made by
FSL MELODIC.
For datasets with a very short TR and a large number of timepoints, this may result
in an unusually high number of components.
By default, dimensionality is limited to a maximum of 200 components.
To override this upper limit one may specify the number of components to be extracted
with --aroma-melodic-dimensionality.
Further details on the implementation are given within the workflow generation
function (:py:func:`~fmriprep.workflows.bold.confounds.init_ica_aroma_wf`).
Note: non-aggressive AROMA denoising is a fundamentally different procedure from its "aggressive" counterpart and cannot be performed only by using a set of noise regressors (a separate GLM with both noise and signal regressors needs to be used). Therefore instead of regressors, fMRIPrep produces non-aggressive denoised 4D NIFTI files in the MNI space:
*space-MNI152NLin6Asym_desc-smoothAROMAnonaggr_bold.nii.gz
Additionally, the MELODIC mix and noise component indices will
be generated, so non-aggressive denoising can be manually performed in the T1w space with fsl_regfilt, e.g.:
fsl_regfilt -i sub-<subject_label>_task-<task_id>_space-T1w_desc-preproc_bold.nii.gz \
-f $(cat sub-<subject_label>_task-<task_id>_AROMAnoiseICs.csv) \
-d sub-<subject_label>_task-<task_id>_desc-MELODIC_mixing.tsv \
-o sub-<subject_label>_task-<task_id>_space-T1w_desc-AROMAnonaggr_bold.nii.gz
Note: The non-steady state volumes are removed for the determination of components in melodic.
Therefore *MELODIC_mixing.tsv may have zero padded rows to account for the volumes not used in
melodic's estimation of components.
A visualization of the AROMA component classification is also included in the HTML reports.
Maps created with maximum intensity projection (glass brain) with a black brain outline. Right hand side of each map: time series (top in seconds), frequency spectrum (bottom in Hertz). Components classified as signal in green; noise in red.
:py:func:`~fmriprep.workflows.bold.t2s.init_bold_t2s_wf`
If multi-echo :abbr:`BOLD (blood-oxygen level-dependent)` data is supplied, this workflow uses the `tedana`_ `T2* workflow`_ to generate an adaptive T2* map and optimally weighted combination of all supplied single echo time series. This optimally combined time series is then carried forward for all subsequent preprocessing steps.
| [Power2017] | Power JD, Plitt M, Kundu P, Bandettini PA, Martin A (2017) Temporal interpolation alters motion in fMRI scans: Magnitudes and consequences for artifact detection. PLOS ONE 12(9): e0182939. doi:10.1371/journal.pone.0182939. |
| [Brett2001] | Brett M, Leff AP, Rorden C, Ashburner J (2001) Spatial Normalization of Brain Images with Focal Lesions Using Cost Function Masking. NeuroImage 14(2) doi:10.006/nimg.2001.0845. |