add a deformation field transform

nibabel/transform/__init__.py

@@ -17,3 +17,4 @@
 """
 from __future__ import absolute_import
 from .base import ImageSpace, Affine
+from .nonlinear import DeformationFieldTransform

nibabel/transform/base.py

@@ -51,9 +51,60 @@ def reference(self):
     def reference(self, image):
         self._reference = ImageSpace(image)
 
-    def resample(self, moving, order=3, dtype=None):
-        '''A virtual method to resample the moving image in the reference space'''
-        raise NotImplementedError
+    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.
+
+        '''
+
+        if output_dtype is None:
+            output_dtype = moving.header.get_data_dtype()
+
+        moving_data = moving.get_data()
+        moved = ndi.geometric_transform(
+            moving_data,
+            mapping=self.map_voxel,
+            output_shape=self.reference.shape,
+            output=output_dtype,
+            order=order,
+            mode=mode,
+            cval=cval,
+            prefilter=prefilter,
+            extra_keywords={'moving': moving},
+        )
+
+        moved_image = moving.__class__(moved, self.reference.affine, moving.header)
+        moved_image.header.set_data_dtype(output_dtype)
+        return moved_image
 
     def map_point(self, coords):
         '''Find the coordinates in moving space corresponding to the
@@ -164,10 +215,10 @@ def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True,
             output_shape=reference.shape, order=order, mode=mode,
             cval=cval, prefilter=prefilter)
 
-        img = moving.__class__(moved, moving.affine, moving.header)
-        img.header.set_data_dtype(output_dtype)
+        moved_image = moving.__class__(moved, reference.affine, moving.header)
+        moved_image.header.set_data_dtype(output_dtype)
 
-        return img
+        return moved_image
 
     def map_point(self, coords, forward=True):
         coords = np.array(coords)

nibabel/transform/nonlinear.py

@@ -0,0 +1,55 @@
+# 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.
+#
+### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
+''' Common interface for transforms '''
+from __future__ import division, print_function, absolute_import
+import numpy as np
+
+from .base import TransformBase
+from ..funcs import four_to_three
+
+
+class DeformationFieldTransform(TransformBase):
+    '''Represents linear transforms on image data'''
+    __slots__ = ['_field']
+
+    def __init__(self, field):
+        self.reference = four_to_three(field)[0]
+        # Set the field components in front
+        self._field = field.get_data()
+
+    def map_voxel(self, index, moving=None):
+        '''
+
+        Examples
+        --------
+
+        >>> import numpy as np
+        >>> import nibabel as nb
+        >>> ref = nb.load('t1_weighted.nii.gz')
+        >>> field = np.zeros(tuple(list(ref.shape) + [3]))
+        >>> field[..., 0] = 4.0
+        >>> fieldimg = nb.Nifti1Image(field, ref.affine, ref.header)
+        >>> xfm = nb.transform.DeformationFieldTransform(fieldimg)
+        >>> new = xfm.resample(ref)
+        >>> new.to_filename('deffield.nii.gz')
+
+        '''
+        vox2ras = self.reference.affine
+        ras2vox = np.linalg.inv(vox2ras if moving is None
+                                else moving.affine)
+
+        index = np.array(index)
+        if index.shape[0] == vox2ras.shape[0] - 1:
+            index = np.append(index, [1.0])
+
+        idx = tuple([int(i) for i in index[:-1]])
+        point = vox2ras.dot(index)[:-1]
+        delta = np.squeeze(self._field[idx + (slice(None), )])
+        newindex = ras2vox.dot(np.append(point + delta, [1]))[:-1]
+        return tuple(newindex)