diff --git a/niworkflows/interfaces/images.py b/niworkflows/interfaces/images.py
index e625ef7fbde..2edcad9c11d 100644
--- a/niworkflows/interfaces/images.py
+++ b/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(
+        traits.Bool(),
+        traits.Float(),
+        default=False,
+        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
 
diff --git a/niworkflows/utils/images.py b/niworkflows/utils/images.py
index 2639fda874a..18693a24fa5 100644
--- a/niworkflows/utils/images.py
+++ b/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
@@ -163,8 +163,6 @@ def resample_by_spacing(in_file, zooms, order=3, clip=True):
     qform, qcode = in_file.get_qform(coded=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,11 +193,18 @@ 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
+        if smooth is True:
+            smooth = np.maximum(0, (pre_zooms / zooms - 1) / 2)
+        data = gaussian_filter(in_file.get_fdata(), smooth)
+    else:
+        data = np.asarray(in_file.dataobj)
+
     # Resample data in the new grid
     resampled = map_coordinates(
         data,
         ijk[:3, :],
-        output=dtype,
         order=order,
         mode="constant",
         cval=0,