function sigma_filter, image, box_width, N_SIGMA=Nsigma, ALL_PIXELS=all, \$ ITERATE=iterate, MONITOR=monitor, \$ KEEP_OUTLIERS=keep, RADIUS=radius, \$ N_CHANGE=nchange, VARIANCE_IMAGE=imvar, DEVIATION_IMAGE=imdev ;+ ; NAME: ; SIGMA_FILTER ; PURPOSE: ; Replace pixels more than a specified pixels deviant from its neighbors ; EXPLANATION: ; Computes the mean and standard deviation of pixels in a box centered at ; each pixel of the image, but excluding the center pixel. If the center ; pixel value exceeds some # of standard deviations from the mean, it is ; replaced by the mean in box. Note option to process pixels on the edges. ; CALLING SEQUENCE: ; Result = sigma_filter( image, box_width, N_sigma=(#), /ALL,/MON ) ; INPUTS: ; image = 2-D image (matrix) ; box_width = width of square filter box, in # pixels (default = 3) ; KEYWORDS: ; N_sigma = # standard deviations to define outliers, floating point, ; recommend > 2, default = 3. For gaussian statistics: ; N_sigma = 1 smooths 35% of pixels, 2 = 5%, 3 = 1%. ; RADIUS = alternative to specify box radius, so box_width = 2*radius+1. ; /ALL_PIXELS causes computation to include edges of image, ; /KEEP causes opposite effect: pixels with values outside of specified ; deviation are not changed, pixels within deviation are smoothed. ; /ITERATE causes sigma_filter to be applied recursively (max = 20 times) ; until no more pixels change (only allowed when N_sigma >= 2). ; /MONITOR prints information about % pixels replaced. ; Optional Outputs: ; N_CHANGE = # of pixels changed (replaced with neighborhood mean). ; VARIANCE = image of pixel neighborhood variances * (N_sigma)^2, ; DEVIATION = image of pixel deviations from neighborhood means, squared. ; CALLS: ; function filter_image( ) ; PROCEDURE: ; Compute mean over moving box-cars using smooth, subtract center values, ; compute variance using smooth on deviations from mean, ; check where pixel deviation from mean is within variance of box, ; replace those pixels in smoothed image (mean) with orignal values, ; return the resulting partial mean image. ; MODIFICATION HISTORY: ; Written, 1991, Frank Varosi and Dan Gezari NASA/GSFC ; F.V.1992, added optional keywords /ITER,/MON,VAR=,DEV=,N_CHANGE=. ; Converted to IDL V5.0 W. Landsman September 1997 ;- if N_elements( radius ) EQ 1 then box_width = 2*radius+1 else begin if N_elements( box_width ) NE 1 then box_width=3 box_width = 2*(fix( box_width )/2) + 1 ;make sure width is odd. endelse if (box_width LT 3) then return,image bw2 = box_width^2 mean=( filter_image( image,SMO=box_width,ALL=all )*bw2 - image )/(bw2-1) if N_elements( Nsigma ) NE 1 then Nsigma=3 if (Nsigma LE 0) then return, mean imdev = (image - mean)^2 fact = float( Nsigma^2 )/(bw2-2) imvar = fact*( filter_image( imdev,SMO=box_width,ALL=all )*bw2 - imdev ) if keyword_set( keep ) then wok = where( imdev GE imvar, nok ) \$ else wok = where( imdev LT imvar, nok ) npix = N_elements( image ) nchange = npix - nok if keyword_set( monitor ) then \$ print, nchange*100./npix, Nsigma, \$ FORM="(F6.2,' % of pixels replaced, N_sigma=',F3.1)" if (nok EQ npix) then return,image if (nok GT 0) then mean[wok] = image[wok] if keyword_set( iterate ) AND (Nsigma GE 2) then begin iterate = iterate+1 if (iterate GT 20) then begin iterate = 1 return,mean endif return, sigma_filter( mean, box_width, N_SIGMA=Nsigma, ALL=all,\$ KEEP=keep, ITER=iterate, MONIT=monitor ) endif return, mean end