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

enter image description here

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 )

[m, n]=size(refIm);
[M,N] = meshgrid(1:m,1:n);
X = [M(:), N(:)];
[ap, bp, rhat]=lsa(angle(r));   %this is my slope a' and b'
[~,w] = max(r(:));
[del_hat2p, del_hat1p] = ind2sub(size(r),w);
R=exp(1j*theta); %E6
theta=atan2(imag(R),real(R)); %E9
[a, b, thetahat]=lsa(theta);
thetapp=theta-rhat;   %equation 14, I guess something amiss here
[app, bpp, thetapphat]=lsa(thetapp);
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);


function [ a, b, hat ] = lsa( theta )

[m, n]=size(theta);
[M,N] = meshgrid(1:m,1:n);
X = [M(:), N(:)];
B=regress(theta(:), X);


Please lighten me on this issue. Thank you so much!


1 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.

  • $\begingroup$ Hey I just saw your answer. I stopped looking at this paper for quite some time. Since now I have the response, would you mind we have a chat? I don't understand what you meant by the actual k. If I'm not wrong it's the frequencies. If so, it's very straight forward for 1D signal with sampling rate, but since this is an image without the sampling rate, would you mind explain how could I obtain it? $\endgroup$ Commented Jan 22, 2018 at 6:23
  • $\begingroup$ I managed to get the k values, however, when you said subtracting the phase due to the integer part (a′,b′), how can I obtain that values? If according to the paper, its says by using the POC, By using POC, I successfully get the values for the integer shift, but not the slopes of the integer shift. The paper is quite confusing. $\endgroup$ Commented Jan 24, 2018 at 13:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.