From f03d54abba834737f6c768f2b8d0139f1df58ef8 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Fri, 28 Apr 2023 08:46:04 +0200 Subject: [PATCH] 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. --- sdcflows/utils/tests/test_tools.py | 74 ++++++++++++++++++++++-------- sdcflows/utils/tools.py | 21 +++++---- 2 files changed, 65 insertions(+), 30 deletions(-) diff --git a/sdcflows/utils/tests/test_tools.py b/sdcflows/utils/tests/test_tools.py index 8b4abd18b8..438f702e2f 100644 --- a/sdcflows/utils/tests/test_tools.py +++ b/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 diff --git a/sdcflows/utils/tools.py b/sdcflows/utils/tools.py index 5a0c67c6b2..1e4c95fd49 100644 --- a/sdcflows/utils/tools.py +++ b/sdcflows/utils/tools.py @@ -43,24 +43,27 @@ 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 @@ -68,12 +71,10 @@ def deoblique_and_zooms(in_reference, oblique, factor=4, padding=1): 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)