;+ ; NAME: ircal_satmask ; PURPOSE: ; Determine saturated regions in IRCAL data and mask to NaN ; Also handles IRCAM and NIRC2 images. ; ; NOTES: ; we define as saturated: ; 1. All pixels > 28000 counts per read ; 2. Pixels > 25000 counts per read which are contiguous with ; other saturated pixels. ; 3. Optionally, pixels remain saturated for some number of exposures ; after they were first saturated. This is to account for persistent ; afterglow seen in heavily saturated regions of the chip. ; ; ; INPUTS: ; images an array of images ; headers an array of FITS headers ; KEYWORDS: ; after= saturated pixels may remain hot due to trapped charge even ; after the star has been moved away. If after=N, then pixels ; will be considered bad for the N exposures following saturation. ; default is 0. ; /maskonly Don't set saturated pixels to NANs, just mark their position ; in the satmask variable ; OUTPUTS: ; satmask pixel mask indicating which pixels are saturated. ; HISTORY: ; Began 2003-11-20 14:49:07 by Marshall Perrin ;- PRO ircal_satmask,images,headers,after=after,satmask=satmask,$ maskonly=maskonly,threshhold=threshhold,$ instrument=instrument,polarimetry=polarimetry if ~(keyword_set(instrument)) then instrument="IRCAL" if instrument eq "IRCAM" then begin ; Leuschner IRCAM if ~(keyword_set(threshhold)) then threshhold =19000 ; is this the right value?? after=0 endif else if instrument eq "NIRC2" then begin ; Keck II NIRC2 ; if ~(keyword_set(threshhold)) then threshhold =14000 ; was 14000, but that seems way low. ; no, really 14000 is correct. why did I think it was 20? Oh, that was ; PRE-SKY-SUBTRACTION, so this is bogus. if ~(keyword_set(threshhold)) then threshhold=19000 ; that's the real saturation threshhold for the chip, or maybe around ; 21000. ; actually sometimes we need 19 for the crab. after=2 endif else begin ; IRCAL ; was 27500. MDP Sept 22, 2004 ;if ~(keyword_set(threshhold)) then threshhold =27000 ; is this the right value?? ;if ~(keyword_set(threshhold)) then threshhold =39000 ; is this the right value?? ;after=0 if ~(keyword_set(threshhold)) then threshhold =24000 ; is this the right value?? endelse sz = size(images) if sz[0] eq 3 then nimgs=sz[3] else nimgs=1 pixelsperimage = sz[1]*sz[2] ; this routine assumes images scaled to IRCAL readout levels; ; i.e. saturation starts a bit under 30,000 if arg_present(satmask) then satmask=bytarr(sz[1],sz[2],sz[3])+1B message,"Masking pixels above "+strc(threshhold)+" per coadd to NaN",/info for i=0L,nimgs-1 do begin ; find what's saturated for each image. im = images[*,*,i] coadds = sxpar(headers[*,i],"NCOADDS") whigh= where(im gt threshhold*coadds,highcount) if highcount le 0 then continue ; go on to next image if this one's fine. ; expand the saturated region to include anything above 25000 moresat = region_grow(im,whigh,thresh=[threshhold*0.888,99000.]*coadds) if (moresat[0] gt -1) then whigh = cmset_op(whigh,'and',moresat) if ~(keyword_set(maskonly)) then images[whigh+pixelsperimage*i] = !values.f_nan if arg_present(satmask) then satmask[whigh+pixelsperimage*i]=0 ;atv,satmask,/bl ; TODO ; what about saturated pixels where the values are driven ; low due to CDS readout? Essentially what we want to do ; maybe is to say any pixels completely surrounded by saturated ; pixels are themselves saturated? ; ; 2004-07-02: A more efficient algorithm than the below (but ; one which I am not going to implement right now until there's ; some actual demand for it) is to just SEARCH2D all the points ; outside of the saturated region. Then whatever remains ; (i.e. rest = allpixels - saturated - outside) is the low core ; which is surrounded by saturated pixels and should be marked ; saturated. ; ; ; ; Here is the old algorithm from November 2003 which I think ; works well enough. It might fail in the case of a saturated ; binary star however. ; this is doable efficiently using HISTOGRAM and the REVERSE_INDICES arrays ; use top and bottom REVERSE_INDICES for each bin? ; Wait, it's not really HISTOGRAM so much as a TOTAL for each row/column ; that we want. mask = bytarr(sz[1],sz[2]) mask[whigh]=1 ; make the saturated region convex in the Y direction xts = total(mask,2) wheresat = where(xts gt 0,wheresatcount) for xi=0L,wheresatcount-1 do begin where2 = where(mask[wheresat[xi],*] eq 1) mask[wheresat[xi],min(where2):max(where2)] = 1 endfor ; WARNING! Do Not Do the following on polarimetry images! ; It'll just make everything between the left and right polarimetry ; regions in NANs. Oops. if ~(keyword_set(polarimetry)) then begin ; make the saturated region convex in the X direction yts = total(mask,1) wheresat = where(yts gt 0,wheresatcount) for yi=0L,wheresatcount-1 do begin where2 = where(mask[*,wheresat[yi]] eq 1) mask[min(where2):max(where2),wheresat[yi]] = 1 endfor endif ; for polarimetry, need to make R and L sides have same masked pixels ; right now just do this to integer pixel accuracy if (keyword_set(polarimetry)) then begin displacement = polarimetry.displacement dx = round(displacement[0]) dy = round(displacement[1]) sz = size(mask) sidemask = mask sidemask[0:SZ[1]/2,*] = 0 ; leave only RHS sidemask = shift(sidemask,-dx,-dy) ; shift to left sidemask2 = mask sidemask2[SZ[1]/2:*,*] = 0 ; leave only LHS sidemask2 = shift(sidemask2,dx,dy) ; shift to right mask2 = mask OR sidemask OR sidemask2 mask=mask2 endif whigh0 = whigh whigh = where(mask eq 1) ; find out the new saturated region. stackpush,wsat,[whigh+pixelsperimage*i] ; add to stack for later use in "after" clause if arg_present(satmask) then satmask[whigh+pixelsperimage*i]=0 endfor ; here, we let the user specify that pixels once saturated remain bad for ; some number of exposures afterwards. if keyword_set(after) and n_elements(wsat) gt 0 then begin for i=1,after do begin if ~(keyword_set(maskonly)) then images[wsat+pixelsperimage*i]=!values.f_NAN if arg_present(satmask) then satmask[wsat+pixelsperimage*i]=0 endfor endif nsat = long(total(satmask eq 0)) message,/info,"Marked as saturated "+strc(nsat)+" pixels" end