Merge branch 'main' into sync_drizzle

.github/workflows/ci.yml

@@ -3,7 +3,7 @@ name: CI
 on:
   push:
     branches:
-      - master
+      - main
       - '*.x'
     tags:
       - "*"

CHANGES.rst

@@ -4,11 +4,28 @@
 Release Notes
 =============
 
-1.15.0 (unreleased)
+
+1.15.1 (Unreleased)
+===================
+
+- Fixed the warning type for the "gaussian", "lanczos2", and "lanczos3" kernels
+  (``DeprecationWarning`` to ``Warning``). [#141]
+
+
+1.15.0 (2024-02-16)
 ===================
 
 - Dropped Python 3.8. [#128]
 
+- Fixed a bug in the pixmap coordinate inversion routine. [#137]
+
+- Deprecated "tophat" kernel which will be remover in the next release. It is
+  not working correctly and should not be used. [#140]
+
+- Added warnings that "gaussian", "lanczos2", and "lanczos3" kernels are not
+  flux conserving. [#140]
+
+
 1.14.4 (2023-11-15)
 ===================
 

README.rst

@@ -5,7 +5,7 @@ drizzle Documentation
     :target: http://www.astropy.org
     :alt: Powered by Astropy Badge
 
-.. image:: https://codecov.io/github/spacetelescope/drizzle/branch/master/graphs/badge.svg
+.. image:: https://codecov.io/github/spacetelescope/drizzle/branch/main/graphs/badge.svg
     :target: https://codecov.io/gh/spacetelescope/drizzle
     :alt: Drizzle's Coverage Status
 

drizzle/dodrizzle.py

@@ -112,8 +112,11 @@ def dodrizzle(insci, input_wcs, inwht,
     kernel: str, optional
         The name of the kernel used to combine the input. The choice of
         kernel controls the distribution of flux over the kernel. The kernel
-        names are: "square", "gaussian", "point", "tophat", "turbo", "lanczos2",
-        and "lanczos3". The square kernel is the default.
+        names are: "square", "turbo", "point", "gaussian", "lanczos2",
+        and "lanczos3".
+
+        .. warning::
+            The "gaussian" and "lanczos2/3" kernels **DO NOT** conserve flux.
 
     fillval: str, optional
         The value a pixel is set to in the output if the input image does

drizzle/drizzle.py

@@ -60,8 +60,11 @@ def __init__(self, infile="", outwcs=None,
         kernel : str, optional
             The name of the kernel used to combine the inputs. The choice of
             kernel controls the distribution of flux over the kernel. The kernel
-            names are: "square", "gaussian", "point", "tophat", "turbo", "lanczos2",
-            and "lanczos3". The square kernel is the default.
+            names are: "square", "turbo", "point", "gaussian", "lanczos2",
+            and "lanczos3".
+
+            .. warning::
+              The "gaussian" and "lanczos2/3" kernels **DO NOT** conserve flux.
 
         fillval : str, otional
             The value a pixel is set to in the output if the input image does

drizzle/tests/test_drizzle.py

@@ -642,3 +642,157 @@ def test_context_planes():
 
     driz.add_image(image, inwcs)
     assert driz.outcon.shape == (2, 10, 10)
+
+
+@pytest.mark.parametrize(
+    'kernel',
+    [
+        'square',
+        'point',
+        'turbo',
+        pytest.param(
+            'lanczos2',
+            marks=pytest.mark.xfail(reason='Not a flux-conserving kernel'),
+        ),
+        pytest.param(
+            'lanczos3',
+            marks=pytest.mark.xfail(reason='Not a flux-conserving kernel'),
+        ),
+        pytest.param(
+            'gaussian',
+            marks=pytest.mark.xfail(reason='Not a flux-conserving kernel'),
+        ),
+    ],
+)
+def test_flux_conservation_nondistorted(kernel):
+    n = 200
+    in_shape = (n, n)
+
+    # input coordinate grid:
+    y, x = np.indices(in_shape, dtype=np.float64)
+
+    # simulate a gaussian "star":
+    fwhm = 2.9
+    x0 = 50.0
+    y0 = 68.0
+    sig = fwhm / (2.0 * np.sqrt(2.0 * np.log(2.0 * fwhm)))
+    sig2 = sig * sig
+    star = np.exp(-0.5 / sig2 * ((x.astype(np.float32) - x0)**2 + (y.astype(np.float32) - y0)**2))
+    in_sci = (star / np.sum(star)).astype(np.float32)  # normalize to 1
+    in_wht = np.ones(in_shape, dtype=np.float32)
+
+    # linear shift:
+    xp = x + 0.5
+    yp = y + 0.2
+
+    pixmap = np.dstack([xp, yp])
+
+    out_shape = (int(yp.max()) + 1, int(xp.max()) + 1)
+    # make sure distorion is not moving flux out of the image towards negative
+    # coordinates (just because of the simple way of how we account for output
+    # image size)
+    assert np.min(xp) > -0.5 and np.min(yp) > -0.5
+
+    out_sci = np.zeros(out_shape, dtype=np.float32)
+    out_ctx = np.zeros(out_shape, dtype=np.int32)
+    out_wht = np.zeros(out_shape, dtype=np.float32)
+
+    cdrizzle.tdriz(
+        in_sci,
+        in_wht,
+        pixmap,
+        out_sci,
+        out_wht,
+        out_ctx,
+        pixfrac=1.0,
+        scale=1.0,
+        kernel=kernel,
+        in_units="cps",
+        expscale=1.0,
+        wtscale=1.0,
+    )
+
+    assert np.allclose(
+        np.sum(out_sci * out_wht),
+        np.sum(in_sci),
+        atol=0.0,
+        rtol=0.0001,
+    )
+
+@pytest.mark.parametrize(
+    'kernel',
+    [
+        'square',
+        'point',
+        'turbo',
+        pytest.param(
+            'lanczos2',
+            marks=pytest.mark.xfail(reason='Not a flux-conserving kernel'),
+        ),
+        pytest.param(
+            'lanczos3',
+            marks=pytest.mark.xfail(reason='Not a flux-conserving kernel'),
+        ),
+        pytest.param(
+            'gaussian',
+            marks=pytest.mark.xfail(reason='Not a flux-conserving kernel'),
+        ),
+    ],
+)
+def test_flux_conservation_distorted(kernel):
+    n = 200
+    in_shape = (n, n)
+
+    # input coordinate grid:
+    y, x = np.indices(in_shape, dtype=np.float64)
+
+    # simulate a gaussian "star":
+    fwhm = 2.9
+    x0 = 50.0
+    y0 = 68.0
+    sig = fwhm / (2.0 * np.sqrt(2.0 * np.log(2.0 * fwhm)))
+    sig2 = sig * sig
+    star = np.exp(-0.5 / sig2 * ((x.astype(np.float32) - x0)**2 + (y.astype(np.float32) - y0)**2))
+    in_sci = (star / np.sum(star)).astype(np.float32)  # normalize to 1
+    in_wht = np.ones(in_shape, dtype=np.float32)
+
+    # linear shift:
+    xp = x + 0.5
+    yp = y + 0.2
+    # add distortion:
+    xp += 1e-4 * x**2 + 1e-5 * x * y
+    yp += 1e-3 * y**2 - 2e-5 * x * y
+
+    pixmap = np.dstack([xp, yp])
+
+    out_shape = (int(yp.max()) + 1, int(xp.max()) + 1)
+    # make sure distorion is not moving (pixels with) flux out of the image
+    # towards negative coordinates (just because of the simple way of how we
+    # account for output image size):
+    assert np.min(xp) > -0.5 and np.min(yp) > -0.5
+
+    out_sci = np.zeros(out_shape, dtype=np.float32)
+    out_ctx = np.zeros(out_shape, dtype=np.int32)
+    out_wht = np.zeros(out_shape, dtype=np.float32)
+
+    cdrizzle.tdriz(
+        in_sci,
+        in_wht,
+        pixmap,
+        out_sci,
+        out_wht,
+        out_ctx,
+        pixfrac=1.0,
+        scale=1.0,
+        kernel=kernel,
+        in_units="cps",
+        expscale=1.0,
+        wtscale=1.0,
+    )
+
+    assert np.allclose(
+        np.sum(out_sci * out_wht),
+        np.sum(in_sci),
+        atol=0.0,
+        rtol=0.0001,
+    )

drizzle/tests/test_overlap_calc.py

@@ -1,5 +1,6 @@
 import pytest
 from math import sqrt
+from itertools import product
 
 import numpy as np
 
@@ -70,6 +71,17 @@ def test_invert_pixmap():
         assert np.allclose(xyin, [xr, yr], atol=0.05)
 
 
+def test_invert_small_pixmap():
+    yin, xin = np.indices((2, 2), dtype=float)
+    pixmap = np.dstack([xin, yin])
+
+    test_coords = list(product(*(2 * [[-0.5, 1.5]])))
+
+    for xr, yr in test_coords:
+        xyin = invert_pixmap(pixmap, [xr, yr], [[-0.5, 1.5], [-0.5, 1.5]])
+        assert np.allclose(xyin, [xr, yr], atol=0.05)
+
+
 def test_poly_intersection_with_self():
     p = [(0, 0), (1, 0), (1, 1), (0, 1)]
 

src/cdrizzleapi.c

@@ -84,6 +84,7 @@ tdriz(PyObject *obj UNUSED_PARAM, PyObject *args, PyObject *keywords)
   struct driz_error_t error;
   struct driz_param_t p;
   integer_t isize[2], psize[2], wsize[2];
+  char warn_msg[96];
 
   driz_log_handle = driz_log_init(driz_log_handle);
   driz_log_message("starting tdriz");
@@ -182,6 +183,22 @@ tdriz(PyObject *obj UNUSED_PARAM, PyObject *args, PyObject *keywords)
     goto _exit;
   }
 
+  if (kernel == kernel_tophat) {
+      if (sprintf(warn_msg,
+            "Kernel '%s' has been deprecated and it will be removed in a future release.",
+            kernel_str) < 1) {
+          strcpy(warn_msg, "Selected kernel has been deprecated and it will be removed in a future release.");
+      }
+      PyErr_WarnEx(PyExc_DeprecationWarning, warn_msg, 1);
+  } else if (kernel == kernel_gaussian || kernel == kernel_lanczos2 || kernel == kernel_lanczos3) {
+      if (sprintf(warn_msg,
+            "Kernel '%s' is not a flux-conserving kernel.",
+            kernel_str) < 1) {
+          strcpy(warn_msg, "Selected kernel is not a flux-conserving kernel.");
+      }
+      PyErr_WarnEx(PyExc_Warning, warn_msg, 1);
+  }
+
   if (pfract <= 0.001){
     printf("kernel reset to POINT due to pfract being set to 0.0...\n");
     kernel_str2enum("point", &kernel, &error);

src/cdrizzlemap.c

@@ -280,8 +280,8 @@ invert_pixmap(struct driz_param_t *par, double xout, double yout, double *xin,
     xmax = ((double)par->xmax) + 0.5;
     ymin = ((double)par->ymin) - 0.5;
     ymax = ((double)par->ymax) + 0.5;
-    dx = xmax;
-    dy = ymax;
+    dx = xmax - xmin;
+    dy = ymax - ymin;
 
     niter = 0;