FIX: Make images "plumb" before running ANTs-SyN (and roll-back afterward)

Because ITK handles registration in "physical" coordinates, when the fixed
image's affine is not aligned with the array's axis (i.e., is not plumb or
aligned to the canonical axis) the transform is obtained with alignment to
one axis, but that is not aligned with the data.

To ensure the PE (phase-encoding) direction of the data array and the
corresponding physical axis are aligned, this PR compensates for potental
rotations.

Resolves #12.

Author: oesteban

sdcflows/interfaces/bspline.py

@@ -241,19 +241,13 @@ def _run_interface(self, runtime):
         # Calculate the physical coordinates of target grid
         targetnii = nb.load(self.inputs.in_target)
         allmask = np.ones_like(targetnii.dataobj, dtype="uint8")
-        # targetaff = targetnii.affine
-        # voxels = np.argwhere(allmask == 1).astype("float32")
-        # points = apply_affine(targetaff.astype("float32"), voxels)
 
         weights = []
         coeffs = []
 
         for cname in self.inputs.in_coeff:
             coeff_nii = nb.load(cname)
             wmat = grid_bspline_weights(targetnii, coeff_nii)
-            # wmat = bspline_weights(
-            #     points, coeff_nii, mem_percent=0.1 if self.inputs.low_mem else None,
-            # )
             weights.append(wmat)
             coeffs.append(coeff_nii.get_fdata(dtype="float32").reshape(-1))
 
@@ -284,6 +278,7 @@ def _run_interface(self, runtime):
         field[..., phaseEncDim] = data.reshape(-1)
         aff = targetnii.affine.copy()
         aff[:3, 3] = 0.0
+        # Multiplying by the affine implicitly applies the voxel size to the shift map
         field = nb.affines.apply_affine(aff, field).reshape(fieldshape)
         warpnii = targetnii.__class__(
             field[:, :, :, np.newaxis, :].astype("float32"), targetnii.affine, None

sdcflows/interfaces/tests/test_utils.py

@@ -3,7 +3,7 @@
 import numpy as np
 import nibabel as nb
 
-from ..utils import Flatten, ConvertWarp
+from ..utils import Flatten, ConvertWarp, Deoblique, Reoblique
 
 
 def test_Flatten(tmpdir):
@@ -43,3 +43,47 @@ def test_ConvertWarp(tmpdir, shape):
     assert nii.header.get_data_dtype() == np.float32
     assert nii.header.get_intent() == ("vector", (), "")
     assert nii.shape == (10, 10, 10, 1, 3)
+
+
+@pytest.mark.parametrize(
+    "angles,oblique",
+    [
+        ((0, 0, 0), False),
+        ((0.9, 0.001, 0.001), True),
+        ((0, 0, 2 * np.pi), False),
+    ],
+)
+def test_Xeoblique(tmpdir, angles, oblique):
+    """Exercise De/Reoblique interfaces."""
+
+    tmpdir.chdir()
+
+    affine = nb.affines.from_matvec(nb.eulerangles.euler2mat(*angles))
+    nb.Nifti1Image(np.zeros((10, 10, 10), dtype="uint8"), affine, None).to_filename(
+        "epi.nii.gz"
+    )
+
+    result = (
+        Deoblique(
+            in_epi="epi.nii.gz",
+            mask_epi="epi.nii.gz",
+            in_anat="epi.nii.gz",
+            mask_anat="epi.nii.gz",
+        )
+        .run()
+        .outputs
+    )
+
+    assert np.allclose(nb.load(result.out_epi).affine, affine) is not oblique
+
+    reoblique = (
+        Reoblique(
+            in_plumb=result.out_epi,
+            in_field=result.out_epi,
+            in_epi="epi.nii.gz",
+        )
+        .run()
+        .outputs
+    )
+
+    assert np.allclose(nb.load(reoblique.out_epi).affine, affine)

sdcflows/interfaces/utils.py

@@ -11,20 +11,30 @@
     OutputMultiObject,
 )
 
+OBLIQUE_THRESHOLD_DEG = 0.5
+
 
 class _FlattenInputSpec(BaseInterfaceInputSpec):
     in_data = InputMultiObject(
-        File(exists=True), mandatory=True, desc="list of input data",
+        File(exists=True),
+        mandatory=True,
+        desc="list of input data",
     )
     in_meta = InputMultiObject(
-        traits.DictStrAny, mandatory=True, desc="list of metadata",
+        traits.DictStrAny,
+        mandatory=True,
+        desc="list of metadata",
     )
     max_trs = traits.Int(50, usedefault=True, desc="only pick first TRs")
 
 
 class _FlattenOutputSpec(TraitedSpec):
     out_list = OutputMultiObject(
-        traits.Tuple(File(exists=True), traits.DictStrAny,), desc="list of output files"
+        traits.Tuple(
+            File(exists=True),
+            traits.DictStrAny,
+        ),
+        desc="list of output files",
     )
     out_data = OutputMultiObject(File(exists=True))
     out_meta = OutputMultiObject(traits.DictStrAny)
@@ -71,6 +81,97 @@ def _run_interface(self, runtime):
         return runtime
 
 
+class _DeobliqueInputSpec(BaseInterfaceInputSpec):
+    in_epi = File(
+        exists=True, mandatory=True, desc="the EPI dataset potentially oblique"
+    )
+    mask_epi = File(
+        exists=True,
+        mandatory=True,
+        desc="a binary mask corresponding to the EPI dataset",
+    )
+    in_anat = File(
+        exists=True,
+        mandatory=True,
+        desc="the corresponging anatomical dataset potentially oblique",
+    )
+    mask_anat = File(
+        exists=True,
+        mandatory=True,
+        desc="a binary mask corresponding to the anatomical dataset",
+    )
+
+
+class _DeobliqueOutputSpec(TraitedSpec):
+    out_epi = File(exists=True, desc="the plumb, EPI dataset")
+    out_anat = File(exists=True, desc="the plumb, anatomical dataset")
+    mask_epi = File(
+        exists=True,
+        mandatory=True,
+        desc="a binary mask corresponding to the EPI dataset",
+    )
+    mask_anat = File(
+        exists=True,
+        mandatory=True,
+        desc="a binary mask corresponding to the anatomical dataset",
+    )
+
+
+class Deoblique(SimpleInterface):
+    """Make a dataset plumb."""
+
+    input_spec = _DeobliqueInputSpec
+    output_spec = _DeobliqueOutputSpec
+
+    def _run_interface(self, runtime):
+        (
+            self._results["out_epi"],
+            self._results["out_anat"],
+            self._results["mask_epi"],
+            self._results["mask_anat"],
+        ) = _deoblique(
+            self.inputs.in_epi,
+            self.inputs.in_anat,
+            self.inputs.mask_epi,
+            self.inputs.mask_anat,
+            newpath=runtime.cwd,
+        )
+        return runtime
+
+
+class _ReobliqueInputSpec(BaseInterfaceInputSpec):
+    in_plumb = File(exists=True, mandatory=True, desc="the plumb EPI image")
+    in_field = File(
+        exists=True,
+        mandatory=True,
+        desc="the plumb field map, extracted from the displacements field estimated by SyN",
+    )
+    in_epi = File(
+        exists=True, mandatory=True, desc="the original, potentially oblique EPI image"
+    )
+
+
+class _ReobliqueOutputSpec(TraitedSpec):
+    out_epi = File(exists=True, desc="the reoblique'd EPI image")
+    out_field = File(exists=True, desc="the reoblique'd EPI image")
+
+
+class Reoblique(SimpleInterface):
+    """Make a dataset plumb."""
+
+    input_spec = _ReobliqueInputSpec
+    output_spec = _ReobliqueOutputSpec
+
+    def _run_interface(self, runtime):
+        (self._results["out_epi"], self._results["out_field"],) = _reoblique(
+            self.inputs.in_epi,
+            self.inputs.in_plumb,
+            self.inputs.in_field,
+            newpath=runtime.cwd,
+        )
+        return runtime
+
+
 def _flatten(inlist, max_trs=50, out_dir=None):
     """
     Split the input EPIs and generate a flattened list with corresponding metadata.
@@ -120,3 +221,63 @@ def _qwarp2ants(in_file, newpath=None):
     data = np.squeeze(nii.get_fdata(dtype="float32"))[..., np.newaxis, :]
     nb.Nifti1Image(data, nii.affine, hdr).to_filename(out_file)
     return out_file
+
+
+def _deoblique(in_epi, in_anat, mask_epi, mask_anat, newpath=None):
+    import numpy as np
+    import nibabel as nb
+    from nipype.utils.filemanip import fname_presuffix
+
+    epinii = nb.load(in_epi)
+    if np.all(nb.affines.obliquity(epinii.affine)) < OBLIQUE_THRESHOLD_DEG:
+        return in_epi, in_anat, mask_epi, mask_anat
+
+    newaff = np.eye(4)
+    newaff[:3, :3] = np.diag(np.array(epinii.header.get_zooms()[:3]))
+    newaff[:3, 3] -= (np.array(epinii.shape[:3]) - 1) * 0.5
+
+    hdr = epinii.header.copy()
+    hdr.set_qform(newaff, code=1)
+    hdr.set_sform(newaff, code=1)
+    newepi = epinii.__class__(epinii.dataobj, newaff, hdr)
+    out_epi = fname_presuffix(in_epi, suffix="_plumb", newpath=newpath)
+    newepi.to_filename(out_epi)
+
+    out_files = [out_epi]
+    for fname in (in_anat, mask_epi, mask_anat):
+        nii = nb.load(fname)
+        hdr = nii.header.copy()
+        hdr.set_qform(newaff, code=1)
+        hdr.set_sform(newaff, code=1)
+        newnii = nii.__class__(np.asanyarray(nii.dataobj), newaff, hdr)
+        out = fname_presuffix(fname, suffix="_plumb", newpath=newpath)
+        newnii.to_filename(out)
+        out_files.append(out)
+
+    return tuple(out_files)
+
+
+def _reoblique(in_epi, in_plumb, in_field, newpath=None):
+    import numpy as np
+    import nibabel as nb
+    from nipype.utils.filemanip import fname_presuffix
+
+    epinii = nb.load(in_epi)
+    if np.all(nb.affines.obliquity(epinii.affine)) < OBLIQUE_THRESHOLD_DEG:
+        return in_plumb, in_field
+
+    out_files = [
+        fname_presuffix(f, suffix="_oriented", newpath=newpath)
+        for f in (in_plumb, in_field)
+    ]
+    plumbnii = nb.load(in_plumb)
+    plumbnii.__class__(plumbnii.dataobj, epinii.affine, epinii.header).to_filename(
+        out_files[0]
+    )
+
+    fmapnii = nb.load(in_field)
+    hdr = fmapnii.header.copy()
+    hdr.set_qform(*epinii.header.get_qform(coded=True))
+    hdr.set_sform(*epinii.header.get_sform(coded=True))
+    fmapnii.__class__(fmapnii.dataobj, epinii.affine, hdr).to_filename(out_files[1])
+    return out_files

sdcflows/workflows/fit/syn.py

@@ -135,6 +135,7 @@ def init_syn_sdc_wf(
     from niworkflows.interfaces.nibabel import Binarize
     from ...interfaces.epi import EPIMask
     from ...utils.misc import front as _pop
+    from ...interfaces.utils import Deoblique, Reoblique
     from ...interfaces.bspline import (
         BSplineApprox,
         DEFAULT_LF_ZOOMS_MM,
@@ -188,6 +189,9 @@ def init_syn_sdc_wf(
     )
     atlas_msk = pe.Node(Binarize(thresh_low=atlas_threshold), name="atlas_msk")
 
+    deoblique = pe.Node(Deoblique(), name="deoblique")
+    reoblique = pe.Node(Reoblique(), name="reoblique")
+
     # SyN Registration Core
     syn = pe.Node(
         Registration(from_file=pkgrf("sdcflows", "data/susceptibility_syn.json")),
@@ -221,24 +225,31 @@ def init_syn_sdc_wf(
                                  (("epi_ref", _pop), "ref_file")]),
         (inputnode, anat2epi, [(("epi_ref", _pop), "reference_image"),
                                ("anat2epi_xfm", "transforms")]),
-        (inputnode, syn, [(("epi_ref", _pop), "moving_image"),
-                          ("epi_mask", "moving_image_masks"),
-                          (("epi_ref", _warp_dir), "restrict_deformation")]),
+        (inputnode, deoblique, [(("epi_ref", _pop), "in_epi"),
+                                ("epi_mask", "mask_epi")]),
+        (inputnode, reoblique, [(("epi_ref", _pop), "in_epi")]),
+        (inputnode, syn, [(("epi_ref", _warp_dir), "restrict_deformation")]),
         (inputnode, unwarp_ref, [(("epi_ref", _pop), "reference_image"),
                                  (("epi_ref", _pop), "input_image")]),
         (inputnode, prior2epi, [(("epi_ref", _pop), "reference_image")]),
         (inputnode, extract_field, [("epi_ref", "epi_meta")]),
         (invert_t1w, anat2epi, [("out_file", "input_image")]),
         (transform_list, prior2epi, [("out", "transforms")]),
         (prior2epi, atlas_msk, [("output_image", "in_file")]),
-        (anat2epi, syn, [("output_image", "fixed_image")]),
-        (atlas_msk, syn, [(("out_mask", _fixed_masks_arg), "fixed_image_masks")]),
+        (anat2epi, deoblique, [("output_image", "in_anat")]),
+        (atlas_msk, deoblique, [("out_mask", "mask_anat")]),
+        (deoblique, syn, [("out_epi", "moving_image"),
+                          ("out_anat", "fixed_image"),
+                          ("mask_epi", "moving_image_masks"),
+                          (("mask_anat", _fixed_masks_arg), "fixed_image_masks")]),
         (syn, extract_field, [("forward_transforms", "in_file")]),
         (syn, unwarp_ref, [("forward_transforms", "transforms")]),
-        (unwarp_ref, epi_mask, [("output_image", "in_file")]),
-        (extract_field, bs_filter, [("out", "in_data")]),
+        (unwarp_ref, reoblique, [("output_image", "in_plumb")]),
+        (reoblique, epi_mask, [("out_epi", "in_file")]),
+        (extract_field, reoblique, [("out", "in_field")]),
+        (reoblique, bs_filter, [("out_field", "in_data")]),
         (epi_mask, bs_filter, [("out_file", "in_mask")]),
-        (unwarp_ref, outputnode, [("output_image", "fmap_ref")]),
+        (reoblique, outputnode, [("out_epi", "fmap_ref")]),
         (epi_mask, outputnode, [("out_file", "fmap_mask")]),
         (bs_filter, outputnode, [
             ("out_extrapolated" if not debug else "out_field", "fmap"),
@@ -308,7 +319,8 @@ def _extract_field(in_file, epi_meta):
         np.squeeze(fieldnii.get_fdata(dtype="float32"))[
             ..., "ijk".index(epi_meta[1]["PhaseEncodingDirection"][0])
         ]
-        / trt * (-1.0 if epi_meta[1]["PhaseEncodingDirection"].endswith("-") else 1.0)
+        / trt
+        * (-1.0 if epi_meta[1]["PhaseEncodingDirection"].endswith("-") else 1.0)
     )
     out_file = fname_presuffix(in_file[0], suffix="_fieldmap")
     nii = nb.Nifti1Image(data, fieldnii.affine, None)