enh: address @matthew-brett's comments

Author: oesteban

nibabel/affines.py

@@ -2,8 +2,6 @@
 # vi: set ft=python sts=4 ts=4 sw=4 et:
 """ Utility routines for working with points and affine transforms
 """
-
-from math import acos, pi as PI
 import numpy as np
 
 from functools import reduce
@@ -299,17 +297,37 @@ def voxel_sizes(affine):
     return np.sqrt(np.sum(top_left ** 2, axis=0))
 
 
-def obliquity(affine, degrees=False):
+def obliquity(affine, all_angles=True):
     r"""
-    Estimate the obliquity an affine's axes represent, in degrees from plumb.
+    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>`_
 
+    Parameters
+    ----------
+    affine : 2D array-like
+        Affine transformation array.  Usually shape (4, 4), but can be any 2D
+        array.
+    all_angles : bool
+        Whether the obliquity of the three axes should be returned, return the
+        *obliquity* of the most deviated axis otherwise.
+
+    Returns
+    -------
+    angles : 1D array-like or float
+        The *obliquity* of each axis with respect to the cardinal axes, in radians.
+        If ``all_angles`` is ``False``, return only that value for the most
+        *oblique* axis.
+
     """
     vs = voxel_sizes(affine)
-    fig_merit = (affine[:-1, :-1] / vs[np.newaxis, ...]).max(axis=1).min()
-    radians = abs(acos(fig_merit))
-    if not degrees:
-        return radians
-    return radians * 180 / PI
+    best_cosines = np.abs((affine[:-1, :-1] / vs).max(axis=1))
+    angles = np.arccos(best_cosines)
+
+    if all_angles:
+        return angles
+
+    return angles.max()

nibabel/tests/test_affines.py

@@ -182,9 +182,13 @@ def test_voxel_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)
-    assert_almost_equal(obliquity(oblique, degrees=True), 5.1569948883)
+    assert_almost_equal(obliquity(aligned), [0.0, 0.0, 0.0])
+    assert_almost_equal(obliquity(oblique) * 180 / pi,
+                        [0.0810285, 5.1569949, 5.1569376])
+    assert_almost_equal(obliquity(oblique, all_angles=False) * 180 / pi,
+                        5.1569949)