Efficient deconvolution methods for astronomical imaging: algorithms and IDL-GPU codes

Abstract

Context. The Richardson-Lucy method is the most popular deconvolution method in astronomy because it preserves the number of counts and the non-negativity of the original object. Regularization is, in general, obtained by an early stopping of Richardson-Lucy iterations. In the case of point-wise objects such as binaries or open star clusters, iterations can be pushed to convergence. However, it is well-known that Richardson-Lucy is an inefficient method. In most cases and, in particular, for low noise levels, acceptable solutions are obtained at the cost of hundreds or thousands of iterations, thus several approaches to accelerating Richardson-Lucy have been proposed. They are mainly based on Richardson-Lucy being a scaled gradient method for the minimization of the Kullback-Leibler divergence, or Csiszár I-divergence, which represents the data-fidelity function in the case of Poisson noise. In this framework, a line search along the descent direction is considered for reducing the number of iterations. Aims. A general optimization method, referred to as the scaled gradient projection method, has been proposed for the constrained minimization of continuously differentiable convex functions. It is applicable to the non-negative minimization of the Kullback-Leibler divergence. If the scaling suggested by Richardson-Lucy is used in this method, then it provides a considerable increase in the efficiency of Richardson-Lucy. Therefore the aim of this paper is to apply the scaled gradient projection method to a number of imaging problems in astronomy such as single image deconvolution, multiple image deconvolution, and boundary effect correction. Methods. Deconvolution methods are proposed by applying the scaled gradient projection method to the minimization of the Kullback-Leibler divergence for the imaging problems mentioned above and the corresponding algorithms are derived and implemented in interactive data language. For all the algorithms, several stopping rules are introduced, including one based on a recently proposed discrepancy principle for Poisson data. To attempt to achieve a further increase in efficiency, we also consider an implementation on graphic processing units. Results. The proposed algorithms are tested on simulated images. The acceleration of scaled gradient projection methods achieved with respect to the corresponding Richardson-Lucy methods strongly depends on both the problem and the specific object to be reconstructed, and in our simulations the improvement achieved ranges from about a factor of 4 to more than 30. Moreover, significant accelerations of up to two orders of magnitude have been observed between the serial and parallel implementations of the algorithms. The codes are available upon request.

On the constrained minimization of smooth Kurdyka—Łojasiewicz functions with the scaled gradient projection method

“…The SGP method is a variable metric forward-backward algorithm [11,12] which has been exploited in the last years for the solution of different real-world inverse problems [5,6,7,16,17,18]. The main difference between SGP and the standard forward-backward schemes is the presence of two independent parameters α k and λ k with a complete different role: while the last one is automatically computed with the Armijo condition (2) to guarantee the sufficient decrease of the objective function, the first one can be chosen to improve the actual convergence rate of the method, exploiting thirty years of literature in numerical optimization [3,13,19,20].…”

Limited-memory scaled gradient projection methods for real-time image deconvolution in microscopy

“…To better explain how SGP improves the RL algorithm, we shortly recall its iteration for problem (4). The reader is referred to [5] for the convergence analysis of the general scheme and to [27,28,29,30,31,32] for examples of SGP applications in different imaging problems. Given an initial feasible x (0) ,…”

Non-parametric PSF estimation from celestial transit solar images using blind deconvolution

“…Therefore, it is important to properly handle the frontiers of our images and avoid both implicit FFT frontier periodizations and the influence of unobserved image values integrated by the filter extension. Inspired by the work of of Almeida & Figueiredo (2013a), we implement their unknown boundary conditions, where instead of expanding the observation using zero padding as in Anconelli et al (2006) and Prato et al (2012), the boundaries are considered unknown in the convolution process and then, by properly selecting the pixels inside the known boundaries, the observed image is obtained. In a nutshell, the unknown boundary conditions method proceeds by considering every convolution in a bigger space obtained by expanding the original one by the radius of the filter in all directions, i.e., by a border of width b.…”

“…Deconvolution is a ubiquitous data processing method that arises in a variety of applications, e.g., biomedical imaging (Jefferies et al 2002), astronomy (Prato et al 2012;Shearer 2013), remote sensing (Dong et al 2011), and video and photo enhancement (Fergus et al 2006). Formally, this problem amounts to recovering a signal x from blurred and corrupted observations y:…”

“…If the observations are corrupted by an additive Gaussian noise, the deconvolution problem can be solved by a regularized least-squares (LS) minimization (Almeida & Almeida 2010). In the presence of Poisson noise, the problem can be formulated using the Kullback-Leibler (KL) divergence as the data fidelity term (Fish et al 1995;Prato et al 2012Prato et al , 2013, which represents a more complex function to minimize. To avoid dealing with the KL divergence, the Poisson corrupted data can also be handled through a Variance Stabilization Transform (VST) (Dupé et al 2009;Shearer et al 2012;Shearer 2013).…”