Skip to content

Commit

Permalink
Add boundary handling modes for adaptive algorithm
Browse files Browse the repository at this point in the history
  • Loading branch information
@svank
svank committed Aug 23, 2022
1 parent 1df0d8e commit 9a0575b
Show file tree
Hide file tree
Showing 5 changed files with 472 additions and 42 deletions.
31 changes: 31 additions & 0 deletions docs/celestial.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -243,6 +243,37 @@ 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.

When, for any one output pixel, the sampling region in the input image
straddles the boundary of the input image or lies entirely outside the input
image, a range of boundary modes can be applied, and this is set with the
``boundary_mode`` option. Allowed values are:

* ``strict`` --- Output pixels will be ``NaN`` if any of their input samples
fall outside the input image.
* ``constant`` --- Samples outside the bounds of the input image are
replaced by a constant value, set with the ``boundary_fill_value`` argument.
Output values become ``NaN`` if there are no valid input samples.
* ``grid-constant`` --- Samples outside the bounds of the input image are
replaced by a constant value, set with the ``boundary_fill_value`` argument.
Output values will be ``boundary_fill_value`` if there are no valid input
samples.
* ``ignore`` --- Samples outside the input image are simply ignored,
contributing neither to the output value nor the sum-of-weights
normalization. If there are no valid input samples, the output value will be
``NaN``.
* ``ignore_threshold`` --- Acts as ``ignore``, unless the total weight that
would have been assigned to the ignored samples exceeds a set fraction of the
total weight across the entire sampling region, set by the
``boundary_ignore_threshold`` argument. In that case, acts as ``strict``.
* ``nearest`` --- Samples outside the input image are replaced by the nearst
in-bounds input pixel.

The input image can also be marked as being cyclic or periodic in the x and/or
y axes with the ``x_cyclic`` and ``y_cyclic`` flags. If these are set, samples
will wrap around to the opposite side of the image, ignoring the
``boundary_mode`` for that axis.


Algorithm Description
---------------------

Expand Down
19 changes: 17 additions & 2 deletions reproject/adaptive/core.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,10 @@ def _reproject_adaptive_2d(array, wcs_in, wcs_out, shape_out,
return_footprint=True, center_jacobian=False,
roundtrip_coords=True, conserve_flux=False,
kernel='Hann', kernel_width=1.3,
sample_region_width=4):
sample_region_width=4,
boundary_mode='ignore', boundary_fill_value=0,
boundary_ignore_threshold=0.5,
x_cyclic=False, y_cyclic=False):
"""
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 Down Expand Up @@ -67,6 +70,14 @@ def _reproject_adaptive_2d(array, wcs_in, wcs_out, shape_out,
sample_region_width : double
The width in pixels of the sample region, used only for the Gaussian
kernel which otherwise has infinite extent.
boundary_mode : str
Boundary handling mode
boundary_fill_value : double
Fill value for 'constant' boundary mode
boundary_ignore_threshold : double
Threshold for 'ignore_threshold' boundary mode, ranging from 0 to 1.
x_cyclic, y_cyclic : bool
Marks in input-image axis as cyclic.
Returns
-------
Expand Down Expand Up @@ -107,7 +118,11 @@ def _reproject_adaptive_2d(array, wcs_in, wcs_out, shape_out,
map_coordinates(array_in, array_out, transformer, out_of_range_nan=True,
center_jacobian=center_jacobian, conserve_flux=conserve_flux,
kernel=kernel, kernel_width=kernel_width,
sample_region_width=sample_region_width)
sample_region_width=sample_region_width,
boundary_mode=boundary_mode,
boundary_fill_value=boundary_fill_value,
boundary_ignore_threshold=boundary_ignore_threshold,
x_cyclic=x_cyclic, y_cyclic=y_cyclic)

if return_footprint:
return array_out, (~np.isnan(array_out)).astype(float)
Expand Down
174 changes: 136 additions & 38 deletions reproject/adaptive/deforest.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -128,47 +128,55 @@ cdef double gaussian_filter(double x, double y, double width) nogil:
@cython.wraparound(False)
@cython.nonecheck(False)
@cython.cdivision(True)
cdef double clip(double x, double vmin, double vmax, int cyclic, int out_of_range_nan) nogil:
cdef double clip(double x, double vmin, double vmax, int cyclic,
int out_of_range_nearest) nogil:
if x < vmin:
if cyclic:
while x < vmin:
x += (vmax-vmin)+1
elif out_of_range_nan:
return nan
else:
elif out_of_range_nearest:
return vmin
else:
return nan
elif x > vmax:
if cyclic:
while x > vmax:
x -= (vmax-vmin)+1
elif out_of_range_nan:
return nan
else:
elif out_of_range_nearest:
return vmax
else:
return x
else:
return nan
return x


@cython.boundscheck(False)
@cython.wraparound(False)
@cython.nonecheck(False)
@cython.cdivision(True)
cdef double sample_array(double[:,:] source, double x, double y, int x_cyclic,
int y_cyclic, int out_of_range_nan) nogil:
x = clip(x, 0, source.shape[1] - 1, x_cyclic, out_of_range_nan)
y = clip(y, 0, source.shape[0] - 1, y_cyclic, out_of_range_nan)
cdef (double, bint) sample_array(double[:,:] source, double x, double y,
int x_cyclic, int y_cyclic, bint out_of_range_nearest) nogil:
x = clip(x, 0, source.shape[1] - 1, x_cyclic, out_of_range_nearest)
y = clip(y, 0, source.shape[0] - 1, y_cyclic, out_of_range_nearest)

if isnan(x) or isnan(y):
return nan
return nan, False

return source[<int> y, <int> x]
return source[<int> y, <int> x], True


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

BOUNDARY_MODES = {}
BOUNDARY_MODES['strict'] = 1
BOUNDARY_MODES['constant'] = 2
BOUNDARY_MODES['grid-constant'] = 3
BOUNDARY_MODES['ignore'] = 4
BOUNDARY_MODES['ignore_threshold'] = 5
BOUNDARY_MODES['nearest'] = 6


@cython.boundscheck(False)
@cython.wraparound(False)
Expand All @@ -178,13 +186,21 @@ def map_coordinates(double[:,:] source, double[:,:] target, Ci, int max_samples_
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, str kernel='Hann', double kernel_width=1.3,
double sample_region_width=4):
double sample_region_width=4, str boundary_mode="ignore",
double boundary_fill_value=0, double boundary_ignore_threshold=0.5):
cdef int kernel_flag
try:
kernel_flag = KERNELS[kernel.lower()]
except KeyError:
raise ValueError("'kernel' must be 'Hann' or 'Gaussian'")

cdef int boundary_flag
try:
boundary_flag = BOUNDARY_MODES[boundary_mode.lower()]
except KeyError:
raise ValueError(
f"boundary_mode '{boundary_mode}' not recognized") from None

cdef np.ndarray[np.float64_t, ndim=3] pixel_target
cdef int delta
if center_jacobian:
Expand Down Expand Up @@ -274,7 +290,8 @@ def map_coordinates(double[:,:] source, double[:,:] target, Ci, int max_samples_
cdef double[:] transformed = np.zeros((2,))
cdef double[:] current_pixel_source = np.zeros((2,))
cdef double[:] current_offset = np.zeros((2,))
cdef double weight_sum = 0.0
cdef double weight_sum
cdef double ignored_weight_sum
cdef double weight
cdef double value
cdef double[:] P1 = np.empty((2,))
Expand All @@ -283,6 +300,7 @@ def map_coordinates(double[:,:] source, double[:,:] target, Ci, int max_samples_
cdef double[:] P4 = np.empty((2,))
cdef double top, bottom, left, right
cdef bint has_sampled_this_row
cdef bint is_good_sample
with nogil:
# Iterate through each pixel in the output image.
for yi in range(target.shape[0]):
Expand Down Expand Up @@ -376,35 +394,100 @@ def map_coordinates(double[:,:] source, double[:,:] target, Ci, int max_samples_
if singularities_nan:
target[yi,xi] = nan
else:
target[yi,xi] = sample_array(source,
pixel_source[yi,xi,0], pixel_source[yi,xi,1],
x_cyclic, y_cyclic, out_of_range_nan)
value, is_good_sample = sample_array(
source, current_pixel_source[0],
current_pixel_source[1], x_cyclic, y_cyclic,
out_of_range_nearest=boundary_flag == 6)
if is_good_sample:
target[yi,xi] = value
elif boundary_flag == 2 or boundary_flag == 3:
target[yi,xi] = boundary_fill_value
else:
target[yi,xi] = nan
continue

top += pixel_source[yi,xi,1]
bottom += pixel_source[yi,xi,1]
right += pixel_source[yi,xi,0]
left += pixel_source[yi,xi,0]

# Clamp sampling region to stay within the source image. (Going
# outside the image plane is handled otherwise, but it's faster
# to just omit those points completely.)
# Draw these points in to the nearest input pixel
bottom = ceil(bottom)
top = floor(top)
left = ceil(left)
right = floor(right)

# Handle the case that the sampling region extends beyond the
# input image boundary. For 'strict' handling, we can set the
# output pixel to NaN right away. For 'ignore' handling, we can
# clamp the region to exclude the out-of-bounds samples. For
# all other boundary modes, we still need to calculate weights
# for each out-of-bounds sample, so we do nothing here.
if not x_cyclic:
right = min(source.shape[1] - 1, right)
left = max(0, left)
if right > source.shape[1] - 1:
if boundary_flag == 1:
target[yi,xi] = nan
continue
if boundary_flag == 4:
right = source.shape[1] - 1
if left < 0:
if boundary_flag == 1:
target[yi,xi] = nan
continue
if boundary_flag == 4:
left = 0
if not y_cyclic:
top = min(source.shape[0] - 1, top)
bottom = max(0, bottom)

target[yi,xi] = 0.0
weight_sum = 0.0
if top > source.shape[0] - 1:
if boundary_flag == 1:
target[yi,xi] = nan
continue
if boundary_flag == 4:
top = source.shape[0] - 1
if bottom < 0:
if boundary_flag == 1:
target[yi,xi] = nan
continue
if boundary_flag == 4:
bottom = 0

# Check whether the sampling region falls entirely outside the
# input image. For strict boundary handling, this is already
# handled above by the partial case. Otherwise, we fill in an
# appropriate value and move along. For some projections, the
# sampling region can become very large when well outside the
# input image, and so this detection becomes an important
# optimization.
if (not x_cyclic and (right < 0 or left > source.shape[1] - 1)
or not y_cyclic
and (top < 0 or bottom > source.shape[0] - 1)):
if boundary_flag == 3:
target[yi,xi] = boundary_fill_value
continue
if (boundary_flag == 2
or boundary_flag == 4
or boundary_flag == 5):
target[yi,xi] = nan
continue
if boundary_flag == 6:
# Just sample one row or column so that we get all of
# the nearest values. Both kernels vary independently
# in x and y, so sampling the full region isn't needed
# when the sampled values are constant in x or in y.
if right < left:
right = left
if top < bottom:
top = bottom

target[yi,xi] = 0
weight_sum = 0
ignored_weight_sum = 0

# Iterate through that bounding box in the input image.
for y in range(<int>ceil(bottom), <int>floor(top)+1):
for y in range(<int> bottom, <int> top+1):
current_pixel_source[1] = y
current_offset[1] = current_pixel_source[1] - pixel_source[yi,xi,1]
has_sampled_this_row = False
for x in range(<int>ceil(left), <int>floor(right)+1):
for x in range(<int> left, <int> right+1):
current_pixel_source[0] = x
current_offset[0] = current_pixel_source[0] - pixel_source[yi,xi,0]
# Find the fractional position of the input location
Expand Down Expand Up @@ -443,16 +526,31 @@ def map_coordinates(double[:,:] source, double[:,:] target, Ci, int max_samples_
continue
has_sampled_this_row = True

value = sample_array(source, current_pixel_source[0],
value, is_good_sample = sample_array(
source, current_pixel_source[0],
current_pixel_source[1], x_cyclic, y_cyclic,
out_of_range_nan)
out_of_range_nearest=(boundary_flag == 6))

if ((boundary_flag == 2 or boundary_flag == 3)
and not is_good_sample):
value = boundary_fill_value
is_good_sample = True

if not isnan(value):
if is_good_sample:
target[yi,xi] += weight * value
weight_sum += weight
target[yi,xi] /= weight_sum
if conserve_flux:
target[yi,xi] *= fabs(det2x2(Ji))
else:
if boundary_flag == 5:
ignored_weight_sum += weight

if (boundary_flag == 5 and
ignored_weight_sum / (ignored_weight_sum + weight_sum)
> boundary_ignore_threshold):
target[yi,xi] = nan
else:
target[yi,xi] /= weight_sum
if conserve_flux:
target[yi,xi] *= fabs(det2x2(Ji))
if progress:
with gil:
sys.stdout.write("\r%d/%d done" % (yi+1, target.shape[0]))
Expand Down
Loading

0 comments on commit 9a0575b

Please sign in to comment.