Merge pull request #74 from astropy/feature/convolution-normalization

Updated the normalization on the convolve fine structure method to be more of a matched filter.

CHANGES.rst

@@ -2,12 +2,15 @@
 ----------------------
 
 - Drop support for Python 3.6.
+- Change to normalization for convolution fine structure method to instead use a matched filter.
 
 1.1.0 (2021-11-19)
 ------------------
 
 - Added the option to add a variance array
 - Added the ability to subtract a background array rather than a single value.
+- To accommodate these changes, we now return the cleaned array in the same units as the user provides, ADU rather than
+  electrons and with the background included.
 
 1.0.5 (2016-08-16)
 ------------------

astroscrappy/astroscrappy.pyx

@@ -126,7 +126,7 @@ def detect_cosmics(indat, inmask=None, inbkg=None, invar=None, float sigclip=4.5
         Method to build the fine structure image:\n
         'median': Use the median filter in the standard LA Cosmic algorithm
         'convolve': Convolve the image with the psf kernel to calculate the
-        fine structure image.
+        fine structure image using a matched filter technique.
         Default: 'median'.
 
     psfmodel : {'gauss', 'gaussx', 'gaussy', 'moffat'}, optional
@@ -259,6 +259,8 @@ def detect_cosmics(indat, inmask=None, inbkg=None, invar=None, float sigclip=4.5
             psfk = moffatkernel(psffwhm, psfbeta, psfsize)
         else:
             raise ValueError('Please choose a supported PSF model.')
+        # Reverse the psf kernel as that is what we use in the convolution throughout
+        psfk[:, :] = psfk[::-1, ::-1]
 
     # Define a cosmic ray mask
     # This is what will be returned at the end
@@ -321,9 +323,19 @@ def detect_cosmics(indat, inmask=None, inbkg=None, invar=None, float sigclip=4.5
         sp = s - sp
         del s
 
+        # Find the candidate cosmic rays
+        goodpix = np.logical_not(mask)
+        cosmics = np.logical_and(sp > sigclip, goodpix)
+
         # Build the fine structure image :
+        # TODO: we can gain a reasonable performance boost if we only calculate the fine structure stuff for cosmic
+        # pixels instead of the whole image
         if fsmode == 'convolve':
-            f = convolve(cleanarr, psfk)
+            # Use a match filter, the (reversed) psf convolved with the image / variance
+            f = convolve(cleanarr / noise / noise, psfk)
+            # Normalize by the sum of the weights to get the flux of the source assuming it was actually a psf.
+            # This will pull the value of pixels that do not look like psfs (e.g. CRs) down.
+            f = f / convolve(1.0 / noise / noise, psfk)
         elif fsmode == 'median':
             if sepmed:
                 f = sepmedfilt5(cleanarr)
@@ -347,9 +359,6 @@ def detect_cosmics(indat, inmask=None, inbkg=None, invar=None, float sigclip=4.5
         del m7
         del noise
 
-        # Find the candidate cosmic rays
-        goodpix = np.logical_not(mask)
-        cosmics = np.logical_and(sp > sigclip, goodpix)
         # S' / F is still not exactly what is in the paper as we have subtracted the "sampling flux" from S
         # via the median filter. This should be an optimization because we have removed the smooth component
         # and are therefore comparing the candidate "CR flux" only.

astroscrappy/tests/fake_data.py

@@ -37,12 +37,14 @@ def make_fake_data():
     # Add sky and sky noise
     imdata += 200
 
+    psf_sigma = 3.5
+
     # Add some fake sources
     for i in range(100):
         x = np.random.uniform(low=0.0, high=1001)
         y = np.random.uniform(low=0.0, high=1001)
         brightness = np.random.uniform(low=1000., high=30000.)
-        imdata += gaussian(imdata.shape, x, y, brightness, 3.5)
+        imdata += gaussian(imdata.shape, x, y, brightness, psf_sigma)
 
     # Add the poisson noise
     imdata = np.float32(np.random.poisson(imdata))

astroscrappy/tests/test_astroscrappy.py

@@ -2,6 +2,7 @@
 
 from ..astroscrappy import detect_cosmics
 from . import fake_data
+import numpy as np
 
 
 def test_main():
@@ -12,3 +13,12 @@ def test_main():
     mask, _clean = detect_cosmics(imdata, readnoise=10., gain=1.0,
                                   sigclip=6, sigfrac=1.0)
     assert (mask == expected_crmask).sum() == (1001 * 1001)
+
+
+def test_with_convolve_fine_structure():
+    imdata, expected_crmask = fake_data.make_fake_data()
+    # Convert from sigma to fwhm. Sigma is taken from the fake data utility.
+    psf_fwhm = 3.5 * 2.0 * np.sqrt(2.0 * np.log(2.0))
+    mask, _clean = detect_cosmics(imdata, readnoise=10., gain=1.0,
+                                  sigclip=6, sigfrac=1.0, fsmode='convolve', psffwhm=psf_fwhm)
+    assert (mask == expected_crmask).sum() == (1001 * 1001)

tox.ini

@@ -68,8 +68,8 @@ extras =
 
 commands =
     pip freeze
-    !cov: pytest --pyargs astroscrappy {toxinidir}/docs {posargs}
-    cov: pytest --pyargs astroscrappy {toxinidir}/docs --cov astroscrappy --cov-config={toxinidir}/setup.cfg {posargs}
+    !cov: pytest --pyargs astroscrappy '{toxinidir}/docs' {posargs}
+    cov: pytest --pyargs astroscrappy '{toxinidir}/docs' --cov astroscrappy --cov-config='{toxinidir}/setup.cfg' {posargs}
 
 [testenv:build_docs]
 changedir = docs