FUNCTION DIcal_IR, imgIn, fitsHdr, modules, VERBOSE=verbose, SILENT=silent, $
	FORCE=force, MAXGAPSIZE=maxgapsize, DN=dn, RAD=rad, REVRAD=revrad, $
	IOVERF=ioverf, REVIOVERF=revioverf, DARKFN=darkfn, $
	USEDARKMODEL=usedarkmodel, COMPRESSFN=compressfn, ADCFN=adcfn, $
	BADPIXSFN=badpixsfn, LINDNFN=lindnfn, FLATFN=flatfn, MTFFN=mtffn, $
	SPECFN=specfn, CONSTMAPFN=constmapfn, SPECMAP=specmap, $
	BADPIXSINTERP=badpixsinterp, MISSINGINTERP=missinginterp, $
	FLAGS=flags, GAINFN=gainfn, HLUTFN=hlutfn, SPIKESIG=spikesig, $
	SPIKEITER=spikeIter,  SPIKEBOX=spikebox, SPIKEMED=spikemed, $
	COMPRESSMETH=compressmeth, MTFALG=mtfalg, MTFPARAM=mtfparam, SNR=snr

;;-----------------------------------------------------------------------------
;; PURPOSE:
;;	Performs the calibration for a HRI IR Deep Impact image.
;;	The following steps are performed:
;;	1. Temps and Voltages in the header are calibrated
;;	2. Decompresses the image if appropriate
;;	3. Flags the saturated pixels
;;	4. Corrects uneven ADC bit weighting
;;	5. Linearize DN values
;;	6. Removes cross talk
;;	7. Subtracts the dark frame
;;	8. Normalizes quadrant gains
;;	9. Applies flat field
;;	10. Flags bad pixels
;;	11. Convert from DN to radiance and/or I/F units
;;	12. Fill in small gaps
;;	13. Remove cosmic ray hits by despiking
;;	14. Remove random gaussian noise
;;	15. MTF and Scattered light correction
;;
;;
;; CALLING SEQUENCE:
;;	RESULT = DIcal_IR(in, fitsHdr, /VERBOSE, /SILENT, /FORCE, $
;;	MAXGAPSIZE=maxgapsize, /RAD, /IOVERF, DARKFN=darkfn, /USEDARKMODEL, $
;;	COMPRESSFN=compressfn, ADCFN=adcfn, BADPIXSFN=badpixsfn, $
;;	FLATFN=flatfn, LINDNFN=lindnfn, MTFFN=mtffn, SPECFN=specfn, $
;;	CONSTMAPFN=constmapfn
;;
;; 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
;;	DARKFN - file name of the FITS image which is a dark frame to 
;;		subtract from the image. If this is not set, a default one 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
;;	FLATFN - filename of the flat field to use
;;	LINDNFN - filename of the linearization polynomial to use
;;	MTFFN - filename of the convolution filter to use
;;	SPECFN - filename of the spectral map to use
;;	CONSTMAPFN - filename of the ST calibration constant map to use
;;	SPECMAP - variable to hold the spectral map for this image
;;	BADPIXSINTERP - if the bad pixels are interpolated over
;;	MISSINGINTERP - if missing data should be interpolated over
;;	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
;;	HLUTFN - filename of the H lambda lookup table 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
;;	COMPRESSMETH - The method used to decompress the image
;;	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 each pixel
;;
;; EXAMPLE:
;;      IDL> imgOut = DIcal_IR(imgIn, fitsHdr, /RAD, /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_linearDN - linearizes the DN values
;;	DIcal_darkFrame - subtracts a dark frame
;;	DIcalIR_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_geom - performs geometric calibration
;;	DIcal_mtf - performs MTF and scattered light corrections
;;	DIcal_RadToIF - convertrs from radiance to I/F units
;;	DIcal_flatfield - applies a flat field
;;	DIcal_gain - normalizes quadrant gain
;;
;; MODIFICATION HISTORY:
;;   2004-06-24  M. Desnoyer    Created
;;
;;-----------------------------------------------------------------------------

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

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
	CASE curStep OF
	0: GOTO, Decompress
	1: GOTO, FlagSat
	2: GOTO, ADCBit
	3: GOTO, linDN
	4: GOTO, XTalk
	5: GOTO, Darks
	6: GOTO, Gain
	7: GOTO, Flat
	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, Output
	ENDCASE
ENDIF

verbose = keyword_set(verbose)

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

;; Make sure a valid output format was chosen
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

;; Convert to floating point numbers
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])

;; Get the smoothed bench temperature
fitsHdr = dical_smootemp(fitsHdr)

;; Calibrate the temperatures and voltages in the header
curStep=0
IF modules.tempvolt THEN BEGIN
	fitsHdr = dical_tempVolt(fitsHdr)
	IF keyword_set(verbose) THEN print, 'Temps and Voltages Calibrated'
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 verbose 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 verbose 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 verbose THEN print, 'Corrected uneven bit weighting'
ENDIF

;; Linearize DN values
linDN:
curStep = 4
IF modules.linDN THEN BEGIN
	out = DIcal_linearDN(out,fitsHdr, FN=lindnfn) 
	IF verbose THEN print, 'Linearized DN values'
ENDIF

;; Remove the crosstalk
XTalk:
curStep=5
IF modules.xTalk THEN BEGIN
	out = DIcal_xtalk(out,fitsHdr, FLAGS=flags)
	IF verbose THEN print, 'Removed cross-talk'
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 = 6
IF modules.dark THEN BEGIN
	out = DIcal_darkFrame(out,fitsHdr, flags, DARKFN=darkfn, USEFIT=usedarkmodel) 
	IF verbose THEN print, 'Subtracted dark frame'
ENDIF

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

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

;; Perform flat fielding
Flat:
curStep = 8
IF modules.flat THEN BEGIN
	out = DIcal_flatfield(out, fitsHdr, FN=flatfn)
	IF verbose THEN print, 'Applied flat-field normalization'
ENDIF

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

;; Convert to I/F units
ChangeUnits:
curStep = 10
dnOut = out
dnHdr = fitsHdr
radOut = dicalir_dntorad(out,fitsHdr, CONSTFN=constmapfn, SPECFN=specfn, $
	SPECMAP=specmap)
radHdr = fitsHdr
revRadOut = radOut
revRadHdr = radHdr
IF ioverf OR revIoverf THEN BEGIN
	iofOut = dical_radtoif(radOut,fitsHdr, SPECFN=specfn, HLUTFN=hlutfn)
	iofHdr = fitsHdr
	revIofOut = iofOut
	revIofHdr = iofHdr
ENDIF
IF verbose THEN print, 'Converted from DN to physical units'

;; Isolate the flag maps for the rest of the pipeline
radFlags = flags
ifFlags = flags

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

;; Remove 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 verbose 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 verbose THEN print, 'Random noise removed'
ENDIF

;; MTF and Scattered light correction. Currently doesn't exist for IR, so force a skip
MTF:
curStep = 14
IF modules.mtf AND 0 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],$
			flags[imgCols[0]+1:imgCols[1]-1,imgRows[0]+1:imgRows[1]-1], $
			radHdr, 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],$
			flags[imgCols[0]+1:imgCols[1]-1,imgRows[0]+1:imgRows[1]-1],$
			iofHdr, FN=mtffn, ALG=mtfalg, PARAM=mtfparam)
	IF verbose THEN print, 'Finished MTF calibration'
ENDIF

Output:
outCnt = dn+rad+revRad+ioverf+revIoverf
out = fltarr(dim[0],dim[1],outCnt)
fitsHdr = strarr((size(fitsHdr, /dimensions))[0],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