Merge pull request #165 from oesteban/enh/syn-oblique-patch

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

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)