Merge branch 'psf-negative-pixels' into 'master'

Remove negative intensities before deconvolution

Closes #95

See merge request LMSAL_HUB/aia_hub/aiapy!107

aiapy/psf/deconvolve.py

@@ -2,6 +2,7 @@
 Deconvolve an AIA image with the channel point spread function
 """
 import copy
+import warnings
 
 import numpy as np
 try:
@@ -11,11 +12,12 @@
     HAS_CUPY = False
 
 from .psf import psf as calculate_psf
+from aiapy.util import AiapyUserWarning
 
 __all__ = ['deconvolve']
 
 
-def deconvolve(smap, psf=None, iterations=25, use_gpu=True):
+def deconvolve(smap, psf=None, iterations=25, clip_negative=True, use_gpu=True):
     """
     Deconvolve an AIA image with the point spread function
 
@@ -38,6 +40,8 @@ def deconvolve(smap, psf=None, iterations=25, use_gpu=True):
         The point spread function. If None, it will be calculated
     iterations : `int`, optional
         Number of iterations in the Richardson-Lucy algorithm
+    clip_negative : `bool`, optional
+        If the image has negative intensity values, set them to zero.
     use_gpu : `bool`, optional
         If True and `~cupy` is installed, do PSF deconvolution on the GPU
         with `~cupy`.
@@ -58,6 +62,12 @@ def deconvolve(smap, psf=None, iterations=25, use_gpu=True):
     """
     # TODO: do we need a check to make sure this is a full-frame image?
     img = smap.data
+    if clip_negative:
+        img = np.where(img < 0, 0, img)
+    if np.any(img < 0):
+        warnings.warn(
+            'Image contains negative intensity values. Consider setting clip_negative to True',
+            AiapyUserWarning)
     if psf is None:
         psf = calculate_psf(smap.wavelength)
     if HAS_CUPY and use_gpu:

aiapy/psf/tests/test_deconvolve.py

@@ -2,10 +2,12 @@
 Test deconvolution
 """
 import pytest
+import numpy as np
 import sunpy.data.test
 import sunpy.map
 
 import aiapy.psf
+from aiapy.util import AiapyUserWarning
 
 
 def test_deconvolve(aia_171_map):
@@ -24,3 +26,17 @@ def test_deconvolve_specify_psf(aia_171_map, psf):
     map_decon = aiapy.psf.deconvolve(aia_171_map, psf=psf, iterations=1)
     assert isinstance(map_decon, sunpy.map.GenericMap)
     assert map_decon.data.shape == aia_171_map.data.shape
+
+
+def test_deconvolve_negative_pixels(aia_171_map, psf):
+    aia_171_map_neg = aia_171_map._new_instance(
+        np.where(aia_171_map.data < 1, -1, aia_171_map.data),
+        aia_171_map.meta,
+    )
+    with pytest.warns(AiapyUserWarning):
+        _ = aiapy.psf.deconvolve(
+            aia_171_map_neg,
+            psf=psf,
+            iterations=1,
+            clip_negative=False,
+        )

changelog/107.feature.rst

@@ -0,0 +1,2 @@
+Add a flag to :func:`aiapy.psf.deconvolve` that sets negative intensity values
+to zero before performing the deconvolution.