From a9caf8f9b3a7a115e247f6f4b459d29fed97b74a Mon Sep 17 00:00:00 2001
From: Leonardo Ayala <leonardo.a027@gmail.com>
Date: Tue, 9 Dec 2025 21:42:03 +0100
Subject: [PATCH] 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(