Merge branch 'to-output-space' into main-master

* to-output-space:
  RF: change to mapped_voxels as vox2out_vox input
  DOC: docstring and comment changes
  NF: function to make affine for selecting slice
  NF: add routine to make output space bounding box
  DOC: reformat affine docstring

Author: matthew-brett

nibabel/affines.py

@@ -15,24 +15,24 @@ def apply_affine(aff, pts):
     coordinate dimension of `pts` should be the last.
 
     For the 3D case, `aff` will be shape (4,4) and `pts` will have final axis
-    length 3 - maybe it will just be N by 3. The return value is the transformed
-    points, in this case::
+    length 3 - maybe it will just be N by 3. The return value is the
+    transformed points, in this case::
 
         res = np.dot(aff[:3,:3], pts.T) + aff[:3,3:4]
         transformed_pts = res.T
 
-    Notice though, that this routine is more general, in that `aff` can have any
-    shape (N,N), and `pts` can have any shape, as long as the last dimension is
-    for the coordinates, and is therefore length N-1.
+    This routine is more general than 3D, in that `aff` can have any shape
+    (N,N), and `pts` can have any shape, as long as the last dimension is for
+    the coordinates, and is therefore length N-1.
 
     Parameters
     ----------
     aff : (N, N) array-like
         Homogenous affine, for 3D points, will be 4 by 4. Contrary to first
         appearance, the affine will be applied on the left of `pts`.
     pts : (..., N-1) array-like
-        Points, where the last dimension contains the coordinates of each point.
-        For 3D, the last dimension will be length 3.
+        Points, where the last dimension contains the coordinates of each
+        point.  For 3D, the last dimension will be length 3.
 
     Returns
     -------

nibabel/spaces.py

@@ -0,0 +1,139 @@
+# emacs: -*- mode: python-mode; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
+#
+#   See COPYING file distributed along with the NiBabel package for the
+#   copyright and license terms.
+#
+### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
+""" Routines to work with spaces
+
+A space is defined by coordinate axes.
+
+A voxel space can be expressed by a shape implying an array, where the axes are
+the axes of the array.
+
+A mapped voxel space (mapped voxels) is either:
+
+* an image, with attributes ``shape`` (the voxel space) and ``affine`` (the
+  mapping), or
+* a length 2 sequence with the same information (shape, affine).
+"""
+
+from itertools import product
+
+import numpy as np
+
+from .affines import apply_affine
+
+
+def vox2out_vox(mapped_voxels, voxel_sizes=None):
+    """ output-aligned shape, affine for input implied by `mapped_voxels`
+
+    The input (voxel) space, and the affine mapping to output space, are given
+    in `mapped_voxels`.
+
+    The output space is implied by the affine, we don't need to know what that
+    is, we just return something with the same (implied) output space.
+
+    Our job is to work out another voxel space where the voxel array axes and
+    the output axes are aligned (top left 3 x 3 of affine is diagonal with all
+    positive entries) and which contains all the voxels of the implied input
+    image at their correct output space positions, once resampled into the
+    output voxel space.
+
+    Parameters
+    ----------
+    mapped_voxels : object or length 2 sequence
+        If object, has attributes ``shape`` giving input voxel shape, and
+        ``affine`` giving mapping of input voxels to output space. If length 2
+        sequence, elements are (shape, affine) with same meaning as above. The
+        affine is a (4, 4) array-like.
+    voxel_sizes : None or sequence
+        Gives the diagonal entries of `output_affine` (except the trailing 1
+        for the homogenous coordinates) (``output_affine == np.diag(voxel_sizes
+        + [1])``). If None, return identity `output_affine`.
+
+    Returns
+    -------
+    output_shape : sequence
+        Shape of output image that has voxel axes aligned to original image
+        output space axes, and encloses all the voxel data from the original
+        image implied by input shape.
+    output_affine : (4, 4) array
+        Affine of output image that has voxel axes aligned to the output axes
+        implied by input affine. Top-left 3 x 3 part of affine is diagonal with
+        all positive entries.  The entries come from `voxel_sizes` if
+        specified, or are all 1.  If the image is < 3D, then the missing
+        dimensions will have a 1 in the matching diagonal.
+    """
+    try:
+        in_shape, in_affine = mapped_voxels.shape, mapped_voxels.affine
+    except AttributeError:
+        in_shape, in_affine = mapped_voxels
+    n_axes = len(in_shape)
+    if n_axes > 3:
+        raise ValueError('This function can only deal with 3D images')
+    if n_axes < 3:
+        in_shape += (1,) * (3 - n_axes)
+    out_vox = np.ones((3,))
+    if not voxel_sizes is None:
+        if not len(voxel_sizes) == n_axes:
+            raise ValueError('voxel sizes length should match shape')
+        if not np.all(np.array(voxel_sizes) > 0):
+            raise ValueError('voxel sizes should all be positive')
+        out_vox[:n_axes] = voxel_sizes
+    in_mn_mx = zip([0, 0, 0], np.array(in_shape) - 1)
+    in_corners = list(product(*in_mn_mx))
+    out_corners = apply_affine(in_affine, in_corners)
+    out_mn = out_corners.min(axis=0)
+    out_mx = out_corners.max(axis=0)
+    out_shape = np.ceil((out_mx - out_mn) / out_vox) + 1
+    out_affine = np.diag(list(out_vox) + [1])
+    out_affine[:3, 3] = out_mn
+    return tuple(int(i) for i in out_shape[:n_axes]), out_affine
+
+
+def slice2volume(index, axis, shape=None):
+    """ Affine expressing selection of a single slice from 3D volume
+
+    Imagine we have taken a slice from an image data array, ``s = data[:, :,
+    index]``.  This function returns the affine to map the array coordinates of
+    ``s`` to the array coordinates of ``data``.
+
+    This can be useful for resampling a single slice from a volume.  For
+    example, to resample slice ``k`` in the space of ``img1`` from the matching
+    spatial voxel values in ``img2``, you might do something like::
+
+        slice_shape = img1.shape[:2]
+        slice_aff = slice2volume(k, 2)
+        whole_aff = np.linalg.inv(img2.affine).dot(img1.affine.dot(slice_aff))
+
+    and then use ``whole_aff`` in ``scipy.ndimage.affine_transform``:
+
+        rzs, trans = to_matvec(whole_aff)
+        data = img2.get_data()
+        new_slice = scipy.ndimage.affine_transform(data, rzs, trans, slice_shape)
+
+    Parameters
+    ----------
+    index : int
+        index of selected slice
+    axis : {0, 1, 2}
+        axis to which `index` applies
+
+    Returns
+    -------
+    slice_aff : shape (4, 3) affine
+        Affine relating input coordinates in a slice to output coordinates in
+        the embedded volume
+    """
+    if index < 0:
+        raise ValueError("Cannot handle negative index")
+    if not 0 <= axis <= 2:
+        raise ValueError("Axis should be between 0 and 2")
+    axes = list(range(4))
+    axes.remove(axis)
+    slice_aff = np.eye(4)[:, axes]
+    slice_aff[axis, -1] = index
+    return slice_aff

nibabel/tests/test_spaces.py

@@ -0,0 +1,113 @@
+""" Tests for spaces module
+"""
+
+import numpy as np
+import numpy.linalg as npl
+
+from ..spaces import vox2out_vox, slice2volume
+from ..affines import apply_affine, from_matvec
+from ..nifti1 import Nifti1Image
+from ..eulerangles import euler2mat
+
+
+from numpy.testing import (assert_almost_equal,
+                           assert_array_equal)
+
+from nose.tools import (assert_true, assert_false, assert_raises,
+                        assert_equal, assert_not_equal)
+
+
+
+def assert_all_in(in_shape, in_affine, out_shape, out_affine):
+    slices = tuple(slice(N) for N in in_shape)
+    n_axes = len(in_shape)
+    in_grid = np.mgrid[slices]
+    in_grid = np.rollaxis(in_grid, 0, n_axes + 1)
+    v2v = npl.inv(out_affine).dot(in_affine)
+    if n_axes < 3: # reduced dimensions case
+        new_v2v = np.eye(n_axes + 1)
+        new_v2v[:n_axes, :n_axes] = v2v[:n_axes, :n_axes]
+        new_v2v[:n_axes, -1] = v2v[:n_axes, -1]
+        v2v = new_v2v
+    out_grid = apply_affine(v2v, in_grid)
+    TINY = 1e-12
+    assert_true(np.all(out_grid > -TINY))
+    assert_true(np.all(out_grid < np.array(out_shape) + TINY))
+
+
+def test_vox2out_vox():
+    # Test world space bounding box
+    # Test basic case, identity, no voxel sizes passed
+    shape, aff = vox2out_vox(((2, 3, 4), np.eye(4)))
+    assert_array_equal(shape, (2, 3, 4))
+    assert_array_equal(aff, np.eye(4))
+    # Some affines as input to the tests
+    trans_123 = [[1, 0, 0, 1], [0, 1, 0, 2], [0, 0, 1, 3], [0, 0, 0, 1]]
+    trans_m123 = [[1, 0, 0, -1], [0, 1, 0, -2], [0, 0, 1, -3], [0, 0, 0, 1]]
+    rot_3 = from_matvec(euler2mat(np.pi / 4), [0, 0, 0])
+    for in_shape, in_aff, vox, out_shape, out_aff in (
+        # Identity
+        ((2, 3, 4), np.eye(4), None, (2, 3, 4), np.eye(4)),
+        # Flip first axis
+        ((2, 3, 4), np.diag([-1, 1, 1, 1]), None,
+         (2, 3, 4), [[1, 0, 0, -1], # axis reversed -> -ve offset
+                     [0, 1, 0, 0],
+                     [0, 0, 1, 0],
+                     [0, 0, 0, 1]]),
+        # zooms for affine > 1 -> larger grid with default 1mm output voxels
+        ((2, 3, 4), np.diag([4, 5, 6, 1]), None,
+         (5, 11, 19), np.eye(4)),
+        # set output voxels to be same size as input. back to original shape
+        ((2, 3, 4), np.diag([4, 5, 6, 1]), (4, 5, 6),
+         (2, 3, 4), np.diag([4, 5, 6, 1])),
+        # Translation preserved in output
+        ((2, 3, 4), trans_123, None,
+         (2, 3, 4), trans_123),
+        ((2, 3, 4), trans_m123, None,
+         (2, 3, 4), trans_m123),
+        # rotation around 3rd axis
+        ((2, 3, 4), rot_3, None,
+         # x diff, y diff now 3 cos pi / 4 == 2.12, ceil to 3, add 1
+         # most negative x now 2 cos pi / 4
+         (4, 4, 4), [[1, 0, 0, -2 * np.cos(np.pi / 4)],
+                     [0, 1, 0, 0],
+                     [0, 0, 1, 0],
+                     [0, 0, 0, 1]]),
+        # Less than 3 axes
+        ((2, 3), np.eye(4), None,
+         (2, 3),  np.eye(4)),
+        ((2,), np.eye(4), None,
+         (2,),  np.eye(4)),
+        # Number of voxel sizes matches length
+        ((2, 3), np.diag([4, 5, 6, 1]), (4, 5),
+         (2, 3), np.diag([4, 5, 1, 1])),
+    ):
+        img = Nifti1Image(np.ones(in_shape), in_aff)
+        for input in ((in_shape, in_aff), img):
+            shape, aff = vox2out_vox(input, vox)
+            assert_all_in(in_shape, in_aff, shape, aff)
+            assert_equal(shape, out_shape)
+            assert_almost_equal(aff, out_aff)
+            assert_true(isinstance(shape, tuple))
+            assert_true(isinstance(shape[0], int))
+    # Enforce number of axes
+    assert_raises(ValueError, vox2out_vox, ((2, 3, 4, 5), np.eye(4)))
+    assert_raises(ValueError, vox2out_vox, ((2, 3, 4, 5, 6), np.eye(4)))
+    # Voxel sizes must be positive
+    assert_raises(ValueError, vox2out_vox, ((2, 3, 4), np.eye(4), [-1, 1, 1]))
+    assert_raises(ValueError, vox2out_vox, ((2, 3, 4), np.eye(4), [1, 0, 1]))
+
+
+def test_slice2volume():
+    # Get affine expressing selection of single slice from volume
+    for axis, def_aff in zip((0, 1, 2), (
+        [[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1]],
+        [[1, 0, 0], [0, 0, 0], [0, 1, 0], [0, 0, 1]],
+        [[1, 0, 0], [0, 1, 0], [0, 0, 0], [0, 0, 1]])):
+        for val in (0, 5, 10):
+            exp_aff = np.array(def_aff)
+            exp_aff[axis, -1] = val
+            assert_array_equal(slice2volume(val, axis), exp_aff)
+    assert_raises(ValueError, slice2volume, -1, 0)
+    assert_raises(ValueError, slice2volume, 0, -1)
+    assert_raises(ValueError, slice2volume, 0, 3)