Add Gaussian smoothing before undersampling
Gaussian smoothing is performed previous to the scaling operation if
zoom factor < 1. This reduces the aliasing introduced by the
undersampling operation.
Choice of sigma according to:
Cardoso et al., "Scale factor point spread function matching:
beyond aliasing in image resampling", MICCAI 2015
@@ -4,7 +4,7 @@ from __future__ import absolute_import, print_function
import warnings
import numpy as np
-import scipy.ndimage
+import scipy.ndimage as ndi
from niftynet.layer.base_layer import RandomisedLayer
@@ -22,7 +22,7 @@ class RandomSpatialScalingLayer(RandomisedLayer):
max_percentage=10.0,
name='random_spatial_scaling'):
super(RandomSpatialScalingLayer, self).__init__(name=name)
- assert min_percentage < max_percentage
+ assert min_percentage <= max_percentage
self._min_percentage = max(min_percentage, -99.9)
self._max_percentage = max_percentage
self._rand_zoom = None
@@ -34,6 +34,18 @@ class RandomSpatialScalingLayer(RandomisedLayer):
size=(spatial_rank,))
self._rand_zoom = (rand_zoom + 100.0) / 100.0
+ def _get_sigma(self, zoom):
+ """
+ Compute optimal standard deviation for Gaussian kernel.
+
+ Cardoso et al., "Scale factor point spread function matching:
+ beyond aliasing in image resampling", MICCAI 2015
+ """
+ k = 1 / zoom
+ variance = (k ** 2 - 1 ** 2) * (2 * np.sqrt(2 * np.log(2))) ** (-2)
+ sigma = np.sqrt(variance)
+ return sigma
+
def _apply_transformation(self, image, interp_order=3):
if interp_order < 0:
return image
@@ -41,19 +53,32 @@ class RandomSpatialScalingLayer(RandomisedLayer):
full_zoom = np.array(self._rand_zoom)
while len(full_zoom) < image.ndim:
full_zoom = np.hstack((full_zoom, [1.0]))
+ smooth = all(full_zoom[:3] < 1) # perform smoothing if undersampling
+ smooth = False
+ if smooth:
+ sigma = self._get_sigma(full_zoom[:3])
if image.ndim == 4:
output = []
for mod in range(image.shape[-1]):
- scaled = scipy.ndimage.zoom(image[..., mod],
- full_zoom[:3],
- order=interp_order)
+ if smooth:
+ original_resolution = ndi.gaussian_filter(
+ image[..., mod], sigma)
+ else:
+ original_resolution = image[..., mod]
+ scaled = ndi.zoom(
+ original_resolution, full_zoom[:3], order=interp_order)
output.append(scaled[..., np.newaxis])
return np.concatenate(output, axis=-1)
- if image.ndim == 3:
- scaled = scipy.ndimage.zoom(image,
- full_zoom[:3],
- order=interp_order)
+ elif image.ndim == 3:
+ if smooth:
+ original_resolution = ndi.gaussian_filter(
+ image, sigma)
+ else:
+ original_resolution = image
+ scaled = ndi.zoom(
+ original_resolution, full_zoom[:3], order=interp_order)
return scaled[..., np.newaxis]
+ else:
raise NotImplementedError('not implemented random scaling')
def layer_op(self, inputs, interp_orders, *args, **kwargs):
