Merge pull request #96 from matthew-brett/add-matvec-functions

NF - add matrix vector functions from nipy

Author: matthew-brett

nibabel/affines.py

@@ -76,6 +76,99 @@ def apply_affine(aff, pts):
     return res.reshape(shape)
 
 
+def to_matvec(transform):
+    """Split a transform into its matrix and vector components.
+
+    The tranformation must be represented in homogeneous coordinates and is
+    split into its rotation matrix and translation vector components.
+
+    Parameters
+    ----------
+    transform : array-like
+        NxM transform matrix in homogeneous coordinates representing an affine
+        transformation from an (N-1)-dimensional space to an (M-1)-dimensional
+        space. An example is a 4x4 transform representing rotations and
+        translations in 3 dimensions. A 4x3 matrix can represent a 2-dimensional
+        plane embedded in 3 dimensional space.
+
+    Returns
+    -------
+    matrix : (N-1, M-1) array
+        Matrix component of `transform`
+    vector : (M-1,) array
+        Vector compoent of `transform`
+
+    See Also
+    --------
+    from_matvec
+
+    Examples
+    --------
+    >>> aff = np.diag([2, 3, 4, 1])
+    >>> aff[:3,3] = [9, 10, 11]
+    >>> to_matvec(aff)
+    (array([[2, 0, 0],
+           [0, 3, 0],
+           [0, 0, 4]]), array([ 9, 10, 11]))
+    """
+    transform = np.asarray(transform)
+    ndimin = transform.shape[0] - 1
+    ndimout = transform.shape[1] - 1
+    matrix = transform[0:ndimin, 0:ndimout]
+    vector = transform[0:ndimin, ndimout]
+    return matrix, vector
+
+
+def from_matvec(matrix, vector=None):
+    """ Combine a matrix and vector into an homogeneous affine
+
+    Combine a rotation / scaling / shearing matrix and translation vector into a
+    transform in homogeneous coordinates.
+
+    Parameters
+    ----------
+    matrix : array-like
+        An NxM array representing the the linear part of the transform.
+        A transform from an M-dimensional space to an N-dimensional space.
+    vector : None or array-like, optional
+        None or an (N,) array representing the translation. None corresponds to
+        an (N,) array of zeros.
+
+    Returns
+    -------
+    xform : array
+        An (N+1, M+1) homogenous transform matrix.
+
+    See Also
+    --------
+    to_matvec
+
+    Examples
+    --------
+    >>> from_matvec(np.diag([2, 3, 4]), [9, 10, 11])
+    array([[ 2,  0,  0,  9],
+           [ 0,  3,  0, 10],
+           [ 0,  0,  4, 11],
+           [ 0,  0,  0,  1]])
+
+    The `vector` argument is optional:
+
+    >>> from_matvec(np.diag([2, 3, 4]))
+    array([[2, 0, 0, 0],
+           [0, 3, 0, 0],
+           [0, 0, 4, 0],
+           [0, 0, 0, 1]])
+    """
+    matrix = np.asarray(matrix)
+    nin, nout = matrix.shape
+    t = np.zeros((nin+1,nout+1), matrix.dtype)
+    t[0:nin, 0:nout] = matrix
+    t[nin, nout] = 1.
+    if not vector is None:
+        t[0:nin, nout] = vector
+    return t
+
+
 def append_diag(aff, steps, starts=()):
     """ Add diagonal elements `steps` and translations `starts` to affine
 

nibabel/tests/test_affines.py

@@ -3,7 +3,8 @@
 
 import numpy as np
 
-from ..affines import apply_affine, append_diag
+from ..affines import (apply_affine, append_diag, to_matvec, from_matvec)
+
 
 from nose.tools import assert_equal, assert_raises
 from numpy.testing import assert_array_equal, assert_almost_equal, \
@@ -65,6 +66,30 @@ def test_apply_affine():
         assert_array_almost_equal(res, exp_res)
 
 
+def test_matrix_vector():
+    for M, N in ((4,4), (5,4), (4, 5)):
+        xform = np.zeros((M, N))
+        xform[:-1,:] = np.random.normal(size=(M-1, N))
+        xform[-1,-1] = 1
+        newmat, newvec = to_matvec(xform)
+        mat = xform[:-1, :-1]
+        vec = xform[:-1, -1]
+        assert_array_equal(newmat, mat)
+        assert_array_equal(newvec, vec)
+        assert_equal(newvec.shape, (M-1,))
+        assert_array_equal(from_matvec(mat, vec), xform)
+        # Check default translation works
+        xform_not = xform[:]
+        xform_not[:-1,:] = 0
+        assert_array_equal(from_matvec(mat), xform)
+        assert_array_equal(from_matvec(mat, None), xform)
+    # Check array-like works
+    newmat, newvec = to_matvec(xform.tolist())
+    assert_array_equal(newmat, mat)
+    assert_array_equal(newvec, vec)
+    assert_array_equal(from_matvec(mat.tolist(), vec.tolist()), xform)
+
+
 def test_append_diag():
     # Routine for appending diagonal elements
     assert_array_equal(append_diag(np.diag([2,3,1]), [1]),