Remove extra interpolation step in adaptive resampling

docs/celestial.rst

@@ -186,16 +186,9 @@ reprojection using the  `DeForest (2004)
 
     >>> from reproject import reproject_adaptive
 
-In addition to the arguments described in :ref:`common`, the order of the
-interpolation to use when sampling values in the input map can be controlled by
-setting the ``order=`` argument to either an integer or a string giving the
-order of the interpolation. Supported strings include:
-
-* ``'nearest-neighbor'``: zeroth order interpolation
-* ``'bilinear'``: first order interpolation
-
-One can also enable a rescaling of output pixel values to conserve flux by
-using the ``conserve_flux`` flag.
+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

reproject/adaptive/core.py

@@ -31,7 +31,7 @@ def __call__(self, pixel_out):
         return pixel_in
 
 
-def _reproject_adaptive_2d(array, wcs_in, wcs_out, shape_out, order=1,
+def _reproject_adaptive_2d(array, wcs_in, wcs_out, shape_out,
                            return_footprint=True, center_jacobian=False,
                            roundtrip_coords=True, conserve_flux=False):
     """
@@ -49,8 +49,6 @@ def _reproject_adaptive_2d(array, wcs_in, wcs_out, shape_out, order=1,
         The output WCS
     shape_out : tuple
         The shape of the output array
-    order : int, optional
-        The order of the interpolation.
     return_footprint : bool
         Whether to return the footprint in addition to the output array.
     center_jacobian : bool
@@ -98,8 +96,7 @@ def _reproject_adaptive_2d(array, wcs_in, wcs_out, shape_out, order=1,
 
     transformer = CoordinateTransformer(wcs_in, wcs_out, roundtrip_coords)
     map_coordinates(array_in, array_out, transformer, out_of_range_nan=True,
-                    order=order, center_jacobian=center_jacobian,
-                    conserve_flux=conserve_flux)
+                    center_jacobian=center_jacobian, conserve_flux=conserve_flux)
 
     if return_footprint:
         return array_out, (~np.isnan(array_out)).astype(float)

reproject/adaptive/deforest.pyx

@@ -152,41 +152,15 @@ cdef double clip(double x, double vmin, double vmax, int cyclic, int out_of_rang
 @cython.wraparound(False)
 @cython.nonecheck(False)
 @cython.cdivision(True)
-cdef double bilinear_interpolation(double[:,:] source, double x, double y, int x_cyclic, int y_cyclic, int out_of_range_nan) nogil:
-
-    x = clip(x, -0.5, source.shape[1] - 0.5, x_cyclic, out_of_range_nan)
-    y = clip(y, -0.5, source.shape[0] - 0.5, y_cyclic, out_of_range_nan)
+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)
 
     if isnan(x) or isnan(y):
         return nan
 
-    cdef int xmin = <int>floor(x)
-    cdef int ymin = <int>floor(y)
-    cdef int xmax = xmin + 1
-    cdef int ymax = ymin + 1
-
-    cdef double fQ11 = source[max(0, ymin), max(0, xmin)]
-    cdef double fQ21 = source[max(0, ymin), min(source.shape[1] - 1, xmax)]
-    cdef double fQ12 = source[min(source.shape[0] - 1, ymax), max(0, xmin)]
-    cdef double fQ22 = source[min(source.shape[0] - 1, ymax), min(source.shape[1] - 1, xmax)]
-
-    return ((fQ11 * (xmax - x) * (ymax - y)
-             + fQ21 * (x - xmin) * (ymax - y)
-             + fQ12 * (xmax - x) * (y - ymin)
-             + fQ22 * (x - xmin) * (y - ymin))
-            * ((xmax - xmin) * (ymax - ymin)))
-
-
-@cython.boundscheck(False)
-@cython.wraparound(False)
-@cython.nonecheck(False)
-@cython.cdivision(True)
-cdef double nearest_neighbour_interpolation(double[:,:] source, double x, double y, int x_cyclic, int y_cyclic, int out_of_range_nan) nogil:
-    y = clip(round(y), 0, source.shape[0]-1, y_cyclic, out_of_range_nan)
-    x = clip(round(x), 0, source.shape[1]-1, x_cyclic, out_of_range_nan)
-    if isnan(y) or isnan(x):
-        return nan
-    return source[<int>y, <int>x]
+    return source[<int> y, <int> x]
 
 
 @cython.boundscheck(False)
@@ -196,7 +170,7 @@ cdef double nearest_neighbour_interpolation(double[:,:] source, double x, double
 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,
-                    int order=0, bint center_jacobian=False):
+                    bint center_jacobian=False):
     cdef np.ndarray[np.float64_t, ndim=3] pixel_target
     cdef int delta
     if center_jacobian:
@@ -210,7 +184,7 @@ def map_coordinates(double[:,:] source, double[:,:] target, Ci, int max_samples_
         # be representing (-1,-1) in the output image.
         delta = -1
 
-    cdef int yi, xi, yoff, xoff
+    cdef int yi, xi, y, x
     for yi in range(pixel_target.shape[0]):
         for xi in range(pixel_target.shape[1]):
             pixel_target[yi,xi,0] = xi + delta
@@ -288,7 +262,7 @@ def map_coordinates(double[:,:] source, double[:,:] target, Ci, int max_samples_
     cdef double[:] current_offset = np.zeros((2,))
     cdef double weight_sum = 0.0
     cdef double weight
-    cdef double interpolated
+    cdef double value
     cdef double[:] P1 = np.empty((2,))
     cdef double[:] P2 = np.empty((2,))
     cdef double[:] P3 = np.empty((2,))
@@ -376,34 +350,37 @@ def map_coordinates(double[:,:] source, double[:,:] target, Ci, int max_samples_
                     if singularities_nan:
                         target[yi,xi] = nan
                     else:
-                        if order == 0:
-                            target[yi,xi] = nearest_neighbour_interpolation(source, pixel_source[yi,xi,0], pixel_source[yi,xi,1], x_cyclic, y_cyclic, out_of_range_nan)
-                        else:
-                            target[yi,xi] = bilinear_interpolation(source, pixel_source[yi,xi,0], pixel_source[yi,xi,1], x_cyclic, y_cyclic, out_of_range_nan)
+                        target[yi,xi] = sample_array(source,
+                                pixel_source[yi,xi,0], pixel_source[yi,xi,1],
+                                x_cyclic, y_cyclic, out_of_range_nan)
                     continue
 
-                # Clamp to the largest offsets that remain within the source
-                # image. (Going out-of-bounds in the image plane would be
-                # handled correctly in the interpolation routines, but skipping
-                # those pixels altogether is faster.)
+                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.)
                 if not x_cyclic:
-                    right = min(source.shape[1] - 0.5 - pixel_source[yi,xi,0], right)
-                    left = max(-0.5 - pixel_source[yi,xi,0], left)
+                    right = min(source.shape[1] - 1, right)
+                    left = max(0, left)
                 if not y_cyclic:
-                    top = min(source.shape[0] - 0.5 - pixel_source[yi,xi,1], top)
-                    bottom = max(-0.5 - pixel_source[yi,xi,1], bottom)
+                    top = min(source.shape[0] - 1, top)
+                    bottom = max(0, bottom)
 
                 target[yi,xi] = 0.0
                 weight_sum = 0.0
 
                 # Iterate through that bounding box in the input image.
-                for yoff in range(<int>ceil(bottom), <int>floor(top)+1):
-                    current_offset[1] = yoff
-                    current_pixel_source[1] = pixel_source[yi,xi,1] + yoff
+                for y in range(<int>ceil(bottom), <int>floor(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 xoff in range(<int>ceil(left), <int>floor(right)+1):
-                        current_offset[0] = xoff
-                        current_pixel_source[0] = pixel_source[yi,xi,0] + xoff
+                    for x in range(<int>ceil(left), <int>floor(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
                         # within the transformed ellipse.
                         transformed[0] = J[0,0] * current_offset[0] + J[0,1] * current_offset[1]
@@ -428,17 +405,14 @@ def map_coordinates(double[:,:] source, double[:,:] target, Ci, int max_samples_
                             if has_sampled_this_row:
                                 break
                             continue
+                        has_sampled_this_row = True
+
+                        value = sample_array(source, current_pixel_source[0],
+                                current_pixel_source[1], x_cyclic, y_cyclic,
+                                out_of_range_nan)
 
-                        # Produce an input-image value to sample. Our output
-                        # pixel doesn't necessarily map to an integer
-                        # coordinate in the input image, and so our input
-                        # samples must be interpolated.
-                        if order == 0:
-                            interpolated = nearest_neighbour_interpolation(source, current_pixel_source[0], current_pixel_source[1], x_cyclic, y_cyclic, out_of_range_nan)
-                        else:
-                            interpolated = bilinear_interpolation(source, current_pixel_source[0], current_pixel_source[1], x_cyclic, y_cyclic, out_of_range_nan)
-                        if not isnan(interpolated):
-                            target[yi,xi] += weight * interpolated
+                        if not isnan(value):
+                            target[yi,xi] += weight * value
                             weight_sum += weight
                 target[yi,xi] /= weight_sum
                 if conserve_flux:

reproject/adaptive/high_level.py

@@ -5,15 +5,10 @@
 
 __all__ = ['reproject_adaptive']
 
-ORDER = {}
-ORDER['nearest-neighbor'] = 0
-ORDER['bilinear'] = 1
-
 
 def reproject_adaptive(input_data, output_projection, shape_out=None, hdu_in=0,
-                       order='bilinear', return_footprint=True,
-                       center_jacobian=False, roundtrip_coords=True,
-                       conserve_flux=False):
+                       return_footprint=True, center_jacobian=False,
+                       roundtrip_coords=True, conserve_flux=False):
     """
     Reproject celestial slices from an 2d array from one WCS to another using
     the DeForest (2004) adaptive resampling algorithm.
@@ -43,15 +38,6 @@ def reproject_adaptive(input_data, output_projection, shape_out=None, hdu_in=0,
     hdu_in : int or str, optional
         If ``input_data`` is a FITS file or an `~astropy.io.fits.HDUList`
         instance, specifies the HDU to use.
-    order : int or str, optional
-        The order of the interpolation. This can be any of the
-        following strings:
-
-            * 'nearest-neighbor'
-            * 'bilinear'
-
-        or an integer. A value of ``0`` indicates nearest neighbor
-        interpolation.
     return_footprint : bool
         Whether to return the footprint in addition to the output array.
     center_jacobian : bool
@@ -109,11 +95,8 @@ def reproject_adaptive(input_data, output_projection, shape_out=None, hdu_in=0,
     array_in, wcs_in = parse_input_data(input_data, hdu_in=hdu_in)
     wcs_out, shape_out = parse_output_projection(output_projection, shape_out=shape_out)
 
-    if isinstance(order, str):
-        order = ORDER[order]
-
     return _reproject_adaptive_2d(array_in, wcs_in, wcs_out, shape_out,
-                                  order=order, return_footprint=return_footprint,
+                                  return_footprint=return_footprint,
                                   center_jacobian=center_jacobian,
                                   roundtrip_coords=roundtrip_coords,
                                   conserve_flux=conserve_flux)

reproject/tests/test_high_level.py

@@ -17,8 +17,7 @@
              'biquadratic',
              'bicubic',
              'flux-conserving',
-             'adaptive-nearest-neighbor',
-             'adaptive-bilinear')
+             'adaptive')
 
 ALL_DTYPES = []
 for endian in ('<', '>'):
@@ -125,8 +124,7 @@ 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,
-                                                 order=projection_type.split('-', 1)[1])
+        data_out, footprint = reproject_adaptive((data_in, header_in), header_out)
     else:
         data_out, footprint = reproject_interp((data_in, header_in), header_out,
                                                order=projection_type)
@@ -166,8 +164,7 @@ 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,
-                                                 order=projection_type.split('-', 1)[1])
+        data_out, footprint = reproject_adaptive((data_in, header_in), header_in)
     else:
         data_out, footprint = reproject_interp((data_in, header_in), header_in,
                                                order=projection_type)