\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 large 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 the 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 active basis models from an image of pavement. 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.  }
\label{fig:pave}
\end{figure}


Figure illustrates the basic idea. Figure displays an image of
pavement. Figure displays a codebook of 3 templates learned from the
pavement image. Each template is represented by an active basis
model. Specifically, each template is a composition of a small number
of Gabor wavelet elements automatically selected from a dictionary of
Gabor wavelets. {\em Each selected Gabor wavelet is depicted by a bar
at the same location, orientation and length as the Gabor wavelet}.
The selected Gabor wavelets are allowed to perturb their locations
and orientations, so that the template becomes deformable, and the
linear basis formed by the selected Gabor wavelet elements become
active.

Figure shows the representation of the original pavement by the
templates in the learned codebook. Here the word ``pavement'' also
serves as a good metaphor because the representation depicted by
figure is essentially a ``pavement'' made of the active basis
templates in the codebook of figure. Specifically, each of the
templates in figure is a spatially translated, scaled, rotated and
deformed version of a template in the codebook in figure. The
bounding box of each template is also shown in figure (the deformed
templates in figure are allowed to have some overlaps with each
other). Clearly, the image can be represented 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 active basis models from an image of pavement. 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.  }
\label{fig:pave}
\end{figure}

Figure 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, the 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 all the Gabor wavelets, in order to encode all
these patches simultaneously.

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.

\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 large 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 investigated 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$ 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 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 range 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 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.

{\em Assumption on orthogonality}. 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 in 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 the coefficients
$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 an 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) = 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[\langle B_{x, s, \alpha}, B_{x_i+\Delta x_{m, i}, s,
\alpha_i+\Delta \alpha_{m, i}} \rangle] > \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}).

\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:
\[
(\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 overall orientation and scale of the
object.

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

In Olshausen-Field model \cite{Olshausen_Field} and its probabilistic
interpretation \cite{Lewicki}, the coefficients of the linear
additive representation are treated separately or 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 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 the 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 the 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_i, s,
\alpha_i + \Delta \alpha_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_i, s, \alpha_i +
\Delta \alpha_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}$,
where $X_m$ is location, $S_m$ is 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$ appear in image representation, it always
means the deformed template.

Then similarly, 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, and
scaled 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
with 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 varying for different $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 scan 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 is essentially the matching
pursuit algorithm with 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 by 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 m, 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 Matching pursuit in both steps}. Both the encoding and
re-learning steps of the learning algorithm are based on 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. It may be possible to replace matching pursuit by more
principled sparse coding algorithms such as basis pursuit
\cite{Chen_basis_pursuit} or equivalently Lasso \cite{lasso}.

{\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 the one by 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 show some examples.

\begin{figure}[htb]
\includegraphics[height=.04\textwidth]{img/template_sorted_brick.png}\\
\includegraphics[height=.16\textwidth]{img/brick.jpg}
\includegraphics[height=.16\textwidth]{img/sketch_image1_brick.png}
\caption{}
\label{fig:brick}
\end{figure}

\begin{figure}[htb]
\includegraphics[height=.04\textwidth]{img/template_sorted_leaf2.png}\\
\includegraphics[height=.15\textwidth]{img/leaf2.jpg}
\includegraphics[height=.15\textwidth]{img/sketch_image1_leaf2.png}
\caption{}
\label{fig:leaf}
\end{figure}

\begin{figure}[htb]
\includegraphics[height=.04\textwidth]{img/template_sorted_peanut.png}\\
\includegraphics[height=.14\textwidth]{img/peanut.jpg}
\includegraphics[height=.14\textwidth]{img/sketch_image1_peanut.png}
\caption{}
\label{fig:peanut}
\end{figure}

\begin{figure}[htb]
\includegraphics[height=.04\textwidth]{img/template_sorted_pave4.png}\\
\includegraphics[height=.14\textwidth]{img/pave4.jpg}
\includegraphics[height=.14\textwidth]{img/sketch_image1_pave4.png}
\includegraphics[height=.14\textwidth]{img/colorsketch_image1_pave4.png}
\caption{}
\label{fig:pave4}
\end{figure}

\begin{figure}[htb]
\includegraphics[height=.04\textwidth]{img/template_sorted_grape2.png}\\
\includegraphics[height=.14\textwidth]{img/grape2.jpg}
\includegraphics[height=.14\textwidth]{img/colorsketch_image1_grape2.png}
\caption{}
\label{fig:grape2}
\end{figure}

\begin{figure}[htb]
\includegraphics[height=.04\textwidth]{img/template_sorted_fence.png}\\
\includegraphics[height=.14\textwidth]{img/fence.jpg}
\includegraphics[height=.14\textwidth]{img/colorsketch_image1_fence.png}
\caption{}
\label{fig:grape2}
\end{figure}

{\em Allowing overlaps}. In our current implementation, we allow
quite some overlaps between the translated and deformed templates.
The overlapped deformed templates often share the same or similar
constituent wavelet elements. Allowing some overlaps may avoid
leaving small pieces of images unexplained.

\subsection{Learning codebooks of partons}

We also apply the learning algorithm to object images that are not
aligned. The learned codebooks of active basis templates correspond
to parts of the objects as well as the backgrounds. These templates
are recuring patterns across multiple images. We may call them
partons, following physicist Feynman.

\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} 
\caption{ }
\label{fig:car}
\end{figure}

\begin{figure}[htb]
\includegraphics[height=.08\textwidth]{img/template_sorted_pretty.png}\\
\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} 
\caption{ }
\label{fig:pretty}
\end{figure}

\begin{figure}[htb]
\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} 
\caption{ }
\label{fig:tea}
\end{figure}

\begin{figure}[htb]
\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{ }
\label{fig:tea}
\end{figure}
\subsection{Using codebooks of partons for classification}

The learned codebooks of partons can be used for 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. In our
experiments, we set $T = 50$. This number is smaller than existing
codebook based methods.


\begin{figure}[htb]
\includegraphics[height=.14\textwidth]{img/template_sorted_bike.png}\\
\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} 
\caption{ }
\label{fig:bike}
\end{figure}


\section{Conclusion}

This paper proposes a method for unsupervised learning a codebook of
active basis templates from training images. Unlike commonly used
codebook learning method, our method gives explicit and meaningful
representations of both the codebooks and the training images. The
learned codebooks appear to be useful for classification.

{\small
\section*{Reproducibility}

Data, code, html and readme files for reproducing all the figures and
tables presented in this paper are contained in the supplementary
materials.


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

\bibitem{Chen_basis_pursuit}
S.~Chen, D.~Donoho, and M.~A. Saunders.
\newblock Atomic decomposition by basis pursuit.
\newblock {\em SIAM Journal on Scientific Computing}, 20:33--61,
1999.

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

\bibitem{Lewicki}
M.~S. Lewicki and B.~A. Olshausen
\newblock Probabilistic framework for the adaptation and comparison
of image codes.
\newblock {\em Journal of the Optical Society of America},
16:1587--1601, 1999.

\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}
B.~A. Olshausen and D.~J. 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}
B.~A. Olshausen P. Sallee and M.~S. Lewicki.
\newblock Learning Sparse image codes using a wavelet pyramid
architecture.
\newblock {\em Advances in Neural Information Processing Systems},
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{lasso}
R.~Tibshirani.
\newblock Regression shrinkage and selection via the lasso.
\newblock {\em Journal of the Royal Statistical Society, B},
58:267--288, 1996.

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

\bibitem{Wu_active_basis_IJCV}
Y.~N. Wu, Z.~Si, H.~Gong, and S.-C. Zhu.
\newblock Learning active basis model for object detection and
recognition.
\newblock {\em International Journal of Computer Vision},
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 Proceedings of European Conference on Computer
Vision}, 2008.

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

\bibitem{Zhu_grammar}
S.~C. Zhu and D.~B. 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}
