Add modes for nan and inf handling to adaptive algo

docs/celestial.rst

@@ -173,7 +173,7 @@ integer or a string giving the order of the interpolation. Supported strings
 include:
 
 * ``'nearest-neighbor'``: zeroth order interpolation
-* ``'bilinear'``: fisst order interpolation
+* ``'bilinear'``: first order interpolation
 * ``'biquadratic'``: second order interpolation
 * ``'bicubic'``: third order interpolation
 
@@ -279,14 +279,25 @@ image, a range of boundary modes can be applied, and this is set with the
   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
+* ``nearest`` --- Samples outside the input image are replaced by the nearest
   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.
 
+This implementation includes several options for handling ``nan`` and ``inf``
+values in the input data, set via the ``bad_value_mode`` argument:
+
+* ``strict`` --- Values of ``nan`` or ``inf`` in the input data are propagated
+  to every output value which samples them.
+* ``ignore`` --- When a sampled input value is ``nan`` or ``inf``, that input
+  pixel is ignored (affected neither the accumulated sum of weighted samples
+  nor the accumulated sum of weights).
+* ``constant`` --- Input values of ``nan`` and ``inf`` are replaced with a
+  constant value, set via the ``bad_fill_value`` argument.
+
 
 Algorithm Description
 ---------------------

reproject/adaptive/core.py

@@ -45,6 +45,8 @@ def _reproject_adaptive_2d(
     boundary_ignore_threshold=0.5,
     x_cyclic=False,
     y_cyclic=False,
+    bad_value_mode="strict",
+    bad_fill_value=0,
 ):
     """
     Reproject celestial slices from an n-d array from one WCS to another
@@ -87,6 +89,10 @@ def _reproject_adaptive_2d(
         Threshold for 'ignore_threshold' boundary mode, ranging from 0 to 1.
     x_cyclic, y_cyclic : bool
         Marks in input-image axis as cyclic.
+    bad_value_mode : str
+        NaN and inf handling mode
+    bad_fill_value : float
+        Fill value for 'constant' bad value mode
 
     Returns
     -------
@@ -167,6 +173,8 @@ def _reproject_adaptive_2d(
         boundary_ignore_threshold=boundary_ignore_threshold,
         x_cyclic=x_cyclic,
         y_cyclic=y_cyclic,
+        bad_value_mode=bad_value_mode,
+        bad_fill_value=bad_fill_value,
     )
 
     array_out.shape = shape_out

reproject/adaptive/deforest.pyx

@@ -43,6 +43,7 @@ cdef double nan = np.nan
 
 cdef extern from "math.h":
     int isnan(double x) nogil
+    int isinf(double x) nogil
 
 
 @cython.boundscheck(False)
@@ -134,6 +135,12 @@ cdef double gaussian_filter(double x, double y, double width) nogil:
 @cython.cdivision(True)
 cdef double clip(double x, double vmin, double vmax, int cyclic,
         int out_of_range_nearest) nogil:
+    """Applies bounary conditions to an intended array coordinate.
+
+    Specifically, if the point is outside the array bounds, this function wraps
+    the coordinate if the boundary is periodic, or clamps to the nearest valid
+    coordinate if desired, or else returns NaN.
+    """
     if x < vmin:
         if cyclic:
             while x < vmin:
@@ -164,6 +171,8 @@ cdef bint sample_array(double[:,:,:] source, double[:] dest,
     y = clip(y, 0, source.shape[1] - 1, y_cyclic, out_of_range_nearest)
 
     if isnan(x) or isnan(y):
+        # Indicates the coordinate is outside the array's bounds and the
+        # boundary-handling mode doesn't provide an alternative coordinate.
         return False
 
     # Cython doesn't like a return type of (double[:], bint), so we put the
@@ -343,6 +352,10 @@ BOUNDARY_MODES['ignore'] = 4
 BOUNDARY_MODES['ignore_threshold'] = 5
 BOUNDARY_MODES['nearest'] = 6
 
+BAD_VALUE_MODES = {}
+BAD_VALUE_MODES['strict'] = 1
+BAD_VALUE_MODES['constant'] = 2
+BAD_VALUE_MODES['ignore'] = 3
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
@@ -355,6 +368,7 @@ def map_coordinates(double[:,:,:] source, double[:,:,:] target, Ci, int max_samp
                     str kernel='gaussian', double kernel_width=1.3,
                     double sample_region_width=4, str boundary_mode="strict",
                     double boundary_fill_value=0, double boundary_ignore_threshold=0.5,
+                    str bad_value_mode="strict", double bad_fill_value=0,
                     ):
     # n.b. the source and target arrays are expected to contain three
     # dimensions---the last two are the image dimensions, while the first
@@ -375,6 +389,13 @@ def map_coordinates(double[:,:,:] source, double[:,:,:] target, Ci, int max_samp
         raise ValueError(
                 f"boundary_mode '{boundary_mode}' not recognized") from None
 
+    cdef int bad_val_flag
+    try:
+        bad_val_flag = BAD_VALUE_MODES[bad_value_mode.lower()]
+    except KeyError:
+        raise ValueError(
+                f"bad_value_mode '{bad_value_mode}' not recognized") from None
+
     cdef np.ndarray[np.float64_t, ndim=3] pixel_target
     cdef int delta
     if center_jacobian:
@@ -477,7 +498,7 @@ def map_coordinates(double[:,:,:] source, double[:,:,:] target, Ci, int max_samp
     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
+    cdef double[:] weight_sum = np.empty(source.shape[0])
     cdef double ignored_weight_sum
     cdef double weight
     cdef double[:] value = np.empty(source.shape[0])
@@ -488,7 +509,7 @@ def map_coordinates(double[:,:,:] source, double[:,:,:] target, Ci, int max_samp
     cdef double top, bottom, left, right
     cdef double determinant
     cdef bint has_sampled_this_row
-    cdef bint is_good_sample
+    cdef bint sample_in_bounds
     with nogil:
         # Iterate through each pixel in the output image.
         for yi in range(target.shape[1]):
@@ -572,12 +593,19 @@ def map_coordinates(double[:,:,:] source, double[:,:,:] target, Ci, int max_samp
                     if singularities_nan:
                         target[:,yi,xi] = nan
                     else:
-                        is_good_sample = sample_array(
+                        sample_in_bounds = sample_array(
                                 source, value, 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
+                        if sample_in_bounds:
+                            for i in range(target.shape[0]):
+                                if bad_val_flag != 1 and (isnan(value[i]) or isinf(value[i])):
+                                    if bad_val_flag == 2:
+                                        target[i,yi,xi] = bad_fill_value
+                                    else:
+                                        target[i,yi,xi] = nan
+                                else:
+                                    target[i,yi,xi] = value[i]
                         elif boundary_flag == 2 or boundary_flag == 3:
                             target[:,yi,xi] = boundary_fill_value
                         else:
@@ -657,7 +685,7 @@ def map_coordinates(double[:,:,:] source, double[:,:,:] target, Ci, int max_samp
                             top = bottom
 
                 target[:,yi,xi] = 0
-                weight_sum = 0
+                weight_sum[:] = 0
                 ignored_weight_sum = 0
 
                 # Iterate through that bounding box in the input image.
@@ -704,35 +732,41 @@ def map_coordinates(double[:,:,:] source, double[:,:,:] target, Ci, int max_samp
                             continue
                         has_sampled_this_row = True
 
-                        is_good_sample = sample_array(
+                        sample_in_bounds = sample_array(
                                 source, value, current_pixel_source[0],
                                 current_pixel_source[1], x_cyclic, y_cyclic,
                                 out_of_range_nearest=(boundary_flag == 6))
 
                         if ((boundary_flag == 2 or boundary_flag == 3)
-                                and not is_good_sample):
+                                and not sample_in_bounds):
                             value[:] = boundary_fill_value
-                            is_good_sample = True
+                            sample_in_bounds = True
 
-                        if is_good_sample:
+                        if sample_in_bounds:
                             for i in range(target.shape[0]):
+                                if bad_val_flag != 1 and (isnan(value[i]) or isinf(value[i])):
+                                    if bad_val_flag == 2:
+                                        value[i] = bad_fill_value
+                                    else:
+                                        # bad_val_flag is 3: 'ignore'
+                                        continue
                                 target[i,yi,xi] += weight * value[i]
-                            weight_sum += weight
+                                weight_sum[i] += weight
                         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:
-                    if conserve_flux:
-                        determinant = fabs(det2x2(Ji))
+                if boundary_flag == 5:
                     for i in range(target.shape[0]):
-                        target[i,yi,xi] /= weight_sum
-                        if conserve_flux:
-                            target[i,yi,xi] *= determinant
+                        if (ignored_weight_sum / (ignored_weight_sum + weight_sum[i])
+                                > boundary_ignore_threshold):
+                            target[i,yi,xi] = nan
+                if conserve_flux:
+                    determinant = fabs(det2x2(Ji))
+                for i in range(target.shape[0]):
+                    target[i,yi,xi] /= weight_sum[i]
+                    if conserve_flux:
+                        target[i,yi,xi] *= determinant
             if progress:
                 with gil:
                     sys.stdout.write("\r%d/%d done" % (yi+1, target.shape[1]))

reproject/adaptive/high_level.py

@@ -30,6 +30,8 @@ def reproject_adaptive(
     boundary_ignore_threshold=0.5,
     x_cyclic=False,
     y_cyclic=False,
+    bad_value_mode="strict",
+    bad_fill_value=0,
 ):
     """
     Reproject a 2D array from one WCS to another using the DeForest (2004)
@@ -169,6 +171,20 @@ def reproject_adaptive(
         Indicates that the x or y axis of the input image should be treated as
         cyclic or periodic. Overrides the boundary mode for that axis, so that
         out-of-bounds samples wrap to the other side of the image.
+    bad_value_mode : str
+        How to handle values of ``nan`` and ``inf`` in the input data. The
+        default is ``strct``. Allowed values are:
+
+            * ``strict`` --- Values of ``nan`` or ``inf`` in the input data are
+              propagated to every output value which samples them.
+            * ``ignore`` --- When a sampled input value is ``nan`` or ``inf``,
+              that input pixel is ignored (affected neither the accumulated sum
+              of weighted samples nor the accumulated sum of weights).
+            * ``constant`` --- Input values of ``nan`` and ``inf`` are replaced
+              with a constant value, set via the ``bad_fill_value`` argument.
+
+    bad_fill_value : double
+        The constant value used by the ``constant`` bad-value mode.
 
     Returns
     -------
@@ -211,6 +227,8 @@ def reproject_adaptive(
             boundary_ignore_threshold=boundary_ignore_threshold,
             x_cyclic=x_cyclic,
             y_cyclic=y_cyclic,
+            bad_value_mode=bad_value_mode,
+            bad_fill_value=bad_fill_value,
         ),
         return_type=return_type,
     )