[ENH] Template-based masking of EPI boldrefs

fmriprep/interfaces/__init__.py

@@ -5,7 +5,7 @@
     ReadSidecarJSON, DerivativesDataSink, BIDSDataGrabber, BIDSFreeSurferDir, BIDSInfo
 )
 from .images import (
-    IntraModalMerge, ValidateImage, TemplateDimensions, Conform
+    IntraModalMerge, ValidateImage, TemplateDimensions, Conform, MatchHeader
 )
 from .freesurfer import (
     StructuralReference, MakeMidthickness, FSInjectBrainExtracted,

fmriprep/interfaces/images.py

@@ -479,6 +479,51 @@ def _run_interface(self, runtime):
         return runtime
 
 
+class MatchHeaderInputSpec(BaseInterfaceInputSpec):
+    reference = File(exists=True, mandatory=True,
+                     desc='NIfTI file with reference header')
+    in_file = File(exists=True, mandatory=True,
+                   desc='NIfTI file which header will be checked')
+
+
+class MatchHeaderOutputSpec(TraitedSpec):
+    out_file = File(exists=True, desc='NIfTI file with fixed header')
+
+
+class MatchHeader(SimpleInterface):
+    input_spec = MatchHeaderInputSpec
+    output_spec = MatchHeaderOutputSpec
+
+    def _run_interface(self, runtime):
+        refhdr = nb.load(self.inputs.reference).header.copy()
+        imgnii = nb.load(self.inputs.in_file)
+        imghdr = imgnii.header.copy()
+
+        imghdr['dim_info'] = refhdr['dim_info']  # dim_info is lost sometimes
+
+        # Set qform
+        qform, qcode = refhdr.get_qform(coded=True)
+        if not np.allclose(qform, imghdr.get_qform()):
+            LOGGER.warning(
+                'q-forms of reference and mask are substantially different')
+        imghdr.set_qform(qform, int(qcode))
+
+        # Set sform
+        sform, scode = refhdr.get_sform(coded=True)
+        if not np.allclose(sform, imghdr.get_sform()):
+            LOGGER.warning(
+                's-forms of reference and mask are substantially different')
+        imghdr.set_sform(sform, int(scode))
+
+        out_file = fname_presuffix(self.inputs.in_file, suffix='_hdr',
+                                   newpath=runtime.cwd)
+
+        imgnii.__class__(imgnii.get_data(), imgnii.affine, imghdr).to_filename(
+            out_file)
+        self._results['out_file'] = out_file
+        return runtime
+
+
 def reorient(in_file, newpath=None):
     """Reorient Nifti files to RAS"""
     out_file = fname_presuffix(in_file, suffix='_ras', newpath=newpath)

fmriprep/workflows/bold/util.py

@@ -26,7 +26,7 @@
 from niworkflows.interfaces.utils import CopyXForm
 
 from ...engine import Workflow
-from ...interfaces import ValidateImage
+from ...interfaces import ValidateImage, MatchHeader
 
 
 DEFAULT_MEMORY_MIN_GB = 0.01
@@ -239,6 +239,10 @@ def init_enhance_and_skullstrip_bold_wf(
         operation='max', kernel_shape='sphere', kernel_size=3.0,
         internal_datatype='char'), name='pre_mask_dilate')
 
+    # Ensure mask's header matches reference's
+    check_hdr = pe.Node(MatchHeader(), name='check_hdr',
+                        run_without_submitting=True)
+
     # Run N4 normally, force num_threads=1 for stability (images are small, no need for >1)
     n4_correct = pe.Node(ants.N4BiasFieldCorrection(dimension=3, copy_header=True),
                          name='n4_correct', n_procs=1)
@@ -322,7 +326,9 @@ def init_enhance_and_skullstrip_bold_wf(
         ])
 
     workflow.connect([
-        (pre_dilate, n4_correct, [('out_file', 'mask_image')]),
+        (inputnode, check_hdr, [('in_file', 'reference')]),
+        (pre_dilate, check_hdr, [('out_file', 'in_file')]),
+        (check_hdr, n4_correct, [('out_file', 'mask_image')]),
         (inputnode, n4_correct, [('in_file', 'input_image')]),
         (inputnode, fixhdr_unifize, [('in_file', 'hdr_file')]),
         (inputnode, fixhdr_skullstrip2, [('in_file', 'hdr_file')]),