cr_reject.pro

basic median filter simulation

    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