From 44b096c9a555b3c97ef19e389a4de671499d0108 Mon Sep 17 00:00:00 2001 From: EmmaRenauld Date: Fri, 26 Apr 2024 09:38:36 -0400 Subject: [PATCH 1/4] Add test to showcase error --- scripts/scil_volume_resample.py | 1 + scripts/tests/test_volume_resample.py | 13 ++++++++++--- 2 files changed, 11 insertions(+), 3 deletions(-) diff --git a/scripts/scil_volume_resample.py b/scripts/scil_volume_resample.py index e65e5b853..1c74e005d 100755 --- a/scripts/scil_volume_resample.py +++ b/scripts/scil_volume_resample.py @@ -68,6 +68,7 @@ def main(): assert_inputs_exist(parser, args.in_image, args.ref) assert_outputs_exist(parser, args, args.out_image) assert_headers_compatible(parser, args.in_image, args.ref) + assert_headers_compatible(parser, args.in_image, args.ref) if args.enforce_dimensions and not args.ref: parser.error("Cannot enforce dimensions without a reference image") diff --git a/scripts/tests/test_volume_resample.py b/scripts/tests/test_volume_resample.py index 947bbd59c..884f2cb7d 100644 --- a/scripts/tests/test_volume_resample.py +++ b/scripts/tests/test_volume_resample.py @@ -10,6 +10,7 @@ # If they already exist, this only takes 5 seconds (check md5sum) fetch_data(get_testing_files_dict(), keys=['others.zip']) tmp_dir = tempfile.TemporaryDirectory() +in_img = os.path.join(SCILPY_HOME, 'others', 'fa.nii.gz') def test_help_option(script_runner): @@ -17,10 +18,16 @@ def test_help_option(script_runner): assert ret.success -def test_execution_others(script_runner, monkeypatch): +def test_execution_given_size(script_runner, monkeypatch): monkeypatch.chdir(os.path.expanduser(tmp_dir.name)) - in_img = os.path.join(SCILPY_HOME, 'others', - 'fa.nii.gz') ret = script_runner.run('scil_volume_resample.py', in_img, 'fa_resample.nii.gz', '--voxel_size', '2') assert ret.success + + +def test_execution_ref(script_runner, monkeypatch): + monkeypatch.chdir(os.path.expanduser(tmp_dir.name)) + ref = os.path.join(SCILPY_HOME, 'others', 'fa_resample.nii.gz') + ret = script_runner.run('scil_volume_resample.py', in_img, + 'fa_resample2.nii.gz', '--ref', ref) + assert ret.success From 9ee46ff11c71bb98e5a65734dea68a5f07b22b22 Mon Sep 17 00:00:00 2001 From: EmmaRenauld Date: Fri, 26 Apr 2024 09:45:34 -0400 Subject: [PATCH 2/4] Fix verification of header --- scilpy/image/volume_operations.py | 7 +++---- scripts/scil_volume_resample.py | 14 +++++++++++--- 2 files changed, 14 insertions(+), 7 deletions(-) diff --git a/scilpy/image/volume_operations.py b/scilpy/image/volume_operations.py index 5f6586f37..2b69af494 100644 --- a/scilpy/image/volume_operations.py +++ b/scilpy/image/volume_operations.py @@ -452,7 +452,7 @@ def _interp_code_to_order(interp_code): return orders[interp_code] -def resample_volume(img, ref=None, res=None, iso_min=False, zoom=None, +def resample_volume(img, ref_img=None, res=None, iso_min=False, zoom=None, interp='lin', enforce_dimensions=False): """ Function to resample a dataset to match the resolution of another @@ -463,7 +463,7 @@ def resample_volume(img, ref=None, res=None, iso_min=False, zoom=None, ---------- img: nib.Nifti1Image Image to resample. - ref: nib.Nifti1Image + ref_img: nib.Nifti1Image, optional Reference volume to resample to. This method is used only if ref is not None. (default: None) res: tuple, shape (3,) or int, optional @@ -492,11 +492,10 @@ def resample_volume(img, ref=None, res=None, iso_min=False, zoom=None, affine = img.affine original_zooms = img.header.get_zooms()[:3] - if ref is not None: + if ref_img is not None: if iso_min or zoom or res: raise ValueError('Please only provide one option amongst ref, res ' ', zoom or iso_min.') - ref_img = nib.load(ref) new_zooms = ref_img.header.get_zooms()[:3] elif res is not None: if iso_min or zoom: diff --git a/scripts/scil_volume_resample.py b/scripts/scil_volume_resample.py index 1c74e005d..8524831fb 100755 --- a/scripts/scil_volume_resample.py +++ b/scripts/scil_volume_resample.py @@ -12,6 +12,7 @@ import logging import nibabel as nib +import numpy as np from scilpy.io.utils import (add_verbose_arg, add_overwrite_arg, assert_inputs_exist, assert_outputs_exist, @@ -67,8 +68,6 @@ def main(): # Checking args assert_inputs_exist(parser, args.in_image, args.ref) assert_outputs_exist(parser, args, args.out_image) - assert_headers_compatible(parser, args.in_image, args.ref) - assert_headers_compatible(parser, args.in_image, args.ref) if args.enforce_dimensions and not args.ref: parser.error("Cannot enforce dimensions without a reference image") @@ -85,8 +84,17 @@ def main(): img = nib.load(args.in_image) + ref_img = None + if args.ref: + ref_img = nib.load(args.ref) + # Must not verify that headers are compatible. But can verify that, at + # least, that the last column of their affine are compatible. + if not np.array_equal(img.affine[:, -1], ref_img.affine[:, -1]): + parser.error("The --ref image should have the same affine as the " + "input image (but with a different sampling).") + # Resampling volume - resampled_img = resample_volume(img, ref=args.ref, res=args.volume_size, + resampled_img = resample_volume(img, ref_img=ref_img, res=args.volume_size, iso_min=args.iso_min, zoom=args.voxel_size, interp=args.interp, enforce_dimensions=args.enforce_dimensions) From 34520872a8b007848d170016b1c3684fb9129407 Mon Sep 17 00:00:00 2001 From: EmmaRenauld Date: Fri, 26 Apr 2024 10:23:02 -0400 Subject: [PATCH 3/4] Add unit tests. Fix nomenclature: res meant image shape --- scilpy/image/tests/test_volume_operations.py | 25 ++++++- scilpy/image/volume_operations.py | 72 ++++++++++---------- scripts/scil_volume_resample.py | 6 +- 3 files changed, 64 insertions(+), 39 deletions(-) diff --git a/scilpy/image/tests/test_volume_operations.py b/scilpy/image/tests/test_volume_operations.py index b33361390..0808f8f2c 100644 --- a/scilpy/image/tests/test_volume_operations.py +++ b/scilpy/image/tests/test_volume_operations.py @@ -130,15 +130,36 @@ def smooth_to_fwhm(): def test_resample_volume(): + # Input image: 6x6x6 (4x4x4 with zeros around) + # affine as np.eye => voxel size 1x1x1 moving3d = np.pad(np.ones((4, 4, 4)), pad_width=1, mode='constant', constant_values=0) + moving3d_img = nib.Nifti1Image(moving3d, np.eye(4)) + + # Ref: 2x2x2, voxel size 3x3x3 ref3d = np.ones((2, 2, 2)) + ref_affine = np.eye(4)*3 + ref_affine[-1, -1] = 1 - moving3d_img = nib.Nifti1Image(moving3d, np.eye(4)) + # 1) Option volume_shape: expecting an output of 2x2x2, which means + # voxel resolution 3x3x3 + resampled_img = resample_volume(moving3d_img, volume_shape=(2, 2, 2), + interp='nn') + assert_equal(resampled_img.get_fdata(), ref3d) + assert resampled_img.affine[0, 0] == 3 - resampled_img = resample_volume(moving3d_img, res=(2, 2, 2), interp='nn') + # 2) Option reference image that is 2x2x2, resolution 3x3x3. + ref_img = nib.Nifti1Image(ref3d, ref_affine) + resampled_img = resample_volume(moving3d_img, ref_img=ref_img, + interp='nn') + assert_equal(resampled_img.get_fdata(), ref3d) + assert resampled_img.affine[0, 0] == 3 + # 3) Option final resolution 3x3x3, should be of shape 2x2x2 + resampled_img = resample_volume(moving3d_img, voxel_res=(3, 3, 3), + interp='nn') assert_equal(resampled_img.get_fdata(), ref3d) + assert resampled_img.affine[0, 0] == 3 def test_normalize_metric_basic(): diff --git a/scilpy/image/volume_operations.py b/scilpy/image/volume_operations.py index 2b69af494..87d5a5668 100644 --- a/scilpy/image/volume_operations.py +++ b/scilpy/image/volume_operations.py @@ -452,12 +452,15 @@ def _interp_code_to_order(interp_code): return orders[interp_code] -def resample_volume(img, ref_img=None, res=None, iso_min=False, zoom=None, +def resample_volume(img, ref_img=None, volume_shape=None, iso_min=False, + voxel_res=None, interp='lin', enforce_dimensions=False): """ - Function to resample a dataset to match the resolution of another - reference dataset or to the resolution specified as in argument. - One of the following options must be chosen: ref, res or iso_min. + Function to resample a dataset to match the resolution of another reference + dataset or to the resolution specified as in argument. + + One (and only one) of the following options must be chosen: + ref, volume_shape, iso_min or voxel_res. Parameters ---------- @@ -466,15 +469,15 @@ def resample_volume(img, ref_img=None, res=None, iso_min=False, zoom=None, ref_img: nib.Nifti1Image, optional Reference volume to resample to. This method is used only if ref is not None. (default: None) - res: tuple, shape (3,) or int, optional - Resolution to resample to. If the value it is set to is Y, it will - resample to an isotropic resolution of Y x Y x Y. This method is used - only if res is not None. (default: None) + volume_shape: tuple, shape (3,) or int, optional + Final shape to resample to. If the value it is set to is Y, it will + resample to an isotropic shape of Y x Y x Y. This method is used + only if volume_shape is not None. (default: None) iso_min: bool, optional - If true, resample the volume to R x R x R with R being the smallest - current voxel dimension. If false, this method is not used. - zoom: tuple, shape (3,) or float, optional - Set the zoom property of the image at the value specified. + If true, resample the volume to R x R x R resolution, with R being the + smallest current voxel dimension. If false, this method is not used. + voxel_res: tuple, shape (3,) or float, optional + Set the zoom property of the image at the specified resolution. interp: str, optional Interpolation mode. 'nn' = nearest neighbour, 'lin' = linear, 'quad' = quadratic, 'cubic' = cubic. (Default: linear) @@ -488,39 +491,38 @@ def resample_volume(img, ref_img=None, res=None, iso_min=False, zoom=None, Resampled image. """ data = np.asanyarray(img.dataobj) - original_res = data.shape + original_shape = data.shape affine = img.affine original_zooms = img.header.get_zooms()[:3] + error_msg = ('Please only provide one option amongst ref_img, ' + 'volume_shape, voxel_res or iso_min.') if ref_img is not None: - if iso_min or zoom or res: - raise ValueError('Please only provide one option amongst ref, res ' - ', zoom or iso_min.') + if iso_min or voxel_res or volume_shape: + raise ValueError(error_msg) new_zooms = ref_img.header.get_zooms()[:3] - elif res is not None: - if iso_min or zoom: - raise ValueError('Please only provide one option amongst ref, res ' - ', zoom or iso_min.') - if len(res) == 1: - res = res * 3 - new_zooms = tuple((o / r) * z for o, r, - z in zip(original_res, res, original_zooms)) + elif volume_shape is not None: + if iso_min or voxel_res: + raise ValueError(error_msg) + if len(volume_shape) == 1: + volume_shape = volume_shape * 3 + + new_zooms = tuple((o / r) * z for o, r, z in + zip(original_shape, volume_shape, original_zooms)) elif iso_min: - if zoom: - raise ValueError('Please only provide one option amongst ref, res ' - ', zoom or iso_min.') + if voxel_res: + raise ValueError(error_msg) min_zoom = min(original_zooms) new_zooms = (min_zoom, min_zoom, min_zoom) - elif zoom: - new_zooms = zoom - if len(zoom) == 1: - new_zooms = zoom * 3 + elif voxel_res: + new_zooms = voxel_res + if len(voxel_res) == 1: + new_zooms = voxel_res * 3 else: raise ValueError("You must choose the resampling method. Either with" - "a reference volume, or a chosen isometric resolution" - ", or an isometric resampling to the smallest current" - " voxel dimension!") + "a reference volume, or a chosen image shape, " + "or chosen resolution, or option iso_min.") interp_choices = ['nn', 'lin', 'quad', 'cubic'] if interp not in interp_choices: @@ -539,7 +541,7 @@ def resample_volume(img, ref_img=None, res=None, iso_min=False, zoom=None, logging.info('Resampled data affine setup: %s', nib.aff2axcodes(affine2)) if enforce_dimensions: - if ref is None: + if ref_img is None: raise ValueError('enforce_dimensions can only be used with the ref' 'method.') else: diff --git a/scripts/scil_volume_resample.py b/scripts/scil_volume_resample.py index 8524831fb..628e879cc 100755 --- a/scripts/scil_volume_resample.py +++ b/scripts/scil_volume_resample.py @@ -94,8 +94,10 @@ def main(): "input image (but with a different sampling).") # Resampling volume - resampled_img = resample_volume(img, ref_img=ref_img, res=args.volume_size, - iso_min=args.iso_min, zoom=args.voxel_size, + resampled_img = resample_volume(img, ref_img=ref_img, + volume_shape=args.volume_size, + iso_min=args.iso_min, + voxel_res=args.voxel_size, interp=args.interp, enforce_dimensions=args.enforce_dimensions) From 78fbd20e15e3a7b3f3245a05df679630f3b6b2b2 Mon Sep 17 00:00:00 2001 From: EmmaRenauld Date: Fri, 26 Apr 2024 15:05:01 -0400 Subject: [PATCH 4/4] Answer PK's comments --- scilpy/image/volume_operations.py | 2 +- scripts/scil_volume_resample.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/scilpy/image/volume_operations.py b/scilpy/image/volume_operations.py index 87d5a5668..db6980a0e 100644 --- a/scilpy/image/volume_operations.py +++ b/scilpy/image/volume_operations.py @@ -507,7 +507,7 @@ def resample_volume(img, ref_img=None, volume_shape=None, iso_min=False, if len(volume_shape) == 1: volume_shape = volume_shape * 3 - new_zooms = tuple((o / r) * z for o, r, z in + new_zooms = tuple((o / v) * z for o, v, z in zip(original_shape, volume_shape, original_zooms)) elif iso_min: diff --git a/scripts/scil_volume_resample.py b/scripts/scil_volume_resample.py index 628e879cc..d8ac113fa 100755 --- a/scripts/scil_volume_resample.py +++ b/scripts/scil_volume_resample.py @@ -88,7 +88,7 @@ def main(): if args.ref: ref_img = nib.load(args.ref) # Must not verify that headers are compatible. But can verify that, at - # least, that the last column of their affine are compatible. + # least, the last columns of their affines are compatible. if not np.array_equal(img.affine[:, -1], ref_img.affine[:, -1]): parser.error("The --ref image should have the same affine as the " "input image (but with a different sampling).")