\chapter{The SIFT Method}
\graphicspath{{figures/}} This chapter is devoted to Lowe's {\em Scale-Invariant Feature Transform} (SIFT
\cite{lowe04:_distin}), a very efficient image comparison method. 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. The method turned out to be also robust to large enough changes in view point angle, which explains its
success.  This method uses as fundamental tool the heat equation or, in other terms, the linear scale space. The heat
equation is used  to simulate all zooms out of both images that have to be compared. Indeed, these images may contain
similar objects taken at different distances. But at least two of the simulated zoomed in images should contain these
objects at the same apparent distance. This is the principal ingredient of the SIFT method, but  other invariance
requirements must be addressed as well.


   Sect. \ref{sec:guide} gives a detailed
  description of the SIFT shape encoding
method.   Sect. \ref{scaleandSIFT} proves mathematically that the SIFT method indeed computes translation, rotation and
scale invariants. This proof is correct under the main assumption that image blur can be assumed to be gaussian, and
that images with a gaussian blur larger than $0.6$ (SIFT takes $0.8$) are approximately (but  accurately) well-sampled
and can therefore be interpolated. Chapter. \ref{blurassumption}  checked the validity of this crucial gaussian blur
assumption.


\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 may then be established by matching the descriptors  of
both images.

  In the SIFT method, 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 we use the parameterization of the
gaussian by its standard deviation $\sigma$,
$$G_\sigma(x,y)=G(\sigma, x, y)=\frac{1}{2\pi\sigma^2} e^{-(x^2+y^2)/2 \sigma^2}.$$

Taking apart all sampling issues and several thresholds whose aim it is to eliminate unreliable features, the whole
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 space Laplacian 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  $\sqrt{\sigma_i^2+\bc^2}$. }

The constant $\bc\simeq 0.8$ is the tentative standard deviation of the 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.eps}
\includegraphics[width=3cm]{figures/pise-nfa001.eps}
 \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}
 Figs \ref{SIFTexample} and
 \ref{fig:sift:scale:inv} show  striking examples of the method scale invariance.
 Lowe claims that 1) his descriptors are invariant with respect to translation, scale and
rotation, and that 2) they provide a robust matching across a substantial range of affine distortions, change in 3D
viewpoint, addition of noise, and change in illumination. In addition, being local, they are robust to occlusion. Thus
they match all requirements for shape recognition algorithms except one: they are not really affine invariant but only
robust to moderate affine distortions.


\section{A Short Guide to SIFT Encoding}
\label{sec:guide}


\begin{figure}
\centering
\includegraphics[width=9cm]{fig1Lowe04.eps}
\caption{Gaussian pyramid for key points extraction (from \cite{lowe04:_distin})} \label{fig1lowe}
\end{figure}

The SIFT encoding algorithm consists of four steps: detection of scale-space extrema (Sect.~\ref{scalespace}), accurate
localization of key points (Sect.~\ref{accurate}), and descriptor construction (Sect.~\ref{descriptor}).


\subsection{Scale-Space Extrema}
\label{scalespace}

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.  We recall that the scale-space representation introduces a smoothing
parameter $\sigma$, the scale, and convolves the image with Gaussian functions of increasing standard deviation
$\sigma$.
  By a classical approximation inspired from psychophysics~\cite{marr82vision}, the
Laplacian of the Gaussian is replaced by a Difference of Gaussians at different scales (DOG). Extrema of the Laplacian
are then replaced by extrema of DOG functions: ${\mathbb D}(\sigma,x,y) =  w(k\sigma,x,y)- w(\sigma,x,y) ,$ where~$k$
is a constant multiplicative factor. Indeed, it is easy to show that ${\mathbb D}(\sigma,x,y)$ is an approximation of
the
 Laplacian:
$${\mathbb D}(\sigma,x,y) \approx (k-1)\sigma^2 (\Delta G_\sigma  \ast  u_0)(x,y).$$

In the terms of David Lowe:

\begin{quotation}
{\em The factor $(k-1)$ in the equation is constant over all scales and  therefore does not influence extrema location.
The approximation error will go to zero as  k goes to 1, but in practice we have found that the approximation has
almost no impact on the stability of extrema detection or localization for even  significant differences in scale, such
as $k=\sqrt{2}$.}
\end{quotation} To be more specific, quoting Lowe again:

\begin{quotation}
\[
{\mathbb  D}(\sigma, x, y)=:(G(k\sigma, x, y) - G(\sigma, x, y)) \ast  u_0(x, y)=  w(k\sigma, x, y) -  w(\sigma, x, y)
\]

{\em The relationship between $D$ and $\sigma^2 \Laplacian G$ can be  understood from the heat diffusion equation
(parameterized in terms of $\sigma$ rather than the more usual $t=\sigma^2$):}
\[
\frac{\partial G}{\partial \sigma} = \sigma \Laplacian G.
\]

{\em From this, we see that $\Laplacian G$ can be computed from the finite  difference approximation to $\partial G /
\partial \sigma$, using the difference  of nearby scales at~$k \sigma$ and~$\sigma$: }
\[
\sigma \Laplacian G = \frac{\partial G}{\partial \sigma} \approx \frac{G(k\sigma, x, y)-G(\sigma, x, y)}{k \sigma -
\sigma}
\]

{\em and therefore, }
\[
 G(k\sigma, x, y)-G(\sigma, x, y) \approx (k-1) \sigma^2 \Laplacian G.
\]

{\em This shows that when the difference-of-Gaussian function has scales differing by a constant factor it already
incorporates the $\sigma^2$ scale normalization required for the scale-invariant Laplacian. }

\end{quotation}
This leads to an efficient computation of local extrema of~$\mathbb  D$ by exploring neighborhoods through a Gaussian
pyramid ; see Figs.~\ref{fig1lowe} and~\ref{fig2lowe}.

\begin{exercise}
Show that the gaussian $G_\sigma$ parameterized by its standard deviation $\sigma$ satisfies as stated by Lowe the
time-dependent heat equation $\frac{\partial G}{\partial \sigma} = \sigma \Laplacian G.$ \finex
\end{exercise}



\begin{figure}
\centering
\includegraphics[width=5cm, height=3cm]{fig2Lowe04.eps}
\caption{Neighborhood for the location of key points (from \cite{lowe04:_distin}). Local extrema are detected by
comparing each sample point in $\mathbb D$ with its eight neighbors at scale $\sigma$ and its nine neighbors in the
scales above and below.} \label{fig2lowe}
\end{figure}


\subsection{Accurate Key Point Detection}
\label{accurate}

In order to achieve sub-pixel accuracy, the interest point  position is slightly corrected thanks to a quadratic
interpolation. Let us call ${\mathbf x}_0=: (\sigma_0, x_0, y_0)$ the current detected point in scale space, which is
known up to the (rough) sampling accuracy in space and scale. Notice that all points ${\mathbf x} =(\sigma, x, y)$ here
are scale-space coordinates. Let us call ${\mathbf x}_1={\mathbf x}_0+ {\mathbf y}$ the real extremum of the DOG
function. The Taylor expansion of~${\mathbb D}$ yields
$${\mathbb D}({\mathbf x}_0+{\mathbf y}) = {\mathbb D}({\mathbf x}_0)+
\left(D {\mathbb D}\right)({\mathbf x}_0)\cdot {\mathbf y} + \frac12\left( D^2
  {\mathbb D}\right)({\mathbf x}_0)({\mathbf y},{\mathbf y})  +o(\|{\mathbf y}\|^2), $$
where ${\mathbb D}$ and its derivatives are evaluated at an interest point and
  ${\mathbf y}$ denotes an
  offset from this point.
Since interest points are extrema of ${\mathbb D}$ in scale space, setting the
  derivative to zero gives:
\begin{equation} \label{Newtonkeypoint}{\mathbf y} = - \left( D^2{\mathbb D}({\mathbf x}_0)\right)^{-1}
 \left(D{\mathbb D}({\bf x}_0)\right),\end{equation}
  which is the sub-pixel correction for a more accurate position
  of the key point of interest.

\begin{exercise} Check that \eqref{Newtonkeypoint} is a point where the gradient of ${\mathbb D}$ vanishes.
\finex
\end{exercise}
%If the offset is larger than 0.5 pixel, then the extremum lies closer
%to another sample point; in this case the sample point is changed.

Since points with low contrast are sensitive to noise, and since points that are poorly localized along an edge are not
reliable, a filtering step is called for. Low contrast points are handled through a simple thresholding step. Edge
points are swept out following the Harris and Stephen's interest points paradigm. Let~$H$ be the following Hessian
matrix:
$$H=\left(
\begin{array}{cc}
{\mathbb D}_{xx} & {\mathbb D}_{xy} \\
{\mathbb D}_{xy} & {\mathbb D}_{yy}
\end{array}
\right).$$ The reliability test is simply to assess whether the ratio between the
 larger eigenvalue and the smaller one is below a threshold~$r$.
This amounts to check: \begin{equation}\label{SIFTthresholdtwo} \frac{\text{Tr}(H)^2}{\text{Det}(H)} <
\frac{(r+1)^2}{r}.
\end{equation}
This rules out standard edge points and puts points of interest at locations which are strong enough
extrema, or saddle points.
\begin{exercise}
Explain why \eqref{SIFTthresholdtwo}  is equivalent imposing that the ratio between the smaller eigenvalue and the
larger eigenvalue of $H$ is smaller than  $r$. These eigenvalues are assumed to have the same sign. Why?
 \finex \end{exercise}


\subsection{Construction of the SIFT descriptor\label{descriptor}}
  In order to extract rotation-invariant patches, an orientation must be assigned to each key point.
Lowe proposes to estimate a semi-local average orientation for each key point. From each sample image~$L_\sigma$,
gradient magnitude and orientation is pre-computed using a $2 \times 2$ scheme. An orientation histogram is assigned to
each key point by accumulating gradient orientations weighted by~1) the corresponding gradient magnitude and by~2) a
Gaussian factor depending on the distance to the considered key point and on the scale. The precision of this histogram
is 10 degrees. Peaks simply correspond to dominant directions of local gradients. Key points are created for each peak
with similar magnitude, and the assigned orientation is refined by local quadratic interpolation of the histogram
values.


\begin{figure}
\centering
\includegraphics[width=9cm]{MuralKeypoints.eps}
\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 covariant to any image
similarity.\label{fig.SIFTkeypoints}}
\end{figure}

Once a scale and an orientation are assigned to each key point, 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}. The next step is
to extract from this patch {\em robust} information. Gradient samples are accumulated into orientation histograms
summarizing the contents over $4 \times 4$ subregions
 surrounding
the key point of interest. Each of the 16 subregions corresponds to a 8-orientations bins histogram, leading to a 128
element feature for each key point (see Fig. \ref{fig7lowe}). Two modifications are made in order to reduce the effects
of illumination changes: histogram values are thresholded to reduce importance of large gradients (in order to deal
with a strong illumination change such as camera saturation), and feature vectors are normalized to unit length (making
them invariant to affine changes in illumination).

\subsection{Final matching}
The outcome is for each image, a few hundreds or thousands SIFT descriptors associated with as many key points. The
descriptors of any image can be compared to the descriptors of any other image, or belonging to a database of
descriptors built up from many images. The only remaining question is to decide when two descriptors match, or not. In
the terms of Lowe again:

\begin{quotation} The best candidate match for each keypoint is found by identifying its nearest neighbor in the database of
keypoints from training images. The nearest neighbor is defined as the keypoint with minimum Euclidean distance for the
invariant descriptor vector. However, many features from an image will not have any correct match in the training
database because they arise from background clutter or were not detected in the training images. Therefore, it would be
useful to have a way to discard features that do not have any good match to the database. A global threshold on
distance to the closest feature does not perform well, as some descriptors are much more discriminative than others. A
more effective measure is obtained by comparing the distance of the closest neighbor to that of the second-closest
neighbor. (...) This measure performs well because correct matches need to have the closest neighbor significantly
closer than the closest incorrect match to achieve reliable matching. For false matches, there will likely be a number
of other false matches within similar distances due to the high dimensionality of the feature space. We can think of
the second-closest match as providing an estimate of the density of false matches within this portion of the feature
space and at the same time identifying specific instances of feature ambiguity. (...) For our object recognition
implementation, we reject all matches in which the distance ratio is greater than 0.8, which eliminates 90$\%$ of the
false matches while discarding less than 5$\%$ of the correct matches.
\end{quotation}

\begin{figure}
\centering
\includegraphics[width=9cm, height=4cm]{fig7Lowe04.eps}
\caption{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. (From \cite{lowe04:_distin}) } \label{fig7lowe}
\end{figure}


  \section{Image acquisition  model underlying SIFT}


\subsection{The camera model}
We always work on the camera CCD plane, whose mesh unit is taken to be 1. We shall always assume that the camera pixels
are indexed by $\relative^2$. The image sampling operator is therefore always $\bS_1$. Our second assumption is that
the digital initial image is well-sampled and obtained by a gaussian kernel. Thus, the digital image is $\bu=\bS_1
G_\delta A u_0$, where $\delta \geq \ft2$, $\ft2\simeq 0.6$ ensures well-sampling (see Chapter \ref{blurassumption}),
and $A$ is a similarity with positive determinant. (In fact Lowe's original paper assumes $\ft2\simeq 0.5$, which
amounts to assume a slight under-sampling of the original image).


\begin{definition}  We model all digital frontal images obtained from a
 given ideal planar object whose frontal infinite
resolution image is $ u_0$ as
\begin{equation}
\label{imageformationmodel}\bu_0=:\bS_1 G_\delta A u_0
\end{equation}
where $\delta\geq \ft2$ and $ A$ is a $A=RH_\lambda\mathcal T$ is the composition of a translation and of a similarity.
\end{definition}
So the possibility of aliasing (under-sampling, $\delta < \ft2$ is discarded). Taking into account the way the digital
image  is blurred and sampled in the SIFT method, we can now
 list the
 SIFT assumptions and formalize the method itself. The description is by far simpler if we  do it without mixing in
 sampling issues. We need not mix them in, since the fact that images are well-sampled at all stages permits
 equivalently to describe all operations with the continuous  images directly, and to deduce afterwards the discrete
 operators on samples. We refer to section \ref{discretesetting} for this passage from continuous to discrete
 operations
 in the well-sampled world.



\subsection{Condensed description of the SIFT method\label{condenseddescription}}
\begin{enumerate}
\item There is an underlying infinite resolution bounded planar image $ u_0$;
\item  The initial digital image is $\bS_1 G_\delta A u_0$
 where $\delta\geq \ft2$,
and $ A=  R H_\lambda\mathcal T$ is the composition of a rotation, a zoom,  and a translation;
\item  the SIFT method computes a sufficient scale-space sampling of $ u(\sigma, \fxx)=
 (G_\sigma  G_\delta A u_0)(\fxx)$, and deduces by the Newton method the accurate location or key points defined as extrema in
 scale-space of the spatial image Laplacian, $\Delta u(\sigma ; \fxx)$;
\item The blurred $ u(\sigma, \cdot)$  image is then re-sampled around each characteristic point with sampling mesh
$\sqrt{\sigma^2 +\ft2^2}$;
\item the directions of the sampling axes are fixed by a dominant direction of the  gradient of $ u(\sigma,\cdot)$ in a
neighborhood, with size proportional to $\sqrt{\sigma^2+\ft2^2}$  around the characteristic  point;
\item the rest of the  operations in the SIFT method is a contrast invariant encoding of the samples around each
characteristic point. It is not needed for the discussion to follow.
\end{enumerate}







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

 In this section, in conformity with the SIFT model of Sect. \ref{condenseddescription}, 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 homothety $ H$, a translation $\mathbf{\mathcal T}$ and a rotation $ R$.
 Thus the digital image is  $\bu=\bS_1  G_\delta  H\mathcal T R  u_0$,
  for some $ H$, $\mathcal T$,  $ R$ as above. Assuming that the image is not aliased boils down, by the experimental
  results of Sect. \ref{blurassumption}, to assuming $\delta\geq 0.8$.


 Consider $\mathbf{\mathcal T}$ an arbitrary  image translation,
   $R$ an arbitrary image rotation,   $H_\lambda$ an
arbitrary image homothety, $G$ an arbitrary gaussian convolution, $D$ the gradient  and $\Laplacian$ the Laplacian, all
applied to continuous images. We say that there is strong commutation of two of these operators if we can exchange the
order of their application to any image. We say that there is weak commutation between two of these operators if we can
exchange their order by changing one of the parameters of one of the operators. For example we have $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, with the following exceptions:
 $R$ and $ G$ commute
strongly, $DH_\lambda=\lambda H_\lambda D$, $\Laplacian H_\lambda=\lambda^2H_\lambda\Laplacian$, and $D$ and
$\Laplacian$ do not commute.
\end{lemma}

\begin{exercise} Check all of the mentioned commutations and give their exact formula. There are six
 kinds of operators: translations,
rotations, homotheties, gaussian convolutions, gradients, and Laplacians. Thus, there are 15 verifications to make.
\finex
\end{exercise}

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

 \begin{proof}
By the weak commutation of translations and rotations with all other operators (Lemma \ref{weakcommutation}), 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 space translations and rotations. The normalization process for each SIFT descriptor
situates the origin at each extremum in turn, thus canceling the translation. 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 $u= G_\beta  H_\lambda u_0$ and $v= G_\delta  H_\mu u_0$ 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+\bc^2}$. 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=: I\bS_1  G_\beta  H_\lambda u_0 =   G_\beta  H_\lambda u_0 \hbox{ and }
v=:I\bS_1  G_\delta H_\mu u_0=  G_\delta H_\mu  u_0.$$ Computing the scale-space of these images amounts to convolve
these images for every $\sigma>0$ with $ G_\sigma,$ which yields, using the commutation relation
\eqref{commutationconvolutionzoomcontinus} and the semigroup property \eqref{continuoussemigroupproperty}:
$$ 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.$$
Set $w(s,\bx)=:  G_s  u_0$. 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 at scale $\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}
\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 from arbitrary distances,
with possibly different camera gaussian blurs,  with an arbitrary camera translation parallel to its focal plane, and
an arbitrary rotation around its optical axe.  Without loss of generality, assume $\lambda\leq \mu$. Then if the camera
blurs are standard ($\beta=\delta=\bc$), each SIFT descriptor of $u=I\bu$ is identical to some SIFT descriptor of
$v=I\bv$. If $\beta\neq \delta$ (or $\beta=\delta \neq \bc$), the SIFT descriptors of $u$ become (quickly) similar to
SIFT descriptors of $v$ 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$.
\end{theorem}

\begin{proof} By the result of Lemma \ref{invariancesiftrotationtranslation}, we can neglect the  effect of translations and
rotations. Therefore assume w.l.o.g. that the images under comparison are as in Lemma
\ref{approximatescaleinvarianceofSIFT}. Assume  a key point $(s_0, \bx_0)$ of $ w$ has scale $s_0\geq
\max(\lambda\beta, \mu\delta).$ This key point has a sampling rate proportional to $s_0$. There is a corresponding key
point $(\sigma_1, \frac{\bx_0}\lambda)$ for $ u$ with sampling rate $\sqrt{\sigma_2^2+\bc^2}$ and a corresponding key
point $(\sigma_2, \frac{\bx_0}\mu)$ with sampling rate $\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
 $ u$ 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 $ u$ and $v$ for $\bx_0$ will be identical  if and only if
$\lambda\sqrt{\sigma_1^2+\bc^2}=\mu\sqrt{\sigma_2^2+\bc^2}$.  Now, we  have
$\lambda\sqrt{\sigma_1^2+\beta^2}=\mu\sqrt{\sigma_2^2+\delta^2}$, which implies
$\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$ are proportional to camera distances to the observed object $ u_0$, they are arbitrary and
generally different. Thus, the only way to ensure \eqref{conditiondescriptorsequal} is to have $\beta=\delta=\bc$,
which means that the blurs of both images (or of both cameras) are ideal and gaussian. In any case, $\beta=\delta=\bc$
does imply that the SIFT descriptors of both images are identical.

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}\simeq \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 28 SIFT key points of a very small digital image $\bu$ are compared to the 86 key
points obtained by zooming in $\bu$ by a 32 factor:  The resulting digital image is the digital image
$\bv=\bS_{\frac{1}{32}} I\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 (22 out of 28) of $\bu$ finds its counterpart in $\bv$. 22 matches are detected between the descriptors as
shown on the right. If we trust Theorem \ref{theo:scale}, all descriptors of $\bu$ should have been retrieved in $\bv$.
This does not fully happen for two reasons.  First, the SIFT method thresholds (not taken into account in the theorem)
eliminate many potential key points. Second, the zero-padding interpolation giving $\bv$ is imperfect near the image
boundaries.

\begin{figure}[htbp]
\begin{center}
\begin{tabular}{ccc}
\includegraphics[width=4cm]{./figures/sift_scale_inv_im1.eps}
&
\includegraphics[width=4cm]{./figures/sift_scale_inv_im2.eps}
&
\includegraphics[width=4cm]{./figures/sift_scale_inv.eps}
\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 28 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 $\bv=\bS_{\frac{1}{32}} I\bu$. It has 86
key points. Right: 22 matches found between $ \bu$ and $\bS_{\frac{1}{32}} I\bu$. } \label{fig:sift:scale:inv}
\end{figure}

\section{Exercises}

\begin{exercise}
The aim of the exercise is to explain why the experiment of Fig. \ref{fig:sift:scale:inv} works, and to illustrate
Theorem \ref{theo:scale}. The digital zoom in by a factor $\lambda$ is nothing but the discrete over-sampling operator
$\SAD_{\frac 1\lambda}$ with sampling step $\frac 1 \lambda$,  defined in \eqref{defdiscretesampling}. Here, $\lambda =
32$.  In the experiment an original digital image $\bu=\bS_1G_\delta u$ is zoomed into $\bv=\SAD_{\frac
1{\lambda}}\bu$.

\skipaline 1) Using the definition of the discrete zoom and the right commutation relations given in this chapter and
in the former one (give their numbers), show that
$$\bv=\bS_{\frac 1\lambda} G_\delta  I\bu=\bS_1 G_{\lambda\delta} H_{\frac 1\lambda}I\bu.$$

\skipaline 2) Is $\bv$ well-sampled if $\bu$ was?

\skipaline 3) By applying carefully Theorem \ref{theo:scale}, assuming that $\delta\simeq \ft2$, discuss why SIFT
manages to match SIFT descriptors of $\bu$ and $\bv$.
 \finex
\end{exercise}

 \section{Comments and references}


  Many variations exist on the computation of 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 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, cao:07}
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}, no known detector 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 actually 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 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). 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 registration~\cite{single-alignment}, 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}.


 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 real invariance of each method when all thresholds have been
eliminated. This real invariance has been proved here for SIFT. The second criterion is the practical validity of the
sampling method used in SIFT, that has been again checked in Chapter \ref{blurassumption}. 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 very 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.






\ifodd\count0
    \clearpage
    \leavevmode
    \thispagestyle{empty}
    \clearpage
\else\relax \fi
