We have compared the algorithm of Richardson-Lucy (RLA) and the Image Space Reconstruction Algorithm (ISRA) for image deconvolution. These iterative algorithms keep the reconstructed image non-negative while maximizing the likelihood for a Poisson Process (RLA) and a Gaussian process (ISRA). We have shown that these algorithms can be analyzed in a similar way when written as modified gradient algorithms. This form permits a better insight of the two algorithms. In particular, it reveals the similar roles played by the variances of the signal in the algorithms.
Numerical results are shown for realistic simulated images corresponding to observations made with a telescope with a small phase aberration. For Poisson noise, no strong difference was found between the results of ISRA and RLA. This later algorithm gives slightly better results than ISRA, as expected, but ISRA's results are very good too. For Gaussian noise, ISRA gives much better results than RLA.
The main drawback of these algorithms is that noise increases when the iteration number is too high. They both tend to create spectral components outside the cut-off frequency of the telescope. Although not shown in the paper, this is also observed for noise-free images. This property derives from the non-linearity of the algorithm due to the non-negativity constraint.
Some form of regularization is therefore necessary. We consider here the elementary method that consists of stopping the iterations before instabilities appear in the solution. To optimally adjust the iteration number, we propose to use a comparison of RLA and ISRA with the Wiener filter. The method is founded on the minimization of the Euclidean distance E[k] between the modulus of the Fourier transform of the solution at a given iteration, and the equivalent quantity computed on the image reconstructed by a Wiener filter. In so doing, we obtain the value of k corresponding, for this criterion, to the optimum regularized non negative solution.
Our assumptions have been confronted to numerical simulations for two images defined with 104 and 106 photons. The results are satisfactory; the comparison ISRA/Wiener or RLA/Wiener gives the exact optimal iteration number (known for simulated data) for the 104 photons image. For the 106 photons image, that comparison leads to an iteration number smaller than the optimum; in that case, the resulting image is slightly over regularized.
The procedure we propose can then be summarized as follows. We first construct a Wiener filter founded on physical considerations concerning the level of additive noise in the data or the number of photons in the image. This is probably the most difficult step in our reconstruction procedure. Then, depending on the statistics of the noise, we use ISRA or RLA with an iteration number that makes E[k] minimum. The use of ISRA or RLA ensures a positive result; the comparison with the Wiener filter guarantees a good signal to noise ratio of the reconstructed image.
However, the simple limitation of the iteration number is not sufficient to obtain a perfect adjustment of all frequency components: a good fit in the low frequency range is always obtained with an excess of power in the high frequencies. The introduction of a penalty term (Tikhonov or entropy) will be necessary in practice to better fit the results of the Wiener filter for all frequencies.
The first tests we have made indicate that a slight improvement of the results is to be expected from that (paper in preparation).
In the present state of our technique, we cannot propose to the reader a software package that would allow an observational astronomer to immediately apply the technique to the images of the telescope. A good knowledge of the Wiener filter is necessary. However, we think that the expected results is worth the required effort in programming and physical understanding of the sources of noise.
Copyright The European Southern Observatory (ESO)