Source code for fmriprep.workflows.anatomical

# -*- 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

"""

import os.path as op

from pkg_resources import resource_filename as pkgr

from niworkflows.nipype.pipeline import engine as pe
from niworkflows.nipype.interfaces import (
    io as nio,
    utility as niu,
    c3,
    freesurfer as fs,
    fsl
)
from niworkflows.nipype.interfaces.ants import BrainExtraction, N4BiasFieldCorrection

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, MakeMidthickness, FSInjectBrainExtracted,
    FSDetectInputs, NormalizeSurf, GiftiNameSource, TemplateDimensions, Conform, Reorient,
    ConcatAffines, RefineBrainMask,
)
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', 't1_aseg', 't1_aparc']), name='outputnode') buffernode = pe.Node(niu.IdentityInterface( fields=['t1_brain', 't1_mask']), name='buffernode') 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) workflow.connect([ (inputnode, anat_template_wf, [('t1w', 'inputnode.t1w')]), (anat_template_wf, skullstrip_ants_wf, [('outputnode.t1_template', 'inputnode.in_file')]), (skullstrip_ants_wf, outputnode, [('outputnode.bias_corrected', 't1_preproc')]), (anat_template_wf, outputnode, [ ('outputnode.template_transforms', 't1_template_transforms')]), (buffernode, outputnode, [('t1_brain', 't1_brain'), ('t1_mask', 't1_mask')]), ]) # 4. Surface reconstruction if freesurfer: surface_recon_wf = init_surface_recon_wf(name='surface_recon_wf', omp_nthreads=omp_nthreads, hires=hires) applyrefined = pe.Node(fsl.ApplyMask(), name='applyrefined') 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'), ('outputnode.out_segs', 'inputnode.ants_segs'), ('outputnode.bias_corrected', 'inputnode.corrected_t1')]), (skullstrip_ants_wf, applyrefined, [ ('outputnode.bias_corrected', 'in_file')]), (surface_recon_wf, applyrefined, [ ('outputnode.out_brainmask', 'mask_file')]), (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'), ('outputnode.out_aseg', 't1_aseg'), ('outputnode.out_aparc', 't1_aparc')]), (applyrefined, buffernode, [('out_file', 't1_brain')]), (surface_recon_wf, buffernode, [ ('outputnode.out_brainmask', 't1_mask')]), ]) else: workflow.connect([ (skullstrip_ants_wf, buffernode, [ ('outputnode.out_file', 't1_brain'), ('outputnode.out_mask', 't1_mask')]), ]) # 5. Segmentation t1_seg = pe.Node(fsl.FAST(segments=True, no_bias=True, probability_maps=True), name='t1_seg', mem_gb=3) workflow.connect([ (buffernode, t1_seg, [('t1_brain', 'in_files')]), (t1_seg, outputnode, [('tissue_class_map', 't1_seg'), ('probability_maps', 't1_tpms')]), ]) # 6. 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='MultiLabel'), name='mni_mask' ) mni_seg = pe.Node( ApplyTransforms(dimension=3, default_value=0, float=True, interpolation='MultiLabel'), name='mni_seg' ) mni_tpms = pe.MapNode( ApplyTransforms(dimension=3, default_value=0, float=True, interpolation='Linear'), iterfield=['input_image'], name='mni_tpms' ) 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')]), (buffernode, t1_2_mni, [('t1_mask', 'moving_mask')]), (buffernode, mni_mask, [('t1_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')]), ]) 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'), ]), ]) if freesurfer: workflow.connect([ (surface_recon_wf, anat_derivatives_wf, [ ('outputnode.out_aseg', 'inputnode.t1_fs_aseg'), ('outputnode.out_aparc', 'inputnode.t1_fs_aparc'), ]), ])
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') 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_template_dimensions, outputnode, [('out_report', 'out_report'), ('t1w_valid_list', 't1w_valid_list')]), ]) if num_t1w == 1: def _get_first(in_list): if isinstance(in_list, (list, tuple)): return in_list[0] return in_list outputnode.inputs.template_transforms = [pkgr('fmriprep', 'data/itkIdentityTransform.txt')] workflow.connect([ (t1_conform, outputnode, [(('out_file', _get_first), 't1_template')]), ]) return workflow # 1. Template (only if several T1w images) # 1a. Correct for bias field: the bias field is an additive factor # in log-transformed intensity units. Therefore, it is not a linear # combination of fields and N4 fails with merged images. # 1b. Align and merge if several T1w images are provided n4_correct = pe.MapNode( N4BiasFieldCorrection(dimension=3, copy_header=True), iterfield='input_image', name='n4_correct', n_procs=1) # n_procs=1 for reproducibility t1_merge = pe.Node( fs.RobustTemplate(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') # 2. 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([ (t1_conform, n4_correct, [('out_file', 'input_image')]), (t1_conform, t1_merge, [ (('out_file', _set_threads, omp_nthreads), 'num_threads'), (('out_file', add_suffix, '_template'), 'out_file')]), (n4_correct, t1_merge, [('output_image', 'in_files')]), (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_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_segs', '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'), ('BrainExtractionSegmentation', 'out_segs'), ('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. This workflow refines the external brainmask using the internal mask implicit the the FreeSurfer's ``aseg.mgz`` segmentation, to reconcile ANTs' and FreeSurfer's brain masks. First, the ``aseg.mgz`` mask from FreeSurfer is refined in two steps, using binary morphological operations: 1. With a binary closing operation the sulci are included into the mask. This results in a smoother brain mask that does not exclude deep, wide sulci. 2. Fill any holes (typically, there could be a hole next to the pineal gland and the corpora quadrigemina if the great cerebral brain is segmented out). Second, the brain mask is grown, including pixels that have a high likelihood to the GM tissue distribution: 3. Dilate and substract the brain mask, defining the region to search for candidate pixels that likely belong to cortical GM. 4. Pixels found in the search region that are labeled as GM by ANTs (during ``antsBrainExtraction.sh``) are directly added to the new mask. 5. Otherwise, estimate GM tissue parameters locally in patches of ``ww`` size, and test the likelihood of the pixel to belong in the GM distribution. This procedure is inspired on mindboggle's solution to the problem: https://github.com/nipy/mindboggle/blob/7f91faaa7664d820fe12ccc52ebaf21d679795e2/mindboggle/guts/segment.py#L1660 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) ants_segs Brain tissue segmentation from ANTS ``antsBrainExtraction.sh`` corrected_t1 INU-corrected, merged T1-weighted 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_brainmask Refined brainmask, derived from FreeSurfer's ``aseg`` volume out_aseg FreeSurfer's aseg segmentation, in native T1w space out_aparc FreeSurfer's aparc+aseg segmentation, in native T1w space 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', 'corrected_t1', 'ants_segs', '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_brainmask', 'out_aseg', 'out_aparc', '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() aseg_to_native_wf = init_segs_to_native_wf() aparc_to_native_wf = init_segs_to_native_wf(segmentation='aparc_aseg') refine = pe.Node(RefineBrainMask(), name='refine') 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')]), # Refine ANTs mask, deriving new mask from FS' aseg (inputnode, refine, [('corrected_t1', 'in_anat'), ('ants_segs', 'in_ants')]), (inputnode, aseg_to_native_wf, [('corrected_t1', 'inputnode.in_file')]), (autorecon_resume_wf, aseg_to_native_wf, [ ('outputnode.subjects_dir', 'inputnode.subjects_dir'), ('outputnode.subject_id', 'inputnode.subject_id')]), (inputnode, aparc_to_native_wf, [('corrected_t1', 'inputnode.in_file')]), (autorecon_resume_wf, aparc_to_native_wf, [ ('outputnode.subjects_dir', 'inputnode.subjects_dir'), ('outputnode.subject_id', 'inputnode.subject_id')]), (aseg_to_native_wf, refine, [('outputnode.out_file', 'in_aseg')]), # 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')]), (refine, outputnode, [('out_file', 'out_brainmask')]), (aseg_to_native_wf, outputnode, [('outputnode.out_file', 'out_aseg')]), (aparc_to_native_wf, outputnode, [('outputnode.out_file', 'out_aparc')]), ])
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_segs_to_native_wf(name='segs_to_native', segmentation='aseg'): """ Get a segmentation from FreeSurfer conformed space into native T1w space .. workflow:: :graph2use: orig :simple_form: yes from fmriprep.workflows.anatomical import init_segs_to_native_wf wf = init_segs_to_native_wf() **Parameters** segmentation The name of a segmentation ('aseg' or 'aparc_aseg' or 'wmparc') **Inputs** in_file Anatomical, merged T1w image after INU correction subjects_dir FreeSurfer SUBJECTS_DIR subject_id FreeSurfer subject ID **Outputs** out_file The selected segmentation, after resampling in native space """ workflow = pe.Workflow(name='%s_%s' % (name, segmentation)) inputnode = pe.Node(niu.IdentityInterface([ 'in_file', 'subjects_dir', 'subject_id']), name='inputnode') outputnode = pe.Node(niu.IdentityInterface(['out_file']), name='outputnode') # Extract the aseg and aparc+aseg outputs fssource = pe.Node(nio.FreeSurferSource(), name='fs_datasource') tonative = pe.Node(fs.Label2Vol(), name='tonative') tonii = pe.Node(fs.MRIConvert(out_type='niigz', resample_type='nearest'), name='tonii') if segmentation.startswith('aparc'): if segmentation == 'aparc_aseg': def _sel(x): return x[0] elif segmentation == 'aparc_a2009s': def _sel(x): return x[1] elif segmentation == 'aparc_dkt': def _sel(x): return x[2] segmentation = (segmentation, _sel) workflow.connect([ (inputnode, fssource, [ ('subjects_dir', 'subjects_dir'), ('subject_id', 'subject_id')]), (inputnode, tonii, [('in_file', 'reslice_like')]), (fssource, tonative, [(segmentation, 'seg_file'), ('rawavg', 'template_file'), ('aseg', 'reg_header')]), (tonative, tonii, [('vol_label_file', 'in_file')]), (tonii, outputnode, [('out_file', 'out_file')]), ]) 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', 't1_fs_aseg', 't1_fs_aparc']), 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: ds_t1_fsaseg = pe.Node( DerivativesDataSink(base_directory=output_dir, suffix='label-aseg_roi'), name='ds_t1_fsaseg', run_without_submitting=True) ds_t1_fsparc = pe.Node( DerivativesDataSink(base_directory=output_dir, suffix='label-aparcaseg_roi'), name='ds_t1_fsparc', run_without_submitting=True) 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')]), (inputnode, ds_t1_fsaseg, [('t1_fs_aseg', 'in_file')]), (inputnode, ds_t1_fsparc, [('t1_fs_aparc', 'in_file')]), (t1_name, ds_t1_fsaseg, [('out', 'source_file')]), (t1_name, ds_t1_fsparc, [('out', 'source_file')]), ]) 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, newpath=None): """Converts labels to masks""" import nibabel as nb import numpy as np from niworkflows.nipype.utils.filemanip import fname_presuffix 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(fname_presuffix( in_file, suffix='_label%03d' % i, newpath=newpath)) nii.__class__(ldata, nii.affine, nii.header).to_filename(out_files[-1]) return out_files