\documentclass[11pt]{article}
\usepackage{epsfig}





\renewcommand{\baselinestretch}{1.2}


\setlength{\topmargin}{-0.5in}
\setlength{\oddsidemargin}{-0.00in}
\setlength{\textwidth}{6.5in}
\setlength{\textheight}{9in}

\newtheorem{theor}{Theorem}
\newtheorem{Prop}{Proposition}


\newcommand{\Local}{./img_eps}

\def\I{{\bf I}}
\def\J{{\bf J}}
\def\B{{\bf B}}
\def\bB{\bar{{\bf B}}}
\def\bC{\bar{{C}}}
\def\S{{\rm SUM}}
\def\M{{\rm MAX}}
\def\OB{\bar{{B}}}
\def\bc{\tilde{c}}
\def\bp{\tilde{p}}
\def\tr{\bar{r}}
\def\rr{\gamma}
\def\E{{\rm E}}
\def\N{{\rm N}}
\def\KL{{\rm KL}}
\def\Var{{\rm Var}}

\begin{document}




\begin{center}
{\Huge Active Basis for Modeling, Learning and Recognizing
Deformable Templates}
\end{center}

\begin{center}
Ying Nian Wu$^1$, Zhangzhang Si$^1$, Haifeng Gong$^{1,2}$, and
Song-Chun Zhu$^{1,2}$\\
\vspace{3mm}
\noindent $^1$ Department of Statistics, University of California, Los Angeles\\
$^2$ Lotus Hill Research Institute, Ezhou, China\\
$\{$ywu, zzsi, hfgong, sczhu$\}$@stat.ucla.edu\\
\vspace{3mm} \small{Revised September 25, 2008}
\end{center}




\begin{center}
{\bf Abstract}
\vspace{1cm}

\parbox{5.8in}
{\setlength{\baselineskip}{0.6cm} This article proposes an active
basis model, a shared sketch algorithm, and a computational
architecture of sum-max maps for representing, learning, and
recognizing deformable templates. In our generative model, a
deformable template is in the form of an active basis, which
consists of a small number of Gabor wavelet elements at selected
locations and orientations. These elements are allowed to slightly
perturb their locations and orientations before they are linearly
combined to generate the observed image. The active basis model, in
particular, the locations and the orientations of the basis
elements, can be learned from training images by the shared sketch
algorithm. The algorithm selects the elements of the active basis
sequentially from a dictionary of Gabor wavelets at a dense
collection of locations and orientations. When an element is
selected at each step, the element is  shared by all the training
images, and the element is perturbed to encode or sketch a nearby
edge segment in each training image. The recognition of the
deformable template from an image can be accomplished by a
computational architecture that alternates the sum maps and the max
maps. The computation of the max maps deforms the active basis to
match the image data, and the computation of the sum maps scores the
template matching by the log-likelihood of the deformed active
basis.}
\end{center}

\vspace{1cm} {\bf Keywords}: Generative model; Object recognition; Shared sketch algorithm; Sum maps and max maps; Wavelet sparse coding.

\vspace{5mm}


{\small This paper is submitted to the International Journal of
Computer Vision (IJCV). A shorter version has been published in ICCV
2007 \cite{Wu_ICCV07}. For reproducibility, the matlab/C codes and
data used in this paper are available at
\underline{www.stat.ucla.edu/$\sim$ywu/ActiveBasis.html}}

\newpage



%%%%%%%%% BODY TEXT
\section{Introduction}

Deformable template \cite{Yuille} is an important element in object recognition \cite{FisEl}. In this article, we propose a generative model, a model-based algorithm, and a computational architecture for representing, learning and recognizing deformable templates.


\subsection{Form of representation}

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=ActiveBasis.eps,height=2in}
\end{tabular}
\end{center}
\caption{Active basis. Each basis element is illustrated by a thin
ellipsoid at a certain location and orientation. The upper half
shows the perturbation of one basis element. By shifting its
location or orientation or both within a limited range, the basis
element (illustrated by a black ellipsoid) can change to other Gabor
wavelet elements (illustrated by the blue ellipsoids).  } \label{fig:ActiveBasis}
\end{figure}

We call our model the active basis model. An active basis consists
of a small number of Gabor wavelet elements at selected locations
and orientations. These elements are allowed to slightly perturb
their locations and orientations before they are linearly combined
to generate the observed image. Figure (\ref{fig:ActiveBasis})
illustrates the basic idea. The lower half of Figure
(\ref{fig:ActiveBasis}) shows an active basis, where each element is
illustrated by a thin ellipsoid at a certain position and with a
certain orientation. The upper half of  Figure (\ref{fig:ActiveBasis}) illustrates the perturbation of one basis
element. Intuitively, each Gabor wavelet element can be considered
a ``stroke.'' The template is formed by a composition of a number of strokes. These strokes can be slightly perturbed, so that the template is deformable.



\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=ex1_0.eps,height=1.5in}
\end{tabular}
\end{center}
\caption{Active basis formed by 60 Gabor wavelet elements. The first
block displays the 60 elements, where each element is represented by
a bar. For each of the other 7 blocks, the left plot is the observed
image, and the right plot displays the 60 Gabor wavelet elements
resulting from locally shifting the 60 elements in the first block
to fit the corresponding observed image. }
\label{fig:ActiveBasisReal}
\end{figure}



Figure (\ref{fig:ActiveBasisReal}) shows a real example. It displays
7 images of cars at the same scale and in the same pose. These
images are defined on a common image lattice, which is the bounding
box of the cars. These images are represented by an active basis
consisting of 60 Gabor wavelet elements, as displayed in the first
block of Figure (\ref{fig:ActiveBasisReal}). Each wavelet element is
represented symbolically by a bar at the same location and with the
same length and orientation as the wavelet element. The length of
each element is about 1/10 of the length of the image patch. These
elements do not have much overlap and are well connected. They form
a common template or an average sketch of the training image
patches. The 60 elements of the active basis in the first block of
Figure (\ref{fig:ActiveBasisReal}) are allowed to locally change
their locations and orientations to code each observed image, as
illustrated by the remaining 7 blocks of Figure
(\ref{fig:ActiveBasisReal}). Within each block, the left plot
displays the observed car image, and the right plot displays the 60
Gabor wavelet elements that are actually used for encoding the
corresponding observed image. They form the deformed active basis
that sketches the observed image.

\subsection{Scheme of learning}

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=SharedSketch.eps,height=3in}
\end{tabular}
\end{center}
\caption{Shared sketch algorithm. A selected element (colored
ellipsoid) is shared by all the training images. For each image, a
perturbed version of the element seeks to sketch a local edge segment near the element by a local maximization operation. The elements of the active basis are selected sequentially according to the Kullback-Leibler divergence between the pooled distribution (colored solid curve) of filter responses and the background distribution (black dotted curve). The divergence can be simplified into a pursuit index, which is the sum of the transformed filter responses. The sum essentially counts the number of edge segments sketched by the perturbed versions of the element.  }
\label{fig:SharedSketch}
\end{figure}

The active basis, in particular, the locations and the orientations
of the basis elements, can be learned from training image patches by
the shared sketch algorithm. The algorithm selects the elements of
the active basis sequentially from a dictionary of Gabor wavelets at
a dense collection of locations and orientations. Figure
(\ref{fig:SharedSketch}) illustrates the selection of three elements
by learning from a sample of training images of cars. When an
element is selected, the element is shared by all the training
images in the sense that a perturbed version of this element is
added to improve the encoding of each image. Specifically, the
element is perturbed to a location and orientation that achieves the local maximum response within a small neighborhood of the selected element, that is, the perturbed version of the selected element seeks to sketch a nearby edge segment in each training image. For instance, when the green element is selected, it is attracted to the nearby edge in each training image. The same is true for the red element and the blue element.

For each element, a distribution of filter responses is pooled over
all the training images at the perturbed locations and orientations.
The elements are selected in an order according to the
Kullback-Leibler divergence between the pooled distribution (solid curve) and a background distribution (dotted curve). The background
distribution is pooled over natural images. With proper  parametrization, the Kullback-Leibler divergence is equivalent to a pursuit index that drives the selection of the elements. This index takes the form of the sum of the transformed filtered responses, summed over all the training images. The transformation is an increasing function that discounts large filter responses. So the pursuit index can be interpreted as a voting of the training images, and the index favors the element whose perturbed versions sketch as many edge segments as possible. After an element is selected, its perturbed version explains away a small part of each training image, and thereby inhibits nearby Gabor wavelet elements from coding the same part of the image. So the selected elements of the active basis are well spaced, and usually form a clear template.

The active basis displayed in Figure (\ref{fig:ActiveBasisReal}) is
learned by the shared sketch algorithm. It is worth noting that for
the last two examples in Figure (\ref{fig:ActiveBasisReal}), the
strong edges in the background are not sketched, because these edges
are not shared by other examples, and such edges are ignored by the
shared sketch algorithm.

\subsection{Architecture of inference}

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=SumMaxmaps.eps,height=4in}
\end{tabular}
\end{center}
\caption{Sum-max maps. The $\S1$ maps are obtained by convolving the
input image with Gabor filters at all the locations and
orientations. The ellipsoids in the $\S1$ maps illustrate the local
filtering or summation operation. The $\M1$ maps are obtained by
applying a local maximization operator to the $\S1$ maps. The arrows
in the $\M1$ maps illustrate the perturbations over which the local
maximization is taken. The $\S2$ maps are computed by applying a
local summation operator to the $\M1$ maps, where the summation is
over the elements of the active basis.  This operation computes the
log-likelihood of the deformed active basis, and can be interpreted
as a shape filter. } \label{fig:SumMaxMaps}
\end{figure}

After learning the active basis from training images, the detection and recognition of the deformable template from a testing image can be accomplished by a computational architecture of sum-max maps. Figure (\ref{fig:SumMaxMaps}) illustrates this architecture. It starts from convolving the image with Gabor filters at all the locations and orientations. The filtered images become the first layer of the sum maps, or $\S1$ maps, because each Gabor filter is a local summation operator. In Figure (\ref{fig:SumMaxMaps}), the thin ellipsoids in the $\S1$ maps illustrate the local filtering or summation operation. Then a layer of max maps, or $\M1$ maps, is computed by applying a local maximization operator to the $\S1$ maps. In Figure (\ref{fig:SumMaxMaps}), the arrows in the $\M1$ maps illustrate that the local maximization is taken over small perturbations of the Gabor wavelets. This local maximization tells us how to deform the active basis to match the image data.

On top of that, a sum map, or $\S2$ map, is computed by applying a local summation operator to the $\M1$ maps. Specifically, we scan the active basis over the whole image lattice, and for each pixel of the $\S2$ map, we compute a weighted sum of the values of the $\M1$ maps, where the summation is over the locations and orientations of the elements of the active basis centered at this pixel. So this is another layer of filtering operation, and can be considered a shape filter. It computes the log-likelihood of the deformed active basis. In Figure (\ref{fig:SumMaxMaps}), the car template in the $\S2$ map illustrates the active basis centered at one pixel. We scan this template over all the pixels to obtain the $\S2$ map, which scores the template matching.

The $\S2$ map is obtained by a local summation operator of fixed shape. However, because the local summation is applied to the $\M1$ maps, shape deformation is automatically accounted for, and the template matching score is invariant to shape deformation.

Besides the log-likelihood scoring for template matching, we also develop a non-probabilistic scoring method based on active correlation between the template and the image.

\subsection{Review of literature}

Our work cultivates key insights from three major theories on
biological and computer vision, namely, Olshausen and Field's theory
of linear sparse coding using Gabor wavelets \cite{Olshausen_Field},
Riesenhuber and Poggio's theory on local maximum pooling of Gabor
filter responses \cite{Riesenhuber_Poggio}, and Viola and Jones's
theory of adaboost learning \cite{adaboost} using Harr wavelets as
weak classifiers \cite{Viola_Johnes}.

Our active basis model is based on Olshausen and Field's linear
representation. Inspired by Riesenhuber and Poggio's operator, we
add local perturbations to the basis elements. Motivated  by Viola
and Jones's  work, we apply the model to images of object shapes, so that the training images share the same set of selected elements, which form a common deformable template \cite{Yuille}. This connects Olshausen and Field's work to shape models such as active contours \cite{Kass_snakes} and active appearance model \cite{Cootes}.

Our algorithmic architecture of sum-max maps is a variation on the
theme of Riesenhuber and Poggio's cortex-like hierarchical structure
of simple and complex cells \cite{Riesenhuber_Poggio}. See also the
recent work of Mutch and Lowe \cite{Lowe} for further improvement.
Our architecture is essentially the same as Riesenhuber and Poggio's
structure up to the $\M1$ maps. However, our operator for computing
the $\S2$ map for template matching is different from the template matching scheme of Riesenhuber and Poggio. Our operator is a local summation over highly selected locations and orientations. Moreover, the computation is not entirely a bottom-up process. After the bottom-up scoring process, a top-down retrieving process traces the locations and orientations of the perturbed elements, which form the deformed template that is matched to the input image. Furthermore, since the local summation is taken over a small number of selected locations and orientations, it may be more efficient than comparing all the intensities of the $\M1$ maps.

Our shared sketch algorithm is similar to the learning scheme of
Viola and Jones's  version of adaboost. However, the algorithm is
guided by a generative model, where the selection of the basis
elements is based on explaining away the image data instead of
fitting the classification boundary. Our learning algorithm can be
considered a parallel version of the matching pursuit algorithm
\cite{Mallat_Zhang} where we select the basis elements to
simultaneously code all the training images. It can also be
considered a variation of the projection pursuit algorithm for
density estimation \cite{Friedman}, where we add to it the local
maximization and local inhibition operations of the matching pursuit
algorithm. The exponential family model we adopt is related to the feature induction schemes of Della Piatra et al. \cite{DellaPietra} and Zhu et al. \cite{Zhu_FRAME}.

Our work is a continuation of our long term search for generative
models and model-based algorithms \cite{Zhu_FRAME, Wu_equivalence,
Zhu_texton, Wu_texture_sketch, GuoZhuWu, GuoZhuWu2}, as well as our attempt to understand these models within a common information-theoretical framework \cite{Wu}. The active basis model can be considered a revision of our previous model on textons \cite{Zhu_texton}. It can also be viewed as an inhomogeneous version of the Markov random field model that we previously developed for textures \cite{Zhu_FRAME}, as we show in \cite{Wu}. More important, the active basis model is a simplest instance of the and-or graph \cite{Zhu_grammar} in the compositional framework \cite{Geman_composition} that we have been studying. The and-or grammar naturally suggests that we can further compose multiple active bases to represent more articulate shapes.





%------------------------------------------------------------------------
\section{Representation, Learning, and Inference}

This section presents the active basis representation, and describes the shared sketch algorithm and the sum-max maps with pseudo-codes.  We leave theoretical underpinnings and justifications to the next section.

\subsection{Gabor wavelets and sparse coding}

\indent{\em A dictionary of Gabor wavelets}.
A Gabor function \cite{Daugman} is of the form
\[
G(x, y) \propto
\exp\{ -[(x/\sigma_x)^2 + (y/\sigma_y)^2]/2\}
e^{i {x}}.
\]
We can translate, rotate, and dilate $G(x, y)$ to obtain a general form of Gabor wavelets:
\[
B_{x, y, s, \alpha}(x', y') =
G({\tilde{x}}/{s}, {\tilde{y}}/{s})/s^2,
\]
where
\[
\tilde{x} =(x'-x)\cos\alpha - (y'-y)\sin\alpha,
\]
\[
\tilde{y} = (x'-x)\sin\alpha + (y'-y)\cos\alpha.
\]
$(x, y)$ is the central position, $s$ is the scale parameter, and $\alpha$ is the orientation. The Gabor wavelets give good fit to the receptive fields of the simple cells in V1 \cite{Daugman}.

The central frequency of $B_{x, y, s, \alpha}$ is $\omega = 1/s$. $B_{x, y, s, \alpha}$  =
($B_{x, y, s, \alpha, 0}$, $B_{x, y, s, \alpha, 1}$),
where $B_{x, y, s, \alpha, 0}$ is the even-symmetric Gabor cosine
component, and $B_{x, y, s, \alpha, 1}$ is the odd-symmetric Gabor
sine component. We always use Gabor wavelets as pairs of cosine and
sine components. We normalize both the Gabor sine and cosine
components to have zero mean and unit $\ell_2$ norm. For each $B_{x, y,
s, \alpha}$, $B_{x, y, s, \alpha, 0}$ and $B_{x, y, s, \alpha, 1}$
are orthogonal to each other.



\indent{\em Operation}.
Let $D$ be the domain of image lattice. The dictionary of Gabor
wavelet elements is ${\rm Dictionary} = \{B_{x, y, s, \alpha},
\forall (x, y, s, \alpha)\}, $ where $(x, y, s, \alpha)$ are densely
sampled: $(x, y) \in D$ with a fine sub-sampling rate (e.g., every 2
pixels), and $\alpha \in \{a\pi/A, a = 0, ..., A-1\}$ (e.g., $A =
15$).

For an image $\I$ defined on domain $D$, the projection coefficient of $\I$ onto $B_{x, y, s, \alpha, \eta}$ or the filter response is
\[
\langle \I, B_{x, y, s,
\alpha, \eta} \rangle = \sum_{x', y'} \I(x', y') B_{x, y, s,
\alpha, \eta}(x', y').
\]
 We write $\langle \I, B_{x, y, s, \alpha}
\rangle$ = ($\langle \I, B_{x, y, s, \alpha, 0} \rangle$, $\langle
\I, B_{x, y, s, \alpha, 1} \rangle$). The local energy is
\[
|\langle \I, B_{x, y, s, \alpha} \rangle|^2 =
 \langle \I, B_{x, y, s, \alpha, 0} \rangle^2
+\langle \I, B_{x, y, s, \alpha, 1}\rangle^2.
 \]
In the computation of $\langle \I, B_{x, y, s, \alpha} \rangle$, $B_{x, y, s, \alpha}$ is a linear filtering operator, which can serve as edge detector or local spectral analyzer.

To make filter responses comparable between different
training images, we need to normalize the images. Let
\begin{eqnarray}
   \sigma^2(s) = \frac{1}{|D|A} \sum_{\alpha} \sum_{(x, y) \in D}
   |\langle \I, B_{x, y, s, \alpha} \rangle|^2, \label{eq:normalize}
\end{eqnarray}
where $|D|$ is the number of pixels in $\I$, and $A$ is the total
number of orientations. For each input image $\I$, we normalize it
to $\I \leftarrow \I /\sigma(s)$.

\indent{\em Representation}.
A deeper perspective is offered by the sparse coding theory of Olshausen and Field \cite{Olshausen_Field}, where $B_{x, y, s, \alpha}$ serves as a representational element. Specifically, for an image $\I$, we can represent it by
\begin{eqnarray}
   \I = \sum_{i=1}^{n} c_i B_i + U, \label{eq:sc}
\end{eqnarray}
where $B_i = B_{x_i, y_i, s_i, \alpha_i}$, $c_i$ are the coefficients, and $U$ is the unexplained residual image. Recall that each $B_i$ is a pair of Gabor cosine and sine components. So $B_i = (B_{i, 0}, B_{i, 1})$. Accordingly, $c_{i} = (c_{i, 0}, c_{i, 1})$ and $c_{i} B_{i} = c_{i, 0} B_{i, 0} + c_{i, 1} B_{i, 1}$. The set of Gabor wavelet elements $\{B_i, i = 1, ..., n\}$ are selected from the dictionary. If the $\{B_i, i=1, ..., n\}$ are orthogonal, i.e., if they do not overlap in spatial domain or frequency domain, then $c_i = \langle \I, B_i \rangle$.

Sparse coding means that for a typical natural image $\I$, we can usually select a small number $n$ of elements from the dictionary, so that a linear combination of these elements can represent $\I$ with small residual $U$. Of course, for different images, we usually select different sets of elements. The wavelet sparse coding representation (\ref{eq:sc}) reduces an image of tens of thousands of pixels to a small number of wavelet elements or strokes. Using the sparse coding principle, Olshausen and Field \cite{Olshausen_Field} were able to learn from natural image patches a dictionary of wavelet elements that closely resemble the properties of the receptive fields of the simple cells in V1.

\subsection{Representation: active basis model}

The sparse coding model (\ref{eq:sc}) is constructed for the whole ensemble of natural images, where for different $\I$, we may represent them with completely different wavelet elements $(B_i, i=1, ..., n)$ with different $n$. In the active basis model, we apply the sparse coding model (\ref{eq:sc}) to image ensembles of various object categories. Then for each category, we require that the images share the same set of wavelet elements $(B_i, i=1, ..., n)$ , and these elements form a common template. However, when we use $(B_i, i=1, ..., n)$  to encode each individual image, we allow the template to slightly deform, by allowing the elements or strokes to perturb their locations and orientations.

Let $\{\I_m, m = 1, ..., M\}$ be a set of training image patches
defined on a common rectangle lattice $D$. We assume that $D$ is the
bounding box of the objects in $\I_m$, and these objects are from
the same category and in the same pose. We shall relax this
assumption later.

Our method is scale specific. We fix $s$ so that the length of
$B_{x, y, s, \alpha}$ (e.g., 17 pixels) is fixed. It is possible to learn multiple templates at multiple scales and then combine them.

The active basis model is a composition of strokes that are perturbable:
\begin{eqnarray}
  &&{\rm Composition:} \;\; \I_m = \sum_{i=1}^{n} c_{m, i} B_{m, i} + U_m,
  \label{eq:ab}\\
  && {\rm Perturbations:} \; \; B_{m, i} \approx B_{i}, \;  i = 1, ..., n, \label{eq:shift}
\end{eqnarray}
where $B_i \in {\rm Dictionary}$, $B_{m, i} \in {\rm Dictionary}$,
$(c_{m, i}, i = 1, ..., n)$ are the coefficients, and $U_m$ is the
unexplained residual image. To define the perturbation $B_{m, i} \approx B_{i}$,
suppose
\begin{eqnarray}
 && B_{i} = B_{x_i, y_i, s, \alpha_i}, \label{eq:ab1} \\
  && B_{m, i} = B_{x_{m, i}, y_{m, i}, s, \alpha_{m, i}},
\end{eqnarray}
then $B_{m, i} \approx B_{i}$ if and only if there exists $(d_{m,
i}, \delta_{m, i})$ such that
\begin{eqnarray}
  &&x_{m, i} = x_i + d_{m, i} \sin\alpha_i, \\
  && y_{m, i} = y_i + d_{m,
  i}\cos\alpha_i, \\
  && \alpha_{m, i} = \alpha_i + \delta_{m, i}, \\
  && d_{m, i} \in [-b_1,
b_1], \; \delta_{m, i} \in [-b_2, b_2]. \label{eq:ab2}
\end{eqnarray}
That is, we allow $B_i$ to shift its location along its normal
direction, and we also allow $B_i$ to shift its orientation. See
Figure  (\ref{fig:ActiveBasis}) for an illustration. We call $(d_{m,
i}, \delta_{m, i})$ the activity or perturbation of $B_i$ in image $\I_m$. $b_1$ and $b_2$ are the bounds for the allowed activities (e.g., $b_1 = 6$ pixels, and $b_2 =
\pi/15$).

In the above notation,  the active basis $\B = (B_i, i = 1, ..., n)$
forms a deformable template. The deformed active basis is $\B_m =
(B_{m, i}, i = 1, ..., n) \approx \B$. See Figure
(\ref{fig:ActiveBasisReal}) for an illustration.

Because we fix the scale $s$ in the linear representation
(\ref{eq:ab}) to (\ref{eq:ab2}), the linear superposition
$\sum_{i=1}^{n} c_{m, i} B_{m, i}$ only explains the frequency band
of $\I_m$ around the frequency $\omega = 1/s$, while leaving the
remaining frequency components to the unexplained $U_m$. $U_m$ can
be further explained by templates at other scales or resolutions.


\subsection{Learning: shared sketch algorithm
\label{sect:la}}

Given the set of training images $\{\I_m, m = 1, ..., M\}$, the shared sketch algorithm sequentially selects $B_i$ and perturbs it to $B_{m, i} \approx B_i$ to sketch each image $\I_m$. The basic idea is to select those $B_i$ so that its perturbed versions $\{B_{m, i}, m = 1, ..., M\}$ sketch as many edge segments as possible in the training images $\{\I_m\}$.

\vspace{.5cm} \noindent{\bf Description of the shared sketch algorithm}

\begin{itemize}
\item[Input]: Training images $\{\I_m, m = 1, ..., M\}$. 

\item[Output]: Template $\B = (B_i, i = 1, ..., n)$, and deformed template $\B_m = (B_{m, i}, i=1, ..., n)$ that is matched to $\I_m$ for $m = 1, ..., M$.
\end{itemize}

\begin{enumerate}
\item
Convolution: For each $m = 1, ..., M$, and for each $B \in {\rm Dictionary}$, compute
$[\I_m, B] = h(|\langle\I_m, B\rangle|^2)$.
Set $i \leftarrow 1$.

\item
Local maximization: For each putative candidate $B_i \in {\rm Dictionary}$, do the following: For each $m = 1, ..., M$, choose the optimal $B_{m, i}$ that maximizes $[\I_m, B_{m, i}]$ among all possible $B_{m, i} \approx B_i$.

\item
Selection: Choose that particular candidate $B_i$ whose  corresponding $\sum_{m=1}^{M} [\I_m, B_{m, i}]$ achieves the maximum among all possible $B_i \in {\rm Dictionary}$. Record this $B_i$ and retrieve the corresponding optimal $B_{m, i} \approx B_i$ for $m = 1, ..., M$.

\item
Non-maximum suppression: For each $m = 1, ..., M$, if $[\I_m, B_{m, i}]>0$, then for every $B \in {\rm Dictionary}$ such that
${\rm corr}(B, B_{m, i}) > \epsilon$, set $[\I_m, B] \leftarrow 0$.

\item
Stop if $i = n$. Otherwise let $i \leftarrow i + 1$, and go back to 2.
\end{enumerate}

In the above description, $h()$ is a monotone increasing (or non-decreasing) transformation that discounts large value of  $|\langle\I_m, B\rangle|^2$. For two Gabor elements $B_1$ and $B_2$, \[
{\rm corr}(B_1, B_2) = \sum_{\eta_1=0}^1\sum_{\eta_2=0}^1 \langle B_{1, \eta_1}, B_{2, \eta_2}\rangle^2
\]
measures their correlation or overlap in spatial and frequency domains. $B_1$ and $B_2$ are orthogonal as long as they do not overlap in either spatial domain or frequency domain.

See Figure (\ref{fig:SharedSketch}) for an illustration of the above algorithm. The algorithm can be considered a parallel version of edge detection. For a putative $B_i$, the local maximization step seeks to sketch a local edge segment in image $\I_m$ by a perturbed version $B_{m, i} \approx B_i$. The selection step seeks to find $B_i$ with the strongest $\sum_{m=1}^{M} [\I_m, B_{m, i}]$, which pools the edge strengths from the training images around $B_i$. After $B_i$ is selected, we retrieve the corresponding $B_{m, i}$, and let $B_{m, i}$ suppress or inhibit nearby overlapping Gabor elements $B$ by setting $[\I_m, B] \leftarrow 0$. So for each image $\I_m$, the selected $(B_{m, i}, i = 1, ..., n)$ are approximately orthogonal to each other.

The algorithm can learn from a single training image. If $M = 1$ and if we forbid perturbations in locations and orientations by setting $b_1 = b_2 = 0$, then the algorithm reduces to usual edge detection.

\indent{\em Transformation of responses}. To understand the transformation $h()$, let us consider a simplified discontinuous one: $h(r) = 1_{r>\xi}$, where $\xi$ is a threshold for edge detection. More specifically, $h(r) = 1$ if $r\geq \xi$, and $h(r) = 0$ otherwise. Then $\sum_{m=1}^{M} h(r_{m, i})$  simply counts the number of detected edge segments in the training images $\{\I_m, m = 1, ..., M\}$. That is, we select $B_i$ and perturb it to $\{B_{m, i}\}$, so that $\{B_{m, i}\}$ sketch as many edge segments as possible.

In this article we entertain the following designs of continuous transformations. The learned templates are not very sensitive to the choice of the transformation.

(1) {Sigmoid transformation.} The transformation is characterized by
a saturation level $\xi$ (e.g., $\xi = 6$),
\begin{eqnarray}
h(r) = {\rm sigmoid}(r) =  \xi \left[ \frac{2}{1 + e^{-2r/\xi}} - 1
\right], \label{eq:sigmoid}
\end{eqnarray}
which increases from 0 to $\xi$.

(2) {Whitening transformation.} Let $q(r)$ be the marginal distribution of $\langle \I, B_{x, y, s, \alpha}\rangle$ where $\I$ is a random sample from natural images. Let $F(t) = q(r>t)$, i.e., the probability that $r>t$ under $q(r)$. The non-linear
whitening transformation is
\begin{eqnarray}
h(r) ={\rm whiten}(r)= -\log F(r). \label{eq:whiten}
\end{eqnarray}

(3) {Thresholding transformation}. A crude but simple approximation to ${\rm whiten}(r)$ is
\begin{eqnarray}
h(r) = {\rm threshold}(r) = \min(r, T), \label{eq:threshold}
\end{eqnarray}
where $T$ is a threshold (e.g., $T= 16$).

\indent{\em Scoring template matching}. Let $\B = (B_i, i = 1, ..., n)$ be the template. For each training image $\I_m$, the template matching is scored by
\begin{eqnarray}
{\rm MATCH}(\I_m, \B) = \sum_{i=1}^{n} \left(\lambda_i [\I_m, B_{m, i}] - \log Z(\lambda_i)\right). \label{eq:match}
\end{eqnarray}
$\lambda_i$ can be calculated directly from $\sum_{m=1}^{M} [\I_m, B_{m, i}]$ in the selection step. $Z()$ is a non-linear function. This template matching score is actually a log-likelihood ratio for an exponential family model, and the weight vector  $\Lambda = (\lambda_i, i=1, ..., n)$ is estimated by maximum likelihood method. See the next section for details.

\indent{\em Active correlation.} We can also use a linear score for template matching:
\begin{eqnarray}
{\rm MATCH}(\I_m, \B) = \sum_{i=1}^{n} \theta_i [\I_m, B_{m, i}]. \label{eq:match1}
\end{eqnarray}
where $h(r) = {\rm whiten}(r)^{1/2}$ or $h(r) = {\rm threshold}(r)^{1/2}$, and $\Theta = (\theta_i, i=1, ..., n)$ is a unit vector, with $\|\Theta\|^2 = \sum_{i=1}^{n} \theta_i^{2} = 1$.

The elements are still selected by the shared sketch algorithm, with the new definition of $h()$. To estimate $\Theta$, we first calculate $\theta_i = \sum_{m=1}^{M} [\I_m, B_{m, i}]/M$, then we normalize $\Theta = (\theta_i, i=1, ..., n)$ to be a unit vector.

The template matching score (\ref{eq:match1}) can be considered the active correlation between the template $\B$ and the image $\I_m$, because $\B$ is deformed to $\B_m = (B_{m, i}, i = 1, ..., n)$ before the inner product is calculated. We may also consider (\ref{eq:match1}) as the inner product between $\I_m$ and the vector $V = \sum_{i=1}^{n} \theta_i B_i$. $V$ is an active vector because $B_i$ can be perturbed to $B_{m, i}$ when we correlate $V$ with $\I_m$.


\subsection{Inference: sum-max maps
\label{sect:recognition}}

After training the active basis model, specifically, after selecting
$\B = (B_i = B_{x_i, y_i, s, \alpha_i}, i = 1, ..., n)$, and
computing weighting vector $\Lambda = (\lambda_i, i = 1, ..., n)$ or $\Theta = (\theta_i, i = 1, ..., n)$, we can use the trained model to detect and sketch the template in a testing image.

Let $\I$ be a testing image, which is larger than the bounding box of the template $\B$. We assume that the bounding box of $\B$ is centered at origin $(x=0, y=0)$. We can scan the template over image $\I$, and at each position $(x, y)$, we fit the active basis model to the image within the bounding box centered at $(x, y)$, and calculate the template matching score according to Equation (\ref{eq:match}) or (\ref{eq:match1}).

The inference algorithm consists of two processes. The first process is a bottom-up scoring process, which calculates $\S1$, $\M1$, $\S2$, $\M2$ scores consecutively. The following are the questions that these scores seek to answer:

$\S1$ maps: Is there an edge segment at this location and orientation?

$\M1$ maps: Is there an edge segment at a {\em nearby} location and orientation?

$\S2$ map: Is there a certain composition of edge segments that form the template at this location?

$\M2$ score: Is there a certain composition within the whole image?

These maps are soft scores, not hard decisions. They are computed in a bottom-up process, $\S1 \rightarrow \M1 \rightarrow \S2 \rightarrow \M2$.

This is to be followed by a top-down retrieving process, which retrieves the central location of the template and then retrieves the locations and orientations of the basis elements of the deformed template. The following are the questions to be answered:

Back to $\M2$ score: If there is a template, where is it?

Back to $\S2$ map: What are the locations and orientations of the elements of the template before deformation?

Back to $\M1$ maps: What are the nearby locations and orientations that these elements are perturbed to?

Back to $\S1$ maps: What are the coefficients of these perturbed elements?

The top-down retrieving process follows the order of $\M2 \rightarrow \S2 \rightarrow \M1 \rightarrow \S1$. The process detects and deforms the template to interpret the observed image.

\vspace{.5cm} \noindent{\bf Pseudo-code for inference algorithm}

\begin{itemize}
\item[Input]: Template $\B = (B_i = B_{x_i, y_i, s, \alpha_i}, i = 1, ..., n)$, $\Lambda = (\lambda_i, i = 1, ..., n)$, and testing image $\I$.

\item[Output]: Location $(\hat{x}, \hat{y})$ of the detected template, and the deformed template $(B_{\hat{x}_i, \hat{y}_i, s, \hat{\alpha}_i}, i = 1, ..., n)$ that is matched to $\I$.
\end{itemize}

\begin{itemize}

\item[Up-1] For all $(x, y) \in D$, and for all $\alpha$, compute the
SUM1 maps:
\begin{eqnarray}
\S1(x, y, s, \alpha) = |\langle \I, B_{x, y, s, \alpha}\rangle|^2.
\label{eq:SUMmap2}
\end{eqnarray}

\item[Up-2] For all $(x, y) \in D$, and for all $\alpha$,  compute the MAX1
maps:
\begin{eqnarray}
\M1(x, y, s, \alpha) = \max_{\tiny{\begin{array}{c}d \in [-b_1,
b_1]\\ \delta \in [-b_2, b_2]\end{array}}} \S1(x + d \sin\alpha, y +
d \cos\alpha, s, \alpha + \delta). \label{eq:MAXmap2}
\end{eqnarray}

\item[Up-3] For all $(x, y) \in D$, and for all $\alpha$,  compute the SUM2
map:
\begin{eqnarray}
\S2(x, y) = \sum_{i=1}^{n} \left[ \lambda_i h(\M1(x+x_i, y+y_i, s, \alpha_i)) - \log Z(\lambda_i)\right].  \label{eq:SUMmap2}
\end{eqnarray}

\item[Up-4]  Compute the MAX2 score:
$\M2 = \max_{x, y} \S2(x, y)$.

\item[Down-4] Retrieve $(\hat{x}, \hat{y})$ that achieves the maximum in the computation of Up-4.

\item[Down-3] Retrieve $(\hat{x}+x_i, \hat{y}+y_i, \alpha_i)$ in the computation of $\M1(x+x_i, y+y_i, s, \alpha_i)$ for $i = 1, ..., n$ in Up-3.

\item[Down-2] Retrieve $\hat{x}_i, \hat{y}_i, \hat{\alpha}_i$ for $i = 1, ..., n$, such that
\begin{eqnarray}
 \M1(\hat{x}+x_i, \hat{y}+y_i, s, \alpha_i) =  \S1(\hat{x}_i, \hat{y}_i, s,  \hat{\alpha}_i)
\end{eqnarray}
in the local maximization operation (\ref{eq:MAXmap2}) of Up-2.

\item[Down-1] Retrieve the coefficients in the computation of $\S1(\hat{x}_i, \hat{y}_i, s, \hat{\alpha}_i)$ for $i = 1, ..., n$ in Up-1
\end{itemize}


Then the Gabor wavelet elements $(B_{\hat{x}_i, \hat{y}_i, s, \hat{\alpha}_i}, i=1, ..., n)$ form the deformed template fitted to image $\I$.

The $\S2$ map in Up-3 scores template matching. The computation of $\S2$ can be considered a shape filter for template matching. Like Gabor filters, it is also a local weighted summation operator. See Figure (\ref{fig:SumMaxMaps}) for an illustration. The shape filter in Up-3 has fixed $(x_i, y_i, \alpha_i, i = 1, ..., n)$. But it is computed on the $\M1$ maps instead of $\S1$ maps, so it is invariant to shape deformation.

For an input image, we can apply the above algorithm at multiple
resolutions of the input image. Then we can choose the resolution
that achieves the maximum $\M2$ score as the optimal resolution.


\subsection{Shared sketch algorithm based on sum-max maps
\label{sect:SS}}

The shared sketch algorithm in Subsection (\ref{sect:la}) can be
expressed more precisely in terms of the sum maps and max maps.

\vspace{.5cm} \noindent{\bf Pseudo-code for shared sketch algorithm}

\begin{itemize}
\item[Input]: Training images $\{\I_m, m = 1, ..., M\}$.

\item[Output]: Template $\B = (B_i = B_{x_i, y_i, s, \alpha_i}, i = 1, ..., n)$, weighting vector $\Lambda = (\lambda_i, i = 1, ..., n)$, and deformed template $\B_m = (B_{m, i} = B_{x_{m, i}, y_{m, i}, s, \alpha_{m, i}}, i=1, ..., n)$ that is matched to $\I_m$ for $m = 1, ..., M$.
\end{itemize}

\begin{enumerate}

\item Convolution: For each $m = 1, ..., M$, for all $(x, y) \in D$, and for all $\alpha$, compute the SUM1 maps:
\begin{eqnarray}
\S1_{m}(x, y, s, \alpha) = |\langle \I_m, B_{x, y, s,
\alpha}\rangle|^2.  \label{eq:SUMmap1}
\end{eqnarray}

\item
Local maximization: For each $m = 1, ..., M$, for all $(x, y) \in D$, and for all $\alpha$, compute the MAX1 maps:
\begin{eqnarray}
\M1_{m}(x, y, s, \alpha) = \max_{\tiny{\begin{array}{c}d \in [-b_1,
b_1]\\ \delta \in [-b_2, b_2]\end{array}}} \S1_{m}(x + d \sin\alpha,
y + d \cos\alpha, s, \alpha + \delta). \label{eq:MAXmap1}
\end{eqnarray}

For each $m = 1, ..., M$, set $\S2_m \leftarrow 0$. Set $i
\leftarrow 1$.

\item
Selection: Find $(x_i, y_i, \alpha_i)$ by maximizing
$\sum_{m=1}^{M} h(\M1_{m}(x, y, s, \alpha))$ over all $(x, y, \alpha)$.

For each  $m = 1, ..., M$, retrieve $({x}_{m, i}, {y}_{m, i}, {\alpha}_{m, i})$ so that
\begin{eqnarray}
r_{m, i} = \M1_{m}(x_i, y_i, s, \alpha_i) = \S1_{m}({x}_{m, i}, {y}_{m, i}, s, {\alpha}_{m, i}) \label{eq:MAXback1}
\end{eqnarray}
in the local maximization computation in (\ref{eq:MAXmap1}). That
is, we perturb $B_i= B_{x_i, y_i, s, \alpha_i}$ to $B_{m, i} = B_{{x}_{m, i}, {y}_{m, i}, s, {\alpha}_{m, i}}$ to fit $\I_m$.

Compute ${\lambda}_i$ from $\sum_{m=1}^{M} h(\M1_{m}(x, y, s, \alpha))$. For each $m = 1, ..., M$, compute
$\S2_m \leftarrow \S2_m + {\lambda}_i h(r_{m, i}) - \log
Z({\lambda}_i)$.


\item
Non-maximum suppression: If $r_{m, i}>0$, then for all those $(x, y, \alpha)$ such that ${\rm corr}(B_{{x}_{m, i}, {y}_{m, i}, s, {\alpha}_{m, i}},$ $B_{x, y, s, \alpha})$ $> \epsilon$, set $ \S1_{m}(x, y, s, \alpha) \leftarrow 0. $


Re-compute the MAX1 maps according to (\ref{eq:MAXmap1}).

\item
Stop if $i = n$. Otherwise let $i \leftarrow i + 1$, and go back to Step 3.
\end{enumerate}

The above algorithm can be easily mapped to computer code. The
following are some remarks on implementing it.

(1) In updating the SUM1 maps and the MAX1 maps in Step 4, we only need to update the parts of the maps that are affected.

(2) The correlation  ${\rm corr}(B_{{x}_{m, i}, {y}_{m, i}, s,
{\alpha}_{m, i}}, B_{x, y, s, \alpha})$ in Step 4 only depends on
$(x_{m, i}-x, y_{m, i}-y, \alpha_{m, i}-\alpha)$. We can store a
correlation function $C(\Delta x, \Delta y, \Delta \alpha) = {\rm
corr}(B_{x+\Delta x, y+ \Delta y, s, \alpha + \Delta \alpha}, B_{x,
y, s, \alpha})$ before we run the algorithm.

(3) After selecting $(x_i, y_i, \alpha_i)$, we need to go back to
retrieve $({x}_{m, i}, {y}_{m, i}, \alpha_{m, i})$. We can record
this information for every $(x, y, \alpha)$ when performing the
local maximization in (\ref{eq:MAXmap1}), and then retrieve the
information for $(x_i, y_i, \alpha_i)$. If we choose not to store
this information beforehand, we can re-do the local maximization in
(\ref{eq:MAXmap1}) for $(x_i, y_i, \alpha_i)$ to retrieve $({x}_{m,
i}, {y}_{m, i}, \alpha_{m, i})$.


The $\S2_m$ score evaluates the matching of $\I_m$ to the learned
template $\B$ according to Equation (\ref{eq:match}). The total
score $\sum_{m=1}^{M} \S2_m$ measures the overall alignment of all
the training images. This alignment score is very useful for
unsupervised learning, where the objects in the training images are of unknown locations, scales, and categories. The alignment score $\sum_{m=1}^{M} \S2_m$ is the criterion that determines these hidden variables.

We would like to point out a subtle difference between the computation of $\S2_m$ in the learning algorithm and the computation of $\S2$ map in the inference algorithm. In the learning algorithm, there is a non-maximum suppression step, where $B_{m, i}$ suppresses nearby overlapping elements. This is necessary for selecting the basis elements. In the inference algorithm, we omit this step for efficiency. This is because the elements selected by the learning algorithm are already well spaced due to the non-maximum suppression in learning, so there is no much need for non-maximum suppression in inference. In this paper, we use the inference algorithm for detection. For classification, we use the learning algorithm where non-maximum suppression is applied. See Section (\ref{sect:DC}) for details.

\section{Theoretical Underpinning}

This section presents theoretical underpinnings of the model and the algorithms presented in the previous section. Readers who are more interested in applications and experiments can jump to the next section.

\subsection{Probability distribution on image intensities \label{sect:model}}

With multiple training images $\{\I_m, m = 1, ..., M\}$ represented
by (\ref{eq:ab}) to (\ref{eq:ab2}), we can pool the probability
distribution of $\{(c_{m, i}, i = 1, ..., n)\}$ as well as the
distribution of $\{U_m\}$ over $m = 1, ..., M$. With these two
distributions, we can obtain the distribution of $\I_m$, or more
specifically, the distribution of $\I_m$ given $\B_m$, $p(\I_m \mid
\B_m)$.  The probability density $p(\I_m \mid B_m)$ can be used for
maximum likelihood learning of $\B$ and $\{\B_m\}$ from training
images. It can also be used for finding $\B_m \approx \B$ and
scoring the template matching in recognition after $\B$ is learned
from training images.

We first simplify the notation using matrices and vectors.  $\I_m$
can be treated as a $|D| \times 1$ column vector, where $|D|$ is the number of pixels. $\B = (B_{i, 0}, B_{i, 1}, i = 1, ..., n)$ can be treated as a $|D| \times 2n$ matrix, where each $B_{i, \eta}$ ($\eta= 0, 1$) is a $|D| \times 1$ vector. Each $\B_{m}$ can be treated as a $|D| \times 2n$ matrix in the same way. We can write $C = (c_{m, 0}, c_{m, 1}, i = 1, ..., n)'$ as a $2n \times 1$ vector. Thus in matrix notation, Equation (\ref{eq:ab}) becomes $\I_m = \B_m C_m + U_m$.

{\em Linear decomposition.} We assume that $\B_m C_m$ is the
projection of $\I_m$ onto the subspace spanned by the column vectors
of $\B_m$, so $C_m = (\B_m'\B_m)^{-1}\B_m' \I_m$. If $\B_m$ is
orthogonal, then $C_m = \B_m' \I_m$. $U_m$ resides in the $|D|-2n$
dimensions that are orthogonal to the columns of $\B_m$. There is no
loss of generality in such an assumption, because if $U_m$ is not
orthogonal to $\B_m$, we can always project $U_m$ onto $\B_m$, and
let $\B_m C_m$ absorb this projection. We can write $U_m = \bB_m
\bC_m$, where $\bB_m$ is a $|D| \times (|D|-2n)$ matrix whose
columns are orthogonal to those of $\B_m$, and $\bC_m$ is a
$(|D|-2n) \times 1$ vector. Thus, $\I_m = \B_m C_m + \bB_m \bC_m$.
There is a one-to-one linear mapping between $\I_m$ and $(C_m,
\bC_m)$. $\bB_m$ and $\bC_m$ can be made implicit in statistical
modeling.

{\em Shape and texture.} Now we are ready to specify the probability
density $p(\I_m \mid \B_m)$. For the linear representation $\I_m =
\B_m C_m + \bB_m \bC_m$,
\begin{eqnarray}
p(\I_m \mid \B_m) = p(C_m, \bC_m)|J_m| = p(C_m) p(\bC_m \mid C_m)
|J_m|, \label{eq:prob}
\end{eqnarray}
where $|J_m|$ is the absolute value of the determinant of the
Jacobian matrix of the linear transformation from $\I_m$ to $(C_m,
\bC_m)$. $p(C_m)$ is the distribution of the coefficients for coding
the foreground shape, and $p(\bC_m \mid C_m)$ is the distribution of
the residual background given the foreground coefficients.

Let $q(\I_m)$ be a reference distribution. We can write $q(\I_m) =
p(C_m) q(\bC_m \mid C_m) |J_m|$ with the same Jacobian $J_m$. We
want to construct $p(\I_m \mid \B_m)$ by modifying $q(\I_m)$.
Specifically, we assume that $p(\bC_m \mid C_m) = q(\bC_m \mid
C_m)$, i.e., the conditional distribution of the residual background
in $p(\I_m)$ is assumed to be the same as that in $q(\I_m)$. Then
\begin{eqnarray}
  p(\I_m \mid \B_m) = q(\I_m) \frac{p(C_m)}{q(C_m)} =
  q(\I_m) \frac{p(c_{m, 1}, ..., c_{m, n})}{q(c_{m, 1}, ..., c_{m,
  n})}, \label{eq:rewrite0}
\end{eqnarray}
where we substitute $p(C_m)$ for $q(C_m)$ to construct a density
$p(\I_m)$ from $q(\I_m)$.

We assume $q(\I_m)$ to be stationary. The following are some choices
of $q(\I_m)$.

(1) White noise distribution. This is the distribution that is often
assumed in linear additive model, and is implicitly assumed in the
least squares criterion for model fitting. Under this reference distribution, $q(c_{m, 1}, ..., c_{m, n})$ is multivariate Gaussian.

(2) The distribution of natural image patches. This is the
distribution that we shall use in this paper. In particular, we make
use of the marginal distribution of filter responses $\langle \I_m,
B_{x, y, s, \alpha} \rangle$ in natural images. It is a heavy tail
distribution that allows occasional strong edges. We do not need to
specify $q(\I_m)$ beyond this marginal distribution.

(3) The Markov random field distribution that matches the marginal
distributions of filter responses $\langle \I_m, B_{x, y, s, \alpha}
\rangle$ in natural images. Such a Markov random field model has
been developed by Zhu and Mumford \cite{Zhu_prior}. It is a more
explicit form of (2).

(4) The Markov random field distribution that matches the marginal
distributions of filter responses of the observed image $\I_m$. Such
a model has been developed by Zhu, Wu, and Mumford \cite{Zhu_FRAME}.
The marginal distributions are pooled from the observed image $\I_m$
over $(x, y) \in D$.

(5) Uniform distribution over the set of images that share certain
marginal statistics pooled over image domain. Such a distribution
can be related to the Markov random field model in (4) as shown in
Wu, Zhu, and Liu \cite{Wu_equivalence}.

The model (\ref{eq:rewrite0}) combines both texture and shape.
$q(\I_m)$ models the background texture, and $\B_m$ and $p(C_m)$
model the foreground shape. The foreground shape pops out from the
background texture, as modeled by the probability ratio
$p(C_m)/q(C_m)$.


{\em Log-likelihood and Kullback-Leiber divergence.} To learn $\B$
and $\{\B_m \approx \B, m = 1, ..., M\}$, we can maximize the
average log-likelihood
\begin{eqnarray}
   \frac{1}{M} \sum_{i=1}^{M} \log \frac{p(\I_m\mid \B_m)}{q(\I_m)} &=&
   \frac{1}{M} \sum_{i=1}^{M} \log \frac{p(c_{m, 1}, ..., c_{m, n})}
   {q(c_{m, 1}, ..., c_{m, n})}. \label{eq:lik}
\end{eqnarray}
The average log-likelihood converges to $\KL(p(C_m)||q(C_m))$ as $M
\rightarrow \infty$, assuming that $p(C_m)$ can be consistently
estimated from the training images. Here $\KL(p||q)$ denotes the
Kullback-Leibler divergence from $p$ to $q$. In order to maximize
the log-likelihood, we want to choose $\B$ and deform it to $\{\B_m
\approx \B\}$ to maximize $\KL(p(C_m)||q(C_m))$, so that we achieve
the maximum contrast between the foreground shape and the background
texture. $\KL(p(C_m)||q(C_m))$ also measures the coding gain
achieved by coding $C_m$ by $p(C_m)$ instead of $q(C_m)$, while
continuing to code the residual background by $q(\bC_m|C_m)$.

It is impossible to select $\B$ and $\{\B_m\}$ all at once. In the
next subsection, we present an algorithm that sequentially pursues
$B_i$ and perturbs it to $\{B_{m, i}\}$.




\subsection{Coupling matching pursuit with projection pursuit}

In this subsection, we describe a shared matching pursuit process for selecting the basis elements $\B = (B_i, i = 1, ..., n)$. The process couples matching pursuit \cite{Mallat_Zhang} with projection pursuit \cite{Friedman}. The shared sketch learning algorithm is an approximation to it.

The matching pursuit is a process that sequentially adds elements
$B_{m, i}, i = 1, ..., n$ to improve the encoding of image $\I_m$.
It has the following form:

\begin{enumerate}
\item For $m = 1, ..., M$, set $U_m \leftarrow \I_m$. Set $i
\leftarrow 1$.

\item For $m = 1, ..., M$, choose $B_{m, i}$. Let $c_{m, i} =
\langle U_m, B_{m, i} \rangle$.

\item For $m = 1, ..., M$, update $U_m \leftarrow U_m - c_{m, i}
B_{m, i}$. Represent $\I_m = c_{m, 1} B_{m, 1} + ...+c_{m, i} B_{m,
i} + U_m$.

\item If $i = n$, stop. Otherwise, set $i \leftarrow i+1$, go back
to Step 2.
\end{enumerate}

We need to add the following three steps to the above matching
pursuit process.

(1) {\em The selection of $B_{m, i}$ given $B_i$}. The original
matching pursuit algorithm selects $B_{m, i} = \arg \max_{B}
|\langle U_m, B \rangle|^2$ in Step 2, where the maximization is
over all $B \in {\rm Dictionary}$, so that $B_{m, i}$ achieves the
best fit to the unexplained residual image $U_m$. In shared matching
pursuit process, however, the $B_{m, i}$ are constrained to be
perturbed versions of a commonly shared $B_i$. Therefore, for each
putative $B_i$, we need to select $B_{m, i} = \arg \max_{B \approx
B_i} |\langle U_m, B \rangle|^2$.

(2) {\em The updating of $p(\I_m)$}. After computing $c_{m, i} =
\langle U_m, B_{m, i} \rangle$ in each iteration $i$, we can pool a
distribution $p_i(c)$ over $\{c_{m, i}, m = 1, ..., M\}$. We can use
such pooled densities $p_1(c)$, ..., $p_n(c)$ to construct the
density $p(\I_m)$.

Specifically, we update $p(\I_m)$ sequentially using projection
pursuit. Let $p_0(\I_m) = q(\I_m)$, i.e., the distribution of
background texture. After selecting $\{B_{m, i}, m = 1, ..., M\}$,
we need to update $p_{i-1}(\I_m)$ to $p_i(\I_m)$. We can apply the
density substitution scheme of projection pursuit, and let
$p_i(\I_m) = p_{i-1}(\I_m) p_i(c_{m, i})/q_{i-1}(c_{m, i})$,
where $q_{i-1}(c)$ is the density
of $c_{m, i} = \langle U_m, B_{m, i} \rangle$ under the current
model $p_{i-1}(\I_m)$. This density substitution scheme is very
similar to the model construction scheme of Equation
(\ref{eq:rewrite0}), except that we use $p_{i-1}(\I_m)$ as the
current background, and we only replace the density of $c_{m, i} =
\langle U_m, B_{m, i} \rangle$ under $p_{i-1}(\I_m)$.
 $c_{m, i} = \langle U_m, B_{m, i}\rangle$ can also be written as
 $c_{m, i} = \langle \I_m, \tilde{B}_{m, i} \rangle$,
where $\tilde{B}_{m, i}$ can be constructed from $B_{m, 1}, ...,
B_{m, i-1}$ and $B_{m, i}$. So $p_{i}(\I_m)$ is a legitimate density
function.

(3) {\em The selection of $B_i$}. We select $B_i$ sequentially by
the maximum likelihood principle. The increase in the average
log-likelihood is $\sum_{m=1}^{M} \log [p_i(\I_m)/p_{i-1}(\I_m)]/M =
\log [p_i(c_{m, i})/q_{i-1}(c_{m, i})]/M \rightarrow
\KL(p_i(c)||q_{i-1}(c))$. So we want to select $B_i$ that achieves
the maximum $\KL(p_i(c)||q_{i-1}(c))$. That is,
$\KL(p_i(c)||q_{i-1}(c))$ is the pursuit index that drives the
selection of $B_i$. Intuitively, this means that we want to select
$B_i$ so that the distribution of the responses of the perturbed
versions $\{B_{m, i} \approx B_i\}$ is most different from what is
predicted by the current model $p_{i-1}(\I_m)$.

With (1), (2), and (3) incorporated into the matching pursuit
process, we will eventually reach the model $p(\I_m) = q(\I_m)
\prod_{i=1}^{n} p_i(c_{m, i})/q_{i-1}(c_{m, i})$. This is an
approximation to the model (\ref{eq:rewrite0}) in the previous
subsection. See Figure (\ref{fig:SharedSketch}) for an illustration of the shared matching pursuit process.

The computational burden in the shared matching pursuit process lies
in the computation of $q_{i-1}(c)$, which requires Monte Carlo sampling from $p_{i-1}(\I_m)$. If we have negative training images from $q(\I_m)$, we can re-weight these negative examples after each iteration, and use these re-weighted examples as samples from $p_{i-1}(\I_m)$.



\subsection{Shared sketch as an approximation}

We can simplify the shared matching pursuit process into a shared
sketch process with the following two approximations.

(1) {\em Non-maximum suppression}. After selecting $B_{m, i}$ and computing $c_{m, i} = \langle U_m, B_{m, i} \rangle$, we need to update
$U_m \leftarrow U_m - c_{m, i} B_{m, i}$, i.e.,
$B_{m, i}$ explains away part of $U_m$ or $\I_m$. This can be considered a soft inhibition. If an element $B$ has a high correlation with $B_{m, i}$, in other words, if $B$ heavily overlaps with $B_{m, i}$ in both spatial domain and frequency domain, then such a redundant $B$ can add little to further explaining $\I_m$, in that after the updating $U_m \leftarrow U_m - c_{m, i} B_{m, i}$, $|\langle U_m, B\rangle|^2$ can be very small. Therefore, we may simply enforce that, for each $\I_m$,  the selected elements of $\{B_{m, i}, i = 1, ..., n\}$ do not overlap with each other, or the selected $\{B_{m, i}, i = 1, ..., n\}$ are orthogonal to each other. Then, after $B_{m, i}$ is selected, we let $B_{m, i}$ suppress any $B$ that overlaps with $B_{m, i}$. For such non-overlapping $(B_{m, i}, i = 1, ..., n)$, $c_{m, i} = \langle U_m, B_{m, i} \rangle = \langle \I_m, B_{m, i} \rangle$. In practice, we allow small correlations between the elements $(B_{m, i}, i = 1, ..., n)$.

(2) {\em Background density}. The current density $p_{i-1}(\I_m)$
results from sequentially updating the densities of $c_{m, 1}$, ...,
$c_{m, i-1}$, starting from $q(\I_m)$. If $B_{m, i}$ has no overlap
with $B_{m, 1}$, ..., $B_{m, i-1}$, then the distribution of $c_{m,
i}$ under $p_{i-1}(\I_m)$, i.e., $q_{i-1}(c)$, can be approximated
by $q(c)$, which is the marginal distribution of $c_{m, i}$ under
$q(\I_m)$. Because $q(\I_m)$ is stationary, $q(c)$ is the same for
all $c_{m, i}, i = 1, ..., n$. Therefore, the pursuit index is
$\KL(p_{i}(c)||q(c))$, where again, $p_i(c)$ is the density pooled
over $\{c_{m, i}, m = 1, ..., M\}$.

If we stop the process after $n$ iterations, then the resulting
model is
\begin{eqnarray}
p(\I_m \mid \B_m) = q(\I_m) \prod_{i=1}^{n} \frac{p_i(c_{m,
i})}{q(c_{m, i})}. \label{eq:ind}
\end{eqnarray}
The pursued model (\ref{eq:ind}) is an approximation to the model
(\ref{eq:rewrite0}). $q(c)$ can be pooled from natural images before
we start the shared sketch process. We do not need negative examples
beyond $q(c)$.


\subsection{Parametrization by exponential family model}

{\em Parametric model.} We can further simplify the Kullback-Leibler divergence by assuming the following exponential family model:
\begin{eqnarray}
    p(c; \lambda) = \frac{1}{Z(\lambda)} \exp\{\lambda h(r)\} q(c),
    \label{eq:exponential}
\end{eqnarray}
where $\lambda>0$ is the parameter, $r = |c|^2$, and $Z(\lambda) =
\int \exp\{\lambda h(r)\}q(c)dc = \E_{q}[\exp\{\lambda h(r)\}]$ is
the normalizing constant. $h(r)$ is an increasing function, so
$p(c; \lambda)$ puts more probability than $q(c)$ on those $c$ with large $r$, so model (\ref{eq:exponential}) is more encouraging towards large $r$. The above model can be justified by the maximum entropy principle \cite{DellaPietra, Zhu_FRAME}.

Let $p(r; \lambda)$ and $q(r)$ be the densities of $r = |c|^2$ under
$p(c; \lambda)$ and $q(c)$ respectively, then $p(c; \lambda)/q(c) =
p(r; \lambda)/q(r) = \exp\{\lambda h(r)\}/Z(\lambda)$. More
specifically, let $c = (c_0, c_1)$, and let $\varphi =
\arctan(c_1/c_0)$ be the local phase, then the conditional
distributions of the local phase $\varphi$ given the local energy
$r$ are the same under both $p(c; \lambda)$ and $q(c)$.

{\em Estimating $p_i$.} We estimate $q(c)$ by pooling a histogram
from natural images. We estimate $p_i(c)$ from $\{c_{m, i} =
\langle\I_m, B_{m, i}\rangle, m = 1, ..., M\}$ by fitting the
density $p(c; \lambda_i)$ to $\{c_{m, i}\}$. Specifically, let us
define the mean parameter $\mu(\lambda) = \E_{\lambda}[h(r)] = \int
h(r) \exp\{\lambda h(r)\} q(r)dr/Z(\lambda) $,  we estimate
$\lambda_i$ by matching $\mu(\lambda_i) = \sum_{m=1}^{M} h(r_{m,
i})/M$, so that $\hat{\lambda}_i = \mu^{-1}(\sum_{m=1}^{M} h(r_{m,
i})/M)$. $\hat{\lambda}_i$ is  the maximum likelihood estimate that
maximizes $\sum_{m=1}^{M} \log [p(c_{m, i}; \lambda_i)/q(c_{m, i})]$
over $\lambda_i$ \cite{DellaPietra}. Thus, we estimate $p_i(c)$ by
$p(c; \hat{\lambda}_i)$.

{\em Selecting $B_i$}. The average log-likelihood $\sum_{m=1}^{M}
\log [p(c_{m, i}; \hat{\lambda}_i)/q(c_{m, i})]/M = \KL(p(c;
\hat{\lambda}_i)||q(c))$.  It is an increasing function of
$\sum_{m=1}^{M} h(r_{m, i})/M$. Therefore, we choose $B_i$ and
perturb it to $\{B_{m, i}\}$ by maximizing the pursuit index
$\sum_{m=1}^{M} h(r_{m, i})$.

{\em Perturbing $B_i$ to $B_{m, i}$. } $p(c; \lambda_i)/q(c)$ is a
monotone increasing function of $r = |c|^2$. This justifies that,
given $B_i$, we should perturb $B_i$ to $B_{m, i}$ to maximize
$|\langle \I_m, B_{m, i} \rangle|^2$, subject to the approximate
non-overlapping constraint. Such $B_{m, i}$ is the maximum
likelihood estimate given $B_i$.

Thus, the estimation of $\lambda_i$, the perturbation of $B_i$ to
$B_{m, i}$, and the selection of $B_i$ all follow the maximum
likelihood principle.

{\em Template matching score}. To score the template matching, we
can compute the log-likelihood ratio
\begin{eqnarray}
\log \frac{p(\I_m \mid \B_m)}{q(\I_m)} &=& \log \prod_{i=1}^{n}
\frac{p_{i}(c_{m, i}; \hat{\lambda}_i)}{q(c_{m, i})} \nonumber \\
&=& \sum_{i=1}^{n} \left[ \hat{\lambda}_i h(r_{m, i}) - \log
Z(\hat{\lambda}_i)\right ], \label{eq:score}
\end{eqnarray}
where $p(\I_m \mid \B_m)$ is the pursued model (\ref{eq:ind}).

Both $\log Z(\lambda)$ and $\mu(\lambda)$ are one-dimensional
functions. We can store their values over a grid of points, and use
nearest neighbor linear interpolation for points in between.

The reader is referred to Wu et al. \cite{Wu} for an information-theoretical perspective of this model, as well as its connection with Markov random fields \cite{Zhu_FRAME} and adaboost \cite{adaboost}. The reader is also referred to Tu \cite{Tu} on a generative model constructed from adaboost.

\subsection{Transformation of filter responses}

The following are explanations why we use the sigmoid and whitening transformations for $h(r)$.

{\em Sigmoid transformation.} The saturation in the sigmoid transformation can be justified by mixture distributions.

Let $p_{\rm on}(r)$ be the density of $r = |\langle \I, B_{x, y, s, \alpha}\rangle|^2$ when the Gabor wavelet $B_{x, y, s, \alpha}$ is on an edge. Let $p_{\rm off}(r)$ be the density of $r$ when the
Gabor wavelet is off the edge. We assume that  $p_{\rm on}(r)$ has a
much longer tail than  $p_{\rm off}(r)$. Let $q(r)$ and $p_i(r)$ be
the densities of $r = |c|^2$ under $q(c)$ and $p_i(c)$ respectively.
It is reasonable to assume that $q(r) = (1- \rho_0) p_{\rm off}(r) +
\rho_0 p_{\rm on}(r)$, and $p_i(r) = (1- \rho_i) p_{\rm off}(r) +
\rho_i p_{\rm on}(r)$, that is, both $q(r)$ and $p_i(r)$ are
mixtures of on-distribution and off-distribution, with $\rho_i >
\rho_0 >0$. As $r \rightarrow \infty$, $\log[p_i(r)/q(r)]
\rightarrow \log(\rho_i/\rho_0) > 0$, i.e., a positive constant. So
we may assume the following log-linear model $\log [p_i(c)/q(c)] =
\log [p_i(r)/q(r)]  = \lambda_i h(r) + {\rm constant}$, where
$\lambda_i > 0$, and $h(r)$ reaches a saturation level as $r
\rightarrow \infty$. This justifies the choice of the sigmoid transformation.

{\em Whitening transformation.} The whitening transformation makes $q(\I_m)$ closer to the white noise distribution, which is a simpler null hypothesis. It also leads to explicit expressions of $\mu(\lambda)$ and $\log Z(\lambda)$.

Let $F(t) = q(r>t)$, i.e., the
probability that $r>t$ under $q(r)$ or $q(\I_m)$.
The reason we call $h(r) = -\log F(r)$ the whitening transformation is that $\Pr(h(r)>t) = \Pr(-\log F(r)>t) = \Pr(F(r)<e^{-t}) = e^{-t}$, i.e., $h(r)$ follows Exponential distribution with unit expectation. This is the distribution of $r$ if $q(\I_m)$ is Gaussian white noise. This is because the local energy $r$ is the sum of the squares of the Gabor sine response and Gabor cosine response, and both of them follow independent Normal distributions if $q(\I_m)$ is Gaussian white noise. So their sum of squares follows a $\chi^2_2$
distribution, which is the Exponential distribution. The
distribution has expectation 1 because we normalize the image to
have unit $\sigma^2(s)$, see Equation (\ref{eq:normalize}).

The whitening transformation changes a long tail distribution $q(r)$
to a short tail Exponential distribution. With the whitening
transformation, under $p(c; \lambda)$ of (\ref{eq:exponential}),
 $h(r) \sim {\rm Exp}(1-\lambda)$,
which is an Exponential distribution with $\mu(\lambda) =
1/(1-\lambda_i)$ and $Z(\lambda) = 1/(1-\lambda_i)$. $\lambda_i$ can
be estimated by $\hat{\lambda}_i = 1 - M/\sum_{m=1}^{M} h(r_{m,
i})$.

\subsection{Active mean vector and active correlation  \label{sect:AC}}

We can replace the log-likelihood score $\log [p(\I_m \mid
\B_m)/q(\I_m)]$ in Equation (\ref{eq:score}) by the correlation
between $\I_m$ and the vector $V_m = \sum_{i=1}^{n} \theta_i B_{m,
i}$, which is defined as
\begin{eqnarray}
    \langle \I_m | V_m \rangle =
    \sum_{i = 1}^{n} \theta_i {\rm whiten}(|\langle \I_m, B_{m,
    i}\rangle|^2)^{1/2}.
    \label{eq:CORR0}
\end{eqnarray}
We assume that $\I_m$ is normalized. The reason we use whitening transformation defined by Equation
(\ref{eq:whiten}) is that after such a transformation, the
distribution of the natural images is closer to  white noise.
Geometrically, the white noise distribution is close to the uniform
distribution over a high dimensional sphere. Image patches (after
normalization and whitening transformation) from the same object
category form a cluster on this sphere. Such a simple picture makes
the concept of correlation geometrically meaningful. The correlation
score (\ref{eq:CORR0}) can be considered the length that $\I_m$
projects on $V_m$. We also filter out the local phase information in
(\ref{eq:CORR0}), because phase is irrelevant for shape. We  call
(\ref{eq:CORR0}) the active correlation between $\I_m$ and the
vector $V = \sum_{i=1}^{n} \theta_i B_i$, because we perturb $V$ to
$V_m$ in order to best correlate with $\I_m$.

For the training images $\{\I_m, m = 1, ..., M\}$, we want to find
the vector $V = \sum_{i=1}^{n} \theta_i B_i$ that best correlates
with $\{\I_m, m = 1, ..., M\}$, by maximizing the sum of the active
correlation scores:
\begin{eqnarray}
\sum_{i=1}^{m} \langle \I_m \mid V_m \approx V\rangle =
\sum_{i=1}^{n} \left[\theta_i   \sum_{m=1}^{M} {\rm whiten}(|\langle
\I_m, B_{m, i}\rangle|^2)^{1/2} \right]. \label{eq:CORR}
\end{eqnarray}
The algorithm for learning $\B = (B_i, i = 1, ..., n)$ and  $\Theta
= (\theta_i, i = 1, ..., n)$ is essentially the same as the shared sketch algorithm in Subsection (\ref{sect:SS}). The resulting $V = \sum_{i=1}^{n} \theta_i B_i$ can be consider the
mean shape of $\{\I_m, m = 1, ..., M\}$. We call it the active mean
vector. Geometrically, $V$ points to the center of the cluster
formed by $\{\I_m, m = 1, ..., M\}$.







\section{Supervised Learning, Detection, and Classification
\label{sect:DC}}

This section applies the learning and inference algorithms to supervised learning, detection, and classification.

\subsection{Learning with given bounding boxes}

In supervised learning, we assume that the training images are defined on the same image lattice which is the bounding box of the objects in these images.

In the experiments in this article, we hand pick the number of basis elements, $n$. In principle, it can be automatically determined by comparing $\sum_{m=1}^{M} h(r_{m, i})/M$ with the average of $h(\M1_m(x, y, s, \alpha))$ in natural images or in the observed image $\I_m$. If the former is no much greater than the latter, we should stop the algorithm. We also hand pick the resize factor of the training images. Of course, in each experiment, the same resize factor is applied to all the training images.

{\em Parameter values.} The following are the parameter values that we used in all the experiments in this paper (unless otherwise stated). Size of Gabor wavelets = $17 \times 17$.
$(x, y)$ is sub-sampled every 2 pixels. The orientation $\alpha$
takes $A = 15$ equally spaced angles in $[0, \pi]$. The
orthogonality tolerance is $\epsilon = .1$. The threshold is $T =
16$ in the threshold transformation (\ref{eq:threshold}). The
saturation level $\xi = 6$ in the sigmoid transformation
(\ref{eq:sigmoid}). The shift along the normal direction $d_{m, i}
\in [-b_1, b_1] = [-6, 6]$ pixels. The shift of orientation
$\delta_{m, i} \in [-b_2, b_2] = \{-1, 0, 1\}\times \pi/15$.

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=ex1_1.eps,height=2in}
\end{tabular}
\end{center}
\caption{Experiment 1.1. The $37$ training images are 82 $\times
164$. The first block displays the learned active basis consisting
of 60 elements. Each element is symbolized by a bar. The rest of the
blocks display the observed images and the corresponding deformed
active bases. The images are displayed in the descending order of
the log-likelihood ratio, which scores the template matching. }
 \label{fig:ex1_1}
\end{figure}

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=ex1_2.eps,height=2in}
\end{tabular}
\end{center}
\caption{Experiment 1.2. The $27$ images are 180 $\times$ 180.
Number of elements is 60.  } \label{fig:ex1_2}
\end{figure}

{\em Experiment 1}. In this experiment, we take $h(r) = {\rm threshold}(r)$, as defined by Equation (\ref{eq:threshold}), so that there is  no need to pool $q(r)$. Other designs of $h(r)$ produce similar results.



In Experiment 1.1, we  apply the shared sketch algorithm to a
training set of $M = 37$ car images. The car images are $82 \times
164$. Figure (\ref{fig:ex1_2}) displays the results, where $n = 60$.
The first block displays the learned active basis $\B = \{B_i, i= 1,
..., n = 60\}$, where each $B_i$ is represented symbolically by a
bar at the same location and with the same length and orientation as
$B_i$. The intensity of the bar that symbolizes $B_i$ is the average
$\sum_{m=1}^{M} h(\M1_{m}(x_i, y_i, s, \alpha_i))/M$. For the
remaining $M$ pairs of plots, the left plot shows $\I_m$, and the
right plot shows $\B_m = (B_{m, i}, i = 1, ..., n)$. The intensity
of the bar that symbolizes $B_{m, i}$ is the squared root of
$h(|\langle \I_m, B_{m, i}\rangle|^2)$. These $M$ examples are
arranged in descending order by the $\S2_m$ scores output by the
algorithm. We can see that all the examples with non-typical poses
are in the lower end.

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=ex1_3.eps,height=1.6in}
\end{tabular}
\end{center}
\caption{Experiment 1.3. The $15$ images are 179 $\times$ 112.
Number of elements is 50.} \label{fig:ex1_3}
\end{figure}

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=ex1_4.eps,height=1in}
\end{tabular}
\end{center}
\caption{Experiment 1.4. The $12$ images are 120 $\times$ 167.
Number of elements is 50.} \label{fig:ex1_4}
\end{figure}

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=ex1_5.eps,height=1in}
\end{tabular}
\end{center}
\caption{Experiment 1.5. The $9$ images are 122 $\times$ 120. Number
of elements is 50.} \label{fig:ex1_5}
\end{figure}

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=ex1_6.eps,height=1in}
\end{tabular}
\end{center}
\caption{Experiment 1.6. The $11$ images are 133 $\times$ 140.
Number of elements is 50.} \label{fig:ex1_6}
\end{figure}

Figures (\ref{fig:ex1_2}) - (\ref{fig:ex1_6}) display more examples,
where the results are obtained by the same algorithm.

{\em Negative experience in Experiment 1.} This experiment requires that the training images are roughly aligned and the objects are in  the same pose. If this is not the case, our method cannot learn clean templates. Also, our method does not do well on objects with strong textures, such as zebras, leopards, tigers, giraffes, etc. The learning algorithm tends to sketch edges in textures.

Later we shall show that our method can be extended so that we can  learn from non-aligned images and find clusters in training images.

\subsection{Detection by inference algorithm}

This section studies the detection task using the inference algorithm based on sum-max maps.

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{cc}
\psfig{figure=ex2_1_original.eps,height=1in}&
\psfig{figure=ex2_1_sketch.eps,height=1in}\\
(a) & (b)
\end{tabular}
\end{center}
\caption{Experiment 2.1. (a) Testing image. The recognition
algorithm is run on 10 resolutions, from 270 $\times$ 360 to 432
$\times$ 536. (b) Superposed with sketch of the 60 elements of the
deformed active basis at the optimal resolution 378 $\times$ 504.
The bounding box of the template is 82 $\times$ 164. }
\label{fig:ex2_1}
\end{figure}

{\em Experiment 2}. Figure (\ref{fig:ex2_1}.a) displays the observed
image in Experiment 2.1. The deformable template is learned in
Experiment 1.1. The bounding box is 82 $\times$ 164. We run the
recognition algorithm on 10 resolutions of the testing image, from
270 $\times$ 360 to 432 $\times$ 536. Figure (\ref{fig:ex2_1}.b)
displays the superposed sketch of $(B_{\hat{x}_i, \hat{y}_i, s,
\hat{\alpha}_i}, i = 1, ..., n = 60)$ at the optimal resolution 378
$\times$ 504. Again, we let $h(r) = {\rm threshold}(r)$, the same as
in Experiment 1.

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{cc}
\psfig{figure=ex2_1_M2.eps,height=1.5in} &
\psfig{figure=ex2_1_S2.eps,height=1.5in}\\ (a) & (b)
\end{tabular}
\end{center}
\caption{Experiment 2.1. (a) $\M2$ scores at resolutions 1 to 10.
(b) $\S2$ map at the optimal resolution 7. } \label{fig:ex2_1_M2}
\end{figure}

Figure (\ref{fig:ex2_1_M2}.a) displays the $\M2$ scores at the 10
resolutions. The maximum is achieved at the 7th resolution. Figure
(\ref{fig:ex2_1_M2}.b) displays the $\S2$ map at this optimal
resolution. Recall that $\S2$ is the log of the likelihood ratio. If
we use the likelihood ratio, then the value at the maximum of the
$\S2$ map far exceeds the values of other pixels.

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=ex2_1_S1.eps,height=1.5in}
\end{tabular}
\end{center}
\caption{Experiment 2.1. Negative square root of $\S1$ maps at the
optimal resolution. There are 15 orientations. }
\label{fig:ex2_1_S1}
\end{figure}

Figure (\ref{fig:ex2_1_S1}) displays $-\S1(x, y, s, \alpha)^{1/2}$.
There are 15 orientations, so there are 15 $\S1$ maps.

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=ex2_1_M1.eps,height=1.5in}
\end{tabular}
\end{center}
\caption{Experiment 2.1. Negative square root of $\M1$ maps at the
optimal resolution. There are 15 orientations. }
\label{fig:ex2_1_M1}
\end{figure}

Figure (\ref{fig:ex2_1_M1}) displays $-\M1(x, y, s, \alpha)^{1/2}$.
There are 15 $\M1$ maps. These $\M1$ maps are blurred versions of
the $\S1$ maps. They themselves are not very meaningful unless they
are combined into the $\S2$ map using the active basis model learned
by Experiment 1.


\begin{figure}[!htb]
\begin{center}
\begin{tabular}{cc}
\psfig{figure=ex2_2_original.eps,height=1.5in}&
\psfig{figure=ex2_2_sketch.eps,height=1.5in}\\ (a) & (b)
\end{tabular}
\end{center}
\caption{Experiment 2.2. (a) Testing image. The algorithm is run on
10 resolutions, from 320 $\times$ 240 to 500 $\times$ 375. (b)
Superposed with sketch of 60 elements of the deformed active basis
at the optimal resolution 400 $\times$ 300. The bounding box of the
template is 180 $\times$ 180. } \label{fig:ex2_2}
\end{figure}

Figure (\ref{fig:ex2_2}.a) displays the observed image in Experiment
2.2. The deformable template is learned in Experiment 1.2. The
bounding box is 180 $\times$ 180. We run the recognition algorithm
on 10 resolutions of the testing image, from 320 $\times$ 240 to 500
$\times$ 375. Figure (\ref{fig:ex2_2}.b) displays the superposed
sketch of $(B_{\hat{x}_i, \hat{y}_i, s, \hat{\alpha}_i}, i = 1, ...,
n = 60)$ at the optimal resolution 400 $\times$ 300. Because the
pose of the bike in the testing image is not entirely the same as
the poses in the training images, the sketched bike in Figure
(\ref{fig:ex2_2}.b) misses part of the front wheel.

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=ex2_3_original.eps,height=1.5in}
\end{tabular}
\end{center}
\caption{Experiment 2.3. Testing image.  } \label{fig:ex2_3}
\end{figure}

Figure (\ref{fig:ex2_3}) displays the observed image in Experiment
2.3. The deformable template is learned in Experiment 1.3. The
bounding box is 179 $\times$ 112. We run the recognition algorithm
on 10 resolutions of the testing image, from 192 $\times$ 220 to 322
$\times$ 368.

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{ccccc}
\psfig{figure=ex2_3_scale1.eps,height=.5in} &
\psfig{figure=ex2_3_scale2.eps,height=.5in} &
\psfig{figure=ex2_3_scale3.eps,height=.5in} &
\psfig{figure=ex2_3_scale4.eps,height=.5in} &
\psfig{figure=ex2_3_scale5.eps,height=.5in} \\
\psfig{figure=ex2_3_scale6.eps,height=.5in} &
\psfig{figure=ex2_3_scale7.eps,height=.5in} &
\psfig{figure=ex2_3_scale8.eps,height=.5in} &
\psfig{figure=ex2_3_scale9.eps,height=.5in} &
\psfig{figure=ex2_3_scale10.eps,height=.5in}
\end{tabular}
\end{center}
\caption{Experiment 2.3. Superposed sketch of 50 elements of the
deformed active basis at each of the 10 resolutions, from 192
$\times$ 220 to 322 $\times$ 368. The bounding box is 179 $\times$
112. } \label{fig:ex2_3S}
\end{figure}

Figure (\ref{fig:ex2_3S}) displays the superposed sketch at each of
the 10 resolutions.

{\em Negative experience in Experiment 2.} Our method can sometimes be distracted by cluttered edges in the background. We may need to combine templates at multiple resolutions to overcome this problem.

\subsection{Classification by learning algorithm}

In this section, we evaluate our method on classification tasks and
compare the ROC curves.

{\em Scoring testing images.} We learn the active basis $\B = (B_i,
i = 1, ..., n)$ and estimate $\Lambda = (\lambda_i, i = 1, ..., n)$
(or $\Theta = (\theta_i, i = 1, ..., n)$ for active correlation)
from the training images. Then for each testing image $\I_m$, we can
compute its score $\S2_{m}$ according to Equation (\ref{eq:score}).
This can be done using the same algorithm as the shared sketch
algorithm in Subsection (\ref{sect:SS}), except that we remove the selection of $(x_i, y_i, \alpha_i)$ and the computation of $\lambda_i$ in the selection step, because $B_i$ and $\lambda_i$ are already available. We obtain the ROC curves based on $\S2_m$ scores.

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=ex3_sigmoid_all.eps,height=2in}
\end{tabular}
\end{center}
\caption{Experiment 3.1. Results obtained by fitting the active
basis model with sigmoid transformation. The saturation level is 6.
The $43$ images are 127 $\times$ 85.  The number of basis elements
$40$. } \label{fig:ex3_sigmoid}
\end{figure}

{\em {Experiment 3}}. Figure (\ref{fig:ex3_sigmoid}) displays the 43
positive training images in Experiment 3.1. It also displays the
corresponding $\B_{m} =(B_{m, i}, i = 1, ..., n)$, $n = 40$. The
results are obtained with $h(r) = {\rm sigmoid}(r)$. The saturation
level $\xi = 6$.  The $q(r)$ is pooled from 157 negative training
images. The bounding boxes of the training images and testing images
are all given. The sizes of the training and testing image patches
are all 127 $\times$ 85.

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=ex3_1.eps,height=2in}
\end{tabular}
\end{center}
\caption{Experiment 3.1. The solid curve is ROC of active basis
model with sigmoid transformation. The dashed curve is ROC of the
adaboost method where the weak classifiers are based on thresholding
the filter responses in $\S1$ maps. } \label{fig:ex_3_1}
\end{figure}

We then test on a separate data set with 88 positives and 474
negatives. The solid curve of Figure (\ref{fig:ex_3_1}) is the ROC
curve for the active basis model with sigmoid transformation. The
AUC (area under curve) is .977.

For comparison, we also display the ROC curve for the adaboost
method. The weak classifiers are obtained by thresholding the filter
responses $\S1_m(x, y, s, \alpha)$ computed in  Step 1(a) of the
shared sketch algorithm in Subsection (\ref{sect:SS}). The threshold
for each weak classifier is selected from a dense grid of points. We
select 80 weak classifiers using the adaboost method. The dashed
curve of Figure (\ref{fig:ex_3_1}) is the ROC curve for the adaboost
method. The AUC is .936.

We also try another type of weak classifiers, which are based on
thresholding the local maxima of filter responses, i.e., $\M1_m(x,
y, s, \alpha)$ obtained in Step 1(b) of the shared sketch algorithm.
We select 40 weak classifiers using the adaboost method. The AUC is
.943. If we increase the number of weak classifiers, the AUC starts
to decrease, suggesting that the classifier starts to overfit.

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{ccc}
\psfig{figure=ex3_sigmoid.eps,height=.4in}&
\psfig{figure=ex3_adaboost.eps,height=.4in}&
\psfig{figure=ex3_sigmoid_tiny.eps,height=.4in}\\
(a) & (b) & (c)
\end{tabular}
\end{center}
\caption{Experiment 3.1. (a) The template learned by the active
basis model. (b) Learned by the adaboost method.  (c) Learned from
tiny training images. } \label{fig:ex_3_2}
\end{figure}

Figure (\ref{fig:ex_3_2}.a) displays the template learned by the
active basis model. Figure  (\ref{fig:ex_3_2}.b) displays the
template learned by the adaboost method, where a weak classifier is
illustrated by a bar that symbolizes the corresponding Gabor
wavelet.

We also try the active basis model with transholding transformation.
The thresholding level is 16. The AUC is .941. The reason that it
does not do as well as sigmoid transformation may be due to the fact
that the thresholding level is too high, so that the strong edges in
the negative examples contribute too much to the log-likelihood
ratio score $\S2_m$. We then try the active basis model with
whitening transformation with $q(r)$ or $F(r)$ pooled from negative
training images. The AUC is .965.

We also try active correlation. The AUC for thresholding
transformation is .971. The AUC for whitening transformation is
.974. The improvement over log-likelihood may be due to the fact
that we take squared root of ${\rm whiten}(r)$ in active
correlation. This discounts strong edges in negative examples.


{\em Combining templates at multiple resolutions.} A low saturation
level in the sigmoid model is helpful to discount strong edges in
the background or in the negative examples. But this may cause
another problem: a large number of weak edges in the cluttered
background may collectively produce high $\S2$ score. One way to
alleviate this problem is to combine templates at multiple
resolutions, because the weak clutter edges cannot survive long in
the scale space \cite{Witkin}, so they will not produce high $\S2$
score at lower resolution.

In addition to the template displayed in Figure
(\ref{fig:ex3_sigmoid}), we learn another template using Gabor
wavelets at a scale twice as large as the original scale. We select
$n=20$ basis elements at this scale. Then in testing, we use the sum
of the $\S2_m$ scores obtained at the two scales. This leads to an
AUC = .982.

{\em Adaptive texture background.} In scoring a testing image
$\I_m$, the $\S2_m$ is computed by $\sum_{i=1}^{n} \log [p_i(c_{m,
i}; \hat{\lambda}_i)/q(c_{m, i})]$, according to Equation
(\ref{eq:score}), where $q(c)$ is pooled from natural images or
negative training images. We can change $q(c)$ to a background
texture model fitted specifically to $\I_m$. Specifically, for each
image $\I_m$, and for each orientation $\alpha$, we can fit a model
$q_m(c; \lambda_\alpha)$ to $\{\langle \I_m, B_{x, y, s,
\alpha}\rangle, \forall (x, y) \in D\}$. The distribution $q_m(c;
\lambda_\alpha)$ is again the exponential family model of the form
specified by Equation (\ref{eq:exponential}). The maximum likelihood
estimate $\hat{\lambda}_{\alpha} = \mu^{-1}(\sum_{(x, y) \in D} h
(|\langle \I_m, B_{x, y, s, \alpha}\rangle|^2)/|D|)$.
 Then we compute the $\S2_m$ score by
 \begin{eqnarray}
\S2_m &=& \sum_{i=1}^{n} \log [p_i(c_{m,i};
\hat{\lambda}_i)/q_m(c_{m,i}; \hat{\lambda}_{\alpha_i})] \nonumber \\
&=& \sum_{i=1}^{n} [(\hat{\lambda}_i h(|\langle \I_m, B_{m,
i}\rangle|^2) - \log Z(\hat{\lambda}_i)) \nonumber \\ && -
(\hat{\lambda}_{\alpha_i} h(|\langle \I_m, B_{m, i}\rangle|^2) -
\log Z(\hat{\lambda}_{\alpha_i}))],
\end{eqnarray}
where $\alpha_i$ is the orientation of $B_{m, i}$. That is, we score
the template matching against the adaptive texture background
summarized by $q_m(c; \hat{\lambda}_{\alpha_i})$. The marginal
distributions pooled over $\{\langle \I_m, B_{x, y, s,
\alpha}\rangle$, $\forall (x, y) \in D\}$ for different $\alpha$ and
$s$ have been used by Zhu, Wu, and Mumford \cite{Zhu_FRAME} for
modeling textures. This new $\S2_m$ score leads to an AUC = .983.

{\em Learning from single training image.} We can also train the
sigmoid model on a single positive training image. We use the 13th
image in Figure (\ref{fig:ex3_sigmoid}). It has a clean background.
We set $b_1 = b_2 = 0$ in learning, that is, we do not allow the
$B_i$ to perturb. Then we set $b_1$ and $b_2$ at their default
values for testing. The AUC is .926.

{\em Learning from tiny images.} We can also train the model on tiny
images. We zoom out the positive training images by a factor of 10,
so the image size is 12 $\times$ 8. In order to learn from such tiny
images, we then zoom in these 12 $\times$ 8 images by a factor of
10, so the image size becomes 120 $\times$ 80. Then we apply the
shared sketch algorithm on these $120 \times 80$ images. The learned
template is shown in Figure (\ref{fig:ex_3_2}.c). Then we test on
the testing images, which also go through the same zooming out and
zooming in operations. The AUC is .957.

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=ex3_sigmoid_all_2.eps,height=2in}
\end{tabular}
\end{center}
\caption{Experiment 3.2.  The $31$ images are 150 $\times$ 120. The
number of basis elements is $40$. } \label{fig:ex3_sigmoid2}
\end{figure}


We also did another experiment on horse images. Figure
(\ref{fig:ex3_sigmoid2}) displays the learning results of Experiment
3.2. There are 31 positive training examples. They are 150 $\times$
120. The number of basis elements is $n = 40$. There are 249
positive testing images. We use the same sets of negative training
and testing images as in Experiment 3.1 (with resizing to match the
size of the positive examples). The AUC scores are as follows:
sigmoid = .987, whiten = .982, whiten with active correlation =
.984, threshold = .981, threshold with active correlation = .985,
combining two templates = .988, adaptive texture background = .989,
single image learning = .980, tiny image learning = .956, adaboost
method by thresholding the $\S1$ maps = .937, adaboost method by
thresholding the $\M1$ maps = .957.

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=ex3_sigmoid_all_3.eps,height=2in}
\end{tabular}
\end{center}
\caption{Experiment 3.3. The $33$ images are 150 $\times$ 100. The
number of basis elements is $40$. } \label{fig:ex3_sigmoid3}
\end{figure}

We did a third experiment on butterfly images. Figure
(\ref{fig:ex3_sigmoid3}) displays the learning results of Experiment
3.3. There are 33 positive training examples. They are 150 $\times$
100. The number of basis elements is $n = 40$. There are 191
positive testing images. We use the same sets of negative training
and testing images as in Experiment 3.1. The AUC for sigmoid model
is .999. The AUC for adaboost by thresholding $\M1$ maps is .994.

Experiment 3 focuses on the situation where the number of positive
training images is relatively small. Compared to adaboost, the
active basis model applies non-maximum suppression instead of re-weighting before selecting the next basis element. The active basis model does not require negative examples (except pooling a marginal histogram from natural images). Moreover, the generative model can be conveniently used for unsupervised learning as we shall show in the following sections.


\subsection{Geometric transformation of template}

 Given a template $\B = (B_{i} = B_{x_i, y_i, s, \alpha_i}, i = 1, ..., n)$, we can transform this template by dilation, rotation, and changing the aspect ratio. This amounts to simple transformations of $(x_i, y_i, \alpha_i, i = 1, ..., n)$.

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{ccc}
\psfig{figure=ex8_1.eps,height=.7in}&
\psfig{figure=ex8_2.eps,height=.7in}&
\psfig{figure=ex8_3.eps,height=.7in}\\
(a) & (b) & (c)
\end{tabular}
\end{center}
\caption{Experiment 4.1. (a) The number of elements is 60. The image size is 252 $\times$ 320. The scale factor is 1.4. The rotation is 1 $\times \pi/15$. The aspect factor is 0.9.  (b) The image size is 200 $\times$ 250. The scale factor is 1.4. The rotation is 1 $\times \pi/15$. The aspect factor is 1. (c) The image size is 248 $\times$ 232. The scale factor is 1.2. The rotation is -1 $\times \pi/15$. The aspect factor is 0.6. } \label{fig:Gt}
\end{figure}

Figure (\ref{fig:Gt}) shows three examples. The bicycle template is learned in Experiment 1, where we use the sigmoid transformation. We then transform it into a collection of templates as different scales, orientations, and aspect ratios. After that, we use these templates to detect the objects by template matching, as in Experiment 2. Finally, we choose the transformed template that gives the best match, and superpose it on the input image.

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{cc}
\psfig{figure=ex8_horse_1.eps,height=.7in}&
\psfig{figure=ex8_horse_2.eps,height=.7in}\\
(a) & (b)
\end{tabular}
\end{center}
\caption{Experiment 4.2. (a) The number of elements is 40. The image size is 138 $\times$ 168. The scale factor is 1. The rotation is 0 $\times \pi/15$. The aspect factor is 0.8. (b) The image size is 192 $\times$ 144. The scale factor is 1. The rotation is 4 $\times \pi/15$. The aspect factor is 1.3. } \label{fig:Gt1}
\end{figure}

Figure (\ref{fig:Gt1}) shows another example with the horse template, learned in Experiment 3.

{\em Negative experience in Experiment 4.} We encountered some difficulty with the bicycle template. When the viewing distance is close, the size of one wheel can be larger than the size of the other wheel, so a single scale factor does not give very good fit. Also, the frontal wheel may turn to a different direction than the back wheel.

The above difficulty suggests that we should better split the bicycle template into two part-templates, and each part-template has its own geometric transformation. We shall study the composition of multiple part-templates later.


\section{Learning from Non-aligned Images}

In this section, we study the problem of learning from images where the objects are of unknown locations and scales.

\subsection{Image alignment}

For the training image patches $\{\I_m, m = 1, ..., M\}$ defined on the same bounding box, such as those in the previous section, we can define their overall alignment by
\begin{eqnarray}
{\rm ALIGN}(\I_m, m = 1, ..., M) =
\sum_{m=1}^{M} {\rm MATCH}(\I_m, \B),
\end{eqnarray}
where $\B$ is the template learned from the image patches, and ${\rm MATCH}(\I_m, \B)$ is the template matching score defined by either (\ref{eq:match}) for log-likelihood or (\ref{eq:match1}) for active correlation. The computation is carried out by the
shared sketch algorithm in Subsection (\ref{sect:SS}), and ${\rm
ALIGN}(\I_m, m = 1, ..., M) = \sum_{m=1}^{M} \S2_m$, where the
$\S2_m$ scores are output by the algorithm.

When the training images $\{\I_m, m = 1, ..., M\}$ are of different sizes, and the objects appear at different locations in the training
images, we need to infer the unknown locations. Let ${\rm box}(x,
y)$ be the rectangular bounding box of the template centered at $(x, y)$. For an image $\I$, let $\I[{\rm box}(x, y)]$ be the image patch cropped from the image $\I$ within ${\rm box}(x, y)$. We want to maximize the alignment score
\begin{eqnarray}
{\rm ALIGN}\left(\I_m[{\rm box}(x_m, y_m)], m = 1, ..., M\right),
\end{eqnarray}
where $(x_m, y_m)$ is the unknown location of the bounding box in
$\I_m$.

The alignment score can be maximized by a greedy algorithm that
iterates the following two steps:

(1) Supervised learning: Given $\{(x_m, y_m), m = 1, ..., M\}$, estimate $(\B, \Lambda)$ from $\{\I_m[{\rm box}(x_m, y_m)]$, $m = 1, ..., M\}$ using the shared sketch algorithm in Subsection (\ref{sect:SS}).

(2) Detection: Given $(\B, \Lambda)$, estimate $(x_m, y_m)$ from each $\I_m$ using the inference algorithm in Subsection
(\ref{sect:recognition}). $(x_m, y_m)$ achieves the maximum of the
$\S2$ map.

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{cccccc}
\psfig{figure=ex5_cat_sym.eps,height=.6in} &
\psfig{figure=ex5_cat1.eps,height=.6in} &
\psfig{figure=ex5_cat2.eps,height=.6in} &
\psfig{figure=ex5_cat3.eps,height=.6in} &
\psfig{figure=ex5_cat4.eps,height=.6in} &
\psfig{figure=ex5_cat5.eps,height=.6in}
\end{tabular}
\end{center}
\caption{Experiment 5a.1. The bounding box of the first image is
given. The size of the bounding box is $136 \times 140$. The number
of elements in the active basis is 60. } \label{fig:ex5_1}
\end{figure}

{\em Experiment 5a:} In this experiment, we initialize the
algorithm by specifying the bounding box for the first training image. Then we estimate $(\B, \Lambda)$ from this single image patch. In learning from the single image patch, we set $b_1 = b_2 = 0$, that is, we do not allow the elements $B_i$ to perturb. After that, we
reset $b_1$ and $b_2$ to their default values, and iterate Step (2)
and Step (1) described above.


\begin{figure}[!htb]
\begin{center}
\begin{tabular}{cccccc}
\psfig{figure=ex5_wolf_sym.eps,height=.5in} &
\psfig{figure=ex5_wolf1.eps,height=.5in} &
\psfig{figure=ex5_wolf2.eps,height=.5in} &
\psfig{figure=ex5_wolf3.eps,height=.5in} &
\psfig{figure=ex5_wolf4.eps,height=.5in} &
\psfig{figure=ex5_wolf5.eps,height=.5in}
\end{tabular}
\end{center}
\caption{Experiment 5a.2. The bounding box is $116 \times 116$.
Number of elements is 50. } \label{fig:ex5_2}
\end{figure}

In Step (2), we search over 9 different resolutions (from 0.8 to 1.2
times the input image size). We crop $\I_m[{\rm box}(x_m, y_m)]$
from the optimal resolution.



\begin{figure}[!htb]
\begin{center}
\begin{tabular}{cccccc}
\psfig{figure=ex5_swan_sym.eps,height=.5in} &
\psfig{figure=ex5_swan1.eps,height=.5in} &
\psfig{figure=ex5_swan2.eps,height=.5in} &
\psfig{figure=ex5_swan3.eps,height=.5in} &
\psfig{figure=ex5_swan4.eps,height=.5in} &
\psfig{figure=ex5_swan5.eps,height=.5in}
\end{tabular}
\end{center}
\caption{Experiment 5a.3.
The bounding box is $115 \times 161$. Number of elements is 50.
} \label{fig:ex5_3}
\end{figure}

We run the algorithm for three iterations. Figures (\ref{fig:ex5_1})
- (\ref{fig:ex5_6}) display some examples.


\begin{figure}[!htb]
\begin{center}
\begin{tabular}{cccccc}
\psfig{figure=ex5_pigeon_sym.eps,height=.5in} &
\psfig{figure=ex5_pigeon1.eps,height=.5in} &
\psfig{figure=ex5_pigeon2.eps,height=.5in} &
\psfig{figure=ex5_pigeon3.eps,height=.5in} &
\psfig{figure=ex5_pigeon4.eps,height=.5in} &
\psfig{figure=ex5_pigeon5.eps,height=.5in}
\end{tabular}
\end{center}
\caption{Experiment 5a.4.
The bounding box is $94 \times 99$. Number of elements is 30.
} \label{fig:ex5_4}
\end{figure}


\begin{figure}[!htb]
\begin{center}
\begin{tabular}{cccccc}
\psfig{figure=ex5_horse_sym.eps,height=.5in} &
\psfig{figure=ex5_horse1.eps,height=.5in} &
\psfig{figure=ex5_horse2.eps,height=.5in} &
\psfig{figure=ex5_horse3.eps,height=.5in} &
\psfig{figure=ex5_horse4.eps,height=.5in} &
\psfig{figure=ex5_horse5.eps,height=.5in}
\end{tabular}
\end{center}
\caption{Experiment 5a.5.
The bounding box is $94 \times 138$. Number of elements is 60.
} \label{fig:ex5_5}
\end{figure}


\begin{figure}[!htb]
\begin{center}
\begin{tabular}{cccccc}
\psfig{figure=ex5_deer_sym.eps,height=.5in} &
\psfig{figure=ex5_deer1.eps,height=.5in} &
\psfig{figure=ex5_deer2.eps,height=.5in} &
\psfig{figure=ex5_deer3.eps,height=.5in} &
\psfig{figure=ex5_deer4.eps,height=.5in} &
\psfig{figure=ex5_deer5.eps,height=.5in}
\end{tabular}
\end{center}
\caption{Experiment 5a.6. The bounding box is $134 \times 148$.
Number of elements is 60. } \label{fig:ex5_6}
\end{figure}

{\em Experiment 5b}. This experiment is about pairwise alignment, that is, we align two images so that they overlap on the common parts. We learn a template from the first image with no activity and without any given bounding box. Then we restore the activity and scan the learned template on the second image using the inference algorithm as in Experiment 2. When the template is partially out of the bound of the second image, we set the responses of those elements that are out of the bound to 0. We scan the template over 7 resolutions of the second image, from .7 to 1.3 times the size of the input image. At the optimal location and resolution, we deform the template and sketch the second image. In this experiment, we use active correlation for learning and inference.


\begin{figure}[!htb]
\begin{center}
\begin{tabular}{cc}
\psfig{figure=ex5b_1_original.eps,width=1.2in}&
\psfig{figure=ex5b_1_find.eps,width=1.2in}\\
(a) & (b)
\end{tabular}
\end{center}
\caption{Experiment 5b.1. The number of elements is 50.} \label{fig:A1}
\end{figure}

Figure (\ref{fig:A1}) displays the pairwise alignment. (a) shows the first image and the template learned from this image. (b) shows the second image, and the sketch of this image by the learned template, where the intensities of the elements are proportional to the square roots of the corresponding responses from the second image. Figure (\ref{fig:A2}) displays the alignment where we switch the order of the two images.

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{cc}
\psfig{figure=ex5b_2_original.eps,width=1.2in}&
\psfig{figure=ex5b_2_find.eps,width=1.2in}\\
(a) & (b)
\end{tabular}
\end{center}
\caption{Experiment 5b.2. The number of elements is 50.} \label{fig:A2}
\end{figure}

Figure (\ref{fig:A3}) displays the result of aligning the first frame of a vide sequence to the 15th frame of the same sequence. The image sequence is cropped from the PETS 2006 benchmark data and is to be used again in Experiment 6.

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{cc}
\psfig{figure=ex5b_5_original.eps,width=1.2in}&
\psfig{figure=ex5b_5_find.eps,width=1.2in}\\
(a) & (b)
\end{tabular}
\end{center}
\caption{Experiment 5b.3. The number of elements is 80.} \label{fig:A3}
\end{figure}

{\em Experiment 5c}. This is a repetition of Experiment 5a, except that we do not assume that the bounding box of the object in the first image is given. We simply start from the template learned from the whole image of the first example. Figure (\ref{fig:A4}) displays two examples. In each example, the first template is learned from the first image, and the template serves as the initialization of the algorithm. The second template is produced after 3 iterations of the algorithm used in Experiment 5a. As in the pairwise alignment, we allow the template to be partially out of the bounds of the images in the detection step.


\begin{figure}[!htb]
\begin{center}
\begin{tabular}{cc|cc}
\psfig{figure=ex5c_wolf_sym0.eps,height=.7in}&
\psfig{figure=ex5c_wolf_sym.eps,height=.7in} &
\psfig{figure=ex5c_horse_sym0.eps,height=.5in}&
\psfig{figure=ex5c_horse_sym.eps,height=.5in}
\end{tabular}
\end{center}
\caption{Experiment 5c. In each example, the first template is the starting template. The second template is learned after 3 iterations. The number of elements of the active basis is 60 in the left example, and 50 in the right example.} \label{fig:A4}
\end{figure}


{\em Negative experience in Experiments 5a-c.} When there are cluttered edges in the background, the detection step may fail to locate the objects. When the objects have large deformations or pose changes, the learned template may not be clean, and may fail to sketch the objects in the training images correctly. In Experiments 5b and 5c, if the objects do not occupy significant portions of the training images, our method may fail to establish correct alignment.

\subsection{Learning part-templates}

The algorithm in Experiment 5a can be used to learn part-templates from training images. In Experiment 5d, we start from a large number of patches cropped from the training images, and for each starting patch, we learn a templates using the same algorithm as in Experiment 5a. Here we use active correlation instead of the log-likelihood for learning and detection.

Then we select the first K templates with the highest alignment scores. We did not perform spatial inhibition between the part-templates. After that, we double the sizes of the input images, and use the same procedure to learn part-templates at a higher resolution.



\begin{figure}[!htb]
\begin{center}
\begin{tabular}{cccc}
\psfig{figure=ex5d_wordcar_sym1.eps,height=.4in}&
\psfig{figure=ex5d_wordcar1_1.eps,height=.4in}&
\psfig{figure=ex5d_wordcar1_2.eps,height=.4in}&
\psfig{figure=ex5d_wordcar1_3.eps,height=.4in}\\
\psfig{figure=ex5d_wordcar_sym2.eps,height=.4in}&
\psfig{figure=ex5d_wordcar2_1.eps,height=.4in}&
\psfig{figure=ex5d_wordcar2_2.eps,height=.4in}&
\psfig{figure=ex5d_wordcar2_3.eps,height=.4in}\\
\psfig{figure=ex5d_wordcar_sym3.eps,height=.4in}&
\psfig{figure=ex5d_wordcar3_1.eps,height=.4in}&
\psfig{figure=ex5d_wordcar3_2.eps,height=.4in}&
\psfig{figure=ex5d_wordcar3_3.eps,height=.4in}
\end{tabular}
\end{center}
\caption{Experiment 5d. The top three part-templates. The size of the bounding box is 100 $\times$ 100. The number of elements is 40. The allowed activity in location is up to 3 pixels. The allowed activity in orientation is up to $\pi/15$. } \label{fig:P1}
\end{figure}

Figure (\ref{fig:P1}) displays the top three part-templates learned from three car images. Because of the large deformations in these three cars, it is impossible to learn a common template for the whole cars, but it is still possible to learn meaningful part templates that correspond to frontal, middle and rear parts of the cars.

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{cccc}
\psfig{figure=ex5d_wordcarL_sym1.eps,height=.4in}&
\psfig{figure=ex5d_wordcarL1_1.eps,height=.4in}&
\psfig{figure=ex5d_wordcarL1_2.eps,height=.4in}&
\psfig{figure=ex5d_wordcarL1_3.eps,height=.4in}\\
\psfig{figure=ex5d_wordcarL_sym2.eps,height=.4in}&
\psfig{figure=ex5d_wordcarL2_1.eps,height=.4in}&
\psfig{figure=ex5d_wordcarL2_2.eps,height=.4in}&
\psfig{figure=ex5d_wordcarL2_3.eps,height=.4in}\\
\end{tabular}
\end{center}
\caption{Experiment 5d. The top two part-templates learned after the sizes of the input images are doubled. The parameters are the same as in Figure (\ref{fig:P1})} \label{fig:P2}
\end{figure}

Figure (\ref{fig:P2}) displays the top two part-templates learned from these three images after we resize these images by a factor of 2.

{\em Negative experience in Experiment 5d.} When the part-template is small relative to the whole objects, the method often fails to establish correct correspondence among the images.

The above difficulty suggests that we should add constraints for more reliable learning of the parts. If the bounding boxes are given as in Experiment 1, we can restrict the ranges of movements of parts in the training images, so that in the detection step, we do not need to search over the whole images. If the bounding boxes are not given as in Experiment 5a, we can simultaneously learn multiple parts while restricting their relative positions, and this is very much like a recursion of Experiment 5a or 5c.



\subsection{Learning moving template from motion sequence}

Our method can also be used to learn a moving deformable template
from a video sequence. Let $(\I_{t}, t = 1, ..., M)$ be a sequence
of frames of an object shape that is moving at a speed $v = (v_x,
v_y)$. We can estimate $v$ and learn a template of the object shape
simultaneously.

At the true speed $v=(v_x, v_y)$, let $\J_{t}^{(v)}(x, y) =
\I_t(x+v_x t, y+v_y t)$, i.e., for frame $t$, we shift the image
lattice back by $vt$, then the object shapes in $\{\J_{t}^{(v)}, t =
1, ..., M\}$ will be well aligned. If we apply the shared sketch
algorithm to $\{\J_t^{(v)}\}$, we shall learn a clean template that
has a high alignment score. We can try all possible $v$, and choose
the $v$ that achieves the maximum alignment score, i.e., we maximize
\begin{eqnarray}
{\rm ALIGN}\left(\J_t^{(v)}, t = 1, ..., M\right) \label{eq:moving}
\end{eqnarray}
over $v$. This is actually a simpler problem than learning from
non-aligned images.

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=ex6_man_score.eps,height=2in}
\end{tabular}
\end{center}
\caption{Experiment 6.1. Alignment scores at different speeds of the
optimal direction. } \label{fig:ex6_1_score}
\end{figure}

In our experiment, we use active correlation to evaluate the
alignment score (\ref{eq:moving}). Before computing this score, we
need to do background subtraction. First, we compute the background
$\S1$: $\S1_0(x, y, s, \alpha) = \sum_{t=1}^{M} \S1_t(x, y, s,
\alpha)/M$. Then we modify $\S1_t(x, y, s, \alpha) \leftarrow
[\S1_t(x, y, s, \alpha) - \S1_0(x, y, s, \alpha)]_{+}$, where
$[r]_{+} = r$ if $r > 0$ and $[r]_+ = 0$ otherwise. For each $v$, we
compute the alignment score using the background subtracted $\S1$
maps.

{\em Experiment 6.} We learn the moving template from a sequence of
19 frames of size 204 $\times$ 258. The image sequence is cropped
from the PETS 2006 benchmark data. We try 5 different directions
$v_x/v_y$, and at each direction, we try 7 different speeds. Figure
(\ref{fig:ex6_1_score}) displays the alignment scores at different
speeds of the optimal direction.

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{ccc|c|ccc}
\psfig{figure=ex6_man_templates1.eps,height=.5in}&
\psfig{figure=ex6_man_templates2.eps,height=.5in}&
\psfig{figure=ex6_man_templates3.eps,height=.5in}&
\psfig{figure=ex6_man_templates4.eps,height=.5in}&
\psfig{figure=ex6_man_templates5.eps,height=.5in}&
\psfig{figure=ex6_man_templates6.eps,height=.5in}&
\psfig{figure=ex6_man_templates7.eps,height=.5in}\\
(1)&(2)&(3)&(4)&(5)&(6)&(7)
\end{tabular}
\end{center}
\caption{Experiment 6.1. Learned templates at different speeds. (4)
is the one at the optimal speed. There are 19 frames of size 204
$\times$ 258. Number of elements is 70.  }
\label{fig:ex6_1_templates}
\end{figure}

Figure (\ref{fig:ex6_1_templates}) displays the templates learned at
different speeds of the optimal direction. Template (4) is learned
at the optimal speed. The number of basis elements is 70.


\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=ex6_man_sequence.eps,height=2in}
\end{tabular}
\end{center}
\caption{Experiment 6.1. Learned template and superposed sketch for
each frame at the optimal speed and direction. }
\label{fig:ex6_1_sequence}
\end{figure}

Figure (\ref{fig:ex6_1_sequence}) displays the learned template and
the superposed sketch for each frame at the optimal speed and
direction.

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=ex6_woman_sequence.eps,height=.9in}
\end{tabular}
\end{center}
\caption{Experiment 6.2. Learned template and superposed sketches at
the optimal speed. The image frames are 180 $\times$ 186. Number of
elements is 80.  } \label{fig:ex6_2_sequence}
\end{figure}

Figure (\ref{fig:ex6_2_sequence}) displays results for another example.

The learned templates can be used for tracking. Compared to fully
unsupervised learning, it provides a more reliable way to learn
templates and their parts.

\section{Clustering and Local Learning}

In this section, we study the problem of clustering, where we need to learn multiple templates from the training sample, which is a mixture of different poses or different categories.

\subsection{EM and K-mean}

{\em Mixture model and EM}. Suppose there are $K$ categories, and each category $k$ can be described by an active basis model $\B^{(k)} = (B_{i}^{(k)}, i = 1, ..., n)$ and $\Lambda^{(k)} =
(\lambda_{i}^{(k)}, i = 1, ..., n)$. Let $\rho^{(k)}$ be the
probability that a training image $\I_m$ comes from cluster $k$, $k
= 1, ..., K$. So $\I_m \sim \sum_{k=1}^{K} \rho^{(k)} p^{(k)}
(\I_{m} \mid \B_{m}^{(k)})$, i.e., a mixture distribution, where
each $p^{(k)} (\I_{m} \mid \B_{m}^{(k)})$ is modeled as in
Subsection (\ref{sect:model}).

We can learn $\{\rho^{(k)}, \B^{(k)}, \Lambda^{(k)}, k = 1, ...,
K\}$ by the EM algorithm. For each image $\I_m$, we define
$(z_{m}^{(k)}, k = 1, ..., K)$ as an indictor vector, where
$z_{m}^{(k)} = 1$ if $\I_m$ comes from cluster $k$, otherwise
$z_{m}^{(k)} = 0$.

{\em E-step.} For each $m = 1, ..., M$ and $k = 1, ..., K$, we
impute
\[
z_{m}^{(k)} = \frac{\rho^{(k)} \exp\{\S2_m^{(k)}\}}{\sum_{k=1}^{K}
\rho^{(k)} \exp\{\S2_m^{(k)}\}}. \label{eq:Estep}
\]
This is a soft classification based on the current models of the
clusters, where each $z_{m}^{(k)}$ becomes a fraction. The
$\S2_m^{(k)}$ scores are obtained in the M-step.

{\em M-step.} For each $k = 1, ..., K$, we learn $\B^{(k)}$ and
$\Lambda^{(k)}$ according to the shared sketch algorithm in
Subsection (\ref{sect:SS}). We only need to make the following
changes to the original version of the learning algorithm.

(1) In Step 2(a), we find $(x_i, y_i, \alpha_i)$ by maximizing
$\sum_{m=1}^{M} {z}_{m}^{(k)}h(\M1(x, y, s, \alpha))$, which is a
weighted sum.

(2) In Step 2(c), we compute $\hat{\lambda}_i$ by
\begin{eqnarray}
\hat{\lambda}_i = \mu^{-1}(\frac{\sum_{m=1}^{M} h(r_{m, i})
{z}_{m}^{(k)}}{\sum_{m=1}^{M} {z}_{m}^{(k)}}),  \label{eq:Mstep}
\end{eqnarray}
that is, we match $\mu(\lambda_i)$ to the weighted average.

(3) After we finish the algorithm,  we attach a superscript $(k)$ to
the resulting $\S2_m$ and $\B$.

We initialize the algorithm by randomly generating
$\{z_{m}^{(k)}\}$, and then iterate the M-step and the E-step. We
stop the algorithm after a few iterations. Then we classify $\I_m$
to the cluster $k_*$ that maximizes $z_{m}^{(k)}$ over all $k = 1,
..., K$.


\begin{figure}[!htb]
\begin{center}
\begin{tabular}{ccc|ccc}
\psfig{figure=ex4_clusterTemplates1.eps,height=.4in}&
\psfig{figure=ex4_clusterTemplates2.eps,height=.4in}&
\psfig{figure=ex4_clusterTemplates3.eps,height=.4in}&
\psfig{figure=ex4_Kmean_clusterTemplates1.eps,height=.4in}&
\psfig{figure=ex4_Kmean_clusterTemplates2.eps,height=.4in}&
\psfig{figure=ex4_Kmean_clusterTemplates3.eps,height=.4in}\\ &(a)& & &(b)&
\end{tabular}
\end{center}
\caption{Experiment 7.1. Learned templates from the mixture of the
three sets of positive training images in Experiment 3. Image size
is 120 $\times$ 150. Number of elements in each template is 40.
Number of iteration is 4. (a) EM. (b) K-mean. } \label{fig:ex4_1}
\end{figure}

{\em Experiment 7.} Figure (\ref{fig:ex4_1}.a) displays the learned
templates $\B^{(k)}, k = 1, ..., 3$ from the mixture of the three
sets of positive training images in Experiment 3. The image size is
120 $\times$ 150. For the head and shoulder images and the butterfly
images, we resize them to be 120 $\times$ 150. The number of
elements in each $\B^{(k)}$ is 40. We run the EM algorithm for 4
iterations. The EM algorithm can easily separate the three clusters.
Only one error is made, where a horse image is classified as a
head-shoulder image.

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{cc|cc}
\psfig{figure=ex4_2_clusterTemplates1.eps,height=.4in}&
\psfig{figure=ex4_2_clusterTemplates2.eps,height=.4in}&
\psfig{figure=ex4_2_Kmean_clusterTemplates1.eps,height=.4in}&
\psfig{figure=ex4_2_Kmean_clusterTemplates2.eps,height=.4in}\\
 (a)& & (b)&
\end{tabular}
\end{center}
\caption{Experiment 7.2. Learned templates for 58 images of horses
facing two different directions. Image size is 120 $\times$ 150.
Number of elements in each template is 50. (a) EM with 4 iterations.
(b) K-mean with 8 iterations.   } \label{fig:ex4_2}
\end{figure}

Figure (\ref{fig:ex4_1}.a) displays the learned templates $\B^{(k)},
k = 1, 2$ from a mixture of 58 images of horses facing different
directions. The image size is 120 $\times$ 150. The number of
elements in each $\B^{(k)}$ is 50. We run the EM algorithm for 4
iterations. The EM algorithm easily separates the two clusters.


{\em K-mean clustering.} We can pose the clustering problem as the
following alignment problem: find $\{(z_{m}^{(k)}, k = 1,
..., K), m = 1, ..., M\}$ to maximize
\begin{eqnarray}
\sum_{k=1}^{K} {\rm ALIGN}\left(\I_m, z_{m}^{(k)} = 1\right),
\label{eq:Kmean}
\end{eqnarray}
where ${\rm ALIGN}\{\I_m, z_{m}^{(k)} = 1\}$ is the alignment score
of the $k$-th cluster. The computation of the alignment score also
produces the template $\B^{(k)}$ for the $k$-th cluster. If we use
active correlation to learn the template for each cluster and score
the multiple alignment within each cluster, then the learned
$(\B^{(k)}, \Theta^{(k)})$ gives us an active mean vector $V^{(k)} =
\sum_{i=1}^{n} \theta_i^{(k)} B_i^{(k)}$ for each cluster. The mean
vector $V^{(k)}$ points to the center of the $k$-th cluster. This
leads to a K-mean algorithm. It is a greedy algorithm that maximizes
(\ref{eq:Kmean}), and it iterates the following two steps:

(1) Given $\{(z_{m}^{(k)}, k = 1, ..., K), m = 1, ..., M\}$,
estimate the mean vector $(\B^{(k)}, \Theta^{(k)})$ from $\{\I_m,
z_{m}^{(k)} = 1\}$ for each $k = 1, ..., K$.

(2) Given $\{\B^{(k)}, \theta^{(k)}, k = 1, ..., K\}$, classify each
image $\I_m$ to a cluster $k_{*}$ that maximizes $\langle\I_m \mid
V_m^{(k)}\rangle$ over all $k = 1, ..., K$. Set $z_{m}^{(k_*)} = 1$,
and set $z_{m}^{(k)} = 0$ for $k \neq k_{*}$.

The implementation of this K-mean algorithm is similar to the EM
algorithm. We only need to make the following modifications.

(1) Change the E-step: we let $z_{m}^{(k_*)} = 1$ if $k_*$ achieves
the maximum of  $\S2_{m}^{(k)}$ among all $k = 1, ..., K$, and we
set the the rest of $z_{m}^{(k)}$ to be 0.

(2) Change the M-step: for each $k = 1, ..., K$, we compute
$\S2_{m}$ and estimate $\B$ and $\Theta$ for each cluster $k$ using
the shared sketch algorithm that maximizes the active correlation.
The algorithm is described in Section (\ref{sect:AC}).

Figure (\ref{fig:ex4_1}.b) and Figure (\ref{fig:ex4_2}.b) display
the learned templates $\B^{(k)}, k = 1, ..., K$ using K-mean
algorithm. We initialize the algorithm with random
$\{z_{m}^{(k)}\}$.

We also did a third experiment where we mix the positive training
examples of head-shoulder images and negative training examples. The
EM and K-mean algorithms can still separate out many of the positive
training examples, although they also mistakenly include some
negative examples into the positive cluster.

{\em Negative experience in Experiment 7.} When the object shapes of different categories are not very different, our method often fails to distinguish them if we start from random clustering.

The above difficulty is not caused by the model or the EM or K-mean iteration, but mainly by the fact that random clustering gives poor initialization.

\subsection{Local learning of representative templates}

The EM and the K-mean clustering methods assume that the number of clusters is given. They also assume that the clusters should be  distinct from each other. The results of the algorithms seem to be severely dependent on initialization. To address these problems, we develop a local learning scheme by modifying the K-mean method in the previous subsection. In this scheme, we learn a local representative template around each input image, by recruiting a small number of nearest neighbors. With the local templates and their nearest neighbors, we can then perform trimming and merging for clustering.

In Experiment 8, we learn local representatives in an ensemble of 123 images of animal heads. Around each image, we learn a template from this image and its $K$ nearest neighbors, using an iterative scheme based on the active correlation. Specifically, we first learn a template for this single image where no activity is allowed. Then using the learned template, we find $K$ nearest neighbors by active correlation, where activity is restored. After that we re-learn the template from this image and its $K$ nearest neighbors. The weight for each of the $K$ nearest neighbors is 1, whereas the weight for this image is $\rho K$. In our experiment, $K = 5$, and $\rho$ = 1/3. Then again, we find the $K$ nearest neighbors, and iterate. After learning all the templates, we trim them to satisfy the constraint that the neighbors of the remaining templates should not overlap (this may be too aggressive). This leaves 16 templates.


\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=ex9_1_represent_trim.eps,height=.8in}
\end{tabular}
\end{center}
\caption{Experiment 8. The 16 representative templates locally learned by active correlation. They are ordered by the total alignment scores. Image size is 100 $\times$ 100. Number of elements is 40. Number of iterations is 3 for learning each template. The allowed activity of location is up to 2 pixels. The allowed activity of orientation is up to $\pi/15$.} \label{fig:L1}
\end{figure}

Figure (\ref{fig:L1}) shows the 16 templates obtained by local learning. They are ordered by the alignment scores computed from the $K$ nearest neighbors.

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=ex9_1_97.eps,height=.4in}\\
\psfig{figure=ex9_1_88.eps,height=.4in}\\
\psfig{figure=ex9_1_16.eps,height=.4in}\\
\psfig{figure=ex9_1_40.eps,height=.4in}
\end{tabular}
\end{center}
\caption{Experiment 8. The top four templates and their neighbors. } \label{fig:L2}
\end{figure}

Figure (\ref{fig:L2}) shows the top four templates. For each template, we also show its five nearest neighbors and their sketches.
\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=ex9_1_112.eps,height=.4in}\\
\psfig{figure=ex9_1_54.eps,height=.4in}\\
\psfig{figure=ex9_1_71.eps,height=.4in}\\
\psfig{figure=ex9_1_38.eps,height=.4in}
\end{tabular}
\end{center}
\caption{Experiment 8. Another four templates and their neighbors. } \label{fig:L3}
\end{figure}

Figure (\ref{fig:L3}) shows another four templates. Although these templates are less clear and clean than the top four templates, they still manage to recruit nearest neighbors that are from identical categories.

We can then use these 16 templates to initialize the K-mean clustering. The results are meaningful. With local learning, we can also perform hierarchical clustering by merging clusters.

Such a local learning scheme identifies local dimensions in the image ensemble. It also furnishes us with a local metric measured by active correlation. It is unclear whether we could use support vector machine \cite{SVM} to select the positive as well as negative templates, and use these support deformable templates for classification.

We have not had much experience with local learning. We tried it on digits, and it produced meaningful results, although there is some mixing between different categories, such as ``4'' or ``9''.


\section{Synthesis by Multi-scale Gabors and DoGs}

In most of the experiments so far, we only analyze the images at a single scale. This is sometimes not enough for detection and classification, where we should learn templates at multiple scales and combine the template matching scores. For instance, Figure (\ref{fig:deer}) displays the template learned from the deer images, where the basis elements are twice the size of the elements in Experiment 1. We can combine these two templates for detection and classification.

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=deerL.eps,height=.9in}
\end{tabular}
\end{center}
\caption{The size of the Gabor wavelet elements are twice the size of the elements in Experiment 1. The number of elements is 15. } \label{fig:deer}
\end{figure}

In computer vision, researchers often distinguish between edges and regions. Actually, these are two relative concepts in the frequency domain. While edges can be captured by high frequency Gabor wavelets, the regions can be encoded by low frequency wavelets, including the difference of Gaussian (DoG) wavelets. To account for both edges and regions, we need to combine Gabor and DoG wavelet elements at multiple frequency bands.

In Experiment 9, we select the wavelet elements of active basis from a dictionary of Gabor wavelets and Difference of Gaussian (DoG) wavelets at different scales. We use the same shared sketch algorithm with sigmoid pursuit index, except that we normalize the filter responses by marginal variance. After selecting the elements and record their responses, we use matching pursuit \cite{Mallat_Zhang}  to reconstruct the images. We need to use matching pursuit for reconstruction because the selected elements are only approximately orthogonal to each other, so the projection coefficients and the reconstruction coefficients are slightly different. The matching pursuit algorithm computes the reconstruction coefficients from the project coefficients.


\begin{figure}[!htb]
\begin{center}
\begin{tabular}{cccc}
\psfig{figure=ex10_1sym1.eps,height=.4in}&
\psfig{figure=ex10_1sym2.eps,height=.4in}&
\psfig{figure=ex10_1sym3.eps,height=.4in}&
\psfig{figure=ex10_1sym4.eps,height=.4in}
\end{tabular}
\end{center}
\caption{Experiment 9.1. The selected Gabor elements (illustrated by bars) at 3 different scales and the selected DoG elements (illustrated by circles, and larger circles are darker than smaller ones). The lengths of the Gabor elements are 35, 25, and 17 pixels respectively. The sizes of the DoG elements are 77 and 55 respectively. The allowed activity of location is 4 pixels for both Gabor and DoG elements.} \label{fig:S1}
\end{figure}

Figure (\ref{fig:S1}) displays the select Gabor and DoG elements. The Gabor elements are illustrated by bars at different sizes, and the DoG elements are illustrated by circles. The radius of a circle is about half of that of the blob represented by the corresponding DoG elements. Larger circles are darker than smaller ones.

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=ex10_1.eps,height=1.2in}
\end{tabular}
\end{center}
\caption{Experiment 9.1. The fist block displays all the 50 selected Gabor and DoG elements. The smaller Gabors are illustrated by darker bars. The remaining blocks display the original images and the corresponding reconstructed images. The image size is 102 $\times$ 100.} \label{fig:S2}
\end{figure}

Figure (\ref{fig:S2}) displays the reconstructed images. The DoG elements are necessary to account for the large regional contrasts.

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=ex10_2.eps,height=1.4in}
\end{tabular}
\end{center}
\caption{Experiment 9.2. The fist block displays all the 50 selected Gabor and DoG elements. The remaining blocks display the original 95 $\times$ 100 images and the corresponding reconstructed images.} \label{fig:S3}
\end{figure}

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=ex10_4.eps,height=1.2in}
\end{tabular}
\end{center}
\caption{Experiment 9.4. The fist block displays all the 40 selected Gabor and DoG elements. The remaining blocks display the 100 $\times$ 70 original images and the corresponding reconstructed images.} \label{fig:S4}
\end{figure}


\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=ex10_5.eps,height=1.in}
\end{tabular}
\end{center}
\caption{Experiment 9.5. The fist block displays all the 40 selected Gabor and DoG elements. The remaining blocks display the original 100 $\times$ 110 images and the corresponding reconstructed images.} \label{fig:S5}
\end{figure}



\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=ex10_7.eps,height=1.in}
\end{tabular}
\end{center}
\caption{Experiment 9.7. The fist block displays all the 50 selected Gabor and DoG elements. The remaining blocks display the original 100 $\times$ 100 images and the corresponding reconstructed images.} \label{fig:S6}
\end{figure}

Figures (\ref{fig:S3}) - (\ref{fig:S6}) display more examples. Ideally, the large Gabor elements gauge the breadths of the edges, while the small Gabor elements gauge the sharpness of the edges. The large DoG elements gauge the sizes of the regions, which are to be contoured by the Gabor elements. It is unclear whether we could use large Gabor and DoG elements to substitute the region concepts, so that we do not have to do image segmentation.


\section{Composing Multiple Part-Templates}

For articulate objects, we need to represent them as compositions of part-templates at different locations and resolutions.

\subsection{Recursive active basis and recursive sum-max maps}

An active basis is a composition of multiple Gabor wavelet elements,
where each element is allowed to shift its location and orientation.
We can further compose multiple active bases, where each active
basis serves as a part-template that is allowed to change its overall location, orientation and scale. We call such a structure a
``recursive active basis,'' which is a template that consists of multiple part-templates. The following experiment illustrates the basic idea.


\begin{figure}[!htb]
\begin{center}
\begin{tabular}{cc}
\psfig{figure=ex7_original.eps,height=1in}&
\psfig{figure=ex7_4.eps,height=1in}\\
(a) & (b)
\end{tabular}
\end{center}
\caption{Experiment 10.1. (a) Input image of 330 $\times$ 496. (b)
Superposed sketch. The bounding box of the front wheel is 112
$\times$ 126. The bounding box of the back wheel is 86 $\times$ 76.
The total number of elements is 60.  } \label{fig:ex7_1}
\end{figure}

Figure (\ref{fig:ex7_1}.a) displays an image of size 330 $\times$
496. Figure (\ref{fig:ex7_1}.b) displays the superposed sketch. The
template is learned in Experiment 1.2 from the bicycle images. See
Figure (\ref{fig:ex1_2}). We split the bicycle template in Figure
(\ref{fig:ex1_2}) horizontally into two part-templates. The bounding
box for the part-template of the front wheel is 112 $\times$ 126. The
bounding box for the part-template of the back wheel is 86 $\times$
76. We allow the two part-templates to locally shift horizontally, so
these two part-templates make up a multi-basis. We then fit this
multi-basis to the tandem bike in Figure (\ref{fig:ex7_1}.a) and
obtain the sketch in Figure (\ref{fig:ex7_1}.b).

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=SumMaxMapsS3.eps,height=5in}
\end{tabular}
\end{center}
\caption{Sum-max maps. A $\S2$ map is computed for each
part-template. For each $\S2$ map, a $\M2$ map is computed by
applying a local maximization operator to the $\S2$ map. Then a
$\S3$ map is computed by summing over the two $\M2$ maps. The $\S3$
map scores the template matching, where the template consists of two
part-templates that are allowed to locally shift their locations. }
\label{fig:SumMaxMapsS3}
\end{figure}

Given the two part-templates, the inference can be accomplished by
alternating the sum maps and max maps as illustrated by Figure
(\ref{fig:SumMaxMapsS3}). Here we have two $\S2$ maps, one for each part-template. On top of each $\S2$ map, there is also a $\M2$ maps. Then on top of the two $\M2$ maps, a $\S3$ maps is computed. After that a $\M3$ score is obtained. These scores are computed by a bottom-up process, and they answer the following questions:

$\S2$ maps: Is there a part-template at this location?

$\M2$ maps: Is there a part-template at a {\em nearby} location?

$\S3$ map: Is there a certain composition of part-templates that form the whole template at this location?

$\M3$ score: Is there a composite template within the whole image?

If there is such a composite template, then a top-down process retrieves the location of the whole template, and then retrieves the locations of the part-templates, and finally the elements of the part-templates.

\vspace{.5cm} \noindent{\bf Inference by recursive sum-max maps}

\begin{itemize}

\item[Up-1]
For $j = 1, 2$, compute $\S2(x, y, j)$ using the inference algorithm of Subsection (\ref{sect:recognition}).

\item[Up-2]
For all $(x, y)$ and $j = 1, 2$, compute
\begin{eqnarray}
\M2(x, y, j) = \max_{\tiny{\begin{array}{l}-b_{x} \leq \Delta x \leq
b_x\\-b_y \leq \Delta y \leq b_y\end{array}}} \S2(x+\Delta x, y +
\Delta y, j), j = 1, 2. \label{eq:M2}
\end{eqnarray}

\item[Up-3]
For all $(x, y)$ and $j = 1, 2$, compute
\begin{eqnarray}
\S3(x, y) = \sum_{j=1}^{2} \M2(x+x_j, y+y_j).
\end{eqnarray}

\item[Up-4]
Compute $\M3 = \max_{x, y} \S3(x, y)$.

\item[Down-4]
Retrieve $(\hat{x}, \hat{y})$ that achieves the maximum in Up-4.

\item[Down-3]
Retrieve $(\hat{x}+x_j, \hat{y}+y_j)$ in Up-3 for $j = 1, 2$.

\item[Down-2]
Retrieve $(\hat{x}_j, \hat{y}_j)$ so that
\begin{eqnarray}
\M2(\hat{x}+x_j, \hat{y}+y_j, j) = \S2(\hat{x}_j, \hat{y}_j, j),
\end{eqnarray}
in the local maximization operation (\ref{eq:M2}) in Up-2.

\item[Down-1] Retrieve the perturbed elements of $j$-th part-template as described the inference algorithm of Subsection (\ref{sect:recognition}).

\end{itemize}


In our experiment, in Step Up-2, we take $b_x = 20$ pixels and $b_y = 4$ pixels. Let $(x_1, y_1)$ and $(x_2, y_2)$ be the central
positions of the bounding boxes for the two part-templates in the
original template learned from regular bicycles. Assume $x_1 < x_2$,
we let $x_1 \leftarrow x_1 - b_x$, $x_2 \leftarrow x_2 + b_x$, and
let the two part-templates shift around the new centers $(x_1, y_1)$
and $(x_2, y_2)$.

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=ex7_Splitscore.eps,height=2in}
\end{tabular}
\end{center}
\caption{Experiment 10.1. $\M3$ scores for different splitting points.
} \label{fig:ex7_3}
\end{figure}

The $\M3$ score in Step 4 measures the template matching, or the
alignments of the two part-templates to the images. This $\M3$ score
can be used to decide where we should split the original bicycle
template. Specifically, we can try different splitting points, and
for each splitting point, we compute the $\M3$ score. Figure
(\ref{fig:ex7_3}) displays the $\M3$  scores for 10 different
splitting points. The result shown in Figure (\ref{fig:ex7_1}.b) is
obtained at the splitting point the achieves the maximum $\M3$
score.


\begin{figure}[!htb]
\begin{center}
\begin{tabular}{cc}
\psfig{figure=ex7_S2_1.eps,height=.7in}&
\psfig{figure=ex7_S2_2.eps,height=.7in}
\end{tabular}
\end{center}
\caption{Experiment 10.1. The two $\S2$ maps for the two part-templates
at the optimal splitting point. } \label{fig:ex7_S2}
\end{figure}

Figure (\ref{fig:ex7_S2}) displays the two $\S2$ maps for the two
part-templates at the optimal splitting point.

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{cc}
\psfig{figure=ex7_M2_1.eps,height=.7in}&
\psfig{figure=ex7_M2_2.eps,height=.7in}
\end{tabular}
\end{center}
\caption{Experiment 10.1. The two $\M2$ maps for the two part-templates
at the optimal splitting point.} \label{fig:ex7_M2}
\end{figure}

Figure (\ref{fig:ex7_M2}) displays the two $\M2$ maps for the two
part-templates at the optimal splitting point.

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=ex7_S3.eps,height=.7in}
\end{tabular}
\end{center}
\caption{$\S3$ map at the optimal splitting point.}
\label{fig:ex7_M3}
\end{figure}

Figure (\ref{fig:ex7_M3}) displays the $\S3$ map at the optimal
splitting point.



\begin{figure}[!htb]
\begin{center}
\begin{tabular}{ccccc}
\psfig{figure=ex7_1.eps,height=.6in}&
\psfig{figure=ex7_2.eps,height=.6in}&
\psfig{figure=ex7_3.eps,height=.6in}&
\psfig{figure=ex7_4.eps,height=.6in}&
\psfig{figure=ex7_5.eps,height=.6in}\\
\psfig{figure=ex7_6.eps,height=.6in}&
\psfig{figure=ex7_7.eps,height=.6in}&
\psfig{figure=ex7_8.eps,height=.6in}&
\psfig{figure=ex7_9.eps,height=.6in}&
\psfig{figure=ex7_10.eps,height=.6in}
\end{tabular}
\end{center}
\caption{Experiment 10.1. Superposed sketches obtained by all the
splitting points. } \label{fig:ex7_split}
\end{figure}

Figure (\ref{fig:ex7_split}) displays superposed sketches obtained
at all the splitting points.

In this experiment, we set the response of a basis element to 0 if it is outside the boundary of the image. So the $\S2$ maps are of the same size as the original image. If the center of a part-template is outside the image boundary, we set the responses of all its basis elements to 0. So the $\S3$ map is also of the same size as the original image.

The recursive active basis can be considered a constellation model \cite{Weber} whose constituent components are active basis. The $\M2$ and $\S2$ maps may have been commonly used in part-based models. Thanks to the work of Riesenhuber and Poggio \cite{Riesenhuber_Poggio}, we are able to extend the SUM and MAX operations down to the image intensities.

\subsection{Account for large deformations}

The recursive active basis and recursive sum-max maps can account for the existence of parts, as illustrated in Experiment 10.1. They can also be used to deal with large deformations.

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{ccc}
\psfig{figure=ex7_horse_original.eps,height=.8in}&
\psfig{figure=ex7_horse_sketch.eps,height=.8in}&
\psfig{figure=ex7_horse_sketch0.eps,height=.8in}\\
(a) & (b) & (c)
\end{tabular}
\end{center}
\caption{Experiment 10.2. (a) The observed image. (b) Superposed with sketch where the horse template is split horizontally into two part-templates. These two part-templates are allowed to move horizontally up to 10 pixels in each direction. (c) Superposed with sketch using the original horse template. In other words, the two part-templates are not allowed to move relative to each other. } \label{fig:M1}
\end{figure}

Figure (\ref{fig:M1}.a) displays the image of horse that we used in Experiment 4, where we change the aspect ratio of the horse template to fit this image. From a 2D point of view, this amounts to a large deformation that cannot be handled by a single-layer active basis model. However, we still can use the same method in Experiment 10.1 to split the original horse template into two part-templates, and allow these two part-templates to move relative to each other. Figure (\ref{fig:M1}.b) displays the result of fitting the recursive active basis at the optimal splitting point. As a comparison, Figure (\ref{fig:M1}.c) displays the result using the original template. The original template does not fit the rear part of the horse very well.

{\em Correlated activities}. At this point, we would like to discuss the modeling of the perturbations or activities of the basis elements in the active basis model. For simplicity, we assume that the perturbations are independent and uniformly distributed within a small range. This leads to the simple inference algorithm based on the sum-max maps. The computational complexity is linear in the size of the image.

The recursive active basis can be considered a scheme to account for correlated activities, where those elements that belong to the same part-template share a common overall movement, in addition to their perturbations relative to the overall location of the part-template.

The method of Experiment 5 may be used to learn part-templates. We shall report our work on this issue elsewhere.

\section{Discussion}

We conclude this article with some comments on our method and a discussion of several related principles.

\subsection{Simplicity} We argue that our method cannot be simpler for the vision tasks that we have studied in this article. The model is no much more complex than a wavelet expansion, except that we add local perturbations to the wavelet elements. The learning algorithm is no much more complex than edge detection, except that we perform it in a parallel fashion on multiple images. The inference algorithm only involves two consecutive operations on top of Gabor filtering. One is a local max filtering and the other is a local sum filtering.

For unsupervised learning, such as learning from non-aligned images and local learning of representative templates, we exploit the fact that our method can learn from single images,  and the iterative learning procedures only involve an extra step of detection or classification in addition to the step of supervised learning. It remains to be seen how reliable it is to initialize from single image learning.

The learning and inference algorithms are very easy to program. There are not many sensitive tuning parameters, although we did tuned a couple of parameters such as the image resize factor and the number of elements. The data and codes for reproducing the experimental results in this paper are available from the reproducibility page. We only include some of the positive results in this paper. To avoid partial reporting, we did summarize our negative experiences in the experiments, although we did not document them. Because of the simplicity of the method, when things go wrong, it is often fairly easy to see where is the problem.

Our work is a simple logical consequence if we assume the validity of Olshausen-Field theory for V1 simple cells \cite{Olshausen_Field}, and the Riesenhuber-Poggio theory for V1 complex cells \cite{Riesenhuber_Poggio}. We, of course, are not in a position to claim any neuroscientific validity of the $\S2$ maps in our method. Our work is also a simple re-tooling of the texton model of Zhu, et al. \cite{Zhu_texton}.

In retrospect, we find that the following three principles provide useful insights.

\subsection{Sparsity} Olshausen and Field \cite{Olshausen_Field} proposed this principle for understanding V1 simple cells, where a typical natural image can be represented by a linear superposition of a small number of Gabor-like wavelet elements at different scales, locations, and orientations, plus a small residual image. The reason for such a sparse representation is that edges are prominent and frequently occurring structures in natural images.

The active basis model can be considered a further sparsification. The reason is as follows. In Olshausen-Field representation, each image is encoded by a set of locations, orientations, and scales of wavelet elements. We can further encode these locations, orientations, and scales by an even smaller number of templates, each being a composition of locations, orientations, and scales. The reason for such a sparser representation is that those templates are prominent and frequently occurring structures in natural images. In Olshausen-Field representation, we need to allow for small residuals in image intensities. Similarly, in this template representation, we need to allow for small residuals in locations, orientations, and scales. Such small residuals become the perturbations in the active basis model, so that the templates are deformable.

\subsection{Compositionality} S. Geman et al.  \cite{Geman_composition} proposed this principle for vision. If we want a compositional representation of image intensities and if we insist on linear representation for simplicity, then it is natural to adopt wavelet representation because the wavelet elements are localized in both spatial and frequency domains. The active basis model follows such a compositional scheme.

Zhu and Mumford \cite{Zhu_grammar} investigated the and-or graph as a recursive compositional scheme for vision, where ``and'' accounts for compositions of constituent elements, while ``or'' accounts for variations in the constituent elements. The active basis model is a
simplest form of an and-or graph, where ``and'' means composition of
wavelet elements, and ``or'' means variations in the locations and
orientations of the elements. The and-or graph is a grammar that can be applied  recursively. The recursive active basis follows such a grammar.

The recursive sum-max maps is a variation of the cortex-like structure proposed by Riesenhuber and Poggio \cite{Riesenhuber_Poggio}. It is a natural hierarchical structure for parsing an image according to the and-or grammar. The sum maps score the and-compositions, and the max maps account for the or-variations. After bottom-up scoring for detection and classification, the top-down retrieving produces the parsing of the image. See also the recent work of L. Zhu et al. \cite{Leo} on a recursive compositional scheme.

\subsection{Invariance} Riesenhuber and Poggio \cite{Riesenhuber_Poggio} proposed this principle for V1 complex cells. While the V1 simple cells capture the essence of the image intensities via Olshausen-Field sparse coding,  the local maximization operation of the V1 complex cells filters out shape deformations, and makes the subsequent  processing invariant to shape deformations. Of course, invariance here is only approximate.

The Riesenhuber-Poggio scheme compares intensities of the $\M1$ maps directly for template matching. We modified their template matching scheme by a weighted sum of the $\M1$ intensities at highly selected locations and orientations. If the locations and orientations of the selected wavelet elements are at the centers of the local perturbations that cause shape deformation, then hopefully, the intensities of the $\M1$ maps of these highly selected locations and orientations are more invariant (and more indicative of the object shapes) than the intensities of other locations and orientations.



\section*{Acknowledgement}
We thank Chuck Fleming for earlier collaboration on some of the
experiments. We thank Alan Yuille, Stefano Soatto, Zhuowen Tu, and Leo Zhu for discussions.  The work is supported by NSF-DMS 0707055, NSF-IIS 0713652, ONR N00014-05-01-0543, and Keck foundation. We acknowledge the use of data sets \cite{Yao_DB} provided by the Lotus Hill Institute, which is supported by a Chinese National 863 grant 2006AA01Z121 and an NSFC grant 60728203.


{\small
\begin{thebibliography}{9}


\bibitem{Cootes}
T. F. Cootes, G. J. Edwards, and C. J. Taylor, ``Active appearance
models,'' {\em IEEE Transactions on Pattern Analysis and Machine
Intelligence}, 23, 681-685, 2001.

\bibitem{SVM}
C. Cortes and V. Vapnik, ``Support vector networks,'' {\em Machine Learning}, 20, 273-297, 1995.

\bibitem{Daugman}
J. Daugman, ``Uncertainty relation for resolution in space, spatial
frequency, and orientation optimized by two-dimensional visual
cortical filters,'' {\em Journal of Optical Society of America}, 2,
1160-1169, 1985.

\bibitem{DellaPietra}
S. Della Pietra,   V. Della Pietra,  and J. Lafferty, ``Inducing
features of random fields,'' {\em IEEE Transactions on Pattern
Analysis and Machine Intelligence}, 19, 380-393, 1997.

\bibitem{FisEl}
M. Fischler and R. Elschlager, ``The representation and matching of
pictorial structures,'' {\it IEEE Transactions on Computers}, C-22,
67-92, 1973.

\bibitem{adaboost}
Y. Freund and  R. E. Schapire, ``A decision-theoretic generalization
of on-line learning and an application to boosting,'' {\em Journal
of Computer and System Sciences}, 55, 119-139, 1997.

\bibitem{Friedman}
J. H.  Friedman, ``Exploratory projection pursuit,''  {\em Journal
of the American Statistical Association}, 82, 249-266, 1987.

\bibitem{Geman_composition}
S. Geman, D. F. Potter, and Z. Chi, ``Composition system,'' {\em
Quarterly of Applied Math}, 60, 707-736, 2002.

\bibitem{GuoZhuWu}
C. Guo, S. C. Zhu, and Y. N. Wu, ``Towards a mathematical theory of
primal sketch and sketchability,'' {\em International Conference on
Computer Vision}, 2, 1228-1235, 2003.

\bibitem{GuoZhuWu2}
C. Guo, S. C. Zhu, and Y. N. Wu, ``Primal sketch: integrating structure and texture,'' {\em Computer Vision and Image Understanding}, 106, 5-19, 2007.


\bibitem{Kass_snakes}
M. Kass, A. Witkin, and D. Terzopoulos, ``Snakes: Active contour
models,'' {\em International Journal of Computer Vision}, 1,
321-331, 1988.

\bibitem{Mallat_Zhang}
S. Mallat and Z. Zhang, ``Matching pursuit in a time-frequency
dictionary,'' {\em IEEE Transactions on Signal Processing}, 41,
3397-415, 1993.

\bibitem{Lowe}
J. Mutch and D. G. Lowe, ``Multiclass object recognition with
sparse, localized features,'' {\em Proceedings of Computer Vision
and Pattern Recognition}, 2006.


\bibitem{Olshausen_Field}
B. A. Olshausen and D. J. Field, ``Emergence of simple-cell
receptive field properties by learning a sparse code for natural
images,'' {\em Nature}, 381, 607-609, 1996.


\bibitem{Riesenhuber_Poggio}
M. Riesenhuber and T. Poggio, ``Hierarchical models of object
recognition in cortex,'' {\em Nature Neuroscience},  2, 1019-1025,
1999.

\bibitem{Tu}
Z. Tu, ``Learning generative models via discriminative approaches,'' {\em Proceedings of IEEE Computer Vision and Pattern Recognition}, 2007.

\bibitem{Viola_Johnes}
P. A. Viola and M. J. Jones, ``Robust real-time face detection,''
{\em International Journal of Computer Vision}, 57, 137-154, 2004.


\bibitem{Weber}
M. Weber, M. Welling and P. Perona, ``Towards automatic discovery of
object categories'', {\em Proceedings of Computer Vision and Pattern
Recognition}, 2000.

\bibitem{Witkin}
A. Witkin, ``Scale-space filtering,'' {\em Proceedings of
International Joint Conference on Artificial Intelligence}, 1983.

\bibitem{Wu}
Y. N. Wu, C. Guo, S. C. Zhu,  ``From information scaling of natural
images to regimes of statistical models,'' {\em Quarterly of Applied
Math}, 66, 81-122, 2008.

\bibitem{Wu_ICCV07}
Y. N. Wu, Z. Si, C. Flemming, and S. C. Zhu, ``Active basis as
deformable templates,'' {\em Proceedings of International Conference
on Computer Vision}, 2007.

\bibitem{Wu_texture_sketch}
Y. N. Wu, S. C. Zhu, and C. Guo,  ``Statistical modeling of texture
sketch,'' \emph{Proceedings of European Conference of Computer
Vision,} 2002.

\bibitem{Wu_equivalence}
Y. N.  Wu, S.C. Zhu, and X.  Liu, ``Equivalence of Julesz ensemble
and FRAME models,'' \emph{International Journal of Computer Vision},
38, 245-261, 2000.


\bibitem{Yao_DB}
Z. Yao, X. Yang, and S. C. Zhu, ``Introduction to a large scale
general purpose groundtruth database: methodology, annotation tools,
and benchmarks,''  {\em 6th International Conference on EMMCVPR},
2007.


\bibitem{Yuille}
A. L. Yuille, P. W. Hallinan, and D. S. Cohen, ``Feature extraction
from faces using deformable templates,'' {\em International Journal
of Computer Vision}, 8, 99-111, 1992.

\bibitem{Leo}
L. Zhu, C. Lin, H. Huang, Y. Chen, and A. Yuille, ``Unsupervised structure learning: hierarchical recursive composition, suspicious coincidence and competitive exclusion,'' \emph{Proceedings of European Conference of Computer Vision,} 2008.

\bibitem{Zhu_texton}
S. C. Zhu, C. E. Guo, Y. Z. Wang, and Z. J. Xu, ``What are
textons,'' {\em International Journal of Computer Vision}, 62,
121-143, 2005.

\bibitem{Zhu_prior}
S. C. Zhu and D. B. Mumford, ``Prior learning and Gibbs
reaction-diffusion,'' {\em IEEE Transactions on Pattern Analysis and
Machine Intelligence}, 19, 1236-1250, 1997.

\bibitem{Zhu_grammar}
S. C. Zhu and D. B. Mumford, ``A stochastic grammar of images,''
{\em Foundations and Trends in Computer Graphics and Vision}, 2,
259-362, 2006.


\bibitem{Zhu_FRAME}
S. C. Zhu, Y. N. Wu, and D. B. Mumford, ``Minimax entropy principle
and its applications in texture modeling,'' {\em Neural
Computation}, 9, 1627-1660, 1997.

\end{thebibliography}
}

\end{document}
