; $Id: leefilt.pro,v 1.3 1996/12/17 23:45:06 lubos Exp $

Function Lee_filter_exact,A,N,SIG
;
; Slow, but accurate Lee filter. Derived from
; Jong-Sen Lee, Optical Engineering 25(5), 636-643 (May 1986)
;
  Z=(AA=FLOAT(A))
  Delta=(T_Y=(T_X=(N-1)/2))
  SZ=size(AA) & N_SAMPLES=SZ[1]
;
  mean=smooth(AA,N,/EDGE_TRUNCATE) & mean2=mean^2.
;
  IF SZ[0] eq 1 THEN BEGIN
     FOR S=T_X,N_SAMPLES-T_X-1 DO $
        Z[S]=TOTAL((AA[S-Delta:S+Delta] - mean[S])^2.)
  ENDIF ELSE BEGIN
     N_LINES=SZ[2]
;
; Compute Variance of Z
;
    FOR L=T_Y, N_LINES-T_Y-1 DO $
      FOR S=T_X,N_SAMPLES-T_X-1 DO $
        Z[S,L]=TOTAL((AA[S-Delta:S+Delta,L-Delta:L+Delta] - mean[S,L])^2.)
  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)
;
  Z=( (TEMPORARY(Z) +mean2) /(1+Sig^2.)) - mean2 
;
; 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 + (AA-mean) * ( Z/(mean2*Sig^2. + Z) ) 
END

Function Lee_filter_fast,A,N,SIG
; Faster, but less accurate Lee filter
; Only recomended when N < 7
;
  AA=FLOAT(A)
;
  mean=SMOOTH(AA,N,/EDGE_TRUNCATE) & mean2=mean^2. 
  AA=AA-mean
;
 VAR_X=((SMOOTH(AA^2.,N,/EDGE_TRUNCATE)+mean2)/(1.+Sig^2.)) - mean2
;
; 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 + AA * (VAR_X / ( mean2 * Sig^2. + VAR_X))
;
END

Function Leefilt, A, N, Sig,EXACT=EXACT
;+
; 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.
; 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.
;-
;
ON_ERROR,2                      ;Return to caller if an error occurs
NP = N_params(0)
IF np lt 3 THEN Sig = 5. else Sig=FLOAT(Sig)	;supply defaults
IF np lt 2 THEN n = 5
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)  ELSE $
	   f=lee_filter_fast(a,2*n+1,Sig) 
	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