BF - Fix io_orientation to process input axes by strength

[leoyala]

This PR introduces changes aiming to make the axis assignment based on image affines more consistent. Instead of iterating over the axis in order range(p), the axis are now iterated over based on the strongest axis according to the rotation-shear matrix (RS).

Alternative implementation:
Add a parameter to io_orientation called sort_by_strength which would activate this new sorting method, while preserving the old behavior. This parameter should then be propagated to as_closest_canonical and other places that might need it. In my opinion this might introduce more complexity and make it harder to maintain.

Notes:
In principle the sorting could be done based on the rotation matrix R instead of RS. However this causes test_io_orientation to fail due to the check when additional columns are added to the affine.

Closes #1449

[leoyala]

After testing the new implementation on a few hundred images that have non-standard orientations (e.g. oblique planes). The previous bug seems to be resolved for most cases, however I still found one image for which loading the image assigns what should be an Axial plane to the Sagittal one.
The affine of this image is:

[
[ 1.94222772e+00, -1.87633261e-01, -2.54019409e-01, 1.09614174e+02],
[ 1.56332731e+00,  3.34649354e-01, -7.97063559e-02, -1.08926643e+02],
[ 2.16188765e+00, -7.34265447e-02,  2.85847723e-01, -8.57995453e+01],
[ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00, 1.00000000e+00]
]

Even more, the old implementation (using range(p)) seems to load this image correctly.
Interestingly, running sflhd on this particular image returns the following, which seems to be the same as the new implementation.

sto_xyz:1	1.942228 -0.187633 -0.254019 109.614174 
sto_xyz:2	1.563327 0.334649 -0.079706 -108.926643 
sto_xyz:3	2.161888 -0.073427 0.285848 -85.799545 
sto_xyz:4	0.000000 0.000000 0.000000 1.000000 
sform_xorient	Left-to-Right
sform_yorient	Posterior-to-Anterior
sform_zorient	Inferior-to-Superior

Loading this image with another tool, like MITK, does seem to load the image in the correct orientation, assigning the axial plane to what the new implementation here assigns the sagittal one.
Any suggestions for improvement are welcomed :)

[effigies]

I think I may be thinking about this differently. The goal of an affine is to accurately describe the orientation of an image as a vector space, and the axis codes are just human conveniences. When your data is highly oblique, there just isn't an axial plane in the data, there's a phase-encoding, frequency-encoding and slice axis. We can assign a "nearest" world axis to each, but the only way to be wrong here is to use a procedure that (a) gives unintuitive results when each axis is dominantly aligned with a different world axis or (b) assigns different labels when the axes are reordered.

Looking at the affine you show, it is currently SAL by the naive implementation, but attempting to rotate to RAS will produce IAR.

Here's the rotation matrix:

array([[ 0.58855388, -0.48034116, -0.65028971],
       [ 0.47373557,  0.85670235, -0.20404827],
       [ 0.6551175 , -0.18797196,  0.73177018]])

What happens if you give MITK the as_closest_canonical'd version (old method)? My guess is that it will be inconsistent.

[codecov]

Codecov Report

❌ Patch coverage is 96.77419% with 1 line in your changes missing coverage. Please review.
✅ Project coverage is 95.42%. Comparing base (0f61bcb) to head (a9caf8f).

Files with missing lines	Patch %	Lines
nibabel/tests/test_orientations.py	96.55%	0 Missing and 1 partial ⚠️

Additional details and impacted files

@@           Coverage Diff           @@
##           master    #1450   +/-   ##
=======================================
  Coverage   95.42%   95.42%           
=======================================
  Files         209      209           
  Lines       29814    29844   +30     
  Branches     4483     4485    +2     
=======================================
+ Hits        28451    28480   +29     
  Misses        930      930           
- Partials      433      434    +1     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[effigies]

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.

[effigies]

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]])

[effigies]

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.

[effigies]

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.