Estimates sub-pixel shift directly in phase region

Question

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!

Answer 1

I am unable to comment on your Matlab code, but eq. (14) seems straightforward to me. You have a phase difference field $\theta$ which depends on the two spatial wavenumber components $k_1$ and $k_2$, and which formes a wrapped plane. The shift you seek is the slope of this plane, expressed as two scalar components $a$ and $b$. Since phase wraps, the slopes are decomposed into integer and fractional parts, $a = a' + a''$ etc. In eq. (14) you obtain the residual phase field by subtracting the phase due to the integer part $(a', b')$. You need the actual $k$ values when you compute $(a'k_1 + b'k_2)$, not just the $k$ indices. Computing the $k$ values is similar to making a frequency axis for a 1D Fourier spectrum analysis.