;------------------------------------------------------------------
; MLSW.PRO
; Weighted Multigrid Least Square Phase Uwrapping
;
; (c)1997 by Andreas Reigber
;------------------------------------------------------------------
; Usage:
; res=mls(phase,weights)
;      phase  = 2D floating point phase array to be unwrapped
;      weight = 2D floating point array containing the weights
;      res  = unwrapped phase
;
; see: M.D. Pritt: "Phase Unwrapping by Means of Multigrid Techniques
; for Interferometric SAR", IEEE Trans. on Geoscience and Remote
; Sensing, Vol. 34, No. 3, pp. 728-738,. 1996
;-----------------------------------------------------------------
; This software has been released under the terms of the GNU Public
; license. See http://www.gnu.org/copyleft/gpl.html for details.
;------------------------------------------------------------------
function mlsw,pha,gew
   sizx = (size(pha))[1]
   sizy = (size(pha))[2]
	dxdyw,pha,gew,dx,dy,wen
	dummy=pha-pha	
	return,wfmg(dummy,dx,dy,wen)
end
;---------------------------------------------------------------------
function relaxw,pphi,dx,dy,wen,nr
   dummy=pphi
   sizx = (size(pphi))(1)
   sizy = (size(pphi))(2) 
   res = pphi  
   for i=1,nr do begin
   	dummy = shift(res - dx,-1,0)
   	dummy(sizx-1,*) = res(sizx-2,*) + dx(sizx-1,*) 
   	res2 =  shift(wen,-1,0) * dummy

   	dummy = shift(res - dy,0,-1)
   	dummy(*,sizy-1) = res(*,sizy-2) + dy(*,sizy-1) 
   	res2 =  res2 + shift(wen,0,-1) * dummy
   	
   	dummy = shift(res,0,1) + shift(res,1,0)
   	dummy(0,*) = dummy(0,*) - res(sizx-1,*) + res(1,*)
   	dummy(*,0) = dummy(*,0) - res(*,sizy-1) + res(*,1)
   	res2 = res2 + wen*(dummy + dx + dy)
   	res  = res2 / (wen + wen + shift(wen,-1,0) + shift(wen,0,-1))
	endfor
return,res
end

;--------------------------------------------------------------------
function wmv,phi_res,dxn,dyn,wen
   sizx = (size(phi_res))(1)
   sizy = (size(phi_res))(2)

   phi_res = relaxw(phi_res,dxn,dyn,wen,2)
   if sizx eq 4 then goto,back2
   dummy = phi_res - shift(phi_res,1,0)
   dummy(0,*) = - dummy(1,*)      
   dxn2 = 2*congrid(smooth(smooth(dxn-dummy,3),3), sizx/2, sizy/2)
   dummy = phi_res - shift(phi_res,0,1)
   dummy(*,0) = - dummy(*,1)         
   dyn2 = 2*congrid(smooth(smooth(dyn-dummy,3),3), sizx/2, sizy/2)
   dyn2(*,0) = -dyn2(*,1)
   dxn2(0,*) = -dxn2(1,*)
   
   dummy=fltarr(sizx,sizy)
   dummy(where(wen gt 0.0025)) = 1
   wen2=rebin(wen,sizx/2,sizy/2)
   dd=where(rebin(dummy,sizx/2,sizy/2) lt 0.6,anz)
   if anz gt 0 then wen2(dd) = 0.0025
   dummy = 0
   dd = 0
	   
	phin2 = fltarr(sizx/2, sizy/2)
	ret = wmv(phin2,dxn2,dyn2,wen2)
	phi_res = phi_res + congrid(ret, sizx, sizy) 
   back2:
 	return,relaxw(phi_res,dxn,dyn,wen,2)
end

;----------------------------------------------------------------------
function wfmg,phin,dxn,dyn,wen
   sizx = (size(phin))(1)
   sizy = (size(phin))(2)
   if sizx eq 4 then goto,back
   dummy = phin - shift(phin,1,0)
   dummy(0,*) = - dummy(1,*)   
   dxn2 = 2*congrid(smooth(smooth(dxn-dummy,3),3), sizx/2, sizy/2)
   dummy = phin - shift(phin,0,1)
   dummy(*,0) = - dummy(*,0)   
   dyn2 = 2*congrid(smooth(smooth(dyn-dummy,3),3), sizx/2, sizy/2)
   dyn2(*,0) = - dyn2(*,1)
   dxn2(0,*) = - dxn2(1,*)
   
   dummy=fltarr(sizx,sizy)
   dummy(where(wen gt 0.0025)) = 1
   wen2=rebin(wen,sizx/2,sizy/2)
   dd=where(rebin(dummy,sizx/2,sizy/2) lt 0.6,anz)
   if anz gt 0 then wen2(dd) = 0.0025
   dummy = 0
   dd = 0
     
	phin2 = fltarr(sizx/2, sizy/2)
	ret = wfmg(phin2,dxn2,dyn2,wen2)
	phin = phin + congrid(ret, sizx, sizy) 
   back:
 	return,wmv(phin,dxn,dyn,wen)
end

;---------------------------------------------------------------------
pro dxdyw,pha,gew,dx,dy,wen
   sizx = (size(pha))(1)
   sizy = (size(pha))(2)
	dx = atan(exp(complex(0,pha - shift(pha,1,0))))
	dy = atan(exp(complex(0,pha - shift(pha,0,1))))
	dx(0,*) = - dx(1,*)
	dy(*,0) = - dy(*,1)	
	wen = gew
	dummy = shift(gew,1,0)
	a = where(dummy lt wen,len)
	if len ne 0 then wen(a) = dummy(a)
	dummy = shift(gew,0,1)
	a = where(dummy lt wen,len)
	if len ne 0 then wen(a) = dummy(a)	
	wen = wen^2
	a = where(wen lt 0.0025,len)
	if len ne 0 then wen(a) = 0.0025	
end