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.
@@ -68,6 +68,9 @@ def reproject_adaptive(input_data, output_projection, shape_out=None, hdu_in=0,
         averaged to produce values at the pixel grid points. This is more
         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
     -------
@@ -89,4 +92,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)