diff --git a/CHANGES.rst b/CHANGES.rst index 88c2a67975c..93218ad5e70 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -29,6 +29,9 @@ New Features - Convolution with un-normalized and un-normalizable kernels is now possible. [#5782] + - Add a new argument, ``normalization_rtol``, to ``convolve_fft``, allowing + the user to specify the relative error tolerance in the normalization of + the convolution kernel. [#5649, #5177] - ``astropy.coordinates`` diff --git a/astropy/convolution/convolve.py b/astropy/convolution/convolve.py index ca8d6db58b0..1d0152d3077 100644 --- a/astropy/convolution/convolve.py +++ b/astropy/convolution/convolve.py @@ -26,7 +26,7 @@ @support_nddata(data='array') def convolve(array, kernel, boundary='fill', fill_value=0., nan_treatment='interpolate', normalize_kernel=True, mask=None, - preserve_nan=False): + preserve_nan=False, normalization_zero_tol=1e-8): ''' Convolve an array with a kernel. @@ -82,6 +82,10 @@ def convolve(array, kernel, boundary='fill', fill_value=0., `None`, no masking will be performed unless ``array`` is a masked array. If ``mask`` is not `None` *and* ``array`` is a masked array, a pixel is masked of it is masked in either ``mask`` *or* ``array.mask``. + normalization_zero_tol: float, optional + The absolute tolerance on whether the kernel is different than zero. + If the kernel sums to zero to within this precision, it cannot be + normalized. Default is "1e-8". Returns ------- @@ -212,7 +216,7 @@ def convolve(array, kernel, boundary='fill', fill_value=0., # explicitly normalize the kernel here, and then scale the image at the # end if normalization was not requested. kernel_sum = kernel_internal.sum() - kernel_sums_to_zero = np.isclose(kernel_sum, 0, atol=1e-8) + kernel_sums_to_zero = np.isclose(kernel_sum, 0, atol=normalization_zero_tol) if (kernel_sum < 1. / MAX_NORMALIZATION or kernel_sums_to_zero) and normalize_kernel: raise Exception("The kernel can't be normalized, because its sum is " @@ -316,11 +320,12 @@ def convolve(array, kernel, boundary='fill', fill_value=0., @support_nddata(data='array') def convolve_fft(array, kernel, boundary='fill', fill_value=0., nan_treatment='interpolate', normalize_kernel=True, + normalization_zero_tol=1e-8, preserve_nan=False, mask=None, crop=True, return_fft=False, fft_pad=None, psf_pad=None, quiet=False, min_wt=0.0, allow_huge=False, fftn=np.fft.fftn, ifftn=np.fft.ifftn, - complex_dtype=np.complex, ): + complex_dtype=np.complex): """ Convolve an ndarray with an nd-kernel. Returns a convolved image with ``shape = array.shape``. Assumes kernel is centered. @@ -379,6 +384,10 @@ def convolve_fft(array, kernel, boundary='fill', fill_value=0., e.g., ``normalize_kernel=np.sum`` means that kernel will be modified to be: ``kernel = kernel / np.sum(kernel)``. If True, defaults to ``normalize_kernel = np.sum``. + normalization_zero_tol: float, optional + The absolute tolerance on whether the kernel is different than zero. + If the kernel sums to zero to within this precision, it cannot be + normalized. Default is "1e-8". preserve_nan : bool After performing convolution, should pixels that were originally NaN again become NaN? @@ -547,7 +556,6 @@ def convolve_fft(array, kernel, boundary='fill', fill_value=0., array[nanmaskarray] = 0 nanmaskkernel = np.isnan(kernel) | np.isinf(kernel) kernel[nanmaskkernel] = 0 - anynan = (np.any(nanmaskarray) or np.any(nanmaskkernel)) if normalize_kernel is True: if kernel.sum() < 1. / MAX_NORMALIZATION: @@ -564,19 +572,17 @@ def convolve_fft(array, kernel, boundary='fill', fill_value=0., normalized_kernel = kernel / kernel_scale else: kernel_scale = kernel.sum() - if np.abs(kernel_scale - 1) < 1e-8: - kernel_scale = 1 - normalized_kernel = kernel - else: - if np.abs(kernel.sum()) < 1e-8: - if nan_treatment == 'interpolate': - raise ValueError('Cannot interpolate NaNs with an unnormalizable kernel') - else: - # the kernel's sum is near-zero, so it can't be scaled - kernel_scale = 1 - normalized_kernel = kernel + if np.abs(kernel_scale) < normalization_zero_tol: + if nan_treatment == 'interpolate': + raise ValueError('Cannot interpolate NaNs with an unnormalizable kernel') else: - normalized_kernel = kernel / kernel_scale + # the kernel's sum is near-zero, so it can't be scaled + kernel_scale = 1 + normalized_kernel = kernel + else: + # the kernel is normalizable; we'll temporarily normalize it + # now and undo the normalization later. + normalized_kernel = kernel / kernel_scale if boundary is None: warnings.warn("The convolve_fft version of boundary=None is " diff --git a/astropy/convolution/tests/test_convolve_fft.py b/astropy/convolution/tests/test_convolve_fft.py index 7b173cd374f..2e5ab6c67d9 100644 --- a/astropy/convolution/tests/test_convolve_fft.py +++ b/astropy/convolution/tests/test_convolve_fft.py @@ -9,6 +9,8 @@ from numpy.testing import assert_array_almost_equal_nulp, assert_allclose from ..convolve import convolve_fft +from ...tests.helper import catch_warnings +from ...utils.exceptions import AstropyUserWarning VALID_DTYPES = [] @@ -307,6 +309,37 @@ def test_normalize_function(self): result = convolve_fft(array, kernel, normalize_kernel=np.max) assert_allclose(result, [3, 6, 5], atol=1E-10) + @pytest.mark.parametrize(option_names, options) + def test_normalization_is_respected(self, boundary, + nan_treatment, + normalize_kernel): + """ + Check that if normalize_kernel is False then the normalization + tolerance is respected. + """ + array = np.array([1, 2, 3]) + # A simple identity kernel to which a non-zero normalization is added. + base_kernel = np.array([1.0]) + + # Use the same normalization error tolerance in all cases. + normalization_rtol = 1e-4 + + # Add the error below to the kernel. + norm_error = [normalization_rtol / 10, normalization_rtol * 10] + + for err in norm_error: + kernel = base_kernel + err + result = convolve_fft(array, kernel, + normalize_kernel=normalize_kernel, + nan_treatment=nan_treatment, + normalization_zero_tol=normalization_rtol) + if normalize_kernel: + # Kernel has been normalized to 1. + assert_allclose(result, array) + else: + # Kernel should not have been normalized... + assert_allclose(result, array * kernel) + class TestConvolve2D(object):