;------------------------------------------------------------------ ; LS.PRO ; ; written by Andreas Reigber ; 5.2.1997 ;------------------------------------------------------------------ ; Fast Least Square Phase Uwrapping (using cosine transform) ; pha = wrapped phase, uw = unwrapped phase return ; ; see: M.D. Pritt and J.S. Shipman: "Least-Square Two-dimensional ; Phase Unwrapping using FFT's", IEEE Trans. on Geoscience and Remote ; Sensing, Vol. 32, No. 3, pp. 706-708, 1994 ;------------------------------------------------------------------ ; This software has been released under the terms of the GNU Public ; license. See http://www.gnu.org/copyleft/gpl.html for details. ;------------------------------------------------------------------ pro ls,pha,uw stt=systime(1) s = size(pha) xsize=s(1) ysize=s(2) uw = atan(exp(complex(0,shift(pha,-1,0)-pha))) uw(xsize-1,*)=atan(exp(complex(0,pha(xsize-2,*)-pha(xsize-1,*)))) x2 = atan(exp(complex(0,pha-shift(pha,1,0)) )) x2(0,*)=atan(exp(complex(0,pha(0,*)-pha(1,*)))) uw = uw - x2 x2 = atan(exp(complex(0,shift(pha,0,-1)-pha) )) x2(*,ysize-1)=atan(exp(complex(0,pha(*,ysize-2)-pha(*,ysize-1)))) uw = uw + x2 x2 = atan(exp(complex(0,pha-shift(pha,0,1)) )) x2(*,0)=atan(exp(complex(0,pha(*,0)-pha(*,1)))) uw = uw - x2 for i=0,xsize-1 do uw(i,*) = float((fft([(uw(i,*))(*),(reverse((uw(i,*))(*)))(1:ysize-2)],-1))(0:ysize-1)) for i=0,ysize-1 do uw(*,i) = float((fft([(uw(*,i))(*),(reverse((uw(*,i))(*)))(1:xsize-2)],-1))(0:xsize-1)) index = findgen(xsize)*!pi/(xsize-1) x2 = cos(index) # (fltarr(ysize) + 1) index = findgen(ysize)*!pi/(ysize-1) x2 = x2 + (fltarr(xsize) + 1) # cos(index) x2(0,0) = x2(0,0) + 0.00001 uw = uw / (2*(x2 - 2)) for i=0,xsize-1 do uw(i,*) = float((fft([(uw(i,*))(*),(reverse((uw(i,*))(*)))(1:ysize-2)],1))(0:ysize-1)) for i=0,ysize-1 do uw(*,i) = float((fft([(uw(*,i))(*),(reverse((uw(*,i))(*)))(1:xsize-2)],1))(0:xsize-1)) end