Skip to content

Commit 253a256

Browse files
committed
enh: address @matthew-brett's comments
1 parent da8da56 commit 253a256

File tree

2 files changed

+21
-11
lines changed

2 files changed

+21
-11
lines changed

nibabel/affines.py

+17-9
Original file line numberDiff line numberDiff line change
@@ -2,8 +2,6 @@
22
# vi: set ft=python sts=4 ts=4 sw=4 et:
33
""" Utility routines for working with points and affine transforms
44
"""
5-
6-
from math import acos, pi as PI
75
import numpy as np
86

97
from functools import reduce
@@ -299,17 +297,27 @@ def voxel_sizes(affine):
299297
return np.sqrt(np.sum(top_left ** 2, axis=0))
300298

301299

302-
def obliquity(affine, degrees=False):
300+
def obliquity(affine):
303301
r"""
304-
Estimate the obliquity an affine's axes represent, in degrees from plumb.
302+
Estimate the *obliquity* an affine's axes represent.
305303
304+
The term *obliquity* is defined here as the rotation of those axes with
305+
respect to the cardinal axes.
306306
This implementation is inspired by `AFNI's implementation
307307
<https://github.com/afni/afni/blob/b6a9f7a21c1f3231ff09efbd861f8975ad48e525/src/thd_coords.c#L660-L698>`_
308308
309+
Parameters
310+
----------
311+
affine : 2D array-like
312+
Affine transformation array. Usually shape (4, 4), but can be any 2D
313+
array.
314+
315+
Returns
316+
-------
317+
angles : 1D array-like
318+
The *obliquity* of each axis with respect to the cardinal axes, in radians.
319+
309320
"""
310321
vs = voxel_sizes(affine)
311-
fig_merit = (affine[:-1, :-1] / vs[np.newaxis, ...]).max(axis=1).min()
312-
radians = abs(acos(fig_merit))
313-
if not degrees:
314-
return radians
315-
return radians * 180 / PI
322+
best_cosines = np.abs((affine[:-1, :-1] / vs).max(axis=1))
323+
return np.arccos(best_cosines)

nibabel/tests/test_affines.py

+4-2
Original file line numberDiff line numberDiff line change
@@ -182,9 +182,11 @@ def test_voxel_sizes():
182182

183183
def test_obliquity():
184184
"""Check the calculation of inclination of an affine axes."""
185+
from math import pi
185186
aligned = np.diag([2.0, 2.0, 2.3, 1.0])
186187
aligned[:-1, -1] = [-10, -10, -7]
187188
R = from_matvec(euler2mat(x=0.09, y=0.001, z=0.001), [0.0, 0.0, 0.0])
188189
oblique = R.dot(aligned)
189-
assert_almost_equal(obliquity(aligned), 0.0)
190-
assert_almost_equal(obliquity(oblique, degrees=True), 5.1569948883)
190+
assert_almost_equal(obliquity(aligned), [0.0, 0.0, 0.0])
191+
assert_almost_equal(obliquity(oblique) * 180 / pi,
192+
[0.0810285, 5.1569949, 5.1569376])

0 commit comments

Comments
 (0)