Skip to content

BF - Fix io_orientation to process input axes by strength #1450

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
leoyala wants to merge 1 commit into
base: master
Choose a base branch
from
Open

BF - Fix io_orientation to process input axes by strength #1450

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 4 additions & 1 deletion nibabel/orientations.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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')
Copy link
Member

@effigies effigies Dec 10, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think RS is the right choice. The whole point of R is to produce a de-sheared rotation matrix.

In the case of the failing test, the RS matrix is:

array([[1.00000000e+00, 1.00000000e+00, 0.00000000e+00, 0.00000000e+00],
       [0.00000000e+00, 2.22044605e-15, 0.00000000e+00, 0.00000000e+00],
       [0.00000000e+00, 0.00000000e+00, 1.00000000e+00, 0.00000000e+00]])

Column 1 has a slightly higher magnitude than col 0, so ends up with a slightly higher maximum cosine in the desheared rotation (R):

array([[7.07106781e-01, 7.07106781e-01, 0.00000000e+00, 0.00000000e+00],
       [7.85046229e-16, 7.85046229e-16, 0.00000000e+00, 0.00000000e+00],
       [0.00000000e+00, 0.00000000e+00, 1.00000000e+00, 0.00000000e+00]])

If we passed in [[1, 1, 0, 0, 0], [0, 0, 0, 0, 0], [0, 0, 1, 0, 0], [0, 0, 0, 0, 1]], then we would have a tie in both the max and argmax of the two columns, so choosing the first column makes sense and the stable argsort will achieve this. In the test case, we have a tie in argmax, but we can break the tie with max, which also solves the problem we were seeing with non-degenerate matrices.

I think the thing to do here is actually to update the test to check three conditions:

  1. Put the positive error in the first column, yielding the current expected result.
  2. Put a negative error in the second column, yielding the current expected result.
  3. Put a positive error in the second column (current test), yielding a swap of the axis labeled nan.

I would be curious to hear @matthew-brett's thoughts here.

for in_ax in in_axes:
col = R[:, in_ax]
if not np.allclose(col, 0):
out_ax = np.argmax(np.abs(col))
Expand Down
52 changes: 52 additions & 0 deletions nibabel/tests/test_orientations.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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],
]
)
Comment on lines +298 to +305
Copy link
Member

@effigies effigies Dec 10, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is the kind of thing I was thinking of:

tricky_affine = <iar_affine>
orig_ornt = io_orientation(tricky_affine)
assert_array_equal(orig_ornt, axcodes2ornt('IAR'))

orig_img = Nifti1Image(np.zeros((2, 3, 4)), tricky_affine)
ras_img = orig_img.as_reoriented(orig_ornt)
ras_ornt = io_orientation(ras_img.affine)

# Verify reorientation swaps first and third axes and the labels follow the axes
assert ras_img.shape == (4, 3, 2)
assert_array_equal(ras_ornt, [[0, 1], [1, 1], [2, 1]])

# Test idempotence
reras_img = ras_img.as_reoriented(ras_ornt)
reras_ornt = io_orientation(reras_img.affine)

assert reras_img.shape == (4, 3, 2)
assert_array_equal(reras_ornt, [[0, 1], [1, 1], [2, 1]])

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)
Comment on lines +311 to +318
Copy link
Member

@effigies effigies Dec 10, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This doesn't really fit in this regression test, as it artificially changes the problem. The columns that were getting flipped were 1 and 3, not 1 and 2. What you're doing here is setting two columns to point mostly superior, then right, then posterior, and one to point mostly superior, then left, then anterior.

This is an interesting edge case, but not one that I think should push us into decent edge-case test, but then we need to agree that a rank-deficient matrix should be considered to be should be considered to have a particular, rather than just explainable, orientation.


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)
Comment on lines +320 to +343
Copy link
Member

@effigies effigies Dec 10, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For a regression test, we just need to have plausible affines and their correct orientations. We do not need to replicate the old function, though. The test should fail if we revert the fix or otherwise change the code to no longer do what's expected, and pass if the fix remains in place.



def test_ornt_transform():
assert_array_equal(
ornt_transform(
Expand Down