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

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 a
factor of 2x 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.

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,54 @@ 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,))
+            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 +485,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/workflows/fit/pepolar.py

@@ -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