Source code for niworkflows.anat.coregistration

# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
#
# Copyright 2021 The NiPreps Developers <nipreps@gmail.com>
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
# We support and encourage derived works from this project, please read
# about our expectations at
#
#     https://www.nipreps.org/community/licensing/
#
"""Workflow for the registration of EPI datasets to anatomical space via reconstructed surfaces."""

from nipype import logging
from nipype.interfaces import utility as niu
from nipype.pipeline import engine as pe

LOGGER = logging.getLogger('workflow')


[docs] def init_bbreg_wf( *, omp_nthreads, debug=False, epi2t1w_init='register', epi2t1w_dof=6, name='bbreg_wf', use_bbr=None, ): """ Build a workflow to run FreeSurfer's ``bbregister``. This workflow uses FreeSurfer's ``bbregister`` to register a EPI image to a T1-weighted structural image. It is a counterpart to :py:func:`~fmriprep.workflows.bold.registration.init_fsl_bbr_wf`, which performs the same task using FSL's FLIRT with a BBR cost function. The ``use_bbr`` option permits a high degree of control over registration. If ``False``, standard, affine coregistration will be performed using FreeSurfer's ``mri_coreg`` tool. If ``True``, ``bbregister`` will be seeded with the initial transform found by ``mri_coreg`` (equivalent to running ``bbregister --init-coreg``). If ``None``, after ``bbregister`` is run, the resulting affine transform will be compared to the initial transform found by ``mri_coreg``. Excessive deviation will result in rejecting the BBR refinement and accepting the original, affine registration. Workflow Graph .. workflow :: :graph2use: orig :simple_form: yes from niworkflows.anat.coregistration import init_bbreg_wf wf = init_bbreg_wf(omp_nthreads=1) Parameters ---------- use_bbr : :obj:`bool` or None Enable/disable boundary-based registration refinement. If ``None``, test BBR result for distortion before accepting. epi2t1w_dof : 6, 9 or 12 Degrees-of-freedom for EPI-T1w registration epi2t1w_init : :obj:`str`, ``"header"`` or ``"register"`` If ``"header"``, use header information for initialization of EPI and T1 images. If ``"register"``, align volumes by their centers. name : :obj:`str`, optional Workflow name (default: ``bbreg_wf``) Inputs ------ in_file Reference EPI image to be registered fsnative2t1w_xfm FSL-style affine matrix translating from FreeSurfer T1.mgz to T1w subjects_dir Sets FreeSurfer's ``$SUBJECTS_DIR`` subject_id FreeSurfer subject ID (must have a corresponding folder in ``$SUBJECTS_DIR``) t1w_brain Unused (see :py:func:`~fmriprep.workflows.bold.registration.init_fsl_bbr_wf`) t1w_dseg Unused (see :py:func:`~fmriprep.workflows.bold.registration.init_fsl_bbr_wf`) Outputs ------- itk_epi_to_t1w Affine transform from the reference EPI to T1w space (ITK format) itk_t1w_to_epi Affine transform from T1w space to EPI space (ITK format) out_report Reportlet for assessing registration quality fallback Boolean indicating whether BBR was rejected (mri_coreg registration returned) """ from ..engine.workflows import LiterateWorkflow as Workflow # See https://github.com/nipreps/fmriprep/issues/768 from ..interfaces.freesurfer import ( PatchedBBRegisterRPT as BBRegisterRPT, ) from ..interfaces.freesurfer import ( PatchedLTAConvert as LTAConvert, ) from ..interfaces.freesurfer import ( PatchedMRICoregRPT as MRICoregRPT, ) from ..interfaces.nitransforms import ConcatenateXFMs workflow = Workflow(name=name) workflow.__desc__ = """\ The EPI reference was then co-registered to the T1w reference using `bbregister` (FreeSurfer) which implements boundary-based registration [@bbr]. Co-registration was configured with {dof} degrees of freedom{reason}. """.format( dof={6: 'six', 9: 'nine', 12: 'twelve'}[epi2t1w_dof], reason='' if epi2t1w_dof == 6 else 'to account for distortions remaining in the EPI reference', ) inputnode = pe.Node( niu.IdentityInterface( [ 'in_file', 'fsnative2t1w_xfm', 'subjects_dir', 'subject_id', # BBRegister 't1w_dseg', # FLIRT BBR 't1w_brain', # FLIRT BBR ] ), name='inputnode', ) outputnode = pe.Node( niu.IdentityInterface(['itk_epi_to_t1w', 'itk_t1w_to_epi', 'out_report', 'fallback']), name='outputnode', ) if epi2t1w_init not in ('register', 'header'): raise ValueError(f'Unknown EPI-T1w initialization option: {epi2t1w_init}') # For now make BBR unconditional - in the future, we can fall back to identity, # but adding the flexibility without testing seems a bit dangerous if epi2t1w_init == 'header': if use_bbr is False: raise ValueError('Cannot disable BBR and use header registration') if use_bbr is None: LOGGER.warning('Initializing BBR with header; affine fallback disabled') use_bbr = True merge_ltas = pe.Node(niu.Merge(2), name='merge_ltas', run_without_submitting=True) concat_xfm = pe.Node(ConcatenateXFMs(inverse=True), name='concat_xfm') # fmt:off workflow.connect([ # Output ITK transforms (inputnode, merge_ltas, [('fsnative2t1w_xfm', 'in2')]), (merge_ltas, concat_xfm, [('out', 'in_xfms')]), (concat_xfm, outputnode, [('out_xfm', 'itk_epi_to_t1w')]), (concat_xfm, outputnode, [('out_inv', 'itk_t1w_to_epi')]), ]) # fmt:on if debug is True: from ..interfaces.nibabel import RegridToZooms downsample = pe.Node(RegridToZooms(zooms=(4.0, 4.0, 4.0), smooth=True), name='downsample') workflow.connect([(inputnode, downsample, [('in_file', 'in_file')])]) mri_coreg = pe.Node( MRICoregRPT( dof=epi2t1w_dof, sep=[4], ftol=0.0001, linmintol=0.01, generate_report=not use_bbr, ), name='mri_coreg', n_procs=omp_nthreads, mem_gb=5, ) # Use mri_coreg if epi2t1w_init == 'register': # fmt:off workflow.connect([ (inputnode, mri_coreg, [('subjects_dir', 'subjects_dir'), ('subject_id', 'subject_id')]), ]) # fmt:on if not debug: workflow.connect(inputnode, 'in_file', mri_coreg, 'source_file') else: workflow.connect(downsample, 'out_file', mri_coreg, 'source_file') # Short-circuit workflow building, use initial registration if use_bbr is False: # fmt:off workflow.connect([ (mri_coreg, outputnode, [('out_report', 'out_report')]), (mri_coreg, merge_ltas, [('out_lta_file', 'in1')]), ]) # fmt:on outputnode.inputs.fallback = True return workflow # Use bbregister bbregister = pe.Node( BBRegisterRPT( dof=epi2t1w_dof, contrast_type='t2', registered_file=True, out_lta_file=True, generate_report=True, ), name='bbregister', mem_gb=12, ) # fmt:off workflow.connect([ (inputnode, bbregister, [('subjects_dir', 'subjects_dir'), ('subject_id', 'subject_id')]), ]) # fmt:on if not debug: workflow.connect(inputnode, 'in_file', bbregister, 'source_file') else: workflow.connect(downsample, 'out_file', bbregister, 'source_file') if epi2t1w_init == 'header': bbregister.inputs.init = 'header' else: workflow.connect([(mri_coreg, bbregister, [('out_lta_file', 'init_reg_file')])]) # Short-circuit workflow building, use boundary-based registration if use_bbr is True: # fmt:off workflow.connect([ (bbregister, outputnode, [('out_report', 'out_report')]), (bbregister, merge_ltas, [('out_lta_file', 'in1')]), ]) # fmt:on outputnode.inputs.fallback = False return workflow # Only reach this point if epi2t1w_init is "register" and use_bbr is None transforms = pe.Node(niu.Merge(2), run_without_submitting=True, name='transforms') reports = pe.Node(niu.Merge(2), run_without_submitting=True, name='reports') lta_ras2ras = pe.MapNode( LTAConvert(out_lta=True), iterfield=['in_lta'], name='lta_ras2ras', mem_gb=2 ) compare_transforms = pe.Node(niu.Function(function=compare_xforms), name='compare_transforms') select_transform = pe.Node(niu.Select(), run_without_submitting=True, name='select_transform') select_report = pe.Node(niu.Select(), run_without_submitting=True, name='select_report') # fmt:off workflow.connect([ (bbregister, transforms, [('out_lta_file', 'in1')]), (mri_coreg, transforms, [('out_lta_file', 'in2')]), # Normalize LTA transforms to RAS2RAS (inputs are VOX2VOX) and compare (transforms, lta_ras2ras, [('out', 'in_lta')]), (lta_ras2ras, compare_transforms, [('out_lta', 'lta_list')]), (compare_transforms, outputnode, [('out', 'fallback')]), # Select output transform (transforms, select_transform, [('out', 'inlist')]), (compare_transforms, select_transform, [('out', 'index')]), (select_transform, merge_ltas, [('out', 'in1')]), # Select output report (bbregister, reports, [('out_report', 'in1')]), (mri_coreg, reports, [('out_report', 'in2')]), (reports, select_report, [('out', 'inlist')]), (compare_transforms, select_report, [('out', 'index')]), (select_report, outputnode, [('out', 'out_report')]), ]) # fmt:on return workflow
[docs] def compare_xforms(lta_list, norm_threshold=15): """ Determine a distance between two affine transformations. Computes a normalized displacement between two affine transforms as the maximum overall displacement of the midpoints of the faces of a cube, when each transform is applied to the cube. This combines displacement resulting from scaling, translation and rotation. Although the norm is in mm, in a scaling context, it is not necessarily equivalent to that distance in translation. We choose a default threshold of 15mm as a rough heuristic. Normalized displacement above 20mm showed clear signs of distortion, while "good" BBR refinements were frequently below 10mm displaced from the rigid transform. The 10-20mm range was more ambiguous, and 15mm chosen as a compromise. This is open to revisiting in either direction. See discussion in `GitHub issue #681 <https://github.com/nipreps/fmriprep/issues/681>`__ and the `underlying implementation <https://github.com/nipy/nipype/blob/56b7c81eedeeae884ba47c80096a5f66bd9f8116/nipype/algorithms/rapidart.py#L108-L159>`__. Parameters ---------- lta_list : :obj:`list` or :obj:`tuple` of :obj:`str` the two given affines in LTA format norm_threshold : :obj:`float` the upper bound limit to the normalized displacement caused by the second transform relative to the first (default: `15`) """ from nipype.algorithms.rapidart import _calc_norm_affine from niworkflows.interfaces.surf import load_transform bbr_affine = load_transform(lta_list[0]) fallback_affine = load_transform(lta_list[1]) norm, _ = _calc_norm_affine([fallback_affine, bbr_affine], use_differences=True) return norm[1] > norm_threshold