An alternating minimization method for blind deconvolution from Poisson data

Blind deconvolution is a particularly challenging inverse problem since information on both the desired target and the acquisition system have to be inferred from the measured data. When the collected data are affected by Poisson noise, this problem is typically addressed by the minimization of the Kullback-Leibler divergence, in which the unknowns are sought in particular feasible sets depending on the a priori information provided by the specific application. If these sets are separated, then the resulting constrained minimization problem can be addressed with an inexact alternating strategy. In this paper we apply this optimization tool to the problem of reconstructing astronomical images from adaptive optics systems, and we show that the proposed approach succeeds in providing very good results in the blind deconvolution of nondense stellar clusters.


Blind deconvolution from Poisson data
Several acquisition systems provide the information as the convolution of an unknown distribution with a hardware-dependent kernel which transforms the target input in the measured output. A common instance of deconvolution problem is the image deblurring, in which the target image is distorted by a point spread function (PSF) whose features depend on the acquisition system and the environmental occurrences [1,2]. If the PSF is not completely known, it has to be estimated within the deconvolution process together with the unknown distribution and the reconstruction problem increases its ill-posedness due to the large number of possible solutions, most of which have no physical meaning. This problem is known in the literature as blind deconvolution, and a key requirement for its solution is the availability of suitable a priori information. In this paper we address the blind deconvolution problem of reconstructing a blurred astronomical image acquired by a ground-based telescope when the PSF is unknown. Due to the Poisson nature of the noise introduced in the recording process of the image, a maximum likelihood approach leads to the reformulation of the imaging problem as the constrained minimization of the Kullback-Leilber (or Csiszár I- [3]) divergence where g is the measured image, h the PSF, f the target image, b a known background radiation, S the set of multi-indices i = (i 1 , i 2 ) and h * f denotes the convolution product. As concerns the feasible sets for f and h, we consider non-negativity and flux conservation for the target image, while for the PSF we impose non-negativity, normalization to 1 and an upper bound deduced by the so-called Strehl ratio (SR), which is the ratio of peak diffraction intensity of an aberrated versus perfect waveform and can be estimated from the telescope and image features (see e.g. [4]). As a result of these choices, the mathematical formulation of the minimization problem becomes where c is the background-subtracted flux of the detected image and s is the upper bound derived from the SR.

Alternating minimization
Problem (2) is a hard task for an optimization method since the objective function is convex if restricted to f or h, but is nonconvex with respect to (f , h) [5], thus admitting multiple stationary points. Due to the separability and the simplicity of the constraints, here we adopt an inexact alternating minimization (IAM) strategy [6,7,8], in which at each (outer) iteration k ≥ 1 the two subproblems h ≤ N h (with N f and N h arbitrarily fixed) of (inner) iterations of the scaled gradient projection (SGP) method [9]. The convergence of the resulting scheme has been proved in Theorem 4.2 of [6]: (3) and (4), respectively, is a stationary point for problem (2). Since SGP has been described in several papers concerning deblurring and denoising problems [9,10,11,12,13,14], here we do not address the steps of the algorithm and its details. A complete and recent description can be found in [15], while the routines in MATLAB and Interactive Data Language (IDL) implementing SGP can be freely downloaded by the website http://www.unife.it/prisma/software. The result in Theorem 2.1 is the usual convergence property proved for a general first-order method in absence of further convexity hypotheses of the objective function. We remark that, thanks to the constraints we adopted, the sequence (f (k) , h (k) ) generated by our method is bounded and therefore the existence of limit points is assured, even if there is no practical way to demonstrate if these points are minimizers.

Numerical experiments
In this section we test the effectiveness of the IAM procedure in an astronomical blind deconvolution problem. In particular, we used simulations of a single mirror of the Large Binocular Telescope (LBT) and its First Light Adaptive Optics (FLAO) system [16], which provides SR up to 0.9 in K-band.

Creating the simulations
We chose two PSFs with SR = 0.81 and 0.60, respectively, while for the astronomical objects we restricted our analysis to the 512×512 images of three star clusters, in which the positions of the stars have been computed following a Gaussian distribution around the center of the image and their magnitudes are randomly distributed around m = 16 with a standard deviation of about 0.4. The number of stars in the two images has been set equal to 50, 100 and 200, respectively. The three targets will be denoted in the following with SC50, SC100 and SC200. We convolved the three objects with the two PSFs, thus obtaining six noise-free images, to which we added a background in K-band, corresponding to about 13.5 mag arcsec −2 . The resulting images have been corrupted with Poisson and additive Gaussian noise (read out noise -RON). In order to avoid saturation of the detector, the images have been obtained by co-adding 30 frames. It follows that the variance of the RON is 30σ 2 , thus corresponding to the RON of 30 frames; we set σ = 10 e − /px. According to the approach proposed in [17], RON compensation is obtained in the deconvolution algorithms by adding the constant 30σ 2 = 3000 to the image and the background. The final images are shown in figure 1.

Initialization and iteration numbers
A key point of the procedure is a good choice of the starting point (f (0) , h (0) ), since different initializations might lead to different limit point due to the nonconvexity of the objective function. In particular, since the first step of the alternating approach is a minimization on the image, a clever initialization for the PSF is crucial to address the sequence in the "basin of attraction" of the desired stationary point. While for f (0) it is rather common to start with a constant image (and several tests with different choices showed very similar results), the initial PSF h (0) has been chosen in two ways: the autocorrelation of the ideal PSF (possibly reiterated until the SR bound is satisfied) and the ideal PSF plus a small constant chosen to satisfy normalization and SR value. As concerns the choice of the number of iterations, for the inner subproblems we used N (k) f = 50 and N (k) h = 1 for all k, while for the outer cycle we performed 3000 iterations. We remark that for the specific content of our images we do not perform an early stopping of the iterations to achieve a regularized solution since the minimizers of the KL divergence are sparse objects (nightsky [18] or checker-board [19] effect). However, we recall that our approach can be easily extended to incorporate regularization terms in the objective function, as a ℓ 1 constraint on the image (inducing sparsity) and/or a Tikhonov term on the PSF (inducing smoothness). In these cases, some Poisson-based criterion for the choice of the parameter(s) would be necessary (see e.g. [20,21]).

Results
The accuracy of the reconstructions has been evaluated with different indices, namely the MARE for the image and the RMSE for the PSF, defined respectively as where q is the number of stars, m i andm i are the reconstructed and the true magnitudes, h andh are the reconstructed and true PSFs and · 2 denotes the Euclidean norm.
The results are shown in table 1 while in figure 2 we reported the reconstructed PSFs for the SC200 dataset. In particular, we can notice that the algorithm is able to recover the stars with a very high accuracy, with values of the MARE which are always below 0.1%. At the same time, both the PSFs are also well recovered, with reconstruction errors varying around a few percent. As concerns the comparison between the two initializations, our results show that the autocorrelation of the ideal PSF seems to provide better reconstructions of the PSF for the higher SR, as attested both by the lower RMSE in table 1 and by the plots in figure 2, in which we can notice a better recovery of the PSF structures. No significant differences are present in all the three tests when SR = 0.60. We have to remark that we always observed a convergent behaviour in the objective function while, in few cases, the errors of the PSF are still decreasing after 3000 iterations, thus indicating that a larger number of iterations could still improve the solutions.

Conclusions and future work
In this paper we present a general optimization approach for the minimization of the KL divergence in the case of blind deconvolution of images corrupted by Poisson noise. The minimization scheme we analyzed exploits the separability of the constraints, providing subsequences which converge to stationary points. Moreover, it is particulary suited for the feasible sets arising in the case studied, made up by box constraints combined with one equality constraint. The application to the reconstruction of star cluster images blurred by LBT-like PSFs seems to confirm that the method is able to reconstruct both the target distribution and the unknown PSF with very high accuracy and paves the way to possible integrations into specific codes developed for accurate photometric and astrometric analysis of star clusters, as the so-called StarFinder [22]. Future work will involve also the generalization of the approach to multiple-image problems as the co-adding of images with different PSFs [23] or the blind deconvolution of multiple images provided by a Fizeau interferometer [24]. Moreover, effective strategies for the reduction of undesired boundary effects [25,26] will be employed, in order to be able to reconstruct effectively clusters with stars positioned near the edges of the image.