From d24ffe8bb3ff3067c82f811676a8256a03be3f76 Mon Sep 17 00:00:00 2001 From: georgewilliamstrong Date: Fri, 27 Oct 2023 14:03:42 +0100 Subject: [PATCH 1/2] add anti-aliasing filter for model interpolation --- stride/problem/data.py | 64 +++++++++++++++++++++++++++++++++--------- 1 file changed, 50 insertions(+), 14 deletions(-) diff --git a/stride/problem/data.py b/stride/problem/data.py index c64576b0..d2613e27 100644 --- a/stride/problem/data.py +++ b/stride/problem/data.py @@ -899,7 +899,7 @@ def stagger_data(self, data, stagger, method='nearest'): return interp - def resample(self, space=None, order=3, prefilter=True, **kwargs): + def resample(self, space=None, **kwargs): """ Resample the internal (non-padded) data given some new space object. @@ -911,7 +911,15 @@ def resample(self, space=None, order=3, prefilter=True, **kwargs): Order of the interplation, default is 3. prefilter : bool, optional Determines if the input array is prefiltered - before interpolation. The default is ``True``. + before interpolation. If downsampling, this defaults to ``False`` as an anti-aliasing filter + will be applied instead. If upsampling, this defaults to ``True``. + anti_alias : bool, optional + Whether a Gaussian filter is applied to smooth the data before interpolation. + The default is ``True``. This is only applied when downsampling. + anti_alias_sigma : float or tuple of floats, optional + Gaussian filter standard deviations used for the anti-aliasing filter. + The default is (d - 1) / 2 where d is the downsampling factor and d > 1. When upsampling, + d < 1, and no anti-aliasing filter is applied. Returns ------- @@ -923,22 +931,18 @@ def resample(self, space=None, order=3, prefilter=True, **kwargs): interp = [] for t in range(data.shape[0]): - interp.append(self.resample_data(data[t], space, - order=order, - prefilter=prefilter)) + interp.append(self.resample_data(data[t], space, **kwargs)) interp = np.stack(interp, axis=0) else: - interp = self.resample_data(self.data, space, - order=order, - prefilter=prefilter) + interp = self.resample_data(self.data, space, **kwargs) self.grid.space = space self._init_shape() self._data = self.pad_data(interp) - def resample_data(self, data, space, order=3, prefilter=True): + def resample_data(self, data, space, **kwargs): """ Resample the data given some new space object. @@ -949,10 +953,18 @@ def resample_data(self, data, space, order=3, prefilter=True): space : Space New space. order : int, optional - Order of the interplation, default is 3. + Order of the interpolation, default is 3. prefilter : bool, optional Determines if the input array is prefiltered - before interpolation. The default is ``True``. + before interpolation. If downsampling, this defaults to ``False`` as an anti-aliasing filter + will be applied instead. If upsampling, this defaults to ``True``. + anti_alias : bool, optional + Whether a Gaussian filter is applied to smooth the data before interpolation. + The default is ``True``. This is only applied when downsampling. + anti_alias_sigma : float or tuple of floats, optional + Gaussian filter standard deviations used for the anti-aliasing filter. + The default is (d - 1) / 2 where d is the downsampling factor and d > 1. When upsampling, + d < 1, and no anti-aliasing filter is applied. Returns ------- @@ -960,11 +972,35 @@ def resample_data(self, data, space, order=3, prefilter=True): Resampled data. """ + order = kwargs.pop('order', 3) + prefilter = kwargs.pop('prefilter', True) + + resampling_factors = np.array([dx_old/dx_new + for dx_old, dx_new in zip(self.space.spacing, space.spacing)]) + + # Anti-aliasing is only required for down-sampling interpolation + if any(factor > 1 for factor in resampling_factors): + anti_alias = kwargs.pop('anti_alias', True) + + if anti_alias: + anti_alias_sigma = kwargs.pop('anti_alias_sigma', None) + + if anti_alias_sigma is not None: + anti_alias_sigma = anti_alias_sigma * np.ones_like(resampling_factors) + + if np.any(anti_alias_sigma < 0): + raise ValueError("Anti-alias standard dev. must be equal to or greater than zero") + + # Estimate anti-alias standard deviations if none provided + else: + anti_alias_sigma = np.maximum(0, (1/resampling_factors - 1) / 2) + + data = scipy.ndimage.gaussian_filter(data, anti_alias_sigma) - resampling_factor = [dx_old/dx_new - for dx_old, dx_new in zip(self.space.spacing, space.spacing)] + # Prefiltering is not necessary if anti-alias filter used + prefilter = False - interp = scipy.ndimage.zoom(data, resampling_factor, + interp = scipy.ndimage.zoom(data, resampling_factors, order=order, prefilter=prefilter) return interp From f8be508817de7d74c30149a9911e0bb2bb80b263 Mon Sep 17 00:00:00 2001 From: georgewilliamstrong Date: Fri, 27 Oct 2023 14:34:25 +0100 Subject: [PATCH 2/2] fix condition --- stride/problem/data.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stride/problem/data.py b/stride/problem/data.py index d2613e27..96fa2973 100644 --- a/stride/problem/data.py +++ b/stride/problem/data.py @@ -979,7 +979,7 @@ def resample_data(self, data, space, **kwargs): for dx_old, dx_new in zip(self.space.spacing, space.spacing)]) # Anti-aliasing is only required for down-sampling interpolation - if any(factor > 1 for factor in resampling_factors): + if any(factor < 1 for factor in resampling_factors): anti_alias = kwargs.pop('anti_alias', True) if anti_alias: