nipy
Saving image after loading as_closest_canonical swaps axis
I have noticed that in some cases (not always) saving an image after being loaded using nib.as_closest_canonical and then re-loading it in the same way causes axis to be transposed. To illustrate this the attached image and the code below can be used.
import nibabel as nib
from pathlib import Path
f = Path("test.nii.gz")
image = nib.as_closest_canonical(nib.load(f))
nib.save(nib.Nifti1Image(image.get_fdata(), image.affine, image.header), "test-canonical.nii.gz")
image_reloaded = nib.as_closest_canonical(nib.load("test-canonical.nii.gz"))
print(image.shape, image_reloaded.shape)
# (1024, 1024, 33) (33, 1024, 1024)
Expected behavior: The axis order should be consistent after saving an image and re-loading it.
Exemplary image: test.nii.gz
Huh, that's some bug.
In []: import nibabel as nb
In []: img = nb.load('Downloads/test.nii.gz')
In []: nb.aff2axcodes(img.affine)
Out[]: ('I', 'A', 'R')
In []: ras = nb.as_closest_canonical(img)
In []: nb.aff2axcodes(ras.affine)
Out[]: ('S', 'A', 'L')
In []: again = nb.as_closest_canonical(ras)
In []: nb.aff2axcodes(again.affine)
Out[]: ('I', 'A', 'R')
This appears to be a pretty weird affine:
In []: np.round(img.affine, 2)
Out[]:
array([[ 2.09, 0.08, 0.11, -16.75],
[ -1.04, 0.16, -0.05, 12.24],
[ -2.33, -0. , 0.12, -11. ],
[ 0. , 0. , 0. , 1. ]])
The first three columns are $i$, $j$, $k$ voxel dimensions, and the first three rows are the cosines with each world dimension $x$, $y$, $z$. Strangely, for both the $i$ and $k$ voxel dimensions, the increment is larger in the $z$ direction than $x$ or $y$. The way we calculate orientations is to just iterate over the voxel dimensions, find the world dimension it most closely matches of the unmatched world dimensions.
For simplicity, we can see the same pattern by pulling the rotation matrix out of the affine:
In []: import transforms3d as t3d
In []: T, R, Z, S = t3d.affines.decompose44(img.affine)
In []: R
Out[]:
array([[ 0.63260427, 0.43813722, 0.63862948],
[-0.3158175 , 0.89885806, -0.30383137],
[-0.70715708, -0.00948534, 0.70699285]])
$i$ points slightly more inferior than right, while $k$ points almost the same amount superior and right. So if we swap this around:
In []: t3d.affines.decompose44(ras.affine)[1]
Out[]:
array([[ 0.63862936, 0.43813722, -0.63260439],
[-0.30383133, 0.89885806, 0.31581754],
[ 0.70699297, -0.00948532, 0.70715696]])
We end up in the same situation. The first dimension corresponds to $z$ (this time in the S direction) and the third dimension is left to claim $x$ (this time in the L direction).
That said, if we run these through fslhd:
fslhd Downloads/test.nii.gz| grep '^s.*_[xyz]'
57: sto_xyz:1 2.087595 0.077025 0.112271 -16.753136
58: sto_xyz:2 -1.042198 0.158019 -0.053414 12.240506
59: sto_xyz:3 -2.333619 -0.001668 0.124289 -10.996345
60: sto_xyz:4 0.000000 0.000000 0.000000 1.000000
61: sform_xorient Superior-to-Inferior
62: sform_yorient Posterior-to-Anterior
63: sform_zorient Left-to-Right
fslhd test-canonical.nii.gz| grep '^s.*_[xyz]'
57: sto_xyz:1 0.112271 0.077025 -2.087595 50.049904
58: sto_xyz:2 -0.053414 0.158019 1.042198 -21.109837
59: sto_xyz:3 0.124289 -0.001668 2.333619 -85.672165
60: sto_xyz:4 0.000000 0.000000 0.000000 1.000000
61: sform_xorient Left-to-Right
62: sform_yorient Posterior-to-Anterior
63: sform_zorient Inferior-to-Superior
So FSL is at least able to do this stably. My guess is that it's working through the orientations by taking them from most to least overlap with a cardinal dimension, rather than just i, j, k, as we do here:
https://github.com/nipy/nibabel/blob/0f61bcbc1c2e4c47291e6ed5ba4c97b24f083002/nibabel/orientations.py#L69-L90
Instead of range(3), we could use np.argsort(np.min(-R**2, axis=0), kind='stable'):
In []: np.argsort(np.min(-R**2, axis=0), kind='stable')
Out[]: array([1, 0, 2])
I do think it would make sense to be able to consistently label a dimension the same axis, regardless of where it shows up in the order.
Would you be interested in submitting a pull request?
That's a huge problem. If you @ nibabel could fix it asap that would be great 🙏
Thanks for the quick reply @effigies, I will try your suggestion in the next days when I find the time :)
PRs are welcome. A fix should come with a regression test that fails prior to the fix and passes after.
If you want to use a small test file (you could just make one with the affected affine and a data shape of (2, 3, 4)), it can go in https://github.com/nipy/nibabel/tree/master/nibabel/tests/data. Otherwise, an image can be created on the fly.
The test should go in tests/test_orientations.py. You can use as_closest_canonical as a reference for how to use those objects.
https://github.com/nipy/nibabel/blob/0f61bcbc1c2e4c47291e6ed5ba4c97b24f083002/nibabel/funcs.py#L180-L209