The FMAPE Algorithm

We use the Bayesian strategy to obtain an iterative algorithm for image restoration. The application of Bayes' theorem to the image restoration problem gives

The most probable image , given data , is obtained by maximizing in Eq. (2) or the product since is constant.

The conditional probability describes the noise in the data and its possible object dependence. It is fully specified in the problem by the likelihood function. As indicated above, we have two processes: the first to form the image on the detector array and the second to read the detector. Taking the whole process into account, the compound likelihood is (Núñez and Llacer 1993)

and its logarithm is

The authors first introduced this compound likelihood for restoration of Poisson data in the presence of read-out noise (Llacer and Núñez 1990). Snyder et al. (1993) use the same likelihood form for CCD cameras. If the process were pure Poisson (no read-out noise), the logarithm of the likelihood would be the classical expression

We use entropy to define the prior probability in a generalization of the concepts originally described by Frieden (1972). Let be the total energy (usually the number of counts or photons) in the object. Assume that there is an intensity increment describing an intensity jump to which we can assign some appropriate physical meaning.

Assume, in addition, that we have prior information regarding the statistical makeup of the object obtained for example, from previous observations or from observations at other wavelengths. If is the prior energy distribution, the prior probability of a photon going to pixel is given by

In that case, the number of ways that the object can occur is

We take the prior probability of the object to be proportional to its multiplicity, as given by Eq. (5). Using Stirling's approximation and taking into account that both and are constants, it is easy to show that the logarithm of the probability is

This expression is the Shannon form of entropy with the inclusion of the parameter and the spatial prior information . This form of the entropy is also called cross-entropy.

The Bayesian function to be maximized is

where is given by

and is a Lagrange multiplier for the conservation of counts. Note that the relative weight of the two terms in the Bayesian function Eq. (6) is controlled by .

To obtain the maximum of Eq. (6), we set and apply the method of successive substitutions, which affords us greater flexibility than other methods and results in rapidly converging algorithms. The Bayesian maximum a posteriori algorithm with entropy prior (FMAPE) is given by the iterative formula

where

In Eq. (7) is the index of the iteration, is a constant to preserve the energy in the form , computed at the end of each iteration; is an arbitrary constant (usually ) to secure positivity in the solution; is a constant to accelerate convergence up to approximately three times (). Constants and do not affect the point to which the algorithm converges. The parameter determines the point of convergence of the solution between the trivial solution and the maximum likelihood solution. is thus the balancing parameter of the Bayesian approach and it is the inverse of the weight of the entropy versus likelihood.

The iterative algorithm in Eqs. (7) and (8) has a number of desirable characteristics: (1) it solves the cases of both pure Poisson data and Poisson data with Gaussian read-out noise, (2) it maintains positivity of the solution, (3) it is easy to implement, (4) it includes case-specific prior information (default map) and flat field corrections, and (5) it removes background and can be accelerated to be faster than the Richardson-Lucy algorithm (Lucy 1974). The main loop (projection and back-projection) of the algorithm is similar in nature to the Expectation Maximization algorithm. The algorithm can be applied to a large number of imaging situations, including CCD and pulse counting cameras both in the presence and in the absence of background.

We can obtain a maximum likelihood algorithm from the general FMAPE by taking the limit when . The algorithm for the maximum likelihood (MLE) case in the presence of read-out noise and background is

where is given by Eq. (8).

In the case of no background and no read-out noise, Eq. (9) becomes

For and disregarding the gain (flat field) corrections (), Eq. (10) is identical to the Richardson-Lucy algorithm.