FUNCTION DIcal_despike, in, fitsHdr, FLAGS=flags, NSIGMA=nsigma, ITERATE=iterate, $
	BOXSIZE=boxsize, MED=med, ALG=alg

;;-----------------------------------------------------------------------------
;; PURPOSE:
;;	Used in the DI Calibration Pipeline.
;;	Runs a despiking routine to remove cosmic ray hits.  
;;
;; CALLING SEQUENCE:
;;	out = DIcal_despike(in, NSIGMA=nsigma, ITERATE=iterate, BOXSIZE=boxsize, /MED)
;;
;; REQUIRED INPUTS:
;;	in - The image to which this pipeline element is going to be applied
;;	fitsHdr - The FITS header of the image
;;
;; OUTPUTS:
;;	RETURN - the image after going through this calibration step
;;
;; OPTIONAL INPUT KEYWORDS:
;;	FLAGS - map of flags to add the despike flags to
;;	NSIGMA - The number of standard deviations a value must be beyond 
;;		to be flaged as spike in algorithm 1 (DEFAULT = 3(mean) 10(median))
;;	ITERATE - Specifies the maximum number of times to iterate algorithm 1. 
;;		Will stop if pixels stop changing
;;	BOXSIZE - How big a box to do statistics for in algorithm 1. Should be an odd number. Default=3
;;	MED - if the statistics used should be median and median deviation instead
;;		of mean and standard deviation. Used for algorithm 1 only.
;;	ALG - Defines which algorithm to use:
;;	  1 -   Spikes are identified by taking the mean and stdev of the pixels in a nxn block 
;;		around each non-edge pixel. If the pixel is above nsigma*stddev, then it is 
;;		a spike. The spike is replaced by the mean of the surrounding good pixels.
;;		Median statistics can also be used. This is the default choice.
;;	  2 -	Use a version of the imgclean algorithm with default parameters. See imgcleanf.pro for details.
;;		Can only be used for VIS images.
;;
;; EXAMPLE:
;;      IDL> imgOut = DIcal_despike(imgIn)
;;
;; PROCEDURES USED (i.e. called directly!):
;;    Functions:
;;    Procedures:
;;
;; MODIFICATION HISTORY:
;;   2004-05-24  M. Desnoyer    Created
;;   2005-05-26  M. Desnoyer   	Added BOXSIZE and MED
;;
;;-----------------------------------------------------------------------------

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

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

IF size(in, /n_dimensions) NE 2 THEN message, 'Invalid Image', /NONAME

@dical_flags

;; Identify where the image is
img = double(in)

dim = size(img, /dimensions)

;; Set defaults
IF NOT keyword_set(NSIGMA) THEN sigma = keyword_set(MED)?10:3 ELSE sigma=nsigma
IF NOT keyword_set(ITERATE) THEN iter = 1 ELSE iter=iterate
IF NOT keyword_set(FLAGS) THEN flags = bytarr(dim[0],dim[1])
IF NOT keyword_set(BOXSIZE) THEN box = 3 ELSE box = boxsize
IF NOT keyword_set(alg) THEN algorithm = 1 ELSE algorithm=alg

;; If the image is IR, force using algorithm 1
;;inst = sxpar(fitsHdr, 'INSTRUME', count=c1)
;;IF (c1 EQ 0) THEN message, 'Invalid FITS Header'
;;IF inst EQ 'HRIIR' THEN algorithm = 1

;; Figure out where the actual image is
POCrows = getPOCrows(fitsHdr,COUNT = cnt, COMPLEMENT = imgRows)
BIAScols = getBIAScols(fitsHdr, COUNT=cnt, COMPLEMENT = imgCols)

;; Do the despiking
CASE algorithm OF
1: BEGIN
	;; Cut of the OC rows/cols
	imgTmp = img[imgCols[0]:imgCols[1],imgRows[0]:imgRows[1]]
	fimg = float(imgTmp)
	flagTmp = flags[imgCols[0]:imgCols[1],imgRows[0]:imgRows[1]]

	;; The sigma filter algorithm
	i=0
	ischange=1
	WHILE (i LT iter AND ischange) DO BEGIN
		IF keyword_set(MED) THEN $
			tmp = median_sigma_filter(fimg, box, n_sigma=sigma) $
		ELSE $
			tmp = sigma_filter(fimg,box,n_sigma=sigma)
		i = i+1
		w = where(tmp NE fimg, cnt)
		ischange = cnt NE 0 
		fimg=tmp
		IF cnt GT 0 THEN BEGIN
			flagTmp[w] = flagTmp[w] OR FLAG_SPIKE
			imgTmp[w] = fimg[w]
		ENDIF
	END

	;; Put the OC rows/cols back on
	img[imgCols[0],imgRows[0]] = imgTmp
	flags[imgCols[0],imgRows[0]] = flagTmp

	;; Update the FITS header
	fxaddpar, fitsHdr, 'DESPIKE', 'T'
	fxaddpar, fitsHdr, 'DESPIKET', sigma
	fxaddpar, fitsHdr, 'DESPIKEI', i
	fxaddpar, fitsHdr, 'DESPIKEB', box
	fxaddpar, fitsHdr, 'DESPIKEM', keyword_set(MED)?'MEDIAN':'MEAN'
END
2: BEGIN
	;; Sergei's imgcleane algorithm
	fimg = img
	imgcleanf, fimg, img, header=fitsHdr

	;; Update the flag map
	w = where(fimg NE img,cnt)
	IF cnt GT 0 THEN flags[w] = flags[w] OR FLAG_SPIKE

	;; Update the FITS header
	fxaddpar, fitsHdr, 'DESPIKE', 'T'
	fxaddpar, fitsHdr, 'DESPIKEM', 'IMGCLEAN'
END
ENDCASE

RETURN, img
END