TEST: Tests using alternate source image

tests/test_scripts/test_fsl_apply_x5.py

@@ -6,6 +6,7 @@ import numpy as np
 import fsl.scripts.fsl_apply_x5 as fsl_apply_x5
 import fsl.data.image           as fslimage
 import fsl.utils.image.resample as resample
+import fsl.utils.image.roi      as roi
 import fsl.transform.x5         as x5
 import fsl.transform.affine     as affine
 import fsl.utils.tempdir        as tempdir
@@ -20,9 +21,13 @@ def _random_affine():
         np.random.random(3))
 
 
-def _random_image(aff):
-    vx, vy, vz = np.random.randint(10, 50, 3)
-    data       = (np.random.random((vx, vy, vz)) - 0.5) * 10
+def _random_image(aff, shape=None):
+
+    if shape is None:
+        shape = np.random.randint(10, 50, 3)
+
+    data = (np.random.random(shape) - 0.5) * 10
+
     return fslimage.Image(data, xform=aff)
 
 
@@ -129,3 +134,118 @@ def test_nonlinear_altref():
 
         assert result.sameSpace(altref)
         assert np.all(result.data == expect)
+
+
+def test_linear_altsrc():
+    with tempdir.tempdir():
+
+        src2ref = _random_affine()
+        src     = _random_image(np.eye(4), (20, 20, 20))
+        ref     = _random_image(src2ref)
+
+        x5.writeLinearX5('xform.x5', src2ref, src, ref)
+
+        src.save('src')
+        ref.save('ref')
+
+        srclo, xf = resample.resample(src, (10, 10, 10))
+        srclo = fslimage.Image(srclo, xform=xf)
+        srchi, xf = resample.resample(src, (40, 40, 40))
+        srchi = fslimage.Image(srchi, xform=xf)
+
+        srcoff = roi.roi(src, [(-10, 10), (-10, 10), (-10, 10)])
+
+        srclo .save('srclo')
+        srchi .save('srchi')
+        srcoff.save('srcoff')
+
+        fsl_apply_x5.main('src    xform.x5 out'   .split())
+        fsl_apply_x5.main('srclo  xform.x5 outlo' .split())
+        fsl_apply_x5.main('srchi  xform.x5 outhi' .split())
+        fsl_apply_x5.main('srcoff xform.x5 outoff'.split())
+
+        out    = fslimage.Image('out')
+        outlo  = fslimage.Image('outlo')
+        outhi  = fslimage.Image('outhi')
+        outoff = fslimage.Image('outoff')
+
+        exp    = resample.resampleToReference(src,    ref, matrix=src2ref)[0]
+        explo  = resample.resampleToReference(srclo,  ref, matrix=src2ref)[0]
+        exphi  = resample.resampleToReference(srchi,  ref, matrix=src2ref)[0]
+        expoff = resample.resampleToReference(srcoff, ref, matrix=src2ref)[0]
+
+        assert out   .sameSpace(ref)
+        assert outlo .sameSpace(ref)
+        assert outhi .sameSpace(ref)
+        assert outoff.sameSpace(ref)
+
+        assert np.all(np.isclose(out   .data, exp))
+        assert np.all(np.isclose(outlo .data, explo))
+        assert np.all(np.isclose(outhi .data, exphi))
+        assert np.all(np.isclose(outoff.data, expoff))
+
+
+def test_nonlinear_altsrc():
+    with tempdir.tempdir():
+
+        src2ref = _random_affine()
+        ref2src = affine.invert(src2ref)
+        src     = _random_image(np.eye(4), (20, 20, 20))
+        ref     = _random_image(src2ref,   (20, 20, 20))
+
+        field   = _affine_field(src, ref, ref2src, 'world', 'world')
+
+        x5.writeNonLinearX5('xform.x5', field)
+
+        src.save('src')
+        ref.save('ref')
+
+        srclo, xf = resample.resample(src, (10, 10, 10), origin='corner')
+        srclo = fslimage.Image(srclo, xform=xf)
+        srchi, xf = resample.resample(src, (40, 40, 40), origin='corner')
+        srchi = fslimage.Image(srchi, xform=xf)
+
+        srcoff = roi.roi(src, [(-10, 10), (-10, 10), (-10, 10)])
+
+        srclo .save('srclo')
+        srchi .save('srchi')
+        srcoff.save('srcoff')
+
+        # The up-sampled source image is subject to
+        # an unavoidable interpolation, as the
+        # deformation field coordinates are affine-
+        # transformed to the space of the alternate
+        # source image. So in this test case we use
+        # nearest-neighbour interp, so that we can
+        # get the same output as a standard linear
+        # resample
+        fsl_apply_x5.main('src    xform.x5 out'             .split())
+        fsl_apply_x5.main('srclo  xform.x5 outlo'           .split())
+        fsl_apply_x5.main('srchi  xform.x5 outhi -i nearest'.split())
+        fsl_apply_x5.main('srcoff xform.x5 outoff'          .split())
+
+        out    = fslimage.Image('out')
+        outlo  = fslimage.Image('outlo')
+        outhi  = fslimage.Image('outhi')
+        outoff = fslimage.Image('outoff')
+
+        exp,    x1 = resample.resampleToReference(src,    ref, matrix=src2ref, mode='constant')
+        explo,  x2 = resample.resampleToReference(srclo,  ref, matrix=src2ref, mode='constant')
+        exphi,  x3 = resample.resampleToReference(srchi,  ref, matrix=src2ref, mode='constant', order=0)
+        expoff, x4 = resample.resampleToReference(srcoff, ref, matrix=src2ref, mode='constant')
+
+        assert out   .sameSpace(ref)
+        assert outlo .sameSpace(ref)
+        assert outhi .sameSpace(ref)
+        assert outoff.sameSpace(ref)
+
+        # We get boundary issues just at the first
+        # voxel, so I'm masking that voxel out
+        for img in (out, outlo, outhi, outoff,
+                    exp, explo, exphi, expoff):
+            img[0, 0, 0] = 0
+
+        assert np.all(np.isclose(out   .data, exp))
+        assert np.all(np.isclose(outlo .data, explo))
+        assert np.all(np.isclose(outhi .data, exphi))
+        assert np.all(np.isclose(outoff.data, expoff))

tests/test_transform/test_nonlinear.py

@@ -7,6 +7,7 @@ import numpy   as np
 
 import fsl.data.image           as fslimage
 import fsl.utils.image.resample as resample
+import fsl.utils.image.roi      as roi
 import fsl.transform.affine     as affine
 import fsl.transform.nonlinear  as nonlinear
 import fsl.transform.fnirt      as fnirt
@@ -362,6 +363,53 @@ def test_applyDeformation():
     assert np.all(np.isclose(expect, result))
 
 
+def test_applyDeformation_altsrc():
+
+    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')
+    result = nonlinear.applyDeformation(
+        altsrc, field, order=1, mode='nearest')
+    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')
+    result = nonlinear.applyDeformation(
+        altsrc, field, order=1, mode='nearest')
+    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')
+    result = nonlinear.applyDeformation(
+        altsrc, field, order=1, mode='nearest')
+    assert np.all(np.isclose(expect, result))
+
 
 def test_applyDeformation_altref():
     src2ref = affine.compose(
@@ -427,8 +475,6 @@ def test_applyDeformation_worldAligned():
     result = nonlinear.applyDeformation(
         src, field, order=1, mode='constant', cval=0)
 
-    fslimage.Image(result, xform=ref.voxToWorldMat).save('result')
-
     expect = expect[1:-1, 1:-1, 1:-1]
     result = result[1:-1, 1:-1, 1:-1]