PRO xrt_clean_ro, image_in, index_in, dark, image_out, grade, model, $ hist=hist, clean_type=clean_type, $ bsub_type=bsub_type, dark_type=dark_type, nsigma=nsigma, $ nmed=nmed, stop_bleed=stop_bleed, verbose=verbose, $ quiet=quiet, qabort=qabort, qstop=qstop ; ========================================================================= ;+ ; PROJECT: ; Solar-B / XRT ; ; NAME: ; XRT_CLEAN_RO ; ; CATEGORY: ; Data calibration ; ; PURPOSE: ; Remove periodic and other "background" signals introduced into ; the data by the read-out electronics. ; ; CALLING SEQUENCE: ; XRT_CLEAN_RO, image_in, index_in, dark, image_out [,grade] ; [,model] ,hist=hist [,clean_type=clean_type] ; [,bsub_type=bsub_type] [,dark_type=dark_type] ; [,nsigma=nsigma] [,nmed=nmed] [,/stop_bleed] ; [,/verbose] [,/quiet] [,qabort=qabort] [,/qstop] ; ; INPUTS: ; IMAGE_IN - [Mandatory] (2-dim float array, [Nx,Ny]) ; The image containing read-out signals. ; INDEX_IN - [Mandatory] index structure associated with IMAGE_IN ; DARK - [Mandatory] (medianed dark) to be used for either ; dark subtraction or background level adjustment. ; Typically generated by . ; (may be non-existant, in which case it is generated) ; GRADE - [Optional] a constant indicating the saturation level ; (default = 2500 DN) ; OR ; a 2-dim byte array, [Nx,Ny] containing ; (at a minimum) ; 0 if pixel is OK ; 1 if pixel is saturated (>2500 DN) ; (note: overwritten on output) ; ; KEYWORDS: ; HIST - [Mandatory] (string array) Contains information ; on type of cleaning done and keyword values ; for later use in updating header information ; CLEAN_TYPE - [Optional] Type of Fourier cleaning to use: ; 0 = prefilter on semi-fixed streaks, then remove ; remaining Fourier streaks, stars, & smudges ; [default] ; 1 = remove semi-fixed streaks only ; 2 = remove Fourier streaks, stars, & smudges ; 3 = don't Fourier clean! ; BSUB_TYPE - [Optional] Type of background subtraction to use: ; 0 = remove Nyquist ringing and large-scale background ; "ramp" (dark) [default] ; 1 = remove Nyquist ringing only ; 2 = remove large-scale background "ramp" (dark) only ; 3 = no (model) background subraction! ; DARK_TYPE - [Optional] Type of dark subtraction to use: ; 0 = Use model dark with median adjusted to median of ; cleaned darks nearby in time (if possible) [default] ; 1 = Use median of cleaned darks nearby in time ; (if possible) ; 2 = Use model darks without median level adjustment ; (Note: if dark_type=1, bsub_type = 3 is used) ; NSIGMA - [Optional] Number of standard deviations ; beyond which an FFT pixel is ; considered bad in "Fourier star" removal. ; Default is 4.5 ; (Note: much below 4 not recommended) ; NMED - [Optional] Number of standard deviations ; above smoothed central minimum background ; to begin shielding (presumed) data from correction ; Default is 3.5 ; (Note: outside range of 2.0 to 4.5 not recommended) ; /STOP_BLEED - [Optional] (Boolean) Set = 1 if pixels corrected for ; saturation "bloom/bleed" at the edge of saturated ; region(s) are to be retained (otherwise they are ; reverted to their input values) ; /QUIET - [Optional] (Boolean) Set for no messages or ; errors. Default is some messages and all errors. ; /VERBOSE - [Optional] (Boolean) Set for all messages, errors, and ; intermediate data listings. Suppresses quiet. ; /QSTOP - [Optional] (Boolean) For debugging. ; ; OUTPUTS: ; IMAGE_OUT - [Mandatory] (2-dim float array, [Nx,Ny]) ; The image with the model of read-out signals removed. ; GRADE - [Optional] a 2-dim byte array (1x1 binning case) or ; a 2-dim double array (other binnings), [Nx,Ny] containing ; the sum of ; 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 ; 16 if pixel is "hot" ; PLUS (optionally, in case of binning, if the /grade_frac ; flag is set) an encoding (in fractional form). ; MODEL - [Optional] (2-dim float array, [Nx,Ny]) ; The model of read-out signals which has been removed ; (subtracted). ; 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. ; ; EXAMPLES: ; ; Basic usage: ; IDL> xrt_clean_ro, image_in, index_in, dark, image_out, hist=hist ; ; You want the the pixel grade and read-out/dark model: ; IDL> xrt_clean_ro, image_in, index_in, dark, image_out, grade, model, $ ; hist=hist ; ; You want to correct saturation bleed-affected pixels, and ; pixel grades: ; IDL> xrt_clean_ro, image_in, index_in, dark, image_out, grade, $ ; hist=hist, stop_bleed=1 ; ; You want pixel grades, but no Fourier cleaning: ; IDL> xrt_clean_ro, image_in, index_in, dark, image_out, grade, $ ; hist=hist, clean_type=3 ; ; You want more aggressive Fourier cleaning, with median dark ; subtraction: ; IDL> xrt_clean_ro, image_in, index_in, dark, image_out, $ ; hist=hist, nsigma=4, nmed=4, dark_type=1 ; ; ; ; COMMON BLOCKS: ; None. ; ; NOTES: ; 1) An XRT image shows a prevelent sinusoidal ringing in x (width) ; at the nyquist frequency, with a peak-to-peak amplitude of ; ~2.6 DN due to different bias levels set in odd vs. even columns. ; There is also a low-level "ski-ramp" in the y (height) ; direction with a shape which can be approximated by an exponential ; decrease (amplitude ~ 4.3 DN, e-folding width ~185 pixels) + a ; weak linear increase on a base of ~42 DN (for 1x1 binning). ; (This is best seen in dark frames). The Fourier power spectrum of ; an XRT dark shows many ; odd features, including "streaks" (narrow ridges of power spanning ; all y, many are semi-fixed in location), "stars" (2-D Voigt ; profiles) and "smudges" (truncated streaks). The latter two ; features are typically variable in y position; all vary ; in amplitude. These features are present (though less clearly ; visible) in the Fourier transform of data as well. This ; procedure attempts to correct for these read-out signals. ; 2) The Fourier filtering can take a while on large images ; (~15s on a 2048x2048 using a MacPro with dual 3GHz Xeon chips). ; Its effects can be dramatic or minimal, but are typically ; confined to low count levels. ; Thus, for certain types of images and certain applications, the ; user may want to skip this part of the cleaning (set clean_type=3). ; See xrt_fourier_vacuuum.pro notes for more details ; on how to use and modify the Fourier cleaning. ; 3) Some Fourier patterns may not be fully removed in images with one ; or more smaller dimensions (<=128). Fourier patterns are not ; removed if a large part of the image (>45-80% depending on ; image size) is saturated. ; ; MODIFICATION HISTORY: ; progver = 'v2007.Feb.14' ; --- (SS, MW) Written by SS. Cleaned up by MW. progver = 'v2007.Mar.06' ; --- (SS, MW) Modified by SS to include ; --- large-scale background removal, reorganize ; --- using function calls . Cleaned up by MW. progver = 'v2007.Apr.23' ; --- (SS, MW) Modified by SS to include a first ; --- stab at Fourier read-out signal removal. ; --- Cleaned up by MW. progver = 'v2007.May.18' ; --- (SS) Modified by SS to include ; --- "patch-over" of saturated regions to ; --- improve Fourier read-out signal ; --- removal, plus a few other minor fixes. progver = 'v2007.Nov.27' ; --- (SS) Added comments/history update about ; --- possible poorer performance for small images. progver = 'v2007.Nov.28' ; --- (SS) Added fix for cases where images are ; --- significantly saturated. progver = 'v2008.May.18' ; --- (SSaar) Added an improved large-scale ; --- background correction (see ) ; --- including corrections for binning, CCD ; --- temperature, and exposure time. progver = 'v2009.Jun.12' ; --- (SSaar) Added an improved dark subtraction ; --- adjusting model results with median of darks ; --- observed nearby in time. Added generation ; --- of initial pixel "grade" map. Added options for ; --- different kinds of dark subtraction. Made ; --- some minor corrections. progver = 'v2009.Jul.08' ; --- (SSaar) Fixed a bug in dark zero level reset progver = 'v2009.Nov.10' ; --- (SSaar) Fixed a bug when i_bleed does not exist ; ; ---- or = -1 progver = 'v2010.Dec.02' ; --- (SSaar) Fixed a bug when bin=4 and improved ; ; --- correction for dark_type=1 progver = 'v2012.Feb.15' ; --- (SSaar) Fixed a bug when i_bleed does ; not exist and clean_type=3 progver = 'v2013.Feb.15' ;--- (SSaar (SAO)) Fixed problem arising when ; ;--- bad images are mislabeled as darks. ; ;--- Simplified the code structure somewhat; ; ;--- improved reuse of already processed median darks ; ;--- (should speed code somewhat). progver = 'v2014.Feb.06' ;--- (JSlavin (SAO)) Pass along quiet flag to ;--- xrt_med_dark ; ;- ; ========================================================================= ;=== Check inputs & initial setup ===================================== ;===== Set some keyword Booleans, plus miscellaneous. q_vb = keyword_set(verbose) q_qt = keyword_set(quiet) and (not q_vb) qstop = keyword_set(qstop) qabort = 0B if not keyword_set(clean_type) then clean_type=0 ; default to full Fourier if not keyword_set(bsub_type) then bsub_type=0 ; default to full background ; subtraction if not keyword_set(dark_type) then dark_type=0 ; default to median adjusted ; dark model if qstop eq 1 then debug=1 else debug=0 if not keyword_set(nsigma) then nsigma=4.5 ; threshold sigma to flag spike if not keyword_set(nmed) then nmed=3.5 ; shield data above nmed*central_median if not keyword_set(stop_bleed) then stop_bleed=0 ; do not replace bleed pixels if n_elements(dark) lt 1 then $ ; calc medianed dark dark = xrt_med_dark(index_in, ddates=dates, darkff=darkff, nday=15, $ quiet=q_qt) ;stop ;===== Thresholds for "small" image size, maximum saturated pixel % sat_level = 2500. if n_elements(grade) eq 0 then grade=sat_level minsize = 128. maxsat = 0.80 ;===== Check inputs. if q_vb then print, 'XRT_CLEAN_RO: Checking input parameters...' ;=== Check number of parameters. if (n_params() ge 3) then begin ; ===== OK endif else begin if (not q_qt) then box_message, 'XRT_CLEAN_RO: Requires >=3 ' $ + 'parameters. Aborting...' qabort = 1B return endelse ;=== Check size and type of input array. dim_check = (size(image_in,/n_dim) eq 2) typ_check = (datatype(image_in) eq 'FLO') or (datatype(image_in) eq 'DOU') if (dim_check and typ_check) then begin ;=== OK endif else begin if not dim_check then begin if (not q_qt) then box_message, 'XRT_CLEAN_RO: Requires a 2D ' $ + 'image array as input. Aborting...' qabort = 1B endif if not typ_check then begin if (not q_qt) then box_message, 'XRT_CLEAN_RO: Requires that ' $ + 'the image is floating type. Aborting...' qabort = 1B endif return endelse ;===== check clean_type if (clean_type lt 0) or (clean_type gt 3) or $ (clean_type ne fix(clean_type)) then begin if (not q_qt) then box_message, 'XRT_CLEAN_RO: Requires that ' $ + 'clean_type = 0, 1, 2, or 3. Defaulting to clean_type = 0' clean_type=0 endif ;===== check bsub_type if (bsub_type lt 0) or (bsub_type gt 3) or $ (bsub_type ne fix(bsub_type)) then begin if (not q_qt) then box_message, 'XRT_CLEAN_RO: Requires that ' $ + 'bsub_type = 0, 1, 2, or 3. Defaulting to bsub_type = 0' bsub_type=0 endif if q_vb then print, ' ...OK' ;===== Get dimensions of image. s = size(image_in) n1 = s(1) n2 = s(2) x1 = findgen(n1) bin = index_in.chip_sum ;===== If Fourier cleaning is used, temporarily smooth over saturated ;===== areas (as needed) and keep record of altered "bloom/bleed" pixels ;===== If number of saturated pixels > 80% (128x128) to 45% (2048x2048) ;===== of the total, default to no Fourier cleaning (clean_type=3) ;===== and add note to History later satflag=0 i_bleed=-1 ; covers case of sat exist but sat_fudge not used if clean_type ne 3 then begin if n_elements(grade) gt 1 then nsat = total(grade eq 1) else $ nsat = total(image_in ge grade) if nsat gt (n1*n2)*maxsat/(n1*n2/minsize^2)^0.1 then satflag=1 if satflag eq 0 then begin if q_vb then print,'Fudging saturated pixels (if any)...' saturation_fudge,image_in,grade,image,i_bleed endif else begin if q_vb then print,'Defaulting to no Fourier cleaning - too' $ + ' many saturated pixels!' clean_type=3 image = image_in endelse endif else image = image_in if q_vb then begin case 1 of (bsub_type eq 1): txtv=' Remove Nyquist ringing only' (bsub_type eq 2): txtv=' Remove large-scale ramp only' (bsub_type eq 3): txtv=' No Nyquist or ramp subtraction' else: txtv=' [Default] Remove Nyquist ringing and large-scale ramp' end print, 'XRT_CLEAN_RO: beginning Nyquist and ramp subtraction; ' print,txtv endif ;print,'ct: ',clean_type,'bt: ',bsub_type,'dt: ',dark_type dflag=(n_elements(dark) gt 1) ; 0 if no valid dark if dflag eq 0 then mdark = 0. else mdark = median(dark) txt='XRT_CLEAN_RO: ' if bsub_type ne 3 then begin ; (if so, then dark_type <> 1) if bsub_type ne 2 then begin ; if not "ramp only" ;====== Remove Nyquist Frequency ringing image_out1 = no_nyquist(image, index_in, fit1) txt=txt+' Remove Nyquist ringing' endif else image_out1 = image if bsub_type ne 1 then begin ; if not "Nyquist only" ;===== Remove Large-Scale component of background ;===== (Pedestal+Read-out+Dark) txt=txt+' Remove large-scale background ramp ' if dark_type eq 0 then begin ; default "optimum" case: subtract model ; with median adjusted to medianed dark, if ; at all possible (else just median or model) if (bin ne 4 or (bin eq 4 and dflag eq 0)) then begin image_out2 = lsback_away(image_out1, index_in, fit2) image_out2 = image_out2 -mdark +median(fit2)*dflag fit2 = fit2 +mdark -median(fit2)*dflag if dflag eq 0 then txt=txt+' via model dark.' else $ txt=txt+' via median adjusted model dark [default].' endif else begin ; use median darks where possible if bin=4 image_out2 = image_out1 - dark fit2 = dark txt=txt+' via medianed dark.' end endif else begin ; dark_type=3, subtract model dark image_out2 = lsback_away(image_out1, index_in, fit2) txt=txt+' via model dark.' end endif endif else begin if dark_type eq 1 then begin ; subtract median dark if at all possible image_out1 = no_nyquist(image, index_in, fit1) if total(mdark) ne 0 then begin image_out2 = image_out1 - dark fit2 = dark txt=txt+' via medianed dark.' endif else begin ; if no good darks, use model image_out2 = lsback_away(image_out1, index_in, fit2) txt=txt+' via model dark [not enough good darks].' end endif else begin txt=txt + ' NO Nyquist or large-scale ramp removal. image_out2 = image end end if bsub_type eq 0 then txt=txt+' [Default]' ;===== construct history update hist=strarr(3) hist(0)=txt if q_vb then begin print, 'XRT_CLEAN_RO: beginning Fourier cleaning ' case 1 of (clean_type eq 1): txt=' Remove semi-fixed streaks only' (clean_type eq 2): txt=' Remove Fourier streaks, stars, & smudges' (clean_type eq 3): txt=' No Fourier cleaning.' else: txt=' [Default] Prefilter on semi-fixed streaks, then ' $ + ' remove remaining Fourier streaks, stars, & smudges' end if satflag eq 1 then txt = txt + ' Too many saturated pixels! ' print, 'clean_type = ',clean_type,' '+txt if nsigma eq 4.5 then txt=' [default] ' else txt='' print,' nsigma = ',nsigma,txt if nmed eq 3.5 then txt=' [default] ' else txt='' print,' nmed = ',nmed,txt endif ;===== Vacuum up naughty read-out features in Fourier space if clean_type ne 3 then begin xrt_fourier_vacuum,image_out2,image_out, $ clean_type=clean_type,nsig=nsig,nmed=nmed,debug=debug ;===== if stop_bleed=0, change "bleed/bloom" boundary pixels to ;===== original values with nyquist and ramp correction only if ((stop_bleed eq 0) and (i_bleed[0] ne -1)) then begin case 1 of (bsub_type eq 0): corr = fit1(i_bleed)+fit2(i_bleed) (bsub_type eq 1): corr = fit1(i_bleed) (bsub_type eq 2): corr = fit2(i_bleed) else: corr = 0. end image_out(i_bleed)=image_in(i_bleed) - corr endif endif else image_out = image_out2 ;===== add to history update case 1 of (clean_type eq 1): txt=' Semi-fixed Fourier streak preclean' (clean_type eq 2): txt=' Fourier cleaning (no streak preclean)' (clean_type eq 3): txt=' No Fourier cleaning. ' else: txt=' Full Fourier preclean/clean [default]' end if satflag eq 1 then txt = txt + ' Too many saturated pixels! ' hist(1)='XRT_CLEAN_RO: clean_type = ' $ + strtrim(string(fix(clean_type)),2) + txt if clean_type ne 3 then begin hist(2)= 'XRT_CLEAN_RO: nmed = ' + strtrim(string(nmed),2) if nmed eq 3.5 then hist(2) = hist(2) + ' [default]' if clean_type ne 1 then begin hist(2)= hist(2) + ' nsigma = ' + strtrim(string(nsigma),2) if nsigma eq 4.5 then hist(2) = hist(2) + ' [default]' endif endif else hist=hist(0:1) if ((n1 le minsize) or (n2 le minsize)) and ((clean_type eq 0) or $ (clean_type eq 2)) then hist = [hist,'XRT_CLEAN_RO: small image '+ $ 'dimension - some Fourier noise patterns may remain uncorrected'] if stop_bleed eq 1 and satflag eq 0 then hist = $ [hist,'XRT_CLEAN_RO: saturation-affected boundary pixels corrected '] hist=['XRT_CLEAN_RO: Running '+progver,hist] ;===== reset saturation level to input value if n_elements(grade) eq 1 then isat=where(image_in ge grade,nsat) $ else isat=where(grade eq 1,nsat) if nsat gt 0 then image_out(isat) = image_in(isat) ;===== generate fit array model = image_in - image_out ;==== generate initial grade array (saturation and bleed only) grade=bytarr(n1,n2) if nsat gt 0 then begin grade(isat) = 1B if n_elements(i_bleed) gt 0 and i_bleed(0) ne -1 then grade(i_bleed) = 2B ; checks i_bleed is there and not = -1 endif END ;======================================================================