Merge pull request #986 from EmmaRenauld/fix_verifications_resample

Fix verifications in volume_resample

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():

scilpy/image/volume_operations.py

@@ -452,29 +452,32 @@ 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, 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
     ----------
     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
-        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,40 +491,38 @@ def resample_volume(img, ref=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]
 
-    if ref 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)
+    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 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 / v) * z for o, v, 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:
@@ -540,7 +541,7 @@ def resample_volume(img, ref=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:

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,7 +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)
 
     if args.enforce_dimensions and not args.ref:
         parser.error("Cannot enforce dimensions without a reference image")
@@ -84,9 +84,20 @@ 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, 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).")
+
     # Resampling volume
-    resampled_img = resample_volume(img, ref=args.ref, 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)
 

scripts/tests/test_volume_resample.py

@@ -10,17 +10,24 @@
 # 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):
     ret = script_runner.run('scil_volume_resample.py', '--help')
     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