# Trajectory PMB Filters for Extended Object Tracking Using Belief Propagation

Yuxuan Xia, Ángel F. García-Fernández, Florian Meyer, Jason L. Williams, Karl
Granström, and Lennart Svensson Y. Xia and L. Svensson are with the Department
of Electrical Engineering, Chalmers University of Technology, Gothenburg,
Sweden. E-mail<EMAIL_ADDRESS>A. F. García-Fernández is with
the Department of Electrical Engineering and Electronics, University of
Liverpool, Liverpool, United Kingdom, and also with the ARIES research centre,
Universidad Antonio de Nebrija, Madrid, Spain. F. Meyer is with the Scripps
Institution of Oceanography and the Department of Electrical and Computer
Engineering, University of California San Diego, La Jolla, CA, USA. J. L.
Williams is with the Commonwealth Scientific and Industrial Research
Organization, Brisbane, Australia. K. Granström is with Zoox Inc., San
Francisco, CA, USA. The work of K. Granström was done when he was with
Chalmers University of Technology.

###### Abstract

In this paper, we propose a Poisson multi-Bernoulli (PMB) filter for extended
object tracking (EOT), which directly estimates the set of object
trajectories, using belief propagation (BP). The proposed filter propagates a
PMB density on the posterior of sets of trajectories through the filtering
recursions over time, where the PMB mixture (PMBM) posterior after the update
step is approximated as a PMB. The efficient PMB approximation relies on
several important theoretical contributions. First, we present a PMBM
conjugate prior on the posterior of sets of trajectories for a generalized
measurement model, in which each object generates an independent set of
measurements. The PMBM density is a conjugate prior in the sense that both the
prediction and the update steps preserve the PMBM form of the density. Second,
we present a factor graph representation of the joint posterior of the PMBM
set of trajectories and association variables for the Poisson spatial
measurement model. Importantly, leveraging the PMBM conjugacy and the factor
graph formulation enables an elegant treatment on undetected objects via a
Poisson point process and efficient inference on sets of trajectories using
BP, where the approximate marginal densities in the PMB approximation can be
obtained without enumeration of different data association hypotheses. To
achieve this, we present a particle-based implementation of the proposed
filter, where smoothed trajectory estimates, if desired, can be obtained via
single-object particle smoothing methods, and its performance for EOT with
ellipsoidal shapes is evaluated in a simulation study.

###### Index Terms:

Multi-object tracking, extended object tracking, random finite sets, sets of
trajectories, factor graph, particle belief propagation.

##  I Introduction

Multi-object tracking (MOT) refers to the process of estimating object
trajectories of interest based on sequences of noisy sensor measurements
obtained from multiple sources [1, 2]. Conventional MOT algorithms are usually
tailored to the point object assumption, where each object is modelled as a
point without spatial extent, and it gives rise to at most one measurement at
each time step. This assumption is, however, unrealistic for modern high-
resolution radar and Lidar sensors, for which it is common that an object
gives rise to multiple measurements per time scan. The tracking of such an
object is called extended object tracking (EOT), and overviews of EOT
literature can be found in [3, 4]. The focus of this paper is on multiple EOT.

Various extended object measurement models exist in the literature, including,
e.g., the Poisson spatial model [5], the set cluster process [6, 7], the Set
of Points on a Rigid Body model [8, 9], and the physics-based model [10].
Among these different measurement models, it is most common to use the Poisson
spatial model, where the set of object detections is modelled by an
inhomogeneous Poisson point process (PPP). The PPP measurement likelihood has
a simple factorization, which avoids explicit associations between
measurements and points on objects, thereby making it convenient to use in
EOT.

The Poisson spatial model has been used in many multiple EOT algorithms, which
can be, in general, summarized into two categories. The first category
contains methods based on random vectors, including, e.g., the multiple
hypothesis tracker (MHT) [11], the probabilistic MHT [12], the joint
(integrated) probabilistic data association filter [13, 14, 15], graph-based
filters [16, 17, 18] using belief propagation (also known as the sum-product
algorithm (SPA)), and the box particle filter [19]. The second category
contains methods based on random finite sets (RFSs) [20], including the
probability hypothesis density (PHD) filter [21, 22, 23, 24, 25, 26], the
cardinalized PHD filter [27], filters based on labelled RFSs [28, 29, 30, 31],
the Poisson multi-Bernoulli mixture (PMBM) filter [32, 33, 34, 35], and its
approximation the Poisson multi-Bernoulli (PMB) filter [36, 37].

For estimating object trajectories, MOT methods based on random vectors link
an object state estimate with a previous estimate or declare the appearance of
a new object. For RFSs-based MOT methods, one approach to estimating
trajectories is to add a unique label to each single-object state such that
each object can be identified over time [25, 28, 29, 30]. This track labelling
procedure may work well in some cases, but it often becomes problematic in
challenging scenarios [38, 39, 40]. A more advantageous approach is to compute
the multi-object posterior on sets of trajectories [40], which captures all
the information about the trajectories. In the context of EOT, this has led to
the development of the trajectory PHD filter [26] and the trajectory PMBM
(TPMBM) filter [34].

Common to all the multiple extended object filters is that the data
association problem, due to the unknown correspondences between measurements
and objects, needs to be efficiently addressed to keep the computational
complexity at a tractable level. A common way is to handle the data
association problem in two stages: first, clustering algorithms are used to
find a set of reasonable measurement partitions; second, for each partition,
the explicit assignment of measurement clusters to objects is either avoided
using the principle approximation [22, 23, 25, 27] or obtained using a 2D
assignment algorithm [28, 29, 30, 32, 34, 35, 36, 37]. These two-step
approaches based on clustering and assignment (C&A) usually work well in
scenarios with well-spaced objects, but they may suffer from a decreased
performance in scenarios where objects move in close proximity. To improve the
tracking performance in such scenarios, [24, 33] use sampling-based methods,
which work directly on the multi-object likelihood and solve the data
association problem in a single step without using clustering. The simulation
results in [33] showed that the sampling-based methods outperform methods
based on C&A in scenarios with both distant and close objects. In addition, a
method combining hierarchical clustering and Gibbs sampling is presented in
[41] for automotive applications.

For both C&A and sampling-based methods, the exhaustive enumeration of
association hypotheses is avoided by truncating hypotheses with negligible
weights. Thus, their performance might be limited in scenarios with high data
association uncertainty. A more scalable data association method for multiple
EOT, which avoids explicit enumeration of local association hypotheses, is to
directly compute the marginal multi-object posterior, where the uncertainty of
unknown data association is marginalized out [14, 15, 42, 43, 44, 16, 17]. The
current state-of-the-art multiple EOT algorithm is the SPA filter proposed in
[17] using particle belief propagation (BP), which is enabled by a factor
graph representation of the joint posterior of the multi-object states and
association variables. Simulation results in [17] show that the SPA filter
outperforms a PMBM filter that uses C&A, in terms of both estimation error and
runtime. However, the multi-object models employed in [17] are subject to two
approximations, as compared to the standard multi-object models [20]. First,
newborn objects are always detected with probability one. Second, object
survival probability is state-independent.

For the standard multi-object models [20] with PPP birth, the posterior
density is of the form PMBM [45] without approximation. The PMBM density has a
compact representation of global hypotheses with probabilistic object
existence in each Bernoulli component and undetected objects represented by a
PPP. The PMBM filtering recursions have been established for both point
objects [46] and extended objects [32], and the resulting PMBM filters have
achieved state-of-the-art performance compared to other RFSs-based
filters111An online course on MOT that explains all details regarding PMBM
filtering and introduces sets of trajectories is available at
https://www.youtube.com/channel/UCa2-fpj6AV8T6JK1uTRuFpw.. Moreover, in the
PMBM filters, the density of undetected objects is propagated over time via a
PPP, whereas in [17] and many other MOT methods using BP [47], the explicit
modelling of undetected objects is ignored, and their PPP intensity becomes
zero in the filtering density due to the modelling assumption. To provide full
trajectory information, the PMBM filtering recursions have been extended to
consider the posterior on sets of trajectories for both point objects [39] and
extended objects [34]. Furthermore, the PMBM filtering recursions are recently
extended in [35] to consider a generalized measurement model, in which each
object generates an independent set of measurements, but the filtering
recursions in [35] have not yet been further generalized to consider posterior
on sets of trajectories.

In this paper, we present the trajectory PMB (TPMB) filters for EOT, where the
marginal association probabilities and Bernoulli densities in the update step
are jointly obtained using BP. By doing so, we leverage 1) the PMBM filtering
recursions on the posterior of sets of trajectories built on the standard
extended object models and 2) the efficient data association method in [17]
using BP. As a comparison, in our previous work on extended object PMB
filtering [37], the PMB approximation is obtained by marginalizing the data
association uncertainties using truncated global hypotheses. Furthermore,
compared to previous works [47, 48, 49] that apply BP only to the data
association variables in RFSs-based methods, this paper presents a complete
pipeline for applying BP to RFSs-based methods. That is, starting from the
closed-form filtering recursion to the factor graph formulation of the joint
posterior of sets of objects and association variables, and the message
passing equations.

In our preliminary work [34], we presented the PMBM filtering recursions for
multiple EOT on posterior of sets of trajectories and an efficient track-
oriented PMBM implementation. In addition, we have proposed multiple EOT using
BP in [17] using a random vector-based derivation. This paper is a significant
extension of [17, 34], and it contains the following contributions:

  * •

We extend the PMBM conjugacy on posterior of sets of trajectories to consider
a generalized measurement model in [35], where each object generates an
independent set of measurements.

  * •

We present a full pipeline for integrating TPMB filtering and BP data
association methods. To do so, we present a factor graph representation of the
PMBM set of trajectories posterior and association variables, on an augmented
trajectory space. Importantly, by leveraging the PMBM filtering recursions, we
obtain a simpler and more general derivation of the factor graph formulation,
with undetected objects represented by a PPP, as compared to the derivation in
[50, Section 1].

  * •

We present the message passing equations, derived using RFSs without the
additional modelling assumptions made in [17], for running BP on the
constructed factor graph. We also present particle implementations for two
TPMB filters: one estimates the set of alive trajectories and the other
estimates the set of all trajectories, which includes both alive and dead
trajectories. In addition, we show how smoothed trajectory estimates, if
desired, can be further obtained using backward simulation [51].

  * •

We evaluate the performance of several implementations of TPMBM and TPMB
filters using different data association methods for EOT, including C&A,
sampling-based method and BP, with ellipsoidal shapes in a simulation study.
The results demonstrate that the TPMB filter using BP has the best trajectory
estimation performance.

The rest of the paper is organized as follows. In Section II, we introduce the
background on multi-object models, sets of trajectories and PMBM. The problem
formulation and PMBM filtering recursions are given in Section III. The factor
graph formulation and the corresponding equations for running loopy BP are
presented in Section IV and Section V, respectively. The particle
implementations of the proposed TPMB filters are provided in Section VI.
Simulation results are presented in Section VII, and conclusions are drawn in
Section VIII.

##  II Background

In this section, we first introduce the notations and the state variables of
interest, and then we give the densities/integrals for trajectories and the
general multi-object modelling. At last, we present the PMBM density
representation and the TPMB approximation. More details on the background can
be found in [40, 52].

For a generic space $D$, the set of finite subsets of $D$ is denoted by ${\cal
F}(D)$, and the cardinality of a set $A\in{\cal F}(D)$ is $|A|$. We use
$\uplus$ to denote the union of sets that are mutually disjoint, $\langle
f,g\rangle$ to denote the inner product $\int f(x)g(x)dx$, and the multi-
object exponential $f^{A}$, for some real-valued function $f$, to denote the
product $\prod_{x\in A}f(x)$ with $f^{\emptyset}=1$ by convention. In
addition, we use $\delta_{x}[\cdot]$ and $\delta_{x}(\cdot)$ to represent the
Kronecker delta and the generalized Dirac delta functions centred at $x$,
respectively [53, 54]. Here, $x$ in $\delta_{x}(\cdot)$ belongs to the space
of $\nu$ Cartesian products of space $\mathcal{X}$, i.e.,
$x\in\mathcal{X}^{\nu}$.

###  II-A State variables

The single-object state $x_{k}\in{\cal X}$ at time step $k$, where ${\cal X}$
is a locally compact, Hausdorff and second-countable (LCHS) space [20, Section
2.2.2], contains information of interest of a single object, including its
kinematic state (e.g., position and velocity) and possibly also its extent
state, which describes the object shape and size. The set of object states at
time step $k$ is an RFS ${\bf
x}_{k}=\\{x_{k}^{1},\dots,x_{k}^{n_{k}}\\}\in{\cal F}({\cal X})$. The
measurements collected by the sensor at time step $k$ is ${\bf
z}_{k}=\\{z_{k}^{1},\dots,z_{k}^{m_{k}}\\}$ with single measurement
$z_{k}^{j}\in\mathbb{R}^{d_{z}}$ for $j\in\\{1,\dots,m_{k}\\}$, and the
sequence of measurement sets up to and including time step $k$ is denoted as
${\bf z}_{1:k}$.

The trajectory of an object is a finite sequence of its states at consecutive
time steps. In this paper, we follow the notational convention in [40, 55, 52,
56] to denote a single trajectory as $X=(t,x^{1:\nu})$ where $t$ is the
initial time step of the trajectory, $\nu$ is its length and
$x^{1:\nu}=(x^{1},\dots,x^{\nu})$ is a finite sequence of object states [40,
55, 52, 56]. The variable $(t,\nu)$ belongs to the set
$I_{(k)}=\\{(t,\nu):0\leq t\leq k,1\leq\nu\leq k-t+1\\}$. A single trajectory
$X$ up to time step $k$ therefore belongs to the space
$T_{(k)}=\uplus_{(t,\nu)\in I_{(k)}}\\{t\\}\times\mathcal{X}^{\nu}$, which is
also LCHS [40, Appendix A]. The set of trajectories up to the current time
step $k$ is denoted as ${\bf X}_{k}\in{\cal F}(T_{(k)})$.

We note that trajectory $X$ is a combination of discrete and continuous
states. Such a hybrid state is not uncommon in MOT: a typical example is the
interacting multiple model [57]. We also note that the dimension of the state
sequence $x^{1:\nu}$ becomes high for long trajectories. This, however, does
not necessarily make the computation of multi-object filtering recursions
complex at each time step, see, e.g., the efficient MOT implementations in
[39, 52].

###  II-B Densities and integrals

Given a real-valued function $\pi(\cdot)$ on the single trajectory space
$T_{(k)}$, its integral is [52]

$\int\pi(X)dX=\sum_{(t,\nu)\in
I_{(k)}}\int\pi\left(t,x^{1:\nu}\right)dx^{1:\nu},$ (1)

which goes through all possible start times, lengths and object states of the
trajectory. Since the single trajectory space $T_{(k)}$ is LCHS, we can apply
Mahler’s FISST theory set integral [20, Section 3.3] to a space of sets of
trajectories ${\cal F}(T_{(k)})$ [58]. Specifically, given a real-valued
function $f(\cdot)$ on the space ${\cal F}(T_{k})$ of finite sets of
trajectories, its set integral is

$\int f({\bf X}_{k})\delta{\bf X}_{k}=\sum_{n=0}^{\infty}\frac{1}{n!}\int
f(\\{X_{1},\dots,X_{n}\\})d(X_{1},\dots,X_{n}).$ (2)

A function $f(\cdot)$ is a multi-trajectory density of a random finite set of
trajectories if $f(\cdot)\geq 0$ and its set integral is one. The reference
measures and measure theoretic integrals for sets of trajectories are defined
in [58].

###  II-C Multi-object modelling

####  II-C1 Multi-object dynamic model

At time step $k$, new objects appear in the sensor’s field of view, following
a PPP with birth intensity $\lambda^{B}_{k}(x_{k})$, independently of any
existing objects. For an existing object $x_{k-1}$, it survives with
probability $p^{S}(x_{k-1})$, and if it survives, its state evolves with a
Markovian transition density $g_{k}(x_{k}|x_{k-1})$, independently of any
other objects. It is also assumed that objects that have disappeared never
reappear.

####  II-C2 Multi-object measurement model

At time step $k$, the set ${\bf z}_{k}$ of measurements is the union of a set
of object-generated measurements and a set of clutter measurements. The
measurements from each object are independent of other objects and of clutter
measurements. The set of clutter measurements is a PPP with Poisson intensity
$\lambda_{k}^{C}(z_{k})=\gamma_{k}^{C}\mu_{k}^{C}(z_{k})$ where
$\gamma_{k}^{C}$ is the Poisson rate and $\mu_{k}^{C}(z_{k})$ is the clutter
spatial distribution.

A single object with state $x_{k}$ generates an independent set ${\bf w}_{k}$
of measurements with density $\ell_{k}({\bf w}_{k}|x_{k})$. Note that in terms
of multi-object conjugate priors, the difference between point and extended
object tracking is only how the single-object measurement likelihood
$\ell_{k}(\cdot|x_{k})$ is defined [32, 28]. For a point object, ${\bf w}_{k}$
is a Bernoulli RFS with at maximum one measurement, whereas for an extended
object, ${\bf w}_{k}$ is commonly modelled as a PPP with state dependent
Poisson rate $\gamma_{k}(x_{k})$ and spatial distribution
$\ell_{k}(\cdot|x_{k})$ [5]. For a set ${\bf w}_{k}$ of measurements generated
by an extended object with state $x_{k}$, its Poisson set density is [20,
Section 21.2.4]

$\ell_{k}({\bf w}_{k}|x_{k})=e^{-\gamma_{k}(x_{k})}\prod_{z_{k}\in{\bf
w}_{k}}\gamma_{k}(x_{k})\ell_{k}(z_{k}|x_{k}),$ (3)

and, accordingly, we have the effective probability of detection, i.e., the
probability that object $x_{k}$ generates at least one measurement

$p(|\mathbf{w}_{k}|>0|x_{k})=1-\ell_{k}(\emptyset|x_{k})=1-e^{-\gamma_{k}(x_{k})},$

where $\emptyset$ denotes an empty set.

###  II-D PMBM conjugate prior

Given the sequence of measurements ${\bf z}_{1:k}$ up to time step $k$ and the
multi-object models in Section II-C for point objects where the set of object-
generated measurement is a Bernoulli RFS, the density
$f_{k|k^{\prime}}(\cdot)$ of the set of trajectories at time step
$k^{\prime}\in\\{k-1,k\\}$ is a PMBM [39, 59] with

$\displaystyle f_{k|k^{\prime}}({\bf X}_{k})$ $\displaystyle=\sum_{{\bf
Y}\uplus{\bf V}={\bf X}_{k}}f^{p}_{k|k^{\prime}}({\bf
Y})f^{mbm}_{k|k^{\prime}}({\bf V}),$ (4) $\displaystyle
f^{p}_{k|k^{\prime}}({\bf Y})$
$\displaystyle=e^{-\left\langle\lambda_{k|k^{\prime}},1\right\rangle}\left[\lambda_{k|k^{\prime}}(\cdot)\right]^{{\bf
Y}},$ (5) $\displaystyle f^{mbm}_{k|k^{\prime}}({\bf V})$
$\displaystyle=\sum_{a\in{\cal
A}_{k|k^{\prime}}}w^{a}_{k|k^{\prime}}\sum_{\uplus_{l=1}^{n_{k|k^{\prime}}}{\bf
X}^{l}={\bf
V}}\prod_{i=1}^{n_{k|k^{\prime}}}f^{i,a^{i}}_{k|k^{\prime}}\left({\bf
X}^{i}\right),$ (6) $\displaystyle f^{i,a^{i}}_{k|k^{\prime}}\left({\bf
X}^{i}\right)$ $\displaystyle=\begin{cases}1-r^{i,a^{i}}_{k|k^{\prime}}&{\bf
X}^{i}=\emptyset\\\
r^{i,a^{i}}_{k|k^{\prime}}f^{i,a^{i}}_{k|k^{\prime}}(X)&{\bf X}^{i}=\\{X\\}\\\
0&\text{otherwise},\end{cases}$ (7)

where the prior density at time step $0$ is given by $f_{0|0}(\cdot)$. The
multi-object density $f_{k|k^{\prime}}(\cdot)$ of the form (4) is shown to be
a multi-object conjugate prior [45, 39], where it can be observed that the PPP
$f^{p}_{k|k^{\prime}}(\cdot)$ represents trajectories of objects that are
hypothesized to exist but have never been detected222The set of undetected
trajectories arises naturally from the standard multi-object models with
Poisson birth, see [45] for further explanations. The PPP of undetected
trajectories is independent of the MBM that represents trajectories detected
at some point., and that the mixture of multi-Bernoulli (MB) components
$f^{mbm}_{k|k^{\prime}}(\cdot)$ represents potential trajectories that have
been detected at least once at some point to time step $k^{\prime}$. In (5),
$\lambda_{k|k^{\prime}}(\cdot)$ denotes the PPP intensity function, whereas in
(6) each MB component describes the distribution of detected trajectories
conditioned on global association hypothesis $a\in{\cal A}_{k|k^{\prime}}$,
and its weight $w^{a}_{k|k^{\prime}}$ represents the probability of the
corresponding data association hypothesis [45]. In (6), there are
$n_{k|k^{\prime}}$ Bernoulli components, indexed by variable $i$, and for each
Bernoulli there are $h^{i}_{k|k^{\prime}}$ local hypotheses. The density of
the $i$-th Bernoulli component with local hypothesis
$a^{i}\in\\{1,\dots,h^{i}_{k|k^{\prime}}\\}$ is given by (7), where
$r^{i,a^{i}}_{k|k^{\prime}}$ is the probability of existence and
$f^{i,a^{i}}_{k|k^{\prime}}(\cdot)$ is the single-trajectory density [59]. A
global hypothesis $a=(a^{1},\dots,a^{n_{k|k^{\prime}}})\in{\cal
A}_{k|k^{\prime}}$ selects a local hypothesis for each Bernoulli component,
where the selected local hypotheses together form a valid data association and
${\cal A}_{k|k^{\prime}}$ is the set of global hypotheses that will be
described in Section III-C. The weight of global hypothesis $a$ satisfies

$w^{a}_{k|k^{\prime}}\propto\prod_{i=1}^{n_{k|k^{\prime}}}w^{i,a^{i}}_{k|k^{\prime}}$
(8)

where $w^{i,a^{i}}_{k|k^{\prime}}$ is the weight of local hypothesis $a_{i}$
for the $i$-th Bernoulli component, and normalization is required to ensure
that $\sum_{a\in{\cal A}_{k|k^{\prime}}}w^{a}_{k|k^{\prime}}=1$.

The explicit modelling of undetected objects is beneficial for automotive
applications, e.g., in traffic scenes where pedestrians walking on side walks
are occluded by parked vehicles along the street. More examples can be found
in [4, 32]. We also note that MOT algorithms based on sets of trajectories,
e.g., the TPMBM filter [39] and the TPMB filter [52], do not produce
interrupted trajectories as this is an impossible event according to the
standard models.

###  II-E Trajectory PMB approximation

The PMB is a special case of PMBM with a single MB component. We can form PMB
approximations using different techniques, see, e.g., [45, 52, 37]. In this
work, we follow [52] where the best PMB approximation is defined by using
Kullback-Leibler divergence (KLD) minimization, on a trajectory space with
auxiliary variables [52]. The KLD on the space of sets of trajectories with
auxiliary variables is an upper bound on the KLD for sets of trajectories
without auxiliary variables [52, Lemma 3].

####  II-E1 PMBM with auxiliary variables

We extend the single trajectory space with an auxiliary variable
$u\in\mathbb{U}_{k|k^{\prime}}=\\{0,1,\dots,n_{k|k^{\prime}}\\}$, such that
$(u,X)\in\mathbb{U}_{k|k^{\prime}}\times T_{(k)}$, and we denote a set of
trajectories with auxiliary variables as ${\widetilde{{\bf X}}}_{k}\in{\cal
F}(\mathbb{U}_{k|k^{\prime}}\times T_{(k)})$. Given $f_{k|k^{\prime}}(\cdot)$
of the form (4), the density $\widetilde{f}_{k|k^{\prime}}(\cdot)$ on the
space ${\cal F}(\mathbb{U}_{k|k^{\prime}}\times T_{(k)})$ can be defined as
[52]

$\widetilde{f}_{k|k^{\prime}}\left(\widetilde{{\bf
X}}_{k}\right)=\widetilde{f}^{p}_{k|k^{\prime}}\left(\widetilde{{\bf
Y}}_{k}\right)\sum_{a\in{\cal
A}_{k|k^{\prime}}}w^{a}_{k|k^{\prime}}\prod_{i=1}^{n_{k|k^{\prime}}}\widetilde{f}^{i,a^{i}}_{k|k^{\prime}}\left(\widetilde{{\bf
X}}^{i}_{k}\right)$ (9)

where, for a given $\widetilde{{\bf X}}_{k}$, $\widetilde{{\bf
Y}}_{k}=\\{(u,X)\in\widetilde{{\bf X}}_{k}:u=0\\}$ and $\widetilde{{\bf
X}}^{i}_{k}=\\{(u,X)\in\widetilde{{\bf X}}_{k}:u=i\\}$, and

$\displaystyle\widetilde{f}^{p}_{k|k^{\prime}}\left(\widetilde{{\bf
Y}}_{k}\right)$
$\displaystyle=e^{-\left\langle\widetilde{\lambda}_{k|k^{\prime}},1\right\rangle}\left[\lambda_{k|k^{\prime}}(\cdot)\right]^{\widetilde{{\bf
Y}}_{k}},$ (10a) $\displaystyle\widetilde{\lambda}_{k|k^{\prime}}(u,X)$
$\displaystyle=\delta_{0}[u]\widetilde{\lambda}_{k|k^{\prime}}(X),$ (10b)
$\displaystyle\widetilde{f}^{i,a^{i}}_{k|k^{\prime}}\left(\widetilde{{\bf
X}}_{k}^{i}\right)$
$\displaystyle=\begin{cases}1-r^{i,a^{i}}_{k|k^{\prime}}&\widetilde{{\bf
X}}_{k}^{i}=\emptyset\\\
r^{i,a^{i}}_{k|k^{\prime}}f^{i,a^{i}}_{k|k^{\prime}}(X)\delta_{i}[u]&\widetilde{{\bf
X}}_{k}^{i}=\\{(u,X)\\}\\\ 0&\text{otherwise}.\end{cases}$ (10c)

Here, the notation $\widetilde{f}_{k|k^{\prime}}(\cdot)$ is adopted to clarify
that it is the distribution on the augmented space incoporating the auxiliary
variables. We also note that the sum over sets in (4) disappears in (9) due to
the use of auxiliary variables as now there is only one possible partition of
$\widetilde{{\bf X}}_{k}$ into $\widetilde{{\bf Y}}_{k}$, $\widetilde{{\bf
X}}^{1}_{k},\dots,\widetilde{{\bf X}}^{n_{k|k^{\prime}}}_{k}$ that provides a
non-zero density.

####  II-E2 PMB approximation

Given a PMBM density $\widetilde{f}_{k|k^{\prime}}(\cdot)$ of the form (9),
the PMB density that minimizes the KLD from
$\widetilde{f}_{k|k^{\prime}}(\cdot)$ has a closed form and is given by [52,
Proposition 2]

$\displaystyle\widetilde{f}^{pmb}_{k|k^{\prime}}(\widetilde{{\bf X}}_{k})$
$\displaystyle=\widetilde{f}^{p}_{k|k^{\prime}}(\widetilde{{\bf
Y}}_{k})\prod_{i=1}^{n_{k|k^{\prime}}}\widetilde{f}^{i}_{k|k^{\prime}}\left(\widetilde{{\bf
X}}^{i}_{k}\right),$ (11a)
$\displaystyle\widetilde{f}^{i}_{k|k^{\prime}}(\widetilde{{\bf X}})$
$\displaystyle=\begin{cases}1-r^{i}_{k|k^{\prime}}&\widetilde{{\bf
X}}=\emptyset\\\
r^{i}_{k|k^{\prime}}f^{i}_{k|k^{\prime}}(X)\delta_{i}[u]&\widetilde{{\bf
X}}=\\{(u,X)\\}\\\ 0&\text{otherwise},\end{cases}$ (11b) $\displaystyle
r^{i}_{k|k^{\prime}}$
$\displaystyle=\sum_{a^{i}=1}^{h^{i}_{k|k^{\prime}}}\overline{w}_{k|k^{\prime}}^{i,a^{i}}r^{i,a^{i}}_{k|k^{\prime}},$
(11c) $\displaystyle f^{i}_{k|k^{\prime}}(X)$
$\displaystyle=\frac{\sum_{a^{i}=1}^{h^{i}_{k|k^{\prime}}}\overline{w}_{k|k^{\prime}}^{i,a^{i}}r^{i,a^{i}}_{k|k^{\prime}}f_{k|k^{\prime}}^{i,a^{i}}(X)}{\sum_{a^{i}=1}^{h^{i}_{k|k^{\prime}}}\overline{w}_{k|k^{\prime}}^{i,a^{i}}r^{i,a^{i}}_{k|k^{\prime}}},$
(11d) $\displaystyle\overline{w}_{k|k^{\prime}}^{i,a^{i}}$
$\displaystyle=\sum_{b\in{\cal
A}_{k|k^{\prime}}:b^{i}=a^{i}}w^{b}_{k|k^{\prime}}.$ (11e)

If we integrate out the auxiliary variables in (11), we obtain the PMB density
without auxiliary variables [52]

${f}^{pmb}_{k|k^{\prime}}({{\bf
X}}_{k})=\sum_{\uplus_{l=1}^{n_{k|k^{\prime}}}{\bf X}^{l}\uplus{\bf Y}={\bf
X}_{k}}{f}^{p}_{k|k^{\prime}}({{\bf
Y}})\prod_{i=1}^{n_{k|k^{\prime}}}{f}^{i}_{k|k^{\prime}}\left({{\bf
X}}^{i}\right)$ (12)

where the $i$-th Bernoulli component $f^{i}_{k|k^{\prime}}({\bf X}^{i})$ is
parameterized by $r^{i}_{k|k^{\prime}}$ in (11c) and $f^{i}_{k|k^{\prime}}(X)$
in (11d).

##  III Problem Formulation and Filtering Recursions

In this section, we present the problem formulations and the PMBM filtering
recursions based on the multi-object models introduced in Section II-C with a
generalized measurement model.

###  III-A Problem formulation

We consider the following two MOT problem formulations:

  1. 1.

The set ${\bf X}_{k}$ of alive trajectories that exist at the current time
step, i.e., $t+\nu-1=k$ for each $(t,x^{1:\nu})\in{\bf X}_{k}$.

  2. 2.

The set ${\bf X}_{k}$ of all trajectories that have existed up to the current
time step $k$, i.e., $t+\nu-1\leq k$ for each $(t,x^{1:\nu})\in{\bf X}_{k}$.

It should be noted that the set of trajectories included in both 1) and 2)
does not depend on the sensor model, but on the object dynamic model in
Section II-C. Therefore, trajectories of objects that move outside the
sensor’s field-of-view, and appear in the field-of-view afterwards will always
belong to 2), and will belong to 1) if the dynamic model indicates they are
alive.

For both problem formulations, the objective is to compute the posterior
density of the set ${\bf X}_{k}$ of trajectories at time step $k$ given the
sequence ${\bf z}_{1:k}$ of measurements up to time step $k$ using Bayesian
recursive filtering.

###  III-B Bayesian models for sets of trajectories

We proceed to give the Bayesian multi-trajectory dynamic and measurement
models for the two types of problem formulations (see also [40, 52]), required
for deriving the Bayesian filtering recursions. The two multi-trajectory
dynamic models extend the multi-object dynamic model in Section II-C, where
each object evolves according to a Markov process, to trajectories.

####  III-B1 Dynamic model for sets of alive trajectories

Given the set ${\bf X}_{k-1}$ of alive trajectories at time step $k-1$, each
$X=(t,x^{1:\nu})\in{\bf X}_{k-1}$ survives with probability
$p^{S}(X)=p^{S}(x^{\nu})$ (i.e., the survival probability depends only on the
final state in the trajectory), and if it survives, its state $(t,x^{1:\nu})$
evolves with a transition density [52, Sec. IV-A-1]

$\displaystyle g_{k}\left(\bar{t},y^{1:\bar{\nu}}|X\right)$
$\displaystyle=\delta_{t}[\bar{t}]\delta_{\nu+1}[\bar{\nu}]\delta_{x^{1:\nu}}\left(y^{1:\bar{\nu}-1}\right)g_{k}\left(y^{\bar{\nu}}|x^{\nu}\right),$
(13)

where the state $y^{\bar{\nu}}$ only depends on $x^{\nu}$ due to the Markov
property.

The set ${\bf X}_{k}$ is the union of the surviving trajectories and the PPP
newborn trajectories with birth intensity

$\lambda^{B}_{k}\left(t,x^{1:\nu}\right)=\delta_{k}[t]\delta_{1}[\nu]\lambda^{B}_{k}\left(x^{\nu}\right).$
(14)

####  III-B2 Dynamic model for sets of all trajectories

Given the set ${\bf X}_{k-1}$ of all trajectories at time step $k-1$, each
$X=(t,x^{1:\nu})\in{\bf X}_{k-1}$ survives with probability $p^{S}(X)=1$ and
evolves with a transition density [52, Sec. IV-A-2]

$\displaystyle
g_{k}\left(\bar{t},y^{1:\bar{\nu}}|X\right)=\delta_{t}[\bar{t}]$
$\displaystyle\times\begin{cases}\delta_{\nu}[\bar{\nu}]\delta_{x^{1:\nu}}\left(y^{1:\bar{\nu}}\right)&\bar{\omega}<k-1\\\
\delta_{\nu}[\bar{\nu}]\delta_{x^{1:\nu}}\left(y^{1:\bar{\nu}}\right)\left(1-p^{S}\left(x^{\nu}\right)\right)&\bar{\omega}=k-1\\\
\delta_{\nu+1}[\bar{\nu}]\delta_{x^{1:\nu}}\left(y^{1:\bar{\nu}-1}\right)p^{S}\left(x^{\nu}\right)g_{k}\left(y^{\bar{\nu}}|x^{\nu}\right)&\bar{\omega}=k\end{cases}$
(15)

where $\bar{\omega}=\bar{t}+\bar{\nu}-1$ and the state $y^{\bar{\nu}}$ only
depends on $x^{\nu}$ due to the Markov property. The birth model is also a PPP
with intensity (14). Note that $p^{S}(X)$ refers to the probability that a
trajectory remains in the considered set of trajectories for both problem
formulations.

####  III-B3 Measurement model for sets of trajectories

The measurement model is the same for sets of alive trajectories and the sets
of all trajectories. Each trajectory $(t,x^{1:\nu})\in{\bf X}_{k}$ generates a
set ${\bf w}_{k}$ of measurements with density $\ell_{k}({\bf
w}_{k}|t,x^{1:\nu})=\ell_{k}({\bf w}_{k}|x^{\nu})$. The clutter model is a PPP
as described in Section II-C2.

###  III-C Filtering recursions

We present the prediction and update for the multi-trajectory dynamic and
measurement models in Section III-B.

######  Proposition 1.

Given the PMBM filtering density on the set of trajectories at time step $k-1$
of the form (4), the predicted density at time step $k$ is a PMBM of the form
(4), with

$\displaystyle\lambda_{k|k-1}(X)$
$\displaystyle=\lambda^{B}_{k}(X)+\left\langle\lambda_{k-1|k-1},g_{k}(X|\cdot)p^{S}(\cdot)\right\rangle,$
(16a) $\displaystyle n_{k|k-1}$ $\displaystyle=n_{k-1|k-1},$ (16b)
$\displaystyle h^{i}_{k|k-1}$ $\displaystyle=h^{i}_{k-1|k-1},$ (16c)
$\displaystyle w^{i,a^{i}}_{k|k-1}$ $\displaystyle=w^{i,a^{i}}_{k-1|k-1},$
(16d) $\displaystyle r^{i,a^{i}}_{k|k-1}$
$\displaystyle=r^{i,a^{i}}_{k-1|k-1}\left\langle
f^{i,a^{i}}_{k-1|k-1},p^{S}\right\rangle,$ (16e) $\displaystyle
f^{i,a^{i}}_{k|k-1}(X)$ $\displaystyle=\frac{\left\langle
f^{i,a^{i}}_{k-1|k-1},g_{k}(X|\cdot)p^{S}(\cdot)\right\rangle}{\left\langle
f^{i,a^{i}}_{k-1|k-1},p^{S}\right\rangle}$ (16f)

where $g_{k}(\cdot|\cdot)$ and $p^{S}(\cdot)$ are chosen depending on the
problem formulation.

The TPMBM prediction given in Proposition 1 is the same as the one in [39] as
it is not affected by the choice of the measurement model.

Before presenting the update step, it is important to define the set of
feasible global hypotheses for a generalized measurement model. We refer to
measurement $z_{k}^{j}$ using the pair $(k,j)$ and the set of all such
measurement pairs up to (and including) time step $k$ is denoted by ${\cal
M}_{k}$. Then, a local hypothesis $a^{i}$ for the $i$-th Bernoulli component
has a set of measurement pairs denoted as ${\cal
M}_{k}^{i,a^{i}}\subseteq{\cal M}_{k}$, and the set ${\cal A}_{k|k^{\prime}}$
of all global hypotheses satisfies

${\cal
A}_{k|k^{\prime}}=\left\\{\left(a^{1},\dots,a^{n_{k|k^{\prime}}}\right):a^{i}\in\left\\{1,\dots,h^{i}_{k|k^{\prime}}\right\\}~{}\forall~{}i,\right.\\\
\left.\biguplus_{i=1}^{n_{k|k^{\prime}}}{\cal M}_{k^{\prime}}^{i,a^{i}}={\cal
M}_{k^{\prime}}\right\\}$ (17)

where the recursive constructions of ${\cal M}_{k^{\prime}}^{i,a^{i}}$ and
${\cal M}_{k^{\prime}}$, with ${\cal M}_{0}=\emptyset$, will be given in
Proposition 2.

We note that every global hypothesis needs to explain the association of every
measurement that has been received so far. The same measurement cannot be
associated to more than one local hypothesis, but more than one measurement
may be associated to the same local hypothesis at the same time step. Each
global hypothesis therefore corresponds to a unique partition of ${\cal
M}_{k^{\prime}}$, and the number of global hypotheses is then given by the
Bell number of $|{\cal M}_{k^{\prime}}|$ [32]. At each time step, each
measurement generates a new Bernoulli component, corresponding to an
undetected object detected for the first time or clutter.

######  Proposition 2.

Given the PMBM predicted density on the set of trajectories at time step $k$
of the form (4), and measurements ${\bf
z}_{k}=\\{z_{k}^{1},\dots,z_{k}^{m_{k}}\\}$, the updated density is a PMBM of
the form (4), with

$\displaystyle n_{k|k}$ $\displaystyle=n_{k|k-1}+m_{k},$ (18)
$\displaystyle{\cal M}_{k}$ $\displaystyle={\cal
M}_{k-1}\cup\left\\{(k,j)|j\in\\{1,\dots,m_{k}\\}\right\\},$ (19)
$\displaystyle\lambda_{k|k}(X)$
$\displaystyle=\ell_{k}(\emptyset|X)\lambda_{k|k-1}(X).$ (20)

For Bernoulli components continuing from previous time steps
$i\in\\{1,\dots,n_{k|k-1}\\}$, a new local hypothesis is included for each
previous local hypothesis and either a misdetection or an update with a non-
empty subset of ${\bf z}_{k}$. The updated number of local hypotheses is
$h^{i}_{k|k}=2^{m_{k}}h^{i}_{k|k-1}$. For missed detection hypotheses,
$i\in\\{1,\dots,n_{k|k-1}\\}$, $a^{i}\in\\{1,\dots,h^{i}_{k|k-1}\\}$, we have

$\displaystyle{\cal M}_{k}^{i,a^{i}}$ $\displaystyle={\cal
M}_{k-1}^{i,a^{i}},$ (21a) $\displaystyle\ell_{k|k}^{i,a^{i},0}$
$\displaystyle=\left\langle
f_{k|k-1}^{i,a^{i}},\ell_{k}(\emptyset|\cdot)\right\rangle,$ (21b)
$\displaystyle w_{k|k}^{i,a^{i}}$
$\displaystyle=w^{i,a^{i}}_{k|k-1}\left(1-r^{i,a^{i}}_{k|k-1}+r^{i,a^{i}}_{k|k-1}\ell_{k|k}^{i,a^{i},0}\right),$
(21c) $\displaystyle r^{i,a^{i}}_{k|k}$
$\displaystyle=\frac{r^{i,a^{i}}_{k|k-1}\ell_{k|k}^{i,a^{i},0}}{1-r^{i,a^{i}}_{k|k-1}+r^{i,a^{i}}_{k|k-1}\ell_{k|k}^{i,a^{i},0}},$
(21d) $\displaystyle f^{i,a^{i}}_{k|k}(X)$
$\displaystyle=\frac{\ell_{k}(\emptyset|X)f^{i,a^{i}}_{k|k-1}(X)}{\ell_{k|k}^{i,a^{i},0}}.$
(21e)

Let ${\bf w}_{k}^{1},\dots,{\bf w}_{k}^{2^{m_{k}}-1}$ be the non-empty subsets
of ${\bf z}_{k}$. For the $i$-th Bernoulli component,
$i\in\\{1,\dots,n_{k|k-1}\\}$, with a local hypothesis
$\tilde{a}^{i}\in\\{1,\dots,h^{i}_{k|k-1}\\}$ in the predicted density, the
new local hypothesis generated by a set ${\bf w}^{j}$ of measurements has
$a^{i}=\tilde{a}^{i}+h^{i}_{k|k-1}j$ with $j\in\\{1,\dots,2^{m_{k}}-1\\}$, and

$\displaystyle{\cal M}_{k}^{i,a^{i}}$ $\displaystyle={\cal
M}_{k-1}^{i,\tilde{a}^{i}}\cup\left\\{(k,p):z_{k}^{p}\in{\bf
w}_{k}^{j}\right\\},$ (22a) $\displaystyle\ell_{k|k}^{i,a^{i},j}$
$\displaystyle=\left\langle f_{k|k-1}^{i,\tilde{a}^{i}},\ell_{k}\left({\bf
w}^{j}_{k}|\cdot\right)\right\rangle,$ (22b) $\displaystyle w_{k|k}^{i,a^{i}}$
$\displaystyle=w_{k|k-1}^{i,\tilde{a}^{i}}r_{k|k-1}^{i,\tilde{a}^{i}}\ell_{k|k}^{i,a^{i},j},$
(22c) $\displaystyle r_{k|k}^{i,a^{i}}$ $\displaystyle=1,$ (22d)
$\displaystyle f_{k|k}^{i,a^{i}}(X)$ $\displaystyle=\frac{\ell_{k}\left({\bf
w}_{k}^{j}|X\right)f_{k|k-1}^{i,\tilde{a}^{i}}(X)}{\ell_{k|k}^{i,a^{i},j}}.$
(22e)

Each new Bernoulli component has a different number of local hypotheses, each
of which is created by a subset of ${\bf z}_{k}$. Let the set ${\cal S}_{i}$
of subsets of measurements associated to the $i$-th new Bernoulli component
($i\in\\{1,\dots,m_{k}\\}$) be recursively built as

${\cal
S}_{i}=\left\\{\left\\{{z}_{k}^{i}\right\\}\right\\}\cup\left(\bigcup_{{\bf
w}\in\cup_{j=1}^{i-1}{\cal S}_{i}}\left\\{\left\\{{z}_{k}^{i}\right\\}\cup{\bf
w}\right\\}\right)$ (23)

with ${\cal S}_{1}=\\{\\{{z}_{k}^{1}\\}\\}$, where $\\{\\{{z}_{k}^{i}\\}\\}$
denotes a set of subsets, containing the single element set $\\{z^{i}_{k}\\}$.
According to (23), the $i$-th new Bernoulli component has $2^{i-1}$ local
hypotheses that are created by non-empty subsets of measurements. Further, let
${\bf w}_{k}^{i,\iota}$ denote the $\iota$-th subset of measurements of the
$i$-th Bernoulli component
($i=n_{k|k-1}+j,j\in\\{1,\dots,m_{k}\\},\iota\in\\{1,\dots,2^{j-1}\\}$). Then,
for the $i$-th Bernoulli component, there are $h^{i}_{k|k}=2^{j-1}+1$ local
hypotheses, one corresponding to a non-existent Bernoulli component
($a^{i}=1$)

${\cal M}_{k}^{i,1}=\emptyset,w_{k|k}^{i,1}=1,r_{k|k}^{i,1}=0,$ (24)

and the others ($a^{i}=\iota+1,\iota\in\\{1,\dots,2^{j-1}\\}$) have

$\displaystyle{\cal M}_{k}^{i,a^{i}}$
$\displaystyle=\left\\{(k,p):z_{k}^{p}\in{\bf w}_{k}^{i,\iota}\right\\},$
(25a) $\displaystyle\ell_{k|k}^{i,\iota}$
$\displaystyle=\left\langle\lambda_{k|k-1},\ell_{k}\left({\bf
w}_{k}^{i,\iota}|\cdot\right)\right\rangle,$ (25b) $\displaystyle
w_{k|k}^{i,a^{i}}$ $\displaystyle=\delta_{1}\left[\left|{\bf
w}_{k}^{i,\iota}\right|\right]\left[\lambda_{k}^{C}\right]^{{\bf
w}_{k}^{i,\iota}}+\ell_{k|k}^{i,\iota},$ (25c) $\displaystyle
r_{k|k}^{i,a^{i}}$
$\displaystyle=\frac{\ell_{k|k}^{i,\iota}}{w_{k|k}^{i,a^{i}}},$ (25d)
$\displaystyle f_{k|k}^{i,a^{i}}(X)$ $\displaystyle=\frac{\ell_{k}\left({\bf
w}_{k}^{i,\iota}|X\right)\lambda_{k|k-1}(X)}{\ell_{k|k}^{i,\iota}}.$ (25e)

Figure 1: An illustration of the local hypothesis structure in Proposition 2
for the case $n_{k|k-1}=0$ and $m_{k}=3$ where the time step $k$ is omitted.
In this case, three Bernoulli components are created, one for each
measurement. The local hypotheses under these components are: 1) $\emptyset$
and $\\{(k,1)\\}$, 2) $\emptyset$, $\\{(k,2)\\}$ and $\\{(k,1),(k,2)\\}$, 3)
$\emptyset$, $\\{(k,3)\\}$, $\\{(k,1),(k,3)\\}$, $\\{(k,2),(k,3)\\}$ and
$\\{(k,1),(k,2),(k,3)\\}$.

Proposition 2 is an extension of the PMBM update for sets of objects [35,
Theorem 1] to sets of trajectories, and a short proof of Proposition 2 is
given in Appendix A.

Compared to the local hypothesis representation in [35, Theorem 1], in
Proposition 2 each measurement, instead of each non-empty subset of
measurements, creates a new Bernoulli component, and the $i$-th new potential
object cannot generate more than $i$ measurements, see Fig. 1 for an example.
These two different local hypothesis representations lead to the same PMBM
posterior, but the one in Proposition 2 facilitates the development of
efficient PMB approximation methods [37].

The established PMBM conjugacy on posterior of sets of trajectories for a
generalized object measurement model paves the way for developing efficient
inference methods on sets of trajectories using more general object
measurement models for both online and offline applications, using, e.g., BP
[47] and Gibbs sampling [60]. A typical example is the extended object
trajectory PMBM filter presented in our preliminary work [34]. In the next
section, we will elaborate on how to derive the factor graph formulation of
the update step by making use of Proposition 2, with the objective to later
develop BP algorithm for online MOT.

##  IV Factor Graph Formulation

The TPMB filters are computationally lighter alternatives to the TPMBM
filters. The prediction of the TPMB filter is a special case of the TPMBM
prediction with only one mixture component. The update of the TPMB filter is
obtained by first performing a Bayesian update, which yields a PMBM
distribution, followed by a PMB approximation. The main challenge of the PMB
approximation is the efficient calculation of the marginal density of set of
trajectories, where the uncertainty on global hypotheses has been marginalized
out. To address this problem, the existing works [35, 37] first prune global
hypotheses with negligible weights using either C&A [34, 32] or sampling-based
methods [33], and then approximately compute the marginal association
probabilities (11e) by enumerating the truncated global hypotheses. Also, we
note that, even if one manages to compute (11e), it may still be
computationally challenging to evaluate (11d) in closed form, see, e.g., [61,
62].

The present work is motivated by an emergent approximation method for multiple
EOT with Poisson spatial model that is both accurate and scalable [17]. The
approximation is based on BP, which jointly computes the marginal density of
object states as well as the marginal association probability of each
measurement, without explicitly enumerating the associations between
measurement clusters and objects. To apply the techniques in [16, 17] to PMB
approximation, we need to represent the joint posterior of set of trajectories
and association variables using a factor graph.

In this section, we give the factor graph representation of the joint
posterior of sets of trajectories and association variables, required for
developing efficient TPMB filters with Poisson spatial model using BP.

###  IV-A Factor graph formulation

Figure 2: Factor graph for the factorization in (3) where variable nodes are
represented using square and factor nodes are represented using circle. The
messages passing in loopy BP via $\widetilde{{\bf X}}_{k}^{1}$ and
$\widetilde{{\bf X}}_{k}^{n_{k|k-1}+m_{k}}$ are also illustrated (in red
lines) where their arguments and iteration indices are omitted for brevity.

Following [17], we use the measurement-oriented association vector
$\beta_{k}=[\beta^{1}_{k},\dots,\beta_{k}^{m_{k}}]^{T}$ where
$\beta_{k}^{j}\in\\{0,1,\dots,n_{k|k}\\}$ for $j\in\\{1,\dots,m_{k}\\}$, to
describe the object-to-measurement associations. Specifically,
$\beta_{k}^{j}=i>0$ if and only if the $j$-th measurement $z_{k}^{j}$ is
associated to the $i$-th potential object (described by the $i$-th Bernoulli
component), whereas $\beta_{k}^{j}=0$ if the $j$-th measurement is a false
alarm and thus not associated to a potential object. According to the local
hypothesis representation for new Bernoulli components in Proposition 2, the
$j$-th measurement $z_{k}^{j}$ cannot be associated to new potential objects
with index $i\in\\{n_{k-1|k}+1,\dots,n_{k-1|k}+j-1\\}$, i.e.,
$\beta_{k}^{j}\in\\{0,1,\dots,n_{k-1|k},n_{k-1|k}+j,\dots,n_{k|k}\\}$.

######  Proposition 3.

Given a predicted TPMB density at time step $k$ of the form (11a), measurement
set ${\bf z}_{k}$, and the PPP measurement likelihood

$\displaystyle\ell_{k}\left({\bf w}_{k}|X\right)$
$\displaystyle=e^{-\gamma_{k}(X)}\prod_{z_{k}\in{\bf
w}_{k}}\gamma_{k}(X)\ell_{k}(z_{k}|X),$ (26a)
$\displaystyle\gamma_{k}\left(t,x^{1:\nu}\right)$
$\displaystyle=\gamma_{k}\left(x^{\nu}\right)\delta_{k}[t+\nu-1],$ (26b)
$\displaystyle\ell_{k}\left(z_{k}|t,x^{1:\nu}\right)$
$\displaystyle=\ell_{k}\left(z_{k}|x^{\nu}\right)\delta_{k}[t+\nu-1],$ (26c)

the joint posterior of the set of trajectories and measurement-oriented
association vector is

$\displaystyle\widetilde{f}_{k|k}\left(\widetilde{{\bf
X}}_{k},\beta_{k}\right)$
$\displaystyle\propto\widetilde{f}^{p}_{k|k}\left(\widetilde{{\bf
Y}}_{k}\right)\prod_{i=1}^{n_{k|k-1}}\Biggl{[}\underline{f}^{i}_{k|k-1}\left(\widetilde{{\bf
X}}^{i}_{k}\right)\prod_{j=1}^{m_{k}}\underline{s}_{k}\left(\widetilde{{\bf
X}}_{k}^{i},\beta_{k}^{j};z_{k}^{j}\right)\Biggr{]}$
$\displaystyle~{}~{}\times\prod_{i=n_{k|k-1}+1}^{n_{k|k}}\Biggl{[}\overline{f}^{i}_{k|k-1}\left(\widetilde{{\bf
X}}^{i}_{k}\right)\overline{s}_{k}\left(\widetilde{{\bf
X}}_{k}^{i},\beta_{k}^{i-n_{k|k-1}};z_{k}^{i-n_{k|k-1}}\right)$
$\displaystyle~{}~{}\times\prod_{j=1}^{i-n_{k|k-1}-1}\underline{s}_{k}\left(\widetilde{{\bf
X}}_{k}^{i},\beta_{k}^{j};z_{k}^{j}\right)\Biggr{]}$ (27)

where $\widetilde{\bf Y}_{k}$ and $\widetilde{\bf X}_{k}^{i}$ with
$i\in\\{1,\dots,n_{k|k}\\}$ are defined as in (9), and

$\displaystyle\underline{f}^{i}_{k|k-1}\left(\widetilde{{\bf X}}\right)$
$\displaystyle=\begin{cases}r^{i}_{k|k-1}f^{i}_{k|k-1}(X)e^{-\gamma_{k}(X)}\delta_{i}[u]&\widetilde{{\bf
X}}=\\{(u,X)\\}\\\ 1-r^{i}_{k|k-1}&\widetilde{{\bf X}}=\emptyset\\\
0&\text{otherwise},\end{cases}$ (28)
$\displaystyle\overline{f}^{i}_{k|k-1}\left(\widetilde{{\bf X}}\right)$
$\displaystyle=\begin{cases}\lambda_{k|k-1}(X)e^{-\gamma_{k}(X)}\delta_{i}[u]&\widetilde{{\bf
X}}=\\{(u,X)\\}\\\ 1&\widetilde{{\bf X}}=\emptyset\\\
0&\text{otherwise},\end{cases}$ (29)
$\displaystyle\underline{s}_{k}\left(\widetilde{{\bf
X}}_{k}^{i},\beta_{k}^{j};z_{k}^{j}\right)$
$\displaystyle=\begin{cases}\frac{\ell_{k}\left(z_{k}^{j}|X\right)\gamma_{k}(X)}{\lambda_{k}^{C}\left(z_{k}^{j}\right)}\delta_{i}[u]&\widetilde{{\bf
X}}_{k}^{i}=\\{(u,X)\\},\beta_{k}^{j}=i\\\ 1&\beta_{k}^{j}\neq i\\\
0&\text{otherwise},\end{cases}$ (30)
$\displaystyle\overline{s}_{k}\left(\widetilde{{\bf
X}}_{k}^{i},\beta_{k}^{j};z_{k}^{j}\right)$
$\displaystyle=\begin{cases}\frac{\ell_{k}\left(z_{k}^{j}|X\right)\gamma_{k}(X)}{\lambda_{k}^{C}\left(z_{k}^{j}\right)}\delta_{i}[u]&\widetilde{{\bf
X}}_{k}^{i}=\\{(u,X)\\},\beta_{k}^{j}=i\\\ 1&\widetilde{{\bf
X}}_{k}^{i}=\emptyset,\beta_{k}^{j}\neq i\\\ 0&\text{otherwise}.\end{cases}$
(31)

The expression of the joint posterior (3) consists of three factors, including
$\widetilde{f}^{p}_{k|k}(\cdot)$ and two products. The factor
$\widetilde{f}^{p}_{k|k}(\cdot)$ describes the posterior density of the set of
undetected trajectories. The first product over $i\in\\{1,\dots,n_{k|k-1}\\}$
describes the joint posterior of the set of existing trajectories and the
measurement-oriented association vector. The second product over
$i\in\\{n_{k|k-1}+1,\dots,n_{k|k}\\}$ describes the joint posterior of the set
of newly detected trajectories and the measurement-oriented association
vector.

Each factor in the two products in (3) can be further factorized as a product
over several smaller factors. In particular, factors (28) and (29) describe
the set $\widetilde{\bf X}^{i}_{k}$ of trajectories prior without relating to
any measurements, respectively for existing and newly detected objects.
Factors (30) and (31) reflect the likelihood of associating a particular
measurement with its corresponding Bernoulli component, respectively for
existing and newly detected objects. We note that an object may have several
associated measurements; the factorized likelihood in (3) results in the
factorized form of (3).

Proposition 3 is proved in Appendix B, and the factor graph corresponding to
(3) is illustrated in Fig. 2. By leveraging Proposition 2 and RFSs, we can
have a concise derivation of the factor graph formulation of the joint
posterior, which is precisely built upon the multi-object models given in
Section II-C and the Poisson spatial measurement model (26) without any
further assumptions, as compared to the long and involved derivation given in
[50, Section 1].

###  IV-B Relation to the factor graph in [17]

The factor graph formulation in [17] is for the joint posterior of set of
object states and measurement-oriented association variables, derived using
random vectors. As a comparison, in this work, we consider inference on sets
of trajectories and the factor graph formulation is derived using RFSs. Also
note that the factor graph in [17] also includes the time evolution of multi-
object dynamics but with the assumption that the object survival probability
$p^{S}$ is state-independent, whereas the factor graph in Fig. 2 only focuses
on the update step. It is also possible to derive a factor graph formulation
that considers the joint posterior of set of trajectories and association
variables until the current time step, but this is not required for the
development of TPMB filters using BP.

Another important difference is that in the factor graph shown in Fig. 2 we
have one additional variable node and factor node for the set of undetected
trajectories, which is a direct result of the explicit modelling of undetected
trajectories in the PMBM conjugate prior. This is not the case in [17], where
it is assumed that newborn objects are detected with probability one such that
the posterior of the set of undetected objects is zero, and that the prior
distribution for newborn objects (cf. (29)) needs to be normalized by
$1-\ell_{k}(\emptyset|\cdot)$ to exclude the case that a newborn object may
generate zero measurement. The difference on the predicted number of
undetected objects at time step $k$ is further elaborated in Appendix C.

Note that both factor graph formulations are not unique. The one in [17]
depends on the mapping between measurements and new potential objects, and a
different mapping may result in a different factorization of the posterior.
The factorization in (3) depends on the local hypothesis representation for
new Bernoulli components since there is a one-to-one mapping between global
hypothesis and measurement-oriented association vector. There is a family of
hypothesis representations of the posterior that only differ in the
representation of new Bernoulli components [37, Theorem 1]. Different
hypothesis representations for the new Bernoulli components yield different
joint posteriors on the space of auxiliary variables, and therefore the
corresponding factorizations are also different.

##  V Loopy Belief Propagation for EOT

The objective is to compute the marginal density of each Bernoulli set
$\widetilde{{\bf X}}^{i}_{k}$ of trajectory with $i\in\\{1,\dots,n_{k|k}\\}$.
Since the factor graph in Fig. 2 has cycles, we apply loopy BP for inference,
where at each time step, messages are computed and processed in parallel [17].
Loopy BP has been applied to solve multiple point object tracking problems
[47]. Since the factor graph formulations are different for point and extended
objects, the message passing equations are also different. In this section, we
present the generic equations for message passing and belief calculation for
the factor graph in Fig. 2. These equations resemble the message passing
equations in [17], but they are derived using RFSs. For notational brevity, we
omit the time index for all the messages.

###  V-A Iterative message passing

The notations for the different messages are given as follows. At iteration
$p\in\\{1,\dots,P\\}$, we denote the message from variable node
$\widetilde{{\bf X}}_{k}^{i}$, $i\in\\{1,\dots,n_{k|k}\\}$, to factor node
$\overline{s}_{k}(\widetilde{{\bf X}}_{k}^{i},\beta_{k}^{j};z_{k}^{j})$ or
$\underline{s}_{k}(\widetilde{{\bf X}}_{k}^{i},\beta_{k}^{j};z_{k}^{j})$,
$j\in\\{1,\dots,m_{k}\\}$, by $\varepsilon_{i,j}^{(p)}(\cdot)$ and the message
in the reverse direction by $\varrho_{j,i}^{(p)}(\cdot)$. In addition, we also
denote the message from factor node $\overline{s}_{k}(\widetilde{{\bf
X}}_{k}^{i},\beta_{k}^{j};z_{k}^{j})$ or $\underline{s}_{k}(\widetilde{{\bf
X}}_{k}^{i},\beta_{k}^{j};z_{k}^{j})$ to variable node $\beta_{k}^{j}$ by
$\vartheta_{i,j}^{(p)}(\cdot)$ and the message in the reverse direction by
$\varphi_{j,i}^{(p)}(\cdot)$.

####  V-A1 Initialization

For $p=1$, we set
$\varepsilon_{i,j}^{(1)}(\cdot)=\underline{f}^{i}_{k|k-1}(\cdot)$ for
$i\in\\{1,\dots,n_{k|k-1}\\}$ and
$\varepsilon_{i,j}^{(1)}(\cdot)=\overline{f}^{i}_{k|k-1}(\cdot)$ for
$i\in\\{n_{k|k-1}+1,\dots,n_{k|k}\\}$, where $n_{k|k}$ is given by (18).
Further, we introduce
$\varepsilon^{(p)}_{i,j}=\int\varepsilon^{(p)}_{i,j}\left(\\{X\\}\right)dX+\varepsilon^{(p)}_{i,j}(\emptyset)$,
where for $p=1$, $i\in\\{1,\dots,n_{k|k-1}\\}$, we have

$\displaystyle\varepsilon^{(1)}_{i,j}$
$\displaystyle=\int\underline{f}^{i}_{k|k-1}\left(\widetilde{{\bf
X}}\right)\delta\widetilde{{\bf X}}$
$\displaystyle=\sum_{u\in\mathbb{U}_{k|k}}\int\underline{f}^{i}_{k|k-1}(\\{(u,X)\\})dX+\underline{f}^{i}_{k|k-1}(\emptyset)$
$\displaystyle=r^{i}_{k|k-1}\int
f^{i}_{k|k-1}(X)e^{-\gamma_{k}(X)}dX+1-r^{i}_{k|k-1},$ (32)

and for $i\in\\{n_{k|k-1}+1,\dots,n_{k|k}\\}$,

$\displaystyle\varepsilon^{(1)}_{i,j}$
$\displaystyle=\int\overline{f}^{i}_{k|k-1}\left(\widetilde{{\bf
X}}\right)\delta\widetilde{{\bf X}}$
$\displaystyle=\sum_{u\in\mathbb{U}_{k|k}}\int\overline{f}_{k|k-1}^{i}(\\{(u,X)\\})dX+\overline{f}_{k|k-1}(\emptyset)$
$\displaystyle=\int\lambda_{k|k-1}(X)e^{-\gamma_{k}(X)}dX+1.$ (33)

By comparing the initialization step with the update of local hypotheses in
Proposition 2, it is not difficult to observe that the initialization of
message $\varepsilon_{i,j}^{(1)}(\cdot)$, $i\in\\{1,\dots,n_{k|k-1}\\}$,
corresponds to the unnormalized misdetection hypothesis density for the $i$-th
Bernoulli component, and that the initialization of message
$\varepsilon_{i,j}^{(1)}(\cdot)$, $i\in\\{n_{k|k-1},\dots,n_{k|k}\\}$,
corresponds to the unnormalized density of the second local hypothesis of the
$i$-th Bernoulli component but with $\ell_{k}(\emptyset|\cdot)$.

####  V-A2 Measurement evaluation

We first describe the messages $\vartheta_{i,j}^{(p)}(\cdot)$, with
$i\in\\{1,\dots,n_{k|k-1}\\}$, $j\in\\{1,\dots,m_{k}\\}$ and
$i\in\\{n_{k|k-1}+1,\dots,n_{k|k}\\}$, $j\in\\{1,\dots,i-n_{k|k-1}-1\\}$, sent
from factor nodes $\underline{s}_{k}(\widetilde{{\bf
X}}_{k}^{i},\beta_{k}^{j};z_{k}^{j})$ to variable nodes $\beta_{k}^{j}$.
Applying the sum-product rule [63, Eq. (6)], we obtain an integral over the
Bernoulli set $\widetilde{{\bf X}}_{k}^{i}$

$\Theta_{i,j}^{(p)}\left(\beta_{k}^{j}\right)=\int\underline{s}_{k}\left(\widetilde{{\bf
X}}_{k}^{i},\beta_{k}^{j};z_{k}^{j}\right)\varepsilon^{(p)}_{i,j}\left(\widetilde{{\bf
X}}_{k}^{i}\right)\delta\widetilde{{\bf X}}_{k}^{i}.$ (34)

We then plug (30) into (34), which yields

$\Theta^{(p)}_{i,j}\left(\beta_{k}^{j}\right)=\begin{cases}\frac{\int\ell_{k}\left(z_{k}^{j}|X\right)\gamma_{k}(X)\varepsilon^{(p)}_{i,j}\left(\\{X\\}\right)dX}{\lambda_{k}^{C}\left(z_{k}^{j}\right)}&\beta_{k}^{j}=i\\\
\varepsilon^{(p)}_{i,j}&\beta_{k}^{j}\neq i.\end{cases}$ (35)

After normalizing (35) by $\varepsilon^{(p)}_{i,j}$, we obtain

$\vartheta^{(p)}_{i,j}\left(\beta_{k}^{j}\right)=\begin{cases}\frac{\int\ell_{k}\left(z_{k}^{j}|X\right)\gamma_{k}(X)\varepsilon^{(p)}_{i,j}\left(\\{X\\}\right)dX}{\lambda_{k}^{C}\left(z_{k}^{j}\right)\varepsilon^{(p)}_{i,j}}&\beta_{k}^{j}=i\\\
1&\beta_{k}^{j}\neq i.\end{cases}$ (36)

Note that multiplying messages by a constant does not alter the resulting
approximate marginals [63], but the above normalization step makes it possible
to perform data association and measurement update more efficiently.

The messages $\vartheta_{i,j}^{(p)}(\cdot)$, with
$i\in\\{n_{k|k-1}+1,\dots,n_{k|k}\\}$, $j=i-n_{k|k-1}$, sent from factor nodes
$\overline{s}_{k}(\widetilde{{\bf X}}_{k}^{i},\beta_{k}^{j};z_{k}^{j})$ to
variable nodes $\beta_{k}^{j}$ can be computed similarly. In particular, we
have

$\vartheta^{(p)}_{i,j}\left(\beta_{k}^{j}\right)=\begin{cases}\frac{\int\ell_{k}\left(z_{k}^{j}|X\right)\gamma_{k}(X)\varepsilon^{(p)}_{i,j}\left(\\{X\\}\right)dX}{\lambda_{k}^{C}\left(z_{k}^{j}\right)\varepsilon^{(p)}_{i,j}(\emptyset)}&\beta_{k}^{j}=i\\\
1&\beta_{k}^{j}\neq i.\end{cases}$ (37)

####  V-A3 Data association

We proceed to describe the messages $\varphi_{j,i}^{(p)}(\cdot)$, with
$i\in\\{1,\dots,n_{k|k-1}\\}$, $j\in\\{1,\dots,m_{k}\\}$ and
$i\in\\{n_{k|k-1}+1,\dots,n_{k|k}\\}$, $j\in\\{1,\dots,i-n_{k|k-1}\\}$, sent
from variable nodes $\beta_{k}^{j}$ to factor nodes
$\underline{s}_{k}(\widetilde{{\bf X}}_{k}^{i},\beta_{k}^{j};z_{k}^{j})$ and
$\overline{s}_{k}(\widetilde{{\bf X}}_{k}^{i},\beta_{k}^{j};z_{k}^{j})$, which
can be expressed as [63, Eq. (5)]

$\varphi_{j,i}^{(p)}\left(\beta_{k}^{j}\right)=\prod_{\begin{subarray}{c}i^{\prime}=1\\\
i^{\prime}\neq
i\end{subarray}}^{n_{k|k}}\vartheta^{(p)}_{i^{\prime},j}\left(\beta_{k}^{j}\right).$
(38)

By plugging (36) and (37) into (38), we obtain

$\varphi_{j,i}^{(p)}\left(\beta_{k}^{j}\right)=\begin{cases}1&\beta_{k}^{j}\in\\{0,i\\}\\\
\vartheta^{(p)}_{i,j}\left(i\right)&\beta_{k}^{j}\in\\{1,\dots,n_{k|k}\\}\setminus\\{i\\}.\end{cases}$
(39)

####  V-A4 Measurement update

Next, the messages $\varrho_{j,i}^{(p)}(\cdot)$, with
$i\in\\{1,\dots,n_{k|k-1}\\}$, $j\in\\{1,\dots,m_{k}\\}$ and
$i\in\\{n_{k|k-1}+1,\dots,n_{k|k}\\}$, $j\in\\{1,\dots,i-n_{k|k-1}-1\\}$, sent
from factor nodes $\underline{s}_{k}(\widetilde{{\bf
X}}_{k}^{i},\beta_{k}^{j};z_{k}^{j})$ to variable nodes $\widetilde{{\bf
X}}_{k}^{i}$ can be computed as [63, Eq. (6)]

$\Upsilon_{j,i}^{(p)}\left(\widetilde{{\bf
X}}_{k}^{i}\right)=\sum_{\beta_{k}^{j}}\underline{s}_{k}\left(\widetilde{{\bf
X}}_{k}^{i},\beta_{k}^{j};z_{k}^{j}\right)\varphi_{j,i}^{(p)}\left(\beta_{k}^{j}\right).$
(40)

By plugging (30) and (39) into (40), we obtain

$\displaystyle\Upsilon_{j,i}^{(p)}\left(\widetilde{{\bf X}}_{k}^{i}\right)$
$\displaystyle=\begin{cases}\left(\frac{\ell_{k}\left(z_{k}^{j}|X\right)\gamma_{k}(X)}{\lambda_{k}^{C}\left(z_{k}^{j}\right)}+\xi_{i,j}^{(p)}\right)\delta_{i}[u]&\widetilde{{\bf
X}}_{k}^{i}=\\{(u,X)\\}\\\ \xi_{i,j}^{(p)}&\widetilde{{\bf
X}}_{k}^{i}=\emptyset\end{cases}$ (41)

where we denote

$\xi_{i,j}^{(p)}=\sum_{\begin{subarray}{c}i^{\prime}=1\\\ i^{\prime}\neq
i\end{subarray}}^{n_{k|k}}\vartheta^{(p)}_{i^{\prime},j}\left(i^{\prime}\right)+1.$
(42)

After normalizing (41) by $\xi_{i,j}^{(p)}$, we obtain

$\displaystyle\varrho_{j,i}^{(p)}\left(\widetilde{{\bf X}}_{k}^{i}\right)$
$\displaystyle=\begin{cases}\left(\frac{\ell_{k}\left(z_{k}^{j}|X\right)\gamma_{k}(X)}{\lambda_{k}^{C}\left(z_{k}^{j}\right)\xi_{i,j}^{(p)}}+1\right)\delta_{i}[u]&\widetilde{{\bf
X}}_{k}^{i}=\\{(u,X)\\}\\\ 1&\widetilde{{\bf
X}}_{k}^{i}=\emptyset.\end{cases}$ (43)

The messages $\varrho_{j,i}^{(p)}(\cdot)$, with
$i\in\\{n_{k|k-1}+1,\dots,n_{k|k}\\}$, $j=i-n_{k|k-1}$, sent from factor nodes
$\overline{s}_{k}(\widetilde{{\bf X}}_{k}^{i},\beta_{k}^{j};z_{k}^{j})$ to
variable nodes $\widetilde{{\bf X}}_{k}^{i}$ can be computed similarly. In
particular, we have

$\varrho_{j,i}^{(p)}\left(\widetilde{{\bf
X}}_{k}^{i}\right)=\begin{cases}\frac{\ell_{k}\left(z_{k}^{j}|X\right)\gamma_{k}(X)}{\lambda_{k}^{C}\left(z_{k}^{j}\right)\xi_{i,j}^{(p)}}\delta_{i}[u]&\widetilde{{\bf
X}}_{k}^{i}=\\{(u,X)\\}\\\ 1&\widetilde{{\bf
X}}_{k}^{i}=\emptyset.\end{cases}$ (44)

####  V-A5 Extrinsic information

At last, the messages $\varepsilon_{i,j}^{(p+1)}(\cdot)$ at iteration $p+1$,
sent from variable nodes $\widetilde{{\bf X}}_{k}^{i}$ to factor nodes
$\underline{s}_{k}(\widetilde{{\bf X}}_{k}^{i},\beta_{k}^{j};z_{k}^{j})$, with
$i\in\\{1,\dots,n_{k|k-1}\\}$, $j\in\\{1,\dots,m_{k}\\}$, can be computed as
[63, Eq. (5)]

$\varepsilon_{i,j}^{(p+1)}\left(\widetilde{{\bf
X}}_{k}^{i}\right)=\underline{f}^{i}_{k|k-1}\left(\widetilde{{\bf
X}}_{k}^{i}\right)\prod_{\begin{subarray}{c}j^{\prime}=1\\\ j^{\prime}\neq
j\end{subarray}}^{m_{k}}\varrho_{j^{\prime},i}^{(p)}\left(\widetilde{{\bf
X}}_{k}^{i}\right),$ (45)

and similarly, for the messages $\varepsilon_{i,j}^{(p+1)}(\cdot)$, sent from
variable nodes $\widetilde{{\bf X}}_{k}^{i}$ to factor nodes
$\underline{s}_{k}(\widetilde{{\bf X}}_{k}^{i},\beta_{k}^{j};z_{k}^{j})$ or
$\overline{s}_{k}(\widetilde{{\bf X}}_{k}^{i},\beta_{k}^{j};z_{k}^{j})$, with
$i\in\\{n_{k|k-1}+1,\dots,n_{k|k}\\}$, $j\in\\{1,\dots,i-n_{k|k-1}\\}$, we
have

$\varepsilon_{i,j}^{(p+1)}\left(\widetilde{{\bf
X}}_{k}^{i}\right)=\overline{f}_{k|k-1}\left(\widetilde{{\bf
X}}_{k}^{i}\right)\prod_{\begin{subarray}{c}j^{\prime}=1\\\ j^{\prime}\neq
j\end{subarray}}^{i-n_{k|k-1}}\varrho_{j^{\prime},i}^{(p)}\left(\widetilde{{\bf
X}}_{k}^{i}\right).$ (46)

###  V-B Belief calculation

After the last iteration $p=P$, we can evaluate the belief
$\widetilde{f}^{i}_{k|k}(\cdot)$ for each Bernoulli component
$i\in\\{1,\dots,n_{k|k}\\}$, which is proportional to the product of all
coming messages [63]. In particular, for $i\in\\{1,\dots,n_{k|k-1}\\}$,

$\widetilde{f}^{i}_{k|k}\left(\widetilde{{\bf
X}}_{k}^{i}\right)\propto\underline{f}^{i}_{k|k-1}\left(\widetilde{{\bf
X}}_{k}^{i}\right)\prod_{\begin{subarray}{c}j=1\end{subarray}}^{m_{k}}\varrho_{j,i}^{(P)}\left(\widetilde{{\bf
X}}_{k}^{i}\right),$ (47)

and for $i\in\\{n_{k|k-1}+1,\dots,n_{k|k}\\}$,

$\widetilde{f}^{i}_{k|k}\left(\widetilde{{\bf
X}}_{k}^{i}\right)\propto\overline{f}^{i}_{k|k-1}\left(\widetilde{{\bf
X}}_{k}^{i}\right)\prod_{\begin{subarray}{c}j=1\end{subarray}}^{i-n_{k|k-1}}\varrho_{j,i}^{(P)}\left(\widetilde{{\bf
X}}_{k}^{i}\right)$ (48)

where normalization is required to ensure that
$\widetilde{f}^{i}_{k|k}(\cdot)$ is a valid set density.

##  VI Particle Implementation of Trajectory PMB Filters Using Belief
Propagation

For general multi-object dynamic and measurement models, the messages in BP
typically cannot be evaluated in closed form. In this section, we present the
particle-based implementation333In theory, it is possible to use a Rao-
Blackwellized particle filter for BP, where we only sample kinematic states
while keeping an analytic representation of object extent state conditioned on
its kinematic state. However, because of the way extrinsic information (45),
(46) is computed, doing so requires an analytic representation for every
possible non-empty subset of the measurements, resulting in an exponential
growth of local hypotheses. of the TPMB filter with Poisson spatial
measurement model for both the set of alive trajectories and the set of all
trajectories. In addition, we discuss aspects that need to be considered in
practical implementations and how the proposed implementation compares to the
one in [17]. We note that only the update step of the proposed TPMB filters
involves BP.

We represent a single-trajectory density/intensity by

$f\left(t,x^{1:\nu}\right)=\sum_{l=1}^{L}w^{(l)}\delta_{t^{(l)}}[t]\delta_{\nu^{(l)}}[\nu]\delta_{\chi^{(l)}}\left(x^{1:\nu}\right)$

where the $l$-th particle has weight $w^{(l)}\geq 0$, and it represents a
single trajectory with start time $t^{(l)}$ and sequence of object states
$\chi^{(l)}\in{\cal X}^{\nu^{(l)}}$, where the length $\nu^{(l)}$ is implicit
in $\chi^{(l)}$ [64]. In addition, we denote the last state of $\chi^{(l)}$ as
$\overline{\chi}^{(l)}$. Furthermore, we note that $f(\cdot)$ is a density if
and only if $\sum_{l=1}^{L}w^{(l)}=1$, and that it can be fully described by
the set of parameters $\\{(w^{(l)},t^{(l)},\chi^{(l)})\\}_{l=1}^{L}$.

For the set of alive trajectories, the single-trajectory density of the $i$-th
Bernoulli component is of the form

$f^{i}_{k|k^{\prime}}\left(t,x^{1:\nu}\right)=\sum_{l=1}^{L^{d}_{k|k^{\prime}}}w^{i,(l)}_{k|k^{\prime}}\delta_{t_{k|k^{\prime}}^{i,(l)}}[t]\delta_{k-t_{k|k^{\prime}}^{i,(l)}+1}[\nu]\delta_{\chi_{k|k^{\prime}}^{i,(l)}}\left(x^{1:\nu}\right),$
(49)

which implies that, if the corresponding trajectory exists, it is alive at
time step $k$ with probability one. The PPP for undetected trajectories has
intensity

$\lambda_{k|k^{\prime}}\left(t,x^{1:\nu}\right)=\sum_{l=1}^{L^{u}_{k|k^{\prime}}}w^{0,(l)}_{k|k^{\prime}}\delta_{t_{k|k^{\prime}}^{0,(l)}}[t]\delta_{k-t_{k|k^{\prime}}^{0,(l)}+1}[\nu]\delta_{\chi_{k|k^{\prime}}^{0,(l)}}\left(x^{1:\nu}\right)$
(50)

where $\sum_{l=1}^{L^{u}_{k|k^{\prime}}}w^{0,(l)}_{k|k^{\prime}}$ gives the
expected number of undetected trajectories.

For the set of all trajectories, the single-trajectory density of the $i$-th
Bernoulli component is of the form

$f^{i}_{k|k^{\prime}}\left(t,x^{1:\nu}\right)=\sum_{l=1}^{L^{d}_{k|k^{\prime}}}w^{i,(l)}_{k|k^{\prime}}\delta_{t_{k|k^{\prime}}^{i,(l)}}[t]\delta_{\nu_{k|k^{\prime}}^{i,(l)}}[\nu]\delta_{\chi_{k|k^{\prime}}^{i,(l)}}\left(x^{1:\nu}\right),$
(51)

and the Poisson intensity for undetected trajectories is of the same form as
(50). In the implementation, we do not account for undetected trajectories
that are not present at the current time step since these are usually not of
practical interest.

###  VI-A Prediction step

####  VI-A1 Set of alive trajectories

Assume that the filtering density for the alive trajectories is a PMB of the
form (12) with $f^{i}_{k-1|k-1}(\cdot)$ and $\lambda_{k-1|k-1}(\cdot)$ given
by (49) and (50). Then, the predicted density is a PMB of the form (12) with

$\displaystyle\lambda_{k|k-1}\left(t,x^{1:\nu}\right)$
$\displaystyle=\sum_{l=1}^{L^{u}_{k-1|k-1}}w^{0,(l)}_{k|k-1}\delta_{t_{k|k-1}^{0,(l)}}[t]$
$\displaystyle~{}~{}~{}\times\delta_{k-t_{k|k-1}^{0,(l)}+1}[\nu]\delta_{\chi_{k|k-1}^{0,(l)}}\left(x^{1:\nu}\right)$
$\displaystyle~{}~{}~{}+\sum_{l=1}^{L^{b}}w_{k}^{b,(l)}\delta_{k}[t]\delta_{1}[\nu]\delta_{\chi^{b,(l)}_{k}}\left(x^{1:\nu}\right),$
(52a) $\displaystyle L^{u}_{k|k-1}$ $\displaystyle=L^{u}_{k-1|k-1}+L^{b},$
(52b) $\displaystyle t_{k|k-1}^{0,(l)}$ $\displaystyle=t_{k-1|k-1}^{0,(l)},$
(52c) $\displaystyle x_{k}^{0,(l)}$ $\displaystyle\sim
g_{k}\left(\cdot|\overline{\chi}^{0,(l)}_{k-1|k-1}\right),$ (52d)
$\displaystyle\chi^{0,(l)}_{k|k-1}$
$\displaystyle=\left(\chi^{0,(l)}_{k-1|k-1},x_{k}^{0,(l)}\right),$ (52e)
$\displaystyle w^{0,(l)}_{k|k-1}$
$\displaystyle=p^{S}\left(\overline{\chi}^{0,(l)}_{k-1|k-1}\right)w^{0,(l)}_{k-1|k-1},$
(52f) $\displaystyle w_{k}^{b,(l)}$
$\displaystyle=\left\langle\lambda^{B}_{k},1\right\rangle/L^{b},$ (52g)
$\displaystyle\chi^{b,(l)}_{k}$
$\displaystyle\sim\lambda_{k}^{B}(\cdot)/\left\langle\lambda^{B}_{k},1\right\rangle$
(52h)

where $L^{b}$ is the number of particles used to represent newborn
trajectories, and we can see that the prediction of undetected trajectories
comprises the prediction of existing particles and the generation of new
particles. As for the predicted Bernoulli component, it is parameterized by

$\displaystyle L^{d}_{k|k-1}$ $\displaystyle=L^{d}_{k-1|k-1},$ (53a)
$\displaystyle r^{i}_{k|k-1}$
$\displaystyle=r^{i}_{k-1|k-1}\sum_{l=1}^{L^{d}_{k|k-1}}p^{S}\left(\overline{\chi}_{k-1|k-1}^{i,(l)}\right)w^{i,(l)}_{k-1|k-1},$
(53b) $\displaystyle t_{k|k-1}^{i,(l)}$ $\displaystyle=t_{k-1|k-1}^{i,(l)},$
(53c) $\displaystyle x_{k}^{i,(l)}$ $\displaystyle\sim
g_{k}\left(\cdot|\overline{\chi}^{i,(l)}_{k-1|k-1}\right),$ (53d)
$\displaystyle\chi_{k|k-1}^{i,(l)}$
$\displaystyle=\left(\chi_{k-1|k-1}^{i,(l)},x_{k}^{i,(l)}\right),$ (53e)
$\displaystyle w^{i,(l)}_{k|k-1}$
$\displaystyle=\frac{p^{S}\left(\overline{\chi}_{k-1|k-1}^{i,(l)}\right)w^{i,(l)}_{k-1|k-1}}{\sum_{l=1}^{L^{d}_{k|k-1}}p^{S}\left(\overline{\chi}_{k-1|k-1}^{i,(l)}\right)w^{i,(l)}_{k-1|k-1}}.$
(53f)

The proposal densities used in the above prediction step for drawing samples
are given by the motion model $g_{k}(\cdot|\cdot)$ and the birth density
$\lambda_{k}^{B}(\cdot)/\left\langle\lambda^{B}_{k},1\right\rangle$. If these
densities are difficult to sample from, then other proposal densities may be
used for importance sampling [65].

####  VI-A2 Set of all trajectories

Assume that the filtering density for the set of all trajectories is a PMB of
the form (12) with $f^{i}_{k-1|k-1}(\cdot)$ and $\lambda_{k-1|k-1}(\cdot)$
given by (51) and (50). Then, the predicted density is a PMB of the form (12)
with $\lambda_{k|k-1}(\cdot)$ given by (52) and $f^{i}_{k|k-1}(\cdot)$
parameterized by

$\displaystyle r^{i}_{k|k-1}$ $\displaystyle=r^{i}_{k-1|k-1},$ (54a)
$\displaystyle L^{d}_{k|k-1}$ $\displaystyle=L^{d}_{k-1|k-1}$
$\displaystyle~{}~{}~{}+\sum_{l=1}^{L^{d}_{k-1|k-1}}\delta_{k-t^{i,(l)}_{k-1|k-1}}\left[\nu_{k-1|k-1}^{i,(l)}\right],$
(54b) $\displaystyle f^{i}_{k|k-1}\left(t,x^{1:\nu}\right)$
$\displaystyle=\sum_{l=1}^{L^{d}_{k-1|k-1}}w^{i,(l)}_{k|k-1}\delta_{t_{k|k-1}^{i,(l)}}[t]$
$\displaystyle~{}~{}~{}\times\delta_{\nu_{k|k-1}^{i,(l)}}[\nu]\delta_{\chi_{k|k-1}^{i,(l)}}\left(x^{1:\nu}\right)$
$\displaystyle~{}~{}~{}+\sum_{l=1}^{L^{d}_{k-1|k-1}}\left[1-p^{S}\left(\overline{\chi}^{i,(l)}_{k-1|k-1}\right)\right]$
$\displaystyle~{}~{}~{}\times
w^{i,(l)}_{k-1|k-1}\delta_{t_{k-1|k-1}^{i,(l)}}[t]$
$\displaystyle~{}~{}~{}\times\delta_{k-t^{i,(l)}_{k-1|k-1}}\left[\nu_{k-1|k-1}^{i,(l)}\right]\delta_{\chi_{k-1|k-1}^{i,(l)}}\left(x^{1:\nu}\right),$
(54c) $\displaystyle t_{k|k-1}^{i,(l)}$ $\displaystyle=t_{k-1|k-1}^{i,(l)},$
(54d) and if $\nu_{k-1|k-1}^{i,(l)}=k-t^{i,(l)}_{k-1|k-1}$, $\displaystyle
x_{k}^{i,(l)}$ $\displaystyle\sim
g_{k}\left(\cdot|\overline{\chi}^{i,(l)}_{k-1|k-1}\right),$ (54e)
$\displaystyle\chi_{k|k-1}^{i,(l)}$
$\displaystyle=\left(\chi_{k-1|k-1}^{i,(l)},x_{k}^{i,(l)}\right),$ (54f)
$\displaystyle w^{i,(l)}_{k|k-1}$
$\displaystyle=p^{S}\left(\overline{\chi}_{k-1|k-1}^{i,(l)}\right)w^{i,(l)}_{k-1|k-1},$
(54g) otherwise $\displaystyle\chi_{k|k-1}^{i,(l)}$
$\displaystyle=\chi_{k-1|k-1}^{i,(l)},$ (54h) $\displaystyle
w^{i,(l)}_{k|k-1}$ $\displaystyle=w^{i,(l)}_{k-1|k-1}.$ (54i)

We note that for each particle used to represent a Bernoulli component, it
remains unchanged if it corresponds to a dead trajectory, and that a copy of
it is created if it corresponds to an alive trajectory. In the latter case,
there is a change in the weight for original particles to account for the
probability that the object dies, and for each copy, its weight is updated
with the survival probability and its state is propagated to the next time
step. The dimensions of particles may be high for long trajectories, but this
does not mean that the computation would be infeasible. In the standard
particle filter implementations, the particle trajectory length also increases
in time, see, e.g., [66, Eq. (40)]. The particle degeneracy problem in
practical implementations is dealt in Section VI-C3. We also note that the
particle implementation of the TPMB prediction step is general, which also
holds for TPMBM and is not limited to EOT.

###  VI-B Update step

The update step is the same for both the set of alive trajectories and the set
of all trajectories, where for each Bernoulli component we only update its
probability of existence and the weights of particles while leaving the states
of particles unchanged. Assume that the predicted density is a PMB of the form
(12) with parameters described in either Section VI-A1 or Section VI-A2, the
updated Poisson intensity $\lambda_{k|k}(\cdot)$ for undetected trajectories
is described by

$\left\\{\left(t^{0,(l)}_{k|k-1},\chi^{0,(l)}_{k|k-1}w_{k|k-1}^{0,(l)}e^{-\gamma_{k}\left(\overline{\chi}_{k|k-1}^{0,(l)}\right)}\right)\right\\}_{l=1}^{L^{u}_{k|k-1}},$

and the implementation of particle BP is given as follows, where we omit the
auxiliary variables in all the messages and beliefs for notational brevity.

####  VI-B1 Initialization

Each message $\varepsilon_{i,j}^{(p)}(\cdot)$ is represented by a scalar and a
set of weighted particles

$\left(\varepsilon_{i,j}^{(p)},\left\\{\left(t_{k|k}^{i,(l)},\chi_{k|k}^{i,(l)},w_{i,j}^{(p,l)}\right)\right\\}_{l=1}^{L_{i}}\right).$

For $p=1$ and $i\in\\{1,\dots,n_{k|k-1}\\}$, we set

$\displaystyle\left\\{\left(t_{k|k}^{i,(l)},\chi_{k|k}^{i,(l)}\right)\right\\}_{l=1}^{L_{i}}=\left\\{\left(t^{i,(l)}_{k|k-1},\chi_{k|k-1}^{i,(l)}\right)\right\\}_{l=1}^{L^{d}_{k|k-1}},$
$\displaystyle w_{i,j}^{(1,l)}$
$\displaystyle=\begin{cases}w^{i,(l)}_{k|k-1}e^{-\gamma_{k}\left(\overline{\chi}_{k|k-1}^{i,(l)}\right)}&k-t^{i,(l)}_{k|k-1}+1=\nu^{i,(l)}_{k|k-1}\\\
w^{i,(l)}_{k|k-1}&\text{otherwise},\end{cases}$
$\displaystyle\varepsilon_{i,j}^{(1)}$
$\displaystyle=r^{i}_{k|k-1}\sum_{l=1}^{L^{d}_{k|k-1}}\delta_{k-t^{i,(l)}_{k|k-1}+1}\left[\nu_{k|k-1}^{i,(l)}\right]w_{k|k-1}^{i,(l)}$
$\displaystyle~{}~{}~{}\times
e^{-\gamma_{k}\left(\overline{\chi}_{k|k-1}^{i,(l)}\right)}+1-r^{i}_{k|k-1},$

and for $i\in\\{n_{k|k-1}+1,\dots,n_{k|k}\\}$, we set

$\displaystyle\left\\{\left(t_{k|k}^{i,(l)},\chi_{k|k}^{i,(l)},w_{i,j}^{(1,l)}\right)\right\\}_{l=1}^{L_{i}}$
$\displaystyle=\left\\{\left(t^{0,(l)}_{k|k-1},\chi_{k|k-1}^{0,(l)},w_{k|k-1}^{0,(l)}e^{-\gamma_{k}\left(\overline{\chi}_{k|k-1}^{0,(l)}\right)}\right)\right\\}_{l=1}^{L^{u}_{k|k-1}},$
$\displaystyle\varepsilon_{i,j}^{(1)}$
$\displaystyle=\sum_{l=1}^{L^{u}_{k|k-1}}w^{0,(l)}_{k|k-1}e^{-\gamma_{k}\left(\overline{\chi}_{k|k-1}^{0,(l)}\right)}+1.$

####  VI-B2 Measurement evaluation

For $i\in\\{1,\dots,n_{k|k-1}\\}$, $j\in\\{1,\dots,m_{k}\\}$ and
$i\in\\{n_{k|k-1}+1,\dots,n_{k|k}\\}$, $j\in\\{1,\dots,i-n_{k|k-1}-1\\}$, we
have

$\displaystyle\vartheta^{(p)}_{i,j}\left(i\right)$
$\displaystyle=\frac{1}{\lambda^{C}_{k}\left(z_{k}^{j}\right)\varepsilon_{i,j}^{(p)}}\sum_{l=1}^{L_{i}}w_{i,j}^{(p,l)}\delta_{k-t^{i,(l)}_{k|k}+1}\left[\nu_{k|k}^{i,(l)}\right]$
$\displaystyle~{}~{}~{}\times\gamma_{k}\left(\overline{\chi}_{k|k}^{i,(l)}\right)\ell_{k}\left(z_{k}^{j}|\overline{\chi}_{k|k}^{i,(l)}\right)$
(55)

As for $i\in\\{n_{k|k-1}+1,\dots,n_{k|k}\\}$, $j=i-n_{k|k-1}$, the message
$\vartheta^{(p)}_{i,j}(i)$ is obtained by replacing $\varepsilon_{i,j}^{(p)}$
with $\varepsilon_{i,j}^{(p)}(\emptyset)$ in (VI-B2). It should be noted that
in (VI-B2) we only consider particles that represent alive trajectories.

####  VI-B3 Data association, measurement update, and extrinsic information

The data association and measurement update steps do not need to be explicitly
implemented. After computing $\vartheta^{(p)}_{i,j}(\cdot)$, we can then
obtain $\xi_{j,i}^{(p)}(\cdot)$ using (42) and (VI-B2). Each message
$\varepsilon_{i,j}^{(p+1)}(\cdot)$ is an unnormalized Bernoulli density,
represented by a scalar and a set of weighted particles

$\left(\varepsilon_{i,j}^{(p+1)},\left\\{\left(t_{k|k}^{i,(l)},\chi_{k|k}^{i,(l)},w_{i,j}^{(p+1,l)}\right)\right\\}_{l=1}^{L_{i}}\right),$

which can be calculated by first computing the particle-based representation
of (43)/(44) and then plugging it into (45)/(46). This yields, for
$i\in\\{1,\dots,n_{k|k-1}\\}$, $j\in\\{1,\dots,m_{k}\\}$,

$\displaystyle\varepsilon_{i,j}^{(p+1)}$
$\displaystyle=r^{i}_{k|k-1}\sum_{l=1}^{L_{i}}w_{i,j}^{(p+1,l)}+1-r^{i}_{k|k-1},$
(56) $\displaystyle w_{i,j}^{(p+1,l)}$
$\displaystyle=w_{i,j}^{(1,l)}\prod_{\begin{subarray}{c}j^{\prime}=1\\\
j^{\prime}\neq
j\end{subarray}}^{m_{k}}\Bigg{[}\frac{\ell_{k}\left(z_{k}^{j^{\prime}}|\overline{\chi}_{k|k}^{i,(l)}\right)\gamma_{k}\left(\overline{\chi}_{k|k}^{i,(l)}\right)}{\lambda_{k}^{C}\left(z_{k}^{j^{\prime}}\right)\xi_{i,j^{\prime}}^{(p)}}$
$\displaystyle~{}~{}~{}\times\delta_{k-t^{i,(l)}_{k|k}+1}\left[\nu_{k|k}^{i,(l)}\right]+1\Bigg{]},$
(57)

for $i\in\\{n_{k|k-1}+1,\dots,n_{k|k}\\}$, $j\in\\{1,\dots,i-n_{k|k-1}\\}$,

$\varepsilon_{i,j}^{(p+1)}=\sum_{l=1}^{L_{i}}w_{i,j}^{(p+1,l)}+1,$ (58)

for $i\in\\{n_{k|k-1}+1,\dots,n_{k|k}\\}$, $j\in\\{1,\dots,i-n_{k|k-1}-1\\}$,

$\displaystyle w_{i,j}^{(p+1,l)}$
$\displaystyle=w_{i,j}^{(1,l)}\frac{\ell_{k}\left(z_{k}^{i-n_{k|k-1}}|\overline{\chi}_{k|k}^{i,(l)}\right)\gamma_{k}\left(\overline{\chi}_{k|k}^{i,(l)}\right)}{\lambda_{k}^{C}\left(z_{k}^{i-n_{k|k-1}}\right)\xi_{i,i-n_{k|k-1}}^{(p)}}$
$\displaystyle~{}~{}~{}\times\prod_{\begin{subarray}{c}j^{\prime}=1\\\
j^{\prime}\neq
j\end{subarray}}^{i-n_{k|k-1}-1}\left[\frac{\ell_{k}\left(z_{k}^{j^{\prime}}|\overline{\chi}_{k|k}^{i,(l)}\right)\gamma_{k}\left(\overline{\chi}_{k|k}^{i,(l)}\right)}{\lambda_{k}^{C}\left(z_{k}^{j^{\prime}}\right)\xi_{i,j^{\prime}}^{(p)}}+1\right],$
(59)

and for $i\in\\{n_{k|k-1}+1,\dots,n_{k|k}\\}$, $j=i-n_{k|k-1}$,

$w_{i,j}^{(p+1,l)}=w_{i,j}^{(1,l)}\prod_{\begin{subarray}{c}j^{\prime}=1\end{subarray}}^{j-1}\left[\frac{\ell_{k}\left(z_{k}^{j^{\prime}}|\overline{\chi}_{k|k}^{i,(l)}\right)\gamma_{k}\left(\overline{\chi}_{k|k}^{i,(l)}\right)}{\lambda_{k}^{C}\left(z_{k}^{j^{\prime}}\right)\xi_{i,j^{\prime}}^{(p)}}+1\right].$
(60)

####  VI-B4 Belief calculation

Each belief $\widetilde{f}^{i}_{k|k}(\cdot)$ is a Bernoulli RFS density,
parameterized by

$\left(r^{i}_{k|k},\left\\{\left(t^{i,(l)}_{k|k},\chi^{i,(l)}_{k|k},w^{i,(l)}_{k|k}\right)\right\\}_{l=1}^{L_{i}}\right),$

which can be calculated similarly as the extrinsic information
$\varepsilon_{i,j}^{(p+1)}(\cdot)$ with the difference that the product in
(57) and (59) now needs to enumerate every element, and that normalization is
required to ensure that $\widetilde{f}^{i}_{k|k}(\cdot)$ is a valid density.

Specifically, we have, for $i\in\\{1,\dots,n_{k|k-1}\\}$,

$\displaystyle w_{i}^{(P,l)}$
$\displaystyle=w_{i,1}^{(1,l)}\prod_{\begin{subarray}{c}j^{\prime}=1\end{subarray}}^{m_{k}}\Bigg{[}\frac{\ell_{k}\left(z_{k}^{j^{\prime}}|\overline{\chi}_{k|k}^{i,(l)}\right)\gamma_{k}\left(\overline{\chi}_{k|k}^{i,(l)}\right)}{\lambda_{k}^{C}\left(z_{k}^{j^{\prime}}\right)\xi_{i,j^{\prime}}^{(P)}}$
$\displaystyle~{}~{}~{}\times\delta_{k-t^{i,(l)}_{k|k}+1}\left[\nu_{k|k}^{i,(l)}\right]+1\Bigg{]},$
(61) $\displaystyle r^{i}_{k|k}$
$\displaystyle=\frac{r^{i}_{k|k-1}\sum_{l=1}^{L_{i}}w_{i}^{(P,l)}}{r^{i}_{k|k-1}\sum_{l=1}^{L_{i}}w_{i}^{(P,l)}+1-r^{i}_{k|k-1}},$
(62)

and for $i\in\\{n_{k|k-1}+1,\dots,n_{k|k}\\}$,

$\displaystyle w_{i}^{(P,l)}$
$\displaystyle=w_{i,1}^{(1,l)}\frac{\ell_{k}\left(z_{k}^{i-n_{k|k-1}}|\overline{\chi}_{k|k}^{i,(l)}\right)\gamma_{k}\left(\overline{\chi}_{k|k}^{i,(l)}\right)}{\lambda_{k}^{C}\left(z_{k}^{i-n_{k|k-1}}\right)\xi_{i,i-n_{k|k-1}}^{(P)}}$
$\displaystyle~{}~{}~{}\times\prod_{\begin{subarray}{c}j^{\prime}=1\end{subarray}}^{i-n_{k|k-1}-1}\left[\frac{\ell_{k}\left(z_{k}^{j^{\prime}}|\overline{\chi}_{k|k}^{i,(l)}\right)\gamma_{k}\left(\overline{\chi}_{k|k}^{i,(l)}\right)}{\lambda_{k}^{C}\left(z_{k}^{j^{\prime}}\right)\xi_{i,j^{\prime}}^{(P)}}+1\right],$
(63) $\displaystyle r^{i}_{k|k}$
$\displaystyle=\frac{\sum_{l=1}^{L_{i}}w_{i}^{(P,l)}}{\sum_{l=1}^{L_{i}}w_{i}^{(P,l)}+1},$
(64)

and the normalized weight is given by

$w^{i,(l)}_{k|k}=\frac{w_{i}^{(P,l)}}{\sum_{l=1}^{L_{i}}w_{i}^{(P,l)}}.$ (65)

###  VI-C Practical considerations

####  VI-C1 Approximations for efficient implementation

For the proposed particle implementation of TPMB, the computational complexity
of its prediction step scales linearly in the number of Bernoulli components
(potential trajectories). As for the update step, its computational complexity
for a fixed number of message passing iterations scales quadratically in the
number of Bernoulli components and measurements [17].

Since the number of Bernoulli components increases with time, we need to prune
Bernoulli components with probability of existence smaller than a threshold.
For the implementation considering the set of all trajectories, the number of
particles in the Bernoulli components quickly increases over time as more
particles are used to represent dead trajectories in the prediction step. To
avoid this, for the $i$-th Bernoulli component, $i\in\\{1,\dots,n_{k|k-1}\\}$,
we first compute the probability mass function of the trajectory end time
$t^{i}_{e}$

$P\left(t^{i}_{e}=t\right)=\sum_{l=1}^{L^{d}_{k|k}}w^{i,(l)}_{k|k}\delta_{t}\left[t^{i,(l)}_{k|k}+\nu_{k|k}^{i,(l)}-1\right],$
(66)

and find those $t^{i}_{e}$ with probability (66) smaller than a threshold.
Then, we discard particles with the corresponding trajectory end times and re-
normalize the particle weights and the probability of existence $r^{i}_{k|k}$.
At last, we note that computational complexity can be further reduced by
censoring of messages and reordering of measurements, see [16] for details.

####  VI-C2 Measurement-driven initialization of newly detected trajectories

In the initialization step of particle BP described in Section VI-B1,
particles describing newly detected trajectories are set as particles
representing undetected trajectories. For non-informative birth densities,
e.g., a single Gaussian with very large covariance and uniform distribution,
it becomes more advantageous to use importance sampling and directly draw
samples from a proposal density related to the measurement likelihood
$\ell_{k}(\cdot|\cdot)$444If $\ell_{k}(z_{k}|x_{k})$ does not depend on all
the elements of $x_{k}$, one can draw samples of elements of $x_{k}$ that do
not depend on $z_{k}$ from a suitable proposal density, e.g., a Gaussian or a
uniform distribution. [48].

Specifically, in the initialization step, for
$i\in\\{n_{k|k-1}+1,\dots,n_{k|k}\\}$ and $l\in\\{1,\dots,L^{u}_{k|k-1}\\}$,
we draw sample $\chi_{k|k}^{i,(l)}$ from a proposal density
$\ell^{\prime}_{k}(\cdot|z_{i-n_{k|k-1}})$, and its corresponding weight
$w_{i,j}^{(1,l)}$ is given by

$w_{i,j}^{(1,l)}\propto\frac{\sum_{l=1}^{L^{u}_{k|k-1}}w^{0,(l)}_{k|k-1}e^{-\gamma_{k}\left(\overline{\chi}_{k|k-1}^{0,(l)}\right)}K(\chi_{k|k}^{i,(l)}-\overline{\chi}_{k|k-1}^{0,(l)})}{\ell^{\prime}_{k}(\chi_{k|k}^{i,(l)}|z_{i-n_{k|k-1}})}$

where $K(\cdot)$ is a user-defined kernel function [48], typically a
multivariate Gaussian with suitable covariance. At last, the weight
$w_{i,j}^{(1,l)}$ needs to be normalized such that its sum is given by
$\sum_{l=1}^{L^{u}_{k|k-1}}w^{0,(l)}_{k|k-1}e^{-\gamma_{k}\left(\overline{\chi}_{k|k-1}^{0,(l)}\right)}$.

In practical scenarios, the birth density is typically uniform or a Gaussian
mixture. For these cases with state-independent Poisson measurement rate
$\gamma_{k}$, it is not necessary to consider particle representation of the
Poisson intensity of undetected trajectories $\lambda_{k|k^{\prime}}(\cdot)$,
and consequently there is no need to draw particles in the prediction step for
undetected trajectories. In particular, when the birth density is uniform,
$\lambda_{k|k^{\prime}}(\cdot)$ can be approximately represented using only a
scalar [48], whereas when the birth density is a Gaussian mixture,
$\lambda_{k|k^{\prime}}(\cdot)$ can be computed using the Gaussian
implementation in [52]. For the more general cases, one can also use a grid of
points to model the Poisson intensity of undetected objects [67].

####  VI-C3 Trajectory estimation

One problem with sequential importance sampling is weight degeneracy, which
can be reduced by resampling particles that represent alive trajectories [65].
Furthermore, resampling can be used to cap the number of particles in single-
trajectory density/intensity representation. We also note that, by keeping the
full sample histories in the prediction step, particle filtering provides an
approximation to the smoothing problem as a by-product. However, the resulting
approximation also tends to be degenerate for alive trajectories at some point
in the past as many particles may share the same history at earlier time steps
[68]. We proceed to describe how to obtain reasonable trajectory estimates in
the proposed TPMB implementation.

We first select Bernoulli components with existence probability above a
threshold, and then from each of them, we obtain a single-trajectory estimate
from a set of weighted particles (i.e., a mixture of Dirac delta single-
trajectory densities) as follows. At time step $k^{\prime}$, given the
particle representation of single-trajectory density of the $i$-th Bernoulli
component,

$\left\\{\left(w_{k^{\prime}|k^{\prime}}^{i,(l)},t_{k^{\prime}|k^{\prime}}^{i,(l)},\chi_{k^{\prime}|k^{\prime}}^{i,(l)}\right)\right\\}_{l=1}^{L^{d}_{k^{\prime}|k^{\prime}}},$

we extract and store the particle representation of the marginal single-object
state at time step $k^{\prime}$,

$\displaystyle\left\\{\left(w_{k^{\prime}|k^{\prime}}^{i,(l)},\overline{\chi}_{k^{\prime}|k^{\prime}}^{i,(l)}\right):l\in\mathbb{L}^{i}_{k^{\prime}}\right\\},$
(67a)
$\displaystyle\mathbb{L}^{i}_{k^{\prime}}=\left\\{l:k^{\prime}-t_{k^{\prime}|k^{\prime}}^{i,(l)}+1=\nu_{k^{\prime}|k^{\prime}}^{i,(l)}\right\\}$
(67b)

where $\nu_{k^{\prime}|k^{\prime}}^{i,(l)}$ is the length of
${\chi}_{k^{\prime}|k^{\prime}}^{i,(l)}$. To extract the trajectory estimate
at time step $k^{\prime}\leq k$, we first compute the probability mass
functions of the trajectory start time $t^{i}_{s}$ and end time $t^{i}_{e}$,
which are given by

$P\left(t^{i}_{s}=t\right)=\sum_{l=1}^{L^{d}_{k|k}}w^{i,(l)}_{k|k}\delta_{t}\left[t^{i,(l)}_{k|k}\right],$
(68)

and (66), respectively. Note that for the implementation considering the set
of current trajectories, we have $P(t^{i}_{e}=k)=1$.

Next, we find the maximum a posteriori estimates $\hat{t}^{i}_{s}$ and
$\hat{t}^{i}_{e}$ of the trajectory start and end times. Then, we can obtain
the object state estimate at time step $k^{\prime}$ with $\hat{t}^{i}_{s}\leq
k^{\prime}\leq\hat{t}^{i}_{e}$ by

$\hat{x}^{i}_{k^{\prime}}=\frac{1}{\sum_{l\in\mathbb{L}_{k^{\prime}}^{i}}w^{i,(l)}_{k^{\prime}|k^{\prime}}}\sum_{l\in{\mathbb{L}_{k^{\prime}}^{i}}}w^{i,(l)}_{k^{\prime}|k^{\prime}}\overline{\chi}_{k^{\prime}|k^{\prime}}^{i,(l)},$
(69)

where
$\\{(w_{k^{\prime}|k^{\prime}}^{i,(l)},\overline{\chi}_{k^{\prime}|k^{\prime}}^{i,(l)}):l\in\mathbb{L}^{i}_{k^{\prime}}\\}$
are pre-stored particles (67) at previous time steps. This can be understood
as using a type of $L$-scan approximation [52] with $L=1$ for the object
states. In the general $L$-scan approximation [55], past states of the
trajectories before the last $L$ time steps are considered independent, and no
smoothing is performed when $L=1$. For Gaussian implementations, a larger $L$
yields more smoothed trajectory estimates. As for particle implementations,
using a large $L$ does not necessarily improve the trajectory estimates due to
the particle history degeneracy problem.

At last, we note that, if desired, individual smoothed trajectory estimate can
be obtained using, e.g., backward simulation [51], which is described in
Appendix D. Backward simulation can be applied either only at the final time
step or in a sliding window in the form of fixed-lag smoothing. Doing so does
not improve the accuracy of the data association results or the object
cardinality estimates. It is also possible to consider joint backward
simulation of the whole multi-trajectory distribution [64], where improved
data associations can be obtained during backward smoothing.

###  VI-D Discussion

The connection between the proposed TPMB filter and the SPA filter [17] can be
understood as follows. The SPA filter may be seen as a TPMB filter for the set
of alive trajectories, where previous object states are marginalized out and
the PPP intensity of undetected objects is set to zero after the update step.
From another perspective, the TPMB filter may be regarded as a derivation of
the SPA filter using RFSs, with undetected objects/trajectories propagated
over time in parallel using the standard PMBM equations for the PPP, which are
equivalent to a zero-measurement PHD filter [69].

There are several important differences between the proposed TPMB filters and
the SPA filter in [17]. First, the set of undetected objects is explicitly
modelled as a PPP in the TPMB filters, and its Poisson intensity is propagated
over time. This becomes advantageous when the state-dependent Poisson
measurement rate $\gamma_{k}(\cdot)$ is small, which could happen in, e.g.,
scenarios with low sensor resolution, occlusion, or long distance between the
sensor and the objects. One such application is EOT with long-range automotive
radar. Second, the object survival probability $p^{S}(\cdot)$ is state-
dependent in the TPMB filters.

Lastly, the TPMB filter directly reports trajectory estimates from local
trajectory hypotheses (trajectory Bernoulli components), and therefore these
trajectory estimates do not have any gaps in between. As a comparison, the SPA
filter in [17] only estimates the set of object states at the current time
step, and post-processing is required to link object state estimates at
different time steps and to bridge possible gaps. In addition, although for
both implementations smoothed object state estimates can be obtained using
particle smoothing techniques, the particle TPMB filters also provide smoothed
estimate of the start time and length of trajectories, in the sense that
estimates at early time steps may be refined at later time steps. This is an
important feature for MOT algorithms based on sets of trajectories.

##  VII Simulation Results

In this section, we present the results from a Monte Carlo simulation with 200
runs where the performance of the following multiple extended object trackers
are compared555MATLAB implementations of TPMB-BP and PMB-BP are available at
https://github.com/yuhsuansia/Trajectory-PMB-EOT-BP. MATLAB implementation of
SPA is available at https://github.com/meyer-ucsd/EOT-TSP-21. MATLAB
implementations of TPMB(M)-CA and TPMB(M)-SO are available at
https://github.com/yuhsuansia/Extended-target-PMBM-tracker.:

  * •

Trajectory PMB filter using BP, referred to as TPMB-BP.

  * •

PMB filter using BP, referred to as PMB-BP.

  * •

The SPA filter [17].

  * •

Trajectory PMBM and PMB filters using clustering and assignment [34], referred
to as TPMBM-CA and TPMB-CA.

  * •

Trajectory PMBM and PMB filters using stochastic optimization [33], referred
to as TPMBM-SO and TPMB-SO.

For all the TPMBM and TPMB filters, we consider the implementations for the
set of all trajectories. We also present PMB-BP, which can be obtained from
TPMB-BP by marginalizing out all the previous states after the prediction step
[59], and its trajectory estimates are obtained by linking object state
estimates that originate from the same first detection (i.e., the same
Bernoulli component). Note that TPMB-BP and PMB-BP have the same filtering
performance in terms of the set of current object states estimate.

###  VII-A Single object model

There are several extended object models available in the literature, see [3]
for an overview. We consider the random matrix model in [70], in which the
object shape is approximated as an ellipse. The random matrix model has been
used in several PMBM implementations [32, 34, 33, 37], thereby making the
comparison easy.

The single-object state $x_{k}=(e_{k},E_{k})$, represented using a tuple [20,
Eq. (21.88)], consists of a kinematic state vector, which describes the two-
dimensional position and velocity of the object, and an extent state $E_{k}$,
which is a $2\times 2$ symmetric positive definite matrix. We assume that the
object moves according to a nearly constant velocity model and its extent
remains unchanged over time. In this case, the object state transition density
is

$g_{k}(x_{k}|x_{k-1})={\cal N}\left(e_{k};Fe_{k-1}+Q\right){\cal
W}\left(E_{k};E_{k-1}/q,q\right),$ $F=I_{2}\otimes\begin{bmatrix}1&T_{s}\\\
0&1\end{bmatrix},\quad
Q=\sigma_{q}^{2}I_{2}\otimes\begin{bmatrix}T_{s}^{3}/3&T_{s}^{2}/2\\\
T_{s}^{2}/2&T_{s}\end{bmatrix}$

where $I_{2}$ is a $2\times 2$ identity matrix, $\otimes$ denotes the
Kronecker product, $T_{s}$ is the sampling period, $\sigma_{q}$ is standard
acceleration deviation, and ${\cal W}\left(E_{k};E_{k-1}/q,q\right)$
represents a Wishart distribution with mean $E_{k}$ and degree of freedom $q$.
The single measurement likelihood is

$\ell_{k}(z_{k}|x_{k})={\cal N}(z_{k};He_{k},\rho E_{k}+R),$
$H=I_{2}\otimes\begin{bmatrix}1&0\end{bmatrix},\quad R=\sigma_{r}^{2}I_{2}$

where $\rho>0$ is a scaling factor and $\sigma_{r}$ is the measurement noise
deviation.

The gamma distribution is the conjugate prior for the Poisson likelihood,
whereas the Gaussian inverse-Wishart (GIW) distribution is the conjugate prior
for the multivariate Gaussian likelihood with unknown mean and covariance.
Therefore, the single-object state density in TPMB(M)-CA and TPMB(M)-SO is a
gamma GIW (GGIW).

###  VII-B Simulation scenario

We consider the same scenario as in [17] where ten object tracks intersect at
the centre of the region of interest of size
$[-150~{}\text{m},150~{}\text{m}]\times[-150~{}\text{m},150~{}\text{m}]$. The
ten objects start moving towards the centre from positions uniformly placed on
a circle of radius $75~{}\text{m}$ around the centre with initial velocity
$10~{}\text{m/s}$, and then they become closely-spaced for some time before
they separate. The extent of each object is obtained by sampling from an
inverse-Wishart distribution with mean $9I_{2}$ and degree of freedom $1000$.
The object survival probability is $p^{S}=0.99$. The following parameters for
the single object model specified in Section VII-A are used: $T_{s}=0.2$,
$\sigma_{q}=1$, $q=1000$, $\rho=1$ and $\sigma_{r}=1$. The PPP clutter is
uniformly distributed in the region of interest with mean $\gamma^{C}=10$.
Three different settings for the mean number of measurements in the Poisson
spatial model are considered: $\gamma\in\\{3,5,7\\}$, and the corresponding
effective probabilities of detection are approximately
$1-\ell(\emptyset|\cdot)\in\\{0.950,0.993,0.999\\}$.

Figure 3: Example realization of true object trajectories. True object extents
at the times objects appear are also shown. The ten objects appear in pairs at
time step $3,6,9,12,15$ and disappear also in pairs at time step
$83,86,89,92,95$, respectively. The true object extents may (partly) overlap
around the centre.

###  VII-C Implementation details

For all the implementations, the Poisson rate of PPP birth is set to 0.01, and
in the birth density, the velocity is Gaussian distributed with zero mean and
covariance $100I_{2}$, and the extent is inverse-Wishart distributed with mean
$9I_{2}$ and degree of freedom $1000$. For particle-based implementations, the
position in the birth density is uniformly distributed in the region of
interest, whereas for GGIW implementations, the position is Gaussian
distributed with zero mean and covariance $150^{2}I_{2}$, and the Poisson
measurement rate is gamma distributed with shape $1000\gamma$ and scale
$1000$.

For implementations using particle BP, each particle describes both kinematic
and extent states, and therefore more particles are needed to guarantee a
reasonable performance as compared to particle-based implementation for point
object tracking [47]. Moreover, new kinematic and extent states are obtained
by sampling the Gaussian and Wishart distributions, respectively. We note that
it is not difficult to generate random Wishart matrices. Empirical results
show that the runtime taken by sampling from the Gaussian and Wishart
distributions is rather marginal.

For PMB-BP and SPA, the number of particles is 2000, whereas for TPMB-BP the
number of particles is adaptive, and we set $L_{b}=2000$. Empirical results
show that using 2000 particles is a good trade-off between runtime and
estimation performance. For all these three implementations, the number of
message passing iterations is set to 3, and the measurement-driven
initialization of newly detected objects, as described in Section VI-C2, is
used, where samples of state elements that do not depend on the measurements
are drawn from the birth density. In addition, we use message censoring and
measurement reordering, as discussed in [16], to facilitate track
initialization and reduce computational complexity. Furthermore, we prune
Bernoulli components with probability of existence smaller than $10^{-3}$. For
TPMB-BP, we further discard particles using (66) with threshold $10^{-4}$.

For GGIW implementations, we prune global hypotheses (MBs) with weight smaller
than $10^{-2}$ and Bernoulli components with probability of existence smaller
than $10^{-3}$. We also prune GGIW components in the PPP intensity of
undetected objects with weight smaller than $10^{-3}$. In addition, we use
ellipsoidal gating with gate size $13.8$ to reduce computational complexity.
For TPMB(M)-CA and TPMB(M)-SO, we consider $L$-scan implementation with $L=1$,
i.e., no smoothing-while-filtering is performed. Furthermore, Bernoulli
components with probability being alive at current time step smaller than
$10^{-4}$ are not updated. For TPMB(M)-CA, we first apply the density-based
spatial clustering of applications with noise (DB-SCAN) using 200 different
distance values equally spaced between $0.1$ and $20$ to obtain a set of
different measurement partitions, and then for each measurement partition and
global hypothesis $a$, we apply Murty’s algorithm666The C++ implementation in
the Tracker Component Library [71] is used. [72] to find the $\lceil
1000w_{k}^{a}\rceil$ best cluster-to-Bernoulli assignments. For TPMB(M)-SO,
the number of iterations in SO at time step $k$ is set to $10m_{k}$.

For all the implementations, the object/trajectory estimates are extracted
from Bernoulli components with probability of existence no smaller than $0.5$.

###  VII-D Performance evaluation

The trajectory estimation performance is evaluated using the linear
programming (LP) metric $d(\cdot,\cdot)$ for sets of trajectories [56] with
parameters: cut-off distance $20$, order $1$, and track switch cost $2$, see
Appendix E for the detailed expression. The LP metric is integrated with the
Gaussian Wasserstein distance (GWD) for performance evaluation of extended
object estimates with ellipsoidal extent [73]. In the simulated scenarios, we
apply the metric at each time step, and normalize it by the time step. This
enables a comparison of how the LP metric (filtering performance) evolves over
time.

We also apply the LP metric [56] at the final time step to measure the
smoothing error. Only the smoothing of kinematic states is performed, which is
given by the backward simulation (with 100 iterations) and Rauch-Tung-Striebel
smoother [65] for particle-based and GGIW implementations, respectively.

###  VII-E Results

(a) $\gamma=3$

(b) $\gamma=5$

(c) $\gamma=7$

Figure 4: Filtering performance with $\gamma\in\\{3,5,7\\}$ in terms of the
normalized LP trajectory metric over time. The line of PMB-BP almost overlaps
the line of SPA.

(a) $\gamma=3$

(b) $\gamma=5$

(c) $\gamma=7$

Figure 5: Smoothing performance (evaluated at the final time step) with
$\gamma\in\\{3,5,7\\}$ in terms of the LP trajectory metric over time. The
line of PMB-BP almost overlaps the line of SPA.

The filtering and smoothing estimation errors in terms of the LP trajectory
metric are presented in Fig. 4 and Fig. 5, respectively. In addition, the
decompositions of the trajectory metrics into costs due to state estimation
error, missed and false detections, as well as track switches, are given in
Table I and Table II. The results show that implementations using particle BP
have the best estimation performance, followed by implementations using SO.

The different implementations have similar estimation performance when objects
are well-spaced. When objects become closely-spaced, implementations using C&A
and SO present increased estimation error, and in particular, implementations
using C&A suffer from false detection error. This is because the clustering
algorithm does not yield reasonable results when objects are in proximity: the
generated measurement partitions contain either many small clusters or a
single big cluster. For TPMB(M)-CA in the simulated scenario, the hypotheses
with many small clusters have higher likelihoods, and therefore they tend to
overestimate the number of objects. This also explains why implementations
using C&A show increased estimation error when the Poisson measurement rate
$\gamma$ increases, which is instead the opposite for implementations using SO
or particle BP that avoid clustering. Compared to implementations using SO,
implementations using particle BP enjoy a better trade-off between estimation
performance and runtime. When the data association uncertainty is high,
implementations using SO need a significant number of iterations in the
sampling step to obtain hypotheses with high likelihoods. However, when the
objects are relatively well-spaced, implementations using SO actually present
better performance. This can be observed in Fig. 5 around time step $15$ to
time step $20$, and the difference is most noticeable in the case $\gamma=7$.

For the three different implementations using particle BP, PMB-BP slightly
outperforms SPA, and the moderate performance advantage of PMB-BP over SPA can
be explained by the fact that PMB-BP also propagates the density of undetected
objects over time. TPMB-BP shows the best trajectory estimation performance
since estimates of trajectory end time at earlier time steps may be improved
at later time steps. This is most obvious for the case $\gamma=3$ where both
PMB-BP and SPA present increased false detection error when objects die. When
$\gamma=3$, the effective probability of detection $0.95$ is smaller than the
object survival probability $0.99$. Therefore, it is likely that the object
death events cannot be immediately reported by the multi-object filter.

The average runtime777MATLAB implementations on a single core of an Intel(R)
Xeon(R) Gold 6226R CPU @ 2.90GHz. of different implementations is given in
Table III. As can be seen, PMB-BP and SPA are the fastest. TPMB-BP is slower
than TPMB-SO, but it has significantly better estimation performance. It has
been observed that a time-consuming part of the current implementation of
TPMB-BP is the appending of particles representing dead trajectories in the
prediction step via struct arrays. Possible improvements may be brought by
using more efficient data structures.

TABLE I: Filtering performance in terms of the normalized LP trajectory metric
(summed over all the time steps) and its decomposition

| $\gamma=3$ | $\gamma=5$ | $\gamma=7$
---|---|---|---
| Total | State | Miss | False | Switch | Total | State | Miss | False | Switch | Total | State | Miss | False | Switch
TPMB-BP | $\underline{1365.1}$ | 1230.4 | 121.2 | 11.0 | 2.6 | $\underline{1039.0}$ | 1009.0 | 20.4 | 8.1 | 1.5 | $\underline{921.7}$ | 890.0 | 3.9 | 26.4 | 1.3
PMB-BP | 1448.4 | 1222.2 | 175.0 | 48.2 | 3.1 | 1091.4 | 995.8 | 75.7 | 18.5 | 1.4 | 962.4 | 886.2 | 21.7 | 53.5 | 1.1
SPA | 1452.1 | 1225.5 | 174.4 | 49.2 | 3.1 | 1091.4 | 996.4 | 75.4 | 18.1 | 1.4 | 964.1 | 886.1 | 23.2 | 54.4 | 1.1
TPMBM-CA | 2128.5 | 1394.2 | 337.4 | 379.7 | 17.2 | 2406.9 | 1258.0 | 99.4 | 1021.3 | 28.2 | 2827.4 | 1173.5 | 61.5 | 1553.9 | 38.5
TPMB-CA | 2100.3 | 1374.9 | 319.0 | 390.1 | 16.3 | 2427.2 | 1236.8 | 86.2 | 1077.0 | 27.2 | 2846.9 | 1164.6 | 55.4 | 1589.3 | 37.6
TPMBM-SO | 1710.8 | 1247.8 | 410.6 | 45.8 | 6.6 | 1241.5 | 1042.9 | 145.8 | 48.0 | 4.9 | 1070.0 | 923.8 | 93.9 | 47.8 | 4.5
TPMB-SO | 1726.1 | 1242.2 | 396.4 | 81.1 | 6.3 | 1296.6 | 1057.6 | 128.7 | 104.8 | 5.4 | 1120.7 | 941.8 | 81.8 | 92.3 | 4.9

TABLE II: Smoothing performance (evaluated at the final time step) in terms of
the LP trajectory metric and its decomposition

| $\gamma=3$ | $\gamma=5$ | $\gamma=7$
---|---|---|---
| Total | State | Miss | False | Switch | Total | State | Miss | False | Switch | Total | State | Miss | False | Switch
TPMB-BP | $\underline{906.1}$ | 832.9 | 65.8 | 6.4 | 1.0 | $\underline{698.9}$ | 681.2 | 11.4 | 5.7 | 0.6 | $\underline{623.7}$ | 606.6 | 1.8 | 14.8 | 0.5
PMB-BP | 1052.9 | 840.2 | 75.1 | 136.4 | 1.2 | 724.2 | 679.7 | 21.0 | 22.9 | 0.6 | 654.6 | 606.5 | 4.8 | 42.9 | 0.5
SPA | 1058.0 | 843.3 | 76.3 | 137.2 | 1.1 | 725.7 | 683.0 | 19.6 | 22.4 | 0.6 | 659.8 | 606.7 | 7.7 | 44.9 | 0.5
TPMBM-CA | 1744.0 | 956.0 | 177.8 | 596.8 | 13.4 | 2464.0 | 900.2 | 56.5 | 1484.7 | 22.5 | 3140.8 | 867.3 | 39.9 | 2201.8 | 31.8
TPMB-CA | 1726.7 | 911.5 | 161.3 | 643.2 | 10.6 | 2483.5 | 859.2 | 47.1 | 1557.5 | 19.5 | 3184.4 | 845.3 | 35.4 | 2275.7 | 28.0
TPMBM-SO | 1117.5 | 820.7 | 212.3 | 79.7 | 4.9 | 810.2 | 671.7 | 83.7 | 51.0 | 3.7 | 723.0 | 602.5 | 60.4 | 56.8 | 3.4
TPMB-SO | 1215.6 | 789.2 | 194.2 | 228.7 | 3.5 | 959.9 | 672.0 | 66.2 | 218.4 | 3.2 | 886.5 | 618.3 | 46.5 | 217.8 | 3.9

TABLE III: Average runtime (in seconds) for a complete Monte Carlo run of different implementations | $\gamma=3$ | $\gamma=5$ | $\gamma=7$
---|---|---|---
TPMB-BP | 88.2 | 86.6 | 107.5
PMB-BP | 26.5 | 42.2 | 60.6
SPA | 26.1 | 42.2 | 60.7
TPMBM-CA | 144.7 | 304.4 | 452.8
TPMB-CA | 103.4 | 236.3 | 375.4
TPMBM-SO | 699.3 | 1330.9 | 2081.0
TPMB-SO | 34.7 | 60.5 | 91.6

##  VIII Conclusion

In this paper, we present a PMBM conjugate prior on the posterior of sets of
trajectories for a generalized measurement model. We also present the factor
graph representation of the joint posterior of the PMBM set of trajectories
and association variables for the Poisson spatial measurement model. Based on
these important theoretical contributions, we present two TPMB filters for
multiple EOT implemented using particle BP: one estimates the set of alive
trajectories, and the other estimates the set of all trajectories. The
proposed implementations show excellent performance advantages when objects
move in proximity, while also providing full trajectory information.

For the future work, it would be interesting to investigate how to extend the
current work to consider multi-scan data associations [74] and tracking co-
existing point and extended objects [35]. Another interesting future work
direction is to study how to integrate the particle flow filter [75], which
does not suffer from the particle path degeneracy problem, into filters based
on sets of trajectories, such that one can directly extract trajectory
estimates from the multi-trajectory posterior density.

## References

  * [1] S. Challa, M. R. Morelande, D. Mušicki, and R. J. Evans, _Fundamentals of object tracking_. Cambridge University Press, 2011.
  * [2] R. L. Streit, R. B. Angle, and M. Efe, _Analytic combinatorics in multiple object tracking_. Springer, 2021.
  * [3] K. Granström, M. Baum, and S. Reuter, “Extended object tracking: Introduction, overview and applications,” _Journal of Advances in Information Fusion_ , vol. 12, no. 2, 2017.
  * [4] K. Granström and M. Baum, “A tutorial on multiple extended object tracking,” _TechRxiv_ , 2022.
  * [5] K. Gilholm and D. Salmond, “Spatial distribution model for tracking extended objects,” _IEE Proceedings-Radar, Sonar and Navigation_ , vol. 152, no. 5, pp. 364–371, 2005.
  * [6] A. Swain and D. Clark, “Extended object filtering using spatial independent cluster processes,” in _13th International Conference on Information Fusion_. IEEE, 2010, pp. 1–8.
  * [7] J. Honer and H. Kaulbersch, “Bayesian extended target tracking with automotive radar using learned spatial distribution models,” in _International Conference on Multisensor Fusion and Integration for Intelligent Systems_. IEEE, 2020, pp. 316–322.
  * [8] M. Buhren and B. Yang, “Simulation of automotive radar target lists using a novel approach of object representation,” in _Intelligent Vehicles Symposium_. IEEE, 2006, pp. 314–319.
  * [9] L. Hammarstrand, L. Svensson, F. Sandblom, and J. Sorstedt, “Extended object tracking using a radar resolution model,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 48, no. 3, pp. 2371–2386, 2012.
  * [10] C. Knill, A. Scheel, and K. Dietmayer, “A direct scattering model for tracking vehicles with high-resolution radars,” in _Intelligent Vehicles Symposium_. IEEE, 2016, pp. 298–303.
  * [11] S. P. Coraluppi and C. A. Carthel, “Multiple-hypothesis tracking for targets producing multiple measurements,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 54, no. 3, pp. 1485–1498, 2018.
  * [12] X. Tang, M. Li, R. Tharmarasa, and T. Kirubarajan, “Seamless tracking of apparent point and extended targets using Gaussian process PMHT,” _IEEE Transactions on Signal Processing_ , vol. 67, no. 18, pp. 4825–4838, 2019.
  * [13] R. Streit, “JPDA intensity filter for tracking multiple extended objects in clutter,” in _Proceedings of International Conference on Information Fusion_. IEEE, 2016, pp. 1477–1484.
  * [14] S. Yang, K. Thormann, and M. Baum, “Linear-time joint probabilistic data association for multiple extended object tracking,” in _10th Sensor Array and Multichannel Signal Processing Workshop_. IEEE, 2018, pp. 6–10.
  * [15] S. Yang, L. M. Wolf, and M. Baum, “Marginal association probabilities for multiple extended objects without enumeration of measurement partitions,” in _23rd International Conference on Information Fusion_ , 2020, pp. 1–8.
  * [16] F. Meyer and J. L. Williams, “Scalable detection and tracking of extended objects,” in _International Conference on Acoustics, Speech and Signal Processing_. IEEE, 2020, pp. 8916–8920.
  * [17] F. Meyer and J. Williams, “Scalable detection and tracking of geometric extended objects,” _IEEE Transactions on Signal Processing_ , vol. 69, no. 6283–6298, 2021.
  * [18] Y. Li, P. Wei, Y. Wei, L. Gao, and H. Zhang, “Loopy sum-product algorithm based joint detection, tracking and classification of extended objects with analytic implementations,” _Signal Processing_ , p. 108520, 2022.
  * [19] A. De Freitas, L. Mihaylova, A. Gning, M. Schikora, M. Ulmke, D. Angelova, and W. Koch, “A box particle filter method for tracking multiple extended objects,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 55, no. 4, pp. 1640–1655, 2018.
  * [20] R. P. Mahler, _Advances in Statistical Multisource-Multitarget Information Fusion_. Artech House Norwood, MA, 2014\.
  * [21] A. Swain and D. Clark, “The PHD filter for extended target tracking with estimable extent shape parameters of varying size,” in _15th International Conference on Information Fusion_. IEEE, 2012, pp. 1111–1118.
  * [22] K. Granström, C. Lundquist, and O. Orguner, “Extended target tracking using a Gaussian-mixture PHD filter,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 48, no. 4, pp. 3268–3286, 2012.
  * [23] K. Granström and U. Orguner, “A PHD filter for tracking multiple extended targets using random matrices,” _IEEE Transactions on Signal Processing_ , vol. 60, no. 11, pp. 5657–5671, 2012.
  * [24] J. Böhler, T. Baur, S. Wirtensohn, and J. Reuter, “Stochastic partitioning for extended object probability hypothesis density filters,” in _Sensor Data Fusion: Trends, Solutions, Applications_. IEEE, 2019, pp. 1–6.
  * [25] S. Yang, F. Teich, and M. Baum, “Network flow labeling for extended target tracking PHD filters,” _IEEE Transactions on Industrial Informatics_ , vol. 15, no. 7, pp. 4164–4171, 2019.
  * [26] J. Sjudin, M. Marcusson, L. Svensson, and L. Hammarstrand, “Extended object tracking using sets of trajectories with a PHD filter,” in _24th International Conference on Information Fusion_. IEEE, 2021, pp. 1–8.
  * [27] C. Lundquist, K. Granström, and U. Orguner, “An extended target CPHD filter and a gamma Gaussian inverse Wishart implementation,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 7, no. 3, pp. 472–483, 2013.
  * [28] M. Beard, S. Reuter, K. Granström, B.-T. Vo, B.-N. Vo, and A. Scheel, “Multiple extended target tracking with labeled random finite sets,” _IEEE Transactions on Signal Processing_ , vol. 64, no. 7, pp. 1638–1653, 2016.
  * [29] T. Hirscher, A. Scheel, S. Reuter, and K. Dietmayer, “Multiple extended object tracking using Gaussian processes,” in _19th International Conference on Information Fusion_. IEEE, 2016, pp. 868–875.
  * [30] A. Daniyan, S. Lambotharan, A. Deligiannis, Y. Gong, and W.-H. Chen, “Bayesian multiple extended target tracking using labeled random finite sets and splines,” _IEEE Transactions on Signal Processing_ , vol. 66, no. 22, pp. 6076–6091, 2018.
  * [31] A. Scheel and K. Dietmayer, “Tracking multiple vehicles using a variational radar model,” _IEEE Transactions on Intelligent Transportation Systems_ , vol. 20, no. 10, pp. 3721–3736, 2018.
  * [32] K. Granström, M. Fatemi, and L. Svensson, “Poisson multi-Bernoulli mixture conjugate prior for multiple extended target filtering,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 56, no. 1, pp. 208–225, 2020.
  * [33] K. Granström, L. Svensson, S. Reuter, Y. Xia, and M. Fatemi, “Likelihood-based data association for extended object tracking using sampling methods,” _IEEE Transactions on Intelligent Vehicles_ , vol. 3, no. 1, pp. 30–45, 2018.
  * [34] Y. Xia, K. Granström, L. Svensson, Á. F. García-Fernández, and J. L. Williams, “Extended target Poisson multi-Bernoulli mixture trackers based on sets of trajectories,” in _22th International Conference on Information Fusion_. IEEE, 2019, pp. 1–8.
  * [35] Á. F. García-Fernández, J. L. Williams, L. Svensson, and Y. Xia, “A Poisson multi-Bernoulli mixture filter for coexisting point and extended targets,” _IEEE Transactions on Signal Processing_ , vol. 69, pp. 2600–2610, 2021.
  * [36] M. Fröhle, K. Granström, and H. Wymeersch, “Decentralized Poisson multi-Bernoulli filtering for vehicle tracking,” _IEEE Access_ , vol. 8, pp. 126 414–126 427, 2020.
  * [37] Y. Xia, K. Granström, L. Svensson, M. Fatemi, Á. F. García-Fernández, and J. L. Williams, “Poisson multi-Bernoulli approximations for multiple extended object filtering,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 58, no. 2, pp. 890–906, 2022.
  * [38] E. H. Aoki, P. K. Mandal, L. Svensson, Y. Boers, and A. Bagchi, “Labeling uncertainty in multitarget tracking,” _IEEE Transactions on Aerospace and Electronic systems_ , vol. 52, no. 3, pp. 1006–1020, 2016.
  * [39] K. Granström, L. Svensson, Y. Xia, J. Williams, and Á. F. García-Femández, “Poisson multi-Bernoulli mixture trackers: Continuity through random finite sets of trajectories,” in _21st International Conference on Information Fusion_. IEEE, 2018, pp. 1–8.
  * [40] Á. F. García-Fernández, L. Svensson, and M. R. Morelande, “Multiple target tracking based on sets of trajectories,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 56, no. 3, pp. 1685–1707, 2019.
  * [41] J. Honer and F. Schmieder, “Gibbs sampling of measurement partitions and associations for extended multi-target tracking,” in _22th International Conference on Information Fusion_. IEEE, 2019, pp. 1–8.
  * [42] F. Meyer and M. Z. Win, “Data association for tracking extended targets,” in _Military Communications Conference_. IEEE, 2019, pp. 337–342.
  * [43] F. Meyer, Z. Liu, and M. Z. Win, “Scalable probabilistic data association with extended objects,” in _IEEE International Conference on Communications Workshops_. IEEE, 2019, pp. 1–6.
  * [44] F. Meyer and M. Z. Win, “Scalable data association for extended object tracking,” _IEEE Transactions on Signal and Information Processing over Networks_ , vol. 6, pp. 491–507, 2020.
  * [45] J. L. Williams, “Marginal multi-Bernoulli filters: RFS derivation of MHT, JIPDA, and association-based member,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 51, no. 3, pp. 1664–1687, 2015.
  * [46] Á. F. García-Fernández, J. L. Williams, K. Granström, and L. Svensson, “Poisson multi-Bernoulli mixture filter: direct derivation and implementation,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 54, no. 4, pp. 1883–1901, 2018.
  * [47] F. Meyer, T. Kropfreiter, J. L. Williams, R. Lau, F. Hlawatsch, P. Braca, and M. Z. Win, “Message passing algorithms for scalable multitarget tracking,” _Proceedings of the IEEE_ , vol. 106, no. 2, pp. 221–259, 2018.
  * [48] T. Kropfreiter, F. Meyer, and F. Hlawatsch, “Sequential Monte Carlo implementation of the track-oriented marginal multi-Bernoulli/Poisson filter,” in _19th International Conference on Information Fusion_. IEEE, 2016, pp. 972–979.
  * [49] ——, “A fast labeled multi-Bernoulli filter using belief propagation,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 56, no. 3, pp. 2478–2488, 2020.
  * [50] F. Meyer and J. L. Williams, “Scalable detection and tracking of geometric extended objects: Supplementary material,” _arXiv preprint arXiv:2103.11279_ , 2021.
  * [51] F. Lindsten and T. B. Schön, “Backward simulation methods for Monte Carlo statistical inference,” _Foundations and Trends® in Machine Learning_ , vol. 6, no. 1, pp. 1–143, 2013\.
  * [52] Á. F. García-Fernández, L. Svensson, J. L. Williams, Y. Xia, and K. Granström, “Trajectory Poisson multi-Bernoulli filters,” _IEEE Transactions on Signal Processing_ , vol. 68, pp. 4933–4945, 2020.
  * [53] R. L. Graham, D. E. Knuth, O. Patashnik, and S. Liu, “Concrete mathematics: a foundation for computer science,” _Computers in Physics_ , vol. 3, no. 5, pp. 106–107, 1989.
  * [54] B.-N. Vo, M. Mallick, Y. Bar-Shalom, S. Coraluppi, R. Osborne, R. Mahler, and B.-t. Vo, “Multitarget tracking,” _Wiley encyclopedia of electrical and electronics engineering_ , 2015.
  * [55] Á. F. García-Fernández and L. Svensson, “Trajectory PHD and CPHD filters,” _IEEE Transactions on Signal Processing_ , vol. 67, no. 22, pp. 5702–5714, 2019.
  * [56] Á. F. García-Fernández, A. S. Rahmathullah, and L. Svensson, “A metric on the space of finite sets of trajectories for evaluation of multi-target tracking algorithms,” _IEEE Transactions on Signal Processing_ , vol. 68, pp. 3917–3928, 2020.
  * [57] H. A. Blom and Y. Bar-Shalom, “The interacting multiple model algorithm for systems with Markovian switching coefficients,” _IEEE Transactions on Automatic Control_ , vol. 33, no. 8, pp. 780–783, 1988.
  * [58] Y. Xia, K. Granström, L. Svensson, Á. F. García-Fernández, and J. L. Williams, “Multi-scan implementation of the trajectory Poisson multi-Bernoulli mixture filter,” _Journal of Advances in Information Fusion_ , vol. 14, no. 2, pp. 213–235, 2019.
  * [59] K. Granström, L. Svensson, Y. Xia, J. Williams, and Á. F. García-Fernández, “Poisson multi-Bernoulli mixtures for sets of trajectories,” _arXiv preprint arXiv:1912.08718_ , 2019.
  * [60] M. Fatemi, K. Granström, L. Svensson, F. J. Ruiz, and L. Hammarstrand, “Poisson multi-Bernoulli mapping using Gibbs sampling,” _IEEE Transactions on Signal Processing_ , vol. 65, no. 11, pp. 2814–2827, 2017.
  * [61] K. Granström and U. Orguner, “On the reduction of Gaussian inverse Wishart mixtures,” in _International Conference on Information Fusion_. IEEE, 2012, pp. 2162–2169.
  * [62] K. Thormann and M. Baum, “Fusion of elliptical extended object estimates parameterized with orientation and axes lengths,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 57, no. 4, pp. 2369–2382, 2021.
  * [63] F. R. Kschischang, B. J. Frey, and H.-A. Loeliger, “Factor graphs and the sum-product algorithm,” _IEEE Transactions on Information Theory_ , vol. 47, no. 2, pp. 498–519, 2001.
  * [64] Y. Xia, L. Svensson, Á. F. García-Fernández, J. L. Williams, D. Svensson, and K. Granström, “Multiple object trajectory estimation using backward simulation,” _IEEE Transactions on Signal Processing_ , vol. 70, pp. 3249–3263, 2022.
  * [65] S. Särkkä, _Bayesian filtering and smoothing_. Cambridge University Press, 2013, no. 3.
  * [66] M. S. Arulampalam, S. Maskell, N. Gordon, and T. Clapp, “A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking,” _IEEE Transactions on Signal Processing_ , vol. 50, no. 2, pp. 174–188, 2002\.
  * [67] P. Bostrom-Rost, D. Axehill, and G. Hendeby, “PMBM filter with partially grid-based birth model with applications in sensor management,” _IEEE Transactions on Aerospace and Electronic Systems_ , 2021.
  * [68] S. Saha, P. K. Mandal, A. Bagchi, Y. Boers, and J. N. Driessen, “Particle based smoothed marginal MAP estimation for general state space models,” _IEEE Transactions on Signal Processing_ , vol. 61, no. 2, pp. 264–273, 2012\.
  * [69] P. R. Horridge and S. Maskell, “Using a probabilistic hypothesis density filter to confirm tracks in a multi-target environment,” in _GI-Jahrestagung_ , 2011.
  * [70] M. Feldmann, D. Franken, and W. Koch, “Tracking of extended objects and group targets using random matrices,” _IEEE Transactions on Signal Processing_ , vol. 59, no. 4, pp. 1409–1420, 2011.
  * [71] D. F. Crouse, “The tracker component library: free routines for rapid prototyping,” _IEEE Aerospace and Electronic Systems Magazine_ , vol. 32, no. 5, pp. 18–27, 2017.
  * [72] ——, “On implementing 2D rectangular assignment algorithms,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 52, no. 4, pp. 1679–1696, 2016.
  * [73] S. Yang, M. Baum, and K. Granström, “Metrics for performance evaluation of elliptic extended object tracking methods,” in _International Conference on Multisensor Fusion and Integration for Intelligent Systems_. IEEE, 2016, pp. 523–528.
  * [74] J. L. Williams and R. A. Lau, “Multiple scan data association by convex variational inference,” _IEEE Transactions on Signal Processing_ , vol. 66, no. 8, pp. 2112–2127, 2018.
  * [75] S. Mori and D. F. Crouse, “A note on particle flow methods for solving bayesian updates,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 58, no. 5, pp. 4796–4802, 2022.

Supplementary Materials

## Appendix A

In this appendix, we present a short proof of Proposition 2.

It has been shown in [35, Theorem 1] that the update of a PMBM prior with
single-object state in a LCHS space and the generalized measurement model in
Section II-C2 is also a PMBM. Furthermore, given that the single-object state
space is LCHS, the single-trajectory space $T_{(k^{\prime})}$ is also LCHS
[40, Appendix A]. This means that finite set statistics [20] can be used on
sets of trajectories. Therefore, Proposition 2 can be understood as an
extension of the PMBM update for sets of objects [35, Theorem 1] to sets of
trajectories.

## Appendix B

In this appendix, we prove Proposition 3. In what follows, we first present
the joint posterior of set of trajectories and global hypothesis augmented
with auxiliary variables, which does not involve the summation over set
partitions (cf. (9)). Then, we proceed to describe how to rewrite the joint
posterior using the object-oriented and measurement-oriented association
variables [16]. At last, we show how to further simplify such joint posterior
to exclude the object-oriented association variables.

###  B-A Joint posterior of trajectories and global hypothesis

Given a predicted TPMB density at time step $k$ of the form (11a) and
measurements ${\bf z}_{k}$, the joint posterior of set of trajectory and
global hypothesis is

$\displaystyle\widetilde{f}_{k|k}(\widetilde{{\bf X}}_{k},a)$
$\displaystyle=\widetilde{f}_{k|k}(\widetilde{{\bf
X}}_{k}|a)w^{a}_{k|k}\psi(a)$
$\displaystyle=w^{a}_{k|k}\widetilde{f}^{p}_{k|k}(\widetilde{{\bf
Y}}_{k})\prod_{i=1}^{n_{k|k}}\widetilde{f}^{i,a^{i}}_{k|k}\left(\widetilde{{\bf
X}}_{k}^{i}\right)\psi(a)$
$\displaystyle\propto\widetilde{f}^{p}_{k|k}(\widetilde{{\bf
Y}}_{k})\prod_{i=1}^{n_{k|k}}w^{i,a^{i}}_{k|k}\widetilde{f}^{i,a^{i}}_{k|k}\left(\widetilde{{\bf
X}}^{i}_{k}\right)\psi(a)$
$\displaystyle=\widetilde{f}^{p}_{k|k}(\widetilde{{\bf
Y}}_{k})\prod_{i=1}^{n_{k|k}}g^{i,a^{i}}_{k|k}\left(\widetilde{{\bf
X}}^{i}_{k}\right)\psi(a)$ (70)

where the constraint on global hypothesis $a$ can be expressed as an indicator
function

$\psi(a)=\begin{cases}1&\biguplus_{i=1}^{n_{k|k}}{\cal M}_{k}^{i,a^{i}}={\cal
M}_{k}\\\ 0&\text{otherwise},\end{cases}$ (71)

and we introduce the unnormalized Bernoulli densities
$g^{i,a^{i}}_{k|k}(\cdot)$, each of which is given by the product of a local
hypothesis density $\widetilde{f}^{i,a^{i}}_{k|k}(\widetilde{{\bf
X}}^{i}_{k})$ and its corresponding weight $w^{i,a^{i}}_{k|k}$. Using their
expressions presented in Proposition 2, we have, for
$i\in\\{1,\dots,n_{k|k-1}\\}$,

$\displaystyle g^{i,1}_{k|k}\left(\widetilde{{\bf X}}\right)$
$\displaystyle=\begin{cases}r^{i}_{k|k-1}\ell_{k}(\emptyset|X)f^{i}_{k|k-1}(X)\delta_{i}[u]&\widetilde{{\bf
X}}=\\{(u,X)\\}\\\ 1-r^{i}_{k|k-1}&\widetilde{{\bf X}}=\emptyset\\\
0&\text{otherwise},\end{cases}$ (72)

for $i\in\\{1,\dots,n_{k|k-1}\\}$, $j\in\\{2,\dots,2^{m_{k}}-1\\}$,
$a^{i}=j+1$,

$\displaystyle g_{k|k}^{i,a^{i}}\left(\widetilde{{\bf X}}\right)$
$\displaystyle=\begin{cases}r_{k|k-1}^{i}\ell_{k}\left({\bf
w}_{k}^{j}|X\right)f^{i}_{k|k-1}(X)\delta_{i}[u]&\widetilde{{\bf
X}}=\\{(u,X)\\}\\\ 0&\text{otherwise},\end{cases}$ (73)

for $i\in\\{n_{k|k-1}+1,\dots,n_{k|k-1}+m_{k}\\}$,

$g^{i,1}_{k|k}\left(\widetilde{{\bf X}}\right)=\begin{cases}1&\widetilde{{\bf
X}}=\emptyset\\\ 0&\text{otherwise},\end{cases}$ (74)

and for $i=n_{k|k-1}+j,j\in\\{1,\dots,m_{k}\\},\iota\in\\{1,\dots,2^{j-1}\\}$,
$a^{i}=\iota+1$

$\displaystyle g^{i,a^{i}}_{k|k}\left(\widetilde{{\bf X}}\right)$
$\displaystyle=\begin{cases}\ell_{k}\left({\bf
w}_{k}^{i,\iota}|X\right)\lambda_{k|k-1}(X)\delta_{i}[u]&\widetilde{{\bf
X}}=\\{(u,X)\\}\\\ \lambda_{k}^{C}\left(z_{k}^{j}\right)&\widetilde{{\bf
X}}=\emptyset,{\bf w}_{k}^{i,\iota}=\left\\{z_{k}^{j}\right\\}\\\
0&\text{otherwise}.\end{cases}$ (75)

###  B-B Data association variables

Until now, we have used variable $a\in{\cal A}_{k|k}$ to represent the global
hypotheses. To obtain a neat factor graph representation of the joint
posterior, it is useful to describe the Bernoulli-to-measurement association
${\cal M}_{k}^{i,a^{i}}$ using association vectors. We first introduce the
binary object-oriented association vector
$\alpha_{k}=[{\alpha^{1}_{k}}^{T},\dots,{\alpha_{k}^{n_{k|k}}}^{T}]^{T}$ where

$\alpha_{k}^{i}=\begin{cases}\left[\alpha_{k}^{i,1},\dots,\alpha_{k}^{i,m_{k}}\right]^{T},&i\in\\{1,\dots,n_{k|k-1}\\}\\\
\left[\alpha_{k}^{i,1},\dots,\alpha_{k}^{i,j}\right]^{T},&i=n_{k|k-1}+j,j\in\\{1,\dots,m_{k}\\}\end{cases}$

and $\alpha_{k}^{i,j}=1$ if and only if the $j$-th measurement $z_{k}^{j}$ is
associated to the $i$-th potential object (described by the $i$-th Bernoulli
component). Note that the object-oriented association vector corresponding to
the $j$-th new Bernoulli component has length $j$. This is a direct result of
the local hypothesis representation for new Bernoulli components specified in
Proposition 2.

For TPMB filters888For the TPMB filter, the global hypothesis does not include
the measurement association history., from global hypothesis $a\in{\cal
A}_{k|k}$ (17) we can represent the data associations at time step $k$ by
object-oriented or measurement-oriented association variables; the latter has
been introduced in Section IV-A. Using a hybrid representation for both
object-oriented and measurement-oriented association variables makes it
possible for developing many scalable MOT algorithms using BP [47]. Later in
Section B-D, we will show that, for EOT with Poisson spatial model, the factor