I try to estimates shift estimation directly in phase region, by following the proposed method in this Sub-pixel Shift Estimation of Image based on the Least Squares Approximation in Phase Region by Fujimoto, Fujisawa and Ikehara (Proceedings of 26th European Signal Processing Conference, EUSIPCO '16, pp. 91–95. IEEE, 2016 PDF). Here is the flow chart described in the paper
If you're still with me, my issue raised when the author try to proceed with the subtraction of the integer shift with the phase difference $θ(k1,k2)$ (Equation (14) to smooth the phase difference. The authors however did mention that they obtained the slope $(a′,b′)$ of the integer shift using conventional phase-only correlation(POC) (Section III, paragraph above Equation (14)).
How can this step can be done? Since the phase difference is of a $60\times60$ matrix (assume the image dimension is so), while the integer shift consists only TWO values. How exactly were the slopes obtained?
Full matlab 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));
[ap, bp, rhat]=lsa(angle(r)); %this is my slope a' and b'
[~,w] = max(r(:));
[del_hat2p, del_hat1p] = ind2sub(size(r),w);
del_hat2p=del_hat2p-1;
del_hat1p=del_hat1p-1;
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-rhat; %equation 14, I guess something amiss here
%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=struct('a',del_hat1, 'b', 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
Please lighten me on this issue. Thank you so much!