fix: revise tests by a human

Author: oesteban

env.yml

@@ -24,6 +24,8 @@ dependencies:
   - nitime=0.10
   - scikit-image=0.22
   - scikit-learn=1.4
+  # SimpleITK, so build doesn't complain about building scikit from sources
+  - simpleitk=2.2
   # Utilities
   - graphviz=9.0
   - pandoc=3.1

nitransforms/tests/test_io.py

@@ -1,6 +1,7 @@
 # emacs: -*- mode: python-mode; py-indent-offset: 4; indent-tabs-mode: nil -*-
 # vi: set ft=python sts=4 ts=4 sw=4 et:
 """I/O test cases."""
+
 import os
 from subprocess import check_call
 from io import StringIO
@@ -15,7 +16,7 @@
 from nibabel.affines import from_matvec
 from scipy.io import loadmat
 from nitransforms.linear import Affine
-from nitransforms import nonlinear as nitnl, linear as nitl
+from nitransforms import nonlinear as nitnl
 from nitransforms.io import (
     afni,
     fsl,
@@ -69,11 +70,13 @@ def test_volume_group_voxel_ordering():
 def test_VG_from_LTA(data_path):
     """Check the affine interpolation from volume geometries."""
     # affine manually clipped after running mri_info on the image
-    oracle = np.loadtxt(StringIO("""\
+    oracle = np.loadtxt(
+        StringIO("""\
 -3.0000   0.0000  -0.0000    91.3027
 -0.0000   2.0575  -2.9111   -25.5251
  0.0000   2.1833   2.7433  -105.0820
- 0.0000   0.0000   0.0000     1.0000"""))
+ 0.0000   0.0000   0.0000     1.0000""")
+    )
 
     lta_text = "\n".join(
         (data_path / "bold-to-t1w.lta").read_text().splitlines()[13:21]
@@ -420,10 +423,17 @@ def test_afni_Displacements():
 
 
 @pytest.mark.parametrize("only_linear", [True, False])
-@pytest.mark.parametrize("h5_path,nxforms", [
-    (_datadir / "affine-antsComposite.h5", 1),
-    (_testdir / "ds-005_sub-01_from-T1w_to-MNI152NLin2009cAsym_mode-image_xfm.h5", 2),
-])
+@pytest.mark.parametrize(
+    "h5_path,nxforms",
+    [
+        (_datadir / "affine-antsComposite.h5", 1),
+        (
+            _testdir
+            / "ds-005_sub-01_from-T1w_to-MNI152NLin2009cAsym_mode-image_xfm.h5",
+            2,
+        ),
+    ],
+)
 def test_itk_h5(tmpdir, only_linear, h5_path, nxforms):
     """Test displacements fields."""
     assert (
@@ -435,7 +445,9 @@ def test_itk_h5(tmpdir, only_linear, h5_path, nxforms):
                 )
             )
         )
-        == nxforms if not only_linear else 1
+        == nxforms
+        if not only_linear
+        else 1
     )
 
     with pytest.raises(TransformFileError):
@@ -466,24 +478,33 @@ def test_regressions(file_type, test_file, data_path):
     file_type.from_filename(data_path / "regressions" / test_file)
 
 
-@pytest.mark.parametrize("parameters", [
-    {"x": 0.1, "y": 0.03, "z": 0.002},
-    {"x": 0.001, "y": 0.3, "z": 0.002},
-    {"x": 0.01, "y": 0.03, "z": 0.2},
-])
+@pytest.mark.parametrize(
+    "parameters",
+    [
+        {"x": 0.1, "y": 0.03, "z": 0.002},
+        {"x": 0.001, "y": 0.3, "z": 0.002},
+        {"x": 0.01, "y": 0.03, "z": 0.2},
+    ],
+)
 @pytest.mark.parametrize("dir_x", (-1, 1))
 @pytest.mark.parametrize("dir_y", (-1, 1))
 @pytest.mark.parametrize("dir_z", (1, -1))
-@pytest.mark.parametrize("swapaxes", [
-    None, (0, 1), (1, 2), (0, 2),
-])
+@pytest.mark.parametrize(
+    "swapaxes",
+    [
+        None,
+        (0, 1),
+        (1, 2),
+        (0, 2),
+    ],
+)
 def test_afni_oblique(tmpdir, parameters, swapaxes, testdata_path, dir_x, dir_y, dir_z):
     tmpdir.chdir()
     img, R = _generate_reoriented(
         testdata_path / "someones_anatomy.nii.gz",
         (dir_x, dir_y, dir_z),
         swapaxes,
-        parameters
+        parameters,
     )
     img.to_filename("orig.nii.gz")
 
@@ -508,9 +529,8 @@ def test_afni_oblique(tmpdir, parameters, swapaxes, testdata_path, dir_x, dir_y,
         "orig.nii.gz",
     )
 
-    diff = (
-        np.asanyarray(img.dataobj, dtype="uint8")
-        - np.asanyarray(nt3drefit.dataobj, dtype="uint8")
+    diff = np.asanyarray(img.dataobj, dtype="uint8") - np.asanyarray(
+        nt3drefit.dataobj, dtype="uint8"
     )
     assert np.sqrt((diff[10:-10, 10:-10, 10:-10] ** 2).mean()) < 0.1
 
@@ -523,14 +543,15 @@ def test_afni_oblique(tmpdir, parameters, swapaxes, testdata_path, dir_x, dir_y,
         "deob_3drefit.nii.gz",
     )
 
-    diff = (
-        np.asanyarray(img.dataobj, dtype="uint8")
-        - np.asanyarray(nt_undo3drefit.dataobj, dtype="uint8")
+    diff = np.asanyarray(img.dataobj, dtype="uint8") - np.asanyarray(
+        nt_undo3drefit.dataobj, dtype="uint8"
     )
     assert np.sqrt((diff[10:-10, 10:-10, 10:-10] ** 2).mean()) < 0.1
 
     # Check the target grid by 3dWarp and the affine & size interpolated by NiTransforms
-    cmd = f"3dWarp -verb -deoblique -NN -prefix {tmpdir}/deob.nii.gz {tmpdir}/orig.nii.gz"
+    cmd = (
+        f"3dWarp -verb -deoblique -NN -prefix {tmpdir}/deob.nii.gz {tmpdir}/orig.nii.gz"
+    )
     assert check_call([cmd], shell=True) == 0
 
     deobnii = nb.load("deob.nii.gz")
@@ -541,11 +562,12 @@ def test_afni_oblique(tmpdir, parameters, swapaxes, testdata_path, dir_x, dir_y,
 
     # Check resampling in deobliqued grid
     ntdeobnii = apply(
-        Affine(np.eye(4), reference=deobnii.__class__(
-            np.zeros(deobshape, dtype="uint8"),
-            deobaff,
-            deobnii.header
-        )),
+        Affine(
+            np.eye(4),
+            reference=deobnii.__class__(
+                np.zeros(deobshape, dtype="uint8"), deobaff, deobnii.header
+            ),
+        ),
         img,
         order=0,
     )
@@ -560,9 +582,8 @@ def test_afni_oblique(tmpdir, parameters, swapaxes, testdata_path, dir_x, dir_y,
     )
     mask = np.asanyarray(ntdeobmask.dataobj, dtype=bool)
 
-    diff = (
-        np.asanyarray(deobnii.dataobj, dtype="uint8")
-        - np.asanyarray(ntdeobnii.dataobj, dtype="uint8")
+    diff = np.asanyarray(deobnii.dataobj, dtype="uint8") - np.asanyarray(
+        ntdeobnii.dataobj, dtype="uint8"
     )
     assert np.sqrt((diff[mask] ** 2).mean()) < 0.1
 
@@ -592,7 +613,7 @@ def _generate_reoriented(path, directions, swapaxes, parameters):
         aff = np.diag((*directions, 1)) @ aff
 
         for ax in range(3):
-            if (directions[ax] == -1):
+            if directions[ax] == -1:
                 aff[ax, 3] = last_xyz[ax]
                 data = np.flip(data, ax)
 
@@ -622,16 +643,15 @@ def test_itk_linear_h5(tmpdir, data_path, testdata_path):
     assert len(h5xfm.xforms) == 1
 
     # File loadable with single affine object
-    itk.ITKLinearTransform.from_filename(
-        data_path / "affine-antsComposite.h5"
-    )
+    itk.ITKLinearTransform.from_filename(data_path / "affine-antsComposite.h5")
 
     with open(data_path / "affine-antsComposite.h5", "rb") as f:
         itk.ITKLinearTransform.from_fileobj(f)
 
     # Exercise only_linear
     itk.ITKCompositeH5.from_filename(
-        testdata_path / "ds-005_sub-01_from-T1w_to-MNI152NLin2009cAsym_mode-image_xfm.h5",
+        testdata_path
+        / "ds-005_sub-01_from-T1w_to-MNI152NLin2009cAsym_mode-image_xfm.h5",
         only_linear=True,
     )
 
@@ -674,55 +694,25 @@ def test_itk_linear_h5(tmpdir, data_path, testdata_path):
     with pytest.raises(TransformIOError):
         itk.ITKLinearTransform.from_filename("test.h5")
 
-# Added tests for h5 orientation bug
-
-
-@pytest.mark.xfail(
-    reason="GH-137/GH-171: displacement field dimension order is wrong",
-    strict=False,
-)
-def test_itk_h5_field_order(tmp_path):
-    """Displacement fields stored in row-major order should fail to round-trip."""
-    shape = (3, 4, 5)
-    vals = np.arange(np.prod(shape), dtype=float).reshape(shape)
-    field = np.stack([vals, vals + 100, vals + 200], axis=0)
-
-    params = field.reshape(-1, order="C")
-    fixed = np.array(list(shape) + [0, 0, 0] + [1, 1, 1] + list(np.eye(3).ravel()), dtype=float)
-    fname = tmp_path / "field.h5"
-    with H5File(fname, "w") as f:
-        grp = f.create_group("TransformGroup")
-        grp.create_group("0")["TransformType"] = np.array([b"CompositeTransform_double_3_3"])
-        g1 = grp.create_group("1")
-        g1["TransformType"] = np.array([b"DisplacementFieldTransform_float_3_3"])
-        g1["TransformFixedParameters"] = fixed
-        g1["TransformParameters"] = params
-
-    img = itk.ITKCompositeH5.from_filename(fname)[0]
-    expected = np.moveaxis(field, 0, -1)
-    expected[..., (0, 1)] *= -1
-    assert np.allclose(img.get_fdata(), expected)
-
 
+# Added tests for h5 orientation bug (#167)
 def _load_composite_testdata(data_path):
     """Return the composite HDF5 and displacement field from regressions."""
     h5file = data_path / "regressions" / "ants_t1_to_mniComposite.h5"
     # Generated using
     # CompositeTransformUtil --disassemble ants_t1_to_mniComposite.h5 \
     #     ants_t1_to_mniComposite
-    warpfile = data_path / "regressions" / (
-        "01_ants_t1_to_mniComposite_DisplacementFieldTransform.nii.gz"
+    warpfile = (
+        data_path
+        / "regressions"
+        / ("01_ants_t1_to_mniComposite_DisplacementFieldTransform.nii.gz")
     )
     if not (h5file.exists() and warpfile.exists()):
         pytest.skip("Composite transform test data not available")
     return h5file, warpfile
 
 
-@pytest.mark.xfail(
-    reason="GH-137/GH-171: displacement field dimension order is wrong",
-    strict=False,
-)
-def test_itk_h5_displacement_mismatch(testdata_path):
+def test_itk_h5_and_displacement_equivalence(testdata_path):
     """Composite displacements should match the standalone field"""
     h5file, warpfile = _load_composite_testdata(testdata_path)
     xforms = itk.ITKCompositeH5.from_filename(h5file)
@@ -733,80 +723,21 @@ def test_itk_h5_displacement_mismatch(testdata_path):
         np.asanyarray(field_h5.dataobj), np.asanyarray(field_img.dataobj)
     )
 
+    points = np.array(
+        [
+            [0.0, 0.0, 0.0],
+            [1.0, 2.0, 3.0],
+            [10.0, -10.0, 5.0],
+            [-5.0, 7.0, -2.0],
+            [12.0, 0.0, -11.0],
+        ]
+    )
 
-def test_itk_h5_transpose_fix(testdata_path):
-    """Check the displacement field orientation explicitly.
-
-    ITK stores displacement fields with the vector dimension leading in
-    Fortran (column-major) order [1]_. Transposing the parameters from the HDF5
-    composite file accordingly should match the standalone displacement image.
-
-    References
-    ----------
-    .. [1] ITK Software Guide. https://itk.org/ItkSoftwareGuide.pdf
-    """
-    h5file, warpfile = _load_composite_testdata(testdata_path)
-
-    with H5File(h5file, "r") as f:
-        group = f["TransformGroup"]["2"]
-        size = group["TransformFixedParameters"][:3].astype(int)
-        params = group["TransformParameters"][:].reshape(*size, 3)
-
-    img = nb.load(warpfile)
-    ref = np.squeeze(np.asanyarray(img.dataobj))
-
-    np.testing.assert_array_equal(params.transpose(2, 1, 0, 3), ref)
-
-
-def test_itk_h5_field_order_fortran(tmp_path):
-    """Verify Fortran-order displacement fields load correctly"""
-    shape = (3, 4, 5)
-    vals = np.arange(np.prod(shape), dtype=float).reshape(shape)
-    field = np.stack([vals, vals + 100, vals + 200], axis=0)
-
-    params = field.reshape(-1, order="F")
-    fixed = np.array(list(shape) + [0, 0, 0] + [1, 1, 1] + list(np.eye(3).ravel()), dtype=float)
-    fname = tmp_path / "field_f.h5"
-    with H5File(fname, "w") as f:
-        grp = f.create_group("TransformGroup")
-        grp.create_group("0")["TransformType"] = np.array([b"CompositeTransform_double_3_3"])
-        g1 = grp.create_group("1")
-        g1["TransformType"] = np.array([b"DisplacementFieldTransform_float_3_3"])
-        g1["TransformFixedParameters"] = fixed
-        g1["TransformParameters"] = params
-
-    img = itk.ITKCompositeH5.from_filename(fname)[0]
-    expected = np.moveaxis(field, 0, -1)
-    expected[..., (0, 1)] *= -1
-    assert np.allclose(img.get_fdata(), expected)
-
-
-@pytest.mark.xfail(strict=False, reason="Results may not match ANTs exactly")
-def test_composite_h5_map_against_ants(testdata_path, tmp_path):
-    """Map points with NiTransforms and compare to ANTs."""
-    h5file = testdata_path / "regressions" / "ants_t1_to_mniComposite.h5"
-    if not h5file.exists():
-        pytest.skip("Composite transform test data not available")
-
-    points = np.array([[0.0, 0.0, 0.0], [1.0, 2.0, 3.0]])
-    csvin = tmp_path / "points.csv"
-    np.savetxt(csvin, points, delimiter=",", header="x,y,z", comments="")
-
-    csvout = tmp_path / "out.csv"
-    cmd = f"antsApplyTransformsToPoints -d 3 -i {csvin} -o {csvout} -t {h5file}"
-    exe = cmd.split()[0]
-    if not shutil.which(exe):
-        pytest.skip(f"Command {exe} not found on host")
-    check_call(cmd, shell=True)
-
-    ants_res = np.genfromtxt(csvout, delimiter=",", names=True)
-    ants_pts = np.vstack([ants_res[n] for n in ("x", "y", "z")]).T
+    # Load the displacements field
+    xfm_h5 = nitnl.DenseFieldTransform(xforms[1])
+    mapped_h5 = xfm_h5.map(points)
 
-    xforms = itk.ITKCompositeH5.from_filename(h5file)
-    chain = (
-        nitl.Affine(xforms[0].to_ras())
-        + nitnl.DenseFieldTransform(xforms[1])
-    )
-    mapped = chain.map(points)
+    xfm_warp = nitnl.DenseFieldTransform(field_img)
+    mapped_warp = xfm_warp.map(points)
 
-    assert np.allclose(mapped, ants_pts, atol=1e-6)
+    np.testing.assert_array_equal(mapped_h5, mapped_warp)