fix: finalize the implementation of a test

after long time scratching my head, it seems mango does not always
interpret orientation correctly. freeview seems to do so though.

Author: oesteban

sdcflows/utils/tests/test_tools.py

@@ -38,43 +38,77 @@ def test_epi_mask(tmpdir, testdata_dir):
 @pytest.mark.parametrize("padding", [0, 1, 4])
 @pytest.mark.parametrize("factor", [1, 4, 0.8])
 @pytest.mark.parametrize("centered", [True, False])
-@pytest.mark.parametrize("rotate", [
-    np.eye(4),
-    # Rotate 90 degrees around x-axis
-    np.array([[0, 1, 0, 0], [0, 0, 1, 0], [1, 0, 0, 0], [0, 0, 0, 1]]),
-])
-def test_deoblique_and_zooms(tmpdir, padding, factor, centered, rotate):
+@pytest.mark.parametrize(
+    "rotate",
+    [
+        np.eye(4),
+        # Rotate 90 degrees around x-axis
+        np.array([[0, 1, 0, 0], [0, 0, 1, 0], [1, 0, 0, 0], [0, 0, 0, 1]]),
+        # Rotate 30 degrees around x-axis
+        nb.affines.from_matvec(
+            nb.eulerangles.euler2mat(0, 0, 30 * np.pi / 180),
+            (0, 0, 0),
+        ),
+        # Rotate 48 degrees around y-axis and translation
+        nb.affines.from_matvec(
+            nb.eulerangles.euler2mat(0, 48 * np.pi / 180, 0),
+            (2.0, 1.2, 0.7),
+        ),
+    ],
+)
+def test_deoblique_and_zooms(tmpdir, padding, factor, centered, rotate, debug=False):
     """Check deoblique and denser."""
     tmpdir.chdir()
 
     # Generate an example reference image
-    ref_data = np.zeros((20, 32, 40), dtype=np.float32)
+    ref_shape = (80, 128, 160) if factor < 1 else (20, 32, 40)
+    ref_data = np.zeros(ref_shape, dtype=np.float32)
     ref_data[1:-2, 10:-11, 1:-2] = 1
     ref_affine = np.diag([2.0, 1.25, 1.0, 1.0])
     ref_zooms = nb.affines.voxel_sizes(ref_affine)
     if centered:
         ref_affine[:3, 3] -= nb.affines.apply_affine(
-            ref_affine, 0.5 * (np.array(ref_data.shape) - 1),
+            ref_affine,
+            0.5 * (np.array(ref_data.shape) - 1),
         )
     ref_img = nb.Nifti1Image(ref_data, ref_affine)
+    ref_img.header.set_qform(ref_affine, 1)
+    ref_img.header.set_sform(ref_affine, 0)
 
     # Generate an example oblique image
-    ob_img = nb.Nifti1Image(ref_data, rotate @ ref_affine)
+    mov_affine = rotate @ ref_affine
+    mov_img = nb.Nifti1Image(ref_data, mov_affine)
+    mov_img.header.set_qform(mov_affine, 1)
+    mov_img.header.set_sform(mov_affine, 0)
 
     # Call function with default parameters
-    out_img = deoblique_and_zooms(ref_img, ob_img, padding=padding, factor=factor)
+    out_img = deoblique_and_zooms(ref_img, mov_img, padding=padding, factor=factor)
 
     # Check output shape and zooms
     assert np.allclose(out_img.header.get_zooms()[:3], ref_zooms / factor)
 
+    # Check idempotency with a lot of tolerance
     ref_resampled = Affine(reference=out_img).apply(ref_img, order=0)
-    resampled = Affine(reference=out_img).apply(ob_img, order=0)
-    ref_img.to_filename("reference.nii.gz")
-    ob_img.to_filename("moving.nii.gz")
-    ref_resampled.to_filename("ref_resampled.nii.gz")
-    resampled.to_filename("resampled.nii.gz")
-    # import pdb; pdb.set_trace()
-    ref_volume = ref_data.sum() * ref_zooms.prod()
-    res_volume = resampled.get_fdata().sum() * np.prod(resampled.header.get_zooms())
-    # 20% of change in volume is too high, must be an error somewhere
-    assert abs(ref_volume - res_volume) < ref_volume * 0.2
+    ref_back = Affine(reference=ref_img).apply(ref_resampled, order=0)
+    resampled = Affine(reference=out_img).apply(mov_img, order=2)
+    if debug:
+        ref_img.to_filename("reference.nii.gz")
+        ref_back.to_filename("reference_back.nii.gz")
+        ref_resampled.to_filename("reference_resampled.nii.gz")
+        mov_img.to_filename("moving.nii.gz")
+        resampled.to_filename("resampled.nii.gz")
+
+    # Allow up to 3% pixels wrong after up(down)sampling and walk back.
+    assert (
+        np.abs(np.clip(ref_back.get_fdata(), 0, 1) - ref_data).sum()
+        < ref_data.size * 0.03
+    )
+
+    vox_vol_out = np.prod(out_img.header.get_zooms())
+    vox_vol_mov = np.prod(mov_img.header.get_zooms())
+    vol_factor = vox_vol_out / vox_vol_mov
+
+    ref_volume = ref_data.sum()
+    res_volume = np.clip(resampled.get_fdata(), 0, 1).sum() * vol_factor
+    # Tolerate up to 2% variation of the volume of the moving image
+    assert abs(1 - res_volume / ref_volume) < 0.02

sdcflows/utils/tools.py

@@ -43,37 +43,38 @@ def deoblique_and_zooms(in_reference, oblique, factor=4, padding=1):
     """
     from itertools import product
     import numpy as np
-    from nibabel.affines import apply_affine
+    from nibabel.affines import apply_affine, rescale_affine
 
     # Reference space metadata
     hdr = in_reference.header.copy()
-    affine = in_reference.affine
+    affine = in_reference.affine.copy()
     ref_shape = np.array(in_reference.shape[:3])
+    ref_zooms = np.array(hdr.get_zooms()[:3])
     _, scode = in_reference.get_sform(coded=True)
     _, qcode = in_reference.get_qform(coded=True)
 
     # Calculate the 8 most extreme coordinates of oblique in in_reference space
     extent = apply_affine(
         oblique.affine,
-        np.array(list(product((0, 1), repeat=3))) * (ref_shape - 1),
+        np.array(list(product((0, 1), repeat=3))) * (np.array(oblique.shape[:3]) - 1),
     )
     extent_ijk = apply_affine(np.linalg.inv(affine), extent)
 
     underflow = np.clip(extent_ijk.min(0) - padding, None, 0).astype(int)
-    overflow = np.ceil(np.clip(extent_ijk.max(0) + padding + 1 - ref_shape, 0, None)).astype(int)
+    overflow = np.ceil(
+        np.clip(extent_ijk.max(0) + padding + 1 - ref_shape, 0, None)
+    ).astype(int)
     if np.any(underflow < 0) or np.any(overflow > 0):
         # Add under/overflow voxels
         ref_shape += np.abs(underflow) + overflow
         # Consistently update origin
         affine[:-1, -1] = apply_affine(affine, underflow)
 
     # Make grid denser
-    if abs(1.0 - factor) > 1e-5:
-        cur_center = apply_affine(affine, 0.5 * (ref_shape - 1))
-        affine[:3, :3] /= factor
-        ref_shape = np.rint(ref_shape * factor)
-        new_center = affine[:3, :3] @ (0.5 * (ref_shape - 1))
-        affine[:3, 3] = cur_center - new_center
+    if abs(1.0 - factor) > 1e-4:
+        new_shape = np.rint(ref_shape * factor)
+        affine = rescale_affine(affine, ref_shape, ref_zooms / factor, new_shape)
+        ref_shape = new_shape
 
     # Generate new reference
     hdr.set_sform(affine, scode)