Estimates sub-pixel shift directly in phase region
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.