ENH: Add a utility to calculate obliquity of affines

[oesteban]

This implementation mimics the implementation of AFNI, and will be useful in #656 .

[codecov]

Codecov Report

Merging #815 into master will increase coverage by <.01%.
The diff coverage is 100%.

@@            Coverage Diff            @@
##           master    #815      +/-   ##
=========================================
+ Coverage    90.1%   90.1%   +<.01%     
=========================================
  Files          96      96              
  Lines       11911   11914       +3     
  Branches     2125    2125              
=========================================
+ Hits        10732   10735       +3     
  Misses        834     834              
  Partials      345     345

Impacted Files	Coverage Δ
nibabel/affines.py	100% <100%> (ø)	⬆️
nibabel/funcs.py	80.28% <0%> (-0.28%)	⬇️
nibabel/spatialimages.py	96.33% <0%> (ø)	⬆️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 1d91b72...891d85c. Read the comment docs.

[matthew-brett]

Thanks for this. Some quick thoughts:

• the test should go into a test function, and it needs some more tests (for not-oblique, for example).
• if the doctest stays, can you find something more obvious as an example? The exponential notation values are hard to read.
• Can you explain the algorithm? I was expecting you to take the dot products of the columns of unit-length, something like this:

R = affine[:-1, :-1] / vs
col_cosines = R @ R.T
upper_cosines = col_cosines[np.triu_indices(3, 1)]
angles = np.arccos(upper_cosines) 
np.any(np.abs(angles - np.pi / 2) > thresh)

[oesteban]

• the test should go into a test function, and it needs some more tests (for not-oblique, for example).

Sure, the latest commit should address this.

• if the doctest stays, can you find something more obvious as an example? The exponential notation values are hard to read.

The doctest did not stay, but now I'm generating the matrix in a more understandable way (first, generate an aligned matrix; second explicitly rotate it).

Can you explain the algorithm? I was expecting you to take the dot products of the columns of unit-length, something like this:

TBH I just transferred AFNI's implementation because it behaves differently when the affines are not alligned when applying transforms, so I needed to identify this situation for AFNI.

Now looking at the code, I believe AFNI's implementation is pretty close to your suggestion since it first projects each column of the affine to the canonical axes, then finds the projection for each axis (i.e. that .max(axis=1)) and finally it computes the arccos only on the smallest of these projected vectors (which will give the largest angle). Your implementation would be calculating the obliquity of the three axes.

I don't mind either way. However, I'd be interested in having an implementation parallel to AFNI's to be able to predict it's behavior. In other words, I can check whether your implementation is equivalent for the largest angle. If not, then rename the function to afni_obliquity or something that notes this particular.

[oesteban]

Comparing implementations

import numpy as np
from math import acos, pi as PI
from nibabel.eulerangles import euler2mat
from nibabel.affines import from_matvec, voxel_sizes
aligned = np.diag([2.0, 2.0, 2.3, 1.0]) 
aligned[:-1, -1] = [-10, -10, -7] 
R = from_matvec(euler2mat(x=0.09, y=0.001, z=0.001), [0.0, 0.0, 0.0]) 
oblique = R.dot(aligned)
vs = voxel_sizes(oblique)
O = oblique[:-1, :-1] / vs[np.newaxis, ...]
fig_merit = O.max(axis=1).min()
afni = abs(acos(fig_merit) * 180 / PI)
col_cosines = O @ O.T
upper_cosines = col_cosines[np.triu_indices(3, 1)]
angles = np.arccos(upper_cosines)
angles_deg = np.abs(angles - np.pi / 2)

>>> afni                                                                                                                                                                                                                                                                                                                                                                                                      
5.156994888317267
>>> angles_deg
array([0., 0., 0.])

I suspect your code is actually checking the orthogonality of the axes.

[effigies]

From a quick glance it looks like they're using the inf-norm of the cosines instead of the 2-norm, so perhaps you could parametrize the norm?

[oesteban]

In https://github.com/nipy/nibabel/pull/815/files#diff-fadc4e2f32f9100f1b5ce2384d00238fR311,

the affine[:-1, :-1] / vs[np.newaxis, ...] is implicitly using the L2-norm for normalization (the L2-norm is calculated in voxel_sizes).
I think the inf-norm is just used to project on axes.

[matthew-brett]

Oh - sorry - yes - my code is definitely checking the orthogonality of axes, and that's because I thought that's what you wanted to test. So now I realize I need to ask - what does 'obliquity' mean, exactly? Some rotation that puts the axes out of alignment with any permutation of the RAS+ axes? In that case R = affine[:-1, :-1] / vs gives the cosines of the angles of the axes with respect to the three canonical axes, right (np.eye(3) @ R)? In which case, M = np.abs(R).max(axis=1) (don't we need the abs?) gives the cosine of the angle corresponding to the closest match of any column of R to np.eye(3), for each of the three canonical axes (rows of np.eye(3)). And then the .min checks for the canonical axis with the worst match to the axes of R (furthest from 1). Then the algorithm is basically correct, and could also be written:

cosines = oblique[:-1, :-1] / vs
best_cosines = np.abs(cosines).max(axis=1)

and so on.

I'd rather stick with radians if possible, we're not using degrees anywhere else, I think.

By the way, I believe the the [np.newaxis, ...] in affine[:-1, :-1] / vs[np.newaxis, ...] is redundant.

[oesteban]

Oh - sorry - yes - my code is definitely checking the orthogonality of axes, and that's because I thought that's what you wanted to test. So now I realize I need to ask - what does 'obliquity' mean, exactly? Some rotation that puts the axes out of alignment with any permutation of the RAS+ axes?

Yes, in AFNI's terms they have an explanation (where they mix up the concept of "obliquity" with the orientation of NIfTI images) here.

I'd rather stick with radians if possible, we're not using degrees anywhere else, I think.

I'll add a degrees keyword argument, defaulting to False, so that we can directly compare to AFNI (after all, the initiative to add this is a result of a need to predict its behavior).

By the way, I believe the the [np.newaxis, ...] in affine[:-1, :-1] / vs[np.newaxis, ...] is redundant.

I preferred to be explicit to make sure the columns were normalized. But yes, both are equivalent. Beyond explicitness, I have no preference - let me know what aligns best with existing code.

[pep8speaks]

Hello @oesteban, Thank you for updating!

Cheers! There are no style issues detected in this Pull Request. 🍻 To test for issues locally, pip install flake8 and then run flake8 nibabel.

Comment last updated at 2019-09-27 15:52:05 UTC

[matthew-brett]

Do you agree with me about abs above?

If you prefer to make the broadcasting rule explicit with vs[np.newaxis, ...], then the eliipsis is confusing - the combination of the not-strictly-necessary new axis, and the ellipsis, looks like you are thinking of more than one dimension. I'd prefer vs[None, :] personally (or just vs).

The degrees keyword seems messy to me, can't we just convert to degrees outside the function to test against AFNI? As in obliquity(affine) * 180 / pi == my_afni_result ?

Can you put more explanation in the docstring aboutt what 'obliquity' is?

[oesteban]

@matthew-brett, did the last commit address all your concerns?

[matthew-brett]

Yes - thanks for your patience.

Can you also add a link to https://sscc.nimh.nih.gov/sscc/dglen/Obliquity in the docstring?

Is it worth checking this calculation for non-orthogonal axes? It will say 'oblique', but I suppose that is what is intended?

[oesteban]

Is it worth checking this calculation for non-orthogonal axes? It will say 'oblique', but I suppose that is what is intended?

Like adding a test where just one of the axes is not orthogonal to the plane where the other two are?.

[oesteban]

Thanks for your patience, anything else I should address here?

[matthew-brett]

OK with me - Chris - OK with you?

[effigies]

I didn't follow this very closely, but I have no objections.