I tried to reproduce the result from this paper, However I'm stuck at the Section III of the paper which involve integer shift. Two images which their phase correlation $\theta(k_1,k_2)$ was calculated at pixel-accuracy(integer shift). I get their phase component as shown(from -$\pi$ to $\pi$). It can be represented by $\theta(k_1,k_2)=\theta_1(k_1,k_2)-\theta_2(k_1,k_2)$. Since the two images shifted in parallel, it can be approximated as $\theta(k_1,k_2) \approx ak_1+bk_2$ where $a$ and $b$ is the slopes obtained through least square method, $k_1$ and $k_2$ are the value $MxN$ of the image. Plotting $\theta(k_1,k_2)$:
And to get phase component with only sub-pixel accuracy(decimal shift), the phase component is subtracted with integer shift. The equation is $\theta''(k_1,k_2)=\theta(k_1,k_2)-(a'k_1+b'k_2)$
However, as you can see from the second figure, the value of $\theta''(k_1,k_2)$ ranges wider than -$\pi$ to $\pi$(as compared to Figure 4 in the paper), which I suspect the slope($a'$,$b'$) was obtained inaccurately. The paragraph above the Equation (14) mentioned that they used conventional POC to obtain the slopes. However I'm not sure did I understood it correctly.
The code I implemented:
function [ output_args ] = phasecorrlsa( refIm, shifIm ) F=(fft2(double(refIm))); G=(fft2(double(shifIm))); [m, n]=size(refIm); [M,N] = meshgrid(1:m,1:n); X = [M(:), N(:)]; R=(F.*conj(G))./abs(F.*conj(G)); r=(ifft2(R)); % surf(r) [apx, bpx, rhat]=lsa(angle(r)); [~,w] = max(r(:)); [del_hat2p, del_hat1p] = ind2sub(size(r),w); del_hat2p=del_hat2p-1; del_hat1p=del_hat1p-1; ap=(del_hat1p*2*pi)/m; %THE SLOPES,obtained through Eq. (12) bp=(del_hat2p*2*pi)/n; theta=angle(F)-angle(G); R=exp(1j*theta); %E6 theta=atan2(imag(R),real(R)); %E9 [a, b, thetahat]=lsa(theta); del_hat1_ts=(m/(2*pi))*a; del_hat2_ts=(n/(2*pi))*b; thetapp=theta-(ap*M+bp*N); figure; surf(theta); figure; surf(thetapp); [app, bpp, thetapphat]=lsa(thetapp); del_hat1=(m/(2*pi))*app+del_hat1p; del_hat2=(n/(2*pi))*bpp+del_hat2p; if del_hat1>n/2, del_hat1=del_hat1p-m; end if del_hat2>m/2, del_hat2=del_hat1p-n; end output_args=[del_hat1 del_hat2]; end function [ a, b, hat ] = lsa( theta ) [m, n]=size(theta); [M,N] = meshgrid(1:m,1:n); X = [M(:), N(:)]; B=regress(theta(:), X); a=B(1); b=B(2); hat=reshape(X*B,m,n); end
Anyone could lighten me what went wrong with my code? Thanks in advance!