Tracking multiple moving objects in images using Markov Chain Monte Carlo

A new Bayesian state and parameter learning algorithm for multiple target tracking models with image observations are proposed. Specifically, a Markov chain Monte Carlo algorithm is designed to sample from the posterior distribution of the unknown time-varying number of targets, their birth, death times and states as well as the model parameters, which constitutes the complete solution to the specific tracking problem we consider. The conventional approach is to pre-process the images to extract point observations and then perform tracking, i.e. infer the target trajectories. We model the image generation process directly to avoid any potential loss of information when extracting point observations using a pre-processing step that is decoupled from the inference algorithm. Numerical examples show that our algorithm has improved tracking performance over commonly used techniques, for both synthetic examples and real florescent microscopy data, especially in the case of dim targets with overlapping illuminated regions.

1. Fab (Fragment antigen binding) is a region on an antibody that binds to antigens.

2. When \(\sigma _h\) is too small, the illuminated region taken into account is smaller than it should be which would cause many more targets than expected. In that case, we should increase \(\sigma _h\).

We thank Kristina Ganzinger and Professor David Klenerman for providing the real data and the code in Weimann et al. (2013) for the comparisons in Sect.  4.2, and Sinan Yıldırım for his careful reading of this paper.

Authors and Affiliations

1. Department of Engineering, University of Cambridge, Cambridge, UK

   Lan Jiang & Sumeetpal S. Singh

Corresponding author

Correspondence to Sumeetpal S. Singh.

This work was supported by the Engineering and Physical Sciences Research Council [Grant Numbers EP/G037590/1, EP/K020153/1].

Appendix

1.1 Appendix 1: Hypothesis testing

The term of (26) to be calculated is \(p( y^r_{t,L(i)} | H_1)\) where \(y^r_{t,L(i)}=\{y^r_{t,j}\}_{j \in L(i)} \) and \(i \in G_t\) is a pixel that is a local maximum of \(y_t^f\). Let (r, c) be its row and column number and

Vector \((\bar{a},\bar{s})\) can be interpreted as the likely intensity \(\bar{a}\) and location \(\bar{s}\) of an undetected target. Below we just write L as the set of pixels under consideration instead of L(i).

Let \(x=(a,s,v)\in \mathbb {R} \times \mathbb {R}^2 \times \mathbb {R}^2\) where as before a denotes intensity, \(s=(s(1),s(2))\) spatial coordinates and \(v=(v(1),v(2))\) spatial velocity. (Recall that a pixel illumination is not a function of velocity.) The aim is to calculate

where \(p(a,s|H_1)\) is either the marginal (or restriction to intensity and spatial position only) of the law of the birth \(\mu _{\psi }\) [see (1)] if proposing the initial state of the birth move, or the pdf \(p(a,s|a_{0:k}, a_{0:k})\) (to be defined below) if extending the target intensity and position trajectory after having created the initial intensity and position. \(p(y^r_{t,L}|a,s)\) is implicitly defined.

Let \(p(y^r_{t,L},a,s|H_1)=p(y^r_{t,L}|a,s) p(a,s|H_1)\). Use the approximation

where \(-D\) is the second-order derivative \(\nabla ^2 \ln p(y^r_{t,L},a,s|H_1)\) evaluated at \((\bar{a},\bar{s})\). Expression (41) is like the Laplace approximation except that the second-order Taylor expansion is computed at \((y^f_{t,i},{\varDelta }r, {\varDelta }c)\) and not the true maximum \(\arg \max _{a,s} p(y^r_{t,L},a,s|H_1)\) to save on the maximisation step, which we find in the numerical examples to be still effective as a component of the birth move. (Moreover, it is a fair simplification for a diffused prior.) Thus,

and the Gaussian distribution in (29) is

Section 4 (numerical examples) assumes a linear and Gaussian model for the targets in both the synthetic and real data examples. The description is now completed by specifying \(p(\bar{a},\bar{s}|H_1)\).

When the birth move is constructing the initial/first state of the target,

and the expression in (26) now simplifies to

Finally, we derive \(p(a,s| H_1)\) in (40) when the birth move is extending the target intensity and position trajectory after having created the initial intensity position pairs \((a_{0}, s_{0}),\ldots ,\) \((a_{k-1}, s_{k-1})\) for some \(k\ge 1\). For a target with a linear Gaussian model (1), the pdf of \(v_{k-1}\) (or \(v_{0:k-1}\)) conditioned on \(s_{0:k-1}\), which is denoted \(p_{\psi }(v_{k-1}|{s}_{0:k-1})\), is a Gaussian. Thus, \(p(a,s| H_1)\) is

and the corresponding expression for \(\rho (y_{t_k,L(i)}^r)\) in (30) is the same as in (44).

1.2 Appendix 2: State proposal

The state proposal of Sect. 3.7 alters the state values of the targets whose trajectories have been partially exchanged. This proposal is defined for the Gaussian model in Sect. 4.

Let \(\left\{ (k,t_{b}^{k},\mathbf {x}^{k})\right\} _{k=1}^{K}\) be the MTT state and assume without loss of generality \(U=\{1,2\}\). The state proposal is decomposed into

where the first term is the probability of selecting (U, t) and is not \(y_{1:n}\) dependent. Using \(y_{1:n}\) and \(\left\{ (k,t_{b}^{k},\mathbf {x}^{k})\right\} _{k=3}^{K}\), generate the residual image as in (23) by subtracting the background intensity and the contribution from all targets except \((t_{b}^{1},\mathbf {x}^{1})\) and \((t_{b}^{2},\mathbf {x}^{2}\)). Let \(y_{t}^{r}=(y_{t,1}^{r},\ldots ,y_{t,m}^{r})\) denote the residual images at time t. The state proposal samples \((\tilde{\mathbf {x}}^{1},\tilde{\mathbf {x}}^{2})\) from the pdf \(q_{s,\theta }(\tilde{\mathbf {x}}^{1},\tilde{\mathbf {x}}^{2}|\hat{\mathbf {x}}^{1:2},\tau _{b}^{1:2},y_{1:n}^{r})\) where \((\hat{\mathbf {x}}^{1},\hat{\mathbf {x}}^{2})\) is defined in (33). Note the dependency on targets \(k>2\) is captured through the residual image.

For brevity, instead of \((\hat{\mathbf {x}}^{1},\hat{\mathbf {x}}^{2})\), we write \(({\mathbf {x}}^{1},{\mathbf {x}}^{2})\). Also, to highlight the intensity, spatial coordinates and velocity components, \(\mathbf {x}^{i}=(a_{0}^{i},s_{0}^{i},v_{0}^{i},\ldots ,a_{l^{i}-1}^{i},s_{l^{i}-1}^{i},v_{l^{i}-1}^{i})\) is written as \((\mathbf {a}^{i},\mathbf {s}^{i},\mathbf {v}^{i})\).

The proposal \(q_{s,\theta }\) does not alter the spatial components, i.e. \(\tilde{\mathbf {s}}^{1}=\mathbf {{s}}^{1}\) and \(\tilde{\mathbf {s}}^{2}=\mathbf {{s}}^{2}\).

The velocities \((\tilde{\mathbf {v}}^{1},\tilde{\mathbf {v}}^{2})\sim p_{\psi }(\tilde{\mathbf {v}}^{1}|\mathbf {{s}}^{1})p_{\psi }(\tilde{\mathbf {v}}^{2}|\mathbf {{s}}^{2})\), i.e. are sampled independently from \(p_{\psi }\), which is the prior pdf of the velocity conditioned on the spatial coordinates values, which is Gaussian.

If targets 1 and 2 exist at time t, then their state values are \(x_{t-\tau _{b}^{1}}^{1}\) and \(x_{t-\tau _{b}^{2}}^{2}\), respectively. We assume (for pixel i) \(y_{t,i}^{r}\sim \mathcal {N}(\cdot |0,\sigma _{r,t}^{2})\) if targets 1 and 2 both do not exist at time t, \(y_{t,i}^{r}\sim \mathcal {N}(\cdot |a_{t-\tau _{b}^{1}}^{1}h_{i}(s_{t-\tau _{b}^{1}}^{1}),\sigma _{r,t}^{2})\) if only target 1 exists and \(y_{t,i}^{r}\sim \mathcal {N}(\cdot |a_{t-\tau _{b}^{1}}^{1}h_{i}(s_{t-\tau _{b}^{1}}^{1})+a_{t-\tau _{b}^{2}}^{2}h_{i}(s_{t-\tau _{b}^{2}}^{2}),\sigma _{r,t}^{2})\) if both exists. The prior probability model [see (1))] for the intensities is independent Gaussians, denoted \(p_{\psi }(\mathbf {a}^{1})p_{\psi }(\mathbf {a}^{2})\). Conditioned on \(y_{1:n}^{r}\), \((\tau _{b}^{1},\mathbf {s}^{1})\) and \((\tau _{b}^{2},\mathbf {s}^{2})\), the posterior pdf for the joint intensities is also a Gaussian, which is denoted by \(q_{s,\theta }(\tilde{\mathbf {a}}^{1},\tilde{\mathbf {a}}^{2} | \mathbf {s}^{1:2},\tau _{b}^{1:2},y_{1:n}^{r})\). Thus

1.3 Appendix 3: CSMC step of Algorithm 1

(Initialisation.) Denote target k’s trajectory by \((t_{b}^{k},\mathbf {x}^{k})\) with \(\mathbf {x}^{k}=(x_{0},\ldots ,x_{l-1})\). Let \(\mathbf {x}'_{1:n}\) denote the MTT state, in the representation of Sect. 2.3, omitting target k. Define target k’s state likelihood to be

[c.f. (11).]

Execute the CSMC algorithm of Whiteley (2010), with backward sampling, for the partially observed Markov model with state law \(\mu _{\psi }(x_{0})\), \(f_{\psi }(x_{t}\vert x_{t-1})\) [see (1)], state likelihood \(g_{t_{b}^{k}+t}(x_{t})\) above and input trajectory \((x_{0},\ldots ,x_{l-1})\) to yield an updated target k trajectory \((x'_{0},\ldots ,x'_{l-1})\). The CSMC implementation used proposes particles from the state law.

1.4 Appendix 4: Gibbs step of Algorithm 1 for the MTT model of Sect. 4

Recall \(K_{t}^{s}\) is a binomial r.v. with success probability \(p_{s}\) and number of trials \(K_{t-1}^{x}\). \(K_t^b\) is a Poisson r.v. with rate \(\lambda _b\). Thus, their posteriors are the Beta and Gamma pdfs:

The updates for \(b_t\) and \(\sigma _{r,t}\) require the following empirical means and variances [c.f. (10], (11)): \(\bar{e}_{t} = m^{-1}\sum _{i=1}^m y_{t,i} - h_i (\mathbf {x}_t)\), \(\beta _1 = \sum _{i=1}^m( y_{t,i} - h_i (\mathbf {x}_t) - \bar{e}_t)^2\), \(\beta _2= \frac{n_0 m}{n_0 + m} (\mu _0 - \bar{e}_{t})^2\),

The Gibbs update for means and variances of the remaining Gaussians of the MTT model of Sect. 4 follows similar update formulae.

Reprints and permissions