;-----------------------------------------------------------------
; lee_pol_ref.pro
;
; 06/2001 by Andreas Reigber
;------------------------------------------------------------------
; Vectorial refined Lee filter of multiplicative image noise 
; using local  statistics. Used for complex polarimetric SAR images. 
; Better preservation of image structures than the 
; original Lee filter by the use of directional windows. 
;
; J.S.Lee et al.: "Speckle Reduction in Multipolarization,
; Multifrequency SAR Imagery", Trans. on Geoscience and Remote
; Sensing, Vol. 29, No. 4, pp. 535-544, 1991
;
; J.S.Lee: "Refined filtering of image noise using local statistics"
; Comp. Graph. Image Process. Vol. 15, pp. 380-389, 1981
;
; Usage:
;    res = lee_pol(image,size)
;   
;    image = 3D complex array containing the polarimetric image
;            the polarizations are interleaved over the third
;            dimension, i.e. image has size (xdim,ydim,3).
;    size  = size of the filter window
;
;    res   = Resulting filtered covariance matrix (xdim,ydim,3,3)
;
; Keywords:
;
;    looks     = Number of looks (default = 1)
;    threshold = Threshold for switching to rectangular boxes
;
;------------------------------------------------------------------
; This software has been released under the terms of the GNU Public
; license. See http://www.gnu.org/copyleft/gpl.html for details.
;------------------------------------------------------------------
function lee_pol_ref,arr,smm,LOOKS=looks,THRESHOLD=threshold
	if not keyword_set(looks) then looks=1
	if not keyword_set(threshold) then threshold=0.5

	siz    = size(arr)	
	anz_rg = siz[1]
	anz_az = siz[2]

	delta = (smm-1)/2
	ol    = delta / 2 + 1
	thr   = threshold
	sig2  = 1.0/looks
	sfak  = 1.0+sig2

	amp   = total(abs(arr),3)
	out   = complexarr(anz_rg,anz_az,3,3)

	fbox  = mbox(smm,anz_rg)
	sbox  = intarr(9,delta^2)
	for i=0,2 do begin
		sbox[i*3,*]   = mbox(delta,anz_rg)+ol*(i-1)*anz_rg-ol
		sbox[i*3+1,*] = mbox(delta,anz_rg)+ol*(i-1)*anz_rg
		sbox[i*3+2,*] = mbox(delta,anz_rg)+ol*(i-1)*anz_rg+ol
	endfor
	dummy = sbox[4,*]
	sbox[4,*] = sbox[8,*]
	sbox[8,*] = dummy
	
	sbox[0,*]   = mbox(delta,anz_rg)-ol
	sbox[1,*]   = mbox(delta,anz_rg)-ol*anz_rg-ol
	sbox[2,*]   = mbox(delta,anz_rg)-ol*anz_rg
	sbox[3,*]   = mbox(delta,anz_rg)-ol*anz_rg+ol
	sbox[4,*]   = mbox(delta,anz_rg)+ol
	sbox[5,*]   = mbox(delta,anz_rg)+ol*anz_rg-ol
	sbox[6,*]   = mbox(delta,anz_rg)+ol*anz_rg
	sbox[7,*]   = mbox(delta,anz_rg)+ol*anz_rg+ol
	sbox[8,*]   = mbox(delta,anz_rg)
	
	cbox = intarr(8,smm*(delta+1))
	for i=0,smm-1 do begin
		cbox[0,i*(delta+1):(i+1)*(delta+1)-1] = indgen(delta+1)+(i-delta)*anz_rg
		cbox[2,i*(delta+1):(i+1)*(delta+1)-1] = indgen(delta+1)*anz_rg+(i-delta)
		cbox[4,i*(delta+1):(i+1)*(delta+1)-1] = -indgen(delta+1)+(i-delta)*anz_rg
		cbox[6,i*(delta+1):(i+1)*(delta+1)-1] = -indgen(delta+1)*anz_rg+(i-delta)  
	endfor
	pos = 0
	for i=0,smm-1 do begin		
		cbox[1,pos:pos+i] = (delta-indgen(i+1))*anz_rg  + (i-delta)
		cbox[3,pos:pos+i] = (delta-indgen(i+1))*anz_rg  - (i-delta)
		cbox[5,pos:pos+i] = (-delta+indgen(i+1))*anz_rg + (i-delta)
		cbox[7,pos:pos+i] = (-delta+indgen(i+1))*anz_rg - (i-delta)
		pos = pos + (i+1)
	endfor
	pol  = [0,anz_rg*anz_az,2*anz_rg*anz_az]
;------------------------------------------------------	
	for i=delta,anz_rg-delta-1 do begin
		for j=delta,anz_az-delta-1 do begin
			pos  = i+j*anz_rg
			stat = total(amp[pos+sbox],2)/9
			grad = max(abs(stat/stat[8]-1.0),dir)
			if grad gt thr then win = reform(cbox[dir,*]) else win=fbox
			win2  = [[win],[win+anz_rg*anz_az],[win+2*anz_rg*anz_az]]
         k3    = arr[pos+win2]
			meanc = (conj(k3) ## transpose(k3))/n_elements(win)
			mean2 = (abs(meanc[0])+abs(meanc[4])+abs(meanc[8]))	
			z     = total((total(abs(k3),2)-sqrt(mean2))^2)/(n_elements(win)-1)
			z     = ((z+mean2)/sfak)-mean2 > 0
			fak   = z/(mean2*sig2+z)
			k3    = arr[pos+pol]
			cov   = adj(k3) ## k3
			out[i,j,*,*]   = meanc + (cov - meanc)*(z/(mean2*sig2+z))
		endfor
	endfor
	return,out
end