ENH: Add slice padding to TOPUP

sdcflows/interfaces/tests/test_utils.py

@@ -25,7 +25,7 @@
 import numpy as np
 import nibabel as nb
 
-from ..utils import Flatten, ConvertWarp, Deoblique, Reoblique
+from ..utils import Flatten, ConvertWarp, Deoblique, Reoblique, PadSlices
 
 
 def test_Flatten(tmpdir):
@@ -106,3 +106,36 @@ def test_Xeoblique(tmpdir, angles, oblique):
     )
 
     assert np.allclose(nb.load(reoblique.out_epi).affine, affine)
+
+
+@pytest.mark.parametrize("in_shape,expected_shape,padded", [
+    ((2,2,2), (2,2,2), False),
+    ((2,2,3), (2,2,4), True),
+    ((3,3,2,2), (3,3,2,2), False),
+    ((3,3,3,2), (3,3,4,2), True),
+])
+def test_pad_slices(tmpdir, in_shape, expected_shape, padded):
+    tmpdir.chdir()
+
+    data = np.random.rand(*in_shape)
+    aff = np.eye(4)
+
+    # RAS
+    img = nb.Nifti1Image(data, aff)
+    img.to_filename("epi-ras.nii.gz")
+    res = PadSlices(in_file="epi-ras.nii.gz").run().outputs
+
+    # LPS
+    newaff = aff.copy()
+    newaff[0, 0] *= -1.0
+    newaff[1, 1] *= -1.0
+    newaff[:2, 3] = aff.dot(np.hstack((np.array(img.shape[:3]) - 1, 1.0)))[:2]
+    img2 = nb.Nifti1Image(np.flip(np.flip(data, 0), 1), newaff)
+    img2.to_filename("epi-lps.nii.gz")
+    res2 = PadSlices(in_file="epi-lps.nii.gz").run().outputs
+
+
+    out_ras = nb.load(res.out_file)
+    out_lps = nb.load(res2.out_file)
+    assert out_ras.shape == out_lps.shape == expected_shape
+    assert res.padded == padded

sdcflows/interfaces/utils.py

@@ -267,6 +267,36 @@ class DenoiseImage(_DenoiseImageBase, _CopyHeaderInterface):
     _copy_header_map = {"output_image": "input_image"}
 
 
+class _PadSlicesInputSpec(BaseInterfaceInputSpec):
+    in_file = File(exists=True, mandatory=True, desc="3D or 4D NIfTI image")
+    axis = traits.Int(
+        2,
+        usedefault=True,
+        desc="The axis through which slices are stacked in the input data"
+    )
+
+
+class _PadSlicesOutputSpec(TraitedSpec):
+    out_file = File(exists=True, desc="The output file with even number of slices")
+    padded = traits.Bool(desc="Indicator if the input image was padded")
+
+
+class PadSlices(SimpleInterface):
+    """
+    Check an image for uneven slices, and add an empty slice if necessary
+
+    This intends to avoid TOPUP's segfault without changing the standard configuration
+    """
+    input_spec = _PadSlicesInputSpec
+    output_spec = _PadSlicesOutputSpec
+
+    def _run_interface(self, runtime):
+        self._results["out_file"], self._results["padded"] = _pad_num_slices(
+            self.inputs.in_file, self.inputs.axis, 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.
@@ -384,3 +414,44 @@ def _reoblique(in_epi, in_plumb, in_field, in_mask=None, newpath=None):
         masknii.__class__(masknii.dataobj, epinii.affine, hdr).to_filename(out_files[2])
 
     return out_files
+
+
+def _pad_num_slices(in_file, ax=2, newpath=None):
+    """
+    Ensure the image has even number of slices to avert TOPUP's segfault.
+
+    Check if image has an even number of slices.
+    If it does, return the image unaltered.
+    Otherwise, return the image with an empty slice added.
+
+    Parameters
+    ----------
+    img : :obj:`str` or :py:class:`~nibabel.spatialimages.SpatialImage`
+        3D or 4D NIfTI image
+    ax : :obj:`int`
+        The axis through which slices are stacked in the input data.
+
+    Returns
+    -------
+    file : :obj:`str`
+        The output file with even number of slices
+    padded : :obj:`bool`
+        Indicator if the input image was padded.
+
+    """
+    import nibabel as nb
+    from nipype.utils.filemanip import fname_presuffix
+    import numpy as np
+
+    img = nb.load(in_file)
+    if img.shape[ax] % 2 == 0:
+        return in_file, False
+
+    pwidth = [(0,0)] * len(img.shape)
+    pwidth[ax] = (0, 1)
+    padded = np.pad(img.dataobj, pwidth)
+    hdr = img.header
+    hdr.set_data_shape(padded.shape)
+    out_file = fname_presuffix(in_file, suffix="_padded", newpath=newpath)
+    img.__class__(padded, img.affine, header=hdr).to_filename(out_file)
+    return out_file, True

sdcflows/workflows/fit/pepolar.py

@@ -102,7 +102,7 @@ def init_topup_wf(
 
     from ...utils.misc import front as _front
     from ...interfaces.epi import GetReadoutTime
-    from ...interfaces.utils import Flatten, UniformGrid
+    from ...interfaces.utils import Flatten, UniformGrid, PadSlices
     from ...interfaces.bspline import TOPUPCoeffReorient
     from ..ancillary import init_brainextraction_wf
 
@@ -138,6 +138,8 @@ def init_topup_wf(
         iterfield=["metadata", "in_file"],
         run_without_submitting=True,
     )
+    pad_blip_slices = pe.Node(PadSlices(), name="pad_blip_slices")
+    pad_ref_slices = pe.Node(PadSlices(), name="pad_ref_slices")
 
     topup = pe.Node(
         TOPUP(
@@ -163,7 +165,8 @@ def init_topup_wf(
         (regrid, concat_blips, [("out_data", "in_files")]),
         (readout_time, topup, [("readout_time", "readout_times"),
                                ("pe_dir_fsl", "encoding_direction")]),
-        (regrid, fix_coeff, [("reference", "fmap_ref")]),
+        (regrid, pad_ref_slices, [("reference", "in_file")]),
+        (pad_ref_slices, fix_coeff, [("out_file", "fmap_ref")]),
         (readout_time, fix_coeff, [(("pe_direction", _front), "pe_dir")]),
         (topup, fix_coeff, [("out_fieldcoef", "in_coeff")]),
         (topup, outputnode, [("out_jacs", "jacobians"),
@@ -180,7 +183,8 @@ def init_topup_wf(
     if not debug:
         # fmt: off
         workflow.connect([
-            (concat_blips, topup, [("out_file", "in_file")]),
+            (concat_blips, pad_blip_slices, [("out_file", "in_file")]),
+            (pad_blip_slices, topup, [("out_file", "in_file")]),
             (topup, ref_average, [("out_corrected", "in_file")]),
             (topup, outputnode, [("out_field", "fmap"),
                                  ("out_warps", "out_warps")]),
@@ -204,7 +208,8 @@ def init_topup_wf(
     # fmt:off
     workflow.connect([
         (concat_blips, realign, [("out_file", "in_file")]),
-        (realign, topup, [("out_file", "in_file")]),
+        (realign, pad_blip_slices, [("out_file", "in_file")]),
+        (pad_blip_slices, topup, [("out_file", "in_file")]),
         (fix_coeff, unwarp, [("out_coeff", "in_coeff")]),
         (realign, split_blips, [("out_file", "in_file")]),
         (split_blips, unwarp, [("out_files", "in_target")]),