enh: full implementation of the resample method

Author: oesteban

nitransforms/base.py

@@ -7,6 +7,7 @@
 #
 ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
 """Common interface for transforms."""
+from pathlib import Path
 import numpy as np
 import h5py
 from nibabel.loadsave import load
@@ -24,15 +25,16 @@ class ImageGrid(object):
 
     def __init__(self, image):
         """Create a gridded sampling reference."""
+        if isinstance(image, (str, Path)):
+            image = load(str(image))
+
         self._affine = image.affine
         self._shape = image.shape
         self._ndim = len(image.shape)
         self._nvox = np.prod(image.shape)  # Do not access data array
         self._ndindex = None
         self._coords = None
         self._inverse = np.linalg.inv(image.affine)
-        if self._ndim not in [2, 3]:
-            raise ValueError('Invalid image space (%d-D)' % self._ndim)
 
     @property
     def affine(self):
@@ -81,13 +83,11 @@ def ndcoords(self):
 
     def ras(self, ijk):
         """Get RAS+ coordinates from input indexes."""
-        ras = self._affine.dot(_as_homogeneous(ijk).T)[:3, ...]
-        return ras.T
+        return _apply_affine(ijk, self._affine, self._ndim)
 
     def index(self, x):
         """Get the image array's indexes corresponding to coordinates."""
-        ijk = self._inverse.dot(_as_homogeneous(x).T)[:3, ...]
-        return ijk.T
+        return _apply_affine(x, self._inverse, self._ndim)
 
     def _to_hdf5(self, group):
         group.attrs['Type'] = 'image'
@@ -173,22 +173,23 @@ def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True,
             moving = load(moving)
 
         moving_data = np.asanyarray(moving.dataobj)
-        if output_dtype is None:
-            output_dtype = moving_data.dtype
+        output_dtype = output_dtype or moving_data.dtype
+        targets = ImageGrid(moving).index(
+            _as_homogeneous(self.map(self.reference.ndcoords.T),
+                            dim=self.reference.ndim))
 
-        moved = ndi.geometric_transform(
+        moved = ndi.map_coordinates(
             moving_data,
-            mapping=self._map_index,
-            output_shape=self.reference.shape,
+            targets.T,
             output=output_dtype,
             order=order,
             mode=mode,
             cval=cval,
             prefilter=prefilter,
-            extra_keywords={'moving': ImageGrid(moving)},
         )
 
-        moved_image = moving.__class__(moved, self.reference.affine, moving.header)
+        moved_image = moving.__class__(moved.reshape(self.reference.shape),
+                                       self.reference.affine, moving.header)
         moved_image.header.set_data_dtype(output_dtype)
         return moved_image
 
@@ -213,10 +214,6 @@ def map(self, x, inverse=False, index=0):
         """
         raise NotImplementedError
 
-    def _map_index(self, ijk, moving):
-        x = self.reference.ras(_as_homogeneous(ijk))
-        return moving.index(self.map(x))
-
     def to_filename(self, filename, fmt='X5'):
         """Store the transform in BIDS-Transforms HDF5 file format (.x5)."""
         with h5py.File(filename, 'w') as out_file:
@@ -232,22 +229,33 @@ def _to_hdf5(self, x5_root):
         raise NotImplementedError
 
 
-def _as_homogeneous(xyz, dtype='float32'):
+def _as_homogeneous(xyz, dtype='float32', dim=3):
     """
     Convert 2D and 3D coordinates into homogeneous coordinates.
 
     Examples
     --------
-    >>> _as_homogeneous((4, 5), dtype='int8').tolist()
+    >>> _as_homogeneous((4, 5), dtype='int8', dim=2).tolist()
     [[4, 5, 1]]
 
     >>> _as_homogeneous((4, 5, 6),dtype='int8').tolist()
     [[4, 5, 6, 1]]
 
+    >>> _as_homogeneous((4, 5, 6, 1),dtype='int8').tolist()
+    [[4, 5, 6, 1]]
+
     >>> _as_homogeneous([(1, 2, 3), (4, 5, 6)]).tolist()
     [[1.0, 2.0, 3.0, 1.0], [4.0, 5.0, 6.0, 1.0]]
 
 
     """
     xyz = np.atleast_2d(np.array(xyz, dtype=dtype))
+    if np.shape(xyz)[-1] == dim + 1:
+        return xyz
+
     return np.hstack((xyz, np.ones((xyz.shape[0], 1), dtype=dtype)))
+
+
+def _apply_affine(x, affine, dim):
+    """Get the image array's indexes corresponding to coordinates."""
+    return affine.dot(_as_homogeneous(x, dim=dim).T)[:dim, ...].T

nitransforms/conftest.py

@@ -35,7 +35,7 @@ def doctest_autoimport(doctest_namespace):
 @pytest.fixture
 def data_path():
     """Return the test data folder."""
-    return os.path.join(os.path.dirname(__file__), 'tests/data')
+    return Path(__file__).parent / 'tests' / 'data'
 
 
 @pytest.fixture

nitransforms/linear.py

@@ -9,13 +9,11 @@
 """Linear transforms."""
 import sys
 import numpy as np
-from scipy import ndimage as ndi
 from pathlib import Path
-import warnings
 
 from nibabel.loadsave import load as loadimg
 from nibabel.affines import from_matvec, voxel_sizes, obliquity
-from .base import TransformBase
+from .base import TransformBase, _as_homogeneous
 from .patched import shape_zoom_affine
 from . import io
 
@@ -74,110 +72,6 @@ def matrix(self):
         """Access the internal representation of this affine."""
         return self._matrix
 
-    def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True,
-                 output_dtype=None):
-        """
-        Resample the moving image in reference space.
-
-        Parameters
-        ----------
-        moving : `spatialimage`
-            The image object containing the data to be resampled in reference
-            space
-        order : int, optional
-            The order of the spline interpolation, default is 3.
-            The order has to be in the range 0-5.
-        mode : {'reflect', 'constant', 'nearest', 'mirror', 'wrap'}, optional
-            Determines how the input image is extended when the resamplings overflows
-            a border. Default is 'constant'.
-        cval : float, optional
-            Constant value for ``mode='constant'``. Default is 0.0.
-        prefilter: bool, optional
-            Determines if the moving image's data array is prefiltered with
-            a spline filter before interpolation. The default is ``True``,
-            which will create a temporary *float64* array of filtered values
-            if *order > 1*. If setting this to ``False``, the output will be
-            slightly blurred if *order > 1*, unless the input is prefiltered,
-            i.e. it is the result of calling the spline filter on the original
-            input.
-
-        Returns
-        -------
-        moved_image : `spatialimage`
-            The moving imaged after resampling to reference space.
-
-        Examples
-        --------
-        >>> a = Affine()
-        >>> a.matrix
-        array([[[1., 0., 0., 0.],
-                [0., 1., 0., 0.],
-                [0., 0., 1., 0.],
-                [0., 0., 0., 1.]]])
-
-        >>> xfm = Affine([[1, 0, 0, 4], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]])
-        >>> ref = nb.load(testfile)
-        >>> xfm.reference = ref
-        >>> refdata = ref.get_fdata()
-        >>> np.allclose(refdata, 0)
-        True
-
-        >>> refdata[5, 5, 5] = 1  # Set a one in the middle voxel
-        >>> moving = nb.Nifti1Image(refdata, ref.affine, ref.header)
-        >>> resampled = xfm.resample(moving, order=0).get_fdata()
-        >>> resampled[1, 5, 5]
-        1.0
-
-        """
-        if output_dtype is None:
-            output_dtype = moving.header.get_data_dtype()
-
-        try:
-            reference = self.reference
-        except ValueError:
-            warnings.warn('No reference space defined, using moving as reference')
-            reference = moving
-
-        nvols = 1
-        if len(moving.shape) > 3:
-            nvols = moving.shape[3]
-
-        movaff = moving.affine
-        movingdata = np.asanyarray(moving.dataobj)
-        if nvols == 1:
-            movingdata = movingdata[..., np.newaxis]
-
-        nmats = self._matrix.shape[0]
-        if nvols != nmats and nmats > 1:
-            raise ValueError("""\
-The moving image contains {0} volumes, while the transform is defined for \
-{1} volumes""".format(nvols, nmats))
-
-        singlemat = None
-        if nmats == 1:
-            singlemat = np.linalg.inv(movaff).dot(self._matrix[0].dot(reference.affine))
-
-        if singlemat is not None and nvols > nmats:
-            warnings.warn('Resampling a 4D volume with a single affine matrix')
-
-        # Compose an index to index affine matrix
-        moved = []
-        for i in range(nvols):
-            i2imat = singlemat if singlemat is not None else np.linalg.inv(
-                movaff).dot(self._matrix[i].dot(reference.affine))
-            mdata = movingdata[..., i]
-            moved += [ndi.affine_transform(
-                mdata, matrix=i2imat[:-1, :-1],
-                offset=i2imat[:-1, -1],
-                output_shape=reference.shape, order=order, mode=mode,
-                cval=cval, prefilter=prefilter)]
-
-        moved_image = moving.__class__(
-            np.squeeze(np.stack(moved, -1)), reference.affine, moving.header)
-        moved_image.header.set_data_dtype(output_dtype)
-
-        return moved_image
-
     def map(self, x, inverse=False, index=0):
         r"""
         Apply :math:`y = f(x)`.
@@ -200,45 +94,17 @@ def map(self, x, inverse=False, index=0):
         --------
         >>> xfm = Affine([[1, 0, 0, 1], [0, 1, 0, 2], [0, 0, 1, 3], [0, 0, 0, 1]])
         >>> xfm.map((0,0,0))
-        array([1, 2, 3])
+        array([[1., 2., 3.]])
 
         >>> xfm.map((0,0,0), inverse=True)
-        array([-1., -2., -3.])
+        array([[-1., -2., -3.]])
 
         """
-        coords = np.array(x)
-        if coords.shape[0] == self._matrix[index].shape[0] - 1:
-            coords = np.append(coords, [1])
+        coords = _as_homogeneous(x, dim=self._matrix[0].shape[0] - 1).T
         affine = self._matrix[index]
-
         if inverse is True:
             affine = np.linalg.inv(self._matrix[index])
-
-        return affine.dot(coords)[:-1]
-
-    def _map_voxel(self, index, nindex=0, moving=None):
-        """Apply ijk' = f_ijk((i, j, k)), equivalent to the above with indexes."""
-        try:
-            reference = self.reference
-        except ValueError:
-            print('Warning: no reference space defined, using moving as reference',
-                  file=sys.stderr)
-            reference = moving
-        else:
-            if moving is None:
-                moving = reference
-        finally:
-            if reference is None:
-                raise ValueError('Reference and moving spaces are both undefined')
-
-        index = np.array(index)
-        if index.shape[0] == self._matrix[nindex].shape[0] - 1:
-            index = np.append(index, [1])
-
-        matrix = reference.affine.dot(
-            self._matrix[nindex].dot(np.linalg.inv(moving.affine))
-        )
-        return tuple(matrix.dot(index)[:-1])
+        return affine.dot(coords).T[..., :-1]
 
     def _to_hdf5(self, x5_root):
         """Serialize this object into the x5 file format."""