From 3cfff0d14f5d36950a2661bb262c52ffabc66bac Mon Sep 17 00:00:00 2001 From: oesteban Date: Thu, 19 Sep 2019 21:31:47 -0700 Subject: [PATCH 1/6] ENH: Add an utility to calculate obliquity of affines This implementation mimics the implementation of AFNI. --- nibabel/affines.py | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/nibabel/affines.py b/nibabel/affines.py index a7d7a4e9b8..03ab8b80e8 100644 --- a/nibabel/affines.py +++ b/nibabel/affines.py @@ -3,6 +3,7 @@ """ Utility routines for working with points and affine transforms """ +from math import acos, pi as PI import numpy as np from functools import reduce @@ -296,3 +297,24 @@ 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, in degrees from plumb. + + This implementation is inspired by `AFNI's implementation + `_ + + >>> affine = np.array([ + ... [2.74999725e+00, -2.74999817e-03, 2.74999954e-03, -7.69847980e+01], + ... [2.98603540e-03, 2.73886840e+00, -2.47165887e-01, -8.36692043e+01], + ... [-2.49170222e-03, 2.47168626e-01, 2.73886865e+00, -8.34056848e+01], + ... [ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.00000000e+00]]) + >>> abs(5.15699 - obliquity(affine)) < 0.0001 + True + + """ + vs = voxel_sizes(affine) + fig_merit = (affine[:-1, :-1] / vs[np.newaxis, ...]).max(axis=1).min() + return abs(acos(fig_merit) * 180 / PI) From c92d5609f37ac2ec18f94606b8bb3075759e27ea Mon Sep 17 00:00:00 2001 From: oesteban Date: Fri, 20 Sep 2019 11:26:42 -0700 Subject: [PATCH 2/6] enh(tests): add not-oblique test, move tests to ``test_affines.py`` Addressing one of @matthew-brett's comments. --- nibabel/affines.py | 8 -------- nibabel/tests/test_affines.py | 12 +++++++++++- 2 files changed, 11 insertions(+), 9 deletions(-) diff --git a/nibabel/affines.py b/nibabel/affines.py index 03ab8b80e8..bc800fd318 100644 --- a/nibabel/affines.py +++ b/nibabel/affines.py @@ -306,14 +306,6 @@ def obliquity(affine): This implementation is inspired by `AFNI's implementation `_ - >>> affine = np.array([ - ... [2.74999725e+00, -2.74999817e-03, 2.74999954e-03, -7.69847980e+01], - ... [2.98603540e-03, 2.73886840e+00, -2.47165887e-01, -8.36692043e+01], - ... [-2.49170222e-03, 2.47168626e-01, 2.73886865e+00, -8.34056848e+01], - ... [ 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.00000000e+00]]) - >>> abs(5.15699 - obliquity(affine)) < 0.0001 - True - """ vs = voxel_sizes(affine) fig_merit = (affine[:-1, :-1] / vs[np.newaxis, ...]).max(axis=1).min() diff --git a/nibabel/tests/test_affines.py b/nibabel/tests/test_affines.py index e66ed46190..143c80312b 100644 --- a/nibabel/tests/test_affines.py +++ b/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,13 @@ 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.""" + 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), 5.1569948883) From da8da564adfccee6e87ec3695b9ce0aa55fa1be9 Mon Sep 17 00:00:00 2001 From: oesteban Date: Sat, 21 Sep 2019 17:46:12 -0700 Subject: [PATCH 3/6] enh: return radians unless degrees=True --- nibabel/affines.py | 7 +++++-- nibabel/tests/test_affines.py | 2 +- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/nibabel/affines.py b/nibabel/affines.py index bc800fd318..462798b6f1 100644 --- a/nibabel/affines.py +++ b/nibabel/affines.py @@ -299,7 +299,7 @@ def voxel_sizes(affine): return np.sqrt(np.sum(top_left ** 2, axis=0)) -def obliquity(affine): +def obliquity(affine, degrees=False): r""" Estimate the obliquity an affine's axes represent, in degrees from plumb. @@ -309,4 +309,7 @@ def obliquity(affine): """ vs = voxel_sizes(affine) fig_merit = (affine[:-1, :-1] / vs[np.newaxis, ...]).max(axis=1).min() - return abs(acos(fig_merit) * 180 / PI) + radians = abs(acos(fig_merit)) + if not degrees: + return radians + return radians * 180 / PI diff --git a/nibabel/tests/test_affines.py b/nibabel/tests/test_affines.py index 143c80312b..db4fbe0630 100644 --- a/nibabel/tests/test_affines.py +++ b/nibabel/tests/test_affines.py @@ -187,4 +187,4 @@ def test_obliquity(): 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), 5.1569948883) + assert_almost_equal(obliquity(oblique, degrees=True), 5.1569948883) From 253a256d18d75664001ef3d9872c25952675b7f4 Mon Sep 17 00:00:00 2001 From: oesteban Date: Mon, 23 Sep 2019 09:31:42 -0700 Subject: [PATCH 4/6] enh: address @matthew-brett's comments --- nibabel/affines.py | 26 +++++++++++++++++--------- nibabel/tests/test_affines.py | 6 ++++-- 2 files changed, 21 insertions(+), 11 deletions(-) diff --git a/nibabel/affines.py b/nibabel/affines.py index 462798b6f1..25c9a4cdbd 100644 --- a/nibabel/affines.py +++ b/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,27 @@ def voxel_sizes(affine): return np.sqrt(np.sum(top_left ** 2, axis=0)) -def obliquity(affine, degrees=False): +def obliquity(affine): 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 `_ + 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) - 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)) + return np.arccos(best_cosines) diff --git a/nibabel/tests/test_affines.py b/nibabel/tests/test_affines.py index db4fbe0630..13b554c5a8 100644 --- a/nibabel/tests/test_affines.py +++ b/nibabel/tests/test_affines.py @@ -182,9 +182,11 @@ 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]) From 04584b6d8e0e50928ee7c5ff09de61471d6e5acb Mon Sep 17 00:00:00 2001 From: oesteban Date: Tue, 24 Sep 2019 10:20:19 -0700 Subject: [PATCH 5/6] doc: add link to AFNI's documentation about *obliquity* [skip ci] --- nibabel/affines.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/nibabel/affines.py b/nibabel/affines.py index 25c9a4cdbd..b09257c7ad 100644 --- a/nibabel/affines.py +++ b/nibabel/affines.py @@ -304,7 +304,9 @@ def obliquity(affine): 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 - `_ + `_. + For further details about *obliquity*, check `AFNI's documentation + _. Parameters ---------- From 891d85c07b41f4908fb3e9ffb9e89c78883a64eb Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Fri, 27 Sep 2019 08:52:01 -0700 Subject: [PATCH 6/6] fix: incorrect operation ordering caught by @matthew-brett --- nibabel/affines.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nibabel/affines.py b/nibabel/affines.py index b09257c7ad..9a37cc9e49 100644 --- a/nibabel/affines.py +++ b/nibabel/affines.py @@ -321,5 +321,5 @@ def obliquity(affine): """ vs = voxel_sizes(affine) - best_cosines = np.abs((affine[:-1, :-1] / vs).max(axis=1)) + best_cosines = np.abs(affine[:-1, :-1] / vs).max(axis=1) return np.arccos(best_cosines)