Merge pull request #154 from oesteban/enh/topup-coeff-conformation

ENH: Massage TOPUP's fieldcoeff files to be compatible with ours

Author: oesteban

sdcflows/conftest.py

@@ -1,4 +1,4 @@
-"""py.test configuration"""
+"""py.test configuration."""
 import os
 from pathlib import Path
 import numpy

sdcflows/interfaces/bspline.py

@@ -317,6 +317,60 @@ def _run_interface(self, runtime):
         return runtime
 
 
+class _TOPUPCoeffReorientInputSpec(BaseInterfaceInputSpec):
+    in_coeff = InputMultiObject(
+        File(exist=True), mandatory=True, desc="input coefficients file(s) from TOPUP"
+    )
+    fmap_ref = File(exists=True, mandatory=True, desc="the fieldmap reference")
+
+
+class _TOPUPCoeffReorientOutputSpec(TraitedSpec):
+    out_coeff = OutputMultiObject(File(exists=True), desc="patched coefficients")
+
+
+class TOPUPCoeffReorient(SimpleInterface):
+    """
+    Revise the orientation of TOPUP-generated B-Spline coefficients.
+
+    TOUP-generated "fieldcoeff" files are just B-Spline fields, where the shape
+    of the field is fixated to be a decimated grid of the original image by an
+    integer factor and added 3 pixels on each dimension.
+    This is one root reason why TOPUP errors (FSL 6) or segfaults (FSL5), when the
+    input image has odd number of voxels along one or more directions.
+
+    These "fieldcoeff" are fixated to be zero-centered, and have "plumb" orientation
+    (as in, aligned with cardinal/imaging axes).
+    The q-form of these NIfTI files is always diagonal, with the decimation factors
+    set on the diagonal (and hence, the voxel zooms).
+    The origin of the q-form is set to the reference image's shape.
+
+    This interface modifies these coefficient files to be fully-fledged NIfTI images
+    aligned with the reference image.
+    Therefore, the s-form header of the coefficients file is updated to match that
+    of the reference file.
+    The s-form header is used because the imaging axes may be oblique.
+
+    The q-form retains the original header and is marked with code 0.
+
+    """
+
+    input_spec = _TOPUPCoeffReorientInputSpec
+    output_spec = _TOPUPCoeffReorientOutputSpec
+
+    def _run_interface(self, runtime):
+        self._results["out_coeff"] = [
+            str(
+                _fix_topup_fieldcoeff(
+                    in_coeff,
+                    self.inputs.fmap_ref,
+                    fname_presuffix(in_coeff, suffix="_fixed", newpath=runtime.cwd),
+                )
+            )
+            for in_coeff in self.inputs.in_coeff
+        ]
+        return runtime
+
+
 def bspline_grid(img, control_zooms_mm=DEFAULT_ZOOMS_MM):
     """Create a :obj:`~nibabel.nifti1.Nifti1Image` embedding the location of control points."""
     if isinstance(img, (str, Path)):
@@ -437,3 +491,36 @@ def _move_coeff(in_coeff, fmap_ref, transform):
         img.__class__(img.dataobj, newaff, img.header).to_filename(out[-1])
 
     return out
+
+
+def _fix_topup_fieldcoeff(in_coeff, fmap_ref, out_file=None):
+    """Read in a coefficients file generated by TOPUP and fix x-form headers."""
+    from pathlib import Path
+    import numpy as np
+    import nibabel as nb
+
+    if out_file is None:
+        out_file = Path("coefficients.nii.gz").absolute()
+
+    coeffnii = nb.load(in_coeff)
+    refnii = nb.load(fmap_ref)
+
+    coeff_shape = np.array(coeffnii.shape[:3])
+    ref_shape = np.array(refnii.shape[:3])
+    factors = coeffnii.header.get_zooms()[:3]
+    if not np.all(coeff_shape == ref_shape // factors + 3):
+        raise ValueError(
+            f"Shape of coefficients file {coeff_shape} does not meet the "
+            f"expectation given the reference's shape {ref_shape}."
+        )
+    newaff = np.eye(4)
+    newaff[:3, :3] = refnii.affine[:3, :3] * factors
+    c_ref = nb.affines.apply_affine(refnii.affine, 0.5 * (ref_shape - 1))
+    c_coeff = nb.affines.apply_affine(newaff, 0.5 * (coeff_shape - 1))
+    newaff[:3, 3] = c_ref - c_coeff
+    header = coeffnii.header.copy()
+    coeffnii.header.set_qform(coeffnii.header.get_qform(coded=False), code=0)
+    coeffnii.header.set_sform(newaff, code=1)
+
+    coeffnii.__class__(coeffnii.dataobj, newaff, header).to_filename(out_file)
+    return out_file

sdcflows/interfaces/tests/test_bspline.py

@@ -4,7 +4,13 @@
 import nibabel as nb
 import pytest
 
-from ..bspline import bspline_grid, Coefficients2Warp, BSplineApprox
+from ..bspline import (
+    bspline_grid,
+    Coefficients2Warp,
+    BSplineApprox,
+    TOPUPCoeffReorient,
+    _fix_topup_fieldcoeff,
+)
 
 
 @pytest.mark.parametrize("testnum", range(100))
@@ -52,3 +58,29 @@ def test_bsplines(tmp_path, testnum):
 
     # Absolute error of the interpolated field is always below 2 Hz
     assert np.all(np.abs(nb.load(test2.outputs.out_error).get_fdata()) < 2)
+
+
+def test_topup_coeffs(tmpdir, testdata_dir):
+    """Check the revision of TOPUP headers."""
+    tmpdir.chdir()
+    result = TOPUPCoeffReorient(
+        in_coeff=str(testdata_dir / "topup-coeff.nii.gz"),
+        fmap_ref=str(testdata_dir / "epi.nii.gz"),
+    ).run()
+
+    nii = nb.load(result.outputs.out_coeff)
+    ctrl = nb.load(testdata_dir / "topup-coeff-fixed.nii.gz")
+    assert np.allclose(nii.affine, ctrl.affine)
+
+    nb.Nifti1Image(nii.get_fdata()[:-1, :-1, :-1], nii.affine, nii.header).to_filename(
+        "failing.nii.gz"
+    )
+
+    with pytest.raises(ValueError):
+        TOPUPCoeffReorient(
+            in_coeff="failing.nii.gz", fmap_ref=str(testdata_dir / "epi.nii.gz"),
+        ).run()
+
+    # Test automatic output file name generation, just for coverage
+    with pytest.raises(ValueError):
+        _fix_topup_fieldcoeff("failing.nii.gz", str(testdata_dir / "epi.nii.gz"))

sdcflows/tests/data/topup-coeff-fixed.nii.gz

@@ -69,6 +69,7 @@ def init_topup_wf(omp_nthreads=1, debug=False, name="pepolar_estimate_wf"):
 
     from ...interfaces.epi import GetReadoutTime, EPIMask
     from ...interfaces.utils import Flatten
+    from ...interfaces.bspline import TOPUPCoeffReorient
 
     workflow = Workflow(name=name)
     workflow.__desc__ = f"""\
@@ -113,6 +114,11 @@ def init_topup_wf(omp_nthreads=1, debug=False, name="pepolar_estimate_wf"):
     )
 
     epi_mask = pe.Node(EPIMask(), name="epi_mask")
+
+    fix_coeff = pe.Node(
+        TOPUPCoeffReorient(), name="fix_coeff", run_without_submitting=True
+    )
+
     # fmt: off
     workflow.connect([
         (inputnode, flatten, [("in_data", "in_data"),
@@ -124,14 +130,16 @@ def init_topup_wf(omp_nthreads=1, debug=False, name="pepolar_estimate_wf"):
         (readout_time, topup, [("readout_time", "readout_times")]),
         (concat_blips, topup, [("out_file", "in_file")]),
         (topup, merge_corrected, [("out_corrected", "in_files")]),
+        (topup, fix_coeff, [("out_fieldcoef", "in_coeff"),
+                            ("out_corrected", "fmap_ref")]),
         (topup, outputnode, [("out_field", "fmap"),
-                             ("out_fieldcoef", "fmap_coeff"),
                              ("out_jacs", "jacobians"),
                              ("out_mats", "xfms"),
                              ("out_warps", "out_warps")]),
         (merge_corrected, epi_mask, [("out_avg", "in_file")]),
         (merge_corrected, outputnode, [("out_avg", "fmap_ref")]),
         (epi_mask, outputnode, [("out_file", "fmap_mask")]),
+        (fix_coeff, outputnode, [("out_coeff", "fmap_coeff")]),
     ])
     # fmt: on