FUNCTION DIcal_VIS, imgIn, fitsHdr, modules, VERBOSE=verbose, SILENT=silent, $
	FORCE=force, MAXGAPSIZE=maxgapsize, DESMEARALG=desmearAlg, RAD=rad, $
	DN=dn, REVRAD=revRad, IOVERF=ioverf, REVIOVERF=revIoverf,  FLATFN=flatfn, $
	DARKFN=darkfn, USEDARKMODEL=usedarkmodel, COMPRESSFN=compressfn, $
	ADCFN=adcfn, BADPIXSFN=badpixsfn, GEOMFN=geomfn, MTFFN=mtffn, $
	VISCONSTFN=visconstfn, FLAGS=flags, GAINFN=gainfn, SPIKESIG=spikesig, $
	SPIKEITER=spikeiter, SPIKEBOX=spikebox, SPIKEMED=spikemed, $
	BADPIXSINTERP=badpixsinterp, MISSINGINTERP=missinginterp, $
	COMPRESSMETH=compressmeth, MTFALG=mtfalg, MTFPARAM=mtfparam, SNR=snr

;;-----------------------------------------------------------------------------
;; PURPOSE:
;;	Performs the calibration for a MRI/HRI/ITS VIS Deep Impact image.
;;	The following steps are performed:
;;	1. Calibrates the temperatures and voltages in the header
;;	2. Decompresses the image if appropriate
;;	3. Flags the saturated pixels
;;	4. Corrects uneven ADC bit weightingk
;;	5. Subtracts the dark frame
;;	6. Remove cross tal
;;	7. Normalizes quadrant gains
;;	8. Flat-field correction
;;	9. Remove transfer smear
;;	10. Flag bad pixels
;;	11. Convert from DN to requested units
;;	12. Fill in small gaps
;;	13. Remove cosmic ray hits via despiking
;;	14. Remove random gaussian noise
;;	15. Deconvolution
;;	16. Geometric calibration
;;
;;
;; CALLING SEQUENCE:
;;	RESULT = DIcal_VIS(in, fitsHdr, /VERBOSE, /SILENT, /FORCE, 
;;		MAXGAPSIZE=size, DESMEARALG=alg, INTERPBADPIXS=interp,
;;		/RAD, /IOVERF, FLATFN=flatfn, DARKFN=darkfn, /USEDARKMODEL)
;;
;; REQUIRED INPUTS:
;;	in - The raw image to calibrate
;;	fitsHdr - The FITS header for the image. The header will be modified
;;		throughout the calibration routine. It will be NxM string
;;		array on output where N is the number of FITS lines and M
;;		is the number of images begin outputed
;;	modules - An activeModules structure defining which modules to apply
;;
;; OUTPUTS:
;;	RETURN - this has the calibrated images. It is a cube whose last 
;;		dimension depends on which of the DN, RAD, REVRAD, IOVERF and
;;		REVIOVERF keywords are set. The images will output in the cube
;;		in the order the keywords are listed above.
;;
;; OPTIONAL INPUT KEYWORDS:
;;	DN - if we should output a calibrated image in DN
;;	RAD - if we should output a calibrated image in radiance units
;;	REVRAD - if we should output a reversibly calibrated image in rad units
;;	IOVERF - if we should output a calibrate image in I/F units
;;	REVIOVERF - if we should ouput a reversibly calibrate image in I/F
;;		units
;;	VERBOSE - flaged if status messages are to be displayed
;;	SILENT - flaged if errors will not be displayed
;;	FORCE - flaged to force calibration to continue even if a step fails.
;;		If this is not set, then if an error occurs in any step,
;;		the calibration fails.
;;	MAXGAPSIZE - the maximum dimension of a gap to fill with interpolated 
;;		data (default = 9999)
;;	DESMEARALG - the algorithm to use for desmearing. 
;;		0 = No Desmearing
;;		1 = Use POC values (default for modes 1-6,9)
;;		2 = Column averaging (default for modes 7,8)
;;	FLATFN - Filename of the FITS image which is a flat field to apply to 
;;		the image. If this is not set, a default flat field is found
;;	DARKFN - Filename of the FITS image which is a dark frame to apply to 
;;		the image. If this is not set, a default dark is found
;;	USEDARKMODEL - If a dark frame is not supplied and this keyword is set
;;		then a pre-fit dark model is used to generate the dark frame
;;	COMPRESSFN - filename of the decompression LUT to use
;;	ADCFN - filename of the bit weighting correction LUT to use
;;	BADPIXSFN - filename of the bad pixel map to use
;;	GEOMFN - filename of the geometric tie points to use
;;	MTFFN - filename of the convolution filter to use
;;	VISCONSTFN - filename of the calibration constant to use
;;	FLAGS - variable to hold an byte of flags for each pixel. This may be an image 
;;		if more than one output is requested. The order is that same as the main
;;		output
;;	GAINFN - filename of the gain map to use
;;	SPIKESIG - Threshold for the despiking routine
;;	SPIKEITER - Maximum number of iterations for the despiking routine
;;	SPIKEBOX - How big a box to use to calculate statistics for despiking
;;	SPIKEMED - Should despiking use Median? Otherwise it will use mean
;;	BADPIXSINTERP - Should bad pixels be interpolated over
;;	MISSINGINTERP - Should missing data be interpolated over
;;	COMPRESSMETH - The method to use for decompression
;;	MTFALG - Algorithm used for deconvolution 
;;		1=constrained least squares, 2=Richardson-Lucy with wavelets
;;	MTFPARAM - Algorithm specific parameter for the deconvolution algorithm
;;	SNR - Output map of the signal to noise ratio of the image
;;
;; EXAMPLE:
;;      IDL> imgOut = DIcal_VIS(imgIn, fitsHdr, /VERBOSE, /FORCE)
;;
;; PROCEDURES USED (i.e. called directly!):
;;	SXPAR - Read a FITS parameter
;;	DIcal_decompress - decompresses the image
;;	DIcal_flagSat - flags the saturated pixels
;;	DIcal_flagGaps - flags gaps in data
;;	DIcal_bitWeight - corrects for uneven bit wighting
;;	DIcal_darkFrame - subtracts a dark frame
;;	DIcal_flatField - performs flat field correction
;;	DIcal_desmear - removes transfer smear
;;	DIcalVIS_DNtoRad - converts the units from DN to radience
;;	DIcal_fixBadPixs - fixes bad pixels
;;	DIcal_despike - removes noise and spikes
;;	DIcal_fillGaps - fills in small gaps
;;	DIcal_denoise - removes gaussian noise
;;	DIcal_geom - performs geometric calibration
;;	DIcal_mtf - performs MTF and scattered light corrections
;;	DIcal_RadtoIF - Converts from radiance to I/F units
;;	DIcal_gain - Normalizes quadrant gain
;;
;; MODIFICATION HISTORY:
;;   2004-05-24  M. Desnoyer    Created
;;
;;-----------------------------------------------------------------------------

doverbose = keyword_set(verbose)

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

IF error NE 0 THEN BEGIN
	fn = SXPAR(fitsHdr, 'FILENAMER', COUNT=c1)
	IF c1 EQ 0 THEN fn=''

	;; Throw an error if an error causes things to stop
	IF not keyword_set(FORCE) THEN BEGIN
		CATCH, /CANCEL
		IF keyword_set(SILENT) THEN $
			message, /NOPRINT
		message, !ERROR_STATE.MSG, /NONAME
	ENDIF

	;; Print a message
	IF not keyword_set(SILENT) THEN $
		PRINT, fn+': '+!ERROR_STATE.MSG
	
	;; Go to the next calibration step if errors don't stop the pipeline
	;; Also flag that that step wasn't applied
	CASE curStep OF
	0: GOTO, Decompress
	1: GOTO, FlagSat
	2: GOTO, ADCBit
	3: GOTO, Darks
	4: GOTO, XTalk
	5: GOTO, Gain
	6: GOTO, Flatfield
	7: GOTO, Desmear
	8: GOTO, FindBad
	9: GOTO, ChangeUnits
	10: BEGIN
		CATCH, /cancel
		message, !ERROR_STATE.MSG, /noname
	END
	11: GOTO, Despike
	12: GOTO, Denoise
	13: GOTO, MTF
	14: GOTO, Geometric
	15: GOTO, Output
	ENDCASE
ENDIF

IF n_params() NE 3 THEN message, 'Invalid number of parameters', /NONAME

;; Check to make sure that RAD and/or IOVERF are flagged
dn = keyword_set(dn)
ioverf = keyword_set(ioverf)
revIoverf = keyword_set(revIoverf)
rad = keyword_set(rad)
revrad = keyword_set(revrad)
IF (NOT (ioverf OR rad OR dn OR revIoverf OR revRad)) THEN $
	message, 'An output type must be specified', /NONAME

out = double(imgIn)

;; Create an empty flag array
flsz = size(flags, /dimensions)
dim = size(imgIn, /dimensions)
IF max(flsz NE dim) THEN flags = bytarr(dim[0], dim[1])

;; Calibrate the temperatures and voltages in the header
curStep=0
IF modules.tempvolt THEN BEGIN
	fitsHdr = dical_tempvolt(fitsHdr)
	IF doverbose THEN print, 'Calibrated temps and voltages'
ENDIF

;; Uncompress the image
Decompress:
curStep = 1
mycmp = sxpar(fitsHdr, 'CMPRESSN')>0
IF (strtrim(sxpar(fitsHdr,'FNCMPRSS'),2) EQ 'C') AND (mycmp EQ 0) THEN $
	out = DIcal_decompress(out,fitsHdr, METHOD=compressmeth, FN=compressfn,$
		FLAGS=flags,NOISE=snr) 
IF doverbose THEN print, 'Image uncompressed'

FlagSat:
curStep = 2

;; Determine where the image, bias cols and POC rows are
POCrows = getPOCrows(fitsHdr,COUNT = cnt, COMPLEMENT = imgRows)
BIAScols = getBIAScols(fitsHdr, COUNT=cnt, COMPLEMENT = imgCols)

;; Flag the saturated pixels
IF modules.flagSat THEN BEGIN
	out = DIcal_flagSat(out,fitsHdr,FLAGS=flags) 
	IF doverbose THEN print, 'Identified saturated pixels'
ENDIF

;; Correct ADC bit weigting
ADCBit:
curStep = 3
IF modules.bitWeight THEN BEGIN
	out = DIcal_bitWeight(out,fitsHdr, FN=adcfn) 
	IF doverbose THEN print, 'Corrected uneven bit weighting'
ENDIF

;; Calculate the noise
snr = DIcal_calcNoise(fitsHdr,out,flags,QUANT=snr)

;; Subtract dark frame
Darks:
;; Identify whether or not to use the dark frame model TODO
curStep = 4
IF modules.dark THEN BEGIN
	out = DIcal_darkFrame(out,fitsHdr, flags, DARKFN=darkfn, USEFIT=usedarkmodel) 
	IF doverbose THEN print, 'Subtracted dark frame'
ENDIF

;; Calculate the signal to noise ratio
snr = abs(out / snr)

;; Remove the crosstalk
XTalk:
curStep=5
IF modules.xTalk THEN BEGIN
	out[imgCols[0],imgRows[0]] = DIcal_xtalk($
		out[imgCols[0]:imgCols[1],imgRows[0]:imgRows[1]],$
		fitsHdr,$
		FLAGS=flags[imgCols[0]:imgCols[1],imgRows[0]:imgRows[1]])
	IF doverbose THEN print, 'Removed cross-talk'
ENDIF

;; Normalize quadrant gains
Gain:
curStep=6
IF modules.gain THEN BEGIN
	out = DIcal_gain(out, fitsHdr, FN=GAINFN)
	IF doverbose THEN print, 'Normalized quadrant gains'
ENDIF

;; Perform the flat field correction
Flatfield:
curStep = 7
IF modules.flat THEN BEGIN
	out = DIcal_flatField(out,fitsHdr, FN=flatfn) 
	IF doverbose THEN print, 'Applied flat-field normalization'
ENDIF

;; Remove the transfer smear
Desmear:
curStep = 8
IF modules.desmear THEN BEGIN
	out = DIcal_desmear(out,fitsHdr, POCRows,imgRows, imgCols, flags, ALG=desmearAlg) 
	IF doverbose THEN print, 'Removed transfer smear'
ENDIF

;; Flag bad pixels
FindBad:
curStep = 9
IF modules.badPixs THEN BEGIN
	out = DICal_fixBadPixs(out,fitsHdr,FN=badpixsfn,FLAGS=flags)
	IF doverbose THEN print, 'Flagged bad pixels'
ENDIF

;; Convert from DN to I/F and/or radiance units
ChangeUnits:
curStep = 10
hdrSize = (size(fitsHdr,/dimensions))[0]
dnOut = out
dnHdr = fitsHdr

;; If we only want dn results, output now
IF NOT (rad OR revRad OR ioverf OR revIoverf) THEN goto, Output

;; Convert to physical units
out = dicalvis_DNtoRadIF(out,fitsHdr,FN=visconstfn, RAD=rad OR revRad, $
	IOVERF=ioverf OR revIoverf)
IF rad OR revRad THEN BEGIN
	radHdr = fitsHdr.rad
	revRadHdr = radHdr
	radOut = out.rad
	revRadOut = radOut
ENDIF
IF ioverf OR revIoverf THEN BEGIN
	iofOut = out.ioverf
	iofHdr = fitsHdr.ioverf
	revIofOut = iofOut
	revIofHdr = iofHdr
ENDIF
IF doverbose THEN print, 'Converted from DN to physical units'

;; Figure out which flag maps to use from here on
radFlags = flags
ifFlags = flags

;; Fill in small gaps using interpolation
FillGaps:
curStep=11
IF modules.fillGaps THEN BEGIN
	IF keyword_set(maxgapsize) EQ 0 THEN maxgapsize = 9999
	IF rad THEN BEGIN
		radOut = DIcal_fillGaps(radOut, radHdr, maxgapsize, radFlags,$
			imgCols, imgRows, BADPIX=badpixsinterp, $
			MISSING=missinginterp)
	ENDIF
	IF ioverf THEN BEGIN
		iofOut = DIcal_fillGaps(iofOut, iofHdr, maxgapsize, ifFlags,$
			imgCols, imgRows, BADPIX=badpixsinterp, $
			MISSING=missinginterp)
	ENDIF
	IF doverbose THEN print, 'Gaps filled'
ENDIF

;; Remove spikes from cosmic ray hits
Despike:
curStep=12
IF modules.despike THEN BEGIN
	IF rad THEN BEGIN
		radOut[imgCols[0],imgRows[0]] = DIcal_despike($
			radOut[imgCols[0]:imgCols[1],imgRows[0]:imgRows[1]],$
			radHdr,$
			FLAGS=tmp,$
			NSIGMA=spikesig, $
			ITERATE=spikeIter, $
			BOXSIZE=spikebox, $
			MED=spikemed)
		radFlags[imgCols[0],imgRows[0]] = tmp OR $
			radFlags[imgCols[0]:imgCols[1],imgRows[0]:imgRows[1]]
	ENDIF
	IF ioverf THEN BEGIN
		iofOut[imgCols[0],imgRows[0]] = DIcal_despike($
			iofOut[imgCols[0]:imgCols[1],imgRows[0]:imgRows[1]], $
			iofHdr, $
			FLAGS=tmp, $
			NSIGMA=spikesig, $
			ITERATE=spikeIter, $
			BOXSIZE=spikebox, $
			MED=spikemed)
		ifFlags[imgCols[0],imgRows[0]] = tmp OR $
			ifFlags[imgCols[0]:imgCols[1],imgRows[0]:imgRows[1]]
	ENDIF
	IF doverbose THEN print, 'Despiked image'
ENDIF

;; Remove random gaussian noise
Denoise:
curStep = 13
IF modules.denoise THEN BEGIN
	IF rad THEN $
		radOut[imgCols[0],imgRows[0]] = DIcal_denoise($
			radOut[imgCols[0]:imgCols[1], imgRows[0]:imgRows[1]], $
			radHdr)
	IF ioverf THEN $
		iofOut[imgCols[0],imgRows[0]] = DIcal_denoise($
			iofOut[imgCols[0]:imgCols[1], imgRows[0]:imgRows[1]], $
			iofHdr)
	IF doverbose THEN print, 'Removed random gaussian noise'
ENDIF

;; MTF and Scattered light correction
MTF:
curStep = 14
IF modules.mtf THEN BEGIN
	IF rad THEN $
		radOut[imgCols[0]+1,imgRows[0]+1] = DIcal_mtf($
			radOut[imgCols[0]+1:imgCols[1]-1,imgRows[0]+1:imgRows[1]-1],$
			radHdr, $
			flags[imgCols[0]+1:imgCols[1]-1,imgRows[0]+1:imgRows[1]-1], $
			FN=mtffn, ALG=mtfalg, PARAM=mtfparam) 
	IF ioverf THEN $
		iofOut[imgCols[0]+1,imgRows[0]+1] = DIcal_mtf($
			iofOut[imgCols[0]+1:imgCols[1]-1,imgRows[0]+1:imgRows[1]-1],$
			iofHdr, $
			flags[imgCols[0]+1:imgCols[1]-1,imgRows[0]+1:imgRows[1]-1],$
			FN=mtffn, ALG=mtfalg, PARAM=mtfparam)
	IF doverbose THEN print, 'Finished Deconvolution'
ENDIF

;; Perform geometric calibration
Geometric:
curStep = 15
IF modules.geom THEN BEGIN
	IF rad THEN $
		radOut = DIcal_geom(radOut,radHdr,FN=geomfn) 
	IF ioverf THEN $
		iofOut = DIcal_geom(iofOut,iofHdr,FN=geomfn)
	IF doverbose THEN print, 'Applied geometric transformation'
ENDIF

Output:
outCnt = dn+rad+revRad+ioverf+revIoverf
out = fltarr(dim[0],dim[1],outCnt,/nozero)
fitsHdr = strarr(hdrSize,outCnt)
flagOut = bytarr(dim[0],dim[1],outCnt,/nozero)
i=0
IF dn THEN BEGIN
	out[*,*,i] = dnOut
	fitsHdr[*,i] = dnHdr
	flagOut[*,*,i] = flags
	i = i+1
ENDIF
IF rad THEN BEGIN
	out[imgCols[0]:imgCols[1],imgRows[0]:imgRows[1],i] = $
		radOut[imgCols[0]:imgCols[1],imgRows[0]:imgRows[1]]
	fitsHdr[*,i] = radHdr
	flagOut[*,*,i] = radFlags
	i = i+1
ENDIF
IF revRad THEN BEGIN
	out[*,*,i] = revRadOut
	fitsHdr[*,i] = revRadHdr
	flagOut[*,*,i] = flags
	i = i+1
ENDIF
IF ioverf THEN BEGIN
	out[imgCols[0]:imgCols[1],imgRows[0]:imgRows[1],i] = $
		iofOut[imgCols[0]:imgCols[1],imgRows[0]:imgRows[1]]
	fitsHdr[*,i] = iofHdr
	flagOut[*,*,i] = ifFlags
	i = i+1
ENDIF
IF revIoverf THEN BEGIN
	out[*,*,i] = revIofOut
	fitsHdr[*,i] = revIofHdr
	flagOut[*,*,i] = flags
ENDIF
flags = flagOut
RETURN, out
END