Deconvolving images with the instrument Point Spread Function (PSF)#
This example demonstrates how to deconvolve an AIA image with the instrument point spread function (PSF).
AIA images are subject to convolution with the instrument point-spread function (PSF) due to effects introduced by the filter mesh of the telescope and the CCD, among others. This has the effect of “blurring” the image. The PSF diffraction pattern may also be particularly noticable during the impulsive phase of a flare where the intensity enhancement is very localized. To remove these artifacts, the PSF must be deconvolved from the image.
First, we’ll use a single level 1 image from the 171 Å channel from
15 March 2019. Note that deconvolution should be performed on level 1 images
only. This is because, as with the level 1 data, the PSF model is defined
on the CCD grid. Once deconvolved, the image can be passed to
(see the Registering and aligning level 1 data example).
Next, we’ll calculate the PSF using
aiapy.psf.psf for the 171 Å channel.
The PSF model accounts for several different effects, including diffraction
from the mesh grating of the filters, charge spreading, and jitter. See
Grigis et al (2012)
for more details. Currently, this only works for
\(4096\times4096\) full frame images.
Note that this will be significantly faster if you have a GPU and the
psf = aiapy.psf.psf(m.wavelength)
We’ll plot just a 500-by-500 pixel section centered on the center pixel. The diffraction “arms” extending from the center pixel can often be seen in flare observations due to the intense, small-scale brightening.
Now that we’ve downloaded our image and computed the PSF, we can deconvolve the image with the PSF using the Richardson-Lucy deconvolution algorithm. Note that passing in the PSF is optional. If you exclude it, it will be calculated automatically. However, when deconvolving many images of the same wavelength, it is most efficient to only calculate the PSF once.
m_deconvolved = aiapy.psf.deconvolve(m, psf=psf)
Let’s compare the convolved and deconvolved images.
norm = ImageNormalize(vmin=0, vmax=1.5e4, stretch=AsinhStretch(0.01)) fig = plt.figure() ax = fig.add_subplot(121, projection=m) m.plot(axes=ax, norm=norm) ax = fig.add_subplot(122, projection=m_deconvolved) m_deconvolved.plot(axes=ax, annotate=False, norm=norm) ax.coords.set_axislabel(" ") ax.coords.set_axislabel(" ") ax.coords.set_ticklabel_visible(visible=False) plt.show()
The differences become a bit more obvious when we zoom in. Note that the deconvolution has the effect of “deblurring” the image.
left_corner = 500 * u.arcsec, -600 * u.arcsec right_corner = 1000 * u.arcsec, -100 * u.arcsec fig = plt.figure() m_sub = m.submap( bottom_left=SkyCoord(*left_corner, frame=m.coordinate_frame), top_right=SkyCoord(*right_corner, frame=m.coordinate_frame), ) ax = fig.add_subplot(121, projection=m_sub) m_sub.plot(axes=ax, norm=norm) m_deconvolved_sub = m_deconvolved.submap( bottom_left=SkyCoord(*left_corner, frame=m_deconvolved.coordinate_frame), top_right=SkyCoord(*right_corner, frame=m_deconvolved.coordinate_frame), ) ax = fig.add_subplot(122, projection=m_deconvolved_sub) m_deconvolved_sub.plot(axes=ax, annotate=False, norm=norm) ax.coords.set_axislabel(" ") ax.coords.set_axislabel(" ") ax.coords.set_ticklabel_visible(visible=False) plt.show()
Total running time of the script: (0 minutes 0.000 seconds)