diff --git a/aiapy/psf/deconvolve.py b/aiapy/psf/deconvolve.py index a81dd2c..869d34c 100644 --- a/aiapy/psf/deconvolve.py +++ b/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: diff --git a/aiapy/psf/tests/test_deconvolve.py b/aiapy/psf/tests/test_deconvolve.py index 4cd02ef..435476c 100644 --- a/aiapy/psf/tests/test_deconvolve.py +++ b/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, + ) diff --git a/changelog/107.feature.rst b/changelog/107.feature.rst new file mode 100644 index 0000000..46ed480 --- /dev/null +++ b/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.