
\def\a{{a}}
\section{Preliminaries on sub-pixel interpolation}\label{sec:theoretical_results}
This section proves a  discrete correlation formula which is faithful to the continuous image interpolates.  Thanks to it, an accurate subpixel matching becomes possible.
Without loss of generality, all considered  images $u$, $u_1$, etc. are defined on a square $[0,\a]^2$ and are supposed to be square integrable. Thus, the Fourier series decomposition applies
\begin{equation}\label{Fourierseries}u(x,y)\!\!=\!\!\! \sum_{k,l\in{\Z}} \tilde u_{k,l}e^{\frac{2i\pi(kx+ly)}\a},\end{equation}
where the $\tilde u_{k,l}$ are  the Fourier series coefficients (or shortly the Fourier coefficients) of $u$. By the classic Fourier series isometry, for any  two square integrable functions $u(\bx)$ and $v(\bx)$ on $[0,a]^2$,

\begin{equation}\label{Fourierseriesisometry}
\int_{[0,a]^2}u(\bx)\overline{v}(\bx)d\bx=a^2\sum_{k,l\in{\Z}}\tilde u_{k,l}\overline{\tilde v_{k,l}}.
\end{equation}

 The digital images  are usually given by their $N^2$ samples $u(\bmm)$ for $\bmm$ in the grid
 $$\Z^1_\a=[0,\a]^2\cap ((\frac \a{2N},\frac \a{2N})+ \frac \a N\Z^2).$$ Similarly,   the  over-sampling grid with four times more samples is denoted by
$$\Z^{\nicefrac {1}{2}}_\a=[0,\a]^2\cap ((\frac \a{4N},\frac \a{4N})+\frac \a {2N}\Z^2).$$

$N$  is always an even integer. In all that follows we shall assume that the images obtained by a stereo vision system are band-limited. This assumption is classical and realistic, the aliasing in good quality CCD cameras being moderate.
  As classical in image processing, under the (forced) $a$-periodicity assumption a
band-limited image  becomes a trigonometric polynomial. This periodicity assumption is not natural, but it only entails a minor drawback, namely a small
distortion near the boundary of the image domain $[0,\a]^2$.  The payoff for the {\it band-limited + periodic assumption} is that the image can be interpolated, and its Fourier coefficients computed from discrete samples.
Indeed, given $N^2$ samples $u_\bmm$ for $\bmm$ in $\Z^1_\a$, there is a unique  trigonometric polynomial in the form \begin{equation}\label{interpolationFFT}u(x,y)\!\!=\!\!\! \sum_{k,l=-N/2}^{N/2-1} \tilde u_{k,l}e^{\frac{2i\pi(kx+ly)}\a}\end{equation}
such that $u(\bmm)=u_{\bmm}$.
We shall call such polynomials {\it $N$-degree trigonometric polynomials}.    The coefficients $\tilde u_{k,l}$ are the {\it Fourier coefficients} of $u$ in the Fourier basis $e^{\frac{2i\pi(kx+ly)}{a}}$, $k,l\in {\Z}$. The map $u_\bmm\to u_{k,l}$ is nothing but the  {\it 2D Discrete Fourier Transform (DFT)}, and
  the map $(u_{\bmm})\to N(\tilde u_{k,l})$ is  an isometry from ${\bf C}^{N^2}$ to itself. The function $u(x,y)$ is therefore usually called the {\it DFT interpolate of the samples $u_\bmm$.}   In consequence, there is an isometry between the set of $N$-degree trigonometric polynomials endowed with the $L^2([0,a]^2)$ norm, and ${\bf C}^{N^2}$ endowed with the usual Euclidean norm:

\begin{equation}\label{isometries}
\int_{[0,a]^2}|u(x,y)|^2=\a^2\sum_{k,l=-N/2}^{N/2-1} |\tilde u_{k,l}|^2=\frac{\a^2}{N^2}\sum_{\bmm\in\Z^1_a} |u(\by+\bmm)|^2,
\end{equation}
where the  $N^2$ samples grid can have an arbitrary origin $\by$.
If $u(\bx)$ and $v(\bx)$ are two degree $N$  trigonometric polynomials, we therefore  also have
\begin{equation}\label{isometries2}
\int_{[0,\a]^2}\!\!\!\!\!u(\bx)\overline{v}(\bx)=\!\a^2\sum_{k,l=-N/2}^{N/2-1} \tilde u_{k,l}\overline{\tilde v_{k,l}}=\!\!\frac{\a^2}{N^2}\sum_{\bmm\in\Z^1_\a} \!\!u(\by+\bmm)\overline{v(\by+\bmm)}.
\end{equation}
where $\overline{v}$ is the complex conjugate of $v$. Taking four  times more samples, it  follows from \eqref{isometries2} that
\begin{equation}\label{isometries3}
\int_{[0,\a]^2}\!\!\!\!\!u(\bx)\overline{v}(\bx)=\!\a^2\sum_{k,l=-N}^{N-1} \tilde u_{k,l}\overline{\tilde v_{k,l}}=\!\frac {\a^2}{4N^2}\sum_{\bmm\in \Z^{\nicefrac {1}{2}}_\a} \!u(\bmm)\overline{v(\bmm)}.
\end{equation}
% AA v7
which is also valid if $u(\bx)$ and $v(\bx)$ are up to $2N$-degree trigonometric polynomials in $\bx$.

%-------------%
This last fact has a first important consequence in block-matching. Consider two images $u_1(\bx)$ and $u_2(\bx)$ on $[0,a]^2$ and a window function $\varphi(\bx)$. Block matching is the search for a value of $\mu$ minimizing the continuous quadratic distance
\begin{equation}\label{definitiondelacible}
e_{\bx_0}(\mu):=\int_{[0,a]^2}\varphi(\bx-\bx_0)\big(u_1(\bx)-u_2(\bx+(\mu,0))\big)^2d\bx.
\end{equation}


\begin{proposition}{\bf (Equality of the discrete and the continuous quadratic distance)}\label{prop:discrete-equals-continuous}
Let  $u_1(\bx)$ and $u_2(\bx)$ be two $N$-degree trigonometric polynomials on $[0,a]^2$ and  let $\varphi(\bx)$ be a window function which we assume to be a $2N$-degree trigonometric polynomial. Then

\begin{equation}
\label{identitedisparitediscretecontinue}e_{\bx_0}(\mu)= e^d_{\bx_0}(\mu),\; \hbox{ where}
\end{equation}
\begin{equation}\label{discretisationdeladistance}
e^d_{\bx_0}(\mu):=\!\frac{\a^2}{4N^2}\sum_{ \bmm \in \Z^{\nicefrac{1}{2}}_\a} \!\varphi(\bmm-\bx_0)\big(u_1(\bmm)-u_2(\bmm+(\mu,0))\big)^2.
\end{equation}
\end{proposition}

The proof follows from \eqref{isometries3}. Indeed,  $\big(u_1(\bx)-u_2(\bx+(\mu,0))\big)^2$  and $\varphi(\bx-\bx_0)$ are both $2N$-degree trigonometric polynomials in $\bx$,
% AA v7
so according to \eqref{isometries3} the discrete scalar product defining $e^d_{\bx_0}(\mu)$ equals the continues scalar product defining $e_{\bx_0}(\mu)$.
% End AA
Thus  {\it the continuous block distance is a finite sum of discrete samples!}


The block distance function $\mu\to e_{\bx_0}(\mu)$, whose minimization is our main objective here,  is also easily  sampled. By  \eqref{discretisationdeladistance}  it is a 2$N$-degree trigonometric polynomial with respect to $\mu$. This proves:

\begin{proposition}[Sub-pixel correlation requires $\times$2 zoom]\label{subpixelcorrelation2zoom}
Let  $u_1(\bx)$ and $u_2(\bx)$ be two $N$-degree trigonometric polynomials. Then the quadratic distance $e^d_{\bx_0}(\mu)$ is well-sampled provided it has at least $2N$ successive samples. Thus the computation of $e^d_{\bx_0}(\mu)$ at half samples $\mu\in\frac{a\Z}2$ (via zero padding) allows the exact reconstruction of $e^d_{\bx_0}(\mu)$ for any real $\mu$ by DFT interpolation.
\end{proposition}

Remark that the last proposition does not require any assumption on the window function $\varphi(\bx)$.
Prop. \ref{subpixelcorrelation2zoom}, which opens the way to rigorous block matching with sub-pixel accuracy, has been noticed (but not implemented) in \cite{Szeliski04}. It is also  used in the MARC method \cite{Giros04} used by the French space agency (CNES).  The above simple  proof of Prop. \ref{subpixelcorrelation2zoom} is new.\\


Sub-pixel block matching  will require to interpolate the noisy images. Thus, following Shannon's classical observation,  the noise itself must also be interpolated as a band-limited function. In the periodic framework it therefore becomes  a trigonometric polynomial.
Assume that $(n_{\bmm}), \bmm \in \Z^1_\a$ are $N^2$ independent   $\mathcal N(0,\sigma^2)$ noise samples.
This amounts to say that $(n_\bmm)$ is a Gaussian vector. Since the DFT is an  isometry, the noise Fourier   coefficients $N(\tilde {n}_{\bk})$ also form a Gaussian vector with diagonal covariance matrix $\sigma^{2} Id$.
By \eqref{isometries2}, the mapping
$(n_\bmm)_{\bmm\in \Z^1_\a}\to (n(\bx+\bmm))_{\bmm\in\Z^1_\a}$ is an isometry  from ${\bf C}^{N^2}$ to itself.
It follows that $n(\bx)$ is $\mathcal N(0,\sigma^2)$ for every $\bx$.


One can also estimate $\Var(n_x(\bx))$, where $n_x(\bx)=\frac{\partial n}{\partial x}(x,y).$
$$\Var(n_x(\bx))=\Var\left(\sum_{k,l=-N/2}^{N/2-1}  \tilde n_{k,l}\frac{2ik\pi}\a e^{2i\frac{k\pi x+l\pi y}\a}\right)$$
$$=\frac{4\pi^2\sigma^2 N}{N^2\a^2}\sum_{k=-N/2}^{N/2-1}k^2 \simeq\frac{4\pi^2\sigma^2}{\a^2N}\frac {N^3}{12}=\frac{\pi^2N^2}{3\a^2}\sigma^2.$$
Since $n(\bx)$ is a normal law, $n(\bx)^2$ is a $\chi^2$  law of order 1. Thus its variance is $2\sigma^4$.
 Finally we shall need to evaluate $\Var(n_1(\bx)n_2(\bx))$, where $n_i$ are two independent  interpolated white noises of the above kind.  Thus $n_1(\bx)n_2(\bx)$ is the product of two normal laws. The expectation of the product is zero and the  variance is therefore $\Var(n_1n_2)= {\bf E }(n_1n_2)^2={\bf E}n_1^2 {\bf E}n_2^2=({\bf E}n^2)^2=\Var(n)^2=\sigma^4$. In summary:

\begin{lemma} \label{noiseinterpolation}  Let $(n_\bmm)_{\bmm\in\Z^1_\a}$ be $N^2$ independent white Gaussian noise samples  with variance $\sigma^2$.  Then the DFT interpolate $n(\bx)$ on $[0,a]^2$ is  $\mathcal N(0,\sigma^2)$ for every $\bx$. If $n_1$ and $n_2$ are two  independent noises like $n$, one  has
\begin{align}\label{variousvariances}
\Var (n^2(\bx))\!\! &=\!\!2\sigma^4,\;\\
\Var((n)_x(\bx))\!\!&\simeq\!\! \frac{\pi^2N^2}{3a^2}\sigma^2,\;\\
\Var (n_1(\bx)n_2(\bx))\!\!&=\!\!\sigma^4.
\end{align}
\end{lemma}
 