(PDF) Tracking Multiple Particles in Fluorescence Time-Lapse Microscopy Images via Probabilistic Data Association

Background: Since 2019, the COVID-19 epidemic has spread across the world, affecting millions of individuals on physical, mental, and socioeconomic levels. Materials and methods: This qualitative study aims to map out the mental health consequences of the COVID-19 pandemic and suggest recovery strategies based on analyzing COVID-19 mental health and psychosocial impact-related comments on social/news media, participant’s observations, and participant action research. The results and suggested recovery strategies were triangulated with literature survey. Results: In COVID sufferers, neuropsychiatric symptoms lasting months were noted while some went on to have a relapse of psychological or neurological conditions. In the general population, common consequences included common mental health disorders, changing patterns of substance abuse particularly alcohol, and increased domestic violence. Family relationships and interactions generally improved but in a minority had strained dynamics. Vulnerable groups included the elderly, women, youth, children, disabled, frontline and health workers, minorities, and severely mentally ill. At the community level, there were signs of collective (mass) hysteria from panic buying to conspiracy theories, public shaming, fake news, and disinformation spreading on social media and mass protests. There were also positive effects such as better understanding of bio-knowledge, interventions and discoveries, online learning, gratitude toward health-care workers, increased community and family cohesion, reduction in road traffic accidents, and drop in other respiratory infections with their resultant morbidity and mortality. Conclusion: Given the widespread mental health consequences of the COVID-19 pandemic, a community-based approach is suggested while treating more severe mental disorders at the primary care or specialist level.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2359541, IEEE Transactions on Medical Imaging SUBMITTED TO IEEE TRANSACTIONS ON MEDICAL IMAGING 1 Tracking Multiple Particles in Fluorescence Time-Lapse Microscopy Images via Probabilistic Data Association William J. Godinez* Member, IEEE, and Karl Rohr Abstract—Tracking subcellular structures as well as viral structures displayed as ‘particles’ in fluorescence microscopy images yields quantitative information on the underlying dynamical processes. We have developed an approach for tracking multiple fluorescent particles based on probabilistic data association. The approach combines a localization scheme that uses a bottom-up strategy based on the spot-enhancing filter as well as a top-down strategy based on an ellipsoidal sampling scheme that uses the Gaussian probability distributions computed by a Kalman filter. The localization scheme yields multiple measurements that are incorporated into the Kalman filter via a combined innovation, where the association probabilities are interpreted as weights calculated using an image likelihood. To track objects in close proximity, we compute the support of each image position relative to the neighboring objects of a tracked object and use this support to re-calculate the weights. To cope with multiple motion models, we integrated the interacting multiple model algorithm. The approach has been successfully applied to synthetic 2D and 3D images as well as to real 2D and 3D microscopy images, and the performance has been quantified. In addition, the approach was successfully applied to the 2D and 3D image data of the recent Particle Tracking Challenge at the IEEE International Symposium on Biomedical Imaging (ISBI) 2012. Index Terms—Biomedical imaging, microscopy images, tracking, virus particles. I. I NTRODUCTION The dynamical behavior of subcellular as well as viral structures is the subject of intensive research. In time-lapse fluorescence microscopy, these structures are tagged with a fluorescent label. Since the optical system of typical microscopes has a limited spatial resolution, these small structures cannot be spatially resolved and instead the structures appear as ‘particles’ in the microscopy images. By tracking these particles we can extract detailed quantitative characterizations of the spatial-temporal processes. Since tracking should be performed for a large number of particles to obtain statistically Copyright (c) 2010 IEEE. Personal use of this material is permitted. However, permission to use this material for any other purposes must be obtained from the IEEE by sending a request to pubs-permissions@ieee.org. Manuscript received April 03, 2014; revised August 08, 2014. Asterisk indicates corresponding author. *W. J. Godinez is with the Department Bioinformatics and Functional Genomics, Biomedical Computer Vision Group, University of Heidelberg, BIOQUANT, IPMB, and DKFZ Heidelberg, Im Neuenheimer Feld 267, 69120 Heidelberg, Germany (email: wgodinez@ieee.org). K. Rohr is with the Department Bioinformatics and Functional Genomics, Biomedical Computer Vision Group, University of Heidelberg, BIOQUANT, IPMB, and DKFZ Heidelberg, Im Neuenheimer Feld 267, 69120 Heidelberg, Germany (email: k.rohr@dkfz.de). sound conclusions, automatic tracking approaches must be robust and efficient. In previous work for tracking fluorescent particles, deterministic approaches consisting of a localization step and a motion correspondence step have been proposed (e.g., [1]–[8]). Other deterministic schemes find minimal cost paths in the volume entailed by space and time (e.g., [9], [10]). Fitting spatialtemporal models is also a common approach (e.g., [11], [12]). While efficient, deterministic approaches ignore the spatial and temporal uncertainties involved in particle localization and position estimation over time. Thus, the approaches typically have problems dealing with challenging tracking scenarios (e.g., low signal-to-noise ratio, high object density). Probabilistic approaches formulated within a Bayesian framework cope with the spatial and temporal uncertainties by defining a posterior distribution on the variables describing the objects given a series of image-derived measurements. For probabilistic schemes that take into account the spatial uncertainty but not the temporal uncertainty see, for example, [13], [14]. For probabilistic schemes that take into account both the spatial and temporal uncertainty, the spatial-temporal posterior can be resolved via a Kalman filter or a particle filter. Approaches based on the Kalman filter (e.g., [15]–[18]) incorporate only the measurements computed by a bottom-up localization algorithm; thus localization and position estimation over time are uncoupled. Also, the Kalman filter typically considers only a single measurement (image position) per object. In contrast, the particle filter [19]–[24] queries directly multiple image positions (represented by random samples) to determine the location of an object (e.g. [25], [26]). The filter performs both localization and position estimation over time and achieves more accurate results (e.g., [21]). Because of the random nature of its top-down localization scheme, the filter uses a relatively large number of samples. This introduces a high computational cost because an image likelihood is evaluated for each sample. While one can improve the efficiency of the particle filter, (e.g., by exploring hierarchically the feature space, e.g., [19], [27], by biasing the samples towards regions with high likelihoods [28], or by marginalization, e.g., [29], [30]), several hundred samples are still required to ensure convergence. Probabilistic approaches have also been proposed for addressing the uncertainty involved in the task of motion correspondence [31]. These include approaches based on the joint probabilistic data association filter (JPDAF, e.g., [32]– [34]), the multiple hypothesis tracking algorithm (MHT, e.g., [35], [36]) as well as approaches sampling correspondences via 0278-0062 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2359541, IEEE Transactions on Medical Imaging SUBMITTED TO IEEE TRANSACTIONS ON MEDICAL IMAGING Monte Carlo schemes (e.g., [37]). These approaches explore thoroughly the space of correspondences over two or multiple time steps, which entails a high computational cost when tracking a very large number of objects. In this work, based on our previous work in [38], we introduce an efficient and robust approach for tracking multiple fluorescent particles based on probabilistic data association (PDA). We propose a localization scheme that combines bottom-up and top-down strategies. The bottom-up strategy relies on the spot-enhancing filter [39] while the top-down strategy generates measurements via an ellipsoidal sampling scheme based on the Gaussian probabilities calculated by the Kalman filter. Since we use an elliptical sampling scheme we refer to our approach as PDAE. The localization scheme generates multiple measurements which are integrated into the Kalman filter using the principle of a combined innovation as proposed in the PDA algorithm (e.g., [40], [41]). Unlike the standard PDA algorithm, our PDAE approach interprets the association probabilities of each measurement as weights relative to the image likelihood of the object. To calculate the weights, we synthesize hypotheses on the possible appearance and location of the object and test the hypotheses against the information found in the images. Thus, rather than using the combined innovation for addressing the correspondence problem, we use the combined innovation as a localization mechanism based on a recognition-by-synthesis scheme. Since the image likelihood considers only one object, an image region may support the measurements of multiple neighboring objects. Thus, multiple filters may converge towards the same image regions. To discourage this, we calculate the support of each image position relative to the neighboring objects of a tracked object and use this relative support to re-calculate the probabilities of each measurement of the tracked object. To incorporate multiple motion models, we combine our PDAE approach with the interacting multiple model (IMM) algorithm (e.g., [15], [42]). Analogous to the particle filter (e.g., [19]), our approach capitalizes on multiple measurements, uses the image data directly, and performs both localization and position estimation over time. Unlike previous linear schemes (e.g., [43], [44] for microtubuli) or radial schemes (e.g., [45] for white blood cells, or [46] for C. elegans worms) for localization in biological imaging, our ellipsoidal sampling scheme for generating measurements takes into consideration the directions along which the main deviations from the predicted position are expected. Also, the size of the ellipsoid is automatically determined. Moreover, our approach can be used for 2D and 3D images. In contrast to previous approaches using the standard PDA algorithm, our approach does not rely exclusively on bottom-up measurements (e.g., [33], [34]) nor does it rely on random sampling (e.g., [32], [41]). In contrast to approaches based on multiple hypothesis tracking (e.g., [35], [36]), the computational cost when dealing with multiple objects is at most quadratic relative to the number of objects (compared to exponential costs). In comparison to approaches based on the IMM algorithm (e.g., [15], [17], [34]), our approach uses a top-down localization strategy, and to compute the likelihood of the motion models, our approach does not assume Gaussian densities. Compared to our previous work 2 in [38], where measurements were generated only along the main axes of the ellipsoids, here measurements are generated along the contours of concentric ellipsoids, which yields better results; additionally we here use a different image likelihood as well as a different scheme for calculating the support relative to the neighboring objects, and we quantitatively evaluated the approach using synthetic images and 3D real image data. Based on synthetic 2D and 3D image sequences we have studied the performance as a function of the image noise and the object density. Our approach has been also successfully applied to 2D and 3D real images displaying HIV-1 particles. To the best of our knowledge, this is the first time that a tracking approach based on the PDA algorithm using a recognitionby-synthesis scheme underpinned by an ellipsoidal sampling scheme is introduced for tracking multiple fluorescent particles in 2D and 3D microscopy image sequences. Our approach generates bottom-up as well as top-down measurements for particle localization and operates in conjunction with either the Kalman filter or the IMM algorithm for position estimation. This paper is organized as follows. In Section II we outline the probabilistic framework underpinning our approach and briefly recall the standard Kalman filter. Section III presents the proposed localization scheme. The approach used for dealing with multiple objects is presented in Section IV. Section V describes the strategy for integrating multiple motion models. The models used for tracking are described in Section VI. Section VII presents the experimental results, and the conclusions are drawn in Section VIII. II. P ROBABILISTIC T RACKING A PPROACH A. Bayesian Framework In our approach, tracking is formulated as a Bayesian sequential estimation problem. Within a one-body state space, it is assumed that a fluorescent particle is represented by a state vector xt and that a noisy measurement yt reflects the true state of xt . At time step t, the aim is to estimate the state xt of a fluorescent particle given a sequence of measurements y1:t . By modeling the temporal evolution using a dynamical model p(xt |xt−1 ) and incorporating measurements derived from the images via a measurement model p(yt |xt ), a Bayesian filter estimates the posterior distribution p(xt |y1:t ) via stochastic propagation and Bayes’ theorem: Z p(xt |y1:t−1 ) = p(xt |xt−1 ) p(xt−1 |y1:t−1 ) dxt−1 (1) p(xt |y1:t ) ∝ p(yt |xt ) p(xt |y1:t−1 ). (2) An estimate of xt can be obtained from the posterior p(xt |y1:t ). For linear and Gaussian models, one can resolve analytically the posterior using a Kalman filter; for non-linear and/or non-parametric models, the particle filter provides a numerical solution. For reasons of efficiency and robustness, we use the Kalman filter with multiple measurements. B. Kalman filter The Kalman filter (e.g., [47], [48]) represents the posterior p(xt |y1:t ) via a Gaussian probability distribution that is 0278-0062 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2359541, IEEE Transactions on Medical Imaging SUBMITTED TO IEEE TRANSACTIONS ON MEDICAL IMAGING 3 parametrized by its mean vector mt and its covariance matrix Pt : p(xt |y1:t ) = N (·; mt , Pt ), (3) where N (·; µ, Σ) represents a Gaussian distribution with mean vector µ and covariance matrix Σ. The Kalman filter also assumes that the dynamical model p(xt |xt−1 ) as well as the measurement model p(yt |xt ) are linear and Gaussian: p(xt |xt−1 ) p(yt |xt ) = = N (xt ; Fxt−1 , Q) N (yt ; Hxt , R), (4) (5) where the transition matrix F and the measurement matrix H are known matrices. The covariance matrices Q and R encode the uncertainty about the prediction and the measurement, respectively. Given a posterior p(xt−1 |y1:t−1 ), the posterior at time step t can be computed by first carrying out a prediction on the state vector as well as on the associated covariance matrix: m̂ = Fmt−1 (6) P̂ = FPt−1 FT + Q. (7) To improve the readability we omit the time step t index in the following. The predicted state vector m̂ can be transformed onto the measurement space to obtain a predicted measurement ŷ: ŷ = Hm̂. (8) The predicted errors as encoded by the predicted covariance matrix P̂ can also be propagated onto the measurement space to obtain a predicted measurement covariance S, which also takes into account the measurement noise process as encoded by R: S = HP̂HT + R. (9) Once a measurement y is derived from the image at time step t, the difference between the actual measurement and the predicted measurement, i.e., the innovation ν, is calculated: ν = y − ŷ. (10) The Kalman gain K is computed as follows: K = P̂HT S−1 . (11) An estimate for the mean vector m is finally calculated by correcting the prediction m̂ with the innovation ν, where the correction is regulated by the Kalman gain K: m = m̂ + Kν. (12) The covariance matrix P can also be calculated by updating the predicted covariance matrix P̂ with the Kalman gain K: P = (I − KH)P̂, where I is the identity matrix. (13) III. M EASUREMENT P ROCESS A. Overview In biological imaging, tracking approaches based on the Kalman filter typically use only a single measurement y in (10) to update each filter (e.g., [15], [16], [35], [49]). In contrast, we consider multiple measurements (e.g., [32], [33]) to update each Kalman filter. These measurements are obtained by considering the image data (i.e., bottom-up localization) as well as by using a localization scheme that takes into account the predicted measurement ŷ (i.e., top-down localization). The measurements are assimilated by the Kalman filter via the combined innovation principle of the probabilistic data association (PDA) algorithm. Below we describe our schemes for bottom-up and top-down localization as well as the scheme for measurement integration via the PDA algorithm. In the following, we describe the measurement process considering only a single particle. The case of multiple particles is considered in Section IV below. B. Bottom-up localization via the Spot-Enhancing Filter For detecting and localizing a particle based on the image data only, any spot detection scheme may be used (see [50], [51] for an evaluation of spot detection schemes in fluorescence microscopy). In our case, we use an approach based on the spot-enhancing filter [39]. Briefly, this approach consists in first applying a Laplacian-of-Gaussian (LoG) filter, which is parametrized by the standard deviations σLoG,xy and σLoG,z , followed by applying a threshold to the filtered image to detect the particle. The threshold is automatically determined based on the mean intensity of the filtered image plus a factor c times the standard deviation of the filtered image intensities. Applying the spot detection scheme to the image ideally yields a single measurement ySEF that corresponds to the tracked particle. Note that this approach does not take into account the predicted measurement ŷ and instead exhaustively examines each position of the image to localize the particle. C. Top-down localization via Ellipsoidal Measurements Instead of exhaustively examining each position, the topdown localization scheme in our approach generates measurements in a neighborhood around the predicted measurement ŷ ∈ Y. In biological imaging (e.g., [15], [30], [33], [52]), typically the measurement space Y includes variables describing the position p of the particle in the image as well as variables describing the appearance (e.g., the peak intensity Imax of the particle). Certainly, one could generate measurements in the neighborhood of ŷ over all dimensions of the space Y. For efficiency reasons, in our approach we only generate position measurements p close to the predicted position p̂ within the position space P defined by the image. The predicted position p̂ can be obtained by projecting the predicted measurement ŷ to P via Φ: p̂ = Φŷ. (14) To define a neighborhood around p̂, first we define a submatrix Sp by transforming S onto P: Sp = ΦSΦT . 0278-0062 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. (15) This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2359541, IEEE Transactions on Medical Imaging SUBMITTED TO IEEE TRANSACTIONS ON MEDICAL IMAGING 4 Our approach generates measurements based on the predicted position distribution N (·; p̂, Sp ). The neighborhood about p̂ is described by the following ellipsoidal validation region Vp̂,Sp (γp ): 2 Vp̂,Sp (γp ) ≡ {p | (p − p̂)T S−1 p (p − p̂) ≤ γp }. (16) Within this region Vp̂,Sp (γp ) one finds position measurements p that entail a Mahalanobis distance to p̂ less than or equal to γp2 . Top-down position measurements within the ellipsoidal region Vp̂,Sp (γp ) can be generated via, for instance, proposing split measurements (e.g., [53]), or performing random sampling (e.g., [32], [41]). In our case we generate measurements by first diagonalizing Sp thus obtaining the semi-axes of the ellipsoidal region Vp̂,Sp (γp ) as: p (17) ri = γp λi ei , where λi and ei are the eigenvalues and eigenvectors of Sp , respectively. The ellipsoid Ep̂,Sp is explicitly described by its centroid p̂ and its semi-axes ri . The measurements are generated by taking positions pj,c along Nc concentric ellipsoidal contours centered at p̂: c pj,c = p̂ + Auj , (18) Nc canon where uj is a position along an ellipsoid E0,S in canonical p position whose semi-axes have the same magnitude as the semi-axes of Ep̂,Sp , A is a rotation matrix describing the orientation of Ep̂,Sp , and c = 1, 2, . . . , Nc is the concentric index. In 2D images, along a single ellipsoidal contour, we take Nj ellipsoidal positions uj that take the following parametric form: ! |r0 | cos 2πj Nj , (19) uj = |r1 | sin 2πj Nj where |r0 | ≥ |r1 |, and j = 1, 2, . . . , Nj . The rotation matrix A is a matrix whose columns are given by the eigenvectors ei : A ≡ e0 e1 , (20) where for the corresponding eigenvalues we have λ0 ≥ λ1 . In 3D images, the measurements pj,k,c are additionally indexed by the parameter k, and are defined as follows: c pj,k,c = p̂ + Auj,k . (21) Nc Here we take Nj Nk ellipsoidal positions uj,k which are defined as follows:   πk |r0 | cos 2πj Nj sin Nk  πk  uj,k =  |r1 | sin 2πj (22) Nj sin Nk  , πk |r2 | cos Nk where |r0 | ≥ |r1 | ≥ |r2 |, j = 1, 2, . . . , Nj , and k = 1, 2, . . . , Nk . The rotation matrix A is given as follows: A ≡ e0 e1 e2 , (23) where λ0 ≥ λ1 ≥ λ2 . Since the measurements are generated along ellipsoidal contours we refer to these measurements as ellipsoidal measurements. Our scheme thus generates Fig. 1. Elliptical measurements (dots) for a 2D anisotropic Gaussian distribution. Brighter intensities for the Gaussian distribution correspond to higher probabilities. Using γp = 3, Nc = 4, and Nj = 16, 64 measurements are generated. Fig. 2. Ellipsoidal measurements (dots) for a 3D anisotropic Gaussian distribution. Brighter intensities for the Gaussian distribution correspond to higher probabilities. Using γp = 3, Nc = 1, Nj = 16, and Nk = 8, 128 measurements are generated. Nc Nj Nk (in 2D, Nk = 1) measurements based on N (·; p̂, Sp ) within the position space P. As an example, Figures 1 and 2 show the ellipsoidal measurements obtained for sample 2D and 3D Gaussian distributions, respectively. The Kalman filter, however, expects measurements y within Y. For this reason we map the position measurements back to Y. Each i-th measurement pi is embedded into Y via the pseudoinverse Φ+ : yi,p = Φ+ pi . (24) The vector yi,p only includes the position information and must be supplemented with the additional variables (e.g., appearance variables) in Y. We extract these variables from the predicted measurement ŷ via a selection matrix Ψ and add them to the embedded vector yi,p to obtain measurement yi : yi = yi,p + Ψŷ. (25) Together with the predicted measurement ŷ the top-down localization strategy yields Nc Nj Nk +1 measurements in total. D. Ellipsoidal Measurements For Bottom-up Measurement The prediction generated at time step t by the Kalman filter does not take into account the image data of time step t (latest image data). Thus the ellipsoidal measurements based on the predicted position distribution N (·; p̂, Sp ) do not 0278-0062 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2359541, IEEE Transactions on Medical Imaging SUBMITTED TO IEEE TRANSACTIONS ON MEDICAL IMAGING 5 take into account the latest image data. The latest bottomup measurement ySEF encodes information about the latest image data. To take into account the latest image data, we generate ellipsoidal measurements based on the measurement distribution N (·; pSEF , Rp ), where pSEF = ΦySEF includes the position variables of the bottom-up measurement ySEF and Rp = ΦRΦT is a submatrix derived from the covariance matrix R that regulates the measurement noise process. For embedding the position measurements into the measurement space, we use the bottom-up measurement ySEF instead of the predicted measurement ŷ in (25). Thus, for N (·; pSEF , Rp ), together with the bottom-up measurement ySEF , we generate Nc Nj Nk + 1 additional measurements. Taking into account the previous Nc Nj Nk + 1 measurements generated for N (·; p̂, Sp ), our measurement process yields up to Nm = 2Nc Nj Nk + 2 measurements in total. Note that if the bottomup measurement ySEF is missing (e.g., due to image noise), N (·; pSEF , Rp ) is undefined and the corresponding ellipsoidal measurements cannot be determined. In that case only half the number of measurements is generated. Thus, regardless of the performance of the bottom-up particle localization scheme, our measurement process supplies the Kalman filter with at least Nm = Nc Nj Nk + 1 measurements at each time step. E. Measurement Integration via Probabilistic Data Association The probabilistic data association algorithm (e.g., [31], [40]) is an approach for solving the problem of motion correspondence, which entails determining one-to-one correspondences between the tracked object and the set of measurements Yt , {y1 , y2 , . . . , yNm } found within the validation region Vŷ,S (γ) (cf. (16)). A standard approach for solving this task is to select the measurement yi that is closest to the predicted measurement ŷ. This local nearest neighbor strategy may lead to incorrect correspondences in difficult tracking scenarios (e.g., low SNR), where the set Yt may miss the true measurement corresponding to the tracked object, and/or may include false measurements. Instead of committing to a single measurement (i.e., hard association), the probabilistic data association algorithm computes an association probability βi to quantify the degree of correspondence between the tracked object and measurement yi . To take into account all measurements within the validation region (i.e., soft association), a combined innovation νcomb is computed: νcomb = Nm X βi ν i (26) i=1 P Nm βi = 1, where β0 is the probability with νi = yi − ŷ and i=0 that none of the measurements corresponds to the tracked object. The combined innovation νcomb is used to perform the update of the prediction in the Kalman filter (see (12)). To compute the association probabilities βi , the standard PDA algorithm assumes that at most one measurement within Yt represents the true measurement corresponding to the tracked object and that all other measurements are false measurements. It is also assumed that the false measurements are independent and identically distributed (i.i.d.) and that they arise from a uniform distribution over the validation region Vŷ,S (γ). Additionally, the number of false measurements is assumed to follow a Poisson distribution. The exact formula for the association probabilities βi as used in the standard PDA algorithm is given in [40] (see also Appendix). In our case, because of our ellipsoidal sampling scheme for localization (cf. Section III-C), some of these fundamental assumptions of the PDA algorithm are not satisfied. For example, taking Yt as the set of measurements generated by our measurement process, the ‘false’ measurements in Yt do not follow a uniform distribution. In addition, since we generate measurements based on either the predicted measurement ŷ or based on the bottom-up measurement ySEF , the two respective sets of ellipsoidal measurements do not follow the same distribution. Also, the number of measurements generated by our scheme is constant and predictable, so the number of false measurements does not follow a Poisson distribution. A final consideration is that in the standard PDA algorithm the association probabilities only take into account the Mahalanobis distance between each measurement yi and the predicted measurement ŷ, i.e., the image data is not directly used in the computation of βi . For these reasons, we adopt a different strategy for computing βi . Instead of treating βi as association probabilities, we interpret them as weights that quantify the probability that the image intensities z conform to the Kalman-based prediction, as represented by the measurements raised with our measurement process. To compute the weights, we query the image likelihood p(z|x), and thus βi , p(z|x = xi ), where xi = Ξ(yi ) and the function Ξ(·) transforms its argument (measurement y) onto the state space X . Given the image intensities z within a region-of-interest (ROI) about the position encoded in x, the image likelihood p(z|x) is defined by the ratio: p(z|x) , po (z|x) . pb (z|b) (27) Here the image likelihood po (z|x) relative to the object is defined as: D(z, g(x))2 po (z|x) ∝ exp − , (28) 2σn2 where g(·) represents the image intensities synthesized via a parametric appearance model g(·), D(·) is the Euclidean distance, and σn describes the expected degree of noise. The image likelihood pb (z|b) relative to the background intensity is given as: D(z, b)2 pb (z|b) ∝ exp − , (29) 2σn2 where each element of b takes the value of the background intensity Ib . Note that a similar image likelihood is also used to compute the weights of samples in tracking approaches based on particle filters (e.g., [19], [21]). Once all weights have been evaluated with the image likelihood p(z|x), the weights βi , i = 1, . . . , Nm are normalized so that they sum to unity. Since in our approach the variables βi are computed using the principle of recognition-by-synthesis, the combined innovation is not used for the problem of motion correspondence (as in 0278-0062 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2359541, IEEE Transactions on Medical Imaging SUBMITTED TO IEEE TRANSACTIONS ON MEDICAL IMAGING the standard PDA algorithm) but rather assimilates multiple measurements into the Kalman filter taking into account the image information directly, i.e., segmentation of a particle is not necessarily required for estimating the position over time. IV. T RACKING M ULTIPLE F LUORESCENT PARTICLES To track multiple objects, we instantiate one Kalman filter per object. The measurement process as outlined in Section III considers only one object and therefore must be adjusted to accommodate the task of tracking multiple fluorescent particles. Concretely, two aspects of the measurement process need additional steps: first, for the bottom-up localization step in Section III-B, we assume that only one bottom-up measurement is generated by the spot detection scheme. This assumption is not satisfied within the context of tracking multiple fluorescent particles, and thus a motion correspondence step is required (see Section IV-A below). Second, we assume that modes in the image likelihood p(z|x) used to evaluate the probabilities β in the measurement integration step in Section III-E are caused by the tracked object only. Within the context of tracking multiple fluorescent particles, such modes may actually be caused by neighboring objects with a similar appearance. To minimize the influence of modes induced by other objects, we calculate the support of an image position relative to the neighboring objects of each tracked object (see Section IV-B). Note that since we use the combined innovation principle of PDA in our measurement process, it would be conceivable that the task of tracking multiple objects could be addressed via the joint probabilistic data association (JPDA) algorithm (e.g., [31], [32], [33]). In the JPDA algorithm, each association probability β is calculated over all possible global association hypotheses between the tracked objects and the measurements. This requires deterministically enumerating all possible global association hypotheses. Since our measurement scheme generates a large number of measurements, the JPDA algorithm would entail a high computational cost. Additionally, in the JPDA algorithm, the same assumptions concerning false measurements are made as in the PDA algorithm. As mentioned above, in our measurement process such fundamental assumptions are not met (cf. Section III-E). The JPDA algorithm also assumes a one-to-one correspondence between the tracked objects and the measurements. However, in our measurement process, the measurements generated with our ellipsoidal sampling scheme are already associated to the corresponding object (filter) and thus we have a many-to-one correspondence. Finally, we do not interpret the weight β as an association probability but rather as a weight computed based on an image likelihood. For these reasons, the JPDA algorithm is not well suited to deal with the task of tracking multiple particles within the context of our measurement process. A. Motion Correspondence When applied to an image displaying multiple fluorescent particles, the spot detection scheme yields a set YSEF , {ySEF,1 , ySEF,2 , . . . , ySEF,Nm,SEF } of Nm,SEF measurements. In our bottom-up localization scheme, we assume that the spot 6 detection scheme yields a single measurement per tracked object. To assign a single measurement to each filter, we find one-to-one correspondences between the set of bottomup measurements YSEF and the set Ŷ , {ŷ1 , ŷ2 , . . . , ŷNobj } of Nobj predicted measurements generated by Nobj Kalman filters. To solve this correspondence problem, we use a global nearest neighbor (GNN) approach [1] based on a graph-theoretical approach for the transportation problem. After finding oneto-one correspondences, predicted measurements (i.e., filters) in Ŷ may remain unmatched. In this case, the measurement process of that filter generates only measurements based on the predicted measurement. Also, if the filter remains unmatched for Nd consecutive time steps, the corresponding trajectory is terminated. Analogously, unmatched measurements in YSEF are used to initialize trajectories. B. Image Support Relative To Neighbors In our measurement process, the probabilities β are calculated by querying the image likelihood p(z|x). This likelihood is constructed under the assumption that only one object is visible in the image, i.e., independently from other objects. Under this assumption, a mode (or peak) in the image likelihood corresponds to the tracked object. Within the context of tracking multiple objects, a peak in the image likelihood may correspond to a neighboring object with a similar appearance. Thus there is a risk that a peak originating from a neighboring object may strongly influence the estimates of the Kalman filter, eventually leading the filter to converge on neighboring objects. To prevent this, the calculation of the probabilities β needs to also take into consideration the probabilities β of the neighboring filters. After calculating the corresponding image likelihoods of all filters, each filter j is associated with m a set of weighted position measurements {pjk = Φykj , βkj }N k=1 . Likewise, the neighboring filters Nb(j) around filter j allow constructing a set of measurements ∪i∈Nb(j) {pik , βki }. By considering all positions p up to pixel accuracy we define a neighborhood map: X X βli , (30) MNb(j) (p) = i i∈Nb(j) l∈Kp where Kpi = {k | pik = p} is the set of indices of measurements in ∪i∈Nb(j) {pik , βki } located at p. The map MNb(j) describes the support of each image position to neighboring filters. The support provided to neighboring objects Nb(j) directly influences the support provided to the j-th tracked object: the higher the support to the neighbors, the lower the support to the j-th tracked object. Thus each position p allows for a relative support Sj (p) = exp(−MNb(j) (p)) for object j. For reasons of robustness, we also take into account neighboring positions p′ around a given position p that entail a distance D(p, p′ ) smaller than a certain value δp . The support Sj (p) at a given position p is defined as:   X (31) MNb(j) (p′ ) log(θ(p, p′ )) Sj (p) = exp  p′ 0278-0062 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2359541, IEEE Transactions on Medical Imaging SUBMITTED TO IEEE TRANSACTIONS ON MEDICAL IMAGING where θ(p, p′ ) describes the spatial dependency between two image positions in general and is defined as follows [54]: D(p, p′ )2 ′ . (32) θ(p, p ) = 1 − exp − σp2 In (31), log(θ(p, p′ )) acts as a weighting factor over the neighboring positions. The parameter σp regulates the strength of the contribution of neighboring positions towards the calculation of the support. Since p′ takes positions over a neighborhood around p, and at some point it takes the same position as p, i.e., p = p′ , we have in such a case that θ(p, p′ ) = 0 and thus log(θ(p, p′ )) = −∞. For numerical reasons, we use a constant value of −5 instead for the logarithm in this case. The probability β̂kj for a measurement ykj of object j is computed as: β̂kj = p(z|x = Ξ(ykj ))Sj (Φykj ). (33) The weights β̂kj are subsequently normalized so that they sum to unity. This process is performed for all objects. Since the probabilities β̂ of the neighbors are eventually different from the original probabilities β that were used to calculate the support, we recalculate the neighborhood map using the new weights β̂ and recalculate the probabilities. We perform this process (i.e., recalculating the neighborhood map and calculating the probabilities) a certain number of times. In practice we found that the values of the weights become relatively stable after five iterations over all objects. If two objects are at the same position, i.e., the prediction for the two objects is the same p̂1 = p̂2 , the resulting support depends for either object to a large extent on the position of the elliptical samples as well as on the behavior of the image likelihood p(z|x) of each object. Often the object with the larger values for the image likelihood gets larger values for the image support. However, if the image likelihood of the two objects is very similar, the scheme also allows, in principle, for the two objects to share the same position. V. M ULTIPLE M OTION M ODELS To better describe the motion of particles that alternate between distinct motion patterns, we use the interacting multiple model (IMM) algorithm (e.g., [42], [55], [15]). In the IMM algorithm, the motion of a particle is described by a finite set Ω of NΩ motion models. One Kalman filter is instantiated per motion model ω ∈ Ω. It is assumed that the dynamical model of a particle is governed by a Markov chain with transition matrix Π. The IMM algorithm recursively computes a model-conditioned posterior density p(xt |ωt , y1:t ) as well as a posterior model probability P (ωt |y1:t ) (for details see, e.g., [55]). In the standard IMM algorithm, a single measurement yt is used to update all NΩ filters. Here, we use the PDAE approach to feed multiple measurements to each filter: each i-th Kalman filter corresponding to the i-th motion model, i ∈ Ω, raises m its own set {yk , βk , ωk = i}N k=1 of ellipsoidal measurements. Note that each measurement is associated with the motion model ωk = i that led to such a measurement. To calculate the likelihood p(yt |ωt = i, y1:t−1 ) of each i-th model, we take the unnormalized values of the weights βk of each i-th filter 7 and normalize them with respect to the sum of all weights of all NΩ filters: βk βk′ = PNΩ P c j∈Υc βj , k ∈ Υi (34) where Υi = {j | βj , ωj = i} indexes the weights of the measurements of the i-th filter. The likelihood p(yt |ωt = i, y1:t−1 ) is given as the sum over its normalized weights βk′ : p(yt |ωt = i, y1:t−1 ) = X βj′ (35) j∈Υi For motion correspondence, we select for each object the filter with the highest predicted probability. This strategy leads to Nobj predicted measurements and so the correspondence scheme described in Section IV-A above is still applicable. For the calculation of the image support (see Section IV-B), we treat all NΩ Nobj filters as independent filters. This allows the NΩ filters of a certain object to compete for image positions among themselves, which helps to select the best motion model among all NΩ motion models relative to the predicted positions of neighboring objects. VI. M ODEL D EFINITIONS In this section, we define models for tracking fluorescent particles. Since in our application we deal with virus particles, we define models suitable for such particles. A. Appearance Model In our approach the intensities of each particle are represented by a 2D or 3D Gaussian function. The Gaussian function is parametrized by the position p = (xp , yp )T or p = (xp , yp , zp )T of the particle in a 2D or 3D image, respectively, by the peak intensity Imax , as well as by the standard deviations σxy (2D case) or σxy , σz (3D case). For position estimation, we consider the velocity vector ṗ = (ẋ, ẏ)T or ṗ = (ẋ, ẏ, ż)T . These parameters are encoded in the state vector x = (x0 , ẋ, y0 , ẏ, Imax , σxy )T (2D case) or x = (x0 , ẋ, y0 , ẏ, z0 , ż, Imax , σxy , σz )T (3D case). For 2D images, the appearance model of a single particle is given by: gGauss2D (x, y; x) = Ib + (x − xp )2 + (y − yp )2 , (Imax − Ib ) exp − 2 2σxy (36) where Ib denotes the background intensity. For 3D images, the appearance model is described by: gGauss3D (x, y, z; x) = Ib + (x − xp )2 + (y − yp )2 (z − zp )2 (Imax − Ib ) exp − . − 2 2σxy 2σz2 (37) 0278-0062 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2359541, IEEE Transactions on Medical Imaging SUBMITTED TO IEEE TRANSACTIONS ON MEDICAL IMAGING 8 B. Motion Model We use a random walk model as well as a directed motion model with constant velocity to describe the motion of a single particle (NΩ = 2). For both models, we assume that the appearance parameters Imax , σxy (2D case) or Imax , σxy , σz (3D case) follow Gaussian random walks with small perturbations. For the random walk model (ω = 1), the dynamical model is given in 2D by F1 = diag(1, 1, 1, 1, 1, 1) and the covariance matrix Q1 = diag(qx , qẋ , qy , qẏ , qImax , qσxy ). In 3D, the dynamical model is defined as F1 = diag(1, 1, 1, 1, 1, 1, 1, 1, 1) and the covariance matrix is Q1 = diag(qx , qẋ , qy , qẏ , qz , qż , qImax , qσxy , qσz ). For the directed motion model (ω = 2), the dynamical model in 2D is defined as:   1 ∆t 0 0 0 0 0 1 0 0 0 0   0 0 1 ∆t 0 0  F2 =  (38) 0 0 0 1 0 0 ,   0 0 0 0 1 0 0 0 0 0 0 1 while the covariance matrix takes the following form:  q ∆t3 q ∆t2  CV CV 0 0 0 0 3 2 2  qCV ∆t  qCV ∆t 0 0 0 0   2   3 2 qCV ∆t qCV ∆t  0 0 0  3 2 2 Q2 =  0 . qCV ∆t  0 qCV ∆t 0 0  0   2  0 0  0 0 0 qImax 0 0 0 0 0 qσxy (39) The dynamical model for directed motion in 3D is defined analogously. C. Measurement Model ′ ′ )T (2D case) , σxy We have measurements y = (x′ , y ′ , Imax ′ ′ ′ ′ ′ ′ T or y = (x , y , z , Imax , σxy , σz ) (3D case). For both models, the measurement matrix is given in 2D by:   1 0 0 0 0 0 0 0 1 0 0 0  (40) H= 0 0 0 0 1 0 . 0 0 0 0 0 1 The uncertainties of the measurements in 2D are defined by R = diag(rx , ry , rImax , rσxy ). In 3D, H and R are defined analogously. The function Ξ(·) is given by: Ξ(y) = H+ y (41) where H+ = HT . The matrix Φ used to project a 2D measurement y onto a 2D position space is defined as: 1 0 0 0 Φ= (42) 0 1 0 0 with the pseudoinverse matrix is given by:  1 Φ = 0 0 Φ+ = ΦT . In 3D, the projection 0 1 0 0 0 0 0 1 0  0 0 0 0 . 0 0 (43) The selection matrix Ψ in the 2D case is given by Ψ = diag(0, 0, 1, 1) while in 3D we have Ψ = diag(0, 0, 0, 1, 1, 1). The transition matrix Π is given by: 0.75 0.36 . (44) Π= 0.25 0.64 The initial probabilities for the motion models are set to T P (ω) = [ 0.6 0.4 ] . We use the Lucas-Kanade optical flow approach [56] to initialize the values of the velocity components. VII. E XPERIMENTAL R ESULTS We have evaluated the proposed PDAE approach using synthetic as well as real image sequences displaying virus particles. A. Experimental Procedures We compare experimentally the proposed PDAE approach with a standard Kalman filter as well as with an approach based on independent particle filters (IPF) [21]. The Kalman filter as well as the IPF use a random walk model in 2D and 3D (see [21] for details). For the PDAE approach we also use a random walk model or a random walk model as well as a directed motion model. The Kalman filter is updated with a single bottom-up measurement as computed by the approach based on the spot-enhancing filter (SEF) (cf. Section III-B) and uses a global nearest neighbor (GNN) scheme for motion correspondence (cf. Section IV-A). For the IPF, the measurement model is defined by the likelihood p(z|x) (27). As in the PDAE approach, the appearance model used by the IPF is also given by a 2D or 3D Gaussian function, cf. (36) and (37), respectively. Within the IPF approach, each independent particle filter uses a bottom-up measurement computed by the SEF approach to generate samples [26], in which case the proposal distribution q(·) of the particle filter is defined by the dynamical model and a Gaussian centered at the bottom-up detection with a covariance matrix. The bottomup measurement is assigned using the same global nearest neighbor scheme for motion correspondence. For numerical stability, we use Ns = 1000 (2D case) or Ns = 8000 (3D case) samples for each independent particle filter. As a comparison, for the PDAE approach, we use Nc = 4 concentric contours and Nj = 16 positions along each elliptical contour. In 3D we additionally evaluate concentric elliptical contours at Nk = 8 positions in z-direction. Together with the ellipsoidal measurements generated using the bottom-up measurement we have Nm = 130 measurements (samples) per filter in total (2D case) or Nm = 1026 measurements (3D case). Table I summarizes the evaluated approaches and their properties. To quantitatively evaluate the performance of the tracking approaches, we computed the tracking accuracy Ptrack : ntrack,correct Ptrack = , (45) ntrack,total where ntrack,correct is the number of correctly computed trajectories and ntrack,total represents the number of true trajectories. Note that Ptrack ∈ [0, 1]. We determine correspondences between true and computed trajectories using a nearest neighbor 0278-0062 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2359541, IEEE Transactions on Medical Imaging SUBMITTED TO IEEE TRANSACTIONS ON MEDICAL IMAGING 9 TABLE I S UMMARY OF THE EVALUATED APPROACHES . Bottom-up detection Motion correspondence Number of measurements Measurement source Search strategy Measurement evaluation Kalman filter IPF PDAE SEF GNN 1 SEF Exhaustive – SEF GNN Ns q(·) Random p(z|x) SEF GNN Nm N (·; p̂, Sp ), N (·; pSEF , Rp ) Ellipsoidal p(z|x) approach that starts at the position of the first time step of a true trajectory and finds the closest position of a computed trajectory at the same time step. The maximum distance is set to 5 pixels. If no such computed position exists, the search is carried out in subsequent time steps until a matching position is found. If this is the case, it is assumed that there is a correspondence between the true and the computed trajectory. Then the number of time steps is counted where the error between the corresponding true and computed positions is below the maximum distance. If this condition does not hold, the computed trajectory no longer matches the true trajectory. The true trajectory may be subsequently matched to another computed trajectory and the count for the number of time steps for which the true trajectory was correctly matched is resumed. Based on this count, for each true trajectory the percentage of correctly computed time steps is calculated, which is given by the number of time steps for which the true trajectory was correctly matched divided by the total number of time steps within the true trajectory. Ideally each true trajectory is matched by a single computed trajectory. If a true trajectory is matched with more than one computed trajectory, the percentage of correctly computed time steps is multiplied by a weight following a Gaussian function: the larger the number of broken trajectories, the lower the weight. The expected number of matching trajectories is 1, so the mean of the Gaussian function is set to 1. The deviation from the expected number of matching trajectories is expected to be low, so the standard deviation of the Gaussian function is set to 1. The number of correctly computed trajectories ntrack,correct is given as the sum over all weighted percentages of correctly computed time steps. A computed trajectory may be matched at most once with a ground truth trajectory (see also [21]). To evaluate the localization accuracy, we also calculated the root mean square error (RMSE) between the computed positions pt,comp and the true positions pt,true of a true trajectory whose positions are indexed in time by T = {. . . , t − 1, t, t + 1, . . .} spanning |T | = Nsteps,true time steps: s X 1 |pt,comp − pt,true |2 (46) RMSE = Nsteps,true t∈T B. Evaluation on Synthetic Images We evaluate the performance of the PDAE approach using synthetic image sequences and consider two experimental scenarios. In the first scenario, we evaluate the localization accuracy of the PDAE approach as a function of the image noise. In the second scenario, we evaluate the tracking accuracy of the approach as a function of the density of the objects. In the synthetic images we render each object using either a 2D Gaussian function (36) or a 3D Gaussian function (37) as the appearance model of the rendered objects. 1) First Synthetic Scenario: In this scenario we evaluate the localization accuracy of the approaches as a function of the image noise. We evaluate two noise models: a Poisson noise model and an EMCCD noise model. a) Poisson Noise Model: Here we assume a Poisson distribution for the noise model. The image noise is reflected by the signal-to-noise ratio (SNR), which is defined as the difference between the peak intensity Imax of an object and the intensity of the background Ib , divided by the standard deviation of the noise level σn ([57]). We set the background intensity to Ib = 10 and vary the peak intensity Imax to explore different SNR levels. We explore the following SNR levels: 11.6, 8.8, 6.5, 4.6, 3.5, 2.8, 2, and 1.3. For each level we generate 30 image sequences. In 2D, each image sequence has 50 time steps. The image dimensions are 64×64 pixels. In 3D, each image sequence has 30 time steps and the image dimensions are 64×64×32 voxels. Each image sequence displays one object. The object is initially positioned at the center of the image and is visible throughout the entire image sequence. The motion of the object is governed by random walk, where the expected magnitude of the displacement over two time steps is set to 1.3 pixels. The true trajectory of the object is the same over all SNR levels. Thus for each SNR level only the (random) noise configuration varies over the corresponding 30 image sequences. For the appearance model, we use a value of σxy = 2 pixels and a value of σz = 1 pixel. The parameters of the tracking approaches are empirically defined as follows. For the detection scheme based on the SEF, we adjust the factor c according to the SNR level and use a standard deviation of σLoG,xy = 1.5 pixels in 2D and a standard deviation of σLoG,z = 1.5 voxels in 3D. For all tracking approaches we use a random motion model. The expected deviations for the position variables are fixed to qx = qy = qz = 4 pixels2 . For the appearance variables, we use the following values: qImax = 1 intensity level unit2 , qσxy = qσz = 0.01 pixels2 . Values for these parameters may be also obtained via, e.g., maximum likelihood estimates (MLE) from sample trajectories. The measurement errors take the following values: rx = ry = rz = 0.01 pixels2 , rImax = 1 intensity level unit2 , qσxy = qσz = 0.36 pixels2 . The threshold for the validation region is set to γp2 = 5.99 in 2D and 0278-0062 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2359541, IEEE Transactions on Medical Imaging SUBMITTED TO IEEE TRANSACTIONS ON MEDICAL IMAGING γp2 = 7.82 in 3D. The expected degree of noise σn is set according to the SNR level. In this scenario we examine the performance of the PDAE approach in comparison to the IPF approach using the same number of samples Ns as the number of measurements Nm in the PDAE approach, viz. Ns = 130 in 2D and Ns = 1026 in 3D. We refer to these approaches as IPF-130 (2D case) and IPF-1026 (3D case), respectively; we refer to the approaches with higher number of samples as IPF-1000 (2D case) and IPF-8000 (3D case). We also examine the performance of the Kalman filter using 2D Gaussian fitting [58] for particle localization. In this scenario (but only in this one), for each approach a single trajectory is initialized throughout the image sequence. The initial position is the true position of the particle at the first time step of the image sequence. This is done to ensure that all approaches have the same initial conditions and also to prevent the influence of erroneous initializations caused by erroneous detections, especially at the lower SNR levels. In this scenario (but only in this one), the trajectory of each approach is terminated at the final time step of the image sequence. This allows quantifying how far away the computed trajectories are from the true positions over time. Thus the number of time steps for the computed trajectories of all approaches is the same as the number of time steps in the ground truth trajectory. In this scenario only, we allow the distances between the true and computed positions to be arbitrarily large. Sample results for all approaches at SNR = 1.3 are shown in Figure 3 and Figure 4 for 2D and 3D images, respectively. It can be seen that the accuracy of the PDAE approach is superior to that of the Kalman filter at this challenging SNR level. The performance of the PDAE approach is quite similar to that of the particle filter using a large number of samples. Figures 5 and 7 show the results of all evaluated approaches in terms of their RMSE as a function of the SNR level. The mean values of the RMSE computed over 30 image sequences corresponding to each SNR level are shown. The standard Kalman filter using Gaussian fitting obtains good results at the lower SNR levels, since it can better discriminate between true and false detections. For both 2D and 3D images, the PDAE approach outperforms generally the standard Kalman filter using either the spotenhancing filter or 2D Gaussian fitting for localization. The bottom-up detection schemes rely solely on the information found in a single image to localize the object. In comparison, the PDAE measurement process is based on a probability distribution computed sequentially over consecutive images by a spatial-temporal filter. Since the PDAE measurement process exploits the spatial-temporal information encoded in the image sequence, it yields improved results in comparison to using bottom-up particle localization schemes only. For the approach based on standard Kalman filters and using the spot-enhancing filter for bottom-up localization, the performance is influenced by the factor c used for calculating the detection threshold (cf. Section III-B). To investigate the influence of this parameter on the tracking performance, we lower the parameter value used originally for calculating the performance in Figure 5 by a certain percentage u = 10, 25, 50%. Figure 6 shows the results. It can be seen that for u = 10% the performance 10 is nearly unchanged. For u = 25% the error for SNR = 2 is slightly lower, however, at lower SNR levels the error is higher because of the larger number of false detections (due to the lower threshold). For u = 50% the performance is significantly worse. Thus c should be chosen conservatively so that the number of false detections is not too large. For both 2D and 3D images, the PDAE approach outperforms the particle filter using the same number of samples (Ns = 130 in 2D and Ns = 1026 in 3D). This indicates that for the IPF a relatively low number of samples is not sufficient to obtain an accurate numerical approximation of the posterior probabilities. The PDAE approach is not subject to such numerical issues since the posterior probabilities are represented analytically via the underlying Kalman filter. Improved results are obtained for the IPF by using a larger number of samples (i.e., Ns = 1000 in 2D or Ns = 8000 in 3D). At the higher SNR levels, the PDAE approach outperforms the IPF with a larger number of samples. A reason for this is that the proposed ellipsoidal sampling scheme explores the position space systematically while the IPF explores this space randomly. The systematic approach yields a comparable performance relative to that of the random scheme with a large number of samples at the lower SNR levels. Since both approaches query the same image likelihood, the computation time of both approaches is linearly dependent on the number of samples (measurements). Thus the computational cost entailed by the PDAE approach is ca. 7.7 times lower than that of the IPF. This is particularly beneficial when tracking a large number of objects. b) EMCCD Noise Model: We have also evaluated the approaches based on synthetic images generated using an EMCCD noise model [59], [60]. Here, the shot noise follows a Poisson distribution while the excess noise follows a Gaussian distribution. The two noise contributions are sequentially embedded into the images. We used the same variance for the two distributions. Then the noise variance σn2 that determines the SNR is given by the sum of the variances of the two distributions. Otherwise, the 2D image sequences are generated analogously as in Section VII-B1a above. To cope with the additional noise contribution, the expected variance σn2 used in the image likelihood (27) is increased according to the sum of the variances of the two noise distributions. Otherwise the parameters of the tracking approaches are the same as those in Section VII-B1a. We compared the PDAE approach with the standard Kalman filter approach as well as the different particle filter approaches. The results in terms of the RMSE as a function of the SNR levels 8.2, 6.3, 4.6, 3.2 are shown in Figure 8. It can be seen that the performance is similar to the results for the synthetic data generated with a Poisson noise model only (cf. Figure 5 and note the different ranges of the SNR levels). Thus, also for EMCCD noise the PDAE approach yields a better result than the other approaches. 2) Second Synthetic Scenario: We evaluate here the performance of the PDAE approach as a function of the object density. The object density dobj is given by the ratio between the number of objects Nobj and the size (area) of the image Simage . Each object occupies a certain area Sobj . The definition 0278-0062 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2359541, IEEE Transactions on Medical Imaging SUBMITTED TO IEEE TRANSACTIONS ON MEDICAL IMAGING Original image t = 34 Original image t = 35 Original image t = 36 Ground truth t = 34 Ground truth t = 35 Ground truth t = 36 Kalman t = 34 Kalman t = 35 Kalman t = 36 Kalman (GaussFit) t = 34 IPF-130 t = 34 Kalman (GaussFit) t = 35 IPF-130 t = 35 11 Original image Ground truth Kalman IPF-1026 IPF-8000 PDAE Kalman (GaussFit) t = 36 IPF-130 t = 36 Fig. 4. Original image, ground truth, and tracking results for the evaluated approaches on synthetic 3D images generated using a Poisson noise model. The SNR is 1.3 and the time step is t = 12. A z-slice (z = 13) of the original volume image as well as the ground truth are shown. Trajectories are rendered as spheres (positions) and sticks (displacement vectors); the cube shows the current position. 1.6 Kalman Kalman (GaussFit) IPF−130 IPF−1000 PDAE 1.4 IPF-1000 t = 35 IPF-1000 t = 36 1.2 RMSE [pixel] IPF-1000 t = 34 1 0.8 0.6 0.4 0.2 0 0 PDAE t = 34 PDAE t = 35 PDAE t = 36 Fig. 3. Original image, ground truth, and tracking results for the evaluated approaches on synthetic images generated using a Poisson noise model. The SNR is 1.3 and the time steps are t = 34, 35, 36. The arrows show the direction of the displacement vector relative to the position at the previous time step. The tip of the arrow indicates the current position. A dot is shown instead of an arrow if the position over two consecutive time steps is the same. Image intensities have been inverted. 2 4 6 8 10 12 SNR Fig. 5. Tracking accuracy (RMSE) as a function of the SNR for 2D synthetic image sequences generated using a Poisson noise model. The mean values computed over 30 image sequences at each SNR level are shown. 0278-0062 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2359541, IEEE Transactions on Medical Imaging SUBMITTED TO IEEE TRANSACTIONS ON MEDICAL IMAGING 12 of dobj , however, does not take into account Sobj . To take into account the area occupied by each object, we calculate the probability of overlap poverlap [13], which is defined as: 3 Kalman u=0% Kalman u=10% Kalman u=25% Kalman u=50% RMSE [pixel] 2.5 poverlap = 2 1.5 1 0.5 0 2 4 6 8 10 12 SNR Fig. 6. Tracking accuracy (RMSE) as a function of the SNR for 2D synthetic image sequences generated using a Poisson noise model. The mean values computed over 30 image sequences at each SNR level are shown for the approach using Kalman filters. The detection parameter c is lowered by a certain percentage u. 3 Kalman IPF−1026 IPF−8000 PDAE RMSE [voxel] 2.5 2 1.5 1 0.5 0 0 2 4 6 8 10 12 SNR Fig. 7. Tracking accuracy (RMSE) as a function of the SNR for 3D synthetic image sequences using a Poisson noise model. The mean values computed over 30 image sequences at each SNR level are shown. 1.6 Kalman IPF−130 IPF−1000 PDAE 1.4 RMSE [pixel] 1.2 1 0.8 0.6 0.4 0.2 0 3 4 5 6 7 8 9 SNR Fig. 8. Tracking accuracy (RMSE) as a function of the SNR for 2D synthetic image sequences generated using an EMCCD noise model. The mean values computed over 30 image sequences at each SNR level are shown. Nobj Sobj . Simage (47) Thus the higher the number of objects, the higher the probability of overlap among objects, and therefore the more tracking errors occur. In this scenario we generate 2D image sequences consisting of 30 time steps. The images have dimensions 100×100 pixels (16-bit) and thus Simage = 10000 pixels2 . We use the following number of objects: Nobj = 10, 20, 30, 40, 50, 75, and 100. We use the same appearance parameters as in the first scenario for rendering all objects. Assuming that each Gaussian object with σxy = 2 pixels can be described by a 2D disk with radius r = 1.5σxy = 3 pixels, the area of each object is given by Sobj = r2 π = 28 pixels2 and thus we explore the following values for the probability of overlap poverlap : 0.03, 0.06, 0.08, 0.11, 0.14, 0.21, 0.28. The SNR is 11.6. The initial image position of each object is random and the motion is governed by random walk. The expected deviation for the displacement takes values of up to 1.3 pixels. For each number of objects Nobj , we generate 30 image sequences. The factor c for the detection scheme based on the SEF is adjusted according to the number of objects, since the image statistics (mean intensity and standard deviation) vary with the number of objects. For all tracking approaches, the expected deviations for the position variables are fixed to qx = qy = qz = 0.56 pixels2 . Trajectories are initialized only at the first time step (only in this scenario). All other parameters are the same as in the first synthetic scenario. Results for all approaches at a probability of overlap poverlap = 0.11 (number of objects Nobj = 40) are presented in Figure 9. The results show that the PDAE approach preserves the identity of neighboring objects relatively well. The performance of the approaches in terms of Ptrack is shown in Figure 10. The performance follows approximately the theoretical detection performance 1 − poverlap [13]. To further examine the performance, we also calculate the Jaccard similarity measure [35] for the positions. The results are shown in Figure 11. For both measures, it can be seen that the PDAE approach consistently outperforms the other approaches. This shows the viability of the image support scheme that computes the support at each image position for tracking multiple objects. We also measured the computation time of each approach. All approaches were implemented in Java within our software ViroTracker [21] which was executed on an Intel Xeon X5650 (6 cores at 2.6 GHz) machine running Linux. Figure 12 shows the computation time as a function of the number of objects (mean computation time over 30 2D image sequences with 30 time steps each.) The approach with the lowest computation time is the Kalman filter. For the IPF approach and the PDAE approach, the computational cost is linearly dependent on the number of samples (measurements) per object. In comparison to the Kalman filter, the PDAE approach is computationally somewhat more expensive (approximately a factor of 3). In comparison to the IPF with 1000 samples, the PDAE approach 0278-0062 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2359541, IEEE Transactions on Medical Imaging SUBMITTED TO IEEE TRANSACTIONS ON MEDICAL IMAGING 13 1 Kalman IPF−1000 PDAE Jaccard similarity 0.95 0.9 0.85 0.8 0.75 0.7 Ground truth 0.65 0 Kalman 0.05 0.1 p 0.15 0.2 0.25 overlap Fig. 11. Tracking accuracy (Jaccard similarity) as a function of the probability of overlap of the objects poverlap . The mean values computed over 30 image sequences for each probability of overlap are shown. 50 45 Kalman IPF−1000 PDAE 40 IPF-1000 PDAE Fig. 9. Ground truth and tracking results for all evaluated approaches on synthetic images. The probability of overlap of the objects is poverlap = 0.11 and the time step is t = 27. The numbers label each trajectory. Insets show the tracking results of the region enclosed with a rectangle at a higher resolution. Time [s] 35 30 25 20 15 10 5 0 0 1 Kalman IPF−1000 PDAE 0.95 40 Nobj 60 80 100 Fig. 12. Computation time (s) as a function of the number of objects. The mean values computed over 30 2D image sequences consisting of 30 time steps for each number of objects are shown. Ptrack 0.9 20 0.85 0.8 0.75 0 0.05 0.1 p 0.15 0.2 0.25 overlap Fig. 10. Tracking accuracy (Ptrack ) as a function of the probability of overlap of the objects poverlap . The mean values computed over 30 image sequences for each probability of overlap are shown. using 130 measurements is about 7.7 times faster, as expected. For example, for Nobj = 40, the PDAE approach entails a mean computation time of 2.5s while the IPF requires about 20s. Thus, in comparison to the Kalman filter, the PDAE approach yields a superior performance at the expense of a somewhat higher computational cost. In comparison to the IPF, the PDAE approach delivers also a better performance in terms of Ptrack at a significantly lower computational cost. C. Evaluation on Real Microscopy Images 1) 2D Images: We have applied our approach to real microscopy 2D image sequences displaying fluorescently labeled HIV-1 particles. The HIV-1 particles are double labeled with GFP and RFP, pseudotyped with the glycoprotein of the vesicular stomatitis virus (VSV-G), and incubated with HeLa cells [61]. A Zeiss Axiovert 200 M microscope with a Roper Scientific Cascade II EM-CCD was used to acquire the images. The recording frequency is 10 Hz. The image sequences consist of 150 up to 400 images (256×256 or 512×512 pixels; 16-bit). We consider nine image sequences. The ground truth for the positions of the particles was determined manually with the commercial package MetaMorph. Between 5 and 43 particles were manually tracked in each image sequence. For the spot detection scheme we use a standard deviation of σLoG,xy = 1.5 pixels. For the adaptive threshold, we use a factor of c = 2.7 for the first four image sequences; for the other image sequences we use c = 3.7. For the PDAE approach we use a random walk model as well as a directed motion model. For the random walk motion model, we use a value of qx = qy = qz = 4 pixels2 . For the directed motion model, we use a value of qCV = 4 pixels2 /frame3 . The noise parameter σn is adjusted based on the measured peak intensity for the particles. Trajectories are initialized and terminated at any time step. All other parameters are set as in the first synthetic scenario (cf. Section VII-B1). In this scenario we also present the results obtained with the PDAE approach using only a random walk (RW) motion model. Additionally 0278-0062 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2359541, IEEE Transactions on Medical Imaging SUBMITTED TO IEEE TRANSACTIONS ON MEDICAL IMAGING we perform a comparison with the ParticleTracker approach [1], the uTrack approach [2], and the multiple-hypothesis tracking (MHT) approach using Kalman filters and multiple motion models [35] as implemented in the Icy software [62]. For all approaches, different parameter values were tested and we apply the values which yielded the best results. For the ParticleTracker approach, we use a radius of 2 pixels, a cutoff of 3 and a percentile of 1 for the first four image sequences; for the other image sequences we set the percentile to 0.3. The linking range is set to 4 frames for the first image sequences and 2 frames for the other image sequences; the displacement is set to 1.73 pixels for the first four image sequences and 4 pixels for the other image sequences. For the uTrack approach we use the scheme based on fitting a mixture of Gaussian models for localization (σ = 1.5 pixels and α = 0.2). This approach uses both a Brownian and a directed motion model for position estimation. Here we set the range for the search radius to 2-4 pixels and set the multiplication factor to 7.82. The gap parameter is set to 12 frames. For the other parameters we used the default values. The MHT approach uses a waveletbased detection scheme for particle detection and we use scale 2 with a coefficient threshold of 110. The motion covariance values for both Brownian and directed motion models are set to 4, 4, and 1. For the other parameters we used the default values. Sample results for all approaches are shown in Figure 13. In this example, the particle alternates between random motion and directed motion. Some approaches yield broken trajectories. In comparison, the PDAE approach using multiple motion models obtains a trajectory without gaps that is very similar to the ground truth trajectory. For all approaches, Table II shows the quantitative results for all nine image sequences in terms of Ptrack . The ParticleTracker approach yields Ptrack = 71% (standard deviation of 11%) while the uTrack approach yields Ptrack = 71% (standard deviation of 12%). For the MHT approach, the mean tracking accuracy is Ptrack = 82% (standard deviation of 9%). The Kalman filter yields Ptrack = 84% (standard deviation of 9%) while the IPF yields a mean accuracy of Ptrack = 85% (standard deviation of 9%). The PDAE approach using only a random walk (RW) motion model yields a mean accuracy of Ptrack = 88% (standard deviation of 8%). The PDAE approach using multiple motion models yields the best results with a mean accuracy of Ptrack = 89% (standard deviation of 8%). We perform a paired-sample t-test to determine the statistical significance of the difference in performance between two approaches. It turns out that the difference in performance between the PDAE approach using multiple motion models and all other approaches is significant (p < 0.05). 2) 3D Images: We have also applied the approaches to a real microscopy 3D image sequence. The image sequence displays budding HIV-1 particles [63]. The particles are labeled with eGFP. Images were acquired with a Nikon TE2000-E spinning disk confocal microscope using an Andor Technology EM-CCD camera. The image sequence consists of 395 time steps. For each time step, a stack of 7 images (slices) with dimensions 255 × 512 pixels was acquired. The ground truth for the 3D positions of the particles was obtained manually with the Manual_Tracking plug-in of ImageJ [64]. In total, 15 14 Ground truth ParticleTracker uTrack MHT Kalman IPF-1000 PDAE (RW) PDAE Fig. 13. Ground truth and tracking results for the evaluated approaches on real 2D microscopy images (time step t = 400). The small rectangles along the trajectories indicate the intermediate positions while the large rectangle indicates the final position. Image intensities have been inverted. particles were manually tracked. For the spot detection scheme we use σLoG,xy = 1.5 voxels, σLoG,z = 1 voxel, and c = 5. Trajectories are initialized and terminated at any time step. All other parameters remain the same as for the 2D image sequences. The performance of all approaches is relatively good. In terms of the tracking accuracy Ptrack , the Kalman filter yields Ptrack =88%, the IPF yields Ptrack =90% while the PDAE approach yields Ptrack =91%. Analogous to the 2D images, the 0278-0062 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2359541, IEEE Transactions on Medical Imaging SUBMITTED TO IEEE TRANSACTIONS ON MEDICAL IMAGING 15 TABLE II R ESULTS FOR 2D REAL IMAGE SEQUENCES IN TERMS OF THE TRACKING ACCURACY PTRACK [%]. T HE MEAN VALUES AND STANDARD DEVIATIONS ARE ALSO SHOWN . Sequence (a) ParticleTracker uTrack MHT Kalman IPF-1000 PDAE (RW) PDAE 1 2 3 4 5 6 7 8 9 73 55 63 67 94 64 82 74 70 75 64 52 66 93 66 85 75 66 91 68 91 78 94 71 87 77 83 96 81 90 82 94 68 86 77 81 96 90 92 82 94 68 87 78 83 97 90 98 88 94 76 87 79 83 98 92 99 89 94 75 88 79 84 Mean Std. Dev. 71 11 71 12 82 9 84 9 85 9 88 8 89 8 (b) steps each. 14 different measures were computed to quantify the tracking performance (e.g., RMSE and Jaccard similarity index for the positions). For the images showing virus particles and receptors we used both a random walk model as well as a directed motion model. For the vesicle images we used a random walk model while for the microtubule images we used a directed motion model. In the challenge, 14 groups worldwide participated with their methods. For the different cases of the image data, the top 3 best-performing methods were identified. It turned out that our approach was quite accurate and yielded the highest number of top 3 occurrences. VIII. D ISCUSSION (c) (d) Fig. 14. Tracking results of the PDAE approach on a real 3D microscopy image sequence (time step t = 368). Individual trajectories are shown. Zslices of the original volume image are displayed. Small spheres along the trajectories represent the intermediate positions while cubes show the current position of the particles. proposed approach performs better than the other approaches. In Figure 14 sample results for individual trajectories obtained with the PDAE approach are presented. In Figure 15 all trajectories computed by the PDAE approach are shown. It can be seen that the approach deals successfully with a large number of particles. D. Evaluation on the Image Data of the ISBI’2012 Particle Tracking Challenge We have also applied the PDAE approach to the image data of the first Particle Tracking Challenge held in conjunction with the IEEE International Symposium on Biomedical Imaging (ISBI) 2012 (see [65]). The challenge comprises four different application scenarios: vesicles, virus particles, receptors, and microtubule tips. Realistic synthetic image sequences were generated for each scenario with different levels of image noise and object densities. For the virus particles, 3D image data was available and for the other scenarios 2D image data. In total, the data consists of 48 image sequences with 100 time We have introduced a new approach for tracking multiple particles in fluorescence microscopy image sequences based on probabilistic data association and an elliptical sampling scheme (PDAE). We have proposed a localization approach based on both a bottom-up strategy using the spot-enhancing filter as well as a top-down strategy using an elliptical sampling scheme for exploring the position space. The multiple measurements generated by our localization approach are integrated via the combined innovation principle of the probabilistic data association algorithm. For calculating the combined innovation, a weight is assigned to each measurement, and this weight is obtained by querying an image likelihood. In our case, the image likelihood considers both the predicted appearance of the object as well as the observed image intensities, and thus our approach uses the image information directly. For tracking multiple objects, we calculate the support that each image position provides to a tracked object relative to the support that the position provides to its neighboring objects. To incorporate multiple motion models, we used the interacting multiple model (IMM) algorithm. In comparison to tracking approaches based on the standard Kalman filter, our approach does not rely solely on a single measurement generated by a spot detection scheme and assigned by a motion correspondence algorithm. Instead, in addition to the measurement generated by a spot detection scheme, our approach steadily supplies measurements to the Kalman filter via the proposed top-down ellipsoidal sampling scheme. Thus, our approach is robust to errors arising from the spot detection scheme and in the correspondence algorithm. Moreover, our approach 0278-0062 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2359541, IEEE Transactions on Medical Imaging SUBMITTED TO IEEE TRANSACTIONS ON MEDICAL IMAGING 16 Fig. 15. Tracking results of the PDAE approach on a real 3D microscopy image sequence (time step t = 368). A z-slice (z = 6) of the original volume image is shown. Small spheres along the trajectories represent the intermediate positions while cubes show the current position of the particles. For visualization purposes, only the 50 previous intermediate positions are shown for each trajectory. Colors are used to distinguish the different trajectories. considers the appearance of individual objects and uses the image information directly, i.e., segmentation of a particle is not necessarily required for estimating the position over time. For the approach based on independent particle filters, a relatively low number of samples (e.g., Ns = 130 in 2D) is not sufficient to obtain an accurate numerical approximation of the posterior probabilities. In comparison, the PDAE approach is not subject to such numerical issues since the posterior probabilities are represented analytically via the underlying Kalman filter, and so the PDAE approach operates well with a relatively low number of samples. Also, the ellipsoidal sampling scheme explores the position space more efficiently in comparison to the random strategy used by the particle filter. In addition to the parameters of the standard Kalman filter (e.g., the covariance matrices determining the noise processes) and of the IMM algorithm (e.g., transition probabilities), the approach is parametrized by the image noise parameter σn as well as by the discretization parameters Nc , Nj , and Nk that determine the number of measurements generated by the ellipsoidal sampling scheme. The noise parameter σn is adapted based on the measured image intensities for the particles. For the discretization parameters we have used fixed values for all our experiments. We also determined the computation time of all approaches relative to the number of objects. The Kalman filter induced the lowest computation time, while the IPF entailed the highest computation time. Our approach was somewhat slower than the Kalman filter (factor of 3) but significantly faster than the particle filter (factor of about 8). The application of our approach to real microscopy image data showed its suitability for tracking HIV-1 particles in 2D as well as 3D real microscopy images. Compared to the Kalman filter and the IPF, the PDAE approach yielded the best tracking accuracy on average. Also in comparison with well-known methods (viz. the ParticleTracker approach, the uTrack approach, and MHT approach), the PDAE approach yielded a higher tracking accuracy. Certainly, the proposed approach has limitations. The approach currently uses the temporal information sequentially. Using the temporal information from the entire image sequence might improve the performance. The performance for larger object overlaps could also be improved by using a more sophisticated algorithm for motion correspondence. Another extension would be to include the image information of multiple channels. We have applied our approach to synthetic images as well to real microscopy images and carried out an experimental comparison with approaches based on the standard Kalman filter as well as based on independent particle filters (IPF). Using synthetic images, we explored different levels of image noise, and characterized the localization performance of the approaches as a function of the SNR. For high SNR levels, our approach outperformed the other approaches. For low SNR levels, our approach performed better than the Kalman filter, and yielded a similar performance as the IPF approach. In addition, by varying the number of objects in synthetic image sequences, we characterized the performance of the approaches as a function of the probability of overlap between objects. It turned out that the performance of the approaches decreased linearly with the probability of overlap. Overall, the proposed PDAE approach outperformed the other approaches. Support of the EU project SysPatho (FP7) and the BMBF project ImmunoQuant (e:Bio) is gratefully acknowledged. We thank Marko Lampe and Barbara Müller (Dept. of Virology, University of Heidelberg) for kindly providing the used 2D real microscopy image data. Sergey Ivanchenko and Don C. Lamb (Dept. of Chemistry and Biochemistry, LMU Munich) have kindly provided the 3D real microscopy image data. ACKNOWLEDGMENTS A PPENDIX Association Probabilities The probabilistic data association algorithm defines the association probabilities βi as follows [40]: ( ai , i = 1, ..., Nm βi = C (48) b i=0 C, 0278-0062 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2359541, IEEE Transactions on Medical Imaging SUBMITTED TO IEEE TRANSACTIONS ON MEDICAL IMAGING 17 where ai , b , C , 1 exp − νiT S−1 νi 2 1 1 − PD PG λF |2πS| 2 PD N m X ai . b+ (49) (50) (51) i=1 Here Nm is the number of measurements, S is the predicted measurement covariance, νi is the innovation due to the i-th measurement, λF represents the expected number of false measurements, PD denotes the detection probability of true targets, and PG is the gate probability that arises from restricting the validation region Vŷ,S (γ) with the threshold γ2. R EFERENCES [1] I. F. 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. [2] K. Jaqaman, D. Loerke, M. Mettlen, H. Kuwata, S. Grinstein, S. L. Schmid, and G. Danuser, “Robust single-particle tracking in live-cell time-lapse sequences,” Nat. Methods, vol. 5, no. 8, pp. 695–702, 2008. [3] C. Collinet, M. Stoter, C. R. Bradshaw, N. Samusik, J. C. Rink, D. Kenski, B. Habermann, F. Buchholz, R. Henschel, M. S. Mueller, W. E. Nagel, E. Fava, Y. Kalaidzidis, and M. Zerial, “Systems survey of endocytosis by multiparametric image analysis,” Nature, vol. 464, no. 7286, pp. 243–249, 2010. [4] S. Jaensch, M. Decker, A. A. Hyman, and E. W. Myers, “Automated tracking and analysis of centrosomes in early Caenorhabditis elegans embryos,” Bioinformatics, vol. 26, no. 12, pp. i13–i20, 2010. [5] A. Matov, M. M. Edvall, G. Yang, and G. Danuser, “Optimal-flow minimum-cost correspondence assignment in particle flow tracking,” Comput. Vision Image Understanding, vol. 115, no. 4, pp. 531–540, 2011. [6] F. Ruhnow, D. Zwicker, and S. Diez, “Tracking single particles and elongated filaments with nanometer precision,” Biophys. J., vol. 100, no. 11, pp. 2820–2828, 2011. [7] K. T. Applegate, S. Besson, A. Matov, M. H. Bagonis, K. Jaqaman, and G. Danuser, “plusTipTracker: Quantitative image analysis software for the measurement of microtubule dynamics,” J. Struct. Biol., vol. 176, no. 2, pp. 168–184, 2011. [8] L. Paavolainen, P. Kankaanpaa, P. Ruusuvuori, G. McNerney, M. Karjalainen, and V. Marjomaki, “Application independent greedy particle tracking method for 3D fluorescence microscopy image series,” in Proc. IEEE Int. Symp. on Biomedical Imaging: From Nano to Macro (ISBI’12), Barcelona, Spain, May 2012, pp. 672–675. [9] 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, pp. 1384–1395, 2005. [10] Q. Xue and M. C. Leake, “A novel multiple particle tracking algorithm for noisy in vivo data by minimal path optimization within the spatiotemporal volume,” in Proc. IEEE Int. Symp. on Biomedical Imaging: From Nano to Macro (ISBI’09), Boston, MA, USA, Jun.-Jul. 2009, pp. 1155–1158. [11] L. Cortes and Y. Amit, “Efficient annotation of vesicle dynamics in video microscopy,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 30, no. 11, pp. 1998–2010, 2008. [12] P. Roudot, C. Kervrann, and F. Waharte, “Lifetime estimation of moving vesicles in frequency-domain fluorescence lifetime imaging microscopy,” in Proc. IEEE Int. Symp. on Biomedical Imaging: From Nano to Macro (ISBI’12), Barcelona, Spain, May 2012, pp. 668–671. [13] A. Serge, N. Bertaux, H. Rigneault, and D. Marguet, “Dynamic multipletarget tracing to probe spatiotemporal cartography of cell membranes,” Nat. Methods, vol. 5, no. 8, pp. 687–694, 2008. [14] D. J. Rolfe, C. I. McLachlan, M. Hirsch, S. R. Needham, C. J. Tynan, S. E. D. Webb, M. L. Martin-Fernandez, and M. P. Hobson, “Automated multidimensional single molecule fluorescence microscopy feature detection and tracking,” Eur. Biophys. J., vol. 40, no. 10, pp. 1167–1186, 2011. [15] A. Genovesio, T. Liedl, V. Emiliani, W. J. Parak, M. Coppey-Moisan, and J.-C. Olivo-Marin, “Multiple particle tracking in 3-D+t microscopy: method and application to the tracking of endocytosed quantum dots,” IEEE Trans. Image. Process., vol. 15, no. 5, pp. 1062–1070, 2006. [16] A. Matov, K. Applegate, P. Kumar, C. Thoma, W. Krek, G. Danuser, and T. Wittmann, “Analysis of microtubule dynamic instability using a plus-end growth marker,” Nat. Methods, vol. 7, no. 9, pp. 761–768, 2010. [17] 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. [18] L. Yang, Z. Qiu, A. Greenaway, and W. Lu, “A new framework for particle detection in low-SNR fluorescence live-cell images and its application for improved particle tracking,” IEEE Trans. Biomed. Eng., vol. 59, no. 7, pp. 2040–2050, 2012. [19] 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. Imaging, vol. 27, no. 6, pp. 789–804, 2008. [20] J. W. Yoon, A. Bruckbauer, W. J. Fitzgerald, and D. Klenerman, “Bayesian inference for improved single molecule fluorescence tracking,” Biophys. J., vol. 94, no. 12, pp. 4932–4947, 2008. [21] W. J. Godinez, M. Lampe, S. Wörz, B. Müller, R. Eils, and K. Rohr, “Deterministic and probabilistic approaches for tracking virus particles in time-lapse fluorescence microscopy image sequences,” Med. Image Anal., vol. 13, no. 2, pp. 325–342, 2009. [22] H. Gribben, P. Miller, J. Zhang, and M. Browne, “Poisson Kalman particle filtering for tracking centrosomes in low-light 3-D confocal image sequences,” in Proc. Int. Machine Vision and Image Processing Conf. (IMVIP’09), Dublin, Ireland, Sep. 2009, pp. 83–88. [23] J. Cardinale, A. Rauch, Y. Barral, G. Székely, and I. F. Sbalzarini, “Bayesian image analysis with on-line confidence estimates and its application to microtubule tracking,” in Proc. IEEE Int. Symp. on Biomedical Imaging: From Nano to Macro (ISBI’09), Boston, MA, USA, Jun.-Jul. 2009, pp. 1091–1094. [24] 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. Imaging, vol. 31, no. 1, pp. 117–130, 2012. [25] J. Vermaak, A. Doucet, and P. Pérez, “Maintaining multi-modality through mixture tracking,” in Proc. Int. Conf. on Computer Vision (ICCV’03), Nice, France, Oct. 2003, pp. 1110–1116. [26] K. Okuma, A. Taleghani, N. de Freitas, J. Little, and D. Lowe, “A boosted particle filter: Multitarget detection and tracking,” in Proc. European Conf. on Computer Vision (ECCV’04), Prague, Czech Republic, May 2004, pp. 28–39. [27] C. Yang, R. Duraiswami, and L. S. Davis, “Fast multiple object tracking via a hierarchical particle filter,” in Proc. Int. Conf. on Computer Vision (ICCV’05), Beijing, China, Oct. 2005, pp. 212–219. [28] Y. Cai, N. de Freitas, and J. J. Little, “Robust visual tracking for multiple targets,” in Proc. European Conf. on Computer Vision (ECCV’06), Graz, Austria, May 2006, pp. 107–118. [29] Z. Khan, T. Balch, and F. Dellaert, “A Rao-Blackwellized particle filter for eigentracking,” in Proc. Computer Vision and Pattern Recognition (CVPR’04), vol. 2, Washington, DC, USA, Jun.-Jul. 2004, pp. 980–986. [30] I. Smal, E. Meijering, K. Draegestein, N. Galjart, I. Grigoriev, A. Akhmanova, M. E. van Royen, A. B. 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. [31] I. J. Cox, “A review of statistical data association for motion correspondence,” Int. J. Comput. Vision, vol. 10, no. 1, pp. 53–66, 1993. [32] I. Smal, W. Niessen, and E. Meijering, “A new detection scheme for multiple object tracking in fluorescence microscopy by joint probabilistic data association filtering,” in Proc. IEEE Int. Symp. on Biomedical Imaging: From Nano to Macro (ISBI’08), Paris, France, May 2008, pp. 264–267. [33] R. A. Kerekes, S. S. Gleason, N. Trivedi, and D. Solecki, “Automated 3-D tracking of centrosomes in sequences of confocal image stacks,” in Proc. Int. Conf. of the IEEE Eng. in Medicine and Biology Society (EMBC’09), Minneapolis, MN, USA, Sep. 2009, pp. 6994–6997. [34] S. H. Rezatofighi, S. Gould, R. Hartley, K. Mele, and W. E. Hughes, “Application of the IMM-JPDA filter to multiple target tracking in total internal reflection fluorescence microscopy images,” in Proc. Conf. on Medical Image Computing Computer-Assisted Intervention (MICCAI’12), Nice, France, Oct. 2012, pp. 357–364. 0278-0062 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2014.2359541, IEEE Transactions on Medical Imaging SUBMITTED TO IEEE TRANSACTIONS ON MEDICAL IMAGING [35] N. Chenouard, I. Bloch, and J.-C. Olivo-Marin, “Multiple hypothesis tracking for cluttered biological image sequences,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 35, no. 11, pp. 2736–2750, 2013. [36] L. Liang, H. Shen, P. D. Camilli, and J. S. Duncan, “Tracking clathrin coated pits with a multiple hypothesis based method,” in Proc. Conf. on Medical Image Computing Computer-Assisted Intervention (MICCAI’10), vol. 2, Beijing, China, Sep. 2010, pp. 315–322. [37] O. Gress and S. Posch, “Trajectory retrieval from Monte Carlo data association samples for tracking in fluorescence microscopy images,” in Proc. IEEE Int. Symp. on Biomedical Imaging: From Nano to Macro (ISBI’12), Barcelona, Spain, May 2012, pp. 374–377. [38] W. J. Godinez, M. Lampe, R. Eils, B. Müller, and K. Rohr, “Tracking multiple particles in fluorescence microscopy images via probabilistic data association,” in Proc. IEEE Int. Symp. on Biomedical Imaging: From Nano to Macro (ISBI’11), Chicago, IL, USA, Mar. 2011, pp. 1925– 1928. [39] D. Sage, F. R. Neumann, F. Hediger, S. M. Gasser, and M. Unser, “Automatic tracking of individual fluorescence particles: Application to the study of chromosome dynamics,” IEEE Trans. Image Process., vol. 14, no. 9, pp. 1372–1382, 2005. [40] T. Kirubarajan and Y. Bar-Shalom, “Probabilistic data association techniques for target tracking in clutter,” Proc. IEEE, vol. 92, no. 3, pp. 536–557, 2004. [41] C. Rasmussen and G. D. Hager, “Probabilistic data association methods for tracking complex visual objects.” IEEE Trans. Pattern. Anal. Mach. Intell., vol. 23, no. 6, pp. 560–576, 2001. [42] H. A. P. Blom, “An efficient filter for abruptly changing systems,” in Proc. IEEE Conf. on Decision and Control, Las Vegas, NV, USA, Dec. 1984, pp. 656–658. [43] M. Sargin, A. Altinok, B. Manjunath, and K. Rose, “Variable length open contour tracking using a deformable trellis,” IEEE Trans. Image. Process., vol. 20, no. 4, pp. 1023–1035, 2011. [44] R. Kidambi, M.-C. Shih, and K. Rose, “Deformable trellises on factor graphs for robust microtubule tracking in clutter,” in Proc. IEEE Int. Symp. on Biomedical Imaging: From Nano to Macro (ISBI’12), Barcelona, Spain, May 2012, pp. 676–679. [45] J. Cui, S. T. Acton, and Z. Lin, “A Monte Carlo approach to rolling leukocyte tracking in vivo,” Med. Image Anal., vol. 10, no. 4, pp. 598– 610, 2006. [46] C. Restif, C. Ibañez-Ventoso, M. Driscoll, and D. N. Metaxas, “Tracking C. elegans swimming for high-throughput phenotyping,” in Proc. IEEE Int. Symp. on Biomedical Imaging: From Nano to Macro (ISBI’11), Chicago, IL, USA, Mar. 2011, pp. 1542–1548. [47] R. E. Kalman, “A new approach to linear filtering and prediction problems,” J. of Basic Eng., vol. 82, pp. 35–45, 1960. [48] B. Anderson and J. Moore, Optimal Filtering. Mineola, NY, USA: Dover Publications, 1979. [49] 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. [50] 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. [51] P. Ruusuvuori, T. Äijö, S. Chowdhury, C. Garmendia-Torres, J. Selinummi, M. Birbaumer, A. M. Dudley, L. Pelkmans, and O. YliHarja, “Evaluation of methods for detection of fluorescence labeled subcellular objects in microscope images,” BMC Bioinf., vol. 11, p. 248, 2010. [52] L. Liang, H. Shen, P. D. Camilli, D. Toomre, and J. S. Duncan, “An expectation maximization based method for subcellular particle tracking using multi-angle TIRF microscopy,” in Proc. Conf. on Medical Image Computing Computer-Assisted Intervention (MICCAI’11), vol. 2, Toronto, Canada, Sep. 2011, pp. 629–636. [53] A. Genovesio and J.-C. Olivo-Marin, “Split and merge data association filter for dense multi-target tracking,” in Proc. Int. Conf. on Pattern Recognition (ICPR’04), vol. 4, Cambridge, UK, Aug. 2004, pp. 677– 680. [54] T. Yu and Y. Wu, “Collaborative tracking of multiple targets,” in Proc. Computer Vision and Pattern Recognition (CVPR’04), vol. 1, Washington D. C., USA, Jun. 2004, pp. 834–841. [55] H. A. P. Blom and Y. Bar-Shalom, “The interacting multiple model algorithm for systems with Markovian switching coefficients,” IEEE Trans. Autom. Control, vol. 33, no. 8, pp. 780–783, 1988. [56] B. D. Lucas and T. Kanade, “An iterative image registration technique with an application to stereo vision,” in Proc. Int. Joint Conf. on Artificial 18 [57] [58] [59] [60] [61] [62] [63] [64] [65] Intelligence (IJCAI’81), Vancouver, BC, Canada, Aug. 1981, pp. 674– 679. M. K. Cheezum, W. F. Walker, and W. H. Guilford, “Quantitative comparison of algorithms for tracking single fluorescent particles,” Biophys. J., vol. 81, pp. 2378–2388, 2001. S. Wörz, P. Sander, M. Pfannmöller, R. Rieker, S. Joos, G. Mechtersheimer, P. Boukamp, P. Lichter, and K. Rohr, “3D geometry-based quantification of colocalizations in multichannel 3D microscopy images of human soft tissue tumors,” IEEE Trans. Med. Imaging, vol. 29, no. 8, pp. 1474 –1484, 2010. L. Liang, “A multiple hypothesis based particle tracking method with application to clathrin mediated endocytosis analysis,” Ph.D. dissertation, Yale University, 2013. L. Liang, H. Shen, Y. Xu, P. D. Camilli, D. Toomre, and J. S. Duncan, “A Bayesian method for 3D estimation of subcellular particle features in multi-angle TIRF microscopy,” in Proc. IEEE Int. Symp. on Biomedical Imaging: From Nano to Macro (ISBI’12), Barcelona, Spain, May 2012, pp. 984–987. M. Lampe, J. A. G. Briggs, T. Endress, B. Glass, S. Riegelsberger, H.-G. Kräusslich, D. C. Lamb, C. Bräuchle, and B. Müller, “Double-labelled HIV-1 particles for study of virus-cell interaction,” Virology, vol. 360, no. 1, pp. 92–104, 2007. F. de Chaumont, S. Dallongeville, N. Chenouard, N. Hervé, S. Pop, T. Provoost, V. Meas-Yedid, P. Pankajakshan, T. Lecomte, Y. Le Montagner, T. Lagache, A. Dufour, and J.-C. C. Olivo-Marin, “Icy: an open bioimage informatics platform for extended reproducible research.” Nat. Methods, vol. 9, no. 7, pp. 690–696, 2012. S. Ivanchenko, W. J. Godinez, M. Lampe, H.-G. Kräusslich, R. Eils, K. Rohr, C. Bräuchle, B. Müller, and D. C. Lamb, “Dynamics of HIV-1 assembly and release,” PLoS Pathog, vol. 5, no. 11, p. e1000652, 2009. M. Abramoff, P. Magelhaes, and S. Ram, “Image Processing with ImageJ,” Biophotonics International, vol. 11, no. 7, pp. 36–42, 2004. N. Chenouard, I. Smal, F. de Chaumont, M. Maska, I. F. Sbalzarini, Y. Gong, J. Cardinale, C. Carthel, S. Coraluppi, M. Winter, A. R. Cohen, W. J. Godinez, K. Rohr, Y. Kalaidzidis, L. Liang, J. Duncan, H. Shen, Y. Xu, K. E. G. Magnusson, J. Jalden, H. M. Blau, P. PaulGilloteaux, P. Roudot, C. Kervrann, F. Waharte, J.-Y. Tinevez, S. L. Shorte, J. Willemse, K. Celler, G. P. van Wezel, H.-W. Dan, Y.-S. Tsai, C. Ortiz de Solorzano, J.-C. Olivo-Marin, and E. Meijering, “Objective comparison of particle tracking methods,” Nat. Methods, vol. 11, no. 3, pp. 281–289, 2014. 0278-0062 (c) 2013 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

RELATED PAPERS

RELATED TOPICS

Log In

Tracking Multiple Particles in Fluorescence Time-Lapse Microscopy Images via Probabilistic Data Association

Tracking Multiple Particles in Fluorescence Time-Lapse Microscopy Images via Probabilistic Data Association

Related Papers

RELATED PAPERS

RELATED TOPICS