NF: 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).

nibabel/affines.py

@@ -254,3 +254,36 @@ 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):
+    """ 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 array (voxel) axis is the distance between the
+    world coordinate of voxel (0, 0, 0) and the world coordinate of voxel (1,
+    0, 0).  For the second axis it is the world coordinate distance of (0, 0,
+    0) and (0, 1, 0) and so on.  The translations (last column of the affine)
+    are common between the world coordinates of (0, 0, 0) and those of (1, 0,
+    0), so the vector between (0, 0, 0) and (1, 0, 0) in world coordinates is
+    the given by (1, 0, 0) transformed by the top left (3, 3) portion of the
+    affine, and this is just the first column of the affine.  To get the
+    distance we take the Euclidean length of the vector.  So, the voxel sizes
+    are given by the Euclidean lengths of the columns 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,26 @@ 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])
+    # Rotations do not change the voxel size
+    for x_rot, y_rot, z_rot in product((0, 0.4), (0, 0.6), (0, 0.8)):
+        rotation = from_matvec(euler2mat(z_rot, y_rot, x_rot))
+        rot_affine = rotation.dot(affine)
+        assert_almost_equal(voxel_sizes(rot_affine), [2, 3, 4])
+        # Translations have no effect
+        rot_affine[:3, 3] = [10, 11, 12]
+        assert_almost_equal(voxel_sizes(rot_affine), [2, 3, 4])
+    # 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)
+    # Does not have to be square
+    aff = np.array([[2, 0, 0, 20],
+                    [0, 3, 0, 30],
+                    [0, 0, 0, 1]])
+    assert_almost_equal(voxel_sizes(aff), [2, 3, 0])