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

;;-----------------------------------------------------------------------------
;; PURPOSE:
;;	Used in the DI Calibration Pipeline.
;;	Runs a despiking routine to remove cosmic ray hits. 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. 
;;
;; 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 (DEFAULT = 3(mean) 10(median))
;;	ITERATE - Specifies the maximum number of times to iterate the
;;		algorithm. Will stop if pixels stop changing
;;	BOXSIZE - How big a box to do statistics for. Should be an odd number. Default=3
;;	MED - if the statistics used should be median and median deviation instead
;;		of mean and standard deviation.
;;
;; 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

;; Do the despiking
i=0
ischange=1
fimg = float(img)
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
		flags[w] = flags[w] OR FLAG_SPIKE
		img[w] = fimg[w]
	ENDIF
END

;; 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'

RETURN, img
END