# Get slopes of integer shift of an image in radians

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!