; $Id: leefilt.pro,v 1.12 2001/11/08 00:08:00 dave Exp $ Function Lee_filter_exact,A,N,SIG, $ DOUBLE=double ; ; Slow, but accurate Lee filter. Derived from ; Jong-Sen Lee, Optical Engineering 25(5), 636-643 (May 1986) ; COMPILE_OPT hidden on_error,2 Delta = (LONG(N)-1L)/2L ;Width of window SZ=size(A) N_SAMPLES=SZ[1] ;# of columns ; mean=smooth(a, N, /EDGE_TRUNCATE) mean2=mean^2 ; IF SZ[0] eq 1 THEN BEGIN ;Vector? z = double ? dblarr(n_samples) : fltarr(n_samples) FOR S=Delta,N_SAMPLES-Delta-1 DO $ Z[S]=TOTAL((A[S-Delta:S+Delta] - mean[S])^2) z = z / (n-1) ;Variance ENDIF ELSE BEGIN ;***** 2D case ***** N_LINES=SZ[2] ; ; Compute Variance of Z ; z = double ? dblarr(n_samples, n_lines) : fltarr(n_samples, n_lines) FOR L=Delta, N_LINES-Delta-1 DO $ FOR S=Delta, N_SAMPLES-Delta-1 DO $ Z[S,L]=TOTAL((A[S-Delta:S+Delta,L-Delta:L+Delta] - mean[S,L])^2) z = z / (n^2-1) ;Variance ENDELSE ; ; Upon starting the next equation, Z = Var(Z). Upon exit, Z = Var(X) ; of equation 19 of Lee, Optical Engineering 25(5), 636-643 (May 1986) ; ; VAR_X = (VAR_Z + Mean^2 )/(Sigma^2 +1) - Mean^2 (19) ; ; Here we constrain to >= 0, because Var(x) can't be negative: Z=((TEMPORARY(Z) +mean2) /(1 + Sig^2)) - mean2 > 0 ; ; return value from equation 21,22 of Lee, 1986. ; K = ( VAR_X/(mean^2 * Sigma^2 + VAR_X) ) (22) ; Filtered_Image = Mean + K * ( Input_Image - Mean) (21) ; return, mean + (A-mean) * ( Z/(mean2*Sig^2 + Z) ) END ; Faster, but less accurate Lee filter. Only recomended when N < 7. Function Lee_filter_fast,A,N,SIG, $ DOUBLE=double COMPILE_OPT hidden on_error,2 Xmean = SMOOTH(a, N, /EDGE_TRUNCATE) ;Ensure floating Xmean2 = Xmean^2 AA = A - Xmean ; Here we VAR_X to be >= 0, because the variance can't be < 0. VAR_X=((SMOOTH(AA^2, N, /EDGE_TRUNCATE)+Xmean2)/(1+Sig^2)) - Xmean2 > 0 ; return value from equation 21,22 of Lee, 1986. ; K = ( VAR_X/(Xmean^2 * Sigma^2 + VAR_X) ) (22) ; Filtered_Image = Xmean + K * ( Input_Image - Xmean) (21) ; RETURN, Xmean + AA * (VAR_X / ( Xmean2 * Sig^2 + VAR_X)) END ;+ ; NAME: ; LEEFILT ; ; PURPOSE: ; Performs the Lee filter algorithm on an image array using a ; box of size 2N+1. This function can also be used on vectors. ; ; CATEGORY: ; E3 Smoothing (image). ; ; CALLING SEQUENCE: ; Result = LEEFILT(A [, N, Sig]) ; ; INPUTS: ; A: The input image array or one-dimensional vector. ; ; OPTIONAL INPUT PARAMETERS: ; N: The size of the filter box is 2N+1. The default value is 5. ; ; Sig: Estimate of the standard deviation. The default is 5. ; If Sig is negative the procedure requests a value to be ; entered, and displays the resulting image or vector. This ; cycle continues until a zero value of Sig is entered. ; ; KEYWORDS: ; ; EXACT: Use this keyword to use a more accurate, but much slower ; Lee filter algorithm. Recommended when N > ~7. ; ; DOUBLE = Set this keyword to force the computations to be done ; in double-precision arithmetic. ; ; OUTPUTS: ; The filtered image or vector is returned. ; ; COMMON BLOCKS: ; None. ; ; SIDE EFFECTS: ; Displays the filtered image in an IDL window using TVSCL or PLOT if ; Sig is negative. ; ; RESTRICTIONS: ; None. ; ; PROCEDURE: ; The LEE (Optical Engineering, Vol 25, No 5, Pg 636-643, May 1986) ; technique smooths additive image noise by generating statistics in ; a local neighborhood and comparing them to the expected values. ; ; MODIFICATION HISTORY: ; Written, 24-Nov-1982, by R. A. Howard, Naval Research Lab, ; Washington, DC 20375 ; Modified, 30-May-1996, SVP. Modifed to match 1986 algorithm & ; Added /EXACT at user suggestion. ; Added PLOT for vector support. ; CT, RSI, May 2000: Added double-precision support. ;- ; Function Leefilt, aIn, nIn, sigIn, $ DOUBLE=double, $ EXACT=EXACT ON_ERROR,2 ;Return to caller if an error occurs tnames = [SIZE(aIn,/TNAME), SIZE(nIn,/TNAME), SIZE(sigIn,/TNAME)] double = (N_ELEMENTS(double) GT 0) ? KEYWORD_SET(double) : $ MAX(tnames EQ 'DOUBLE') NP = N_params(0) IF np lt 3 THEN sigIn = 5 ;supply defaults sig = double ? DOUBLE(sigIn) : FLOAT(sigIn) IF np lt 2 THEN nIn = 5 dim = SIZE(aIn,/DIMENSIONS) n = FLOOR(nIn) IF ((n LT 1) OR (n GE MIN(dim)/2)) THEN MESSAGE, $ 'Filter width 2N+1 must be > 2 and smaller than array dimensions.' a = double ? DOUBLE(aIn) : FLOAT(aIn) pl = sig LE 0. ;true if interactive mode loop: IF pl THEN read,'Type in Sigma (0 to quit) : ',sig IF sig EQ 0 THEN goto,endp IF keyword_set(EXACT) THEN $ f=lee_filter_exact(a,2*n+1,Sig, DOUBLE=double) $ ELSE f=lee_filter_fast(a,2*n+1,Sig, DOUBLE=double) IF pl THEN BEGIN IF (size(f))[0] eq 1 then $ PLOT,f,xtitle='Element of A',ytitle='Value of A',$ title='Lee Filtered A, N ='+string(N,FORMAT='(I3)')+', Sigma = '+$ string(Sig,FORMAT='(F5.2)') else TVSCL,f GOTO , loop ENDIF endp: RETURN,f END