Add antialias filtering by recursive Gaussian to resize

SimpleITK/utilities/resize.py

@@ -2,6 +2,13 @@
 from typing import Sequence
 from typing import Union
 from math import ceil
+import collections
+
+
+def _issequence(obj):
+    if isinstance(obj, basestring):
+        return False
+    return isinstance(obj, collections.Sequence)
 
 
 def resize(
@@ -10,8 +17,9 @@ def resize(
     isotropic: bool = True,
     fill: bool = True,
     interpolator=sitk.sitkLinear,
-    outside_pixel_value: Union[int, float] = 0,
-    use_nearest_extrapolator=False,
+    outside_pixel_value: float = 0.0,
+    use_nearest_extrapolator: bool = False,
+    anti_aliasing_sigma: Union[None, float, Sequence[float]] = None,
 ) -> sitk.Image:
     """
     Resize an image to an arbitrary size while retaining the original image's spacial location.
@@ -30,8 +38,12 @@ def resize(
     new_size.
     :param outside_pixel_value: Value used for padding.
     :param interpolator: Interpolator used for resampling.
-    :param use_nearest_extrapolator: If True, use a nearest neighbor extrapolator when resampling, instead of
+    :param use_nearest_extrapolator: If True, use a nearest neighbor for extrapolation when resampling, overridding the
     constant fill value.
+    :param anti_aliasing_sigma: If zero no antialiasing is performed. If a scalar, it is used as the sigma value in
+     physical units for all axes. If None or a sequence, the sigma value for each axis is calculated as
+     $sigma = (new_spacing - old_spacing) / 2$ in physical units. Gaussian smoothing is performed prior to resampling
+     for antialiasing.
     :return: A SimpleITK image with desired size.
     """
     new_spacing = [
@@ -59,6 +71,28 @@ def resize(
 
     new_origin = image.TransformContinuousIndexToPhysicalPoint(new_origin_cidx)
 
+    input_pixel_type = image.GetPixelID()
+
+    if anti_aliasing_sigma is None:
+        # (s-1)/2.0 is the standard deviation of the Gaussian kernel in index space, where s downsample factor defined
+        # by nspc/ospc.
+        anti_aliasing_sigma = [
+            max((nspc - ospc) / 2.0, 0.0)
+            for ospc, nspc in zip(image.GetSpacing(), new_spacing)
+        ]
+    elif not _issequence(anti_aliasing_sigma):
+        anti_aliasing_sigma = [anti_aliasing_sigma] * image.GetDimension()
+
+    if any([s < 0.0 for s in anti_aliasing_sigma]):
+        raise ValueError("anti_aliasing_sigma must be positive, or None.")
+
+    if all([s > 0.0 for s in anti_aliasing_sigma]):
+        image = sitk.SmoothingRecursiveGaussian(image, anti_aliasing_sigma)
+    else:
+        for d, s in enumerate(anti_aliasing_sigma):
+            if s > 0.0:
+                image = sitk.RecursiveGaussian(image, sigma=s, direction=d)
+
     return sitk.Resample(
         image,
         size=new_size,
@@ -68,4 +102,5 @@ def resize(
         defaultPixelValue=outside_pixel_value,
         interpolator=interpolator,
         useNearestNeighborExtrapolator=use_nearest_extrapolator,
+        outputPixelType=input_pixel_type,
     )