FIX: Double-check dtypes within tests and increase RMSE tolerance

nitransforms/conftest.py

@@ -7,6 +7,7 @@
 import tempfile
 
 _data = None
+_brainmask = None
 _testdir = Path(os.getenv("TEST_DATA_HOME", "~/.nitransforms/testdata")).expanduser()
 
 
@@ -51,29 +52,47 @@ def get_testdata():
     if _data is not None:
         return _data
 
-    img = nb.load(_testdir / "someones_anatomy.nii.gz")
-    imgaff = img.affine
+    return _reorient(_testdir / "someones_anatomy.nii.gz")
+
+
+@pytest.fixture
+def get_testmask():
+    """Generate data in the requested orientation."""
+    global _brainmask
 
+    if _brainmask is not None:
+        return _brainmask
+
+    return _reorient(_testdir / "someones_anatomy_brainmask.nii.gz")
+
+
+def _reorient(path):
+    """Reorient the input NIfTI file."""
+    img = nb.load(path)
+    imgaff = img.affine
     _data = {"RAS": img}
     newaff = imgaff.copy()
     newaff[0, 0] *= -1.0
     newaff[0, 3] = imgaff.dot(np.hstack((np.array(img.shape[:3]) - 1, 1.0)))[0]
-    _data["LAS"] = nb.Nifti1Image(np.flip(img.get_fdata(), 0), newaff, img.header)
+    _data["LAS"] = nb.Nifti1Image(np.flip(np.asanyarray(img.dataobj), 0), newaff, img.header)
+    _data["LAS"].set_data_dtype(img.get_data_dtype())
     newaff = imgaff.copy()
     newaff[0, 0] *= -1.0
     newaff[1, 1] *= -1.0
     newaff[:2, 3] = imgaff.dot(np.hstack((np.array(img.shape[:3]) - 1, 1.0)))[:2]
     _data["LPS"] = nb.Nifti1Image(
-        np.flip(np.flip(img.get_fdata(), 0), 1), newaff, img.header
+        np.flip(np.flip(np.asanyarray(img.dataobj), 0), 1), newaff, img.header
     )
+    _data["LPS"].set_data_dtype(img.get_data_dtype())
     A = nb.volumeutils.shape_zoom_affine(
         img.shape, img.header.get_zooms(), x_flip=False
     )
     R = nb.affines.from_matvec(nb.eulerangles.euler2mat(x=0.09, y=0.001, z=0.001))
     newaff = R.dot(A)
-    oblique_img = nb.Nifti1Image(img.get_fdata(), newaff, img.header)
+    oblique_img = nb.Nifti1Image(np.asanyarray(img.dataobj), newaff, img.header)
     oblique_img.header.set_qform(newaff, 1)
     oblique_img.header.set_sform(newaff, 1)
     _data["oblique"] = oblique_img
+    _data["oblique"].set_data_dtype(img.get_data_dtype())
 
     return _data

nitransforms/tests/test_linear.py

@@ -14,23 +14,23 @@
 from .. import linear as nitl
 from .utils import assert_affines_by_filename
 
-TESTS_BORDER_TOLERANCE = 0.05
+RMSE_TOL = 0.1
 APPLY_LINEAR_CMD = {
     "fsl": """\
 flirt -setbackground 0 -interp nearestneighbour -in {moving} -ref {reference} \
--applyxfm -init {transform} -out resampled.nii.gz\
+-applyxfm -init {transform} -out {resampled}\
 """.format,
     "itk": """\
 antsApplyTransforms -d 3 -r {reference} -i {moving} \
--o resampled.nii.gz -n NearestNeighbor -t {transform} --float\
+-o {resampled} -n NearestNeighbor -t {transform} --float\
 """.format,
     "afni": """\
 3dAllineate -base {reference} -input {moving} \
--prefix resampled.nii.gz -1Dmatrix_apply {transform} -final NN\
+-prefix {resampled} -1Dmatrix_apply {transform} -final NN\
 """.format,
     "fs": """\
 mri_vol2vol --mov {moving} --targ {reference} --lta {transform} \
---o resampled.nii.gz --nearest""".format,
+--o {resampled} --nearest""".format,
 }
 
 
@@ -121,13 +121,15 @@ def test_linear_save(tmpdir, data_path, get_testdata, image_orientation, sw_tool
     assert_affines_by_filename(xfm_fname1, xfm_fname2)
 
 
-@pytest.mark.parametrize("image_orientation", ["RAS", "LAS", "LPS",])  # 'oblique',
+@pytest.mark.parametrize("image_orientation", ["RAS", "LAS", "LPS", ])  # 'oblique',
 @pytest.mark.parametrize("sw_tool", ["itk", "fsl", "afni", "fs"])
-def test_apply_linear_transform(tmpdir, get_testdata, image_orientation, sw_tool):
+def test_apply_linear_transform(tmpdir, get_testdata, get_testmask, image_orientation, sw_tool):
     """Check implementation of exporting affines to formats."""
     tmpdir.chdir()
 
     img = get_testdata[image_orientation]
+    msk = get_testmask[image_orientation]
+
     # Generate test transform
     T = from_matvec(euler2mat(x=0.9, y=0.001, z=0.001), [4.0, 2.0, -1.0])
     xfm = nitl.Affine(T)
@@ -140,33 +142,60 @@ def test_apply_linear_transform(tmpdir, get_testdata, image_orientation, sw_tool
         ext = ".lta"
 
     img.to_filename("img.nii.gz")
+    msk.to_filename("mask.nii.gz")
+
+    # Write out transform file (software-dependent)
     xfm_fname = "M.%s%s" % (sw_tool, ext)
     xfm.to_filename(xfm_fname, fmt=sw_tool)
 
     cmd = APPLY_LINEAR_CMD[sw_tool](
         transform=os.path.abspath(xfm_fname),
-        reference=os.path.abspath("img.nii.gz"),
-        moving=os.path.abspath("img.nii.gz"),
+        reference=os.path.abspath("mask.nii.gz"),
+        moving=os.path.abspath("mask.nii.gz"),
+        resampled=os.path.abspath("resampled_brainmask.nii.gz"),
     )
 
     # 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))
 
+    # resample mask
+    exit_code = check_call([cmd], shell=True)
+    assert exit_code == 0
+    sw_moved_mask = nb.load("resampled_brainmask.nii.gz")
+    nt_moved_mask = xfm.apply(msk, order=0)
+    nt_moved_mask.set_data_dtype(msk.get_data_dtype())
+    diff = np.asanyarray(sw_moved_mask.dataobj) - np.asanyarray(nt_moved_mask.dataobj)
+
+    assert np.sqrt((diff ** 2).mean()) < RMSE_TOL
+
+    cmd = APPLY_LINEAR_CMD[sw_tool](
+        transform=os.path.abspath(xfm_fname),
+        reference=os.path.abspath("img.nii.gz"),
+        moving=os.path.abspath("img.nii.gz"),
+        resampled=os.path.abspath("resampled.nii.gz"),
+    )
+
+    brainmask = np.asanyarray(nt_moved_mask.dataobj, dtype=bool)
+
     exit_code = check_call([cmd], shell=True)
     assert exit_code == 0
     sw_moved = nb.load("resampled.nii.gz")
+    sw_moved.set_data_dtype(img.get_data_dtype())
 
     nt_moved = xfm.apply(img, order=0)
-    diff = sw_moved.get_fdata() - nt_moved.get_fdata()
+    diff = (sw_moved.get_fdata() - nt_moved.get_fdata())
+    diff[~brainmask] = 0.0
+    diff[np.abs(diff) < 1e-3] = 0
+
     # A certain tolerance is necessary because of resampling at borders
-    assert (np.abs(diff) > 1e-3).sum() / diff.size < TESTS_BORDER_TOLERANCE
+    assert np.sqrt((diff ** 2).mean()) < RMSE_TOL
 
     nt_moved = xfm.apply("img.nii.gz", 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
+    assert np.sqrt((diff[brainmask] ** 2).mean()) < RMSE_TOL
 
 
 def test_Affine_to_x5(tmpdir, testdata_path):