# A Quantitative Analysis on Two RFS-Based Filtering Methods for Multicell Tracking

Authors: Yayun Ren; Benlian Xu
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Multiobject filters developed from the theory of random finite sets (RFS) have recently become well-known methods for solving multiobject tracking problem. In this paper, we present two RFS-based filtering methods, Gaussian mixture probability hypothesis density (GM-PHD) filter and multi-Bernoulli filter, to quantitatively analyze their performance on tracking multiple cells in a series of low-contrast image sequences. The GM-PHD filter, under linear Gaussian assumptions on the cell dynamics and birth process, applies the PHD recursion to propagate the posterior intensity in an analytic form, while the multi-Bernoulli filter estimates the multitarget posterior density through propagating the parameters of a multi-Bernoulli RFS that approximates the posterior density of multitarget RFS. Numerous performance comparisons between the two RFS-based methods are carried out on two real cell images sequences and demonstrate that both yield satisfactory results that are in good agreement with manual tracking method.


---

## Body

## 1. Introduction

In previous related research in the field of medicine and biology, people often use manual visual inspection to observe living cells behaviors, such as the moving velocity and the density of cell population. Manual analysis methods always are limited by many challenges including large volume of image data and tedious manual operations. Therefore, it is of great significance to find an automatic and reliable way to track multiple cells. Over the last decades, many automatic tracking methods have been developed with rapid development of processing and computer vision technologies. Furthermore, numerous applications have been found for analyzing time-lapse cell microscopy imagery including fluorescent image [1, 2] and nonstaining cell image [3].In the field of cell tracking, automatic and accurate tracking needs to overcome problems that mainly stem from two factors: cell factors such as population size, mitosis, and shape variability and imaging factors such as level of noise, image resolution ratio, and image contrast. The main objective of cell motion analysis is to determine the number of cells and the state of each cell, and visual multicell tracking methods can be divided into two categories, deterministic methods and stochastic methods. Deterministic methods usually handle the detections and tracking tasks separately [4–7], but tracking may fail under problematic imaging conditions such as large cell density, cell division events, or segmentation errors.Stochastic methods mostly fall in the category of Bayesian framework, and tracking based on Bayesian probabilistic framework methods has also been proposed in recent years and is more robust to low resolution and signal-to-noise (SNR) scenarios than other tracking methods. Successful instances include multi-cell tracking in low contrast image sequences where the number of cells varies over time, cell shape deformation, and partial overlap [8, 9].Recently, Bayesian filtering based on Random Finite Set (RFS) theory has been proposed to solve multitarget tracking problems [10]. In the RFS formulation, the sets of targets and observations are modeled as random finite sets. This allows the problem of dynamically estimating multiple targets in presence of clutter and association uncertainty to be cast as a single-object filtering problem in a Bayesian filtering framework, where the object is the random set [10]. It is well known that this type of filter gives the optimal solution to multitarget tracking problems but is computationally intractable due to set integral calculation involved. The probabilistic hypothesis density (PHD) filter proposed by Mahler [10] overcomes this obstacle and significantly reduces the computational complexity. It is noted that the PHD recursion still involves integrals and does not have a closed-form solution. The implementation of the recursion is usually based on approximations for single-object distributions which resort to either sequential Monte Carlo method yielding sequential Monte-Carlo PHD (SMC-PHD) filter [11, 12] or the Gaussian mixture PHD (GM-PHD) filter [13, 14]. There are some reports on multitarget tracking in visual data using PHD filter. Examples include tracking of pedestrians in image sequences [15], face tracking [16], and vehicle tracking [17]. In these works, the SMC-PHD employs additional heuristics for partitioning particles to extract final state estimates and for target initialization.The GM-PHD filter is first applied to track multiple cells in [13] which focuses on the study of cells genealogy and find that the filter could track the lineage and motion of cells relatively well. The GM-PHD filter also exhibits good performance in maintaining continuous tracks for the cells in motion, even in low contrast image [18]. A new closed-form recursion for the LGJMS-PHD filter has been presented in [19] by incorporating state-dependent transition probability and the spawned transition probability. The filter has noticeably lower processing time in the cases of numerous cells and noisy detections. Recently, a sequential Monte Carlo (SMC) implementation of multi-Bernoulli filter is presented in [8], called “track-before-detect” technique. This method by passes the detection module and makes use of the spatiotemporal information directly extracted from poor image sequences. It demonstrates very good performance in cell tracking application. Stochastic methods generally make better use of the spatiotemporal information than deterministic methods and yield more robust tracking outcomes especially for poor cell image data.These cell tracking methods introduced above have their own properties and show different performance during implementation. This paper first quantitatively analyzes and compares the performance of GM-PHD and multi-Bernoulli filters [8, 9] in tracking multiple cells in a series of low-contrast image sequences. The cell states including cell number, position, velocity, and acceleration are estimated by both tracking methods, and the differences of obtained results are quantitatively analyzed by performance measures such as RMSE and OSPA distance. Furthermore, the detection and association errors and processing time are further discussed to compare the overall performance of the two tracking methods.This paper is organized as follows. In Section2, we briefly review the random finite set model in multitarget tracking problem. In Section 3, we present the principles and adapt the formulation of GM-PHD filter and multi-Bernoulli filter for multi-cell tracking. In Section 4, experimental results on real cell image sequences are presented to quantitatively analyze and compare the performance of the two methods. Finally, conclusions are summarized in Section 5.


## 2. Background on RFS Model

In a single target scenario, the state and measurement at any time are two vectors of possibly different dimensions, and these vectors evolve in time, with a constant dimension of each vector. However, in a multitarget scenario, some targets may disappear, new targets may appear, and the surviving targets may evolve to new states according to their individual dynamics. Obviously, the number and the state of targets and measurements may vary with time. In other words, the dimensions of the multitarget state and multitarget measurement evolve in time, and there is no ordering for the elements of multitarget state and measurement. Thus, such uncertainty in a multitarget system is characterized by modeling the multitarget state and multitarget measurement as random finite sets (RFS)Ξ
k and Σ
k on the (single target) state and observation spaces E
s and E
O, respectively. The RFS Ξ
k encapsulates all aspects of multitarget motion such as the time-varying number of targets, individual target motion, target birth, spawning, and target interactions. Given a realization X
k
-
1 of the RFS Ξ
k
-
1 at time k
-
1, the multitarget state at time k is modeled by the RFS
(1)
Ξ
k
=
S
k
(
Χ
k
-
1
)
∪
B
k
(
Χ
k
-
1
)
∪
Γ
k
,
where S
k
(
Χ
k
-
1
) denotes the RFS of targets that have survived from Χ
k
-
1 according to the dynamics at time k, B
k
(
Χ
k
-
1
) denotes the RFS of targets spawned from Χ
k
-
1 at time k, and Γ
k is the RFS of targets that appear spontaneously at time k. Similarly, the RFS Σ
k encapsulates all sensor characteristics such as measurement noise, sensor field of view, and false alarms. Given a realization X
k of Ξ
k at time k, the multitarget measurement is modeled by the RFS
(2)
Σ
k
=
Θ
k
(
X
k
)
∪
C
k
(
X
k
)
,
where Θ
k
(
X
k
) denotes the RFS of measurements generated by X
k and C
k
(
X
k
) denotes the RFS of clutter or false alarms.The problem of multitarget filtering concerns the estimation of the multitarget stateX
k at time k given the set Z
1
:
k
≡
(
Z
1
,
…
,
Z
k
) of all multitarget measurements up to time k. Let f
k
∣
k
-
1
(
·
∣
·
) be the multitarget transition density based on the dynamic model, and let g
k
(
·
∣
·
) be the multitarget likelihood. Supposing that the multitarget posterior density p
k
-
1
∣
k
-
1
(
X
k
-
1
∣
Z
1
:
k
-
1
) is known, the posterior density at time k can be calculated recursively. The optimal multitarget Bayes filter based on RFS theory can be described by the recursion, as in [10],
(3)
p
k
∣
k
-
1
(
X
k
∣
Z
1
:
k
-
1
)
=
∫
f
k
∣
k
-
1
(
X
k
∣
W
)
p
k
-
1
∣
k
-
1
(
W
∣
Z
1
:
k
-
1
)
δ
W
,
(4)
p
k
∣
k
(
X
k
∣
Z
1
:
k
)
=
p
k
∣
k
-
1
(
X
k
∣
Z
1
:
k
-
1
)
g
k
(
Z
k
∣
X
k
)
∫
p
k
∣
k
-
1
(
W
∣
Z
1
:
k
-
1
)
g
k
(
Z
k
∣
W
)
δ
W
,
where p
k
∣
k
-
1
(
X
k
∣
Z
1
:
k
-
1
) is the predicted value of multitarget posterior density at time k and δ is an appropriate reference measure on the state space.


## 3. Methods

This section provides the adaptions of the GM-PHD filter and multi-Bernoulli filter for solving multi-cell tracking problem, and the corresponding implementation issues are discussed in details as well.


### 3.1. The GM-PHD Filter for Multi-Cell Tracking

The optimal multitarget Bayes filtering recursions involve multiple integrals on multitarget state space and they are computationally intractable. The probability hypothesis density (PHD) filter, proposed by Mahler [20], helps to overcome this problem. The first statistical moment of the probability density function p
k
∣
k
(
X
k
∣
Z
1
:
k
), or PHD, is denoted by υ
k
∣
k
(
x
), and it has a useful property that ∫
S
υ
k
∣
k
(
x
)
d
x is equal to the expected number of targets in S. Hence, the PHD is not a probability density and does not necessarily sum to 1 but may be interpreted as a density of expected targets. The PHD filter recursion can be formulated as the prediction equation and the update equation and represented by
(5)
υ
k
∣
k
-
1
(
x
)
=
γ
k
(
x
)
+
∫
β
k
∣
k
-
1
(
x
∣
ζ
)
υ
k
-
1
(
ζ
)
+
∫
p
S
,
k
(
ζ
)
f
k
∣
k
-
1
(
x
∣
ζ
)
,
υ
k
(
x
)
=
[
1
-
p
D
,
k
(
x
)
]
υ
k
∣
k
-
1
(
x
)
+
∑
z
∈
Z
k
p
D
,
k
(
x
)
g
k
(
z
∣
x
)
υ
k
∣
k
-
1
(
x
)
Κ
k
(
z
)
+
∫
p
D
,
k
(
ξ
)
g
k
(
z
∣
ξ
)
υ
k
∣
k
-
1
(
ξ
)
d
ξ
,
where p
S
,
k
(
·
) is the probability that a target still exists at time k, p
D
,
k
(
·
) is the probability of detection at time k, and γ
k
(
·
), β
k
∣
k
-
1
(
·
), and K
k
(
·
) denote the intensities of birth, spawn and clutter, respectively. It is noted that the PHD recursion does not admit closed-form solutions in general, and numerical integration suffers from the curse of dimensionality. The Gaussian mixture probability hypothesis density (GM-PHD) filter [14] provides a closed-form solution to the PHD filter under the assumption of linear Gaussian multitarget dynamic and measurement models and approximates the PHD and other intensities with Gaussian mixture.In this section, we will adapt the original GM-PHD filter to solve multi-cell tracking problem, and many challenges are involved. These challenges include varying cell population densities due to cells leaving or entering the field of view and complex cellular topologies such as shape deformation, close contact, partial overlap, and uneven motions of cells. To present the GM-PHD filtering method for solving multi-cell tracking problem, we consider the following assumptions first.(
A
.
1
)
Each cell evolves and generates observations independently of one another.
(
A
.
2
)
Clutter is Poisson and independent of cell-originated measurements.
(
A
.
3
)
The survival and detection probabilities of each cell are state independent; that is,(6)
p
S
,
k
(
x
)
=
p
S
,
k
,
p
D
,
k
(
x
)
=
p
D
,
k
.
(
A
.
4
)
The predicted multi-cell state RFS governed byυ
k
∣
k
-
1 is Poisson.
(
A
.
5
)
Each cell follows a linear Gaussian dynamical model and the cell-originated measurement model has a linear Gaussian form; that is,(7)
f
k
∣
k
-
1
(
x
∣
ζ
)
=
N
(
x
;
F
k
-
1
ζ
,
Q
k
-
1
)
,
g
k
(
z
∣
x
)
=
N
(
z
;
H
k
x
,
R
k
)
,
whereN
(
·
;
m
,
P
) denotes a Gaussian density with mean m and covariance P, F
k
-
1 is the state transition matrix, H
k is the measurement matrix, and Q
k
-
1 and R
k are the process noise covariance and the observation noise covariance, respectively.
(
A
.
6
)
The intensities of cell birth and spawn RFS are in the form of Gaussian mixtures as follows:(8)
γ
k
(
x
)
=
∑
i
=
1
J
γ
,
k
ω
γ
,
k
(
i
)
N
(
x
;
m
γ
,
k
(
i
)
,
P
γ
,
k
(
i
)
)
,
(9)
β
k
∣
k
-
1
(
x
k
∣
ζ
)
=
∑
j
=
1
J
β
,
k
ω
β
,
k
(
j
)
N
(
x
k
;
F
β
,
k
-
1
(
j
)
ζ
+
d
β
,
k
-
1
(
j
)
,
Q
β
,
k
-
1
(
j
)
)
,where J
γ
,
k, ω
γ
,
k
(
i
), m
γ
,
k
(
i
), and P
γ
,
k
(
i
) are given model parameters that determine the shape of the cell birth intensity. J
β
,
k, ω
β
,
k
(
j
), F
β
,
k
-
1
(
j
), d
β
,
k
-
1
(
j
), and Q
β
,
k
-
1
(
j
) determine the shape of the spawning intensity of a cell with previous state ζ. For simplicity, we do not consider the cell mitosis, and the corresponding spawn term β
k
∣
k
-
1
(
x
k
∣
ζ
) is omitted here.Prediction Step. Suppose that assumptions (A.1)–(A.6) hold and that the posterior intensity of multiple cells at time k
-
1 is a Gaussian mixture of the form
(10)
υ
k
-
1
(
x
)
=
∑
i
=
1
J
k
-
1
ω
k
-
1
(
i
)
N
(
x
;
m
k
-
1
(
i
)
,
P
k
-
1
(
i
)
)
.
Then, the predicted intensity for time k is also a Gaussian mixture [14] and is given by
(11)
υ
k
∣
k
-
1
(
x
)
=
γ
k
(
x
)
+
υ
S
,
k
∣
k
-
1
(
x
)
,
where γ
k
(
x
) is given in (8) due to the spontaneous births such as new cells entering the field of view and υ
S
,
k
∣
k
-
1
(
x
) is the Gaussian mixture derived from cells that survived from the previous time step,
(12)
υ
S
,
k
∣
k
-
1
(
x
)
=
p
S
,
k
∑
j
=
1
J
k
-
1
ω
k
-
1
(
j
)
N
(
x
;
m
S
,
k
∣
k
-
1
(
j
)
,
P
S
,
k
∣
k
-
1
(
j
)
)
,
m
S
,
k
∣
k
-
1
(
j
)
=
F
k
-
1
m
k
-
1
(
j
)
,
P
S
,
k
∣
k
-
1
(
j
)
=
Q
k
-
1
+
F
k
-
1
P
k
-
1
(
j
)
F
k
-
1
T
.Update Step. Suppose that assumptions (A.1)–(A.6) hold and that the predicted intensity of multi-cell at time k is a Gaussian mixture of the form
(13)
υ
k
∣
k
-
1
(
x
)
=
∑
i
=
1
J
k
∣
k
-
1
ω
k
∣
k
-
1
(
i
)
N
(
x
;
m
k
∣
k
-
1
(
i
)
,
P
k
∣
k
-
1
(
i
)
)
.
Then, the posterior intensity at time k is also a Gaussian mixture, which consists of a misdetection term and a detection term and is given by
(14)
υ
k
(
x
)
=
(
1
-
p
D
,
k
)
υ
k
∣
k
-
1
(
x
)
+
∑
z
∈
Z
k
υ
D
,
k
(
x
;
z
)
,
where
(15)
υ
D
,
k
(
x
;
z
)
=
∑
j
=
1
J
k
∣
k
-
1
ω
k
(
j
)
(
z
)
N
(
x
;
m
k
∣
k
(
j
)
(
z
)
,
P
k
∣
k
(
j
)
)
,
ω
k
(
j
)
(
z
)
=
p
D
,
k
ω
k
∣
k
-
1
(
j
)
q
k
(
j
)
(
z
)
K
k
(
z
)
+
p
D
,
k
∑
i
=
1
J
k
∣
k
-
1
ω
k
∣
k
-
1
(
i
)
q
k
(
i
)
(
z
)
,
q
k
(
j
)
(
z
)
=
N
(
z
;
H
k
m
k
∣
k
-
1
(
j
)
,
R
k
+
H
k
P
k
∣
k
-
1
(
j
)
H
k
T
)
,
m
k
∣
k
-
1
(
j
)
(
z
)
=
m
k
∣
k
-
1
(
j
)
+
K
k
(
j
)
(
z
-
H
k
m
k
∣
k
-
1
(
j
)
)
,
P
k
∣
k
(
j
)
=
[
I
-
K
k
(
j
)
H
k
]
P
k
∣
k
-
1
(
j
)
,
K
k
(
j
)
=
P
k
∣
k
-
1
(
j
)
H
k
T
(
H
k
P
k
∣
k
-
1
(
j
)
H
k
T
+
R
k
)
-
1
.
Given the Gaussian mixture intensities υ
k
∣
k
-
1
(
x
) and υ
k
(
x
), the corresponding expected number of cells N
^
k
∣
k
-
1 and N
^
k can be obtained by summing up the appropriate weights. Specifically, the mean of the predicted number of cells N
^
k
∣
k
-
1 is obtained by adding the mean number of birth cells and the mean number of surviving cells based on (11)
(16)
N
^
k
∣
k
-
1
=
∑
i
=
1
J
γ
,
k
ω
γ
,
k
(
i
)
+
p
S
,
k
N
^
k
-
1
.
Similarly, the mean of the updated number of cells N
^
k is the sum of the number of measurements and the mean number of undetected cells based on (14)
(17)
N
^
k
=
(
1
-
p
D
,
k
)
N
^
k
∣
k
-
1
+
∑
z
∈
Z
k
∑
j
=
1
J
k
∣
k
-
1
ω
k
(
j
)
(
z
)
.
It can be observed that the Gaussian mixture PHD filter suffers from computation problems due to the increasing number of Gaussian components with the evolution of the process over time. It needs pruning and merging procedures to reduce the amount of calculation and improve the efficiency of tracking. For completeness, we summarize the main steps of the GM-PHD filter in Appendix A.Multicell State Extraction. In the GM-PHD filter, unlike the traditional PHD filter where the clustering process is required to extract multiobject state and usually leads to unreliable results for those objects in the proximity, the state and the number of cells can be extracted according to the Gaussian mixture representation of the posterior intensity. A good approximation to the Gaussian mixture posterior intensity
(18)
υ
k
(
x
)
=
∑
i
=
1
J
k
ω
k
(
i
)
N
(
x
;
m
k
(
i
)
,
P
k
(
i
)
)
can be obtained after pruning and merging steps. The Gaussians with weak weights can be ignored, while the means of the Gaussians with greater weights than threshold T
E (ω
k
(
i
)
>
T
E) are selected as the output cells, and we can get the position and velocity of cell centroid in our experiments.


### 3.2. The Multi-Bernoulli Filter for Multi-Cell Tracking

The multi-Bernoulli filter recursion, proposed by Mahler [10], is a direct approximation to multitarget posterior density through propagating the parameters of a multi-Bernoulli RFS. When applied to the field of cell tracking, the multi-Bernoulli filter recursion updates a finite but time-varying number of hypothesized cell tracks, each characterized by the probability of existence and the probability density of the current hypothesized cell state.In addition to assumptions (A.1)–(A.3) in Section3.1, we further assume that(
A
.
7
)
cell births follow a multi-Bernoulli RFS independent of cell survivals.Prediction Step. Suppose that the assumptions presented above hold and that the posterior multi-cell density at time k
-
1 is a multi-Bernoulli of the form
(19)
π
k
-
1
=
{
(
r
k
-
1
(
i
)
,
p
k
-
1
(
i
)
(
·
)
)
}
i
=
1
M
k
-
1
,
where r
(
i
) and p
(
i
)
(
·
) denote the existence probability and the probability density of the ith independent Bernoulli RFS, respectively. Then, the predicted multi-cell intensity for time k is also a multi-Bernoulli [10] and is given by
(20)
π
k
∣
k
-
1
=
{
(
r
P
,
k
∣
k
-
1
(
i
)
,
p
P
,
k
∣
k
-
1
(
i
)
(
·
)
)
}
i
=
1
M
k
-
1
∪
{
(
r
Γ
,
k
(
i
)
,
p
Γ
,
k
(
i
)
(
·
)
)
}
i
=
1
M
Γ
,
k
,
where {
(
r
Γ
,
k
(
i
)
,
p
Γ
,
k
(
i
)
(
·
)
)
}
i
=
1
M
Γ
,
k are the multi-Bernoulli parameters of the RFS of spontaneously born cells and the parameters {
(
r
P
,
k
∣
k
-
1
(
i
)
,
p
P
,
k
∣
k
-
1
(
i
)
(
·
)
)
}
i
=
1
M
k
-
1 are given by
(21)
r
k
∣
k
-
1
(
i
)
=
r
k
-
1
(
i
)
〈
p
k
-
1
(
i
)
(
·
)
,
p
S
,
k
(
·
)
〉
,
p
k
∣
k
-
1
(
i
)
(
x
)
=
〈
f
k
∣
k
-
1
(
x
∣
·
)
,
p
k
-
1
(
i
)
(
·
)
p
S
,
k
(
·
)
〉
〈
p
k
-
1
(
i
)
(
·
)
,
p
S
,
k
(
·
)
〉
,
where 〈
f
1
(
·
)
,
f
2
(
·
)
〉 denotes the standard inner product ∫
f
1
(
x
)
,
f
2
(
x
)
d
x and f
k
∣
k
-
1
(
·
∣
ζ
) is the transition density at time k given previous state ζ. In essence, the multi-Bernoulli parameter set for the predicted multi-cell density π
k
∣
k
-
1 is formed by the union of the multi-Bernoulli parameter sets for the surviving cells (the first term in (20)) and cell births (the second term in (20)). The total number of predicted cells is M
k
∣
k
-
1
=
M
k
-
1
+
M
Γ
,
k.Update Step. Suppose that the assumptions mentioned above hold and that the predicted multi-cell density for k is a multi-Bernoulli of the form
(22)
π
k
∣
k
-
1
=
{
(
r
k
∣
k
-
1
(
i
)
,
p
k
∣
k
-
1
(
i
)
(
·
)
)
}
i
=
1
M
k
∣
k
-
1
.
Then, given by Bayes rule (4), the posterior multi-cell density can be approximated by a multi-Bernoulli with the parameters
(23)
r
k
(
i
)
=
r
k
∣
k
-
1
(
i
)
〈
p
k
∣
k
-
1
(
i
)
(
·
)
,
g
k
(
·
)
〉
1
-
r
k
∣
k
-
1
(
i
)
+
r
k
∣
k
-
1
(
i
)
〈
p
k
∣
k
-
1
(
i
)
(
·
)
,
g
k
(
·
)
〉
,
p
k
(
i
)
(
·
)
=
p
k
∣
k
-
1
(
i
)
(
·
)
g
k
(
·
)
〈
p
k
∣
k
-
1
(
i
)
(
·
)
,
g
k
(
·
)
〉
,
where g
k
(
·
) is likelihood function, which will be discussed in Section 3.3.The sequential Monte Carlo (SMC) implementation for the multi-Bernoulli filter is briefly presented in AppendixB. Analogous to the standard particle filter, degeneracy is inevitable [21]. We resample particles for each hypothesized track after the update step to reduce the effect of degeneracy. Furthermore, the number of multi-Bernoulli parameters will grow fast as time step evolves. It needs the pruning procedure to reduce the computational load. To reduce the number of particles, pruning of hypothesized tracks is performed by discarding those with existence probabilities below a threshold T
P (we used 0.01) at each time step. In order to estimate the correct number and states of cells, the hypothesized cells with substantial overlap should be merged.Multicell State Extraction. According to the parameters of multi-Bernoulli RFS representing the multi-cell state after pruning and merging steps, the number of cells and their states can be estimated via finding the existence probabilities that are larger than a threshold (we used 0.5 in our experiments). In an alternative approach, the estimated number of cells is the cardinality mean or mode, and the cell individual state can be calculated by the weighted average of the particles of the corresponding density.


### 3.3. Implementation Issues

#### 3.3.1. Issues in GM-PHD Filter

As implicated in (14), the measurement values from each frame are required to compute the means, covariance, and weights of υ
k
(
x
) using (15). To calculate the measurement values, we propose a stable and robust multi-cell detection method in our studied data images. First, we use the grey level transformation to enhance the quality of the original image in order to facilitate separating cells from the image background. After the image preprocessing, we employ the image binarization processing. Since cell regions are darker than the background in our experiments, image inversion should be operated on the binary image. In order to make the contour of connected domain closer to the real cell, mathematical morphology processing, that is, close operation, is employed to compute the area, centroid, length, and width of each connected component. Finally, the output cell measurements can be calculated by ignoring the connected components with meaningless area and length-breadth ratio. As an example, Figure 1 shows the image processing results for a given real cell image using our proposed method and it can be observed that cell measurements, that is, the centroids of three connected components in Figure 1(f), can be effectively and accurately detected.Cell detection results. (a) Original cell image. (b) Grey level transformation. (c) Image binarization. (d) Image inversion. (e) Closed operation. (f) Output measurement values.
(a)
(b)
(c)
(d)
(e)
(f)Our experiments focus on maintaining continuous tracks for multi-cell in motion. Therefore, we assume that the cells follow linear Gaussian dynamics and measurement models as (7), respectively, and each state of cell is represented by
(24)
x
k
=
[
x
k
,
x
˙
k
,
y
k
,
y
˙
k
]
T
,
where (
x
k
,
y
k
) denotes the position of cell centroid and (
x
˙
k
,
y
˙
k
) denotes the corresponding cell velocity at time k. Furthermore, the motion and observation process can be modeled as
(25)
x
k
+
1
=
F
k
x
k
+
u
k
,
z
k
+
1
=
H
k
+
1
x
k
+
1
+
v
k
+
1
,
where we have
(26)
F
k
=
[
1
Δ
0
0
0
1
0
0
0
0
1
Δ
0
0
0
1
]
,
H
k
=
[
1
0
0
0
0
0
1
0
]
and Δ is the sampling period. Furthermore, u
k
~
N
(
0
,
Q
k
) and v
k
~
N
(
0
,
R
k
) are the process noise vector and the observation noise vector, respectively, and the covariance matrices Q
k and R
k of the process and measurement noise are given by
(27)
Q
k
=
[
B
0
2
0
2
B
]
σ
v
2
,
R
k
=
σ
w
2
I
2
,
where B
=
[
Δ
3
/
3
Δ
2
/
2
Δ
2
/
2
Δ
], 0
n and I
n denote the n
×
n zero and identity matrices, respectively, and σ
v and σ
w are the standard deviations of the process and measurement noise, respectively. The taken values of these parameters for multi-cell tracking are determined empirically according to the dynamics of cells (we set σ
v
=
1 pixel and σ
w
=
1.5 pixels in our experiments).


#### 3.3.2. Issues in Multi-Bernoulli Filter

The multi-Bernoulli filter, also called a “track-before-detect” method in [8], by passed the detection module and exploited the spatiotemporal information directly from the image sequence to estimate the multiobject state. The key step is to model the measurement model and to design the likelihood function g
k
(
·
∣
·
) in (4). Let X
=
{
x
1
,
…
,
x
n
} denote the multi-cell state vectors, and let y
k
=
[
y
1
,
k
,
y
2
,
k
,
…
,
y
m
,
k
] denote the image observation comprising an array of m pixels (or bins) value y
i
,
k of the ith pixel, which can be a real number or a vector depending on the application. If each cell or state is represented by a rectangular blob, then the corresponding histogram v
i based on the cell state x
i can be calculated (
i
=
1
,
…
,
n
). Based on the histogram difference comparison between the interested blob region and the cell training database, the similarity between the candidate and real cell can be estimated by (28) using total T
k cell training vectors {
v
j
*
}
j
=
1
T
k. We have
(28)
g
k
(
y
k
∣
x
i
)
=
(
e
-
(
1
-
d
(
v
i
)
)
γ
)
ξ
,
(29)
d
(
v
i
)
=
1
T
k
∑
j
=
1
T
k
∑
l
=
1
C
min
⁡
⁡
(
v
i
(
l
)
,
v
j
*
(
l
)
)
,
where γ and ξ are the adjustment coefficients designed for achieving better likelihood discrimination, v
i
(
l
) denotes the value of the lth element of the color histogram vector v
i, and C is the total number of elements in a histogram vector and takes the values of 256 × 3 for an RGB image.In addition, it is observed that the shape and brightness of cells might change over time. In order to reduce the negative effect of this interference, the cell training template should be updated in time. Specifically, assume that the cell templatet
k has total T
k vectors {
v
j
*
}
j
=
1
T
k at time step k, and if the similarity between the potential and template is greater than a threshold T
0 (i.e., g
k
(
t
k
∣
x
^
)
>
T
0), the potential x
^ is added to the cell template resulting as t
k
+
1
=
{
t
k
,
x
^
}. Obviously, the dimension of the cell template gradually increases with the update process, and we set T
max
⁡
=
10 to limit the dimension of cell template in our experiments for the sake of easy implementation and less computational burden.


#### 3.3.3. Label Management

The common approach of label management for the GM-PHD filter and multi-Bernoulli filter is briefly described as follows. Each Gaussian component or Bernoulli particle at timek simply inherits the label of its predecessor (in prediction step) if it belongs to one of the first predicted M
k
-
1 Gaussian components or Bernoulli particles; otherwise, it belongs to a newly born target and is assigned a new label, and the labels remain unchanged during the update step of the filter [10, 22]. However, this method works best when the targets are well separated. In fact, cell movement is random and changeable over time. When two cells occlude each other, they are easily merged in one frame and separated in another frame. As a result, one of them would be detected by the filter as a newly born cell, and would be assigned a new label as a new track. To solve the label management problem, we present a method developed from a graph theoretical framework in which the previously missed cells are taken into account besides the most recent cells, for their association with currently detected cells [10]. This label management scheme has a memory to keep the labels of missed cells for later use in case they reappear. The recorded labels will only be removed from memory if they are missed for several consecutive frames. This method employs a three-stage constrained nearest neighbor graph search to find candidate tracks in a bipartite graph. The search is constrained in the sense that an edge is chosen only if its cost is less than a certain limit [8].Initialization. The initialization takes place in the first frame. The estimates of locations of cells are read and labeled and play the role of previous cells for the next frame. After initialization, an iterative labeling contained three-stage constrained nearest neighbor graph search and memory update is repeated for each of the next frames.Stage 1. Some of the current estimated cells are matched to some of the estimates from the previous iteration.Stage 2. The remainder of the current estimates is matched to some previous cells that have been missed at least once.Stage 3. The remainder of the current estimates is checked if they are close enough to the boundaries of the image, and if validated, they are assigned new labels.Memory Update. The previous cells that have been missed for too long are removed from the memory, and other identity data is updated for the next iteration.


## 3.1. The GM-PHD Filter for Multi-Cell Tracking

The optimal multitarget Bayes filtering recursions involve multiple integrals on multitarget state space and they are computationally intractable. The probability hypothesis density (PHD) filter, proposed by Mahler [20], helps to overcome this problem. The first statistical moment of the probability density function p
k
∣
k
(
X
k
∣
Z
1
:
k
), or PHD, is denoted by υ
k
∣
k
(
x
), and it has a useful property that ∫
S
υ
k
∣
k
(
x
)
d
x is equal to the expected number of targets in S. Hence, the PHD is not a probability density and does not necessarily sum to 1 but may be interpreted as a density of expected targets. The PHD filter recursion can be formulated as the prediction equation and the update equation and represented by
(5)
υ
k
∣
k
-
1
(
x
)
=
γ
k
(
x
)
+
∫
β
k
∣
k
-
1
(
x
∣
ζ
)
υ
k
-
1
(
ζ
)
+
∫
p
S
,
k
(
ζ
)
f
k
∣
k
-
1
(
x
∣
ζ
)
,
υ
k
(
x
)
=
[
1
-
p
D
,
k
(
x
)
]
υ
k
∣
k
-
1
(
x
)
+
∑
z
∈
Z
k
p
D
,
k
(
x
)
g
k
(
z
∣
x
)
υ
k
∣
k
-
1
(
x
)
Κ
k
(
z
)
+
∫
p
D
,
k
(
ξ
)
g
k
(
z
∣
ξ
)
υ
k
∣
k
-
1
(
ξ
)
d
ξ
,
where p
S
,
k
(
·
) is the probability that a target still exists at time k, p
D
,
k
(
·
) is the probability of detection at time k, and γ
k
(
·
), β
k
∣
k
-
1
(
·
), and K
k
(
·
) denote the intensities of birth, spawn and clutter, respectively. It is noted that the PHD recursion does not admit closed-form solutions in general, and numerical integration suffers from the curse of dimensionality. The Gaussian mixture probability hypothesis density (GM-PHD) filter [14] provides a closed-form solution to the PHD filter under the assumption of linear Gaussian multitarget dynamic and measurement models and approximates the PHD and other intensities with Gaussian mixture.In this section, we will adapt the original GM-PHD filter to solve multi-cell tracking problem, and many challenges are involved. These challenges include varying cell population densities due to cells leaving or entering the field of view and complex cellular topologies such as shape deformation, close contact, partial overlap, and uneven motions of cells. To present the GM-PHD filtering method for solving multi-cell tracking problem, we consider the following assumptions first.(
A
.
1
)
Each cell evolves and generates observations independently of one another.
(
A
.
2
)
Clutter is Poisson and independent of cell-originated measurements.
(
A
.
3
)
The survival and detection probabilities of each cell are state independent; that is,(6)
p
S
,
k
(
x
)
=
p
S
,
k
,
p
D
,
k
(
x
)
=
p
D
,
k
.
(
A
.
4
)
The predicted multi-cell state RFS governed byυ
k
∣
k
-
1 is Poisson.
(
A
.
5
)
Each cell follows a linear Gaussian dynamical model and the cell-originated measurement model has a linear Gaussian form; that is,(7)
f
k
∣
k
-
1
(
x
∣
ζ
)
=
N
(
x
;
F
k
-
1
ζ
,
Q
k
-
1
)
,
g
k
(
z
∣
x
)
=
N
(
z
;
H
k
x
,
R
k
)
,
whereN
(
·
;
m
,
P
) denotes a Gaussian density with mean m and covariance P, F
k
-
1 is the state transition matrix, H
k is the measurement matrix, and Q
k
-
1 and R
k are the process noise covariance and the observation noise covariance, respectively.
(
A
.
6
)
The intensities of cell birth and spawn RFS are in the form of Gaussian mixtures as follows:(8)
γ
k
(
x
)
=
∑
i
=
1
J
γ
,
k
ω
γ
,
k
(
i
)
N
(
x
;
m
γ
,
k
(
i
)
,
P
γ
,
k
(
i
)
)
,
(9)
β
k
∣
k
-
1
(
x
k
∣
ζ
)
=
∑
j
=
1
J
β
,
k
ω
β
,
k
(
j
)
N
(
x
k
;
F
β
,
k
-
1
(
j
)
ζ
+
d
β
,
k
-
1
(
j
)
,
Q
β
,
k
-
1
(
j
)
)
,where J
γ
,
k, ω
γ
,
k
(
i
), m
γ
,
k
(
i
), and P
γ
,
k
(
i
) are given model parameters that determine the shape of the cell birth intensity. J
β
,
k, ω
β
,
k
(
j
), F
β
,
k
-
1
(
j
), d
β
,
k
-
1
(
j
), and Q
β
,
k
-
1
(
j
) determine the shape of the spawning intensity of a cell with previous state ζ. For simplicity, we do not consider the cell mitosis, and the corresponding spawn term β
k
∣
k
-
1
(
x
k
∣
ζ
) is omitted here.Prediction Step. Suppose that assumptions (A.1)–(A.6) hold and that the posterior intensity of multiple cells at time k
-
1 is a Gaussian mixture of the form
(10)
υ
k
-
1
(
x
)
=
∑
i
=
1
J
k
-
1
ω
k
-
1
(
i
)
N
(
x
;
m
k
-
1
(
i
)
,
P
k
-
1
(
i
)
)
.
Then, the predicted intensity for time k is also a Gaussian mixture [14] and is given by
(11)
υ
k
∣
k
-
1
(
x
)
=
γ
k
(
x
)
+
υ
S
,
k
∣
k
-
1
(
x
)
,
where γ
k
(
x
) is given in (8) due to the spontaneous births such as new cells entering the field of view and υ
S
,
k
∣
k
-
1
(
x
) is the Gaussian mixture derived from cells that survived from the previous time step,
(12)
υ
S
,
k
∣
k
-
1
(
x
)
=
p
S
,
k
∑
j
=
1
J
k
-
1
ω
k
-
1
(
j
)
N
(
x
;
m
S
,
k
∣
k
-
1
(
j
)
,
P
S
,
k
∣
k
-
1
(
j
)
)
,
m
S
,
k
∣
k
-
1
(
j
)
=
F
k
-
1
m
k
-
1
(
j
)
,
P
S
,
k
∣
k
-
1
(
j
)
=
Q
k
-
1
+
F
k
-
1
P
k
-
1
(
j
)
F
k
-
1
T
.Update Step. Suppose that assumptions (A.1)–(A.6) hold and that the predicted intensity of multi-cell at time k is a Gaussian mixture of the form
(13)
υ
k
∣
k
-
1
(
x
)
=
∑
i
=
1
J
k
∣
k
-
1
ω
k
∣
k
-
1
(
i
)
N
(
x
;
m
k
∣
k
-
1
(
i
)
,
P
k
∣
k
-
1
(
i
)
)
.
Then, the posterior intensity at time k is also a Gaussian mixture, which consists of a misdetection term and a detection term and is given by
(14)
υ
k
(
x
)
=
(
1
-
p
D
,
k
)
υ
k
∣
k
-
1
(
x
)
+
∑
z
∈
Z
k
υ
D
,
k
(
x
;
z
)
,
where
(15)
υ
D
,
k
(
x
;
z
)
=
∑
j
=
1
J
k
∣
k
-
1
ω
k
(
j
)
(
z
)
N
(
x
;
m
k
∣
k
(
j
)
(
z
)
,
P
k
∣
k
(
j
)
)
,
ω
k
(
j
)
(
z
)
=
p
D
,
k
ω
k
∣
k
-
1
(
j
)
q
k
(
j
)
(
z
)
K
k
(
z
)
+
p
D
,
k
∑
i
=
1
J
k
∣
k
-
1
ω
k
∣
k
-
1
(
i
)
q
k
(
i
)
(
z
)
,
q
k
(
j
)
(
z
)
=
N
(
z
;
H
k
m
k
∣
k
-
1
(
j
)
,
R
k
+
H
k
P
k
∣
k
-
1
(
j
)
H
k
T
)
,
m
k
∣
k
-
1
(
j
)
(
z
)
=
m
k
∣
k
-
1
(
j
)
+
K
k
(
j
)
(
z
-
H
k
m
k
∣
k
-
1
(
j
)
)
,
P
k
∣
k
(
j
)
=
[
I
-
K
k
(
j
)
H
k
]
P
k
∣
k
-
1
(
j
)
,
K
k
(
j
)
=
P
k
∣
k
-
1
(
j
)
H
k
T
(
H
k
P
k
∣
k
-
1
(
j
)
H
k
T
+
R
k
)
-
1
.
Given the Gaussian mixture intensities υ
k
∣
k
-
1
(
x
) and υ
k
(
x
), the corresponding expected number of cells N
^
k
∣
k
-
1 and N
^
k can be obtained by summing up the appropriate weights. Specifically, the mean of the predicted number of cells N
^
k
∣
k
-
1 is obtained by adding the mean number of birth cells and the mean number of surviving cells based on (11)
(16)
N
^
k
∣
k
-
1
=
∑
i
=
1
J
γ
,
k
ω
γ
,
k
(
i
)
+
p
S
,
k
N
^
k
-
1
.
Similarly, the mean of the updated number of cells N
^
k is the sum of the number of measurements and the mean number of undetected cells based on (14)
(17)
N
^
k
=
(
1
-
p
D
,
k
)
N
^
k
∣
k
-
1
+
∑
z
∈
Z
k
∑
j
=
1
J
k
∣
k
-
1
ω
k
(
j
)
(
z
)
.
It can be observed that the Gaussian mixture PHD filter suffers from computation problems due to the increasing number of Gaussian components with the evolution of the process over time. It needs pruning and merging procedures to reduce the amount of calculation and improve the efficiency of tracking. For completeness, we summarize the main steps of the GM-PHD filter in Appendix A.Multicell State Extraction. In the GM-PHD filter, unlike the traditional PHD filter where the clustering process is required to extract multiobject state and usually leads to unreliable results for those objects in the proximity, the state and the number of cells can be extracted according to the Gaussian mixture representation of the posterior intensity. A good approximation to the Gaussian mixture posterior intensity
(18)
υ
k
(
x
)
=
∑
i
=
1
J
k
ω
k
(
i
)
N
(
x
;
m
k
(
i
)
,
P
k
(
i
)
)
can be obtained after pruning and merging steps. The Gaussians with weak weights can be ignored, while the means of the Gaussians with greater weights than threshold T
E (ω
k
(
i
)
>
T
E) are selected as the output cells, and we can get the position and velocity of cell centroid in our experiments.


## 3.2. The Multi-Bernoulli Filter for Multi-Cell Tracking

The multi-Bernoulli filter recursion, proposed by Mahler [10], is a direct approximation to multitarget posterior density through propagating the parameters of a multi-Bernoulli RFS. When applied to the field of cell tracking, the multi-Bernoulli filter recursion updates a finite but time-varying number of hypothesized cell tracks, each characterized by the probability of existence and the probability density of the current hypothesized cell state.In addition to assumptions (A.1)–(A.3) in Section3.1, we further assume that(
A
.
7
)
cell births follow a multi-Bernoulli RFS independent of cell survivals.Prediction Step. Suppose that the assumptions presented above hold and that the posterior multi-cell density at time k
-
1 is a multi-Bernoulli of the form
(19)
π
k
-
1
=
{
(
r
k
-
1
(
i
)
,
p
k
-
1
(
i
)
(
·
)
)
}
i
=
1
M
k
-
1
,
where r
(
i
) and p
(
i
)
(
·
) denote the existence probability and the probability density of the ith independent Bernoulli RFS, respectively. Then, the predicted multi-cell intensity for time k is also a multi-Bernoulli [10] and is given by
(20)
π
k
∣
k
-
1
=
{
(
r
P
,
k
∣
k
-
1
(
i
)
,
p
P
,
k
∣
k
-
1
(
i
)
(
·
)
)
}
i
=
1
M
k
-
1
∪
{
(
r
Γ
,
k
(
i
)
,
p
Γ
,
k
(
i
)
(
·
)
)
}
i
=
1
M
Γ
,
k
,
where {
(
r
Γ
,
k
(
i
)
,
p
Γ
,
k
(
i
)
(
·
)
)
}
i
=
1
M
Γ
,
k are the multi-Bernoulli parameters of the RFS of spontaneously born cells and the parameters {
(
r
P
,
k
∣
k
-
1
(
i
)
,
p
P
,
k
∣
k
-
1
(
i
)
(
·
)
)
}
i
=
1
M
k
-
1 are given by
(21)
r
k
∣
k
-
1
(
i
)
=
r
k
-
1
(
i
)
〈
p
k
-
1
(
i
)
(
·
)
,
p
S
,
k
(
·
)
〉
,
p
k
∣
k
-
1
(
i
)
(
x
)
=
〈
f
k
∣
k
-
1
(
x
∣
·
)
,
p
k
-
1
(
i
)
(
·
)
p
S
,
k
(
·
)
〉
〈
p
k
-
1
(
i
)
(
·
)
,
p
S
,
k
(
·
)
〉
,
where 〈
f
1
(
·
)
,
f
2
(
·
)
〉 denotes the standard inner product ∫
f
1
(
x
)
,
f
2
(
x
)
d
x and f
k
∣
k
-
1
(
·
∣
ζ
) is the transition density at time k given previous state ζ. In essence, the multi-Bernoulli parameter set for the predicted multi-cell density π
k
∣
k
-
1 is formed by the union of the multi-Bernoulli parameter sets for the surviving cells (the first term in (20)) and cell births (the second term in (20)). The total number of predicted cells is M
k
∣
k
-
1
=
M
k
-
1
+
M
Γ
,
k.Update Step. Suppose that the assumptions mentioned above hold and that the predicted multi-cell density for k is a multi-Bernoulli of the form
(22)
π
k
∣
k
-
1
=
{
(
r
k
∣
k
-
1
(
i
)
,
p
k
∣
k
-
1
(
i
)
(
·
)
)
}
i
=
1
M
k
∣
k
-
1
.
Then, given by Bayes rule (4), the posterior multi-cell density can be approximated by a multi-Bernoulli with the parameters
(23)
r
k
(
i
)
=
r
k
∣
k
-
1
(
i
)
〈
p
k
∣
k
-
1
(
i
)
(
·
)
,
g
k
(
·
)
〉
1
-
r
k
∣
k
-
1
(
i
)
+
r
k
∣
k
-
1
(
i
)
〈
p
k
∣
k
-
1
(
i
)
(
·
)
,
g
k
(
·
)
〉
,
p
k
(
i
)
(
·
)
=
p
k
∣
k
-
1
(
i
)
(
·
)
g
k
(
·
)
〈
p
k
∣
k
-
1
(
i
)
(
·
)
,
g
k
(
·
)
〉
,
where g
k
(
·
) is likelihood function, which will be discussed in Section 3.3.The sequential Monte Carlo (SMC) implementation for the multi-Bernoulli filter is briefly presented in AppendixB. Analogous to the standard particle filter, degeneracy is inevitable [21]. We resample particles for each hypothesized track after the update step to reduce the effect of degeneracy. Furthermore, the number of multi-Bernoulli parameters will grow fast as time step evolves. It needs the pruning procedure to reduce the computational load. To reduce the number of particles, pruning of hypothesized tracks is performed by discarding those with existence probabilities below a threshold T
P (we used 0.01) at each time step. In order to estimate the correct number and states of cells, the hypothesized cells with substantial overlap should be merged.Multicell State Extraction. According to the parameters of multi-Bernoulli RFS representing the multi-cell state after pruning and merging steps, the number of cells and their states can be estimated via finding the existence probabilities that are larger than a threshold (we used 0.5 in our experiments). In an alternative approach, the estimated number of cells is the cardinality mean or mode, and the cell individual state can be calculated by the weighted average of the particles of the corresponding density.


## 3.3. Implementation Issues

### 3.3.1. Issues in GM-PHD Filter

As implicated in (14), the measurement values from each frame are required to compute the means, covariance, and weights of υ
k
(
x
) using (15). To calculate the measurement values, we propose a stable and robust multi-cell detection method in our studied data images. First, we use the grey level transformation to enhance the quality of the original image in order to facilitate separating cells from the image background. After the image preprocessing, we employ the image binarization processing. Since cell regions are darker than the background in our experiments, image inversion should be operated on the binary image. In order to make the contour of connected domain closer to the real cell, mathematical morphology processing, that is, close operation, is employed to compute the area, centroid, length, and width of each connected component. Finally, the output cell measurements can be calculated by ignoring the connected components with meaningless area and length-breadth ratio. As an example, Figure 1 shows the image processing results for a given real cell image using our proposed method and it can be observed that cell measurements, that is, the centroids of three connected components in Figure 1(f), can be effectively and accurately detected.Cell detection results. (a) Original cell image. (b) Grey level transformation. (c) Image binarization. (d) Image inversion. (e) Closed operation. (f) Output measurement values.
(a)
(b)
(c)
(d)
(e)
(f)Our experiments focus on maintaining continuous tracks for multi-cell in motion. Therefore, we assume that the cells follow linear Gaussian dynamics and measurement models as (7), respectively, and each state of cell is represented by
(24)
x
k
=
[
x
k
,
x
˙
k
,
y
k
,
y
˙
k
]
T
,
where (
x
k
,
y
k
) denotes the position of cell centroid and (
x
˙
k
,
y
˙
k
) denotes the corresponding cell velocity at time k. Furthermore, the motion and observation process can be modeled as
(25)
x
k
+
1
=
F
k
x
k
+
u
k
,
z
k
+
1
=
H
k
+
1
x
k
+
1
+
v
k
+
1
,
where we have
(26)
F
k
=
[
1
Δ
0
0
0
1
0
0
0
0
1
Δ
0
0
0
1
]
,
H
k
=
[
1
0
0
0
0
0
1
0
]
and Δ is the sampling period. Furthermore, u
k
~
N
(
0
,
Q
k
) and v
k
~
N
(
0
,
R
k
) are the process noise vector and the observation noise vector, respectively, and the covariance matrices Q
k and R
k of the process and measurement noise are given by
(27)
Q
k
=
[
B
0
2
0
2
B
]
σ
v
2
,
R
k
=
σ
w
2
I
2
,
where B
=
[
Δ
3
/
3
Δ
2
/
2
Δ
2
/
2
Δ
], 0
n and I
n denote the n
×
n zero and identity matrices, respectively, and σ
v and σ
w are the standard deviations of the process and measurement noise, respectively. The taken values of these parameters for multi-cell tracking are determined empirically according to the dynamics of cells (we set σ
v
=
1 pixel and σ
w
=
1.5 pixels in our experiments).


### 3.3.2. Issues in Multi-Bernoulli Filter

The multi-Bernoulli filter, also called a “track-before-detect” method in [8], by passed the detection module and exploited the spatiotemporal information directly from the image sequence to estimate the multiobject state. The key step is to model the measurement model and to design the likelihood function g
k
(
·
∣
·
) in (4). Let X
=
{
x
1
,
…
,
x
n
} denote the multi-cell state vectors, and let y
k
=
[
y
1
,
k
,
y
2
,
k
,
…
,
y
m
,
k
] denote the image observation comprising an array of m pixels (or bins) value y
i
,
k of the ith pixel, which can be a real number or a vector depending on the application. If each cell or state is represented by a rectangular blob, then the corresponding histogram v
i based on the cell state x
i can be calculated (
i
=
1
,
…
,
n
). Based on the histogram difference comparison between the interested blob region and the cell training database, the similarity between the candidate and real cell can be estimated by (28) using total T
k cell training vectors {
v
j
*
}
j
=
1
T
k. We have
(28)
g
k
(
y
k
∣
x
i
)
=
(
e
-
(
1
-
d
(
v
i
)
)
γ
)
ξ
,
(29)
d
(
v
i
)
=
1
T
k
∑
j
=
1
T
k
∑
l
=
1
C
min
⁡
⁡
(
v
i
(
l
)
,
v
j
*
(
l
)
)
,
where γ and ξ are the adjustment coefficients designed for achieving better likelihood discrimination, v
i
(
l
) denotes the value of the lth element of the color histogram vector v
i, and C is the total number of elements in a histogram vector and takes the values of 256 × 3 for an RGB image.In addition, it is observed that the shape and brightness of cells might change over time. In order to reduce the negative effect of this interference, the cell training template should be updated in time. Specifically, assume that the cell templatet
k has total T
k vectors {
v
j
*
}
j
=
1
T
k at time step k, and if the similarity between the potential and template is greater than a threshold T
0 (i.e., g
k
(
t
k
∣
x
^
)
>
T
0), the potential x
^ is added to the cell template resulting as t
k
+
1
=
{
t
k
,
x
^
}. Obviously, the dimension of the cell template gradually increases with the update process, and we set T
max
⁡
=
10 to limit the dimension of cell template in our experiments for the sake of easy implementation and less computational burden.


### 3.3.3. Label Management

The common approach of label management for the GM-PHD filter and multi-Bernoulli filter is briefly described as follows. Each Gaussian component or Bernoulli particle at timek simply inherits the label of its predecessor (in prediction step) if it belongs to one of the first predicted M
k
-
1 Gaussian components or Bernoulli particles; otherwise, it belongs to a newly born target and is assigned a new label, and the labels remain unchanged during the update step of the filter [10, 22]. However, this method works best when the targets are well separated. In fact, cell movement is random and changeable over time. When two cells occlude each other, they are easily merged in one frame and separated in another frame. As a result, one of them would be detected by the filter as a newly born cell, and would be assigned a new label as a new track. To solve the label management problem, we present a method developed from a graph theoretical framework in which the previously missed cells are taken into account besides the most recent cells, for their association with currently detected cells [10]. This label management scheme has a memory to keep the labels of missed cells for later use in case they reappear. The recorded labels will only be removed from memory if they are missed for several consecutive frames. This method employs a three-stage constrained nearest neighbor graph search to find candidate tracks in a bipartite graph. The search is constrained in the sense that an edge is chosen only if its cost is less than a certain limit [8].Initialization. The initialization takes place in the first frame. The estimates of locations of cells are read and labeled and play the role of previous cells for the next frame. After initialization, an iterative labeling contained three-stage constrained nearest neighbor graph search and memory update is repeated for each of the next frames.Stage 1. Some of the current estimated cells are matched to some of the estimates from the previous iteration.Stage 2. The remainder of the current estimates is matched to some previous cells that have been missed at least once.Stage 3. The remainder of the current estimates is checked if they are close enough to the boundaries of the image, and if validated, they are assigned new labels.Memory Update. The previous cells that have been missed for too long are removed from the memory, and other identity data is updated for the next iteration.


## 3.3.1. Issues in GM-PHD Filter

As implicated in (14), the measurement values from each frame are required to compute the means, covariance, and weights of υ
k
(
x
) using (15). To calculate the measurement values, we propose a stable and robust multi-cell detection method in our studied data images. First, we use the grey level transformation to enhance the quality of the original image in order to facilitate separating cells from the image background. After the image preprocessing, we employ the image binarization processing. Since cell regions are darker than the background in our experiments, image inversion should be operated on the binary image. In order to make the contour of connected domain closer to the real cell, mathematical morphology processing, that is, close operation, is employed to compute the area, centroid, length, and width of each connected component. Finally, the output cell measurements can be calculated by ignoring the connected components with meaningless area and length-breadth ratio. As an example, Figure 1 shows the image processing results for a given real cell image using our proposed method and it can be observed that cell measurements, that is, the centroids of three connected components in Figure 1(f), can be effectively and accurately detected.Cell detection results. (a) Original cell image. (b) Grey level transformation. (c) Image binarization. (d) Image inversion. (e) Closed operation. (f) Output measurement values.
(a)
(b)
(c)
(d)
(e)
(f)Our experiments focus on maintaining continuous tracks for multi-cell in motion. Therefore, we assume that the cells follow linear Gaussian dynamics and measurement models as (7), respectively, and each state of cell is represented by
(24)
x
k
=
[
x
k
,
x
˙
k
,
y
k
,
y
˙
k
]
T
,
where (
x
k
,
y
k
) denotes the position of cell centroid and (
x
˙
k
,
y
˙
k
) denotes the corresponding cell velocity at time k. Furthermore, the motion and observation process can be modeled as
(25)
x
k
+
1
=
F
k
x
k
+
u
k
,
z
k
+
1
=
H
k
+
1
x
k
+
1
+
v
k
+
1
,
where we have
(26)
F
k
=
[
1
Δ
0
0
0
1
0
0
0
0
1
Δ
0
0
0
1
]
,
H
k
=
[
1
0
0
0
0
0
1
0
]
and Δ is the sampling period. Furthermore, u
k
~
N
(
0
,
Q
k
) and v
k
~
N
(
0
,
R
k
) are the process noise vector and the observation noise vector, respectively, and the covariance matrices Q
k and R
k of the process and measurement noise are given by
(27)
Q
k
=
[
B
0
2
0
2
B
]
σ
v
2
,
R
k
=
σ
w
2
I
2
,
where B
=
[
Δ
3
/
3
Δ
2
/
2
Δ
2
/
2
Δ
], 0
n and I
n denote the n
×
n zero and identity matrices, respectively, and σ
v and σ
w are the standard deviations of the process and measurement noise, respectively. The taken values of these parameters for multi-cell tracking are determined empirically according to the dynamics of cells (we set σ
v
=
1 pixel and σ
w
=
1.5 pixels in our experiments).


## 3.3.2. Issues in Multi-Bernoulli Filter

The multi-Bernoulli filter, also called a “track-before-detect” method in [8], by passed the detection module and exploited the spatiotemporal information directly from the image sequence to estimate the multiobject state. The key step is to model the measurement model and to design the likelihood function g
k
(
·
∣
·
) in (4). Let X
=
{
x
1
,
…
,
x
n
} denote the multi-cell state vectors, and let y
k
=
[
y
1
,
k
,
y
2
,
k
,
…
,
y
m
,
k
] denote the image observation comprising an array of m pixels (or bins) value y
i
,
k of the ith pixel, which can be a real number or a vector depending on the application. If each cell or state is represented by a rectangular blob, then the corresponding histogram v
i based on the cell state x
i can be calculated (
i
=
1
,
…
,
n
). Based on the histogram difference comparison between the interested blob region and the cell training database, the similarity between the candidate and real cell can be estimated by (28) using total T
k cell training vectors {
v
j
*
}
j
=
1
T
k. We have
(28)
g
k
(
y
k
∣
x
i
)
=
(
e
-
(
1
-
d
(
v
i
)
)
γ
)
ξ
,
(29)
d
(
v
i
)
=
1
T
k
∑
j
=
1
T
k
∑
l
=
1
C
min
⁡
⁡
(
v
i
(
l
)
,
v
j
*
(
l
)
)
,
where γ and ξ are the adjustment coefficients designed for achieving better likelihood discrimination, v
i
(
l
) denotes the value of the lth element of the color histogram vector v
i, and C is the total number of elements in a histogram vector and takes the values of 256 × 3 for an RGB image.In addition, it is observed that the shape and brightness of cells might change over time. In order to reduce the negative effect of this interference, the cell training template should be updated in time. Specifically, assume that the cell templatet
k has total T
k vectors {
v
j
*
}
j
=
1
T
k at time step k, and if the similarity between the potential and template is greater than a threshold T
0 (i.e., g
k
(
t
k
∣
x
^
)
>
T
0), the potential x
^ is added to the cell template resulting as t
k
+
1
=
{
t
k
,
x
^
}. Obviously, the dimension of the cell template gradually increases with the update process, and we set T
max
⁡
=
10 to limit the dimension of cell template in our experiments for the sake of easy implementation and less computational burden.


## 3.3.3. Label Management

The common approach of label management for the GM-PHD filter and multi-Bernoulli filter is briefly described as follows. Each Gaussian component or Bernoulli particle at timek simply inherits the label of its predecessor (in prediction step) if it belongs to one of the first predicted M
k
-
1 Gaussian components or Bernoulli particles; otherwise, it belongs to a newly born target and is assigned a new label, and the labels remain unchanged during the update step of the filter [10, 22]. However, this method works best when the targets are well separated. In fact, cell movement is random and changeable over time. When two cells occlude each other, they are easily merged in one frame and separated in another frame. As a result, one of them would be detected by the filter as a newly born cell, and would be assigned a new label as a new track. To solve the label management problem, we present a method developed from a graph theoretical framework in which the previously missed cells are taken into account besides the most recent cells, for their association with currently detected cells [10]. This label management scheme has a memory to keep the labels of missed cells for later use in case they reappear. The recorded labels will only be removed from memory if they are missed for several consecutive frames. This method employs a three-stage constrained nearest neighbor graph search to find candidate tracks in a bipartite graph. The search is constrained in the sense that an edge is chosen only if its cost is less than a certain limit [8].Initialization. The initialization takes place in the first frame. The estimates of locations of cells are read and labeled and play the role of previous cells for the next frame. After initialization, an iterative labeling contained three-stage constrained nearest neighbor graph search and memory update is repeated for each of the next frames.Stage 1. Some of the current estimated cells are matched to some of the estimates from the previous iteration.Stage 2. The remainder of the current estimates is matched to some previous cells that have been missed at least once.Stage 3. The remainder of the current estimates is checked if they are close enough to the boundaries of the image, and if validated, they are assigned new labels.Memory Update. The previous cells that have been missed for too long are removed from the memory, and other identity data is updated for the next iteration.


## 4. Experiments and Discussions

In this section we will apply the above presented two RFS-based tracking methods on two real cell image sequences and quantitatively analyze and compare their performance in different scenarios. First, we present several multicell-related measures to evaluate the tracking performance in Section4.1. Then the analysis and comparison of the results of two experiments are presented in Section 4.2.


### 4.1. Performance Measures

(1) Root mean square error (RMSE): due to lack of ground truth for real data, we evaluate the accuracy performance of tracking by comparison with manual tracking results using a traditional quantitative performance measure, that is, the root mean square error (RMSE):(30)
RMSE
k
=
1
IR
∑
i
=
1
IR
RMSE
i
,
k
2
with
(31)
RMSE
i
,
k
2
=
1
M
∑
m
=
1
M
∥
x
m
,
k
-
x
^
m
,
k
∥
2
,
where IR is the total number of successfully tracking results (we used IR
=
10 in our experiments), M is the total number of cells at step k, and x
m
,
k and x
^
m
,
k are defined as the true and estimate positions of the cell m at time step k, respectively.(2) Optimal subpattern assignment (OSPA): besides the RMSE measure, in multiobject tracking system, we usually use an alternative metric OSPA designed under the basis of Wasserstein distance (WD), proposed in [23]. It overcomes some drawbacks by introducing a parameter which controls the relative emphasis of localizations and cardinality errors, and can accommodate cardinality differences in a mathematically consistent and physically meaningful sense. The definition of OSPA metric is briefly described as follows. Let d
(
c
)
(
x
,
y
)
≜
min
⁡
(
c
,
d
(
x
,
y
)
) denote the distance between x and y with a nonnegative cut-off parameter c (we used c
=
100 in our experiments); Ψ
k is the set of permutations on {
1,2
,
…
,
k
} for any k
∈
ℕ
=
{
1,2
,
…
}. For 1
≤
p
<
∞, m
,
n
≥
0 and any two finite subsets X
=
{
x
1
,
x
2
,
…
,
x
m
} and Y
=
{
y
1
,
y
2
,
…
,
y
n
}, if m
≤
n, the OSPA is then defined as
(32)
d
-
p
(
c
)
(
X
,
Y
)
≜
(
1
n
(
min
⁡
ψ
∈
Ψ
n
∑
i
=
1
m
d
(
c
)
(
x
i
,
y
ψ
(
i
)
)
p

+

c
p
(
n
-
m
)
)
)
1
/
p

,
and d
-
p
(
c
)
(
X
,
Y
)
=
d
-
p
(
c
)
(
Y
,
X
) otherwise. In addition, if p
→
∞, the OSPA metric is
(33)
d
-
∞
(
c
)
(
X
,
Y
)
≜
{
min
⁡
ψ
∈
Ψ
n
max
⁡
1
≤
i
≤
n
d
-
(
c
)
(
x
i
,
y
ψ
(
i
)
)
if
m
=
n
,
c
otherwise
.
Note that if m
=
n
=
0 we set the distance to 0 directly for both cases.(3) Detection and association errors: there are some other measures to evaluate the performance of cell detection and association such as false negative rate (FNR), false alarm rate (FAR), label switching rate (LSR), and lost tracks ratio (LTR). FNR and FAR quantify the detection performance, defined as the total number of cells that are missed and the total number of nonexisting cells that are detected, respectively, both normalized over the accumulative total number of true cells over all frames. LSR and LTR quantify the association performance. LSR is the number of label switching events normalized over total number of ground truth tracks crossing events, while LTR is the number of tracks lost for more than 50% of their lifetime normalized over total number of ground truth tracks.


### 4.2. Experiments on Real Cell Image Sequences

We evaluate the performances of the two RFS-based methods on two real low-contrast multi-cell image sequences including different cell dynamics, cell morphology (shape) variation, closely moving cells, and varying number of cells in different frames. In all cases, we assume a constant survival probabilityp
S
=
0.99 and a constant detection probability p
D
=
0.98. For the GM-PHD filter, we consider a predefined model for birth Gaussian components γ
k, that is, we divide the cell image into nine independent parts and Gaussian components are initialized with the location being uniformly distributed within each part and with a random velocity located within [
-
5
pixels
,
5
pixels
], and the variance P
γ
,
k
=
diag
⁡
(
2
pixels
,
0.2
pixel
/
s
,
2
pixels
,
0.2
pixel
/
s
) and the weight of each Gaussian component is a random number. Moreover, the number of Gaussian components of each part is a random number between 100 and 200; the standard deviation of the Gaussian observation model R
k
=
diag
⁡
(
2
pixels
,
2
pixels
) and the intensity of clutter is set to 1
0
-
6. For the multi-Bernoulli filter, we consider a predefined model for birth Bernoulli terms denoted by known parameters {
(
r
Γ
,
k
(
i
)
,
p
Γ
,
k
(
i
)
(
·
)
)
}
i
=
1
M
Γ
,
k where the density p
Γ
,
k
(
i
)
(
·
) is represented by the particles {
(
w
Γ
,
k
(
i
,
j
)
,
x
Γ
,
k
(
i
,
j
)
)
}
j
=
1
L
Γ; in our experiments, we assume that one cell appears in each of nine parts of the image with the location of the particle being uniformly distributed within the part (i.e., M
Γ
,
k
=
9) and with a random velocity located within [
-
5
pixels
,
5
pixels
], and set r
Γ
(
1
)
=
⋯
=
r
Γ
(
9
)
=
0.1 and L
Γ is a random number between 100 and 200; furthermore, the process noise is a Gaussian variable with zero mean and variance P
k
=
diag
⁡
(
2
pixels
,
0.2
pixel
/
s
,
2
pixels
,
0.2
pixel
/
s
); the initial cell template t
1 is determined according to cells in the first frame and keeps updated using the method discussed in Section 3.3.Case
1. The sequences are 128
×
128 pixels RGB image with 35 frames, and the sampling period of image is Δ
=
4 s. Each cell is indicated by the rectangle 10
×
10 pixels, as shown in Figure 2, and it can be observed that the upper cell’s displacement is relatively smaller and moves with nearly round shape, while the lower cell has a larger displacement and moves with irregular shape across the entire image and cell 3 partially enters the image from the lower boundary at frame 35. Figure 3 gives the position estimates of each cell in x and y directions, respectively, which are compared with the results of manual tracking method. Figure 4 plots the velocity and acceleration estimates of cell 2 using different methods. Figures 2 to 4 show that the GM-PHD filtering method and multi-Bernoulli filtering method can both automatically track each cell; moreover, the velocity and acceleration of cell can be roughly estimated. Figure 5 plots the RMSE curves (IR
=
10) and shows that the average errors of the GM-PHD filter are 0.20 pixel and 0.16 pixel in x and y directions, respectively, while 1.15 pixels and 1.31 pixels are computed for the multi-Bernoulli filtering method. Obviously, both methods have high accuracy, but the RMSE of the multi-Bernoulli filtering method is more than six times that of the former.Tracking partial results of original RGB image sequences, respectively, frame 9, 19, 22, 30, and 35.
(a)
Partial original RGB image sequences
(b)
Partial results of GM-PHD filter
(c)
Partial results of multi-Bernoulli filterFigure 3
The corresponding positions of each cell inx and y directions per frame using different methods.Cell instant velocity (a) and acceleration (b) estimates (cell 2 in Figure2).
(a)
(b)Figure 5
RMSE in cell position estimation using the GM-PHD filter and multi-Bernoulli filter.Figure6 shows the estimated average number of cells in 10 independent Monte-Carlo runs against the true cell number along with the corresponding OSPA distance. It can be observed that the OSPA distance becomes larger for penalizing such discrepancy when the estimated number of cell does not coincide with the true one (such as the OSPA distance varying from 0.56 to 70.72 in frame 7 using the GM-PHD filter), and the OSPA distance remains a small value if the estimated number of cell is equal to the true one. In frames 7 and 17, the estimated number by the GM-PHD filtering method is always larger than the true numbers due to clutters. And the errors of estimated number by the multi-Bernoulli filtering method often happen at frame 35 and between frame 19 and frame 27. The average OSPA distance of each filter is 7.50 and 15.03, respectively. Through comparative analysis, it can be observed that the average OSPA distance of the GM-PHD filtering method is only half that of the multi-Bernoulli filtering method; that is, the performance using the GM-PHD filter is increased by 100.4%.Cell number estimate and OSPA distance of two methods. ((a) GM-PHD filter and (b) multi-Bernoulli filter).
(a)
(b)We record all false negative reports, false alarm reports, label switching reports, and tracks lost reports in each frame over 100 Monte-Carlo simulations, and the averaged FNR, FAR, LSR, and LTR are computed as presented in Table1. We find there is only a small amount of false negative reports in frame 1 or frame 2 when using the GM-PHD filter, and the multi-Bernoulli filter is insensitive to cell 3 in Figure 2. By comparing the results, the GM-PHD filtering method performs better than the multi-Bernoulli filtering method. Finally, the total time of treatment is an important performance index for multi-cell tracking. We calculate the running time of each module such as detection module, prediction module, update module, plot module, and others (including pruning, merging, and label management). As shown in Figure 7, the average total time of treatment for the GM-PHD filter and the multi-Bernoulli filter is 38.6118 s and 30.0957 s, respectively (averaged over 10 independent runs), and it is noted that most of computation has been spent in the update module.Table 1
Performance comparison results.
Methods
FAR (%)
FNR (%)
LSR (%)
LTR (%)
GM-PHD
0
2.70
0
0
Multi-Bernoulli
1.13
5.55
7.00
0.50Figure 7
The average total time of treatment to different methods (10 independent runs).Case
2. In this case, the sequences are 201
×
201pixels RGB image with 20 frames. Each cell is indicated by the rectangle 15
×
15pixels, and the sampling period of image is Δ
=
3 s. Figure 8 presents an example of successfully tracking results using the GM-PHD filtering method and the multi-Bernoulli filtering method, respectively. It can be observed that two methods both track at most 7 cells over the scene and one cell disappears in frame 12 and three cells (i.e., cell 5, cell 6, and cell 7) sequentially enter the scene from the boundary of frames 2, 12, and 18, respectively. Figure 9 plots the corresponding position estimate of each cell in x and y directions per frame and indicates that two methods have different sensitivity to some cells which are partially entering or leaving the image such as cell 5 and cell 6. Figure 10 plots one cell’s instant velocity and acceleration estimates by different filters. It can be observed that the GM-PHD filter constantly corrects the velocity to make cell close to the true speed. For the multi-Bernoulli filter, although it can make the speed close to the true roughly, there is a bigger fluctuation than the GM-PHD filter.Tracking partial results of original RGB image sequences, respectively, frame 2, 8, 10, 12, and 18.
(a)
Partial original RGB image sequences
(b)
Partial results of GM-PHD filter
(c)
Partial results of multi-Bernoulli filterFigure 9
The corresponding positions estimate of each cell inx and y directions per frame using different methods.Cell instant velocity (a) and acceleration (b) estimates (cell 4 in Figure7).
(a)
(b)Figure11 plots the RMSE curves in cell position estimations. Although the RMSE using the GM-PHD filter is unexpectedly larger in frames 8 and 11 due to the unsatisfactory measurement of detection, the accuracy of the GM-PHD filter is globally superior to the multi-Bernoulli filter, specifically improved by 343.5%. Figure 12 shows the estimated average number of cells in 10 independent runs.Figure 11
RMSE in cell position estimations using different methods.Cell number estimate and OSPA distance of two methods. ((a) The GM-PHD filter and (b) the multi-Bernoulli filter.)
(a)
(b)Monte-Carlo runs against the true cell number along with the corresponding OSPA distance. It can be observed that in frames 10 and 11 the estimated cell numbers using two methods are not consistent with the true number due to cell 6 looming in the image boundary and result in large OSPA distance in these frames. The average OSPA distance of each filter is 6.27 and 21.89, respectively, and we can deduce that the performance of the GM-PHD filter is more than 3 times higher than the multi-Bernoulli filter. Furthermore, the averaged FNR, FAR, LSR, and LTR once again prove the superior performance of the GM-PHD filter in Table2. Figure 13 records the average total time of treatment, 70.9262 s and 32.6588 s, respectively, and shows that the spending time of the GM-PHD filter is more than twice the spending time of the multi-Bernoulli filter.Table 2
Performance comparison results.
Methods
FAR (%)
FNR (%)
LSR (%)
LTR (%)
GM-PHD
0
2.58
0
0
Multi-Bernoulli
0.10
5.26
2.00
1.14Figure 13
The average total time of treatment to different methods (10 independent runs).


## 4.1. Performance Measures

(1) Root mean square error (RMSE): due to lack of ground truth for real data, we evaluate the accuracy performance of tracking by comparison with manual tracking results using a traditional quantitative performance measure, that is, the root mean square error (RMSE):(30)
RMSE
k
=
1
IR
∑
i
=
1
IR
RMSE
i
,
k
2
with
(31)
RMSE
i
,
k
2
=
1
M
∑
m
=
1
M
∥
x
m
,
k
-
x
^
m
,
k
∥
2
,
where IR is the total number of successfully tracking results (we used IR
=
10 in our experiments), M is the total number of cells at step k, and x
m
,
k and x
^
m
,
k are defined as the true and estimate positions of the cell m at time step k, respectively.(2) Optimal subpattern assignment (OSPA): besides the RMSE measure, in multiobject tracking system, we usually use an alternative metric OSPA designed under the basis of Wasserstein distance (WD), proposed in [23]. It overcomes some drawbacks by introducing a parameter which controls the relative emphasis of localizations and cardinality errors, and can accommodate cardinality differences in a mathematically consistent and physically meaningful sense. The definition of OSPA metric is briefly described as follows. Let d
(
c
)
(
x
,
y
)
≜
min
⁡
(
c
,
d
(
x
,
y
)
) denote the distance between x and y with a nonnegative cut-off parameter c (we used c
=
100 in our experiments); Ψ
k is the set of permutations on {
1,2
,
…
,
k
} for any k
∈
ℕ
=
{
1,2
,
…
}. For 1
≤
p
<
∞, m
,
n
≥
0 and any two finite subsets X
=
{
x
1
,
x
2
,
…
,
x
m
} and Y
=
{
y
1
,
y
2
,
…
,
y
n
}, if m
≤
n, the OSPA is then defined as
(32)
d
-
p
(
c
)
(
X
,
Y
)
≜
(
1
n
(
min
⁡
ψ
∈
Ψ
n
∑
i
=
1
m
d
(
c
)
(
x
i
,
y
ψ
(
i
)
)
p

+

c
p
(
n
-
m
)
)
)
1
/
p

,
and d
-
p
(
c
)
(
X
,
Y
)
=
d
-
p
(
c
)
(
Y
,
X
) otherwise. In addition, if p
→
∞, the OSPA metric is
(33)
d
-
∞
(
c
)
(
X
,
Y
)
≜
{
min
⁡
ψ
∈
Ψ
n
max
⁡
1
≤
i
≤
n
d
-
(
c
)
(
x
i
,
y
ψ
(
i
)
)
if
m
=
n
,
c
otherwise
.
Note that if m
=
n
=
0 we set the distance to 0 directly for both cases.(3) Detection and association errors: there are some other measures to evaluate the performance of cell detection and association such as false negative rate (FNR), false alarm rate (FAR), label switching rate (LSR), and lost tracks ratio (LTR). FNR and FAR quantify the detection performance, defined as the total number of cells that are missed and the total number of nonexisting cells that are detected, respectively, both normalized over the accumulative total number of true cells over all frames. LSR and LTR quantify the association performance. LSR is the number of label switching events normalized over total number of ground truth tracks crossing events, while LTR is the number of tracks lost for more than 50% of their lifetime normalized over total number of ground truth tracks.


## 4.2. Experiments on Real Cell Image Sequences

We evaluate the performances of the two RFS-based methods on two real low-contrast multi-cell image sequences including different cell dynamics, cell morphology (shape) variation, closely moving cells, and varying number of cells in different frames. In all cases, we assume a constant survival probabilityp
S
=
0.99 and a constant detection probability p
D
=
0.98. For the GM-PHD filter, we consider a predefined model for birth Gaussian components γ
k, that is, we divide the cell image into nine independent parts and Gaussian components are initialized with the location being uniformly distributed within each part and with a random velocity located within [
-
5
pixels
,
5
pixels
], and the variance P
γ
,
k
=
diag
⁡
(
2
pixels
,
0.2
pixel
/
s
,
2
pixels
,
0.2
pixel
/
s
) and the weight of each Gaussian component is a random number. Moreover, the number of Gaussian components of each part is a random number between 100 and 200; the standard deviation of the Gaussian observation model R
k
=
diag
⁡
(
2
pixels
,
2
pixels
) and the intensity of clutter is set to 1
0
-
6. For the multi-Bernoulli filter, we consider a predefined model for birth Bernoulli terms denoted by known parameters {
(
r
Γ
,
k
(
i
)
,
p
Γ
,
k
(
i
)
(
·
)
)
}
i
=
1
M
Γ
,
k where the density p
Γ
,
k
(
i
)
(
·
) is represented by the particles {
(
w
Γ
,
k
(
i
,
j
)
,
x
Γ
,
k
(
i
,
j
)
)
}
j
=
1
L
Γ; in our experiments, we assume that one cell appears in each of nine parts of the image with the location of the particle being uniformly distributed within the part (i.e., M
Γ
,
k
=
9) and with a random velocity located within [
-
5
pixels
,
5
pixels
], and set r
Γ
(
1
)
=
⋯
=
r
Γ
(
9
)
=
0.1 and L
Γ is a random number between 100 and 200; furthermore, the process noise is a Gaussian variable with zero mean and variance P
k
=
diag
⁡
(
2
pixels
,
0.2
pixel
/
s
,
2
pixels
,
0.2
pixel
/
s
); the initial cell template t
1 is determined according to cells in the first frame and keeps updated using the method discussed in Section 3.3.Case
1. The sequences are 128
×
128 pixels RGB image with 35 frames, and the sampling period of image is Δ
=
4 s. Each cell is indicated by the rectangle 10
×
10 pixels, as shown in Figure 2, and it can be observed that the upper cell’s displacement is relatively smaller and moves with nearly round shape, while the lower cell has a larger displacement and moves with irregular shape across the entire image and cell 3 partially enters the image from the lower boundary at frame 35. Figure 3 gives the position estimates of each cell in x and y directions, respectively, which are compared with the results of manual tracking method. Figure 4 plots the velocity and acceleration estimates of cell 2 using different methods. Figures 2 to 4 show that the GM-PHD filtering method and multi-Bernoulli filtering method can both automatically track each cell; moreover, the velocity and acceleration of cell can be roughly estimated. Figure 5 plots the RMSE curves (IR
=
10) and shows that the average errors of the GM-PHD filter are 0.20 pixel and 0.16 pixel in x and y directions, respectively, while 1.15 pixels and 1.31 pixels are computed for the multi-Bernoulli filtering method. Obviously, both methods have high accuracy, but the RMSE of the multi-Bernoulli filtering method is more than six times that of the former.Tracking partial results of original RGB image sequences, respectively, frame 9, 19, 22, 30, and 35.
(a)
Partial original RGB image sequences
(b)
Partial results of GM-PHD filter
(c)
Partial results of multi-Bernoulli filterFigure 3
The corresponding positions of each cell inx and y directions per frame using different methods.Cell instant velocity (a) and acceleration (b) estimates (cell 2 in Figure2).
(a)
(b)Figure 5
RMSE in cell position estimation using the GM-PHD filter and multi-Bernoulli filter.Figure6 shows the estimated average number of cells in 10 independent Monte-Carlo runs against the true cell number along with the corresponding OSPA distance. It can be observed that the OSPA distance becomes larger for penalizing such discrepancy when the estimated number of cell does not coincide with the true one (such as the OSPA distance varying from 0.56 to 70.72 in frame 7 using the GM-PHD filter), and the OSPA distance remains a small value if the estimated number of cell is equal to the true one. In frames 7 and 17, the estimated number by the GM-PHD filtering method is always larger than the true numbers due to clutters. And the errors of estimated number by the multi-Bernoulli filtering method often happen at frame 35 and between frame 19 and frame 27. The average OSPA distance of each filter is 7.50 and 15.03, respectively. Through comparative analysis, it can be observed that the average OSPA distance of the GM-PHD filtering method is only half that of the multi-Bernoulli filtering method; that is, the performance using the GM-PHD filter is increased by 100.4%.Cell number estimate and OSPA distance of two methods. ((a) GM-PHD filter and (b) multi-Bernoulli filter).
(a)
(b)We record all false negative reports, false alarm reports, label switching reports, and tracks lost reports in each frame over 100 Monte-Carlo simulations, and the averaged FNR, FAR, LSR, and LTR are computed as presented in Table1. We find there is only a small amount of false negative reports in frame 1 or frame 2 when using the GM-PHD filter, and the multi-Bernoulli filter is insensitive to cell 3 in Figure 2. By comparing the results, the GM-PHD filtering method performs better than the multi-Bernoulli filtering method. Finally, the total time of treatment is an important performance index for multi-cell tracking. We calculate the running time of each module such as detection module, prediction module, update module, plot module, and others (including pruning, merging, and label management). As shown in Figure 7, the average total time of treatment for the GM-PHD filter and the multi-Bernoulli filter is 38.6118 s and 30.0957 s, respectively (averaged over 10 independent runs), and it is noted that most of computation has been spent in the update module.Table 1
Performance comparison results.
Methods
FAR (%)
FNR (%)
LSR (%)
LTR (%)
GM-PHD
0
2.70
0
0
Multi-Bernoulli
1.13
5.55
7.00
0.50Figure 7
The average total time of treatment to different methods (10 independent runs).Case
2. In this case, the sequences are 201
×
201pixels RGB image with 20 frames. Each cell is indicated by the rectangle 15
×
15pixels, and the sampling period of image is Δ
=
3 s. Figure 8 presents an example of successfully tracking results using the GM-PHD filtering method and the multi-Bernoulli filtering method, respectively. It can be observed that two methods both track at most 7 cells over the scene and one cell disappears in frame 12 and three cells (i.e., cell 5, cell 6, and cell 7) sequentially enter the scene from the boundary of frames 2, 12, and 18, respectively. Figure 9 plots the corresponding position estimate of each cell in x and y directions per frame and indicates that two methods have different sensitivity to some cells which are partially entering or leaving the image such as cell 5 and cell 6. Figure 10 plots one cell’s instant velocity and acceleration estimates by different filters. It can be observed that the GM-PHD filter constantly corrects the velocity to make cell close to the true speed. For the multi-Bernoulli filter, although it can make the speed close to the true roughly, there is a bigger fluctuation than the GM-PHD filter.Tracking partial results of original RGB image sequences, respectively, frame 2, 8, 10, 12, and 18.
(a)
Partial original RGB image sequences
(b)
Partial results of GM-PHD filter
(c)
Partial results of multi-Bernoulli filterFigure 9
The corresponding positions estimate of each cell inx and y directions per frame using different methods.Cell instant velocity (a) and acceleration (b) estimates (cell 4 in Figure7).
(a)
(b)Figure11 plots the RMSE curves in cell position estimations. Although the RMSE using the GM-PHD filter is unexpectedly larger in frames 8 and 11 due to the unsatisfactory measurement of detection, the accuracy of the GM-PHD filter is globally superior to the multi-Bernoulli filter, specifically improved by 343.5%. Figure 12 shows the estimated average number of cells in 10 independent runs.Figure 11
RMSE in cell position estimations using different methods.Cell number estimate and OSPA distance of two methods. ((a) The GM-PHD filter and (b) the multi-Bernoulli filter.)
(a)
(b)Monte-Carlo runs against the true cell number along with the corresponding OSPA distance. It can be observed that in frames 10 and 11 the estimated cell numbers using two methods are not consistent with the true number due to cell 6 looming in the image boundary and result in large OSPA distance in these frames. The average OSPA distance of each filter is 6.27 and 21.89, respectively, and we can deduce that the performance of the GM-PHD filter is more than 3 times higher than the multi-Bernoulli filter. Furthermore, the averaged FNR, FAR, LSR, and LTR once again prove the superior performance of the GM-PHD filter in Table2. Figure 13 records the average total time of treatment, 70.9262 s and 32.6588 s, respectively, and shows that the spending time of the GM-PHD filter is more than twice the spending time of the multi-Bernoulli filter.Table 2
Performance comparison results.
Methods
FAR (%)
FNR (%)
LSR (%)
LTR (%)
GM-PHD
0
2.58
0
0
Multi-Bernoulli
0.10
5.26
2.00
1.14Figure 13
The average total time of treatment to different methods (10 independent runs).


## 5. Conclusions

In this paper, we investigate the multi-cell tracking performances using two RFS-based filtering methods, that is, Gaussian mixture probability hypothesis density filter and multi-Bernoulli filter. Numerous performance comparisons between the two filtering methods have been carried on. According to the comparisons of our proposed performance measures, the GM-PHD filtering method has higher accuracy than the multi-Bernoulli filtering method. Specifically, the GM-PHD filter outperforms the multi-Bernoulli filter in terms of RMSE and OSPA distance, improved by 463.4% and 174.8%, respectively. Furthermore, the FAR and FNR (detection error) and LSR and LTR (tracking error) using the GM-PHD filtering method are reduced by 1.33%, 2.77%, 4.50% and 0.82% (averaged over two cases), respectively. It is noteworthy that the multi-Bernoulli filtering method has an advantage in the processing time, especially in Case 2 where the running speed is more than twice as fast as the GM-PHD. Specifically, for the GM-PHD filtering method, the average total time of Cases 1 and 2 is 38.6118 s and 70.9262 s, respectively, while for the multi-Bernoulli filtering method they are 30.0957 s and 32.6588 s, respectively. The main cause is that the number of cells in Case 2 (up to 7 cells) is bigger than Case 1 (up to 3 cells); thus it needs to spend more time in the update module. However, it is noted that the multi-Bernoulli filter has noticeably lower processing time in the case of tracking numerous cells. Finally, factors affecting the performance of the GM-PHD filter include detection module (affected the tracking accuracy) and cell number (affected the processing time), but the main factor affecting the performance of the multi-Bernoulli filter is the design of the likelihood function.



---

*Source: 495765-2014-01-22.xml*