FIX: Revise the TOPUP workflow and related utilities

docs/requirements.txt

@@ -4,7 +4,7 @@ matplotlib >= 2.2.0
 nibabel
 nipype >= 1.5.1
 traits < 6.4
-niworkflows ~= 1.6.0
+niworkflows ~= 1.6.3
 numpy
 packaging
 pydot >= 1.2.3

min-requirements.txt

@@ -2,10 +2,11 @@
 attrs==20.1.0
 nibabel==3.0.1
 nipype==1.5.1
-niworkflows==1.4.2
+traits<6.4
+niworkflows==1.6.3
 nitransforms==21.0.0
 numpy==1.21.0
-pybids==0.12.1
+pybids==0.15.1
 scikit-image==0.18
 scipy==1.6.0
-templateflow==0.6
+templateflow

requirements.txt

@@ -2,10 +2,11 @@
 attrs>=20.1.0
 nibabel>=3.0.1
 nipype<2.0,>=1.5.1
-niworkflows<1.6,>=1.4.2
+traits<6.4
+niworkflows~=1.6.3
 nitransforms>=21.0.0
 numpy>=1.21.0
-pybids>=0.12.1
+pybids>=0.15.1
 scikit-image>=0.18
 scipy>=1.6.0
-templateflow>=0.6
+templateflow

sdcflows/data/flirtsch/b02b0.cnf

@@ -1,21 +1,21 @@
 # Resolution (knot-spacing) of warps in mm
---warpres=20,16,14,12,10,6,4,4,4
+--warpres=14,12,10,6,4,4,4
 # Subsampling level (a value of 2 indicates that a 2x2x2 neighbourhood is collapsed to 1 voxel)
---subsamp=2,2,2,2,2,1,1,1,1
+--subsamp=2,2,2,1,1,1,1
 # FWHM of gaussian smoothing
---fwhm=8,6,4,3,3,2,1,0,0
+--fwhm=6,4,3,2,1,0,0
 # Maximum number of iterations
---miter=5,5,5,5,5,10,10,20,20
+--miter=5,5,5,10,10,20,20
 # Relative weight of regularisation
---lambda=0.005,0.001,0.0001,0.000015,0.000005,0.0000005,0.00000005,0.0000000005,0.00000000001
+--lambda=0.001,0.0001,0.000005,0.0000005,0.00000005,0.0000000005,0.00000000001
 # If set to 1 lambda is multiplied by the current average squared difference
 --ssqlambda=1
 # Regularisation model
 --regmod=bending_energy
 # If set to 1 movements are estimated along with the field
---estmov=1,1,1,1,1,0,0,0,0
+--estmov=0
 # 0=Levenberg-Marquardt, 1=Scaled Conjugate Gradient
---minmet=0,0,0,0,0,1,1,1,1
+--minmet=0,0,0,1,1,1,1
 # Quadratic or cubic splines
 --splineorder=3
 # Precision for calculation and storage of Hessian

sdcflows/interfaces/bspline.py

@@ -21,6 +21,7 @@
 #     https://www.nipreps.org/community/licensing/
 #
 """Filtering of :math:`B_0` field mappings with B-Splines."""
+from itertools import product
 from pathlib import Path
 import numpy as np
 import nibabel as nb
@@ -369,6 +370,11 @@ class _TOPUPCoeffReorientInputSpec(BaseInterfaceInputSpec):
         File(exist=True), mandatory=True, desc="input coefficients file(s) from TOPUP"
     )
     fmap_ref = File(exists=True, mandatory=True, desc="the fieldmap reference")
+    pe_dir = traits.Enum(
+        *["".join(p) for p in product("ijkxyz", ("", "-"))],
+        mandatory=True,
+        desc="phase encoding direction"
+    )
 
 
 class _TOPUPCoeffReorientOutputSpec(TraitedSpec):
@@ -390,6 +396,8 @@ class TOPUPCoeffReorient(SimpleInterface):
     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.
+    All the director cosines of the output coefficients will be positive.
+    In other words, the output orientation is either RAS, ARS, ASR, SAR, or SRA.
 
     This interface modifies these coefficient files to be fully-fledged NIfTI images
     aligned with the reference image.
@@ -410,6 +418,7 @@ def _run_interface(self, runtime):
                 _fix_topup_fieldcoeff(
                     in_coeff,
                     self.inputs.fmap_ref,
+                    self.inputs.pe_dir,
                     out_file=fname_presuffix(
                         in_coeff, suffix="_fixed", newpath=runtime.cwd
                     ),
@@ -448,36 +457,63 @@ def bspline_grid(img, control_zooms_mm=DEFAULT_ZOOMS_MM):
     return img.__class__(np.zeros(bs_shape, dtype="float32"), bs_affine)
 
 
-def _fix_topup_fieldcoeff(in_coeff, fmap_ref, out_file=None):
+def _fix_topup_fieldcoeff(in_coeff, fmap_ref, pe_dir, 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
+    from sdcflows.utils.tools import ensure_positive_cosines
 
     if out_file is None:
         out_file = Path("coefficients.nii.gz").absolute()
 
-    coeffnii = nb.load(in_coeff)
+    # Load reference image
     refnii = nb.load(fmap_ref)
 
+    # Load coefficients
+    coeffnii = nb.load(in_coeff)
+    # Coefficients - flip LR and overwrite coeffnii variable
+    # Internal data orientation of FSL is LAS, so coefficients will be LR flipped,
+    # and because the affine does not encode orientation (factors instead), this flip
+    # always is implicit.
+    # If the PE direction is x/i, the flip in the axis direction causes that the
+    # fieldmap estimation must also be inverted in direction (multiply by -1.0)
+    reverse_pe = -1.0 if "i" in pe_dir.replace("x", "i") else 1.0
+    lr_axis = "".join(nb.aff2axcodes(coeffnii.affine)).index("R")
+    coeffnii = coeffnii.__class__(
+        reverse_pe * np.flip(np.asanyarray(coeffnii.dataobj), axis=lr_axis),
+        coeffnii.affine,
+        coeffnii.header,
+    )
+    # Ensure reference has positive director cosines
+    refnii, ref_axcodes = ensure_positive_cosines(refnii)
+
+    # Get matrix of B-Spline control knots
     coeff_shape = np.array(coeffnii.shape[:3])
+    # Get factors (w.r.t. reference's pixel sizes) to calculate separation btw control points
     factors = np.array(coeffnii.header.get_zooms()[:3])
+    # Shape checking
     ref_shape = np.array(refnii.shape[:3])
     exp_shape = ref_shape // factors + 3 * (factors > 1)
     if not np.all(coeff_shape == exp_shape):
         raise ValueError(
             f"Shape of coefficients file {coeff_shape} does not meet the "
             f"expected shape {exp_shape} (toupup factors are {factors})."
         )
+
+    # Contextualize the control points in space with a proper NIfTI affine
     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
+
+    # Edit coefficient's header
     header = coeffnii.header.copy()
-    coeffnii.header.set_qform(coeffnii.header.get_qform(coded=False), code=0)
-    coeffnii.header.set_sform(newaff, code=1)
+    header.set_qform(newaff, code=1)
+    header.set_sform(newaff, code=1)
 
+    # Write out fixed (generalized) coefficients
     coeffnii.__class__(coeffnii.dataobj, newaff, header).to_filename(out_file)
     return out_file
 
@@ -498,13 +534,15 @@ def _chunks(inlist, chunksize):
         yield str(p)
 
 
-def _b0_resampler(data, coeffs, pe, ro, hmc_xfm=None, unwarp=None, newpath=None):
+def _b0_resampler(in_file, coeffs, pe, ro, hmc_xfm=None, unwarp=None, newpath=None):
     """Outsource the resampler into a separate callable function to allow parallelization."""
     from functools import partial
+    from niworkflows.interfaces.nibabel import reorient_image
+    from sdcflows.utils.tools import ensure_positive_cosines
 
     # Prepare output names
     filename = partial(fname_presuffix, newpath=newpath)
-    retval = [filename(data, suffix=s) for s in ("_unwarped", "_xfm", "_field")]
+    retval = [filename(in_file, suffix=s) for s in ("_unwarped", "_xfm", "_field")]
 
     if unwarp is None:
         from sdcflows.transform import B0FieldTransform
@@ -520,12 +558,23 @@ def _b0_resampler(data, coeffs, pe, ro, hmc_xfm=None, unwarp=None, newpath=None)
 
         unwarp.xfm = Affine(XFMLoader.from_filename(hmc_xfm).to_ras())
 
-    if unwarp.fit(data):
-        unwarp.shifts.to_filename(retval[2])
+    # Reorient input to match that of the coefficients, i.e., to have positive director cosines
+    reoriented_img, axcodes = ensure_positive_cosines(nb.load(in_file))
+
+    if unwarp.fit(reoriented_img):
+        unwarp.mapped.to_filename(retval[2])
     else:
         retval[2] = None
 
-    unwarp.apply(nb.load(data), ro_time=ro, pe_dir=pe).to_filename(retval[0])
+    # Unwarp the reoriented image, and restore original orientation
+    unwarped_img = reorient_image(
+        unwarp.apply(reoriented_img, ro_time=ro, pe_dir=pe),
+        axcodes,
+    )
+    # Write out to disk
+    unwarped_img.to_filename(retval[0])
+
+    # Store the corresponding spatial transformation
     unwarp.to_displacements(ro_time=ro, pe_dir=pe).to_filename(retval[1])
 
     return retval

sdcflows/interfaces/epi.py

@@ -24,11 +24,13 @@
 
 from nipype.interfaces.base import (
     BaseInterfaceInputSpec,
-    TraitedSpec,
     File,
+    InputMultiObject,
+    OutputMultiObject,
+    SimpleInterface,
+    TraitedSpec,
     isdefined,
     traits,
-    SimpleInterface,
 )
 
 
@@ -64,3 +66,76 @@ def _run_interface(self, runtime):
             .replace("k", "z")
         )
         return runtime
+
+
+class _SortPEBlipsInputSpec(BaseInterfaceInputSpec):
+    in_data = InputMultiObject(
+        File(exists=True),
+        mandatory=True,
+        desc="list of input data",
+    )
+    pe_dirs_fsl = InputMultiObject(
+        traits.Enum("x", "x-", "y", "y-", "z", "z-"),
+        mandatory=True,
+        desc="list of PE directions, in FSL's conventions",
+    )
+    readout_times = InputMultiObject(
+        traits.Float,
+        mandatory=True,
+        desc="list of total readout times"
+    )
+
+
+class _SortPEBlipsOutputSpec(TraitedSpec):
+    out_data = OutputMultiObject(
+        File(),
+        desc="list of input data",
+    )
+    pe_dirs = OutputMultiObject(
+        traits.Enum("i", "i-", "j", "j-", "k", "k-"),
+        desc="list of PE directions, in BIDS's conventions",
+    )
+    pe_dirs_fsl = OutputMultiObject(
+        traits.Enum("x", "x-", "y", "y-", "z", "z-"),
+        desc="list of PE directions, in FSL's conventions",
+    )
+    readout_times = OutputMultiObject(
+        traits.Float,
+        desc="list of total readout times"
+    )
+
+
+class SortPEBlips(SimpleInterface):
+    """Sort PE blips so they are consistently fed into TOPUP."""
+
+    input_spec = _SortPEBlipsInputSpec
+    output_spec = _SortPEBlipsOutputSpec
+
+    def _run_interface(self, runtime):
+        # Put sign first
+        blips = [
+            f"+{pe[0]}" if len(pe) == 1 else f"-{pe[0]}"
+            for pe in self.inputs.pe_dirs_fsl
+        ]
+        sorted_inputs = sorted(zip(
+            blips,
+            self.inputs.readout_times,
+            self.inputs.in_data,
+        ))
+
+        (
+            self._results["pe_dirs_fsl"],
+            self._results["readout_times"],
+            self._results["out_data"],
+        ) = zip(*sorted_inputs)
+
+        # Put sign back last
+        self._results["pe_dirs_fsl"] = [
+            pe[1] if pe.startswith("+") else f"{pe[1]}-"
+            for pe in self._results["pe_dirs_fsl"]
+        ]
+        self._results["pe_dirs"] = [
+            pe.replace("x", "i").replace("y", "j").replace("z", "k")
+            for pe in self._results["pe_dirs_fsl"]
+        ]
+        return runtime

sdcflows/interfaces/tests/test_bspline.py

@@ -88,6 +88,7 @@ def test_topup_coeffs(tmpdir, testdata_dir):
     result = TOPUPCoeffReorient(
         in_coeff=str(testdata_dir / "topup-coeff.nii.gz"),
         fmap_ref=str(testdata_dir / "epi.nii.gz"),
+        pe_dir="j",
     ).run()
 
     nii = nb.load(result.outputs.out_coeff)
@@ -102,11 +103,12 @@ def test_topup_coeffs(tmpdir, testdata_dir):
         TOPUPCoeffReorient(
             in_coeff="failing.nii.gz",
             fmap_ref=str(testdata_dir / "epi.nii.gz"),
+            pe_dir="j",
         ).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"))
+        _fix_topup_fieldcoeff("failing.nii.gz", str(testdata_dir / "epi.nii.gz"), "i")
 
 
 @pytest.mark.skipif(os.getenv("GITHUB_ACTIONS") == "true", reason="this is GH Actions")

sdcflows/interfaces/tests/test_epi.py

@@ -0,0 +1,45 @@
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+#
+# Copyright 2022 The NiPreps Developers <nipreps@gmail.com>
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+# We support and encourage derived works from this project, please read
+# about our expectations at
+#
+#     https://www.nipreps.org/community/licensing/
+#
+"""Test EPI interfaces."""
+from pathlib import Path
+from ..epi import SortPEBlips
+
+
+def test_sort_pe_blips(tmpdir):
+
+    tmpdir.chdir()
+
+    input_comb = [("x-", 0.08), ("x-", 0.04), ("y-", 0.05), ("y", 0.05), ("x", 0.05)]
+
+    fnames = []
+    for i in range(len(input_comb)):
+        fnames.append(f"file{i}.nii")
+        Path(fnames[-1]).write_text("")
+
+    result = SortPEBlips(
+        in_data=fnames,
+        pe_dirs_fsl=[pe for pe, _ in input_comb],
+        readout_times=[trt for _, trt in input_comb],
+    ).run()
+
+    assert result.outputs.out_data == [f"file{i}.nii" for i in (4, 3, 1, 0, 2)]