diff --git a/fmriprep/data/reports-spec-func.yml b/fmriprep/data/reports-spec-func.yml index 04acea2c..cf932014 100644 --- a/fmriprep/data/reports-spec-func.yml +++ b/fmriprep/data/reports-spec-func.yml @@ -104,6 +104,21 @@ sections: anatomical white matter mask, which appears as a red contour. static: false subtitle: Alignment of functional and anatomical MRI data (coregistration) + - button: + type: button + class: btn btn-primary + data-toggle: collapse + data-target: '#flippedcoreg' + text: 'Alignment of functional and flipped anatomical MRI data' + - div: + class: collapse + id: 'flippedcoreg' + reportlets: + - bids: {datatype: figures, desc: flippedcoreg, suffix: bold} + caption: This panel shows the alignment of the reference EPI (BOLD) image to the + left-right flipped anatomical (T1-weighted) image. + static: false + subtitle: Left-right flip check, alignment of functional and flipped anatomical MRI data - bids: {datatype: figures, desc: rois, suffix: bold} caption: Brain mask calculated on the BOLD signal (red contour), along with the regions of interest (ROIs) used for the estimation of physiological and movement diff --git a/fmriprep/data/reports-spec.yml b/fmriprep/data/reports-spec.yml index 4ec02bda..17bc8caf 100644 --- a/fmriprep/data/reports-spec.yml +++ b/fmriprep/data/reports-spec.yml @@ -130,6 +130,21 @@ sections: anatomical white matter mask, which appears as a red contour. static: false subtitle: Alignment of functional and anatomical MRI data (coregistration) + - button: + type: button + class: btn btn-primary + data-toggle: collapse + data-target: '#flippedcoreg' + text: 'Alignment of functional and flipped anatomical MRI data' + - div: + class: collapse + id: 'flippedcoreg' + reportlets: + - bids: {datatype: figures, desc: flippedcoreg, suffix: bold} + caption: This panel shows the alignment of the reference EPI (BOLD) image to the + left-right flipped anatomical (T1-weighted) image. + static: false + subtitle: Left-right flip check, alignment of functional and flipped anatomical MRI data - bids: {datatype: figures, desc: rois, suffix: bold} caption: Brain mask calculated on the BOLD signal (red contour), along with the regions of interest (ROIs) used for the estimation of physiological and movement diff --git a/fmriprep/interfaces/reports.py b/fmriprep/interfaces/reports.py index 8f26c374..76d91a36 100644 --- a/fmriprep/interfaces/reports.py +++ b/fmriprep/interfaces/reports.py @@ -67,6 +67,17 @@ \t\t\t
  • Slice timing correction: {stc}
  • \t\t\t
  • Susceptibility distortion correction: {sdc}
  • \t\t\t
  • Registration: {registration}
  • +\t\t\t
  • Left-right flip check warning: {lr_flip_warning}
  • +\t\t +\t\t\t +\t\t\t\t +\t\t\t\t +\t\t\t +\t\t\t +\t\t\t\t +\t\t\t\t +\t\t\t +\t\t
    Original Registration CostFlipped Registration Cost
    {cost_original}{cost_flipped}
    \t\t\t
  • Non-steady-state volumes: {dummy_scan_desc}
  • \t\t \t\t @@ -218,6 +229,12 @@ class FunctionalSummaryInputSpec(TraitedSpec): desc='Whether to initialize registration with the "header"' ' or by centering the volumes ("t1w" or "t2w")', ) + flip_info = traits.Dict( + traits.Enum('lr_flip_warning', 'cost_original', 'cost_flipped'), + traits.Either(traits.Bool(), traits.Float()), + desc='Left-right flip check warning and registration costs', + mandatory=True, + ) tr = traits.Float(desc='Repetition time', mandatory=True) dummy_scans = traits.Either(traits.Int(), None, desc='number of dummy scans specified by user') algo_dummy_scans = traits.Int(desc='number of dummy scans determined by algorithm') @@ -279,6 +296,12 @@ def _generate_segment(self): if n_echos > 2: multiecho = f'Multi-echo EPI sequence: {n_echos} echoes.' + lr_flip_warning = ( + 'LR flip detected' + if self.inputs.flip_info.get('lr_flip_warning', False) + else 'none' + ) + return FUNCTIONAL_TEMPLATE.format( pedir=pedir, stc=stc, @@ -288,6 +311,9 @@ def _generate_segment(self): dummy_scan_desc=dummy_scan_msg, multiecho=multiecho, ornt=self.inputs.orientation, + lr_flip_warning=lr_flip_warning, + cost_original=self.input.flip_info.get('cost_original', None), + cost_flipped=self.input.flip_info.get('cost_flipped', None), ) @@ -369,3 +395,40 @@ def get_world_pedir(ornt, pe_direction): f'Orientation: {ornt}; PE dir: {pe_direction}' ) return 'Could not be determined - assuming Anterior-Posterior' + + +class _CheckFlipInputSpec(BaseInterfaceInputSpec): + cost_original = File( + exists=True, + mandatory=True, + desc='cost associated with registration of BOLD to original T1w images', + ) + cost_flipped = File( + exists=True, + mandatory=True, + desc='cost associated with registration of BOLD to the flipped T1w images', + ) + + +class _CheckFlipOutputSpec(TraitedSpec): + flip_info = traits.Dict( + traits.Enum('warning', 'cost_original', 'cost_flipped'), + traits.Either(traits.Bool(), traits.Float()), + desc='Left-right flip check warning and registration costs', + mandatory=True, + ) + + +class CheckFlip(SimpleInterface): + """Check for a LR flip by comparing registration cost functions.""" + + input_spec = _CheckFlipInputSpec + output_spec = _CheckFlipOutputSpec + + def _run_interface(self, runtime): + self._results['flip_info'] = { + 'warning': self.inputs.cost_flipped < self.inputs.cost_original, + 'cost_original': self.inputs.cost_original, + 'cost_flipped': self.inputs.cost_flipped, + } + return runtime diff --git a/fmriprep/workflows/bold/fit.py b/fmriprep/workflows/bold/fit.py index f66d49eb..3161c332 100644 --- a/fmriprep/workflows/bold/fit.py +++ b/fmriprep/workflows/bold/fit.py @@ -289,6 +289,9 @@ def init_bold_fit_wf( 'boldref2fmap_xfm', 'movpar_file', 'rmsd_file', + # LR flip check + 'fboldref2anat_xfm', + 'flipped_boldref', ], ), name='outputnode', @@ -373,6 +376,8 @@ def init_bold_fit_wf( ('coreg_boldref', 'inputnode.coreg_boldref'), ('bold_mask', 'inputnode.bold_mask'), ('boldref2anat_xfm', 'inputnode.boldref2anat_xfm'), + ('fboldref2anat_xfm', 'inputnode.fboldref2anat_xfm'), + ('flipped_boldref', 'inputnode.flipped_boldref'), ]), (summary, func_fit_reports_wf, [('out_report', 'inputnode.summary_report')]), ]) @@ -632,7 +637,11 @@ def init_bold_fit_wf( (regref_buffer, ds_boldreg_wf, [('boldref', 'inputnode.source_files')]), (bold_reg_wf, ds_boldreg_wf, [('outputnode.itk_bold_to_t1', 'inputnode.xform')]), (ds_boldreg_wf, outputnode, [('outputnode.xform', 'boldref2anat_xfm')]), - (bold_reg_wf, summary, [('outputnode.fallback', 'fallback')]), + (bold_reg_wf, outputnode, [ + ('outputnode.itk_fbold_to_t1', 'fboldref2anat_xfm'), + ('outputnode.flipped_boldref', 'flipped_boldref')]), + (bold_reg_wf, summary, [('outputnode.fallback', 'fallback'), + ('outputnode.flip_info', 'flip_info')]), ]) # fmt:on else: diff --git a/fmriprep/workflows/bold/outputs.py b/fmriprep/workflows/bold/outputs.py index 4f3450b1..17f1e0db 100644 --- a/fmriprep/workflows/bold/outputs.py +++ b/fmriprep/workflows/bold/outputs.py @@ -203,6 +203,9 @@ def init_func_fit_reports_wf( 't1w_dseg', 'fieldmap', 'fmap_ref', + # LR flip check + 'fboldref2anat_xfm', + 'flipped_boldref', # May be missing 'subject_id', 'subjects_dir', @@ -267,6 +270,29 @@ def init_func_fit_reports_wf( mem_gb=1, ) + # LR flip check + t1w_flipped_boldref = pe.Node( + ApplyTransforms( + dimension=3, + default_value=0, + float=True, + invert_transform_flags=[True], + interpolation='LanczosWindowedSinc', + ), + name='t1w_flipped_boldref', + mem_gb=1, + ) + flipped_boldref_wm = pe.Node( + ApplyTransforms( + dimension=3, + default_value=0, + invert_transform_flags=[True], + interpolation='NearestNeighbor', + ), + name='boldref_wm', + mem_gb=1, + ) + # fmt:off workflow.connect([ (inputnode, ds_summary, [ @@ -288,6 +314,18 @@ def init_func_fit_reports_wf( ('boldref2anat_xfm', 'transforms'), ]), (t1w_wm, boldref_wm, [('out', 'input_image')]), + # LR flip check + (inputnode, t1w_flipped_boldref, [ + ('t1w_preproc', 'input_image'), + ('flipped_boldref', 'reference_image'), + ('fboldref2anat_xfm', 'transforms'), + ]), + (inputnode, flipped_boldref_wm, [ + ('flipped_boldref', 'reference_image'), + ('fboldref2anat_xfm', 'transforms'), + ]), + (t1w_wm, flipped_boldref_wm, [('out', 'input_image')]), + ]) # fmt:on @@ -409,6 +447,27 @@ def init_func_fit_reports_wf( name='ds_epi_t1_report', ) + flipped_epi_t1_report = pe.Node( + SimpleBeforeAfter( + before_label='T1w', + after_label='EPI', + dismiss_affine=True, + ), + name='flipped_epi_t1_report', + mem_gb=0.1, + ) + + ds_flipped_epi_t1_report = pe.Node( + DerivativesDataSink( + base_directory=output_dir, + desc='flippedcoreg', + suffix='bold', + datatype='figures', + dismiss_entities=dismiss_echo(), + ), + name='ds_flipped_epi_t1_report', + ) + # fmt:off workflow.connect([ (inputnode, epi_t1_report, [('coreg_boldref', 'after')]), @@ -416,6 +475,11 @@ def init_func_fit_reports_wf( (boldref_wm, epi_t1_report, [('output_image', 'wm_seg')]), (inputnode, ds_epi_t1_report, [('source_file', 'source_file')]), (epi_t1_report, ds_epi_t1_report, [('out_report', 'in_file')]), + (inputnode, flipped_epi_t1_report, [('flipped_boldref', 'after')]), + (t1w_flipped_boldref, flipped_epi_t1_report, [('output_image', 'before')]), + (flipped_boldref_wm, flipped_epi_t1_report, [('output_image', 'wm_seg')]), + (inputnode, ds_flipped_epi_t1_report, [('source_file', 'source_file')]), + (flipped_epi_t1_report, ds_flipped_epi_t1_report, [('out_report', 'in_file')]), ]) # fmt:on diff --git a/fmriprep/workflows/bold/registration.py b/fmriprep/workflows/bold/registration.py index 3f7b5e1d..a53070ce 100644 --- a/fmriprep/workflows/bold/registration.py +++ b/fmriprep/workflows/bold/registration.py @@ -127,6 +127,8 @@ def init_bold_reg_wf( Affine transform from T1 space to BOLD space (ITK format) fallback Boolean indicating whether BBR was rejected (mri_coreg registration returned) + flip_info + Information regarding whether a left-right flip was detected See Also -------- @@ -153,7 +155,16 @@ def init_bold_reg_wf( ) outputnode = pe.Node( - niu.IdentityInterface(fields=['itk_bold_to_t1', 'itk_t1_to_bold', 'fallback']), + niu.IdentityInterface( + fields=[ + 'itk_bold_to_t1', + 'itk_t1_to_bold', + 'fallback', + 'flip_info', + 'itk_fbold_to_t1', + 'flipped_boldref', + ] + ), name='outputnode', ) @@ -187,6 +198,9 @@ def init_bold_reg_wf( ('outputnode.itk_bold_to_t1', 'itk_bold_to_t1'), ('outputnode.itk_t1_to_bold', 'itk_t1_to_bold'), ('outputnode.fallback', 'fallback'), + ('outputnode.flip_info', 'flip_info'), + ('outputnode.itk_fbold_to_t1', 'itk_fbold_to_t1'), + ('outputnode.flipped_boldref', 'flipped_boldref'), ]), ]) # fmt:skip @@ -267,13 +281,16 @@ def init_bbreg_wf( Affine transform from T1 space to BOLD space (ITK format) fallback Boolean indicating whether BBR was rejected (mri_coreg registration returned) + flip_info + Information regarding whether a left-right flip was detected """ from nipype.interfaces.freesurfer import BBRegister from niworkflows.engine.workflows import LiterateWorkflow as Workflow - from niworkflows.interfaces.nitransforms import ConcatenateXFMs + from niworkflows.interfaces.morphology import AxisFlip from fmriprep.interfaces.patches import FreeSurferSource, MRICoreg + from fmriprep.interfaces.reports import CheckFlip workflow = Workflow(name=name) workflow.__desc__ = """\ @@ -308,7 +325,16 @@ def init_bbreg_wf( name='inputnode', ) outputnode = pe.Node( - niu.IdentityInterface(['itk_bold_to_t1', 'itk_t1_to_bold', 'fallback']), + niu.IdentityInterface( + [ + 'itk_bold_to_t1', + 'itk_t1_to_bold', + 'fallback', + 'flip_info', + 'itk_fbold_to_t1', + 'flipped_boldref', + ] + ), name='outputnode', ) @@ -347,24 +373,23 @@ def init_bbreg_wf( if bold2anat_init == 'header': bbregister.inputs.init = 'header' - transforms = pe.Node(niu.Merge(2), run_without_submitting=True, name='transforms') - # In cases where Merge(2) only has `in1` or `in2` defined - # output list will just contain a single element - select_transform = pe.Node( - niu.Select(index=0), run_without_submitting=True, name='select_transform' + lr_flip = pe.Node(AxisFlip(axis=0), name='flip') + bbregister_flipped = pe.Node( + BBRegister( + dof=bold2anat_dof, + contrast_type='t2', + out_lta_file=True, + ), + name='bbregister', + mem_gb=12, ) - merge_ltas = pe.Node(niu.Merge(2), name='merge_ltas', run_without_submitting=True) - concat_xfm = pe.Node(ConcatenateXFMs(inverse=True), name='concat_xfm') + if bold2anat_init == 'header': + bbregister_flipped.inputs.init = 'header' - workflow.connect([ - (inputnode, merge_ltas, [('fsnative2t1w_xfm', 'in2')]), - # Wire up the co-registration alternatives - (transforms, select_transform, [('out', 'inlist')]), - (select_transform, merge_ltas, [('out', 'in1')]), - (merge_ltas, concat_xfm, [('out', 'in_xfms')]), - (concat_xfm, outputnode, [('out_xfm', 'itk_bold_to_t1')]), - (concat_xfm, outputnode, [('out_inv', 'itk_t1_to_bold')]), - ]) # fmt:skip + check_flip = pe.Node(CheckFlip(), name='check_flip') + + transform_wf = _transform_handling_wf(use_bbr=use_bbr) + transform_flipped_wf = _transform_handling_wf(use_bbr=use_bbr) # Do not initialize with header, use mri_coreg if bold2anat_init != 'header': @@ -372,7 +397,8 @@ def init_bbreg_wf( (inputnode, mri_coreg, [('subjects_dir', 'subjects_dir'), ('subject_id', 'subject_id'), ('in_file', 'source_file')]), - (mri_coreg, transforms, [('out_lta_file', 'in2')]), + (mri_coreg, transform_wf, [('out_lta_file', 'inputnode.in2')]), + (mri_coreg, transform_flipped_wf, [('out_lta_file', 'inputnode.in2')]), ]) # fmt:skip if use_t2w: @@ -390,15 +416,86 @@ def init_bbreg_wf( # Otherwise bbregister will also be used workflow.connect(mri_coreg, 'out_lta_file', bbregister, 'init_reg_file') + workflow.connect(mri_coreg, 'out_lta_file', bbregister_flipped, 'init_reg_file') # Use bbregister workflow.connect([ (inputnode, bbregister, [('subjects_dir', 'subjects_dir'), ('subject_id', 'subject_id'), ('in_file', 'source_file')]), - (bbregister, transforms, [('out_lta_file', 'in1')]), + (inputnode, bbregister_flipped, [('subjects_dir', 'subjects_dir'), + ('subject_id', 'subject_id')]), + (inputnode, lr_flip), [('in_file', 'in_file')], + (lr_flip, bbregister_flipped), [('out_file', 'source_file')], + (lr_flip, outputnode), [('out_file', 'flipped_boldref')], + (bbregister, check_flip), [('min_cost_file', 'cost_original')], + (bbregister, transform_wf, [('out_lta_file', 'inputnode.in1')]), + (bbregister_flipped, check_flip), [('min_cost_file', 'cost_flipped')], + (bbregister_flipped, transform_flipped_wf, [('out_lta_file', 'inputnode.in1')]), + (check_flip, outputnode, [('flip_info', 'flip_info')]), + (transform_wf, outputnode, [('outputnode.itk_bold_to_t1', 'itk_bold_to_t1'), + ('outputnode.itk_t1_to_bold', 'itk_t1_to_bold'), + ('outputnode.fallback', 'fallback')]), + (transform_flipped_wf, outputnode, [('outputnode.itk_bold_to_t1', 'itk_fbold_to_t1')]), ]) # fmt:skip + return workflow + + +def _transform_handling_wf(use_bbr: bool, name: str = 'transform_handling_wf'): + """ + Wire up the co-registration alternatives + + Parameters + ---------- + use_bbr : :obj:`bool` or None + Enable/disable boundary-based registration refinement. + If ``None``, test BBR result for distortion before accepting. + """ + from niworkflows.engine.workflows import LiterateWorkflow as Workflow + from niworkflows.interfaces.nitransforms import ConcatenateXFMs + + workflow = Workflow(name=name) + inputnode = pe.Node( + niu.IdentityInterface( + [ + 'in1', + 'in2', + ] + ), + name='inputnode', + ) + outputnode = pe.Node( + niu.IdentityInterface( + [ + 'itk_bold_to_t1', + 'itk_t1_to_bold', + 'fallback', + ] + ), + name='outputnode', + ) + + transforms = pe.Node(niu.Merge(2), run_without_submitting=True, name='transforms') + # In cases where Merge(2) only has `in1` or `in2` defined + # output list will just contain a single element + select_transform = pe.Node( + niu.Select(index=0), run_without_submitting=True, name='select_transform' + ) + merge_ltas = pe.Node(niu.Merge(2), name='merge_ltas', run_without_submitting=True) + concat_xfm = pe.Node(ConcatenateXFMs(inverse=True), name='concat_xfm') + + workflow.connect( + [ + (inputnode, transforms, [('in1', 'in1'), ('in2', 'in2')]), + (transforms, select_transform, [('out', 'inlist')]), + (select_transform, merge_ltas, [('out', 'in1')]), + (merge_ltas, concat_xfm, [('out', 'in_xfms')]), + (concat_xfm, outputnode, [('out_xfm', 'itk_bold_to_t1')]), + (concat_xfm, outputnode, [('out_inv', 'itk_t1_to_bold')]), + ] + ) + # Short-circuit workflow building, use boundary-based registration if use_bbr is True: outputnode.inputs.fallback = False @@ -531,7 +628,15 @@ def init_fsl_bbr_wf( name='inputnode', ) outputnode = pe.Node( - niu.IdentityInterface(['itk_bold_to_t1', 'itk_t1_to_bold', 'fallback']), + niu.IdentityInterface( + [ + 'itk_bold_to_t1', + 'itk_t1_to_bold', + 'fallback', + 'flip_info', + 'itk_fbold_to_t1', + ] + ), name='outputnode', ) @@ -576,6 +681,18 @@ def init_fsl_bbr_wf( name='fsl2itk_inv', mem_gb=DEFAULT_MEMORY_MIN_GB, ) + + def flip_detection_not_implemented_yet(): + return { + 'warning': 'Flip detection not implemented yet for alignment with FSL flirt', + 'cost_original': '', + 'cost_flipped': '', + } + + check_flip = pe.Node( + niu.Function(function=flip_detection_not_implemented_yet, output_names=['flip_info']), + name='check_flip', + ) # fmt:off workflow.connect([ (inputnode, mask_t1w_brain, [('t1w_preproc', 'in_file'), @@ -590,6 +707,7 @@ def init_fsl_bbr_wf( (invt_bbr, fsl2itk_inv, [('out_file', 'transform_file')]), (fsl2itk_fwd, outputnode, [('itk_transform', 'itk_bold_to_t1')]), (fsl2itk_inv, outputnode, [('itk_transform', 'itk_t1_to_bold')]), + (check_flip, outputnode, [('flip_info', 'flip_info')]), ]) # fmt:on diff --git a/pyproject.toml b/pyproject.toml index f707950a..2b3a3235 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -25,7 +25,7 @@ dependencies = [ "nireports @ git+https://github.com/nipreps/nireports.git@main", "nitime", "nitransforms >= 21.0.0", - "niworkflows @ git+https://github.com/nipreps/niworkflows.git@master", + "niworkflows @ git+https://github.com/celprov/niworkflows.git@enh/flip_image", "numpy >= 1.22", "packaging", "pandas",