Merge pull request #128 from matthew-brett/fix-orientation-hangs

BF: fix nasty bug for orientations and 0 columns

Getting orientations for columns of all zeros in an affine resulted in
numpy svd hanging. I realized in the mean time that the io_orientation
function was not returning the correct number of rows. Extend tests,
fix, improve docstrings.

Author: matthew-brett

nibabel/funcs.py

@@ -10,7 +10,7 @@
 ''' Processor functions for images '''
 import numpy as np
 
-from .orientations import (io_orientation, orientation_affine, flip_axis,
+from .orientations import (io_orientation, inv_ornt_aff, flip_axis,
                            apply_orientation, OrientationError)
 from .loadsave import load
 
@@ -189,7 +189,7 @@ def as_closest_canonical(img, enforce_diag=False):
             raise OrientationError('Transformed affine is not diagonal')
         return img
     shape = img.shape
-    t_aff = orientation_affine(ornt, shape)
+    t_aff = inv_ornt_aff(ornt, shape)
     out_aff = np.dot(aff, t_aff)
     # check if we are going to end up with something diagonal
     if enforce_diag and not _aff_is_diag(aff):

nibabel/orientations.py

@@ -19,42 +19,43 @@ class OrientationError(Exception):
 def io_orientation(affine, tol=None):
     ''' Orientation of input axes in terms of output axes for `affine`
 
-    Valid for an affine transformation from ``m`` dimensions to ``n``
-    dimensions (``affine.shape == (n+1, m+1)``).
+    Valid for an affine transformation from ``p`` dimensions to ``q``
+    dimensions (``affine.shape == (q + 1, p + 1)``).
 
     The calculated orientations can be used to transform associated
-    arrays to best match the output orientations. If ``n`` > ``m``, then
+    arrays to best match the output orientations. If ``p`` > ``q``, then
     some of the output axes should be considered dropped in this
     orientation.
 
     Parameters
     ----------
-    affine : (n+1,m+1) ndarray-like
-       Transformation affine from ``m`` inputs to ``n`` outputs.
-       Usually this will be a shape (4,4) matrix, transforming 3 inputs
-       to 3 outputs, but the code also handles the more general case
+    affine : (q+1, p+1) ndarray-like
+       Transformation affine from ``p`` inputs to ``q`` outputs.  Usually this
+       will be a shape (4,4) matrix, transforming 3 inputs to 3 outputs, but the
+       code also handles the more general case
     tol : {None, float}, optional
-       threshold below which SVD values of the affine are considered
-       zero. If `tol` is None, and ``S`` is an array with singular
-       values for `affine`, and ``eps`` is the epsilon value for
-       datatype of ``S``, then `tol` set to ``S.max() * eps``.
+       threshold below which SVD values of the affine are considered zero. If
+       `tol` is None, and ``S`` is an array with singular values for `affine`,
+       and ``eps`` is the epsilon value for datatype of ``S``, then `tol` set to
+       ``S.max() * eps``.
 
     Returns
     -------
-    orientations : (n,2) ndarray
-       one row per input axis, where the first value in each row is the
-       closest corresponding output axis, and the second value is 1 if
-       the input axis is in the same direction as the corresponding
-       output axis, , and -1 if it is in the opposite direction, to the
-       corresponding output axis.  If a row is [np.nan, np.nan], which
-       can happen when n > m, then this row should be considered
-       dropped.
+    orientations : (p, 2) ndarray
+       one row per input axis, where the first value in each row is the closest
+       corresponding output axis. The second value in each row is 1 if the input
+       axis is in the same direction as the corresponding output axis and -1 if
+       it is in the opposite direction.  If a row is [np.nan, np.nan], which can
+       happen when p > q, then this row should be considered dropped.
     '''
     affine = np.asarray(affine)
-    n, m = affine.shape[0]-1, affine.shape[1]-1
-    # extract the underlying rotation matrix
-    RZS = affine[:n, :m]
+    q, p = affine.shape[0]-1, affine.shape[1]-1
+    # extract the underlying rotation, zoom, shear matrix
+    RZS = affine[:q, :p]
     zooms = np.sqrt(np.sum(RZS * RZS, axis=0))
+    # Zooms can be zero, in which case all elements in the column are zero, and
+    # we can leave them as they are
+    zooms[zooms == 0] = 1
     RS = RZS / zooms
     # Transform below is polar decomposition, returning the closest
     # shearless matrix R to RS
@@ -66,13 +67,13 @@ def io_orientation(affine, tol=None):
     R = np.dot(P[:, keep], Qs[keep])
     # the matrix R is such that np.dot(R,R.T) is projection onto the
     # columns of P[:,keep] and np.dot(R.T,R) is projection onto the rows
-    # of Qs[keep].  R (== np.dot(R, np.eye(m))) gives rotation of the
+    # of Qs[keep].  R (== np.dot(R, np.eye(p))) gives rotation of the
     # unit input vectors to output coordinates.  Therefore, the row
     # index of abs max R[:,N], is the output axis changing most as input
     # axis N changes.  In case there are ties, we choose the axes
     # iteratively, removing used axes from consideration as we go
-    ornt = np.ones((n, 2), dtype=np.int8) * np.nan
-    for in_ax in range(m):
+    ornt = np.ones((p, 2), dtype=np.int8) * np.nan
+    for in_ax in range(p):
         col = R[:, in_ax]
         if not np.alltrue(np.equal(col, 0)):
             out_ax = np.argmax(np.abs(col))
@@ -88,51 +89,6 @@ def io_orientation(affine, tol=None):
     return ornt
 
 
-def _ornt_to_affine(orientations):
-    ''' Create affine transformation matrix determined by orientations.
-
-    This transformation will simply flip, transpose, and possibly drop some
-    coordinates.
-
-    Parameters
-    ----------
-    orientations : (n,2) ndarray
-       one row per input axis, where the first value in each row is the
-       closest corresponding output axis, and the second value is 1 if
-       the input axis is in the same direction as the corresponding
-       output axis, and -1 if it is in the opposite direction, to the
-       corresponding output axis.  If a row has first entry np.nan, then
-       this axis is dropped from the output.
-
-    Returns
-    -------
-    affine : (m+1,n+1) ndarray
-       matrix representing flipping / dropping axes. m is equal to the
-       number of rows of orientations[:,0] that are not np.nan
-    '''
-    ornt = np.asarray(orientations)
-    n = ornt.shape[0]
-    keep = ~np.isnan(ornt[:, 1])
-    # These are the input coordinate axes that do have a matching output
-    # column in the orientation.  That is, if the 2nd row is [np.nan,
-    # np.nan] then the orientation indicates that no output axes of an
-    # affine with this orientation matches the 2nd input coordinate
-    # axis.  This would happen if the 2nd row of the affine that
-    # generated ornt was [0,0,0,*]
-    axes_kept = np.arange(n)[keep]
-    m = keep.sum(0)
-    # the matrix P represents the affine transform impled by ornt. If
-    # all entries of ornt are not np.nan, then P is square otherwise it
-    # has more columns than rows indicating some coordinates were
-    # dropped
-    P = np.zeros((m + 1, n + 1))
-    P[-1, -1] = 1
-    for idx in range(m):
-        axs, flip = ornt[axes_kept[idx]]
-        P[idx, axs] = flip
-    return P
-
-
 def apply_orientation(arr, ornt):
     ''' Apply transformations implied by `ornt` to the first
     n axes of the array `arr`
@@ -164,7 +120,6 @@ def apply_orientation(arr, ornt):
     if np.any(np.isnan(ornt[:, 0])):
         raise OrientationError('Cannot drop coordinates when '
                                'applying orientation to data')
-    shape = t_arr.shape
     # apply ornt transformations
     for ax, flip in enumerate(ornt[:, 1]):
         if flip == -1:
@@ -176,8 +131,8 @@ def apply_orientation(arr, ornt):
     return t_arr
 
 
-def orientation_affine(ornt, shape):
-    ''' Affine transform resulting from transforms implied in `ornt`
+def inv_ornt_aff(ornt, shape):
+    ''' Affine transform reversing transforms implied in `ornt`
 
     Imagine you have an array ``arr`` of shape `shape`, and you apply the
     transforms implied by `ornt` (more below), to get ``tarr``.
@@ -187,41 +142,55 @@ def orientation_affine(ornt, shape):
 
     Parameters
     ----------
-    ornt : (n,2) ndarray
-       orientation transform. ``ornt[N,1]` is flip of axis N of the
-       array implied by `shape`, where 1 means no flip and -1 means
-       flip.  For example, if ``N==0`` and ``ornt[0,1] == -1``, and
-       there's an array ``arr`` of shape `shape`, the flip would
-       correspond to the effect of ``np.flipud(arr)``.  ``ornt[:,0]`` is
-       the transpose that needs to be done to the implied array, as in
-       ``arr.transpose(ornt[:,0])``
-
-    shape : length n sequence
+    ornt : (p, 2) ndarray
+       orientation transform. ``ornt[P, 1]` is flip of axis N of the array
+       implied by `shape`, where 1 means no flip and -1 means flip.  For
+       example, if ``P==0`` and ``ornt[0, 1] == -1``, and there's an array
+       ``arr`` of shape `shape`, the flip would correspond to the effect of
+       ``np.flipud(arr)``.  ``ornt[:,0]`` gives us the (reverse of the)
+       transpose that has been done to ``arr``.  If there are any NaNs in
+       `ornt`, we raise an ``OrientationError`` (see notes)
+    shape : length p sequence
        shape of array you may transform with `ornt`
 
     Returns
     -------
-    transformed_affine : (n+1,n+1) ndarray
-       An array ``arr`` (shape `shape`) might be transformed according
-       to `ornt`, resulting in a transformed array ``tarr``.
-       `transformed_affine` is the transform that takes you from array
-       coordinates in ``tarr`` to array coordinates in ``arr``.
+    transform_affine : (p + 1, p + 1) ndarray
+       An array ``arr`` (shape `shape`) might be transformed according to
+       `ornt`, resulting in a transformed array ``tarr``.  `transformed_affine`
+       is the transform that takes you from array coordinates in ``tarr`` to
+       array coordinates in ``arr``.
+
+    Notes
+    -----
+    If a row in `ornt` contains NaN, this means that the input row does not
+    influence the output space, and is thus effectively dropped from the output
+    space.  In that case one ``tarr`` coordinate maps to many ``arr``
+    coordinates, we can't invert the transform, and we raise an error
     '''
     ornt = np.asarray(ornt)
-    n = ornt.shape[0]
-    shape = np.array(shape)[:n]
+    if np.any(np.isnan(ornt)):
+        raise OrientationError("We cannot invert orientation transform")
+    p = ornt.shape[0]
+    shape = np.array(shape)[:p]
     # ornt implies a flip, followed by a transpose.   We need the affine
     # that inverts these.  Thus we need the affine that first undoes the
     # effect of the transpose, then undoes the effects of the flip.
     # ornt indicates the transpose that has occurred to get the current
-    # ordering, relative to canonical, so we just use that
-    undo_reorder = np.eye(n + 1)[list(ornt[:, 0]) + [n], :]
+    # ordering, relative to canonical, so we just use that.
+    # undo_reorder is a row permutatation matrix
+    undo_reorder = np.eye(p + 1)[list(ornt[:, 0]) + [p], :]
     undo_flip = np.diag(list(ornt[:, 1]) + [1.0])
     center_trans = -(shape - 1) / 2.0
-    undo_flip[:n, n] = (ornt[:, 1] * center_trans) - center_trans
+    undo_flip[:p, p] = (ornt[:, 1] * center_trans) - center_trans
     return np.dot(undo_flip, undo_reorder)
 
 
+@np.deprecate_with_doc("Please use inv_ornt_aff instead")
+def orientation_affine(ornt, shape):
+    return inv_ornt_aff(ornt, shape)
+
+
 def flip_axis(arr, axis=0):
     ''' Flip contents of `axis` in array `arr`
 

nibabel/tests/test_orientations.py

@@ -14,10 +14,11 @@
 
 from numpy.testing import assert_array_equal, assert_array_almost_equal
 
-from ..orientations import (io_orientation, orientation_affine,
-                            flip_axis, _ornt_to_affine,
-                            apply_orientation, OrientationError,
-                            ornt2axcodes, aff2axcodes)
+from ..orientations import (io_orientation, inv_ornt_aff, flip_axis,
+                            apply_orientation, OrientationError, ornt2axcodes,
+                            aff2axcodes)
+
+from ..affines import from_matvec, to_matvec
 
 
 IN_ARRS = [np.eye(4),
@@ -147,23 +148,34 @@ def test_flip_axis():
 
 
 def test_io_orientation():
-    shape = (2,3,4)
-    for in_arr, out_ornt in zip(IN_ARRS, OUT_ORNTS):
-        ornt = io_orientation(in_arr)
-        assert_array_equal(ornt, out_ornt)
-        taff = orientation_affine(ornt, shape)
-        assert_true(same_transform(taff, ornt, shape))
-        for axno in range(3):
-            arr = in_arr.copy()
-            ex_ornt = out_ornt.copy()
-            # flip the input axis in affine
-            arr[:,axno] *= -1
-            # check that result shows flip
-            ex_ornt[axno, 1] *= -1
-            ornt = io_orientation(arr)
-            assert_array_equal(ornt, ex_ornt)
-            taff = orientation_affine(ornt, shape)
+    for shape in ((2,3,4), (20, 15, 7)):
+        for in_arr, out_ornt in zip(IN_ARRS, OUT_ORNTS):
+            ornt = io_orientation(in_arr)
+            assert_array_equal(ornt, out_ornt)
+            taff = inv_ornt_aff(ornt, shape)
             assert_true(same_transform(taff, ornt, shape))
+            for axno in range(3):
+                arr = in_arr.copy()
+                ex_ornt = out_ornt.copy()
+                # flip the input axis in affine
+                arr[:,axno] *= -1
+                # check that result shows flip
+                ex_ornt[axno, 1] *= -1
+                ornt = io_orientation(arr)
+                assert_array_equal(ornt, ex_ornt)
+                taff = inv_ornt_aff(ornt, shape)
+                assert_true(same_transform(taff, ornt, shape))
+    # Test nasty hang for zero columns
+    rzs = np.c_[np.diag([2, 3, 4, 5]), np.zeros((4,3))]
+    arr = from_matvec(rzs, [15,16,17,18])
+    ornt = io_orientation(arr)
+    assert_array_equal(ornt, [[0, 1],
+                              [1, 1],
+                              [2, 1],
+                              [3, 1],
+                              [np.nan, np.nan],
+                              [np.nan, np.nan],
+                              [np.nan, np.nan]])
 
 
 def test_ornt2axcodes():
@@ -205,56 +217,8 @@ def test_aff2axcodes():
                  ('B', 'R', 'U'))
 
 
-def test_drop_coord():
-    # given a 5x4 affine from slicing an fmri,
-    # the orientations code should easily reorder and drop the t
-    # axis
-
-    # this affine has output coordinates '-y','z','x' and is at t=16
-    sliced_fmri_affine = np.array([[0,-1,0,3],
-                                   [0,0,2,5],
-                                   [3,0,0,4],
-                                   [0,0,0,16],
-                                   [0,0,0,1]])
-    ornt = io_orientation(sliced_fmri_affine)
-    affine_that_drops_t_reorders_and_flips = _ornt_to_affine(ornt)
-    final_affine = np.dot(affine_that_drops_t_reorders_and_flips, 
-                          sliced_fmri_affine)
-    # the output will be diagonal
-    # with the 'y' row having been flipped and the 't' row dropped
-    assert_array_equal(final_affine,
-                       np.array([[3,0,0,4],
-                                 [0,1,0,-3],
-                                 [0,0,2,5],
-                                 [0,0,0,1]]))
-
-
-def test_ornt_to_affine():
-    # this orientation indicates that the first output
-    # axis of the affine is closest to the vector [0,0,-1],
-    # the last is closest to [1,0,0] and 
-    # that the y coordinate ([0,1,0]) is dropped
-    ornt = [[2,-1],
-            [np.nan,np.nan],
-            [0,1]]
-    # the reordering/flipping is represented by an affine that 
-    # takes the 3rd output coordinate and maps it to the
-    # first, takes the 3rd, maps it to first and flips it
-    A = np.array([[0,0,-1,0],
-                  [1,0,0,0],
-                  [0,0,0,1]])
-    assert_array_equal(A, _ornt_to_affine(ornt))
-    # a more complicated example. only the 1st, 3rd and 6th
-    # coordinates appear in the output
-    ornt = [[3,-1],
-            [np.nan,np.nan],
-            [0,1],
-            [np.nan,np.nan],
-            [np.nan,np.nan],
-            [1,-1]]
-    B = np.array([[0,0,0,-1,0,0,0],
-                  [1,0,0,0,0,0,0],
-                  [0,-1,0,0,0,0,0],
-                  [0,0,0,0,0,0,1]])
-    assert_array_equal(B, _ornt_to_affine(ornt))
-
+def test_inv_ornt_aff():
+    # Extra tests for inv_ornt_aff routines (also tested in
+    # io_orientations test)
+    assert_raises(OrientationError, inv_ornt_aff,
+                  [[0, 1], [1, -1], [np.nan, np.nan]], (3, 4, 5))