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