Merge pull request #413 from matthew-brett/add-voxel-size

MRG: function to return voxel size from affine

Voxel sizes are the Euclidean lengths of the columns of the affine
(excluding the last column, and assuming zeros in the last row).

Author: matthew-brett

nibabel/affines.py

@@ -254,3 +254,45 @@ def dot_reduce(*args):
         ...  arg[N-2].dot(arg[N-1])))...``
     """
     return reduce(lambda x, y: np.dot(y, x), args[::-1])
+
+
+def voxel_sizes(affine):
+    r""" Return voxel size for each input axis given `affine`
+
+    The `affine` is the mapping between array (voxel) coordinates and mm
+    (world) coordinates.
+
+    The voxel size for the first voxel (array) axis is the distance moved in
+    world coordinates when moving one unit along the first voxel (array) axis.
+    This is the distance between the world coordinate of voxel (0, 0, 0) and
+    the world coordinate of voxel (1, 0, 0).  The world coordinate vector of
+    voxel coordinate vector (0, 0, 0) is given by ``v0 = affine.dot((0, 0, 0,
+    1)[:3]``.  The world coordinate vector of voxel vector (1, 0, 0) is
+    ``v1_ax1 = affine.dot((1, 0, 0, 1))[:3]``.  The final 1 in the voxel
+    vectors and the ``[:3]`` at the end are because the affine works on
+    homogenous coodinates.  The translations part of the affine is ``trans =
+    affine[:3, 3]``, and the rotations, zooms and shearing part of the affine
+    is ``rzs = affine[:3, :3]``. Because of the final 1 in the input voxel
+    vector, ``v0 == rzs.dot((0, 0, 0)) + trans``, and ``v1_ax1 == rzs.dot((1,
+    0, 0)) + trans``, and the difference vector is ``rzs.dot((0, 0, 0)) -
+    rzs.dot((1, 0, 0)) == rzs.dot((1, 0, 0)) == rzs[:, 0]``.  The distance
+    vectors in world coordinates between (0, 0, 0) and (1, 0, 0), (0, 1, 0),
+    (0, 0, 1) are given by ``rzs.dot(np.eye(3)) = rzs``.  The voxel sizes are
+    the Euclidean lengths of the distance vectors.  So, the voxel sizes are
+    the Euclidean lengths of the columns of the affine (excluding the last row
+    and column of the affine).
+
+    Parameters
+    ----------
+    affine : 2D array-like
+        Affine transformation array.  Usually shape (4, 4), but can be any 2D
+        array.
+
+    Returns
+    -------
+    vox_sizes : 1D array
+        Voxel sizes for each input axis of affine.  Usually 1D array length 3,
+        but in general has length (N-1) where input `affine` is shape (M, N).
+    """
+    top_left = affine[:-1, :-1]
+    return np.sqrt(np.sum(top_left ** 2, axis=0))

nibabel/tests/test_affines.py

@@ -1,10 +1,13 @@
 # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
 # vi: set ft=python sts=4 ts=4 sw=4 et:
 
+from itertools import product
+
 import numpy as np
 
+from ..eulerangles import euler2mat
 from ..affines import (AffineError, apply_affine, append_diag, to_matvec,
-                       from_matvec, dot_reduce)
+                       from_matvec, dot_reduce, voxel_sizes)
 
 
 from nose.tools import assert_equal, assert_raises
@@ -144,3 +147,34 @@ def test_dot_reduce():
                        np.dot(mat2, np.dot(vec, mat)))
     assert_array_equal(dot_reduce(mat, vec, mat2, ),
                        np.dot(mat, np.dot(vec, mat2)))
+
+
+def test_voxel_sizes():
+    affine = np.diag([2, 3, 4, 1])
+    assert_almost_equal(voxel_sizes(affine), [2, 3, 4])
+    # Some example rotations
+    rotations = []
+    for x_rot, y_rot, z_rot in product((0, 0.4), (0, 0.6), (0, 0.8)):
+        rotations.append(euler2mat(z_rot, y_rot, x_rot))
+    # Works on any size of array
+    for n in range(2, 10):
+        vox_sizes = np.arange(n) + 4.1
+        aff = np.diag(list(vox_sizes) + [1])
+        assert_almost_equal(voxel_sizes(aff), vox_sizes)
+        # Translations make no difference
+        aff[:-1, -1] = np.arange(n) + 10
+        assert_almost_equal(voxel_sizes(aff), vox_sizes)
+        # Does not have to be square
+        new_row = np.vstack((np.zeros(n + 1), aff))
+        assert_almost_equal(voxel_sizes(new_row), vox_sizes)
+        new_col = np.c_[np.zeros(n + 1), aff]
+        assert_almost_equal(voxel_sizes(new_col),
+                            [0] + list(vox_sizes))
+        if n < 3:
+            continue
+        # Rotations do not change the voxel size
+        for rotation in rotations:
+            rot_affine = np.eye(n + 1)
+            rot_affine[:3, :3] = rotation
+            full_aff = rot_affine.dot(aff)
+            assert_almost_equal(voxel_sizes(full_aff), vox_sizes)