From 960da0bed55096824b532a2dbd0ccf67dc57e0fb Mon Sep 17 00:00:00 2001 From: Sam Van Kooten Date: Wed, 12 Jan 2022 16:51:18 -0700 Subject: [PATCH] Add Gaussian filter kernel for adaptive resampling --- docs/celestial.rst | 54 ++++++++++--- reproject/adaptive/core.py | 15 +++- reproject/adaptive/deforest.pyx | 125 ++++++++++++++++++----------- reproject/adaptive/high_level.py | 33 +++++--- reproject/tests/test_high_level.py | 15 +++- 5 files changed, 166 insertions(+), 76 deletions(-) diff --git a/docs/celestial.rst b/docs/celestial.rst index e4a927ad5..cbe3ad7de 100644 --- a/docs/celestial.rst +++ b/docs/celestial.rst @@ -181,18 +181,46 @@ Adaptive resampling =================== The :func:`~reproject.reproject_adaptive` function can be used to carry out -reprojection using the `DeForest (2004) +anti-aliased reprojection using the `DeForest (2004) `_ algorithm:: >>> from reproject import reproject_adaptive -In addition to the arguments described in :ref:`common`, one can enable a -rescaling of output pixel values to conserve flux by using the -``conserve_flux`` flag. - -Additionally, one can control the calculation of the Jacobian used in this -algorithm with the ``center_jacobian`` flag. The Jacobian matrix represents how -the corresponding input-image coordinate varies as you move between output +Options +------- + +In addition to the arguments described in :ref:`common`, one can use the +options described below. + +A rescaling of output pixel values to conserve flux can be enabled with the +``conserve_flux`` flag. (Flux conservation is enhanced with a Gaussian +kernel---see below.) + +The kernel used for interpolation and averaging can be controlled with a set of +options. The ``kernel`` argument can be set to 'hann' or 'gaussian' to set the +function being used. The Hann window is the default, and the Gaussian window +improves anti-aliasing and photometric accuracy (or flux conservation, when the +flux-conserving mode is enabled) at the cost of blurring the output image by a +few pixels. The ``kernel_width`` argument sets the width of the Gaussian +kernel, in pixels, and is ignored for the Hann window. This width is measured +between the Gaussian's :math:`\pm 1 \sigma` points. The default value is 1.3 +for the Gaussian, chosen to minimize blurring without compromising accuracy. +Lower values may introduce photometric errors or leave input pixels +under-sampled, while larger values may improve anti-aliasing behavior but will +increase blurring of the output image. Since the Gaussian function has infinite +extent, it must be truncated. This is done by sampling within a region of +finite size. The width in pixels of the sampling region is determined by the +coordinate transform and scaled by the ``sample_region_width`` option, and this +scaling represents a trade-off between accuracy and computation speed. The +default value of 4 represents a reasonable choice, with errors in extreme cases +typically limited to less than one percent, while a value of 5 typically reduces +extreme errors to a fraction of a percent. (The ``sample_region_width`` option +has no effect for the Hann window, as that window does not have infinite +extent.) + +One can control the calculation of the Jacobian used in this +algorithm with the ``center_jacobian`` flag. The Jacobian matrix represents +how the corresponding input-image coordinate varies as you move between output pixels (or d(input image coordinate) / d(output image coordinate)), and serves as a local linearization of the coordinate transformation. When this flag is ``True``, the Jacobian is calculated at pixel grid points by calculating the @@ -208,16 +236,18 @@ points. This is more efficient, and the loss of accuracy is extremely small for transformations that vary smoothly between pixels. The default (``False``) is to use the faster option. +Algorithm Description +--------------------- + Broadly speaking, the algorithm works by approximating the footprint of each output pixel by an elliptical shape in the input image which is stretched and rotated by the transformation (as described by the Jacobian mentioned above), then finding the weighted average of samples inside -that ellipse, where the weight is 1 at the center of the ellipse, and 0 at the -side, and the shape of the weight function is given by an analytical -distribution (currently we use a Hann function). +that ellipse, where the shape of the weighting function is given by an analytical +distribution (Hann and Gaussian functions are supported in this implementation). To illustrate the benefits of this method, we consider a simple case -where the reprojection includes a large change in resoluton. We choose +where the reprojection includes a large change in resolution. We choose to use an artificial data example to better illustrate the differences: .. plot:: diff --git a/reproject/adaptive/core.py b/reproject/adaptive/core.py index 92b52d97c..93a27c6ad 100644 --- a/reproject/adaptive/core.py +++ b/reproject/adaptive/core.py @@ -33,7 +33,9 @@ def __call__(self, pixel_out): def _reproject_adaptive_2d(array, wcs_in, wcs_out, shape_out, return_footprint=True, center_jacobian=False, - roundtrip_coords=True, conserve_flux=False): + roundtrip_coords=True, conserve_flux=False, + kernel='Hann', kernel_width=1.3, + sample_region_width=4): """ Reproject celestial slices from an n-d array from one WCS to another using the DeForest (2004) algorithm [1]_, and assuming all other dimensions @@ -58,6 +60,13 @@ def _reproject_adaptive_2d(array, wcs_in, wcs_out, shape_out, directions. conserve_flux : bool Whether to rescale output pixel values so flux is conserved + kernel : str + The averaging kernel to use. + kernel_width : double + The width of the kernel in pixels. Applies only to the Gaussian kernel. + sample_region_width : double + The width in pixels of the sample region, used only for the Gaussian + kernel which otherwise has infinite extent. Returns ------- @@ -96,7 +105,9 @@ def _reproject_adaptive_2d(array, wcs_in, wcs_out, shape_out, transformer = CoordinateTransformer(wcs_in, wcs_out, roundtrip_coords) map_coordinates(array_in, array_out, transformer, out_of_range_nan=True, - center_jacobian=center_jacobian, conserve_flux=conserve_flux) + center_jacobian=center_jacobian, conserve_flux=conserve_flux, + kernel=kernel, kernel_width=kernel_width, + sample_region_width=sample_region_width) if return_footprint: return array_out, (~np.isnan(array_out)).astype(float) diff --git a/reproject/adaptive/deforest.pyx b/reproject/adaptive/deforest.pyx index f35504d99..41ab048b0 100644 --- a/reproject/adaptive/deforest.pyx +++ b/reproject/adaptive/deforest.pyx @@ -119,8 +119,8 @@ cdef double hanning_filter(double x, double y) nogil: @cython.wraparound(False) @cython.nonecheck(False) @cython.cdivision(True) -cdef double gaussian_filter(double x, double y) nogil: - return exp(-(x*x+y*y) * 1.386294) +cdef double gaussian_filter(double x, double y, double width) nogil: + return exp(-(x*x+y*y) / (width*width) * 2) @cython.boundscheck(False) @@ -163,6 +163,12 @@ cdef double sample_array(double[:,:] source, double x, double y, int x_cyclic, return source[ y, x] +KERNELS = {} +KERNELS['hann'] = 0 +KERNELS['hanning'] = KERNELS['hann'] +KERNELS['gaussian'] = 1 + + @cython.boundscheck(False) @cython.wraparound(False) @cython.nonecheck(False) @@ -170,7 +176,14 @@ cdef double sample_array(double[:,:] source, double x, double y, int x_cyclic, def map_coordinates(double[:,:] source, double[:,:] target, Ci, int max_samples_width=-1, int conserve_flux=False, int progress=False, int singularities_nan=False, int x_cyclic=False, int y_cyclic=False, int out_of_range_nan=False, - bint center_jacobian=False): + bint center_jacobian=False, str kernel='Hann', double kernel_width=1.3, + double sample_region_width=4): + cdef int kernel_flag + try: + kernel_flag = KERNELS[kernel.lower()] + except KeyError: + raise ValueError("'kernel' must be 'Hann' or 'Gaussian'") + cdef np.ndarray[np.float64_t, ndim=3] pixel_target cdef int delta if center_jacobian: @@ -303,48 +316,60 @@ def map_coordinates(double[:,:] source, double[:,:] target, Ci, int max_samples_ s[1] = 1.0/s[1] svd2x2_compose(V, s, U, J) - # We'll need to sample some number of input images to set this - # output pixel. Later on, we'll compute weights to assign to - # each input pixel with a Hanning window, and that window will - # assign weights of zero outside some range. Right now, we'll - # determine a search region within the input image---a bounding - # box around those pixels that will be assigned non-zero - # weights. + # We'll need to sample some number of input image pixels to set + # this output pixel. Later on, we'll compute weights to assign + # to each input pixel, and they will be at or near zero outside + # some range. Right now, we'll determine a search region within + # the input image---a bounding box around those pixels that + # will be assigned non-zero weights. # - # We do that by identifying the locations in the input image of - # the corners of a square region centered around the output - # pixel (using the local linearization of the transformation). - # Those transformed coordinates will set our bounding box. - # - # The output-plane region we're transforming is twice the width - # of a pixel---it runs to the centers of the neighboring - # pixels, rather than the edges of those pixels. When we use - # the Hann window as our filter function, having that window - # stretch to the neighboring pixel centers ensures that, at - # every point, the sum of the overlapping Hann windows is 1, - # and therefore that every input-image pixel is fully - # distributed into some combination of output pixels (in the - # limit of a Jacobian that is constant across all output - # pixels). - - # Transform the corners of the output-plane region to the input - # plane. - P1[0] = - 1 * Ji_padded[0, 0] + 1 * Ji_padded[0, 1] - P1[1] = - 1 * Ji_padded[1, 0] + 1 * Ji_padded[1, 1] - P2[0] = + 1 * Ji_padded[0, 0] + 1 * Ji_padded[0, 1] - P2[1] = + 1 * Ji_padded[1, 0] + 1 * Ji_padded[1, 1] - P3[0] = - 1 * Ji_padded[0, 0] - 1 * Ji_padded[0, 1] - P3[1] = - 1 * Ji_padded[1, 0] - 1 * Ji_padded[1, 1] - P4[0] = + 1 * Ji_padded[0, 0] - 1 * Ji_padded[0, 1] - P4[1] = + 1 * Ji_padded[1, 0] - 1 * Ji_padded[1, 1] - - # Find a bounding box around the transformed coordinates. - # (Check all four points at each step, since sometimes negative - # Jacobian values will mirror the transformed pixel.) - top = max(P1[1], P2[1], P3[1], P4[1]) - bottom = min(P1[1], P2[1], P3[1], P4[1]) - right = max(P1[0], P2[0], P3[0], P4[0]) - left = min(P1[0], P2[0], P3[0], P4[0]) + # We do that by defining a square region in the output plane + # centered on the output pixel, and transforming its corners to + # the input plane (using the local linearization of the + # transformation). Those transformed coordinates will set our + # bounding box. + if kernel_flag == 0: + # The Hann window is zero outside +/-1, so + # that's how far we need to go. + # + # The Hann window width is twice the width of a pixel---it + # runs to the centers of the neighboring pixels, rather + # than the edges of those pixels. This ensures that, at + # every point, the sum of the overlapping Hann windows is + # 1, and therefore that every input-image pixel is fully + # distributed into some combination of output pixels (in + # the limit of a Jacobian that is constant across all + # output pixels). + P1[0] = - 1 * Ji_padded[0, 0] + 1 * Ji_padded[0, 1] + P1[1] = - 1 * Ji_padded[1, 0] + 1 * Ji_padded[1, 1] + P2[0] = + 1 * Ji_padded[0, 0] + 1 * Ji_padded[0, 1] + P2[1] = + 1 * Ji_padded[1, 0] + 1 * Ji_padded[1, 1] + P3[0] = - 1 * Ji_padded[0, 0] - 1 * Ji_padded[0, 1] + P3[1] = - 1 * Ji_padded[1, 0] - 1 * Ji_padded[1, 1] + P4[0] = + 1 * Ji_padded[0, 0] - 1 * Ji_padded[0, 1] + P4[1] = + 1 * Ji_padded[1, 0] - 1 * Ji_padded[1, 1] + + # Find a bounding box around the transformed coordinates. + # (Check all four points at each step, in case a negative + # Jacobian value is mirroring the transformed pixel.) + top = max(P1[1], P2[1], P3[1], P4[1]) + bottom = min(P1[1], P2[1], P3[1], P4[1]) + right = max(P1[0], P2[0], P3[0], P4[0]) + left = min(P1[0], P2[0], P3[0], P4[0]) + elif kernel_flag == 1: + # The Gaussian window is non-zero everywhere, but it's + # close to zero almost everywhere. Sampling the whole input + # image isn't tractable, so we truncate and sample only + # within a certain region. + # n.b. `s` currently contains the reciprocal of the + # singular values + top = sample_region_width / (2 * min(s[0], s[1])) + bottom = -top + right = top + left = -right + else: + with gil: + raise ValueError("Invalid kernel type") if max_samples_width > 0 and max(right-left, top-bottom) > max_samples_width: if singularities_nan: @@ -388,7 +413,17 @@ def map_coordinates(double[:,:] source, double[:,:] target, Ci, int max_samples_ # Compute an averaging weight to be assigned to this # input location. - weight = hanning_filter(transformed[0], transformed[1]) + if kernel_flag == 0: + weight = hanning_filter( + transformed[0], transformed[1]) + elif kernel_flag == 1: + weight = gaussian_filter( + transformed[0], + transformed[1], + kernel_width) + else: + with gil: + raise ValueError("Invalid kernel type") if weight == 0: # As we move along each row in the image, we'll # first be seeing input-plane pixels that don't map diff --git a/reproject/adaptive/high_level.py b/reproject/adaptive/high_level.py index 7fa3e3e59..44076f0fc 100644 --- a/reproject/adaptive/high_level.py +++ b/reproject/adaptive/high_level.py @@ -8,7 +8,8 @@ def reproject_adaptive(input_data, output_projection, shape_out=None, hdu_in=0, return_footprint=True, center_jacobian=False, - roundtrip_coords=True, conserve_flux=False): + roundtrip_coords=True, conserve_flux=False, + kernel='Hann', kernel_width=1.3, sample_region_width=4): """ Reproject celestial slices from an 2d array from one WCS to another using the DeForest (2004) adaptive resampling algorithm. @@ -61,24 +62,28 @@ def reproject_adaptive(input_data, output_projection, shape_out=None, hdu_in=0, directions. conserve_flux : bool Whether to rescale output pixel values so flux is conserved. - kernel: str + kernel : str The averaging kernel to use. Allowed values are 'Hann' and 'Gaussian'. Case-insensitive. The Gaussian kernel produces better photometric accuracy at the cost of some blurring (on the scale of a few pixels). - kernel_width: double - The width of the kernel in pixels, measuring to the edge of the Hann - window or to +/- 1 sigma for the Gaussian window. If negative, a - default appropriate to the chosen kernel is used (2 for Hann, 1.5 for - Gaussian). Reducing this width may introduce photometric errors, while - increasing it will increase blurring of the output image. - sample_region_width: double + kernel_width : double + The width of the kernel in pixels, measuring to +/- 1 sigma for the + Gaussian window. Does not apply to the Hann window. Reducing this width + may introduce photometric errors or leave input pixels under-sampled, + while increasing it may improve the degree of anti-aliasing but will + increase blurring of the output image. If this width is changed from + the default, a proportional change should be made to the value of + sample_region_width to maintain an equivalent degree of photometric + accuracy. + sample_region_width : double The width in pixels of the output-image region which, when transformed to the input plane, defines the region to be sampled for each output pixel. Used only for the Gaussian kernel, which otherwise has infinite extent. This value sets a trade-off between accuracy and computation - time, with better accuracy at higher values. The default value of 4 may - yield errors of up to a few percent in some cases, while a value of 5 - should ensure errors are much less than a percent. + time, with better accuracy at higher values. The default value of 4, + with the default kernel width, should limit the most extreme errors to + less than one percent. Higher values will offer even more photometric + accuracy. Returns ------- @@ -99,4 +104,6 @@ def reproject_adaptive(input_data, output_projection, shape_out=None, hdu_in=0, return_footprint=return_footprint, center_jacobian=center_jacobian, roundtrip_coords=roundtrip_coords, - conserve_flux=conserve_flux) + conserve_flux=conserve_flux, + kernel=kernel, kernel_width=kernel_width, + sample_region_width=sample_region_width) diff --git a/reproject/tests/test_high_level.py b/reproject/tests/test_high_level.py index 9a7d195ef..e0f7453e5 100644 --- a/reproject/tests/test_high_level.py +++ b/reproject/tests/test_high_level.py @@ -17,7 +17,8 @@ 'biquadratic', 'bicubic', 'flux-conserving', - 'adaptive') + 'adaptive-hann', + 'adaptive-gaussian') ALL_DTYPES = [] for endian in ('<', '>'): @@ -124,7 +125,8 @@ def test_surface_brightness(projection_type, dtype): if projection_type == 'flux-conserving': data_out, footprint = reproject_exact((data_in, header_in), header_out) elif projection_type.startswith('adaptive'): - data_out, footprint = reproject_adaptive((data_in, header_in), header_out) + data_out, footprint = reproject_adaptive((data_in, header_in), header_out, + kernel=projection_type.split('-', 1)[1]) else: data_out, footprint = reproject_interp((data_in, header_in), header_out, order=projection_type) @@ -164,7 +166,8 @@ def test_identity_projection(projection_type): if projection_type == 'flux-conserving': data_out, footprint = reproject_exact((data_in, header_in), header_in) elif projection_type.startswith('adaptive'): - data_out, footprint = reproject_adaptive((data_in, header_in), header_in) + data_out, footprint = reproject_adaptive((data_in, header_in), header_in, + kernel=projection_type.split('-', 1)[1]) else: data_out, footprint = reproject_interp((data_in, header_in), header_in, order=projection_type) @@ -173,4 +176,8 @@ def test_identity_projection(projection_type): # and the footprint values to be ~ones. expected_footprint = np.ones((header_in['NAXIS2'], header_in['NAXIS1'])) np.testing.assert_allclose(footprint, expected_footprint) - np.testing.assert_allclose(data_in, data_out, rtol=1e-6) + if projection_type != 'adaptive-gaussian': + np.testing.assert_allclose(data_in, data_out, rtol=1e-6) + else: + # The Gaussian kernel in the adaptive algorithm blurs the image + assert np.all(data_in != data_out)