FUNCTION xrt_despike, image_in, sens=sens, spike_map=ss, hist=hist, $ verbose=verbose, quiet=quiet, qabort=qabort, $ qstop=qstop ; ========================================================================= ;+ ; PROJECT: ; Solar-B / XRT ; ; NAME: ; XRT_DESPIKE ; ; 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_DESPIKE(image_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. ; ; KEYWORDS: ; ; SENS - [Optional] (number scalar) ; Sensitivity to spots. If a pixel has a value greater ; than (SENS * Mean(neighbors)), then it is identified ; as a spike. Values of 1.1 to 1.30 are reasonable. ; Default value = 1.1. ; /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_despike(image_in) ; ; Find out which pixels were altered: ; IDL> image_out = xrt_despike(image_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 average neighbor. (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.1. 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) 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 = 'v2007.May.18' ; --- (P.Jibben, M.Weber) Written. Jibben derived ; the cleaning algorithms based upon the ; ANA software and "The Handbook of ; Astronomical Image Processing" by Berry & ; Burnell. Weber implemented a faster version ; of the code. ; ;- ; ========================================================================= ;=== 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_DESPIKE' prognul=' ' def_sens = 1.1 ;=== 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 1) 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 avg for nearest neighbors. kernel = [[1,1,1],[1,0,1],[1,1,1]] / 8. img_avg = convol(image_in, kernel) ;=== Compare pix/avg versus threshold, ;=== for pixels not near the edge. ss = where( (image_in[2:nx-3,2:ny-3] / img_avg[2:nx-3,2:ny-3]) $ ge sens, 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) ;=== Replace spikes with medians. image_out = image_in image_out[ss] = spike_med ;=== Else if didn't find any spikes. endif else image_out = image_in 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 ;======================================================================