I Introduction
The ability to accurately monitor cellular and subcellular structures in their native biological environment has enormous potential in addressing open questions in cell biology. In various applications, one of the key steps for understanding biological phenomena is to assess the motion of these structures. Recent developments in timelapse cell microscopy imaging systems have had a great impact on the analysis of these dynamics. However, visual inspection of sequences acquired by these imaging techniques requires manual tracking of many tiny structures in numerous noisy images. Thus, automated multitarget tracking methods have been extensively used in different biological applications in the last decade [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25].
Despite significant technical advances made in automatically tracking moving objects, tracking microscopic structures, known as particles [26], remains a challenging task due to the complex nature of biological sequences. The microscopic sequences are usually populated with similar tiny structures having intricate motion patterns and sophisticated interactions with other structures. Moreover, the structures may enter or disappear from the field of view or be occluded by other objects. In addition in some imaging techniques, i.e. fluorescence microscopy imaging, the sequences are contaminated with a high degree of noise. Under these conditions, detection methods usually fail to detect all particles and generate many spurious measurements (clutter) [27, 14]. To be successfully applied in many biological applications, multitarget tracking methods should be able to track an unknown and timevarying number of similar particles in the presence of clutter and detection uncertainty.
To this end, many particle tracking approaches have been proposed in literature [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24]. Some of the most popular particle tracking approaches are based on detection followed by a deterministic linking procedure. Here, each particle is separately detected in each frame. Then, a deterministic solution, e.g. an optimization technique, links the corresponding targets between frames [3, 4, 5, 6, 7]. The performance of these algorithms is often sensitive to the detection algorithm and may degrade in the presence of complex target dynamics and highly cluttered detections resulting from very noisy sequences.
Bayesian filtering approaches are another class of tracking algorithms that have become popular for particle tracking applications in recent years [10, 11, 12, 13, 14, 15, 24, 16, 17, 18, 19, 20, 21, 22, 23]. These approaches better deal with tracking multiple particles in cluttered measurements by incorporating prior knowledge of target dynamics and measurement models [14, 15, 19, 17, 16, 18]. However, knowledge of the clutter rate and detection profile of the chosen detection method, are of critical importance in Bayesian multitarget tracking.
Most existing multitarget tracking solutions assume known and fixed detection and false alarm parameters. However, these parameters cannot be computed in many practical applications and worse, it is not even known whether they are timeinvariant. In biological imaging techniques, the noise characteristic and the background intensity of sequences may change during the acquisition process which make detection profile and clutter rate timevariant. For instance, injection of a stimulus such as insulin into pancreatic beta cells increases the noise level and overall intensity of the sequences acquired by total internal reflection fluorescence microscopy (TIRFM) [28] (Fig. 1). Thus, the ability of multitarget trackers to accommodate unknown clutter and detection parameters is crucial in timelapse cell microscopy since mismatches in these model parameters inevitably result in erroneous tracking outputs.
The key contribution of this paper is to propose an effective solution to the problem of tracking multiple maneuvering particles in unknown and timevarying false alarm and detection rates. To the best of our knowledge, this practical problem has not been discussed in the biological signal processing literature so far. To address this challenging problem, we use the recent generation of Bayesian filters based on random finite set theory. This framework provides an elegant mathematical formulation for multitarget systems and allows us to deal with the aforementioned complexities. Our proposed approach is based on the Cardinalized Probability Hypothesis Density (CPHD) filter for unknown clutter rate and detection profile known as the CPHD filter [29]. Clutter rate and detection probabilities are estimated by the CPHD filter bootstrapped onto a CPHD filter that outputs target estimates. This bootstrap idea was inspired by [30] which requires known and uniform probability of detection. Our proposed solution can accommodate unknown and nonuniform probability of detection. To cope with maneuvering motion of subcellular structures [19, 18, 21, 17, 16, 20], we also propose the multiple model implementation of the aforementioned filters in this paper. We will show that the proposed method is able to deal with tracking maneuvering structures in the presence of unknown and timevarying clutter rate and detection profile.
Ii Background
The Bayesian estimation paradigm deals with the problem of inferring the state of a target from a sequence of noisy measurements. The state vector contains all relevant dynamic information about the target while the measurements contain what can be directly observed from the sequences. The notion of Bayes optimality is fundamental to the Bayesian estimation paradigm and is welldefined for singletarget tracking
[31]. While Bayes optimal tracking techniques such as Kalman and particle filters are formulated for a singletarget, they can be algorithmically extended to track a timevarying number of targets by combining them with data association [32, 33, 34]. However, the notion of Bayes optimality does not carry over. Nonetheless this approach to multitarget tracking has been used in a wide range of applications, including particle tracking [13, 21, 22, 14, 19, 18, 11, 12, 23, 20].The random finite set (RFS) approach, introduced by Mahler generalizes the notion of Bayes optimality to multitarget system using RFS theory [31]. The key difference with other approaches is that the collection of target states at any given time is treated as a setvalued multitarget state. This representation provides an elegant formulation for multitarget system which avoids the need for data association. The RFS approach has generated substantial interest in recent years with the development of the Probability Hypothesis Density (PHD) filter [35], Cardinalized PHD (CPHD) filter [36], multiBernoulli filters [37, 38], labeled RFS filter [39] and a host of many applications [40] including cell and particle tracking [16, 17, 15, 24].
Iia RFS for multitarget tracking
In the Bayesian estimation paradigm, the state and measurement are treated as realizations of random variables which are suitable for single target and single measurement system. However, in a multitarget system, the collection of states and measurements can be naturally represented as finite sets.
Let and be the states of all targets and all measurements at time , respectively. Over time, some of these targets may disappear, new targets may appear, and the surviving targets evolve to new states. Moreover, due to poor detector performance, only some targets are detected at each time step and many measurements are spurious detections (clutter). Thus, we can conveniently represent the targets and measurements at each time slice with two finite sets as
(1)  
Consequently, the concept of a random finite set (RFS) is required to cast the multitarget estimation problem in the Bayesian framework. Intuitively, an RFS is a finitesetvalued random variable. What distinguishes an RFS from a random vector is that the number of constituent variables is random and the variables themselves are random, distinct and unordered. Mahler’s Finite Set Statistics (FISST) provides powerful yet practical mathematical tools for dealing with RFSs [35, 31], based on the notion of integration and density that is consistent with point process theory [41]. The centerpiece of the RFS approach is the so called Bayes multitarget filter, a generalization of the singletarget (optimal) Bayes filter to accommodate multiple targets based on multitarget dynamic and measurement models.
IiA1 MultiTarget Dynamical Model
Given a multitarget state at time , each target either continues to exist at time with probability and moves to a new state with probability density , or dies with probability and takes on the value . Thus, given a state at time , its behavior at time is modeled by a Bernoulli RFS, , which has FISST density^{1}^{1}1For more detail on RFS density and probability density please see [31] given by
(2)  
Assuming that individual targets move independently, each of the above Bernoulli RFSs are mutually independent conditioned on . Thus the survival or death of all existing targets from time to time is modeled by the multiBernoulli RFS
(3) 
The appearance of new targets at time is modeled by an RFS of spontaneous births which is usually specified as an i.i.d. cluster RFS with intensity function and cardinality distribution ^{2}^{2}2An i.i.d. cluster RFS can be statistically represented by its intensity and cardinality distribution. The cardinality distribution,
, provides a discrete probability distribution over the number of the constituent random variables,
, while the intensity distribution at the point in the state space is interpreted as the instantaneous expected number of the variables that exist at that point [31].. Consequently, the RFS multitarget state at time is given by the union(4) 
The likelihood of the multitarget state , at time , is described by the multitarget transition density which incorporates target dynamics, births and deaths [31].
IiA2 MultiTarget Measurement Model
Given a multitarget state at time , each target , at time , is either detected with probability and generates an observation with likelihood , or missed with probability and generates the value , i.e. each target generates a Bernoulli RFS, , which has FISST density given by
(5) 
The set of measurements generated by , at time , is modeled by the multiBernoulli RFS
(6) 
In addition, the sensor receives a set of false/spurious measurements or clutter, modeled by an RFS , which is usually specified as an i.i.d. cluster RFS with intensity function and cardinality distribution [31]. Consequently, at time , the multitarget measurement generated by the multitarget state is formed by the union
(7) 
The likelihood of the multitarget measurement , given the multitarget state, at time , is described by the multitarget measurement likelihood which incorporates detection uncertainty and clutter.
IiA3 Bayes MultiTarget Filter
The multitarget posterior density contains all statistical information about the multitarget state at time and can be recursively computed by the aforementioned multitarget transition density and measurement likelihood using the following prediction and update equations.
(8)  
(9) 
where is the set integral [31].
This filter generalizes the notion of Bayes optimality to multitarget system using RFS theory. Further details can be found in [31].
IiB The PHD filters
Although the above framework provides an elegant Bayesian formulation of the multitarget filtering problem, it can be computationally expensive for applications with large number of targets [35]. In [35], Mahler proposed to propagate the probability hypothesis density (PHD), or posterior intensity distribution, of the targets
which is the first statistical moment of the probability density function
. The primary weakness of the PHD filter is a loss of higher order cardinality information which causes erratic estimates in the number of targets especially when the density of targets is high [36, 42, 16].Mahler [36] subsequently proposed the Cardinalized PHD (CPHD) filter which jointly propagates the posterior intensity function, , and cardinality distribution, . The propagation of the cardinality distribution makes the filter more robust in the estimation of the number of targets [36].
As with all Bayesian filtering frameworks, this filter also requires some predefined models such as the single target Markov transition density for the dynamic model, measurement likelihood and probabilities of survival and detection . Moreover, the models for new born targets and false alarms (clutter) are described by the birth cardinality distribution and the birth intensity function , and the clutter cardinality distribution and the clutter intensity function , respectively [36]. Obviously, this filter requires knowledge of the detection probability and clutter distributions similar to other Bayesian filters.
In [29], it was shown that the CPHD filter can also be reformulated such that it jointly estimates clutter distributions and detection probability while tracking. We refer to this version of the filter as the CPHD filter. However, the CPHD cannot naturally perform as well as the CPHD filter when exact knowledge of detection probability and clutter density is available.
Iii The Bootstrap CPHD Filter
To propose an effective multitarget particle tracker for practical applications with unknown and timevarying clutter rate and detection profile, we borrow the bootstrap idea proposed in [30]. However, in addition to the clutter rate, we calculate the detection probability, using the CPHD filter, and feed these to a robust multitarget tracking filter such as the CPHD filter (Fig. 2). We use the CPHD filter for target estimation because it is a good tradeoff between accuracy and speed. Nonetheless, it can be replaced by any multitarget filter that requires knowledge of false alarm and detection rate. In this paper, we assume that and the clutter distributions, and , are the meta parameters to be estimated for the CPHD filter in each time frame. However, we will show that the estimation of the mean of the clutter rate, , is sufficient to uniquely define the and .
To deal with maneuvering dynamics, we assume that target motion can be characterized by multiple dynamic models. Thus, here, we propose the multiple model (or jumpMarkov) implementation of both the CPHD and CPHD filters and use these filters as the tracker and the parameter estimator, respectively.
Iiia Multiple Model CPHD filter
Complex maneuvering motions can be often characterized by multiple simpler dynamic models. The idea of tracking a target with multiple motion models is to augment the state of the targets by the index of the model [43].
For notational convenience, we remove the time index in our formulation throughout the paper. However, the random variables and distributions are generally timeindexed. All random variables at time and at time are simply denoted by and , respectively.
Let and denote the kinematic state space and the multiple models state space, respectively. Define the augmented state space as , where denotes the Cartesian product. The underscore notation is used for the augmented state space so that denotes for the kinematic state and for the model state. The integral of a function is given by
(10) 
The measurement likelihood and probabilities of survival and detection for the augmented state vector can be respectively written as
(11) 
(12) 
(13) 
Also, for the augmented transition density and birth intensity, we have the following factored forms,
(14) 
(15) 
where is the probability transition from model to and is the probability of birth for model [43]. As described in Section II, the filter also requires the knowledge of the birth cardinality distribution .
Clutter is usually assumed to be Poisson RFS and independent from the target state [31, 33, 32, 34] with
(16) 
(17) 
where
is a known probability distribution, e.g. uniform distribution, representing how the false alarms are distributed over the measurement space and
. Therefore, given , knowledge of the parameter is sufficient to define the clutter distribution.Substituting the aforementioned augmented model terms (Eqs. 1117) into the conventional CPHD equations [42] and using Eq. 10 yields the recursive equations for the multiple model or jump Markov CPHD filter (see Appendix A).
Fig. 3 depicts a simple graphical representation of the proposed multiple model CPHD filter. At each time step , the intensity and cardinality distributions are predicted and updated from the intensity and cardinality distributions in the previous time step and using the augmented model terms (Eqs. 1117). The intensity distribution provides information on the state of the targets while the number of targets can be estimated using the cardinality distribution. Specifically, the number of targets is estimated using the posterior mode [42]. Note that this filter requires the clutter rate and the detection probability as prior knowledge.
IiiB Multiple Model CPHD filter
We borrow a key idea from Mahler et al. [29] to develop a multiple model CPHD filter with unknown clutter rate and detection profile. To estimate clutter, the idea is to model it as a random finite set of false targets which is statistically independent from the set of actual targets. These false targets are defined by their own models such as birth, death, survival and detection probabilities and transition model. Therefore, the multitarget state is composed of a disjoint combination of these two finite sets including actual targets and clutter. This hybrid multitarget state allows us to track both actual and false targets simultaneously. To deal with multiple model dynamics and unknown detection probability, the state of targets are augmented by the index of the models and the unknown detection probability. Finally this hybrid and augmented state is estimated from the sequence of finite sets of measurements generated by both real targets and clutter while accommodating multiple model dynamics and estimating the detection probability and clutter rate.
Let denote the state space for actual targets, denote the state space for clutter generators and , and denote the state space for the unknown detection probability, and the multiple models for actual targets and clutter, respectively. Define the hybrid and augmented state space as
(18) 
where denotes a disjoint union. The double dot notation is used throughout to denote a function or variable defined on the hybrid state space and the underscore notation is used for the augmented state space. Therefore, represents a hybrid and augmented state such that and denotes the actual and clutter states comprising the kinematic state and augmented components, respectively. The integral of a function is given by
(19) 
In theory, there may exist multiple model clutter generators. However, in most tracking applications, it is assumed that the clutter is uniformly distributed with Poisson cardinality [31, 33, 32, 34]
. As a result, multiple model clutter generators with uniformPoisson distribution can be substituted by a single model clutter generator with uniform distribution and higher Poisson mean. Therefore, we continue our derivations with the single model clutter generators by eliminating the variable
from the equations. Deriving the equations for applications with multiple clutter generators and nonuniform distribution is straightforward.In addition, since clutter generators are identical and the false targets do not follow any motion pattern, it is reasonable to ignore any functional dependence on the state of a clutter generator , [29]. As a result, the probabilities of survival and detection and the measurement likelihood for the hybrid and augmented state vector can be respectively written as
(20) 
(21) 
(22) 
where is clutter density which is often assumed to be uniform distribution in the measurement space. As the survival probabilities and the measurement likelihoods are independent from the detection probabilities, the terms and are removed from the equations. Obviously, the detection probabilities are only dependent on and .
For the hybrid and augmented birth intensity and transition density, we have the following factored forms,
(23) 
(24) 
where is the transition density for the detection probability and due to the independence of false alarms from the clutter state . Furthermore, the cardinality distribution of birth for the hybrid space is [29].
By substituting the hybrid and augmented state space model terms (Eqs. 20–24) into the conventional CPHD equations [42] and using Eq. 19, the recursive equations for this filter can be calculated (see Appendix B).
Fig. 4 shows a simple graphical representation of the proposed multiple model CPHD filter. At each time step , the target’s intensity , the clutter intensity and the hybrid cardinality distribution are predicted and updated from the target’s intensity , the clutter intensity and the hybrid cardinality distribution in the previous time step and using the hybrid model terms (Eqs. 20–24).
The posterior cardinality provides information on the total number of real and clutter targets. Therefore, the number of actual targets cannot be estimated using the posterior mode , since includes both real targets and clutter. In this filter, the posterior mean is used for estimation of the number of real targets [29]. Therefore, in contrast with the CPHD filter, this filter cannot benefit from propagation of the cardinality distribution for the estimated number of the targets. Similar to the PHD filter, this leads the erratic estimates of the number of targets which worsens the tracking result. To this end, we use this filter only as the estimator in our bootstrap filter. The estimated detection and mean clutter are respectively calculated as follows [29].
(25)  
where is the inner product operator, and .
Iv Experimental results
To evaluate the performance of our proposed filters, we apply them to multitarget tracking in 2D Total Internal Reflection Fluorescence Microscopy (TIRFM) sequences. TIRFM is an imaging technique that enables visualization of subcellular structures such as vesicles that are on or close to the plasma membrane of cells [28]. The vesicles are very tiny subcellular structures and are seen in TIRFM sequences as small particles (bright spots) moving with varying dynamics. These small particles appear and disappear from the field of view and can be occluded by other structures. Due to limitations in the TIRFM acquisition process, the sequences are contaminated with a high level of noise. We consider the case where the main characteristics of the sequences such as noise level and the background intensity gradually increase over time (Fig. 1).
In this specific application, we are interested in the overall motion of the vesicles before and after injection of insulin, not motion of a single object. In this case, the primary concern is not how well a tracking approach maintains the identity of a tracked object over time. Instead, we are mainly concerned with how well a tracker avoids false and missed tracks.
We compare the results of our bootstrap filter against those of several popular stateoftheart particle tracking methods including two reliable deterministic linking techniques, PTracker [6] and UTracker [4], and two robust detection based traditional Bayesian filtering approaches, MHT [14] and IMMJPDA [19].^{3}^{3}3The performance of the first three trackers, PTracker, UTracker and MHT, has been evaluated in a recent particle tracking challenge and their results have been reported in [26]. According to this article, they can be assumed as three top performing methods for similar applications. Our results are also compared against the result of a multiple model PHD filter [16]. Finally, the results of our other derived filters such as the multiplemodel CPHD (MMCPHD) and CPHD (MMCPHD) filters are also reported here in order to show the efficiency of our bootstrap idea.
Iva Setup and Implementation Details
To fairly evaluate the performance of all tracking algorithms, the same detection lists were provided for all competing tracking methods. We chose the detector proposed in [44] as it performs reliably in our synthetic and real sequences.
Having ground truth, we can accurately calculate the clutter (false positive) and the detection probability (true positive rate) of the chosen detector for each time frame. True positive and false positive were calculated based on the optimal subpattern assignment between the ground truth and detected point sets. In this case, if the distance between two optimally assigned points is less than a predefined distance ( pixel), the detected point is counted as true positive, otherwise it is a false positive.
For the MHT, IMMJPDA, MMPHD and MMCPHD filters, the average number of false detections per frame and the mean value of the detection probability (the optimal value for these parameters) were used as the predefined clutter rate and detection probability. However, accurate knowledge of these values is not possible in many practical applications and thus, the reported results for these filters are optimistic. For the MMCPHD filter and similarly for our bootstrap (BMMCPHD) filter, the clutter rate and detection probability are adaptively estimated using our proposed framework.
To ensure the validity of our experiments, for the independently implemented tracking methods such as MHT, PTracker and UTracker, we attempted to either estimate their parameters from the ground truth or to find the values that resulted in their best performance. Moreover, the multiple motion model implementation of all competing tracking methods was used. For the IMMJPDA and different PHD and CPHD filters, we chose identical state and measurement vectors and dynamic models. We modeled the state of each particle by its position, , and velocity, . The measurements vector also contained the estimated position of the particles as . To model maneuvering motion of particles, two linear dynamics including random walk and small acceleration motion models were used [19, 18].
Since the target dynamics and measurement models in this application can be properly expressed by multiple linear and Gaussian terms, we used the Gaussian mixture implementation [42] and the BetaGaussian mixture approach [29] for analytical implementations of the MMCPHD and MMCPHD filters, respectively. The birth intensity distribution
for all the PHD and CPHD filters was set as a Gaussian distribution centered at the image with a very high standard deviation
[16].As previously discussed, the PHD and CPHD filters propagate the intensity distribution of all targets in each time frame . By calculating the estimated number of targets , the state of all targets in each time frame can be easily extracted using [45, 41, 42]^{4}^{4}4For example, in the Gaussian mixture implementations, the mean of Gaussian terms with the highest weights in can be used as the estimation of the state vector at each time frame [45].. However, in this filtering framework, the temporal correspondences between the estimated states is not considered. In other words, the output of the PHD and CPHD filters is a set of the estimated states for each time frame without considering the identity of trajectories. Thus, the dynamics of an individual target cannot be evaluated. To deal with this, some authors combine the PHD filter with a track management technique to maintain the identity of tracks [46, 47]. To avoid any computational burden due to a track management step, we simply use our tag propagation scheme [16], which only propagates the identity of the intensity distributions for all the PHD and CPHD filters. As it only considers the previous time step to propagate the identities [16], this is not the most reliable approach for identitytotrack assignment. However, since we are interested in the overall, not individual, motion of the vesicles, we used this approach in this application.
IvB Evaluation Metric
To qualitatively assess the performance of these tracking methods, we used a recent and popular metric based on optimal subpattern assignment (OSPA) [48] followed by an extension of this metric for multitarget tracking applications, known as OSPAT [49].
IvB1 OSPA Metric
This metric measures the distance between two sets of points, which can be the set of estimated tracks and the ground truth tracks at each time frame, and represents different aspects of multitarget tracking performance such as track accuracy, track truncation and missed or false tracks by a single value.
For two arbitrary sets and , where , the metric of order is defined as
(26) 
where is all feasible sets of permutations on and
(27) 
where is the cutoff parameter. For , is calculated, instead.
In Eq. 26, the first term, known as the location error, measures track accuracy while the second term, the so called cardinality error, represents the error for missed or false tracks. In this metric, the parameter
controls the sensitivity to outlier estimates that are distant from the true targets. Moreover, the metric penalizes the cardinality error by the cutoff parameter
as a higher value for this parameter penalizes the false and missing targets more. At a first glance, the selection of the parameters and seems to be critical for the performance of different methods. However, it was proved that the relation between the performance of different methods is independent of the parameters and and different values for these parameters only change the scale of the error, and do not affect the method ranking [48].IvB2 OspaT
Although the OSPA metric measures the multitarget tracking errors such as track accuracy and missed or false tracks, it does not evaluate some other errors such as inconsistent labeling (label switching between the targets in crossing cases) and incorrect label initiation in the case of track truncation. An extension of the OSPA metric, known as OSPAT^{5}^{5}5Contrary to the OSPA, the OSPAT is not a metric in the formal mathematical sense. However, it can be used as a performance measure., considers the aforementioned errors by measuring the distance between two sets of labeled points.
Let us assume that we have two labeled sets such that and . In a multitarget tracking problem, these would be the ground truth and estimated state of the targets with their labels in each time frame. In this performance measure, the label correspondence between the estimated and ground truth sets is first performed based on an optimal global assignment using the trajectories’ temporal information [49]. Then, in Eq. 26 is substituted by the following distance.
(28)  
where if and are the corresponding labels, and otherwise. The parameter should be between and and controls the penalty assigned to the labeling error. The case assigns no penalty, and assigns the maximum penalty. In this paper, we reported the results with .
IvC Evaluation on synthetic data
To quantitatively evaluate the tracking algorithms, they were first evaluated using realistic synthetic movies generated by the framework proposed in [50]. The sequences simulated using this framework reflect the difficulties existing in real TIRFM sequences while providing accurate ground truth. In these sequences, the aim is to mimic the effect of the injection of a stimulus such as insulin into a pancreatic cell. Fig. 5 shows two frames of the synthetic sequences using this framework which are comparable with the real TIRFM sequences shown in Fig. 1.
Each synthetic sequence was simulated with spatial resolution of nm/pixel and temporal resolution fps and consists of timevarying number of targets (on average particles per frame) moving through time frames inside a cell membrane (an estimated background) that is extracted from real TIRFM sequences with effective region pixels. The spots were generated in different sizes, pixels ( nm). The dynamics of the targets were modeled using random walk and linear movement. Furthermore, targets can switch between these two dynamics. The number of intersecting and touching spots in each frame was counted according to the Rayleigh resolution [51]. The averaged percentage of intersecting and touching spots per frame in these sequences is equal to .
Due to the 3D motion of the structures, their intensity changes according to the TIRFM exponential equation [50]. Therefore, they may either temporarily or permanently disappear from the frames. New born targets may also gradually appear from the background. The sequences are contaminated with Poisson noise and the main characteristics of the sequences such as background intensity and noise level gradually change based on a model extracted from real sequences. Due to spatiotemporally varying noise level and backgrounds and dynamic intensity of spots, the signaltonoise ratio (SNR)^{6}^{6}6SNR value is used for representing the level of noise and is calculated using definition of Smal et al. [22] which is the difference in intensity between the object and the background, divided by the standard deviation of the object noise. of an object cannot be constant. Instead, it varies between and .
Similar to the real TIRFM data, the spots in the synthetic sequences fade over time as noise level and background intensity gradually increases. From a biological perspective, it is important to assess the dynamics of vesicles after injection of the stimulus which leads the escalation of noise level and background intensity. Therefore, the threshold in the detection method needs to be set low to ensure consistent detections of the objects in order to avoid early track termination. To this end, we chose a value for the threshold such that the averaged detection probability for all sequences is about . However, this scenario noticeably increases the clutter rate and its variation over time. Fig. 6 shows that the mean clutter rate estimated using the proposed MMCPHD filter can appropriately follow the quick changes in the ground truth clutter rate.
In this experiment, the performance of the trackers was evaluated using these timevarying and highly cluttered detections and their results are reported in Table I. In this Table, the errors are averaged over the number of frames in all synthetic sequences. The results show that our bootstrap filter benefits from a reliable estimator, the MMCPHD filter, which accurately estimates the detection probability and clutter rate. Consequently compared to the other tracker, its tracking results contain fewer false and missing tracks, which decreases the cardinality error significantly. Similarly, this error is also low in the MMCPHD filter. The other Bayesian trackers cannot benefit from these estimations as their parameters are fixed. Therefore, they have higher cardinality error. In comparison with traditional Bayesian trackers such as the IMMJPDA and MHT trackers, the RFS filters have relatively better cardinality error as they properly incorporate birth, death and clutter models in their formulations. However, the higher difference between the OSPA and OSPAT errors of the all RFS filters represent that their results include more identity switch errors compared to the other trackers due to the simple tag propagation scheme used.
The performance of the Ptracker as one of the deterministic trackers seems to be sensitive to cluttered measurements. Therefore, it has the highest cardinality error of the trackers. In contrast, the results show that another deterministic tracking scheme, Utracker, can robustly track the particles while dealing with highly cluttered detections.
The PTracker has the lowest location error of all trackers. Although this error reflects the accuracy of the trackers in tracking particles, its lower value can be also an artifact of the tracker’s poor performance. The trackers with higher cardinality error have relatively lower location error and vice versa.
Method  Location  Cardinality  OSPA  OSPAT 

BMMCPHD  
MMCPHD  
MMCPHD  
MMPHD [16]  
IMMJPDA [19]  
MHT [14]  
PTracker [6]  
UTracker [4] 
The method rankings may change in other scenarios, e.g. where the clutter rate is low and its variation over time can be ignored. In our application, the clutter rate and its variation over time can be decreased by increasing the detection threshold value. For example, we chose a value such that the averaged clutter rate and detection probability are respectively about and . In this case, the detected points from a target are inconsistent and many faint particles are not detected after initial time frames. This case is undesirable for our application due to many truncated tracks with short life time. Nevertheless, we report the tracker’s OSPA and OSPAT errors in this case in Table II in order to compare their different performance in the low and high clutter rate scenarios. Note that both false tracks in high clutter scenario and missed or truncated tracks in low clutter scenario increase the OSPA and OSPAT errors.
Method  Location  Cardinality  OSPA  OSPAT 

BMMCPHD  
MMCPHD  
MMCPHD  
MMPHD [16]  
IMMJPDA [19]  
MHT [14]  
PTracker [6]  
UTracker [4] 
The results show that the trackers such as the IMMJPDA filter and PTracker which have the worst performance in highly cluttered measurements perform better for a low clutter rate with long missed detections compared to other trackers. In contrast, the UTracker cannot perform reliably in this case as it has the highest cardinality error. Estimating declining trend in the detection probability (Fig. 7) helps our bootstrap tracker to still have lower OSPA error compared to the other trackers, but relatively high OSPAT error due to the simple tag propagation scheme.
In a nutshell, the performance of the trackers varies according to the detector reliability. Our bootstrap approach has superior performance for highly cluttered detections with the timevarying rates.
IvD Evaluation on real data
The tracking methods were also tested on a real TIRFM sequence with spatial and temporal resolution of nm/pixel and fps respectively. The image sequences were acquired from a pancreatic beta cell injected by insulin during the acquisition.
In order to prepare a ground truth, an independent expert manually tracked all visible structures ( trajectories from spots) in a part of the cell within time frames using the freely available software MTrackJ [52]. The annotated trajectories were doublechecked by another biologist to ensure reliability of the ground truth.
In order to detect and track all structures, especially after injection of insulin, the threshold in the detection method was set low, which noticeably increases clutter rate and its variation over time (Fig. 8). The results of the trackers for the real data are reported in Table III. In this Table, the errors are averaged over the number of time steps in the real image sequence.
The higher clutter rate along with significantly higher number of time frames and numerous faint structures existing in real sequences cause higher values of the OSPA and OSPAT errors of the all trackers for real sequences compared to the synthetic data. In addition, the results confirm our arguments about the performance of the trackers in the synthetic sequences for the same high clutter scenario. The ability of our BMMCPHD to properly track true targets while accurately estimating the clutter rate and detection probability (Fig. 8) leads to the lowest cardinality and also OSPA errors.
Since there were always faint structures in real sequences that even experienced annotators were unable to detect or determine whether they are real or false targets, the reliability of the manual ground truth cannot be completely guaranteed. In this case, the quantitative results may be biased to a specific method. To maximize the validity of our comparison, the results of the tracking were also visually assessed by the experts. This assessment also showed that our bootstrap filter can better detect and track the real vesicles, especially faint ones, while avoiding tracking false targets.
In Fig. 9 (a), the tracking results of several crossing particles with maneuvering motions using the proposed bootstrap filter are shown. The results show that it can properly deal with the maneuvering motion of the targets. However, the results are not error free and include some missing and false tracks and switching labeling errors. Fig. 9 (b) also demonstrates the performance of our proposed bootstrap filter in tracking two faint particles moving through the synthetic and real sequences with timevarying background and noise level.
V Conclusion
The estimation of clutter rate and detection probability helps improve tracking of the targets in particle tracking applications where measurements include many spurious and missed detections with unknown and timevarying parameters. The RFS framework provides a principled solution to deal with the estimation of these parameters which was not possible in previous approaches. The CPHD filter is one of the filters derived based on RFS theory which is able to estimate these parameters. However, the filter cannot naturally perform as well as the CPHD filter with known parameters. To this end, in this paper we proposed a bootstrap filter by combination of the CPHD and CPHD filters. To accommodate the maneuvering motion, we also proposed a multiple model implementation of the filters. Therefore, the clutter rate and detection probability are estimated by the multiple model CPHD filter bootstrapped onto a multiple model CPHD filter that outputs target estimates.
The proposed approach was evaluated on a challenging particle tracking application where vesicles move in noisy sequences of total internal reflection fluorescence microscopy while the noise characteristic and background intensity of the sequences change during the acquisition process. In this application which the clutter rate and detection probability are timevarying, we demonstrated that our bootstrap filter is able to better track the real targets and more reliably avoid false tracks compared to the other RFS filters such as the MMPHD, MMCPHD and MMCPHD as well as other stateoftheart particle tracking methods such as the IMMJPDA, MHT, PTracker and UTracker in both synthetic and real sequences. However, in the applications where the motion of each individual targets is required, the tag propagation scheme may need to be changed to a track management algorithm, e.g., that proposed in [47, 46] or more generally, the tracker in the bootstrap filter can be replaced by any multitarget tracker that requires knowledge of false alarm and detection rate and performs reliably in the label assignments. In the RFS concept, a principled solution to the track labeling problem using labeled RFS [39] has been recently proposed which can be also used for this purpose.
Obviously, our approach may not be a good choice for some applications. For example, in the cases where the particles can be easily detected or the rates for false alarms and missed detections are negligible, known or timeinvariant, there is no point to adaptively estimate the clutter rate and the detection profile. In addition, to estimate these timevarying rates, the MMCPHD requires some probabilistic priors describing how the clutter rate and the detection probability change over time [29]. Therefore, this approach requires more parameters compared to other trackers. Moreover, the bootstrap method filters twice in each frame in order to estimate the state of the targets. Consequently, its processing time is more than each of its component filters, MMCPHD and MMCPHD. Our nonoptimized MATLAB code was run on an ordinary PC (Intel Core Quad, GHz CPU, GB RAM). The average CPU processing time per frame per target in the highest considered clutter rate (around detected positions per frame) for the multiple model CPHD, CPHD and our bootstrap filter are about ms, ms and ms respectively.
Due to linearity of the dynamic and measurement models in our application, we only used the multiple model Gaussian linear implementation of all filters. The recursive equations derived in the appendices can be used for both Gaussian linear and nonGaussian nonlinear system models. In order to apply the proposed framework to applications with nonlinear nonGaussian system models, the sequential Monte Carlo (SMC) implementation of the filters is required. This implementation will also allows us to further evaluate the performance of our bootstrap filter and compare it against the SMCPHD [41], SMCCPHD [31] and the traditional particle [34] filters.
Appendix
For completeness, we now derive the recursive equations for both the multiple model CPHD and CPHD filters in this section. The following notations are used in throughout this Appendix. The binomial and permutation coefficients are denoted by and , respectively. is the inner product operation between two continuous or two discrete functions. The elementary symmetric function of order defined for a finite set of real numbers is denoted by
(29) 
where is the cardinality of a set and [42].
A Multiple Model CPHD Recursions
Prediction step: Suppose at time , the posterior cardinality distribution and posterior intensity are known. The predicted cardinality distribution and predicted intensity are calculated by
(30) 
(31)  
where
(32) 
Note that in the multiple model approach, the inner product function operates on both the kinematic state and the model, i.e. .
Update step: If at time , the predicted cardinality distribution , predicted intensity and set of measurement are given, the updated cardinality distribution and updated intensity are calculated by
(33) 
(34) 
where
(35)  
where
(36) 
and
(37) 
B Multiple Model CPHD Recursions
Prediction step: Suppose at time , the hybrid cardinality distributions and the posterior intensity distribution for actual targets and clutter generators are given. The hybrid predicted cardinality distribution and predicted intensity for actual targets and clutter generators are calculated by
(38) 
where
(39) 
(40) 
(41) 
Update step: If at time , the predicted intensity for actual targets , the predicted intensity for clutter generators , the predicted hybrid cardinality distribution and set of measurement are all given and the function defined as follows:
(42) 
where
(43) 
where and .
The updated cardinality distribution and the updated intensity distribution for actual targets and clutter generators are given as follows
(44) 
(45) 
(46) 
Acknowledgment
The authors would like to thank Dr. William E. Hughes and the research team from the Garvan Institute of Medical Research, NSW, Australia, for providing the real TIRFM sequences and their manual ground truth. The authors also thank Dr. Anthony Dick from the University of Adelaide whose comments and suggestions improved presentation of the paper.
References
 [1] A. Matov, M. M. Edvall, G. Yang, and G. Danuser, “Optimalflow minimumcost correspondence assignment in particle flow tracking,” Comput. Vis. Image Und., vol. 115, no. 4, pp. 531–540, 2011.
 [2] S. Bonneau, M. Dahan, and L. D. Cohen, “Single quantum dot tracking based on perceptual grouping using minimal paths in a spatiotemporal volume,” IEEE Trans. Image Process., vol. 14, no. 9, pp. 1384–1395, 2005.
 [3] G. Mashanov and J. Molloy, “Automatic detection of single fluorophores in live cells,” Biophys. J., vol. 92, no. 6, pp. 2199–2211, 2007.
 [4] K. Jaqaman, D. Loerke, M. Mettlen, H. Kuwata, S. Grinstein, S. L. Schmid, and G. Danuser, “Robust singleparticle tracking in livecell timelapse sequences,” Nat. methods, vol. 5, no. 8, pp. 695–702, 2008.
 [5] M. Dewan, M. Ahmad, and M. Swamy, “Tracking biological cells in timelapse microscopy: an adaptive technique combining motion and topological features,” IEEE Trans. Biomed. Eng., vol. 58, no. 6, pp. 1637–1647, 2011.
 [6] I. Sbalzarini and P. Koumoutsakos, “Feature point tracking and trajectory analysis for video imaging in cell biology,” J. Struct. Biol., vol. 151, no. 2, pp. 182–195, 2005.
 [7] D. Padfield, J. Rittscher, and B. Roysam, “Coupled minimumcost flow cell tracking for highthroughput quantitative analysis,” Med. Image Anal., vol. 15, no. 4, pp. 650–668, 2011.
 [8] E. Meijering, I. Smal, and G. Danuser, “Tracking in molecular bioimaging,” IEEE Signal Proc. Mag., vol. 23, no. 3, pp. 46–53, 2006.
 [9] D. House, M. Walker, Z. Wu, J. Wong, and M. Betke, “Tracking of cell populations to understand their spatiotemporal behavior in response to physical stimuli,” in Proc. IEEE Conf. Comput. Vis. Mach. Intell. (CVPR 2009), 2009, pp. 186–193.

[10]
L. Liang, H. Shen, P. De Camilli, D. Toomre, and J. Duncan, “An expectation maximization based method for subcellular particle tracking using multiangle TIRF microscopy,” in
Medical Imag. Comput. ComputerAssist. Interv. (MICCAI 2011), 2011, pp. 629–636.  [11] N. H. Nguyen, S. Keller, E. Norris, T. T. Huynh, M. G. Clemens, and M. C. Shin, “Tracking colliding cells in vivo microscopy,” IEEE Trans. Biomed. Eng., vol. 58, no. 8, pp. 2391–2400, 2011.
 [12] A. Genovesio, T. Liedl, V. Emiliani, W. Parak, M. CoppeyMoisan, and J. OlivoMarin, “Multiple particle tracking in 3D+t microscopy: Method and application to the tracking of endocytosed quantum dots,” IEEE Trans. Image Process., vol. 15, no. 5, pp. 1062–1070, 2006.
 [13] L. Yang, Z. Qiu, A. Greenaway, and W. Lu, “A new framework for particle detection in lowSNR fluorescence livecell images and its application for improved particle tracking,” IEEE Trans. Biomed. Eng., vol. 59, no. 7, pp. 2040–2050, 2012.
 [14] N. Chenouard, I. Bloch, and J.C. OlivoMarin, “Multiple hypothesis tracking for cluttered biological image sequences,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 35, no. 11, pp. 2736–2750, 2013.
 [15] R. Juang, A. Levchenko, and P. Burlina, “Tracking cell motion using GMPHD,” in Proc. IEEE Int. Symp. Biomed. Imag. (ISBI 2009), 2009, pp. 1154–1157.
 [16] S. H. Rezatofighi, S. Gould, B.N. Vo, K. Mele, W. E. Hughes, and R. Hartley, “A multiple model probability hypothesis density tracker for timelapse cell microscopy sequences,” in Inform. Process. Medical Imag. (IPMI 2013), 2013, pp. 110–122.
 [17] T. Wood, C. Yates, D. Wilkinson, and G. Rosser, “Simplified multitarget tracking using the PHD filter for microscopic video data,” IEEE Trans. Circ. Syst. Vid., vol. 22, no. 5, pp. 702–713, 2012.
 [18] L. Feng, Y. Xu, Y. Yang, and X. Zheng, “Multiple dense particle tracking in fluorescence microscopy images based on multidimensional assignment,” J. Struct. Biol., vol. 173, no. 2, pp. 219–228, 2011.
 [19] S. H. Rezatofighi, S. Gould, R. Hartley, K. Mele, and W. E. Hughes, “Application of the IMMJPDA filter to multiple target tracking in total internal reflection fluorescence microscopy images,” in Medical Imag. Comput. ComputerAssist. Interv. (MICCAI 2012), 2012, pp. 357–364.
 [20] K. Li, E. D. Miller, M. Chen, T. Kanade, L. E. Weiss, and P. G. Campbell, “Cell population tracking and lineage construction with spatiotemporal context,” Med. Image Anal., vol. 12, no. 5, pp. 546–566, 2008.
 [21] I. Smal, E. Meijering, K. Draegestein, N. Galjart, I. Grigoriev, A. Akhmanova, M. Van Royen, A. Houtsmuller, and W. Niessen, “Multiple object tracking in molecular bioimaging by RaoBlackwellized marginal particle filtering,” Med. Image Anal., vol. 12, no. 6, pp. 764–777, 2008.
 [22] I. Smal, K. Draegestein, N. Galjart, W. Niessen, and E. Meijering, “Particle filtering for multiple object tracking in dynamic fluorescence microscopy images: Application to microtubule growth analysis,” IEEE Trans. Med. Imag., vol. 27, no. 6, pp. 789–804, 2008.
 [23] W. Godinez, M. Lampe, S. Wörz, B. Müller, R. Eils, and K. Rohr, “Deterministic and probabilistic approaches for tracking virus particles in timelapse fluorescence microscopy image sequences,” Med. Image Anal., vol. 13, no. 2, pp. 325–342, 2009.
 [24] R. Hoseinnezhad, B.N. Vo, B.T. Vo, and D. Suter, “Visual tracking of numerous targets via multiBernoulli filtering of image data,” Pattern Recogn., vol. 45, no. 10, pp. 3625–3635, 2012.
 [25] L. Yuan, Y. F. Zheng, J. Zhu, L. Wang, and A. Brown, “Object tracking with particle filtering in fluorescence microscopy images: Application to the motion of neurofilaments in axons,” IEEE Trans. Med. Imag., vol. 31, no. 1, pp. 117–130, 2012.
 [26] N. Chenouard, I. Smal, F. De Chaumont, M. Maška, I. F. Sbalzarini, Y. Gong, J. Cardinale, C. Carthel, S. Coraluppi, M. Winter et al., “Objective comparison of particle tracking methods,” Nat. methods, vol. 11, pp. 281––289, 2014.
 [27] I. Smal, M. Loog, W. Niessen, and E. Meijering, “Quantitative comparison of spot detection methods in fluorescence microscopy,” IEEE Trans. Med. Imaging, vol. 29, no. 2, pp. 282–301, 2010.
 [28] J. Burchfield, J. Lopez, K. Mele, P. Vallotton, and W. Hughes, “Exocytotic vesicle behaviour assessed by TIRFM,” Traffic, vol. 11, pp. 429–439, 2010.
 [29] R. P. Mahler, B.T. Vo, and B.N. Vo, “CPHD filtering with unknown clutter rate and detection profile,” IEEE Trans. Signal Process., vol. 59, no. 8, pp. 3497–3513, 2011.
 [30] M. Beard, B.T. Vo, and B.N. Vo, “Multitarget filtering with unknown clutter density using a bootstrap GMCPHD filter,” IEEE Signal Proc. Let., vol. 20, no. 4, pp. 323–326, 2013.
 [31] R. P. Mahler, Statistical multisourcemultitarget information fusion. Artech House Boston, 2007, vol. 685.
 [32] Y. BarShalom, Tracking and data association. Academic Press Professional, Inc., 1987.
 [33] S. S. Blackman, Multipletarget tracking with radar applications. Artech House, Inc.,, 1986.
 [34] C. Hue, J.P. Le Cadre, and P. Pérez, “Tracking multiple objects with particle filtering,” IEEE Trans. Aerosp. Electron. Syst., vol. 38, no. 3, pp. 791–812, 2002.
 [35] R. Mahler, “Multitarget Bayes filtering via firstorder multitarget moments,” IEEE Trans. Aerosp. Electron. Syst., vol. 39, no. 4, pp. 1152–1178, 2003.
 [36] ——, “PHD filters of higher order in target number,” IEEE Trans. Aerosp. Electron. Syst., vol. 43, no. 4, pp. 1523–1543, 2007.
 [37] B.T. Vo, B.N. Vo, and A. Cantoni, “The cardinality balanced multitarget multiBernoulli filter and its implementations,” IEEE Trans. Signal Process., vol. 57, no. 2, pp. 409–423, 2009.
 [38] B.N. Vo, B.T. Vo, N.T. Pham, and D. Suter, “Joint detection and estimation of multiple objects from image observations,” IEEE Trans. Signal Process., vol. 58, no. 10, pp. 5129–5141, 2010.

[39]
B.T. Vo and B.N. Vo, “Labeled random finite sets and multiobject conjugate priors,”
IEEE Trans. Signal Process., vol. 61, no. 13, pp. 3460–3475, 2013.  [40] R. P. Mahler, Advances in Statistical MultisourceMultitarget Information Fusion. Artech House Boston, 2014.
 [41] B.N. Vo, S. Singh, and A. Doucet, “Sequential Monte Carlo methods for multitarget filtering with random finite sets,” IEEE Trans. Aerosp. Electron. Syst., vol. 41, no. 4, pp. 1224–1245, 2005.
 [42] B.T. Vo, B.N. Vo, and A. Cantoni, “Analytic implementations of the cardinalized probability hypothesis density filter,” IEEE Trans. Signal Process., vol. 55, no. 7, pp. 3553–3567, 2007.
 [43] S. Pasha, B.N. Vo, H. Tuan, and W. Ma, “A Gaussian mixture PHD filter for jump Markov system models,” IEEE Trans. Aerosp. Electron. Syst., vol. 45, no. 3, pp. 919–936, 2009.
 [44] S. H. Rezatofighi, R. Hartley, and W. E. Hughes, “A new approach for spot detection in total internal reflection fluorescence microscopy,” in Proc. IEEE Int. Symp. Biomed. Imag. (ISBI 2012), 2012, pp. 860–863.
 [45] B.N. Vo and W. Ma, “The Gaussian mixture probability hypothesis density filter,” IEEE Trans. Signal Process., vol. 54, no. 11, pp. 4091–4104, 2006.
 [46] L. Lin, Y. BarShalom, and T. Kirubarajan, “Track labeling and PHD filter for multitarget tracking,” IEEE Trans. Aerosp. Electron. Syst., vol. 42, no. 3, pp. 778–795, 2006.
 [47] K. Panta, B.N. Vo, and S. Singh, “Novel data association schemes for the probability hypothesis density filter,” IEEE Trans. Aerosp. Electron. Syst., vol. 43, no. 2, pp. 556–570, 2007.
 [48] D. Schuhmacher, B.T. Vo, and B.N. Vo, “A consistent metric for performance evaluation of multiobject filters,” IEEE Trans. Signal Process., vol. 56, no. 8, pp. 3447–3457, 2008.
 [49] B. Ristic, B.N. Vo, D. Clark, and B.T. Vo, “A metric for performance evaluation of multitarget tracking algorithms,” IEEE Trans. Signal Process., vol. 59, no. 7, pp. 3452–3457, 2011.
 [50] S. H. Rezatofighi, W. T. E. Pitkeathly, S. Gould, R. Hartley, K. Mele, W. E. Hughes, and J. G. Burchfield, “A framework for generating realistic synthetic sequences of total internal reflection fluorescence microscopy images,” in Proc. IEEE Int. Symp. Biomed. Imag. (ISBI 2013), 2013.
 [51] L. Rayleigh, “Investigations in optics, with special reference to the spectroscope,” Philos. Mag., vol. 8, no. 49, pp. 261–274, 1879.
 [52] E. Meijering, “MTrackJ: A Java program for manual object tracking,” [Online].Available:http://www.imagescience.org/meijering/software/mtrackj/.
Comments
There are no comments yet.