Table of Contents -- Index -- PS reprint -- PDF reprint

Astronomical Data Analysis Software and Systems VII
ASP Conference Series, Vol. 145, 1998
Editors: R. Albrecht, R. N. Hook and H. A. Bushouse

Positive Iterative Deconvolution in Comparison to Richardson-Lucy Like Algorithms

Matthias Pruksch and Frank Fleischmann
OES-Optische und elektronische Systeme GmbH, Dr.Neumeyer-Str.240, D-91349 Egloffstein, F.R. Germany

 

Abstract:

Positive iterative deconvolution is an algorithm that applies non-linear constraints, conserves energy, and delivers stable results at high noise-levels. This is also true for Richardson-Lucy like algorithms which follow a statistical approach to the deconvolution problem.

In two-dimensional computer experiments, star-like and planet-like objects are convolved with band-limited point-spread functions. Photon noise and read-out noise are applied to these images as well as to the point-spread functions being used for deconvolution. Why Richardson-Lucy like algorithms favor star-like objects and the difference in computational efforts are discussed.

         

1. Introduction

The increasing availability of computing power fuels the rising interest in iterative algorithms to reconstruct an unknown object distribution O(x), that was blurred by a linear system's point-spread function P(x). The measured image distribution I(x) is then known to be

 

I(x) = P(x) * O(x),

with * denoting convolution.

In the case of band-limited point-spread functions or point-spread functions with incomplete coverage of the Fourier domain (interferometry), information is lost and therefore, deconvolution is not possible. Instead of one unique solution, a space of distributions solves (1) (Lannes, Roques & Casanove 1987). This raises the question of how the algorithms choose their reconstruction out of the space of possible solutions.

2. Algorithms

Richardson-Lucy like algorithms use a statistical model for image formation and are based on the Bayes formula  

$$p(P*O\vert I) = \frac{p(I\vert P*O)\ p(P*O)}{p(I)},$$ (2)

where p(P*O|I) denotes the probability of an event at (P*O)(x), if an event at I(x) occured. Müller (1997) optimized (2) for O(x) by functional variation. For the case of setting the probability p(O(x)) to Poisson statistic the optimization leads to the algorithm invented by Richardson (1972) and Lucy (1974)

$$\hat O_{i+1}(x) = \hat O_i(x) \left[P(-x)* \frac{I(x)}{P(x)*\hat O_i(x)} \right],$$ (3)

with $\hat O_i(x)$ denoting the iterated object distribution for reconstruction. Setting the probability distribution to have Gauss statistics (Müller 1997) leads to

$$\hat O_{i+1}(x) = \hat O_i(x) \frac{I(x)* P(-x)}{[P(x)*\hat O_i(x)]*P(-x)}.$$ (4)

Both algorithms intrinsically apply the positivity constraint and conserve energy.

Positive iterative deconvolution expands the deconvolution to a converging sum. In between iteration-steps the positivity constraint is applied such that the equation

$$I(x) = \hat O_i(x)*P(x)+R_i(x)$$ (5)

is always satisfied, with R_i(x) denoting the residual of the sum. (Pruksch & Fleischmann 1997)

3. Simulation

In order to test the algorithms for astronomical objects, the simulated object distribution O(x) consists of one planet-like object, and five star-like objects as shown in Figure 1 on the left. The planet-like object in the center is a picture of minor planet Ida . The star-cluster on the lower left has brightness ratios of 1:2:3:4 (right:left:upper:center). The star in the upper right is ten times brighter than the weakest star in the cluster and allows to determine the point-spread function.

  

Figure 1: Simulated object on the left, and point-spread function that simulates aberrations by atmospheric turbulence of a ground-based telescope on the right (short exposure).

All pictures consist of 256 by 256 pixels and show Ida in full contrast (gamma correction 0.7). The smaller picture on the right side of each picture shows phase (no gamma correction) and modulus (gamma correction 0.25) of the corresponding spectrum. The spectrum is shown for positive horizontal frequencies only, since the values of the negative frequencies are the complex conjugate of the positive ones. The origin of the phase is located on the left center and the origin of the modulus in the center of the small picture.

Different band-limited point-spread functions have been applied to the object, but only one is discussed and illustrated in Figure 1 on the right. The point-spread function simulates wavefront aberrations of approximately $\pm 3\pi$. This is typical for a ground-based telescope that suffers from atmospheric turbulence. The band-limit is visible in the spectrum, since any value beyond the circle is zero (grey for phase and black for modulus).

  

Figure 2: Convolved image with noise consistent to 10⁸ photons per image on the left, and reconstruction by Richardson-Lucy algorithm on the right.

The peaks in the point-spread function are best seen in Figure 2 on the left which shows the convolved image I(x). In order to avoid spurious information due to round-off errors, the point-spread function was calculated and convolved with the object in the Fourier domain. The contributions beyond the band-limit in Figure 2 on the left originate from noise that was added to the image after convolution.

4. Results

The number of iterations needed by the algorithms to deconvolve the images was in the range of 20 to 100. The computational effort per iteration is dominated by the number of Fourier transforms. The algorithm by Müller is the fastest by performing only two Fourier transforms per iteration, followed by positive iterative deconvolution needing three and Richardson-Lucy with four. To measure convergence, another two Fourier transforms are needed for the Richardson-Lucy like algorithms. If convergence is tested for every iteration-step then positive iterative deconvolution is the fastest algorithm.

All algorithms resolve the objects, as shown by the reconstructions by Richardson-Lucy in Figure 2 on the right, the algorithm by Müller in Figure 3 on the left and positive iterative deconvolution in Figure 3 on the right.

  

Figure 3: Reconstruction by Müller algorithm on the left, and reconstruction by positive iterative deconvolution on the right.

The Richardson-Lucy like algorithms seem to perform better in the reconstruction of stars, whereas positive iterative deconvolution shows more details for the planet-like object. This is due to the fact that positive iterative deconvolution uses only the positivity constraint, whereas the other algorithms allow statistical deviations. This is clearly seen by comparing the phase distribution of the reconstructions with that of the original object. Since the phase is the most important information, positive iterative deconvolution delivers in this sense the best solution.

5. Summary

Richardson-Lucy like algorithms lead to satisfactory reconstructions. Especially in situations with bad signal-to-noise ratio, the results are very stable. Due to the statistical process of optimization, a different phase distribution is reconstructed, which is favors star-like objects.

Positive iterative deconvolution reconstructs the original data as far as it is available in the measured image and adds information only consistent to that data and the positivity constraint. Therefore, it shows equal quality for all objects. Since positive iterative deconvolution is also the fastest algorithm it is the algorithm of choice, if the reconstructions have to be consistent with measurements.

Acknowledgments:

The work reported here has received financial support from the OES-Optische und elektronische Systeme GmbH.

References:

Lannes, A., Roques, S., & Casanove, M. J. 1987, J.Mod.Opt., 34, 161

Lucy, L. B. 1974, AJ, 79, 745

Müller, R. 1997, private communication

Pruksch, M. & Fleischmann, F. 1997, accepted by Comput.Phys.

Richardson, W. H. 1972, J.Opt.Soc.Am., 62, 55

© Copyright 1998 Astronomical Society of the Pacific, 390 Ashton Avenue, San Francisco, California 94112, USA