RF: Pull out pedir reorientation into a tested routine

Author: effigies

sdcflows/interfaces/bspline.py

@@ -41,6 +41,7 @@
 )
 
 from sdcflows.transform import grid_bspline_weights
+from sdcflows.utils.tools import reorient_pedir
 
 
 LOW_MEM_BLOCK_SIZE = 1000
@@ -522,43 +523,27 @@ def _fix_topup_fieldcoeff(in_coeff, fmap_ref, pe_dir, out_file=None):
     if out_file is None:
         out_file = Path("coefficients.nii.gz").absolute()
 
-    # Load reference image
     refnii = nb.load(fmap_ref)
-
-    # Load coefficients
     coeffnii = nb.load(in_coeff)
 
-    # Orientations are [[axis, flip], [axis, flip], [axis, flip]] arrays,
-    # where RAS is "no change"
-    # e.g. PSL == [[1, -1], [2, 1], [0, -1]]
+    # Coefficients are in LAS space, and we will convert to RAS.
+    # Reorient the reference image and phase-encoding direction to RAS
     ref_ornt = nb.io_orientation(refnii.affine)
-    coeff_ornt = nb.io_orientation(coeffnii.affine)
-    # Get ref_ornt relative to coeff_ornt instead of relative to RAS
-    ref2coeff = nb.orientations.ornt_transform(ref_ornt, coeff_ornt).astype(int)
-
-    # Find the axis index and flip of PE direction in reference space
-    ref_pe_axis = "ijk".find(pe_dir[0])
-    if ref_pe_axis == -1:
-        ref_pe_axis = "xyz".index(pe_dir[0])
-
-    # Transform to coefficient space
-    coeff_pe_axis = "ijk"[ref2coeff[ref_pe_axis, 0]]
+    refnii_ras = refnii.as_reoriented(ref_ornt)
+    coeff_pe_dir = reorient_pedir(pe_dir, ref_ornt)
 
     # 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 coeff_pe_axis == "i" else 1.0
-    lr_axis = np.nonzero(coeff_ornt[:, 0] == 0)[0]
+    reverse_pe = -1.0 if coeff_pe_dir[0] == "i" else 1.0
     coeffnii = coeffnii.__class__(
-        reverse_pe * np.flip(np.asanyarray(coeffnii.dataobj), axis=lr_axis),
+        reverse_pe * np.flip(np.asanyarray(coeffnii.dataobj), axis=0),
         coeffnii.affine,
         coeffnii.header,
     )
-    # Ensure reference has positive director cosines
-    refnii_ras = nb.as_closest_canonical(refnii)
 
     # Get matrix of B-Spline control knots
     coeff_shape = np.array(coeffnii.shape[:3])

sdcflows/interfaces/utils.py

@@ -38,6 +38,8 @@
 )
 from nipype.interfaces.mixins import CopyHeaderInterface as _CopyHeaderInterface
 
+from sdcflows.utils.tools import reorient_pedir
+
 OBLIQUE_THRESHOLD_DEG = 0.5
 
 
@@ -202,17 +204,7 @@ def _run_interface(self, runtime):
 
         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 ""))
+        pe_dirs = [reorient_pedir(pe_dir, img2target) for pe_dir in self.inputs.pe_dir]
 
         self._results = dict(
             out_file=fname_presuffix(

sdcflows/utils/tools.py

@@ -154,3 +154,59 @@ def brain_masker(in_file, out_file=None, padding=5):
     img.__class__(data, img.affine, img.header).to_filename(out_brain)
 
     return str(out_brain), str(out_probseg), str(out_mask)
+
+
+def reorient_pedir(pe_dir, source_ornt=None, target_ornt=None):
+    """Return updated PhaseEncodingDirection to account for an image data array rotation
+
+    This function can accept one or two orientations.
+    Passing only a source orientation will rotate into RAS:
+
+    >>> from nibabel.orientations import axcodes2ornt
+    >>> reorient_pedir("j", axcodes2ornt("RAS"))
+    'j'
+    >>> reorient_pedir("i", axcodes2ornt("PSL"))
+    'j-'
+
+    Passing only a target_ornt will rotate from RAS:
+
+    >>> reorient_pedir("j", target_ornt=axcodes2ornt("RAS"))
+    'j'
+    >>> reorient_pedir("i", target_ornt=axcodes2ornt("PSL"))
+    'k-'
+
+    Passing both will rotate from source to target.
+
+    >>> reorient_pedir("j", axcodes2ornt("LPS"), axcodes2ornt("AIR"))
+    'i-'
+
+    Passing a transform orientation as source_ornt will perform the transform,
+    and passing it as ``target_ornt`` will invert the transform:
+
+    >>> from nibabel.orientations import ornt_transform
+    >>> xfm = ornt_transform(axcodes2ornt("LPS"), axcodes2ornt("AIR"))
+    >>> reorient_pedir("j", xfm)
+    'i-'
+    >>> reorient_pedir("j", target_ornt=xfm)
+    'k-'
+    """
+    from nibabel.orientations import ornt_transform
+
+    if source_ornt is None:
+        source_ornt = [[0, 1], [1, 1], [2, 1]]
+    if target_ornt is None:
+        target_ornt = [[0, 1], [1, 1], [2, 1]]
+
+    xfm = ornt_transform(source_ornt, target_ornt).astype(int)  # shape: (3, 2)
+
+    directions = "ijk" if pe_dir[0] in "ijk" else "xyz"
+
+    source_axis = directions.index(pe_dir[0])
+    source_flip = pe_dir[1:] == "-"
+
+    axis_xfm = xfm[source_axis, :]  # shape: (2,)
+
+    target_axis = directions[axis_xfm[0]]
+    target_flip = source_flip ^ (axis_xfm[1] == -1)
+
+    return f"{target_axis}-" if target_flip else target_axis