\documentclass[10pt,journal,letterpaper,compsoc]{IEEEtran}
\usepackage{bbm}
\usepackage{amsfonts}
\usepackage[usenames]{color}
\usepackage{cite}
\usepackage{psfrag}
\usepackage{subfigure}
%\usepackage{url}
%\usepackage{stfloats}
\usepackage{amsmath}

\usepackage{times}
\usepackage{epsfig}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{cancel}
\usepackage{algorithm}
\usepackage{algorithmic}




\hyphenation{op-tical net-works semi-conduc-tor}



\graphicspath{{figures/}}
\begin{document}

%
% paper title
% can use linebreaks \\ within to get better formatting as desired
\title{Particle Filter for Unordered Observations}

\author{Longin~Jan~Latecki, %~\IEEEmembership{Member,~IEEE,}
        Nagesh~Adluru, %~\IEEEmembership{Student~Member,~IEEE}
        and~Xingwei~Yang%~\IEEEmembership{Student~Member,~IEEE}% <-this % stops a space
\IEEEcompsocitemizethanks{\IEEEcompsocthanksitem L. J. Latecki and X. Yang are with the Dept.
of Computer and Information Sciences, Temple University, Philadelphia.\protect\\
% note need leading \protect in front of \\ to get a newline within \thanks as
% \\ is fragile and will error, could use \hfil\break instead.
E-mails: latecki@temple.edu, xingwei.yang@temple.edu
\IEEEcompsocthanksitem Nagesh Adluru is with ?.}% <-this % stops a space
\thanks{}}





% The paper headers
\markboth{Journal of \LaTeX\ Class Files,~Vol.~6, No.~1, January~2007}%
{Shell \MakeLowercase{\textit{et al.}}: Bare Demo of IEEEtran.cls for Computer Society Journals}
% The only time the second header will appear is for the odd numbered pages
% after the title page when using the twoside option.
%
% *** Note that you probably will NOT want to include the author's ***
% *** name in the headers of peer review papers.                   ***
% You can use \ifCLASSOPTIONpeerreview for conditional compilation here if
% you desire.



% The publisher's ID mark at the bottom of the page is less important with
% Computer Society journal papers as those publications place the marks
% outside of the main text columns and, therefore, unlike regular IEEE
% journals, the available text space is not reduced by their presence.
% If you want to put a publisher's ID mark on the page you can do it like
% this:
%\IEEEpubid{0000--0000/00\$00.00~\copyright~2007 IEEE}
% or like this to get the Computer Society new two part style.
%\IEEEpubid{\makebox[\columnwidth]{\hfill 0000--0000/00/\$00.00~\copyright~2007 IEEE}%
%\hspace{\columnsep}\makebox[\columnwidth]{Published by the IEEE Computer Society\hfill}}
% Remember, if you use this you must call \IEEEpubidadjcol in the second
% column for its text to clear the IEEEpubid mark (Computer Society jorunal
% papers don't need this extra clearance.)




% for Computer Society papers, we must declare the abstract and index terms
% PRIOR to the title within the \IEEEcompsoctitleabstractindextext IEEEtran
% command as these need to go into the title area created by \maketitle.
\IEEEcompsoctitleabstractindextext{%
\begin{abstract}
Usually in the Particle Filter (PF) framework it is assumed that the observations arrive sequentially,
e.g., the observations are naturally ordered by their time stamps in the tracking scenario.
Based on this assumption, the posterior density over the corresponding hidden states is estimated.

In this paper we relax the assumption of having
ordered observations and extend the PF framework to simultaneously estimate the
posterior density while exploring different orders of observations.
This significantly broadens the scope of applications of PF methods.
We demonstrate that the proposed framework can be applied to solve a broad spectrum of challenging problems,
e.g., building jigsaw puzzles, detecting and localizing target objects in images,
and inferring a globally consistent map given a set of partial maps in the task of multi robot map merging.
The provided experimental results demonstrate that the proposed framework
yields comparable results to state-of-the-art algorithms for object detection and localization in images.
\end{abstract}

\begin{keywords}
Particle Filter, Object Detection and Recognition, Contour Grouping.
\end{keywords}}


% make the title area
\maketitle


% To allow for easy dual compilation without having to reenter the
% abstract/keywords data, the \IEEEcompsoctitleabstractindextext text will
% not be used in maketitle, but will appear (i.e., to be "transported")
% here as \IEEEdisplaynotcompsoctitleabstractindextext when compsoc mode
% is not selected <OR> if conference mode is selected - because compsoc
% conference papers position the abstract like regular (non-compsoc)
% papers do!
\IEEEdisplaynotcompsoctitleabstractindextext
% \IEEEdisplaynotcompsoctitleabstractindextext has no effect when using
% compsoc under a non-conference mode.


% For peer review papers, you can put extra information on the cover
% page as needed:
% \ifCLASSOPTIONpeerreview
% \begin{center} \bfseries EDICS Category: 3-BBND \end{center}
% \fi
%
% For peerreview papers, this IEEEtran command inserts a page break and
% creates the second title. It will be ignored for other modes.
%\IEEEpeerreviewmaketitle




\section{Introduction} \label{sec:Intro}

We illustrate the key ideas of the proposed approach with an example of a jigsaw puzzle,
where the goal is to estimate the poses (location and the orientation) of all puzzle pieces so that
they fit according to their shape constraints and form a meaningful image.
When a puzzle is created, one usually starts with an image on a paper or wooden board, e.g., see Fig.~\ref{fig:puzzle}(a).
Then the board is cut into jigsaw puzzle pieces.
In our scenario, the player is only presented the jigsaw puzzle pieces in a random arrangement
(e.g., see Fig.~\ref{fig:puzzle}(b)) and is not shown the original image.
Thus, the player has never seen the original image.
To be more precise, the player is given $m$ puzzle pieces described by a set of $m$ observations $Z=\{ z_1, \dots, z_m \}$.
Each observation $z_t$ describes piece $t$ and is given by a vector of features
that represents the shape and the partial image depicted on piece $t$.
The state of piece $t$, which describes its is pose, is represented by a random variable (RV) $X_t$.
Each value $x_t$ of $X_t$ represents 2D coordinates and the orientation of piece $t$.
Our goal is to determine the state vector $(x_1, \dots, x_m )$ that solves the given puzzle,
i.e., the relative poses of the puzzle pieces that yield the original image.
In the statistical framework such a vector of poses maximizes the posterior distribution $p(X_1, \dots, X_m | Z)$.
In other words, $(x_1, \dots, x_m )$ is a vector of most likely poses of the puzzle pieces given the observed features $Z$.
The posterior distribution is a function of how well adjacent pieces fit together with respect to their shape and their images,
which can be measured by similarity of their observed and derived features.
Of course, the goal is to build the original image, which however in unknown to the player.
Since the whole original image is not observed, $Z$ is a set of partial observations that describe local features of each puzzle piece.
Therefore, the posterior distribution, which we want to infer, is a function the set of partial observations $Z$,
i.e., we use the observations $Z$ to evaluate the likelihood of state vectors $(x_1, \dots, x_m )$.
Then, usually, the state vector that maximizes the posterior distribution is selected as the solution.
We observe that the initial poses are irrelevant, since the player is initially given a random configuration of the puzzle pieces.
We also observe that the posterior distribution $p(X_1, \dots, X_m | Z)$ usually is very complicated and has many local maxima.
This is particulary the case when the local information of the puzzle pieces is not very descriptive.
Therefore, it is often important to estimate the whole posterior distribution.

\begin{figure*}[bt]
\centering
\begin{tabular}{c c c c c c}
(a) & \includegraphics[width=0.18\linewidth]{puzzle2.eps} &
(b) & \includegraphics[width=0.3\linewidth]{puzzle1.eps} \vspace{0.3in} &
(c) & \includegraphics[width=0.28\linewidth]{puzzle3num.eps}
\end{tabular}
\vspace{-16pt}
\caption{\label{fig:puzzle}
The goal is to build the original image (a) given the jigsaw puzzle pieces (b).
The original image is not known, thus, it needs to be estimated given the observations shown in (b).
Even if the puzzle pieces were numbered (c), following their numbering may not be helpful in solving the puzzle.}
\end{figure*}


It is possible to apply the classical particle filter (PF) framework as stochastic optimization to solve this problem by utilizing
a fix order of observations $Z=\{ z_1, \dots, z_m \}$, e.g., by following a given numbering of the puzzle pieces.
However, by doing so, we would have selected an arbitrary order,
and the puzzle construction may fail because of the selected order.
Would we have selected a different order, the localization task could have been successful.
Moreover, the observations $Z=\{ z_1, \dots, z_m \}$ are given simultaneously at the same time.
Consequently, there is no reason to favor any particular order without utilizing this fact.

The classical particle filter (PF) framework has been developed for sequential state estimation like tracking or robot localization.
There, the observations arrive sequentially and are indexed by their time stamps.
In the proposed framework we extend PF to handle the situations where we have unordered set of observations that are given simultaneously
like in our example of jigsaw puzzle.
One of our key observations is the fact that it is possible to extend the importance sampling from the proposal distribution
so that the particle resampling automatically determines most informative orders of observations
allowing us to simultaneously consider different orders of observation in the process estimating the posterior distribution.
If our goal is maximizing the posterior distribution, we can view the order of observations that led to the global maximum as a most informative.
Of course, this order does not have to be unique, since it is possible that other orders of observations lead to the same global maximum.



\section{Basic Facts about Particle Filters}\label{sec:background}

In this section we summarize basic facts about Particle Filters (PFs).
They will be utilized in the next section when we introduce the proposed framework.

Our goal is to compute a posterior distribution $p(X_1, \dots, X_m \ | \ Z )$,
where $(X_1, \dots, X_m)$ is a vector of random variables (RVs)
and $Z=\{ z_1, \dots, z_m \}$ is a set of observations.
This will allow us to find value assignments $X_t=\widehat{x}_t$ for $t=1, \dots, m$
to RVs that maximize this posterior:
\begin{align}
& \widehat{x}_{1:m} = \underset{x_{1:m}}{\operatorname{argmax}} \ p(X_1=x_{1}, \dots, X_m=x_{m} \ | \ Z),
\label{eqG2}
\end{align}
where $x_{1:m}= (x_1, \dots, x_m) \in \mathcal{X}^m$ is a state space vector.
As it is commonly the case, we will abbreviate
$$p(X_1=x_{1}, \dots, X_m=x_{m} \ | \ Z)=p(x_{1:m} \ | \ Z).$$
We will achieve our goal by approximating the posterior distribution with a finite number of samples
in the framework of Bayesian Importance Sampling (BIS).
Since it is usually difficult to draw samples from the pdf $p(x_{1:m}|Z)$,
we will draw samples $x_{1:m}^{(i)} \sim q(x_{1:m}|Z)$ for $i=1, \dots, N$
from a proposal pdf $q$, form which samples are easily generated.
Then approximation to the density $p$ is given by
\begin{align}\label{eq:PFn1}
p(x_{1:m}|Z)& \approx \sum_{i=1}^N w^{(i)} \delta(x_{1:m} - x_{1:m}^{(i)}),
\end{align}
where $\delta$ is the Dirac delta function and
\begin{align}\label{eq:PFn2}
w^{(i)} = \frac{p(x_{1:m}^{(i)}|Z)}{q(x_{1:m}^{(i)}|Z)}
\end{align}
are normalized weights (so that they sum to one).

Since it is still hard to draw samples from $q$ due to high dimensionality of $x_{1:m}$,
Sequential Importance Sampling (SIS) is usually utilized.
In the classical PF approaches, samples are generated recursively
following the order of dimensions in the vector of RVs $X=( X_1, \dots, X_t )$:
\begin{align}\label{eq:sampling}
x_{t}^{(i)} \sim q_t(x|x_{1:t-1}, Z) = q_t(x|x_{1:t-1}, z_{1:t})
\end{align}
for $t = 1, \dots m$, and the particles are built sequentially $x_{1:t}^{(i)} =(x_{1:t-1}^{(i)}, x_{t}^{(i)})$ for $i=1, \dots, N$.
The subscript $t$ in $q_t$ indicates from which dimension of the state space the samples are generated,
i.e., from which RV.
We use this notation to stress that $q_t$ is a pdf from which we sample.
(For $t = 1$, $x_{1:t-1}=x_{1:0}$ denotes an empty vector, thus, we sample $x_{1}^{(i)} \sim q_1(x|Z)$.)
Since $q$ factorizes as
\begin{align}\label{eq:factor1}
& q(x_{1:m} | Z) = q(x_1 | Z) \prod_{t=2}^m q(x_t | x_{1:t-1}, Z),
\end{align}
we obtain that $x_{1:m}^{(i)} \sim q(x_{1:m} | Z)$.
In over words, by sampling recursively $x_{t}^{(i)}$ from each dimension $t$ according to \eqref{eq:sampling}
and appending the samples to a vector $x_{1:m}^{(i)}$ we obtain a sample from $q(x_{1:m} | Z)$

Since at a given iteration we have a partial sample vector $x_{1:t}^{(i)}$ for $t<m$,
we also need an evaluation procedure of this partial sample vector.
For this we observe that the weights can be recursively updated according to
\begin{align}\label{eq:PFn3}
w(x_{1:t}^{(i)}) = \frac{p(z_t|x_{1:t}^{(i)},z_{1:t-1}) p(x_{t}^{(i)}|x_{1:t-1}^{(i)})}
{q(x_{t}^{(i)}|x_{1:t-1}^{(i)},z_{1:t})} \ w(x_{1:t-1}^{(i)}) .
\end{align}
The key observation is that when $t=m$, the weight $w(x_{1:m}^{(i)})$ of particle $(i)$ is equal to $w^{(i)}$
(defined in \eqref{eq:PFn2}). The derivation of this fact can be found in \cite{ThrunBook05}.
Now we are able to formulate the SIS theorem in our setting: \\

\noindent {\bf Theorem 1. }
By sampling particles according to \eqref{eq:sampling} and weighting them according \eqref{eq:PFn3},
we obtain a set of weighted samples form $p(x_{1:m}|Z)$ once $t=m$.
Consequently, we can approximate $p(x_{1:m}|Z)$ with any precision
if the number of particles $N$ is sufficiently large. Thus, we can write
\begin{align}\label{eq:PFSIS}
p(x_{1:m}|Z)& \approx \sum_{i=1}^N w(x_{1:t}^{(i)}) \delta(x_{1:m} - x_{1:m}^{(i)}).
\end{align}


The equation (\ref{eq:PFn3}) can be simplified by making a common assumption that
$q(x_{t}^{(i)}|x_{1:t-1}^{(i)},z_{1:t}) = p(x_{t}^{(i)}|x_{1:t-1}^{(i)})$,
i.e., the proposal distribution depends only on the partial state vector and does not depends on the observations.
This assumption significantly simplifies the recursive weight update formula
\begin{align}\label{eq:PFn4}
w(x_{1:t}^{(i)}) = w(x_{1:t-1}^{(i)}) p(z_t|x_{1:t}^{(i)},z_{1:t-1}),
\end{align}
and, consequently, the samples are generated from
\begin{align}\label{eq:sampling2}
x_{t}^{(i)} \sim p_t(x|x_{1:t-1}^{(i)}),
\end{align}
i.e., the current observations $z_{1:t-1}$ have no direct influence on the samples.
Analog to \eqref{eq:sampling} we added the index $t$ in \eqref{eq:sampling2}
to indicate the dimension of the state space from which the samples are generated.
The influence of observations on the samples is introduced when particles are resampled according to their weights, which is a common step in PF algorithms,
since particle weights depend on observations.

In order to motivate the proposed approach, we derive now
the recursive weight update in \eqref{eq:PFn3} in a different way than
is usually reported in the literature, e.g., in \cite{ThrunBook05}:
\begin{align}\label{eq:PFn3.5}
&w(x_{1:t}^{(i)}) = \frac{p(x_{1:t}^{(i)}|Z)}{q(x_{1:t}^{(i)}|Z)} =
\frac{ p(x^{(i)}_{t}|x_{1:t-1}^{(i)},Z) \ p(x_{1:t-1}^{(i)}|Z) }
{ q(x^{(i)}_{t}|x_{1:t-1}^{(i)},Z) \ q(x_{1:t-1}^{(i)}|Z) } \nonumber \\
&=\frac{ p(x^{(i)}_{t}|x_{1:t-1}^{(i)},Z) }{ q(x^{(i)}_{t}|x_{1:t-1}^{(i)},Z) } \ w(x_{1:t-1}^{(i)}) \\
&= \frac{ p(Z|x_{1:t}^{(i)}) } { p(Z|x_{1:t-1}^{(i)})  }  \
\frac{ p(x^{(i)}_{t}|x_{1:t-1}^{(i)}) } {  q(x^{(i)}_{t}|x_{1:t-1}^{(i)},Z) } \ w(x_{1:t-1}^{(i)}) \label{eq:PFn3.6}
\end{align}
Eq.~\eqref{eq:PFn3.6} follows from \eqref{eq:PFn3.5} by
Bayes rule interchanging $x^{(i)}_{t}$ and $Z$ in $p(x^{(i)}_{t}|x_{1:t-1}^{(i)},Z)$.


Since both formulas \eqref{eq:PFn3} and \eqref{eq:PFn3.6} are equal to  $w(x_{1:t}^{(i)})$,
by setting \eqref{eq:PFn3} = \eqref{eq:PFn3.6},
we obtain
\begin{align}\label{eq:PFn3.8}
& p(z_t|x_{1:t}^{(i)},z_{1:t-1}) = \frac{p(Z | x_{1:t}^{(i)})}{p(Z |x_{1:t-1}^{(i)})}
\end{align}
Equation \eqref{eq:PFn3.8} provides an intuitive justification for the reclusive weight update in \eqref{eq:PFn4}.
The weight of a particle $(i)$ increases, if adding the new state $x_{t}^{(i)}$ to the state vector $x_{1:t-1}^{(i)}$
increases the probability of the observations $Z$,
i.e., if $p(Z | x_{1:t}^{(i)}) > p(Z |x_{1:t-1}^{(i)})$. \\


Now we outline a {\bf standard PF algorithm}.
For every time step $t=1, \dots, m$ and for every particle $i =
1, \dots, N$ execute the following three steps:

\noindent 1) \textbf{Importance sampling / proposal:}
Sample followers of particle $(i)$ according to \eqref{eq:sampling2} (a special case of \eqref{eq:sampling}) and set
$x_{1:t}^{(i)} = (x_{1:t-1}^{(i)}, x_{t}^{(i)})$.

\noindent 2) \textbf{Importance weighting / evaluation:}
An importance weight is assigned
to each particle $x_{1:t}^{(i)}$ according to \eqref{eq:PFn4} (a special case of \eqref{eq:PFn3}).

\noindent 3) \textbf{Resampling:}
Sample with replacement $N$ new particles form the current set of particles
$$ \{ x_{1:t}^{(i)} | i=1, \dots, N \}$$
according to their weights.
We obtain a set of new particles $x_{1:t}^{(i)}$ for $i=1, \dots, N$, all with weights of $1/N$.
This procedure is called Sampling Importance Resampling (SIR) \cite{ThrunBook05}.
It is an important part of any PF algorithm, since resampling prevents weight degeneration of particles.





A common assumption underlying the equations \eqref{eq:sampling} and \eqref{eq:PFn3} is that there exist two known functions
$f$ and $h$ such that
\begin{align}\label{eq:PFa1}
x_{t} &= f(x_{1:t-1}) + u_t\\
z_{t} &= h(x_{t}) + v_t \label{eq:PFa2}
\end{align}
where both $u_t$ and $v_t$ are mutually independent and identically distributed sequences with known probability density functions (pdfs). Often they are assumed to be Gaussians that model state prediction noise $u_t$
and observation noise $v_t$.
(Under the Markov assumption, \eqref{eq:PFa1} is simplified to $x_{t} = f(x_{t-1}) + u_t$,
i.e., that we can determine $x_{t}$ if we know the previous state $x_{t-1}$.)




\section{Unordered Observations}\label{sec:new}

In the proposed approach, the order of the observations is not predetermined,
in particular, we do not follow the order of indices of the observations in $Z$.
Our key idea is to extend the PF framework to examine all possible orders of observations
and to follow the most informative orders.
This way we are able to utilize the the most informative observations first.
Intuitively, it makes sense, for example,
if the first puzzle piece has shape very similar to many other puzzle pieces
and the second puzzle pieces has a very distinctive shape,
then our approach will first process the second puzzle piece, since it is more informative.


We stress that the standard SIS in Eq.~\ref{eq:sampling}
and particle evaluation in Eq.~\ref{eq:PFn3} utilize the sequential order of the RVs
reflected in the order of dimensions in the state space $(x_1, \dots, x_m )$.
In many applications, this order is determined naturally by the time stamp of the observations,
e.g., a single robot is collecting laser measurements at consecutive time points,
in which case $x_t$ denotes the robot pose at time $t$.
The goal of our work is to extend SIS to applications in which there is no natural order of observations like the case of multi robot map merging.

The key idea of the proposed approach is not to utilize the fix order of the dimensions,
but instead explore all possible orders of the dimensions $(x_{i_1}, \dots, x_{i_m})$ (or equivalently RVs)
so that the corresponding sequence of observations $Z=(z_{i_1}, \dots, z_{i_m})$ is most informative.
To achieve this we modify the first step of the PF algorithm so that
the importance sampling is performed for every dimension not yet represented by the current particle.
Thus, we sample follower states $x_{i_{t}}$ of particle $(x^{(i)}_{i_1}, \dots, x^{(i)}_{i_{t-1}})$
from each dimension that is not represented in $(x^{(i)}_{i_1}, \dots, x^{(i)}_{i_{t-1}})$.
For example, our puzzle in Fig.~\ref{fig:puzzle}(c) has the total of six puzzle pieces
and the current state vector of particle $(i)$ is
$(x^{(i)}_{i_{1}},  x^{(i)}_{i_{2}}) = (x^{(i)}_{5},  x^{(i)}_{3})$, where $i_{1}=5, i_{2}=3$,
meaning that we have already placed pieces 5 and 3,
Then in the next step we sample the follower states from dimensions representing the poses of the remaining puzzle pieces,
which are pieces 1, 2, 4, 6 in our example.



To simplify the notation we replace the double indexing of the state variables $x^{(i)}_{i_s}$
with a bijection (onto and one-to-one function)
$<\cdot>: \{1, \dots, m \} \to \{1, \dots, m \}$ and use the shorthand notation
$<1:t>$ to denote $(<1>, <2>, \dots, <t>)$.
For example, if $(i_{1},  i_{2}, i_{3})=(2,3,1)$, then
$$<1:3>=(<1>, <2>,<3>)=(i_{1},  i_{2}, i_{3})=(2,3,1)$$
and we denote $x^{(i)}_{<1:3>} = x^{(i)}_{i_{1:3}} = (x^{(i)}_{i_{1}},  x^{(i)}_{i_{2}}, x^{(i)}_{i_{3}})$.


In order to introduce the extended sampling method,
we need to revise the underlying assumption \eqref{eq:PFa1} to the existence of a function
\begin{align}\label{eq:PFa3}
x_{s} &= f(x_{<1:t-1>}) + u_{s,<1:t-1>}  \ \ \mathrm{for} \ \  s \in \mathcal{N} (x_{<1:t-1>}),
\end{align}
where $\mathcal{N} (x_{<1:t-1>})$ denotes the neighborhood of the set of states $x_{<1:t-1>}$.
The neighborhood is application dependent with a simplest example being
$\mathcal{N} (x_{1:t-1}) = \{ t \}$, in which case we obtain \eqref{eq:PFa1} as a special case of \ref{eq:PFa3}.

In our jigsaw puzzle example, the neighborhood of the pose states of the two puzzle pieces 3 and 5 can be composed of all the remaining pieces, e.g.,
$\mathcal{N} (x_{<1,2>}) = \mathcal{N} (x_{3,5}) = \{1, 2, 4, 5 \}$.
Then function $f$ allows to predict the poses of the puzzle pieces $s \in \mathcal{N} (x_{<1,2>})$.
Consequently, \eqref{eq:PFa3} means that we can estimate the relative pose of the remaing puzzle pieces
for a given configuration of already positioned pieces.





We stress that now the sequence of states $x_{<1:t-1>}$ visited before time $t$ is not a sequence of consecutive numbers
$(1, \dots, t-1)$ but any subsequence $(i_1, \dots, i_{t-1})$ formed by $t-1$ different numbers in
$\{1, \dots, m \}$.
Due to noise factor $u_{s,<1:t-1>}$ in \eqref{eq:PFa1}, each value $x_s$ is an estimate or prediction of a true unknown value.



\section{Particle Filter for Unordered Observations}\label{sec:PF}
We are ready now to precise formulate the proposed importance sampling.
At iteration $t$, the assumption \eqref{eq:PFa3} allows us to sample
\begin{align}\label{eq:sampling2.5}
x_{s}^{(i)} \sim p_s(x|x_{<1:t-1>}^{(i)}) \ \ \mathrm{for \ each} \ \  s \in \mathcal{N} (<1:t-1>),
\end{align}
where the index $s$ at the posterior pdf $p_s$ indicates that we sample from RV in dimension $s$.
Thus, we generate at least one sample for each dimension $s \in \mathcal{N} (<1:t-1>)$.

For example, if the neighborhood of $<1:t-1>$ is composed of all dimensions not present in $<1:t-1>$, i.e.,
\begin{align}\label{eq:neighb1}
\mathcal{N}(<1:t-1>) = \{1, \dots, m \} \setminus <1:t-1>
\end{align}
and we generate exactly one sample from each
dimension (i.e., RV) $s \in \mathcal{N} (<1:t-1>)$, then
at iteration $t$ particle $(i)$ has $m-t+1$ followers.
Each follower is a sample from a different dimension of the state vector (i.e., a different RV).

Let us assume that we have total of $m=5$ pieces in our jigsaw puzzle and after iteration 3,
we determined the poses of puzzle pieces $<1:3>=(3, 1, 5)$.
We then consider the two remaining pieces as followers of $<1:3>=(3, 1, 5)$ at iteration 4.
The poses of pieces 2 and 4 are determined with reference to the poses of pieces $3, 1, 5$,
i.e., by sampling
$x_{2}^{(i)} \sim p_2(x|x_{<1:3>}^{(i)})$ and $x_{4}^{(i)} \sim p_4(x|x_{<1:3>}^{(i)})$.
Of course, this process is repeated for each particle (i) for $i=1, \dots, N$,
where $N$ is the number of particles.



In contrast, in the standard application of rule \eqref{eq:sampling2},
at each iteration $t$ particle $(i)$ has one follower.
Even when sometimes each particle $(i)$ has many followers, all followers are in the same dimension (samples form the same RV),
which means that we only determine possible locations of say puzzle piece 2
and do not consider locations of puzzle piece 4 for particle (i),
since a strict order of the state dimensions is followed in the classical setting.

We do not make any Markov assumption in \eqref{eq:sampling2.5},
i.e., the new state $x^{(i)}_{s}$ is dependent on all previous states $x_{<1:t-1>}^{(i)}$ for each particle $(i)$.




We derive now the recursive weight update formula for the unordered set of static observations $Z$.
This derivation is analog to our derivation in \eqref{eq:PFn3.6}.
For every $t$ from 2 to $m$, we have
\begin{align}\label{EQ:pf3}
&w(x_{<1:t>}^{(i)})
=\frac{ p(x_{<1:t>}^{(i)}  \ | \ Z) }{ q(x_{<1:t>}^{(i)}  \ | \ Z) } \nonumber \\
&=\frac{ p(x^{(i)}_{<t>}|x_{<1:t-1>}^{(i)},Z) \ p(x_{<1:t-1>}^{(i)}|Z) }
{ q(x^{(i)}_{<t>}|x_{<1:t-1>}^{(i)},Z) \ q(x_{<1:t-1>}^{(i)}|Z) } \nonumber \\
&=\frac{ p(x^{(i)}_{<t>}|x_{<1:t-1>}^{(i)},Z) }{ q(x^{(i)}_{<t>}|x_{<1:t-1>}^{(i)},Z) } \ w(x_{<1:t-1>}^{(i)}) \nonumber \\
&= \frac{ p(Z|x_{<1:t>}^{(i)}) p(x^{(i)}_{<t>}|x_{<1:t-1>}^{(i)}) }
{ p(Z|x_{<1:t-1>}^{(i)}) q(x^{(i)}_{<t>}|x_{<1:t-1>}^{(i)},Z) } \ w(x_{<1:t-1>}^{(i)})
\end{align}
The only difference to \eqref{eq:PFn3.6} is that we replace the fix order of states $1:t-1$ with a dynamic order $<1:t-1>$
that may be different for each particle $(i)$.
(In particular, to obtain the last equation, we apply Bayes rule to decompose
$p(x^{(i)}_{<t>}|x_{<1:t-1>}^{(i)},Z)$ that interchanges $x^{(i)}_{<t>}$ and $Z$.)

As it is often the case in PF applications, we assume that
$q(x^{(i)}_{<t>}|x_{<1:t-1>}^{(i)},Z) = p(x^{(i)}_{<t>}|x_{<1:t-1>}^{(i)})$.
Using this simple exploration based proposal the weight recursion in
\eqref{EQ:pf3} becomes:
\begin{align}\label{eqn:recur-wgt}
&w(x_{<1:t>}^{(i)})=
&=w(x_{<1:t-1>}^{(i)}) \ \frac{p(Z|x_{<1:t>}^{(i)})}{p(Z|x_{<1:t-1>}^{(i)})}
\end{align}

Our next derivation step differs fundamentally from the standard PF framework.
By recursive substitution of weights in \eqref{eqn:recur-wgt},
i.e., by applying \eqref{eqn:recur-wgt} to $w(x_{<1:t-1>}^{(i)}), w(x_{<1:t-2>}^{(i)}), \dots, w(x_{<1:2>}^{(i)})$,
we obtain
\begin{align}\label{eqn:ind-wgt}
&w(x_{<1:t>}^{(i)})
=w(x_{<1:t-2>}^{(i)}) \ \frac{\cancel{p(Z|x_{<1:t-1>}^{(i)})}}{p(Z|x_{<1:t-2>}^{(i)})} \ \frac{p(Z|x_{<1:t>}^{(i)})}{\cancel{p(Z|x_{<1:t-1>}^{(i)})}} \nonumber \\
&= \dots =w(x_{<1>}^{(i)}) \ \frac{p(Z|x_{<1:t>}^{(i)})}{p(Z|x_{<1>}^{(i)})}
\end{align}
Under the assumption that all particles have the same initial weight $w(x_{<1>}^{(i)})$
and the same initial observation probability $p(Z|x_{<1>}^{(i)})$ for $i = 1, \dots, N$, we obtain
\begin{align}\label{eqn:partweight}
&w(x_{<1:t>}^{(i)})= p(Z|x_{<1:t>}^{(i)})
\end{align}
Since for $t=m$ we have $w(x_{<1:m>}^{(i)})=w(x_{1:m}^{(i)})$,
we obtain that the weights computed by the recursive formulas \eqref{eqn:recur-wgt} -
\eqref{eqn:partweight} are equal to the weights in Eq.~\eqref{eq:PFn2}.

For comparison, the corresponding weight update in the standard PF
framework \cite{ThrunBook05} is
\begin{align}\label{eq:PFstandard}
w(x_{1:t}^{(i)})&=
w(x_{1:t-1}^{(i)}) \ p(z_t | x_{1:t-1}^{(i)}, x_t),
\end{align}
where $z_t$ denotes the new observations obtained at time $t$.

A distinctive feature of our derivations is the fact that we use the whole set of observations $Z$.
In particular, this means that we do not decompose $Z$ into partial sequences following some order of observations.
Consequently, our weight update formula \eqref{eqn:partweight} evaluates how likely
the whole set of observations $Z$ is conditioned on the current state vector $x_{<1:t>}^{(i)}$.


We outline now the proposed {\bf PF for unordered observations (PFUO)} algorithm.
For every time step $t=1, \dots, m$ and for every particle $i =
1, \dots, N$ execute the following three steps:

\noindent 1) \textbf{Importance sampling / proposal:}
Sample followers $x_{s}^{(i)}$ of particle $(i)$ from each dimension $s \in \mathcal{N} (<1:t-1>)$ according to \eqref{eq:sampling2.5},
which we repeat here for completeness,
\begin{align}
x_{s}^{(i)} \sim p_s(x|x_{<1:t-1>}^{(i)})
\end{align}
and set $x_{<1:t-1>,s}^{(i)} = (x_{<1:t-1>}^{(i)}, x_{s}^{(i)})$.
%We can think of each extended particle $x_{<1:t-1>,s}^{(i)}$ as an instance of particle $x_{<1:t>}^{(i)}$
%obtained by appending the state value $x_{<t>}^{(i)}$ of dimension $<t>=s$ to the state vector.
We need the second index $s$, since in our framework,
particle $x_{<1:t-1>}^{(i)}$ usually has more than one follower and $s$ indexes the followers of $x_{<1:t-1>}^{(i)}$.
This means that the single particle $x_{<1:t-1>}^{(i)}$ is multiplied and extended to several follower particles
$x_{<1:t-1>,s}^{(i)}$.
If we generate one sample for each dimension $s \in \mathcal{N} (<1:t-1>)$, then
we the number of followers $x_{<1:t-1>,s}^{(i)}$ is equal to the number of dimensions in $\mathcal{N} (<1:t-1>)$.

%
%In the first iteration ($t=1$) we generate $m$ samples
%\begin{align}\label{eq:sampling2.7}
%x_{<1>}^{(i)}=x_{s}^{(i)} \sim p_s(x) \ \ \mathrm{for} \ \  s \in \{1, \dots, m \}.
%\end{align}
%As is the case in the standard PF framework, often the initial distribution $p_s(x)$ is assumed to be uniform.

\noindent 2) \textbf{Importance weighting/evaluation:}
An individual importance weight is assigned
to each particle follower $x_{<1:t-1>,s}^{(i)}$ according to Eq.~\ref{eqn:partweight},
where we substitute  $x_{<1:t-1>,s}^{(i)}$ for $x_{<1:t>}^{(i)}$, i.e.,
\begin{align}\label{eqn:partweightv2}
&w(x_{<1:t-1>,s}^{(i)})= p(Z|x_{<1:t-1>,s}^{(i)})
\end{align}

\noindent 3) \textbf{Resampling:}
Sample with replacement $N$ new particles form the current set of particles
\begin{align}  \label{eqn:set1}
& \{ x_{<1:t-1>,s}^{(i)} |  \ i=1, \dots, N, \ s \in \mathcal{N} (<1:t-1>^{(i)}) \}.
\end{align}
according to the weights.
We obtain a set of new particles $x_{<1:t>}^{(i)}$ for $i=1, \dots, N$, all with weights of $1/N$.
This is the standard Sampling Importance Resampling (SIR) step \cite{ThrunBook05} as in the classical PF framework,
but the set of particles that is resampled is different.
We recall that the index of dimensions $<1:t-1>$ may be different for each particle $(i)$,
which is why we use $<1:t-1>^{(i)}$ in \eqref{eqn:set1}. \\



\noindent \textbf{Theorem 2.} \
The PFUO computes an approximation of Eq.~\ref{eq:PFn1},
i.e., for $t=m$ the particles $x_{<1:t>}^{(i)}$ with $i=1, \dots, N$
provide an approximation to the posterior $p(x_{1:m}|Z)$ for sufficiently large $N$, i.e.,
\begin{align}\label{eq:PFn1.1}
p(x_{1:m}|Z)& \approx \sum_{i=1}^N w(x_{1:m}^{(i)}) \delta(x_{1:m} - x_{1:m}^{(i)}).
\end{align}
Before we prove this theorem, we need to make an important remark.
Since $t=m$, all dimensions of considered RVs are present in each particle $x_{<1:t>}^{(i)}$.
Therefore, we can simplify the notation and write $x_{<1:t>}^{(i)}= x_{1:m}^{(i)}$.
We stress again that the particular order of dimensions $<1:t>$ is extremely important in our approach
in order to derive the state vector $x_{<1:t>}^{(i)}$ of particle $(i)$.
However, once we have this state vector,
we can simply present the state vector in the original order of dimensions.\\

\noindent  \textbf{Proof of Theorem 2.} \
Our proof is based on Eq.~\eqref{eq:PFn1.1} with weights defined in Eq.~\eqref{eq:PFn2}.
We observe that the weights computed by the recursive formulas
\eqref{eqn:recur-wgt} - \eqref{eqn:partweight} are equal to the weights in Eq.~\eqref{eq:PFn2},
since for $t=m$ we have $w(x_{<1:m>}^{(i)})=w(x_{1:m}^{(i)})$.

The Sampling Importance Resampling (SIR) replaces weighted particles
with $N$ particles with the weight equal to $1/N$, which provides an approximation to the same target pdf.
This proves the theorem. \\





The fact that we can consider more than one follower of each particle and
reduce the number of followers by resampling is known in the PF literature and
is referred to as prior boosting \cite{gordon93novel,carpenter99building}.
It is used to capture multi-modal likelihood regions.
We stress that the resampling in our framework plays an additional and a very crucial role.
It selects the the most informative random variables (i.e., state space dimensions) as followers of particles.
Since the weights of $x_{<1:t-1>,s}^{(i)}$ are determined by the observations $Z$,
and the resampling uses the weights to selects new particles $x_{<1:t>}^{(i)}$
by adding not yet considered dimensions,
the resampling determines the order of RVs, i.e., the bijection $<t>$ for $t=1, \dots m$.
Consequently, the order of RVs is heavily determined by observations $Z$,
and this order may be different for each particle $(i)$.
This is in strong contrast to the classical PF, where observations $Z$ have no influence on the order of RVs,
which is fixed.


In many applications, the weight \eqref{eqn:partweight} of each particle is based on the evaluation
of the similarity of the predicted observations $h(x_{<1:t-1>,s}^{(i)})$,
i.e., the observations derived from the current state vector $x_{<1:t-1>,s}^{(i)}$,
to the real observations in set $Z$.
If we adopt the standard PF assumption \eqref{eq:PFa2} that underlies particle weight computation,
we can only relate the predicted observations $h(x_{<1:t-1>,s}^{(i)})$
to the corresponding real observations $z_{<1:t-1>,s}$.
Then we utilize a special case of \eqref{eqn:partweightv2} to compute the particle weights
\begin{align}\label{eqn:partweight2}
&w(x_{<1:t-1>,s}^{(i)})= p(z_{<1:t-1>,s} | x_{<1:t-1>,s}^{(i)}).
\end{align}
This means that we only use the subset of observations that corresponds to the current dimensions of the state vector to evaluate particle $(i)$.
For instance, if the current state vector of particle $(i)$ is $x_{<1>,s}^{(i)} = x_{5, 3}^{(i)}$,
meaning that we have positioned the puzzle piece $<1>=5$ and the piece $s=3$
as shown in our puzzle example in Fig.~\ref{fig:puzzle}(c),
then the weight of particle $(i)$ depends on how likely the corresponding observations
$z_{<1>,s} = z_{5,3} =(z_5, z_3)$ are conditioned on the state vector $x_{<1>,s}^{(i)} = (x_5^{(i)}, x_3^{(i)})$.
In our jigsaw puzzle example, this means that we evaluate the likelihood of
the appearance features $(z_5, z_3)$ of puzzle pieces 5 and 3 conditioned on their poses $(x_5^{(i)}, x_3^{(i)})$.
This can be computed by evaluating the fit of the puzzle pieces positioned according to the poses $(x_5^{(i)}, x_3^{(i)})$.
Since the pieces 5 and 3 match well, the particle shown in Fig.~\ref{fig:puzzle}(c) would obtain a high weight.

However, in our framework it is also possible to generalize the assumption \eqref{eq:PFa2}
so that we always use the whole set of observations $Z$ to compute the particle weights,
i.e., we directly use formula \eqref{eqn:partweight}.
This may increase the discriminative power of the particle evaluation.
For instance, if we use the appearance of all six puzzle pieces $Z=\{z_1, z_2, z_3, z_4, z_5, z_6 \}$
to compute the weight of particle $x_{<1>,s}^{(i)} = (x_5^{(i)}, x_3^{(i)})$,
we can additionally consider how distinctive the appearance features $(z_5, z_3)$ are in the context of all feature vectors in $Z$.
In this case, the value of the weight depends on two factors,
\begin{itemize}
\item how descriptive a given vector of observations $z_{<1:t-1>,s}$ is in the context of all observations $Z$ and
\item how good the prediction $h(x_{<1:t-1>,s}^{(i)})$ of observations $z_{<1:t-1>,s}$ is.
\end{itemize}




\section{Looping Particle Filter for Unordered Observations}\label{sec:PF}

In the PF algorithm introduced in Section \ref{sec:PF}, the number of iterations $t$ is equal to the dimensionality $m$ of the state space.
In our jigsaw puzzle example with six puzzle pieces, this means that $t=6$.
Consequently a large number of particles $N$ is needed to well approximate the posterior distribution.
According to Theorem 2, such a number of particles always exist, but it may be hard in practice to set it properly.
Therefore, we propose in this section a modification of the importance sampling / proposal step that makes number of iterations
independent from the dimensionality of the state space.


In order to obtain a second version of the proposed algorithm, called {\bf Looping PF for unordered observations (L-PFUO)} algorithm,
we only need to slightly modify the neighborhood in the importance sampling / proposal step in \eqref{eq:sampling2.5}.
For iteration steps $t \leq m$ we keep the user defined neighborhood structure, e.g., as defined in \eqref{eq:neighb1},
but for $t >m$ we set the neighborhood structure to be
\begin{align}\label{eq:neighb2}
\mathcal{N}(<1:t-1>) = \{1, \dots, m \}.
\end{align}
We observe that for $t >m$ our L-PFUO becomes a variant of Gibbs sampler,
since we loop over the state space and replace values in one dimension by sampling from the conditional proposal distribution
(with values of the other dimensions being fixed).
However, our PF steps (2) and (3) are very different from the existing Gibbs sampler approaches.
While there exist Gibbs sampler approaches that weight samples \cite{}, none of them performs resampling as in our step (3).

In L-PFUO, the length of each state vector is equal to $t$ if $t <m$ and is equal to $m$ for every iteration $t \geq m$.
Thus, L-PFUO can be viewed as a mixture of PF and Gibbs sampler approaches.
For $t <m$ we follow the PF framework but with unordered observations
and incrementally build the state vector until we reach its full length $m$.
This makes L-PFUO significantly more efficient than Gibbs sampler,
which must be initialized with the full length state vectors and therefore requires a lengthy burn in phase
(i.e., the first $N$ samples are discarded, where $N$ can be very large).
Once $t = m$, L-PFUO switches to looping over the state space.
L-PFUO loops following the order determined by the resampling according to the particle weights,
which allows us to generate more samples from the most informative dimensions of the state space.
In contrast, most Gibbs samplers loop in a predefined deterministic order or in a completely random order \cite{}.

As we stated at the beginning of this section, L-PFUO algorithm allows us to estimate the posterior distribution
with a significantly smaller number of particles than PFUO.
This means that we trade the smaller number of particles $N$ for the larger number or iterations $t$.
While PFUO has a clearly defined stop condition $t=m$, we need to specify the stop condition for L-PFUO algorithm.
We simply iterate L-PFUO until the posterior distribution stabilizes,
i.e., the change in the posterior distribution is below a specified threshold.
This simple stop condition indicates the convergence of L-PFUO algorithm to a stable solution.
In contrast, there does not seem to exist any condition on the number of particles $N$ that would indicate reaching a stable solution of PFUO.





\section{Multi Robot Map Merging}

In this section, we describe an illustrative application of the proposed framework to the problem of merging partial maps acquired by multiple robots. We use map merging as an illustrative instance of the jigsaw puzzle problem where the individual puzzle pieces are the local maps which have to be put together into a full global map. The goal becomes to estimate the poses of the individual maps such that all the individual maps align into a consistent global map.

\subsection{System Description}
We have a team of $m$ robots that built $m$ local maps, $Z=\{ z_1, \dots, z_m \}$.
Each map is represented as a 2D occupancy grid with three types of cells representing free space, obstacles, and unknown space.
Some maps have partial overlap, i.e., they have common parts, while some do not have any overlap.
%which could be obtained by using standard single robot SLAM techniques (\cite{ThrunBook05}). 
Let $X=\{x_1,\dots,x_m\}$ be the poses of the local maps in a global reference frame that need to be estimated. 
Without getting into full details we explain the two most important components needed in importance weighting and proposal distribution viz. $p(Z|x_{<1:t-1>,s}^{(i)})$ and $p_s(x|x_{<1:t-1>}^{(i)})$. 
More details on the system can be found in \cite{adluru_et_al:icpr08}.
An example result is shown in Fig.~\ref{fig:vslam-results}.

\begin{figure*}[htb]
\begin{center}
\begin{tabular}{cc}
\includegraphics[width=0.17\textwidth]{figures/3robotsadaptive.eps}&
\includegraphics[width=0.26\textwidth]{figures/5robotsadaptive.eps}\\
\includegraphics[width=0.34\textwidth]{figures/6robotsadaptive.eps}&
\includegraphics[width=0.34\textwidth]{figures/10robotsadaptive.eps}
\end{tabular}
\end{center}
\caption{Merged maps of $3,5,6$ and $10$ robots. The map indices are
placed at the local origins of the individual maps. The sequence of
the maps merged is shown in the respective titles.}
%\label{fig:vslam-corners}
\label{fig:vslam-results}
\end{figure*}

\subsection{Importance weights}\label{sec:scan-match}
The importance weights depend on how well the partial maps fit together.
We denote with $M_{<1:t-1>,s}^{(i)}$ a partially assembled map of particle $(i)$ obtained by putting the local maps
$z_{<1:t-1>,s}$ at poses $x^{(i)}_{<1:t-1>,s}$ in the global coordinate system.
We define
\begin{equation}
 p(Z|x_{<1:t-1>,s}^{(i)})=\Psi(M_{<1:t-1>,s}) \cdot \gamma(M_{<1:t-1>,s})
\end{equation}
where $\Psi$ is computed following the description in \cite{birkmerge06}. 
The basic idea is to treat the images as arrays of pixels (occupancy grid) and reward matching
pixels based on their values and relative positions. The arrays are
scanned for each class of pixels (obstacles, free space, and unknown). 
For each pixel, the distance to
its closest similarly valued pixel in the other image is computed
and its score is taken to be proportional to a Gaussian transform of
this distance. Fortunately the distances can be computed in linear
time in the number of pixels in the images using "distance-maps" \cite{birkmerge06}.

$\gamma$ is designed in such a way that its magnitude is
made proportional to the merge inconsistency. The inconsistency of a merge is defined as the mismatch in the
perception of the robot between two differently calculated
positions. The mismatch in robot
perception is based on the number of disagreeing pixels in the two
images. Thus its computation involves similar steps as those for
$\Psi$.  The larger the number of matching pixels ($c_1=c_2$), the
larger $\Psi$ is.  In contrast to this, the smaller the number of
mismatched pixels ($c_1 \neq c_2$), the larger $\gamma$ is.

<Why do we need $\gamma$?>


\subsection{Proposal distribution}\label{sec:motion-prop}
The discrete probability distribution function (pdf), $p_s(x|x_{<1:t-1>}^{(i)})$, is simulated using ``structure registration`` among local maps. The initial distribution is uniform and can pick any of the local maps, which becomes the reference coordinate (global) frame in which the puzzle is completed. The poses of rest of the puzzles are determined in reference to this frame. Pose of each frame can be imagined at the center of each of the puzzle.

Similar structures among local maps are extracted and registered using closed form solutions like
in \cite{horn88orthonormal,horn87quaternion,arun87svd}. The
correspondences are obtained by shape matching between similar
structures. There could be several similar structures between local maps
$z_t$ and $x_u$. As a consequence, our
proposal distribution is usually multi-modal with peaks around the poses
predicted using structure registration. 
%Fig. \ref{fig:virt-motion} shows a sample pose update process with two modes. 
We note that each possible pose update is actually a
transformation of the local map, $z_s$, into the global frame of $M_{<1:t-1>,s}$.

%\begin{figure*}[htb]
%\begin{center}
%\begin{tabular}{ccc}
%\includegraphics[width=0.3\textwidth]{figures/virt_motion11.eps}&
%\includegraphics[width=0.22\textwidth]{figures/virt_motion2.eps}&
%\includegraphics[width=0.22\textwidth]{figures/virt_motion3.eps}
%\end{tabular}
%\end{center}
%\caption{Left: The pose of a local map, $Z[\varphi_{t-1}]$ in the global frame of $m_{t-1}$ is shown
%with an arrow inside the circle. The new local map, $Z[\varphi_t]$
%is shown in its local frame. The dotted red lines are the
%trajectories of the individual robots. Middle: One possible
%structure registration updates the pose of $Z[\varphi_t]$. Right: Another
%possible structure registration leads to a different pose update. The blue line indicates that both the puzzles are in the same frame.} \label{fig:virt-motion}
%\end{figure*}



\subsection{Experimental Results}
We applied our technique to the data collected by
National Institute of Standards and Technology (NIST) in a "maze"
type of environment\footnote{We would like to thank Raj Madhavan
(NIST) for providing the dataset.}. We are provided only with the
local maps without any initial estimate of the relative poses. For
feature extraction we fit lines to the data points and determine
corners. For typical indoor local maps, corners are sufficiently
distinct feature descriptors to represent the environment structure.
On average there were about $10$ corners in a
local map and there were about $5$ matching corner pairs between two
local maps.
Fig. \ref{fig:vslam-results} shows maps of best particles at iterations $3,5,6$ and $10$,
which corresponds to global maps merged from local maps of $3,5,6$ and $10$ robots.
As can be seen we could successfully merge maps of $10$ robots.
As virtual trajectory we list the state dimensions, which are the indices of the robots, used to build the global maps at corresponding iterations.
We observe that at different iterations the order is different,
which nicely illustrates the fact that the proposed framework explores different orders of particles.




\section{Object Detection and Recognition by Groping Contour Fragments}

\subsection{System Description}
In this section we present an application of the proposed framework to object detection and recognition in images.
Two detection examples are illustrated in Fig~\ref{fig:ParticleEx}.
The white curves show edge fragments obtained by edge detection.
The contour segments that form detected horses in these images are marked with colors,
where the colors indicate the order of detecting and localizing these segments.
We observe that the two orders are different, which nicely illustrates the advantages of the proposed framework,
e.g., the best particle for the image shown in (b) begins with the state dimension representing the position of the horse head,
which makes sense intuitively, since the head contour can be clearly identified in (b).
The situation is different in (a), where the horse head is positioned in the cluttered background.
The first dimension of the best particle in (a) represents the the state of a front leg.

\begin{figure}[tbh]
\centering
\begin{tabular}{c c}
    \includegraphics[width=3.8cm]{ParticleXingwei1.eps}  & \includegraphics[width=3.8cm]{ParticleXingwei2.eps}\\
    (a) & (b)
    \end{tabular}
\caption{Examples of two different inferred orders of detected contour parts.
Colors represent the order, which is 1=red, 2=cyan, 3=blue, 4=green, 5=yellow, and 6=black. }
\label{fig:ParticleEx}
\end{figure}

Our detection is model driven, where a model is defined by the contour of an example object of a given class
that can be hand drawn or extracted from an example image, e.g., see the apple logo in Fig.~\ref{fig:bundle.model}.
Since edge fragments in real images may be broken and may have missing parts,
we decompose our model contour into possibly overlapping model contour segments $S = \{s_1, \dots, s_k \}$,
which are grouped into part bundles.
An example bundle decomposition is shown in Fig.~\ref{fig:bundle.model}.
The main constraint for the bundle design is to ensure that a rough shape sketch obtained by selecting one part form each bundle
still resembles the model contour.
A cognitive motivation behind our bundle decomposition scheme
is that an object can be recognized even if some parts of it are
missing, as can be observed in Fig. \ref{fig:PartMiss}. There are
several reasons why parts of objects can be missing in real images:
occlusion, missing or weak edge information, and failures of edge detection algorithms.
The selection of parts and their grouping into bundles was designed manually.
We have one part bundle model per shape class.
%but parts in each bundle may come for more than one object in a given shape class.


\begin{figure}[tbh]
\begin{center}
\includegraphics[width=.7\linewidth]{Applelogo_Bundle_nummark.eps}
\caption{The contour model of the apple and the corresponding part
bundles. The contour is shown in the center.
The 11 contour fragments are decomposed into four part-bundles.}
\label{fig:bundle.model}
\vspace{-5mm}
\end{center}
\end{figure}



Formally, $\mathcal{B} = \{B_k\}_{k=1}^{m}$, where $B_k \subset S$ and $m\leq k$,
is a part bundle decomposition of $S$ if and only if
$\bigcup \mathcal{B}=S$ and $ B_i\bigcap B_j=\emptyset$ for $i,j=1,...,m$ and $i \neq j$.
For example, the part bundle decomposition of the apple logo in Fig.~\ref{fig:bundle.model}
is composed of 11 contour segments grouped into four bundles
i.e.,$k=11$ and $m=4$. The four bundles are represented by four RVs.
In general, if a model is composed of $m$ part bundles, our goal is to estimate the states of $m$ RVs $X_1, \dots, X_m$.
The values of each RV $X_t$ range over a 4D parameter space
$x_t=(s_t, h_t, v_t, \alpha_t)$, where $s_t$ is an index of a selected segment from the part bundle $B_t$,
$(h_t, v_t)$ describe the horizontal and vertical coordinates of the centroid of this segment on the image,
and $\alpha_t$ describes its orientation.
(We assume that the scale is known, but it is also possible to estimate it in our framework).

Our shape model $\mathcal{M}$ of a given object class is composed of the part bundle $\mathcal{B}$
and a set of contours of positive instances of objects in this class.
The instances represent shape constraints of the object class.
We measure the similarity of deformed parts on the image (due to rotation $\alpha_t$) to the set of the instance contours
in order to ensure that the deformed model fragments still resemble the shape of object in this class.


\begin{figure}[tbh]
\begin{center}
\includegraphics[width=0.8\linewidth]{show_part_missing.eps}
%\vspace{-3mm}
\caption{Parts of the objects are missing both due to missing edges
and due to broken edge links.}
\vspace{-7mm}
\label{fig:PartMiss}
\end{center}
\end{figure}


The two phases of the proposed algorithm L-PFUO have a natural interpretation in this framework.
Given a particle $x_{1:t-1}^{(i)}$,
at steps $t=2, \dots, m$ we construct a rough shape sketch of the model shape in the edge image.
Then for $t>m$ we loop over the RVs and refine the placement of the model segments on the image
with possibly selecting different segments form the same bundle.
Intuitively this means that the particles first trace a rough shape before adjusting the shape details.
Of course, we can achieve the same effect with our first algorithm PFUO if we simply have a larger number of particles.

\subsection{Importance Weights}
The weight of each particle is computed based on the shape similarity of the contour constructed from the selected model segments
and the edges in the image. We compute it as a Gaussian of the Oriented chamfer matching (OCD)~\cite{shotton_1}.
However, we also need to consider the deformation of the model introduced by the rotation of the model segments on the image.
We measure this deformation as similarity to known positive model instances.
To be more precise, in this application we have a static set of observations
$Z=(I,\mathcal{M})$, where $I$ is a given edge image and $\mathcal{M}$ is a model of a given shape class,
since both $I$ and $\mathcal{M}$ are known and fixed.
Since and $I$ and $\mathcal{M}$ can be viewed as independent conditioned on $x_{<1:t-1>,s}^{(i)}$
we obtain:
\begin{equation}
p(Z|x_{<1:t-1>,s}^{(i)})=p(I|x_{<1:t-1>,s}^{(i)}) p(\mathcal{M}|x_{<1:t-1>,s}^{(i)}). \label{eq:ci1}
\end{equation}
The first factor defines the probability that the model parts positioned on the image according to the parameters of particle $x_{<1:t-1>,s}^{(i)}$
led to the edges around them in the edge map $I$ and we define it as
\begin{align}
& p(I|x_{<1:t-1>,a}^{(i)})=\exp(-\beta \cdot OCD_{I}(x_{<1:t-1>,s}^{(i)})),
\label{eq: unary}
\end{align}
where $OCD_{I}$ returns the Oriented Chamfer distance.
Consequently, $OCD_{I}$ measures how well the constructed partial model matches to the edge map $I$.

The second factor in \eqref{eq:ci1} allows us to ensure that
the rotated model parts on the image still resemble objects in the given shape class.
It can be expressed as similarity of the model parts on the image (deformed by rotations $\alpha_t$)
to the set of the instance contours, which is also measured with OCD.

\subsection{Proposal Distribution}
In order to execute the derived PF algorithm, it remains to define the proposal distribution.
The initial proposal distribution $p_s(x)$ is simply the probability of finding a model part from part bundle $B_s$ at each location in the image $I$.
It is obtained as a Gaussian of the $OCM_{I}$ computed in the sliding window fashion at every location in $I$.
In our application, the conditional proposal distribution $p_s(x|x_{<1:t-1>}^{(i)})$
is the probability of finding a model part from part bundle $B_s$ around the locations
of the previously positioned parts according to the parameters of particle $x_{<1:t-1>}^{(i)}$,
i.e., the location of a model part from part bundle $B_s$ is constrained by the already positioned model parts.
It is also computed as a Gaussian of the $OCM_{I}$.
While the initial proposal distribution is computed at every image location,
the conditional proposal distribution is only computed at regions of interest determined by the previously placed model parts.


\subsection{Experimental Results}
This section is based on our experimental evaluation reported in \cite{YangLateckiECCV2010}.
Specific implementation details of our shape similarity and model construction can be found in \cite{YangLateckiECCV2010}.

We have tested our algorithm on three widely used data sets: the extended Weizmann Horses \cite{WeizmanHorse,ShottonPAMI07}, the ETHZ shapes \cite{FerrariIJCV09} and the TU Darmstadt Database \cite{Leibe04}. During the testing for Weizmann Horses, only 12 automatically selected horse silhouettes with one hand decomposed horse are used to learn the shape model. All the other images are used for testing. The edge maps for this dataset are obtained by Canny edge detector.
We also test our method on the class of giraffe in ETHZ shape dataset \cite{FerrariIJCV09}.
The reason why we only select the category giraffes from ETHZ is that our model learning method can only transfer between objects with similar structure and giraffe is the only object in ETHZ having similar structure to horse.
Only one hand decomposed horse and 6 automatically selected giraffe silhouettes are used to learn the giraffe model.
Further, we work on the cow dataset the TU Darmstadt Database \cite{Leibe04},
since cows have similar structure with the above two classes.
It contains 111 images. Only one hand decomposed horse and 6 automatically selected cow silhouettes are used to learn the cow model. The edge maps for this dataset are obtained by Canny edge detector.

To adapt to large scale variance, we generate multiple models by resizing the original ones to 5 to 8 scales,
and choose as the final result from the best score in all the scales.
We not only report our results on the commonly used bounding box intersection,
but also the accuracy of our boundary localization.

We first evaluate the ability of the proposed approach to localize objects in cluttered images
using bounding-box intersection, which is widely used in traditional object detection task.
We adopt the strict standards of PASCAL Challenge criterion: a detection is counted as correct only
if the intersection-over-union ratio with the ground-truth bounding-box is greater than $50\%$.


Fig.~\ref{fig:PR_giraff} reports precision-recall (P/R) curve and detection rate vs false positive per image (DR/FPPI) curve
for the class Giraffes in ETHZ dataset \cite{FerrariIJCV09}.
In P/R, we compare to Lu et al. \cite{ChengEnICCV09}, Zhu et al. \cite{QihuiECCV08}, Ommer and Malik \cite{OmmerICCV2009} and Ferrari et al. \cite{FerrariIJCV09}, whose results are quoted from \cite{ChengEnICCV09}.
In DR/FPPI, as Ferrari et al. \cite{FerrariIJCV09,FerrariPAMI08}, Ommer and Malik \cite{OmmerICCV2009} and Lu et al. \cite{ChengEnICCV09} provide their results, we compare to them.
As Ravishankar et al. \cite{MittalECCV08} do not give their curves, we do not compare to them in Fig. \ref{fig:PR_giraff}.
According to the curves, we are better than Lu et al. \cite{ChengEnICCV09}, Ommer and Malik \cite{OmmerICCV2009}, Ferrari et al. \cite{FerrariIJCV09,FerrariPAMI08} and perform equally well as Zhu et al. \cite{QihuiECCV08}.
The performance of the proposed method illustrates its ability
to cope with substantial nonrigid deformations, which are present in the class Giraffes.
This is demonstrated by our example results in Fig.~\ref{fig: Total_Ex}(a).

\begin{figure*}[tbh]
\begin{center}
\begin{tabular}{c c}
\includegraphics*[width=5.5cm]{curves_0.5_PR_v2.eps}&\includegraphics*[width=6cm]{fppi_compare.eps}%{fppi_Giraffes_50_overlap_compare.eps}
\end{tabular}
\end{center}
\caption{Precision-recall curve and detection rate (DR) vs false positive per image (FPPI) curve for the class Giraffes in ETHZ dataset.}
\label{fig:PR_giraff}
\end{figure*}
%\vspace{-3mm}


\begin{figure*}[tbh]
\begin{center}
\begin{tabular}{c c}
\includegraphics*[width=8cm]{giraffes_image.eps}&
\includegraphics*[width=8cm]{horse_cow_image.eps}\\
(a) & (b)
\end{tabular}
\end{center}
\caption{Examples of detection results for Giraffes, horses and cows.}
\label{fig: Total_Ex}
\end{figure*}


Table \ref{Table:Cow_Horse} compares our detection rate to \cite{AlanYuillePAMI09,ShottonPAMI07}
on Weizman Hores and TU Darmstadt Cows.
The detection rate on horses is estimated from the DR/FPPI curve in \cite{ShottonPAMI07}.
The DR/FPPI curve for cows is not available in \cite{ShottonPAMI07}.
The method in \cite{AlanYuillePAMI09} is also matching based, while \cite{ShottonPAMI07} is a classification method.
Some examples of our horse and cow detection results are shown in Fig. \ref{fig: Total_Ex}(b).
The detection precision/recall area under curve (AUC) is a standard performance measure on the Weizmann Horses dataset.
The AUC for our approach is $79.84\%$, which is comparable to the result $80.32\%$ in Xiang et al.~\cite{XiangICCV09}.
We compare to them as they also use the explicit shape model and matching based method for object detection.
The AUC of classification based methods \cite{ShottonPAMI07,RandomForCVPR09} is $84.98\%$ and $96\%$, respectively.
We observe that classification based methods are bounding box classifiers
and utilize significantly more information than matching based methods as ours.
This explains why our detection rate and AUC is lower than \cite{ShottonPAMI07,RandomForCVPR09}.


\setlength{\tabcolsep}{2pt}
\begin{table}[h]
\begin{center}
\caption{Detection rate.}\label{Table:Cow_Horse}
\begin{tabular}{|c|c|c|c|}
\hline
& Our method & Zhu et al.~\cite{AlanYuillePAMI09}& Shotton et al.~\cite{ShottonPAMI07}\\
\hline
Horses&$93.97\%$&$86.0\%$&$95.20\%$\\
\hline
Cows&$90.38\%$&$88.6\%$&N/A\\
\hline
\end{tabular}
\end{center}
\end{table}

The proposed approach can not only succeed in extensive cluttered images,
but also handles the problem of large range of scales and intra-class variability.
This is demonstrated by several examples in Fig.~\ref{fig: Total_Ex}.
The images in the bottom right of Fig.~\ref{fig: Total_Ex}(a) with red rectangles are the ones we fail to detect.
The images of horses in Fig.~\ref{fig: Total_Ex}(b) with red rectangles are false positives
in the negative images provided by Shotton et. al.~\cite{ShottonPAMI07} to complement the Weizmann horse dataset.
They show that the false positives in the negative set are caused by really very cluttered edges
or by the structure of edges happening to match to the model very well.
Interestingly, the rightmost false positive of horses is due to a camel, whose shape is very similar to that of a horse.



The statistical inference framework presented in this paper offers one important advantage compared to texture based and classification methods like \cite{HOGcvpr05,RandomForCVPR09,MarPrizeICCV09}.
It can localize object boundaries, rather than just bounding-boxes.

In order to quantify how accurately the output shapes match to true boundaries, we use the coverage and precision measures defined in \cite{FerrariIJCV09}.
Coverage is the percentage of points from ground-truth boundaries closer than a threshold $t$ to the output shapes of the proposed approach.
Reversely, precision is the percentage of points from output shapes closer than $t$ to any point of ground-truth boundaries.
As in \cite{FerrariIJCV09} $t$ is set to $4\%$ of the diagonal of the ground-truth bounding box.
The measures are complementary.
Coverage captures how much of the object boundary has been recovered by the algorithm, whereas precision reports how much of the algorithm's output lies on the object boundaries.
These measurements seem to be more suitable for evaluating shape based approaches than
bounding-box evaluation, which cannot quantify the accuracy of contour localization.
It is possible to have bounding-box intersection larger than 0.5
without having correctly identified the ground-truth object boundaries.
Two examples of horse detection are shown in Fig. \ref{fig: Total_Ex}(b) with green rectangles.

The first two columns of Table \ref{Table:Pre_Cov_Gir} show coverage and precision averaged over all images of the class giraffes in ETHZ dataset in comparison to the results in \cite{FerrariIJCV09}.
We measure the coverage and precision for the correct detections at $0.4$ FPPI, following \cite{FerrariIJCV09}.
The coverage of the proposed approach is over $11\%$ better than \cite{FerrariIJCV09}, which shows that our approach can efficiently recover the true boundary of objects.
The precision is a little lower than \cite{FerrariIJCV09}.
More importantly, the detection rate at our $0.4$ FPPI is $86.75\%$.
However, even for $20\%$ bounding box intersection, the detection rate at $0.4$ FPPI in \cite{FerrariIJCV09} is only around $60\%$ ,
which is significantly lower than our result.



For horses and cows, the coverage and precision are obtained over all correct detections.
The third column of Table. \ref{Table:Pre_Cov_Gir} shows the coverage and precision of the proposed method on the Weizmann Horse dataset.
As the edges are significantly worse than the ones provided for the giraffes,
both measures are worse than the results on giraffes.
The coverage and precision results for cow are shown in the fourth column of Table. \ref{Table:Pre_Cov_Gir}.
Due to less intra-shape variance, the precision is $92.02\%$, which is much higher than giraffes and horses.
However, the coverage is only $73.86\%$.
The main reason for the difference between these two values is that our model has a gap,
since we removed the contour part representing the horse tail from the horse contour used for part decomposition.
Thus, even if the model and object match perfectly, the coverage score cannot be perfect (see examples in Fig. \ref{fig: Total_Ex}).

\setlength{\tabcolsep}{2pt}
\begin{table}[h]
\begin{center}
\caption{Accuracy of the boundary localization.}\label{Table:Pre_Cov_Gir}
\begin{tabular}{|c|c|c|c|c|}
\hline
&Ours& Results in \cite{FerrariIJCV09}& Ours&Ours\\
& on giraffes& on giraffes & on horses& on cows\\
\hline
Coverage&$79.4\%$&$68.5\%$&$77.5\%$&$73.86\%$\\
\hline
Precision&$74.6\%$&$77.3\%$&$61.7\%$&$92.02\%$\\
\hline
\end{tabular}
\end{center}
\end{table}



\section{Related Work}\label{sec:RelatedWork}

Particle filters (PF) are also known as sequential Monte Carlo methods (SMC) for model estimation based on simulation.
There is large number of articles published on PF and we refer to two excellent books \cite{Doucet2001,Liu2001} for an overview.
To our best knowledge, the proposed PF framework with unordered observations is novel and has not been considered before by other authors.
PF can be viewed as a powerful inference framework that is utilized in many applications.
One of the leading examples is the progress in robot localization and mapping based on PF \cite{ThrunBook05}.
Particle filter (PF) has been used for object detection previously \cite{ParticleIJCV2001,ParticleICCV1999}.
They mainly utilize PF in the classical tracking/filtering scenario with
a pre-defined order of states and observations.
In contrast, our method explores different orders of states and observations on the fly and let particles follow the most informative orders, 
which is theoretically quite different from the traditional PF framework.

The proposed PF framework is also related to Gibbs sampling.
The algorithm for Gibbs sampling was first introduced in the seminal paper by Geman and Geman \cite{GemanGemanPAMI1984}.
Since then there have been many sampling algorithms like Gibbs sampler,
e.g., Hot Coupling \cite{hot_coupling_nips2005}, Tree sampling, Swendsen-Wang sampling etc.
But most of them assume restrictive conditional independencies.
Recently Hamze et. al. proposed a very generic importance sampling method called Large Flip Importance Sampling (LFIS) to sample from the posterior \cite{hamze_freitas_uai2007}. 
The main motivation for their approach comes from N-Fold Way \cite{bortz_etal:jcp1975} and Tabu search \cite{glover:jc1989}, 
where they use heuristics to improve the sampling of the exponential state space using memory and heuristics to design good moves in the state space.
Since the moves are no-longer MCMC in the traditional sense they introduce importance weights 
to the distinct states visited by $N$ copies of the sampler.
In this paper we combine the strengths of both PF and Gibbs sampler based approaches.
In particular, our approach can be viewed as improved Gibbs sampler that employs the PF weighting scheme as well as the PF resampling scheme.


A preliminarily version of the proposed PF inference framework has also been applied to object detection and recognition in \cite{ChengEnICCV09},
where the shape similarity was based on direct matching of model contour parts to image contour fragments.
The presented PF inference framework has been first applied to object detection and recognition in \cite{YangLateckiECCV2010}.
This article focused on theoretical presentation of the proposed PF inference framework
and provides derivation of all key formulas as well as proves their theoretical properties.


%\section{Rao-Blackwellized Particle Filter <we probably do not need this section>}
%The following derivation is inspired by Rao-Blackwellized particle filter often used in SLAM \cite{Grisetti07}.
%Our goal is to derive the posterior $ p(x_{1:m}, \mu |Z)$ not only over the set of stats $x_{1:m}$ but also over possible shape models $\mu$.
%To make this estimation possible, we use the following factorization
%\begin{align}\label{eq:RB1}
%& p(x_{1:t}, \mu |Z) = p(\mu | x_{1:t}, Z)  \ p(x_{1:t} | Z)
%\end{align}
%for $t=1, \dots, m$.
%This factorization allows us to first estimate only the set of states
%and then compute the shape model.
%The key assumption here is that the shape model is a deterministic function of the set of states.
%Hence, the posterior over models $p(\mu | x_{1:t}, Z)$ can be computed analytically
%since $x_{1:t}, Z$ are known.
%This technique is often referred to as Rao-Blackwellization.
%
%By substituting Eq.~\eqref{eq:RB1} to Eq.~\eqref{eq:PFn1.1}, we obtain
%\begin{align}\label{eq:RB2}
%& p(x_{1:m}, \mu |Z) \approx p(\mu | x_{1:m}, Z)  \sum_{i=1}^N w(x_{1:m}^{(i)}) \delta(x_{1:m} - x_{1:m}^{(i)}) \nonumber \\
%&=\sum_{i=1}^N p(\mu | x_{1:m}, Z) w(x_{1:m}^{(i)}) \delta(x_{1:m} - x_{1:m}^{(i)})
%\end{align}
%By denoting $ u(x_{1:m}^{(i)}) = p(\mu | x_{1:m}^{(i)}, Z) w(x_{1:m}^{(i)})$, we obtain
%\begin{align}\label{eq:RB2}
%& p(x_{1:m}, \mu |Z) \approx \sum_{i=1}^N u(x_{1:m}^{(i)}) \delta(x_{1:m} - x_{1:m}^{(i)}),.
%\end{align}
%Consequently, if we modify the weight \eqref{eqn:partweight} of each particle to
%\begin{align}\label{eqn:RBpartweight}
%&u(x_{<1:t-1>}^{(i)}) = p(\mu | x_{<1:t-1>}^{(i)}, Z) p(Z|x_{<1:t-1>}^{(i)}),
%\end{align}
%we obtain a PF algorithm that estimates the posterior $p(x_{1:t}, \mu |Z)$.






\section*{Acknowledgments}
The work has been supported in part by the NSF Grants IIS-0812118, BCS-0924164, the AFOSR Grant FA9550-09-1-0207,
and the DOE Award 71498-001-09.

\bibliographystyle{IEEEtran}
\bibliography{Bib/RobotMapping/EMrobot,Bib/RobotMapping/FFS,Bib/Shape/shape,Bib/Shape/contours,Bib/BMI/bmi,Bib/BMI/mrfs,Bib/Shape/seg,Bib/ModelLearning2010}



\end{document} 