\documentclass[10pt,twocolumn,letterpaper]{article}

\usepackage{cvpr}
\usepackage{times}
\usepackage{epsfig}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{amssymb}

\def\I{{\bf I}}
\def\B{{\bf B}}
\def\R{{\bf R}}
\def\corr{{\rm corr}}
\def\E{{\rm E}}

% Include other packages here, before hyperref.

% If you comment hyperref and then uncomment it, you should delete
% egpaper.aux before re-running latex.  (Or just hit 'q' on the first latex
% run, let it finish, and you should be clear).
\usepackage[pagebackref=true,breaklinks=true,letterpaper=true,colorlinks,bookmarks=false]{hyperref}


% \cvprfinalcopy % *** Uncomment this line for the final submission

\def\cvprPaperID{556} % *** Enter the CVPR Paper ID here
\def\httilde{\mbox{\tt\raisebox{-.5ex}{\symbol{126}}}}

% Pages are numbered in submission mode, and unnumbered in camera-ready
\ifcvprfinal\pagestyle{empty}\fi
\begin{document}

%%%%%%%%% TITLE
\title{Learning Codebooks of Deformable Templates for Sparse Image
Representation}

\author{First Author\\
Institution1\\
Institution1 address\\
{\tt\small firstauthor@i1.org}
% For a paper whose authors are all at the same institution,
% omit the following lines up until the closing ``}''.
% Additional authors and addresses can be added with ``\and'',
% just like the second author.
% To save space, use either the email address or home page, not both
\and
Second Author\\
Institution2\\
First line of institution2 address\\
{\small\url{http://www.author.org/~second}} }

\maketitle
% \thispagestyle{empty}

%%%%%%%%% ABSTRACT
\begin{abstract}
This paper proposes a method for unsupervised learning of codebooks
of deformable templates for sparse representations of images of
various objects and textures. Each deformable template is represented
by an active basis model, which is a composition of a small number of
Gabor wavelets automatically selected from a dictionary of such
wavelets. The selected wavelets are allowed to perturb their
locations and orientations so that the active basis template can
deform to encode the observed images. For a given set of training
images, our method is to learn a codebook of such active basis
templates, so that each training image can be represented by a small
number of templates automatically selected from the codebook. The
learning algorithm iterates the following two steps. (1) Image
encoding by a template matching pursuit algorithm. (2) Codebook
re-learning by a shared matching pursuit algorithm. Our experiments
show that the method is capable of learning meaningful codebooks from
images of textures and objects, and the codebooks give meaningful
representations of these images. Our experiments also show that the
learned codebooks can be used as features for classification.
\end{abstract}

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

Learning codebooks or dictionaries of representational elements for
images of natural scenes or images of various types of textures and
objects is one of the most fundamental problems in both computer and
biological vision. This paper proposes a method for unsupervised
learning of codebooks of deformable templates that give sparse
representations of images.

\subsection{Explicit codebook and representation}

Our method is based on an explicit generative model of image
intensities, where each image is represented by a small number of
templates automatically selected from a codebook of templates. Each
template in the codebook is further encoded by a small number of
Gabor wavelet elements automatically selected from a dictionary of
Gabor wavelets. The codebook is learned from training images without
any labelling or annotation. We also require that the number of
templates in the codebook to be as small as possible.

\begin{figure}[htb]
 \includegraphics[height=.045\textwidth]{img/template_sorted_pave2.png}\\
 \includegraphics[height=.17\textwidth]{img/pave2.jpg}
\includegraphics[height=.17\textwidth]{img/sketch_image1_pave2.png}\\
\caption{Learning a codebook of 3 active basis templates from an image of pavement. {\em Each black bar represents a Gabor wavelet at the
same location, orientation and length}. The size of the image is 578 by 434 pixels. Each active basis template (i.e. entry in the codebook) is of size 60 by 60 pixels, with no more than 10 selected wavelets. The input image is paved by 43 deformed templates selected from the codebook, with their bounding boxes shown in green color.}
\label{fig:pave}
\end{figure}


Figure \ref{fig:pave} illustrates the basic idea. The first row displays a codebook of 3 templates learned from a
pavement image on the left of the second row. Each template in the codebook is represented by an active basis
model. Specifically, each template is a composition of a small number
of Gabor wavelets automatically selected from a dictionary of
Gabor wavelets. {\em Each selected Gabor wavelet is depicted by a bar
that has the same location, orientation and length as the Gabor wavelet}.
The selected Gabor wavelets are allowed to perturb their locations
and orientations, so the template becomes deformable, and the
linear basis formed by the selected Gabor wavelet elements becomes
active. The sketch on the right of the second row displays the representation of the original pavement image using the
templates in the codebook. Here the word ``pavement'' also
serves as a good metaphor because the representation is essentially a ``pavement'' made of the active basis
templates. Specifically, each of the
templates in the sketch is a {\em spatially translated, scaled, rotated and
deformed version} of a template in the codebook. The
bounding box of each template in the sketch is also shown in green color (the deformed
templates are allowed to have some overlaps with each
other). As shown in the caption, the image can be paved by a relatively small
number of templates, and such a representation is said to be sparse.
It is sparser than wavelet representation because the coding
elements are compositions or groups of wavelets.

\begin{figure}[htb]
\includegraphics[height=.07\textwidth]{img/template_sorted_cat.png}\\
\includegraphics[height=.05\textwidth]{img/cat4.jpg}
\includegraphics[height=.05\textwidth]{img/sketch_image4.png} \hspace{1mm}
\includegraphics[height=.05\textwidth]{img/cat12.jpg}
\includegraphics[height=.05\textwidth]{img/sketch_image12.png} \hspace{1mm}
\includegraphics[height=.05\textwidth]{img/cat24.jpg}
\includegraphics[height=.05\textwidth]{img/sketch_image24.png} \hspace{1mm}
\includegraphics[height=.05\textwidth]{img/cat25.jpg}
\includegraphics[height=.05\textwidth]{img/sketch_image25.png}
\caption{Learning a codebook of 20 active basis templates from 25 images of cats (not assumed to be aligned). The input images are resized to have roughly 22,500 pixels. Each active basis template (i.e. entry in the codebook) is of size 72 by 72 pixels, with no more than 12 selected wavelets. On average about 8 templates are used to sketch each image. }
\label{fig:cat}
\end{figure}

Figure \ref{fig:cat} shows another example, where a codebook of active basis
templates is learned from a set of training images of cats.

\subsection{Unsupervised learning algorithm}

For a given set of training images without any labeling or
annotation, our learning algorithm iterates the following two steps:

{\em Step 1: Image encoding by template matching pursuit:} Given the
current codebook of templates, for each training image, a template
matching pursuit algorithm sequentially selects the templates from
the codebook to pave the whole image.

{\em Step 2: Codebook re-learning by shared matching pursuit:} Given
the current encoding, for each template, a shared matching pursuit
algorithm re-learns the template from all the image patches that are
currently covered by this template. The algorithm sequentially
selects the constituent Gabor wavelet elements of the template from
the dictionary of Gabor wavelets, in order to encode all
these patches simultaneously.

Our experiments show that our method is capable of learning
meaningful codebooks from images of textures and objects, and the
codebooks give meaningful representations of these images. Our
experiments also show that the learned codebooks can be used as
features for classification.

\subsection{Past work}

Our generative model is built on the wavelet sparse coding model of
Olshausen and Field \cite{Olshausen_Field} and the active basis model
of Wu et al. \cite{Wu_active_basis_IJCV}. In Olshausen-Field model,
an image is represented by a small number of wavelet elements
selected from a dictionary. In our model, an image is
represented by a smaller number of groups of wavelet elements, and
these groups of wavelets exhibit recurring compositional patterns
that can be represented by active basis models.

The part of our work on textons can be considered a continuation of
the work Zhu et al. \cite{Zhu_texton}. The learned codebooks can
serve as Or-nodes in the And-Or graph studied by Zhu and Mumford
\cite{Zhu_grammar}. The learning method is also related to L. Zhu et
al. \cite{LeoZhu}.

\section{Olshausen-Field model}

Olshausen and Field \cite{Olshausen_Field} propose that the role of
simple V1 cells is to provide sparse representations of natural
images. This section reviews the Olshausen-Field model in order to
set the stage and fix the notation.

\subsection{Overcompleteness and sparsity} \label{sect:OF}

Let $\{\I_m, m = 1, ..., M\}$ be a set of training image patches
(e.g. $12 \times 12$), Olshausen-Field model seeks to represent these
images by
\begin{eqnarray}
    \I_m = \sum_{i=1}^{N} c_{m, i} B_i + U_m, \label{eq:OF1}
\end{eqnarray}
where $(B_i, i = 1, ..., N)$ is a dictionary or codebook of basis
elements of the same dimensionality as $\I_m$, $c_{m, i}$ are the
coefficients, and $U_m$ is the unexplained residual image. $N$ is
often assumed to be greater than the dimensionality of $\I_m$ (e.g.
$N = 2 \times 12 \times 12$), so the dictionary is said to be
overcomplete. On the other hand, the number of coefficients $(c_{m,
i}, i = 1, ..., N)$ that are non-zero or significantly different from
zero is assumed to be very small for each image $\I_m$. The
dictionary of $(B_i, i = 1, ..., N)$ can be learned automatically from the training
images $\{\I_m\}$ by imposing a sparsity regularization function.

\subsection{Self-similarity and geometric attributes}

One may also assume that the dictionary of the basis elements are
translated, rotated and dilated version of one another, as in
Olshausen et al. \cite{Olshausen_wavelet}, so that each $B_i$ can be
written as $B_{x, s, \alpha}$, where $x$ is the location (a
two-dimensional vector), $s$ is the scale, and $\alpha$ is the
orientation. We call such a dictionary {self-similar}, and we call
$(x, s, \alpha)$ the geometric attribute of $B_{x, s, \alpha}$.

Model (\ref{eq:OF1}) then becomes
\begin{eqnarray}
\I_m = \sum_{x, s, \alpha} c_{m, x, s, \alpha} B_{x, s, \alpha} +
U_m, \label{eq:OF2}
\end{eqnarray}
where $B_{x, s, \alpha}$ are translated, rotated and dilated copies
of a {\em single} basis element, e.g. $B = B_{x=0, s=1, \alpha =0}$,
and $(x, s, \alpha)$ are properly discretized. $B$ can be learned
from training images $\{\I_m\}$.

{\em Assumption on wavelets in this paper}. From now on, we assume
that the dictionary of wavelets is self-similar, and $(B_{x, s,
\alpha}, \forall (x, s, \alpha))$ is already available. It can either
be learned or designed. In the following, we assume that $B_{x, s,
\alpha}$ is a Gabor wavelet, and we also assume that $B_{x, s,
\alpha}$ is normalized to have unit $L_2$ norm so that $|B_{x, y,
\alpha}|^2 = 1$. $B_{x, s, \alpha}$ may also be a pair of Gabor sine
and cosine wavelets, so that for each Gabor wavelet $B$, $B = (B_0,
B_1)$. The corresponding coefficient $c = (c_0, c_1)$, and $cB = c_0
B_0 + c_1 B_1$. For projection $\langle \I, B\rangle = (\langle \I,
B_0\rangle, \langle \I, B_1\rangle)$, and $|\langle \I, B\rangle|^2 =
\langle \I, B_0\rangle^2 + \langle \I, B_1\rangle^2$.

\subsection{spatial point process}

Given the dictionary $(B_{x, s, \alpha}, \forall (x, s, \alpha))$,
the encoding of an image $\I_m$ amounts to inferring $(c_{m, x, s,
\alpha}, \forall (x, s, \alpha))$ in (\ref{eq:OF2}) under the
sparsity constraint, which means that only a small number of $(c_{m,
x, s, \alpha})$ are non-zero. That is, we seek to encode $\I_m$ by
\begin{eqnarray}
\I_m = \sum_{i=1}^{n} c_{m, i} B_{x_{m,i}, s_{m,i}, \alpha_{m,i}} +
U_m, \label{eq:OF3}
\end{eqnarray}
where $n \ll N$ is a small number, and $(x_{m, i}, s_{m, i}, \alpha_{m, i},
i = 1, ..., n)$ are the geometric attributes of the selected wavelet
elements whose coefficients $(c_{m, i})$ are non-zero. $(x_{m, i},
s_{m, i}, \alpha_{m, i}, i = 1, ..., n)$ form a spatial point
process. (we continue to use $i$ to index the basis elements, but
here $i$ only runs through the $n$ selected basis elements instead of
all the $N$ basis elements as in (\ref{eq:OF1})).

\section{Active basis model}

The active basis model was proposed by Wu et al.
\cite{Wu_active_basis_IJCV} for modeling deformable compositional
patterns of the selected wavelet elements.

\subsection{A single template for aligned image patches}
Suppose we have a set of training image patches $\{\I_m, m = 1, ...,
M\}$. This time they are defined on the same bounding box. The
objects in these images come from the same category. They appear at
the same location, scale and orientation, and in the same pose within
the bounding box. The active basis model is of the following form
\begin{eqnarray}
\I_m = \sum_{i=1}^{n} c_{m, i} B_{x_i + \Delta x_{m, i}, s, \alpha_i
+ \Delta \alpha_{m, i}} + U_m,
\end{eqnarray}
where $\B = (B_{x_i, s, \alpha_i}, i = 1, ..., n)$ form the original
template. Here we assume that the scale $s$ is fixed and given. $\B_m
= (B_{x_i + \Delta x_{m, i}, s, \alpha_i + \Delta \alpha_{m, i}}, i =
1, ..., n)$ is the deformed template for encoding $\I_m$, where
$(\Delta x_{m, i}, \Delta \alpha_{m, i})$ are the perturbations in
the location and orientation respectively, in order to account for
deformation. Both $\Delta x_{m, i}$ and  $\Delta \alpha_{m, i}$ are
assumed to vary within limited ranges (e.g. $\Delta x_{m, i} \in [-3,
3]$ pixels, and $\Delta \alpha_{m, i} \in [-\pi/16, \pi/16]$).

\subsection{Shared matching pursuit: local maximum pooling and
arg-max explaining-away inhibition} \label{sect:SMP0}

The selection of $B_{x_i, s, \alpha_i}$ and the inference of its
perturbed versions $B_{x_i + \Delta x_{m, i}, s, \alpha_i + \Delta
\alpha_{m, i}}$ can be accomplished by a shared matching pursuit
algorithm that greedily minimizes $\sum_{m=1}^{M} |\I_m -
\sum_{i=1}^{n} c_{m, i} B_{x_i + \Delta x_{m, i}, s, \alpha_i +
\Delta \alpha_{m, i}}|^2$ (recall that the wavelet elements are
normalized to have unit $L_2$ norm).

[0] Initialize $i \leftarrow 0$. For $m = 1, ..., M$, initialize the
residual image $U_m \leftarrow \I_m$.

[1] $i \leftarrow i + 1$. Select the next element by $(x_i, \alpha_i)
= \arg\max_{x, \alpha} \sum_{m=1}^{M} \max_{\Delta x, \Delta \alpha}
|\langle U_m, B_{x+\Delta x, s, \alpha+\Delta \alpha}\rangle|^2$,
where $\max_{\Delta x, \Delta \alpha}$ is a local maximum pooling
within the small ranges of $\Delta x_{m, i}$ and  $\Delta \alpha_{m,
i}$ .

[2] For $m = 1, ..., M$, given $(x_i, \alpha_i)$, infer the
perturbations in location and orientation by retrieving the arg-max
in the local maximum pooling of step [1]: $(\Delta x_{m, i}, \Delta
\alpha_{m, i}) = \arg\max_{\Delta x, \Delta \alpha} |\langle U_m,
B_{x_i+\Delta x, s, \alpha_i+\Delta \alpha}\rangle|^2. $ Let $c_{m,
i} \leftarrow \langle U_m, B_{x_i+\Delta x_{m, i}, s, \alpha_i+\Delta
\alpha_{m, i}}\rangle$, and update the residual image $U_m
\leftarrow$ $U_m - c_{m, i} B_{x_i+\Delta x_{m, i}, s,
\alpha_i+\Delta \alpha_{m, i}}$.

[3] Stop if $i = n$, else go back to step [1].

The above algorithm is a generalization of the matching pursuit
algorithm of Mallat and Zhang \cite{Mallat_Zhang}. In step [1], the $\max_{\Delta x, \Delta \alpha}$ is the local
maximum pooling, proposed by Riesenhuber and Poggio
\cite{Riesenhuber_Poggio} as a function of V1 complex cells. The
selected $B_{x_i, s, \alpha_i}$ is supposed to encode all the images
$\{\I_m\}$ simultaneously, subject to local perturbations. That is,
$B_{x_i, s, \alpha_i}$ is shared by all the images. After the
selection of the shared wavelet element $B_{x_i, s, \alpha_i}$, we
infer its perturbations by retrieving the arg-max of the local
maximum pooling. This arg-max wavelet element $B_{x_i+\Delta x_{m,
i}, s, \alpha_i+\Delta \alpha_{m, i}}$ then explains away a small
part from $U_m$, thereby implicitly inhibits nearby wavelet elements
from being selected in the future.

{\em Assumption on orthogonality in this paper}. Because of the arg-max
explaining-away inhibition, the wavelet elements in each deformed
template $\B_m = (B_{x_i+\Delta x_{m, i}, s, \alpha_i+\Delta
\alpha_{m, i}}, i = 1, ..., n)$ usually have little overlaps with
each other. So from now on, we shall assume that these wavelet
elements are orthogonal to each other, so that the coefficient $c_{m,
i} = \langle \I_m, B_{x_i+\Delta x_{m, i}, s, \alpha_i+\Delta
\alpha_{m, i}}\rangle$. We write $C_m = (c_{m, i}, i = 1, ..., n)$.
In practice, we allow small overlaps between the elements of $\B_m$.

\subsection{Statistical modeling: foreground pops out from natural
image background}

The above algorithm implicitly assumes that the residual $U_m$ is
Gausian white noise. This assumption is wrong. A better assumption is
to assume that $U_m$ follows the same distribution as that of natural
images. Specifically, the distribution of $\I_m$ given the deformed
template $\B_m = (B_{x_i + \Delta x_{m, i}, s, \alpha_i + \Delta
\alpha_{m, i}}, i = 1, ..., n)$, i.e., $p(\I_m \mid \B_m)$, is
obtained by modifying the distribution of natural images $q(\I_m)$ in
such a way that we only change the distribution of 
$C_m = (c_{m, i}, i = 1, ..., n)$ from $q(C_m)$ to $p(C_m)$, while
leaving the conditional distribution of $U_m$ given $C_m$ unchanged.
Here $p(C_m)$ and $q(C_m)$ are the distributions of $C_m$ under
$p(\I_m \mid \B_m)$ and $q(\I_m)$ respective. We should find $\B$ and
perturb it to $\{\B_m\}$ so that $p(C_m)$ and $q(C_m)$ have the
maximum contrast in terms of the Kullback-Leibler divergence. Thus
the model is in the form of foreground $p(C_m)$ popping out from
background $q(\I_m)$. Specifically, $p(\I_m \mid \B_m) = q(\I_m)
p(C_m)/q(C_m)$. Such a density substitution scheme was first used in
projection pursuit density estimation
\cite{Friedman_projection_pursuit}.

For computational simplicity, we further assume that $(c_{m, i}, i =
1, ..., n)$ are independent given $\B_m$, under both $p$ and $q$, so
we have
\[
p(\I_m \mid \B_m) = q(\I_m) \prod_{i=1}^{n} \frac{p_i(c_{m,
i})}{q(c_{m, i})},
\]
where $q(c)$ is assumed to be the same for $i = 1, ..., n$ because
$q(\I_m)$ is stationary. $q(c)$ can be pooled from natural images in
the form of a heavy-tailed histogram of Gabor filter responses.

For parametric modeling, we assume the following exponential family
model,
\begin{eqnarray}
p_i(c) = \frac{1}{Z(\lambda_i)} \exp\{\lambda_i h(|c|^2) \} q(c),
\label{eq:EXP}
\end{eqnarray}
where $h(r)$ is a function of the response $r = |c|^2$ that saturates
for large $r$. Specifically, $h(r) = \xi [{2}/{(1 + e^{-2r/\xi})} -
1]$. $h(r)$ behaves like $h(r) \approx r$ for small $r$, but $h(r)
\rightarrow \xi$ (e.g., $\xi = 6$) as $r \rightarrow \infty$. The
reason we assume this saturation function is that in the natural
image background, there are also occasional (albeit less frequent)
edges that give equally large responses $r$ as those in foreground,
so the probability ratio $p_i(c)/q(c)$ should go to a positive
constant instead of 0 for large $|c|^2$. $Z(\lambda) = \int
\exp\{\lambda h(r)\}q(c)dc = \E_{q}[\exp\{\lambda h(r)\}]$ is the
normalizing constant. $\mu(\lambda) = \E_{\lambda}[h(r)]$ is the mean
parameter.

\subsection{Shared matching
pursuit revised: saturation and hard inhibition} \label{sect:SMP1}

We can revise the shared matching pursuit in section
(\ref{sect:SMP0}) in order to maximize the log-likelihood instead of
minimize the squared loss as in section (\ref{sect:SMP0}). The
algorithm is as follows.

[0] Initialize $i \leftarrow 0$. For $m = 1, ..., M$, initialize the
response maps $R_m(x, \alpha) \leftarrow \langle \I_m, B_{x, s,
\alpha}\rangle$ for all $(x, \alpha)$.

[1] $i \leftarrow i + 1$. Select the next wavelet element by finding
\[
(x_i, \alpha_i) = \arg\max_{x, \alpha} \sum_{m=1}^{M} \max_{\Delta x,
\Delta \alpha} h(|R_m(x+\Delta x, \alpha+\Delta \alpha)|^2),
\]
where $\max_{\Delta x, \Delta \alpha}$ is again local maximum
pooling.

[2] For $m = 1, ..., M$, given $(x_i, \alpha_i)$, infer the
perturbations by retrieving the arg-max in the local maximum pooling
of step [1]:
\[
(\Delta x_{m, i}, \Delta \alpha_{m, i}) = \arg\max_{\Delta x, \Delta
\alpha} |R_m(x_i+\Delta x, \alpha_i+\Delta \alpha)|^2.
\]
Let $c_{m, i} \leftarrow R_m(x_i+\Delta x_{m, i}, \alpha_i+\Delta
\alpha_{m, i})$, and update $R_m(x, \alpha) \leftarrow$ 0 if
$\corr[B_{x, s, \alpha}, B_{x_i+\Delta x_{m, i}, s,
\alpha_i+\Delta \alpha_{m, i}}] > \epsilon$.  Then estimate
$\hat{\lambda}_i = \mu^{-1}(\sum_{m=1}^{M} h(|c_{m, i}|^2)/M)$.

[3] Stop if $i = n$, else go back to 1.

There are two modifications to the original shared matching pursuit
in section (\ref{sect:SMP0}). (1) In step [1], we apply the
saturation function $h()$ to the response. This is justified by
maximum likelihood based on the exponential family model
(\ref{eq:EXP}). Intuitively, it means we discount very large
responses, because there can also be occasional (though less
frequent) large responses caused by the edges in natural image
background. (2) In step [2], the arg-max element $B_{x_i+\Delta x_{m,
i}, s, \alpha_i+\Delta \alpha_{m, i}}$ directly inhibits nearby
elements whose correlation with it is greater than a tolerance
$\epsilon$, instead of explaining away from $U_m$ and inhibiting
nearby elements indirectly. The correlation is defined as the square
of the inner product. This hard inhibition is to approximately
enforce the orthogonality assumption in section (\ref{sect:SMP0}). 
{\em $n$ may be adaptively chosen by setting a threshold on $\hat{\lambda}_i$.  }

\subsection{Template matching and shape filter}
\label{sect:shapefilter}

After learning the template from training images $\{\I_m\}$, for a
testing image $\I$ that contains an instance of the object, we can
detect it by scanning the template over the whole image and compute
the template matching score.

[1] For every pixel $X$, compute the log-likelihood $l(X)$, which
serves as the template matching score at putative location $X$:
\begin{eqnarray}
\sum_{i=1}^{n} \big[ \lambda_i \max_{\Delta x, \Delta \alpha}
h(|\langle \I, B_{X+x_i+\Delta x, s, \alpha_i + \Delta
\alpha}\rangle|^2)  - \log Z(\lambda_i) \big]. \label{eq:llr}
\end{eqnarray}

[2] Find maximum likelihood $\hat{X} = \arg\max_{X} l(X)$. For $i =
1, ..., n$, inferring perturbations by retrieving the arg-max in the
local maximum pooling in step [1]:
\[
(\Delta x_i, \Delta \alpha_i) = \arg\max_{\Delta x, \Delta \alpha}
|\langle \I, B_{\hat{X}+x_i+\Delta x, s, \alpha_i + \Delta
\alpha}\rangle|^2.
\]

[3] Return the location $\hat{X}$, and the translated and deformed
template $(B_{\hat{X}+x_i+\Delta x_i, s, \alpha_i + \Delta \alpha_i},
i = 1, ..., n)$.

$l(\hat{X})$ can also be used for classification.

{\em Shape filter}. The template $\B$ can be viewed as a shape
filter, and the template matching score $l(X)$ is the filter response
at location $X$. We may even write $\B*\I(X) = l(X)$. Step [2] finds
the object and deforms the template. Step [3] not only returns a
bounding box, but also sketches the detected object by the deformed
template.

{\em Rotation and multi-resolution}. We can rotate the template and
scan the template over multiple resolutions of the original image, to
account for uncertainties of the orientation and scale of the
object.

\section{Unsupervised learning of codebooks of active basis
templates}

In Olshausen-Field model (\ref{eq:OF1}), the coefficients are assumed to 
be independent for simplicity. 
A natural question is how to correct this
assumption, or specifically, how to add another layer of model on all the
coefficients?

Since we argue that the Olshausen-Field model (\ref{eq:OF1}) with a
dictionary of self-similar wavelets as in (\ref{eq:OF2}) is
essentially a spatial point process as explicated in (\ref{eq:OF3}),
the model on top of all the coefficients should focus on the
geometric patterns formed by the wavelet elements with non-zero
coefficients. In particular, we may search for recurring
compositional patterns of the selected wavelet elements. These
compositional patterns can be modeled by active basis models.

In this section, we shall specify our representation where each
representational element is an active basis template. Then we shall
describe the algorithm for learning codebooks of active basis
templates from training images.

\subsection{Representation}

In this section, we strive to write down our model in a form that is
analogous to the Olshausen-Field model, by using compactified notation.

{\em Compactified notation}. As the first step of this exercise of
compactification, let us slightly generalize the active basis model
by assuming that the template may appear at location $X_m$ in image
$\I_m$, then
\begin{eqnarray}
\I_m &=& \sum_{i=1}^{n} c_{m, i} B_{X_m + x_i + \Delta x_{m, i}, s,
\alpha_i + \Delta \alpha_{m, i}} + U_m \nonumber \\
&=& C_m \B_{X_m}+U_m, \label{eq:pack1}
\end{eqnarray}
where $\B = (B_{x_i, s, \alpha_i}, i = 1, ..., n)$ is the original
template, $\B_{X_m} = (B_{X_m + x_i + \Delta x_{m, i}, s, \alpha_i +
\Delta \alpha_{m, i}}, i = 1, ..., n)$ is the deformed template spatially
translated to $X_m$. Equation (\ref{eq:pack1}) is written in what we
call the compactified notation.

For each image $\I$ and each $X$, we can also define
\begin{eqnarray}
\langle \I, \B_X \rangle = \sum_{i=1}^{n} \big[ \lambda_i
\max_{\Delta x, \Delta \alpha} h(|\langle \I, B_{X+x_i+\Delta x, s,
\alpha_i + \Delta \alpha}\rangle|^2) \nonumber\\ - \log Z(\lambda_i)
\big],  \label{eq:llr2}
\end{eqnarray}
which is $l(X) = \B*\I(X)$ in section (\ref{sect:shapefilter}).

As the next step of this compactification exercise, in addition to
spatial translation and deformation, we can also rotate and scale
template. So a more general model is $\I_m = C_m \B_{X_m, S_m, A_m}+U_m$,
where $X_m$ is the location, $S_m$ is the scale, and $A_m$ is the
orientation. The scaling of the template is implemented by changing
the resolution of the original image. We adopt the convention that
whenever the notation $\B$ appears in image representation, it always
means the deformed template.

We can also define $\langle \I, \B_{X, S, A}\rangle$
in a similar way as in equation (\ref{eq:llr2}).

{\em Compactified representation}. Now suppose we have a codebook of
$T$ templates, $\{\B^{(t)}, t = 1, ..., T\}$. Then we can represent
an image by $K$ templates that are spatially translated, rotated, 
scaled and deformed versions of these $T$ templates in the codebook.
\begin{eqnarray}
\I_m = \sum_{k=1}^{K} C_{m, k} \B^{(t_{m, k})}_{X_{m, k}, S_{m, k},
A_{m, k}} + U_m, \label{eq:MODEL}
\end{eqnarray}
where each $\B^{(t_k)}_{X_{m, k}, S_{m, k}, A_{m, k}}$ is obtained by
translating the template $\B^{(t_k)}$ in the codebook to location
$X_{m, k}$, scale it to scale $S_{m, k}$, rotate it to orientation
$A_{m, k}$, and deform it to match $\I_m$. Again, the scaling of a
template can be implemented by changing the resolution of the image.

{\em Packing and unpacking}. The above representation is in analogy
to  equation (\ref{eq:OF3}) in section (\ref{sect:OF}), which we
copy here: $\I_m = \sum_{i=1}^{n} c_{m, i} B_{x_{m, i}, s_{m, i},
\alpha_{m, i}} + U_m$. The difference is that each $\B^{(t_k)}_{X_{m,
k}, S_{m, k}, A_{m, k}}$ is itself a group of wavelet elements that
follow a certain composition pattern $t_k$. Because of such grouping or packing, the number of templates $K$ needed to code
$\I_m$ is expected to be smaller than $n$. Specifically, if each
template is a group of $g$ wavelet elements, then $n = Kg$ (if there
is no overlaps between templates). In other words, we can unpack
model (\ref{eq:MODEL}) into the wavelet expansion model
(\ref{eq:OF3}). The reason that it is advantageous to pack the
wavelets into groups is that these groups exhibit $T$ types of
recurring compositional patterns to be discovered automatically by
the learning algorithm.

\subsection{Learning algorithm} \label{sect:learning}

We want to learn the codebook $\{\B^{(t)}, t = 1, ..., T\}$ from
training images $\{\I_m, m = 1, ..., M\}$, and at the same time,
infer the representation for each $\I_m$ in terms of the model
(\ref{eq:MODEL}). The objective function is the log-likelihood
\begin{eqnarray}
\sum_{m=1}^{M} \sum_{k = 1}^{K} \langle \I_m, \B^{(t_{m, k})}_{X_{m,
k}, S_{m, k}, A_{m, k}} \rangle, \label{eq:LLK}
\end{eqnarray}
which is the sum of the log-likelihood of all the selected templates (recall that we should not use the
simple squared loss, which implicitly and incorrectly assumes that
$U_m$ is white noise, see section (\ref{sect:SMP1})).
$K$ is assumed to be small and  may vary for different $\I_m$.

The learning algorithm is guided by the above objective function. It iterates an
image encoding step and a codebook revision step.

{\em Step 1: Image encoding by template matching pursuit}. Suppose we
are given the current codebook $\{\B^{(t)}, t = 1, ..., T\}$. Then
for each $\I_m$, $m = 1, ..., M$, the template matching pursuit
algorithm seeks to represent it by sequentially selecting a small
number of templates from the codebook in order to pave $\I_m$.

[0] Initialize the response maps $\R_m^{(t)}(X, S, A) \leftarrow
\langle \I_m, \B^{(t)}_{X, S, A}\rangle$, for all $(X, S, A, t)$.
This can be accomplished by first rotating the template to
orientation $A$, and then scanning the rotated template over the image
zoomed to the resolution corresponds to scale $S$. The larger the $S$
is, the smaller the resolution is. Let $k \leftarrow 1$.

[1] Select the translated, rotated, scaled and deformed template by
finding the global maximum of the response maps: $(X_{m, k}, S_{m,
k}, A_{m, k}, t_{m, k}) = \arg\max_{X, S, A, t} \R_m^{(t)}(X, S, A)$.

[2] Let the selected arg-max template inhibits overlapping candidate
templates, i.e., for all $(X, S, A, t)$, if corr$[\B^{(t)}_{X, S,
A}$, $\B^{(t_{m, k})}_{X_{m, k}, S_{m, k}, A_{m, k}}]$ $> \epsilon$,
then set the response $\R_m^{(t)}(X, S, A) \leftarrow -\infty$. $k
\leftarrow k+1$.

[3] Stop if all $\R_m^{(t)}(X, S, A, t) \leq 0$. Otherwise go to step
[1].

This template matching pursuit algorithm implements a hard inhibition. The corr$[\B_1$, $\B_2]$ here
is simply defined as the overlap between the bounding boxes of the
two templates $\B_1$ and $\B_2$. It can also be defined in a more
refined manner in terms of the inner products between the constituent
elements of $\B_1$ and $\B_2$. More rigorously, we could even update
the residual image by $U_m \leftarrow U_m - C_m\B^{(t_{m,
k})}_{X_{m,k}, S_{m, k}, A_{m, k}}$ as in the original version of
matching pursuit. But the current simplified version is more
efficient and works well enough.

{\em Step 2: Codebook re-learning by shared matching pursuit}. For
each $t = 1, ..., T$ in the codebook, the shared matching pursuit
algorithm re-learns $B^{(t)}$ from all the image patches that are
covered by $\B^{(t)}$ by selecting the constituent Gabor wavelet
elements of $\B^{(t)}$. This consists of the following two
operations:

[1] Image patch cropping. For each $\I_m$, go through all the
selected templates $\{\B^{(t_{m, k})}_{X_{m, k}, S_{m, k}, A_{m, k}},
\forall k\}$ that pave $\I_m$. If $t_{m, k} = t$, then crop the
image patch of $\I_m$ (at the resolution that corresponds to $S_{m,
k}$) covered by the bounding box of the template $\B^{(t_{m,
k})}_{X_{m, k}, S_{m, k}, A_{m, k}}$.

[2] Template re-learning. Re-learn template $\B^{(t)}$ from all the
image patches covered by $\B^{(t)}$ that are cropped in step [1], and
with their bounding boxes aligned.  The learning is accomplished by
the shared matching pursuit algorithm of section (\ref{sect:SMP1}),
which sequentially selects Gabor wavelet elements to encode all these
image patches simultaneously.

In practice, we do not have to deal with these image patches
directly. We only need to crop the maps of Gabor filter responses,
and feed them into the shared sketch algorithm.

{\em Polarization or specialization}. The learning algorithm starts
from a codebook of templates learned from randomly cropped image
patches. As a result, the initial templates are rather random and
meaningless, and the differences among them are small. However, as
the algorithm proceeds, the small differences among the initial
templates quickly start a polarizing or specializing process, where
the templates become more and more different, and they specialize in
coding different types of image patches. One may start the algorithm
multiple times and select the one that achieves the maximum of the
log-likelihood (\ref{eq:LLK}).

{\em Generalized matching pursuit in both steps}. Both the encoding and
re-learning steps of the learning algorithm are generalizations of matching
pursuit. In fact, the whole learning algorithm can be viewed as an
encoding algorithm, which seeks to automatically discover the
recurring compositional patterns that are otherwise overlooked by the
plain matching pursuit algorithm. The re-learning of each template
can be viewed as encoding multiple image patches by a single
template. 

{\em No early decision on wavelet representation} For learning the
codebooks, it is tempting to first apply the plain matching pursuit
algorithm to each training image, and then search for recurring
compositional patterns in the selected wavelets produced by the plain
matching pursuit. The problem with such a two-step sequential scheme
is that the wavelet representation produced by the plain matching
pursuit is an early decision or early commitment. Presumably, there
may be many other wavelet representations that are equally good or
comparable in terms of sparsity, but they may be much more regular in
terms of forming recurring compositional patterns. So we have to
obtain the wavelet representation and discover the compositional
patterns simultaneously or by an iterative algorithm. Not making
early decision on wavelet expansion or edge detection is a key
difference between our learning algorithm and those of Zhu et al.
\cite{Zhu_texton} and L. Zhu et al. \cite{LeoZhu}.

{\em Biological plausibility}. The Olshausen-Field model is a model
for simple V1 cells. The local max pooling of Riesenhuber and Poggio
\cite{Riesenhuber_Poggio} and the arg-max retrieval and inhibition of
the active basis model may be related to complex V1 cells. The
codebooks of active basis templates may be related to V2 cells.

\section{Experiments}

This section presents experimental results that show that our
learning algorithm is capable of learning explicit and meaningful
codebooks that give explicit and meaningful representations of
training images.

\subsection{Learning codebooks of textons}

We first apply the learning algorithm to texture images, where we
learn codebooks of active basis templates that can be used to pave
the texture images. The learned templates are recurring patterns
within the texture images, and we may call them textons, following
Julesz. Figure \ref{fig:textons} shows some examples. {\em Some members 
in the codebook may disappear if they fail to be selected in the encoding step.}

\begin{figure}[htb]
\begin{tabular}{cccc}
  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
 \includegraphics[height=.04\textwidth]{img/template_sorted_brick.png} & \hspace{-53mm}\includegraphics[height=.1\textwidth]{img/brick.jpg} & \hspace{8mm}\includegraphics[height=.1\textwidth]{img/sketch_image1_brick.png}
\end{tabular}\\
\begin{tabular}{cccc}
  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
   \includegraphics[height=.04\textwidth]{img/template_sorted_leaf2.png} & \hspace{-53mm}\includegraphics[height=.1\textwidth]{img/leaf2.jpg} & \includegraphics[height=.1\textwidth]{img/sketch_image1_leaf2.png}
\end{tabular}\\
\begin{tabular}{cccc}
  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
   \includegraphics[height=.04\textwidth]{img/template_sorted_peanut.png} & \hspace{-53mm}\includegraphics[height=.1\textwidth]{img/peanut.jpg} & \includegraphics[height=.1\textwidth]{img/sketch_image1_peanut.png}
\end{tabular}\\
\begin{tabular}{cccc}
  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
   \includegraphics[height=.04\textwidth]{img/template_sorted_pave4.png} & \hspace{-53mm}\includegraphics[height=.1\textwidth]{img/pave4.jpg} & \hspace{11mm}\includegraphics[height=.1\textwidth]{img/sketch_image1_pave4.png}
\end{tabular}\\
\begin{tabular}{cccc}
  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
   \includegraphics[height=.04\textwidth]{img/template_sorted_grape2.png} & \hspace{-53mm}\includegraphics[height=.1\textwidth]{img/grape2.jpg} & \includegraphics[height=.1\textwidth]{img/colorsketch_image1_grape2.png}
\end{tabular}\\
\begin{tabular}{cccc}
   \includegraphics[height=.04\textwidth]{img/template_sorted_fence.png} & \hspace{-53mm}\includegraphics[height=.1\textwidth]{img/fence.jpg} & \includegraphics[height=.1\textwidth]{img/colorsketch_image1_fence.png}
\end{tabular}\\
\caption{Learning codebooks of texons. For each texture image, we learn a codebook of 3 textons (for the wall example only 2 survived).
The sizes of the images: bricks (291 by 271), leaves (509 by 382), peanuts (289 by 217), fabric (275 by 275), grapes (555 by 417), and fence (554 by 418). Each active basis template (i.e. entry in the texton codebooks) is of size 60 by 60 pixels, with no more than 12 selected wavelets. A small number of translated and deformed templates are used to pave each input image: bricks (43), leaves (94), peanuts (40), fabric (29), grapes (83), and fence (105). For the last two examples, we use different colors for different textons. }
\label{fig:textons}
\end{figure}

{\em Allowing overlaps}. In our current implementation, we allow
quite some overlaps between the translated and deformed templates to avoid 
that situations where small pieces of images fall through the cracks and unpaved. 

\subsection{Learning codebooks of partons}

We also apply the learning algorithm to object images {\em without assuming 
image alignment}. The learned codebooks of active basis templates correspond
to parts of the objects as well as the backgrounds. These templates
are recurring patterns across multiple images. We may call them
partons, following the physicist Feynman. Figure \ref{fig:partons} shows
some examples.


\begin{figure}[htb]
\includegraphics[height=.07\textwidth]{img/template_sorted_car.png}  \\
\includegraphics[height=.05\textwidth]{img/car9.jpg}
\includegraphics[height=.05\textwidth]{img/sketch_image9.png} \hspace{1mm}
\includegraphics[height=.05\textwidth]{img/car33.jpg}
\includegraphics[height=.05\textwidth]{img/sketch_image33.png} \\
\includegraphics[height=.05\textwidth]{img/car48.jpg}
\includegraphics[height=.05\textwidth]{img/sketch_image48.png} \hspace{1mm}
\includegraphics[height=.05\textwidth]{img/car16.jpg}
\includegraphics[height=.05\textwidth]{img/sketch_image16.png}\\

%-----
\includegraphics[height=.08\textwidth]{img/template_sorted_pretty.png} \hspace{-7mm}\\
\includegraphics[height=.07\textwidth]{img/pretty1.jpg}
\includegraphics[height=.07\textwidth]{img/sketch_image1_pretty.png} \hspace{1mm}
\includegraphics[height=.07\textwidth]{img/pretty4.jpg}
\includegraphics[height=.07\textwidth]{img/sketch_image4_pretty.png} \hspace{1mm}
\includegraphics[height=.07\textwidth]{img/pretty8.jpg}
\includegraphics[height=.07\textwidth]{img/sketch_image8_pretty.png}\\

%-------
\includegraphics[height=.07\textwidth]{img/template_sorted_tea.png}\\
\includegraphics[height=.05\textwidth]{img/tea19.jpg}
\includegraphics[height=.05\textwidth]{img/sketch_image19_tea.png} \hspace{1mm}
\includegraphics[height=.05\textwidth]{img/tea53.jpg}
\includegraphics[height=.05\textwidth]{img/sketch_image53_tea.png} \hspace{1mm}
\includegraphics[height=.05\textwidth]{img/tea60.jpg}
\includegraphics[height=.05\textwidth]{img/sketch_image60_tea.png} \hspace{1mm}
\includegraphics[height=.05\textwidth]{img/tea23.jpg}
\includegraphics[height=.05\textwidth]{img/sketch_image23_tea.png}\\

%-------
\includegraphics[height=.07\textwidth]{img/template_sorted_horse.png}\\
\includegraphics[height=.06\textwidth]{img/horse2.jpg}
\includegraphics[height=.06\textwidth]{img/sketch_image2_horse.png} \hspace{1mm}
\includegraphics[height=.06\textwidth]{img/horse12.jpg}
\includegraphics[height=.06\textwidth]{img/sketch_image12_horse.png} \hspace{1mm}
\includegraphics[height=.06\textwidth]{img/horse10.jpg}
\includegraphics[height=.06\textwidth]{img/sketch_image10_horse.png}\\

\caption{Learning codebooks of partons. For each data set, we learn a codebook of $20$ partons (for the flower example 19 survived).
The number of training images: car (50), flower (9), teapot (61), and horse (50).
Each active basis template (i.e. entry in the codebook) is of size 60 by 60 pixels, with no more than 12 selected wavelets. All images are resized to around 22,500 pixels. }
\label{fig:partons}
\end{figure}


\begin{table}[htb]
\begin{center}
%% increase table row spacing, adjust to taste
\renewcommand{\arraystretch}{1.3}
\caption{Classification accuracies on face and bicycle datasets. } \label{tab:classification_accuracy}
% The array package and the MDW tools package offers better commands
%% for making tables than plain LaTeX2e's tabular which is used here.
\begin{tabular}{|c|c|c|c|}
\hline
dataset&method&codebook size& accuracy\\
\hline\hline
Face&\begin{tabular}{c}
       % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
       k-means+sift \\
       k-means+sift \\
       our approach \\
     \end{tabular}
& \begin{tabular}{c}
       % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
       $50$ \\
       $500$ \\
       $20$ \\
     \end{tabular}  &\begin{tabular}{c}
       % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
       $84.15\%$ \\
       $89.31\%$ \\
       $99.56\%$ \\
     \end{tabular}\\
\hline
Bicycle&\begin{tabular}{c}
       % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
       k-means+sift \\
       k-means+sift \\
       our approach \\
     \end{tabular}
&
\begin{tabular}{c}
       % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
       $50$ \\
       $500$ \\
       $50$ \\
     \end{tabular}
&
\begin{tabular}{c}
       % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
       $53.95\%$ \\
       $76.69\%$ \\
       $81.23\%$ 
     \end{tabular}
\\
\hline
\end{tabular}
\end{center} 
\end{table}



\subsection{Using codebooks of partons for classification}

The learned codebooks of partons can be used for image classification.
Specifically let $\{\B^{(t)}, t = 1, ..., T\}$ be a codebook learned
from positive training images, then for each testing image $\I$, we
scan rotated versions of $\B^{(t)}$ over multiple resolutions of
$\I$, and then take the global maximum of the template matching score
as defined in section (\ref{sect:shapefilter}). We threshold the
maximum score below at 0. This gives us $T$ scores. We then feed
these $T$ scores to SVM \cite{svm} for classification.

We report experimental results on two commonly used datasets: Caltech-101 face \footnote{downloaded from  http://www.robots.ox.ac.uk/~vgg/data3.html} and Graz02 bicycle \footnote{downloaded from http: //www.emt.tugraz.at/~pinz/data/GRAZ\_02/}. The Caltech-101 \cite{FeiFei} contains $450$ face images and $900$ background images. Among them, we randomly select $225$ face images for training, and all the other images for testing. The Graz02 dataset has $365$ bicycle images and $380$ background images. We train on $100$ randomly selected bicycle images, and test on all the remaining bicycle and background images. For each data set, the accuracies are computed by averaging over 10 repetitions. We compared our method with k-means + sift \cite{Kmean,SIFT} under identical experiment settings. The latter is perhaps the most commonly used approach for image classification. Table \ref{tab:classification_accuracy} presents the experimental results. Our approach achieves a higher classification accuracy than k-means + sift with much smaller codebook size. Figure \ref{fig:bikeface} gives the codebooks, each obtained from one repetition.

\begin{figure}[htb]
\includegraphics[height=.09\textwidth]{img/template_sorted_bike2.jpg}\\
\includegraphics[height=.05\textwidth]{img/bike25.jpg}
\includegraphics[height=.05\textwidth]{img/sketch_image25_bike.png} \hspace{1mm}
\includegraphics[height=.05\textwidth]{img/bike70.jpg}
\includegraphics[height=.05\textwidth]{img/sketch_image70_bike.png} \hspace{1mm}
\includegraphics[height=.05\textwidth]{img/bike83.jpg}
\includegraphics[height=.05\textwidth]{img/sketch_image83_bike.png}\\
\vspace{1mm}\\
\includegraphics[height=.09\textwidth]{img/template_iter10face.png}
%\includegraphics[height=.07\textwidth]{img/template_sorted_face.png}\\
%\includegraphics[height=.05\textwidth]{img/face12.jpg}
%\includegraphics[height=.05\textwidth]{img/sketch_image12_face.png} 
%\includegraphics[height=.05\textwidth]{img/face44.jpg}
%\includegraphics[height=.05\textwidth]{img/sketch_image44_face.png} 
%\includegraphics[height=.05\textwidth]{img/face77.jpg}
%\includegraphics[height=.05\textwidth]{img/sketch_image77_face.png}
\caption{For bicycles we learn a codebook of 50 partons (48 survived, each parton has no more than 15 selected wavelets, size is 72 by 72 pixels) from 100 training images (each about 40,000 pixels). For faces, we learn a codebook of $20$ templates (each template has no more than 40 selected wavelets, size is 120 by 120 pixels) from 225 images (each about 40,000 pixels).}
\label{fig:bikeface}
\end{figure}

\section{Conclusion}

Although our work makes use of active basis model, our model is novel, 
and our learning algorithm is much more versatile and powerful than that 
used in the active basis model \cite{Wu_active_basis_IJCV}. 
\\\\
{\bf Reproducibility}: Data and code can be found in the supplementary materials.

{\small
\bibliographystyle{abbrv}
\begin{thebibliography}{99}

\bibitem{Kmean}
G Csurka, C Dance, L Fan, J Willamowski, and C Bray.
\newblock Visual categorization with bags of keypoints. 
\newblock {\em Workshop of ECCV}, 2004.

\bibitem{FeiFei}
L Fei-Fei, R Fergus, and P Perona. 
\newblock Learning generative visual models from few training examples: an incremental
Bayesian approach tested on 101 object categories. 
\newblock{\em CVPR Workshop}, 2004.

\bibitem{Friedman_projection_pursuit}
JH Friedman.
\newblock Exploratory projection pursuit.
\newblock {\em Journal of the American Statistical Association},
82:249--266, 1987.

\bibitem{SIFT}
D Lowe. 
\newblock Distinctive image features from scale-invariant
keypoints. 
\newblock{\em IJCV}, 2004.

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

\bibitem{Olshausen_Field}
BA Olshausen and DJ Field.
\newblock Emergence of simple-cell receptive field properties by
learning a sparse code for natural images.
\newblock {\em Nature}, 381:607--609, 1996.

\bibitem{Olshausen_wavelet}
BA Olshausen, P Sallee and MS Lewicki.
\newblock Learning Sparse image codes using a wavelet pyramid
architecture.
\newblock {\em NIPS},
13: 887-893, 2001.

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

\bibitem{svm}
VN Vapnik.
\newblock {\em The nature of statistical learning theory}. 2000.

\bibitem{Wu_active_basis_IJCV}
YN Wu, Z Si, H Gong, and SC Zhu.
\newblock Learning active basis model for object detection and
recognition.
\newblock {\em IJCV},
90:198--235, 2010.

\bibitem{LeoZhu}
L Zhu, C Lin, H Huang, Y Chen, and A Yuille.
\newblock Unsupervised structure learning: hierarchical recursive
composition, suspicious coincidence and competitive exclusion.
\newblock In {\em ECCV}, 2008.

\bibitem{Zhu_texton}
SC Zhu, C Guo, Y Wang, and Z Xu.
\newblock What are textons?
\newblock {\em IJCV},
62:121--143, 2005

\bibitem{Zhu_grammar}
SC Zhu and DB Mumford.
\newblock A stochastic grammar of images.
\newblock {\em Foundations and Trends in Computer Graphics and
Vision}, 2:259--362, 2006.


\end{thebibliography}

}


\end{document}
