;+ ;NAME: ; IMSCRUB ;PURPOSE: ; Remove spikes, lines, cosmic ray hits, etc. from an otherwise ; smooth image. Bad pixels are identified by contrast between ; the image and a median-smoothed image. Then the bad pixels are ; replaced from a smoothed image created by ck_convol (which handles ; bad pixels more intelligently than IDL's convol) ;CALLING SEQUENCE: ; result = imscrub(image [, thresh=thresh] [, medwidth=medwidth] $ ; [, image_marked=image_marked) ;INPUT PARAMETERS: ; image --- 2d float or double array. If there are negative pixels, ; imscrub may not behave well. ;OPTIONAL KEYWORD INPUTS: ; thresh --- intensity-normalized threshold for declaring a ; bad pixel. If a pixel differs from the median smoothed ; image by more than thresh * median, then ; the pixel is bad. If the signal is in counts, and the ; noise is mainly Poisson, then thresh is roughly the ; number of standard deviations. Default = 0.5. ; medwidth --- width of the median neighborhood. Default = 5 pixels. ; expand_bad --- if set, expand the width of the bad pixel regions ; by as many pixels as the keyword is set to. ;OPTIONAL KEYWORD OUTPUTS: ; image_marked --- copy of image with bad pixels set to NaN. ;MODIFICATION HISTORY: ; 2006-Jun-19 C. Kankelborg ;- function imscrub, image, thresh=thresh, medwidth=medwidth, $ image_marked=image_marked, expand_bad=expand_bad isize = size(image) Nx = isize[1] Ny = isize[2] NaN = !values.f_nan ;used for marking bad data. maxit = 16 ;maximum number of smoothing iterations kernel = [[1,2,1],[2,4,2],[1,2,1]]/16.0 ;fairly minimal smoothing kernel if not keyword_set(thresh) then thresh = 0.5 if not keyword_set(medwidth) then medwidth = 5 ;First, create median smoothed image: med = median(image, medwidth) ;Find and mark bad pixels. badpix = float(abs(image - med) gt thresh*med) ;Expand bad pixel map if desired. if keyword_set(expand_bad) then begin for i=1,expand_bad do begin badpix = convol(badpix, kernel, /edge_truncate) endfor endif ss = where(badpix) ss_good = where(~badpix) image_marked = image if (ss[0] eq -1) then begin message,'No bad pixels found.',/informational return, image endif else begin print,'Found ',n_elements(ss),' bad pixels.' endelse image_marked[ss] = NaN ;Replace bad pixels as appropriate smoothed = image_marked for i = 1L,maxit do begin print,' Smoothing iteration ',i smoothed = ck_convol(smoothed, kernel, $ method = 'redeem_taint', /edge_truncate) smoothed[ss_good] = image[ss_good] ;no need to smooth good pixels! if min(finite(smoothed[ss])) then break endfor result = image_marked result[ss] = smoothed[ss] return,result end