Support 32-bit input/output

Author: svank

reproject/adaptive/core.py

@@ -117,8 +117,10 @@ def _reproject_adaptive_2d(
     "broadcast" across those images.
     """
 
-    # Make sure image is floating point
-    array_in = np.asarray(array, dtype=float)
+    # Make sure image is floating point, but preserve the floating-point precision if it's floating point already
+    array_in = np.asarray(array)
+    if not np.issubdtype(array.dtype, np.floating):
+        array_in = np.asarray(array, dtype=float)
     shape_out = tuple(shape_out)
 
     # Check dimensionality of WCS and shape_out

reproject/adaptive/deforest.pyx

@@ -33,6 +33,7 @@ import numpy as np
 
 cimport cython
 cimport numpy as np
+from cython cimport floating
 from libc.math cimport atan2, ceil, cos, exp, fabs, floor, round, sin, sqrt
 from libc.stdlib cimport qsort
 
@@ -155,7 +156,7 @@ cdef double clip(double x, double vmin, double vmax, int cyclic,
 @cython.wraparound(False)
 @cython.nonecheck(False)
 @cython.cdivision(True)
-cdef bint sample_array(const double[:,:,:] source, double[:] dest,
+cdef bint sample_array(const floating[:,:,:] source, floating[:] dest,
         double x, double y, int x_cyclic, int y_cyclic,
         bint out_of_range_nearest) noexcept nogil:
     x = clip(x, 0, source.shape[2] - 1, x_cyclic, out_of_range_nearest)
@@ -352,7 +353,7 @@ BAD_VALUE_MODES['ignore'] = 3
 @cython.wraparound(False)
 @cython.nonecheck(False)
 @cython.cdivision(True)
-def map_coordinates(const double[:,:,:] source, double[:,:,:] target, Ci, int max_samples_width=-1,
+def map_coordinates(const floating[:,:,:] source, floating[:,:,:] 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,
                     bint center_jacobian=False, bint despiked_jacobian=False,
@@ -492,7 +493,7 @@ def map_coordinates(const double[:,:,:] source, double[:,:,:] target, Ci, int ma
     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])
+    cdef floating[:] value = np.empty_like(source, shape=source.shape[0])
     cdef double[:] P1 = np.empty((2,))
     cdef double[:] P2 = np.empty((2,))
     cdef double[:] P3 = np.empty((2,))

reproject/adaptive/tests/test_core.py

@@ -62,6 +62,43 @@ def test_reproject_adaptive_2d(wcsapi, center_jacobian, roundtrip_coords):
     return array_footprint_to_hdulist(array_out, footprint_out, header_out)
 
 
+def test_reproject_adaptive_dtypes():
+    # Set up initial array with pattern
+    data_in = np.zeros((256, 256))
+    data_in[::20, :] = 1
+    data_in[:, ::20] = 1
+    data_in[10::20, 10::20] = 1
+
+    # Define a simple input WCS
+    wcs_in = WCS(naxis=2)
+    wcs_in.wcs.crpix = 128.5, 128.5
+    wcs_in.wcs.cdelt = -0.01, 0.01
+
+    # Define a lower resolution output WCS
+    wcs_out = WCS(naxis=2)
+    wcs_out.wcs.crpix = 30.5, 30.5
+    wcs_out.wcs.cdelt = -0.0427, 0.0427
+
+    header_out = wcs_out.to_header()
+
+    array_out, footprint_out = reproject_adaptive(
+        (data_in, wcs_in),
+        wcs_out,
+        shape_out=(60, 60),
+    )
+
+    array_out_32bit = np.empty_like(array_out, dtype=np.float32)
+    _, footprint_out = reproject_adaptive(
+        (data_in.astype(np.float32), wcs_in),
+        wcs_out,
+        shape_out=(60, 60),
+        output_array=array_out_32bit,
+    )
+
+    # Check that surface brightness is conserved in the unrotated case
+    assert_allclose(array_out, array_out_32bit, rtol=0.1)
+
+
 @pytest.mark.parametrize("axis", ("x", "y"))
 @pytest.mark.parametrize("center_jacobian", (True, False))
 def test_reproject_adaptive_despike_jacobian(axis, center_jacobian):