From a9caf8f9b3a7a115e247f6f4b459d29fed97b74a Mon Sep 17 00:00:00 2001 From: Leonardo Ayala Date: Tue, 9 Dec 2025 21:42:03 +0100 Subject: [PATCH 1/2] BF - Fix io_orientation to process input axes by strength, ensuring consistent labeling. Add regression test for handling competing axes in orientation determination. --- nibabel/orientations.py | 5 ++- nibabel/tests/test_orientations.py | 52 ++++++++++++++++++++++++++++++ 2 files changed, 56 insertions(+), 1 deletion(-) diff --git a/nibabel/orientations.py b/nibabel/orientations.py index f1cdd228b..75f0890da 100644 --- a/nibabel/orientations.py +++ b/nibabel/orientations.py @@ -75,7 +75,10 @@ def io_orientation(affine, tol=None): # axis N changes. In case there are ties, we choose the axes # iteratively, removing used axes from consideration as we go ornt = np.ones((p, 2), dtype=np.int8) * np.nan - for in_ax in range(p): + # Process input axes from strongest to weakest (stable on ties) so a given + # dimension is labeled consistently regardless of the original order. + in_axes = np.argsort(np.min(-(RS**2), axis=0), kind='stable') + for in_ax in in_axes: col = R[:, in_ax] if not np.allclose(col, 0): out_ax = np.argmax(np.abs(col)) diff --git a/nibabel/tests/test_orientations.py b/nibabel/tests/test_orientations.py index e7c32d786..b7ea73ed0 100644 --- a/nibabel/tests/test_orientations.py +++ b/nibabel/tests/test_orientations.py @@ -13,6 +13,7 @@ from numpy.testing import assert_array_equal from ..affines import from_matvec, to_matvec +from ..nifti1 import Nifti1Image from ..orientations import ( OrientationError, aff2axcodes, @@ -291,6 +292,57 @@ def test_io_orientation(): ) +def test_io_orientation_column_strength_regression(): + # Build a small image using the real-world affine that motivated the + # stronger column ordering. + affine = np.array( + [ + [1.12271041e-01, 7.70245194e-02, -2.08759499e00, 5.00499039e01], + [-5.34135476e-02, 1.58019245e-01, 1.04219818e00, -2.11098356e01], + [1.24289364e-01, -1.66752085e-03, 2.33361936e00, -8.56721640e01], + [0.00000000e00, 0.00000000e00, 0.00000000e00, 1.00000000e00], + ] + ) + img = Nifti1Image(np.zeros((2, 3, 4), dtype=np.float32), affine) + + # Sanity check: current orientation for the provided affine. + assert_array_equal(io_orientation(img.affine), [[0, 1], [1, 1], [2, 1]]) + + # Duplicate the first column to make two axes compete for the same output + # axis. The fixed code (ordering by RS strengths) keeps the strongest axis + # and drops the duplicate; the buggy SVD-ordered variant would pick the + # wrong column. + dup_affine = affine.copy() + dup_affine[:3, 1] = dup_affine[:3, 0] + expected = np.array([[0, 1], [1, -1], [2, 1]], dtype=np.int8) + assert_array_equal(io_orientation(dup_affine), expected) + + def buggy_io_orientation(aff): + # Replicates the pre-fix iteration order (range(p)) that could flip + # assignments when columns compete for the same output axis. + q, p = aff.shape[0] - 1, aff.shape[1] - 1 + rzs = aff[:q, :p] + zooms = np.sqrt(np.sum(rzs * rzs, axis=0)) + zooms[zooms == 0] = 1 + rs = rzs / zooms + P, S, Qs = np.linalg.svd(rs, full_matrices=False) + tol = S.max() * max(rs.shape) * np.finfo(S.dtype).eps + keep = S > tol + R = np.dot(P[:, keep], Qs[keep]) + ornt = np.ones((p, 2), dtype=np.int8) * np.nan + for in_ax in range(p): + col = R[:, in_ax] + if not np.allclose(col, 0): + out_ax = np.argmax(np.abs(col)) + ornt[in_ax, 0] = out_ax + ornt[in_ax, 1] = -1 if col[out_ax] < 0 else 1 + R[out_ax, :] = 0 + return ornt + + # check that the buggy orientation is not the expected orientation + assert not np.array_equal(buggy_io_orientation(dup_affine), expected) + + def test_ornt_transform(): assert_array_equal( ornt_transform( From 927be9130b42ff636f665ece106e8be6cc44208d Mon Sep 17 00:00:00 2001 From: Leonardo Ayala Date: Fri, 19 Dec 2025 23:34:07 +0100 Subject: [PATCH 2/2] BF - Refactors io_orientation to use rotation matrix for axes sorting. Updates regression test to verify orientation handling. --- nibabel/orientations.py | 2 +- nibabel/tests/test_orientations.py | 62 +++++++++++------------------- 2 files changed, 24 insertions(+), 40 deletions(-) diff --git a/nibabel/orientations.py b/nibabel/orientations.py index 75f0890da..f19908782 100644 --- a/nibabel/orientations.py +++ b/nibabel/orientations.py @@ -77,7 +77,7 @@ def io_orientation(affine, tol=None): ornt = np.ones((p, 2), dtype=np.int8) * np.nan # Process input axes from strongest to weakest (stable on ties) so a given # dimension is labeled consistently regardless of the original order. - in_axes = np.argsort(np.min(-(RS**2), axis=0), kind='stable') + in_axes = np.argsort(np.min(-(R**2), axis=0), kind='stable') for in_ax in in_axes: col = R[:, in_ax] if not np.allclose(col, 0): diff --git a/nibabel/tests/test_orientations.py b/nibabel/tests/test_orientations.py index b7ea73ed0..7b4dfb495 100644 --- a/nibabel/tests/test_orientations.py +++ b/nibabel/tests/test_orientations.py @@ -297,50 +297,34 @@ def test_io_orientation_column_strength_regression(): # stronger column ordering. affine = np.array( [ - [1.12271041e-01, 7.70245194e-02, -2.08759499e00, 5.00499039e01], - [-5.34135476e-02, 1.58019245e-01, 1.04219818e00, -2.11098356e01], - [1.24289364e-01, -1.66752085e-03, 2.33361936e00, -8.56721640e01], + [2.08759499e00, 7.70245194e-02, 1.12271041e-01, -1.67531357e01], + [-1.04219818e00, 1.58019245e-01, -5.34135476e-02, 1.22405062e01], + [-2.33361936e00, -1.66752085e-03, 1.24289364e-01, -1.09963446e01], [0.00000000e00, 0.00000000e00, 0.00000000e00, 1.00000000e00], ] ) + + # create image from afffine and reorient img = Nifti1Image(np.zeros((2, 3, 4), dtype=np.float32), affine) + print(aff2axcodes(img.affine)) + + # check that orientation is as expected + img_ornt = io_orientation(affine) + assert_array_equal(img_ornt, axcodes2ornt('IAR')) + + # reorient image to RAS + img_ras = img.as_reoriented(img_ornt) + img_ras_ornt = io_orientation(img_ras.affine) + + # Verify reorientation swaps first and third axes and the labels follow the axes + assert_array_equal(io_orientation(img.affine), [[2, -1], [1, 1], [0, 1]]) + + # Test idempotence + img_ras_reornt = img_ras.as_reoriented(img_ras_ornt) + img_ras_reornt_ornt = io_orientation(img_ras_reornt.affine) - # Sanity check: current orientation for the provided affine. - assert_array_equal(io_orientation(img.affine), [[0, 1], [1, 1], [2, 1]]) - - # Duplicate the first column to make two axes compete for the same output - # axis. The fixed code (ordering by RS strengths) keeps the strongest axis - # and drops the duplicate; the buggy SVD-ordered variant would pick the - # wrong column. - dup_affine = affine.copy() - dup_affine[:3, 1] = dup_affine[:3, 0] - expected = np.array([[0, 1], [1, -1], [2, 1]], dtype=np.int8) - assert_array_equal(io_orientation(dup_affine), expected) - - def buggy_io_orientation(aff): - # Replicates the pre-fix iteration order (range(p)) that could flip - # assignments when columns compete for the same output axis. - q, p = aff.shape[0] - 1, aff.shape[1] - 1 - rzs = aff[:q, :p] - zooms = np.sqrt(np.sum(rzs * rzs, axis=0)) - zooms[zooms == 0] = 1 - rs = rzs / zooms - P, S, Qs = np.linalg.svd(rs, full_matrices=False) - tol = S.max() * max(rs.shape) * np.finfo(S.dtype).eps - keep = S > tol - R = np.dot(P[:, keep], Qs[keep]) - ornt = np.ones((p, 2), dtype=np.int8) * np.nan - for in_ax in range(p): - col = R[:, in_ax] - if not np.allclose(col, 0): - out_ax = np.argmax(np.abs(col)) - ornt[in_ax, 0] = out_ax - ornt[in_ax, 1] = -1 if col[out_ax] < 0 else 1 - R[out_ax, :] = 0 - return ornt - - # check that the buggy orientation is not the expected orientation - assert not np.array_equal(buggy_io_orientation(dup_affine), expected) + assert img_ras_reornt.shape == (4, 3, 2) + assert_array_equal(img_ras_reornt_ornt, [[0, 1], [1, 1], [2, 1]]) def test_ornt_transform():