Merge pull request #356 from effigies/enh/validate-phase12

ENH: Check registration of magnitude1/magnitude2 images and update headers

Author: oesteban

sdcflows/interfaces/fmap.py

@@ -21,8 +21,10 @@
 #     https://www.nipreps.org/community/licensing/
 #
 """Interfaces to deal with the various types of fieldmap sources."""
+import os
 import numpy as np
 import nibabel as nb
+import nitransforms as nt
 from nipype.utils.filemanip import fname_presuffix
 from nipype import logging
 from nipype.interfaces.base import (
@@ -31,7 +33,10 @@
     File,
     traits,
     SimpleInterface,
+    InputMultiObject,
+    OutputMultiObject,
 )
+from nipype.interfaces import freesurfer as fs
 
 LOGGER = logging.getLogger("nipype.interface")
 
@@ -201,3 +206,185 @@ def _run_interface(self, runtime):
 
         fmapnii.to_filename(self._results["out_file"])
         return runtime
+
+
+class _CheckRegisterInputSpec(TraitedSpec):
+    mag_files = InputMultiObject(
+        File(exists=True),
+        mandatory=True,
+        minlen=1,
+        maxlen=2,
+        desc="Magnitude image(s) to verify registration",
+    )
+    fmap_files = InputMultiObject(
+        File(exists=True),
+        mandatory=True,
+        minlen=1,
+        maxlen=2,
+        desc="Phase(diff) or fieldmap image(s) to update affines",
+    )
+    rot_thresh = traits.Float(
+        0.02, usedefault=True, mandatory=True, desc="rotation threshold in radians"
+    )
+    trans_thresh = traits.Float(
+        1., usedefault=True, mandatory=True, desc="translation threshold in mm"
+    )
+
+
+class _CheckRegisterOutputSpec(TraitedSpec):
+    mag_files = OutputMultiObject(
+        File,
+        desc="Magnitude image(s) verified to be in register, with consistent affines",
+    )
+    fmap_files = OutputMultiObject(
+        File,
+        desc="Fieldmap files with consistent affines",
+    )
+
+
+class CheckRegister(SimpleInterface):
+    """Check registration of one or more images and paired files
+
+    Use cases:
+
+    Phase1 + Phase2:
+
+        >>> fmap_files = ['*_phase1.nii.gz', '*_phase2.nii.gz']
+        >>> mag_files = ['*_magnitude1.nii.gz', '*_magnitude2.nii.gz']
+
+    Phasediff (2 magnitude):
+
+        >>> fmap_files = ['*_phasediff.nii.gz']
+        >>> mag_files = ['*_magnitude1.nii.gz', '*_magnitude2.nii.gz']
+
+    Phasediff (1 magnitude):
+
+        >>> fmap_files = ['*_phasediff.nii.gz']
+        >>> mag_files = ['*_magnitude1.nii.gz']
+
+    Fieldmap (1 magnitude):
+
+        >>> fmap_files = ['*_fieldmap.nii.gz']
+        >>> mag_files = ['*_magnitude.nii.gz']
+
+    In general, we expect all files to have the same affine and shape,
+    and this will be a pass-through interface. If there are two magnitude
+    files where the affine differs, then we will register them and inspect
+    the registration parameters for evidence of significantly different FoV.
+    The default values of 0.02rad and 1mm both correspond to shifts of 1mm
+    (assuming 50mm radius brain).
+
+    This may run ``mri_robust_register``, so should not be run without
+    submitting.
+    """
+
+    input_spec = _CheckRegisterInputSpec
+    output_spec = _CheckRegisterOutputSpec
+
+    def _run_interface(self, runtime):
+        mag_files = self.inputs.mag_files
+        fmap_files = self.inputs.fmap_files
+
+        mag_imgs = [nb.load(fname) for fname in mag_files]
+        fmap_imgs = [nb.load(fname) for fname in fmap_files]
+
+        # Baseline check: paired magnitude/phase maps are basically the same
+        for mag, fmap in zip(mag_imgs, fmap_imgs):
+            msg = _check_gross_geometry(mag, fmap)
+            if msg is not None:
+                LOGGER.critical(msg)
+                raise ValueError(msg)
+
+        # Verify images are in register before conforming affines
+        if len(mag_files) == 2:
+            msg = _check_gross_geometry(mag_imgs[0], mag_imgs[1])
+            if msg is not None:
+                LOGGER.critical(msg)
+                raise ValueError(msg)
+
+            # If affines match, do not attempt to register
+            # Treat this as an assertion by the scanner or curator that they are aligned
+            if not np.allclose(mag_imgs[0].affine, mag_imgs[1].affine):
+                reg = fs.RobustRegister(
+                    target_file=mag_files[0],
+                    source_file=mag_files[1],
+                    auto_sens=True,
+                )
+                result = reg.run()
+                lta = nt.io.lta.FSLinearTransform.from_filename(result.outputs.out_reg_file)
+                mat, vec = nb.affines.to_matvec(lta.to_ras())
+                angles = np.abs(nb.eulerangles.mat2euler(mat))
+                rot_thresh, trans_thresh = (
+                    self.inputs.rot_thresh,
+                    self.inputs.trans_thresh,
+                )
+
+                if np.any(angles > rot_thresh) or np.any(vec > trans_thresh):
+                    LOGGER.critical(
+                        "Magnitude files {mag_files} are not in register with rotation "
+                        f"threshold {self.inputs.rot_thresh} and translation threshold "
+                        f"{self.inputs.trans_thresh}. Please manually verify images "
+                        "are in register and update the image headers before running SDC."
+                    )
+                    raise ValueError(
+                        "Magnitude 1/2 orientation mismatch too big to ignore."
+                    )
+
+        # Probably redundant, but we could hit this error
+        # with phase1/magnitude1 + wonky phase2 + no magnitude2
+        if len(fmap_files) == 2:
+            msg = _check_gross_geometry(fmap_imgs[0], fmap_imgs[1])
+            if msg is not None:
+                LOGGER.critical(msg)
+                raise ValueError(msg)
+
+        # Check/copy affines
+        out_mags = [_conform_img(mag, mag_imgs[0], runtime.cwd) for mag in mag_imgs]
+        out_fmaps = [_conform_img(fmap, mag_imgs[0], runtime.cwd) for fmap in fmap_imgs]
+
+        self._results = {
+            "mag_files": out_mags,
+            "fmap_files": out_fmaps,
+        }
+
+        return runtime
+
+
+def _conform_img(
+    img: nb.spatialimages.SpatialImage,
+    target_img: nb.spatialimages.SpatialImage,
+    cwd: str,
+) -> str:
+    """Return path to image matching target_img geometry
+
+    Copy target_affine to a new image if necessary.
+    """
+    if np.allclose(img.affine, target_img.affine):
+        return img.get_filename()
+
+    basename = os.basename(img.get_filename())
+    out_file = os.path.join(cwd, basename)
+
+    LOGGER.info(f"Copying affine to {basename}")
+    new_img = img.__class__(img.dataobj, target_img.affine, img.header)
+    new_img.to_filename(out_file)
+
+    return out_file
+
+
+def _check_gross_geometry(
+    img1: nb.spatialimages.SpatialImage,
+    img2: nb.spatialimages.SpatialImage,
+):
+    if img1.shape[:3] != img2.shape[:3]:
+        return (
+            "Images have shape mismatch: "
+            f"{img1.get_filename()} {img1.shape}, "
+            f"{img2.get_filename()} {img2.shape}"
+        )
+    if nb.aff2axcodes(img1.affine) != nb.aff2axcodes(img2.affine):
+        return (
+            "Images have orientation mismatch: "
+            f"{img1.get_filename()} {''.join(nb.aff2axcodes(img1.affine))}, "
+            f"{img2.get_filename()} {''.join(nb.aff2axcodes(img2.affine))}"
+        )

sdcflows/workflows/fit/fieldmap.py

@@ -90,6 +90,7 @@ def init_fmap_wf(omp_nthreads=1, sloppy=False, debug=False, mode="phasediff", na
         DEFAULT_HF_ZOOMS_MM,
         DEFAULT_ZOOMS_MM,
     )
+    from ...interfaces.fmap import CheckRegister
 
     workflow = Workflow(name=name)
     inputnode = pe.Node(niu.IdentityInterface(fields=INPUT_FIELDS), name="inputnode")
@@ -98,6 +99,18 @@ def init_fmap_wf(omp_nthreads=1, sloppy=False, debug=False, mode="phasediff", na
         name="outputnode",
     )
 
+    def _unzip(fmap_spec):
+        return fmap_spec
+
+    unzip = pe.MapNode(
+        niu.Function(function=_unzip, output_names=["fmap_file", "meta"]),
+        run_without_submitting=True,
+        iterfield=["fmap_spec"],
+        name="unzip",
+    )
+
+    check_register = pe.Node(CheckRegister(), name='check_register')
+
     magnitude_wf = init_magnitude_wf(omp_nthreads=omp_nthreads)
     bs_filter = pe.Node(BSplineApprox(), name="bs_filter")
     bs_filter.interface._always_run = debug
@@ -108,7 +121,10 @@ def init_fmap_wf(omp_nthreads=1, sloppy=False, debug=False, mode="phasediff", na
 
     # fmt: off
     workflow.connect([
-        (inputnode, magnitude_wf, [("magnitude", "inputnode.magnitude")]),
+        (inputnode, unzip, [("fieldmap", "fmap_spec")]),
+        (inputnode, check_register, [("magnitude", "mag_files")]),
+        (unzip, check_register, [("fmap_file", "fmap_files")]),
+        (check_register, magnitude_wf, [("mag_files", "inputnode.magnitude")]),
         (magnitude_wf, bs_filter, [("outputnode.fmap_mask", "in_mask")]),
         (magnitude_wf, outputnode, [
             ("outputnode.fmap_mask", "fmap_mask"),
@@ -131,7 +147,8 @@ def init_fmap_wf(omp_nthreads=1, sloppy=False, debug=False, mode="phasediff", na
 
         # fmt: off
         workflow.connect([
-            (inputnode, phdiff_wf, [("fieldmap", "inputnode.phase")]),
+            (unzip, phdiff_wf, [("meta", "inputnode.phase_meta")]),
+            (check_register, phdiff_wf, [("fmap_files", "inputnode.phase")]),
             (magnitude_wf, phdiff_wf, [
                 ("outputnode.fmap_ref", "inputnode.magnitude"),
                 ("outputnode.fmap_mask", "inputnode.mask"),
@@ -160,7 +177,7 @@ def init_fmap_wf(omp_nthreads=1, sloppy=False, debug=False, mode="phasediff", na
         # fmt: off
         workflow.connect([
             (inputnode, units, [(("fieldmap", _get_units), "units")]),
-            (inputnode, fmapmrg, [(("fieldmap", _get_file), "in_files")]),
+            (check_register, fmapmrg, [("fmap_files", "in_files")]),
             (fmapmrg, units, [("out_avg", "in_file")]),
             (units, bs_filter, [("out_file", "in_data")]),
         ])
@@ -298,24 +315,15 @@ def init_phdiff_wf(omp_nthreads, debug=False, name="phdiff_wf"):
 """
 
     inputnode = pe.Node(
-        niu.IdentityInterface(fields=["magnitude", "phase", "mask"]), name="inputnode"
+        niu.IdentityInterface(fields=["magnitude", "phase", "phase_meta", "mask"]),
+        name="inputnode",
     )
 
     outputnode = pe.Node(
         niu.IdentityInterface(fields=["fieldmap"]),
         name="outputnode",
     )
 
-    def _split(phase):
-        return phase
-
-    split = pe.MapNode(  # We cannot use an inline connection function with MapNode
-        niu.Function(function=_split, output_names=["map_file", "meta"]),
-        iterfield=["phase"],
-        run_without_submitting=True,
-        name="split",
-    )
-
     # phase diff -> radians
     phmap2rads = pe.MapNode(
         PhaseMap2rads(),
@@ -334,15 +342,14 @@ def _split(phase):
 
     # fmt: off
     workflow.connect([
-        (inputnode, split, [("phase", "phase")]),
+        (inputnode, phmap2rads, [("phase", "in_file")]),
+        (inputnode, calc_phdiff, [("phase_meta", "in_meta")]),
         (inputnode, prelude, [("magnitude", "magnitude_file"),
                               ("mask", "mask_file")]),
-        (split, phmap2rads, [("map_file", "in_file")]),
         (phmap2rads, calc_phdiff, [("out_file", "in_phases")]),
-        (split, calc_phdiff, [("meta", "in_meta")]),
         (calc_phdiff, prelude, [("phase_diff", "phase_file")]),
-        (prelude, compfmap, [("unwrapped_phase_file", "in_file")]),
         (calc_phdiff, compfmap, [("metadata", "metadata")]),
+        (prelude, compfmap, [("unwrapped_phase_file", "in_file")]),
         (compfmap, outputnode, [("out_file", "fieldmap")]),
     ])
     # fmt: on