;------------------------------------------------------------------
; MLS.PRO
; Unweighted Multigrid Least Square Phase Uwrapping
;
; (c)1997 by Andreas Reigber
;------------------------------------------------------------------
; Usage:
; res=mls(pha)
;      pha  = 2D floating point phase array to be unwrapped
;      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 mls,pha
   sizx = (size(pha))[1]
   sizy = (size(pha))[2]
	dxdy,pha,dx,dy
	dummy=pha-pha	
	return,ufmg(dummy,dx,dy)
end
;---------------------------------------------------------------------

function relax,pphi,dx,dy,nr
   dummy=pphi
   sizx = (size(pphi))(1)
   sizy = (size(pphi))(2) 
   res = pphi  
   for i=1,nr do begin
		res2 = shift(res,0,1)+shift(res,0,-1)+shift(res,1,0)+shift(res,-1,0)
		res2(0,*) = res2(0,*) - res(sizx-1,*) + res(1,*)
		res2(sizx-1,*) = res2(sizx-1,*) - res(0,*) + res(sizx-2,*)
		res2(*,0) = res2(*,0) - res(*,sizy-1) + res(*,1)
		res2(*,sizy-1) = res2(*,sizy-1) - res(*,0) + res(*,sizy-2)
		dummy = shift(dx,-1,0) - dx
		dummy(sizx-1,*) = -2 * dx(sizx-1,*)
		res = res2 - dummy
		dummy = shift(dy,0,-1) - dy
		dummy(*,sizy-1) = -2 * dy(*,sizy-1)
		res = (res - dummy)/4
	endfor
return,res
end

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

   phi_res = relax(phi_res,dxn,dyn,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,*)
	phin2 = fltarr(sizx/2, sizy/2)
	ret = umv(phin2,dxn2,dyn2)
	phi_res = phi_res + congrid(ret, sizx, sizy) 
   back2:
 	return,relax(phi_res,dxn,dyn,2)
end

;----------------------------------------------------------------------
function ufmg,phin,dxn,dyn
   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,*)
	phin2 = fltarr(sizx/2, sizy/2)
	ret = ufmg(phin2,dxn2,dyn2)
	phin = phin + congrid(ret, sizx, sizy) 
   back:
 	return,umv(phin,dxn,dyn)
end

;---------------------------------------------------------------------
pro dxdy,pha,dx,dy
   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)
end