;+ ; 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 ; Use COMPLEMENT keyword to WHERE() W. Landsman Jan 2008 ;- pro Max_Likelihood, data, psf, deconv, Re_conv, FT_PSF=psf_ft, NO_FT=noft, \$ GAUSSIAN=gaussian, \$ POSITIVITY_EPS=epsilon, \$ UNDERFLOW_ZERO=under compile_opt idl2 if N_elements( deconv ) NE N_elements( data ) then begin deconv = data deconv[*] = total( data )/N_elements( data ) Re_conv = 0 endif if N_elements( under ) NE 1 then under = 1.e-22 if N_elements( epsilon ) NE 1 then epsilon = -1 if N_elements( Re_conv ) NE N_elements( deconv ) then \$ Re_conv = convolve( positivity( deconv, EPS=epsilon ), psf, \$ FT_PSF=psf_ft, NO_FT=noft ) if keyword_set( gaussian ) then begin deconv = deconv + convolve( data - Re_conv, psf, /CORREL, \$ FT_PSF=psf_ft, NO_FT=noft ) endif else begin wp = where( Re_conv GT under, npos, \$ ncomplement=nneg,complement=wz) if (npos GT 0) then Re_conv[wp] = ( data[wp]/Re_conv[wp] ) > 0 if (nneg GT 0) then Re_conv[wz] = 1. deconv = deconv * convolve( Re_conv, psf, FT_PSF=psf_ft, \$ /CORREL, NO_FT=noft ) endelse if N_params() GE 4 then \$ Re_conv = convolve( positivity( deconv, EPS=epsilon ), psf, \$ FT_PSF = psf_ft, NO_FT = noft ) end