Condensation | conditional density propagation for visual tracking Michael Isard and Andrew Blake, Department of Engineering Science, University of Oxford, Oxford OX1 3PJ, UK. Int. J. Computer Vision, in press (1998). ABSTRACT The problem of tracking curves in dense visual clutter is challenging. Kalman ltering is inadequate because it is based on Gaussian densities which, being unimodal, cannot represent simultaneous alternative hypotheses. The Condensation algorithm uses \factored sampling", previously applied to the interpretation of static images, in which the probability distribution of possible interpretations is represented by a randomly generated set. Condensation uses learned dynamical models, together with visual observations, to propagate the random set over time. The result is highly robust tracking of agile motion. Notwithstanding the use of stochastic methods, the algorithm runs in near real-time. Contents 1 2 3 4 5 6 7 8 A B C Tracking curves in clutter Discrete-time propagation of state density Factored sampling The Condensation algorithm Stochastic dynamical models for curve motion Observation model Applying the Condensation algorithm to video-streams Conclusions Non-linear ltering Derivation of the sampling rule Asymptotic correctness of the Condensation Algorithm 2 3 6 8 10 13 17 26 31 33 34 Isard and Blake. IJCV, in press, (1998). 2 1 Tracking curves in clutter The purpose of this paper1 is to establish a stochastic framework for tracking curves in visual clutter, using a sampling algorithm. The approach is rooted in ideas from statistics, control theory and computer vision. The problem is to track outlines and features of foreground objects, modelled as curves, as they move in substantial clutter, and to do it at, or close to, video frame-rate. This is challenging because elements in the background clutter may mimic parts of foreground features. In the most severe case of camouage, the background may consist of objects similar to the foreground object, for instance when a person is moving past a crowd. Our approach aims to dissolve the resulting ambiguity by applying probabilistic models of object shape and motion to analyse the video-stream. The degree of generality of these models is pitched carefully: suciently specic for eective disambiguation but suciently general to be broadly applicable over entire classes of foreground objects. 1.1 Modelling shape and motion Eective methods have arisen in computer vision for modelling shape and motion. When suitable geometric models of a moving object are available, they can be matched eectively to image data, though usually at considerable computational cost (Hogg, 1983; Lowe, 1991; Sullivan, 1992; Huttenlocher et al., 1993). Once an object has been located approximately, tracking it in subsequent images becomes more ecient computationally (Lowe, 1992), especially if motion is modelled as well as shape (Gennery, 1992; Harris, 1992). One important facility is the modelling of curve segments which interact with images (Fischler and Elschlager, 1973; Yuille and Hallinan, 1992) or image sequences (Kass et al., 1987; Dickmanns and Graefe, 1988). This is more general than modelling entire objects but more clutter-resistant than applying signalprocessing to low-level corners or edges. The methods to be discussed here have been applied at this level, to segments of parametric B-spline curves (Bartels et al., 1987) tracking over image sequences (Menet et al., 1990; Cipolla and Blake, 1990). The B-spline curves could, in theory, be parameterised by their control points. In practice this allows too many degrees of freedom for stable tracking and it is necessary to restrict the curve to a low-dimensional parameter x, for example over an ane space (Koenderink and Van Doorn, 1991; Ullman and Basri, 1991; Blake et al., 1993), or more generally allowing a \shape-space" of non-rigid motion (Cootes et al., 1993). Finally, prior probability densities can be dened over the curves (Cootes et al., 1993) represented by appropriate parameter vectors x, and also over their motions (Terzopoulos and Metaxas, 1991; Blake et al., 1993), and this constitutes a powerful facility for tracking. Reasonable defaults can be chosen for those densities. However, it is obviously more satisfactory to measure or estimate them from data-sequences (x1 ; x2 ; : : :). Algorithms to do this, assuming Gaussian densities, are known in the control-theory literature (Goodwin and Sin, 1984) and have been applied in computer vision (Blake and Isard, 1994; Baumberg and Hogg, 1995). Given the prior, and an observation density that characterises the statistical variability of image data z given a curve state x, a posterior distribution can, in principle, be estimated for xt given zt at successive times t. This paper has appeared in short form (Isard and Blake, 1996) as joint winner of the prize of the European Conference on Computer Vision, 1996. 1 Isard and Blake. IJCV, in press, (1998). 3 1.2 Kalman lters and data-association Spatio-temporal estimation, the tracking of shape and position over time, has been dealt with thoroughly by Kalman ltering, in the relatively clutter-free case in which p(xt ) can satisfactorily be modelled as Gaussian (Dickmanns and Graefe, 1988; Harris, 1992; Gennery, 1992; Rehg and Kanade, 1994; Matthies et al., 1989) and can be applied to curves (Terzopoulos and Szeliski, 1992; Blake et al., 1993). These solutions work relatively poorly in clutter which causes the density for xt to be multi-modal and therefore non-Gaussian. With simple, discrete features such as points or corners combinatorial data-association methods can be eective with clutter but combinatorial methods to do not apply naturally to curves. There remains a need for an appropriately general probabilistic mechanism to handle multi-modal density functions. 1.3 Temporal propagation of conditional densities The Kalman lter as a recursive linear estimator is a special case, applying only to Gaussian densities, of a more general probability density propagation process. In continuous time this can be described in terms of diusion, governed by a \Fokker-Planck" equation (Astrom, 1970), in which the density for xt drifts and spreads under the action of a stochastic model of its dynamics. In the simple Gaussian case, the diusion is purely linear and the density function evolves as a Gaussian pulse that translates, spreads and is reinforced, remaining Gaussian throughout, as in gure 1, a process that is described analytically and exactly by the Kalman lter. The random component of the dynamical model leads to spreading | increasing uncertainty | while the deterministic component causes the density function to drift bodily. The eect of an external observation zt is to superimpose a reactive eect on the diusion in which the density tends to peak in the vicinity of observations. In clutter, there are typically several competing observations and these tend to encourage a non-Gaussian state-density (gure 2). The Condensation algorithm is designed to address this more general situation. It has the striking property that, generality notwithstanding, it is a considerably simpler algorithm than the Kalman lter. Moreover, despite its use of random sampling which is often thought to be computationally inecient, the Condensation algorithm runs in near real-time. This is because tracking over time maintains relatively tight distributions for shape at successive time-steps, and particularly so given the availability of accurate, learned models of shape and motion. 2 Discrete-time propagation of state density For computational purposes, the propagation process must be set out in terms of discrete time t. The state of the modelled object at time t is denoted xt and its history is Xt = fx1 ; : : : ; xt g. Similarly the set of image features at time t is zt with history Zt = fz1 ; : : : ; zt g. Note that no functional assumptions (linearity, Gaussianity, unimodality) are made about densities in the general treatment, though particular choices will be made in due course in order to demonstrate the approach. Isard and Blake. IJCV, in press, (1998). 4 deterministic drift p(x) p(x) x x stochastic diffusion p(x) z p(x) x x reactive effect of measurement Figure 1: Kalman lter as density propagation. In the case of Gaussian prior, process and observation densities, and assuming linear dynamics, the propagation process of gure 2 reduces to a diusing Gaussian state density, represented completely by its evolving (multivariate) mean and variance | precisely what a Kalman lter computes. 2.1 Stochastic dynamics A somewhat general assumption is made for the probabilistic framework that the object dynamics form a temporal Markov chain so that p(xt jXt?1 ) = p(xt jxt?1 ) (1) | the new state is conditioned directly only on the immediately preceding state, independent of the earlier history. This still allows quite general dynamics, including stochastic dierence equations of arbitrary order; we use second order models and details are given later. The dynamics are entirely determined therefore by the form of the conditional density p(xt jxt?1 ). For instance, p(xt jxt?1 ) / exp ? 21 (xt ? xt?1 ? 1)2 represents a one-dimensional random walk (discrete diusion) whose step length is a standard normal variate, superimposed on a rightward drift at unit speed. Of course, for realistic problems, the state x is multi-dimensional and the density is more complex (and, in the applications presented later, learned from training sequences). Isard and Blake. IJCV, in press, (1998). 5 deterministic drift p(x) p(x) x x stochastic diffusion p(x) z p(x) x x reactive effect of measurement Figure 2: Probability density propagation. Propagation is depicted here as it occurs over a discrete time-step. There are three phases: drift due to the deterministic component of object dynamics; diusion due to the random component; reactive reinforcement due to observations. 2.2 Measurement Observations zt are assumed to be independent, both mutually and with respect to the dynamical process. This is expressed probabilistically as follows: p(Zt?1 ; xt jXt?1 ) = p(xt jXt?1 ) tY ?1 i=1 p(zi jxi): (2) Note that integrating over xt implies the mutual conditional independence of observations: p(Zt jXt ) = t Y i=1 p(zi jxi ): (3) The observation process is therefore dened by specifying the conditional density p(zt jxt ) at each time t, and later, in computational examples, we take this to be a time-independent function p(zjx). Suce it to say for now that, in clutter, the observation density is multi-modal. Details will be given in section 6 Isard and Blake. IJCV, in press, (1998). 6 2.3 Propagation Given a continuous-valued Markov chain with independent observations, the conditional statedensity pt at time t is dened by pt(xt ) p(xt jZt ): This represents all information about the state at time t that is deducible from the entire datastream up to that time. The rule for propagation of state density over time is: p(xt jZt ) = kt p(zt jxt )p(xt jZt?1 ); where p(xt jZt?1 ) = Z xt?1 p(xt jxt?1 )p(xt?1 jZt?1 ) (4) (5) and kt is a normalisation constant that does not depend on xt . The validity of the rule is proved in the appendix. The propagation rule (4) should be interpreted simply as the equivalent of the Bayes' rule (6) for inferring posterior state density from data, for the time-varying case. The eective prior p(xt jZt?1 ) is actually a prediction taken from the posterior p(xt?1 jZt?1 ) from the previous time-step, onto which is superimposed one time-step from the dynamical model (Fokker-Planck drift plus diusion as in gure 2), which is expressed in (5). Multiplication in (4) by the observation density p(zt jxt ) in the Bayesian manner then applies the reactive eect expected from observations. Because the observation density is non-Gaussian, the evolving state density p(xt jZt ) is also generally non-Gaussian. The problem now is how to apply a nonlinear lter to evaluate the state density over time, without incurring excessive computational load. Inevitably this means approximating. Numerous approaches, including \multiple hypothesis tracking", have been proposed but prove unsuitable for use with curves as opposed to discrete features | details are given in the appendix. In this paper we propose a sampling approach which is described in the following two sections. 3 Factored sampling This section describes rst the factored sampling algorithm dealing with non-Gaussian observations in single images. Then factored sampling is extended in the following section to deal with temporal image sequences. A standard problem in statistical pattern recognition is to nd an object parameterised as x with prior p(x), using data z from a single image. The posterior density p(xjz) represents all the knowledge about x that is deducible from the data. It can be evaluated in principle by applying Bayes' rule (Papoulis, 1990) to obtain p(xjz) = kp(zjx)p(x) (6) where k is a normalisation constant that is independent of x. In cases where p(zjx) is suciently complex that p(xjz) cannot be evaluated simply in closed form, iterative sampling techniques can be used (Geman and Geman, 1984; Ripley and Sutherland, 1990; Grenander et al., 1991; Storvik, 1994). The factored sampling algorithm (Grenander et al., 1991) generates a random variate x from a distribution p~(x) that approximates the posterior p(xjz). First a sample-set Isard and Blake. IJCV, in press, (1998). 7 fs ; : : : ; s N g is generated from the prior density p(x) and then an index i 2 f1; : : : ; N g is (1) ( ) chosen with probability i , where and (i) i = PNpz (s )(j) j =1 pz (s ) pz (x) = p(zjx); the conditional observation density. The value x0 = xi chosen in this fashion has a distribution which approximates the posterior p(xjz) increasingly accurately as N increases (gure 3). Probability posterior density weighted sample State x Figure 3: Factored sampling. A set of points s(i) , the centres of the blobs in the gure, is sampled randomly from a prior density p(x). Each sample is assigned a weight i (depicted by blob area) in proportion to the value of the observation density p(zjx = s(i) ). The weighted point-set then serves as a representation of the posterior density p(xjz), suitable for sampling. The one-dimensional case illustrated here extends naturally to the practical case that the density is dened over several position and shape variables. Note that posterior mean properties E [g(x)jz] can be generated directly from the samples fs n g by weighting with pz (x) to give: ( ) E [g(x)jz] PN (n) )pz (s(n) ) : nP =1 g (s N p (s(n) ) n=1 z (7) For example, the mean can be estimated using g(x) = x (illustrated in gure 4) and the variance using g(x) = xxT . In the case that p(x) is a spatial Gauss-Markov process, Gibbs sampling from p(x) has been used to generate the random variates fs(1) ; : : : ; s(N ) g. Otherwise, for low-dimensional parameterisations as in this paper, standard, direct methods can be used for Gaussians2 (Press et al., 1988). Note that, in the case that the density p(zjx) is normal, the mean obtained by factored sampling is consistent with an estimate obtained more conventionally, and eciently, from linear least squares estimation. For multi-modal distributions which Note: the presence of clutter causes p(z x) to be non-Gaussian, but the prior p(x) may still happily be Gaussian, and that is what will be assumed in our experiments. 2 j Isard and Blake. IJCV, in press, (1998). 8 cannot be approximated as normal, so that linear estimators are unusable, estimates of mean x by factored sampling continue to apply. a) b) Figure 4: Sample-set representation of shape distributions The sample-set representation of probability distributions, illustrated in one dimension in gure 3, is illustrated here (a) as it applies to the distribution of a multi-dimensional curve parameter x. Each sample s(n) is shown as a curve (of varying position and shape) with a thickness proportional to the weight n. The weighted mean of the sample set (b) serves as an estimator of the distribution mean. 4 The Condensation algorithm The Condensation algorithm is based on factored sampling but extended to apply iteratively to successive images in a sequence. The same sampling strategy has been developed elsewhere (Gordon et al., 1993; Kitagawa, 1996), presented as developments of Monte-Carlo methods. Jump-diusion tracking (Miller et al., 1995) may also be related to the approach described here. Given that the process at each time-step is a self-contained iteration of factored sampling, the output of an iteration will be a weighted, time-stamped sample-set, denoted fs(tn) ; n = 1; : : : ; N g with weights t(n) , representing approximately the conditional state-density p(xt jZt ) at time t. How is this sample-set obtained? Clearly the process must begin with a prior density and the eective prior for time-step t should be p(xt jZt?1 ). This prior is of course multi-modal in general and no functional representation of it is available. It is derived from the sample set representation f(s(t?n)1 ; t(?n)1 ); n = 1; : : : ; N g of p(xt?1 jZt?1 ), the output from the previous time-step, to which prediction (5) must then be applied. The iterative process as applied to sample-sets, depicted in gure 5, mirrors the continuous diusion process in gure 2. At the top of the diagram, the output from time-step t ? 1 is the weighted sample-set f(s(t?n)1 ; t(?n)1 ); n = 1; : : : ; N g. The aim is to maintain, at successive time-steps, sample sets of xed size N , so that the algorithm can be guaranteed to run within a given computational resource. The rst operation therefore is to sample (with replacement) Isard and Blake. IJCV, in press, (1998). 9 (n) (n) sk−1 ,πk−1 drift diffuse observation density (n) sk measure (n) (n) sk , πk Figure 5: One time-step in the Condensation algorithm. Each of the three steps | driftdiuse-measure | of the probabilistic propagation process of gure 2 is represented by steps in the Condensation algorithm. Isard and Blake. IJCV, in press, (1998). 10 N times from the set fs(t?n)1 g, choosing a given element with probability t(?n)1 . Some elements, especially those with high weights, may be chosen several times, leading to identical copies of elements in the new set. Others with relatively low weights may not be chosen at all. Each element chosen from the new set is now subjected to the predictive steps. First, an element undergoes drift and, since this is deterministic, identical elements in the new set undergo the same drift. This is apparent in the gure. The second predictive step, diusion, is random and identical elements now split because each undergoes its own independent Brownian motion step. At this stage, the sample set fs(tn) g for the new time-step has been generated but, as yet, without its weights; it is approximately a fair random sample from the eective prior density p(xt jZt?1 ) for time-step t. Finally, the observation step from factored sampling is applied, generating weights from the observation density p(zt jxt ) to obtain the sample-set representation f(s(tn) ; t(n) )g of state-density for time t. Figure 6 gives a synopsis of the algorithm. Note the use of cumulative weights c(t?j )1 (constructed in step 3) to achieve ecient sampling in step 1. After any time-step, it is possible to \report" on the current state, for example by evaluating some moment of the state density as shown. One of the striking properties of the Condensation algorithm is its simplicity, compared with the Kalman lter, despite its generality. Largely this is due to the absence of the Riccati equation which appears in the Kalman lter for the propagation of covariance. The Riccati equation is relatively complex computationally but is not required in the Condensation algorithm which instead deals with variability by sampling, involving the repeated computation of a relatively simple propagation formula. 5 Stochastic dynamical models for curve motion In order to apply the Condensation algorithm, which is general, to tracking curves in imagestreams, specic probability densities must be established both for the dynamics of the object and for the observation process. In the examples described here, x is a linear parameterisation of the curve and allowed transformations of the curve are represented by linear transformations of x. The Condensation algorithm itself does not demand necessarily a linear parameterisation though linearity is an attraction for another reason | the availability of algorithms to learn object dynamics. The algorithm could also be used, in principle, with non-linear parameterised kinematics | for instance representing an articulated hand in terms of joint angles (Rehg and Kanade, 1994). 5.1 Linear parameterisations of splines for tracking We represent the state of a tracked object following methods established for tracking using a Kalman lter (Blake et al., 1995). Objects are modelled as a curve (or set of curves), typically though not necessarily the occluding contour, and represented at time t by a parameterised image curve r(s; t). The parameterisation is in terms of B-splines, so r(s; t) = (B(s) Qx(t); B(s) Qy (t)) for 0 s L (8) where B(s) is a vector (B (s); : : : ; BNB (s))T of B-spline basis functions, Qx and Qy are vectors 1 of B-spline control point coordinates and L is the number of spans. It usually desirable (Blake Isard and Blake. IJCV, in press, (1998). 11 Iterate From the \old" sample-set fs(t?n)1 ; t(?n)1 ; c(t?n)1 ; n = 1; : : : ; N g at time-step t ? 1, construct a \new" sample-set fs(tn) ; t(n) ; c(tn) g; n = 1; : : : ; N for time t. Construct the nth of N new samples as follows: 1. Select a sample s0t (n) as follows: (a) generate a random number r 2 [0; 1], uniformly distributed. (b) nd, by binary subdivision, the smallest j for which c(t?j )1 r (c) set s0t (n) = s(t?j )1 2. Predict by sampling from p(xt jxt?1 = s0 (tn) ) to choose each s(tn) . For instance, in the case that the dynamics are governed by a linear stochastic dierential equation, the new sample value may be generated as: st(n) = As0 (tn) + B wt(n) where wt(n) is a vector of standard normal random variates, and BB T is the process noise covariance | see section 5. 3. Measure and weight the new position in terms of the measured features zt : t(n) = p(zt jxt = s(tn) ) P then normalise so that n t(n) = 1 and store together with cumulative probability as (s(tn) ; t(n) ; c(tn) ) where c(0) t = 0; (n) ct = c(tn?1) + t(n) (n = 1; : : : ; N ): Once the N samples have been constructed: estimate, if desired, moments of the tracked position at time-step t as E [f (xt )] = N X n=1 t(n) f st(n) obtaining, for instance, a mean position using f (x) = x. Figure 6: The Condensation algorithm. Isard and Blake. IJCV, in press, (1998). 12 et al., 1993) to restrict the conguration of the spline to a shape-space of vectors X dened by ! ! Qx = W X + Qxy ; Qy Q (9) where the matrix W is a matrix of rank NX considerably lower than the 2NB degrees of freedom of the unconstrained spline. Typically the shape-space may allow ane deformations of the template shape Q, or more generally a space of rigid and non-rigid deformations. The space is constructed by applying an appropriate combination of three methods to build a W -matrix: 1. determining analytically combinations of contours derived from one or more views (Ullman and Basri, 1991; Koenderink and Van Doorn, 1991; Blake et al., 1993), a method that is usable both for ane spaces and for certain classes of articulated object; 2. capturing sequences of key frames of the object in dierent poses (Blake et al., 1995); 3. performing principal components analysis on a set of outlines of the deforming object (Cootes et al., 1993; Baumberg and Hogg, 1994) to derive a small set of representative contours. 5.2 Dynamical model Exploiting earlier work on dynamical modelling (Blake et al., 1993; Blake et al., 1995), object dynamics are modelled as a 2nd order process, conveniently represented in discrete time t as a 2nd order linear dierence equation: xt ? x = A(xt? ? x ) + B wt 1 (10) where wt are independent vectors of independent standard normal variables, the state-vector xt = XXt?t 1 ! ; (11) and where x is the mean value of the state and A; B are matrices representing the deterministic and stochastic components of the dynamical model respectively. The system is a set of damped oscillators, whose modes, natural frequencies and damping constants are determined by A, driven by random accelerations coupled into the dynamics via B from the noise term B w. While it is possible to set sensible defaults for A, x and B , it is more satisfactory and eective to estimate them from input data taken while the object performs typical motions. Methods for doing this via Maximum Likelihood Estimation are essential to the work described here and are described fully elsewhere (Blake et al., 1995; Reynard et al., 1996). The dynamical model can be re-expressed in such a way as to make quite clear that it is a temporal Markov chain: (12) p(xt jxt?1 ) / exp ? 21 kB ?1 ((xt ? x) ? A(xt?1 ? x))k2 where k : : : k is the Euclidean norm. It is therefore clear that the learned dynamical models are appropriate for use in the Condensation algorithm. Isard and Blake. IJCV, in press, (1998). 13 5.3 Initial conditions Initial conditions for tracking can be determined by specifying the prior density p(x0 ), and if this is Gaussian, direct sampling can be used to initialise the Condensation algorithm. Alternatively it is possible simply to allow the density p(xt ) to settle to a steady state p(x1 ), in the absence of object measurements. Provided the learned dynamics are stable (free of undamped oscillations) a unique steady state exists. Furthermore, if p(x0 ) is Gaussian, p(x1 ) is Gaussian with parameters that can be computed by iterating the Riccati equation (Gelb, 1974). At this point the density function represents an envelope of possible congurations of the object, as learned during the training phase. (Background clutter, if present, will modify and bias this envelope to some extent.) Then, as soon as the foreground object arrives and is measured, the density p(xt ) begins to evolve appropriately. 6 Observation model The observation process dened by p(zt jxt ) is assumed here to be stationary in time (though the Condensation algorithm does not necessarily demand this) so a static function p(zjx) needs to be specied. As yet we have no capability to estimate it from data, though that would be ideal, so some reasonable assumptions must be made. First a measurement model for one-dimensional data with clutter is suggested. Then an extension is proposed for two-dimensional observations that is also used later in computational experiments. 6.1 One-dimensional observations in clutter In one dimension, observations reduce to a set of scalar positions fz = (z1 ; z2 ; : : : ; zM )g and the observation density has the form p(zjx) where x is one-dimensional position. The multiplicity of measurements reects the presence of clutter so either one of the events m = ftrue measurement is zm g; m = 1; : : : ; M P occurs, or else the target object is not visible with probability q = 1 ? m P (m ). Such reasoning about clutter and false alarms is commonly used in target tracking (Bar-Shalom and Fortmann, 1988). Now the observation density can be expressed as p(zjx) = qp(zjclutter) + M X m=1 p(zjx; m )P (m ): A reasonable functional form for this can be obtained by making some specic assumptions: that3 P (m ) = p; 8 m, that the clutter is a Poisson process along the line with spatial density and that any true target measurement is unbiased and normally distributed with standard deviation . This leads to 2 X exp ? 2m2 (13) p(zjx) / 1 + p 1 2 m where = q and m = zm ? x, and is illustrated in gure 7. Peaks in the density function There could be some benet in allowing the P (m ) to vary with m to reect varying degrees of feature-anity, based on contrast, colour or orientation. 3 Isard and Blake. Measured features IJCV, in press, (1998). z1 z2 z3 14 z4 z5 z6 x p(z | x) x Figure 7: One-dimensional observation model. A probabilistic observation model allowing for clutter and the possibility of missing the target altogether is specied here as a conditional density p(zjx). correspond to measured features and the state density will tend to be reinforced in the Condensation algorithm at such points. The background level reects the possibility that the true target has not been detected at all. The eect on tracking behaviour is to provide for the possibility of \tunneling": a good hypothesis should survive a transitory failure of observations due, for example, to occlusion of the tracked object. The parameters (units of distance) and (units of inverse distance) must be chosen, though in principle they could be estimated from data by observing measurement error and both the density of clutter and probability of non-detection q. Considerable economy can be applied, in practice, in the evaluation of the observation density. Given a hypothesised position x in the \observation" step (gure 6) it is not necessary to attend to all features z1 ; : : : ; zM . Any m for which p1 exp ? 2m2 1 2 2 can be neglected and this sets a search window around the position x outside which measurements can be ignored. For practical values of the constants the search window will have a width of a few . In practice the clutter is suciently sparse and is suciently small that the search window rarely contains more than one feature. Note that the density p(zjx) represents the information about x given a xed number M of measurements. Potentially, the event M that there are M measurements, regardless of the actual values of those measurements, also constitutes information about x. However, we can reasonably assume here that P ( M jx) = P ( M ); Isard and Blake. IJCV, in press, (1998). 15 for instance because x is assumed to lie always within the image window. In that case, by Bayes' theorem, p(xj M ) = p(x) | the event M provides no additional information about the position x. (If x is allowed also to fall outside the image window then the event M is informative: a value of M well above the mean value for the background clutter enhances the probability that x lies within the window.) 6.2 Two-dimensional observations In a two-dimensional image, the set of observations z is, in principle, the entire set of features visible in the image. However, an important aspect of earlier systems in achieving real-time performance (Lowe, 1992; Harris, 1992; Blake et al., 1993) has been the restriction of measurement to a sparse set of lines normal to the tracked curve. These two apparently conicting ideas can be resolved as follows. The observation density p(zjx) in two dimensions describes the distribution of a (linearly) parameterised image curve z(s), given a hypothetical shape in the form of a curve r(s); 0 s 1, represented by a shape parameter x. The two-dimensional density be derived as an extension of the one-dimensional case. It is assumed that a mapping g(s) is known that associates each point z(s) on the image curve with a point r(g(s)) on the shape. In practice this mapping is set up by tracing normals from the curve r. Note that g(s) is not necessarily injective because z(s) includes clutter as well as foreground features. Next the one-dimensional density (13) is approximated in a more amenable form that neglects the possibility of more than one feature lying inside the search interval: (14) p(zjx) / exp ? 21 2 f (1; ) where f ( ; ) = min( 2 ; 2 ); p p = 2 log(1= 2) is a spatial scale constant, and 1 is the m with smallest magnitude, representing the feature lying closest to the hypothesised position x. A natural extension to two dimensions is then ZL (15) p(zjx) = Z exp ? 21r f (z1 (s) ? r(s); ) ds 0 in which r is a variance constant and z1 (s) is the closest associated feature to r(s): z (s) = z(s0) where s0 = arg s02min jr(s) ? z(s0)j: g?1 s 1 ( ) Note that the constant of proportionality (\partition function") Z (x) is an unknown function. We make the assumption that the variation of Z with x is slow compared with the other term in (15) so that Z can be treated as constant. It remains to establish whether this assumption is justied. The observation density (15) can be computed via a discrete approximation, the simplest being: ) ( M X 1 (16) p(zjx) / exp ? 2rM f (z1 (sm) ? r(sm); ) ; m=1 p where sm = m=M . This is simply the product of one-dimensional densities (14) with = rM , evaluated independently along M curve normals as in gure 8. Isard and Blake. IJCV, in press, (1998). 16 Figure 8: Observation process. The thick line is a hypothesised shape, represented as a parametric spline curve. The spines are curve normals along which high-contrast features (white crosses) are sought. Isard and Blake. IJCV, in press, (1998). 17 7 Applying the Condensation algorithm to video-streams Four examples are shown here of the practical ecacy of the Condensation algorithm. Movie (MPEG) versions of some results are available on the web at http://www.robots.ox.ac.uk/~ab/. 7.1 Tracking a multi-modal distribution The ability of the Condensation algorithm to represent multi-modal distributions was tested using a 70 frame (2.8 second) sequence of a cluttered room containing three people each facing the camera (gure 9). One of the people moves from right to left, in front of the other two. The Figure 9: Tracking three people in a cluttered room. The rst frame of a sequence in which one gure moves from right to left in front of two stationary gures. shape-space for tracking is built from a hand-drawn template of head and shoulders (gure 8) which is then allowed to deform via planar ane transformations . A Kalman lter contourtracker (Blake et al., 1993) with default motion parameters is able to track a single moving person just well enough to obtain a sequence of outline curves that is usable as training data. Given the high level of clutter, adequate performance with the Kalman lter is obtained here by means of background modelling (Rowe and Blake, 1996), a statistical form of background subtraction, which eectively removes clutter from the image data before it is tracked. It transpires, for this particular training set, that the learned motions comprise primarily horizontal translation, with vertical translation and horizontal and vertical shear present to a lesser degree. The learned shape and motion model can now be installed as p(xt jxt?1 ) in the Condensation algorithm which is now run on a test sequence but without the benet of background modelling, so that the background clutter is now visible to the tracker. Figure 10 shows how the state-density evolves as tracking progresses. Initialisation is performed simply by iterating the stochastic model, in the absence of measurements, to its steady state and it can be seen that this corresponds, at time 0, to a roughly Gaussian distribution, as expected. The distribution rapidly Isard and Blake. IJCV, in press, (1998). 18 0 ms 400 ms 800 ms 1200 ms 1600 ms 2000 ms 2400 ms Time 2800 ms Figure 10: Tracking with multi-modal state-density. An approximate depiction of the (2) state-density is shown, computed by smoothing the distribution of point masses s(1) t ; st ; : : : in the Condensation algorithm. The density is, of course, multi-dimensional; its projection onto the horizontal translation axis is shown here. The initial distribution is roughly Gaussian but this rapidly evolves to acquire peaks corresponding to each of the three people in the scene. The right-most peak drifts leftwards, following the moving person, coalescing with and separating from the other two peaks as it moves. Having specied a tracker for one person we eectively have, for free, a multi-person tracker, owing to the innate ability of the Condensation algorithm to maintain multiple hypotheses. Isard and Blake. IJCV, in press, (1998). 19 collapses down to three peaks which are then maintained appropriately even during temporary occlusion. Although the tracker was designed to track just one person, the Condensation algorithm allows the tracking of all three, for free; the ability to represent multi-modal distributions eectively provides multiple hypothesis capability. Tracking is based on frame rate (40 ms) sampling in this experiment and distributions are plotted in the gure for alternate frames. The experiment was run using a distribution of N = 1000 samples per time-step. 7.2 Tracking rapid motions through clutter The ability to track more agile motion, still against clutter, was tested using a 500 eld (10 second) sequence of a girl dancing vigorously to a Scottish reel. The shape-space for tracking was planar ane, based on a hand-drawn template curve for the head outline. The training sequence consisted of dancing against a largely uncluttered background, tracked by a Kalman lter contour-tracker with default dynamics to record 140 elds (2.8 seconds) of tracked head positions, the most that could be tracked before losing lock. Those 140 elds were sucient to learn a bootstrap motion model which then allowed the Kalman lter to track the training data for 800 elds (16 seconds) before loss of lock. The motion model obtained from these 800 elds was used in experiments with the Condensation tracker and applied to the test data, now including clutter. Figure 11 shows some stills from the test sequence, with a trail of preceding head positions to indicate motion. The motion is primarily translation, with some horizontal shear apparent as the dancer turns her head. Representing the state density with N = 100 samples at each time-step proves just sucient for successful tracking. As in the previous example, a prior density can be computed as the steady state of the motion model and, in this case, that yields a prior for position that spreads across most of the image area, as might be expected given the range of the dance. Such a broad distribution cannot eectively be represented by just N = 100 samples. One alternative is to increase N in the early stages of tracking, and this is done in a later experiment. Alternatively, the prior can be based on a narrower distribution whose centre is positioned by hand over the object at time 0, and that is what was done here. Observation parameters were = 24; = 7 with M = 18 normals. Figure 12 shows the motion of the centroid of the estimated head position as tracked both by the Condensation algorithm and by a Kalman lter using the same motion model. The Condensation tracker correctly estimated head position throughout the sequence, but after about 40 elds (0.80 s), the Kalman lter was distracted by clutter, never to recover. Given that there is only one moving person in this experiment, unlike the previous one in which there were three, it might seem that a unimodal representation of the state density should suce. This is emphatically not the case. The facility to represent multiple modes is crucial to robustness as gure 13 illustrates. The gure shows how the distribution becomes misaligned (at 900ms), reacting to the distracting form of the computer screen. After a further 20ms the distribution splits into two distinct peaks, one corresponding to clutter (the screen), one to the dancer's head. At this point the clutter peak actually has the higher posterior probability | a unimodal tracker, for instance a Kalman lter, would almost certainly discard the lower peak, rendering it unable to recover. The Condensation algorithm however, capable as it is of carrying several hypotheses simultaneously, does recover rapidly as the clutter peak decays for lack of conrmatory observation, leaving just one peak corresponding to the dancer at 960 ms. Isard and Blake. IJCV, in press, (1998). 20 eld 91 (1820 ms) eld 121 (2420 ms) eld 221 (4420 ms) eld 265 (5300 ms) Figure 11: Tracking agile motion in clutter. The test sequence consists of 500 elds (10 seconds) of agile dance against a cluttered background. The dancer's head is tracked through the sequence. Several representative elds are shown here, each with a trail of successive mean tracked head positions at intervals of 40 ms. The Condensation algorithm used N = 100 samples per time-step to obtain these results. Isard and Blake. IJCV, in press, (1998). 21 Y Y X Time = 10 s X Time = 10 s Condensation tracker Kalman filter tracker Figure 12: The Condensation tracker succeeds where a Kalman lter fails. The estimated centroid for the sequence shown in gure 11 is plotted against time for the entire 500 eld sequence, as tracked rst by the Condensation tracker, then by a comparable Kalman lter tracker. The Condensation algorithm correctly estimates the head position throughout the sequence. The Kalman lter tracks briey, but is soon distracted by clutter and never recovers. 7.3 Tracking an articulated object The preceding sequences show motion taking place in ane shape-spaces of just 6 dimensions. High dimensionality is one of the factors, in addition to agility and clutter, that makes tracking hard (Blake et al., 1993). In order to investigate tracking performance in higher dimensions, we used a 500 eld (10 second) test sequence of a hand translating, rotating, and exing its ngers independently, over a highly cluttered desk scene (gure 14). Figure 15 shows just how severe the clutter problem is | the hand is immersed in a dense eld of edges. A model of shape and motion model was learned from training sequences of hand motion against a plain background, tracked by Kalman lter (using signed edges to help to disambiguate nger boundaries). The procedure comprised several stages, creative assembly of methods from the available \toolkit" for learning (Blake et al., 1995). 1. Shape-space was constructed from 6 templates drawn around the hand with the palm in a xed orientation and with the ngers and thumb in various congurations. The 6 templates combined linearly to form a 5-dimensional space of deformations which were then added to the space of translations to form a 7 dimensional shape-space. 2. Default dynamics in the shape-space above were adequate to track a clutter-free training sequence of 600 frames in which the palm of the hand maintained an approximately xed attitude. 3. Principal components analysis: the sequence of 600 hand outlines was replicated with each hand contour rotated through 90 degrees and the sequences concatenated to give a Isard and Blake. IJCV, in press, (1998). 22 eld 45 (900 ms) eld 46 (920 ms) eld 47 (940 ms) eld 48 (960 ms) Figure 13: Recovering from tracking failure. Detail from 4 consecutive elds of the sequence illustrated in gure 11. Each sample from the distribution is plotted on the image, with intensity scaled to indicate its posterior probability. (Most of the N = 100 samples have too low a probability to be visible in this display.) In eld 45, the distribution is misaligned, and has begun to diverge. In elds 46 and 47 it has split into two distinct peaks, the larger attracted to background clutter, but converges back onto the dancer in eld 48. Isard and Blake. IJCV, in press, (1998). 23 Figure 14: A hand moving over a cluttered desk. Field 0 of a 500 eld (10 second) sequence in which the hand translates, rotates, and the ngers and thumb ex independently. Edge detector Condensation Figure 15: Severe clutter. Detail of one eld (gure 14) from the test-sequence shows the high level of potential ambiguity. Output from a directional Gaussian edge detector shows that there are many clutter edges present as potential distractors. Isard and Blake. IJCV, in press, (1998). 24 sequence of 1200 deformations. Projecting out the translational component of motion, the application of Principal Component Analysis (PCA) to the sequence of residual deformations of the 1200 contours established a 10-dimensional space that accounted almost entirely for deformation. This was then combined with the translational space to form a 12-dimensional shape-space that accounted both for the exing of ngers and thumb and also for rotations of the palm. 4. Bootstrapping: a Kalman lter with default dynamics in the 12-dimensional shape-space was sucient to track a training sequence of 800 elds of the hand translating, rotating, and exing ngers and thumb slowly. This was used to learn a model of motion. 5. Re-learning: that motion model was installed in a Kalman lter used to track another, faster training-sequence of 800 elds. This allowed a model for more agile motion to be learned, which was then used in experiments with the Condensation tracker. Figure 16: Tracking a exing hand across a cluttered desk. Representative stills from a 500 eld (10 second) sequence show a hand moving over a highly cluttered desk scene. The ngers and thumb ex independently, and the hand translates and rotates. Here the Condensation algorithm uses N = 1500 samples per time-step initially, dropping gradually over 4 elds to N = 500 for the tracking of the remainder of the sequence. The mean conguration of the contour is displayed. Figure 16 shows detail of a series of images from a tracked, 500 eld test-sequence. The initial state density was simply the steady state of the motion model, obtained by allowing the lter to iterate in the absence of observations. Tracker initialisation was facilitated by using more samples per time-step (N = 1500) at time t = 0, falling gradually to 500 over the rst 4 elds. The rest of the sequence was tracked using N = 500. As with the previous example of the dancer, clutter can distract the tracker but the ability to represent multi-modal state density means that tracking can recover. 7.4 Tracking a camouaged object Finally, we tested the ability of the algorithm to track rapid motion against background distraction in the extreme case that background objects actually mimiced the tracked object. A 12 second (600 eld) sequence showed a bush blowing in the wind, the task being to track one particular leaf. A template was drawn by hand around a still of one chosen leaf and allowed to Isard and Blake. IJCV, in press, (1998). 25 undergo ane deformations during tracking. Given that a clutter-free training sequence was not available, the motion model was again learned by means of a bootstrap procedure. A tracker with default dynamics proved capable of tracking the rst 150 elds of a training sequence before losing the leaf, and those tracked positions allowed a rst approximation to the model to be learned. Installing that in a Condensation tracker, the entire sequence could be tracked, though with occasional misalignments. Finally a third learned model was sucient to track accurately the entire 12-second training sequence. Despite occasional violent gusts of wind and temporary obscuration by another leaf, the Condensation algorithm successfully followed the object, as gure 17 shows. In fact, tracking is accurate enough using N = 1200 samples to 1.46 s 2.66 s 5.54 s 7.30 s Figure 17: Tracking with camouage. The aim is to track a single camouaged moving leaf in this 12-second sequence of a bush blowing in the wind. Despite the heavy clutter of distractors which actually mimic the foreground object, and occasional violent gusts of wind, the chosen foreground leaf is successfully tracked throughout the sequence. Representative stills depict mean contour congurations, with preceding tracked leaf positions plotted at 40 ms intervals to indicate motion. separate the foreground leaf from the background reliably, an eect which can otherwise only be achieved using \blue-screening". Having obtained the model iteratively as above, independent Isard and Blake. IJCV, in press, (1998). 26 test sequences could be tracked without further training. With N = 1200 samples per time-step the tracker runs at 6.5 Hz on a SGI Indy SC4400 200MHz workstation. Reducing this to N = 200 increases processing speed to video frame-rate (25 Hz), at the cost of occasional misalignments in the mean conguration of the contour. Observation parameters were = 8; = 3 with M = 21 normals. 8 Conclusions Tracking in clutter is hard because of the essential multi-modality of the conditional observation density p(zjx). In the case of curves multiple-hypothesis tracking is inapplicable and a new approach is needed. The Condensation algorithm is a fusion of the statistical factored sampling algorithm for static, non-Gaussian problems with a stochastic model for object motion. The result is an algorithm for tracking rigid and non-rigid motion which has been demonstrated to be far more eective in clutter than comparable Kalman lters. Performance of the Condensation algorithm improves as the sample size parameter N increases; formally computational complexity is O(N log N ), although this can be made O(N ) with a minor modication to the sampling procedure. Impressive results have been demonstrated for models with between 6 and 12 degrees of freedom, even when N is as low as 100{200. Performance in several cases was improved still further with an increased value N 1000. In a 6-dimensional shape-space, the system currently runs with N = 100 in real-time (50Hz) on a desk-top graphics workstation (SGI Indy R4400SC, 200 MHz). The new approach raises a number of questions. First, alternative observation models could be explored in order to make greater use of image intensity variations, though without sacricing too much in the way of photometric invariance. It is to be hoped in the interests of eciency that, as happens with the search window in the edge-based case, computational attention could be concentrated in a band around the hypothesised curve without signicant loss of accuracy in the model. Such a model would have echoes of correlation matching but of course without the exhaustive search characteristic of correlation matchers which is quite infeasible in more than two or three dimensions. Secondly, the availability of general state densities suggests the need for more general representations of those densities. When the density is approximately unimodal, rst and second moments may be adequate to convey the likely states, but in the multi-modal case, as for example when several people are tracked simultaneously, the mean conguration is not a particularly useful statistic | it meaninglessly combines the congurations of the three people. An alternative is to attempt to develop a mode nder capable of pin-pointing several modes when present. More generally there is a need for \operators" to interrogate densities: for instance, an operator to nd a person moving to the right, or to nd the tallest person. Perhaps such operators could be formulated as hypothesis tests applied to sample sets. A third question concerns the random sampling scheme and its eciency. Factored sampling can be inecient as the modes of p(zjx) become narrow. One approach is \importance sampling" (Ripley, 1987) in which a heuristically chosen distribution, approximating p(zjx), is used to concentrate random sampling around modes. However, this has the drawback that the prior p(x) must be repeatedly evaluated whereas, in temporal propagation, the prior (prediction) p(xt jzt?1 ) cannot be evaluated pointwise, only sampled. Fourthly and nally, it is striking that the density propagation equation (4) in the Condensation algorithm is a continuous form of the propagation rule of the \forward algorithm" Isard and Blake. IJCV, in press, (1998). 27 for Hidden Markov Models (HMMs) (Rabiner and Bing-Hwang, 1993). The integral over continuous states in (5) becomes a summation over discrete states in the HMM, with p(xt jxt?1 ) represented by a transition matrix. This suggests a natural opportunity to combine the two so that mixed discrete/continuous states could be propagated over time. This would allow switching between multiple models, for instance walk-trot-canter-gallop, each model represented by a stochastic dierential equation, with transitions governed by a discrete conditional probability matrix. It seems likely that such a system could be executed as a Condensation tracker. A further challenge is to develop a learning algorithm for mixed dynamical models. Acknowledgements The authors would like to acknowledge the support of the EPSRC. They are also grateful for discussions with Roger Brockett, John Kent, Ian Reid, David Mumford, Brian Ripley, David Reynard, Simon Rowe, Andrew Wildenberg and Andrew Zisserman, and for experimental assistance from Sarah Blake. Isard and Blake. IJCV, in press, (1998). 28 References Anderson and Moore (1979). Optimal ltering. Prentice Hall. Astrom, K. J. (1970). Introduction to stochastic control theory. Academic Press. Bar-Shalom, Y. and Fortmann, T. (1988). Tracking and Data Association. Academic Press. Bartels, R., Beatty, J., and Barsky, B. (1987). An Introduction to Splines for use in Computer Graphics and Geometric Modeling. Morgan Kaufmann. Baumberg, A. and Hogg, D. (1994). Learning exible models from image sequences. In Eklundh, J.-O., editor, Proc. 3rd European Conference on Computer Vision, 299{308. Springer{ Verlag. Baumberg, A. and Hogg, D. (1995). Generating spatiotemporal models from examples. In Proc. of the British Machine Vision Conference, 2, 413{422. Blake, A., Curwen, R., and Zisserman, A. (1993). A framework for spatio-temporal control in the tracking of visual contours. Int. Journal of Computer Vision, 11, 2, 127{145. Blake, A. and Isard, M. (1994). 3D position, attitude and shape input using video tracking of hands and lips. In Proc. Siggraph, 185{192. ACM. Blake, A., Isard, M., and Reynard, D. (1995). Learning to track the visual motion of contours. Articial Intelligence, 78, 101{134. Bucy, R. (1969). Bayes theorem and digital realizations for non-linear lters. J. Astronautical Sciences, 17, 2, 80{94. Cipolla, R. and Blake, A. (1990). The dynamic analysis of apparent contours. In Proc. 3rd Int. Conf. on Computer Vision, 616{625. Cootes, T., Taylor, C., Lanitis, A., Cooper, D., and Graham, J. (1993). Building and using exible models incorporating grey-level information. In Proc. 4th Int. Conf. on Computer Vision, 242{246. Dickmanns, E. and Graefe, V. (1988). Applications of dynamic monocular machine vision. Machine Vision and Applications, 1, 241{261. Fischler, M. and Bolles, R. (1981). Random sample consensus: a paradigm for model tting with application to image analysis and automated cartography. Commun. Assoc. Comp. Mach., 24, 381{95. Fischler, M. A. and Elschlager, R. A. (1973). The representation and matching of pictorial structures. IEEE. Trans. Computers, C-22, 1. Gelb, A., editor (1974). Applied Optimal Estimation. MIT Press, Cambridge, MA. Geman, S. and Geman, D. (1984). Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images. IEEE Trans. Pattern Analysis and Machine Intelligence, 6, 6, 721{ 741. Isard and Blake. IJCV, in press, (1998). 29 Gennery, D. (1992). Visual tracking of known three-dimensional objects. Int. Journal of Computer Vision, 7:3, 243{270. Goodwin, C. and Sin, K. (1984). Adaptive ltering prediction and control. Prentice-Hall. Gordon, N., Salmond, D., and Smith, A. (1993). Novel approach to nonlinear/non-gaussian bayesian state estimation. IEE Proc. F, 140, 2, 107{113. Grenander, U., Chow, Y., and Keenan, D. M. (1991). HANDS. A Pattern Theoretical Study of Biological Shapes. Springer-Verlag. New York. Hager, G. (1990). Sensor fusion and planning: a computational approach. Kluwer Academic Publishers. Harris, C. (1992). Tracking with rigid models. In Blake, A. and Yuille, A., editors, Active Vision, 59{74. MIT. Hogg, D. (1983). Model-based vision: a program to see a walking person. Image and Vision Computing, 1, 1, 5{20. Huttenlocher, D., Noh, J., and Rucklidge, W. (1993). Tracking non-rigid objects in complex scenes. In Proc. 4th Int. Conf. on Computer Vision, 93{101. Isard, M. and Blake, A. (1996). Visual tracking by stochastic propagation of conditional density. In Proc. 4th European Conf. on Computer Vision, 343{356, Cambridge, England. Kass, M., Witkin, A., and Terzopoulos, D. (1987). Snakes: Active contour models. In Proc. 1st Int. Conf. on Computer Vision, 259{268. Kitagawa, G. (1996). Monte-carlo lter and smoother for non-Gaussian nonlinear state space models. J. Computational and Graphical Statistics, 5, 1, 1{25. Koenderink, J. and Van Doorn, A. (1991). Ane structure from motion. J. Optical Soc. of America A., 8, 2, 337{385. Lowe, D. (1991). Fitting parameterised 3D models to images. IEEE Trans. Pattern Analysis and Machine Intelligence, 13:5, 441{450. Lowe, D. (1992). Robust model-based motion tracking through the integration of search and estimation. Int. Journal of Computer Vision, 8, 2, 113{122. Matthies, L. H., Kanade, T., and Szeliski, R. (1989). Kalman lter-based algorithms for estimating depth from image sequences. Int. Journal of Computer Vision, 3, 209{236. Menet, S., Saint-Marc, P., and Medioni, G. (1990). B-snakes: implementation and application to stereo. In Proceedings DARPA, 720{726. Miller, M., Srivasta, A., and Grenander, U. (1995). Conditional-mean estimation via jumpdiusion processes in multiple target tracking/recognition. IEEE Transactions on Signal Processing, 43, 11, 2678{2690. Papoulis, A. (1990). Probability and Statistics. Prentice-Hall. Isard and Blake. IJCV, in press, (1998). 30 Press, W., Teukolsky, S., Vetterling, W., and Flannery, B. (1988). Numerical Recipes in C. Cambridge University Press. Rabiner, L. and Bing-Hwang, J. (1993). Fundamentals of speech recognition. Prentice-Hall. Rao, B. (1992). Data association methods for tracking systems. In Blake, A. and Yuille, A., editors, Active Vision, 91{105. MIT. Rao, C. (1973). Linear Statistical Inference and Its Applications. John Wiley and Sons, New York. Rehg, J. and Kanade, T. (1994). Visual tracking of high dof articulated structures: an application to human hand tracking. In Eklundh, J.-O., editor, Proc. 3rd European Conference on Computer Vision, 35{46. Springer{Verlag. Reynard, D., Wildenberg, A., Blake, A., and Marchant, J. (1996). Learning dynamics of complex motions from image sequences. In Proc. 4th European Conf. on Computer Vision, 357{368, Cambridge, England. Ripley, B. (1987). Stochastic simulation. New York: Wiley. Ripley, B. and Sutherland, A. (1990). Finding spiral structures in images of galaxies. Phil. Trans. R. Soc. Lond. A., 332, 1627, 477{485. Rowe, S. and Blake, A. (1996). Statistical feature modelling for active contours. In Proc. 4th European Conf. on Computer Vision, 560{569, Cambridge, England. Sorenson, H. W. and Alspach, D. L. (1971). Recursive bayesian estimation using gaussian sums. Automatica, 7, 465{479. Storvik, G. (1994). A Bayesian approach to dynamic contours through stochastic sampling and simulated annealing. IEEE Trans. Pattern Analysis and Machine Intelligence, 16, 10, 976{986. Sullivan, G. (1992). Visual interpretation of known objects in constrained scenes. Phil. Trans. R. Soc. Lond. B., 337, 361{370. Terzopoulos, D. and Metaxas, D. (1991). Dynamic 3D models with local and global deformations: deformable superquadrics. IEEE Trans. Pattern Analysis and Machine Intelligence, 13, 7. Terzopoulos, D. and Szeliski, R. (1992). Tracking with Kalman snakes. In Blake, A. and Yuille, A., editors, Active Vision, 3{20. MIT. Ullman, S. and Basri, R. (1991). Recognition by linear combinations of models. IEEE Trans. Pattern Analysis and Machine Intelligence, 13, 10, 992{1006. Yuille, A. and Hallinan, P. (1992). Deformable templates. In Blake, A. and Yuille, A., editors, Active Vision, 20{38. MIT. Isard and Blake. IJCV, in press, (1998). 31 A Non-linear ltering There are four distinct probability distributions represented in a non-linear Bayesian lter. Three of them form part of the problem specication and the fourth constitutes the solution. The three specied distributions are: 1. the prior density p(x) for the state x 2. the process density p(xt jxt?1 ) that describes the stochastic dynamics 3. the observation density p(zjx) and the lter evolves over time to generate, as the solution at each time-step, the state-density pt(x) where pt (xt ) p(xt jZt ). Only when all of the three specied distributions are Gaussian is the state-density pt also Gaussian. Otherwise, for non-Gaussian pt , it is possible to use one of a number of approximate lters, depending on which of the specied densities it is that is non-Gaussian. A.1 Non-Gaussian prior density The case that the prior density is non-Gaussian is the simplest to deal with provided only that it can adequately be represented (or approximated) as an additive Gaussian mixture: p0 (x) = M X m=1 w(m) G(x; (m) ; P0(m) ): In that case, provided that other specied densities are Gaussian, the state density can also be represented as a corresponding mixture pt (x) = M X m=1 w(m) G(x; t(m) ; Pt(m) ) in which the means t(m) and variances Pt(m) vary over time but the weights w(m) are xed. Each of the M mixture components evolves as an independent Gaussian so that, in fact, the state density is just a sum of densities from M independent linear Kalman lters. A.2 Non-Gaussian process density Non-Gaussian state densities can arise from the nature of the process either because the dynamics are driven by non-Gaussian process noise, or, more generally, because the deterministic dynamics are non-linear. One approach to ltering is then to approximate the dynamics by Taylor expansion as a linear process with time-varying coecients and proceed as for linear Kalman lters. This generates a Gaussian representation of the evolving state-density which may be a good approximation depending on the nature of the non-linearity. This is the basis of the \Extended Kalman Filter" (EKF) (Gelb, 1974; Bar-Shalom and Fortmann, 1988). Alternatively, one can attempt a mixture representation, as earlier, but now allowing the weights w(m) also to vary over time. Unfortunately, even allowing dynamic re-weighting (Sorenson and Alspach, 1971) does not produce exact solutions for pt (x), because the individual Gaussian components do not remain Gaussian over time. For example, consider the case in which the process Isard and Blake. IJCV, in press, (1998). 32 density p(xt jxt?1 ) is itself an additive mixture of k > 1 Gaussian components. According to the Bayesian propagation equation (5) each component of pt splits into k separate components in the transition from time n to time n + 1; the total number of components in pt grows exponentially as kt . Clearly pt must be approximated at each time-step to prune back the number of components (Anderson and Moore, 1979) within some resource-limited bound M . Eectively there are Mk full Kalman lters running at each time-step, each bringing the computational expense of its own Riccati equation step. Clearly the success of this approach depends on how well the densities pt and p(xt jxt?1 ) can be approximated with a modest number Mk of components. A.3 Non-Gaussian observation density In the case of visual tracking in clutter, non-linearity of the tracking lter arises, as we have seen, because the observation density p(zjx) is non-Gaussian and, furthermore, is multi-modal so that it cannot be well approximated by a single Gaussian. Each of the methods just mentioned for handling non-Gaussian process density, the EKF and Gaussian mixtures, are relevant also to non-Gaussian observation density but continue to have the same drawbacks. Note that, in the case of Gaussian mixtures, the number of mixture components again proliferates at each time-step of (4), albeit via a dierent mechanism involving products of Gaussians rather than convolutions. Even this assumes that the observation density can be approximated as a mixture but in clutter this becomes rather inecient, requiring at least one component per visible feature. There is an additional class of techniques which applies to this case when the non-Gaussian state density arises from clutter of a particular sort. In the simplest case, one of a nite set of measurements zt = fzt;1 ; : : : ; zt;k g at time n is to be associated with the state xt at time t. Heuristic mechanisms such as the validation gate and the probabilistic data-association lter (PDAF) (Bar-Shalom and Fortmann, 1988) attempt to deal with the ambiguity of association. Alternatively it can, in principle, be dealt with exactly by \multiple hypothesis ltering" but with computational cost that grows exponentially over time and which is therefore ruled out in practice. The \RANSAC" algorithm (Fischler and Bolles, 1981) deals probabilistically with multiple observations but the observations have to be discrete, and there is no mechanism for temporal propagation. More complex methods including the Joint PDAF (JPDAF) (Bar-Shalom and Fortmann, 1988; Rao, 1992) address the more dicult problem of associating not simply single features but subsequences of Zt with the state. However, these methods rely on the existence of discrete features. In contour tracking the features are continuous curves and so are not naturally amenable to discrete association. A.4 Direct integration Finally, one very general approach to nonlinear ltering must be mentioned. This is simply to integrate (5) directly, using a suitable numerical representation of the state density such as nite elements. This in essence is what (Bucy, 1969) proposed and more recently (Hager, 1990) investigated with respect to robotics applications. It is usable in one or two dimensions but, complexity being exponential in the dimension, is altogether infeasible for problems of dimension around 6{20, typical of the tracking problems dealt with here. The Condensation algorithm is designed to oer a viable alternative. Isard and Blake. IJCV, in press, (1998). 33 B Derivation of the sampling rule The correctness of the sampling rule (4) on page 6 is proved by rst deriving two lemmas from the independence assumption (2). (This is similar to the derivation found in (Bar-Shalom and Fortmann, 1988), except that our independence assumptions are explicitly specied.) Lemma 1 p(zt jXt ; Zt?1 ) = p(zt jxt ): Proof: p(Zt jXt ) = p(zt ; Zt?1 jXt ) = p(zt jZt?1 ; Xt )p(Zt?1 jXt ) = p(zt jZt?1 ; Xt ) tY ?1 i=1 p(zijxi ): (Taking (3) at time t and integrating w.r.t. zt yields the reduction of the second term in line 2.) Now, using (3) again gives the result. Lemma 2 Proof: from (2) so p(xt jXt?1 ; Zt?1 ) = p(xt jxt?1 ): p(xt ; Zt?1 jXt?1 ) = p(xt jXt?1 )p(Zt?1 jXt?1 ) p(xt jZt?1 ; Xt?1 ) = p(xt jXt?1 ) = p(xt jxt?1 ); using the Markov assumption (1). Derivation of the propagation formula: consider p(Xt jZt ) = p(zt jXt ;pZ(zt?jZ)p(X)t jZt? ) t t? = kt p(zt jXt ; Zt? )p(Xt jZt? ) = kt p(zt jxt )p(Xt jZt? ) using lemma 1: Now integrating w.r.t. Xt? gives p(xt jZt ) = kt p(zt jxt )p(xt jZt? ): 1 1 1 1 1 1 1 1 The last term can be expanded: p(xt jZt?1 ) = = = which is the required result. Z p(xt jXt?1 ; Zt?1 )p(Xt?1 jZt?1 ) ZXt?1 Z p(xt jxt?1 )p(Xt?1 jZt?1 ) Zxt?1 Xt?2 p(xt jxt?1 )p(xt?1 jZt?1 ) xt?1 using lemma 2 Isard and Blake. IJCV, in press, (1998). 34 C Asymptotic correctness of the Condensation Algorithm The Condensation algorithm is validated here by a probabilistic argument showing that the sample-set representation of conditional density is correct, asymptotically, as the size N of the sample set at each time-step gets large. The argument is based on the one by Grenander et al. to justify their factored sampling algorithm for interpretation of static images. They use the standard probabilistic tool of \weak convergence" (Rao, 1973) and the \weak law of large numbers" to show that a posterior distribution inferred by factored sampling can be made arbitrarily accurate by choosing N suciently large. No formal indication is given as to how large N should be for a given level of accuracy, something which is determined in practice by experimentation. In the proof that follows, the correctness proof for factored sampling of a static image is made inductive so that it can be applied to successive images in a sequence. This would be sucient to apply several independent images to the estimation of a static underlying object. A further generalisation takes account of the predictive step (step 2 of the Condensation algorithm) that deals with the dynamics of an object in motion. C.1 Factored sampling The asymptotic correctness of the factored sampling algorithm (section 3) is expressed in a theorem of Grenander et al (1991): Theorem 3 (Factored sampling) If p0pz is an (absolutely continuous) density function (with a suitable normalisation constant) then for any given value x p~(x) ! p0 (x)pz (x); weakly; as N ! 1 | pointwise, weak convergence of the density function to the required posterior. (Recall p~ is the density function of the random variate x generated by factored sampling, as dened in section 3.) The proof of the theorem was given by Grenander et al.. C.2 Dynamic extension of factored sampling The rst step in the extension for dynamic problems is to state a corollary of the theorem above that generalises it slightly to the case where the prior is not known exactly but has itself been simulated approximately. Cor. 4 (Weak factored sampling) The sequence s1; : : : ; sN is now generated by sampling from a density ps chosen such that ps(x) ! p0 (x); weakly; as N ! 1; where convergence is uniform with respect to x. Provided pz is bounded, the random variate x0 generated from the sn as before has a density function p~ for which p~(x) ! p0 (x)pz (x) weakly; as N ! 1 and convergence is uniform with respect to x. The proof of this corollary is straightforward. Isard and Blake. IJCV, in press, (1998). 35 C.3 Propagation of approximated state density First note that the samples s(tn) generated by the algorithm can themselves be regarded as random variables. Using the corollary it is possible to establish that asymptotically the probability density of any given s(tn) converges to the desired probability density p(xt jZt?1 ) . From now on the limit symbol ! is used to denote weak, uniform convergence of density functions as N ! 1. The correctness result is expressed in the theorem below. We rst require a normalisation assumption for the process density, that Z p(xt jxt?1 ) dxt?1 is bounded4 : (17) Theorem 5 (Weak propagation) Each sample stn ; n = 1; : : : ; N at time t is drawn from a ( ) distribution with density p~t such that p~t(xt ) ! p(xt jZt?1 ): Proof The proof is inductive. Suppose the result holds for p~t?1 ; then after step 1 of the algorithm in gure 6, by the corollary, and observing that the sampling probabilities are t(?n)1 / p(zt?1 jxt?1 = s(t?n)1 ); each s0 (t?n)1 has a density p0t?1 such that p0t?1 ! t?1 p(xt?1 jZt?2 )p(zt?1 jxt?1 ) where t?1 is a normalisation constant so that p0t?1 ! p(xt?1 jZt?1 ): In step 2 of the algorithm the random dynamical step is applied to s0 (tn) to give s(tn) with density p00 such that p00 (xt ) = = ! Z Z Z p(xt jxt?1 = s0 (tn) ) p(s0 (tn) ) ds0 (tn) p(xt jxt?1 )p0 (xt?1 ) dxt?1 p(xt jxt?1 )p(xt?1 jZt?1 ) dxt?1 (making use of (17)) = p(xt jZt?1 ) and this is the required density function for s(tn) , establishing the inductive step as required. Finally the ground instance is straightforward. Initial samples s0 (1n) are drawn in step 1 from the prior p0 so that, after step 2 of the algorithm, the s(1n) are sampled predictions for time t = 1 from a density p~1 such that p~1 (x1 ) = p(x1 ) p(x1 jZ0 ): it. 4 This assumption is not restrictive in practice but is a little inelegant and perhaps there is a way to do without Isard and Blake. IJCV, in press, (1998). (Z0 is an empty set) so certainly p~1 (x1 ) ! p(x1 jZ0 ) 36 as required. Note that convergence has not been proved to be uniform in t. For a given xed t, there is convergence as N ! 1 but nothing is said about the limit t ! 1. In practice this could mean that at later times t larger values of N may be required, though that could depend also on other factors such as the nature of the dynamical model.

© Copyright 2020