ENH: Add a utility to calculate obliquity of affines

nibabel/affines.py

@@ -2,7 +2,6 @@
 # vi: set ft=python sts=4 ts=4 sw=4 et:
 """ Utility routines for working with points and affine transforms
 """
-
 import numpy as np
 
 from functools import reduce
@@ -296,3 +295,31 @@ def voxel_sizes(affine):
     """
     top_left = affine[:-1, :-1]
     return np.sqrt(np.sum(top_left ** 2, axis=0))
+
+
+def obliquity(affine):
+    r"""
+    Estimate the *obliquity* an affine's axes represent.
+
+    The term *obliquity* is defined here as the rotation of those axes with
+    respect to the cardinal axes.
+    This implementation is inspired by `AFNI's implementation
+    <https://github.com/afni/afni/blob/b6a9f7a21c1f3231ff09efbd861f8975ad48e525/src/thd_coords.c#L660-L698>`_.
+    For further details about *obliquity*, check `AFNI's documentation
+    <https://sscc.nimh.nih.gov/sscc/dglen/Obliquity>_.
+
+    Parameters
+    ----------
+    affine : 2D array-like
+        Affine transformation array.  Usually shape (4, 4), but can be any 2D
+        array.
+
+    Returns
+    -------
+    angles : 1D array-like
+        The *obliquity* of each axis with respect to the cardinal axes, in radians.
+
+    """
+    vs = voxel_sizes(affine)
+    best_cosines = np.abs(affine[:-1, :-1] / vs).max(axis=1)
+    return np.arccos(best_cosines)

nibabel/tests/test_affines.py

@@ -7,7 +7,7 @@
 
 from ..eulerangles import euler2mat
 from ..affines import (AffineError, apply_affine, append_diag, to_matvec,
-                       from_matvec, dot_reduce, voxel_sizes)
+                       from_matvec, dot_reduce, voxel_sizes, obliquity)
 
 
 from nose.tools import assert_equal, assert_raises
@@ -178,3 +178,15 @@ def test_voxel_sizes():
             rot_affine[:3, :3] = rotation
             full_aff = rot_affine.dot(aff)
             assert_almost_equal(voxel_sizes(full_aff), vox_sizes)
+
+
+def test_obliquity():
+    """Check the calculation of inclination of an affine axes."""
+    from math import pi
+    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)
+    assert_almost_equal(obliquity(aligned), [0.0, 0.0, 0.0])
+    assert_almost_equal(obliquity(oblique) * 180 / pi,
+                        [0.0810285, 5.1569949, 5.1569376])