enh: replace ANTs' ResampleImageBySpacing with RegridToZooms

Author: oesteban

niworkflows/anat/ants.py

@@ -18,19 +18,19 @@
     N4BiasFieldCorrection,
     ThresholdImage,
 )
-from nipype.interfaces.ants.utils import AI, ResampleImageBySpacing
+from nipype.interfaces.ants.utils import AI
 
 from ..utils.misc import get_template_specs
 from ..utils.connections import pop_file as _pop
 
 # niworkflows
-from ..interfaces.nibabel import ApplyMask
 from ..interfaces.ants import ImageMath
 from ..interfaces.fixes import (
     FixHeaderRegistration as Registration,
     FixHeaderApplyTransforms as ApplyTransforms,
 )
-from ..interfaces.nibabel import Binarize
+from ..interfaces.images import RegridToZooms
+from ..interfaces.nibabel import ApplyMask, Binarize
 
 
 ATROPOS_MODELS = {
@@ -221,15 +221,9 @@ def init_brain_extraction_wf(
         iterfield=["input_image"],
     )
 
-    res_tmpl = pe.Node(
-        ResampleImageBySpacing(out_spacing=(4, 4, 4), apply_smoothing=True),
-        name="res_tmpl",
-    )
-    res_tmpl.inputs.input_image = tpl_target_path
-    res_target = pe.Node(
-        ResampleImageBySpacing(out_spacing=(4, 4, 4), apply_smoothing=True),
-        name="res_target",
-    )
+    res_tmpl = pe.Node(RegridToZooms(in_file=tpl_target_path, zooms=(4, 4, 4), smooth=True),
+                       name="res_tmpl")
+    res_target = pe.Node(RegridToZooms(zooms=(4, 4, 4), smooth=True), name="res_target")
 
     lap_tmpl = pe.Node(ImageMath(operation="Laplacian", op2="1.5 1"), name="lap_tmpl")
     lap_tmpl.inputs.op1 = tpl_target_path
@@ -337,9 +331,9 @@ def init_brain_extraction_wf(
         (inputnode, norm, [("in_mask", fixed_mask_trait)]),
         (inputnode, map_brainmask, [(("in_files", _pop), "reference_image")]),
         (trunc, inu_n4, [("output_image", "input_image")]),
-        (inu_n4, res_target, [(("output_image", _pop), "input_image")]),
-        (res_tmpl, init_aff, [("output_image", "fixed_image")]),
-        (res_target, init_aff, [("output_image", "moving_image")]),
+        (inu_n4, res_target, [(("output_image", _pop), "in_file")]),
+        (res_tmpl, init_aff, [("out_file", "fixed_image")]),
+        (res_target, init_aff, [("out_file", "moving_image")]),
         (init_aff, norm, [("output_transform", "initial_moving_transform")]),
         (norm, map_brainmask, [
             ("reverse_transforms", "transforms"),

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

@@ -127,7 +127,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
@@ -143,7 +143,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
@@ -174,6 +173,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,