Holmes (1992) has derived an iterative algorithm for finding the maximum likelihood solution for both the unblurred image and the PSF in the presence of Poisson statistics. This problem is usually referred to as blind deconvolution. A maximum likelihood approach to blind deconvolution is certain to fail unless one can place very strong constraints on the properties of the PSF. For example, if one applies the Holmes iteration to HST images with no constraints on the PSF, the solution found is that the unblurred image looks exactly like the data and the PSF is a delta-function. (A perfect fit to the data, but ludicrously far from the truth!) If one applies a band-limit constraint on the PSF based on the size of the HST aperture, the result is little better: the PSF solution looks like a diffraction-limited PSF rather than the true spherically aberrated PSF.
If one has either a star within the image or a separate PSF observation, one can do much better by forcing the selected PSF objects to be point sources in the restored image. One simple approach is to start with an initial guess for the image that has a delta-function at the position of the star and is zero in the other nearby pixels. Then it is really possible to get estimates of both the unblurred image and the PSF from the Holmes iteration.
This may look to be of only academic interest, but in fact it may be the solution to an important problem for HST images. Most HST deconvolutions are carried out using observed PSFs, because theoretically computed PSFs (using the Tiny TIM software, for example) are usually not a very good match to the observations. The noise in observed PSFs presents a problem, however: when using a noisy PSF, the restoration algorithm ought to account for the fact that both the image and the PSF are noisy. The Richardson-Lucy method (and other commonly used image restoration methods) incorrectly assume that the PSF is perfectly known. The blind deconvolution approach allows one to construct an algorithm that takes as input a noisy image and a noisy PSF observation and to construct as output estimates of both the PSF and the unblurred image.
Early experiments indicate that this method may be especially useful in
crowded stellar fields where it is difficult to find isolated PSF
stars. In that case one can start with a very noisy PSF star (either
on the edge of the crowded region or from a separate observation) and
then can improve the PSF estimate using the overlapping PSFs of stars
in the crowded region (Fig. 4). The final PSF is
considerably more accurate than the original noisy PSF, which in turn leads to better restorations of the unblurred image.