Add option to disable roundtrip-checking of transforms

docs/celestial.rst

@@ -102,6 +102,19 @@ also be specified, and be given the shape of the output image using the Numpy
 :class:`~astropy.io.fits.Header`, does not contain information about image
 size).
 
+For the interpolation and adaptive algorithms, an optional third argument,
+``roundtrip_coords`` is accepted. By default, after coordinates are transformed
+from the output plane to the input plane, the input-plane coordinates are
+transformed back to the output plane to ensure that the transformation is
+defined in both directions. This doubles the amount of
+coordinate-transformation work to be done. In speed-critical situations, where
+it is known that the coordinate transformation is defined everywhere, this
+extra work can be disabled by setting ``roundtrip_coords=False``. (Note that
+this is not a question of whether each output pixel maps to an existing *pixel*
+in the input image and vice-versa, but whether it maps to a valid *coordinate*
+in the coordinate system of the input image---regardless of whether that
+coordinate falls within the bounds of the input image.)
+
 As an example, we start off by opening a FITS file using Astropy::
 
     >>> from astropy.io import fits

reproject/adaptive/core.py

@@ -5,27 +5,35 @@
 import numpy as np
 
 from .deforest import map_coordinates
-from ..wcs_utils import efficient_pixel_to_pixel_with_roundtrip
+from ..wcs_utils import (efficient_pixel_to_pixel_with_roundtrip,
+        efficient_pixel_to_pixel)
 
 
 __all__ = ['_reproject_adaptive_2d']
 
 
 class CoordinateTransformer:
 
-    def __init__(self, wcs_in, wcs_out):
+    def __init__(self, wcs_in, wcs_out, roundtrip_coords):
         self.wcs_in = wcs_in
         self.wcs_out = wcs_out
+        self.roundtrip_coords = roundtrip_coords
 
     def __call__(self, pixel_out):
         pixel_out = pixel_out[:, :, 0], pixel_out[:, :, 1]
-        pixel_in = efficient_pixel_to_pixel_with_roundtrip(self.wcs_out, self.wcs_in, *pixel_out)
+        if self.roundtrip_coords:
+            pixel_in = efficient_pixel_to_pixel_with_roundtrip(
+                    self.wcs_out, self.wcs_in, *pixel_out)
+        else:
+            pixel_in = efficient_pixel_to_pixel(
+                    self.wcs_out, self.wcs_in, *pixel_out)
         pixel_in = np.array(pixel_in).transpose().swapaxes(0, 1)
         return pixel_in
 
 
 def _reproject_adaptive_2d(array, wcs_in, wcs_out, shape_out, order=1,
-                           return_footprint=True, center_jacobian=False):
+                           return_footprint=True, center_jacobian=False,
+                           roundtrip_coords=True):
     """
     Reproject celestial slices from an n-d array from one WCS to another
     using the DeForest (2004) algorithm [1]_, and assuming all other dimensions
@@ -47,6 +55,9 @@ def _reproject_adaptive_2d(array, wcs_in, wcs_out, shape_out, order=1,
         Whether to return the footprint in addition to the output array.
     center_jacobian : bool
         Whether to compute centered Jacobians
+    roundtrip_coords : bool
+        Whether to veryfiy that coordinate transformations are defined in both
+        directions.
 
     Returns
     -------
@@ -83,7 +94,7 @@ def _reproject_adaptive_2d(array, wcs_in, wcs_out, shape_out, order=1,
     # Create output array
     array_out = np.zeros(shape_out)
 
-    transformer = CoordinateTransformer(wcs_in, wcs_out)
+    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)
 

reproject/adaptive/high_level.py

@@ -12,7 +12,7 @@
 
 def reproject_adaptive(input_data, output_projection, shape_out=None, hdu_in=0,
                        order='bilinear', return_footprint=True,
-                       center_jacobian=False):
+                       center_jacobian=False, roundtrip_coords=True):
     """
     Reproject celestial slices from an 2d array from one WCS to another using
     the DeForest (2004) adaptive resampling algorithm.
@@ -69,6 +69,9 @@ def reproject_adaptive(input_data, output_projection, shape_out=None, hdu_in=0,
         efficient, and the loss of accuracy is extremely small for
         transformations that vary smoothly between pixels. Defaults to
         ``False``.
+    roundtrip_coords : bool
+        Whether to verify that coordinate transformations are defined in both
+        directions.
 
     Returns
     -------
@@ -90,4 +93,5 @@ def reproject_adaptive(input_data, output_projection, shape_out=None, hdu_in=0,
 
     return _reproject_adaptive_2d(array_in, wcs_in, wcs_out, shape_out,
                                   order=order, return_footprint=return_footprint,
-                                  center_jacobian=center_jacobian)
+                                  center_jacobian=center_jacobian,
+                                  roundtrip_coords=roundtrip_coords)

reproject/adaptive/tests/test_core.py

@@ -22,7 +22,8 @@ def as_high_level_wcs(wcs):
 @pytest.mark.array_compare(single_reference=True)
 @pytest.mark.parametrize('wcsapi', (False, True))
 @pytest.mark.parametrize('center_jacobian', (False, True))
-def test_reproject_adaptive_2d(wcsapi, center_jacobian):
+@pytest.mark.parametrize('roundtrip_coords', (False, True))
+def test_reproject_adaptive_2d(wcsapi, center_jacobian, roundtrip_coords):
 
     # Set up initial array with pattern
     data_in = np.zeros((256, 256))
@@ -48,7 +49,8 @@ def test_reproject_adaptive_2d(wcsapi, center_jacobian):
 
     array_out, footprint_out = reproject_adaptive(
             (data_in, wcs_in), wcs_out, shape_out=(60, 60),
-            center_jacobian=center_jacobian)
+            center_jacobian=center_jacobian,
+            roundtrip_coords=roundtrip_coords)
 
     # Check that surface brightness is conserved in the unrotated case
     assert_allclose(np.nansum(data_in), np.nansum(array_out) * (256 / 60) ** 2, rtol=0.1)
@@ -64,7 +66,8 @@ def test_reproject_adaptive_2d(wcsapi, center_jacobian):
 
 @pytest.mark.array_compare(single_reference=True)
 @pytest.mark.parametrize('center_jacobian', (False, True))
-def test_reproject_adaptive_2d_rotated(center_jacobian):
+@pytest.mark.parametrize('roundtrip_coords', (False, True))
+def test_reproject_adaptive_2d_rotated(center_jacobian, roundtrip_coords):
 
     # Set up initial array with pattern
     data_in = np.zeros((256, 256))
@@ -87,7 +90,8 @@ def test_reproject_adaptive_2d_rotated(center_jacobian):
 
     array_out, footprint_out = reproject_adaptive(
             (data_in, wcs_in), wcs_out, shape_out=(60, 60),
-            center_jacobian=center_jacobian)
+            center_jacobian=center_jacobian,
+            roundtrip_coords=roundtrip_coords)
 
     # ASTROPY_LT_40: astropy v4.0 introduced new default header keywords,
     # once we support only astropy 4.0 and later we can update the reference
@@ -98,7 +102,8 @@ def test_reproject_adaptive_2d_rotated(center_jacobian):
     return array_footprint_to_hdulist(array_out, footprint_out, header_out)
 
 
-def test_reproject_adaptive_high_aliasing_potential():
+@pytest.mark.parametrize('roundtrip_coords', (False, True))
+def test_reproject_adaptive_high_aliasing_potential(roundtrip_coords):
     # Generate sample data with vertical stripes alternating with every column
     data_in = np.arange(40*40).reshape((40, 40))
     data_in = (data_in) % 2
@@ -118,7 +123,8 @@ def test_reproject_adaptive_high_aliasing_potential():
 
     array_out = reproject_adaptive((data_in, wcs_in),
                                    wcs_out, shape_out=(11, 6),
-                                   return_footprint=False)
+                                   return_footprint=False,
+                                   roundtrip_coords=roundtrip_coords)
 
     # The CDELT1 value in wcs_out produces a down-sampling by a factor of two
     # along the output x axis. With the input image containing vertical lines
@@ -138,7 +144,8 @@ def test_reproject_adaptive_high_aliasing_potential():
                       [np.sin(angle), np.cos(angle)]]
     array_out = reproject_adaptive((data_in, wcs_in),
                                    wcs_out, shape_out=(11, 6),
-                                   return_footprint=False)
+                                   return_footprint=False,
+                                   roundtrip_coords=roundtrip_coords)
     np.testing.assert_allclose(array_out, 0.5)
 
     # But if we add a 90-degree rotation to the input coordinates, then when
@@ -151,7 +158,8 @@ def test_reproject_adaptive_high_aliasing_potential():
                      [np.sin(angle), np.cos(angle)]]
     array_out = reproject_adaptive((data_in, wcs_in),
                                    wcs_out, shape_out=(11, 6),
-                                   return_footprint=False)
+                                   return_footprint=False,
+                                   roundtrip_coords=roundtrip_coords)
 
     # Generate the expected pattern of alternating stripes
     data_ref = np.arange(array_out.shape[1]) % 2
@@ -162,7 +170,8 @@ def test_reproject_adaptive_high_aliasing_potential():
     wcs_out.wcs.pc = [[1, 0], [0, 1]]
     array_out = reproject_adaptive((data_in, wcs_in),
                                    wcs_out, shape_out=(11, 6),
-                                   return_footprint=False)
+                                   return_footprint=False,
+                                   roundtrip_coords=roundtrip_coords)
     data_ref = np.arange(array_out.shape[0]) % 2
     data_ref = np.vstack([data_ref] * array_out.shape[1]).T
     np.testing.assert_allclose(array_out, data_ref)

reproject/interpolation/core.py

@@ -4,7 +4,8 @@
 from astropy.wcs import WCS
 
 from ..array_utils import map_coordinates
-from ..wcs_utils import efficient_pixel_to_pixel_with_roundtrip, has_celestial
+from ..wcs_utils import (efficient_pixel_to_pixel_with_roundtrip,
+        efficient_pixel_to_pixel, has_celestial)
 
 
 def _validate_wcs(wcs_in, wcs_out, shape_out):
@@ -56,7 +57,7 @@ def _validate_array_out(array_out, array, shape_out):
 
 
 def _reproject_full(array, wcs_in, wcs_out, shape_out, order=1, array_out=None,
-                    return_footprint=True):
+                    return_footprint=True, roundtrip_coords=True):
     """
     Reproject n-dimensional data to a new projection using interpolation.
 
@@ -80,7 +81,12 @@ def _reproject_full(array, wcs_in, wcs_out, shape_out, order=1, array_out=None,
                             indexing='ij', sparse=False, copy=False)
     pixel_out = [p.ravel() for p in pixel_out]
     # For each pixel in the ouput array, get the pixel value in the input WCS
-    pixel_in = efficient_pixel_to_pixel_with_roundtrip(wcs_out, wcs_in, *pixel_out[::-1])[::-1]
+    if roundtrip_coords:
+        pixel_in = efficient_pixel_to_pixel_with_roundtrip(
+                wcs_out, wcs_in, *pixel_out[::-1])[::-1]
+    else:
+        pixel_in = efficient_pixel_to_pixel(
+                wcs_out, wcs_in, *pixel_out[::-1])[::-1]
     pixel_in = np.array(pixel_in)
 
     if array_out is not None:

reproject/interpolation/high_level.py

@@ -17,7 +17,8 @@
 @deprecated_renamed_argument('independent_celestial_slices', None, since='0.6')
 def reproject_interp(input_data, output_projection, shape_out=None, hdu_in=0,
                      order='bilinear', independent_celestial_slices=False,
-                     output_array=None, return_footprint=True):
+                     output_array=None, return_footprint=True,
+                     roundtrip_coords=True):
     """
     Reproject data to a new projection using interpolation (this is typically
     the fastest way to reproject an image).
@@ -68,6 +69,9 @@ def reproject_interp(input_data, output_projection, shape_out=None, hdu_in=0,
         extremely large files.
     return_footprint : bool
         Whether to return the footprint in addition to the output array.
+    roundtrip_coords : bool
+        Whether to verify that coordinate transformations are defined in both
+        directions.
 
     Returns
     -------
@@ -86,5 +90,6 @@ def reproject_interp(input_data, output_projection, shape_out=None, hdu_in=0,
     if isinstance(order, str):
         order = ORDER[order]
 
-    return _reproject_full(array_in, wcs_in, wcs_out, shape_out=shape_out, order=order,
-                           array_out=output_array, return_footprint=return_footprint)
+    return _reproject_full(array_in, wcs_in, wcs_out, shape_out=shape_out,
+            order=order, array_out=output_array,
+            return_footprint=return_footprint, roundtrip_coords=roundtrip_coords)