Merge branch 'enh/nonlinear' into 'master'

Enh/nonlinear

See merge request fsl/fslpy!211

CHANGELOG.rst

@@ -15,6 +15,8 @@ Added
 * New ``firstDot`` option to the :func:`.path.getExt`,
   :func:`.path.removeExt`, and :func:`.path.splitExt`, functions, offering
   rudimentary support for double-barrelled filenames.
+* The :func:`.nonlinear.applyDeformation` function now accepts a ``premat``
+  affine, which is applied to the input image before the deformation field.
 
 
 Changed

fsl/transform/nonlinear.py

@@ -679,7 +679,13 @@ def convertDeformationSpace(field, from_, to):
         defType=field.deformationType)
 
 
-def applyDeformation(image, field, ref=None, order=1, mode=None, cval=None):
+def applyDeformation(image,
+                     field,
+                     ref=None,
+                     order=1,
+                     mode=None,
+                     cval=None,
+                     premat=None):
     """Applies a :class:`DeformationField` to an :class:`.Image`.
 
     The image is transformed into the space of the field's reference image
@@ -691,25 +697,30 @@ def applyDeformation(image, field, ref=None, order=1, mode=None, cval=None):
     the input image. It is therefore assumed that an alternate ``ref`` is
     aligned in world coordinates with the field's actual reference image.
 
-    :arg image: :class:`.Image` to be transformed
+    :arg image:  :class:`.Image` to be transformed
 
-    :arg field: :class:`DeformationField` to use
+    :arg field:  :class:`DeformationField` to use
 
-    :arg ref:   Alternate reference image - if not provided, ``field.ref``
-                is used
+    :arg ref:    Alternate reference image - if not provided, ``field.ref``
+                 is used
 
-    :arg order: Spline interpolation order, passed through to the
-                ``scipy.ndimage.affine_transform`` function - ``0``
-                corresponds to nearest neighbour interpolation, ``1``
-                (the default) to linear interpolation, and ``3`` to
-                cubic interpolation.
+    :arg order:  Spline interpolation order, passed through to the
+                 ``scipy.ndimage.affine_transform`` function - ``0``
+                 corresponds to nearest neighbour interpolation, ``1``
+                 (the default) to linear interpolation, and ``3`` to
+                 cubic interpolation.
 
-    :arg mode:  How to handle regions which are outside of the image FOV.
-                Defaults to `''nearest'``.
+    :arg mode:   How to handle regions which are outside of the image FOV.
+                 Defaults to `''nearest'``.
 
-    :arg cval:  Constant value to use when ``mode='constant'``.
+    :arg cval:   Constant value to use when ``mode='constant'``.
 
-    :return:    ``numpy.array`` containing the transformed image data.
+    :arg premat: Optional affine transform which can be used if ``image`` is
+                 not in the same space as ``field.src``. Assumed to transform
+                 from ``image`` **voxel** coordinates into ``field.src``
+                 **voxel** coordinates.
+
+    :return:     ``numpy.array`` containing the transformed image data.
     """
 
     if order is None: order = 1
@@ -763,12 +774,15 @@ def applyDeformation(image, field, ref=None, order=1, mode=None, cval=None):
     # We assume world-world alignment
     # between the original source
     # and the image to be resampled
-    if not image.sameSpace(src):
-        post  = affine.concat(image.getAffine('world', 'voxel'),
-                              src  .getAffine('voxel', 'world'))
+    if (premat is not None) or (not image.sameSpace(src)):
+        if premat is None:
+            premat = affine.concat(image.getAffine('world', 'voxel'),
+                                   src  .getAffine('voxel', 'world'))
+        else:
+            premat = affine.invert(premat)
         shape = field.shape
         field = field.reshape((-1, 3))
-        field = affine.transform(field, post)
+        field = affine.transform(field, premat)
         field = field.reshape(shape)
 
     field = field.transpose((3, 0, 1, 2))

tests/test_transform/test_nonlinear.py

@@ -411,6 +411,60 @@ def test_applyDeformation_altsrc():
     assert np.all(np.isclose(expect, result))
 
 
+def test_applyDeformation_premat():
+
+    src2ref = affine.compose(
+        np.random.randint(2, 5, 3),
+        np.random.randint(1, 10, 3),
+        [0, 0, 0])
+    ref2src = affine.invert(src2ref)
+
+    srcdata = np.random.randint(1, 65536, (10, 10, 10))
+    refdata = np.random.randint(1, 65536, (10, 10, 10))
+
+    src   = fslimage.Image(srcdata)
+    ref   = fslimage.Image(refdata, xform=src2ref)
+    field = _affine_field(src, ref, ref2src, 'world', 'world')
+
+    # First try a down-sampled version
+    # of the original source image
+    altsrc, xf = resample.resample(src, (5, 5, 5), origin='corner')
+    altsrc     = fslimage.Image(altsrc, xform=xf, header=src.header)
+    expect, xf = resample.resampleToReference(
+        altsrc, ref, matrix=src2ref, order=1, mode='nearest')
+    premat = affine.concat(src   .getAffine('world', 'voxel'),
+                           altsrc.getAffine('voxel', 'world'))
+    result = nonlinear.applyDeformation(
+        altsrc, field, order=1, mode='nearest', premat=premat)
+    assert np.all(np.isclose(expect, result))
+
+    # Now try a down-sampled ROI
+    # of the original source image
+    altsrc     = roi.roi(src, [(2, 9), (2, 9), (2, 9)])
+    altsrc, xf = resample.resample(altsrc, (4, 4, 4))
+    altsrc     = fslimage.Image(altsrc, xform=xf, header=src.header)
+    expect, xf = resample.resampleToReference(
+        altsrc, ref, matrix=src2ref, order=1, mode='nearest')
+    premat = affine.concat(src   .getAffine('world', 'voxel'),
+                           altsrc.getAffine('voxel', 'world'))
+    result = nonlinear.applyDeformation(
+        altsrc, field, order=1, mode='nearest', premat=premat)
+    assert np.all(np.isclose(expect, result))
+
+    # down-sampled and offset ROI
+    # of the original source image
+    altsrc     = roi.roi(src, [(-5, 8), (-5, 8), (-5, 8)])
+    altsrc, xf = resample.resample(altsrc, (6, 6, 6))
+    altsrc     = fslimage.Image(altsrc, xform=xf, header=src.header)
+    expect, xf = resample.resampleToReference(
+        altsrc, ref, matrix=src2ref, order=1, mode='nearest')
+    premat = affine.concat(src   .getAffine('world', 'voxel'),
+                           altsrc.getAffine('voxel', 'world'))
+    result = nonlinear.applyDeformation(
+        altsrc, field, order=1, mode='nearest', premat=premat)
+    assert np.all(np.isclose(expect, result))
+
+
 def test_applyDeformation_altref():
     src2ref = affine.compose(
         np.random.randint(2, 5, 3),