fix: address offset with factors

Author: oesteban

sdcflows/utils/tests/test_tools.py

@@ -36,31 +36,45 @@ def test_epi_mask(tmpdir, testdata_dir):
 
 
 @pytest.mark.parametrize("padding", [0, 1, 4])
-@pytest.mark.parametrize("factor", [1, 4, 0.5])
-def test_deoblique_and_zooms(padding, factor):
+@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):
     """Check deoblique and denser."""
+    tmpdir.chdir()
 
     # Generate an example reference image
-    ref_data = np.zeros((20, 30, 40), dtype=np.float32)
-    ref_affine = np.diag([1.0, 1.2, 0.8, 1.0])
+    ref_data = np.zeros((20, 32, 40), 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_img = nb.Nifti1Image(ref_data, ref_affine)
-    ref_zooms = np.array(ref_img.header.get_zooms()[:3])
 
     # Generate an example oblique image
-    ob_data = np.ones_like(ref_data)
-
-    # Rotate 90 degrees around x-axis
-    rotate = np.array([[0, 1, 0, 0], [0, 0, 1, 0], [1, 0, 0, 0], [0, 0, 0, 1]])
-    ob_img = nb.Nifti1Image(ob_data, rotate @ ref_affine)
+    ob_img = nb.Nifti1Image(ref_data, rotate @ ref_affine)
 
     # Call function with default parameters
     out_img = deoblique_and_zooms(ref_img, ob_img, padding=padding, factor=factor)
 
     # Check output shape and zooms
     assert np.allclose(out_img.header.get_zooms()[:3], ref_zooms / factor)
 
+    ref_resampled = Affine(reference=out_img).apply(ref_img, order=0)
     resampled = Affine(reference=out_img).apply(ob_img, order=0)
-    ref_volume = ref_data.size * ref_zooms.prod()
+    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

sdcflows/utils/tools.py

@@ -43,13 +43,12 @@ 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, rescale_affine
+    from nibabel.affines import apply_affine
 
     # Reference space metadata
     hdr = in_reference.header.copy()
     affine = in_reference.affine
     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)
 
@@ -60,9 +59,6 @@ def deoblique_and_zooms(in_reference, oblique, factor=4, padding=1):
     )
     extent_ijk = apply_affine(np.linalg.inv(affine), extent)
 
-    if isinstance(padding, int):
-        padding = tuple([padding] * 3)
-
     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)
     if np.any(underflow < 0) or np.any(overflow > 0):
@@ -73,9 +69,11 @@ def deoblique_and_zooms(in_reference, oblique, factor=4, padding=1):
 
     # Make grid denser
     if abs(1.0 - factor) > 1e-5:
-        new_shape = np.rint(ref_shape * factor)
-        affine = rescale_affine(affine, ref_shape, ref_zooms / factor, new_shape)
-        ref_shape = new_shape
+        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
 
     # Generate new reference
     hdr.set_sform(affine, scode)