;+ ; NAME: ; XRT_TUP_CONTAM_M ; PURPOSE: ; Touch-up contaminants feature on the XRT CCD and make the images ; less untolerable to see... ; CALLING SEQUENCE: ; XRT_TUP_CONTAM, dat_in, dat_out, spotmap, gband=gband, dc=dc ; INPUT: ; dat_in - Image data cube (or a single image). ; OPTIONAL (KEYWORD) INPUT: ; spotmap - Array identifying location of the contaminants ; feature on the CCD. ; OUTPUT: ; RETURNED - Array with touch-up made. ; SIDE EFFECT: ; None. ; CATEGORY: ; Hinode XRT data analysis. ; HISTORY: ; Version 1 2007/08/10 (Fri) T.Sakao written. For G-Band images. ; 2 2007/08/12 (Sun) Working version. ; 3 2007/08/14 (Tue) Worked also for X-ray images (Al/Mesh). ; 4 2007/08/15 (Wed) MEDIAN used instead of Average for ; Ring Data to avoid effect of ; possible spike pixels. ; 5 2007/08/24 (Fri) M.Weber. Replaced for-do loop over ; blobs with array operations. ; Changed name from to ; . ; 6 2007/08/30 (Thu) M.Weber. Implemented D.McKenzie's fix ; for missing edge pixels. ;- ; ; FUNCTION XRT_TUP_CONTAM_M, dat_in, spotmap st_time = systime(1) sz = SIZE(dat_in) ;; nloop = number of images in datacube. IF sz(0) EQ 2 THEN nloop = 1 $ ELSE nloop = sz(3) nx = sz(1) ; X pixels ny = sz(2) ; Y pixels ;; This next little bit with damemap and incl_border ;; is because doesn't count edge pixels, ;; so we add a fake edge, and then remove it after we ;; we get res and incl_border. damemap = fltarr(nx+2,ny+2) damemap[1:nx,1:ny] = spotmap incl_border = $ damemap + SHIFT(damemap,1,0) + SHIFT(damemap,-1,0) + SHIFT(damemap,0,1) $ + SHIFT(damemap,0,-1) nz = WHERE(incl_border GT 0) incl_border(nz) = 1 res = LABEL_REGION(incl_border,/all) res = res[1:nx,1:ny] incl_border = incl_border[1:nx,1:ny] ;; Map of spot rings. rings = incl_border - spotmap ;; Rings as labeled regions. res_rings = long(res * rings) ;; Number of spots. n = MAX(res) ;; Maximum ring size (in numbers of pixels). mrs = max((histogram(res_rings))(lindgen(n)+1)) ;; === This section prepares addresses for the spot rings, ;; and requires a loop over spots. However, it only runs once, ;; before the data are introduced. ;; ss_rings will hold array addresses for each ring. ;; The 1st dimension of ss_rings spans the number of ;; pixels associated with the ring for a given spot, ;; but may be padded with zeros. The 2nd dimension spans ;; the number of spots, but has a padded "0" entry, so ;; that the 2nd dimension index matches the spot ID. ss_rings = lonarr(mrs,n+1) - 1L array_idxs = lindgen(nx,ny) ss_sort = sort(res_rings) array_idxs = array_idxs[ss_sort] res_rings2 = res_rings[ss_sort] ss_uniq = uniq(res_rings2) starts = shift(ss_uniq,1)+1 starts[0] = 0 lengths = ss_uniq - shift(ss_uniq,1) lengths[0] = starts[1] for ispot = 1,n do begin ss_rings[0:lengths[ispot]-1,ispot] = $ array_idxs[starts[ispot]:starts[ispot]+lengths[ispot]-1] endfor mask1 = (ss_rings ne -1L) ;; === This section handles the datacube. dat_out = dat_in FOR j = 0, nloop-1 DO BEGIN image = dat_in[*,*,j] ;; Some valid ring pixels will correspond to data le 0, ;; and we don't want to count those. data1 = (image[ss_rings] * mask1) > 0 ;; q_data = each column has ones or zeros. One = valid ring ;; data value gt 0. Sum over ones for that spot says how ;; many ring values gt 0. q_data = data1 gt 0 npix = total(q_data,1) ;; npix = 1D array; [N_spots + 1] ;; Do sorting on whole set of blobs, all at once. ;; We assume this is Level-0 data, so every pixel value ;; is within the range [0,4095]. We separate spots by adding ;; (Spot_ID * 5000) to the data values for that ring. That way, ;; we can sort the data values for each ring, but we do not ;; mix the rings. data1 = temporary(data1) + (lindgen(mrs,n+1) / mrs * 5000) data2 = data1[sort(data1)] ss_med = ((lindgen(n)+2) * mrs) - 1 - long((npix[1:*]-1)/2) med_vals = data2[ss_med] - ((lindgen(n)+1)*5000) ;; Make med_vals's indices match up again with the region labels. med_vals = [0, med_vals] ;; Payoff. dat_out[*,*,j] = image*(1-spotmap) + med_vals[res]*spotmap ENDFOR en_time = systime(1) runtime = en_time - st_time return, dat_out END ;;=========================================================================== PRO xrt_tup_contam, ind_in, dat_in, ind_out, dat_out, spotmaps=spotmaps, $ runtime=runtime, batch=batch, verbose=verbose, $ quiet=quiet, qabort=qabort, threshold = threshold ; ========================================================================= ;+ ; ; PROJECT: ; Solar-B / XRT ; ; NAME: ; ; XRT_TUP_CONTAM ; ; CATEGORY: ; ; Image calibration ; ; PURPOSE: ; ; Replace CCD contamination "spots" with neighbor median patching, ; and update the metadata. ; ; CALLING SEQUENCE: ; ; XRT_TUP_CONTAM, ind_in, dat_in, ind_out, dat_out ; [,spotmaps=spotmaps] [,runtime=runtime] [,batch=batch] ; [,verbose=verbose] [,quiet=quiet] [,qabort=qabort] ; [,threshold = threshold] ; ; ; INPUTS: ; ; IND_IN - [Mandatory] (structure array, [Nimg]) ; This is the INDEX corresponding to the DATA ; that is to be processed. ; DAT_IN - [Mandatory] (number array, [Nx,Ny,Nimg]) ; The datacube to be cleaned. It may have ; a mixture of ROI shapes and CHIP_SUM levels ; (see Note #2). ; THRESHOLD - [Optional] value from 0 to 1 indicating level of ; fractional spottedness a binned pixel must have ; before it is deemed a spot and corrected. ; Default = 1 (1x1) or 0.5 (other binnings). ; ; KEYWORDS: ; ; /BATCH - [Optional] (Boolean) ; This switch is only for XRT institutions ; running large pipeline batch jobs. It avoids ; difficulties with running in batch, but greatly ; slows the program, so should not be used by ; end-users. ; /VERBOSE - [Optional] (Boolean) If set, print out extra ; information. Overrides "/quiet" (see Note #2). ; /QUIET - [Optional] (Boolean) If set, suppress messages ; (see Note #2). ; ; OUTPUTS: ; ; IND_OUT - [Mandatory] (structure array, [Nimg]) ; This is the modified INDEX corresponding to the ; DATA that has been processed. The HISTORY tag ; will have an added line from this routine. ; DAT_OUT - [Mandatory] (number array, [Nx,Ny,Nimg]) ; The spot-corrected datacube. ; SPOTMAPS - [Optional] (integer array, [Nx,Ny,Nimg]) ; A datacube matched to DAT_OUT. Values are: ; 0: pixel was not spot-corrected. ; 1: pixel was spot-corrected; ; RUNTIME - [Optional] (number scalar) ; The runtime for this routine, in seconds. ; 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 call to get correct a datacube and indices: ; IDL> xrt_tup_contam, ind_in, dat_in, ind_out, dat_out ; ; Get maps of which pixels were patched for each image: ; IDL> xrt_tup_contam, ind_in, dat_in, ind_out, dat_out, spotmaps=spotmaps ; ; COMMON BLOCKS: ; ; xrt_spotblock ; ; NOTES: ; ; 1) This routine is basically a wrapper for . ; ; 2) This routine is designed to handle a mixture of ROI shapes and ; levels of on-chip summing. However, each smaller image must ; be in the lowest addresses of its slice of the datacube. ; For example, if a [2048,2048,10] datacube has 9x images 2048x2048, ; and has 1x image 512x512 as the last image, then the small image ; should be at [0:511,0:511,9]. If instead it is located at ; [768:1279,768:1279,9], then the program will produce incorrect ; results. ; ; 3) There are three levels of verbosity. ; a) "verbose" = highest priority. All errors and messages are ; displayed. ("if q_vb") ; b) "quiet" = lower priority. No errors or messages are ; displayed. ("if q_qt") ; c) neither = lowest priority. All errors and some messages are ; displayed. ("if not q_qt") ; ; 4) Only certain sets of data need to be spot-corrected. Here are ; the criteria: ; - FW1 = any of these {Open, Al-poly, C-poly, Be-thin}; AND ; - FW2 = any of these {Open, Al-mesh, Ti-poly, G-band}; AND ; - Observation date > xrt_spottime_st ("2007-JUL-23", currently); AND ; - Observation date < xrt_spottime_en ("2008-FEB-01", currently); AND ; - Image is NOT a Dark; AND ; - Image has not been previously processed by . ; ; 5) Corrections for 4x4 and (especially) 8x8 binning are not as ; effective, due to spot blending/fractionization ; ; MODIFICATION HISTORY: ; progver = 'v2007.Aug.31' ;--- (MW) Written. Designed as a wrapper to ; . progver = 'v2007.Sep.03' ;--- (MW) Updated the selection criteria to include ; Be-thin, filter combinations, and HISTORY. ; Also added spotmaster "ver" and "src" to ; common block. progver = 'v2008.Feb.20' ;--- (MW) Updated to allow for existence of ; multiple spotmaps through time. progver = 'v2008.Nov.21' ;--- (MW) Changed path separators from ; forward slashes to path_sep(). progver = 'v2009.Jun.29' ;--- (RK) Only the pixel occupied with spots is ; defined as the spot in a binned image. progver = 'v2009.Jul.29' ;--- (KR) Fixed typo in rebinned spotmap progver = 'v2013.Apr.03' ;--- (SS) Allowed variable threshold for designation ; of a spot in binned environment, modified ; history to include threshold. Added note 5 ; progver = 'v2025.Jul.11' ;--- (LRG) Changed hardcoding of file location in batch mode ;- ; ========================================================================= ;=== Initial setup ==================================================== if (not keyword_set(batch)) then begin common xrt_spotblock, xrt_spotmaster, xrt_spottime_st, $ xrt_spottime_en, xrt_spotmstr_ver, xrt_spotmstr_src endif st_time = systime(1) ps = path_sep() ;=== Set Booleans which control print statements. ;=== Keyword "verbose" overrules "quiet". q_vb = keyword_set(verbose) q_qt = keyword_set(quiet) and (not q_vb) ;=== Initialize program constants. prognam = 'XRT_TUP_CONTAM' prognul = ' ' qabort = 0B qstop = keyword_set(qstop) q_thr = keyword_set(threshold) ;=== If the common block hasn't yet been defined, ;=== then define it. This should reduce the number of times ;=== that the "geny" reference file must be read. if n_elements(xrt_spotmaster) eq 0 then begin file_prefix = get_logenv('$SSW_XRT') + ps + 'idl' $ + ps + 'util' + ps restgenx, file=file_prefix+'xrt_spotmasters', xrt_spotmaster, $ xrt_spottime_st, xrt_spottime_en, xrt_spotmstr_ver, $ xrt_spotmstr_src endif ;=== Check inputs ===================================================== ;=== Check for 4 mandatory parameters. if q_vb then print, prognam+': Checking input parameters...' if q_vb then print, prognul+' ...OK' if (n_params() ne 4) then begin if (not q_qt) then box_message, prognam+': Requires 4 parameters. ' $ + 'Aborting...' qabort = 1B return endif if q_vb then print, prognul+' ...OK' ;===== Verify that index0 has necessary tags. if q_vb then print, prognam+': Checking whether this data has ' $ + 'necessary XRT header tags...' if not required_tags(ind_in, 'HISTORY,EC_FW1,EC_FW2,DATE_OBS,' + $ 'EC_IMTYP,RPOS_COL,RPOS_ROW,RSIZ_COL,RSIZ_ROW,' + $ 'CHIP_SUM', missing=missing) then begin if (not q_qt) then box_message, [prognam+': This input header ' $ + 'does not have necessary XRT tags:', $ prognul+' '+missing, $ prognul+' Aborting...'] qabort = 1B return endif if q_vb then print, prognul+' ' $ + ' ...OK' ;=== Main section ===================================================== ;=== Prepare the output variables. ind_out = ind_in dat_out = dat_in spotmaps = byte(dat_in) * 0b hist_val = strarr(n_elements(ind_in)) + 'No correction necessary.' ;=== Select only those images which need to be patched. ;=== The first block identifies by HISTORY tag those images ;=== which haven't been checked previously. unfixed = bytarr(n_elements(ind_in)) + 1 dummy1 = grep(prognam, ind_in.history, index=ss_prev) if (ss_prev[0] ne -1) then begin ss_prev = ss_prev / n_elements(ind_in[0].history) ss_prev = (ss_prev[uniq(ss_prev)]) unfixed[ss_prev] = 0 ;; Discount previously fixed indexes. endif ;=== Since there are multiple maps, we do all ;=== the corrections for a map, in a row, before moving ;=== on to the next map. nmaps = n_elements(xrt_spottime_st) for imap = 0,(nmaps-1) do begin ;=== This second block checks various selection criteria. ;=== See Note #4 to see what the criteria are in English. ss_fix = where((ind_in.ec_fw1 le 3) and (ind_in.ec_fw2 le 3) and $ (anytim(ind_in.date_obs) ge anytim(xrt_spottime_st[imap])) and $ (anytim(ind_in.date_obs) lt anytim(xrt_spottime_en[imap])) and $ (ind_in.ec_imtyp eq 0) and (unfixed), cnt_filt) ;=== Now do corrections if there are any to be made. if (cnt_filt gt 0) then begin ;=== Want to identify images by their combination of ;=== ROI size (xsiz, ysiz), ROI position (xpos, ypos), ;=== and level of on-chip summing (chsm). ;=== Want to sort and group by unique combinations. ;=== To do this, we make a single number key that combines ;=== all 5 parameters. key1 = ind_in[ss_fix].rpos_col/8d $ + ind_in[ss_fix].rpos_row*256d/8 $ + ind_in[ss_fix].rsiz_col*(256d^2)/64 $ + ind_in[ss_fix].rsiz_row*(256d^2 * 33)/64 $ + ind_in[ss_fix].chip_sum*(256d^2 * 33d^2) ;=== Unique values of key1: key_vals = key1[uniq(key1, sort(key1))] n_keys = n_elements(key_vals) ;=== Loop over unique key values. for ikey = 0,(n_keys-1) do begin ss_key = where(key1 eq key_vals[ikey], cnt_sskey) xpos = ind_in[ss_fix[ss_key[0]]].rpos_col ypos = ind_in[ss_fix[ss_key[0]]].rpos_row xsiz = ind_in[ss_fix[ss_key[0]]].rsiz_col ysiz = ind_in[ss_fix[ss_key[0]]].rsiz_row chsm = ind_in[ss_fix[ss_key[0]]].chip_sum if q_thr eq 0 then threshold = 0.5 + 0.5*(chsm eq 1) ; default thresh=1 ;=== Now make a spotmap for this ROI and CHIP_SUM from ;=== the master spotmap. sptmp = xrt_spotmaster[xpos:xpos+xsiz-1,ypos:ypos+ysiz-1,imap] ;=== This next bit is to apply chip_sum rebinning to sptmp. if (chsm ne 1) then begin ; if image is binned ; ; sptmp=rebin(temporary(sptmp)*chsm^2,xsiz/chsm,xsiz/chsm) ne 0 ; This is equivalent with 'v2008.Nov.21' ; sptmp=rebin(temporary(sptmp)*chsm^2,xsiz/chsm,ysiz/chsm) eq chsm^2 ; Only the pixel occupied with spots is defined as the spot in a binned image. ; sptmp0=rebin(sptmp*chsm^2,xsiz/chsm,ysiz/chsm) ne 0 ; output map has ; ANY % spots >0 sptmp=rebin(temporary(sptmp)*chsm^2,xsiz/chsm,ysiz/chsm) ge $ threshold*chsm^2 ; For the purposes of cosmetic correction tho, ; only the pixel occupied >=thresh fraction of ; spots is defined as the spot in a binned image endif else sptmp0 = sptmp ; else set spotmap out = used ;=== Pass to subroutine that does spot corrections. dat_out[0:xsiz/chsm-1,0:ysiz/chsm-1,ss_fix[ss_key]] = $ xrt_tup_contam_m(dat_in[0:xsiz/chsm-1,0:ysiz/chsm-1,ss_fix[ss_key]], $ sptmp) ;=== Update the output variable spotmaps. for ii = 0,(cnt_sskey-1) do $ spotmaps[0:xsiz/chsm-1,0:ysiz/chsm-1,ss_fix[ss_key[ii]]] = sptmp0 ;=== Make note for the image histories. if chsm ne 1 then hend = ', threshold = ' + string(threshold) + '.' $ else hend = '.' hist_val[ss_fix[ss_key]] = 'Spots corrected with neighborhood ' + $ 'medians. Spot mastermap version = ' + xrt_spotmstr_ver[imap] + hend endfor ;; ikey loop endif ;; cnt_filt endfor ;; imap loop ;=== Finish =========================================================== ;=== Now write the correction notes to the image HISTORY tag. hist_prefix = prognam + ' ' + progver + ': (' + systim() + ') ' hist_entry = hist_prefix + hist_val update_history, ind_out, hist_entry, /noroutine, /mode en_time = systime(1) runtime = en_time - st_time END ;======================================================================