#!/usr/bin/env python
# -*- coding: utf-8 -*-
# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
"""
Anatomical reference preprocessing workflows
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
.. autofunction:: init_anat_preproc_wf
.. autofunction:: init_skullstrip_ants_wf
Surface preprocessing
+++++++++++++++++++++
``fmriprep`` uses FreeSurfer_ to reconstruct surfaces from T1w/T2w
structural images.
.. autofunction:: init_surface_recon_wf
.. autofunction:: init_autorecon_resume_wf
.. autofunction:: init_gifti_surface_wf
"""
# Originally coded by Craig Moodie. Refactored by the CRN Developers.
import os.path as op
from niworkflows.nipype.interfaces import c3, freesurfer as fs
from niworkflows.nipype.interfaces import utility as niu
from niworkflows.nipype.interfaces import io as nio
from niworkflows.nipype.interfaces import fsl
from niworkflows.nipype.interfaces.ants import BrainExtraction
from niworkflows.nipype.pipeline import engine as pe
from niworkflows.interfaces.registration import RobustMNINormalizationRPT
import niworkflows.data as nid
from niworkflows.interfaces.masks import ROIsPlot
from niworkflows.interfaces.segmentation import ReconAllRPT
from niworkflows.interfaces.fixes import FixHeaderApplyTransforms as ApplyTransforms
from ..interfaces import (
DerivativesDataSink, StructuralReference, MakeMidthickness, FSInjectBrainExtracted,
FSDetectInputs, NormalizeSurf, GiftiNameSource, TemplateDimensions, Conform, Reorient,
ConcatAffines
)
from ..utils.misc import fix_multi_T1w_source_name, add_suffix
# pylint: disable=R0914
[docs]def init_anat_preproc_wf(skull_strip_template, output_spaces, template, debug,
freesurfer, longitudinal, omp_nthreads, hires, reportlets_dir,
output_dir, num_t1w,
name='anat_preproc_wf'):
r"""
This workflow controls the anatomical preprocessing stages of FMRIPREP.
This includes:
- Creation of a structural template
- Skull-stripping and bias correction
- Tissue segmentation
- Normalization
- Surface reconstruction with FreeSurfer
.. workflow::
:graph2use: orig
:simple_form: yes
from fmriprep.workflows.anatomical import init_anat_preproc_wf
wf = init_anat_preproc_wf(omp_nthreads=1,
reportlets_dir='.',
output_dir='.',
template='MNI152NLin2009cAsym',
output_spaces=['T1w', 'fsnative',
'template', 'fsaverage5'],
skull_strip_template='OASIS',
freesurfer=True,
longitudinal=False,
debug=False,
hires=True,
num_t1w=1)
**Parameters**
skull_strip_template : str
Name of ANTs skull-stripping template ('OASIS' or 'NKI')
output_spaces : list
List of output spaces functional images are to be resampled to.
Some pipeline components will only be instantiated for some output spaces.
Valid spaces:
- T1w
- template
- fsnative
- fsaverage (or other pre-existing FreeSurfer templates)
template : str
Name of template targeted by `'template'` output space
debug : bool
Enable debugging outputs
freesurfer : bool
Enable FreeSurfer surface reconstruction (may increase runtime)
longitudinal : bool
Create unbiased structural template, regardless of number of inputs
(may increase runtime)
omp_nthreads : int
Maximum number of threads an individual process may use
hires : bool
Enable sub-millimeter preprocessing in FreeSurfer
reportlets_dir : str
Directory in which to save reportlets
output_dir : str
Directory in which to save derivatives
name : str, optional
Workflow name (default: anat_preproc_wf)
**Inputs**
t1w
List of T1-weighted structural images
t2w
List of T2-weighted structural images
subjects_dir
FreeSurfer SUBJECTS_DIR
**Outputs**
t1_preproc
Bias-corrected structural template, defining T1w space
t1_brain
Skull-stripped ``t1_preproc``
t1_mask
Mask of the skull-stripped template image
t1_seg
Segmentation of preprocessed structural image, including
gray-matter (GM), white-matter (WM) and cerebrospinal fluid (CSF)
t1_tpms
List of tissue probability maps in T1w space
t1_2_mni
T1w template, normalized to MNI space
t1_2_mni_forward_transform
ANTs-compatible affine-and-warp transform file
t1_2_mni_reverse_transform
ANTs-compatible affine-and-warp transform file (inverse)
mni_mask
Mask of skull-stripped template, in MNI space
mni_seg
Segmentation, resampled into MNI space
mni_tpms
List of tissue probability maps in MNI space
subjects_dir
FreeSurfer SUBJECTS_DIR
subject_id
FreeSurfer subject ID
t1_2_fsnative_forward_transform
LTA-style affine matrix translating from T1w to FreeSurfer-conformed subject space
t1_2_fsnative_reverse_transform
LTA-style affine matrix translating from FreeSurfer-conformed subject space to T1w
surfaces
GIFTI surfaces (gray/white boundary, midthickness, pial, inflated)
**Subworkflows**
* :py:func:`~fmriprep.workflows.anatomical.init_skullstrip_ants_wf`
* :py:func:`~fmriprep.workflows.anatomical.init_surface_recon_wf`
"""
workflow = pe.Workflow(name=name)
inputnode = pe.Node(
niu.IdentityInterface(fields=['t1w', 't2w', 'subjects_dir', 'subject_id']),
name='inputnode')
outputnode = pe.Node(niu.IdentityInterface(
fields=['t1_preproc', 't1_brain', 't1_mask', 't1_seg', 't1_tpms',
't1_2_mni', 't1_2_mni_forward_transform', 't1_2_mni_reverse_transform',
'mni_mask', 'mni_seg', 'mni_tpms',
'template_transforms',
'subjects_dir', 'subject_id', 't1_2_fsnative_forward_transform',
't1_2_fsnative_reverse_transform', 'surfaces']),
name='outputnode')
anat_template_wf = init_anat_template_wf(longitudinal=longitudinal, omp_nthreads=omp_nthreads,
num_t1w=num_t1w)
# 3. Skull-stripping
# Bias field correction is handled in skull strip workflows.
skullstrip_ants_wf = init_skullstrip_ants_wf(name='skullstrip_ants_wf',
skull_strip_template=skull_strip_template,
debug=debug,
omp_nthreads=omp_nthreads)
# 4. Segmentation
t1_seg = pe.Node(fsl.FAST(segments=True, no_bias=True, probability_maps=True),
name='t1_seg', mem_gb=3)
# 5. Spatial normalization (T1w to MNI registration)
t1_2_mni = pe.Node(
RobustMNINormalizationRPT(
float=True,
generate_report=True,
flavor='testing' if debug else 'precise',
),
name='t1_2_mni',
n_procs=omp_nthreads,
mem_gb=2
)
# Resample the brain mask and the tissue probability maps into mni space
mni_mask = pe.Node(
ApplyTransforms(dimension=3, default_value=0, float=True,
interpolation='NearestNeighbor'),
name='mni_mask'
)
mni_seg = pe.Node(
ApplyTransforms(dimension=3, default_value=0, float=True,
interpolation='NearestNeighbor'),
name='mni_seg'
)
mni_tpms = pe.MapNode(
ApplyTransforms(dimension=3, default_value=0, float=True,
interpolation='Linear'),
iterfield=['input_image'],
name='mni_tpms'
)
workflow.connect([
(inputnode, anat_template_wf, [('t1w', 'inputnode.t1w')]),
(anat_template_wf, skullstrip_ants_wf, [('outputnode.t1_template', 'inputnode.in_file')]),
(skullstrip_ants_wf, t1_seg, [('outputnode.out_file', 'in_files')]),
(skullstrip_ants_wf, outputnode, [('outputnode.bias_corrected', 't1_preproc'),
('outputnode.out_file', 't1_brain'),
('outputnode.out_mask', 't1_mask')]),
(t1_seg, outputnode, [('tissue_class_map', 't1_seg'),
('probability_maps', 't1_tpms')]),
(anat_template_wf, outputnode, [
('outputnode.template_transforms', 't1_template_transforms')]),
])
if 'template' in output_spaces:
template_str = nid.TEMPLATE_MAP[template]
ref_img = op.join(nid.get_dataset(template_str), '1mm_T1.nii.gz')
t1_2_mni.inputs.template = template_str
mni_mask.inputs.reference_image = ref_img
mni_seg.inputs.reference_image = ref_img
mni_tpms.inputs.reference_image = ref_img
workflow.connect([
(skullstrip_ants_wf, t1_2_mni, [('outputnode.bias_corrected', 'moving_image')]),
(skullstrip_ants_wf, t1_2_mni, [('outputnode.out_mask', 'moving_mask')]),
(skullstrip_ants_wf, mni_mask, [('outputnode.out_mask', 'input_image')]),
(t1_2_mni, mni_mask, [('composite_transform', 'transforms')]),
(t1_seg, mni_seg, [('tissue_class_map', 'input_image')]),
(t1_2_mni, mni_seg, [('composite_transform', 'transforms')]),
(t1_seg, mni_tpms, [('probability_maps', 'input_image')]),
(t1_2_mni, mni_tpms, [('composite_transform', 'transforms')]),
(t1_2_mni, outputnode, [
('warped_image', 't1_2_mni'),
('composite_transform', 't1_2_mni_forward_transform'),
('inverse_composite_transform', 't1_2_mni_reverse_transform')]),
(mni_mask, outputnode, [('output_image', 'mni_mask')]),
(mni_seg, outputnode, [('output_image', 'mni_seg')]),
(mni_tpms, outputnode, [('output_image', 'mni_tpms')]),
])
# 6. FreeSurfer reconstruction
if freesurfer:
surface_recon_wf = init_surface_recon_wf(name='surface_recon_wf',
omp_nthreads=omp_nthreads, hires=hires)
workflow.connect([
(inputnode, surface_recon_wf, [
('t2w', 'inputnode.t2w'),
('subjects_dir', 'inputnode.subjects_dir'),
('subject_id', 'inputnode.subject_id')]),
(anat_template_wf, surface_recon_wf, [('outputnode.t1_template', 'inputnode.t1w')]),
(skullstrip_ants_wf, surface_recon_wf, [
('outputnode.out_file', 'inputnode.skullstripped_t1')]),
(surface_recon_wf, outputnode, [
('outputnode.subjects_dir', 'subjects_dir'),
('outputnode.subject_id', 'subject_id'),
('outputnode.t1_2_fsnative_forward_transform', 't1_2_fsnative_forward_transform'),
('outputnode.t1_2_fsnative_reverse_transform', 't1_2_fsnative_reverse_transform'),
('outputnode.surfaces', 'surfaces')]),
])
seg2msks = pe.Node(niu.Function(function=_seg2msks), name='seg2msks')
seg_rpt = pe.Node(ROIsPlot(colors=['r', 'magenta', 'b', 'g']), name='seg_rpt')
anat_reports_wf = init_anat_reports_wf(
reportlets_dir=reportlets_dir, output_spaces=output_spaces, template=template,
freesurfer=freesurfer)
workflow.connect([
(inputnode, anat_reports_wf, [
(('t1w', fix_multi_T1w_source_name), 'inputnode.source_file')]),
(anat_template_wf, anat_reports_wf, [
('outputnode.out_report', 'inputnode.t1_conform_report')]),
(anat_template_wf, seg_rpt, [
('outputnode.t1_template', 'in_file')]),
(t1_seg, seg2msks, [('tissue_class_map', 'in_file')]),
(seg2msks, seg_rpt, [('out', 'in_rois')]),
(outputnode, seg_rpt, [('t1_mask', 'in_mask')]),
(seg_rpt, anat_reports_wf, [('out_report', 'inputnode.seg_report')]),
])
if freesurfer:
workflow.connect([
(surface_recon_wf, anat_reports_wf, [
('outputnode.out_report', 'inputnode.recon_report')])
])
if 'template' in output_spaces:
workflow.connect([
(t1_2_mni, anat_reports_wf, [('out_report', 'inputnode.t1_2_mni_report')]),
])
anat_derivatives_wf = init_anat_derivatives_wf(output_dir=output_dir,
output_spaces=output_spaces,
template=template,
freesurfer=freesurfer)
workflow.connect([
(anat_template_wf, anat_derivatives_wf, [
('outputnode.t1w_valid_list', 'inputnode.source_files')]),
(outputnode, anat_derivatives_wf, [
('t1_template_transforms', 'inputnode.t1_template_transforms'),
('t1_preproc', 'inputnode.t1_preproc'),
('t1_mask', 'inputnode.t1_mask'),
('t1_seg', 'inputnode.t1_seg'),
('t1_tpms', 'inputnode.t1_tpms'),
('t1_2_mni_forward_transform', 'inputnode.t1_2_mni_forward_transform'),
('t1_2_mni_reverse_transform', 'inputnode.t1_2_mni_reverse_transform'),
('t1_2_mni', 'inputnode.t1_2_mni'),
('mni_mask', 'inputnode.mni_mask'),
('mni_seg', 'inputnode.mni_seg'),
('mni_tpms', 'inputnode.mni_tpms'),
('t1_2_fsnative_forward_transform', 'inputnode.t1_2_fsnative_forward_transform'),
('surfaces', 'inputnode.surfaces'),
]),
])
return workflow
def init_anat_template_wf(longitudinal, omp_nthreads, num_t1w, name='anat_template_wf'):
r"""
This workflow generates a canonically oriented structural template from
input T1w images.
.. workflow::
:graph2use: orig
:simple_form: yes
from fmriprep.workflows.anatomical import init_anat_template_wf
wf = init_anat_template_wf(longitudinal=False, omp_nthreads=1, num_t1w=1)
**Parameters**
longitudinal : bool
Create unbiased structural template, regardless of number of inputs
(may increase runtime)
omp_nthreads : int
Maximum number of threads an individual process may use
name : str, optional
Workflow name (default: anat_template_wf)
**Inputs**
t1w
List of T1-weighted structural images
**Outputs**
t1_template
Structural template, defining T1w space
template_transforms
List of affine transforms from ``t1_template`` to original T1w images
out_report
Conformation report
"""
workflow = pe.Workflow(name=name)
inputnode = pe.Node(niu.IdentityInterface(fields=['t1w']), name='inputnode')
outputnode = pe.Node(niu.IdentityInterface(
fields=['t1_template', 't1w_valid_list', 'template_transforms', 'out_report']),
name='outputnode')
# 0. Reorient T1w image(s) to RAS and resample to common voxel space
t1_template_dimensions = pe.Node(TemplateDimensions(), name='t1_template_dimensions')
t1_conform = pe.MapNode(Conform(), iterfield='in_file', name='t1_conform')
# 1. Align and merge if several T1w images are provided
t1_merge = pe.Node(
# StructuralReference is fs.RobustTemplate if > 1 volume, copying otherwise
StructuralReference(auto_detect_sensitivity=True,
initial_timepoint=1, # For deterministic behavior
intensity_scaling=True, # 7-DOF (rigid + intensity)
subsample_threshold=200,
fixed_timepoint=not longitudinal,
no_iteration=not longitudinal,
transform_outputs=True,
),
mem_gb=2 * num_t1w - 1,
name='t1_merge')
# Reorient template to RAS, if needed (mri_robust_template may set to LIA)
t1_reorient = pe.Node(Reorient(), name='t1_reorient')
lta_to_fsl = pe.MapNode(fs.utils.LTAConvert(out_fsl=True), iterfield=['in_lta'],
name='lta_to_fsl')
concat_affines = pe.MapNode(
ConcatAffines(3, invert=True), iterfield=['mat_AtoB', 'mat_BtoC'],
name='concat_affines', run_without_submitting=True)
fsl_to_itk = pe.MapNode(c3.C3dAffineTool(fsl2ras=True, itk_transform=True),
iterfield=['transform_file', 'source_file'], name='fsl_to_itk')
def set_threads(in_list, maximum):
return min(len(in_list), maximum)
workflow.connect([
(inputnode, t1_template_dimensions, [('t1w', 't1w_list')]),
(t1_template_dimensions, t1_conform, [
('t1w_valid_list', 'in_file'),
('target_zooms', 'target_zooms'),
('target_shape', 'target_shape')]),
(t1_conform, t1_merge, [
('out_file', 'in_files'),
(('out_file', set_threads, omp_nthreads), 'num_threads'),
(('out_file', add_suffix, '_template'), 'out_file')]),
(t1_merge, t1_reorient, [('out_file', 'in_file')]),
# Combine orientation and template transforms
(t1_merge, lta_to_fsl, [('transform_outputs', 'in_lta')]),
(t1_conform, concat_affines, [('transform', 'mat_AtoB')]),
(lta_to_fsl, concat_affines, [('out_fsl', 'mat_BtoC')]),
(t1_reorient, concat_affines, [('transform', 'mat_CtoD')]),
(t1_template_dimensions, fsl_to_itk, [('t1w_valid_list', 'source_file')]),
(t1_reorient, fsl_to_itk, [('out_file', 'reference_file')]),
(concat_affines, fsl_to_itk, [('out_mat', 'transform_file')]),
# Output
(t1_template_dimensions, outputnode, [('out_report', 'out_report'),
('t1w_valid_list', 't1w_valid_list')]),
(t1_reorient, outputnode, [('out_file', 't1_template')]),
(fsl_to_itk, outputnode, [('itk_transform', 'template_transforms')]),
])
return workflow
[docs]def init_skullstrip_ants_wf(skull_strip_template, debug, omp_nthreads, name='skullstrip_ants_wf'):
r"""
This workflow performs skull-stripping using ANTs' ``BrainExtraction.sh``
.. workflow::
:graph2use: orig
:simple_form: yes
from fmriprep.workflows.anatomical import init_skullstrip_ants_wf
wf = init_skullstrip_ants_wf(skull_strip_template='OASIS', debug=False, omp_nthreads=1)
**Parameters**
skull_strip_template : str
Name of ANTs skull-stripping template ('OASIS' or 'NKI')
debug : bool
Enable debugging outputs
omp_nthreads : int
Maximum number of threads an individual process may use
**Inputs**
in_file
T1-weighted structural image to skull-strip
**Outputs**
bias_corrected
Bias-corrected ``in_file``, before skull-stripping
out_file
Skull-stripped ``in_file``
out_mask
Binary mask of the skull-stripped ``in_file``
out_report
Reportlet visualizing quality of skull-stripping
"""
if skull_strip_template == 'OASIS':
from niworkflows.data import get_ants_oasis_template_ras
template_dir = get_ants_oasis_template_ras()
brain_template = op.join(template_dir, 'T_template0.nii.gz')
brain_probability_mask = op.join(
template_dir, 'T_template0_BrainCerebellumProbabilityMask.nii.gz')
extraction_registration_mask = op.join(
template_dir, 'T_template0_BrainCerebellumRegistrationMask.nii.gz')
elif skull_strip_template == 'NKI':
from niworkflows.data.getters import get_ants_nki_template_ras
template_dir = get_ants_nki_template_ras()
brain_template = op.join(template_dir, 'T_template.nii.gz')
brain_probability_mask = op.join(
template_dir, 'T_template_BrainCerebellumProbabilityMask.nii.gz')
extraction_registration_mask = op.join(
template_dir, 'T_template_BrainCerebellumExtractionMask.nii.gz')
else:
raise ValueError("Unknown skull-stripping template; select from {OASIS, NKI}")
workflow = pe.Workflow(name=name)
inputnode = pe.Node(niu.IdentityInterface(fields=['in_file']),
name='inputnode')
outputnode = pe.Node(niu.IdentityInterface(
fields=['bias_corrected', 'out_file', 'out_mask', 'out_report']), name='outputnode')
t1_skull_strip = pe.Node(
BrainExtraction(dimension=3, use_floatingpoint_precision=1, debug=debug,
keep_temporary_files=1),
name='t1_skull_strip', n_procs=omp_nthreads)
t1_skull_strip.inputs.brain_template = brain_template
t1_skull_strip.inputs.brain_probability_mask = brain_probability_mask
t1_skull_strip.inputs.extraction_registration_mask = extraction_registration_mask
workflow.connect([
(inputnode, t1_skull_strip, [('in_file', 'anatomical_image')]),
(t1_skull_strip, outputnode, [('BrainExtractionMask', 'out_mask'),
('BrainExtractionBrain', 'out_file'),
('N4Corrected0', 'bias_corrected')])
])
return workflow
[docs]def init_surface_recon_wf(omp_nthreads, hires, name='surface_recon_wf'):
r"""
This workflow reconstructs anatomical surfaces using FreeSurfer's ``recon-all``.
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.
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 an externally computed skull-stripping mask.
The final phase resumes reconstruction, using the T2w image to assist
in finding the pial surface, if available.
See :py:func:`~fmriprep.workflows.anatomical.init_autorecon_resume_wf` for details.
Memory annotations for FreeSurfer are based off `their documentation
<https://surfer.nmr.mgh.harvard.edu/fswiki/SystemRequirements>`_.
They specify an allocation of 4GB per subject. Here we define 5GB
to have a certain margin.
.. workflow::
:graph2use: orig
:simple_form: yes
from fmriprep.workflows.anatomical import init_surface_recon_wf
wf = init_surface_recon_wf(omp_nthreads=1, hires=True)
**Parameters**
omp_nthreads : int
Maximum number of threads an individual process may use
hires : bool
Enable sub-millimeter preprocessing in FreeSurfer
**Inputs**
t1w
List of T1-weighted structural images
t2w
List of T2-weighted structural images (only first used)
skullstripped_t1
Skull-stripped T1-weighted image (or mask of image)
subjects_dir
FreeSurfer SUBJECTS_DIR
subject_id
FreeSurfer subject ID
**Outputs**
subjects_dir
FreeSurfer SUBJECTS_DIR
subject_id
FreeSurfer subject ID
t1_2_fsnative_forward_transform
LTA-style affine matrix translating from T1w to FreeSurfer-conformed subject space
t1_2_fsnative_reverse_transform
LTA-style affine matrix translating from FreeSurfer-conformed subject space to T1w
surfaces
GIFTI surfaces for gray/white matter boundary, pial surface,
midthickness (or graymid) surface, and inflated surfaces
out_report
Reportlet visualizing quality of surface alignment
**Subworkflows**
* :py:func:`~fmriprep.workflows.anatomical.init_autorecon_resume_wf`
* :py:func:`~fmriprep.workflows.anatomical.init_gifti_surface_wf`
"""
workflow = pe.Workflow(name=name)
inputnode = pe.Node(
niu.IdentityInterface(
fields=['t1w', 't2w', 'skullstripped_t1', 'subjects_dir', 'subject_id']),
name='inputnode')
outputnode = pe.Node(
niu.IdentityInterface(
fields=['subjects_dir', 'subject_id', 't1_2_fsnative_forward_transform',
't1_2_fsnative_reverse_transform', 'surfaces',
'out_report']),
name='outputnode')
recon_config = pe.Node(FSDetectInputs(hires_enabled=hires), name='recon_config')
autorecon1 = pe.Node(
fs.ReconAll(directive='autorecon1', flags='-noskullstrip', openmp=omp_nthreads),
name='autorecon1', n_procs=omp_nthreads, mem_gb=5)
autorecon1.interface._can_resume = False
skull_strip_extern = pe.Node(FSInjectBrainExtracted(), name='skull_strip_extern')
fsnative_2_t1_xfm = pe.Node(fs.RobustRegister(auto_sens=True, est_int_scale=True),
name='fsnative_2_t1_xfm')
t1_2_fsnative_xfm = pe.Node(fs.utils.LTAConvert(out_lta=True, invert=True),
name='t1_2_fsnative_xfm')
autorecon_resume_wf = init_autorecon_resume_wf(omp_nthreads=omp_nthreads)
gifti_surface_wf = init_gifti_surface_wf()
workflow.connect([
# Configuration
(inputnode, recon_config, [('t1w', 't1w_list'),
('t2w', 't2w_list')]),
# Passing subjects_dir / subject_id enforces serial order
(inputnode, autorecon1, [('subjects_dir', 'subjects_dir'),
('subject_id', 'subject_id')]),
(autorecon1, skull_strip_extern, [('subjects_dir', 'subjects_dir'),
('subject_id', 'subject_id')]),
(skull_strip_extern, autorecon_resume_wf, [('subjects_dir', 'inputnode.subjects_dir'),
('subject_id', 'inputnode.subject_id')]),
(autorecon_resume_wf, gifti_surface_wf, [
('outputnode.subjects_dir', 'inputnode.subjects_dir'),
('outputnode.subject_id', 'inputnode.subject_id')]),
# Reconstruction phases
(inputnode, autorecon1, [('t1w', 'T1_files')]),
(recon_config, autorecon1, [('t2w', 'T2_file'),
('hires', 'hires'),
# First run only (recon-all saves expert options)
('mris_inflate', 'mris_inflate')]),
(inputnode, skull_strip_extern, [('skullstripped_t1', 'in_brain')]),
(recon_config, autorecon_resume_wf, [('use_t2w', 'inputnode.use_T2')]),
# Construct transform from FreeSurfer conformed image to FMRIPREP
# reoriented image
(inputnode, fsnative_2_t1_xfm, [('t1w', 'target_file')]),
(autorecon1, fsnative_2_t1_xfm, [('T1', 'source_file')]),
(fsnative_2_t1_xfm, gifti_surface_wf, [
('out_reg_file', 'inputnode.t1_2_fsnative_reverse_transform')]),
(fsnative_2_t1_xfm, t1_2_fsnative_xfm, [('out_reg_file', 'in_lta')]),
# Output
(autorecon_resume_wf, outputnode, [('outputnode.subjects_dir', 'subjects_dir'),
('outputnode.subject_id', 'subject_id'),
('outputnode.out_report', 'out_report')]),
(gifti_surface_wf, outputnode, [('outputnode.surfaces', 'surfaces')]),
(t1_2_fsnative_xfm, outputnode, [('out_lta', 't1_2_fsnative_forward_transform')]),
(fsnative_2_t1_xfm, outputnode, [('out_reg_file', 't1_2_fsnative_reverse_transform')]),
])
return workflow
[docs]def init_autorecon_resume_wf(omp_nthreads, name='autorecon_resume_wf'):
r"""
This workflow resumes recon-all execution, assuming the `-autorecon1` stage
has been completed.
In order to utilize resources efficiently, this is broken down into five
sub-stages; after the first stage, the second and third stages may be run
simultaneously, and the fourth and fifth stages may be run simultaneously,
if resources permit::
$ recon-all -sd <output dir>/freesurfer -subjid sub-<subject_label> \
-autorecon2-volonly
$ recon-all -sd <output dir>/freesurfer -subjid sub-<subject_label> \
-autorecon-hemi lh \
-noparcstats -nocortparc2 -noparcstats2 -nocortparc3 \
-noparcstats3 -nopctsurfcon -nohyporelabel -noaparc2aseg \
-noapas2aseg -nosegstats -nowmparc -nobalabels
$ recon-all -sd <output dir>/freesurfer -subjid sub-<subject_label> \
-autorecon-hemi rh \
-noparcstats -nocortparc2 -noparcstats2 -nocortparc3 \
-noparcstats3 -nopctsurfcon -nohyporelabel -noaparc2aseg \
-noapas2aseg -nosegstats -nowmparc -nobalabels
$ recon-all -sd <output dir>/freesurfer -subjid sub-<subject_label> \
-autorecon3 -hemi lh -T2pial
$ recon-all -sd <output dir>/freesurfer -subjid sub-<subject_label> \
-autorecon3 -hemi rh -T2pial
The excluded steps in the second and third stages (``-no<option>``) are not
fully hemisphere independent, and are therefore postponed to the final two
stages.
.. workflow::
:graph2use: orig
:simple_form: yes
from fmriprep.workflows.anatomical import init_autorecon_resume_wf
wf = init_autorecon_resume_wf(omp_nthreads=1)
**Inputs**
subjects_dir
FreeSurfer SUBJECTS_DIR
subject_id
FreeSurfer subject ID
use_T2
Refine pial surface using T2w images
**Outputs**
subjects_dir
FreeSurfer SUBJECTS_DIR
subject_id
FreeSurfer subject ID
out_report
Reportlet visualizing quality of surface alignment
"""
workflow = pe.Workflow(name=name)
inputnode = pe.Node(
niu.IdentityInterface(
fields=['subjects_dir', 'subject_id', 'use_T2']),
name='inputnode')
outputnode = pe.Node(
niu.IdentityInterface(
fields=['subjects_dir', 'subject_id', 'out_report']),
name='outputnode')
autorecon2_vol = pe.Node(
fs.ReconAll(directive='autorecon2-volonly', openmp=omp_nthreads),
n_procs=omp_nthreads, mem_gb=5, name='autorecon2_vol')
autorecon_surfs = pe.MapNode(
fs.ReconAll(
directive='autorecon-hemi',
flags=['-noparcstats', '-nocortparc2', '-noparcstats2',
'-nocortparc3', '-noparcstats3', '-nopctsurfcon',
'-nohyporelabel', '-noaparc2aseg', '-noapas2aseg',
'-nosegstats', '-nowmparc', '-nobalabels'],
openmp=omp_nthreads),
iterfield='hemi', n_procs=omp_nthreads, mem_gb=5,
name='autorecon_surfs')
autorecon_surfs.inputs.hemi = ['lh', 'rh']
autorecon3 = pe.MapNode(
fs.ReconAll(directive='autorecon3', openmp=omp_nthreads),
iterfield='hemi', n_procs=omp_nthreads, mem_gb=5,
name='autorecon3')
autorecon3.inputs.hemi = ['lh', 'rh']
# Only generate the report once; should be nothing to do
recon_report = pe.Node(
ReconAllRPT(directive='autorecon3', generate_report=True),
name='recon_report', mem_gb=5)
def _dedup(in_list):
vals = set(in_list)
if len(vals) > 1:
raise ValueError(
"Non-identical values can't be deduplicated:\n{!r}".format(in_list))
return vals.pop()
workflow.connect([
(inputnode, autorecon3, [('use_T2', 'use_T2')]),
(inputnode, autorecon2_vol, [('subjects_dir', 'subjects_dir'),
('subject_id', 'subject_id')]),
(autorecon2_vol, autorecon_surfs, [('subjects_dir', 'subjects_dir'),
('subject_id', 'subject_id')]),
(autorecon_surfs, autorecon3, [(('subjects_dir', _dedup), 'subjects_dir'),
(('subject_id', _dedup), 'subject_id')]),
(autorecon3, outputnode, [(('subjects_dir', _dedup), 'subjects_dir'),
(('subject_id', _dedup), 'subject_id')]),
(autorecon3, recon_report, [(('subjects_dir', _dedup), 'subjects_dir'),
(('subject_id', _dedup), 'subject_id')]),
(recon_report, outputnode, [('out_report', 'out_report')]),
])
return workflow
[docs]def init_gifti_surface_wf(name='gifti_surface_wf'):
r"""
This workflow prepares GIFTI surfaces from a FreeSurfer subjects directory
If midthickness (or graymid) surfaces do not exist, they are generated and
saved to the subject directory as ``lh/rh.midthickness``.
These, along with the gray/white matter boundary (``lh/rh.smoothwm``), pial
sufaces (``lh/rh.pial``) and inflated surfaces (``lh/rh.inflated``) are
converted to GIFTI files.
Additionally, the vertex coordinates are :py:class:`recentered
<fmriprep.interfaces.NormalizeSurf>` to align with native T1w space.
.. workflow::
:graph2use: orig
:simple_form: yes
from fmriprep.workflows.anatomical import init_gifti_surface_wf
wf = init_gifti_surface_wf()
**Inputs**
subjects_dir
FreeSurfer SUBJECTS_DIR
subject_id
FreeSurfer subject ID
t1_2_fsnative_reverse_transform
LTA formatted affine transform file (inverse)
**Outputs**
surfaces
GIFTI surfaces for gray/white matter boundary, pial surface,
midthickness (or graymid) surface, and inflated surfaces
"""
workflow = pe.Workflow(name=name)
inputnode = pe.Node(niu.IdentityInterface(['subjects_dir', 'subject_id',
't1_2_fsnative_reverse_transform']),
name='inputnode')
outputnode = pe.Node(niu.IdentityInterface(['surfaces']), name='outputnode')
get_surfaces = pe.Node(nio.FreeSurferSource(), name='get_surfaces')
midthickness = pe.MapNode(
MakeMidthickness(thickness=True, distance=0.5, out_name='midthickness'),
iterfield='in_file',
name='midthickness')
save_midthickness = pe.Node(nio.DataSink(parameterization=False),
name='save_midthickness')
surface_list = pe.Node(niu.Merge(4, ravel_inputs=True),
name='surface_list', run_without_submitting=True)
fs_2_gii = pe.MapNode(fs.MRIsConvert(out_datatype='gii'),
iterfield='in_file', name='fs_2_gii')
fix_surfs = pe.MapNode(NormalizeSurf(), iterfield='in_file', name='fix_surfs')
workflow.connect([
(inputnode, get_surfaces, [('subjects_dir', 'subjects_dir'),
('subject_id', 'subject_id')]),
(inputnode, save_midthickness, [('subjects_dir', 'base_directory'),
('subject_id', 'container')]),
# Generate midthickness surfaces and save to FreeSurfer derivatives
(get_surfaces, midthickness, [('smoothwm', 'in_file'),
('graymid', 'graymid')]),
(midthickness, save_midthickness, [('out_file', 'surf.@graymid')]),
# Produce valid GIFTI surface files (dense mesh)
(get_surfaces, surface_list, [('smoothwm', 'in1'),
('pial', 'in2'),
('inflated', 'in3')]),
(save_midthickness, surface_list, [('out_file', 'in4')]),
(surface_list, fs_2_gii, [('out', 'in_file')]),
(fs_2_gii, fix_surfs, [('converted', 'in_file')]),
(inputnode, fix_surfs, [('t1_2_fsnative_reverse_transform', 'transform_file')]),
(fix_surfs, outputnode, [('out_file', 'surfaces')]),
])
return workflow
def init_anat_reports_wf(reportlets_dir, output_spaces,
template, freesurfer, name='anat_reports_wf'):
"""
Set up a battery of datasinks to store reports in the right location
"""
workflow = pe.Workflow(name=name)
inputnode = pe.Node(
niu.IdentityInterface(
fields=['source_file', 't1_conform_report', 'seg_report',
't1_2_mni_report', 'recon_report']),
name='inputnode')
ds_t1_conform_report = pe.Node(
DerivativesDataSink(base_directory=reportlets_dir, suffix='conform'),
name='ds_t1_conform_report', run_without_submitting=True)
ds_t1_2_mni_report = pe.Node(
DerivativesDataSink(base_directory=reportlets_dir, suffix='t1_2_mni'),
name='ds_t1_2_mni_report', run_without_submitting=True)
ds_t1_seg_mask_report = pe.Node(
DerivativesDataSink(base_directory=reportlets_dir, suffix='seg_brainmask'),
name='ds_t1_seg_mask_report', run_without_submitting=True)
ds_recon_report = pe.Node(
DerivativesDataSink(base_directory=reportlets_dir, suffix='reconall'),
name='ds_recon_report', run_without_submitting=True)
workflow.connect([
(inputnode, ds_t1_conform_report, [('source_file', 'source_file'),
('t1_conform_report', 'in_file')]),
(inputnode, ds_t1_seg_mask_report, [('source_file', 'source_file'),
('seg_report', 'in_file')]),
])
if freesurfer:
workflow.connect([
(inputnode, ds_recon_report, [('source_file', 'source_file'),
('recon_report', 'in_file')])
])
if 'template' in output_spaces:
workflow.connect([
(inputnode, ds_t1_2_mni_report, [('source_file', 'source_file'),
('t1_2_mni_report', 'in_file')])
])
return workflow
def init_anat_derivatives_wf(output_dir, output_spaces, template, freesurfer,
name='anat_derivatives_wf'):
"""
Set up a battery of datasinks to store derivatives in the right location
"""
workflow = pe.Workflow(name=name)
inputnode = pe.Node(
niu.IdentityInterface(
fields=['source_files', 't1_template_transforms',
't1_preproc', 't1_mask', 't1_seg', 't1_tpms',
't1_2_mni_forward_transform', 't1_2_mni_reverse_transform',
't1_2_mni', 'mni_mask', 'mni_seg', 'mni_tpms',
't1_2_fsnative_forward_transform', 'surfaces']),
name='inputnode')
t1_name = pe.Node(niu.Function(function=fix_multi_T1w_source_name), name='t1_name')
ds_t1_preproc = pe.Node(
DerivativesDataSink(base_directory=output_dir, suffix='preproc'),
name='ds_t1_preproc', run_without_submitting=True)
ds_t1_mask = pe.Node(
DerivativesDataSink(base_directory=output_dir, suffix='brainmask'),
name='ds_t1_mask', run_without_submitting=True)
ds_t1_seg = pe.Node(
DerivativesDataSink(base_directory=output_dir, suffix='dtissue'),
name='ds_t1_seg', run_without_submitting=True)
ds_t1_tpms = pe.Node(
DerivativesDataSink(base_directory=output_dir,
suffix='class-{extra_value}_probtissue'),
name='ds_t1_tpms', run_without_submitting=True)
ds_t1_tpms.inputs.extra_values = ['CSF', 'GM', 'WM']
suffix_fmt = 'space-{}_{}'.format
ds_t1_mni = pe.Node(
DerivativesDataSink(base_directory=output_dir,
suffix=suffix_fmt(template, 'preproc')),
name='ds_t1_mni', run_without_submitting=True)
ds_mni_mask = pe.Node(
DerivativesDataSink(base_directory=output_dir,
suffix=suffix_fmt(template, 'brainmask')),
name='ds_mni_mask', run_without_submitting=True)
ds_mni_seg = pe.Node(
DerivativesDataSink(base_directory=output_dir,
suffix=suffix_fmt(template, 'dtissue')),
name='ds_mni_seg', run_without_submitting=True)
ds_mni_tpms = pe.Node(
DerivativesDataSink(base_directory=output_dir,
suffix=suffix_fmt(template, 'class-{extra_value}_probtissue')),
name='ds_mni_tpms', run_without_submitting=True)
ds_mni_tpms.inputs.extra_values = ['CSF', 'GM', 'WM']
# Transforms
suffix_fmt = 'space-{}_target-{}_{}'.format
ds_t1_mni_inv_warp = pe.Node(
DerivativesDataSink(base_directory=output_dir,
suffix=suffix_fmt(template, 'T1w', 'warp')),
name='ds_t1_mni_inv_warp', run_without_submitting=True)
ds_t1_template_transforms = pe.MapNode(
DerivativesDataSink(base_directory=output_dir, suffix=suffix_fmt('orig', 'T1w', 'affine')),
iterfield=['source_file', 'in_file'],
name='ds_t1_template_transforms', run_without_submitting=True)
suffix_fmt = 'target-{}_{}'.format
ds_t1_mni_warp = pe.Node(
DerivativesDataSink(base_directory=output_dir, suffix=suffix_fmt(template, 'warp')),
name='ds_t1_mni_warp', run_without_submitting=True)
lta_2_itk = pe.Node(fs.utils.LTAConvert(out_itk=True), name='lta_2_itk')
ds_t1_fsnative = pe.Node(
DerivativesDataSink(base_directory=output_dir, suffix=suffix_fmt('fsnative', 'affine')),
name='ds_t1_fsnative', run_without_submitting=True)
name_surfs = pe.MapNode(GiftiNameSource(pattern=r'(?P<LR>[lr])h.(?P<surf>.+)_converted.gii',
template='{surf}.{LR}.surf'),
iterfield='in_file',
name='name_surfs',
run_without_submitting=True)
ds_surfs = pe.MapNode(
DerivativesDataSink(base_directory=output_dir),
iterfield=['in_file', 'suffix'], name='ds_surfs', run_without_submitting=True)
workflow.connect([
(inputnode, t1_name, [('source_files', 'in_files')]),
(inputnode, ds_t1_template_transforms, [('source_files', 'source_file'),
('t1_template_transforms', 'in_file')]),
(inputnode, ds_t1_preproc, [('t1_preproc', 'in_file')]),
(inputnode, ds_t1_mask, [('t1_mask', 'in_file')]),
(inputnode, ds_t1_seg, [('t1_seg', 'in_file')]),
(inputnode, ds_t1_tpms, [('t1_tpms', 'in_file')]),
(t1_name, ds_t1_preproc, [('out', 'source_file')]),
(t1_name, ds_t1_mask, [('out', 'source_file')]),
(t1_name, ds_t1_seg, [('out', 'source_file')]),
(t1_name, ds_t1_tpms, [('out', 'source_file')]),
])
if freesurfer:
workflow.connect([
(inputnode, lta_2_itk, [('t1_2_fsnative_forward_transform', 'in_lta')]),
(t1_name, ds_t1_fsnative, [('out', 'source_file')]),
(lta_2_itk, ds_t1_fsnative, [('out_itk', 'in_file')]),
(inputnode, name_surfs, [('surfaces', 'in_file')]),
(inputnode, ds_surfs, [('surfaces', 'in_file')]),
(t1_name, ds_surfs, [('out', 'source_file')]),
(name_surfs, ds_surfs, [('out_name', 'suffix')]),
])
if 'template' in output_spaces:
workflow.connect([
(inputnode, ds_t1_mni_warp, [('t1_2_mni_forward_transform', 'in_file')]),
(inputnode, ds_t1_mni_inv_warp, [('t1_2_mni_reverse_transform', 'in_file')]),
(inputnode, ds_t1_mni, [('t1_2_mni', 'in_file')]),
(inputnode, ds_mni_mask, [('mni_mask', 'in_file')]),
(inputnode, ds_mni_seg, [('mni_seg', 'in_file')]),
(inputnode, ds_mni_tpms, [('mni_tpms', 'in_file')]),
(t1_name, ds_t1_mni_warp, [('out', 'source_file')]),
(t1_name, ds_t1_mni_inv_warp, [('out', 'source_file')]),
(t1_name, ds_t1_mni, [('out', 'source_file')]),
(t1_name, ds_mni_mask, [('out', 'source_file')]),
(t1_name, ds_mni_seg, [('out', 'source_file')]),
(t1_name, ds_mni_tpms, [('out', 'source_file')]),
])
return workflow
def _seg2msks(in_file):
"""Converts labels to masks"""
import nibabel as nb
import numpy as np
from os import path as op
nii = nb.load(in_file)
labels = nii.get_data()
out_files = []
for i in range(1, 4):
ldata = np.zeros_like(labels)
ldata[labels == i] = 1
out_files.append(op.abspath('label%d.nii.gz' % i))
nii.__class__(ldata, nii.affine, nii.header).to_filename(out_files[-1])
return out_files