\documentclass{article}
\usepackage{amsmath, amsthm, amssymb}
\usepackage{latexsym}
\usepackage{doc}
\usepackage{exscale}
\usepackage{fontenc}
\usepackage{graphicx}
\usepackage{epsfig}


%\usepackage{multirow} % the packet is in the current directory of report.tex

%--------------------------------------------------------------
% Page Layout
%--------------------------------------------------------------
\setlength{\textwidth}{15cm}
\setlength{\columnsep}{0.8cm}
%\setlength{\parindent}{4ex}
\setlength{\parskip}{0.1cm}
\setlength{\oddsidemargin}{0.5cm}
\setlength{\evensidemargin}{0.5cm}
\setlength{\topmargin}{2.5cm}
% \setlength{\textheight}{25.5cm}\author{Guoshen YU}
\newtheorem{comment}{Comment}
%\newtheorem{assumption}{SIFT assumption}
\newtheorem{question}{QUESTION}
\newtheorem{theorem}{Theorem}
\newtheorem{corollary}{Corollary}
\newtheorem{definition}{Definition}
\newtheorem{lemma}{Lemma}
\newtheorem{remark}{Remark}
\newtheorem{Guoshen}{Guoshen:}
\newtheorem{Jean-Michel}{Jean-Michel}
\def\enddemo{\hfill $\circ$ \newline\newline}


\newtheorem{proposition}{Proposition}
\newtheorem{assumption}{Assumption}
\def\Laplacian{\Delta}
\def\bA{{\bf A}}
\def\bB{{\bf B}}
\def\bE{{\bf E}}
\def\bS{{\bf S}}
\def\bL{{\bf L}}
\def\bG{{\bf G}}
\def\bR{{\bf R}}
\def\bH{{\bf H}}
\def\bT{{\bf T}}
\def\bI{{\bf I}}
\def\RR{\mathbb R}

\def\bu{{\bf u}}
\def\bv{{\bf v}}
\def\SA{{\bf S}}
\def\bx{{\bf x}}
\def\by{{\bf y}}
\def\bw{{\bf w}}
\def\relative{\mathbb Z} \pagestyle{headings}
\def\bc{{\mathbf c}}
\def\bB{{\mathbf B}}
\def\bq{{\mathbf q}}
\def\bd{{\bf d}}
\def\bn{{\bf n}}
\def\bm{{\bf m}}
\def\bk{{\bf k}}
\def\SAD{\mathcal S}
\def\ft2{{\bf c}}
\def\sft3{{\bf d}}

\newcommand{\nada}[1]{}

\def\real{\bf R}
\newcommand{\Var}{\operatorname{Var}}
\def\indic{{\bf 1}}
\def\relative{{\bf Z}}



\newcommand{\sinc}{{\ensuremath{\operatorname{sinc}}}}
\newcommand{\half}{{\ensuremath{\nicefrac{1}{2}}}}
\newcommand{\R}{\mathbb{R}} %real numbers
\newcommand{\Z}{\mathbb{Z}} %integer numbers
\newcommand{\w}{{\varphi_{q}}} % correlation window around x_0
\newcommand{\scalar}[2]{{\left\langle\,#1,\,#2\right\rangle}} % scalar product
\newcommand{\scalarX}[3]{{\scalar{#1}{#2}_{#3}}} % weighted scalar product
\newcommand{\scalarphi}[2]{\scalarX{#1}{#2}{\w}} % \varphi_x0 weighted scalar product
\newcommand{\normX}[2]{\| #1 \|_{#2}}
\newcommand{\normphi}[1]{\| #1 \|_\w}
\newcommand{\Proba}{\mathbb{P}}

\begin{document}


\title{Is SIFT Scale Invariant?}
\author{Jean-Michel Morel \and  Guoshen Yu}

\graphicspath{{figures/}}


\maketitle
\begin{abstract}
This note is devoted to a mathematical
exploration of whether  Lowe's {\em Scale-Invariant Feature Transform} (SIFT)
\cite{lowe04:_distin}, a very successful image matching method, is
similarity invariant as claimed. It is proved that the method is scale invariant only if the initial image blurs is exactly guessed. Yet, even a large error on the initial blur is quickly attenuated by this multiscale method, when the scale of analysis increases. In consequence, its
scale invariance is almost perfect. The mathematical arguments are
given under the assumption that the Gaussian smoothing  performed by SIFT gives
an aliasing free sampling of the image evolution. The validity of this main assumption is
confirmed by a rigorous experimental procedure, and by a mathematical
proof. These results explain why SIFT outperforms all other image
feature extraction methods when it comes to scale invariance.
\end{abstract}

%\tableofcontents


\section{Introduction}
\label{sec:intro}


Image comparison is a fundamental step in many computer vision and image processing applications. A typical image
matching method  first detects points of interest, then selects a region around each point, and finally associates with
each region a descriptor. Correspondences between two images can then be established by matching the descriptors  of
both images. The images under comparison may have been taken under arbitrary viewpoints. Therefore the invariance to the viewpoint is crucial in
image comparison.  Many variations exist on the computation of invariant interest points,
following the pioneering work of Harris and Stephens~\cite{harris88combined}. The Harris-Laplace and Hessian-Laplace
region detectors~\cite{MikolajczykS01, MikolajczykS04} are invariant to rotation and scale changes. Some moment-based
region detectors~\cite{DBLP:conf/eccv/LindebergG94, baumberg00:_reliab} including the Harris-Affine and Hessian-Affine
region detectors~\cite{MikolajczykS02, MikolajczykS04}, an edge-based region detector~\cite{TuytelaarsG04}, an
intensity- based region detector~\cite{TuytelaarsG04}, an entropy-based region detector~\cite{KadirZB04}, and two
independently developed level line-based region detectors MSER (``maximally stable extremal
region'')~\cite{matas02mser} and LLD (``level line descriptor'')~\cite{ muse03unsupervised, muse04:_acontrario}
are designed to be invariant to affine transformations. These two methods stem from the  Monasse image registration
method \cite{monasse99contrast} that used well contrasted extremal regions to register images. MSER is the most
efficient one and has shown better performance than other affine invariant detectors~\cite{mikolajczyk05}. However, as
pointed out in ~\cite{lowe04:_distin}, none of these detectors is actually fully affine invariant: All of them start with
initial feature scales and locations selected in a non-affine invariant manner.  The difficulty comes from the scale
change from an image to another: This change of scale is in fact an under-sampling, which means that the images differ
by a blur.

In his milestone paper~\cite{lowe04:_distin}, Lowe has addressed this central problem and has proposed the so called
scale-invariant feature transform (SIFT) descriptor,
that is claimed to be invariant to image translations and rotations,  to scale
changes (blur), and robust to illumination changes. It is also
surprisingly robust to large enough orientation changes of the
viewpoint (up to 60 degrees when one of the images is in frontal view). Based on the scale-space
theory~\cite{lindeberg94scalespace}, the SIFT procedure simulates
all Gaussian blurs and normalizes local patches around  scale
covariant image key points that are Laplacian extrema. A number of
SIFT variants and extensions, including
PCA-SIFT~\cite{DBLP:conf/cvpr/KeS04} and gradient
location-orientation histogram
(GLOH)~\cite{mikolajczyk05:_performance}, that claim to have better
robustness and distinctiveness with scaled-down complexity have been
developed ever since~\cite{1273738, 1180760}. Demonstrated to be
superior to other descriptors~\cite{MikolajczykS03,
mikolajczyk05:_performance}, SIFT has been popularly applied for
scene recognition~\cite{1178710, DBLP:conf/nips/MoreelsP04,
DBLP:conf/psivt/Ruiz-del-SolarLD07, veloso-learning,
DBLP:conf/clor/GordonL06, 1291311} and
detection~\cite{DBLP:conf/scia/FritzSKP05,
DBLP:conf/sspr/NegreTGHLC06}, robot
localization~\cite{bennewitz-metric, and-scalable, joseph-building}, image
retrieval~\cite{DBLP:conf/civr/HareL04}, motion
tracking~\cite{1018380, 1029650}, 3D modeling and
reconstruction~\cite{riggi06, 1231035}, building
panoramas~\cite{1141966, brown03recognising}, or photo
management~\cite{1242816, 1180737, 1160945}. Recently, a variant of SIFT named
affine-SIFT (ASIFT) has been mathematically proved to be fully affine invariant and has
shown to give excellent performance even
under very large viewpoint changes~\cite{morel2009asift}.

The initial goal of the SIFT method is to compare two images (or two
image parts)  that can be deduced from each other (or from a common
one) by a rotation, a translation,  and a zoom. In this method, following a classical
paradigm, stable points of interest are supposed to lie at extrema
of the Laplacian of the image in the image scale-space
representation.  The scale-space representation introduces a
smoothing parameter $\sigma$.  Images $u_0$  are smoothed at
several scales to obtain $w(\sigma, x,y) := ( G_\sigma \ast
u_0)(x,y)$, where
$$G_\sigma(x,y)=\frac{1}{2\pi\sigma^2} e^{-(x^2+y^2)/2 \sigma^2}$$
is the 2D-Gaussian function with integral 1 and standard deviation $\sigma$.  The  notation~$\ast$ stands for the space
2-D convolution in $(x,y)$.  The description of the SIFT method involves sampling issues, which we shall discuss later.

Taking apart all sampling issues and several thresholds whose aim it is to eliminate unreliable features, the whole SIFT
method can be summarized in one single sentence:

\medskip \noindent{\it {\bf One sentence description} The SIFT method computes scale-space extrema $(\sigma_i, x_i,y_i)$
of the  Laplacian in $(x,y)$ of $w(\sigma, x,y)$, and then
samples for each
one of these extrema a square image  patch whose origin is $(x_i,y_i)$, whose $x$-direction is one of the dominant
gradients around $(x_i,y_i)$, and whose sampling rate is proportional to (and usually smaller than) $\sqrt{\sigma_i^2+\bc^2}$. }

The constant $\bc\approx 0.5$ is the
tentative standard deviation of
the initial image blur. The resulting samples of the digital patch at scale
$\sigma_i$ are encoded by their gradient direction, which is
invariant under nondecreasing contrast changes. This accounts for
the robustness of the method to illumination changes. In addition,
only local histograms of the direction of the gradient are kept,
which accounts for the robustness of the final descriptor to changes
of view angle (see Fig. \ref{fig7lowe}).

\begin{figure}
\centering
%\includegraphics[width=3cm]{figures/pise-lowe7.jpg}
\includegraphics[width=3cm]{figures/pise-nfa001.jpg}
\caption{A result of the SIFT method, using an outliers elimination method  \cite{ rabin07acontrario}.
Pairs of matching points are connected by segments.}
\label{SIFTexample}
\end{figure}

The  goal of this  paper is to give the mathematical formalism
to examine whether the method indeed is scale invariant, and to discuss whether its main
assumption, that images are  well-sampled under Gaussian blur, does not entail significant errors. We shall not propose a new variant or an
extension of the SIFT method; on the contrary we intend to
demonstrate
that no other method will ever improve more than marginally the SIFT scale invariance (see Figs.~\ref{SIFTexample} and
\ref{fig:sift:scale:inv} for striking examples). To the best of our knowledge,
and in spite of the more
than ten thousand papers quoting and using SIFT, the analysis presented here does not seem to have been done previously.

The paper is organized as follows. In Section \ref{operators}, a
simple formalism is introduced to obtain a condensed
description of the SIFT shape encoding
method. Using this formalism Section~\ref{scaleandSIFT} proves
mathematically that the SIFT method indeed computes translation,
rotation and scale invariants. This proof is correct under the main
assumption that the image initial blur is Gaussian, and that
images with a Gaussian blur larger than $0.8$ can be accurately retrieved by interpolation from their samples.
Section~\ref{blurassumption} checks the validity of this crucial
well sampling assumption through an experimental procedure and
mathematical proofs.


\section{Image operators formalizing SIFT\label{operators}}

The analysis of the scale invariance is much easier on the continuous images
whose samples form the digital image. The Shannon-Whittaker interpolation permits a perfect reconstruction of a continuous image from a discrete one, when the continuous image has been
well-sampled~\cite{shannon48mathematical}. Under the assumption that
a Gaussian filtering can deliver well-sampled images up to a negligible error, (the validity  of this assumption  will be
confirmed in Section~\ref{blurassumption}), this section
gives a formalized description of the SIFT procedure.


We denote by  $u(\bx)$ a continuous  image defined for every $\bx=(x,y)\in \RR^2$.
All  \textit{continuous image} operators, including the
sampling operator itself, will be written in  capital letters $A$, $B$. Their composition will be, for the sake of simplicity, written as a mere juxtaposition $AB$.
For any similarity transform $A$ its application on
$u$ is defined as $A u(\bx):= u(A\bx)$. For instance
$H_\lambda u(\bx):=u(\lambda\bx)$ denotes an expansion of
$u$ by a factor $\lambda^{-1}$. In the same way if $R$ is a
rotation, $R u := u(R\bx)$ is the image rotation by $R^{-1}$.


\subsection{Sampling and interpolation\label{samplingandinterpolation}}

Digital (discrete) images are only defined for  $\bk=(k,l)\in\relative^2$ and are denoted in bold by $\bu(\bk)$.
For example the  $\delta$-sampled image
$\bu=\SA_\delta u$ is  defined on $\relative^2$ by \begin{equation} (\SA_\delta u)(k,l):=u(k\delta,
l\delta); \label{defsampling}\end{equation} Conversely,  the
Shannon interpolate of a digital image is defined as follows
\cite{gasquet1999faa}. Let $\bu$ be a digital image, defined on
$\relative^2$ and such that $\sum_{\bk\in \relative^2}|\bu(\bk)|^2<\infty$
and $\sum_{\bk\in \relative^2}|u(\bk)|<\infty.$ These
conditions are for example satisfied if the digital image has a
finite number of non-zero samples. We call Shannon
interpolate $\bI\bu$ of $\bu$ the only $L^2(\RR^2)$ function $u$ which coincides with $\bu$ on the samples $\bk$ and with spectrum  support contained in
$(-\pi,\pi)^2$.   $\bI\bu$ is defined by the Shannon-Whittaker formula
\begin{equation}
\label{eqn:shannon:interpolation}
\bI\bu(x,y):=\sum_{(k,l)\in \relative^2} \bu(k,l)\hbox{sinc}
(x-k)\hbox{sinc} (y-l),
\end{equation}
where $\hbox{sinc }x:=\frac{\sin \pi
x}{\pi  x}.$
The Shannon interpolation  has the fundamental property
$\SA_1\bI\bu=\bu$. Conversely, if $u$ is $L^2$ and  band-limited in
$(-\pi,\pi)^2$, then
\begin{equation}\label{interpolationrecovery}\bI\SA_1 u = u.\end{equation}
In that ideal situation we  say that $u$ is {\it band-limited}. We shall also say that the digital image $\bu=\bS_1u$ is
{\it well-sampled} if it was obtained from a band-limited image $u$, and therefore permits to go back to $u$.

While the Shannon interpolation~\eqref{eqn:shannon:interpolation}  allows a perfect
reconstruction of a continuous image from a discrete one when the image is
well-sampled, this interpolation $\bI$ is unfortunately impractical since it assumes an infinite number
of available samples. Instead, the DFT interpolation $\bI_d$ that relies only on a limited number of
samples in the discrete image rectangle is predominant in image processing.  The DFT interpolation $\bI_d$
will be defined in Section~\ref{subsec:shannon:vs:DFT} where its approximation to
the Shannon interpolation $\bI$ will be checked.
In the following we will thus proceed in this dual framework. For invariance proofs involving
translations, rotations and  zooms, the Shannon interpolation is necessary  because we need to reason in a space
of functions on $\real^2$ which is invariant by these geometric transforms. When it comes to
 image interpolation operations on the limited image domain, $\bI_d$ will be invoked.  It is well documented that the difference between these two interpolations causes negligible errors, but we shall anyway confirm this fact here.

\subsection{The Gaussian filtering}
The Gaussian convolution that implements the scale-space plays a key role in the SIFT procedure.
$G_\sigma$ will denote  the  convolution operator on $\RR^2$  with  the
Gaussian kernel $G_\sigma(x_1,x_2)=\frac 1{2\pi\sigma^2}e^{-{(x_1^2+x_2^2)}/{2\sigma^2}}$, and the Gaussian kernel itself.
Thus we simply write $G_\sigma u(x,y):= (G_\sigma \ast u)(x,y).$
$G_\sigma$ satisfies the semigroup property
\begin{equation}
\label{continuoussemigroupproperty}G_{\sigma} G_{\beta}= G_{\sqrt{\sigma^2+ \beta^2}}.\end{equation} The proof of
the next formula is a mere change of variables in the integral defining the convolution.
\begin{equation}\label{commutationconvolutionzoomcontinus}G_\sigma  H_\gamma u=
H_\gamma G_{\sigma\gamma}  u.\end{equation}
The discrete Gaussian convolution applied to a digital image
is defined as a digital operator by
\begin{equation}\label{defGdelta} \bG_\delta \bu= : \bS_1 G_\delta \bI
\bu.\end{equation}
This  discrete convolution is nothing but the continuous Gaussian convolution applied to the underlying  continuous image.
This definition maintains the Gaussian semi-group property used repeatedly in SIFT,
\begin{equation}\label{discretesemigroupproperty}\bG_{\delta} \bG_{\beta}= \bG_{\sqrt{\delta^2+ \beta^2}}.\end{equation}
Indeed, using  twice \eqref{defGdelta} and once \eqref{continuoussemigroupproperty} and \eqref{interpolationrecovery},
$$\bG_{\delta} \bG_{\beta} \bu=\SA_1G_{\delta}   \bI\SA_1G_{\beta}  \bI\bu=\SA_1G_{\delta} G_{\beta}
\bI\bu=\SA_1G_{\sqrt{\delta^2+ \beta^2}}\bI\bu=\bG_{\sqrt{\delta^2+ \beta^2}}\bu.$$
(The same formulas are adapted without alteration when replacing $\bI$ by $\bI_d$.)

The SIFT method makes the following assumption, whose validity will be confirmed both experimentally
and mathematically in Section~\ref{blurassumption}.
\begin{assumption}
\label{proposition:gaussian:well:sampling}
For every $\sigma$  larger than 0.8 and every Shannon interpolated digital image $u_0$, the Gaussian blurred image $G_\sigma
u_0$ satisfies the Shannon inversion formula up to a negligible error, namely $I\bS_1G_\sigma u_0 \simeq G_\sigma u_0$.
\end{assumption}


\subsection{Formalized Scale Invariant Features Transform\label{condenseddescription}}
The Assumption~\ref{proposition:gaussian:well:sampling} that
the Gaussian pre-filtering leads to a nearly aliasing-free subsampling
allows a perfect reconstruction of continuous images from discrete ones with
Shannon-Whittaker interpolation. The main steps of the SIFT method can therefore
be formalized in a continuous setting as follows.

\begin{enumerate}
\item {\bf Geometry:} there is an underlying infinite resolution  planar image $u_0(\bx)$ that has undergone a
similarity $A u_0$
(modeling the composition of a rotation, a translation, and a homothety) before sampling.
\item {\bf Sampling and blur:} the camera blur is assumed to be a Gaussian with
standard deviation $\bc$.
%The typical value of $\bc$ will be fixed thereafter. In Lowe's paper, $\bc$ belongs to $[0.5, 0.8]$.
The initial digital image is therefore $\bu=\bS_1G_{\bc}Au_0$;
\item  {\bf Sampled scale space:} the SIFT method computes enough samples of the scale space function
$u(\sigma,\cdot)=G_\sigma G_{\bc}Au_0$ to detect accurately ``key points''
$(\sigma, \bx)$, defined as scale and space local extrema of
$\Delta u(\sigma,\cdot)$.
\item  {\bf Covariant resampling:} a $32\times 32$ grid  centered at the key point is used to sample $u(\sigma, \cdot)$  around each key point $(\sigma, \bx)$. The  grid mesh is
proportional to
$\sqrt{\sigma^2 + \bc^2}$. The directions of the sampling grid axes are fixed by a dominant direction of
$\nabla u(\sigma,\cdot)$ in a
neighborhood of the key point, whose size is also proportional to the key point scale $\sigma$. This yields for each key point a rotation, translation and scale invariant square sampled subimage samples in which the
 four parameters of $A$ have been eliminated (see Fig. \ref{fig.SIFTkeypoints});
\item {\bf Illumination invariance:} the final SIFT descriptors keep mainly the orientation of the samples gradient, to gain
invariance with respect to light conditions.
\end{enumerate}


\begin{figure}
\centering
\includegraphics[width=9cm, height=4cm]{fig7Lowe04.jpg}
\caption{\small Each key-point is associated {\it a square image patch whose size is proportional to the scale and
whose side direction is given by the assigned direction}.  Example of a $2 \times 2$ descriptor array  of orientation
histograms (right) computed from an $8 \times 8$ set of samples (left). The orientation histograms are quantized into 8
directions and the length of each arrow corresponds to the magnitude of the histogram entry.} \label{fig7lowe}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=7cm]{MuralKeypoints.jpg}
\caption{SIFT key points. The arrow starting point, length and the orientation signify respectively the key point
position, scale, and  dominant  orientation.  These features are claimed to be covariant to any image
similarity.\label{fig.SIFTkeypoints}}
\end{figure}
Steps 1 to 5 are the main steps of the method. We have omitted all details that are not relevant in the discussion to
follow. Let them be mentioned briefly. The Laplacian extrema are kept only if they are larger than a fixed threshold
that eliminates small features mainly due to noise. This threshold is not scale invariant. The ratio of the eigenvalues
of the Hessian of the Laplacian must be close enough to 1 to ensure a good key point localization. (Typically, straight
edge points have only one large Hessian eigenvalue, are poorly localized, and are therefore ruled out by this second
threshold, which is scale invariant.)

The SIFT method assumes that the initial image satisfies $\bc=0.5$
(meaning that it is the result of a convolution  with a Gaussian with standard deviation
$\bc$). This implies a slight under-sampling which is compensated by
a complementary Gaussian blur applied to the image. Following
Assumption~\ref{proposition:gaussian:well:sampling} it increases
the initial blur to 0.8. In accordance with this choice, a
$2 \times 2$ subsampling in the SIFT scale-space computations can be
made only when  a $2\times 0.8=1.6$ Gaussian blur has been reached.

Postponing the verification of the SIFT main Assumption~\ref{proposition:gaussian:well:sampling}
to Section~\ref{blurassumption}, the next Section proves that SIFT has an almost
perfect scale invariance, which is the main result of the present paper.

\section{Scale and SIFT: consistency of the method}
\label{scaleandSIFT}

Let $\mathbf{\mathcal T}$, $R$, $H$ and $G$ be respectively an arbitrary  image translation,
an arbitrary image rotation,  an
arbitrary image homothety, and an arbitrary Gaussian
convolution, all applied to continuous images. We say that there is
strong commutation if we can exchange the order of application of
two of these operators. We say that there is weak commutation
between two of these operators if we have (e.g.)
$R\mathbf{\mathcal T}=\mathbf{\mathcal T'}R$,
meaning that given $R$ and $\mathbf{\mathcal T}$ there is
$\mathbf{\mathcal T'}$ such that the former relation occurs. The
next lemma is straightforward.
\begin{lemma}\label{weakcommutation}
All of the aforementioned operators weakly commute. In addition,
$R$ and $G$ commute strongly.
\end{lemma}

In this Section, in conformity with the SIFT model of
Section~\ref{operators}, the digital image is a frontal view of an
infinite resolution ideal image $u_0$. In that case,
$ A= H\mathcal T R$ is  the composition of a rotation $R$, a
translation $\mathbf{\mathcal T}$ and a homothety $H$.
Thus the digital image is  $\bu=\bS_1 G_\delta H\mathcal TR u_0$,
for some $H$, $\mathcal T$,  $R$. Assuming that the image is not aliased boils down, by the experimental
results of Section~\ref{blurassumption}, to assuming $\delta\geq 0.8$.
The next Lemma shows that SIFT is rotation- and translation-invariant.

\begin{lemma} \label{invariancesiftrotationtranslation} For any rotation $R$
 and any translation $\mathcal T$, the SIFT descriptors of
$\bS_1G_\delta H\mathcal TRu_0$
are identical to those of $\bS_1G_\delta H u_0.$
\end{lemma}

\begin{proof}
Using the weak commutation of translations and rotations with all other operators (Lemma \ref{weakcommutation}), it is
easily checked that the SIFT method is rotation and translation invariant: The SIFT descriptors of a rotated or
translated image are identical to those of the original. Indeed, the set of scale space  Laplacian extrema is covariant
to translations and rotations. Then the normalization process for each SIFT descriptor situates the origin at each
extremum in turn, thus canceling the translation, and the local sampling grid defining  the SIFT patch has axes given
by peaks in its gradient direction histogram. Such peaks are translation invariant and rotation covariant. Thus, the
normalization of the direction also cancels the rotation.
\end{proof}


\begin{lemma}\label{approximatescaleinvarianceofSIFT}
Let $\bu$ and $\bv$ be two digital images that are frontal snapshots of
the same continuous flat image $u_0$, $\bu= \bS_1 G_\beta
H_\lambda u_0 \hbox{ and } \bv:=\bS_1 G_\delta H_\mu u_0,$
taken at different distances, with different Gaussian blurs and
possibly different sampling rates. Let $w(\sigma, \bx):=(G_\sigma
u_0)(\bx)$   denote the scale space of $u_0$. Then the scale spaces of
$\bu$ and $\bv$ are
$$u(\sigma, \bx)=w(\lambda\sqrt{\sigma^2+\beta^2}, \lambda\bx)\; \hbox{ and }\;
v(\sigma, \bx)=w(\mu\sqrt{\sigma^2+\delta^2}, \mu\bx).$$ If $(s_0,\bx_0)$ is a key point of $w$ satisfying
$s_0\geq \max (\lambda\beta, \mu\delta)$,  then it corresponds to a key point of $u$ at the scale $\sigma_1$ such
that $\lambda\sqrt{\sigma_1^2+\beta^2}=s_0$, whose SIFT descriptor is sampled with mesh $\sqrt{\sigma_1^2+\bc^2}$,
where $\bc$ is the tentative standard deviation of
the initial image blur as described in Section~\ref{condenseddescription}. In the
same way $(s_0,\bx_0)$ corresponds to a key point of $v$ at scale $\sigma_2$ such that
$s_0=\mu\sqrt{\sigma_2^2+\delta^2}$, whose SIFT descriptor is sampled with mesh $\sqrt{\sigma_2^2+\bc^2}$.
\end{lemma}

\begin{proof}
The interpolated initial images are by \eqref{interpolationrecovery}
$$u:= \bI\bS_1 G_\beta H_\lambda u_0 =  \bG_\beta H_\lambda u_0 \hbox{ and }
v:=\bI\bS_1 G_\delta H_\mu u_0= G_\delta H_\mu u_0.$$
Computing the scale-space of these images amounts to convolve them
for every $\sigma>0$ with $G_\sigma,$ which yields, using the Gaussian
semigroup property~\eqref{continuoussemigroupproperty} and
the commutation relation \eqref{commutationconvolutionzoomcontinus}:
$$u(\sigma, \cdot )= G_\sigma  G_\beta H_\lambda u_0=
G_{\sqrt{\sigma^2+\beta^2}} H_\lambda u_0= H_\lambda  G_{\lambda\sqrt{\sigma^2+\beta^2}} u_0.$$ By
the same calculation, this function is compared by SIFT with
$$v(\sigma,\cdot)=
H_\mu G_{\mu\sqrt{\sigma^2+\delta^2}} u_0$$.
Let us set $w(s,\bx):= (G_s u_0)(\bx)$. Then the scale spaces compared
by SIFT are
$$u(\sigma, \bx)=w(\lambda\sqrt{\sigma^2+\beta^2}, \lambda\bx)\; \hbox{ and }\;
v(\sigma, \bx)=w(\mu\sqrt{\sigma^2+\delta^2}, \mu\bx).$$ Let us consider an extremal point $(s_0,\bx_0)$ of the
Laplacian of the scale space function $w$. If $s_0\ \geq \max(\lambda\beta, \mu\delta)$, an extremal point occurs at
scales $\sigma_1$ for (the Laplacian of) $u(\sigma,\bx)$ and $\sigma_2$ for (the Laplacian of) $v(\sigma, \bx)$
satisfying
\begin{equation}\label{scaletransform}
s_0=\lambda\sqrt{\sigma_1^2+\beta^2}=\mu\sqrt{\sigma_2^2+\delta^2}.\end{equation} We recall that each SIFT descriptor
at a key point $(\sigma_1,\bx_1)$ is computed from space samples of $\bx\to u(\sigma, \bx)$. The origin of the local
grid is $\bx_1$, the intrinsic axes are fixed by one of the dominant directions of the gradient of
$u(\sigma_1,\cdot)$ around $\bx_1$, in a circular neighborhood whose size is proportional to $\sigma_1$.  The SIFT
descriptor sampling rate around the key point is proportional to $\sqrt{\sigma_1^2+\bc^2}$ for $u(\sigma_1,\bx)$,
and to $\sqrt{\sigma_2^2+\bc^2}$ for $u(\sigma_2, \bx)$, as described in Section~\ref{condenseddescription}.
\end{proof}


\begin{theorem} \label{theo:scale} Let $\bu$ and $\bv$ be two digital images that are frontal
snapshots of the same continuous flat image $u_0$, $\bu= \bS_1
G_\beta H_\lambda\mathcal T R u_0 \hbox{ and } \bv:=\bS_1
G_\delta H_\mu u_0,$ taken at different distances, with
different Gaussian blurs and possibly different sampling rates, and
up to a camera translation and rotation around its optical axis.
Without loss of generality, assume $\lambda\leq \mu$. Then if  the
initial blurs are identical for both images (if $\beta=\delta=\bc$), then each SIFT descriptor of
$\bu$ is identical to a SIFT descriptor of $\bv$. If $\beta\neq \delta$
(or $\beta=\delta \neq \bc$), the SIFT descriptors of $\bu$ and $\bv$
become (quickly) similar when their scales grow, namely as soon as
$\frac{\sigma_1}{\max(\bc, \beta)}\gg 1$ and $\frac{\sigma_2}{\max
(\bc, \delta)} \gg 1$,  where $\sigma_1$ and
$\sigma_2$ are respectively the scales of the key points in the two
images.
\end{theorem}

\begin{proof} By the result of Lemma \ref{invariancesiftrotationtranslation}, we can neglect the  effect of translations and
rotations. Therefore assume without loss of generality that the images under comparison are as in Lemma
\ref{approximatescaleinvarianceofSIFT}. Consider a key point $(s_0, \bx_0)$ of $w$ with scale $s_0\geq
\max(\lambda\beta, \mu\delta).$   Following Lemma~\ref{approximatescaleinvarianceofSIFT}, there is a corresponding key
point $(\sigma_1, \frac{\bx_0}\lambda)$ for $\bu$ whose sampling rate is fixed by the method to $\sqrt{\sigma_1^2+\bc^2}$ and a corresponding key
point $(\sigma_2, \frac{\bx_0}\mu)$ whose sampling rate is fixed by the method to $\sqrt{\sigma_2^2+\bc^2}$ for $\bv$. To have a common reference
for these sampling rates, it is convenient to refer to the corresponding sampling rates for $w(s_0,\bx)$, which are
$\lambda\sqrt{\sigma_1^2+\bc^2}$ for the SIFT descriptors of
$\bu$ at scale $\sigma_1$, and
$\mu\sqrt{\sigma_2^2+\bc^2}$ for the descriptors of $\bv$ at scale $\sigma_2$.
Thus the SIFT descriptors of $\bu$ and $\bv$ for $\bx_0$ will be identical  if and only if
$\lambda\sqrt{\sigma_1^2+\bc^2}=\mu\sqrt{\sigma_2^2+\bc^2}$.  Since we  have
$\lambda\sqrt{\sigma_1^2+\beta^2}=\mu\sqrt{\sigma_2^2+\delta^2}$, the SIFT descriptors of $\bu$ and $\bv$ are identical if and only if
\begin{equation}\label{implicationapproximative}
\lambda\sqrt{\sigma_1^2+\beta^2}=\mu\sqrt{\sigma_2^2+\delta^2}\Rightarrow
\lambda\sqrt{\sigma_1^2+\bc^2}=\mu\sqrt{\sigma_2^2+\bc^2}.
\end{equation}
In other terms
$\lambda\sqrt{\sigma_1^2+\bc^2}=\mu\sqrt{\sigma_2^2+\bc^2}$ if and only if
\begin{equation}\label{conditiondescriptorsequal}\lambda^2\beta^2-\mu^2\delta^2=(\lambda^2-\mu^2)\bc^2.\end{equation}
Since $\lambda$ and $\mu$ correspond to camera distances to the observed object $u_0$, their values are arbitrary.
Thus in general the only way to get \eqref{conditiondescriptorsequal} is to have $\beta=\delta=\bc$,  which  means that
the blurs of both images have been guessed correctly.

The second statement is straightforward: if $\sigma_1$ and $\sigma_2$
are large enough with respect to $\beta$, $\delta$ and $\bc$, the
relation
$\lambda\sqrt{\sigma_1^2+\beta^2}=\mu\sqrt{\sigma_2^2+\delta^2}$,
implies $\lambda\sqrt{\sigma_1^2+\bc^2}\approx
\mu\sqrt{\sigma_2^2+\bc^2}$.
\end{proof}


The almost perfect scale invariance of SIFT stated in Theorem \ref{theo:scale} is illustrated by the striking example
of Fig.~\ref{fig:sift:scale:inv}. The 25 SIFT key points of a very small image $\bu$ are compared to the 60 key points
obtained by zooming in $\bu$ by a 32 factor:  The resulting digital image is $\bv:=\bS_{\frac{1}{32}} \bI_d \bu$, again obtained by
zero-padding. For better observability, both images  are displayed with the same size by enlarging the pixels of $\bu$.
Almost each
key point (18 out of 25) of $\bu$ finds its counterpart in $\bv$. 18 matches are detected between the descriptors as shown
on the right.
Let us check how this extreme example is covered by Theorem
\ref{theo:scale}. We compare an initial image $\bu=\bS_1 G_\delta \bI_d\bu_0$ (with $\delta=\bc$) with its zoomed in version $\bv=\bS_{\frac 1{32}} G_\delta  \bI_d \bu_0$. But  we have by~\eqref{commutationconvolutionzoomcontinus} $$\bv = \bS_{\frac 1{32}} G_\delta\bI_d \bu_0=\bS_1 H_{\frac 1{32}}G_\delta\bI_d\bu_0=
\bS_1 G_{32\delta} H_{\frac 1{32}}\bI_d\bu_0.$$
Here the numerical application of the relations in the above proof give:
We want (\ref{implicationapproximative}) to hold approximately, where $\mu=1$, $\lambda = \frac{1}{32}$, $\beta = 32\delta$.
Thus we want
$\frac 1{32}\sqrt{\sigma_1^2+(32\delta)^2}= \sqrt{\sigma_2^2+\delta^2}$
to imply $\frac 1{32}\sqrt{\sigma_1^2+\bc^2}\approx
\sqrt{\sigma_2^2+\bc^2}$ which means
$\sqrt{(\frac{\sigma_1}{32})^2+\bc^2}= \sqrt{\sigma_2^2+\bc^2}$
to imply $\sqrt{(\frac {\sigma_1}{32})^2+(\frac{\bc}{32})^2}\approx
\sqrt{\sigma_2^2+\bc^2}$. This is true only if $\sigma_1$ is significantly larger than $32$, which is true, since $\sigma_1$ is the scale of the SIFT descriptors in the image $\bv$, which has been zoomed in by a 32 factor.

\begin{figure}[htbp]
\begin{center}
\begin{tabular}{ccc}
\includegraphics[width=4.17cm]{./figures/sift_keys1.jpg}
&
\includegraphics[width=4cm]{./figures/sift_keys2.png}
&
\includegraphics[width=6cm]{./figures/sift_zoomx32_18matches.jpg}
\end{tabular}
\end{center}
\caption{Scale invariance of SIFT, an illustration of Theorem \ref{theo:scale}. Left:  a very small digital image $\bu$
with its 25 key points. For the conventions to represent key points and matches, see the comments in Fig.
\ref{fig.SIFTkeypoints}. Middle: this image is over sampled by a 32 factor to $\bS_{\frac{1}{32}} \bI_d \bu$. It has 60 key
points. Right: 18 matches found between $\bu$ and $\bS_{\frac{1}{32}} \bI_d \bu$. A zoom of the small image $\bu$ on the up-left corner
is shown in the bottom left. It can be observed that all the matches are correct.} \label{fig:sift:scale:inv}
\end{figure}




By the second part of Theorem~\ref{theo:scale}  the reliability of
the SIFT matching increases with scale. This fact is illustrated in
Fig.~\ref{fig:sift:scale:approx}. Starting from a high resolution
image $\bu_0$, two images $\bu$ and $\bv$ are obtained by simulated zoom
out, $\bu=\bS_1 G_\beta H_\lambda\bI_d\bu_0=\bS_\lambda G_{\lambda
\beta} \bI_d \bu_0$ and $\bv=\bS_\mu G_{\mu \delta} \bI_d \bu_0,$ with
$\lambda=2$, $\mu=4$, $\beta=\delta=0.8$. Pairs of SIFT descriptors
of $\bu$ and $\bv$ in correspondence, established by a SIFT matching,
are compared using an Euclidean distance $\mathbf{d}$.   The scale
rate $\sigma_1 / \sigma_2$ as well as the distance $d$ between the
matched key points are plotted against $\sigma_2$ in
Fig.~\ref{fig:sift:scale:approx}. That $\sigma_1 / \sigma_2 \approx
2$ for all key points confirms that the SIFT matching process is
reliable. As stated by the theorem, the rate $\sigma_1 / \sigma_2$
goes to $\mu / \lambda=2$ when $\sigma_2$ increases, and the
distance $\mathbf{d}$ goes down. However, as also apparent in the numerical result, when the scale is small
($\sigma_2 < 1$), $\sigma_1 / \sigma_2$ is very different from 2 and
$\mathbf{d}$ is large.

\begin{figure}[htbp]
\begin{center}
\begin{tabular}{ccc}
\includegraphics[width=8cm]{./figures/graf1.jpg} &
\includegraphics[width=4cm]{./figures/sift_scale_approx_dist_im1.jpg} &
\includegraphics[width=2cm]{./figures/sift_scale_approx_dist_im2.jpg} \\
\end{tabular}
\begin{tabular}{c}
\includegraphics[width=10cm]{./figures/sift_scale_approx_dist2.jpg} \\
\end{tabular}
\end{center}
\caption{Top (from left to right): $\bu_0$, $\bu$, $\bv$. Middle: Rate
of scales $\sigma_1 / \sigma_2$ of matched key points in $\bu$ and $\bv$
against $\sigma_2$. Bottom: Distance between matched descriptors of
$\bu$ and $\bv$ against $\sigma_2$.} \label{fig:sift:scale:approx}
\end{figure}

\section{The right Gaussian blur to achieve well-sampling}
\label{blurassumption}

This section first checks the close relation between the impractical Shannon
interpolation, with which we had proved the scale invariance of SIFT, and the
DFT interpolation that is always applied in practice for image resampling.  The rest of the section
checks the SIFT main Assumption~\ref{proposition:gaussian:well:sampling}.
Assumption~\ref{proposition:gaussian:well:sampling} is  that
a sufficiently broad Gaussian blur may lead to well-sampled images.
In principle, a Gaussian blur cannot lead to well-sampled images, because it is not {\it
stricto sensu} band limited. Thus, after Gaussian blur, the image is no more well sampled in the sense
defined above. However, it will shown that Assumption~\ref{proposition:gaussian:well:sampling}
is approximatively true. Section~\ref{subsec:gaussian:experiment}
defines an experimental procedure which checks that a Gaussian blur works.
This procedure will fix the minimum standard deviation of the Gaussian blur ensuring well-sampling,
up to a minor error. This error will be formally computed and shown to be very small in
Section~\ref{subsec:gaussian:math}.


\subsection{Shannon interpolation vs DFT interpolation}
\label{subsec:shannon:vs:DFT}

The Shannon interpolation, assuming infinitely many samples, is impractical. The image processing interpolation is slightly different. A digital image is usually given by its samples on a rectangle (for  simplicity we take it to be a square $\{0,\dots N-1\}^2$). These samples are (implicitly) extended by $N$-periodicity into a periodic image on $\real^2$, and it is  assumed that the resulting underlying image is band-limited with spectrum in $[-\pi, \pi]^2$. Needless to say, the periodicity assumption is enforced and contradicts in some extent the band-limited assumption. Even if the samples came from a really band-limited image, the loss of the samples outside the sampling square is irreversible. Thus, the band limited image and periodic image which is finally interpolated from the samples is different from the Shannon interpolate, but easier to compute. Indeed, given
a digital image $\bu(\bk)$ with $\bk=(k,l)\in \{0,\dots, N-1\}^2$  its  band-limited $N$-periodic interpolate is nothing but the trigonometric polynomial
\begin{equation}\label{dftinterpolated}
\bI_d u(\bx) =\!\!\! \sum_{\bm=(m,n)\in[-N/2,N/2-1]^2}  \tilde
u_\bm e^{\frac{2i\pi\bm.\bx}N},
\end{equation}
%\textbf{(Guoshen: I think it should be
%\begin{equation}\label{dftinterpolated}
%\bI_d \bu(\bx) =\!\!\! \sum_{\bm=(m,n)\in[-N/2,N/2-1]^2}  \tilde
%u_\bm e^{\frac{2i\pi\bm.\bx}N},
%\end{equation}
%)}
where $\tilde u_\bm$ with $\bm = (m,n)$ are the Discrete Fourier
Transform (DFT) coefficients of the $N^2$ samples $u(\bk)$.

%This is a somewhat confusing situation in image processing: we have the Shannon impractical interpolation $\bI$ and a practical discrete one,  $\bI_d$. In both cases the underlying image is supposed to be band-limited. Nevertheless, since only the image samples on a square are known, applying $\bI$ forces one to fix values for the unknown samples outside the square. Thus, $\bI$ is impractical but, it remains the correct tool to make formal proofs on images. On the other hand, $\bI_d$ is no less correct than $\bI$ and permits fast computations. It is therefore predominant in image processing. We shall therefore make our proofs in this dual framework. For sub-sampling issues, the images will be interpolated by $\bI_d$. For invariance proofs involving translations, rotations and zooms, the Shannon interpolation is necessary  because we need to reason in a space of functions on $\real^2$ which is invariant by these geometric transforms.

It is not the object of this paper to investigate in depth the interpolation error caused by ignoring the samples outside the image domain. But, still, a fast numerical experiment was made to get some hints about this {\it DFT interpolation error}. The classic digital test image Lena was inserted in a larger digital image $\bu_z$ by surrounding its samples with a frame of black  (zero) samples. The width of the frame evolved from $z=0$ to $z=50$.  For each image $\bu_z$, a $2 \times 2$ zoom-in by zero-padding~\footnote{The zero-padding permits to compute the values on a finer grid of the digital image with~\eqref{dftinterpolated}. It proceeds as follows: a) apply the $N\times N$ DFT  to the samples $\bu(k,l)$, $(k,l)\in \{0,\dots, N-1\}^2$ to get the discrete Fourier coefficients, $\tilde u_{\bm}$, $\bm\in \{-N/2, N/2-1\}^2$. b) complete these DFT coefficients to get a matrix $\tilde u_\bm$, $\bm \in \{-N, N-1\}^2$, the new DFT coefficients being zeros; c) Apply the  inverse $
2N\times 2N$  DFT to get $4N^2$ samples $v_\bk$, $\bk\in \{0,\dots, 2N\}^2$. Then $v_{2\bk}=u_\bk$ for $\bk \in \{0,\dots, N-1\}^2$.} was performed.  Let us call $\bv_z$ the digital image which is the restriction of the zoomed in image to the new samples (excluding the known samples, and also all samples in the added frame).

To check the impact of the exterior samples on the interpolation inside, the root mean square error
$\textrm{RMS}(\bv_z, \bv_{z+1})$  was computed. As plotted in Fig.~\ref{bandwidth_vs_error}-(a), $\textrm{RMS}(\bv_z, \bv_{z+1})$ converges to zero when the width of the black frame added outside the image increases. Assuming therefore $\bv_{50}$ to be very close to the final Shannon interpolate (which fixes all exterior samples to zero), Fig.~\ref{bandwidth_vs_error}-(b) shows  the root mean square error $\textrm{RMS}(\bv_z, \bv_{50})$ as a function of $z$. The initial error (about 2 for an image ranging from 0 to 255) is not negligible at all.
However, for $z=20$ this error becomes negligible. This simple experiment (which should be completed by a theoretical study in the spirit of Section~\ref{blurassumption}) indicates  that  {\it the perturbations of image interpolation due to the ignorance of exterior samples only affect significantly the image samples near the boundary. The samples at distance larger than 20 to the image boundary can be indifferently computed by DFT or Shannon interpolation. } This justifies our dual use of the  interpolation operators $\bI$ and $\bI_d$: They accurately coincide away from the image boundary.

\begin{figure}
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=7cm]{figures/bandwidth_vs_error.jpg} &
 \includegraphics[width=7cm]{figures/borderwidth_vs_error.jpg} \\
\textbf{(a)} & \textbf{(b)}\\
\end{tabular}
\end{center}
\caption{\textbf{(a).} The root mean square error  $\textrm{RMS}(\bv_z, \bv_{z+1})$ between successive interpolation. \textbf{(b).}
The root mean square error  $\textrm{RMS}(\bv_z, \bv_{50})$ between the current and the last interpolations.}
\label{bandwidth_vs_error}
\end{figure}

\subsection{Experimental Verification of Assumption~\ref{proposition:gaussian:well:sampling}}
\label{subsec:gaussian:experiment}

\begin{figure}[]
\begin{center}
\begin{tabular}{cc}
\hspace{-3ex} \includegraphics[width=6cm]{./figures/Lena1.jpg} &
\includegraphics[width=6cm]{./figures/GaussSize_Lena_errorl2_relative.jpg}  \vspace{-1ex}
\\
\hspace{-3ex}\includegraphics[width=6cm]{./figures/GaussSize_Lena_u1_c0p4.jpg}
&
\includegraphics[width=6cm]{./figures/GaussSize_Lena_u2_c0p4.jpg}  \vspace{-1ex}
\\
\hspace{-3ex}\includegraphics[width=6cm]{./figures/GaussSize_Lena_u1_c0p8.jpg}
&
\includegraphics[width=6cm]{./figures/GaussSize_Lena_u2_c0p8.jpg}  \vspace{-1ex}
\end{tabular}
\vspace{-4ex}
\end{center}
\caption{Top left: the classic test image $\bu$. Top right: relative error $\epsilon(\bu_1, \bu_2) = \|\bu_1-\bu_2\|_2 / \|\bu_1\|_2$ vs
$\ft2_r(4)$, the error introduced by a $2 \times 2$ sampling
after a Gaussian blur with $\sigma = \ft2_r(2)$.  Middle (from left to right): $\bu_1$ and $\bu_2$ (zoomed)
with $\ft2_r(2)=0.8$,  $\epsilon(\bu_1, \bu_2)=1.1e-02$. The aliasing error is large and conspicuous ringing artifacts appear
on the right image (refer to the electronic copy for better visibility). Bottom (from left to
right): $\bu_1$ and $\bu_2$ (zoomed) with $ \ft2_r(2)=1.6$.
$\epsilon(\bu_1, \bu_2)=9.2e-04$. These images are almost identical, which confirms that a Gaussian blur with standard deviation 1.6 is  sufficient before $2 \times 2$-sub-sampling.} \label{fig:guassian:MSELena}
\end{figure}

To understand the right Gaussian blur to achieve a well-sampled image, we shall
distinguish two types of  blur. The \textit{absolute} blur
with standard deviation $\ft2_a$ is the one that must be applied to
an ideal infinite resolution (blur free) image to create an
approximately band-limited image before 1-sampling.   The
\textit{relative} blur
$\sigma=\ft2_r(t)$ is the one that must be applied to a well-sampled image before a sub-sampling by a factor of $t$.
In the case of Gaussian blur, because of the semi-group formula
\eqref{continuoussemigroupproperty}, the relation between the
absolute and relative blur is
$$t^2 \ft2_a^2 = \ft2_r^2(t) + \ft2_a^2,$$ which yields
\begin{equation}
\label{eqn:cr:ca}
\ft2_r(t) = \ft2_a \sqrt{t^2-1}.
\end{equation}
In consequence, if $t \gg 1$, then $\ft2_r(t)
\approx \ft2_a t$.

Two experiments have been designed to calculate the anti-aliasing
absolute Gaussian blur $\ft2_a$ ensuring that an image is
approximately well-sampled. The first experiment finds the appropriate
relative blur $\bc_r(t)$ to achieve well-sampled images and then calculates
the absolute blur $\bc_a$ using~\eqref{eqn:cr:ca}.
$\bc_r(t)$. It compares for
several values of $\ft2_r(t)$ the digital images
$$\bu_1 :=\bG_{\ft2_r(t)} \bu=\bS_1 G_{\ft2_r(t)} \bI_d\bu\;\; \hbox{ and }\;\;
\bu_2 :=(\bS_{1/t}\bI_d) \bS_t
G_{\ft2_r(t)} \bI_d\bu,$$ where $\bu$ is an initial digital image that is
well-sampled, $\bS_{t}$ is a $t$-sub-sampling operator ($t$ is an integer), $\bS_{1/t}$ a $t$-over-sampling operator, and
$\bI_d$  the DFT  interpolation operator. The discrete convolution by a Gaussian is defined in \eqref{defGdelta}.
Since $t$ is an integer, the  $t$-sub-sampling is
trivial. The DFT over-sampling
$\bS_{1/t}\bI_d$ with an integer zoom factor $t$ is
computed by zero-padding as explained in footnote 1.
If the anti-aliasing filter size $\ft2_r(t)$ is too
small, $\bu_1$ and $\bu_2$ can be very different. The right value of  $\ft2_r(t)$ should be the
smallest value permitting that, systematically, $\bu_1 \approx \bu_2$.

Fig.~\ref{fig:guassian:MSELena} shows $\bu_1$ and $\bu_2$ with $t=2$ and
plots their relative error $\epsilon(\bu_1, \bu_2) = \|\bu_1-\bu_2\|_2 / \|\bu_1\|_2$. An anti-aliasing
filter with $\ft2_r(2)=0.8$ is clearly not broad enough:  $\bu_2$
presents strong ringing artifacts. The ringing artifact is instead
hardly noticeable with $\ft2_r(2)=1.6$. The value $\ft2_r(2)\approx
1.6$ is a good visual candidate, and this choice is confirmed by the
curve showing that $\epsilon(\bu_1, \bu_2)$ decays rapidly until $\ft2_r(2)$
gets close to $1.6$, and is stable and small thereafter. By
\eqref{eqn:cr:ca}, this value of $\ft2_r$ yields $\ft2_a \approx 0.8$.
This value has been confirmed by experiments on ten digital images.
A doubt can be cast on this experiment, however: Its result slightly
depends on the assumption that the initial blur on $\bu$ is equal to
$\ft2_a$.

In a second experiment, $\ft2_a$ has been evaluated directly by using a binary image $\bu_0$, thus containing a minimal
blur. As illustrated in Fig.~\ref{fig:guassian:MSELena:binary}, $\bu_0$ is obtained by binarizing the classic test image Lena
(Fig.~\ref{fig:guassian:MSELena}), the threshold being the image median value. Since $\bu_0$ is now close to be blur-free, we can compare
for several values of $\ft2_a$ and for $t=2$, which is large enough, the digital images
$$\bu_1 :=\bG_{t\ft2_a} \bu=\bS_1 G_{t\ft2_a} \bI_d\bu\;\; \hbox{ and }\;\;
\bu_2 := (\bS_{1/t}\bI_d) \bS_tG_{t\ft2_a} \bI_d\bu,$$ As shown in
Fig.~\ref{fig:guassian:MSELena:binary}, $\ft2_a=0.8$ is the smallest
value ensuring no visual ringing in $\bu_2$. Under this value, for
example for $\ft2_a=0.4$, clear ringing artifacts are present in
$\bu_2$. That $\ft2_a=0.8$ is the correct value is confirmed by the
$\epsilon(\bu_1, \bu_2)$ curve showing that the relative error decays
rapidly until $\ft2_a$ goes down to $0.8$, and is stable  and small
thereafter. The result, confirmed in ten experiments with different
initial images, is consistent with the value obtained in the first
experimental setting.

Fig.~\ref{fig:guassian:MSEnoise} illustrates the same experiment on
a Gaussian white noise image that has constant spectrum energy over
all frequencies. This is not at all an unrealistic image, since textured images can have parts with a flat spectrum  which is not so much different from white noise.   The critical value $\ft2_a=0.8$ is reconfirmed. Under this
value, for example for $\ft2_a=0.4$, $u_2$ is visually very different
from $u_1$, while the two very close (numerically and visually as well) with $\ft2_a=0.8$. The relative error for a white noise, about $4\%$ cannot be considered negligible. But it is a worst case, no digital image having such a flat spectrum. The  experimental value of this error will be retrieved by a formal calculation, and shown to be independent of the image size.


\begin{figure}[]
\begin{center}
\begin{tabular}{cc}
\hspace{-3ex}
\includegraphics[width=6cm]{./figures/LenaBinary1.jpg} &
\includegraphics[width=6cm]{./figures/GaussSize_LenaBinary_errorl2_relative.jpg} \vspace{-2ex}
\\
\hspace{-3ex}
\includegraphics[width=6cm]{./figures/GaussSize_LenaBinary_u1_c0p4.jpg}
&
\includegraphics[width=6cm]{./figures/GaussSize_LenaBinary_u2_c0p4.jpg} \vspace{-2ex}
\\
\hspace{-3ex}
\includegraphics[width=6cm]{./figures/GaussSize_LenaBinary_u1_c0p8.jpg}
&
\includegraphics[width=6cm]{./figures/GaussSize_LenaBinary_u2_c0p8.jpg}
\end{tabular}
\vspace{-4ex}
\end{center}
\caption{Top left: $\bu$ Binarized Lena (gray-levels 50 and 0). Top
right: $\epsilon(\bu_1, \bu_2) = \|\bu_1-\bu_2\|_2 / \|\bu_1\|_2$ vs $\ft2_a$, the error introduced by a $2 \times 2$ sampling
after a Gaussian blur with $\sigma = 2 \ft2_a$. Middle (from left to right):
$\bu_1$ and $\bu_2$ (zoomed) with $\ft2_a=0.4$. $\epsilon(\bu_1, \bu_2) = 4.7e-02$.
The aliasing error is large and conspicuous ringing artifacts appear
on the right image (refer to the electronic copy for better visibility).
Bottom (from left to right): $\bu_1$ and $\bu_2$ (zoomed) with $
\ft2_a=0.8$. $\epsilon(\bu_1, \bu_2) = 2.5e-03$. These last two images are almost identical, while strong aliasing artifacts are visible with $\ft2_a=0.4$.}
\label{fig:guassian:MSELena:binary}
\end{figure}


\begin{figure}[]
\begin{center}
\begin{tabular}{cc}
\hspace{-3ex}
\includegraphics[width=6cm]{./figures/WhiteNoise.jpg} &
\includegraphics[width=6cm]{./figures/GaussSize_WhiteNoise_errorl2_relative.jpg}  \vspace{-1ex}
\\
\hspace{-3ex}
\includegraphics[width=6cm]{./figures/GaussSize_WhiteNoise_u1_c0p4.jpg}
&
\includegraphics[width=6cm]{./figures/GaussSize_WhiteNoise_u2_c0p4.jpg} \vspace{-1ex}
\\
\hspace{-3ex}
\includegraphics[width=6cm]{./figures/GaussSize_WhiteNoise_u1_c0p8.jpg}
&
\includegraphics[width=6cm]{./figures/GaussSize_WhiteNoise_u2_c0p8.jpg} \vspace{-1ex}
\end{tabular}
\vspace{-4ex}
\end{center}
\caption{Top left: $\bu$ Gaussian white noise. Top
right: $\epsilon(\bu_1, \bu_2) = \|\bu_1-\bu_2\|_2 / \|\bu_1\|_2$ vs $\ft2_a$, the error introduced by a $2 \times 2$ sampling
after a Gaussian blur with $\sigma = 2 \ft2_a$. Middle (from left to right):
$\bu_1$ and $\bu_2$ (zoomed) with $\ft2_a=0.4$. $\epsilon(\bu_1, \bu_2) = 0.549$ and the images are also visually very dissimilar.
Bottom (from left to right): $\bu_1$ and $\bu_2$ (zoomed) with $
\ft2_a=0.8$. $\epsilon(\bu_1, \bu_2) = 0.0389$, and the images are visually identical.} \label{fig:guassian:MSEnoise}
\end{figure}

\subsection{Mathematical confirmation of Assumption~\ref{proposition:gaussian:well:sampling}}
\label{subsec:gaussian:math}

This section  confirms by exact computations the experimental estimates made on
simulated Gaussian white noise.  The spectrum of natural digital
images has a spectrum usually decaying much faster than a white noise. Thus,
the aliasing error can be safely estimated on white noise, which has
a flat spectrum. Sampling operations shown to be safe for white noise will therefore also be safe for
natural images. We will first calculate the aliasing error in the general case, and then
show that after an anti-aliasing Gaussian filtering with a relative blur $\sigma=2 \times 0.8=1.6$, the
aliasing error created by a $2 \times 2$ subsampling on a Gaussian white noise is negligible.

Consider a digital image $\bu(\bk)$ with $\bk=(k,l)\in \{0,\dots, N-1\}^2$ and its DFT interpolate
\begin{equation}
u(\bx)=(\bI_d\bu)(\bx) =\!\!\! \sum_{\bm=(m,n)\in[-N/2,N/2-1]^2}  \tilde
u_\bm e^{\frac{2i\pi\bm.\bx}N},\;\; \bx\in \real^2,
\end{equation}
where $\tilde u_\bm$ with $\bm = (m,n)$ are the Discrete Fourier
Transform (DFT) coefficients of the $N^2$ samples $u(\bk)$, and the Fourier Series coefficients of $u$ as well.   Recall the classic isometries between samples, DFT coefficients, and the $L^2$ norm of the image:
\begin{equation}\label{Fourierisometries}
||u||^2_{L^2([0,N]^2)}= ||\bu||^2_{l^2(\{0,\dots, N-1\}^2)}= N^2 ||\tilde u||^2_{l^2(\{-N/2,\dots, N/2-1\}^2)}.
\end{equation}
The scale space convolves
$u$ with a Gaussian $G_\sigma(\bx)=\frac
1{2\pi\sigma^2}e^{-\frac{\bx^2}{2\sigma^2}}.$ By computing its DFT, an
easy calculation shows that the result of the convolution of $u$
with $G_\sigma$ is
\begin{equation}
v(\bx):= (G_\sigma\ast u)(\bx)=\sum_{\bm\in[-\frac N2,\frac
N2-1]^2} \tilde {u}_{\bm} \hat
{G}_\sigma \left(\frac{2\bm\pi}N\right)e^{\frac{2i\pi\bm.\bx}N}
\end{equation}
with
$\hat G_\sigma(\xi)=e^{\frac{-\sigma^2\xi^2}2}$.
Indeed, $$(G_\sigma\ast e^{\frac{2i\pi\bm.\bx}N})(\bx)=\int_{\R^2}G_\sigma(\by)e^{\frac{2i\pi\bm.(\bx-\by)}N}d\by=
e^{\frac{2i\pi\bm.\bx}N}
\int_{\R^2}G_\sigma(\by)e^{-\frac{2i\pi\bm.\by}N}d\by= e^{\frac{2i\pi\bm.\bx}N}\hat{G}_\sigma\left(\frac{2\pi\bm}{N}\right).$$
The convolved digital image is
\begin{equation}
\label{eqn:DFT:v} v(\bk)=\sum_{\bm \in[-\frac N2,\frac N2-1]^2}
\tilde v_\bm e^{\frac{2i\pi\bm.\bk}N}.
\end{equation} Subsampling this image by a factor of $2 \times 2$ and
 interpolating the remaining
$N^2/4$ samples with  zero-padding interpolation yields a new trigonometric polynomial
\begin{equation}
\label{eqn:DFT:w} w(\bx)=\sum_{\bm\in [-N/4, N/4-1]^2} \tilde
w_\bm e^{\frac{2i\pi\bm.\bx}N}.
\end{equation}
Let us compute $\tilde w_\bm$.
\begin{lemma}
\label{sousech2-2D} Let $\bw=\bS_2\bv$ be the $\frac N2\times \frac N2$ image obtained by a $2 \times 2$ subsampling of a digital $N\times N$ image $\bv$.   Then $\tilde w _{m,n}$, the DFT of  $\bw$, satisfies for $m,n=-\frac N4,...,\frac N4-1$,
\begin{equation}\label{aliasing2Dparsousechantillonnage}
\tilde w_{m,n}=  \sum_{(\varepsilon_1,\varepsilon_2)\in\{0,1,-1\}}\tilde v_{m+\varepsilon_1\frac N2,
n+\varepsilon_2\frac N2}.
\end{equation}
\end{lemma}
Notice that $\tilde v_{m,n}$ is  supported  by $[-\frac N2, \frac N2-1]^2$ and that $\tilde w_{m,n}$ is supported by $[-\frac N4, \frac N4-1]^2$. For each $(m,n)$, only at most  four coefficients out of the nine $v_{m+\varepsilon_1\frac N2,
n+\varepsilon_2\frac N2} $ are non-zero, because if for example $m,n\leq 0$, $\tilde v_{m-\frac N2,n}$, $\tilde v_{m-\frac N2,n-\frac N2}$, $\tilde v_{m,n-\frac N2}$, $\tilde v_{m+\frac N2,n-\frac N2}$,  and $\tilde v_{m-\frac N2,n+\frac N2}$ are all null. To make apparent the four non-zero coefficients, we take the convention to $N\times N$ periodize  $(m,n)\to\tilde v_{m,n}$. With this convention, (\ref{aliasing2Dparsousechantillonnage}) rewrites
\begin{equation}\label{aliasing2Dparsousechantillonnagebis}
\forall (m,n)\in \left[-\frac N4, \frac N4-1\right]^2, \;\tilde w_{m,n}= \tilde v_{m,n}+ \tilde v_{m+ \frac N2,n+\frac N2}+ \tilde v_{m+ \frac N2,n}+ \tilde v_{m,n+\frac N2}.
\end{equation}
The digital image $\bw$ is twice smaller than $\bv$. To compare it with  $\bv$ we  must go back to the underlying trigonometric polynomials $v(\bx)$ and $w(\bx)$ such that  $\bw=\SA_2 w$ and
$\bv=\SA_1 v$, and we must compare their corresponding samples. Thus set $\bw^1=\SA_1 w$. The image $\bw^1$ is obtained by zero padding from $\tilde\bw$ by setting $\tilde w^1_{m,n}=0$ if $(m,n)\in \{-\frac N2,..., \frac N2-1\}^2\setminus \{-\frac N4, \dots,\frac N4-1\}^2$ and $\tilde w^1_{m,n}=\tilde w_{m,n}$ if $(m,n)\in \{-\frac N4,\dots, \frac N4-1\}^2$.
Then using~\eqref{eqn:DFT:v},~\eqref{eqn:DFT:w} and~\eqref{Fourierisometries}, we have
\begin{equation}
\label{eqn:error} ||v-w||_{L^2([0,N]^2)}^2 =||\bv-\bw_1||_{l^2(\{0,..., N-1\}^2)}^2=
N^2||\tilde v-\tilde w||_{l^2(\{-\frac N2,\dots,\frac N2-1\})}^2=e_h(v) + e_a(v).
\end{equation}
Using the above $N$-periodicity convention for the Fourier coefficients $\tilde v(m,n)$, we have proved
the following proposition.
%
\begin{proposition}
\label{prop:aliasing:error}
Let $u$ be a digital image and $v:=
G_\sigma \ast u$ be an image obtained after a Gaussian convolution with standard deviation
$\sigma$. Then the variance of the error incurred in doing a $2 \times 2$-sub-sampling followed by a zero-padding interpolation is
\begin{equation}
\label{eqn:error:downsample} ||v-w||_2^2 = e_h(v) + e_a(v),
\end{equation}
where
{\scriptsize
\begin{eqnarray}
\label{eqn:eh} e_h(v) & = &   N^2\left[\sum_{(m,n)\in[-N/2,
N/2-1]^2\setminus[-N/4, N/4-1]^2}|\tilde v(m,n)|^2 \right] \\
&=&   N^2\left[\sum_{(m,n)\in[-N/4, N/4-1]^2}\left(|\tilde
v({m,n+N/2})|^2 + |\tilde v({m+N/2,n})|^2 + |\tilde
v({m+N/2,n+N/2})|^2\right) \right]
\end{eqnarray}}
is the error due to the elimination of the high frequencies and
{\small
\begin{equation} e_a(v) =  {N^2}\left[\sum_{(m,n)\in[-N/4,
N/4-1]^2}|\tilde v({m,n+N/2}) + \tilde v({m+N/2,n}) + \tilde
v({m+N/2,n+N/2})|^2 \right]
\end{equation}
}
is the error introduced by the spectrum aliasing, with
\begin{equation}
\label{eqn:DFT:gauss} \tilde v(\bm) = {\tilde {u}({\bm})}\hat
G_\sigma\left(\frac{2\bm\pi}N\right).
\end{equation} \end{proposition}
%
Following Shannon~\cite{shannon48mathematical} we call white noise image with variance $\tau^2$ a continuous image
$$v(\bx)=\sum_{\bm\in [-N/2, N/2-1]^2} \tilde v_\bm e^{\frac{2i\bm\bx}N}$$
whose samples $v(k,l)$, $(k,l)\in \{0, \dots, N-1\}^2$ are i.i.d. $\mathcal N(0,\tau^2)$ Gaussian variables.
Using Proposition~\ref{prop:aliasing:error}, the next Corollary shows that with a Gaussian anti-aliasing filtering with
$\sigma = 2 \times 0.8 = 1.6$, the aliasing error created by a $2 \times 2$ subsampling on a Gaussian white noise image is very small.

\begin{corollary} \label{SNRBruit}Let $u$ be a Gaussian white noise image and
$v:=u \ast  G_\sigma$ be an image obtained after a Gaussian convolution with standard deviation
$\sigma$. Then the relative error incurred in doing a $2 \times 2$ subsampling followed by a zero-padding $2 \times 2$ interpolation is
\begin{equation}
\label{eqn:error:noise}
\epsilon(\sigma):= \frac{\sqrt{{\bf E}||v-w||^2}}{\sqrt{{\bf E}||v||^2}}  =
\sqrt{2 \frac{\sum_{\bm\in[-N/2, N/2-1]^2\setminus[-N/4, N/4-1]^2}
e^{-4\pi^2\sigma^2 \frac{|\bm|^2}{N^2}}}
{\sum_{\bm\in[-N/2, N/2-1]^2}
e^{-4\pi^2\sigma^2 \frac{|\bm|^2}{N^2}}}};
\end{equation}
\begin{equation}\label{approxerreur}
\epsilon(\sigma)\approx  \sqrt{  2\frac{\int_{[-\frac 12, \frac 12]^2\setminus[-\frac 14, \frac 14]^2} e^{-4\pi^2\sigma^2\bx^2}d\bx}{\int_{[-\frac 12, \frac 12]^2 } e^{-4\pi^2\sigma^2\bx^2}d\bx}}.
\end{equation}
In particular, $\epsilon(1.6)\approx 0.0389$.
\end{corollary}
The aim of \eqref{approxerreur} is to  point out that for white noise, the relative error incurred by a Gaussian sub-sampling is independent from the size $N$  of the image.

\begin{proof} The DFT coefficients $\tilde{u}(\bm)$ of a Gaussian
white noise image $u$ with standard deviation $\tau$ are independent Gaussian variables with standard
deviation $\frac{\tau}N$.   Using~\eqref{eqn:error:downsample} and~\eqref{eqn:DFT:gauss}, we therefore have

$\bE ||v-w||^2=\bE e_h(v) + \bE e_a(v)$ and
$$\bE e_a(v)=N^2\!\!\!\!\!\!\!\!\!\!\sum_{(m,n)\in \{-N/4, N/4-1\}^2}\bE \left|\tilde v(m,n+\frac N2)+\tilde v(m+\frac N2, n)+\tilde v(m+\frac N2, n+\frac N2)\right|^2$$
$$=N^2\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\sum_{(m,n)\in \{-N/4, N/4-1\}^2}\!\!\!\!\!\!\!\!\!\!\tau^2\left(\hat G_\sigma\left(\frac{2\pi}N(m, n\!+\!\frac N2)\right)^2\!+\!
\hat G_\sigma\left(\frac{2\pi}N(m\!+\!\frac N2, n)\right)^2\!+\!\hat G_\sigma\left(\frac{2\pi}N(m\!+\!\frac N2, n\!+\!\frac N2)\right)^2\right).$$
Using \eqref{eqn:eh} we deduce that
$$\bE ||v-w||^2=\bE e_h(v) + \bE e_a(v)=2\bE e_a(v)= 2 N^2 \tau^2 \sum_{\bm\in\{-N/2, \dots N/2-1\}^2\setminus\{-N/4,\dots, N/4-1\}^2} \hat G_\sigma\left(\frac{2\pi\bm}N\right)^2$$
$$= 2 N^2 \tau^2 \sum_{\bm\in\{-N/2, \dots N/2-1\}^2\setminus\{-N/4,\dots, N/4-1\}^2} e^{-4\pi^2\sigma^2\frac{|\bm|^2}{N^2}}$$
On the other hand,
$$\bE ||v||^2=\bE \sum_{(k,l)\in\{0,\dots, N-1\}^2}|v(k,l)|^2=N^2 \bE||\tilde v||^2_{l^2(\{-N/2,\dots,N/2-1\}^2)}=$$
$$=N^2\sum_{\bm\in \{-N/2,\dots,N/2-1\}^2}\bE|\tilde u_\bm|^2\hat G_\sigma(\frac{2\bm \pi}{N})^2=N^2\tau^2 \sum_{\bm\in\{-N/2,\dots,N/2-1\}^2} e^{-4\pi^2\sigma^2\frac{|\bm|^2}{N^2}}.$$
Thus
$$\frac{\sqrt{{\bf E}||v-w||^2}}{\sqrt{{\bf E}||v||^2}}  = \sqrt{2 \frac{\sum_{(m,n)\in\{-N/2, \dots N/2-1\}^2\setminus\{-N/4,\dots, N/4-1\}^2}
e^{-4\pi^2\sigma^2\frac{|\bm|^2}{N^2}}}{\sum_{(m,n)\in\{-N/2, \dots N/2-1\}^2}  e^{-4\pi^2\sigma^2\frac{|\bm|^2}{N^2}}}}$$
$$
\approx \sqrt{  2 \frac{\int_{[-\frac 12, \frac 12]^2\setminus[-\frac 14, \frac 14]^2} e^{-4\pi^2\sigma^2\bx^2}d\bx}{\int_{[-\frac 12, \frac 12]^2 } e^{-4\pi^2\sigma^2\bx^2}d\bx}}.$$
The last equivalent is computed by noticing that the discrete sums are Riemann sums for the integrals on the right.
For $\sigma=1.6$, this value  is approximately $0.0389$.
\end{proof}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The error $\epsilon(1.6)\approx 0.0389$ calculated from~\eqref{eqn:error:noise}
was numerically verified with the experiments on Gaussian white noise. It is actually not negligible,  but must be considered as a worst case which is never
approximated by real images.
Indeed, the power spectra of natural
images having a much faster decay than Gaussian white noise, the value
$\epsilon(1.6)$ is smaller on natural images: For Lena, Peppers,
Baboon and Barbara as shown in Figure~\ref{fig:standard:images},
$\epsilon(1.6)$ is respectively $9.2e-04$, $ 1.2e-03$, $1.9e-03$,
$2.0e-03$.  
\begin{figure}[htb]
\begin{center}
\begin{tabular}{cccc}
\includegraphics[width=3cm]{./figures/Lena.jpg} &
\includegraphics[width=3cm]{./figures/Peppers.jpg} &
\includegraphics[width=3cm]{./figures/Baboon.jpg} &
\includegraphics[width=3cm]{./figures/Barbara.jpg} \vspace{1ex}\\
\textbf{Lena} & \textbf{Peppers} & \textbf{Baboon} & \textbf{Barbara}
\end{tabular}
\end{center}
\caption{Standard test images. From left to right: Lena, Peppers, Baboon, and Barbara.} \label{fig:standard:images}
\end{figure}
\section{Conclusion} Our overall conclusion is that no substantial improvement of the SIFT method can be ever hoped,
as far as translation, rotation and scale invariance are concerned. As pointed out by several benchmarks, the
robustness and repeatability of the SIFT descriptors  outperforms other methods. However, such benchmarks mix three
very different criteria that, in our opinion, should have been discussed separately. The first one is the formal
invariance of each method when all thresholds have been eliminated. This formal invariance has been proved here for SIFT when the  initial blurs of both images are equal to a known value $\bc$, and it has been proved to be approximately true even with images having undergone very different blurs, like in the surprising experiment of fig. \ref{fig:sift:scale:approx}.
The second criterion is the practical validity of the sampling method used in SIFT, that has been again checked in the
present note. The last criterion is the clever fixing of several thresholds in the SIFT method ensuring robustness,
repeatability, and a low false alarm rate. This one has been extensively tested and confirmed in previous benchmark
papers (see also the recent and complete report \cite{Fevrier07}). We think, however, that the success of SIFT in
these benchmarks is primarily due to its full scale invariance.

\bibliographystyle{plain}
%\bibliography{these_biblio_final}
\bibliography{biblio_SIFT_consistency}

\end{document}
