diff --git a/nibabel/orientations.py b/nibabel/orientations.py
index f1cdd228b..f19908782 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(-(R**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..7b4dfb495 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,41 @@ 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(
+        [
+            [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)
+
+    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():
     assert_array_equal(
         ornt_transform(