PRO xrt_pixel_grade, index_in, grade, image, ngrade, grade_frac=grade_frac, $ unprocessed=unprocessed, interpret=interpret, qabort=qabort, quiet=quiet ; ========================================================================= ;+ ; PROJECT: ; Solar-B / XRT ; ; NAME: ; XRT_PIXEL_GRADE ; ; CATEGORY: ; Data calibration ; ; PURPOSE: ; Mode 1 [default]: ; Identify "grade" of XRT pixels in an Nx x Ny array; i.e., whether it is ; affected by saturation, saturation bloom/bleed, contamination spots, ; dust/weakness (low QE), or is "hot", and optionally (for binned data) ; the fraction of the pixel so affected. ; Mode 2 [set interpret=1]: ; Interpret a grade array by generating an Nx x Ny x 5 array which ; gives the locations (and optionally, in the binned case when ; grade_frac=1, the fractional contribution to the pixel) of ; saturation, bloom/bleed, contamination spot, dust/weakness, or ; hot pixels . ; ; CALLING SEQUENCE: ; XRT_PIXEL_GRADE, index_in ,grade ,image [,ngrade] [,grade_frac=grade_frac] $ ; [,unprocessed=unprocessed] [,interpret=interpret] [,qabort=qabort] $ ; [,quiet=quiet] ; ; INPUTS: ; INDEX_IN - [Mandatory] index structure ; associated with IMAGE_IN and GRADE ; GRADE - [Mandatory] ; (Mode 1) a constant indicating the saturation level ; (default = 2500 DN) ; OR ; (Mode 1) a 2-dim array, [Nx,Ny] containing ; 0 if pixel is OK ; 1 if pixel is saturated (>2500 DN) ; Overwritten on output ; OR ; (Mode 2) a 2-dim pixel grade array, [Nx,Ny] ; IMAGE - [Mandatory (mode 1), Optional (mode 2)] ; (2-dim float array, [Nx,Ny]) ; The unprocessed (raw) image to be analyzed (unprocessed=1) ; Overwritten on output if interpret=1. ; ; KEYWORDS: ; /GRADE_FRAC - [Optional] (Boolean) Set = 1 if ; info on the fractions of binned pixels of various ; grades is desired (see GRADE below). ; [default = 0] ; /UNPROCESSED - [Optional] (Boolean) Set = 1 if IMAGE input has not ; been processed in any way. ; [default = 0] ; /INTERPRET - [Optional, mode 1; Mandatory, mode 2] (Boolean) Set = 1 ; to interpret an existing GRADE [Nx,Ny] array (mode 2), ; or a newly generated GRADE array (optionally, in mode 1), ; into an output IMAGE [Nx,Ny,5] containing the locations ; of various pixel grades (see GRADE below). Non-zero ; IMAGE(*,*,0) = saturated pixels, IMAGE(*,*,2) = spot ; pixels, etc. ; If GRADE_FRAC=1, the value indicates the fractional ; effect of the grade in a binned pixel (see GRADE). ; [default = 0] ; ; ; OUTPUTS: ; GRADE - [Mandatory (Mode 1)] a 2-D byte array (1x1 binning case) or ; a 2-D double array (other binnings), [Nx,Ny] containing ; in a compact form the pixel "grade", which is the sum of ; -1 if pixel is an over-scan pixel (calibration images only) ; 0 if pixel is OK ; 1 if pixel is saturated (default >2500 DN) ; 2 if pixel is a victim of saturation "bloom/bleed" ; 4 if pixel is in a contamination spot ; 8 if pixel is in a dust speck or weak (low QE) ; 16 if pixel is "hot" ; 32 if pixel is a "growth" of a dust speck (possibly ; due to conatimination adhering to dust specks) ; PLUS optionally, in case of binning (if the /grade_frac ; flag is set), an encoding in fractional form ; of the number of 1x1 pixels affected within the binned ; pixel. Note that this is relevant only for spot, dust/weak ; and hot grades. For each of these grade types, two ; digits to the right of the decimal place are reserved ; for the number of 1x1 pixels within the binned pixel ; which are of the given type. The digits are assigned in ; the same order as listed above, eg, for a 4x4 binned ; pixel which had 1x1 pixels of which 7 were OK, 3 were ; spot, 2 were dust/weak, and 4 were hot, the grade would be ; 2+8+16 = 26, plus 0.030204, or 26.030204 ; This info is potentially useful for thresholding pixels ; by the degree of "damage" - eg, if only one 1x1 pixel ; in an 8x8 pixel is affected by a contamination spot, ; (ie, only 1/64th spotted) the user may wish to treat it as ; "OK" for purposes of their further analysis. ; IMAGE - [Optional (mode 1), Mandatory (mode 2)] Use as an output ; requires that /interpret flag is set. ; An [Nx,Ny,5] bytearr (grade_frac=0 or unset) or double ; array (grade_frac=1) containing saturation, bloom/bleed, ; spot, dust/weak, and hot pixel maps (in order); with ; fractional contribution to the pixel from each type (if ; grade_frac=1). See /interpret. ; NGRADE - [Optional] Integer array of number of pixels affected ; by given grade ; QABORT - [Optional] (Boolean) Indicates that the program ; exited gracefully without completing. (Might be ; useful for calling programs.) ; 0: Program ran to completion. ; 1: Program aborted before completion. ; QUIET - [Optional] (Boolean) Lets program run without printing ; ; ; EXAMPLES: ; ; Mode 1: generate a GRADE array from IMAGE data with ; associated index array INDEX_IN: ; IDL> grade=2500. ; the saturation level in DN; also works if grade does not ; ; exist or if grade is an Boolean array indicating the ; ; position of saturated pixels ; IDL> xrt_pixel_grade, index_in, grade, image, unprocessed=1 ; ; ; Mode 1: Generate GRADE with fractional occupation of grades for IMAGE ; with binned pixels ; IDL> xrt_pixel_grade, index_in, grade, image, grade_frac=1, unprocessed=1 ; ; ; Mode 1: Generate GRADE with fractional occupation of grades for binned ; pixels, and interpreted output array (IMAGE); Note: in this case ; input GRADE array should already be a Boolean saturated pixel map, ; or GRADE = 2500, as in the first example. Input IMAGE, if used ; (ie, with grade = 2500), will be overwritten on output ; IDL> xrt_pixel_grade, index_in, grade, image, grade_frac=1, interpret=1, $ ; unprocessed=1 ; IDL> sat = image(*,*,0) ; location of saturated pixels (0 or 1) ; IDL> bleed = image(*,*,1) ; location of bloom/bleed pixels (0 or 1) ; IDL> spot = image(*,*,2) ; location & fractional contrib. of spot pixels ; IDL> dust = image(*,*,3) ; " " " " dust/weak " " ; IDL> hot = image(*,*,4) ; " " " " hot " " ; IDL> dustgrow = image(*,*,5) ; estimated location of dust growth pixels (0-1) ; ; ; ; Mode 2: Generate interpreted grade array (IMAGE) from input ; GRADE array ; IDL> xrt_pixel_grade, index_in, grade, image, interpret=1 ; ; ; Mode 2: Generate interpreted grade array (IMAGE) including ; fractional occupation from input GRADE array (original ; data array must have been binned) ; IDL> xrt_pixel_grade, index_in, grade, image, interpret=1, grade_frac=1 ; ; ; COMMON BLOCKS: ; None. ; ; NOTES: ; 1) Current dust/weak and hot arrays are provisional; the hot array ; in particular is probably time variable. The current one based on ; the joint hot pixels in two sets of darks (27/07/07 and 02/11/09). ; The dust/weak map is modified from a preliminary analysis ; by Kano-san. They will be updated (and time variability, if any, ; included) as better maps become available. ; ; 2) Contamination spots pose problems mostly for thinner filters ; (Be-thin, C-poly, Ti-poly, Al-poly, Al-mesh and combinations ; thereof) and have virtually no effect for thicker filters. ; The spot locations given here indicate the maximum extent of their ; effect (i.e., in Al-mesh or G-band); thicker filters show ; progressively weaker effects, with some spots having little or ; no effect. Since the contaminant is primarily a long wavelength ; absorber, the spots are often less prominent (fractionally) when ; illuminated by hotter (active region/flare) plasma. ; ; 3) Note that not all "hot" pixels are equally "problematic" and ; their degree of deviation appears to be temperature and exposure ; time related (higher CCD temperatures and longer exposures ; cause more "hot" pixels). It is best to view a pixel ; designated "hot" with caution; depending on the circumstances ; though, it may still be useable. ; 4) Output for calibration size images (2176x2112) has values of -1 ; outside of the 2048x2048 non-overscanned part of the array starting ; at [x=64,y=0] ; 5) The area near dust specks started to accumulate a semi-opaque ; contaminant, roughly starting after the ; May 2012 light leak event. The contaminant is similar, but possibly ; not identical to that producing the contamination spots. It is ; partly removed by periodic CCD bakeouts, growing back between ; them. This code now contains an approximate model for this ; time dependent contaminant. The exact pixels covered by this ; "dust growth" is less certain, so fractions are not computed. ; The dust growth model will likely be improved over time. ; ; ; MODIFICATION HISTORY: ; progver = 'v2009.Jun.14' ; --- (SSaar) Vers. 0 written. progver = 'v2009.Jul.10' ; --- (SSaar) Vers. 1 written. Round-off errors ; & decoding bugs corrected. progver = 'v2009.Sep.12' ; --- (SSaar) Changed dust map to modified Kano-san ; version. progver = 'v2009.Dec.08' ; --- (SSaar) Changed variable name to avoid ; confusion with function of same name. progver = 'v2010.Aug.19' ; --- (SSaar) Rewote parts of the header to clarify ; operation of program. progver = 'v2010.Dec.06' ; --- (SSaar) Fixed processing for 2176x2112 images, ; updated array positioning keywords used, ; added note 4 progver = 'v2011.Aug.15' ; --- (SSaar) 2nd update to array position keywords, ; fixed a bug with cal frame processing progver = 'v2013.Feb.10' ; --- (SSaar) Updated to 2 epoch spot cotamination, ; added qabort, quiet keywords, removed goto, ; added check for prepped input progver = 'v2013.Apr.01' ; --- (SSaar) Changed processing of 2176x2112 images; ; now gives "skirt" pixels a -1 value. Changed ; note (4) to reflect this. Corrected bloom ; determination. progver = 'v2013.Oct.31' ; --- (SSaar) Added computation of dust ; growth pixels and Note 5. ; ;- ; ========================================================================= ;===== Set some keyword Booleans, plus miscellaneous. qabort=0 sat_level=2500. if keyword_set(interpret) then flagi=1 else flagi=0 if keyword_set(grade_frac) then flagf=1 else flagf=0 if keyword_set(unprocessed) then flagu=1 else flagu=0 if keyword_set(quiet) then flagq=1 else flagq=0 if flagu eq 1 then begin sizi=size(image) if sizi(0) lt 2 then begin ; if unprocesed flag set but no input image, if flagq eq 0 then $ print,'I need an unprocessed input image to analyze ... returning.' ; complain & return abort=1 return endif endif if n_elements(image) gt 0 and flagu eq 1 then flagim=1 else flagim=0 path0=getenv('SSW_XRT')+'/idl/util/' ;path0='' ; ***** change to above when ready for submit ;===== Warning messages & eject for improper input; plus mode setting hist0 = get_history(index_in,'XRT_PREP',found=found) ; is data prepped? if found then begin if flagq eq 0 then print,'Data must not be already prepped... returning.' abort=1 ; if so, abort & return return endif sizg=size(grade) if flagi eq 1 then begin ; if interpreting .... ; print,'Warning: input image will be overwritten with output' + $ ; ' interpreted grade array' if sizg(0) lt 2 then begin ; ...and no input grade array... if flagim eq 0 then begin ; plus no input image, if flagq eq 0 then $ print,'I need a GRADE array or unprocessed image to interpret' + $ ' ... returning.' ; complain & return qabort=1 return endif else mode=1 ; ...and input image but no grade array, mode=1 endif else mode=2 ; ...and if input grade array, mode=2 endif else mode=1 ; if not interpreting, mode=1 index=index_in bin=index.chip_sum n1=index.naxis1 n2=index.naxis2 if index.naxis1 gt 2048 then flagc=1 else flagc=0 ; flag cal images pr=index.pos_row pc=index.pos_col if n_elements(bin) gt 1 then begin if flagq eq 0 then print,'Warning... multple images. Im only doing ' + $ 'the first one.' n1=n1(0) n2=n2(0) bin=bin(0) y00=pr(0) x00=pc(0) endif y00=pr(0) x00=pc(0) n1a=n1<2048 n2a=n2<2048 ;y00=index(0).p1col ;x00=index(0).p1row time=anytim(index(0).date_obs) ; time of observation if mode eq 1 then begin ; make a grade array blank=intarr(n1,n2) ; initialize array flagb=0 if sizg(0) lt 2 then begin ; if grade not 2-D if sizg(0) eq 1 then sat_level=grade ; if sat_level externally set grade=image gt sat_level ; find saturated pixels flagb=1 ; set flag to look for bloom pix endif else begin if total(grade ge 4) ne 0 then begin if flagq eq 0 then begin print,'Warning: this grade array has more than saturation ' print,'information; I''m not going to mess with it.' endif qabort=1 return endif if total(fix(grade) eq 1) gt 0 and total(fix(grade) eq 2) eq 0 then $ flagb=1 ; if n_sat>0 but no bloom, set bloom flag ; total(grade ge eq 1) gt 0 ? endelse ;stop ;* print,'flagb:',flagb if flagb eq 1 then begin ; if blooms undefined or not yet done mask_edge=blank+1 mask_edge(*,[0,n2-1])=0 ; mask to kill wrap-around bloom=mask_edge*(shift(grade,0,1) +shift(grade,0,-1)) ; find bloom pixels bloom = bloom*(grade eq 0)<1 ; ensure no double counts grade = grade + bloom*2 ; add to grade endif else bloom = fix(grade) eq 2 nsat=total(grade eq 1) ;stop ; =============== read in spot, dust & hot pixel files , full-frame restgenx, file=path0+'xrt_spotmasters',spot00,tst,tend nep=(size(spot00))[3] iep=-1 for j=0,nep-1 do begin if (time gt anytim(tst(j))) and (time le anytim(tend(j))) then iep=j ; find spot epoch endfor if iep ne -1 then spot = reform(spot00(*,*,iep)) else spot =intarr(2048,2048) ; extract spot array files=['dust_pix1.dat','hot_pix0.dat'] len=[1035,377] for j=0,1 do begin ; read dust,hot arrays dat=intarr(2048,2048) openr,1,path0+files(j) ; note: hot pixels seem worse arr=intarr(2,len(j)) ; in longer exposures readf,1,arr close,1 xx=reform(arr(0,*)) ; read in x & y coords of features yy=reform(arr(1,*)) dat(xx,yy)=1 ; construct integer full-frame image if j eq 0 then dust=dat else hot = dat endfor spotx=spot(x00:x00+n1a*bin-1,y00:y00+n2a*bin-1) ; extract FOV (unbinned size) dustx=dust(x00:x00+n1a*bin-1,y00:y00+n2a*bin-1) ; extract FOV (unbinned size) hotx=hot(x00:x00+n1a*bin-1,y00:y00+n2a*bin-1) ; extract FOV (unbinned size) hotrb=rebin(hotx*1.,n1a,n2a) ; rebin dustrb=rebin(dustx*1.,n1a,n2a) ; rebin spotrb=rebin(spotx*1.,n1a,n2a) ; rebin graderbi = (hotrb gt 0)*16 + (dustrb gt 0)*8 + (spotrb gt 0)*4 ; rebinned ; integer grade for hot, ; dust, & spot tdustgrow='2012-08-10T00:00:00.00Z' ; dust growth starts after eclipse ; season #1 post light leak tdgrow=anytim(tdustgrow) ; start of "dust growth" era tmx=anytim(index.date_obs) ; time of obs if tmx ge tdgrow then begin ; if in "dust growth" era... xrt_dustgrow,index,dustg,dust ; estimate growth dustgx=dustg(x00:x00+n1a*bin-1,y00:y00+n2a*bin-1) ; extract unbinned FOV dustgrb=rebin(dustgx*1.,n1a,n2a) ; rebin graderbi = graderbi + (dustgrb gt 0)*(dustrb eq 0)*32 ; add dust growth endif else dustgrb=0. if flagc eq 1 then begin gradei = blank - 1 ; array of -1s for overscan gradei(64,0) = byte(grade(64:64+2047,0:2047) + graderbi) endif else gradei = byte(grade + graderbi) ; add to spot,bloom ngrade=[nsat,total(bloom),total(spotrb gt 0),total(dustrb gt 0), $ total(hotrb gt 0),total(dustgrb gt 0)] ; # pix in each grade if flagf eq 1 then begin ; if grade fracs wanted if flagc eq 0 then begin ; (and not a cal image), encode grade = gradei + $ ; encode # of 1x1 spotrb*(spotrb lt 1 and spotrb gt 0)/1d2*bin^2 + $ ; spot dustrb*(dustrb lt 1 and dustrb gt 0)/1d4*bin^2 + $ ; dust & hot hotrb*(hotrb lt 1 and hotrb gt 0)/1d6*bin^2 ; per binned pixel endif else begin ; (if cal image) grade0 = blank -1. ; -1s for overscan grade0(64,0) = grade(64:64+2047,0:2047) + $ spotrb*(spotrb lt 1 and spotrb gt 0)/1d2*bin^2 + $ ; spot dustrb*(dustrb lt 1 and dustrb gt 0)/1d4*bin^2 + $ ; dust & hot pix hotrb*(hotrb lt 1 and hotrb gt 0)/1d6*bin^2 ; per binned pixel grade = grade0 end endif else grade = gradei ; dont want fracs? use int grade ;stop endif if mode eq 2 or (mode eq 1 and flagi eq 1) then begin ; if mode = 2 or ; interpret=1 (interpret grade) image=fltarr(n1,n2,6) ; construct output image dustg = (grade ge 32) ; decode dust growth grade0 = grade - 32.*dustg hot = grade0 ge 16 ; decode hot grade0 = grade - 16.*hot dust = grade0 ge 8 ; decode integer dust grade0 = grade0 - 8.*dust spot = grade0 ge 4 ; decode integer spot grade0 = grade0 - 4.*spot bloom = grade0 ge 2 ; decode saturation bloom/bleed grade0 = grade0 - 2.*bloom sat = grade0 ge 1 ; decode saturated ngrade=[total(sat),total(bloom),total(spot),total(dust),total(hot), $ total(dustg)] ; # pix ;stop del=1d-8 ; small bump to correct rounding errors if flagf eq 1 then begin ; if grade fractions are desired fgrade=grade-fix(grade+del) fspot = fix(fgrade*100.+del)*1d0/bin^2 ; decode spot fraction spot = spot*(abs(fspot) le del) + fspot ; add whole & spot frac to spot fgrade = fgrade*100. - fix(fgrade*100.+del) >0 fdust = fix(fgrade*100.+del)*1d0/bin^2 ; decode dust fraction dust = dust*(abs(fdust) le del) + fdust ; add whole & dust frac to dust fgrade = fgrade*100. - fix(fgrade*100.+del) >0 fhot = fix(fgrade*100.+del)*1d0/bin^2 ; decode hot fraction hot = hot*(abs(fhot) le del) + fhot ; add whole & hot frac to hot ;stop endif image[0,0,0]=sat ; insert into image image[0,0,1]=bloom image[0,0,2]=spot image[0,0,3]=dust image[0,0,4]=hot image[0,0,5]=dustg endif return END ;====================================================================== ; =========================================================================