ENH: Add smooth input to RegridToZooms

Resampling (up- or down-) typically requires a gussian smoothing to be
run before, so that outliers are not introduced (esp. when
downsampling).

This PR adds this option, relying on scipy.

Author: oesteban

niworkflows/interfaces/images.py

@@ -42,6 +42,13 @@ class _RegridToZoomsInputSpec(BaseInterfaceInputSpec):
         usedefault=True,
         desc="clip the data array within the original image's range",
     )
+    smooth = traits.Either(
+        False,
+        traits.Bool(),
+        traits.Float(),
+        usedefault=True,
+        desc="apply gaussian smoothing before resampling"
+    )
 
 
 class _RegridToZoomsOutputSpec(TraitedSpec):
@@ -65,6 +72,7 @@ def _run_interface(self, runtime):
             self.inputs.zooms,
             order=self.inputs.order,
             clip=self.inputs.clip,
+            smooth=self.inputs.smooth,
         ).to_filename(self._results["out_file"])
         return runtime
 

niworkflows/utils/images.py

@@ -148,7 +148,7 @@ def dseg_label(in_seg, label, newpath=None):
     return out_file
 
 
-def resample_by_spacing(in_file, zooms, order=3, clip=True):
+def resample_by_spacing(in_file, zooms, order=3, clip=True, smooth=False):
     """Regrid the input image to match the new zooms."""
     from pathlib import Path
     import numpy as np
@@ -164,7 +164,6 @@ def resample_by_spacing(in_file, zooms, order=3, clip=True):
 
     hdr = in_file.header.copy()
     dtype = hdr.get_data_dtype()
-    data = np.asanyarray(in_file.dataobj)
     zooms = np.array(zooms)
 
     # Calculate the factors to normalize voxel size to the specific zooms
@@ -195,6 +194,13 @@ def resample_by_spacing(in_file, zooms, order=3, clip=True):
         new_card.dot(np.vstack((new_grid, np.ones((1, new_grid.shape[1])))))
     )
 
+    if smooth:
+        from scipy.ndimage import gaussian_filter
+        data = gaussian_filter(in_file.get_fdata(),
+                               2 if smooth is True else smooth).astype(dtype)
+    else:
+        data = np.asanyarray(in_file.dataobj)
+
     # Resample data in the new grid
     resampled = map_coordinates(
         data,