FUNCTION pocMean, in, flags, N
;;-----------------------------------------------------------------------------
;; PURPOSE:
;;	Calculates the resistant mean of POC rows. 
;;
;; CALLING SEQUENCE:
;;	mean = pocMean(img, flags, N)
;;
;; REQUIRED INPUTS:
;;	in - The values to take the mean of in a column by column fashion
;;	flags - flags specifying where the good and bad data are
;;	N - Number of standard deviations beyond which to drop data
;;
;; OUTPUTS:
;;	RETURN - a row of resistant means. If data is invalid, returns a 0
;;
;; OPTIONAL INPUT KEYWORDS:
;;
;; EXAMPLE:
;;      IDL> mean = pocMean(img, flags, N)
;;
;; PROCEDURES USED (i.e. called directly!):
;;
;; MODIFICATION HISTORY:
;;   2005-05-10  M. Desnoyer    Created
;;
;;-----------------------------------------------------------------------------

img = double(in)

;; Figure out the size of the data
dim = size(img, /dimensions)

;; Replace all invalid data with 0s
badLoc = where(flags NE 0,badCnt)
IF badCnt GT 0 THEN img[badLoc] = 0

;; Figure out how many valid values in each column
valCnt = total(flags EQ 0,2)

;; Calculate the median and absolute median deviations
med = dblarr(dim[0],dim[1])
medDev = med
FOR i=0, dim[0]-1 DO BEGIN
	IF valCnt[i] GT 0 THEN BEGIN
		w = where(flags[i,*] EQ 0)
		med[i,*] = median(img[i,w],/even)
		absDev = abs(img[i,w]-med[i,0])
		IF max(absDev) EQ 0 THEN absDev[*] = 1
		medDev[i,*] = median(absDev,/even)/0.6745
		IF medDev[i,0] LT 1.0E-24 THEN medDev[i,*] = avg(absDev)/.8
	ENDIF ELSE BEGIN
		med[i,*] = 0
		medDev[i,*] = 0
	ENDELSE
END

;; Throw out median outliers and update teh valid count in each column
outlierLoc = where(abs(img-med) GT abs(N*medDev) AND flags EQ 0,outlierCnt)
IF outlierCnt GT 0 THEN BEGIN
	outlierMap = bytarr(dim[0],dim[1])
	outlierMap[outlierLoc] = 1
	img[outlierLoc] = 0
	valCnt = valCnt - total(outlierMap,2)
ENDIF
		

;; Calculate mean and standard deviations
mn = rebin(total(img,2)/valCnt,dim[0],dim[1])
dev = sqrt(total((img - mn)^2,2)/(valCnt-1))
w = where(dev EQ 0,zcnt)
IF zcnt GT 0 THEN dev[w] = 1
dev = rebin(dev,dim[0],dim[1])

;; Throw out outliers and update the valid count in each column
outlierLoc = where(abs(img-mn) GE abs(N*dev) AND flags EQ 0,outlierCnt)
IF outlierCnt GT 0 THEN BEGIN
	outlierMap = bytarr(dim[0],dim[1])
	outlierMap[outlierLoc] = 1
	img[outlierLoc] = 0
	valCnt = valCnt - total(outlierMap,2)
ENDIF

;; Calculate the final resistant mean with zeros for bad means
out = dblarr(dim[0])
w = where(valCnt GT 0,cnt, complement=badLoc, ncomplement=badCnt)
IF badCnt GT 0 THEN out[badLoc] = 0
IF cnt GT 0 THEN out[w] = total(img[w,*],2)/valCnt[w]

;; Clear any math errors that might have occured
tmp = check_math()

RETURN, out
END