cr_reject.pro

basic median filter simulation

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