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





\renewcommand{\baselinestretch}{1.2}


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

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


\newcommand{\Local}{./img_eps}

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

\begin{document}




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

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




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

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

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

\vspace{5mm}


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

\newpage

\tableofcontents

\newpage

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


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


\subsection{The model}

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

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



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



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

\subsection{The learning algorithm}

\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=SharedSketch.eps,height=3in}
\end{tabular}
\end{center}
\caption{Shared sketch algorithm. A selected element (colored
ellipsoid) is shared by all the training images. For each image, a
perturbed version of the element sketches a local edge segment near
the element. The elements of the active basis are selected
sequentially according to the Kullback-Leibler divergence of the
pooled distribution (colored solid curve) of filter responses from
the background distribution (black dotted curve). }
\label{fig:SharedSketch}
\end{figure}

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

For each element, a distribution of filter responses is pooled over
all the training images at the perturbed locations and orientations.
The elements are selected in an order according to the
Kullback-Leibler divergence of the pooled distribution (solid curve)
from a background distribution (dotted curve). The background
distribution is  pooled over natural images. The Kullback-Leibler
divergence can be simplified into a pursuit index that drives the
selection of the elements. After an element is selected, its
perturbed version explains away a small part of each training image,
and thereby inhibits nearby Gabor wavelet elements from coding the
same part of the image. So the selected elements of the active basis
are well spaced, and usually form a clear template.


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

\subsection{The computational architecture}

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

The recognition of the deformable template from an image can be
accomplished by a computational architecture of sum-max maps. Figure
(\ref{fig:SumMaxMaps}) illustrates this architecture. It starts from
convolving the image with Gabor filters at all the locations and
orientations. The filtered images become the first layer of the sum
maps, or $\S1$ maps, because each Gabor filter is a local summation
operator. In Figure (\ref{fig:SumMaxMaps}), the thin ellipsoids in
the $\S1$ maps illustrate the local filtering or summation
operation. Then a layer of max maps, or $\M1$ maps, is computed by
applying a local maximization operator to the $\S1$ maps. In Figure
(\ref{fig:SumMaxMaps}), the arrows in the $\M1$ maps illustrate that
the local maximization is taken over small perturbations of the
Gabor wavelets. This local maximization tells us how to deform the
active basis to match the image data. After that, a sum map, or
$\S2$ map, is computed by applying a local summation operator to the
$\M2$ maps. Specifically, we scan the active basis over the whole
image lattice, and for each pixel of the $\S2$ map, we compute a
weighted sum of the values of the $\M1$ maps. So this is another
layer of filtering operation, and can be considered a shape filter.
The weighted sum is over the elements of the active basis. It
computes the log-likelihood of the deformed active basis, because
the sum operator is applied to the $\M1$ maps. In Figure
(\ref{fig:SumMaxMaps}), the car template in the $\S2$ map
illustrates the active basis centered at one pixel.  We scan this
template over all the pixels to obtain the $\S2$ map, which scores
the template matching.

\subsection{Past work}

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

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

Our algorithmic architecture of sum-max maps is a variation on the
theme of Riesenhuber and Poggio's cortex-like hierarchical structure
of  simple and complex cells \cite{Riesenhuber_Poggio}. See also the
recent work of Mutch and Lowe \cite{Lowe} for further improvement.
Our architecture is essentially the same as Riesenhuber and Poggio's
structure up to the $\M1$ maps. However, our operator for computing
the $\S2$ map is different from that of Riesenhuber and Poggio. Our
operator is guided by a generative model based on Olshausen and
Field's representation. Specifically, our operator computes the
log-likelihood of the deformed active basis. It is a local summation
operator that can be interpreted as a shape filter.

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

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




%------------------------------------------------------------------------
\section{Active basis model}

\subsection{Gabor wavelets}

A Gabor function \cite{Daugman} is of the form $G(x, y) \propto
\exp\{ -[(x/\sigma_x)^2 + (y/\sigma_y)^2]/2\} e^{i {x}}$. We can
translate, rotate, and dilate $G(x, y)$ to obtain a general form of
Gabor wavelets: $B_{x, y, s, \alpha}(x', y') = G({\tilde{x}}/{s},
{\tilde{y}}/{s})/s^2, $ where $\tilde{x} =(x'-x)\cos\alpha -
(y'-y)\sin\alpha$, $\tilde{y} = (x'-x)\sin\alpha +
(y'-y)\cos\alpha$. $(x, y)$ is the central position, $s$ is the
scale parameter, and $\alpha$ is the orientation. The central
frequency of $B_{x, y, s, \alpha}$ is $\omega = 1/s$. $B_{x, y, s,
\alpha}$  = ($B_{x, y, s, \alpha, 0}$, $B_{x, y, s, \alpha, 1}$),
where $B_{x, y, s, \alpha, 0}$ is the even-symmetric Gabor cosine
component, and $B_{x, y, s, \alpha, 1}$ is the odd-symmetric Gabor
sine component. We always use Gabor wavelets as pairs of cosine and
sine components. We normalize both the Gabor sine and cosine
components to have zero mean and unit $l_2$ norm. For each $B_{x, y,
s, \alpha}$, $B_{x, y, s, \alpha, 0}$ and $B_{x, y, s, \alpha, 1}$
are orthogonal to each other.

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

For image $\I$ and Gabor wavelet $B_{x, y, s, \alpha, \eta}$, $\eta
= 0, 1$, the projection coefficient of $\I$ onto $B_{x, y, s,
\alpha, \eta}$ or the filter response is $\langle \I, B_{x, y, s,
\alpha, \eta} \rangle$ $= \sum_{x', y'} \I(x', y') B_{x, y, s,
\alpha, \eta}(x', y')$. We write $\langle \I, B_{x, y, s, \alpha}
\rangle$ = ($\langle \I, B_{x, y, s, \alpha, 0} \rangle$, $\langle
\I, B_{x, y, s, \alpha, 1} \rangle$). The local energy is $|\langle
\I, B_{x, y, s, \alpha} \rangle|^2$ = $\langle \I, B_{x, y, s,
\alpha, 0} \rangle^2$$+\langle \I, B_{x, y, s, \alpha, 1}
\rangle^2$.


\subsection{Active basis representation}

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

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

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

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

Recall that each $B_i$ or $B_{m, i}$ is a pair of Gabor cosine and
sine components. So $B_i = (B_{i, 0}, B_{i, 1})$ and $B_{m, i} =
(B_{m, i, 0}, B_{m, i, 1})$. Accordingly, $c_{m, i} = (c_{m, i, 0},
c_{m, i, 1})$ and $c_{m, i} B_{m, i} = c_{m, i, 0} B_{m, i, 0} +
c_{m, i, 1} B_{m, i, 1}$.

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

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


\subsection{Statistical formulation \label{sect:model}}

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

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

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

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

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

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

(1) White noise distribution. This is the distribution that is often
assumed in linear additive model, and is implicitly assumed in the
least squares criterion for modeling fitting.

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

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

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

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

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


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

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

\section{Shared sketch algorithm}

In this section, we first describe a shared matching pursuit process
that combines matching pursuit \cite{Mallat_Zhang} and projection
pursuit \cite{Friedman} . Then we present the shared sketch process
which is a simplified and approximate version of the shared matching
pursuit process. After that, we present the shared sketch algorithm.
Finally, we engage in a discussion on designing the pursuit index
and scoring the template matching.

\subsection{Shared matching pursuit process}

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

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

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

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

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

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

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

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

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

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

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

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



\subsection{Shared sketch process}

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

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

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

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

In practice, we enforce approximate orthogonality. Two Gabor
elements $B_1$ and $B_2$ are approximately orthogonal if ${\rm
corr}(B_1, B_2) = \sum_{\eta_1=0}^1\sum_{\eta_2=0}^1 \langle B_{1,
\eta_1}, B_{2, \eta_2}\rangle^2 < \epsilon$, (e.g., $\epsilon =
.1$).


\subsection{The algorithm  \label{sect:la}}

We can replace the pursuit index $\KL(p_i(c)||q(c))$ by
$\sum_{m=1}^{M} h(r_{m, i})$, where $r_{m, i} = |c_{m, i}|^2$ is the
local energy, and $h(r)\geq 0$ is a monotone increasing function.

To get an intuitive idea of this index, let us consider a simplified
example where $h(r) = 1_{r>\xi}$, where $\xi$ is a threshold for
edge detection. More specifically, $h(r) = 1$ if $r\geq \xi$, and
$h(r) = 0$ otherwise. Then $\sum_{m=1}^{M} h(r_{m, i})$  simply
counts the number of detected edge segments in the training images
$\{\I_m, m = 1, ..., M\}$. That is, we select $B_i$ and perturb it
to $\{B_{m, i}\}$, so that $\{B_{m, i}\}$ sketch as many edge
segments as possible. For more general $h(r)$, the pursuit index
$\sum_{m=1}^{M} h(r_{m, i})$ prefers large values of $\{r_{m, i}\}$.

Before presenting our designs of $h(r)$, we first present the shared
sketch algorithm, which is exceedingly simple, only slightly more
complex than regular edge detection.

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

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

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

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

See Figure (\ref{fig:SharedSketch}) for an illustration. In Step 2
of the $i$-th iteration, if $B_{m, i}$ is used to sketch $\I_m$,
then $B_{m, i}$ inhibits nearby Gabor elements $B$ from encoding
$\I_m$ by setting $[\I_m, B] \leftarrow 0$. In Step 1 of the $i$-th
iteration, the optimization of $B_{m, i}$ will avoid those Gabor
elements that have been inhibited by $B_{m, 1}$, ..., $B_{m, i-1}$.
This enforces the approximate non-overlapping constraint.

\subsection{Pursuit index and template matching score}

{\em Parametric model.} The pursuit index $\sum_{m=1}^{M} h(r_{m,
i})$ comes from the following exponential family model:
\begin{eqnarray}
    p(c; \lambda) = \frac{1}{Z(\lambda)} \exp\{\lambda h(r)\} q(c),
    \label{eq:exponential}
\end{eqnarray}
where $\lambda>0$ is the parameter, $r = |c|^2$, and $Z(\lambda) =
\int \exp\{\lambda h(r)\}q(c)dc = \E_{q}[\exp\{\lambda h(r)\}]$ is
the normalizing constant. Because $h(r)$ is an increasing function
and $\lambda>0$, $p(c; \lambda)$ puts more probability than $q(c)$
on those $c$ with large $r$, so model (\ref{eq:exponential}) is more
encouraging towards large $r$.

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

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

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

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

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

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

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

The following are two designs of $h(r)$.

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

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

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

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

A crude but simple approximation to ${\rm whiten}(r)$ is
\begin{eqnarray}
h(r) = {\rm threshold}(r) = \min(r, T), \label{eq:threshold}
\end{eqnarray}
where $T$ is a threshold (e.g., $T= 16$).With this approximation,
there is no need to estimate the marginal density $q(r)$ or the tail
probability $F(r)$.



\section{Sum maps and max maps}

\subsection{Learning algorithm using sum maps and max maps
\label{sect:SS}}

The learning algorithm in Subsection (\ref{sect:la}) can be
expressed more precisely in terms of the sum maps and max maps. See
Figure (\ref{fig:SumMaxMaps}) for an illustration of the sum maps
and max maps.

\vspace{.5cm}

\noindent{\bf Algorithm 1: Shared sketch algorithm}

\begin{enumerate}

\item Initialize $\S1$ maps and $\M1$ maps.

\begin{enumerate}

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

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

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

\end{enumerate}

\item Select the next basis element and its perturbed versions.

\begin{enumerate}

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


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

\item  Compute $\hat{\lambda}_i = \mu^{-1}(\sum_{m=1}^{M} h(r_{m,
i})/M)$.

\item For each $m = 1, ..., M$, compute
$\S2_m \leftarrow \S2_m + \hat{\lambda}_i h(r_{m, i}) - \log
Z(\hat{\lambda}_i)$.
\end{enumerate}

\item Inhibition on $\S1$ maps and updating $\M1$ maps.
\begin{enumerate}

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


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

\end{enumerate}

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

\vspace{.5cm}

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

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

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

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


The $\S2_m$ score evaluates the matching of $\I_m$ to the learned
template $\B$ according to Equation (\ref{eq:score}). The total
score $\sum_{m=1}^{M} \S2_m$ measures the overall alignment of all
the training images. This alignment score is very useful for
unsupervised learning.

In this paper, we hand pick the number of basis elements, $n$, in
our experiments. In principle, it can be automatically selected by
comparing $\sum_{m=1}^{M} h(r_{m, i})/M$ with the average of
$h(\M1_m(x, y, s, \alpha))$ in natural images or in the observed
image $\I_m$. If the former is no much greater than the latter, we
should stop the algorithm.


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

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

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

The following are the parameter values that we used in all the
experiments in this paper (unless otherwise stated).

{\em Parameter values.} Size of Gabor wavelets = $17 \times 17$.
$(x, y)$ is sub-sampled every 2 pixels. The orientation $\alpha$
takes $A = 15$ equally spaced angles in $[0, \pi]$. The
orthogonality tolerance is $\epsilon = .1$. The threshold is $T =
16$ in the threshold transformation (\ref{eq:threshold}). The
saturation level $\xi = 6$ in the sigmoid transformation
(\ref{eq:sigmoid}). The shift along the normal direction $d_{m, i}
\in [-b_1, b_1] = [-6, 6]$ pixels. The shift of orientation
$\delta_{m, i} \in [-b_2, b_2] = \{-1, 0, 1\}\times \pi/15$. We
tuned two parameters in our experiments. One is the number of
elements, $n$, which is usually between 30 and 80. The other
parameter is the resolution of the images at which we train the
model.

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

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

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

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

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

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

There are two shortcomings in the method in Experiment 1. One is
that the bounding boxes of the objects need to be given. In
Experiment 5, we shall show that this requirement can be relaxed.
The other problem is that our method does not do well on object
shapes with strong textures, such as zebras, leopards, tigers,
giraffes, etc.  This is because we take $q(\I_m)$ to be the
distribution of natural images. We need stronger $q(\I_m)$ to model
the texture in $\I_m$.



\subsection{Recognition algorithm using sum maps and max maps
\label{sect:recognition}}

After training the active basis model, specifically, after selecting
$\B = (B_i = B_{x_i, y_i, s, \alpha_i}, i = 1, ..., n)$, and
computing $\Lambda = (\lambda_i, i = 1, ..., n)$ (when there is no
danger of confusion, we use $\lambda_i$ instead of $\hat{\lambda}_i$
for notational simplicity), we can use the trained model to detect
and sketch the template in a testing image. Let $\I$ be a testing
image, which is larger than the bounding box of the template $\B$.
We assume that the bounding box of $\B$ is centered at origin $(x=0,
y=0)$. We can scan the template over image $\I$, and at each
position $(x, y)$, we fit the active basis model to the image within
the bounding box centered at $(x, y)$, and calculate the
log-likelihood ratio according to Equation (\ref{eq:score}) to score
the template matching. We call the map of the log-likelihood ratio
scores the $\S2$ map. The recognition algorithm is as follows.

\vspace{.5cm} \noindent{\bf Algorithm 2: Recognition algorithm}

\begin{enumerate}

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

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

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

\item Compute the MAX2 score:
$\M2 = \max_{x, y} \S2(x, y)$. Let $(\hat{x}, \hat{y})$ be the
position that achieves the maximum.

\item
For $i = 1, ..., n$, retrieve $\hat{x}_i, \hat{y}_i, \hat{\alpha}_i$
such that
\begin{eqnarray}
 \M1(\hat{x}+x_i, \hat{y}+y_i, s, \alpha_i) =  \S1(\hat{x}_i, \hat{y}_i, s,
 \hat{\alpha}_i)
\end{eqnarray}
in the local maximization operation (\ref{eq:MAXmap2}) in Step 2.
Then sketch the Gabor wavelet $B_{\hat{x}_i, \hat{y}_i, s,
\hat{\alpha}_i}$ by a bar.

\end{enumerate}
\vspace{.5cm}

The recognition algorithm is a lot simpler than the shared sketch
algorithm in Subsection (\ref{sect:SS}). It only involves the
computations of $\S1$, $\M1$, $\S2$ maps, and a $\M2$ score that
finds the best match to the template. Steps 1, 2, and 3 can be
accomplished by massively parallel computing. Step 3 computes the
log-likelihood ratio score at each $(x, y)$. This step can be
considered a shape filter for template matching. Like Gabor filters,
it is also a local weighted summation operator. See Figure
(\ref{fig:SumMaxMaps}) for an illustration.

For simplicity, we ignore the inhibition step in Step 3(a) of the
shared sketch algorithm in Subsection (\ref{sect:SS}). The elements
of the learned active basis $\B = (B_i, i = 1, ..., n)$ are usually
well spaced, because of the local inhibition in Step 3(a) of the
shared sketch algorithm. Therefore, when we deform the learned
active basis $\B$ to fit the testing image, the perturbed elements
$(B_{\hat{x}_i, \hat{y}_i, s, \hat{\alpha}_i}, i = 1, ..., n)$
usually do not have much overlap, and there is no much need for
inhibition.

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

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

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

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

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

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

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

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

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


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

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

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

Figure (\ref{fig:ex2_2_M2}.a) displays the $\M2$ scores of the 10
resolutions, where the maximum is achieved at the 5th resolution.
Figure (\ref{fig:ex2_2_M2}.b) displays the $\S2$ map at this optimal
resolution.

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

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

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

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

{\em Partial sums.} The elements of $\B = (B_i, i = 1, ..., n)$ are
selected sequentially in the shared sketch algorithm. In computing
the $\S2$ map in the recognition algorithm, we can sequentially
update $\S2(x, y)$ by $\S2_{i}(x, y) = \S2_{i-1}(x, y) + \lambda_i
h(\M1(x+x_i, y+y_i, s, \alpha_i)) - \log Z(\lambda_i)$. At each step
$i$, we may decide to stop updating those positions $(x, y)$ where
the partial sum $\S2_i(x, y)$ is too low. That is, we should abandon
the computation on the hopeless or unpromising pixels. This will
improve the efficiency of the algorithm.

The shortcoming of Experiment 2 is that we assume there is one
object in the testing image that matches the learned template. In
future work, we shall extend our method to deal with the situations
where there are multiple instances or zero instance of the template.



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

The deformable template $\B = (B_i, i = 1, ..., n)$ is parameterized
by $\Lambda = (\lambda_i, i = 1, ..., n)$. The log-likelihood ratio
score $\log [p(\I_m \mid \B_m)/q(\I_m)]$ in Equation
(\ref{eq:score}) is non-linear in $\lambda_i$, because of the $\log
Z(\lambda_i)$ term. This motivates us to introduce a simpler linear
score without explicit probabilistic assumptions.

We parameterize the deformable template  $\B = (B_i, i = 1, ..., n)$
by $\Theta = (\theta_i, i = 1, ..., n)$, where $\theta_i > 0$, and
$\Theta$ is normalized to have unit norm:  $\|\Theta\|^2 =
\sum_{i=1}^{n} \theta_i^2 = 1$. Then we have a vector $V =
\sum_{i=1}^{n} \theta_i B_i$. Because $\B$ is an active basis that
can be perturbed to $\B_{m} \approx \B$ for image $\I_m$, $V =
\sum_{i=1}^{n} \theta_i B_i$ is an active vector that can be
perturbed to $V_m = \sum_{i=1}^{n} \theta_i B_{m, i} \approx V$.
$V_m$ is approximately a  unit vector because $(B_{m, i}, i = 1,
..., n)$ are approximately orthogonal to each other.

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

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

The algorithm for learning $\B = (B_i, i = 1, ..., n)$ and  $\Theta
= (\theta_i, i = 1, ..., n)$ is very similar to the shared sketch
algorithm in Subsection (\ref{sect:SS}). The following are the
modifications.

(1) In Step 2(a), change the function $h(r)$ to ${\rm
whiten}(r)^{1/2}$.

(2) Change Step 2(c) to: Compute $\theta_i = \sum_{m=1}^{M} {\rm
whiten}(r_{m, i})^{1/2}/M$.

(3) Change Step 2(d) to: For $m = 1, ..., M$, compute $\S2_{m}
\leftarrow \S2_{m} + \theta_i {\rm whiten}(r_{m, i})^{1/2}$.

(4) At the end of the algorithm, compute $\|\Theta\|^2$, let
$\theta_i \leftarrow \theta_i/\|\Theta\|$, and $\S2_{m} \leftarrow
\S2_m/\|\Theta\|$.


We can also use active correlation for detection. The algorithm is
similar to the recognition algorithm in Subsection
(\ref{sect:recognition}). The only modification is to change Step 3
to:
\begin{eqnarray}
\S2(x, y) = \sum_{i=1}^{n} \theta_i {\rm whiten}(\M1(x+x_i, y+y_i,
s, \alpha))^{1/2}. \label{eq:ActiveCorrelation}
\end{eqnarray}

We can replace the whitening transformation (\ref{eq:whiten}) by the
thresholding transformation (\ref{eq:threshold}), which is simpler.

We obtain similar results in Experiments 1 and 2 using active
correlation.

\section{Evaluation on classification tasks}

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

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

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

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




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

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

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

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

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

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

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

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


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

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

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

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

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

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


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

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

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

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


\section{Unsupervised learning}

\subsection{EM and multiple alignment}

In this section, we consider two unsupervised learning problems. The
first problem is clustering. The second problem is concerned with
the situation where the objects can appear at different locations
and scales in the training images.

Because the active basis model is a generative model, the above two
unsupervised learning problems can still be treated in the maximum
likelihood framework, where the unknown categories or locations can
be incorporated into the active basis model as hidden variables. The
maximum likelihood computation can be accomplished by the EM
algorithm, where the E-step imputes the hidden variables, and the
M-step learns the active basis model given the hidden variables.

Unsupervised learning using active basis can also be considered more
specifically as a problem of multiple alignment. Here by multiple
alignment we mean the alignment of multiple images, which includes
but is not limited to the pairwise alignment of two images. For the
clustering problem, we want to separate the training images into
different clusters, so that images within the same cluster are well
aligned. In order to learn from training images where the objects
may appear at different locations, we want to identify the bounding
box for each image, so that the image patches within the bounding
boxes are well aligned.

The alignment can be scored by either the log-likelihood or the
active correlation. For image patches $\{\I_m, m = 1, ..., M\}$, we
use the notation ${\rm ALIGN}(\I_m, m = 1, ..., M)$ to denote the
total log-likelihood ratio maximized over $\B$, $\{\B_{m}\}$, and
$\Lambda$, or the total active correlation maximized over $\B$,
$\{\B_{m}\}$, and $\Theta$. The maximization is carried out by the
shared sketch algorithm in Subsection (\ref{sect:SS}), and ${\rm
ALIGN}(\I_m, m = 1, ..., M) = \sum_{m=1}^{M} \S2_m$, where the
$\S2_m$ scores are output by the algorithm. When we compute ${\rm
ALIGN}(\I_m, m = 1, ..., M)$ by the shared sketch algorithm, we also
learn the template $\B$, as well as the deformed template $\B_m$ for
each image $\I_m$.


\subsection{Clustering: EM and K-mean}

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

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

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

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

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

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

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

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


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

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

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

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


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

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

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

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

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

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

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

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

The shortcoming of Experiment 4 is that our method is less
successful when the object shapes of different categories are not
very different, for instance, our method cannot separate horses and
cows very well.

\subsection{Learning from non-aligned training images}

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

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

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

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

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

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


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

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



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

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


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


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


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

The shortcoming of Experiment 5 is that it is still not completely
unsupervised in that it requires the bounding box of one training
image. We believe this shortcoming can be overcome. We may crop many
different patches in the training images, and run the algorithm
using each patch as a starting point. Then we choose the best
multiple alignment based on the total of the $\S2_m$ scores. One
issue in this scheme is how to do it in an efficient manner.

\section{Learning moving template from motion sequence}

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

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

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

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

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

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

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


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

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


\begin{figure}[!htb]
\begin{center}
\begin{tabular}{c}
\psfig{figure=ex6_woman_score.eps,height=2in}
\end{tabular}
\end{center}
\caption{Experiment 6.2. Alignment scores at 10 different speeds at
the horizontal direction. } \label{fig:ex6_2_score}
\end{figure}

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

Figures (\ref{fig:ex6_2_score}) and (\ref{fig:ex6_2_sequence})
display results for another example.

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

\section{Composing multiple templates}

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


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

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

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

Given the two sub-templates, the recognition can be accomplished by
alternating the sum maps and max maps as illustrated by Figure
(\ref{fig:SumMaxMapsS3}). Because we have two sub-templates, we
obtain two $\S2$ maps, one for each sub-template. For each $\S2$
map, we apply a local maximization operator to obtain a $\M2$ map.
After that, we sum over the two $\M2$ maps to obtain the $\S3$ map.
This map scores the goodness of fit of the multi-basis.

The algorithm is described as follows.

\vspace{.5cm} \noindent{\bf Algorithm 3: Multi-basis algorithm}

\begin{enumerate}

\item
For $j = 1, 2$, compute $\S2(x, y, j)$ according to the recognition
algorithm of Subsection (\ref{sect:recognition}).

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

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

\item Let $\M3 = \max_{x, y} \S3(x, y)$. Suppose $(\hat{x},
\hat{y})$ achieves the maximum. For $j = 1, 2$, retrieve
$(\hat{x}_j, \hat{y}_j)$ so that
\begin{eqnarray}
\M2(\hat{x}+x_j, \hat{y}+y_j, j) = \S2(\hat{x}_j, \hat{y}_j, j),
\end{eqnarray}
in the local maximization operation (\ref{eq:M2}) in Step 2. Sketch
$j$-th sub-template as described in Step 5 of the recognition
algorithm of Subsection (\ref{sect:recognition}).

\end{enumerate}

\vspace{.5cm}

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

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

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


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

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

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

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

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

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


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

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

This experiment only serves as an illustration of multi-basis and
recursive sum-max maps. It also offers a simple example where we may
extract the sub-templates or parts based on the alignment score. We
are now working on learning multi-basis from non-aligned training
images. The method in Experiment 5 can be useful for extracting
sub-templates or highly alignable parts.

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


{\em Speed on current computer}. We ran our experiments on current
PC (Core 2 Duo 2.4GHz). In learning from the 27 regular bicycle
images, it took 3.6 seconds to filter the 27 images of $180 \times
180$, and it took 5 seconds to learn the 60 basis elements of the
bike. In analyzing the 330 $\times$ 496 image in Figure
(\ref{fig:ex7_1}.a), it took 2.1 seconds to compute $\S1$ maps, and
2 seconds for $\M1$ maps. The two $\S2$ maps took 0.4 and 0.7
seconds respectively. The two $\M2$ maps took 0.9 and 0.9 seconds
respectively.


\section{Discussion}

\indent{\em Active basis as and-or graph}. The active basis is a
simplest form of an and-or graph. The and-or graph is a useful
scheme for representing recursive compositional structures
\cite{Zhu_grammar, Geman_composition}. In an and-or graph, ``and''
means composition of parts, and ``or'' means variations of parts,
including the variations of their locations, scales, and
orientations. In active basis, ``and'' means composition of Gabor
wavelet elements, and ``or'' means variations in the locations and
orientations of the Gabor elements.

The and-or graph is a grammar that can be applied  recursively,
where parts can be further decomposed into sub-parts. A multi-basis
follows such a grammar. It is composed of multiple active bases, and
each active basis is composed of multiple Gabor elements.

The sum-max maps are natural hierarchical structures  for parsing an
image according to the and-or grammar. The sum maps score the
and-composition, and the max maps account for the or-variations.

{\em Bottom-up scoring and top-down sketching.} Our computation of
the sum-max maps in the recognition algorithm should not be
classified as a bottom-up process. For instance, in Experiment 7,
after a tandem bike is detected according to the scores in the $\S3$
map, we then go back to the $\M2$ maps to identify the actual
locations of the two wheels, which achieve local maxima of the $\S2$
maps. After that, we then go back to the $\M1$ maps to identify the
actual locations and orientations of the Gabor elements that make up
each of the two wheels. These locations and orientations achieve the
local maxima of the $\S1$ maps. This is a top-down sketching
process, which is triggered after the bottom-up scoring process
declares that an object is detected.

If we are to apply our method to real situations of object
recognition, we need a large number of sum maps and max maps,
because there are a large number of objects and their parts. It is
important to develop efficient scheduling scheme so that we do not
have to compute all these maps, and for those maps that we do need
to compute, we do not have to compute these maps in completion.



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


{\small
\begin{thebibliography}{9}


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


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

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

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

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

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

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

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

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

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

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


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


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

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


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

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

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

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

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

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


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


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

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

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

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


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

\end{thebibliography}
}

\end{document}
