From 1f70c7bd953ece3741dd89ad936b66f49098366a Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 27 Apr 2023 14:55:55 +0200 Subject: [PATCH 1/6] enh: function that calculates the new reference grid --- sdcflows/utils/tools.py | 65 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 65 insertions(+) diff --git a/sdcflows/utils/tools.py b/sdcflows/utils/tools.py index a40b3c5a32..34444a43d7 100644 --- a/sdcflows/utils/tools.py +++ b/sdcflows/utils/tools.py @@ -23,6 +23,71 @@ """Image processing tools.""" +def deoblique_and_zooms(in_reference, oblique, factor=4, padding=1): + """ + Generate a sampling reference aligned with in_reference fully covering oblique. + + Parameters + ---------- + in_reference : :obj:`~nibabel.spatialimages.SpatialImage` + The sampling reference. + oblique : :obj:`~nibabel.spatialimages.SpatialImage` + The rotated coordinate system whose extent should be fully covered by the + generated reference. + factor : :obj:`int` + A factor to increase the resolution of the generated reference. + padding : :obj:`int` + Number of additional voxels around the most extreme positions of the projection of + oblique on to the reference. + + """ + from itertools import product + import numpy as np + from nibabel.affines import apply_affine, rescale_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) + + # 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), + ) + 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): + # 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: + 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) + hdr.set_qform(affine, qcode) + + return in_reference.__class__( + np.zeros(ref_shape.astype(int), dtype=hdr.get_data_dtype()), + affine, + hdr, + ) + + def resample_to_zooms(in_file, zooms, order=3, prefilter=True): """Resample the input data to a new grid with the requested zooms.""" from pathlib import Path From ebc926db77536a2fee4f48aa03c85f9f17cf6d62 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 27 Apr 2023 15:40:19 +0200 Subject: [PATCH 2/6] wip: add a test --- sdcflows/utils/tests/test_tools.py | 37 +++++++++++++++++++++++++++++- 1 file changed, 36 insertions(+), 1 deletion(-) diff --git a/sdcflows/utils/tests/test_tools.py b/sdcflows/utils/tests/test_tools.py index 3b433faf23..6da6043d15 100644 --- a/sdcflows/utils/tests/test_tools.py +++ b/sdcflows/utils/tests/test_tools.py @@ -22,8 +22,10 @@ # """Test EPI manipulation routines.""" import numpy as np +import pytest import nibabel as nb -from ..tools import brain_masker +from nitransforms.linear import Affine +from sdcflows.utils.tools import brain_masker, deoblique_and_zooms def test_epi_mask(tmpdir, testdata_dir): @@ -31,3 +33,36 @@ def test_epi_mask(tmpdir, testdata_dir): tmpdir.chdir() mask = brain_masker(testdata_dir / "epi.nii.gz")[-1] assert abs(np.asanyarray(nb.load(mask).dataobj).sum() - 166476) < 10 + + +@pytest.mark.parametrize("padding", [0, 1, 4]) +@pytest.mark.parametrize("factor", [1, 4, 0.5]) +def test_deoblique_and_zooms(padding, factor): + """Check deoblique and denser.""" + + # Generate an example reference image + ref_data = np.zeros((20, 30, 40), dtype=np.float32) + ref_affine = np.eye(4) + ref_affine[:3, :3] = np.diag([1.0, 1.2, 0.8]) # Set zooms to (2, 3, 4) + 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 = np.eye(4) + + # Rotate 90 degrees around x-axis + rotate[:3, :3] = np.array([[0, 1, 0], [0, 0, 1], [1, 0, 0]]) + ob_img = nb.Nifti1Image(ob_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) + + resampled = Affine(reference=out_img).apply(ob_img, order=0) + ref_volume = ref_data.size * 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 From fb142d58bd3d33984108de531434a1a4927ba416 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 27 Apr 2023 17:13:34 +0200 Subject: [PATCH 3/6] fix: test minor items Co-authored-by: Chris Markiewicz --- sdcflows/utils/tests/test_tools.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/sdcflows/utils/tests/test_tools.py b/sdcflows/utils/tests/test_tools.py index 6da6043d15..0817eb9374 100644 --- a/sdcflows/utils/tests/test_tools.py +++ b/sdcflows/utils/tests/test_tools.py @@ -42,17 +42,15 @@ def test_deoblique_and_zooms(padding, factor): # Generate an example reference image ref_data = np.zeros((20, 30, 40), dtype=np.float32) - ref_affine = np.eye(4) - ref_affine[:3, :3] = np.diag([1.0, 1.2, 0.8]) # Set zooms to (2, 3, 4) + ref_affine = np.diag([1.0, 1.2, 0.8, 1.0]) 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 = np.eye(4) # Rotate 90 degrees around x-axis - rotate[:3, :3] = np.array([[0, 1, 0], [0, 0, 1], [1, 0, 0]]) + 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) # Call function with default parameters From ac802663a98dd7228e926109b91f0ac359be3cf1 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 27 Apr 2023 19:40:16 +0200 Subject: [PATCH 4/6] fix: address offset with factors --- sdcflows/utils/tests/test_tools.py | 36 +++++++++++++++++++++--------- sdcflows/utils/tools.py | 14 +++++------- 2 files changed, 31 insertions(+), 19 deletions(-) diff --git a/sdcflows/utils/tests/test_tools.py b/sdcflows/utils/tests/test_tools.py index 0817eb9374..8b4abd18b8 100644 --- a/sdcflows/utils/tests/test_tools.py +++ b/sdcflows/utils/tests/test_tools.py @@ -36,22 +36,30 @@ 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) @@ -59,8 +67,14 @@ def test_deoblique_and_zooms(padding, 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 diff --git a/sdcflows/utils/tools.py b/sdcflows/utils/tools.py index 34444a43d7..5a0c67c6b2 100644 --- a/sdcflows/utils/tools.py +++ b/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) From f03d54abba834737f6c768f2b8d0139f1df58ef8 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Fri, 28 Apr 2023 08:46:04 +0200 Subject: [PATCH 5/6] 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) From f01cee13425aba2542fee85011453c039477650d Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Tue, 2 May 2023 09:28:57 +0200 Subject: [PATCH 6/6] Apply suggestions from code review Co-authored-by: Mathias Goncalves Co-authored-by: Chris Markiewicz --- sdcflows/utils/tools.py | 24 ++++++++++++++++-------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/sdcflows/utils/tools.py b/sdcflows/utils/tools.py index 1e4c95fd49..99200504e0 100644 --- a/sdcflows/utils/tools.py +++ b/sdcflows/utils/tools.py @@ -21,9 +21,16 @@ # https://www.nipreps.org/community/licensing/ # """Image processing tools.""" +import nibabel as nb -def deoblique_and_zooms(in_reference, oblique, factor=4, padding=1): +def deoblique_and_zooms( + in_reference: nb.spatialimages.SpatialImage, + oblique: nb.spatialimages.SpatialImage, + factor: int = 4, + padding: int = 1, + factor_tol: float = 1e-4, +): """ Generate a sampling reference aligned with in_reference fully covering oblique. @@ -39,6 +46,8 @@ def deoblique_and_zooms(in_reference, oblique, factor=4, padding=1): padding : :obj:`int` Number of additional voxels around the most extreme positions of the projection of oblique on to the reference. + factor_tol : :obj:`float` + Absolute tolerance to determine whether factor is one. """ from itertools import product @@ -54,11 +63,10 @@ def deoblique_and_zooms(in_reference, oblique, factor=4, padding=1): _, 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))) * (np.array(oblique.shape[:3]) - 1), + corners = np.array(list(product((0, 1), repeat=3))) * ( + np.array(oblique.shape[:3]) - 1 ) - extent_ijk = apply_affine(np.linalg.inv(affine), extent) + extent_ijk = apply_affine(np.linalg.inv(affine) @ oblique.affine, corners) underflow = np.clip(extent_ijk.min(0) - padding, None, 0).astype(int) overflow = np.ceil( @@ -66,12 +74,12 @@ def deoblique_and_zooms(in_reference, oblique, factor=4, padding=1): ).astype(int) if np.any(underflow < 0) or np.any(overflow > 0): # Add under/overflow voxels - ref_shape += np.abs(underflow) + overflow + ref_shape += overflow - underflow # Consistently update origin affine[:-1, -1] = apply_affine(affine, underflow) # Make grid denser - if abs(1.0 - factor) > 1e-4: + if abs(1.0 - factor) > factor_tol: new_shape = np.rint(ref_shape * factor) affine = rescale_affine(affine, ref_shape, ref_zooms / factor, new_shape) ref_shape = new_shape @@ -81,7 +89,7 @@ def deoblique_and_zooms(in_reference, oblique, factor=4, padding=1): hdr.set_qform(affine, qcode) return in_reference.__class__( - np.zeros(ref_shape.astype(int), dtype=hdr.get_data_dtype()), + nb.fileslice.strided_scalar(ref_shape.astype(int)), affine, hdr, )