\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.21\linewidth]{Lighthouse.eps} &
(b) & \includegraphics[width=0.3\linewidth]{Lighthouse_6piece.eps} &
(c) & \includegraphics[width=0.21\linewidth]{Lighthouse_2piece.eps}
\end{tabular}
%\vspace{-5pt}
\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).
The empty squares in (c) form possible locations for the puzzle pieces in (b).}
\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 $x_{1:m}=(x_1, \dots, x_m )$ that maximizes the posterior $p(x_{1: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, e.g., the number in Fig.~\ref{fig:puzzle}(b).
The value of the state $x_t$ of puzzle piece $t$ is a location of an empty square in the square grid,
e.g., the value of $x_t$ is the index of an empty square in the square gird shown in Fig.~\ref{fig:puzzle}(c).
Our goal is to determine the state vector $x_{1:m}$ that
maximizes the posterior probability $p(x_{1:m} | Z)$.
Since the original image is not provided,
this probability is determined based on pairwise appearance consistency of the local puzzle images,
i.e., the posterior distribution is a function of how well adjacent pieces fit together once they are placed on the grid.
In other words, a vector of grid locations $x_{1:m}$ maximizes $p(x_{1:m} | Z)$ if
the puzzle pieces placed at these locations form the most consistent image.
We observe that the posterior distribution $p(x_{1: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 summarize basic facts about Particle Filters (PFs).
They will be utilized in the following section when we introduce the proposed framework.
Given is a sequence of observations $Z=\{ z_1, \dots, z_m \}$,
i.e., the observations are ordered.
Our goal is to maximize the posterior distribution $p(x_{1:m} \ | \ Z )$,
i.e., to find the values $\widehat{x}_{t}$ of states $x_t$ such that
\begin{align}
& \widehat{x}_{1:m} = \underset{x_{1:m}}{\operatorname{argmax}} \ p(x_{1:m} \ | \ Z),
\label{eqG2}
\end{align}
where $x_{1:m}= (x_1, \dots, x_m) \in \mathcal{X}^m$ is a state space vector.
We observe that each state $x_t$ corresponds to observation $z_t$ for $t=1, \dots, m$.

This goal can be achieved 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).
The state vector $x_{1:m}^{(i)}$ with the largest weight $w^{(i)}$ is then taken as the solution of \eqref{eqG2}.

Since it is still computationally intractable 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 state vector $x_{1:m}= (x_1, \dots, x_m)$:
\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 vector the samples are generated.
We use this notation to stress that $q_t$ is a probability density function (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_1(x_1 | Z) \prod_{t=2}^m q_t(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 \cite{ThrunBook05}:
\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_t(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}).
Consequently, we have arrived at the SIS theorem: \\

\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$.
Hence 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}


Weight equation \eqref{eq:PFn3} can be simplified by making a common assumption that
$q_t(x_{t}^{(i)}|x_{1:t-1}^{(i)},z_{1:t}) = p(x_{t}^{(i)}|x_{1:t-1}^{(i)})$,
i.e., we take as the proposal distribution the conditional pdf of the state at time $t$ conditioned on the current state vector
$x_{1:t-1}^{(i)}$.
This assumption simplifies the recursive weight update formula 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}
and implies that the samples are generated from
\begin{align}\label{eq:sampling2}
x_{t}^{(i)} \sim p_t(x|x_{1:t-1}^{(i)}).
\end{align}
Analog to \eqref{eq:sampling} we added the index $t$ to $p$ in \eqref{eq:sampling2}
to indicate the dimension of the state space from which the samples are generated.

Now we summarize the derived {\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 $N$ 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$, and renormalize their weights to sum to one.
This procedure is a variant of 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$.
%Intuitively, $f$ allows us to predict the value of state $x_{t}$
%if we know the value of the previous states $x_{1: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 the most likely 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.






In order to formally define the proposed sampling rule, w need to explicitly represent different orders of states
with a permutation
$\sigma: \{1, \dots, m \} \to \{1, \dots, m \}$ and use the shorthand notation
$\sigma(1:t)$ to denote $(\sigma(1), \sigma(2), \dots, \sigma(t))$ for $t \leq m$.
For each particle $(i)$ we may have a different permeation $\sigma^{(i)}$,
but we will drop the index $(i)$ from $\sigma^{(i)}$ in expressions like $x^{(i)}_{\sigma(1:t)}$,
since the state vector already carries the particle index $(i)$.
With reference to our toy example, Fig.~\ref{fig:puzzle}(c) shows the configuration at time $t=2$,
where puzzle pieces numbered 3 and 1 in Fig.~\ref{fig:puzzle}(b) are placed at locations 1 and 2, correspondingly.
Hence $\sigma^{(i)}(1:2)=(3,1)$ and $x^{(i)}_{\sigma(1:2)} = (1, 2)$.
We stress that now the sequence of states $x_{\sigma(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 \}$.


We are ready now to precise formulate the proposed importance sampling.
At each iteration $t \leq m$, for each particle $(i)$ and for each $s \in \overline{ \sigma^{(i)}(1:t-1) }$, we sample
\begin{align}\label{eq:sampling2.5}
x_{s}^{(i)} \sim p_s(x|x_{\sigma(1:t-1)}^{(i)}),
\end{align}
where $\overline{\sigma^{(i)}(1:t-1)} = (1, \dots, m ) \setminus \sigma^{(i)}(1:t-1)$,
i.e., the indices in $1:m$ that are not present in $\sigma^{(i)}(1:t-1)$ for $t \leq m$,
and the subscript $s$ at the posterior pdf $p_s$ indicates that we sample from values for state $s$.
We generate at least one sample for each state $s \in \overline{\sigma^{(i)}(1:t-1)}$.
This means that the single particle $x_{\sigma(1:t-1)}^{(i)}$ is multiplied and extended to several follower particles
$x_{\sigma(1:t-1),s}^{(i)}$.
Consequently,
at iteration $t<m$ particle $(i)$ has $m-t+1$ followers.
Each follower is a sample from a different dimension of the state vector (i.e., represents a location of a different puzzle piece).
Going back to our toy puzzle example, we recall that
the current state vector of particle $(i)$ in Fig.~\ref{fig:puzzle}(c) at time $t=2$ is
$x^{(i)}_{\sigma(1:2)} = (1, 2)$, where $\sigma^{(i)}(1:2)=(3,1)$.
For sampling at time $t=3$, we have $\overline{\sigma^{(i)}(1:t-1)} = (2, 4, 5, 6)$.
Consequently, we sample four followers of particle $(i)$ in \eqref{eq:sampling2.5},
one for each state $s=2, 4, 5, 6$, where $x_{2}^{(i)}$ is the sampled location of puzzle piece 2, $x_{4}^{(i)}$ is the sampled location of puzzle piece 4, and so on.


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 samples form the same state,
since there is a unique state and only one state at time $t$,
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_{\sigma(1:t-1)}^{(i)}$ for each particle $(i)$.



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

We outline in this section the proposed {\bf PF for unordered observations (PFUO)} algorithm.
The key change is in the importance sampling step in comparison to the standard PF algorithm presented in Section \ref{sec:background}.
The other two steps are adjusted to accommodate this key change.
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 \overline{\sigma^{(i)}(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_{\sigma(1:t-1)}^{(i)}),
\end{align}
and set $x_{\sigma(1:t)}^{(i,s)} = (x_{\sigma(1:t-1)}^{(i)}, x_{s}^{(i)})$ and $\sigma^{(i,s)}(t)=s$,
which means that $\sigma^{(i,s)}(1:t) = (\sigma(1:t-1), s)$.
As stated before, we drop the superscript $(i,s)$ in $x_{\sigma(1:t)}^{(i,s)}$,
since it is already present as the particle index.


\noindent 2) \textbf{Importance weighting/evaluation:}
An individual importance weight is assigned
to each follower particle $x_{\sigma(1:t)}^{(i,s)}$ according to
\begin{align}\label{eqn:partweightv2}
& w(x_{\sigma(1:t)}^{(i,s)}) = w(x_{\sigma(1:t-1)}^{(i)}) p(z_s|x_{\sigma(1:t)}^{(i,s)},z_{\sigma^{(i)}(1:t-1)}),
\end{align}

\noindent 3) \textbf{Resampling:}
Sample with replacement $N$ new particles form the current set of $N \times (m-t+1)$ particles
\begin{align}  \label{eqn:set1}
& \{ x_{\sigma(1:t)}^{(i,s)} |  \ i=1, \dots, N, \ s \in \overline{\sigma^{(i)}(1:t-1)} \}.
\end{align}
according to the weights.
We obtain a set of new particles $x_{\sigma(1:t)}^{(i)}$ for $i=1, \dots, N$,
i.e., we assign new indices $i=1, \dots, N$ to the new particles.
We also renormalize their weights to sum to one.
This is a variant of 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 observe that the particle weight evaluation in \eqref{eqn:partweightv2} is analog to \eqref{eq:PFn4} in that
the conditional probability of observation $z_s$ is a function of
two corresponding sequences of observations and states plus the state $x_s$.
The only difference is that the sequences may not be in the original index order but in the order determined by the permutation $\sigma^{(i)}$
applied to $1:t-1$.



When $t=m$ all state dimensions are present in the state vector $x_{\sigma(1:m)}^{(i)}$ for each particle $(i)$.
Hence we can reorder the sequence of state dimensions $\sigma^{(i)}(1:m)$ to form the original order $1:m$,
i.e., we apply the inverse of $\sigma^{(i)}$ to $\sigma^{(i)}(1:m)$.
We denote the so obtained state vector with $x_{1:m}^{(i)}$,
it is composed form the same state values as $x_{\sigma(1:m)}^{(i)}$ but sorted according to the original state indices $1:m$.
We set $w(x_{1:m}^{(i)}) =w(x_{\sigma(1:m)}^{(i)})$, i.e., the particle weight remains invariant to the sorting.
The following theorem guarantees that PFUO computes an approximation to \eqref{eq:PFn1}. \\

\noindent \textbf{Theorem 2.} \
For $t=m$ the particles $x_{1:m}^{(i)}$ with $i=1, \dots, N$
provide an approximation to the posterior distribution $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}
\noindent  \textbf{Proof of Theorem 2.} \
Let us fix a single particle index $(i)$.
We have that
\begin{align}\label{eq:PFperm1}
& w(x_{1:m}^{(i)}) =w(x_{\sigma(1:m)}^{(i)}) = w^{(i)} = \frac{p(x_{\sigma(1:m)}^{(i)}|Z)}{q(x_{\sigma(1:m)}^{(i)}|Z)}.
\end{align}
To see this, let us take the state dimension order $\sigma^{(i)}(1:m)$ as the initial order of state dimensions.
It will then not change during the computation of particle $(i)$, i.e., the recursive weight update formula
simply follows the standard weight update of PF.
Consequently, $x_{\sigma^{(i)}(1:m)}^{(i)}$ is a weighted sample from $p(x_{\sigma^{(i)}(1:m)}^{(i)}|Z)$.

Now $p(x_{\sigma^{(i)}(1:m)}^{(i)}|Z)$ is the joint distribution of $m$ random variables, and as such the order of the random variables is not relevant.
(This follows from the fact that a joint probability is defined as the probability of the intersection of the sets
representing events corresponding to the value assignments of the random variables,
and set intersection is independent of the order of sets.)
Consequently, we have
\begin{align}\label{eq:PFperm2}
p(x_{\sigma^{(i)}(1:m)}^{(i)}|Z) = p(x_{1:m}^{(i)}|Z).
\end{align}
The equation \eqref{eq:PFperm2} is true for every particle $(i)$, although we may have different permutations $\sigma^{(i)}(1:m)$ for different particles.
Consequently, $x_{1:m}^{(i)}$ is a weighted sample from $p(x_{1:m}|Z)$ for every $i = 1, \dots, N$, which proves the theorem.

We stress again that the particular order of dimensions $\sigma^{(i)}(1:m)$ is extremely important for the proposed PFUO.
However, once we have this state vector,
we can simply present the state vector in the original order of dimensions.

The Sampling Importance Resampling (SIR) replaces current particles
with $N$ new particles with the weights normalized to sum to one, 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 orders of states.
Since the weights of $w(x_{\sigma(1:t)}^{(i,s)})$ are determined by the corresponding order of observations $z_{\sigma^{(i)}(1:t-1)}$,
and the resampling uses the weights to selects new particles $x_{\sigma(1:t)}^{(i)}$,
the resampling determines the order of state dimensions.
Consequently, the order of state dimensions is heavily determined by their corresponding observations,
and this order may be different for each particle $(i)$.
This is in strong contrast to the classical PF, where observations are considered only in one order $Z$.


In order to utilize the derived PF algorithm to solve the jigsaw puzzle problem,
we need to define the proposal pdf \eqref{eq:unsamp} and the conditional pdf of a new observation \eqref{eqn:partweightv2}.
Both are defined in the next section.


\section{Implementation Details} \label{sec:impdet}
Given is an rectangular grid of $m$ empty squares $G=\{g_1, \dots, g_m\}$,
as illustrated in Fig.~\ref{fig:puzzle}(c) for $m=6$,
where each $g_i$ is the index of a square in the grid.
In order to solve an image jigsaw puzzle with the puzzle pieces $P= \{1, \dots, m \}$,
we need to assign locations on $S$ to the puzzle pieces in $P$.
The observation associated with each puzzle piece is the local image depicted on it,
i.e., $z_i$ is a $K \times K \times 3$ matrix of pixel color values for $i \in P$
and $Z=\{z_1, \dots, z_m \}$ is the set of observations.

A state value assignment at time $t \leq m$ of our PF algorithm
is given by a state vector $x_{\sigma(1:t)} =(x_{\sigma(1)}, \dots, x_{\sigma(t)})$,
where each $\sigma(i) \in P$ is an index of a puzzle piece for $j=1, \dots, t$ 
and $x_{\sigma(j)} \in S$
meaning that puzzle piece $\sigma(i)$ is placed on the grid square with index $x_{\sigma(i)}$.
The vector of corresponding observations $z_{\sigma(1:t)} =(z_{\sigma(1)}, \dots, z_{\sigma(t)})$
represents the local images of the puzzle pieces.
In this section we drop the particle index $(i)$, since all definitions apply to every particle.

We define now an affinity matrix $A$ representing the compatibility of the local images on the puzzle pieces.
It is a 3D matrix with the third dimension being an adjacency type, 
since two puzzle pieces can be adjacent in four different ways:
left/right, right/left, top/bottom, and bottom/top, which we denote with LR, RL, TB, and BT.
Hence the dimension of $A$ are $m \times m \times 4$.

In order to be able to compare our experimental results to the results in \cite{PuzzleMIT2010}
we define $A$ following the definitions in \cite{PuzzleMIT2010}.
They first define the dissimilarity-based compatibility $D$.
$D$ measures dissimilarity between puzzle piece images $z_j$, $z_i$ by summing the squared $LAB$ color difference along the boundary,
e.g., the left/right (LR) dissimilarity is defined as
\begin{equation} 
D(j, i, LR) = \sum_{k=1}^K\sum_{c = 1}^3 (z_j(k,u,c) - z_i(k,v,c))^2,
\end{equation}
where $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 weight of the LR connection is given by
\begin{equation} \label{eq:binpot}
A(j,i, LR) =  \exp(- \frac{D(j,i,LR)}{2 \sigma^2}),
\end{equation}
where $\sigma$ is adaptively set as the difference between the smallest and the second smallest $D$ values
between puzzle piece $i$ and all other pieces in $P$,
see \cite{PuzzleMIT2010} for more details.


Given is a vector $x_{\sigma(1:t)}$ of state values (i.e., grid locations) at time $t$.
Let $\mathcal{C}(\sigma(j)) \subset \{ LR, RL, TB, BT \}$ denotes a set of possible connections of puzzle piece $\sigma(j)$ for $j=1, \dots, t$,
e.g., if the left and top sides of $\sigma(j)$ are adjacent to the other puzzle pieces $\sigma(1:t)$
and the right and bottom sides are adjacent to empty grid squares,
then $\mathcal{C}(\sigma(j)) \subset \{ RL, BT \}$.
We can extend this definition to the whole state vector by 
$$\mathcal{C}(\sigma(1:t)) = \bigcup \{ \mathcal{C}(\sigma(j)) \ | \ j=1, \dots, t \}.$$
The proposal distribution $p_s(\cdot |x_{\sigma(1:t-1)}) : G \to \mathbb{R}$ is assigns
a probability of placing puzzle piece $s$ to each grid square $x$.
$p_s(x |x_{\sigma(1:t-1)}) =0$ if $x$ is occupied or is not adjacent to any square in $\sigma(1:t-1)$.
Now say $x$ is adjacent and is to the right of grid square $x_{\sigma(j)}$ for some $j=1, \dots, t$.
Then
\begin{align} \label{eq:propP1}
p_s(x |x_{\sigma(1:t-1)}) \propto A(s,\sigma(j),RL).
\end{align}
Intuitively, this probability is proportional to the LR similarity between
puzzle pieces $s$ and $\sigma(j)$.

Now let $x_s$ be a sample from \eqref{eq:unsamp} at time $t$, 
and as above $x_s$ is adjacent and is to the right of grid square $x_{\sigma(j)}$ for some $j=1, \dots, t$.
Then
\begin{align}\label{eq:partweightv21}
& p(z_s|x_{\sigma(1:t)},z_{\sigma(1:t-1)}) \propto A(s,\sigma(j),RL),
\end{align}
where we recall that $\sigma(t)=s$ and $\sigma(1:t)=(\sigma(1:t-1), \sigma(t))$.
The definitions for other adjacency types $LR, TB, BT$ are analogous.

To summarize the proposal distribution is a function of how well puzzle piece $s$ fits to the already placed pieces
and assigns the probability of placing $s$ to all grid squares,
while in the evaluation we already know the grid location of puzzle piece $s$ as well as its adjacent square.
We then use this information to compute the evaluation probability according to $A$.
Hence, both the proposal and evaluation of a given particle are functions of how well adjacent pieces fit together
following the order in which the pieces has been added.
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}.
The comparison is presented in the next section.

For a given image jigsaw puzzle with $m$ pieces, the time complexity of the proposed inference framework is
$O(m^2 N)$, where $N$ is the number of particles.


\section{Experimental Results} \label{sec: Ex}

We compare the image jigsaw puzzle solutions obtained by the proposed PF inference framework
to the solutions of the loopy believe propagation used in \cite{PuzzleMIT2010} under identical settings.
We ran the software released by the authors of \cite{PuzzleMIT2010} to obtain their results.
The results are compared on the dataset provided in \cite{PuzzleMIT2010}, which we call MIT Dataset.
It is composed of 20 images.
In addition, we also consider an extended dataset composed of 40 images, i.e., we added 20 images.
As we will see below the results of both methods on the original and extended datasets are comparable.

The experimental results in \cite{PuzzleMIT2010} are conducted in two different settings:
with and without any prior on the target image layout.
In \cite{PuzzleMIT2008} the prior 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.
We focus on the results without any prior of the image layout.
Consequently, we focus on 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 as defined in Section \ref{sec:impdet}.

In the probabilistic framework in \cite{PuzzleMIT2010}, each grid location is assigned a puzzle piece.
In our PF framework, it is more natural to assign a grid location to a puzzle piece.
The solutions of both methods are equivalent, since a final puzzle solution is a set of $m$ pairs composed of
(puzzle piece, grid location), where $m$ is the number of pairs.
We call such pairs the solution pairs.


Three types of evaluation methods are used in \cite{PuzzleMIT2010}. 
Each focuses on different aspects of the quality of the obtained puzzle solutions.
The most natural and strict one is \textbf{Direct Comparison}.
It simply computes the percentage of correctly placed puzzle pieces,
i.e., for a puzzle with $m$ pieces, Direct Comparison is the number of correct solution pairs divided by the total number of solution pairs $m$.

The second measure is \textbf{Cluster Comparison}. It is less strict than Direct Comparison.
It tolerates an assignment error as long as the puzzle piece 
is assigned to a location that belongs to a similar puzzle piece.
The puzzle pieces are first clusters into clusters of similar pieces 
and only the cluster labels are used in the count of the number of correct solution pairs.

Moreover, due to lack of prior knowledge of target image, the reconstructed image may shift compared to the ground truth image. 
Therefore, the third measure called \textbf{Neighbor Comparison} 
is used to evaluate the label consistency of adjacent puzzle pieces independent of their grid location,
i.e., the location of two adjacent puzzle pieces is considered correct
if two puzzle pieces are left-right neighbors in the ground truth image and they are also left-right neighbors in he inferred image
and similarly for the other three possible adjacency relations.
Neighbor Comparison is the fraction of correct adjacent puzzle pieces. 
This measurement do not penalize the accuracy as long as the adjacent patches in original image are adjacent in the reconstructed image.


The results on the MIT Dataset are shown in \ref{Table: ExMIT108} and on the extended dataset in Table \ref{Table: Ex2MIT108}.
We stress that the proposed PF inference framework outperforms the loopy believe propagation in all three performance measures.
Moreover, the reconstruction accuracy of the original images by our algorithm is three times better 
according to the most natural measure of Direct Comparison.

In order to demonstrate that the considered image jigsaw puzzle problem is also very challenging to humans,
we show some example results in Fig. \ref{fig:Comp}.
There we show the original images, but we stress again in our experimental evaluation, the original images are not provided.
Fig. \ref{fig:Comp} also demonstrates that the reconstructed images obtained by the proposed algorithm compare 
very favorably to the results of \cite{PuzzleMIT2010}.

Both methods are initialized with one anchor patch, i.e., with one correct (puzzle piece, grid location) pair.
We always assign a correct image patch to the top left corner of the image.
In all our experiments we divide each test image into 108 square patches resulting in $m=108$ puzzle pieces.
Of course, the more more patches, the smaller are the images on the puzzle pieces, thus, larger number of patches makes the problem harder.
Since the performance of \cite{PuzzleMIT2010} on 108 puzzle pieces is already very low,
in particular, as measured by Direct Comparison,
a larger number of pieces does not seem to lead to any meaningful results.
%This is also confirmed by the results reported in Fig. 8(a) in \cite{PuzzleMIT2010}.
%Since no statistics is reported in \cite{PuzzleMIT2010} for a system without any image evidence,
%we focus on the results for a system with noisy evidence,
%which are definitely not worse than the results without any image evidence.
%There the value of Direct Comparison with one anchor patch for 221 jigsaw pieces is below $20\%$.


\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.}
\label{Table: ExMIT108}
\end{table}



\begin{table}
\begin{center}
\begin{tabular}{|c|c|c|}
\hline
 & \cite{PuzzleMIT2010} & Our algorithm \\
 \hline
 Direct Comparison &0.2137& 0.7097\\
 \hline
 Neighbor Comparison &0.6458 & 0.8770\\
 \hline
 Cluster Comparison &0.4500 & 0.8018\\
 \hline
\end{tabular}
\end{center}
\caption{Experimental results on the extended dataset.}
\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 images in the first row are the original images.
 The jigsaw puzzle solution of Cho et. al \cite{PuzzleMIT2010} are shown in the second row. 
 The solutions of our algorithm are in the third row.}
\label{fig:Comp}
\end{figure*}


In order to demonstrate the dynamic of the proposed PF inference,
we show reconstructed images of the best particle at different times (iterations) in Fig. \ref{fig:PF_Iter}.
The puzzle pieces of the best particles at different iterations are different.
This illustrates the facts that different dimensions of the state vector are considered by different particles
and the orders of considered dimensions are different. 



\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 images of the best particle at different iterations.}
\label{fig:PF_Iter}
\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}
We introduce a novel inference framework for solving labeling or assignment problems expressed as dynamic state estimation.
Our key contribution is an extension of the PF framework to work with unordered observations.
Intuitively, this means that the weighted particles explore different orders of label assignment hypothesis.
We prove that so obtained weighted particles represent samples form the target, posterior distribution of label assignments.
We evaluate the labeling performance on a problem of image jigsaw puzzles, which is very challenging even for humans.
As the experimental results demonstrate, the proposed PF inference significantly outperforms the loopy belief propagation.

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

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


\end{document}
