🐛 Invert MEDIC displacement in data space; drop redundant flip_sign

tests/test_utilities.py

@@ -235,19 +235,22 @@ def _zero_3vec_field(shape=(8, 8, 8)) -> nib.Nifti1Image:
     return nib.Nifti1Image(np.zeros((*shape, 3), dtype=np.float32), affine)
 
 
-def test_invert_displacement_field_zero():
+@pytest.mark.parametrize("ignore_rotation", [False, True])
+def test_invert_displacement_field_zero(ignore_rotation):
     """Inverting a zero field gives back zero (within numerical noise) and
-    preserves the 3-channel last axis."""
+    preserves the 3-channel last axis. Covers both the physical-space and the
+    ignore_rotation (voxel-grid frame) inversion paths."""
     field = _zero_3vec_field()
-    inverted = invert_displacement_field(field)
+    inverted = invert_displacement_field(field, ignore_rotation=ignore_rotation)
     assert inverted.shape == field.shape
     assert_allclose(inverted.get_fdata(), 0.0, atol=1e-5)
 
 
-def test_invert_displacement_maps_zero():
+@pytest.mark.parametrize("ignore_rotation", [False, True])
+def test_invert_displacement_maps_zero(ignore_rotation):
     affine = np.diag([2.0, 2.0, 2.0, 1.0])
     dmap = nib.Nifti1Image(np.zeros((6, 6, 6, 2), dtype=np.float32), affine)
-    inverted = invert_displacement_maps(dmap, axis="y")
+    inverted = invert_displacement_maps(dmap, axis="y", ignore_rotation=ignore_rotation)
     assert inverted.shape == dmap.shape
     assert_allclose(inverted.get_fdata(), 0.0, atol=1e-5)
 

warpkit/distortion.py

@@ -126,12 +126,14 @@ def medic(
 
     # invert displacement maps (these are in undistorted space)
     displacement_maps = invert_displacement_maps(
-        inv_displacement_maps, phase_encoding_direction
+        inv_displacement_maps, phase_encoding_direction, ignore_rotation=True
     )
 
-    # convert correction maps back to undistorted space field map
+    # convert correction maps back to undistorted space field map. No
+    # flip_sign here: the correlation check below sets the sign relative to
+    # the native field map, which makes an up-front flip redundant.
     field_maps = displacement_maps_to_field_maps(
-        displacement_maps, total_readout_time, phase_encoding_direction, flip_sign=True
+        displacement_maps, total_readout_time, phase_encoding_direction
     )
 
     # check if we need to flip sign of field maps (this is done by comparing sign of correlation between

warpkit/utilities.py

@@ -493,7 +493,10 @@ def get_ras_orient_transform(
 
 
 def invert_displacement_maps(
-    displacement_maps: nib.Nifti1Image, axis: str = "y", verbose: bool = False
+    displacement_maps: nib.Nifti1Image,
+    axis: str = "y",
+    verbose: bool = False,
+    ignore_rotation: bool = False,
 ) -> nib.Nifti1Image:
     """Invert displacement maps
 
@@ -505,6 +508,13 @@ def invert_displacement_maps(
         Axis displacement maps are along, by default "y"
     verbose : bool, optional
         Print debugging information, by default False
+    ignore_rotation : bool, optional
+        Invert in data (voxel) space rather than physical space by passing an
+        identity rotation to the inverter instead of the affine's rotation, by
+        default False. The displacement map is composed along a voxel axis, so
+        for oblique acquisitions the affine's rotation makes ITK invert in a
+        physical frame that does not align with that axis; setting this keeps
+        the inversion in the grid frame the map is actually defined in.
 
     Returns
     -------
@@ -527,6 +537,11 @@ def invert_displacement_maps(
     # split affine into components
     translations, rotations, zooms, _ = decompose44(displacement_maps_ras.affine)
 
+    # invert in data space: drop the affine's rotation so ITK inverts in the
+    # voxel grid frame the map is composed along, not the physical orientation
+    if ignore_rotation:
+        rotations = np.eye(3)
+
     # invert maps
     new_data = np.zeros(data.shape, dtype=np.float32)
     logging.info("Inverting displacement maps...")
@@ -554,7 +569,9 @@ def invert_displacement_maps(
 
 
 def invert_displacement_field(
-    displacement_field: nib.Nifti1Image, verbose: bool = False
+    displacement_field: nib.Nifti1Image,
+    verbose: bool = False,
+    ignore_rotation: bool = False,
 ) -> nib.Nifti1Image:
     """Invert displacement field
 
@@ -564,6 +581,10 @@ def invert_displacement_field(
         Displacement field data in mm
     verbose : bool, optional
         Print debugging information, by default False
+    ignore_rotation : bool, optional
+        Invert in data (voxel) space rather than physical space by passing an
+        identity rotation to the inverter instead of the affine's rotation, by
+        default False. See :func:`invert_displacement_maps` for details.
 
     Returns
     -------
@@ -583,6 +604,11 @@ def invert_displacement_field(
     # split affine into components
     translations, rotations, zooms, _ = decompose44(displacement_field_ras.affine)
 
+    # invert in data space: drop the affine's rotation so ITK inverts in the
+    # voxel grid frame rather than the physical orientation
+    if ignore_rotation:
+        rotations = np.eye(3)
+
     # Pad spatial dims only — leaving the 3-channel axis at size 3 — so we
     # can avoid edge effects of the inverse without inflating the channel
     # count. (Earlier versions padded all 4 axes which produced a 5-channel