FORWARD_FUNCTION DIcal_mtf, getMTFfilter
;;-----------------------------------------------------------------------------
;; PURPOSE:
;;	Used in the DI Calibration Pipeline.
;;	Applies a deconvolution to compensate for smearing and out-of-focus images
;;
;; CALLING SEQUENCE:
;;	out = DIcal_mtf(in, fitsHdr, flags, FN=fn, ALG=alg, PARAM=param)
;;
;; REQUIRED INPUTS:
;;	in - The image to which this pipeline element is going to be applied
;;	fitsHdr - The FITS header for the image
;;	flags - the flags specifying info about the pixels
;;
;; OUTPUTS:
;;	RETURN - the image after going through this calibration step
;;
;; OPTIONAL INPUT KEYWORDS:
;;	FN - The filename containing the psf. The psf is
;;		stored as an image in FITS format. If this is not specified,
;;		it is found on a database.
;;	ALG - deconvolution algorithm to use.
;;		1 = Don Lindler's least squares solution (default because of speed)
;;		2 = Richardson Lucy with wavelet noise suppression
;;	PARAM - parameter specific to the algorithm being used. For the algorithm:
;;		1: This represents the reciprocal lagrange multiplier of the trial (gammat)
;;		2: TODO
;;
;; EXAMPLE:
;;      IDL> imgOut = DIcal_mtf(imgIn, fitsHdr, flags)
;;
;; PROCEDURES USED (i.e. called directly!):
;;	SXPAR - Gets FITS header items
;;
;; MODIFICATION HISTORY:
;;   2004-05-24  M. Desnoyer    Created
;;   2005-04-12  M. Desnoyer    Reordered functions
;;   2005-04-26  M. Desnoyer	Added ALG and PARAM
;;
;;-----------------------------------------------------------------------------

;;-----------------------------------------------------------------------------
;; PURPOSE:
;;	Retrieves the convolution filter for MTF and scattered light correction
;;	The filter is in the spatial domain. This is in other words, a PSF
;;
;; CALLING SEQUENCE:
;;	out = getMTFfilter(fitsHdr, inst, date, filt, FN=fn)
;;
;; REQUIRED INPUTS:
;;	fitsHdr - The FITS header for the image
;;	inst - the instrument that took the image
;;	date - the date the image was taken
;;	filt - number of the filter used
;;
;; OUTPUTS:
;;	RETURN - the MTF and scattered light correction filter in the spatial
;;		domain
;;
;; OPTIONAL INPUT KEYWORDS:
;;	FN - Filename of the convolution filter to use. If this is not 
;;		specified, a database is used.
;;
;; EXAMPLE:
;;      IDL> filt = getMTFfilter(hdr, inst, date, filt)
;;
;; PROCEDURES USED (i.e. called directly!):
;;
;;
;; MODIFICATION HISTORY:
;;   2004-07-16  M. Desnoyer    Created
;;   2005-05-27  M. Desnoyer    Added filter to the search
;;
;;-----------------------------------------------------------------------------
FUNCTION getMTFfilter, fitsHdr, inst, date, filt, FN=fn

ON_ERROR, 2

;; Get the filename
IF not keyword_set(fn) THEN BEGIN
	serv = getSQLserver()
	db = 'di_calib'
	tbl = 'PSF'
	sel = ['Filepath']
	cond = 'Instrument = "'+inst+'" AND Date <= "'+date+$
		'" AND Filter = '+filt+' order by Date desc, Version desc limit 1'

	webfn = di_sqlquery(serv, db, tbl, sel, cond, dim=dim)

	IF min(dim) EQ 0 THEN message, 'Could not Find Calibration File', /NONAME

	fn1 = getLocFn(webFn, subdir='PSF')
ENDIF ELSE fn1=fn

;; Open the filter
filt =  readFITS(fn1[0], calHdr, /SILENT)

IF (size(filt, /n_dimensions)) EQ 0 THEN $
	message, !ERROR_STATE.MSG, /NONAME

;; Change the fits header
fxaddpar, fitsHdr, 'DECONPSF', $
	strupcase(strmid(fn1[0], strpos(fn1[0], path_sep(), /reverse_search)+1))

RETURN, filt

END

;;-----------------------------------------------------------------------------
;; Main function
;;-----------------------------------------------------------------------------
FUNCTION DIcal_mtf, in, fitsHdr, flags, FN=fn, ALG=alg, PARAM=param

;; Do error handling 
;CATCH, error
error=0

IF error NE 0 THEN BEGIN
	CATCH, /CANCEL
	message, 'Deconvolution Correction - ' + !ERROR_STATE.MSG, /NONAME
ENDIF

;; Default inputs
IF NOT keyword_set(ALG) THEN algUsed = 1 ELSE algUsed=alg
IF NOT keyword_set(param) THEN BEGIN
	CASE algUsed OF
		1: algPar = 0.001	;; Gammat
		2: algPar = 100
		ELSE: message, 'Invalid algorithm chosen'
	ENDCASE
ENDIF ELSE algPar = param

@dical_flags

out = double(in)

;; Get the date
date = SXPAR(fitsHdr, 'OBSDATE', count=c1)

;; Get the instrument
inst = SXPAR(fitsHdr, 'INSTRUME', count=c2)

;; Get the filter
filt = SXPAR(fitsHdr, 'IMGH036', count=c3)

;; Make sure everything was in the FITS header
IF (c1 EQ 0) OR (c2 EQ 0) OR (c3 EQ 0)$
	THEN message, 'Invalid FITS header', /NONAME

;; Convert the hall sensor reading to the real filter number
filt = imgh036ToFilter(filt, inst, date)

;; Replace all the gaps with 0s
IF n_params() EQ 3 THEN BEGIN
	gapLoc = where((flags AND FLAG_GAP) NE 0, gapCnt)
	IF gapCnt GT 0 THEN out[gapLoc] = 0
ENDIF

;; Get the psf
psf = getMTFfilter(fitsHdr, inst, date, filt, FN=fn)
psf = psf/total(psf)

;; If the psf has a value of 1 in it, then it is a trivial psf and we shouldn't waste
;; time trying to deconvolve
IF max(psf) GE 1.0 THEN BEGIN
	fxaddpar, fitsHdr, 'DECONPSF', ''
	return, in
ENDIF
;; Pad the image to make its dimensions a power of two and thus make 
;; the fourier transforms faster
inDim = size(in,/dimensions)
dim = 2^ceil(alog(inDim)/alog(2))
img = dblarr(dim[0],dim[1])

offset = dim/2-inDim/2
;; Center the image
IF min(offset) GT 0 THEN BEGIN
	img[offset[0],offset[1]] = out
	;; Pad around the image by repeating the edge pixels
	n = offset+inDim
	img[offset[0],0] = rebin(img(offset[0]:n[0]-1,offset[1]),inDim[0],offset[1])
	img[offset[0],n[1]] = rebin(img(offset[0]:n[0]-1,n[1]-1),inDim[0],dim[1]-n[1])
	img[0,0] = rebin(img(offset[0],*),offset[0],dim[1])
	img[n[0],0] = rebin(img(n[0]-1,*),dim[0]-n[0],dim[1])
ENDIF ELSE img = out

;; Force the image to be positive to help out the deconvolution routine
img = img>0

;; Perform the algorithm
CASE algUsed OF
1:BEGIN	;;Constrained Least Squares
	conlsq2, img, psf, decon, trial=img, gammat=algPar
END
2:BEGIN	;;Richardson-Lucy with wavelet denoising

	;; Hard code these values for now because I don't understand them well enough
	gain = 28.5
	constant_noise = 0.77
	biasCols = getBiasCols(fitsHdr)>0
	con_var = constant_noise^2*gain+median(in[biasCols,*])

	;; Force positivity after noise
	img = img>(-con_var)

	;; Make the psf bigger
	psf_big = dblarr(dim[0],dim[1])
	psf_sz = size(psf,/dimensions)
	psf_big[0,0] = psf
	psf_big = shift(psf_big,-psf_sz[0]/2,-psf_sz[1]/2)
	fft_psf = fft(psf_big,-1)

	decon = replicate(mean(img),dim[0],dim[1])
	
	;; Actually call the routine
	dhwaveiter, decon, 1, img, gain, con_var, fft_psf, algPar, /init, speedup=2, dotv=0
	
END
ELSE:	message, 'Invalid algorithm chosen'
ENDCASE

out = decon[offset[0]:n[0]-1,offset[1]:n[1]-1]

;; Update the FITS header
fxaddpar, fitsHdr, 'DECON', 'T'
CASE algUsed OF
	1: algName = 'Constrained Least Squares'
	2: algName = 'Richardson-Lucy with Wavelet Denoising'
	ELSE: message, 'Should not get here'
ENDCASE
fxaddpar, fitsHdr, 'DECONALG', algName
fxaddpar, fitsHdr, 'DECONV', algPar

RETURN, out
END