MAINT: Add new test to check DisplacementsFieldTransforms

Author: oesteban

nitransforms/__init__.py

@@ -17,7 +17,7 @@
    transform
 """
 from .linear import Affine
-from .nonlinear import DeformationFieldTransform
+from .nonlinear import DisplacementsFieldTransform
 
 
-__all__ = ['Affine', 'DeformationFieldTransform']
+__all__ = ['Affine', 'DisplacementsFieldTransform']

nitransforms/base.py

@@ -9,6 +9,7 @@
 """Common interface for transforms."""
 import numpy as np
 import h5py
+from nibabel.loadsave import load
 
 from scipy import ndimage as ndi
 
@@ -168,6 +169,9 @@ def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True,
             The moving imaged after resampling to reference space.
 
         """
+        if isinstance(moving, str):
+            moving = load(moving)
+
         moving_data = np.asanyarray(moving.dataobj)
         if output_dtype is None:
             output_dtype = moving_data.dtype

nitransforms/nonlinear.py

@@ -9,7 +9,7 @@
 """Nonlinear transforms."""
 from warnings import warn
 import numpy as np
-from scipy import ndimage as ndi
+# from scipy import ndimage as ndi
 # from gridbspline.maths import cubic
 
 from .base import ImageGrid, TransformBase
@@ -18,14 +18,14 @@
 # vbspl = np.vectorize(cubic)
 
 
-class DeformationFieldTransform(TransformBase):
+class DisplacementsFieldTransform(TransformBase):
     """Represents a dense field of displacements (one vector per voxel)."""
 
     __slots__ = ['_field']
 
     def __init__(self, field, reference=None):
         """Create a dense deformation field transform."""
-        super(DeformationFieldTransform, self).__init__()
+        super(DisplacementsFieldTransform, self).__init__()
         self._field = np.asanyarray(field.dataobj)
 
         ndim = self._field.ndim - 1
@@ -38,7 +38,8 @@ def __init__(self, field, reference=None):
         # displacement vector per voxel in output (reference)
         # space
         if reference is None:
-            reference = four_to_three(field)[0]
+            reference = field.__class__(np.zeros(self._field.shape[:-1]),
+                                        field.affine, field.header)
         elif reference.shape[:ndim] != field.shape[:ndim]:
             raise ValueError(
                 'Reference ({}) and field ({}) must have the same '
@@ -71,7 +72,7 @@ def map(self, x, inverse=False, index=0):
         >>> field = np.zeros((10, 10, 10, 3))
         >>> field[..., 0] = 4.0
         >>> fieldimg = nb.Nifti1Image(field, np.diag([2., 2., 2., 1.]))
-        >>> xfm = DeformationFieldTransform(fieldimg)
+        >>> xfm = DisplacementsFieldTransform(fieldimg)
         >>> xfm([4.0, 4.0, 4.0]).tolist()
         [[8.0, 4.0, 4.0]]
 

nitransforms/tests/data/ds-005_sub-01_T1w.nii.gz

@@ -1,6 +1,6 @@
 # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
 # vi: set ft=python sts=4 ts=4 sw=4 et:
-"""Tests of the transform module."""
+"""Tests of linear transforms."""
 import os
 import pytest
 import numpy as np

nitransforms/tests/data/ds-005_sub-01_from-OASIS_to-T1_warp.nii.gz

@@ -0,0 +1,54 @@
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+"""Tests of nonlinear transforms."""
+import os
+import shutil
+from subprocess import check_call
+import pytest
+
+import numpy as np
+import nibabel as nb
+from ..nonlinear import DisplacementsFieldTransform
+
+TESTS_BORDER_TOLERANCE = 0.05
+APPLY_NONLINEAR_CMD = {
+    'itk': """\
+antsApplyTransforms -d 3 -r {reference} -i {moving} \
+-o resampled.nii.gz -n NearestNeighbor -t {transform} --float\
+""".format,
+}
+
+
+@pytest.mark.parametrize('sw_tool', ['itk'])
+def test_displacements_field(tmp_path, data_path, sw_tool):
+    os.chdir(str(tmp_path))
+    img_fname = os.path.join(data_path, 'tpl-OASIS30ANTs_T1w.nii.gz')
+    xfm_fname = os.path.join(
+        data_path, 'ds-005_sub-01_from-OASIS_to-T1_warp.nii.gz')
+    ants_warp = nb.load(xfm_fname)
+    field = nb.Nifti1Image(
+        np.squeeze(np.asanyarray(ants_warp.dataobj)),
+        ants_warp.affine, ants_warp.header
+    )
+
+    xfm = DisplacementsFieldTransform(field)
+
+    # Then apply the transform and cross-check with software
+    cmd = APPLY_NONLINEAR_CMD[sw_tool](
+        transform=os.path.abspath(xfm_fname),
+        reference=img_fname,
+        moving=img_fname)
+
+    # skip test if command is not available on host
+    exe = cmd.split(" ", 1)[0]
+    if not shutil.which(exe):
+        pytest.skip("Command {} not found on host".format(exe))
+
+    exit_code = check_call([cmd], shell=True)
+    assert exit_code == 0
+    sw_moved = nb.load('resampled.nii.gz')
+
+    nt_moved = xfm.resample(img_fname, order=0)
+    diff = sw_moved.get_fdata() - nt_moved.get_fdata()
+    # A certain tolerance is necessary because of resampling at borders
+    assert (np.abs(diff) > 1e-3).sum() / diff.size < TESTS_BORDER_TOLERANCE