From 15d1aa0dce56a7ebeb7c0805ec27ba8a9d660f4e Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Wed, 12 Aug 2020 01:59:09 -0700 Subject: [PATCH 1/5] 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. --- niworkflows/interfaces/images.py | 8 ++++++++ niworkflows/utils/images.py | 10 ++++++++-- 2 files changed, 16 insertions(+), 2 deletions(-) diff --git a/niworkflows/interfaces/images.py b/niworkflows/interfaces/images.py index e625ef7fbde..4535bc4cadc 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( + 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 diff --git a/niworkflows/utils/images.py b/niworkflows/utils/images.py index 2639fda874a..64480b09fc8 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 @@ -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, From cf6d6a0ff247587fbd218a61d9f4ace609659256 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Wed, 12 Aug 2020 15:26:54 +0200 Subject: [PATCH 2/5] Update niworkflows/interfaces/images.py --- niworkflows/interfaces/images.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/niworkflows/interfaces/images.py b/niworkflows/interfaces/images.py index 4535bc4cadc..2edcad9c11d 100644 --- a/niworkflows/interfaces/images.py +++ b/niworkflows/interfaces/images.py @@ -43,9 +43,9 @@ class _RegridToZoomsInputSpec(BaseInterfaceInputSpec): desc="clip the data array within the original image's range", ) smooth = traits.Either( - False, traits.Bool(), traits.Float(), + default=False, usedefault=True, desc="apply gaussian smoothing before resampling" ) From f803bdd0c889c0bf520297f8d83ca3e0e75c3030 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 13 Aug 2020 06:01:10 -0700 Subject: [PATCH 3/5] fix: address code review comments --- niworkflows/utils/images.py | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/niworkflows/utils/images.py b/niworkflows/utils/images.py index 64480b09fc8..225186fc7ed 100644 --- a/niworkflows/utils/images.py +++ b/niworkflows/utils/images.py @@ -195,11 +195,16 @@ def resample_by_spacing(in_file, zooms, order=3, clip=True, smooth=False): ) if smooth: + from numbers import Integral from scipy.ndimage import gaussian_filter - data = gaussian_filter(in_file.get_fdata(), - 2 if smooth is True else smooth).astype(dtype) + if smooth is True: + smooth = np.maximum(0, (pre_zooms / zooms - 1) / 2) + data = gaussian_filter(in_file.get_fdata(), smooth) + if issubclass(dtype, Integral): + data = np.round(data) + data = data.astype(dtype) else: - data = np.asanyarray(in_file.dataobj) + data = np.asarray(in_file.dataobj, dtype=dtype) # Resample data in the new grid resampled = map_coordinates( From 5a2fdd17de123165ed484905be972f634403b698 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 13 Aug 2020 08:52:43 -0700 Subject: [PATCH 4/5] enh: eliminate casting of smoothed data let ``map_coordinates`` do the trick. --- niworkflows/utils/images.py | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/niworkflows/utils/images.py b/niworkflows/utils/images.py index 225186fc7ed..a26402b45a0 100644 --- a/niworkflows/utils/images.py +++ b/niworkflows/utils/images.py @@ -163,7 +163,6 @@ def resample_by_spacing(in_file, zooms, order=3, clip=True, smooth=False): qform, qcode = in_file.get_qform(coded=True) hdr = in_file.header.copy() - dtype = hdr.get_data_dtype() zooms = np.array(zooms) # Calculate the factors to normalize voxel size to the specific zooms @@ -195,22 +194,18 @@ def resample_by_spacing(in_file, zooms, order=3, clip=True, smooth=False): ) if smooth: - from numbers import Integral 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) - if issubclass(dtype, Integral): - data = np.round(data) - data = data.astype(dtype) else: - data = np.asarray(in_file.dataobj, dtype=dtype) + data = np.asarray(in_file.dataobj) # Resample data in the new grid resampled = map_coordinates( data, ijk[:3, :], - output=dtype, + output=hdr.get_data_dtype(), order=order, mode="constant", cval=0, From f52c7fbfbb391cacb3ada8e800eb3bd97aee510c Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Thu, 13 Aug 2020 18:43:14 +0200 Subject: [PATCH 5/5] Update niworkflows/utils/images.py Co-authored-by: Chris Markiewicz --- niworkflows/utils/images.py | 1 - 1 file changed, 1 deletion(-) diff --git a/niworkflows/utils/images.py b/niworkflows/utils/images.py index a26402b45a0..18693a24fa5 100644 --- a/niworkflows/utils/images.py +++ b/niworkflows/utils/images.py @@ -205,7 +205,6 @@ def resample_by_spacing(in_file, zooms, order=3, clip=True, smooth=False): resampled = map_coordinates( data, ijk[:3, :], - output=hdr.get_data_dtype(), order=order, mode="constant", cval=0,