;+

function sxt_image_deconvolve,index,data,unc_data,satpix, $
                              psf=psf,nodecon=nodecon, $
                              delta=delta,moffat=moffat, $
                              slevel=slevel,spower=spower, $
                              no_shadow=no_shadow, $
                              log10t=log10t,totalunc=totalunc, $
                              fftguess=fftguess,scat_fraction=scat_fraction

;NAME:
;     SXT_IMAGE_DECONVOLVE
;PURPOSE:
;     Remove PSF from SXT data
;CATEGORY:
;CALLING SEQUENCE:
;     out = sxt_image_deconvolve(index,data,unc_data,satpix)
;INPUTS:
;     index = SXT index
;     data = SXT data cube
;     unc_data = output from sxt_prep
;     satpix = output from sxt_prep
;OPTIONAL INPUT PARAMETERS:
;KEYWORD PARAMETERS
;     /fftguess = Use the direct FFT inversion as the initial guess
;     /moffat = Use only the moffat core, no scattering wings 
;               (passed to sxt_psf)
;     /delta = Do not use the moffat core, only the scattering wings 
;                (passed to sxt_psf)
;     log10t = Log10(temperature) passed to sxt_psf (def = log(3.e6)).
;     scat_fraction, slevel, spower = passed to sxt_psf.pro.
;     /totalunc = indicates the the uncertainty array passed in should be
;                 used as is and not combined with the sxt_dn_unc program.
;                 If the uncertainty array is taken straight from sxt_prep,
;                 do *not* use this keyword.
;     /nodecon = skips the deconvolution.  You really don't want to use
;                this, it is only for testing.
;     psf = returns with the PSF used in the deconvolution.
;OUTPUTS:
;     out = corrected data cube
;COMMON BLOCKS:
;SIDE EFFECTS:
;RESTRICTIONS:
;     The data input must be corrected with sxt_prep first!
;PROCEDURE:
;     Uses image_deconvolve.pro (a simple MEM algorithm)
;MODIFICATION HISTORY:
;     T. Metcalf  2000-10-04
;                 2001-01-15 Fixed minor indexing bug.
;                 2001-01-16 Added check for sigma>1.0
;                 2001-04-24 Fixed to work with the new sxt_psf
;-

nimage = n_elements(index)
out = float(data)

for iimage = 0L,nimage-1L do begin

   res = 2^gt_res(index[iimage])
   guard = 64/res

   ; Get PSF
   nx = n_elements(data(*,0,iimage))+guard
   ny = n_elements(data(0,*,iimage))+guard
   xy = double(gt_center(index[iimage])) ; force PSF calculation to be double
   filter = gt_filtb(index[iimage])
   psfdim = 1024 < (nx*res) < (ny*res)
   psfsmall = sxt_psf(filter,xy[0],xy[1],dim=psfdim, $
                      delta=delta,moffat=moffat, $
                      slevel=slevel, spower=spower, no_shadow=no_shadow, $
                      res=0,log10t=log10t,scat_fraction=scat_fraction)
   ; Rebin PSF to correct resolution
   ; rebin uses neighbor averaging, which is the right thing to do.
   psfsmall = rebin(psfsmall,psfdim/res,psfdim/res)  
   pnx = n_elements(psfsmall(*,0))
   pny = n_elements(psfsmall(0,*))
   nx2 = floor(nx/2.)
   ny2 = floor(ny/2.)
   pnx2 = floor(pnx/2.)
   pny2 = floor(pny/2.)
   psf = dblarr(nx,ny)
   psf[nx2-pnx2,ny2-pny2] = psfsmall   ; Center PSF
   psf = float(psf/total(psf))  ; normalize PSF

   ; Estimate the error
   if keyword_set(log10t) then te00 = log10t $
   else te00 = alog10(3.e6)  ; Using fixed value makes this approximate
   ph_per_dn = sxt_flux(Te00,index[iimage],/photons,/nover) / $
               sxt_flux(Te00,index[iimage],gain_ccd=gain_ccd,/nover)
   ;;convert = ph_per_dn*gt_expdur(index,/orig)/gt_expdur(index)
   ;;sigma = sqrt((data*convert)>1.)/convert
   if keyword_set(totalunc) then sigma = unc_data[*,*,iimage] $
   else $
      sigma = float(sxt_dn_unc(index[iimage],data[*,*,iimage], $
                               unc_data[*,*,iimage]))

   sssat = where(satpix[*,*,iimage])
   if sssat[0] GE 0 then mask = sssat else delvarx,mask

   if n_elements(mask) GT 0 then begin
      ; Value is unknown so give it a large error
      if mask[0] GE 0 then begin
         dtemp = data[*,*,iimage]
         sigma[mask] = dtemp[mask]*2.0 
         delvarx,dtemp
      endif
   endif

   ; Bigger areas that include the guard band
   old = fltarr(nx,ny)
   unc = fltarr(nx,ny)
   mas = intarr(nx,ny)
   old[*] = 0.
   old[0,0] = data[*,*,iimage]
   unc[*] = 1000.              ; for the guard band
   unc[0,0] = sigma
   mas[*] = 1                  ; Mask out the guard band
   mas[0,0] = satpix[*,*,iimage]
   mask = where(mas)

   ; Do the deconvolution
   print,!stime


   if NOT keyword_set(nodecon) then begin
      if NOT keyword_set(positive_definite) then begin
         pedestal = 0. > (5.-min(data[*,*,iimage]))
         message,/info,strcompress('Pedestal='+string(pedestal))
      endif else pedestal = 0.0
      bad = where((old+pedestal) LE 0.0,nbad)
      if nbad GT 0 then begin
         message,/info,strcompress('Found ' + string(nbad) + $
                                   ' pixels LE 0.  Setting to 0.01')
      endif
      ; Use photons rather than DN with usepoisson=1 to get correct Poisson 
      ; statistics in the image.  The convert back to DN.  Choose one of
      ; the following two lines:

      ;ph_per_dn = 1 & usepoisson=0   ; For Guassian statistics
      usepoisson = 1                  ; For Poisson statistics

      new = image_deconvolve(((old+pedestal)>0.01)*ph_per_dn, $
                             psf,(unc>1.)*ph_per_dn,uselog=1,mask=mask, $
                             usepoisson=usepoisson,flux=1,usepsfcorr=1, $
                             fftguess=fftguess)/ph_per_dn $
            - pedestal
      print,!stime   
      out[*,*,iimage] = new[0:nx-guard-1,0:ny-guard-1]
   endif else message,/info,'Skipping deconvolution'

endfor

return,out

end