idlastro / FITS I/O: MAX_LIKELIHOOD

[Source code]

NAME
MAX_LIKELIHOOD
PURPOSE
Maximum likelihood deconvolution of an image or a spectrum.
EXPLANATION
Deconvolution of an observed image (or spectrum) given the 
instrument point spread response function (spatially invariant psf).
Performs iteration based on the Maximum Likelihood solution for
the restoration of a blurred image (or spectrum) with additive noise.
Maximum Likelihood formulation can assume Poisson noise statistics
or Gaussian additive noise, yielding two types of iteration.
CALLING SEQUENCE
for i=1,Niter do Max_Likelihood, data, psf, deconv, FT_PSF=psf_ft
INPUTS PARAMETERS
data = observed image or spectrum, should be mostly positive,
                        with mean sky (background) near zero.
psf = Point Spread Function of the observing instrument,
                        (response to a point source, must sum to unity).
INPUT/OUTPUT PARAMETERS
deconv = as input: the result of previous call to Max_Likelihood,
        (initial guess on first call, default = average of data),
        as output: result of one more iteration by Max_Likelihood.
Re_conv = (optional) the current deconv image reconvolved with PSF
        for use in next iteration and to check convergence.
OPTIONAL INPUT KEYWORDS
/GAUSSIAN causes max-likelihood iteration for Gaussian additive noise
          to be used,  otherwise the default is Poisson statistics.
  FT_PSF = passes (out/in) the Fourier transform of the PSF,
          so that it can be reused for the next time procedure is called,
/NO_FT overrides the use of FFT, using the IDL function convol() instead.
  POSITIVITY_EPS = value of epsilon passed to function positivity,
                  default = -1 which means no action (identity).
  UNDERFLOW_ZERO = cutoff to consider as zero, if numbers less than this.
EXTERNAL CALLS
function convolve( image, psf ) for convolutions using FFT or otherwise.
function positivity( image, EPS= ) to make image positive.
METHOD
Maximum Likelihood solution is a fixed point of an iterative eq.
(derived by setting partial derivatives of Log(Likelihood) to zero).
Poisson noise case was derived by Richardson(1972) & Lucy(1974).
Gaussian noise case is similar with subtraction instead of division.
NOTES
WARNING: The Poisson case may not conserve flux for an odd image size.  
This behavior is being investigated.
HISTORY
written: Frank Varosi at NASA/GSFC, 1992.
F.V. 1993, added optional arg. Re_conv (to avoid doing it twice).
Converted to IDL V5.0   W. Landsman   September 1997
se COMPLEMENT keyword to WHERE()   W. Landsman  Jan 2008