% Questo file .tex è compilabile su Kile per Fedora 10 con latex+dvi2ps+ps2pdf  oppure su Mac OS X 10.6 Snow Leopard con comando LaTeX e script TeX + DVI 
% per produrre il .pdf finale (solo immagini in eps).

\documentclass[12pt,journal,compsoc]{pcks/IEEEtran}
\usepackage[nocompress]{cite}
\usepackage[utf8]{inputenc}
\setlength{\parindent}{0pt}
\usepackage[T1]{fontenc}

\usepackage{indentfirst}
\usepackage{float}
\usepackage{pcks/floatflt}
\usepackage{verbatim}
\usepackage{textcomp}
\usepackage[final]{listings}

\usepackage{multirow}

\usepackage[dvips]{graphicx}
\usepackage[dvipsnames]{xcolor}
\usepackage{epsfig}
\usepackage{psfrag}
\usepackage[caption=false,font=normalsize,labelfont=sf,textfont=sf]{subfig}
\usepackage{fixltx2e}
\usepackage[dvips]{thumbpdf}

\usepackage{amssymb}
\usepackage[cmex10]{amsmath}
\interdisplaylinepenalty=2500
\usepackage{amsbsy}
\usepackage{amsfonts}
\usepackage{pcks/faktor}
\usepackage{pcks/chemarrow}
\usepackage{mathtools}
\usepackage{verbatim}
\usepackage{eucal}
\usepackage{mdwmath}
\usepackage{mdwtab}
\usepackage{eqparbox}
\usepackage{marvosym}

\usepackage{pdflscape}

\usepackage{url}
\newcommand\ALZhyperrefoptions{bookmarks=true,bookmarksnumbered=true,
pdfpagemode={UseOutlines},plainpages=false,pdfpagelabels=true,
,colorlinks=true,
linkcolor={black},citecolor={black},pagecolor={black},filecolor=red,urlcolor={black}
}
\usepackage[\ALZhyperrefoptions,breaklinks=true,dvips,hyperindex]{hyperref}
\usepackage{breakurl}
\hypersetup{
pdftitle={Camera Simulator and Multi Target Tracking},
pdfsubject={Camera Simulator and Multi Target Tracking},
pdfauthor={D'Elia, Tavian Zugno},
%pdfpagemode=FullScreen
pdfkeywords={Camera, Simulator, Multi, Target, Tracking, Stable, Marriage, Data, Association, Motion, Capture},
pdfstartview=FitH,backref=true,pagebackref=true,
}


\usepackage{pcks/algorithmicx/algorithmicx}
\usepackage{algorithm}
\usepackage{pcks/algorithmicx/algpseudocode}
\newcommand{\INPUT}{\State \textbf{input:} }
\newcommand{\OUTPUT}{\State \textbf{output:} }

% precisazione sillabazioni
\hyphenation{tracking} 
\hyphenation{zu-gno}


\usepackage{pgf,pgfsys,pgffor}
\usepackage{tikz}
\usetikzlibrary{arrows,decorations.pathmorphing,backgrounds,fit,mindmap,calc,through,scopes,fadings,positioning,automata,calendar,shapes,er,matrix,folding,patterns,petri,plothandlers,plotmarks,shadows,topaths,through,trees}


\numberwithin{equation}{section}
 \numberwithin{figure}{section} % (per le figure)
 
\newtheorem{thm}[equation]{Theorem} % per avere la stessa denominazione per teoremi / formule / assiomi
\newtheorem{prop}[equation]{Proposition}
\newtheorem{lemma}[equation]{Lemma} % metti equation quindi al posto di thm
\newtheorem{cor}[equation]{Corollary}


\newcommand{\ud}{\mathrm{d}} % differenziale per integrali

\newcommand{\trasp}{^{\mathrm{T}}} % T trasposta dritta

% \DeclareMathSizes{12}{12}{8}{16} % cambia la grandezza dei font nelle equazioni
        

%  qui si usa mcode.sty per includere i file .m di MatLab, ottimizzato all'uopo
% per i colori del bordo si vedano il pdf relativo
\usepackage[numbered,framed]{pcks/mcode}

% per avere più sottosezioni, qui 5
\usepackage{pcks/subsections}
\setcounter{secnumdepth}{5}
\setcounter{tocdepth}{5}


% nuovi comandi per argmin ed argmax
\newcommand{\argmax}[1]{\underset{#1}{\operatorname{argmax}}\;}
\newcommand{\argmin}[1]{\underset{#1}{\operatorname{argmax}}\;}






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%   WARNING:  

% ``no in outer par mode'' OR ``Infinite glue shrinkage found in box being split'' are errors due to the fact that figures/algorithms/biography cannot be divided in 2 coloums.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



\begin{document}
\title{Simulation and Multi-Target Tracking}

\author{Edoardo~D'Elia,
        Carlo~Tavian,
        and~Alberto~Zugno,% <-this % stops a space
\IEEEcompsocitemizethanks{
\IEEEcompsocthanksitem Edoardo~D'Elia, ID no. 603406, E-mail: deliaedo@dei.unipd.it
\IEEEcompsocthanksitem Carlo~Tavian, ID no. 607367, E-mail: tavianca@dei.unipd.it
\IEEEcompsocthanksitem Alberto~Zugno, ID no. 607378, E-mail: zugnoalb@dei.unipd.it}% <-this % stops a spac,e
\thanks{Students at the Department of Information Engineering, University of Padua. 
Report presented February 15 2010.}}


%\IEEEpubid{\makebox[\columnwidth]{\hfill 0000--0000/00/\$00.00~\copyright~2007 IEEE}%
%\hspace{\columnsep}\makebox[\columnwidth]{Published by the IEEE Computer Society\hfill}}

%\IEEEspecialpapernotice{(Invited Paper)}

\IEEEcompsoctitleabstractindextext{%


\begin{IEEEkeywords}
Camera Simulator, Multiple Target Tracking, SMP, MCMCDA, Data Association, Motion Capture
\end{IEEEkeywords}}

\maketitle



%\title{\textbf{Simulation and Multitarget Tracking}}
%\author{
%\IEEEauthorblockN{}
%  \IEEEauthorblockA{
%    
%  }

%\and

%\IEEEauthorblockN{}
%  \IEEEauthorblockA{
%    matricola:  \\
%    
%  }

%\and

%\IEEEauthorblockN{}
%  \IEEEauthorblockA{
%    matricola:  \\
%    
%  }
%}

%\maketitle

\begin{abstract}
Our purpose is to create a camera simulator and a multiple target
tracking system for motion capture applications.
Starting from literature and state of the art techniques, we focused upon the
Markov Chain Monte Carlo (MCMC) method to obtain a MAP estimation of marker's 
tracks. In particular, this method is based on the well known
Metropolis-Hastings algorithm. This method led to very good results (even with
false alarms and missing measurements), although it required long computation times.
Searching for innovative solutions, we implemented also an algorithm to resolve
the Stable Marriage Problem, and we adapted it to the perform data association.
This implementation is very fast and results accurate, as drawbacks false alarms and missing measurements are not taken into consideration.
For both approaches with provide ideas for possible improvements to get shorter
computation times (MCMCDA) or to take into account disturbance factors (SMP).
\end{abstract}
\section{Introduction}

\subsection{Goals}

\subsubsection{Multi-Target Tracking}

Multi-target tracking is a well known problem in literature. Given a varying number of targets $k_t$ and measurements $n_t$, the purpose of multi-target tracking is to associate measurements with the targets which have generated them. There are many possible disturbance factors, such as false alarms, corrupted or missing reports, non-constant sampling time, etc.

There are a lot of fields in engineering where multiple-target tracking is widely used, such as motion capture, surveillance systems, computer vision and sensor network. Our interest is focused on motion capture area. In our motion capture application we have to track a variable number of markers that are indiscernible, so the only information we have about them is their position at some sampling time. Given that, our principal problem become the data association: we must distinguish markers resting only upon their dynamics; whereas other applications like computer vision can rely upon some descriptor of tracked objects like color, shape, etc. that makes data association problem more reasonable. Motion capture is not the only scenario where data association is the principal problem: for example, in a radar tracking of unidentified vehicles (surveillance) we don't have any data other than their position.
A few number of algorithms were developed for this purpose; the most widely used are JPDA~\footnote{Joint Probabilistic Data Association} for a fixed number of target and the MHT~\footnote{Multiple Hypothesis Testing} for a variable number of targets. These algorithm are used since the Eighties and are generally slow. A new faster method is \emph{Markov Chain Monte Carlo Data Association}~\cite{MCMCDA}. This one outperforms MHT and JPDA in efficiency and accuracy under many conditions. Since a lot of applications in these fields require to be real-time we have implemented the MCMCDA algorithm and we have tested it on some synthetic data. We have also implemented a data association algorithm based on the SMP~\footnote{Stable Marriage Problem}; this is notably faster than the previous one but as a drawback, it works only with fixed number of targets and without false alarms or undetected measurements (unless we introduce some improvements that we will discuss later).

\subsubsection*{Camera Simulator}

In motion capture scenarios a camera simulator should be useful to test various algorithms and to generate synthetic data. In literature there are a bunch of camera models, but only the pin-hole camera model is widely used. We came up with an effective camera simulator which has been used also by other groups which participates in same project.

\subsection{Results}

We managed to implement a MCMCDA algorithm to solve the data association problem. This implementation has turned out quite effective but very slow whereupon not applicable to a real-time application with many targets. We have also adapted the Stable Marriage to Data Association problem and implemented an algorithm that is very effective under some conditions and may be applied in real-time applications.

\newpage
\section{Literature and State of the Art}

Multiple-target tracking (MTT) is an essential component of surveillance-related systems. Its application areas include missile defense, air traffic control, and --- in our case --- motion capture techniques. A general formulation of the problem assumes an unknown and variable number of targets that moves with continuity in a given region. In the single-sensor version, the states of these targets are sampled by the sensor and the noisy measurements are provided to the tracking system. The detection probability is less than one then the targets could be undetected at some sampling time. In addition, there are false reports of possible targets or clutter measurements that arise independently from targets of interest. A primary task of an MTT system is data association, i.e. partitioning the measurements into disjoint sets each generated from a single source (target or clutter). Secondary goal the states estimation based on measurements originating from interesting targets. The data association problem could be formulated in several ways. In the single scan data association, the raw measurements are processed one scan at a time and the target states get updated accordingly. Alternatively several sets of measurements could be collected and processed together in batch mode –- this is the multi scan data association.

Several methods now exist to handle the data association problem. These may roughly grouped into two categories: Bayesian and non-Bayesian. Among the Bayesian methods, there is the well known Joint Probabilistic Data Association Filter (JPDA)~\cite{B-S}, which is a single scan filter where the existent targets' states get updated on the basis of the latest set of measurements (scan). Data association is handled by summing over the probabilities of all feasible partitions, where the targets can't share a measurement and each target could be the source of at most one measurement per scan. Another well known approach is the multiple hypothesis tracker (MHT)~\cite{MHT}, in which each hypothesis make a target-to-measurement match and, little by little observations come up, a new set of hypothesis is formulated increasing the number of the previous ones. The hypothesis with the highest posterior is returned as a solution. MHT is capable of initiating and terminating tracks. However, the number of hypothesis involved in computation grows exponentially over the time. Thus, in order to overcome this computational complexity, a certain number of pruning and clustering methods must be used at expense of optimality.

The non-Bayesian approaches are characterized by hard measurement-to-track association, such that some cost function is minimized. Then the problem may be reformulated as an integer linear programming problem or, to be more precise, as a multidimensional assignment problem, which is NP-hard. Therefore, for the multi scan data association, one should invoke some approximations schemes such as Lagrangian relaxation techniques that relax some constraints and solve the relaxed problem trough linear programming methods~\cite{LagRelax}.

Another option to solve the multi scan data association problem is by using stochastic search methods. In~\cite{MCMCDA} the problem was solved by applying the Markov Chain Monte Carlo (MCMC) method to obtain the partition with maximum posterior. Using the Metropolis-Hastings algorithm, the authors proposed a set of moves to modify the measurements partition, such that sampling from the posterior distribution was possible after few thousands of moves. They showed a remarkable performance of the algorithm in comparison to the MHT method in terms of accuracy of the solution and running time.

\subsection{Joint Probabilistic Data Association}

\begin{figure}
  \centering
  \includegraphics[width=.45\textwidth]{graphs/edo_tracking.pdf}
  \label{fig:matter}
  \caption{Two targets with a common measurement.}
\end{figure}

The Joint Probabilistic Data Association is an extension of the PDA method, which dealt with a single target in clutter, to the situation where there is a \emph{known} number of targets in clutter. When there are several targets in the same neighborhood, measurements from one target can fall in the validation region of a neighboring target (see \figurename~\ref{fig:matter}). This can happen over several sampling times and acts as ``persistent interference''. Since the PDA algorithm models all the incorrect measurements as ``random interference'', with independent uniform spatial distributions, its performance can degrade significantly when the existence of a neighboring target gives rise to interference that is not correctly modeled.

Similarly to the PDA algorithm, the joint probabilistic data association (JPDA) algorithm computes the probabilities of association of only the latest set of measurements $Y_k$ to the various targets --- this is a non-backscan approach. Each target has a dynamic and a measurement model. The models for the various targets do not have to be the same. The state estimation is done for each target as in PDA, i.e. a Kalman filtering and a measurement validation, but the measurement-to-target association probabilities are computed in a different way --- jointly across the targets.

The key to the JPDA algorithm is the evaluation of the conditional probabilities of the following joint association events pertaining to the current time\footnote{In order to reduce notational complexity, the time index $k$ will be omitted from $\theta, \boldsymbol{\theta}, \omega, \boldsymbol{\Omega}$, etc. whenever it can be done without causing confusion.} $k$:
\begin{equation*}
  \boldsymbol{\theta} = \bigcap_{j=1}^{m_k} \theta_{jt_j},
\end{equation*}
where
\begin{multline*}
  \theta_{jt} = \{ \text{measurement} \: j \: \text{originated from target} \: t \}, \\
  j=1,\dots,m_k; \quad t=0,1,\dots,T
\end{multline*}
and $t_j$ is the index of the target to which measurement $j$ is associated in the event under consideration.

The derivation of the probabilities of joint events is done using Bayes' rule as follows. The conditional probability of the joint association event $\boldsymbol{\theta}(k)$ at time $k$ (with the time argument now indicated explicitly) is
\begin{equation*}
  \begin{split}
    P[\boldsymbol{\theta}(k)|Y^k] &= P[\boldsymbol{\theta}(k)|Y(k),Y^{k-1}] \\
    & = \frac{1}{c}P[Y(k)|\boldsymbol{\theta}(k),Y^{k-1}]P[\boldsymbol{\theta}(k)|Y^{k-1}] \\
    &= \frac{1}{c}P[Y(k)|\boldsymbol{\theta}(k),Y^{k-1}]P[\boldsymbol{\theta}(k)]
  \end{split}
\end{equation*}
where $c$ is a normalizing constant. Note that the conditioning on $m_k$ is implicit in the vector $\boldsymbol{\theta}(k)$ and that the irrelevant prior conditioning term has been omitted in the last line of the equation. The PDF on the right-hand side is
\begin{equation*}
  P[Y(k)|\boldsymbol{\theta}(k),Y^{k-1}] = \prod_{j=1}^{m_k}P[\boldsymbol{y}_j(k)|\theta_{jt_j}(k),Y^{k-1}].
\end{equation*}
The conditional PDF of a measurement given its origin is assumed to be
\begin{multline*}
  P[y_j(k)|\theta_{jt_j}(k),Y^{k-1}] = \\
  = 
  \begin{cases}
    \mathcal{N}_{t_j}(y_j(k)), & \text{if}\: \tau\{\boldsymbol{\theta}(k)\} = 1 \\
    V^{-1}, & \text{if}\: \tau\{\boldsymbol{\theta}(k)\} = 0,
  \end{cases}
\end{multline*}
where a measurement associated with target $t_j$ has the Gaussian PDF
\begin{equation*}
  \mathcal{N}_{t_j}(y_j(k)) = \mathcal{N}_{t_j}(y_j(k);\hat{y}_j^t(k|k-1),S_j^t(k))
\end{equation*}
and $\hat{y}_j^t$ denotes the predicted measurement for target $t_j$, with associated innovation covariance $S_j^t$. Measurements not associated with any target are assumed uniformly distributed in the surveillance region of volume $V$. In addition the marginal association probabilities are obtained from the joint probabilities by summing over all the joint events in which the marginal event of interest occurs.

\subsection{Track-Splitting Filter}

In this procedure, following initialization, the track is ``split'' into separate hypothesis tracks at time $k = 1$, one for every measurement at that falls in the validation region around the location $\hat{y}(1|0)$ where the measurement is expected, the track is split. Thus for each measurement a separate updated state and covariance are computed via the Kalman filter equations and propagated forward to yield another validation region at $k = 2$. For each new validation region at $k = 2$, the procedure is then repeated.

Since the number of branches into which the track is split can grow exponentially, the likelihood function of each split track is computed and the unlikely ones are discarded. Denote the $l^{th}$ sequence of measurements up to time $k$ as
\begin{equation*}
  Y^{k,l} = \{y_{i_{1,l}}(1),\dots,y_{i_{k,l}}(k)\}
\end{equation*}
where $y_i(j)$ is the $i^{th}$ measurement at time $j$. The likelihood function of this sequence being a track, i.e., of the event that its components originated from the same target, denoted
\begin{equation*}
  \theta^{k,l} := \{Y^{k,l} \text{is a correct track}\}
\end{equation*}
is the joint probability density function
\begin{equation*}
  \Lambda(\theta^{k,l}) = P[Y^{k,l}|\theta^{k,l}] = P[y_{i_{1,l}}(1),\dots,y_{i_{k,l}}(k)|\theta^{k,l}].
\end{equation*}
The track initialization information (initial estimate) is subsumed. The joint PDF can be expressed as
\begin{equation*}
  \Lambda(\theta^{k,l})  = \prod_{j=1}^k P[y_{i_{j,l}}(j)|Y^{j-1},\theta^{k,l}].
\end{equation*}
Under the linear Gaussian assumptions,
\begin{align*}
  P[y(j)|Y^{j-1},\theta^{k,l}] &= \mathcal{N}(y(j);\hat{y}(j|j-1),S(j)) \\
  &= \mathcal{N}(e(j);0,S(j)),
\end{align*}
where the last form follows from the definition of the innovations $e(j)$ and the particular form of the Gaussian density. The subscripts on the measurements $y$ have been omitted to avoid making $e$ and $S$ very complicated. From the above equations it follows
\begin{align*}
  \Lambda(\theta^{k,l}) &= \prod_{j=1}^k{|2\pi S(j)|}^{-\frac{1}{2}}\exp\{-\frac{1}{2}e^T(j)S^{-1}(j)e(j)\} \\
  &= c_k\exp\{-\frac{1}{2}e^T(j)S^{-1}(j)e(j)\}.
\end{align*}
Note that this likelihood function assumes implicitly that the target detection probability is unity.

The modified log-likelihood function, without explicitly indicating the index $l$, is
\begin{equation*}
  \lambda(k) = -2\log\left[\frac{\Lambda(\theta^{k,l})}{c_k}\right] = \sum_{j=1}^ke^T(j)S^{-1}(j)e(j),
\end{equation*}
which can be also computed recursively as follows:
\begin{equation*}
  \lambda(k) = \lambda(k-1) + e^T(k)S^{-1}(k)e(k).
\end{equation*}
The last term above has a chi-square density with $n_y$ degrees of freedom and, since the innovations are independent, the log-likelihood function at time $k$ is chi-square distributed with $kn_y$ degrees of freedom. Note that it is a measure of the ``goodness of fit'' of the measurements.

The statistical test for accepting a track is that the log-likelihood function satisfies
\begin{equation*}
  \lambda(k) \leq \delta,
\end{equation*}
where the threshold $\delta$ follows from chi-square with $kn_y$ degrees of freedom,
\begin{equation*}
  P[\chi_{kn_y}^2 > \delta] = \alpha
\end{equation*}
In the above, $\alpha$ is the probability that a true track will be rejected and it is taken, typically, as $1\%$.

\subsection{Multi Hypothesis Tracker}

The Multi Hypothesis Tracker (MHT) is a method for calculating the probabilities of various data association hypothesis.

In addition to the above data association capabilities, the algorithm include the desirable characteristic of multiple-scan correlation, clustering and recursiveness. Multiple-scan correlation is the capability to use later measurements to aid in prior correlations (associations) of measurements with targets. This feature is usually found in batch-processing or track-splitting algorithms. Clustering is the process of dividing the whole set of targets and measurements into independent groups (or clusters). Instead of solving a big problem, a number of smaller problems are solved in parallel. Finally, it is desirable for an algorithm to be recursive so that all the previous data do not have to be reprocessed whenever a new data is received.

The MHT take into account of some ``sensor trouble'' by using the detection and false alarm statistics, the expected density of unknown targets, and the accuracy of the target estimates.

\begin{figure}
  \centering
  \includegraphics[width=.45\textwidth,height=.55\textheight]{graphs/schema_mht.pdf}
  \label{MHT}
  \caption{Flow diagram of multiple-target tracking algorithm.}
\end{figure}

A flow diagram of the tracking algorithm is shown in \figurename~\ref{MHT}. Most of the processing is done within the four subroutines shown in the figure. The \textsc{cluster} subroutine associates measurements with previous clusters. If two or more previously independent clusters are associated because of a measurement, then the two clusters are combined into a ``super cluster''. A new cluster is formed for any measurement not associated with a prior cluster. As part of initialization program, previously known targets form their own individual clusters.

The \textsc{hypgen} subroutine forms new data-association hypothesis for the set of measurements associated with each cluster. The probability is calculated and target estimates are updated for each hypothesis of each cluster.

Both the \textsc{cluster} and \textsc{hypgen} subroutines use the \textsc{reduce} subroutine for eliminating unlikely hypothesis or combining hypothesis with similar target estimates. Once the set of hypothesis is simplified by this procedure, uniquely associated measurements are eliminated from the hypothesis matrix by the \textsc{mash} subroutine. Tentative targets become confirmed targets if they were the unique origin of the eliminated measurement.


\subsection{Linear Assignment Problem}

Assignment problems deal with the question how to assign $n$ items (e.g. jobs) to $n$ machines (or workers) in the best possible way. They consist of two components: the assignment as underlying combinatorial structure and an objective function modeling the ``best way''.

Mathematically an assignment is nothing else than a bijective mapping of a finite set into itself, i.e. a permutation. Assignments can be modeled and visualized in different ways: in a permutation matrix  or as an adjacency matrix of a bipartite graph $G = (V, W)$, where the vertex sets $V$ and $W$ have $n$ vertices, i.e. $|V| = |W| = n$.

The set of all assignments of $n$ has $n!$ elements. We can describe this set by the following equations called assignment constraints
\begin{align*}
  &\sum_{i=1}^nx_{ij} = 1 & &\forall j = 1, \dots, n \\
  &\sum_{j=1}^nx_{ij} = 1 & &\forall i = 1, \dots, n \\
  &x_{ij} \in \{0,1\} & &\forall i,j = 1, \dots, n.
\end{align*}
Let $A$ be the coefficient matrix of the system of equations. Matrix $A$ is \emph{totally unimodular}, i.e. every square submatrix of $A$ has a determinant of value $+1$, $-1$ or $0$.

Multi-dimensional (sometimes referred as \emph{multi-index}) assignment problems (MAP) are natural extensions of the linear assignment problem. They have been considered for the first time by Pierskalla~\cite{LAP:Pier}. The most prominent representatives of this class are axial and planar $3$-dimensional assignment problems. The general formulation of the MAP is
\begin{align*}
  & \min & & \sum_{i_1=1}^n \dots \sum_{i_d=1}^n c_{i_1\dots i_d}x_{i_1\dots i_d} \\
  & \text{s.t.} & & \sum_{i_2=1}^n \dots \sum_{i_d=1}^n x_{i_1\dots i_d} = 1, \quad i_1 = 1, \dots, n \\
  & & & \sum_{i_1=1}^n \sum_{i_3=1}^n \dots \sum_{i_d=1}^n x_{i_1\dots i_d} = 1, \quad i_2 = 1, \dots, n \\
  & & & \quad \vdots \\
  & & & \sum_{i_1=1}^n \dots \sum_{i_{d-1}=1}^n x_{i_1\dots i_d} = 1, \quad i_d = 1, \dots, n \\
  & & & x_{i_1 \dots i_d} \in \{0,1\} \quad \text{for} \quad 1 \leq i_1, i_2, \dots, i_d \leq n,  
\end{align*}
with $n^d$ coefficient $c_{i_1 \dots i_d}$.

In terms of graphs a multidimensional assignment problem can be described as follows: let a complete $d$-partite graph $G = (V_1,V_2,\dots,V_d;E)$ with vertex sets $V_i$, $|V_i| = n$, $i = 1,2,\dots,d$, and edge set $E$ be given. A subset $X$ of $V = \bigcup_i^d V_i$ is a clique, if it meets every set $V_i$ in exactly one vertex. A $d$-dimensional assignment is a partition of $V$ into $n$ pairwise disjoint cliques. If $c$ is a real valued cost function defined on the set of cliques of $G = (V_1,V_2,\dots,V_d;E)$, the $d$-dimensional assignment problem asks for a $d$-dimensional assignment of minimum cost.

Multidimensional assignment problems in their general form have recently found some applications as a means to solve data association problems. More specifically, the central problem in any multi-target tracking and multi-sensor surveillance is the data association problem of partitioning the observations into tracks and false alarms in real time. General classes of these problems can be formulated as multidimensional assignment problems.


\subsection{MCMCDA}

This is the principal method we followed in our work. The algorithm is based on~\cite{MCMCDA} where the authors propose an effective way to solve the MTT problem. The idea is to use a Monte Carlo method to ``work out'' the formidable Bayesian estimation matter. To achieve this purpose they used a Metropolis-Hastings algorithm to approximate the optimal solution. More precisely a Markov chain is used for sampling instead of enumerating over all the possible association.

The Monte Carlo Markov Chain Data Association (MCMCDA) algorithm presented in this report is the on-line multiple scan version that includes detection failure, false alarms, and track initiation and termination. In addition, this scheme of filtering is particularly suitable for situation where the number of targets is unknown and changes over the time.

The solution space $\Omega$ contains association histories over multiple steps, as well as considering all possible numbers of targets at each step. The procedure features efficient mechanisms to search over this large solution space in addition to birth and death moves to add or remove tracks.

For further detail see section~\ref{Alberto}.

\subsection{Stable Marriage Problem}

Imagine you are a matchmaker, with one hundred female clients, and one hundred male clients. Each of the women has given you a complete list of the hundred men, ordered by her preference: her first choice, second choice, and so on. Each of the men has given you a list of the women, ranked similarly. It is your job to arrange one hundred happy marriages.

In this problem, we have a set of n men and n women. Each person has their own preference list of the persons they want to marry. Our job is to determine an assignment where each man is married to one and only one woman (monogamous and heterosexual).

Each man has a list of women ordered by his preference and each woman has a similarly ranked list. Every man is on every woman's list and every woman is on every man's list. The goal is to have a set of stable marriages between the men and the women.

When given two married pairs, (\Male,~\FEMALE) and (\MaleMale,~\Female), if man \Male\ prefers another woman \Female\ more than his current wife \FEMALE\ and woman \Female\ prefers \Male\ more than her current husband \MaleMale, then (\Male,~\FEMALE) is called a dissatisfied pair. The marriage is said to be a stable marriage if there are no dissatisfied pairs.

Now, our problem is to ``marry'' targets with measurements in a stable (and satisfied) way. In this report we discuss this method with more detail in section~\ref{Carlo}.

\newpage
\section{Camera simulator}

\subsection{Ideal Camera Model}

In order to build a camera simulator we have to define a suitable mathematical model for a generic camera. All camera's aim is to represent 3D points of real world in a 2D surface (the image plane): in order to do this a projection
is performed.

The most widely used camera model is the pin-hole camera, which is presented schematically in the figure~\ref{fig:schema_pinhole}.
\begin{figure}[tbp]
  \centering
  \includegraphics[width=0.45\textwidth]{graphs/pinhole2Dscheme.png}
  \caption{Pin-hole schematics}
  \label{fig:schema_pinhole}
\end{figure}
In the pin-hole model we suppose that every ray of light that impress the ccd (or film) and form the camera image has passed through a point called optical center. 
In a real camera, the fact that all the ray of light pass through the same point is clearly an approximation; nevertheless for small lens this is a very good approximation that can be improved caring of some kind of distortions.

From the knowledge of the focal length and the coordinates $X=[X_c;Y_c;Z_c]^T$ of the world point (relative to a reference frame centered at the optical center with z-axis being the optical axis) we can calculate his projection coordinates ($x=[x,y]^T$) on the image plane. 

\begin{figure}[tbp]
  \centering
  \includegraphics[width=0.45\textwidth]{graphs/updown.pdf}
  \caption{Real pin-hole camera upside-down effect}
  \label{fig:updown}
\end{figure}

In a real camera, the image should be upside-down with respect to the object (Figure~\ref{fig:updown}) so, for simplicity, we consider the image plane as if it would be located before the optical center.
By the property of similar triangles we have that
$$x=\frac{fX_c}{Zc};\:y=\frac{fY_c}{Zc}$$,
where the position of point $(X_c,Y_c,Z_C)$ is referred to the camera's coordinate system.
To express the world's point in to camera's coordinate system we must know the position and orientation of the camera (extrinsic parameters) and perform a rototranslation (figure~\ref{fig:rototranslation}).

\begin{figure}[tbp]
  \centering
  \includegraphics[width=0.45\textwidth]{graphs/rototranslation.png}
  \caption{Rototranslation from world to camera frame}
  \label{fig:rototranslation}
\end{figure}

Named $R$ the rotation matrix and $T$ the translation vector of the camera-centered coordinate with respect to the world coordinates
we have:
$$\left[\begin{array}{c}
X_{c}\\
Y_{c}\\
Z_{c}\end{array}\right]=R\cdot\left[\begin{array}{c}
X\\
Y\\
Z\end{array}\right]+T$$

\begin{figure}[tbp]
  \centering
  \includegraphics[width=0.45\textwidth]{graphs/ccd.png}
  \caption{Image reference set translation. Adjust the position of the optical center projection to refer to the image reference set that has the origin in it's upper left corner.}
  \label{fig:ccd}
\end{figure}

Since all the cameras used today are digital (i.e. the image is represented by a matrix of finite dimensional pixel), it becomes handy to express the position of the projected point using the pixel as measure unit as shown in figure~\ref{fig:ccd}. Moreover, in a digital
image, the origin of the coordinate system is on the top-left corner: $(x,y)$ coordinates are related with $(u,v)$ coordinate
by the following equation:

$$\begin{array}{c}
u=u_{0}+\frac{{K_{u}fX_{c}}}{Z_{c}}\\
v=v_{0}+\frac{{K_{v}fY_{c}}}{Z_{c}}\end{array}$$

Where $K_u$, $K_v$ are the pixel width and height, and $u_{0}$ $v_{0}$ are the offset of the principal point\footnote{Principal point: the projection of the
optical center on the image plane.}
respect the $(u,v)$ coordinate system ($K_u$, $K_v$, $u_{0}$, $v_{0}$ are called intrinsic parameters).

If we express the points position in homogeneous coordinates all these equations for n different points can be rearranged in matrix products:


$$
\left[\begin{array}{cccc}
su_{1} & su_{2} & \cdots & su_{n}\\
sv_{1} & sv_{2} & \cdots & sv_{n}\\
s & s & \cdots & s\end{array}\right]=K\cdot F\cdot RT\cdot P_{3dpts}=
$$
$$=\left[\begin{array}{ccc}
k_{u} & 0 & u_{0}\\
0 & k_{v} & v_{0}\\
0 & 0 & 1\end{array}\right]\left[\begin{array}{cccc}
f & 0 & 0 & 0\\
0 & f & 0 & 0\\
0 & 0 & 1 & 0\end{array}\right]\cdot
$$
%%
$$
\cdot \left[\begin{array}{cccc}
r_{11} & r_{12} & r_{13} & T_{x}\\
r_{21} & r_{22} & r_{23} & T_{y}\\
r_{31} & r_{32} & r_{33} & T_{z}\\
0 & 0 & 0 & 1\end{array}\right]
\left[\begin{array}{cccc}
x_{1} & x_{2} & \cdots & x_{n}\\
y_{1} & y_{2} & \cdots & y_{n}\\
z_{1} & z_{2} & \cdots & z_{n}\\
1 & 1 & \cdots & 1\end{array}\right]
$$

Where s is a scale factor. The matrix obtained by the product of first two matrix in the previous equation is called intrinsic parameters matrix.

\subsection{Distortions}

If we want a complete description of the camera we are going to simulate, it could be useful to take into account some kind of distortion.

\begin{figure}[tbp]
  \centering
  \includegraphics[width=0.45\textwidth]{graphs/skew.pdf}
  \caption{Skew Distortion}
  \label{fig:skew}
\end{figure}

The simpler distortion we can implement is the skew factor (figure~\ref{fig:skew}): introducing the parameter $s_w$ in the intrinsic parameters matrix
replacing $K$ matrix with $K'$:
$$K'=\left[\begin{array}{ccc}
k_{u} & s_w & u_{0}\\
0 & k_{v} & v_{0}\\
0 & 0 & 1\end{array}\right]
$$

Another kind of aberration that typically affects cameras is the radial distortion (Figure\ref{fig:radial}). This one is not linear so we can't model it within the previous matrix equations. The simplest effective model for such a distortion is:
$$
\begin{array}{c}
x=x_{d}(1+a_{1}r^{2}+a_{2}r^{4})\\
y=y_{d}(1+a_{1}r^{2}+a_{2}r^{4})\end{array}$$

\begin{figure}[tbp]
  \centering
  \includegraphics[width=0.45\textwidth]{graphs/radial.pdf}
  \caption{Two type of radial distortions: barrel and pincushion}
  \label{fig:radial}
\end{figure}

Depending on the sign of $a_{1}$ we have two opposite kind of distortion: barrel or pincushion.
Barrel distortion~\footnote{In barrel distortion, image magnification decreases with distance from the optical axis. The apparent effect is that of an image which has been mapped around a sphere (or barrel).} (illustrated in figure~\ref{fig:radial}) typically will have a positive term for $a_{1}$ where as pincushion distortion~\footnote{In pincushion distortion, image magnification increases with the distance from the optical axis. The visible effect is that lines that do not go through the center of the image are bowed inwards, towards the center of the image, like a pincushion.} will have a negative value.
This kind of distortion should be taken into account in the case of pan-tilt-zoom cameras, because appears most visibly when the widest angle (shortest focal length) is selected either with a fixed or a zoom lens.

In addition to radial distortions and skew, there are other types of lens distortion which affect the quality and geometry of an image and are not radially symmetric. This type of distortion (tangential) is not taken into account, since it is usually small compared with that due to radial lens distortion at the typical image resolutions of digital cameras.

\subsection{Visibility and Visualization}
With the previous described model all points are projected on the image plane, but not all are really visible because:
\begin{itemize}
\item the projection can be outside the image plane since it has limited dimensions,
\item the point is behind image plane
\item the point is too far from the camera
\end{itemize}
The simulator checks, while calculating the position of 3D points respect the camera frame, the distance of any point from the camera and marks as not visible all the points that are behind or too far from it.
Similarly, after the projections on an infinite image plane, it marks as not visible all the points that lie outside a finite area of it.

To understand a plot of a complex structure of points could be useful to see the points of different colors and a line structure that links the points. To do this the simulator provides a script to plot the projected points that requires in input a connection matrix that points out which point and which connections we must plot and specify also their colors.

\subsection{Possible Improvements}

In the actual version just geometric points and their connections can be plotted.
A useful improvement could be to consider limited polygonal surface, starting with very simple one (triangles and rectangles). 
The problem in this case is that this surfaces can cover the points behind them and the simulator must mark them as not visible. The simplest algorithm is to project the objects in decreasing order of distance, and mark as not visible old objects as they are covered by the newer. This procedure (Painter's algorithm) is going to fail when there are cyclically overlapping polygons so that is impossible to determine which is the nearest one.

In real applications the most used algorithm is \emph{Z-Buffering} (algorithm~\ref{zdepth}; used for example in OpenGL and DirectX) 
\begin{algorithm}
  \caption{Z-depth algorithm}
  \label{zdepth}
  \begin{algorithmic}[1]
    \For {polygon $P$ in the polygon list}
      \For {pixel$(x,y)$ that intersects $P$}
	\State Calculate z-depth of $P$ at $(x,y)$
	\If {z-depth $<$ z-buffer$[x,y]$}
	   \State z-buffer$(x,y)$=z-depth
	   \State mark previous $(x,y)$ point as not visible
	\EndIf
      \EndFor
    \EndFor\\
    \Return 
  \end{algorithmic}
\end{algorithm}

The size of this algorithm depends mostly on the resolution of the image because it must calculate the depth of every single pixel and is typically implemented in hardware. There is also a distributed approach called \emph{W-Buffering}.
In our case, since we are interested only on the visibility of geometric points (which represent the markers) we could apply this algorithm to the subset of pixel of the camera where a point has been projected. This result in a considerably reduction of complexity, with the drawback that we aren't able to represent correctly the polygon, which can therefor act only as ``invisible curtains'' that hide the points (figure~\ref{fig:zbuf}).

\begin{figure}[tbp]
  \centering
  \includegraphics[width=0.45\textwidth]{graphs/schema_zbuf.pdf}
  \caption{Hidden areas and points. The polygons are not visible on the image plane, but they are used to create some surfaces that hide the points whose projection rays cross their inside.}
  \label{fig:zbuf}
\end{figure}

% \begin{figure}[tbp]
%   \centering
%   \includegraphics[width=0.45\textwidth]{radial.png}
%   \caption{Area hidden by polygon}
% \end{figure}
% This can be done assigning to every polygon a surface normal. The sense of this vector is chosen following a conventional rule based on the order of vertices: for example, the sense is pointing towards the observator if this one see the vertices in clockwise order.


\newpage
\textcolor{white}{\_}\\
\newpage
\section{MCMCDA}\label{Alberto}
Follow the trajectory of a single particle in clear space is an extreme simple problem. There is no error variance, and no other track to consider. On the opposite side, tracking and labelling multiple objects is generally a NP-Hard problem (as for \emph{JPDA}).\\
There are some algorithms that solve this problem with a simplified approximated implementation. This reduces the problem to a polynomial computational order. Among various methods, at the state of the art the best one is known as \emph{Markov Chain Monte Carlo Data Association}  \cite{MCMCB}.
% spiegazione markov chain
\subsection{Model Formulation}
We suppose that there is an hidden Markov chain $\mathcal{M}$ \emph{suitable for modelling track labelling}. Observation time is a sliding window of $T\in \mathbb{Z}^{+}$ instants, in a region $\mathcal{R}$ with volume $V$. If the current instant is called $H$, the sliding window is the set $\left\{H-T,\dots,H\right\}$.  Let $K$ be the unknown number of existing or existed target at the present instant $H$. At each time, a target disappears with probability $p_{z}$ and remains alive with probability $1-p_{z}$ (a Bernoulli distribution). Otherwise, the number of new arising object at each time is patterned after a Poisson distribution with parameter $\lambda_{b}V$, where $\lambda_{b}$ is the birth rate.\\
The (noisy) Markov chain of the trajectory of each object is, including observing the measurements:
\textcolor{white}{\_}
\begin{equation}
\left\{ \begin{array}{l}
x_{t+1}^{q}=F^{k}\left(x_{t}^{q}\right)+w_{t}^{q}\\

y_{t}^{j} = \left\{ \begin{array}{l}
H^{j}\left(x_{t}^{k}\right)+v_{t}^{j}\\
u_{t}\\
\end{array} \right.\\

\end{array} \right.
\end{equation}
where the $j^{\textrm{th}}$ observation is from target $x^{k}_{t}$ and $u_{t}\sim \textrm{Unif}\left(\mathcal{R}\right)$ stands for a random process for false alarms.
The multi-target tracking and labelling problem is to reconstruct the exact target associations and their trajectories, from measurements. It can be modelled as
\textcolor{white}{\_}
\begin{equation}
\left[ \begin{array}{l}
\tau_{k}=\tau_{\left[H-T,\dots,H\right]_{k}}=\textrm{Alg}\Big(y_{\left[H-T,\dots,H\right]}\Big)\\ \\
\tau_{k_{a}}\cap\tau_{k_{b}}=\emptyset \qquad \forall a\neq b \in \{0,\dots,K\}\\ \\
\omega=\left\{\tau_{0},\tau_{1},\dots,\tau_{K}\right\}\\

\end{array} \right.
\end{equation}
where $\omega$ is the proposed partition of trajectories, taken from $\boldsymbol{\Omega}$, the space of all possible collections of partitions of any target paths. Obviously, $\textrm{Alg}(Y)$ symbolizes the selected algorithm in order to compute the most likely tracks from the measurements.\\
$\tau_{0}$ is a special array containing all measurements labelled as \emph{false alarms}, invalid data.\\
The set of all measurements at a certain frame $t$ is
\textcolor{white}{\_}
\begin{equation}
Y_{t}=\left\{y_{t}^{j}\quad | \quad j=1,\dots,n_{t}\right\}
\end{equation}
where $n_{t}$ is the number of observations taken at that instant.
The set of all measurements up to the current frame is
\textcolor{white}{\_}
\begin{equation}
Y_{t}=\left\{Y_{t} \quad | \quad t=1,\dots,H\right\}
\end{equation}
even though implementations cited in this paper use only a subset, taken from a window of frames:
\textcolor{white}{\_}
\begin{equation}
Y_{W}=\left\{Y_{t} \quad | \quad t=H-T,\dots,H\right\}
\end{equation}
Finally, given the observations, for a certain distribution $\omega$ there is an \emph{a posteriori} probability $P[\omega | Y ]$.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Online Multi-Scan Markov Chain Monte Carlo Data Association algorithm}
\emph{Markov Chain Monte Carlo} (aka \emph{MCMC}), as in its Metropolis–Hastings implementation, is a method to generate repeatedly indirect sequences of random samples from an unknown probability distribution, in order to approximate it. Obviously, as assumption is that there are no deterministic ways for extract the exact results.\\
The standing idea under the \emph{MCMCDA} algorithm is that, at a given time $t$, previous detected tracks $\tau_{i}^{t-1}$ and measurements up to $Y_{j}^{t}$ are part of an hidden (asymptotic) Markov Chain, the cited unknown distribution for the method. States of this chain are the true targets $x_{k}^{t}$ (as well as their trajectory), but all we know are the (noisy) observations $Y_{j}^{t}$. The problem is to associate each measurement to the right target, and to recognize the \emph{false alarms}.\\
% 			Cenni SINGLEscan
MCMCDA is based upon JPDA. In a best case, number of target is known and fixed. Therefore, there is no need to look after all the history of previous detected tracks $\tau_{i}^{t-1}$ and past measurements. The \emph{single-scan} version si made for that, it is reduced to a simple filtering scheme.\\
% 			Spiegazione MULTIscan
Generally, number of targets is unknown and may change over time. So, it is possible that previous associations may not be correct in light of new evidence: e.g. a track $\tau_{i}$ might be an union of two tracks, otherwise unused measurements $Y_{j}^{t-h}$ can produce a trajectory for a new target (appeared at same point in the past).\\
% 			ONLINE
Unfortunately, the computational complexity of the multi-scan MCMCDA algorithm increases with the time and the number of measurements considered  \cite{MCMCDA}. Otherwise, old associations and observations are less influential then newer  \cite{MCMCDA}. Therefore, it is logical to assume a sliding window of time $[t_{\textrm{curr}}-t_{\textrm{win}}+1, \dots, t_{\textrm{curr}}]$, where $t_{\textrm{curr}}$ is the present instant and $t_{\textrm{win}}$ the size of the window. Only subsets of $Y$ and $\omega$ are used, therefore generally this version of the algorithm returns a suboptimal solution, but in a reasonable time. This implementation is called \emph{online multi-scan MCMCDA}.\\
% 			Descrizione base dell'implementazione dell'algoritmo
The essence of the algorithm is proposing a random move (chosen following the \emph{proposal distribution} $\omega'$) for a track from past associations and up to current measurements $\omega$. Possible moves are: birth, death, split, merge, extension, reduction, update and switch. If the move is feasible, it is randomly considered as acceptable. Given a number $u\sim \textrm{Unif}\left[0,1\right]$, the acceptance pseudo-probability $\textrm{\emph{Acceptance}}(\omega,\omega')$ must be greater than $u$. Typically, $\textrm{\emph{Acceptance}}(\omega,\omega')$ may assume values next to $0$ or is equal to $1$. Finally, the new track distribution is accepted if its posterior conditional probability is greater then the initial's one.\\
This procedure is iterated $N_{mc}$ times, that stands for the number of iterations in order to reach the asymptotic distribution of the hidden Markov chain.
The pseudocode of the algorithm is illustrated in code block \ref{algmulti}.

\begin{algorithm} 
\caption{MCMCDA algorithm} \label{algmulti}

\begin{algorithmic}[1]
\INPUT $Y$, noisy 3D points (up to current frame)
\INPUT $\omega$, labelling \& tracking history of all $\tau_{i}$ (up to previous frame)
\State $\hat{\omega}:=\omega$
\State
 	
\For{$n=1,\dots ,N_{mc}$} 
	\State $\omega' = \textrm{\emph{Proposal Distribution}}(\hat{\omega},Y)$   %%%%%%%%%%%%%%%%%%%%%% CHANGED
	\State $u\sim \textrm{Unif}\left[0,1\right]$
	\If{$u<\textrm{\emph{Acceptance}}(\omega,\omega')$}
		\State $\omega:=\omega'$
	\EndIf

	\State
	\If{$P[\omega | Y ] > P[\hat{\omega} | Y ]$}
		\State $\hat{\omega}:=\omega$
	\EndIf	
\EndFor

\State
\OUTPUT $\hat{\omega}$ 

\end{algorithmic}
\end{algorithm}


% 			Spiegazioni implementazione della proposal

\subsection*{Proposal Distribution}
In order to constructing a Markov Chain with the desired properties, a random move $m$ is chosen for the $\textrm{\emph{Proposal Distribution}}(\hat{\omega},Y)$, according to a probability distribution $\xi_{K,H}(m)$. If no other condition is given, $\xi_{K,H}(m)$ is modeled as (discrete) uniformly distributed. At the first frame and iteration, when no $\tau_{k}$ is present, the only possible move is the birth of a track. After the first iteration ($H=1$) there are only two eventualities: either a track can be deleted either another can originate.\\
When there is only a single target, no merging or switching between tracks is possible.\\
If the chosen move is feasible, it is considered as acceptable for a new partition $\omega'$, otherwise the initial $\hat{\omega}$ remains unchanged. Possible moves are: birth, death, split, merge, extension, reduction, track update and track switch. The moves are graphically illustrated in figure \ref{fig:proposaldistro}\\

 In \cite{MCMCDA} and \cite{MCMCB} the $\textrm{\emph{Proposal Distribution}}(\cdot,Y)$ has as input $\omega$, not $\hat{\omega}$. This is an error: the proposed $\omega'$ would be computed upon a move that might pass the acceptance check, but maybe based on trajectories less probable then $\hat{\omega}$ (as $\hat{\omega}$ is the most probable partition). Every iteration, a wronger partition might be extracted and approved for the next cycle. Imposing $\hat{\omega}$ as input to the function, only better partitions can be extracted and approved. \\
 As showed in simulation results section, this fix improves the efficiency of the algorithm from $\approx15\%$ of target correctly detected and tracked up to $\approx99.9\%$ of them.

\subsubsection*{Feasible Moves}
\textbf{1) Birth Move}. We choose randomly an instant in the sliding window, $t_{1}\sim \textrm{Unif}\left\{1,\dots,T-1\right\}$. Then we choose at random a second instant slightly afterwards, $H-t_{1}+d_{1}$, choosing $d_{1}$ from discrete distribution $d_{1}\sim \textrm{Unif}\left\{1,\dots,\bar{d}\right\}$.
A birth move is feasible if there are ``free'' (not already associated to a certain $\tau_{k}$) data at $H-t_{1}$ that are, in turn, linkable to other free data at $H-t_{1}+d_{1}$. A link is possible if a path between data at different times can exist, given a maximum velocity-per-frame $\bar{v}$\footnote{In \cite{MCMCbay} the set of all linkable data is modelled as mathematical structure, called \emph{the neighborhood tree of measurements}.}.\\
If the move is feasible, a new track can originate, and the number of all tracks increase: $K_{1}:=K+1$. Selecting randomly from linkable data at each time $H-t_{1}$, $H-t_{1}+d_{1}$, $H-t_{1}+d_{1}+d_{2}$, etc, the track is recursively extended. Adding measurements is a process that may terminate at a certain frame, with probability $p_{z}$. The proposed partition is $\omega'=\omega\cup \{\tau_{K_{1}}\}$\footnote{As implicit assumption, $\omega'=\omega\cup \{\tau_{K_{1}}\}$ includes the proposing of $\tau_{0}'=\tau_{0} \backslash \{\tau_{K_{1}}\}$. Equally for other moves.}\\
\textbf{2) Death Move}. Selected a not-already-deleted track $k_{d}$ from $\textrm{Unif}\left\{1,\dots,K\right\}$, it will be erased from $\omega$. The partition will be reassigned as $\omega'=\omega\backslash \{\tau_{k_{d}}\}$.\\
\textbf{3) Split Move}. We select an alive track $k_{s}$ from $\textrm{Unif}\left\{1,\dots,K\right\}$ with at least 4 instants contained, so $|\tau_{k_{s}}|\geq 4 $. Then, we randomly choose a time $t_{r}\sim \textrm{Unif}\left\{t_{2},\dots,t_{|\tau_{k_{s}}|-2}\right\}$. A new track originates as $\tau_{k_{s_{2}}}=\left\{\tau_{k_{s_{2}}}(t_{r+1}),\dots,\tau_{k_{s_{2}}}(t_{|\tau_{k_{s}}|-2})\right\}$, and the old one becomes $\tau_{k_{s_{1}}}=\left\{\tau_{k_{s_{1}}}(t_{1}),\dots,\tau_{k_{s_{1}}}(t_{r})\right\}$. The global number of tracks increases, $K_{1}=K+1$, and the $K_{1}^{\textrm{th}}$ is $\tau_{k_{s_{2}}}$ obviously. The proposed partition will be $\omega'=\big( \omega \backslash \{\tau_{k_{s}}\} \big)  \cup \{\tau_{k_{s_{1}}}\} \cup \{\tau_{k_{s_{2}}}\}$.  \\
\textbf{4) Merge Move}. For a merge move we must select (uniformly at random) a couple of 2 tracks $(\tau_{k_{s_{1}}},\tau_{k_{s_{2}}})$, the first must end before the latter and they must be \emph{linkable}, as defined above. Then, they will be combined in a single track $\tau_{k_{s}}=\tau_{k_{s_{1}}}\cup \tau_{k_{s_{2}}} $. The proposed partition will be $\omega'=\big( \omega \backslash ( \{\tau_{k_{s_{1}}}\}, \{\tau_{k_{s_{2}}}\}  ) \big)\cup \{\tau_{k_{s}}\}$. \\
\textbf{5) Extension Move}. Selected an alive track $k_{e}$ from $\textrm{Unif}\left\{1,\dots,K\right\}$, we assign free linkable measurements after the last instant of the track $t_{|\tau_{k_{e}}|}$ as in the birth move. The proposed partition will be $\omega'= \omega \cup \big( \{\tau_{k_{s}}\} \cup \{\tau_{k_{s}}(t_{|\tau_{k_{e}}|+1}),\dots,\tau_{k_{s}}(t_{|\tau_{k_{e}}|'})\}\big)$. \\
\textbf{6) Reduction Move}. Select an alive track $k_{rd}$ from $\textrm{Unif}\left\{1,\dots,K\right\}$ and an instant $t_{r}\sim \textrm{Unif}\left\{t_{2},\dots,t_{|\tau_{k_{rd}}|-1}\right\}$. The track will be shortened to $\tau_{k_{rd'}}=\left\{\tau_{k_{rd}}(t_{1}),\dots,\tau_{k_{rd}}(t_{r})\right\}$. Measurements after $t_{r}$ are freed. The proposed partition becomes $\omega'=\big( \omega \backslash \{\tau_{k_{rd}}\} \big)  \cup \{\tau_{k_{rd'}}\}$.\\
\textbf{7) Track Update Move}. Selected an alive track $k_{u}$ from $\textrm{Unif}\left\{1,\dots,K\right\}$, it will be shortened as in the reduction move, from the beginning up to the chosen instant $t_{r}$. Then, we reassign free linkable measurements after that time. The proposed partition becomes $\omega'=\big( \omega \backslash \{\tau_{k_{u}}\} \big)  \cup \{\tau_{k_{u'}}\}$.\\
\textbf{8) Track Switch Move}. For switching tracks, at first we must (uniformly at random) select a pair of already associated measurements of two distinct tracks in different instants. We name this pair $\big(\tau_{k_{s_{1}}}(t_{p}),\tau_{k_{s_{2}}}(t_{q})\big)$. The pair must satisfy the \emph{linking bond}: $\tau_{k_{s_{1}}}(t_{p})$ should be able to reach $\tau_{k_{s_{2}}}(t_{q})$ and its following associations in a time proportional to $\bar{v}$ and $\tau_{k_{s_{2}}}(t_{q})$ should be able to reach the following associations to $\tau_{k_{s_{1}}}(t_{p})$. This is the essence of the switching. Chosen the pair, associations after $t_{p}$ and $t_{q}$ are (respectively) exchanged. Switched tracks will become
\textcolor{white}{\_}
\begin{equation*}
\begin{array}{c}
\tau_{k_{1'}}=\left\{\tau_{k_{1}}(t_{1}),\dots,\tau_{k_{1}}(t_{p}),\tau_{k_{2}}(t_{q+1}),\dots,\tau_{k_{2}}(t_{|\tau_{k_{2}}|})\right\}\\ \\
\tau_{k_{2'}}=\left\{\tau_{k_{2}}(t_{1}),\dots,\tau_{k_{2}}(t_{q}),\tau_{k_{1}}(t_{p+1}),\dots,\tau_{k_{1}}(t_{|\tau_{k_{1}}|})\right\}\\
\end{array} 
\end{equation*}
Therefore, the proposed partition will become $\omega'=\big( \omega \backslash ( \{\tau_{k_{1}}\}, \{\tau_{k_{2}}\}  ) \big)\cup \{\tau_{k_{1'}}\}\cup \{\tau_{k_{2'}}\}$.\\

% figura 
\begin{landscape}
\begin{figure*}[p]
\begin{center}
    \includegraphics[width=17.5cm,keepaspectratio,trim=1cm 8.3cm 1cm 1cm,angle=-90]{graphs/proposal.pdf}% trim = left bottom right top
    \caption{Graphical illustration of feasible moves for a proposal partition.}
        \label{fig:proposaldistro}
    \end{center}
\end{figure*}
 \end{landscape} 
 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 			Spiegazioni implementazione della accettazione
\subsection*{Acceptance Probability}
We now describe the implementation of the acceptance probability mentioned in the previous sections.
The acceptance probability is:
\begin{equation}
A(\omega,\omega')=min\left(1, \frac{\pi(\omega')q(\omega',\omega)}{\pi(\omega)q(\omega,\omega')}\right)
\label{eq:acceptance}
\end{equation}
That is, if we have a feasible solution $\omega$ and we propose a new solution $\omega'$ (chosen at random following the proposal distribution $q(\omega,\omega')$), we accept the move and switch to new solution with probability $A(\omega,\omega')$.
In our case the distribution $\pi(\omega)$ is represented by the posterior $P(\omega|Y_{1:t})$ This distribution take into account all the aspect of a feasible solution $\omega$:

\begin{multline*}
P(\omega|T_{1:T}) = \frac{1}{Z_0}\prod_{\tau\epsilon\omega\backslash\{\tau_{0}\}}\prod_{i=2}^{|\tau|}\mathcal{N}(\tau(t_{i});\hat{y_{t_{i}}}(\tau),B_{t_{i}})\times\\
 \times \prod_{t=1}^{T}p_{z}^{z_t}(1-p_z)^{m_{t-1}z_t}\cdot P_d^{d_t}(1-P_d)^{u_t}\lambda_b^{a_t}\lambda_f^{f_t}
\end{multline*}

where $m_t$ denote the number of target at time $t$, $a_t$ the new targets, $z_t$ the terminated ones, $d_t$ the number of detections, $u_t$ the number of undetected targets, $f_t$ the false alarms and $n_t$ the total number of measurements (including false alarms).
The value $d_t$ is simply the length of the measurements vector; to determine $z_t$ we have to count how many tracks that at time $t-1$ where associated to a measurements don't have any associations at time $t$; to determine $a_t$ we must count how many new tracks have been created (i.e how many measurements are associated with tracks that have been created at time $t$). Therefore we have:
\begin{eqnarray*}
m_{t} & = & m_{t-1} + a_{t} - z_{t}\\
u_{t} & = & m_{t}-d_{t}\\
f_{t} & = & n_{t}-d_{t}\\
\end{eqnarray*}
Moreover, we have to evaluate the Gaussian $\mathcal{N}(\hat{y_{t_{i}}}(\tau),B_{t_{i}})$ at point $\tau(t_{i})$ for each track $\tau$ and for each instant $t_i$.
This means that, fixed a feasible solution $\omega$, we have to create a Kalman filter for each track starting from the instant where the track begins until the current instant and evaluate the Gaussian at each step. This is computationally the heavier section of the whole algorithm because it must run the Kalman filter on the entire story for each track of each proposed feasible solution $\omega$.
Since in equation~\ref{eq:acceptance} appears only the ratio $\frac{\pi(\omega')}{\pi(\omega)}$ we don't have to compute the normalization constant $Z_0$. 

To calculate the value of the proposal distribution, $q(\cdot,\cdot)$, we must take into account the distribution $\xi_{K,H}(\cdot)$ evaluating it at the point $m$ which represents the kind of move we have proposed.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection*{Possible Improvements}
In order to optimizing the MCMCDA algorithm for applying it to marker tracking for motion capture, some enhancements can be used.

\subsubsection*{Shaping $\xi_{K,H}(m)$}
Motion capture hardware sets are typically optimized to minimize the probability of track's disappearing, $p_{z}$.  For most scenes recorded, number of actors (and therefore markers) normally remains constant for the entire length of the scene (supposing off-screen view as cut from scene shot). So, birth moves after the first frame are not very probable, like death ones. Besides, an extension move is much more likely than the others. To generate a balanced Markov chain, moves such as reduction, track update and track switch have to be almost probable as extension ones. Split and merge ones will be considered less probable.\\
For example, the distribution for choosing moves after the first frame can be reshaped as
\textcolor{white}{\_}\\
\begin{equation}
\xi_{K,H}(m) =\left\{ \begin{array}{lcl}
0.033 & m=1 & P\left[\textrm{birth}\right]=1/30\\
0.066 & m=2 & P\left[\textrm{death}\right]=1/30\\   
0.133 & m=3 & P\left[\textrm{split}\right]=2/30\\
0.2 & m=4 & P\left[\textrm{merge}\right]=2/30\\
0.5 & m=5 & P\left[\textrm{ext.}\right]=9/30\\
 0.666 & m=6 & P\left[\textrm{red.}\right]=5/30\\
 0.833  & m=7 & P\left[\textrm{t. switch}\right]=5/30\\
 1 & m=8 & P\left[\textrm{t. update}\right]=5/30\\
\end{array} \right.
\end{equation}
At the first frame $\xi_{K,1}(m)$ remains as supposed in a general situation.


\subsubsection*{Two-by-Two Track Distance Bond Pruning}
\begin{figure*}[h]
\begin{center}
    \includegraphics[width=17.5cm,keepaspectratio,trim=2.4cm 5cm 3cm 7cm]{graphs/distance_bond_check.pdf}% trim = left bottom right top
    \caption{Example of two-by-two track distance bond pruning. Tracks $\tau_{a}$ and $\tau_{b}$ respect the bond. Tracks $\tau_{a}$ and $\tau_{c}$ are never at the same distance, and $\tau_{b}$ and $\tau_{c}$ don't maintain a constant distance for all frames in the sliding window. Besides, $\tau_{c}$ has a missing measurement.  Only $\tau_{c}$ is not labeled as ``sure'' in the MCMCDA algorithm.}
        \label{fig:chdist}
    \end{center}
\end{figure*}

\begin{figure*}[!h] % figure is 2 pages afterwards
\vspace{-25pt}
\centerline{\subfloat[$\textrm{\emph{Proposal Distribution}}(\omega,Y)$]{\includegraphics[width=0.55\textwidth,keepaspectratio]{graphs/MCMCDAw_w.pdf}%
\label{30case} }
\hfil
\subfloat[$\textrm{\emph{Proposal Distribution}}(\hat{\omega},Y)$]{\includegraphics[width=0.55\textwidth,keepaspectratio]{graphs/MCMCDAw_hat.pdf}%
\label{80case}}}
\caption{Simulation results of a rotating and translating cube (depicted in figure \ref{mcmcdacube}), computed with the original algorithm (a) and with the fixed version (b), with $N_{mc}=30$, no false alarms and no missing data. Brown tracks are labelled as $\tau_{k}$ with $k>8$, others are referred to points from $1$ to $8$ of the cube.}
\label{fig_prop}
\end{figure*}

In motion capture applications, markers are fixed to actor's body, that is a structure with joints. Hence, there are available shape bonds and many measurement are related with each others. An efficient way to capitalize on them is introducing a preliminary check on relative distances. At each frame we compute the distance of each different pair of measurements. This \emph{pre-check} can be easily parallelized. We suppose to use an wide sliding window. Suppose that, when reached the steady state (and therefore a full sliding window is available) we have extracted some tracks with MCMCDA algorithm. If a pair of tracks maintain the same distance from each other during each frame of the sliding window (even rotating), they are considered both as correctly detected target trajectories. The check is made two by two for all tracks. If the distance bond is confirmed for 3 tracks, two by two, their prediction have even less error. With four tracks bonded, a point in 3D space is settled. Only for the current frame $H$, they will be pruned from Metropolis–Hastings method iterations, labeling them as ``sure''. They can be extended, but not deleted, split, reduced, switched or updated. Tracks that never pass the test  can be considered as \emph{free trajectories}, \emph{trajectories with missing measurements} or \emph{false alarms}.\\
This shape bond speeds up the achievement of the stationary distribution of the Markov chain, as some tracks (states of the chain) have already reached their final state.\\
A random sequence $U_{n}$ converges almost surely if  $P\big[\lim_{n\to\infty}U_{n} = U\big]=1$. Supposing the distance between two tracks at time t in the sliding window as a random variable $d_{t}=\textrm{distance}(\tau_{1},\tau_{2})$, a sequence $D_{t}$ can be modeled as
\begin{equation}
\begin{array}{c}\displaystyle
D_{t}=P\big[d_{t}=d \quad \forall t \in [H-T+1,T]\big]=\\ \displaystyle
=\prod_{t=2}^{T}P\big[d_{t}=d | d_{t-1}=d\big]
\end{array}
\end{equation}
but, knowing that $d_{t}=d$ at each frame on a large enough sliding window $[H-T+1,T]$, we can approximate $P\big[\lim_{n\to\infty}D_{n}=1\big]=P\big[\lim_{n\to T}D_{n}=1 \big]=1$
that converges with probability $1$. Convergence almost surely implies convergence in probability and in distribution, so the MCMCDA algorithm might exclude these tracks from computation.


\subsubsection*{Volume or Distance Clusters Splitting up}
A further refinement to the MCMCDA method might be dividing measurements into groups, assigning to each one a CPU.\\
As total volume of tracking is known, the original algorithm can be parallelized. Splitting up in smaller portion of space, we can consider as input only a restriction of $\omega$ and $Y$, such as $\omega_{v}$ and $Y_{v}$. An instance of MCMCDA runs for each portion of space.\\
This approach has some limitations. If a region is far from the scene, its parallelized CPU might be unused. Otherwise, some regions have an high density of markers, so their associated CPU are put under much more stress than the others.\\
A more efficient approach is dividing measurements (and consequently their related track $\tau_{k}$) into clusters, based on relative distances between 3D data. Set a maximum distance $D_{\textrm{max}}$, a cluster rises as set of all near measurements in a certain region of space. An instance of MCMCDA runs for each cluster. As said above, computing distances can be easily parallelized. This enhancement is also scalable. Furthermore, any cluster brings to a strongly connected graph, or can be modeled as a small world network with greatest achievable local clustering coefficient (one). This can be capitalized by distributed algorithm.\\
This approach has some limitations as well. A cluster may be composed by not many measurements, otherwise some other might have a high number of data (such as the ``cluster of an hand''). Clusters might also not remain the same during time.  \\
Finally, regrouping tracks (finding the right continuation) of divided measurements remains a common problem for splitting up methods.

\newpage
\subsection*{Simulations and results}
As shown in \ref{fig_prop}, only $\hat{\omega}$ (as input) is able to get the algorithm to the stationary distribution of the Markov chain, even with few iterations on perfect data.\\
In Figure \ref{fig_sim} some simulations computed on the fixed version and improved distribution $\xi_{K,H}(m) $ are illustrated, with optimized \textsc{matlab}\textcopyright  code run on a pc with 4 Gb of RAM and Intel Core 2 Duo processor.. As $N_{mc}$ is increasing in order, efficiency of the method can be increased. Also time elapsed raises. If disturbance raises, more iterations have to be done. Missing data don't influence on trajectories only at higher number of cycles-per-frame.

\begin{figure*}[!h]
\vspace{-25pt}
\centerline{\subfloat[Original rotating and translating cube]{\includegraphics[width=0.55\textwidth,keepaspectratio,trim=2cm 9cm 3cm 14cm]{graphs/MCMCDAcube.pdf}%
\label{mcmcdacube} }
\hfil
\subfloat[$N_{mc}=30$]{\includegraphics[width=0.55\textwidth,keepaspectratio]{graphs/MCMCDA10f30iter.pdf}%
\label{30case}}}
\label{fig_sim}
\end{figure*}
\begin{figure*}[h]
\ContinuedFloat
\vspace{-25pt}
\centerline{\subfloat[$N_{mc}=80$]{\includegraphics[width=0.55\textwidth,keepaspectratio]{graphs/MCMCDA10f80iter.pdf}%
\label{80case} }
\hfil
\subfloat[$N_{mc}=1000$]{\includegraphics[width=0.55\textwidth,keepaspectratio]{graphs/MCMCDA10f1000iter.pdf}%
\label{1000case}}}
\caption{Original moving object and tracking simulation results computed by MCMCDA algorithm, with $10\%$ false alarms and $10\%$ missing data in $100$ frames, increasing $N_{mc}$. Other parameters: $\bar{v}=0.15$ cm/frame, $\bar{d}=2$.  Brown tracks represent trajectories labeled above $8^{\textrm{th}}$ target, therefore incorrectly associated. }
\label{fig_sim}
\vspace{-20pt}
\end{figure*}

\begin{figure*}[!h]
\centerline{\begin{tabular}{|c|c|c|} \hline
\multicolumn{3}{|c|}{Simulation results} \\ \hline
\multirow{3}{*}{$N_{mc}=30$} & wrong associations & $\sim36\%$ \\
& generated tracks & $12$ \\
& missing labelled marker & $\sim10\%$ \\
& runtime & 40.70 minutes \\ \hline
\multirow{3}{*}{$N_{mc}=80$} & wrong associations & $\sim12\%$ \\
& generated tracks & $9$ \\
& missing labelled marker & $\sim6\%$ \\
& runtime & 102.30 minutes \\ \hline
\multirow{3}{*}{$N_{mc}=1000$} & wrong associations & $\sim0.1\%$ \\
& generated tracks & $8$ \\
& missing labelled marker & $\sim0.1\%$ \\ 
& runtime & 364.75 minutes \\ \hline
\end{tabular}}
\end{figure*}


\newpage
\textcolor{white}{\_}\\
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\newpage
\textcolor{white}{\_}\\
\newpage
\section{Stable Marriage Problem}\label{Carlo}

In a data labeling problem we have to associate measurements with entities which have generated them. We thought that the Stable Marriage Problem could became helpful because it provides a way to bind the elements of two different sets in the best way, accordingly to some index that has to be defined.

\subsection{SMP Definition}

Let $M=\{A,B,C,...\}$ be a set of males and $F=\{a,b,c,...\}$ a set of females with the same cardinality of $M$. Each male have a ranked list of preference of each female, in descending order. Each female has a similarly ranked list, so that in every female list all man are present and viceversa.

A marriage is a match between one element of $M$ and one of $F$ (e.g. $A-b$).
A pair $X-y$ is a dissatisfied one if in the solution $S$ exists two marriages, $X-z$ and $W-y$ such that $X$ prefer $y$ more than is current partner $z$ and woman $y$ prefers $X$ more than her current partner $W$.
A set of marriages $M$ is called stable if there are no dissatisfied pair.

\begin{figure}[tbp]
  \centering
  \includegraphics[width=0.45\textwidth]{graphs/marriage.pdf}
  \caption{Example of two SMP solutions. The first is a stable set of marriages, the second has some dissatisfied pairs: $A$ prefers $a$ more then is current mate and so do $a$, hence they should be paired in a correct SMP solution. }
  \label{fig:smp}
\end{figure}


\subsection{SMP Algorithm}

A simple and deterministic algorithm to determine a stable marriage is shown in algorithm~\ref{smp_proposal}
\begin{algorithm}
  \caption{Proposal algorithm}
  \label{smp_proposal}
  \begin{algorithmic}[1]
%     \REQUIRE $X = \{x_1, \dots, x_n\}$, $k$
%     \ENSURE  $M = \{\mu_1, \dots, \mu_k\}$, $A$
%     \WHILE{$M \leftarrow M_0$ (e.g. scelta casuale in $X$)}
%     \{$x_j \in X$}
%     \STATE $a_{ij} = {\arg\min}_{i \in \{1,\dots,n\}}\|x_j - \mu_i\|$
%     \ENDFOR
    \While{$\exists$ unpaired man $X$}
    \State pick first woman $w$ on $X$'s list
    \State remove $w$ from $X$'s list
    \If{$w$ is married}
	\If{$w$ prefer $X$ to her mate $Y$}
	  \State set $X-w$ as married
	  \State set $Y$ as unpaired
	\Else
	  \State $X$ remains unpaired
	\EndIf
    \Else
	\State set $X-w$ as married
    \EndIf
    \EndWhile\\
    \Return
  \end{algorithmic}
\end{algorithm}

Each man make a marriage proposal to the woman he preference, and then remove the woman from his list (so that he is not going to make a new proposal to the same woman). Each woman decide to accept the proposal if she is unengaged or if she prefer the new man rather then his actual mate. The algorithm terminates when all males are engaged.

\subsection{Algorithm Analysis}

We have to prove that the algorithm terminates with all men (and women) paired and that the solution is a stable marriage.
Every time that a man make a proposal to a woman, this one is removed from his list. A man whose proposal has been refused will continue making proposals to other women until the end of his list, so that every woman receives at least one proposal. The woman will always accept her first proposal so no woman can be left unmarried.
In the worst case the proposal-algorithm makes $O(n^2)$ proposal.
To show that it is a stable marriage, let’s assume we have a dissatisfied pair, $X-b$, where in the marriage
they are paired as $X-a$ and $Y-b$. Since $X$ prefers woman $b$ over his current partner $a$, then he must have
proposed to $b$ before $a$. Woman $b$ either rejected him or accepted him, but dropped him for another better
man than $X$. Thus, $b$ must prefer $Y$ to $X$, contradicting our assumption that $b$ is dissatisfied, therefore the solution it is a
stable marriage set.
It can be shown that, although the worst case of the algorithm is $O(n^2)$, the expected (average) case is $O(n\lg n)$ and deviation is small from the expected value (for the complete proof see \cite{SMP}).

\subsection{Optimality}
We now define a criterion for optimality. A marriage between a man A and a woman B is \emph{feasible} if there exists a stable pairing in which A and B are married.
A pairing is \emph{male-optimal} if every man is married with his highest ranked feasible partner.
We can prove that the proposal algorithm~\ref{smp_proposal} is a male-optimal one.

Suppose that the solution were not male-optimal. Let be $M$ one man that has just been rejected by his optimal woman $W$. Since the women always accept their first proposal, then $W$ must have been paired width a men $Z$ that she prefers more then $M$. $M$ has been just rejected and the proposal are made following the preferences order, that means that $W$ is in the same or higher rank in the $Z$'s preferences list.
Moreover, since $W$ is the optimal partner for $M$, there exist a solution $S$ where $M$ is paired with $W$~\footnote{this solution would be achieved if $M$ had made his proposal before $Z$}. But in this solution $S$ we would have that $W$ prefers $Z$ to $M$ and $Z$ prefers $W$ to his mate in $S$. Thus solution $S$ is unstable and this is a contradiction because we have proven that the proposal algorithm provides stable marriages.
The solution could be female-optimal if we switch the roles of men and women so that are the women the ones to make the proposal to men.

\subsection{Application to the Data Labeling}

In the data labeling we have two sets (the observations and the targets) we want to put in correspondence.
If we make the assumption that:
\begin{itemize}
  \item the number of marker is constant
  \item there aren't false alarms
  \item all the marker are detected
\end{itemize}
then the number of the targets is equal to the number of observations (measurements) at every time. We can consider the targets set as the males set and the observations set as the females set. We need just to impose a rule to express the preferences both for the males and the females.
Let be $M={t_1, t_2, ... t_n}$ the set of targets, and $F={y_1, y_2, ... y_n}$ the set of new measurements a time T. We calculate the prediction (and variance) for time $T$ for each $t_k$ based on his past history ($[T_0,...,T]$) with a Kalman filter.
Then, fixed $t_k$, we evaluate the Gaussian distribution $\mathcal{N}(\hat{t}_k,\Sigma_{t_k})$ on every $y_i$: sorting these values decreasing we obtain the males (targets) preference list.
For the female preferences, we simply sort the males by descending order of distance from the female item (figure \ref{fig:crit}).

Since we have proven that the SMP algorithm is male-optimal, we rest assured that, for every single frame, the solution is optimal for the male criterion in the sense of optimality we have previously defined.

\begin{figure}[tbp]
  \centering
  \includegraphics[width=0.40\textwidth]{graphs/criteri.pdf}
  \caption{Creating preferences list: males (top) and females (bottom)}
  \label{fig:crit}
\end{figure}

\subsection{Simulation Result}

We ran some tests to investigate the behavior of this algorithm under different conditions.
The first test we ran was with eight points moving with variable speed and direction in a 3D space. The dynamical model of the system we applied to the Kalman filter was a \emph{Discrete White Noise Acceleration} (DWNA)~\footnote{We choose DWNA because it's a simple model that we think could describe well the marker's dynamic in a wide range of real application in the motion capture scenario. Clearly, to get the best results, the model of the Kalman filter must be adapted to the application.}. Our \textsc{matlab}\textcopyright implementation of SMP took ~1.1 seconds to complete the associations over 100 instants; The associations have been always correct.
We next run a test with 16 points arranged on the vertexes of two cubes. One of them was still, while the other pass through it while translating and rotating. Also in this case the associations were all correct although at some time the vertexes of the two cubes were lying upon each other.
Other more complicated tests show the goodness of this approach.
We tried with 8 points following a random walk. Note that in this case would have required to change the Kalman filter used to set the preferences list, since DWNA is not a good model for points that follow a random walk. The associations were all correct even when the points came close to each other. As a test designed specifically to put in trouble our algorithm, we set up 8 points that followed a ``particular'' random walk, being attracted towards axis origin. Even in this case the associations were very accurate: on the total frames, 82\% were correct, 11\% frames had one swap between a pair of marker~\footnote{Note that there can't be only one marker mistaken, since if it is wrong, it must have been taken the right place of another one} and the rest have more than two wrong associations. See table~\ref{tbl:percentage} for details.
$$
\left|\begin{array}{cc}
\#\text{wrong associations} & \text{percentage}\\
0\text{ (all correct)} & \sim82\%\\
2 & \sim11\%\\
3 & \sim2.5\%\\
4 & \sim2.5\%\\
5 & \sim0.5\%\\
6 & \sim0.5\%\\
7 & \sim0.5\%\\
8 & \sim0.5\%\end{array}\right|.
\label{tbl:percentage}
$$

\begin{figure}[tbp]
  \centering
  \includegraphics[width=0.45\textwidth]{graphs/SMP_simulazione1.pdf}
  \caption{Time evolution of 8 tracked markers lying on the vertexes of a cube tranlating and rotating}
  \label{fig:smp_simu1}
\end{figure}

\begin{figure}[tbp]
  \centering
  \includegraphics[width=0.45\textwidth]{graphs/SMP_simulazione2.pdf}
  \caption{Time evolution of 16 tracked markers lying on the vertexes of a cube tranlating and rotating}
  \label{fig:smp_simu2}
\end{figure}

\begin{figure}[tbp]
  \centering
  \includegraphics[width=0.45\textwidth]{graphs/SMP_simulazione3_RW.pdf}
  \caption{Time evolution of 8 tracked markers attracted towards origin while following a random walk}
  \label{fig:smp_simu2}
\end{figure}


% : in the case of two targets following trajectories lying near each other for some instants we found that they were frequently mistaken, similarly points which follow random walk path get no correct associations when they come near to each other. Note that in this case would have required to change the Kalman filter used to set the preferences list, since DWNA is not a good model for points that follow a random walk.

\subsection*{ Possible Improvements}

\subsubsection*{Variable Number of Targets, False Alarms and Undetection}

The assumptions made on number of targets, false alarms and detection of targets are very restrictive in a real application.
It's interesting to analyze solutions that allow to relax these assumptions.

If the number of marker can increase, we simply start some more Kalman filters when we find that the number of observation has increased. Less simple is the opposite case: if the number of marker can decrease we have to eliminate some tracks, that means to choose which of the previous time target is less likely to have a correspondent measurement at present time.
In this case we need to find a rule to delete the less probable measure. One way could be to to analyze the females (measures) preference list and keep only the first $K$~\footnote{with K number of targets (males)} measurements. A work-around that remits the choose of which track to eliminate to the SMP algorithm is to create $n=Nwoman-Nman$ factitious women such that they most certainly are the lower ranked ones for every man~\footnote{For example, we could place them very far away from every target; or we can just fill the males list adding factitious women on the end of the list until we reach the length of $n$ elements.}. Doing this, we can delete the men (tracks) which are associated by the SMP algorithm to the factitious women.

If we have false alarms or marker that can be undetected, this approach isn't applicable because we have a different number of males and females individual. There are some algorithm that extend the SMP concept to these cases; for example we could adopt the hospitals (or residents) problem. This is variant of SMP where the role of the female individuals is played by the hospitals, that can accept more then one proposal. Clearly, we must provide some other rules to determine which one of the individuals (measures) assigned to the same hospital (target) is the most probable or a way to merge the different measure into a single one.

\subsubsection*{Past History}

For each step our implementation takes into account the whole history for the filtering, but only the last measure to decide the associations between the two sets.
A way to improve performance could be to redefine the rules for the preferences to take into account all the measure of an interval.
% Consider, for example, a sliding window of two instants. We should:
% -Calculate the preferences ............. Discorso della difficoltà computazionele n!


\newpage
\section{Conclusions}

Our work lead to a \textsc{matlab}\textcopyright\ implementation of MCMCDA that result too slow for a real-time application. Even if implemented in a more efficient programming language (like C or C++ that are at least one order of magnitude faster) the processing takes too much time. On a quality level, our implementation lead some good results, recognizing correctly the targets in our simulations, that where limited to a few (max 80) instants due to the complexity of the algorithm. The sore point of algorithm is that it has to run a Kalman filter on the whole history (or on the whole sliding window) for each proposed solution and it has to propose a very high number of solutions to have a good approximation of the MCMCDA. 
To propose a solution applicable to real time application we have adapted the stable marriage problem. Our implementation of this algorithm has results in a fast and effective method for some cases, i.e. where we have constant number of targets and no false alarms or undetected targets. This method have wide leeways to improve quality of associations~\footnote{changing the way individuals choose their preferences} and, considering it's variants, can be extended to the case of a non constant number of target and with the presence of false alarms and undetected measurements.
\newpage
\appendices

\section{Markov Chain Monte Carlo method}

Markov Chain Monte Carlo method is a general method, its aim is extract a probability distribution as stationary distribution of a ad-hoc constructed Markov chain. Normally it is possible to construct a suitable chain for the purpose of the algorithm. This stationary distribution is considered as acquired after a great number of iterations, e.g. $50000$ for approximately $400$ targets \cite{MCMCDA}.

Its applications are, for example, longitudinal studies (like computing the right association between some targets and measurements), numerically calculating multi-dimensional integrals, radiocarbon dating, time trends for disease incidence and mortality  \cite{MCMC}.

In order to minimize the number of iterations, various random walk algorithms has been used, e.g. \emph{Metropolis–Hastings}, \emph{Gibbs sampling}, \emph{Slice sampling}, \emph{Multiple-try Metropolis}.

John Von Neumann and Stanislaw Ulam elaborated the original idea after 1944, conceived by Enrico Fermi in the 1930s.

\section{Metropolis–Hastings algorithm}
The Metropolis–Hastings algorithm is an algorithm (used for example by the MCMC method) for obtaining a sequence of random samples from a probability distribution for which direct sampling is difficult. It acts like a random walk at each step, using a proposal density and a method for rejecting proposed moves.
The Gibbs sampling algorithm is a special case of the Metropolis–Hastings algorithm which is usually faster and easier to use but is less generally applicable. The algorithm was named after Nicholas Constantine Metropolis. 


\newpage
\begin{thebibliography}{99}

\bibitem{LAP:Pier} Pierskalla W.P., \emph{The multidimensional assignment problem}. Operations Research 16, 1968, 422–431.

\bibitem{B-S} Bar-Shalom Y. and Fortmann T.E., \emph{Tracking and Data Association} Academic Press, San Diego, 1988.

\bibitem{MHT} Reid D.B., \emph{An algorithm for tracking multiple targets} IEEE Trans. on Automatic Control, vol. AC-24, no. 6, pp. 843–854, December 1979.

\bibitem{LagRelax} Deb S., Yeddanapudi M., Pattipati K. and Bar-Shalom Y., \emph{A generalized s-d assignment algorithm for multisensor-multitarget state estimation} IEEE Trans. on Aerospace and Electronic Systems, vol. 33, no. 2, pp. 523–538, April 1997

\bibitem{MCMCDA} Oh S., Russel S. and Sastry S., \emph{Markov chain monte carlo data association for general multiple-target tracking problems} IEEE Conf. on Decision and Control, 2004.

\bibitem{MCMC} Spiegelhalter D. J., Richardson S., W. R. Gilks, \emph{Markov chain Monte Carlo in practice} Chapman \& Hall 2004.
 
\bibitem{MCMCB} Oh S., \emph{Bayesian Formulation of Data Association and Markov Chain Monte Carlo Data Association} IEEE Conf. on Decision and Control, 2004.

\bibitem{SMP} Hunt W. \emph{The Stable Marriage Problem}

\bibitem{3DV} Ma Y., Soatto S., Koseckà J., Sastry S.S. \emph{An Initaion to 3-D Vision} Springer

\end{thebibliography}

\end{document}
