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

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

\usepackage{color}

\graphicspath{{figures/}}
% Include other packages here, before hyperref.

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


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

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

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

%%%%%%%%% TITLE
\title{Particle Filter with Unordered Observations for Solving Image Jigsaw Puzzles}

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

\maketitle
% \thispagestyle{empty}

%%%%%%%%% ABSTRACT
\begin{abstract}
  We deal with an image jigsaw puzzle problem, which is defined as reconstructing an image from a set of square and non-overlapping image patches.
  It is known that this is a NP-complete problem, and it is also challenging for humans,
  since in the considered framework the original image is not given.
  Recently a graphical model has been proposed to solve this problem, where each patch location is a node and each image patch is a label.
  The target label probability function is then maximized using loopy belief propagation.
  In this paper we use the same target function and exactly the same pairwise potentials.
  Our main contribution is a novel inference approach in the sampling framework of Particle Filter (PF).
  Usually in the 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 the jigsaw puzzle problem all observations (puzzle pieces) are given at once without any particular order.
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 the PF inference.
  Our experimental results demonstrate that the proposed inference framework significantly outperforms loopy belief propagation
  in solving the image jigsaw puzzle problem.
  In particular, the extended PF inference triples the accuracy of the label assignment of the loopy belief propagation.
\end{abstract}

%%%%%%%%% BODY TEXT
\section{Introduction and Problem Formulation}
As shown in \cite{PuzzleDemaine1988} the jigsaw puzzle problem is
NP-complete if the pairwise affinity among jigsaw pieces is unreliable.
Following \cite{PuzzleMIT2010},
we focus on reconstructing the original image from square and non-overlapping patches. 
This type of puzzles does not contain the shape information of individual pieces,
 which is quite important to determine the pairwise affinities among them. 
 This makes the problem more challenging, since it is more difficult to evaluate pairwise affinities among puzzles.
Thus, only the image content information of the puzzle pieces is available.
This is different from most of the previous approaches \cite{PuzzleKimia2001,PuzzleGoldberg2002,PuzzleRadack1982,PuzzleWolfson1988},
where the shape of the puzzle pieces is utilized.



\begin{figure*}[t]
\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} {\color{red} change this example to 6 squares.}
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 as in (c), following their numbering may not be helpful in solving the puzzle.}
\end{figure*}


Now we briefly review the PF inference.
We begin with a classical tracking example.
A robot is moving around and taking images at discrete time intervals.
The images form a sequence of observations $Z=( z_1, \dots, z_m )$, where $z_t$ is an image taken at time $t$.
With each observation $z_t$ there is associated a hidden state $x_t$.
In our example the value of $x_t$ is the robot pose (its 2D position plus orientation).
The goal of PF inferences, is to derive the most likely sequence of the hidden states,
i.e., to find a state vector $\textbf{x}=(x_1, \dots, x_m )$ that maximizes the posterior $P(x_1, \dots, x_m | Z)$.
We observe that here the observations are ordered following their time stamps.
In PF inference, this order is utilized to sequentially infer the values of states $x_t$ for $t=1, \dots, m$.
Now imagine that the robot's clock broke and the time stamps are random.
Thus, we are given a set of observations $Z=\{ z_1, \dots, z_m \}$, they are indexed but their index is irrelevant.
Of course, we can still associate state $x_t$ with observation $z_t$, but the set of observations is not ordered,
and consequently, there state vectors $x_t$ are not ordered.
Thus, we deal with unordered observations.
This is exactly the scenario of the image jigsaw puzzle problem, see Fig.~\ref{fig:puzzle}.
We are given $m$ square puzzle pieces described by a set of $m$ observations $Z=\{ z_1, \dots, z_m \}$.
Each observation $z_t$ describes part of the original image on piece $t$ and is given by a vector of features,
which are the color values of the pixels on piece $t$ in our experimental results.
The puzzle pieces are numbered with index $t$, but their numbering is random.
Our goal is to determine the state vector $\textbf{x}=(x_1, \dots, x_m )$ that
maximizes the posterior probability $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$.
A pose of a puzzle piece is simply its 2D location in a grid of empty squares,
i.e., the value of each $x_t$ is the index or label of an empty square in the square gird.
Since there is no particular order of the puzzle pieces,
the posterior distribution is a function of how well adjacent pieces fit together once the pieces a placed on the grid,
i.e., assigned the state values.
We 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 image information of the puzzle pieces is not very descriptive.

Our main contribution is a new PF inference framework that works in this scenario.
In the proposed framework we extend PF to handle the situations where we have unordered set of observations that are given simultaneously.
One of our key ideas 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.
This allows us to simultaneously consider different orders of observation in the process of 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.


In our experimental results, we compare the solutions of \eqref{Eq: MitTarget} obtained by the proposed inference framework
to the solutions of the loopy believe propagation used in \cite{PuzzleMIT2010} under identical settings
on the dataset provided in \cite{PuzzleMIT2010}.
In particular, we use exactly the same dissimilarity-based compatibility $D$ as defined above.
The proposed PF inference significantly outperforms the loopy believe propagation in all evaluation measures.
The main of these measures is the accuracy of the label assignment,
where the difference is most significant.
The accuracy of the loopy believe propagation is $24\%$ while the accuracy of the proposed PF inference is $69\%$.





The classical PF framework has been developed for sequential state estimation like tracking \cite{PFTrackingECCV2004,PFTrackingCVPR2005}
or robot localization \cite{Thrun02d,Fox00Par}.
There, the observations arrive sequentially and are indexed by their time stamps, as our tracking example illustrates.
It is possible to apply the classical PF framework as stochastic optimization to solve this problem by utilizing
a fix order of observations.
However, by doing so, we would have selected an arbitrary order, and the puzzle construction may fail because of the selected order.
If we would have selected a different order, the labeling task could have been successful.
Moreover, the observations are given simultaneously at the same time.
Consequently, there is no reason to favor any particular order without utilizing this fact.




The rest of the paper is organized as follows. 
Section \ref{sec: RelatedWork} describes several closely related approaches. 
The details of the proposed framework are explained in Section \ref{sec: framework}. 
Section \ref{sec: Ex} shows and evaluates the experimental results not only the dataset from \cite{PuzzleMIT2010}, 
but also an extended dataset.

%
%
%\begin{figure}[t]
%\begin{center}
%\begin{tabular}{c c}
%\includegraphics*[width=4cm]{Greedy_Results_20_2.eps}& \includegraphics*[width=4cm]{Iter_107.eps}\\
%\end{tabular}
%\end{center}
%\caption{{\color{red} show first the puzzle pieces, then random, then the MIT  and our results.}
%The left image is the result of greedy algorithm; The right image is the result of our algorithm. The top left patch (red) is the starting patch.}
%\label{fig: GreedyCom}
%\end{figure}
%




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


In this section we outline a {\bf standard PF algorithm}.
It optimizes the posterior distribution $p(x_1, \dots, x_m | Z)$ following the order of observations $Z=( z_1, \dots, z_m )$.
For every time step $t=1, \dots, m$ and for every particle $i =
1, \dots, N$ it executes the following three steps:

\noindent 1) \textbf{Importance sampling / proposal:}
Sample followers of particle $(i)$ according to 
\begin{align}\label{eq:sampling2}
x_{t}^{(i)} \sim p_t(x|x_{1:t-1}^{(i)}),
\end{align}
and set
$x_{1:t}^{(i)} = (x_{1:t-1}^{(i)}, x_{t}^{(i)})$, where
$x_{1:t}$ is a short hand for vector of states $(x_1, \dots, x_t)$.

\noindent 2) \textbf{Importance weighting / evaluation:}
An importance weight is assigned
to each particle $x_{1:t}^{(i)}$ according to
\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}

\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:sampling2} and \eqref{eq:PFn4} 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$.
Intuitively, $f$ allows us to predict the value of state $x_{t}$ 
if we know the value of the previous state $x_{t-1}$,
while $h$ allows us to evaluate the accuracy of the prediction of the value of state $x_{t}$ according to the observations obtained at time $t$.




\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 a local image very similar to many other puzzle pieces
and the second puzzle pieces has a very distinctive local image that matches only a few other pieces,
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:sampling2}
and particle evaluation in Eq.~\ref{eq:PFn4} utilize the sequential order of the states $x_{1:m}=(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.

The key idea of the proposed approach is not to utilize the fix order of the states,
but instead explore all possible orders of the states $(x_{i_1}, \dots, x_{i_m})$
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 the jigsaw puzzle problem, the neighborhood is composed of adjacent squares in the square grid.
%Consequently, \eqref{eq:PFa3} means that we can estimate the relative pose of the remaining puzzle pieces
%around the already positioned pieces, i.e., position them adjacent to one of the 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 state 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 state
dimension $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=6$ pieces in our jigsaw puzzle and after iteration 4,
we determined the poses of puzzle pieces $<1:4>=(3, 1, 5, 6)$.
We then consider the two remaining pieces as followers of $<1:4>$ at iteration 5.
The poses of pieces 2 and 4 are determined with reference to the poses of pieces $3, 1, 5, 6$,
i.e., by sampling
$x_{2}^{(i)} \sim p_2(x|x_{<1:4>}^{(i)})$ and $x_{4}^{(i)} \sim p_4(x|x_{<1:4>}^{(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 state),
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}
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}
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} \label{eq:unsamp}
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 1.} \
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 1.} \
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{Implementation Details} \label{sec:impdet}
In order to implement the proposed PF inference framework, we need to define the proposal distribution \eqref{eq:unsamp} and 
the importance weighting \eqref{eqn:partweightv2}.
Given is an rectangular grid of $m$ empty squares $G=\{g_1, \dots, g_m\}$, 
where each $g_i$ is the index of the unique location on the grid.
$G$ is the set of locations assigned to the puzzle pieces $P= \{1, \dots, m \}$.
A partial label assignment is given by a state vector $x_{<1:t>} =(x_{<1>}, \dots, x_{<t>})$,
where each $<i> \in P$ is an index of a puzzle piece and  $x_{<i>} \in G$ is the location of piece $<i>$;
meaning that piece $<i>$ is placed at location $x_{<i>}$.

On the grid $G$ we have a standard neighborhood structure denoted $GN$,
where $GN(g)$ is the four  squares adjacent to square $g$, i.e., the left/right/top/bottom squares, 
if $g$ is no a boundary square. The neighborhood is smaller for boundary squares in the grid.
$GN$ naturally extends to a set of grid locations $H \subset G$, i.e., $GN(H) = \cup \{GN(g) | g \in H \}$.

The proposal distribution over $G$ is given by the probability 
of putting a puzzle piece $s \in P \setminus <1:t-1>$ at any location $x$ on grid $G$ defined as 
\begin{equation}
\label{eq:propdef}
p_s(x|x_{<1:t-1>}) \propto \left\{ \begin{array}{rl}
 1 &\mbox{ if $x \in GN(x_{<1:t-1>})$} \\
  0 &\mbox{ otherwise}
       \end{array} \right.
\end{equation}
Thus, location $x$ has a positive probability if and only if it is adjacent to one of the occupied locations $x_{<1:t-1>}$.


Now we define the cognitional distribution required for the weight definition in \eqref{eqn:partweightv2}:
\begin{align}\label{eq:pz1}
&p(Z|x_{<1:t>}) \propto p(x_{<1:t>} | Z) p(x_{<1:t>}),
\end{align}
where $p(x_{<1:t>})$ is defined below and we use the following factorization 
\begin{align}\label{eq:pz2}
& p(x_{<1:t>} | Z) = p(x_{<1>} | Z) \prod_{l=2}^t p(x_{<l>} | x_{<1:l-1>}, Z)
\end{align}
In order to use the same evaluation function as in \cite{PuzzleMIT2010} in our framework,
we assume that $p(x_{<l>} | x_{<1:l-1>}, Z) = p(x_{<l>} | x_{<l-1>}, Z)$, which is the standard Markov assumption.
Consequently, we obtain
\begin{align}\label{eq:pz3}
& p(x_{<1:t>} | Z) = p(x_{<1>} | Z) \prod_{l=2}^t p(x_{<l>} | x_{<l-1>}, Z)
\end{align} 

We use the dissimilarity-based compatibility $D$ from \cite{PuzzleMIT2010} to define $p(x_{<k>} | x_{<k-1>}, Z)$.
Let $<k>=j$ and let $<k-1>=i$ for $i \neq j$.
There are four types of $D$ determined by the spatial adjacency of $x_j$ and $x_i$ in the square grid $G$,
which can be $x_j$ is to the left/right/top/bottom of $x_i$.
$D$ measures dissimilarity between puzzle piece observations $z_j$, $z_i$ by summing the squared $LAB$ color difference along the boundary,
e.g., the left-right (LR) dissimilarity between $x_j$ and $x_i$ is based on color differences of adjacent pixels
\begin{equation}
D_{LR}(x_j, x_i) = \sum_{k=1}^K\sum_{c = 1}^3 (z_j(k,u,c) - z_i(k,v,c))^2
\end{equation}
if $x_j$ is to the left of $x_i$ and zero otherwise.
where each observation $z_t$ is represented by a local image matrix of puzzle piece $t$ of size $K \times K \times 3$,
$u$ indexes the last column of $z_j$ and $v$ indexes the first column of $z_i$.
$D(x_j,x_i)$ is not a distance, since it is not symmetric, i.e., $D(x_j,x_i) \neq D(x_i,x_j)$.
Finally, the probability $p_{i,j}(x_j|x_i,Z)$ is obtained by
\begin{equation} \label{eq:binpot}
p(x_j|x_i, Z) = p(x_j|x_i, z_j, z_i) \propto \exp(- \frac{D(x_j,x_i)}{2 \sigma^2})
\end{equation}
if $x_j$ is adjacent to $x_i$ and zero otherwise.
Here $\sigma$ is adaptively set as the difference between the smallest and the second smallest $D(x_j,x_i)$ among all $x_j$,
see \cite{PuzzleMIT2010} for more details.
The equality $p(x_j|x_i, Z) = p(x_j|x_i, z_j, z_i)$ is based on the assumption that 
conditioned on $x_i$ the probability of placing $x_j$ next to $x_i$ only depends on the local images $z_j, z_i$.

Finally $p(x_{<1:t>})$ represents constraints on the spatial arrangement of the locations of current puzzle pieces on the square grid. 



Hence, the posterior distribution is a function of how well adjacent pieces fit together as defined in \eqref{eq:binpot}.




\section{The Image Jigsaw Puzzle Problem as Labeling in a Graphical Model} \label{sec:MIT}
In order to compare the proposed PF inference to the loopy believe propagation for solving the image jigsaw puzzle,
we use exactly the same functions to define the probabilities as in \cite{PuzzleMIT2010}.
First we briefly present the problem formulation introduced in \cite{PuzzleMIT2010},
where the image jigsaw puzzle problem is expressed as a labeling problem in a graphical model.
Each square grid location is a node in the graph and each puzzle piece (patch) is a label.
This formulation flips the previously considered label assignment problem.
We indexed each state variable $x_t$ with a puzzle piece number, and assigned to $x_t$ the location in the square gird.
In \cite{PuzzleMIT2010}, the $x_t$ state variable is indexed with the square grid location,
and it is assigned a number of a puzzle piece as a label.
Obviously, the two formulations are equivalent, and from now on we follow the formulation in \cite{PuzzleMIT2010}.

In order to solve a given image jigsaw puzzle,
we seek a state vector $\textbf{x} = (x_1, \dots, x_m)$ that maximizes the probability
\begin{equation}
P(\textbf{x} | Z) = \frac{1}{\zeta}\prod_{i = 1} ^{m} \prod_{j\in \mathcal{N}(i)} p_{i,j}(x_j|x_i, Z) E(\textbf{x})
\label{Eq: MitTarget}
\end{equation}
where the value of $x_i$ represents the puzzle piece assigned to location $i$,
i.e., a puzzle piece placed at the grid location $i$,
and $p_{i,j}(x_j|x_i)$ is the probability of placing a patch $x_j$
on the square grid in the neighborhood $\mathcal{N}(i)$ of another patch $x_i$.
Hence the neighborhood $\mathcal{N}$ from Section \ref{sec:new} is instantiated as the square grid neighbors
of the already located puzzle pieces.
Since patches should not be used more than once, $E(\textbf{x})$ is used to softly enforce this constraint.
Finally, $\zeta$ is a normalizing constant.

We observe that \eqref{Eq: MitTarget} is a special case of equation (1) in \cite{PuzzleMIT2010}.
It is obtained by dropping
the unary potentials $p(x_i)$ and the probability term that evaluates the probability of label assignment $x_i$
with respect to some knowledge of the whole image layout.
In \cite{} {\color{red} ref [4] in the MIT paper}, the knowledge of the image layout is given by a low resolution version of the original image.
\cite{PuzzleMIT2010} weakens this assumption to a statistics of the possible image layout.
In contrast, we do not use any knowledge of the image layout,
i.e., we do not have any any prior knowledge about the layout of the target image.
Consequently, we consider a harder problem, since we only use the pairwise relation between the image patches.
It is given by pair-wise compatibilities of located puzzle pieces.

The authors of \cite{PuzzleMIT2010} do not explicitly condition the probability in \eqref{Eq: MitTarget} on the observations $Z$,
but this condition is implicitly present in their assumptions.


\section{Experimental Results} \label{sec: Ex}
% the experimental results on MIT dataset and the extended one
We conduct experiments not only on the MIT dataset with 20 images, but also the extended 30 images we downloaded from internet. The performance evaluation methods are introduced in Section \ref{sec: PerformEval}. The experimental results on MIT dataset are shown in Section \ref{sec: ExMIT} and we show the results on extended dataset in Section \ref{sec: ExExtended}. To fairly compare to \cite{PuzzleMIT2010}, we also assign one patch to its groundtruth location in their algorithm. Moreover, since the performance of \cite{PuzzleMIT2010} on 108 pieces of puzzles is already very low, we do not need compare results with them on more number of pieces.

\subsection{Performance Evaluation Methods}\label{sec: PerformEval}
There are three types of evaluation methods \cite{PuzzleMIT2010}, each focuses on different aspects. $x_i$ represents the patch assigned by the proposed algorithm at location $i$. $GT(i)$ is the ground truth patch for location $i$.

The most straightforward and strict one is \textbf{Direct Comparison} (see Eq. \ref{Eq: DC}), where the inferred labels are compared to the ground truth labels directly. The accuracy measures the extent of correct localization of patches.
\begin{equation}
DC = \frac{1}{N} \sum_{i = 1}^N (x_i = GT(i))
\label{Eq: DC}
\end{equation}

However, since the differences between small patches are not descriptive enough, similar patches may be assigned to wrong locations. To measure this kind of situation, the inferred patch labels are mapped to clusters they belong to, and are compared to the ground truth cluster label, which is called \textbf{Cluster Comparison}(see Eq. \ref{Eq: CC}). The measurement tolerates the assignment error as long as the assigned patch belongs to the same cluster as the ground truth path.

\begin{equation}
CC = \frac{1}{N} \sum_{i = 1}^N (Cluster(x_i) = Cluster(GT(i)))
\label{Eq: CC}
\end{equation}

Moreover, due to lack of prior knowledge of target image, the reconstructed image may shift compared to the ground truth. For each assigned patch label, we compute the fraction of four neighbor nodes that the algorithm assigned the correct patch. The reconstruction accuracy is the average fraction of correct neighbor labels. This measurement do not penalize the accuracy as long as the adjacent patches in original image are adjacent in the reconstructed image.

\subsection{MIT Dataset}\label{sec: ExMIT}
MIT Dataset contains 20 images and we evaluate with two different number of patches 108. The results are shown at Table \ref{Table: ExMIT108}. The process of our algorithm to reconstruct the images is illustrated in Fig. \ref{fig:PF_Iter}. It is obvious that the order of Random variables is different for different images, which is automatically determined by the algorithm. The experimental results have demonstrated the excellent performance of our algorithm.

\begin{table}
\begin{center}
\begin{tabular}{|c|c|c|}
\hline
 & \cite{PuzzleMIT2010} & Our algorithm \\
 \hline
 Direct Comparison &0.2366& 0.6921\\
 \hline
 Neighbor Comparison &0.6628 & 0.8620\\
 \hline
 Cluster Comparison & 0.4657& 0.7810\\
 \hline
\end{tabular}
\end{center}
\caption{Experimental Results on MIT Dataset with 108 pieces of jigsaws}
\label{Table: ExMIT108}
\end{table}


\begin{figure*}[t]
\begin{center}
\begin{tabular}{|c|c|c|c|c|c|}
\hline
$Iter=1$ & $Iter=10$ & $Iter=20$ & $Iter=30$ & $Iter=40$ & $Iter=50$ \\
\hline
\includegraphics*[width=2.5cm]{Iter_1.eps}& \includegraphics*[width=2.5cm]{Iter_10.eps} & \includegraphics*[width=2.5cm]{Iter_20.eps}& \includegraphics*[width=2.5cm]{Iter_30.eps}& \includegraphics*[width=2.5cm]{Iter_40.eps}& \includegraphics*[width=2.5cm]{Iter_50.eps}\\
\hline
$Iter=60$ & $Iter=70$ & $Iter=80$ & $Iter=90$ & $Iter=100$ & $Iter=107$\\
\hline
\includegraphics*[width=2.5cm]{Iter_60.eps}& \includegraphics*[width=2.5cm]{Iter_70.eps}& \includegraphics*[width=2.5cm]{Iter_80.eps}& \includegraphics*[width=2.5cm]{Iter_90.eps}& \includegraphics*[width=2.5cm]{Iter_100.eps}& \includegraphics*[width=2.5cm]{Iter_107.eps}\\
\hline
\end{tabular}
\end{center}

\begin{center}
\begin{tabular}{|c|c|c|c|c|c|}
\hline
$Iter=1$ & $Iter=10$ & $Iter=20$ & $Iter=30$ & $Iter=40$ & $Iter=50$ \\
\hline
\includegraphics*[width=2.5cm]{Iter_2_1.eps}& \includegraphics*[width=2.5cm]{Iter_2_10.eps} & \includegraphics*[width=2.5cm]{Iter_2_20.eps}& \includegraphics*[width=2.5cm]{Iter_2_30.eps}& \includegraphics*[width=2.5cm]{Iter_2_40.eps}& \includegraphics*[width=2.5cm]{Iter_2_50.eps}\\
\hline
$Iter=60$ & $Iter=70$ & $Iter=80$ & $Iter=90$ & $Iter=100$ & $Iter=107$\\
\hline
\includegraphics*[width=2.5cm]{Iter_2_60.eps}& \includegraphics*[width=2.5cm]{Iter_2_70.eps}& \includegraphics*[width=2.5cm]{Iter_2_80.eps}& \includegraphics*[width=2.5cm]{Iter_2_90.eps}& \includegraphics*[width=2.5cm]{Iter_2_100.eps}& \includegraphics*[width=2.5cm]{Iter_2_107.eps}\\
\hline
\end{tabular}
\end{center}

\caption{The reconstructed image of the best particle at different iterations. There are 108 pieces of patches for each image and the top left patch (red patch) is given as start point. The two examples demonstrate the process of reconstructing image by particle filter.}
\label{fig:PF_Iter}
\end{figure*}

\subsection{Extended Dataset} \label{sec: ExExtended}
We randomly select 30 images from Google image with classes from: football, sky, airplane and so on. This will demonstrate that the algorithm is able to perform well on general cases. The results are shown at Table \ref{Table: Ex2MIT108}. Some Candidate experimental results and comparison to \cite{PuzzleMIT2010} are shown in Fig. \ref{fig:Comp}. The experimental results have demonstrated the excellent performance of our algorithm.

\begin{table}
\begin{center}
\begin{tabular}{|c|c|c|}
\hline
 & \cite{PuzzleMIT2010} & Our algorithm \\
 \hline
 Direct Comparison &0.4034& 0.7234\\
 \hline
 Neighbor Comparison &0.7129 & 0.8865\\
 \hline
 Cluster Comparison & 0.5841& 0.8210\\
 \hline
\end{tabular}
\end{center}
\caption{Experimental Results on Extended Dataset with 108 pieces of jigsaws}
\label{Table: Ex2MIT108}
\end{table}


\begin{figure*}[t]
\begin{center}
\begin{tabular}{c c c c c}

\includegraphics*[width=3cm]{Original_5.eps}& \includegraphics*[width=3cm]{Original_9.eps}& \includegraphics*[width=3cm]{Original_13.eps}& \includegraphics*[width=3cm]{Original_14.eps}& \includegraphics*[width=3cm]{Original_18.eps}\\

\includegraphics*[width=3cm]{MIT_Results_4_1.eps}& \includegraphics*[width=3cm]{MIT_Results_8_1.eps} & \includegraphics*[width=3cm]{MIT_Results_12_1.eps}& \includegraphics*[width=3cm]{MIT_Results_13_1.eps}& \includegraphics*[width=3cm]{MIT_Results_17_1.eps}\\

\includegraphics*[width=3cm]{Our_Results_5_1.eps}& \includegraphics*[width=3cm]{Our_Results_9_1.eps}& \includegraphics*[width=3cm]{Our_Results_13_1.eps}& \includegraphics*[width=3cm]{Our_Results_14_1.eps}& \includegraphics*[width=3cm]{Our_Results_18_1.eps}\\

\end{tabular}
\end{center}
\caption{The comparison between Our algorithm and \cite{PuzzleMIT2010} on the extended Dataset. The images in the first row are the original image. The results of Cho et. al \cite{PuzzleMIT2010} are shown in the second row. The images in the third row are the results of our algorithm. It is clear that the proposed method obtains better results.}
\label{fig:Comp}
\end{figure*}



\section{Related Work} \label{sec: RelatedWork}
The first work on Jigsaw Puzzle Problem is \cite{PuzzleFreeman1964}. Since shape is an important clue for accurate pairwise relation, many methods \cite{PuzzleKimia2001,PuzzleGoldberg2002,PuzzleRadack1982,PuzzleWolfson1988} focus on matching distinct shapes among jigsaw pieces to solve the problem. The pairwise relations among jigsaw pieces are measured by the fitness of shapes. It is straightforward to consider both the shape and image content \cite{PuzzleICIP2006,PuzzlePRL2008,PuzzlePRL2003}. Most methods solve the problem with a greedy algorithm and report results on just one or few images. Recently, Cho et. al \cite{PuzzleMIT2010} propose to use loopy belief propagation to solve the problem and evaluate on a test dataset.

{\color{red} describe PF approaches}



\section{Conclusion} \label{sec: Conclusion}

%-------------------------------------------------------------------------

{\small
\bibliographystyle{ieee}
\bibliography{JigsawPuzzle}
}


\end{document}
