diff --git a/tests/test_utilities.py b/tests/test_utilities.py index fc375b7..93670d2 100644 --- a/tests/test_utilities.py +++ b/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) diff --git a/warpkit/distortion.py b/warpkit/distortion.py index bd3de22..f225f79 100644 --- a/warpkit/distortion.py +++ b/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 diff --git a/warpkit/utilities.py b/warpkit/utilities.py index 8b4c4ab..8a13067 100644 --- a/warpkit/utilities.py +++ b/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