PRO CDS_CLEAN, DATA, FILL_MISSING=FILL_MISSING, MINSIG=MINSIG, $ ERRMSG=ERRMSG, NOFILL=NOFILL ;+ ; Project : SOHO - CDS ; ; Name : CDS_CLEAN ; ; Purpose : Removes cosmic rays from CDS data ; ; Category : Class3, Analysis ; ; Explanation : This routine removes cosmic rays and missing pixels from NIS ; data structures. ; ; Cosmic rays are recognized as pixels which are very high ; compared to neighboring exposures. Thus, this procedure can ; only be used on data which contain three or more exposures. ; Examples include images built up from rastering the slit ; across the sun, or from a "sit-and-stare" type program. ; ; The algorithm used to recognize cosmic rays is far from ; perfect. Use it at your own risk. ; ; Syntax : CDS_CLEAN, DATA ; ; Examples : A = READCDSFITS('s3300r00.fits') ; CDS_CLEAN, A ; ; Inputs : DATA = CDS data structure to be cleaned. ; ; Opt. Inputs : None. ; ; Outputs : The input data structure is returned with the cosmic rays ; removed. ; ; Opt. Outputs: None. ; ; Keywords : FILL_MISSING = If set, then the routine will attempt to fill in ; missing pixels as well as removing cosmic rays. ; ; NOFILL = If set, then cosmic rays are only marked as missing, ; and are not filled in. When /NOFILL is used, then ; the value of FILL_MISSING is ignored. ; ; MINSIG = Minimum value to use when estimating the standard ; deviation in the pixels. ; ; Calls : AVERAGE, SIG_ARRAY, GET_IM_KEYWORD ; ; Common : None. ; ; Restrictions: None. ; ; Side effects: Computationally intensive. ; ; Prev. Hist. : Based on CDS_CLEAN_IMAGE. ; ; History : Version 1, 09-Aug-1996, William Thompson, GSFC ; Version 2, 12-Aug-1996, William Thompson, GSFC ; Only average vertically. ; Version 3, 17-Sep-1996, William Thompson, GSFC ; Fix bug, use ST_WINDATA ; Version 4, 25-Nov-1996, William Thompson, GSFC ; Added keyword NOFILL ; ; Contact : WTHOMPSON ;- ; ; ON_ERROR, 2 ; IF N_PARAMS() NE 1 THEN BEGIN MESSAGE = 'Syntax: CDS_CLEAN, DATA' GOTO, HANDLE_ERROR ENDIF ; ; Make sure that it is VDS data we're dealing with. ; IF DATA.HEADER.DETECTOR NE 'NIS' THEN BEGIN MESSAGE = 'Cannot clean GIS data' GOTO, HANDLE_ERROR ENDIF ; ; Determine the number of windows, and step through them. ; NWIN = N_ELEMENTS(DATA.DETDESC) FOR IWIN = 0,NWIN-1 DO $ IF DATA.DETDESC(IWIN).IXSTOP(0) GE 0 THEN BEGIN ; ; Extract the data array from the input structure. Save the dimensions for ; later. ; ARRAY = DETDATA(DATA, IWIN+1) DESC = DETDESC(DATA, IWIN+1) SA = SIZE(ARRAY) ; ; Determine the dimension which refers to multiple exposures. ; IF SA(0) EQ 4 THEN DIMENSION = 4 ELSE DIMENSION = 2 ; ; If either the FILL_MISSING or NOFILL keyword was set, then use the MISSING ; value from the data structure. Otherwise, use the smallest value in the ; array, minus 1. ; IF KEYWORD_SET(FILL_MISSING) OR KEYWORD_SET(NOFILL) THEN $ MISSING = DESC.MISSING ELSE MISSING = MIN(ARRAY) - 1 ; ; Get the number of elements in the requested dimension, and the number of ; elements in all the dimensions both before and after. ; IF DIMENSION GT SA(0) THEN BEGIN MESSAGE = 'ARRAY doesn''t have ' + TRIM(DIMENSION) + $ ' dimensions' GOTO, HANDLE_ERROR ENDIF ND = SA(DIMENSION) IF ND LT 3 THEN BEGIN MESSAGE = 'Not enough exposures to clean data' GOTO, HANDLE_ERROR ENDIF IF DIMENSION GT 1 THEN NB = PRODUCT(SA(1:DIMENSION-1)) ELSE NB = 1 IF DIMENSION LT SA(0) THEN NA = PRODUCT(SA(1+DIMENSION:SA(0))) $ ELSE NA = 1 ; ; Create a temporary working array. ; NB = LONG(NB) NA = LONG(NA) TEMP = REFORM(1.*ARRAY,NB,ND,NA) GOOD = TEMP NE MISSING BAD = TEMP EQ MISSING W_BAD = WHERE(BAD, N_BAD) IF N_BAD GT 0 THEN TEMP(W_BAD) = 0 ; ; Form the boxcar total and total square of the image over three vertical ; pixels at each point. (Note that some wrap-around occurs at the boundaries ; when DIMENSION=4. ; TOT = TEMP ASQ = TEMP^2 NAV = GOOD ; IF DIMENSION EQ 2 THEN BEGIN TOT(*,*,0:NA-2) = TOT(*,*,0:NA-2) + TEMP(*,*,1:*) ASQ(*,*,0:NA-2) = ASQ(*,*,0:NA-2) + TEMP(*,*,1:*)^2 NAV(*,*,0:NA-2) = NAV(*,*,0:NA-2) + GOOD(*,*,1:*) TOT(*,*,1:*) = TOT(*,*,1:*) + TEMP(*,*,0:NA-2) ASQ(*,*,1:*) = ASQ(*,*,1:*) + TEMP(*,*,0:NA-2)^2 NAV(*,*,1:*) = NAV(*,*,1:*) + GOOD(*,*,0:NA-2) END ELSE BEGIN TOT(0:NB-2,*,*) = TOT(0:NB-2,*,*) + TEMP(1:*,*,*) ASQ(0:NB-2,*,*) = ASQ(0:NB-2,*,*) + TEMP(1:*,*,*)^2 NAV(0:NB-2,*,*) = NAV(0:NB-2,*,*) + GOOD(1:*,*,*) TOT(1:*,*,*) = TOT(1:*,*,*) + TEMP(0:NB-2,*,*) ASQ(1:*,*,*) = ASQ(1:*,*,*) + TEMP(0:NB-2,*,*)^2 NAV(1:*,*,*) = NAV(1:*,*,*) + GOOD(0:NB-2,*,*) ENDELSE ; ; Step through each column in the image, and see if it's more than 3 sigma ; above its (wrapped around) neighbors. ; TAV = TOT TAV(*,1:ND-2,*) = (TOT(*,0:ND-3,*) + TOT(*,2:*,*)) / $ ((NAV(*,0:ND-3,*) + NAV(*,2:*,*)) > 1) TAV(*,0,*) = (TOT(*,ND-1,*) + TOT(*,1,*)) / $ ((NAV(*,ND-1,*) + NAV(*,1,*)) > 1) TAV(*,ND-1,*) = (TOT(*,ND-2,*) + TOT(*,0,*)) / $ ((NAV(*,ND-2,*) + NAV(*,0,*)) > 1) ; TSQ = ASQ TSQ(*,1:ND-2,*) = (ASQ(*,0:ND-3,*) + ASQ(*,2:*,*)) / $ ((NAV(*,0:ND-3,*) + NAV(*,2:*,*)) > 1) TSQ(*,0,*) = (ASQ(*,ND-1,*) + ASQ(*,1,*)) / $ ((NAV(*,ND-1,*) + NAV(*,1,*)) > 1) TSQ(*,ND-1,*) = (ASQ(*,ND-2,*) + ASQ(*,0,*)) / $ ((NAV(*,ND-2,*) + NAV(*,0,*)) > 1) ; SIG = SQRT((TSQ - TAV^2) > 0) IF N_ELEMENTS(MINSIG) EQ 1 THEN SIG = SIG > MINSIG W = WHERE((TEMP GT TAV+3*SIG) AND (SIG GT 0), N_CR) IF N_CR GT 0 THEN BAD(W) = 1 ; ; Set all the cosmic rays found to the missing pixel flag value. ; W_BAD = WHERE(BAD, N_BAD) IF N_BAD GT 0 THEN BEGIN TEMP(W_BAD) = MISSING ; ; Fill in the missing pixels. ; IF NOT KEYWORD_SET(NOFILL) THEN FOR IA = 0, NA-1 DO BEGIN FOR IB = 0, NB-1 DO BEGIN TEMP2 = REFORM(TEMP(IB,*,IA)) FILL_MISSING, TEMP2, MISSING TEMP(IB,*,IA) = TEMP2 ENDFOR ENDFOR ENDIF ; ; Put the result back into the input structure. ; ARRAY = REFORM(TEMPORARY(TEMP), SA(1:SA(0))) ST_WINDATA, DATA, IWIN, ARRAY, /NO_COPY ENDIF ;IWIN valid GOTO, FINISH ; ; Error handling point. ; HANDLE_ERROR: IF N_ELEMENTS(ERRMSG) NE 0 THEN ERRMSG = 'CDS_CLEAN: ' + MESSAGE $ ELSE MESSAGE, MESSAGE ; ; Exit point. ; FINISH: RETURN END