FUNCTION xrt_despike2, image_in, index_in, sens=sens, spike_map=ss, $ hist=hist, verbose=verbose, quiet=quiet, qabort=qabort, $ qstop=qstop ; ========================================================================= ;+ ; PROJECT: ; Solar-B / XRT ; ; NAME: ; XRT_DESPIKE2 ; ; CATEGORY: ; Data calibration ; ; PURPOSE: ; To smooth out "bad" pixels due to cosmic ray hits. Generally, ; this is intended as a cosmetic fix, not as a scientific ; recalibration of the spiked pixels. However, unspiked pixels ; are unaffected, and the SPIKE_MAP can be used to identify the ; corrected pixels. See "Notes" for more information. ; ; CALLING SEQUENCE: ; Result = XRT_DESPIKE2(image_in ,index_in [,sens=sens] ; [,spike_map=spike_map] ; [,hist=hist] [,/verbose] [,/quiet] [,/qstop] ; [,qabort=qabort] ) ; ; INPUTS: ; ; IMAGE_IN - [Mandatory] (2-dim number array, [Nx,Ny]) ; The input image array. ; INDEX_IN - [Mandatory] (index structure corresponding to IMAGE_IN. ; ; KEYWORDS: ; ; SENS - [Optional] (number scalar) ; Fractional sensitivity to spikes. If a pixel has a value ; > (SENS * Median(neighbors)), then it is identified ; as a spike. Values from 1.3 to 2.0 are reasonable. ; Default value = 1.4. ; /VERBOSE - [Optional] (Boolean) If set, print out extra ; information. Overrides "/quiet" (see Note #3). ; /QUIET - [Optional] (Boolean) If set, suppress messages ; (see Note #3). ; /QSTOP - [Optional] (Boolean) For debugging. ; ; OUTPUTS: ; ; Return - [Mandatory] (2-dim number array, [Nx,Ny]) ; The output image array, with the spikes reduced. ; A spike value is replaced with the median of the 16 ; pixels marking a 5x5 square around the spike. ; SPIKE_MAP - [Optional] (Long array, [Nspikes]) ; The addresses of the despiked pixels. These are the ; only pixels altered by this routine. ; HIST - [Optional] (String scalar) ; Brief report on how many spiked pixels were found. ; Called "hist" because it fits the format for an ; index.HISTORY entry. (This keyword is used by ; , for example.) ; 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> image_out = xrt_despike2(image_in, index_in) ; ; Find out which pixels were altered: ; IDL> image_out = xrt_despike2(image_in, index_in, spike_map=spike_map) ; IDL> help, image_in[spike_map] ; ; ; COMMON BLOCKS: ; None. ; ; NOTES: ; ; 1) Spiked pixels are identified if they are too much larger ; than their median neighbor by several measures (absolute, ratio ; and fractionally; See SENS.) Only such pixels are ; altered by this routine. They are replaced with the median ; of a slightly-removed ring of pixels. (See "Return".) ; ; 2) Multiple passes may be required to affect heavy image ; spiking. Generally, one pass does fine for lightly spiked ; images, at SENS=1.4. A second pass should handle most of the ; heavily spiked images. Picky people may need to experiment a ; bit with the number of passes and the SENS level. ; ; 3) Caution should be taken with more highly binned data (4x4 and up) ; as small real features can be mis-identified as spikes (e.g., ; flare kernels, XRBPs) ; ; 4) 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") ; ; ; MODIFICATION HISTORY: ; progver = 'v2014.Jan.28' ; --- (S. Saar) Written. Based on , ; adding more tests to reduce false detections, ; and over-corrections. Much testing by P. Jibben. ; ;- ; ========================================================================= ;=== Initial setup ======================================= ;=== Set Booleans which control print statements. ;=== Keyword "verbose" overrules "quiet". q_vb = keyword_set(verbose) q_qt = keyword_set(quiet) and (q_vb eq 0) ;=== Initialize program constants. prognam='XRT_DESPIKE2' prognul=' ' bin=index_in.chip_sum def_sens = 1.4 ; minimum fractional enhancement min_lev = 5.*bin ; minimum enhancement hist0 = get_history(index_in, 'XRT_PREP', found=found) bak=0. case 1 of ; set min ratio & minimal background if needed (bin eq 1): begin if found eq 0 then bak=41. min_rat=3.5 end (bin eq 2): begin if found eq 0 then bak=144. min_rat=2.5 end (bin eq 4): begin if found eq 0 then bak=348. min_rat=2.5 ; 2.5? end else: begin if found eq 0 then bak=770. min_rat=2.25 ; 2.25? end endcase ;=== Announce version of the program. if q_vb then box_message, prognam+': Running ' + progver + '.' if q_vb then print, prognam+': Performing initial setup...' ;=== Set some keyword Booleans. qstop = keyword_set(qstop) qabort = 0B if q_vb then print, prognam+' ...OK' ;=== Check inputs ======================================== if q_vb then print, prognam+': Checking inputs...' ;=== Check for number of inputs. if (n_params() ne 2) then begin if (not q_qt) then box_message, prognam+': Incorrect number of ' $ + 'parameters. Aborting.' qabort = 1B return, 0 endif ;=== Check for image dimensions. imsiz = size(image_in) if (imsiz[0] ne 2) then begin if (not q_qt) then box_message, prognam+': Input IMAGE_IN must ' $ + 'be a 2D image. Aborting.' qabort = 1B return, 0 endif ;=== Check SENS value. if keyword_set(sens) then begin if (n_elements(sens) ne 1) then begin if (not q_qt) then box_message, [prognam+': Input SENS must ' $ + 'be a scalar.', prognul+' Using the default value = ' $ + strcompress(def_sens,/rem) + '.'] sens = def_sens endif if (not is_number(sens)) then begin if (not q_qt) then box_message, [prognam+': Input SENS must ' $ + 'be a number.', prognul+' Using the default value = ' $ + strcompress(def_sens,/rem) + '.'] sens = def_sens endif if (sens lt 1.0) then begin if (not q_qt) then box_message, [prognam+': Input SENS must ' $ + 'be a number > 1.0.', prognul+' Using the default value = ' $ + strcompress(def_sens,/rem) + '.'] sens = def_sens endif endif else begin sens = def_sens endelse if q_vb then print, prognul+' ...OK' ;=== Find and remove spikes ============================== if q_vb then print, prognam+': Fixing spikes...' nx = n_elements(image_in[*,0]) ny = n_elements(image_in[0,*]) ;=== Find median for nearest neighbors. img_med = median(image_in*1.,5) ; kernel = [[1,1,1],[1,0,1],[1,1,1]] / 8. ; img_avg = convol(image_in*1., kernel) ;=== Compare pix/avg versus threshold & be sure enhancement > min_rate ;=== for pixels not near the edge. im_inx = image_in[2:nx-3,2:ny-3] im_mx = img_med[2:nx-3,2:ny-3] ss =where((im_inx - bak)/((im_mx - bak)>1.0) gt min_rat and $ (im_inx - im_mx)/(im_mx-bak+41.) gt (sens-1.) and $ (im_inx - im_mx) ge min_lev, count) if q_vb then print, prognam+': Found '+strcompress(count,/rem) $ +' pixels to be fixed.' ;=== Found some spikes? Then proceed. if (count ge 1) then begin ;=== Have to convert ss to be relative to entire image. aa = bytarr(nx-4,ny-4) aa[ss] = 1 bb = bytarr(nx,ny) bb[2:nx-3,2:ny-3] = aa ss = where(bb eq 1, count) ;=== Prepare an "address" filter for +2 pixels around spikes. filter = [-2*nx-2, -2*nx-1, -2*nx, -2*nx+1, -2*nx+2, -nx-2, -nx+2, $ -2, 2, nx-2, nx+2, 2*nx-2, 2*nx-1, 2*nx, 2*nx+1, 2*nx+2 ] ;=== spike_ring is the filter applied to each spike address. spike_ring = lonarr(16,count) for ii = 0,15 do spike_ring[ii,*] = ss + filter[ii] ;=== For each spike, find median of its ring. ;=== This bit lets us do them all simultaneously. ; spike_med = median(image_in[spike_ring], dim=1, /even) spike_med = median(img_med[spike_ring], dim=1, /even) ;=== Replace spikes with medians. image_out = image_in image_out[ss] = spike_med ;=== Else if didn't find any spikes. endif else begin image_out = image_in endelse if q_vb then print, prognul+' ...OK' ;=== Finish ============================================== hist = prognam+': Despiked '+strcompress(count,/rem)+' pixels, at a ' $ + 'sensitivity of '+strcompress(sens,/rem)+'.' return, image_out END ;======================================================================