STARRED: a two-channel deconvolution method with Starlet regularization

The spatial resolution of astronomical images is limited by atmospheric turbulence and diffraction in the telescope optics, resulting in blurred images. This makes it difficult to accurately measure the brightness of blended objects because the contributions from adjacent objects are mixed in a time-variable manner due to changes in the atmospheric conditions. However, this effect can be corrected by characterizing the Point Spread Function (PSF), which describes how a point source is blurred on a detector. This function can be estimated from the stars in the field of view, which provides a natural sampling of the PSF across the entire field of view. Once the PSF is estimated, it can be removed from the data through the so-called deconvolution process, leading to images of improved spatial resolution. The deconvolution operation is an ill-posed inverse problem due to noise and pixelization of the data. To solve this problem, regularization is necessary to guarantee the robustness of the solution. Regularization can take the form of a sparse prior, meaning that the recovered solution can be represented with only a few basis eigenvectors. STARRED is a Python package developed in the context of the COSMOGRAIL collaboration and applies to a vast variety of astronomical problems. It proposes to use an isotropic wavelet basis, called Starlets, to regularize the solution of the deconvolution problem. This family of wavelets has been shown to be well-suited to represent astronomical objects. STARRED provides two modules to first reconstruct the PSF, and then perform the deconvolution. It is based on two key concepts: i) the image is reconstructed in two separate channels, one for the point sources and one for the extended sources, and ii) the code relies on the deliberate choice of not completely removing the effect of the PSF, but rather bringing the image to a higher resolution.


Summary
The spatial resolution of astronomical images is limited by atmospheric turbulence and diffraction in the telescope optics, resulting in blurred images.This makes it difficult to accurately measure the brightness of blended objects because the contributions from adjacent objects are mixed in a time-variable manner due to changes in the atmospheric conditions.However, this effect can be corrected by characterizing the Point Spread Function (PSF), which describes how a point source is blurred on a detector.This function can be estimated from the stars in the field of view, which provides a natural sampling of the PSF across the entire field of view.
Once the PSF is estimated, it can be removed from the data through the so-called deconvolution process, leading to images of improved spatial resolution.The deconvolution operation is an ill-posed inverse problem due to noise and pixelization of the data.To solve this problem, regularization is necessary to guarantee the robustness of the solution.Regularization can take the form of a sparse prior, meaning that the recovered solution can be represented with only a few basis eigenvectors.
STARRED is a Python package developed in the context of the COSMOGRAIL collaboration and applies to a vast variety of astronomical problems.It proposes to use an isotropic wavelet basis, called Starlets (Starck et al., 2015), to regularize the solution of the deconvolution problem.This family of wavelets has been shown to be well-suited to represent astronomical objects.STARRED provides two modules to first reconstruct the PSF, and then perform the deconvolution.It is based on two key concepts: i) the image is reconstructed in two separate channels, one for the point sources and one for the extended sources, and ii) the code relies on the deliberate choice of not completely removing the effect of the PSF, but rather bringing the image to a higher resolution with a known Gaussian PSF.This last point allows us to suppress the deconvolution artifacts, which occur when attempting to deconvolve to an infinite resolution, as most of other techniques do.Finally, STARRED uses JAX automatic differentiation to ensure gradient-informed optimization of this high-dimension optimization problem (Bradbury et al., 2018).

Statement of need
Image deconvolution is a widespread problem, not only in astronomy, but also in other scientific fields such as medical imaging and microscopy.In astronomy, STARRED can be applied to a wide range of cases, such as extracting light curves from lensed quasars, extracting light curves of Cepeheids in crowded fields, deblending supernovae and their host galaxies, among many others.The software is not limited to single images, but can also process time series data.The most popular early image deconvolution methods include Högbom (1974) and Lucy & Walsh (2003).The latter, the Richardson-Lucy algorithm, was updated to decompose images into point source and extended source channels (Becker et al., 2003).An extension to more than two channels was proposed in Bontekoe et al. (1994).These works have been further developed and enhanced under Bayesian frameworks (Selig & Enßlin, 2015) and even using neural networks with a 1-channel approach (Akhaury et al., 2022;Sureau et al., 2020).
What all these algorithms have in common is that they intend to fully correct for the PSF, which leads to artifacts due to the poor representation of frequencies aliased outside the frequency domain allowed by the pixel size of the deconvolved image.In other words, should this effect be completely eliminated, then it would be necessary to sample the deconvolved image at infinitely small intervals.This drawback is overcome with MCS (Magain et al., 1998) since a narrower version of the original PSF is generated to satisfy Nyquist-Shannon sampling theorem, and thus its outputs are not prone to artifacts.
STARRED is based on the principles of the Fortran code MCS and Firedec (Cantale et al., 2016).However, STARRED has several important improvements over these two methods: i) it uses a family of isotropic wavelets well-suited to represent astronomical objects, Starlets, ii) it is auto-differentiable, facilitating the use of optimizers based on gradient descent, iii) it is faster than previous codes (220 CPU seconds versus 260 with MCS for the case of Figure 1) and GPU scalable, and iv) it can handle time series by processing all data together in a multi-epoch deconvolution.The use of a second channel for the point sources allows us to explicitly account for the highest frequencies, which cannot be sparsely represented with Starlets.
The scalability, speed of the code and its ability to handle time series are especially relevant for processing the large volume of data expected from future large sky time-domain surveys.In addition, STARRED is much more robust to local minima than MCS or Firedec since it uses gradient descent with momentum algorithms and because of a well-justified regularization, i.e., Starlets, where galaxies are known to be sparse.

Functionality
The main functionality of STARRED is built around deconvolution and PSF generation classes.The starred.psf module receives as input several stars of the field of interest and optionally the noise maps of the observations.It is capable of computing the convolution kernel that relates the original image to an image with a known Gaussian PSF with a Full Width at Half Maximum (FWHM) of 2 pixels.We call it the narrow PSF.It is expressed as a sum of an analytic function with circular symmetry and a Starlet regularized grid of pixels.
The starred.deconvolution module features all the necessary classes to deconvolve an image or a time series.It is necessary to provide the observations, the retrieved narrow PSF and, optionally, the noise maps.This step decomposes the deconvolved image into two channels of point sources and extended sources.The latter are represented by a Starlet regularized grid of pixels.During this process, the image is decomposed into several scales, each of which captures particular frequency features.It is then reconstructed from only the highest coefficients of each scale, which contain the signal but not the noise.This procedure, referred as soft thresholding, allows us to remove the noise from both the reconstructed PSF and the deconvolved image.
STARRED is implemented in JAX for a fast computation of all the derivatives of the free parameters with automatic differentiation.It relies on several gradient descent algorithms from the optax and scipy Python packages.
When processing time series, this module allows us to perform a joint deconvolution of all the images, benefiting from the combined signal-to-noise ratio of all frames but allowing to measure the photometry of the point sources in each individual image.This module can therefore extract the light curves of blended point sources even superposed on a background of any arbitrary shape.
Finally, the starred.plotsmodule can be used to generate web-friendly images with the appropriate cuts.
retain copyright and release the work under a Creative Commons Attribution 4.0 International License (CC BY 4.0).
This is particularly useful in the context of large time-domain sky surveys like Rubin-LSST(Vera C. Rubin Observatory LSST Solar System Science Collaboration et al., 2021).STARRED can jointly deconvolve multiple observation epochs based on their respective estimated PSFs, thereby removing epoch-to-epoch PSF variations before obtaining photometric measurements.An example of such multi-epoch deconvolution is shown in Figure1.Here, the position of the point sources and the extended source channel are constrained across all epochs, whereas the amplitude of the point sources can vary thoughout the time series.

Figure 1 :
Figure 1: Multi-epoch deconvolution of the DES2038-4008 lensed quasar, SDSS-r filter.From left to right, the panels show: a stacked image of all 31 epochs used in the deconvolution (32x32 pixels); STARRED multi-epoch deconvolution showing only the background (second panel, 64x64 pixels) and also including the point sources (third panel, 64x64 pixels) with a runtime of 17 seconds (GPU: RTX 3060); a deep, high-resolution observation (Shajib et al., 2018) from the Hubble Space Telescope, F160W filter (Credit: TDCOSMO Collaboration, GO-15320, PI: Treu).