PRO cr_reject, input_cube, rd_noise_dn, dark_dn, gain, mult_noise, $ combined_image, combined_noise, combined_npix, $ MASK_CUBE=mask_cube, NOISE_CUBE=noise_cube, $ NSIG=nsig, MEDIAN_LOOP=median_loop, MEAN_LOOP=mean_loop, $ MINIMUM_LOOP=minimum_loop, INIT_MED=init_med, $ INIT_MIN=init_min, INIT_MEAN=init_mean, EXPTIME=exptime,$ BIAS=bias, VERBOSE=verbose, $ TRACKING_SET=tracking_set, DILATION=dilation, DFACTOR=dfactor, $ NOSKYADJUST=noskyadjust,NOCLEARMASK=noclearmask, $ XMEDSKY=xmedsky, RESTORE_SKY=restore_sky, $ SKYVALS=skyvals, NULL_VALUE=null_value, $ INPUT_MASK=input_mask, WEIGHTING=weighting, SKYBOX=skybox ;+ ; NAME: ; CR_REJECT ; ; PURPOSE: ; General, iterative cosmic ray rejection using two or more input images. ; ; EXPLANATION: ; Uses a noise model input by the user, rather than ; determining noise empirically from the images themselves. ; ; The image returned has the combined exposure time of all the input ; images, unless the bias flag is set, in which case the mean is ; returned. This image is computed by summation (or taking mean) ; regardless of loop and initialization options (see below). ; ; CALLING SEQUENCE: ; cr_reject, input_cube, rd_noise_dn, dark_dn, gain, mult_noise, $ ; combined_image, combined_npix, combined_noise ; ; MODIFIED ARGUMENT: ; input_cube - Cube in which each plane is an input image. ; If the noise model is used (rd_noise_dn, dark_dn, ; gain), then input_cube must be in units of DN. ; If an input noise cube is supplied (rd_noise_dn ; <0), then the units of input_cube and noise_cube ; merely need to be consistent. ; ; This array is used as a buffer and its contents ; are not guaranteed on output (although for now, ; weighting=0 with /restore_sky should give you back ; your input unaltered). ; ; INPUT ARGUMENTS: ; rd_noise_dn - Read noise per pixel. Units are DN. ; If negative, then the user supplies an error cube ; via the keyword noise_cube. In the latter case, ; mult_noise still applies, since it is basically a fudge. ; dark_dn - Dark rate in DN per pixel per s. This can be a scalar, ; or it can be a dark image divided by the exposure ; time. ; gain - Electrons per DN. ; mult_noise - Coefficient for multiplicative noise term -- helps ; account for differing PSFs or subpixel image shifts. ; ; INPUT KEYWORDS: ; exptime - If the images have different exposure times, pass ; them in a vector. You can leave this off for ; frames with the same exposure time, but dark counts ; won't be treated correctly. ; verbose - If set, lots of output. ; nsig - Rejection limit in units of pixel-to-pixel noise ; (sigma) on each input image. Can be specified as ; an array, in which case the dimension gives the ; maximum number of iterations to run. (Default = ; [8, 6, 4] ; dilation - With dfactor, provides functionality similar to the ; expansion of the CR with pfactor and radius in STSDAS ; crrej. Dilate gives the size of the border to be ; tested around each initially detected CR pixel. ; E.g., dilate=1 searches a 9 X 9 area centered on the ; original pixel. If dfactor is set, the default is 1. ; dfactor - See dilation. This parameter is equivalent to pfactor ; in STSDAS crrej. The current threshold for rejection ; is multiplied by this factor when doing the search ; with the dilated mask. If dilation is set, the default ; for this parameter is 0.5. ; bias - Set if combining biases (divides through by number ; of images at end, since exposure time is 0). ; tracking_set - Subscripts of pixels to be followed through the ; computation. ; noskyadjust - Sky not to be subtracted before rejection tests. Default ; is to do the subtraction. ; xmedsky - Flag. If set, the sky is computed as a 1-d array ; which is a column-by-column median. This is intended ; for STIS slitless spectra. If sky adjustment is ; disabled, this keyword has no effect. ; input_mask - Mask cube input by the user. Should be byte data ; because it's boolean. 1 means use the pixel, ; and 0 means reject the pixel - these rejections ; are in addition to those done by the CR rejection ; algorithm as such. ; ; The following keywords control how the current guess at a CR-free ; "check image" is recomputed on each iteration: ; ; median_loop - If set, the check image for each iteration is ; the pixel-by-pixel median. THE MEAN IS ; RETURNED in combined_image even if you set ; this option. (Default is mean_loop.) ; minimum_loop - If set, the check image for each iteration is ; the pixel-by-pixel minimum. THE MEAN IS ; RETURNED in combined_image even if you set ; this option. (Default is mean_loop.) ; mean_loop - If set, the check image for each iteration is ; the pixel-by-pixel mean. (Same as the default.) ; noclearmask - By default, the mask of CR flags is reset before ; every iteration, and a pixel that has been ; rejected has a chance to get back in the game ; if the average migrates toward its value. If this ; keyword is set, then any rejected pixel stays ; rejected in subsequent iterations. Note that what ; stsdas.hst_calib.wfpc.crrej does is the same ; as the default. For this procedure, the default ; was NOT to clear the flags, until 20 Oct. 1997. ; restore_sky - Flag. If set, the routine adds the sky back into ; input_cube before returning. Works only if ; weighting=0. ; null_value - Value to be used for output pixels to which no ; input pixels contribute. Default=0 ; weighting - Selects weighting scheme in final image ; combination: ; 0 (default) - Poissonian weighting - co-add ; detected DN from non-CR pixels. (Pixel-by- ; pixel scaling up to total exposure time, ; for pixels in stack where some rejected.) ; Equivalent to exptime weighting of rates. ; 1 or more - Sky and read noise weighting of rates. ; Computed as weighted average of DN rates, ; with total exp time multiplied back in ; afterward. ; ; In all cases, the image is returned as a sum in ; DN with the total exposure time of the image ; stack, and with total sky added back in. ; ; The following keywords allow the initial guess at a CR-free "check ; image" to be of a different kind from the iterative guesses: ; ; init_med - If set, the initial check image is ; the pixel-by-pixel median. (Not permitted if ; input_cube has fewer than 3 planes; default is minimum.) ; init_mean - If set, the initial check image is ; the pixel-by-pixel mean. (Default is minimum.) ; init_min - If set, the initial check image is ; the pixel-by-pixel minimum. (Same as the default.) ; ; OUTPUT ARGUMENTS:: ; combined_image - Mean image with CRs removed. ; combined_npix - Byte (or integer) image of same dimensions as ; combined_image, with each element containing ; the number of non-CR stacked pixels that ; went into the result. ; combined_noise - Noise in combined image according to noise model ; or supplied noise cube. ; ; OUTPUT KEYWORDS: ; mask_cube - CR masks for each input image. 1 means ; good pixel; 0 means CR pixel. ; skyvals - Sky value array. For an image cube with N planes, ; this array is fltarr(N) if the sky is a scalar per ; image plane, and fltarr(XDIM, N), is the XMEDSKY ; is set. ; ; INPUT/OUTPUT KEYWORD: ; noise_cube - Estimated noise in each pixel of input_cube as ; returned (if rd_noise_dn ge 0), or input noise ; per pixel of image_cube (if rd_noise_dn lt 0). ; skybox - X0, X1, Y0, Y1 bounds of image section used ; to compute sky. If supplied by user, this ; region is used. If not supplied, the ; image bounds are returned. This parameter might ; be used (for instance) if the imaging area ; doesn't include the whole chip. ; ; COMMON BLOCKS: none ; ; SIDE EFFECTS: none ; ; METHOD: ; ; COMPARISON WITH STSDAS ; ; Cr_reject emulates the crrej routine in stsdas.hst_calib.wfpc. ; The two routines have been verified to give identical results ; (except for some pixels along the edge of the image) under the ; following conditions: ; ; no sky adjustment ; no dilation of CRs ; initialization of trial image with minimum ; taking mean for each trial image after first (no choice ; in crrej) ; ; Dilation introduces a difference between crrej and this routine ; around the very edge of the image, because the IDL mask ; manipulation routines don't handle the edge the same way as crrej ; does. Away from the edge, crrej and cr_reject are the same with ; respect to dilation. ; ; The main difference between crrej and cr_reject is in the sky ; computation. Cr_reject does a DAOPHOT I sky computation on a ; subset of pixels grabbed from the image, whereas crrej searches ; for a histogram mode. ; ; REMARKS ON USAGE ; ; The default is that the initial guess at a CR-free image is the ; pixel-by-pixel minimum of all the input images. The pixels ; cut from each component image are the ones more than nsig(0) sigma ; from this minimum image. The next iteration is based on the ; mean of the cleaned-up component images, and the cut is taken ; at nsig(1) sigma. The next iteration is also based on the mean with ; the cut taken at nsig(2) sigma. ; ; The user can specify an arbitrary sequence of sigma cuts, e.g., ; nsig=[6,2] or nsig=[10,9,8,7]. The user can also specify that ; the initial guess is the median (/init_med) rather than the ; minimum, or even the mean. The iterated cleaned_up images after ; the first guess can be computed as the mean or the median ; (/mean_loop or /median_loop). The minimum_loop option is also ; specified, but this is a trivial case, and you wouldn't want ; to use it except perhaps for testing. ; ; The routine takes into account exposure time if you want it to, ; i.e., if the pieces of the CR-split aren't exactly the same. ; For bias frames (exposure time 0), set /bias to return the mean ; rather than the total of the cleaned-up component images. ; ; The crrej pfactor and radius to propagate the detected CRs ; outward from their initial locations have been implemented ; in slightly different form using the IDL DILATE function. ; ; It is possible to end up with output pixels to which no valid ; input pixels contribute. These end up with the value ; NULL_VALUE, and the corresponding pixels of combined_npix are ; also returned as 0. This result can occur when the pixel is ; very noisy across the whole image stack, i.e., if all the ; values are, at any step of the process, far from the stack ; average at that position even after rejecting the real ; outliers. Because pixels are flagged symmetrically N sigma ; above and below the current combined image (see code), all ; the pixels at a given position can end up getting flagged. ; (At least, that's how I think it happens.) ; ; MODIFICATION HISTORY: ; 5 Mar. 1997 - Written. R. S. Hill ; 14 Mar. 1997 - Changed to masking approach to keep better ; statistics and return CR-affected pixels to user. ; Option to track subset of pixels added. ; Dilation of initially detected CRs added. ; Other small changes. RSH ; 17 Mar. 1997 - Arglist and treatment of exposure times fiddled ; to mesh better with stis_cr. RSH ; 25 Mar. 1997 - Fixed bug if dilation finds nothing. RSH ; 4 Apr. 1997 - Changed name to cr_reject. RSH ; 15 Apr. 1997 - Restyled with emacs, nothing else done. RSH ; 18 Jun. 1997 - Input noise cube allowed. RSH ; 19 Jun. 1997 - Multiplicative noise deleted from final error. RSH ; 20 Jun. 1997 - Fixed error in using input noise cube. RSH ; 12 July 1997 - Sky adjustment option. RSH ; 27 Aug. 1997 - Dilation kernel made round, not square, and ; floating-point radius allowed. RSH ; 16 Sep. 1997 - Clearmask added. Intermediate as well as final ; mean is exptime weighted. RSH ; 17 Sep. 1997 - Redundant zeroes around dilation kernel trimmed. RSH ; 1 Oct. 1997 - Bugfix in preceding. RSH ; 21 Oct. 1997 - Clearmask changed to noclearmask. Bug in robust ; array division fixed (misplaced parens). Sky as ; a function of X (option). RSH ; 30 Jan. 1998 - Restore_sky keyword added. RSH ; 5 Feb. 1998 - Quick help corrected and updated. RSH ; 6 Feb. 1998 - Fixed bug in execution sequence for tracking_set ; option. RSH ; 18 Mar. 1998 - Eliminated confusing maxiter spec. Added ; null_value keyword. RSH ; 15 May 1998 - Input_mask keyword. RSH ; 27 May 1998 - Initialization of minimum image corrected. NRC, RSH ; 9 June 1998 - Input mask cube processing corrected. RSH ; 21 Sep. 1998 - Weighting keyword added. RSH ; 7 Oct. 1998 - Fixed bug in input_mask processing (introduced ; in preceding update). Input_mask passed to ; skyadj_cube. RSH ; 5 Mar. 1999 - Force init_min for 2 planes. RSH ; 1 Oct. 1999 - Make sure weighting=1 not given with noise cube. RSH ; 1 Dec. 1999 - Corrections to doc; restore_sky needs weighting=0. RSH ; 17 Mar. 2000 - SKYBOX added. RSH ;- on_error,0 IF n_params(0) LT 6 THEN BEGIN print,'CALLING SEQUENCE: cr_reject, input_cube, rd_noise_dn, $' print,' dark_dn, gain, mult_noise, combined_image, combined_noise, $' print,' combined_npix' print,'KEYWORD PARAMETERS: nsig, exptime, bias, verbose,' print,' tracking_set, median_loop, mean_loop, minimum_loop, ' print,' init_med, init_mean, init_min,' print,' mask_cube, noise_cube, dilation, dfactor, noclearmask, ' print,' noskyadjust, xmedsky, restore_sky, skyvals, null_value' print,' input_mask, weighting, skybox' return ENDIF verbose = keyword_set(verbose) xmed = keyword_set(xmedsky) track = n_elements(tracking_set) GT 0 sz = size(input_cube) IF sz[0] NE 3 THEN BEGIN print,'CR_REJECT: Input cube must have 3 dimensions.' return ENDIF IF n_elements(input_mask) GT 0 THEN BEGIN szinpm = size(input_mask) wsz = where(szinpm[0:3] NE sz[0:3], cwsz) IF cwsz GT 0 THEN BEGIN print,'CR_REJECT: INPUT_MASK must be same size as IMAGE_CUBE.' return ENDIF ELSE BEGIN IF verbose THEN print,'CR_REJECT: Using INPUT_MASK.' ENDELSE use_input_mask = 1b ENDIF ELSE BEGIN use_input_mask = 0b ENDELSE xdim = sz[1] ydim = sz[2] nimg = sz[3] npix = xdim*ydim usemedian = keyword_set(median_loop) usemean = keyword_set(mean_loop) usemin = keyword_set(minimum_loop) IF (usemean + usemedian + usemin) GT 1 THEN BEGIN print,'CR_REJECT: Specify only one of MEDIAN_LOOP, MEAN_LOOP' $ + ', or MINIMUM_LOOP' return ENDIF IF (usemean + usemedian + usemin) EQ 0 THEN BEGIN usemean = 1 ENDIF inimed = keyword_set(init_med) inimean = keyword_set(init_mean) inimin = keyword_set(init_min) IF (inimean + inimed + inimin) GT 1 THEN BEGIN print,'CR_REJECT: Specify only one of INIT_MED, INIT_MEAN,' $ + ' or INIT_MIN.' return ENDIF IF (inimean + inimed + inimin) EQ 0 THEN BEGIN inimin = 1 ENDIF IF nimg LT 3 AND inimed THEN BEGIN inimed = 0 inimin = 1 IF verbose THEN BEGIN print,'CR_REJECT: INIT_MED only permitted for 3 or more ' $ + 'images.' print,' Forcing INIT_MIN.' ENDIF ENDIF ; ; Accumulation mode for bad pixels. ; IF keyword_set(noclearmask) THEN BEGIN clearmask = 0b IF verbose THEN print,'CR_REJECT: CR flags accumulate strictly.' ENDIF ELSE BEGIN clearmask = 1b IF verbose THEN print,'CR_REJECT: CR flags cleared between iterations.' ENDELSE ; ; Default iterations. ; IF (n_elements(nsig) LT 1) THEN BEGIN nsig = [8, 6, 4] ENDIF sig_limit = nsig maxiter = n_elements(nsig) IF n_elements(null_value) LT 1 THEN null_value=0 IF verbose THEN BEGIN print,'CR_REJECT: Iteration spec: ' print,' nsig = ',nsig print,' maxiter = ',maxiter print,' null_value = ',null_value ENDIF ; IF n_elements(exptime) NE 0 THEN BEGIN IF n_elements(exptime) NE nimg THEN BEGIN print,'CR_REJECT: EXPTIME must have one element per input image.' return ENDIF zexp = 0b FOR i=0,nimg-1 DO zexp = zexp OR (exptime[i] LE 0.0) IF zexp THEN BEGIN save_expt = exptime exptime = make_array(nimg, value=1.0) IF verbose THEN print, $ 'CR_REJECT: All exposure times <= 0.' ENDIF ENDIF ELSE BEGIN zexp = 1b save_expt = make_array(nimg, value=0.0) exptime = make_array(nimg, value=1.0) ENDELSE etot = total(exptime) IF n_elements(weighting) GT 0 THEN BEGIN wgt = weighting wgt = round(wgt) IF wgt ne 0 and wgt ne 1 THEN BEGIN print, 'CR_REJECT: Weighting must be 0 or 1' print,' Executing return' return ENDIF ENDIF ELSE BEGIN wgt = 0 ENDELSE IF verbose THEN BEGIN print,'CR_REJECT: gain = ',gain IF n_elements(dark_dn) EQ 1 THEN BEGIN print,' dark rate = ',dark_dn ENDIF ELSE BEGIN print,' dark image supplied ' ENDELSE print,' read noise = ',rd_noise_dn print,' multiplicative noise coefficient = ',mult_noise print,' number of images = ',nimg print,' exposure times: ' print,exptime print,' total exposure time = ',etot CASE wgt OF 0: print,' flux to be co-added' 1: print,' weighting of rate by sky and read noise' ENDCASE ENDIF ; ; Process dilation specs ; IF keyword_set(dilation) OR keyword_set(dfactor) THEN BEGIN do_dilation = 1b IF n_elements(dilation) LT 1 THEN dilation = 1 IF n_elements(dfactor) LT 1 THEN dfactor = 0.5 IF (dilation LE 0) OR (dfactor LE 0) THEN BEGIN print,'CR_REJECT: Dilation specs not valid: ' print,' dilation = ',dilation print,' dfactor = ',dfactor return ENDIF kdim = 1 + 2*floor(dilation+1.e-4) kernel = make_array(kdim, kdim, value=1b) half_kern = fix(kdim/2) wkz = where(shift(dist(kdim),half_kern,half_kern) $ GT (dilation+0.0001), ckz) IF ckz GT 0 THEN kernel[wkz] = 0b IF verbose THEN BEGIN print,'CR_REJECT: Dilation will be done. Specs:' print,' dilation = ',dilation print,' dfactor = ',dfactor print,' kernel = ' print,kernel ENDIF ENDIF ELSE BEGIN do_dilation = 0b IF verbose THEN print,'CR_REJECT: Mask dilation will not be done.' ENDELSE IF verbose THEN print,'CR_REJECT: Initializing noise and mask cube.' IF rd_noise_dn GE 0 THEN BEGIN IF verbose THEN print,'CR_REJECT: Noise cube computed.' supplied = 0b noise_cube = 0.0*input_cube FOR i=0, nimg-1 DO BEGIN noise_cube[0,0,i] = sqrt((rd_noise_dn^2 $ + ((input_cube[*,*,i] $ + dark_dn*exptime[i])>0)/gain) > 0.0) ENDFOR ENDIF ELSE BEGIN IF verbose THEN print,'CR_REJECT: Noise cube supplied.' supplied = 1b IF wgt EQ 1 THEN BEGIN print, 'CR_REJECT: WEIGHTING=1 incompatible with supplying ', $ 'noise cube.' print, ' Executing return.' return ENDIF ENDELSE ; ; Mask flags CR with zeroes ; mask_cube = make_array(xdim, ydim, nimg, value=1B) IF nimg LE 255 THEN ivalue=byte(nimg) ELSE ivalue=fix(nimg) combined_npix = make_array(xdim, ydim, value=ivalue) IF keyword_set(noskyadjust) THEN BEGIN skyvals = fltarr(nimg) totsky = 0 ENDIF ELSE BEGIN IF verbose THEN print,'CR_REJECT: Sky adjustment being made.' skyadj_cube, input_cube, skyvals, totsky, $ verbose=verbose, xmedsky=xmed, input_mask=input_mask, $ region=skybox ENDELSE IF verbose THEN print,'CR_REJECT: Scaling by exposure time.' FOR i=0,nimg-1 DO BEGIN input_cube[0,0,i] = input_cube[*,*,i]/exptime[i] noise_cube[0,0,i] = noise_cube[*,*,i]/exptime[i] ENDFOR ; ; Initialization of main loop. ; ncut_tot = lonarr(nimg) cr_subs = lonarr(npix) IF inimin OR usemin THEN flagval = max(input_cube)+1 IF inimed THEN BEGIN IF verbose THEN print,'CR_REJECT: Initializing with median.' IF use_input_mask THEN BEGIN medarr,input_cube,combined_image,input_mask ENDIF ELSE BEGIN medarr,input_cube,combined_image ENDELSE ENDIF ELSE IF inimean THEN BEGIN IF verbose THEN print,'CR_REJECT: Initializing with mean.' IF use_input_mask THEN BEGIN tm = total(input_mask,3) > 1e-6 combined_image = total(input_cube*input_mask,3)/tm wz = where(temporary(tm) le 0.001, cwz) IF cwz GT 0 THEN $ combined_image[temporary(wz)] = 0 ENDIF ELSE BEGIN combined_image = total(input_cube,3)/nimg ENDELSE ENDIF ELSE IF inimin THEN BEGIN IF verbose THEN print,'CR_REJECT: Initializing with minimum.' IF use_input_mask THEN BEGIN combined_image = make_array(xdim,ydim,value=flagval) FOR i=0, nimg-1 DO BEGIN indx = where(input_mask[*,*,i] gt 0, cindx) IF cindx GT 0 THEN $ combined_image[indx] = $ (combined_image < input_cube[*,*,i])[indx] ENDFOR wf = where(combined_image EQ flagval, cf) IF cf GT 0 THEN combined_image[wf] = null_value ENDIF ELSE BEGIN combined_image = input_cube[*,*,0] FOR i=1, nimg-1 DO BEGIN combined_image = (combined_image < input_cube[*,*,i]) ENDFOR ENDELSE ENDIF ELSE BEGIN print,'CR_REJECT: Logic error in program initializing check image.' return ENDELSE ; ; ---------------- MAIN CR REJECTION LOOP. ------------------ ; iter=0 main_loop: iter=iter+1 IF clearmask THEN mask_cube[*]=1b IF track THEN BEGIN print,'CR_REJECT: Tracking. Iter = ',strtrim(iter,2) print,' Combined_image: ' print,combined_image[tracking_set] FOR i = 0, nimg-1 DO BEGIN print,' Image ', strtrim(i,2), ':' print,(input_cube[*,*,i])[tracking_set] print,' Noise ', strtrim(i,2), ':' print,(noise_cube[*,*,i])[tracking_set] print,' Mask ', strtrim(i,2), ':' print,(mask_cube[*,*,i])[tracking_set] ENDFOR ENDIF IF verbose THEN BEGIN print,'CR_REJECT: Beginning iteration number ',strtrim(iter,2) print,' Sigma limit = ',sig_limit[iter-1] ENDIF FOR i=0, nimg-1 DO BEGIN skyarray = fltarr(xdim, ydim) IF xmed THEN BEGIN FOR jl = 0,ydim-1 DO skyarray[0,jl] = skyvals[*,i] ENDIF ELSE BEGIN skyarray[*] = skyvals[i] ENDELSE model_image = $ (temporary(skyarray) + (combined_image + dark_dn)*exptime[i])>0 IF supplied THEN BEGIN current_var = noise_cube[*,*,i]^2 $ + ((mult_noise*temporary(model_image))/exptime[i])^2 ENDIF ELSE BEGIN current_var = (rd_noise_dn^2 + model_image/gain $ + (mult_noise*temporary(model_image))^2) $ / (exptime[i]^2) ENDELSE IF track THEN BEGIN print,'CR_REJECT: Tracking. Iter = ',strtrim(iter,2), $ ' Image = ',strtrim(i,2) print,' Current_var: ' print,current_var[tracking_set] ENDIF testnoise = sig_limit[iter-1] * sqrt(temporary(current_var)) IF track THEN BEGIN print,' Testnoise: ' print,testnoise[tracking_set] ENDIF ; ; Absolute value used so that if you remove too much, at least you ; won't introduce a new bias. ; cr_subs[0] = $ where(abs(input_cube[*,*,i] - combined_image) $ GT testnoise, count) IF count GT 0 THEN BEGIN mask_cube[i*npix + cr_subs[0:count-1]] $ = replicate(0b,count) ENDIF IF verbose THEN print,'CR_REJECT: ',strtrim(count,2), $ ' pixels flagged in image ',strtrim(i,2) ; ; Dilation of mask ; count2 = 0 IF do_dilation THEN BEGIN tempw = where(dilate(1b-mask_cube[*,*,i], kernel),dct) IF dct GT 0 THEN BEGIN ic1 = input_cube[npix*i + tempw] tn1 = testnoise[tempw] cmi = combined_image[tempw] tewsub = where(abs(temporary(ic1) $ - temporary(cmi)) $ GT (dfactor*temporary(tn1)), count2) cr_subs[0] = (temporary(tempw))[temporary(tewsub)>0] IF count2 GT 0 THEN BEGIN mask_cube[i*npix + cr_subs[0:count2-1]] $ = replicate(0b,count2) ENDIF ENDIF IF verbose THEN print,'CR_REJECT: Mask dilation performed. ', $ strtrim(count2,2), ' pixels flagged in image ',strtrim(i,2) ENDIF ENDFOR FOR i=0, nimg-1 DO BEGIN cr_subs[0] = where(1b-mask_cube[*,*,i],count) ; IF verbose THEN print,'CR_REJECT: ',strtrim(count,2), $ ; ' accumulated flags in image ',strtrim(i,2) ; IF count GT 0 THEN BEGIN ; input_cube(i*npix + cr_subs(0:count-1)) $ ; = combined_image(cr_subs(0:count-1)) ; noise_cube(i*npix + cr_subs(0:count-1)) $ ; = sqrt(current_var(cr_subs(0:count-1))) ; ENDIF ENDFOR IF use_input_mask THEN BEGIN combined_npix[0,0] = total((mask_cube AND input_mask),3) ENDIF ELSE BEGIN combined_npix[0,0] = total(mask_cube,3) ENDELSE ; ; Loop termination condition. ; IF (iter GE maxiter) THEN GOTO,end_main_loop IF usemedian THEN BEGIN IF verbose THEN print,'CR_REJECT: Taking median.' IF use_input_mask THEN BEGIN medarr,input_cube,combined_image,mask_cube AND input_mask ENDIF ELSE BEGIN medarr,input_cube,combined_image,mask_cube ENDELSE ENDIF ELSE IF usemean THEN BEGIN IF verbose THEN print,'CR_REJECT: Taking mean.' IF use_input_mask THEN BEGIN maskprod = input_mask[*,*,0] AND mask_cube[*,*,0] combined_image = input_cube[*,*,0]*maskprod*exptime[0] combined_expt = temporary(maskprod)*exptime[0] IF nimg GT 1 THEN BEGIN FOR i=1,nimg-1 DO BEGIN maskprod = input_mask[*,*,i] AND mask_cube[*,*,i] combined_image = combined_image $ + input_cube[*,*,i]*maskprod*exptime[i] combined_expt = combined_expt $ + temporary(maskprod)*exptime[i] ENDFOR ENDIF wexpt0 = where(combined_expt LE 0,cexpt0) combined_image = combined_image / (combined_expt>1e-6) IF cexpt0 GT 0 THEN combined_image[wexpt0] = 0 ENDIF ELSE BEGIN combined_image = input_cube[*,*,0]*mask_cube[*,*,0]*exptime[0] combined_expt = mask_cube[*,*,0]*exptime[0] IF nimg GT 1 THEN BEGIN FOR i=1,nimg-1 DO BEGIN combined_image = combined_image $ + input_cube[*,*,i]*mask_cube[*,*,i]*exptime[i] combined_expt = combined_expt $ + mask_cube[*,*,i]*exptime[i] ENDFOR ENDIF wexpt0 = where(combined_expt LE 0,cexpt0) combined_image = combined_image / (combined_expt>1e-6) IF cexpt0 GT 0 THEN combined_image[wexpt0] = 0 ENDELSE ENDIF ELSE IF usemin THEN BEGIN IF verbose THEN print,'CR_REJECT: Taking minimum.' IF use_input_mask THEN BEGIN combined_image[*] = flagval FOR i=0, nimg-1 DO BEGIN indx = where((input_mask[*,*,i] $ AND mask_cube[*,*,i]) gt 0, cindx) IF cindx GT 0 THEN $ combined_image[indx] = $ (combined_image < input_cube[*,*,i])[indx] ENDFOR wf = where(combined_image EQ flagval, cf) IF cf GT 0 THEN combined_image[wf] = null_value ENDIF ELSE BEGIN combined_image = input_cube[*,*,0] FOR i=1, nimg-1 DO BEGIN combined_image = (combined_image < input_cube[*,*,i]) ENDFOR ENDELSE IF use_input_mask THEN BEGIn combined_image = input_cube[*,*,0]*input_mask[*,*,0] FOR i=1, nimg-1 DO BEGIN combined_image = (combined_image < input_cube[*,*,i] $ *input_mask[*,*,i]) ENDFOR ENDIF ELSE BEGIN combined_image = input_cube[*,*,0] FOR i=1, nimg-1 DO BEGIN combined_image = (combined_image < input_cube[*,*,i]) ENDFOR ENDELSE ENDIF ELSE BEGIN print,'CR_REJECT: Logic error in program recomputing check image.' return ENDELSE GOTO,main_loop END_main_loop: ; ; End of CR rejection loop. ; IF verbose THEN BEGIN FOR i=0,nimg-1 DO BEGIN wdummy = where(1b-mask_cube[*,*,i],count) ncut_tot[i] = count ENDFOR print,'CR_REJECT: Total pixels changed: ' print,ncut_tot ENDIF IF track THEN BEGIN print,'CR_REJECT: Tracking. After loop exit.' print,' Combined_image: ' print,combined_image[tracking_set] ; print,' Current_var: ' ; print,current_var[tracking_set] FOR i = 0, nimg-1 DO BEGIN print,' Image ', strtrim(i,2), ':' print,(input_cube[*,*,i])[tracking_set] print,' Noise ', strtrim(i,2), ':' print,(noise_cube[*,*,i])[tracking_set] print,' Mask ', strtrim(i,2), ':' print,(mask_cube[*,*,i])[tracking_set] ENDFOR ENDIF ; ; Compute weights according to scheme chosen ; xrepl = make_array(dim=xdim,value=1) yrepl = make_array(dim=ydim,value=1) IF wgt EQ 0 THEN BEGIN wgts = xrepl # exptime ENDIF ELSE BEGIN IF xmed THEN skytmp = skyvals>1e-6 ELSE skytmp = xrepl # (skyvals>1e-6) exp2tmp = xrepl # (exptime^2) sky_rate_var = temporary(skytmp)/gain/exp2tmp ron_rate_var = rd_noise_dn^2/temporary(exp2tmp) wgts = 1.0/(temporary(sky_rate_var) + temporary(ron_rate_var)) ENDELSE ; ; Do the final co-addition ; wgt_coeff = fltarr(xdim, ydim) FOR i=0,nimg-1 DO BEGIN plane_wgts = wgts[*,i] # yrepl input_cube[0,0,i] = input_cube[*,*,i]*plane_wgts noise_cube[0,0,i] = noise_cube[*,*,i]*plane_wgts IF use_input_mask THEN BEGIN mcim = (mask_cube[*,*,i] AND input_mask[*,*,i]) ENDIF ELSE BEGIN mcim = mask_cube[*,*,i] ENDELSE wgt_coeff[0,0] = wgt_coeff + temporary(mcim) * temporary(plane_wgts) ENDFOR wh0 = where(combined_npix EQ 0,c0) wgt_coeff = etot/(wgt_coeff > 1.0e-8) IF c0 GT 0 THEN wgt_coeff[wh0] = 0.0 IF verbose THEN BEGIN IF c0 GT 0 THEN $ print,'CR_REJECT: ',strtrim(c0,2),' pixels rejected on all inputs.' ENDIF IF use_input_mask THEN BEGIN IF xmed THEN BEGIN combined_image = wgt_coeff * total(input_cube $ * (mask_cube AND input_mask),3) $ + totsky#yrepl ENDIF ELSE BEGIN combined_image = wgt_coeff * total(input_cube $ * (mask_cube AND input_mask),3) $ + totsky ENDELSE combined_noise = wgt_coeff * sqrt(total((noise_cube $ * (mask_cube AND input_mask))^2,3)) ENDIF ELSE BEGIN IF xmed THEN BEGIN combined_image = wgt_coeff * total(input_cube*mask_cube,3) $ + totsky#yrepl ENDIF ELSE BEGIN combined_image = wgt_coeff * total(input_cube*mask_cube,3) $ + totsky ENDELSE combined_noise = wgt_coeff * sqrt(total((noise_cube*mask_cube)^2,3)) ENDELSE IF keyword_set(bias) THEN BEGIN print,'CR_REJECT: Bias flag set -- returning mean instead of total.' combined_image = combined_image/nimg combined_noise = combined_noise/nimg ENDIF IF c0 GT 0 THEN combined_image[wh0] = null_value IF keyword_set(restore_sky) THEN BEGIN IF wgt EQ 0 THEN BEGIN IF verbose THEN print,'CR_REJECT: Adding sky back into data cube' IF xmed THEN BEGIN FOR i=0,nimg-1 DO BEGIN FOR j=0, ydim-1 DO input_cube[0,j,i] = input_cube[*,j,i] $ + skyvals[*,i] ENDFOR ENDIF ELSE BEGIN FOR i=0,nimg-1 DO $ input_cube[0,0,i] = input_cube[*,*,i] + skyvals[i] ENDELSE ENDIF ELSE BEGIN print, 'CR_REJECT: /RESTORE_SKY ignored because weighting spec ' $ + 'not zero.' ENDELSE ENDIF IF zexp THEN exptime = save_expt return END