RF: Reorient PE dir along with image pre-TOPUP

Author: effigies

sdcflows/interfaces/utils.py

@@ -21,6 +21,7 @@
 #     https://www.nipreps.org/community/licensing/
 #
 """Utilities."""
+from itertools import product
 from nipype.interfaces.base import (
     BaseInterfaceInputSpec,
     TraitedSpec,
@@ -147,6 +148,84 @@ def _run_interface(self, runtime):
         return runtime
 
 
+class _ReorientImageAndMetadataInputSpec(TraitedSpec):
+    in_file = File(exists=True, mandatory=True, desc="Input 3- or 4D image")
+    target_orientation = traits.Str(
+        desc="Axis codes of coordinate system to reorient to"
+    )
+    pe_dir = InputMultiObject(
+        traits.Enum(
+            *["".join(p) for p in product("ijkxyz", ("", "-"))],
+            mandatory=True,
+            desc="Phase encoding direction",
+        )
+    )
+
+
+class _ReorientImageAndMetadataOutputSpec(TraitedSpec):
+    out_file = File(desc="Reoriented image")
+    pe_dir = OutputMultiObject(
+        traits.Enum(
+            *["".join(p) for p in product("ijkxyz", ("", "-"))],
+            desc="Phase encoding direction in reoriented image",
+        )
+    )
+
+
+class ReorientImageAndMetadata(SimpleInterface):
+    input_spec = _ReorientImageAndMetadataInputSpec
+    output_spec = _ReorientImageAndMetadataOutputSpec
+
+    def _run_interface(self, runtime):
+        import numpy as np
+        import nibabel as nb
+        from nipype.utils.filemanip import fname_presuffix
+
+        target = self.inputs.target_orientation.upper()
+        if not all(code in "RASLPI" for code in target):
+            raise ValueError(
+                "Invalid orientation code {self.inputs.target_orientation}"
+            )
+
+        img = nb.load(self.inputs.in_file)
+        img2target = nb.orientations.ornt_transform(
+            nb.io_orientation(img.affine),
+            nb.orientations.axcodes2ornt(self.inputs.target_orientation),
+        ).astype(int)
+
+        # Identity transform
+        if np.array_equal(img2target, [[0, 1], [1, 1], [2, 1]]):
+            self._results = dict(
+                out_file=self.inputs.in_file,
+                pe_dir=self.inputs.pe_dir,
+            )
+
+        reoriented = img.as_reoriented(img2target)
+
+        pe_dirs = []
+        for pe_dir in self.inputs.pe_dir:
+            directions = "ijk" if pe_dir[0] in "ijk" else "xyz"
+
+            orig_axis = directions.index(pe_dir[0])
+            orig_flip = pe_dir[1:] == "-"
+
+            reoriented_axis = directions[img2target[orig_axis, 0]]
+            reoriented_flip = orig_flip ^ (img2target[orig_axis, 1] == -1)
+
+            pe_dirs.append(reoriented_axis + ("-" if reoriented_flip else ""))
+
+        self._results = dict(
+            out_file=fname_presuffix(
+                self.inputs.in_file, suffix=target, newpath=runtime.cwd
+            ),
+            pe_dir=pe_dirs,
+        )
+
+        reoriented.to_filename(self._results["out_file"])
+
+        return runtime
+
+
 class _ConvertWarpInputSpec(BaseInterfaceInputSpec):
     in_file = File(exists=True, mandatory=True, desc="output of 3dQwarp")
 

sdcflows/workflows/fit/pepolar.py

@@ -92,7 +92,7 @@ def init_topup_wf(
 
     from ...utils.misc import front as _front
     from ...interfaces.epi import GetReadoutTime, SortPEBlips
-    from ...interfaces.utils import UniformGrid, PadSlices
+    from ...interfaces.utils import UniformGrid, PadSlices, ReorientImageAndMetadata
     from ...interfaces.bspline import TOPUPCoeffReorient
     from ..ancillary import init_brainextraction_wf
 
@@ -145,7 +145,7 @@ def init_topup_wf(
     setwise_avg = pe.Node(RobustAverage(num_threads=omp_nthreads), name="setwise_avg")
     # The core of the implementation
     # Feed the input images in LAS orientation, so FSL does not run funky reorientations
-    to_las = pe.Node(ReorientImage(target_orientation="LAS"), name="to_las")
+    to_las = pe.Node(ReorientImageAndMetadata(target_orientation="LAS"), name="to_las")
     topup = pe.Node(
         TOPUP(
             config=_pkg_fname("sdcflows", f"data/flirtsch/b02b0{'_quick' * sloppy}.cnf")
@@ -172,16 +172,19 @@ def init_topup_wf(
         (regrid, sort_pe_blips, [("out_data", "in_data")]),
         (readout_time, sort_pe_blips, [("readout_time", "readout_times"),
                                        ("pe_dir_fsl", "pe_dirs_fsl")]),
-        (sort_pe_blips, topup, [("readout_times", "readout_times"),
-                                ("pe_dirs_fsl", "encoding_direction")]),
-        (sort_pe_blips, fix_coeff, [(("pe_dirs", _front), "pe_dir")]),
+        (sort_pe_blips, topup, [("readout_times", "readout_times")]),
         (setwise_avg, fix_coeff, [("out_file", "fmap_ref")]),
         (sort_pe_blips, concat_blips, [("out_data", "in_files")]),
         (concat_blips, pad_blip_slices, [("out_file", "in_file")]),
         (pad_blip_slices, setwise_avg, [("out_file", "in_file")]),
         (setwise_avg, to_las, [("out_hmc_volumes", "in_file")]),
-        (to_las, topup, [("out_file", "in_file")]),
+        (sort_pe_blips, to_las, [("pe_dirs_fsl", "pe_dir")]),
+        (to_las, topup, [
+            ("out_file", "in_file"),
+            ("pe_dir", "encoding_direction"),
+        ]),
         (topup, fix_coeff, [("out_fieldcoef", "in_coeff")]),
+        (to_las, fix_coeff, [(("pe_dir", _front), "pe_dir")]),
         (topup, outputnode, [("out_jacs", "jacobians"),
                              ("out_mats", "xfms")]),
         (ref_average, brainextraction_wf, [("out_file", "inputnode.in_file")]),