# Three-dimensional Tracking of a Large Number of High Dynamic Objects from
Multiple Views using Current Statistical Model

Nianhao Xie
College of Aerospace Science and Engineering
National University of Defense Technology
Changsha, China
<EMAIL_ADDRESS>

###### Abstract

Three-dimensional tracking of multiple objects from multiple views has a wide
range of applications, especially in the study of bio-cluster behavior which
requires precise trajectories of research objects. However, there are
significant temporal-spatial association uncertainties when the objects are
similar to each other, frequently maneuver, and cluster in large numbers.
Aiming at such a multi-view multi-object 3D tracking scenario, a current
statistical model based Kalman particle filter (CSKPF) method is proposed
following the Bayesian tracking-while-reconstruction framework. The CSKPF
algorithm predicts the objects’ states and estimates the objects’ state
covariance by the current statistical model to importance particle sampling
efficiency, and suppresses the measurement noise by the Kalman filter. The
simulation experiments prove that the CSKPF method can improve the tracking
integrity, continuity, and precision compared with the existing constant
velocity based particle filter (CVPF) method. The real experiment on fruitfly
clusters also confirms the effectiveness of the CSKPF method.

_K_ eywords multiple object tracking $\cdot$ 3D reconstruction $\cdot$ current
statistical model $\cdot$ bio-cluster $\cdot$ Drosophila melanogaster

##  1 Introduction

Multi-object tracking (MOT) are widely applied on traffic control, sports
analysis, battlefield reconnaissance, etc. In the tracking of a biol-cluster,
e.g. fish school, birds flock and insect swarm [1, 2, 3] , there are more
objects and higher motion kinematics compared with the traditional MOT
scenario, which makes the general MOT methods not applicable. Moreover, 3D
trajectories of biol-cluster are needed which requires images from
synchronized multiple cameras to do object 3D reconstruction at the same time.

The 3D multi-view multi-object tracking faces two major challenges compared
with the single-object tracking [4] or static objects 3D reconstruction [5] .
One is temporal association problem, i.e. how to determine the detection111For
the sake of clarity, the concepts of _measurement_ , _detection_ , and
_obervation_ in this paper are firstly be claimed. On an image, the foreground
pixels of objects form many _measurements_. The measurement of a certain
object is its _detection_ and the state reconstructed by the detection is its
_observation_. association of each object between the current frame and the
previous frames on a certain view. The other one is the spatial association
problem, i.e. how to determine the detection association of each object on
different views in a certain frame. Muti-view multi-object tracking is
essentially a problem of assigning temporal-spatial association of detections.
The temporal association decides the 2D tracking while the spatial association
decides the 3D reconstruction.

For the temporal association problem, the associated detection on the current
frame to the previous frame can be constrained to a smaller region using
motion continuity constraints. For the spatial association problem, the
associated detection on one view to another view can be constrained on the
epipolar. Assuming that the association credibility between two measurements
can be defined, there are statistical-based methods that assign the nearest
measurement to the target detection, such as joint probability data
association filter (JPDAF) and multiple hypothesis tracking (MHT), which is
effective when the number of objects is not large and the measurement
distinction is not small [6] . The methods based on linear allocation
establish a cost matrix between the target detection and the candidate
measurement, and then the Hungarian algorithm is used to solve the optimal
detection-measurement allocation problem[7] .

The order of reconstruction and tracking is important. According to the
relationship between the temporal association and spatial association, the
existing multi-view multi-object tracking method can be divided into three
categories: a) Firstly, temporal associations are assigned to generate 2D
trajectories on each view, and then the spatial association of these 2D
trajectories is assigned to reconstruct 3D trajectories, named tracking-to-
reconstruction method [8] . This type of algorithm cannot effectively use the
3D motion prediction of the object to deal with the occlusion problem, which
often leads to trajectory interruptions and switch. b) Firstly, the spatial
association is assigned to reconstruct the 3D points in each frame, and then
the temporal association of these 3D points is assigned to link 3D
trajectories, named reconstruction-to-tracking method [9] . This type of
method often results in many illusory 3D trajectories when the measurement is
not discriminative enough, because the spatial association will generate many
"ghost" 3D points. c) Assuming the initial 3D position of each object is
known, firstly, the 3D position in the next frame is predicted by the
kinematic model and reprojected to each view by spatial association
constraints. The prediction can be confirmed if there is overlap between the
reprojected detection and the measurements on every view, which makes the
temporal association natural satisfies the spatial association constraint.
This type of method, named tracking-while-reconstruction method [10] , can
effectively eliminate phantoms and is insensitive to occlusion, which can
obtain continuous and reliable trajectories.

In many cases, especially for bio-clusters, it is difficult to obtain an
object’s kinematic model. Chen et al. [10] used the assumption of uniform
linear motion to predict the object’s position which is reasonable when the
frame rate is large enough. Wang et al. [11] used the LSTM network to learn
the kinematic model of a single fruitfly, which makes the position prediction
more reasonable but largely increases calculation burden. Considering that the
individual’s motion is a variable acceleration motion with no sudden
acceleration, in this paper, the current statistical model (CSM) [12] is used
to predict not only the position but also the velocity of the object. We
establish a tracking-while-reconstructing framework based on Bayesian
inference and solve the temporal-spatial association problem by particle
filter method, where the object’s states and covariance are predicted by the
current statistical model. After the detections are assigned to each object,
the state observations are determined. The Kalman filter (KF) is used to
suppress measurement noise and improve tracking precision. The proposed
method, named current statistical model based Kalman particle filter (CSKPF)
tracking algorithm, can not only build the trajectories with higher integrity,
continuity, and precision but also can estimate the target velocity during
tracking, which is essential to analyze the clusters behaviors.

This paper is arranged as follows: Section 1 introduces the significance and
challenges of multi-view multi-object tracking, especially the 3D tracking of
a large number of objects with high kinematics, and summarizes existing
temporal-spatial association methods and tracking-reconstruction frameworks.
Section 2 extends the Bayesian inference framework from single-view single-
object tracking to the multi-view multi-object tracking-while-reconstructing.
Then, the existing constant velocity based particle filter (CVPF) method is
explain in the above framework, which leads to our CSKPF method. Section 3
first defines the evaluation index for multi-view multi-object tracking, and
then the tracking performance of CVPF and CSKPF methods are compared by
simulated and real data, which proves the improvement of the proposed method.
Finally, the paper is summarized in section 4.

##  2 Method

In this section, the Bayesian framework for single-view single-object tracking
is reviewed and the reason why it cannot be directly applied to multi-view
multi-object tracking is analyzed. Then, a multi-view multi-object Bayesian
inference framework to deal with the above problems is built. Finally, the
existing CVPF method and the proposed CSKPF method based on this framework are
explained.

###  2.1 The Bayesian single-view single-object tracking framework

The discrete kinematic model and obervation model of the target can be
formulated as

$\mathbf{x}_{t}=\mathbf{f}\left(\mathbf{x}_{t-1},\mathbf{v}_{t-1}\right),\mathbf{y}_{t}=\mathbf{h}\left(\mathbf{x}_{t},\mathbf{n}_{t}\right),$
(1)

where $\mathbf{x}_{t},\mathbf{y}_{t}$ are the state and obervation of the
target at moment $t$, and $\mathbf{v}_{t},\mathbf{n}_{t}$ are the process
noise and obervation noise, and $\mathbf{f}$ and $\mathbf{h}$ are the state
transition function and measurement function, respectively.

From the Bayesian point of view, the object tracking is to estimate the
current state of $\mathbf{x}_{t}$ based on historical obervation
$\mathbf{y}_{1:t}$. Assuming that the state transition is a first-order Markov
process, and the _posterior_ state probability
$p\left(\mathbf{x}_{t-1}|\mathbf{y}_{1:t-1}\right)$ at time $t-1$ is known,
then the _priori_ state at time $t$ can be _predicted_ by the state transition
model $p\left(\mathbf{x}_{t}|\mathbf{x}_{t-1}\right)$ before getting the
obervation $\mathbf{y}_{t}$ at moment $t$, i.e.

$p\left(\mathbf{x}_{t}|\mathbf{y}_{1:t-1}\right)=\int
p\left(\mathbf{x}_{t}|\mathbf{x}_{t-1}\right)p\left(\mathbf{x}_{t-1}|\mathbf{y}_{1:t-1}\right)\,d\mathbf{x}_{t-1}.$
(2)

After the obervation $\mathbf{y}_{t}$ is known, the _priori_ state can be
_corrected_ to _posterior_ state

$p\left(\mathbf{x}_{t}|\mathbf{y}_{1:t}\right)\propto
p\left(\mathbf{y}_{t}|\mathbf{x}_{t}\right)p\left(\mathbf{x}_{t}|\mathbf{y}_{1:t-1}\right),$
(3)

where $p\left(\mathbf{y}_{t}|\mathbf{x}_{t}\right)$ means the probability that
an obervation $\mathbf{y}_{t}$ occurs when the state is $\mathbf{x}_{t}$.

It seems that by generating as many trackers as there are objects and making
each tracker run independently, the multi-view multi-object tracking problem
can be solved by the above single-view single-object framework. However, on
the one hand, the target detection cannot be determined before the temporal
association of detection; on the other hand, the target observation cannot be
reconstructed before the spatial association of detection. In brief,
$\mathbf{y}_{t}$ is unknown because the temporal-spatial association is
pendent, which causes the _correct_ step of Bayesian inference
unimplementable. An improved Bayesian framework is designed in next subsection
to determine the temporal-spatial association which makes the single-object
Bayesian implementable in a multi-view multi-object tracking scenario.

###  2.2 A Bayesian multi-view multi-object tracking-while-reconstruction
framework

As mention above, the framework treat each object independently. Without loss
of generality, as shown in Figure 1, only one tracker’s temporal-spatial
association Bayesian inference will be explained. The horizontal axis is the
temporal axis driven by the kinematic model, and the left half-plane
represents the previous moment $t-1$, and the right half-plane is the current
moment $t$. The vertical axis is the spatial axis determined by the camera
projection model, the upper half-plane is 3D space, and the lower half-plane
represents the 2D space of all views. Supposed that the state
$\mathbf{x}_{t-1}$ has been tracked according to the historical measurements.
Then the _prior_ state $\mathbf{x}_{t|t-1}$ can be predicted based on the
kinematic model $p\left(\mathbf{x}_{t}|\mathbf{x}_{t-1}\right)$. Finally, the
projection of $\mathbf{x}_{t|t-1}$ on each view can be calculated by the
camera projection model. The state transfers from the third quadrant through
the second and first quadrants to the fourth quadrant is constrained by the
temporal association constraints and spatial association constraints.

Figure 1: A Bayesian multi-view multi-object tracking-while-reconstruction
framework. The filled circles are 3D state and the dashed circles are their
projection on each view. The solid circles are the 2D measurements, which are
colored red/yellow if they are overlapping/non-overlapping with the state
projection.

In the fourth quadrant, the temporal-spatial association can be determined by
the position and appearence similarity between the target and candidate
detection, which means the target observation $\mathbf{y}_{t}$ is known. To be
specific, in the fourth quadrant, if the projection
$\mathbf{P}^{v}(\mathbf{x}_{t|t-1})$ has overlap with measurement
$\chi_{v,t,k_{v}}$ on view $v$, then the candidate
$\mathbf{z}_{t,k}=\langle\chi_{1,t,k_{1}},\cdots,\chi_{V,t,k_{V}}\rangle,k=1,\cdots,K$
is a potential temporal-spatial association. All potential candidates
constitute the _candidate association group_
$\mathbf{z}_{t}=\\{\mathbf{z}_{t,k}\\}_{k=1}^{K}$. The credibility of the
candidate association is defined as

$w(\mathbf{z}_{t,k})=p\left(\mathbf{z}_{t,k}|\mathbf{x}_{t|t-1}\right)=\prod_{v=1}^{V}{\frac{1}{\left(\sqrt{2\pi}\sigma_{v}\right)^{-\frac{V}{2}}}\exp\left({\tau_{v}\left(\mathbf{x}_{t|t-1},\chi_{v,t,k_{v}}\right)-1}\right)},$
(4)

where $\tau_{v}\left(\mathbf{x}_{t|t-1},\chi_{v,t,k_{v}}\right)$ is the
appearence similarity between the projection
$\mathbf{P}^{v}(\mathbf{x}_{t|t-1})$ and measurement $\chi_{v,t,k_{v}}$. Thus,
the target observation is the most reliable candidate association

$\mathbf{y}_{t}=\mathop{\arg\max}\limits_{k}{w(\mathbf{z}_{t,k})}.$ (5)

According to the Bayesian inference formula (3), the corrected target state is

$\mathbf{x}_{t}=\mathop{\arg\max}\limits_{\hat{\mathbf{x}}_{t}}{p\left(\hat{\mathbf{x}}_{t}|\mathbf{y}_{1:t}\right)}.$
(6)

The proposed framework can track muliple objects while reconstruction. The
tracker for each object runs independently in parallel and one measurement can
be associated with more than one object, which improves tracking efficiency
and avoids tracking interruption caused by image occlusion in some views.

###  2.3 The constant velocity based particle filter (CVPF) method

It is diffcult to solve the above tracking framework because there is an
integral term in Equation (2). Therefore, the Monte Carlo sampling method is
used to solve the Bayesian inference problem [13], i.e. particle filter
method. Specifically, Chen et al. [10] proposed the CVPF method (as shown in
Algorithm 1) based on constant velocity kinematic model and sampling
importance resampling (SIR).

Algorithm 1 CVPF algorithm

Parameter: number of particles $N_{p}$, particle variance $\sigma^{2}$
Input: the measurements $\mathcal{B}_{t,v}$ of each frame $t$ and each view
$v$
Output: 3D trajectories of all objects
Initialization:
Calculate temporal-temporal association among
$\mathcal{B}_{1,v},\mathcal{B}_{2,v},v=1,\cdots,V$ using method of exhaustion
and obtain $N_{2}$ matching pairs. All associated measurements on the second
frame are put into the associated measurement set $\mathcal{M}_{2,v}$. Create
$N_{2}$ new trackers and put them into the active tracker set
$\mathcal{T}_{2}$ and initialize the inactive tracker set $\mathcal{L}_{2}$ as
empty set.

1:for $t$ = $3,\cdots,T$ do

2: Initialize $\mathcal{M}_{t,v}=\emptyset$

3: for $\tau$ in $\mathcal{T}_{t-1}$ do

4:  _Predict_ the state $\mathbf{x}_{t|t-1}$ of the tracker $\tau$ from
$\mathbf{x}_{t-1}$ using Equation (7).

5:  _Sample_ $N_{p}$ particles
$\mathbf{x}^{i}_{t|t-1}\sim\mathcal{N}\left(\mathbf{x}_{t|t-1},\sigma^{2}\right)$
and renew their weights by $1/N_{p}$.

6: for $i=1,\cdots,N_{p}$ do

7: Calculate $\mathbf{z}^{i}_{t|t-1}$ and their credibilities for
$\mathbf{x}_{t|t-1}^{i}$ using Equation (4).

8: Determine the particle weights $w_{i}$ using Equation (9).

9: end for

10: if $\forall i,w_{i}=0$ then

11: Move the tracker $\tau$ from the active tracker set $\mathcal{T}_{t-1}$ to
inactive $\mathcal{L}_{t-1}$.

12: else

13:  _Estimate_ the state $\hat{\mathbf{x}}_{t}$ by Equation (10).

14: Calculate $\mathbf{z}_{t}$ its credibilities for $\hat{\mathbf{x}}_{t}$
using Equation (4).

15: Put the most likely association $\mathbf{z}_{t}^{*}$ into the association
measurements set $\mathcal{M}_{t,v}$.

16: Reconstruct the most likely obervation $\mathbf{y}_{t}$ by
$\mathbf{z}_{t}^{*}$ using stereovision method.

17: $\mathbf{x}_{t}\leftarrow\hat{\mathbf{x}}_{t}$.

18: end if

19: end for

20: Do temporal-spatial association among the unassociated measurements
$\mathcal{B}_{t-1,v}/\mathcal{M}_{t-1,v}$,
$\mathcal{B}_{t,v}/\mathcal{M}_{t,v}$ using method of exhaustion and
initialize new trackers $\mathcal{R}_{t}$.

21: Update the tracker set by
$\mathcal{T}_{t}\leftarrow\mathcal{T}_{t-1}\bigcup\mathcal{R}_{t}$.

22:end for

23:Parse the trajectories of all objects from the historical states of
trackers $\mathcal{T}_{T}$ and $\mathcal{L}_{T}$.

####  2.3.1 State prediction

Denote the target state as
$\mathbf{x}=[x,\dot{x},y,\dot{y},z,\dot{z}]^{\top}$, the precited state based
on constant velocity hypothesis is

$\mathbf{x}_{t|t-1}=\mathbf{F}\mathbf{x}_{t},$ (7)

where $\mathbf{F}$ is the state transition matrix

${\bf{F}}=\left[{\begin{array}[]{*{20}{c}}{{{\bf{F}}_{1}}}&{\bf{0}}&{\bf{0}}\\\
{\bf{0}}&{{{\bf{F}}_{1}}}&{\bf{0}}\\\
{\bf{0}}&{\bf{0}}&{{{\bf{F}}_{1}}}\end{array}}\right],~{}\mathrm{and}~{}~{}\bf{F}_{1}=\left[{\begin{array}[]{*{20}{c}}1&\Delta
t\\\ 0&1\end{array}}\right],$ (8)

and $\Delta t$ is the sampling interval.

####  2.3.2 Particle sampling

Suppose there is $N_{p}$ particles, the importance probability density
function of the particle $\mathbf{x}_{t|t-1}^{i},i=1,\cdots,N_{p}$, is
$q\left({{\bf{x}}_{t}^{i}|{\bf{x}}_{t-1}^{i},{{\bf{y}}_{1:t}}}\right)=p\left({{\bf{x}}_{t}^{i}|{\bf{x}}_{t-1}^{i}}\right)$
according to the SIR filter method. The CVPF method assumes that
$p\left({{\bf{x}}_{t}^{i}|{\bf{x}}_{t-1}^{i}}\right)$ obeys Gaussian
distribution with mean value $\mathbf{x}_{t|t-1}$ and variance $\sigma^{2}$.
By sampling $N_{p}$ particles from the the importance probability density
function
$\mathbf{x}_{t|t-1}^{i}\sim\mathcal{N}(\mathbf{x}_{t|t-1},\sigma^{2})$, and
renewing the particle weights as $w_{i}=1/N_{p}$, the candidate association
group $\mathbf{z}_{t|t-1}^{i}$ of $i$-th particle and its credibility can be
calculated according to the tracking framework.

####  2.3.3 State correlation

According to Equation (4), the particle weights are updated to

$w_{i}=\mathop{\max}\limits_{k}{w(\mathbf{z}_{{t|t-1},k}^{i})},$ (9)

where $\mathbf{z}_{{t|t-1},k}^{i}$ is the $k$-th candidate association of the
particle $\mathbf{x}_{t|t-1}^{i}$. Thus, the estimated target state is

${\hat{\mathbf{x}}_{t}}=\sum\limits_{i=1}^{n}{{w_{i}}{\bf{x}}_{t|t-1}^{i}}.$
(10)

The CVPF method takes the ${\hat{\mathbf{x}}_{t}}$ as the state correction
${\mathbf{x}_{t}}$. By reprojecting the state ${\hat{\mathbf{x}}_{t}}$ to each
view, the most likely association $\mathbf{z}_{t}^{*}$ can be determined
according to the Equation (5) and the observation $\mathbf{y}_{t}$ can be
calculated by stereovision method. Finally, new trackers can be initialized
from the unassociated measurement.

###  2.4 Currect statistical model based Kalman particle filter (CSKPF) method

The constant velocity hypothesis of CVPF method makes it be unsuitable for
maneuvering target. Moreover, the particle covariance is always constant
$\sigma^{2}$, which does not reflect the state prediction uncertainty. The
proposed CSKPF method (as shown in Algorithm 2) will improve the CVPF method
by solving the two problems.

Algorithm 2 CSKPF algorithm

Parameter: number of particles $N_{p}$
observation noise matrix $\mathbf{R}$, initial process noise $\mathbf{Q}$,
initial covariance $\mathbf{P}$
warm-up time $T_{h}$, reciprocal of maneuver time $\alpha_{i}$, maximum
possible acceleration $a_{\mathrm{max},i}$
Input: the measurements $\mathcal{B}_{t,v}$ of each frame $t$ and each view
$v$
Output: 3D trajectories of all targets
Initialization: same as the CVPF method.

1:for $t$ = $3,\cdots,T_{h}$ do

2: $\divideontimes$ Execute CVPF and record the acceleration as the current
average acceleration $\bar{\bf{a}}$.

3:end for

4:for $t$ = $T_{h}+1,T_{h}+2,\cdots,T$ do

5: Initialize $\mathcal{M}_{t,v}=\emptyset$

6: for $\tau$ in $\mathcal{T}_{t-1}$ do

7: * _Predict_ the state $\mathbf{x}_{t|t-1}$ of the tracker $\tau$ from
$\mathbf{x}_{t-1}$ using the Equation (11).

8: $\divideontimes$ _Predict_ the covariance $\mathbf{P}_{t|t-1}$ of the
tracker $\tau$ using the Equation (12).

9: * _Sample_ $N_{p}$ particles
$\mathbf{x}_{t|t-1}\sim\mathcal{N}\left(\mathbf{x}_{t|t-1},\mathbf{P}_{t|t-1}\right)$
and renew their weights by $1/N_{p}$.

10: for $i=1,\cdots,N_{p}$ do

11: Calculate $\mathbf{z}^{i}_{t|t-1}$ and their credibilities for
$\mathbf{x}_{t|t-1}^{i}$ using Equation (4).

12: Determine the particle weights $w_{i}$ using Equation (9).

13: end for

14: if $\forall i,w_{i}=0$ then

15: Move the tracker $\tau$ from the active tracker set $\mathcal{T}_{t-1}$ to
inactive $\mathcal{L}_{t-1}$.

16: else

17:  _Estimate_ the state $\hat{\mathbf{x}}_{t}$ by Equation (10).

18: Calculate $\mathbf{z}_{t}$ its credibilities for $\hat{\mathbf{x}}_{t}$
using Equation (4).

19: Put the most likely association $\mathbf{z}_{t}^{*}$ into the association
measurements set $\mathcal{M}_{t,v}$.

20: Reconstruct the most likely obervation $\mathbf{y}_{t}$ by
$\mathbf{z}_{t}^{*}$ using stereovision method.

21: $\divideontimes$ _Correct_ the state $\mathbf{x}_{t}$ by Equation (18) and
$\mathbf{P}_{t}$ by Equation (19).

22: end if

23: end for

24: Do temporal-spatial association among the unassociated measurements
$\mathcal{B}_{t-1,v}/\mathcal{M}_{t-1,v}$,
$\mathcal{B}_{t,v}/\mathcal{M}_{t,v}$ using method of exhaustion and
initialize new trackers $\mathcal{R}_{t}$.

25: Update the tracker set by
$\mathcal{T}_{t}\leftarrow\mathcal{T}_{t-1}\bigcup\mathcal{R}_{t}$.

26:end for

27:Parse the trajectories of all objects from the historical states of
trackers $\mathcal{T}_{T}$ and $\mathcal{L}_{T}$.

Note: * marks the revised steps, and $\divideontimes$ marks additional steps
campared to CVPF method.

####  2.4.1 State prediction

The object motion is modeled by currect statistical model which assumes that
the object acceleration is bounded and obeys modified Rayleigh distribution to
the average acceleration in a past period of time. Denote the target state as
$\mathbf{x}=[x,\dot{x},\ddot{x},y,\dot{y},\ddot{y},z,\dot{z},\ddot{z}]^{\top}$
and predict the target state by

$\mathbf{x}_{t|t-1}=\mathbf{G}\mathbf{x}_{t}+\mathbf{U}\bar{\mathbf{a}},$ (11)

and estimate its covariance by

$\mathbf{P}_{t|t-1}=\mathbf{G}\mathbf{P}_{t-1}\mathbf{G}^{\top}+\mathbf{Q}_{t},$
(12)

where $\bar{\mathbf{a}}$ is the current average acceleration and

${\bf{G}}=\left[{\begin{array}[]{*{20}{c}}{{{\bf{G}}_{1}}}&{\bf{0}}&{\bf{0}}\\\
{\bf{0}}&{{{\bf{G}}_{1}}}&{\bf{0}}\\\
{\bf{0}}&{\bf{0}}&{{{\bf{G}}_{1}}}\end{array}}\right],~{}\mathrm{and}~{}\bf{G}_{1}=\left[{\begin{array}[]{*{20}{c}}1&{\Delta
t}&{\frac{1}{{{\alpha^{2}}}}\left({-1+\alpha\Delta t+{e^{-\alpha\Delta
t}}}\right)}\\\ 0&1&{\frac{1}{\alpha}\left({1-{e^{-\alpha\Delta
t}}}\right)}\\\ 0&0&{{e^{-\alpha\Delta t}}}\end{array}}\right]$ (13)

is the state transition matrix,

${\bf{U}}=\left[{\begin{array}[]{*{20}{c}}{{{\bf{u}}_{1}}}&{\bf{0}}&{\bf{0}}\\\
{\bf{0}}&{{{\bf{u}}_{1}}}&{\bf{0}}\\\
{\bf{0}}&{\bf{0}}&{{{\bf{u}}_{1}}}\end{array}}\right],~{}\mathrm{and}~{}{{\bf{u}}_{1}}=\left[{\begin{array}[]{*{20}{c}}{\frac{1}{{{\alpha^{2}}}}\left({-\alpha\Delta
t+\frac{{{\alpha^{2}}{{\left({\Delta t}\right)}^{2}}}}{2}+1-{e^{-alpha\Delta
t}}}\right)}\\\ {\frac{1}{\alpha}\left({\alpha\Delta t-1+{e^{-\alpha\Delta
t}}}\right)}\\\ {1-{e^{-\alpha\Delta t}}}\end{array}}\right]$ (14)

is the input control matrix. The the process noise estimation matrix

${\bf{Q}}=\left[{\begin{array}[]{*{20}{l}}{{{\bf{Q}}_{1}}}&{\bf{0}}&{\bf{0}}\\\
{\bf{0}}&{{{\bf{Q}}_{2}}}&{\bf{0}}\\\
{\bf{0}}&{\bf{0}}&{{{\bf{Q}}_{3}}}\end{array}}\right],~{}\mathrm{and}~{}\begin{array}[]{l}{{\bf{Q}}_{i}}=2\alpha_{i}\sigma_{i}^{2}\left[{\begin{array}[]{*{20}{c}}{{q_{11}(\alpha_{i})}}&{{q_{12}(\alpha_{i})}}&{{q_{13}(\alpha_{i})}}\\\
{{q_{12}(\alpha_{i})}}&{{q_{22}(\alpha_{i})}}&{{q_{23}(\alpha_{i})}}\\\
{{q_{13}(\alpha_{i})}}&{{q_{23}(\alpha_{i})}}&{{q_{33}(\alpha_{i})}}\end{array}}\right],i=1,2,3\end{array}\\\
$ (15)

keeps changing over time, where

$\sigma_{i}^{2}=\left\\{{\begin{array}[]{*{20}{c}}{{\frac{4-\pi}{\pi}}{{[{a_{{\rm{max}},i}}-\bar{a}_{i}(t)]}^{2}}{\rm{}}}\\\
{{\frac{4-\pi}{\pi}}{{[-{a_{{\rm{max}},i}}-\bar{a}_{i}(t)]}^{2}}}\end{array}}\right.\begin{array}[]{*{20}{l}}{}&{\bar{a}_{i}(t)>0}\\\
{}&{\bar{a}_{i}(t)<0}\end{array},\\\ $ (16)

and

$\begin{array}[]{l}q_{11}=\frac{1}{{2{\alpha^{5}}}}\left[{1-{e^{-2\alpha\Delta
t}}+2\alpha\Delta
t+\frac{2}{3}{\alpha^{3}}\Delta{t^{3}}-2{\alpha^{2}}\Delta{t^{2}}-4\alpha\Delta
t{e^{-\alpha\Delta t}}}\right]\\\
q_{12}=\frac{1}{{2{\alpha^{4}}}}\left[{{e^{-2\alpha\Delta
t}}+1-2{e^{-\alpha\Delta t}}+2\alpha\Delta t{e^{-\alpha\Delta
t}}-2\alpha\Delta t+{\alpha^{2}}\Delta{t^{2}}}\right]\\\
q_{13}=\frac{1}{{2{\alpha^{3}}}}\left[{1-{e^{-2\alpha\Delta t}}-2\alpha\Delta
t{e^{-\alpha\Delta t}}}\right]\\\
q_{22}=\frac{1}{{2{\alpha^{3}}}}\left[{4{e^{-\alpha\Delta
t}}-3-{e^{-2\alpha\Delta t}}+2\alpha\Delta t}\right]\\\
q_{23}=\frac{1}{{2{\alpha^{2}}}}\left[{{e^{-2\alpha\Delta
t}}+1-2{e^{-\alpha\Delta t}}}\right]\\\
q_{33}=\frac{1}{{2\alpha}}\left[{1-{e^{-2\alpha\Delta t}}}\right].\end{array}$
(17)

The currect statistical model parameters $\alpha_{i}$ is the reciprocal of the
maneuvering time constant, and $a_{{\rm{max}},i}$ is the maximum possible
acceleration.

####  2.4.2 Particle sampling

The particle sampling method of CSKPF method is almost the same with that of
CVPF method, except that the CSKFP method assumes that
$p\left({{\bf{x}}_{t}^{i}|{\bf{x}}_{t-1}^{i}}\right)$ obeys Gaussian
distribution with mean value $\mathbf{x}_{t|t-1}$ _but_ variance
$\mathbf{P}_{t|t-1})$, which means the particles are intensive when the
covariance is small and sparse when the covariance is large. In this way, the
particle sampling is more efficient.

####  2.4.3 State correction

Similar to the CVPF method, the state estimation $\hat{\mathbf{x}}_{t}$ and
its corresponding observation $\mathbf{y}_{t}$ can be obtained, but the CSKPF
loop does not end here. The state estimation and its covariance is then
corrected by observation $\mathbf{y}_{t}$ using Kalman filter, i.e.

$\mathbf{x}_{t}=\hat{\mathbf{x}}_{t}+\mathbf{K}_{t}\left(\mathbf{y}_{t}-\mathbf{H}\hat{\mathbf{x}}_{t}\right),$
(18)
$\mathbf{P}_{t}=\mathbf{P}_{t|t-1}-\mathbf{K}_{t}\mathbf{H}\mathbf{P}_{t|t-1},$
(19)

where

$\mathbf{K}_{t}=\mathbf{P}_{t|t-1}\mathbf{H}^{\top}\left(\mathbf{H}\mathbf{P}_{t|t-1}\mathbf{H}^{\top}+\mathbf{R}\right)^{-1}$
(20)

is the Kalman gain matrix, $\mathbf{H}$ is the observation matrix and
$\mathbf{R}$ is the observation noise matrix.

####  2.4.4 Notes

It should be noted that the current statistical model assumes that the
acceleration obeys the modified Rayleigh distribution whose mean value is the
prediction of the current acceleration. Therefore, the CSKPF method requires a
period of "warm-up" time to estimate the mean acceleration. In practice, at
first $T_{h}$ moment, a tracker uses both CVPF and CSKPF method to predict and
correct the state, but only accept the CVPF output. The mean acceleration is
updated at every CSKPF tracker step so that the estimation of acceleration
converges to a reasonable value after $T_{h}$ frames. When $t\geq T_{h}$, only
CSKPF runs.

As can be seen, the CSKPF can not only build the trajectories of objects but
also estimate their velocity and acceleration at the same time, which are
useful in motion analysis.

##  3 Experiments

In this section, the CVPF and CSKPF method are compared based on the proposed
evaluation index on simulated data and verified on the fruitfly data in the
real world.

###  3.1 Evaluation index

The tracking accuracy index MOTA and the tracking precision index MOTP is
widely used to evaluate multiple object tracking (MOT) method[14]. Based on
MOTA and MOTP, we propose more intuitively index which customizes for 3D
trajectories construction, i.e., tracking integrity, continuity and precision.

Suppose that the series
$\mathcal{P}^{i}=\\{\mathbf{p}_{t_{i}}^{i},\cdots,\mathbf{p}_{t}^{i},\cdots,\mathbf{p}_{T_{i}}^{i}\\}$
is the $i$-th true trajectory which starts at $t_{i}$ and end at $T_{i}$. The
$j$-th tracked trajectory is
$\mathcal{Q}^{j}=\\{\mathbf{q}_{t_{i}}^{j},\cdots,\mathbf{q}_{t}^{j},\cdots,\mathbf{q}_{T_{i}}^{j}\\}$.
Note that at some moment $t$, $\mathbf{q}_{t}^{j}$ may not exist due to miss
tracking. Define the position discrepancy between $\mathcal{P}^{i}$ and
$\mathcal{Q}^{j}$ at moment $t$ by

${d_{ij}}\left(t\right)=\left\\{{\begin{array}[]{*{20}{c}}{\left|{{\bf{p}}_{t}^{i}-{\bf{q}}_{t}^{j}}\right|}&{\mathrm{if}~{}{\bf{q}}_{t}^{j}{\rm{~{}exist}}}\\\
{{d_{0}}}&{{\rm{otherwize}}}\end{array}}\right.,i=1,\cdots,N,j=1,\cdots,M,$
(21)

where $d_{0}$ is the threshold distance. Then, the discrepancy between
$\mathcal{P}^{i}$ and $\mathcal{Q}^{j}$ among the whole tracking is

$d_{ij}=\frac{1}{(T_{i}-t_{i})}\sum_{t=t_{i}}^{T_{i}}d_{ij}(t).$ (22)

Integrity. For the $i$-th true trajectory, more integral tracking means it can
be tracked as long as possible. Therefore, define the tracking integrity as
possibility that the true trajectories tracked, i.e.,

$\mathrm{Integrity}=\frac{\sum_{i=1}^{N}{\sum_{t=t_{i}}^{T_{i}}{\delta\left(k_{t}^{i}\right)}}}{\sum_{i=1}^{N}{(T_{i}-t_{i})}},$
(23)

where $\delta\left(\cdot\right)$ is Dirac function.

Continuity. For the $i$-th true trajectory, more continuous tracking means its
matching trajectories are as less as possible. Using
IDSW$\left(\mathbf{k}^{i}\right)$ to count the switch times of matching series
$\mathbf{k}^{i}$, the tracking continuity can be represented by

$\mathrm{Continuity}=1-\frac{\sum_{i=1}^{N}{\mathrm{IDSW}\left(\mathbf{k}_{i}\right)}}{\sum_{i=1}^{N}{(T_{i}-t_{i})}}.$
(24)

Precision. For the $i$-th true trajectory, the tracking precision is defined
by the average discrepancy at $k_{t}^{i}\neq 0$ between true and tracked
trajectories, i.e.,

$\mathrm{Precision}=\frac{\sum_{i=1}^{N}{\sum_{t=t_{i}}^{T_{i}}{d_{ik_{t}^{i}(t)}}}}{\sum_{i=1}^{N}{\sum_{t=t_{i}}^{T_{i}}{\sum_{j=1}^{M}}}{b_{t,ij}}}.$
(25)

###  3.2 Simulation experiment

####  3.2.1 Generate true trajectories

The simulation data is generated in the coordinate system shown in Figure 3.
There are $N$ objects whose positions are randomly initialized in the box of
$x_{0},y_{0},z_{0}\in[-20,20]$. The $i$-th object’s speed is
$V_{i}(t)=6+2\sin\left(\frac{2\pi}{5}t+\alpha_{i,0}\right)$, and its heading
angle is
$\xi_{i}(t)=\frac{1}{2}A_{\xi,i}\left(1+\cos\left(\frac{\pi}{10}t+\xi_{i,0}\right)\right)$
and the climbing angle is
$\gamma_{i}(t)=\frac{1}{4}A_{\gamma,i}\cos\left(\frac{\pi}{10}t+\gamma_{i,0}\right)$
at moment $t$, where
$\alpha_{i,0},\xi_{i,0},\gamma_{i,0},\in[0,2\pi),A_{\xi,i},A_{\gamma,i}\in[-1,1]$
are randomly initialized. Then, the velocity at moment $t$ is

$\mathbf{v}_{i}(t)=V_{i}(t)\begin{bmatrix}\cos\gamma_{i}\cos\xi_{i}&\cos\gamma_{i}\sin\xi_{i}&\sin\gamma_{i}\end{bmatrix}^{\top}$
(26)

The true trajectories are integrated from their velocities.

Figure 2: Coordinate system and defination of heading angle $\xi$ and climbing
angle $\gamma$.

Figure 3: Example true value trajectories for $N=60$.

####  3.2.2 Render noised measurements

Suppose that the object entity is a ball with a radius of 0.5. As shown in
Figure 3, two cameras are orthogonally placed, and appropriate camera
parameters are selected to ensure that all objects are in the view field of
the two cameras. Gaussian noise is added to each object’s projection to
simulate the measurement error. An example is shown in Figure 3 and Figure 4.

Figure 4: Randered images in (a) view 1 and (b) view 2 for $N=60$ at moment
$t=1$. The blue regions are box and the black dots are objects.

####  3.2.3 Result

The performance of CVPF and CSKPF are compared when there are
$N=$1,20,40,60,80,100, 120, 140, 160 objects and repeats 50 times for each
$N$. The particle number of the two methods $N_{p}=100$. The particle
covariance $\sigma^{2}=0.09$ in CVPF method and $\alpha=5,a_{\mathrm{max}}=5$
in CSKPF method. The time step $\Delta t=0.1$s, and the tracking lasts 5
seconds.

Figure 5: The tracking (a) integrity, (b) continuity, (c) precision of CVPF
and CSKPF method in simulated data. The solid line represents the average
index of 50 experiments, and the corresponding shade surrounds one standard
deviation.

As shown in Figure 5, when there is only one object in the scene, the two
methods both successfully build integral and continuous trajectory because
there is no temporal-spatial association uncertainty for one object.
Nevertheless, the CSKPF method has better precision than the CVPF algorithm.
With the increase of objects, the temporal-spatial association uncertainties
increases. As a result, the integrity and continuity of the two methods both
show a downward trend, and the tracking errors increase. However, the
integrity and continuity deterioration of the CSKPF method is more moderate.
When there are 160 objects, the CSKPF integrity and continuity can still reach
about 0.85 and 0.995, respectively, while that of CVPF has dropped below 0.7
and 0.99. In terms of precision, the tracking errors of the two methods
increase linearly, but that of CSKPF is always lower than CVPF.

In summary, the simulation experiments prove that the proposed CSKPF method
can improve tracking integrity, continuity, and precision.

###  3.3 Fruitflies’ experiment

The tracking of the biological swarm is a challenging scenario for the multi-
view multi-object tracking method. The CVPF and CSKPF algorithm are further
compared on the flight data of fruitflies (_Drosophila melanogaster_) recorded
in a laboratory environment.

###  3.4 Data acquisition

Figure 6: image acquisition system.

As shown in the Figure 6, the fruitflies fly in a $400\times 400\times 400$mm
box glued by transparent acrylic board. The box is placed between infrared
high-speed digital cameras (PCO Dimax HS4 2000x2000@100FPS, Zeiss Biogon T*
35mm f/2 ZM) and the infrared light tabula which are covered with a diffuser
to produce soft and flicker-free light. The two calibrated cameras are
synchronized and orthogonally placed about 900 mm away from the center of the
box.

The measurements are expressed as pixel blobs
$\chi_{v,t,k},v=1,2,t=1,\cdots,T,k=1,\cdots,K$ segmented by Gaussian
background modeling method, whcih is the same with [10]. Other detection
method [15] can be used, but it is not discussed here.

####  3.4.1 Particle weight

Define the geometrical overlap ratio of the particle $\mathbf{x}^{i}_{t|t-1}$
and the measurement $\chi_{v,t,k_{v}}$ by

$\eta\left(\mathbf{x}_{t|t-1}^{i},\chi_{v,t,k_{v}}\right)=\frac{|\mathbf{P}^{v}\Upsilon(\mathbf{x}^{i}_{t|t-1})\bigcap\chi_{v,t,k_{v}}|}{|\chi^{i}_{v,t,k_{v}}|},v=1,2$
(27)

where the $|\cdot|$ operator counts the number of pixels, $\mathbf{P}^{v}$ is
the projection matrix of the camera $v$, $\Upsilon(\mathbf{x}^{i}_{t|t-1})$ is
a ball with the center $\mathbf{x}^{i}_{t|t-1}$. (The diameter of the ball is
set to 3mm because the average body length of fruitfly is about 3 mm). Denote
the most likely association at the moment of $t-1$ is
$\langle\chi_{1,t-1,k_{1}^{*}},\chi_{2,t-1,k_{2}^{*}}\rangle$, the credibility
of the $k$-th association ${\langle\chi_{1,t,k_{1}},\chi_{2,t,k_{2}}\rangle}$
on view $v$ is

$\tau_{v,k}\left(\mathbf{x}_{t|t-1}^{i}\right)=\eta\left(\mathbf{x}_{t|t-1}^{i},\chi_{v,t,k_{v}}\right)\mathrm{NCC}\left(\chi_{v,t,k_{v}},\chi_{v,t-1,k_{v}^{*}}\right)$
(28)

where NCC$(\cdot,\cdot)$ is the normalized covariance of the target appearence
and the measurement. The credibility of particle $\mathbf{x}_{t|t-1}^{i}$ on
view $v$ is $\tau_{v}=\mathop{\arg\max}\limits_{k}{\tau_{v,k}}$, and the
weight $w_{i}$ of the particle $\mathbf{x}^{i}_{t|t-1}$ can be updated by
Equation (4).

####  3.4.2 Result

The performance of CVPF and CSKPF are compared on the 100 frames data. The
particle number of the two methods $N_{p}=300$. The particle covariance
$\sigma^{2}=4$mm in CVPF method and $\alpha=1,a_{\mathrm{max}}=0.1$ in CSKPF
method. The time step $\Delta t=0.01$s.

Figure 7: Trajectories of Fruitfly reconstructed by (a) CVPF and (b) CVKPF
method.

As shown in Figure 7, the trajectories of CSKPF are smoother than that of
CVPF. Moreover, by counting the trajectories that are longer than $x$ in
Figure 9, it can be seen that there are more CSKPF trajectories than CVPF at
$30<x<50$, which means CSKPF can build a more integral trajectory. It should
be noted that at $50<x<60$, the integrity of the CVPF is slightly higher
because CVPF regards the trajectory that rebounded by the wall of the box as
one trajectory while the CSKPF does not accept these kinds of unsmooth
trajectories. In practice, only continuous trajectories are concerned. Thus,
the integrity of the CSKPF method is better than CVPF.

Figure 8: Trajectory integrity.

Figure 9: A typical trajectory. The triangle mark the starting point of the
trajectory, and the circles mark the end point of the trajectory.

One typical trajectory is shown in Figure 9 and it position and velocity curve
are shown in Figure 10 and Figure 11. Obviously, the CVPF trajectory is
unsmooth and even temporally deviates from the object. The deviation of CVPF
occurs around 30 frames according to the position curve. By reporjecting the
trajectory between the $25\sim 36$ frame to each view in Figure 12. It is
obvious that the deviation is caused by another fruitfly that flew towards the
target fruitfly in view 2 at 30 frames. The image of the two fruitflies is too
close for CVPF to distinguish while CSKPF successes thanks to better motion
prediction which is proved by the smoother and more reasonable velocity
prediction. On the other hand, better velocity estimation is important for
motion analysis.

In summary, the experiment of fruitfly flight proves that the proposed CSKPF
method can more effectively track a large amount of high kinematic objects
whose motion model is unknown.

Figure 10: Position curve in (a) x, (b) y, (c) z axis of the typical
trajectory.

Figure 11: Velocity curve in (a) x, (b) y, (c) z axis of the typical
trajectory.

Figure 12: Reprojection of the typical trajectory in (a) view 1 and (b) view
2.

##  4 Conclusion

Aiming at the multi-view multi-object tracking problem, this paper proposes
the CSKPF method which improves the existing CVPF method by introducing the
current statistical model to predict the motion, using the state covariance
matrix to resample particles and using the Kalman filter to suppress the
measurement noise. Both simulation experiments and real-world fruitfly
experiments prove the advantages of CSKPF over the baseline CVPF. The proposed
method is also suitable for tracking other highly kinematic objects even the
kinematic model is unknown. In future research, the missed assumption can be
introduced to deal with the trajectory interruption caused by the missed
detection, and other fast algorithms to solve the Bayesian inference problem
like unscented Kalman filter can be applied to improve the calculation
efficiency.

## Acknowledgments

None

## References

  * [1] Alessandro Attanasi, Andrea Cavagna, Lorenzo Del Castello, Irene Giardina, Stefania Melillo, Leonardo Parisi, Oliver Pohl, Bruno Rossaro, Edward Shen, Edmondo Silvestri, and Massimiliano Viale.  Collective behaviour without collective order in wild swarms of midges.  PLoS Computational Biology, 10(7):e1003697, July 2014.
  * [2] Julia Watzek, Mark E. Hauber, Katharine M. Jack, Julie R. Murrell, Stacey R. Tecot, and Sarah F. Brosnan.  Modelling collective decision-making: Insights into collective anti-predator behaviors from an agent-based approach.  Behavioural Processes, 193, 2021.
  * [3] Delia S. Shelton, Sierra G. Shelton, Danita K. Daniel, Manickam Raja, Anuradha Bhat, Robyn L. Tanguay, Dennis M. Higgs, and Emília P. Martins.  Collective Behavior in Wild Zebrafish.  Zebrafish, 17(4):243–252, 2020.
  * [4] Jae Woo Kim, Seong Joon Bae, Seongjin Park, and Do Hyung Kim.  Three-dimensional object tracking with occlusion using plenoptic image sequences.  Optical Engineering, 57(6):061604, 2018.
  * [5] JinLong Shi, Danping Zou, Suqin Bai, Qiang Qian, and Linbin Pang.  Reconstruction of dense three-dimensional shapes for outdoor scenes from an image sequence.  Optical Engineering, 52(12):123104, 2013.
  * [6] I.J. Cox and S.L. Hingorani.  An efficient implementation and evaluation of reid’s multiple hypothesis tracking algorithm for visual tracking.  In Proceedings of 12th International Conference on Pattern Recognition, volume 1, pages 437–442, Jerusalem, Israel, 1994. IEEE Comput. Soc. Press.
  * [7] Hao Du, Danping Zou, and Yan Qiu Chen.  Relative epipolar motion of tracked features for correspondence in binocular stereo.  In 2007 IEEE 11th International Conference on Computer Vision, pages 1–8, Rio de Janeiro, Brazil, 2007. IEEE.
  * [8] Hai Shan Wu, Qi Zhao, Danping Zou, and Yan Qiu Chen.  Automated 3d trajectory measuring of large numbers of moving particles.  Optics Express, 19(8):7646, 2011.
  * [9] Reza Ardekani, Anurag Biyani, Justin E. Dalton, Julia B. Saltz, Michelle N. Arbeitman, John Tower, Sergey Nuzhdin, and Simon Tavaré.  Three-dimensional tracking and behaviour monitoring of multiple fruit flies.  Journal of The Royal Society Interface, 10(78):20120547, 2013.
  * [10] Xi En Cheng, Zhi-Ming Qian, Shuo Hong Wang, Nan Jiang, Aike Guo, and Yan Qiu Chen.  A novel method for tracking individuals of fruit fly swarms flying in a laboratory flight arena.  PLOS ONE, 10(6), June 2015.
  * [11] Shuo Hong Wang, Hai Feng Su, Xi En Cheng, Ye Liu, Aike Quo, and Yan Qiu Chen.  Tracking the 3D position and orientation of flying swarms with learned kinematic pattern using LSTM network.  In 2017 IEEE International Conference on Multimedia and Expo (ICME), pages 1225–1230, Hong Kong, Hong Kong, 2017. IEEE.
  * [12] Hongren Zhou, Zhongliang Jing, and Peide Wang.  Maneuvering Target Tracking.  National Defence Industry Press, Beijing, 1991.
  * [13] James V. Candy.  Bayesian Signal Processing: Classical, Modern, and Particle Filtering Methods.  John Wiley & Sons, Inc., Hoboken, NJ, USA, 2016.
  * [14] Laura Leal-Taixe, Anton Milan, Ian Reid, Stefan Roth, and Konrad Schindler.  Motchallenge 2015: Towards a benchmark for multi-target tracking.  arXiv:1504.01942v1, 2015.
  * [15] Jingneng Fu, Honggang Wei, Hui Zhang, and Xiaodong Gao.  Three-dimensional pipeline hough transform for small target detection.  Optical Engineering, 60(2):1 – 15, 2021.