FUNCTION DIcal_desmear, in, fitsHdr, pocRows, imgRows, imgCols, flags, $
	ALG=alg

;;-----------------------------------------------------------------------------
;; PURPOSE:
;;	Used in the DI Calibration Pipeline.
;;	Removes the transfer smear from the image. Transfer smear occurs 
;;	because the image is loaded and removed from the CCD line by line.
;;	Thus, a row receives light from unwanted locations as it is transfered
;;	across the CCD. This smearing effect can be removed assuming that the
;;	camera is held steady and thus the light hitting the detector is 
;;	constant. The amount of charge accumulated in the transfer can be
;;	calculated two different ways. The first way is to use the POC lines,
;;	which have only been exposed when transfer occurs. However, the POC
;;	lines are not reliable in modes 7 and 8 and thus another option is
;;	availible. This option assumes that the smear can be approximated 
;;	by the average of the rows on the edge. This assumption holds as long
;;	as there is a dark area near the edge. 
;;
;; CALLING SEQUENCE:
;;	out = DIcal_desmear, in, fitsHdr, ALG=alg, DIVTHRESH=divthresh
;;
;; REQUIRED INPUTS:
;;	in - The image to which this pipeline element is going to be applied
;;	fitsHdr - The FITS header for the image
;;	pocRows - list of rows in the image that contain POC values
;;	imgRows - the bounds on the rows of the frame that are actually part of the image
;;	imgCols - the bounds on the columns of the frame that are actually part of the 
;;		real image. 
;;	flags - flags specifying where the good and bad data are
;;
;; OUTPUTS:
;;	RETURN - the image after going through this calibration step
;;
;; OPTIONAL INPUT KEYWORDS:
;;	alg - The desmear algorithim to apply.
;;		0 - No desmearing occurs
;;		1 - Smear is calculated using the POC values (default for modes
;;			1-6,9)
;;		2 - Smear is calculated using column averages (default for modes
;;			7,8)
;; EXAMPLE:
;;      IDL> out = DIcal_desmear(imgIn, fitsHdr, ALG=1)
;;
;; PROCEDURES USED (i.e. called directly!):
;;	SXPAR - reads values from a FITS header
;;	resistant_mean - Calculates a mean resistant to outliers
;;	dical_denoise - used to smooth image to find edges
;;
;; MODIFICATION HISTORY:
;;   2004-05-24  M. Desnoyer    Created
;;   2005-05-10  M. Desnoyer    Created new resistant mean approach that handles missing
;;				data for POC averaging
;;   2005-06-08  M. Desnoyer    Fixed column averaging routine
;;
;;-----------------------------------------------------------------------------

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

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

;; Includes
@dical_flags

;; Get the integration time
intT = SXPAR(fitsHdr, 'INTTIME', COUNT = c1)

;; Get the image mode
mode = SXPAR(fitsHdr, 'IMGMODE', COUNT = c2)

;; Make sure the FITS header contains the parameters
IF (c1 EQ 0) or (c2 EQ 0) THEN $
	message, 'Invalid FITS header', /NONAME

;; Default which algorithm to use
IF n_elements(alg) EQ 0 THEN BEGIN
	IF (mode EQ 7) OR (mode EQ 8) THEN usealg = 2 ELSE usealg=1
ENDIF ELSE usealg=alg

img = double(in)

;; Get the dimensions of the picture
dim = size(img,/dimensions)

;; Default the smear values to 0 (alg 1)
kTop = dblarr(dim[0])
kBot = kTop
divFact = 1	;; Factor to divide the result by

;; Get the smear constant
CASE usealg OF

;; Smear correction using POC values
1: BEGIN
	pocCnt = (size(pocRows,/DIMENSIONS))[0]	

	;; Make sure there are POC values
	IF pocCnt EQ 0 THEN BEGIN ;; Abort
		usealg=0
		BREAK 
	ENDIF


	;; Calculate the resistant means
	tmean = pocMean(img[*,pocRows[0:pocCnt/2-1]],flags[*,pocRows[0:pocCnt/2-1]],3.0)
	bmean = pocMean(img[*,pocRows[pocCnt/2:*]],flags[*,pocRows[pocCnt/2:*]],3.0)

	;; Calculate the constant to subtract
	kTop = tmean/4.0
	kBot = bmean/4.0
   END

;; No smear correction
0: BREAK

;; Column Averaging
2: BEGIN
	;; This algorithm determines the transfer smear by summing up the dark
	;; subtracted values along each column and weighting that sum by the
	;; transfer smear time. This value is then subtracted column by column.
	;; This approximation works so long as the unseen parts of the frame are 
	;; dark

	;; Sum along the columns in the top and bottom half
	rTop = total(img[*,imgRows[0]+1:dim[1]/2-1],2)
	rBot = total(img[*,dim[1]/2:imgRows[1]-1],2)

	;; Approximate the unknown signal by that around the edges
	edgeRows = imgRows[0] > 2	;; Number of rows to use for approximation
	fillNeed = 512-(dim[1]/2-imgRows[0]-1)	;; How many estimated rows are required to fill full frame
	rTop = rTop + pocMean(img[*,imgRows[0]+1:imgRows[0]+edgeRows+1],$
		flags[*,imgRows[0]+1:imgRows[0]+edgeRows+1],3.0)*fillNeed
	rBot = rBot + pocMean(img[*,imgRows[1]-1-edgeRows:imgRows[1]-1],$
		flags[*,imgRows[1]-1-edgeRows:imgRows[1]-1],3.0)*fillNeed

	;; Determine the transfer time per row
	transT = 32.0/3

	;; Calculate the constant to subtract on the top and bottom
	kTop = (transT * rTop) / (intT*1000 + 512*transT)
	kBot = (transT * rBot) / (intT*1000 + 512*transT)

	;; Set the division factor
	divFact = 1 - transT/(intT*1000)
	
   END

;; Didn't choose a valid algorithm
ELSE: message, 'Invalid algorithm', /NONAME
ENDCASE

;; Generate the smear frame to subtract
smear = [ [rebin([[kTop],[kTop]],dim[0],dim[1]/2)] $
          , [rebin([[kBot],[kBot]],dim[0],dim[1]/2)] ]

;; Subtract the smear frame except for where the bias columns are
img[imgCols[0],0] = (img[imgCols[0]:imgCols[1],*]-smear[imgCols[0]:imgCols[1],*])/divFact

;; Update the FITS header
IF usealg EQ 1 THEN usealg = 'POC ROWS' $
ELSE IF usealg EQ 2 THEN usealg = 'COLUMN AVERAGING' $
ELSE usealg = ''
fxaddpar, fitsHdr, 'SMEAR', (usealg EQ '')?'F':'T'
fxaddpar, fitsHdr, 'SMEARV', usealg

RETURN, img
END