Skip to content

Commit

Permalink
Add Gaussian filter kernel for adaptive resampling
Browse files Browse the repository at this point in the history
  • Loading branch information
@svank
svank committed May 21, 2022
1 parent a5aee3f commit 960da0b
Show file tree
Hide file tree
Showing 5 changed files with 166 additions and 76 deletions.
54 changes: 42 additions & 12 deletions docs/celestial.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
<https://doi.org/10.1023/B:SOLA.0000021743.24248.b0>`_ 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
Expand All @@ -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::
Expand Down
15 changes: 13 additions & 2 deletions reproject/adaptive/core.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
-------
Expand Down Expand Up @@ -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)
Expand Down
125 changes: 80 additions & 45 deletions reproject/adaptive/deforest.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -163,14 +163,27 @@ cdef double sample_array(double[:,:] source, double x, double y, int x_cyclic,
return source[<int> y, <int> x]


KERNELS = {}
KERNELS['hann'] = 0
KERNELS['hanning'] = KERNELS['hann']
KERNELS['gaussian'] = 1


@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
@cython.cdivision(True)
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:
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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
Expand Down
33 changes: 20 additions & 13 deletions reproject/adaptive/high_level.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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
-------
Expand All @@ -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)
Loading

0 comments on commit 960da0b

Please sign in to comment.