Multi-Target Tracking and Occlusion Handling with Learned Variational Bayesian Clusters and a Social Force Model

This paper considers the problem of multiple human target tracking in a sequence of video data. A solution is proposed which is able to deal with the challenges of a varying number of targets, interactions and when every target gives rise to multiple measurements. The developed novel algorithm comprises variational Bayesian clustering combined with a social force model, integrated within a particle filter with an enhanced prediction step. It performs measurement-to-target association by automatically detecting the measurement relevance. The performance of the developed algorithm is evaluated over several sequences from publicly available data sets: AV16.3, CAVIAR and PETS2006, which demonstrates that the proposed algorithm successfully initializes and tracks a variable number of targets in the presence of complex occlusions. A comparison with state-of-the-art techniques due to Khan et al., Laet et al. and Czyz et al. shows improved tracking performance.


I. INTRODUCTION
A. Motivation T HIS paper presents a robust multiple tracking algorithm which can be used to track a variable number of people moving in a room or enclosed environment. Hence, the estimation of the unknown number of targets and their states is considered in each video frame. Multi-target tracking has many applications, such as surveillance, intelligent transportation, behavior analysis and human computer interfaces [1]- [3]. It is a challenging problem and a wealth of research has been undertaken to provide more efficient solutions [1], [4], [5]. Two of the major challenges associated with visual tracking of multiple targets are: 1) mitigating occlusions and 2) handling a variable number of targets. Therefore, the focus of this paper is to present a technique to robustly handle these challenges.
Overcoming occlusions is a difficult task because during an occlusion it is not possible to receive measurements originating from the occluded objects. Significant research effort has been made to address this problem, e.g. [1], [4], [6], [7].
Many 2-D tracking systems rely on appearance models of people. For instance [8] applies the kernel density approach while [9] is based on color histograms, gradients and texture models to track human objects. These techniques use template matching which is not robust to occlusions because the occluded parts of targets cannot be matched. They perform an exhaustive search for the desired target in the whole video frame which obviously requires much processing time. These techniques are also sensitive to illumination changes.
Occlusions can be handled with the help of efficient association of available data to the targets. Most of the existing multiple target tracking algorithms assume that the targets generate one measurement at a time [1], [10]. In the case of human tracking from video and many other tracking applications, targets may generate more than one measurement [4], [11]. To exploit multiple measurements, most algorithms add a preprocessing step which involves extracting features from the measurements [9]. This preprocessing step solves the problem to some extent but it results in information loss due to the dimensionality reduction and hence generally degrades the tracking results. Another problem with many of these algorithms is that they adopt the hard assignment technique [9], [12] wherein likelihood models are employed to calculate the probability of the measurements.
Recently, in [4] an approach relying on clustering and a joint probabilistic data association filter (JPDAF) was proposed to overcome occlusions. Rather than extracting features from the measurements, this approach groups the measurements into clusters and then assigns clusters to respective targets. This approach is attractive in the sense that it prevents information loss but it fails to provide a robust solution to mitigate the occlusion problem. This is because it only utilizes the location information of targets and hence the identity of individual targets is not maintained over the tracking period. As a result, the tracker may confuse two targets in close interaction, eventually causing a tracker switching problem. Other advances of such filtering approaches to multiple target tracking include multiple detection JPDAF (MD-JPDAF) [13] and interacting multiple model JPDAF (IMM-JPDAF) [14].
Interaction models have been proposed in the literature to mitigate the occlusion problem [1], [15]. The interaction model presented in [1] exploits a Markov random field approach to penalize the predicted state of particles which may cause occlusions. This approach works well for tracking multiple targets where two or more of them do not occupy the same space, but it does not address tracking failures caused by intertarget occlusions during target crossovers. Dynamic models can be divided broadly into macroscopic and microscopic models [16]. Macroscopic models focus on the dynamics of a collection of targets. Microscopic models deal with the dynamics of individual targets by taking into account the behavior of every single target and how they react to the movement of other targets and static obstacles.
A representative of microscopic models is the social force model [17]. The social force model can be used to predict the motion of every individual target, i.e. the driving forces which guide the target towards a certain goal and the repulsive forces from other targets and static obstacles. A modified social force model is used in [18] and in [19] for modeling the movements of people in non-observed areas.
Many multi-target tracking algorithms in the literature only consider the case when the number of targets is known and fixed [3], [12], [15], [20]. In [1] a reversible jump Markov chain Monte Carlo (RJMCMC) sampling technique is described but in the experimental results a very strong assumption is made that the targets (ants) are restricted to enter or leave from a very small region (nest site) in a video frame. Random finite set (RFS) theory techniques have been proposed for tracking multiple targets, e.g. [21], [22], in particular the target states are considered as an RFS.
An appealing approach where the number of targets is estimated is proposed in [4] and is based on the evaluated number of clusters. Non-rigid bodies such as humans can however produce multiple clusters per target. Therefore, the number of clusters does not always remain equal to the number of targets. Hence, calculating the number of targets on the basis of the number of clusters can be inaccurate in the case of human target tracking.

C. Main Contributions
The main contributions of this work are: 1. An algorithm that provides a solution to multiple target tracking and complex inter-target interactions and occlusions. Dealing with occlusions is based on social forces between targets. 2. An improved data association technique which clusters the measurements and then uses their locations and features for accurate target identification. 3. A new technique based on the estimated positions of targets, size and location of the clusters, geometry of the monitored area, along with a death and birth concept to handle robustly the variable number of targets. The framework proposed in [4] uses the JPDAF together with the variational Bayesian clustering technique for data association and a simple dynamic model for state transition. Our proposed technique differs from that of [4] in three main aspects: 1) in the dynamic model which is based on social forces between targets, 2) a new features based data association technique, and 3) a new technique for estimation of number of targets which is based on the information from clusters, position of existing targets and geometry of monitored area location.
The remainder of the paper is organized as follows: Section II describes the proposed algorithm, Section III presents the proposed data association algorithm, Section IV explains the clustering process and Section V describes the proposed technique for estimating the number of targets. Extensive experimental validation is given in Section VI. Section VII contains discussion and finally, conclusions are presented in Section VIII.

A. Bayesian estimation
The goal of the multi-target tracking process is to track the state of N unknown targets. The state of each target is represented as x i k , (i = 1, . . . , N ), which is a column vector containing the position and velocity information of the i th target at time k. The joint state of all the targets is constructed as the concatenation of the individual target states Within the Bayesian framework, the tracking problem consists of estimating the belief of the state vector x k given the measurement vector y 1:k . The objective is to sequentially estimate the posterior probability distribution p(x k |y 1:k ) at every time step k. The posterior state distribution p(x k |y 1:k ) can be estimated in two steps. First, the prediction step [23] is where p(x k |x k−1 ) is the state transition probability. After the arrival of the latest measurements, the update step becomes p(x k |y 1:k ) = p(y k |x k )p(x k |y 1:k−1 ) p(y k |y 1:k−1 ) , where p(y k |x k ) is the measurement likelihood function. The two step tracking process yields a closed form expression only for linear and Gaussian dynamic models [10]. Particle filtering [23] is a popular solution for suboptimal estimation of the posterior distribution p(x k |y 1:k ), especially in the nonlinear and non-Gaussian case. The posterior state distribution is estimated by a set of random samples x s k , s = 1, . . . , N s and their associated weights w s k at time k where N s is the total number of particles and δ(·) denotes a multivariate Dirac delta function. We apply particle filtering within the JPDAF framework to estimate the states of targets. The JPDAF recursively updates the marginal posterior distribution p(x i k |y 1:k ) for each target. In our work, given the state vector x i k−1 of target i, the next set of particles at time step k is predicted by the social force dynamic model (equations (6)-(9)) which is described later in Section II-C. The state transition model is represented by equation (9). The N s particles and their corresponding weights can approximate the marginal posterior at time step k.
To assign weights to each sample we follow the likelihood model along with the data association framework which is explained in Section III.

B. Measurement Model
The algorithm considers L number of pixels in a single video frame captured by one camera as input measurements. From the sequentially arriving video frames, the silhouette of the moving targets (foreground) is obtained by background subtraction [24]. At time k, the measurement vector is: The state x k at time k is estimated by using the measurement vector where e m,k is the measurement noise vector. After the background subtraction [24], the measurements at time k are: which are assumed to correspond to the moving targets (please see Section VI for details). The measurement function h k (·) is represented by some features of the video frames, e.g. the color histogram. The measurement vectorỹ k contains only a set of foreground pixels. After background subtraction the number of pixels with non-zero intensity values is reduced to M , where M << L. The number of foreground pixels M is variable over time.
The background subtracted silhouette regions are clustered with the variational Bayesian algorithm described in Section IV. During the clustering process each data point (pixel)ỹ j k contains only the coordinates of the pixel. During the data association stage we use the red, green, blue (RGB) color information contained in each pixel to extract color features of a cluster.

C. The Social Force Model
The system dynamics model is defined by the following equation x where ξ k−1 is the system noise vector. Under the Markovian assumption our state transition model for the i th target is represented as where F i is the social force being applied on target i.
The motion of every target can be influenced by a number of factors, e.g. by the presence of other targets and static obstacles. The interactions between them are modeled here by means of a social force model [17], the key idea of which is to calculate attraction and repulsion forces which represent the social behavior of humans. The so-called "social forces" depend on the distances between the targets.
Following Newtonian dynamics, a change of states stems from the existence of exterior forces. Given a target i, the total number of targets that influence target i at time k, is N i . The overall force F i applied on target i is the sum of the forces exerted by all the neighboring targets, F i = j∈Ni F j . Most of the existing force based models [17] consider only repulsive forces between targets to avoid collisions. However, in reality there can be many other types of social forces between targets due to different social behavior of targets, for instance attraction, spreading and following [25].
Suppose a target i at time step k has N i neighbors, therefore, it would have N i links with other targets l = 1 . . . N i . Only those targets are considered as neighbors of target i which are within a threshold distanced from i. We consider that a social behavior over link l is given by ϕ l ∈ {1 . . . Γ} where Γ is the total number of social behaviors. An interaction mode vector ϕ i = [ϕ 1 , ϕ 2 , . . . , ϕ Ni ] for target i is a vector comprising the different social behaviors over all the links. The total number of interaction modes is S = Γ Ni .
A force due to one interaction mode is calculated as a sum of forces over all the social links where F(ϕ i ) is the force due to social mode ϕ i and φ l is the force over the link l. In our work we consider three social forces: repulsion, attraction and non-interaction. These broadly encompass all possible human motions. People naturally tend to avoid collisions between each other. This is reflected by repulsion forces. Attraction accounts for the behaviors when people approach each other for the intention of meeting, this behavior is usually ignored in the existing force based models. The non-interactions mode represents the behavior of independent motion of every person.
The repulsive force applied by target j on target i over the link l is calculated as [25] φ l − = f r exp when d l ij is less than an empirically found thresholdd, defined as 3m in Table I in Section VI-E; where b which is set equal tod is the boundary in which a target j has its influence of force on target i. The Euclidean distance between targets i and j over link l is defined as d l ij , r l ij = r l i + r l j is the sum of radii of influence of targets i and j over link l and f r is the magnitude of the repulsive force. The unit vector from target j to target i is represented as u ji .
When targets approach each other, repulsive forces act, to avoid collisions. The repulsive force therefore increases as described by equation (7). Similarly, when targets move apart the repulsion decreases because it is less likely for a collision to occur, and therefore they tend to come closer to each other. Attraction is a phenomenon opposite to repulsion. When targets move apart, they tend to attract each other for the intension of meeting and when they move closer, the repulsion increases which means they are less attracted. This attraction phenomenon is described by equation (8). However, if targets are at a distance more than a threshold, it can be assumed they do not influence each other.
An attractive force applied by target j on target i over the link l is calculated as [25] φ l when d l ij is less than an empirically found thresholdd, defined as 3m in Table I in Section VI-E, where f a is the magnitude of the attractive force and force due to non-interactions is considered to be zero.
To represent the state x k of targets at every time step k, a pixel coordinate system is used. Since, our state transition equation (9) below is based on the social force model, to predict the new particles x s k , position and velocity values contained within the state vector x k−1 and the particles at time step k−1 are converted to the image ground coordinate system (a process known as homography [26]). After predicting the new set of N s particles for state x k these particles are converted back to pixel coordinates. For conversion between the two coordinate systems the four point homography technique is implemented, as used in [27]. The position and velocity at time step k in pixel coordinates are represented as p k and v k , respectively, whereas, the position and velocity in ground coordinates systems are represented asp k andṽ k respectively. At every time step a new state is predicted with respect to interaction mode ϕ i according to the following model [25] where ∆t is the time interval between two frames, m represents mass of the target and ξ k is the system noise vector. The social force model gives more accurate predictions than the commonly used constant velocity model [28] and helps tracking objects in the presence of close interaction. This cannot be done with a simple constant velocity motion model.
The implementation of the force model is further explained in Algorithm 1 at the end of the next section.

III. DATA ASSOCIATION AND LIKELIHOOD
To deal with the data association problem we represent a measurement to target (M → T ) association hypothesis by the parameter ψ k , whereas Ψ k represents the set of all hypotheses. At this stage we assume that one target can generate one measurement. However, later in this section we relax this assumption and develop a data association technique for associating multiple measurements (clusters of measurements) to one target. If we assume that measurements are independent of each other then the likelihood model conditioned on the state x k of targets and measurement to target association hypothesis ψ k can be represented as in [10] (equation (11)) where J o and J are, respectively, measurement indices sets, corresponding to the clutter measurements and measurements from the targets to be tracked, for convenance of notation the time dependence is not shown on these quantities. The terms p c (·) and p t (·) are characterizing, respectively, the clutter likelihood and the measurement likelihood association to each target i whose state is represented as x i,j k . If we assume that the clutter likelihood p c (·) is uniform over the volume γ of the measurement space and f k represents the total number of clutter measurements, then the probability p(y k |x k , ψ k ) can be represented as The form of p t (y j k |x i,j k ) is given in Section III-C.

A. The JPDAF Framework
The standard JPDAF algorithm estimates the marginal distribution of each target by following the prediction and update process described in equations (1) and (2). In the development below, we exploit the framework from [10]. The prediction step for each target independently becomes where we assume that the prediction at time step k is based on the measurements only at k − 1. The JPDAF defines the likelihood of measurements with respect to target i as in [10] (equation (25)) where A i,j k is the association probability that the i th target is associated with the j th measurement and A i,o k is the association probability that the target i remains undetected; for notational convenience, in equation (13) and the remainder of the paper we drop the subscript t on the target likelihood.With prediction and likelihood functions defined by equations (12) and (13), respectively, the JPDAF estimates the posterior probability of the state of target i p(x i k |y 1:k ) ∝ p(y k |x i k )p(x i k |y 1:k−1 ).
The association probability A i,j k can be defined as in [10] (equation (27)) where Ψ i,j k represents the set of all those hypotheses which associate the j th measurement to the i th target. The probability p(ψ k |y 1:k ) can be represented as equation (28) in [10] p(ψ k |y 1: where p(ψ k |x k ) is the probability of assignment ψ k given the current state of the objects. We assume that all assignments have equal prior probability and hence p(ψ k |x k ) can be approximated by a constant. If we define a normalizing constant π ensuring that p(ψ k |y 1:k ) integrates to one and define the predictive likelihood p(y j k |y 1:k−1 ) as then equation (16) becomes

B. Particle Filtering Based Approach
The standard JPDAF assumes a Gaussian marginal filtering distribution of the individual targets [10], [28]. In this paper we adopt a particle filtering based approach which represents the state of every target with the help of samples. Therefore, we modify equation (18) as where N s is the total number of particles which are required to estimate the filtering distribution over every individual target andŵ i,s k represents the associated predictive weights as described in [10]. We can modify equation (19) as Once the probability p(ψ k |y 1:k ) is calculated according to (20), we can substitute it in equation (15) to calculate the association probability. We can further use this association probability to calculate the measurement likelihood from equation (13). After drawing particles from a suitably constructed proposal distribution [23], the weights associated with particles are calculated,ŵ i,s k = p(y k |x i k ). The particles and their weights can then be used to approximate the posterior distribution for individual targets where x i,s k is the s th particle for target i and w i,s k is the associated weight. For estimation of a state we predict N ϕ particles with respect to every interaction mode. There are S interaction modes, therefore, N s = S × N ϕ represents the total number of samples.

C. Algorithm to Associate Multiple Measurements to a Target
In the case of tracking people in video, every person generates multiple measurements. To avoid any information loss we relax the assumption that every person can generate a single measurement. Here we present a data association technique for associating multiple measurements to every target. To achieve this, the algorithm groups the foreground pixels into clusters with the help of a variational Bayesian clustering technique, and then instead of associating one measurement to one target we associate clusters to every target. The clustering process returns κ clusters, where the number of clusters is not predefined or fixed. Clusters at discrete time k are regions represented as Z 1 k , Z q k , . . . , Z κ k , where Z q k is the q th cluster which contains certain measurements, i.e. a number of vectorsỹ j k of foreground pixels. The complete clustering process is described in Section IV. The clustering process aims to assign measurements to the respective clusters and gives as output a correspondence matrix indicates, at discrete time k, to which cluster, the measurement vectorỹ j k corresponds. All but one of the elements of b j k is zero, if e.g.ỹ j k belongs to the q th cluster then the correspondence vector will be, b j k = [0, . . . , 0, 1, 0, . . . , 0], which shows only the q th element of b j k is nonzero. Please note that from this section onwards all equations are written with respect to the measurement vector y k that contains foreground pixels.
We modify the measurement to target (M → T ) association probability A i,j k to calculate the cluster to target (Z → T ) association probability A i,q k , where q represents the cluster index where Ψ i,q k represents the set of all those hypotheses which associate the q th cluster to the i th target. We modify equation (20) to define the probability p(ψ k |ỹ 1:k ) as whereỹ j k : b j,q k = 1 represents only those measurements which are associated with cluster q. The measurement likelihood defined in equation (13) is modified as follows where q : b j,q k = 1 represents that the j th measurement is associated with the q th cluster and A i,o k represents the probability that the cluster is not associated to the i th target and A i,(q:b j,q k =1) k is the probability that the cluster is associated to the i th target. To limit the number of hypotheses when the number of targets increases, we have adopted Murty's algorithm [29]. This algorithm returns the k-best hypotheses. The elements of the cost matrix for Murty's algorithm are calculated with the particles x i,s where cost i q represents the cost of assigning the q th cluster to the i th target.
The JPDAF [30] data association technique proposed in [4] performs well when two targets do not undergo a full occlusion. However, this data association technique sometimes fails to associate measurements correctly, especially when two or more targets partially occlude each other or when targets come out of full occlusion. This is due to the fact that this data association is performed on the basis of location information without any information about the features of targets. This can result in missed detections and wrong identifications.
In order to cope with these challenges under occlusions, a data association technique is proposed which assigns clusters to targets by exploiting both location and features of clusters. First, we define the term ỹ j where µ s = x i,s k is the mean vector and Σ s is the fixed and known diagonal covariance matrix. The association probability p(Z q k |x i,s k ) in equation (26) is calculated by extracting features from the clusters at time step k and comparing them with a reference feature (template) of target i. The reference feature for target i is extracted when it first appears in the monitored area. To improve the computational efficiency, the probability p(Z q k |x i,s k ) is calculated only when targets approach each other. This association probability is defined as where σ 2 is the measurement noise variance, H i ref and H q k are histograms of features of the i th target and q th cluster, respectively and is the distance between them. One possible option is to use the Bhattacharyya distance [31]. The parameter ϑ is a predefined threshold, and the probability of occlusion is defined as where d i min is the minimum Euclidean distance among a set of distances which is calculated for target i from all other targets at time k, and ∆ c is the distance constant which limits the influence of target n on target i. A higher value of minimum distance d i min results in a smaller value of probability of occlusions and when this probability drops below threshold ϑ the probability p(Z q k |x i,s k ) becomes unity.

IV. VARIATIONAL BAYESIAN CLUSTERING -FROM MEASUREMENTS TO CLUSTERS
The clustering process aims to subdivide the foreground image into regions corresponding to the different targets. A variational Bayesian clustering technique is developed to group the M foreground pixels into κ clusters. Each cluster at time index k is represented by its center µ q k , where q = 1, . . . , κ. A binary indicator variable b j,q k ∈ {0, 1} represents to which of the q th clusters data pointỹ j k is assigned; for example, when the data pointỹ j k is assigned to the q th cluster then b j,q k = 1, and b j,l k = 0 for l = q. We assume that the foreground pixel locations are modeled by a mixture of Gaussian distributions. The clustering task can be viewed as fitting mixtures of Gaussian distributions [32] to the foreground measurements. Every cluster Z q k of foreground pixel locations is assumed to be modeled by a Gaussian distribution, N (ỹ j k |µ q k , Σ q k ) which is one component of the mixture of Gaussian distributions. Each cluster has a mean Algorithm 1 The social force model Input: N s particles and the state of targets at time step k − 1 Output: N s particles for time step k 1: Convert the state at k − 1 from pixel coordinates system to real world ground coordinates. 2: Use ground coordinates to calculate distances between targets. 3: Based on distances create links between targets. 4: for i = 1, ..., N (where N is the number of targets) do 5: Create S = Γ Ni social modes 6: Calculate forces due to every social mode by using equations (6), (7) and (8). 7: for s = 1 : Γ Ni : N s (where N s is the number of particles) do 8: Convert the particle x i,s k−1 to the ground plane. 9: Add random noise to get x i,s k−1 + ξ k−1 . 10: for s ϕ i = 1, ..., Γ Ni do 11: Predict particle x i,s+sϕ i −1 k−1 by using particle the social mode ϕ i in equation (9). 12: to the pixel coordinates. 13: end for 14: end for 15: end for vector µ q k and a covariance matrix Σ q k . Hence, the probability of a data pointỹ j k can be represented as where k represents the probability of selecting the q th component of the mixture which is the probability of assigning the j th measurement to the q th cluster. If we assume that all the measurements are independent, then the log likelihood function becomes [32] The data pointỹ j k belongs to only one cluster at a time and this association is characterized by a correspondence variable vector b j k . Therefore, it can be represented as C q k = p(b j,q k = 1) where the mixing coefficients must satisfy the following conditions: 0 < C q k ≤ 1 and κ q=1 C q k = 1. Since only one element of b j k is equal to 1, the probability of the mixing coefficient can also be written as Similarly, we assume that p(ỹ j k |b j,q Since only one element of b j k is equal to 1, we can write again As a result of the clustering process we will obtain estimates of the probabilistic weight vector C k , the mean vectors µ k , the covariance matrices Σ k and correspondence matrix B k . The estimation problem can be simplified by introducing hidden (latent) variables. We consider the correspondence variables B k as latent variables and suppose that they are known. Then {ỹ k , B k } represents the complete dataset. The log likelihood function for this complete data set becomes ln p(ỹ k , B k |C k , µ k , Σ k ). Now the new set of unknown parameters contains B k , C k , µ k and Σ k which can be estimated by the proposed adaptive variational Bayesian clustering algorithm. If, for simplicity, these parameters are represented as a single parameter Θ k = (B k , C k , µ k , Σ k ), then the desired distribution can be represented as p(Θ k |ỹ k ). The log marginal distribution p(ỹ k ) can be decomposed as described in [32] (Section 9.4 and equations (10.2)-(10.4)) We define the distribution Q(Θ k ) as an approximation of the desired distribution. Then the objective of the variational Bayesian method is to optimize this distribution, by minimizing the Kullback Leibler (KL) divergence The symbol · · represents the KL divergence, and Minimizing the KL divergence is equivalent to maximizing the lower bound L(q). The maximum of this lower bound can be achieved when the approximate distribution Q(Θ k ) is exactly equal to the desired posterior distribution p(Θ k |ỹ k ).
The joint distribution p(Θ k ,ỹ k ) can be decomposed as [32] p We assume that the variational distribution Q(Θ k ) = Q(B k , C k , µ k , Σ k ) can be factorized between latent variable B k and parameters C k , µ k and Σ k .
A similar assumption is made in [32], please see equation (10.42). Optimization of the variational distribution can therefore be represented as By considering equation (10.54) from [32] the distribution Q * (C k , µ k , Σ k ) can further be decomposed into where (·) * represents the optimum distribution. Therefore, the optimum distribution can be written as This shows that the optimum over the joint distribution is equivalent to obtaining Q * (B k ), Q * (C k ) and Q * (µ k , Σ k ). Therefore, the optimum distributions over Θ k can be evaluated by optimizing L(q) with respect to all parameters one-by-one. A general form of optimization can be written as in equation (10.9) of [32] ln where E d =c [.] represents the mathematical expectation with respect to the distributions over all the elements in represents the optimum approximate distribution over the c th component of Θ k . We next evaluate the optimum distributions Q * (B k |C k ), Q * (C k ) and Q * (µ q k , Σ q k ) by using equation (41). 1) Optimum distribution over B k : According to equation (41), the optimum distribution over B k can be written as Probabilities p(B k |C k ) and p(ỹ k |B k , µ k , Σ k ) can be defined by using equations (32) and (33) respectively and By using equations (42), (43) and (44), the optimum distribution over the correspondence variable B k becomes where r j,q k is the responsibility that component q takes to explain the measurementỹ j k . The derivation of equation (45) and r j,q k can be found in Appendix-A. 2) Optimum distribution over C k : Before evaluating the optimum distributions Q * (C k ) and Q * (µ k , Σ k ) we need to first define their priors. The Dirichlet distribution is chosen as a prior over the mixing coefficient C k where Dir(·) denotes the Dirichlet distribution, and α • is an effective prior number of observations associated with each component of the mixture. Using equation (41) we can write where const. represents a constant. The optimum distributions Q * (C k ) over C k can be calculated by using equations (44), (46) and (47), which becomes where α k = [α 1 k , . . . , α κ k ] and one of its components α q k can be defined as and r j,q k is the responsibility that component q takes to explain the measurementỹ j k . The derivation of (48) can be found in Appendix-B.
3) Optimum distribution over µ k and Σ k : The prior over the mean µ k and the covariance matrix Σ k is defined by the independent Gaussian-Wishart distribution where m • , β • , W • and υ • are the prior parameters. By using equations (44), (46) and (47), decomposition of equation (37) and following the steps explained in Appendix-C, the distribution becomes where m q , β q , W q and υ q are defined from where,ȳ and The variational Bayesian technique operates in two steps to optimize the posterior distributions of unknown variables and parameters. In the first step it calculates the responsibilities r j,q k using equation (85) and in the second step it uses these responsibilities to optimize the distributions by using equations (45), (48) and (52). These steps are repeated until some convergence criterion is met. In our work we monitor the lower bound L(q) after every iteration to test the convergence. When the algorithm converges, the value of the lower bound does not change more than a small amount. The clustering algorithm is further summarized in Algorithm 2 given in Section V.
One of the important advantages of variational Bayesian clustering is that it can automatically determine the number of clusters by using the measurements. This can be achieved if we set the parameter α • less then 1. This helps to obtain a solution which contains a minimum number of clusters to represent the data [4].
The position and shape of clusters are defined using the parameters m • and W • , respectively. A possible choice for selecting these priors is given in Section VI. This stage returns the minimum possible number of clusters and their associated measurements which are defined by the correspondence matrix B k . In the next stage of the algorithm these clusters are assigned to the targets by using a data association technique as described after Section III-C.

V. VARIABLE NUMBER OF TARGETS
Many multi-target tracking algorithms assume that the number of targets is fixed and known [3], [12], [15], [20]. In [33] the JPDAF is extended to a variable number of targets. In other works, including [4] the posterior probability of the number of targets N is estimated given the number of clusters κ at each discrete time step k where p(κ k |N k ) is the probability of the κ k clusters given N k targets. Here we deal with changeable shapes and calculate the number of targets in a different way. The number of targets is estimated based on the: 1) location µ q k of clusters at time step k, 2) size of clusters at time step k, and 3) position vector x k−1 of the targets at the previous time step k − 1.
The number of targets in the first video frame is determined based on the variational Bayesian clustering. At every following video frame the number N k of targets is evaluated, by the following stochastic relationship, similarly to [34] The variable k ∈ {−1, 0, 1} is defined as where p(death) corresponds to the probability of decrementing the number of targets, similarly, p(birth) represents the probability of incrementing number of targets and T hr is an appropriately chosen threshold. We assume that at a given time only one target can enter or leave. In the monitored area people can only enter or leave through a specifically known region in a video frame. This known region is called the red region. This red region is modeled by a Gaussian distribution with its center µ r and covariance matrix Σ r . Multiple Gaussian distributions can be used if the red region was disjoint.

A. Probability of death
The probability of death of a target depends on two things: 1) the probability p(Z r k ) of finding a cluster Z r k in the red region and 2) the probability p(Existing T arget|Z r k , x k−1 ) that the cluster found in the red region is due to an existing target (denoted by the variable Existing T arget). This means that when a target is leaving the monitored area, a cluster is found in the red region, i.e. at the entrance or exit point, and this cluster will be due to an existing target. Therefore, probability of death is modeled as p(death) = p(Z r k )p(Existing T arget|Z r k , x k−1 ). (62) In order to estimate p(Z r k ), we calculate the probability of the measurements (foreground pixels) from Note that we are not considering measurements of one cluster only, because a target in the red region may generate more than one cluster. Therefore, measurements which have probability p(ỹ j k |µ r , Σ r ) greater than a set threshold value are considered as measurements found in the red region and all these measurements are therefore grouped to form cluster Z r k . Since we assume that maximum one target can be in the red region, then, all the measurements of Z r k are considered to be generated by one target. For an N p number of measurements in cluster Z r k , the probability of finding a cluster in the red region depends on the number of pixels in Z r k and is calculated as where ∆ p is a pixel constant which limits the range of pixels to be considered for calculating p(Z r k ). The second probability in equation (62) is the probability that cluster Z r k is generated by an existing target. This probability depends on the location of existing targets and the location of cluster Z r k . The probability p(Existing T arget|Z r k , x k−1 ) is calculated as where ∆ d is a constant which can be chosen experimentally.
To calculate the distanced min , we chose the minimum distance among the distances which are found from the centroid of Z r k , centroid(Z r k ) to each of the existing targets at time step k − 1. Note that thisd min is different from the d min calculated in equation (28).
where · denotes the Euclidean distance operator. Note that the probability Z r k depends on the distance between the centroid of Z r k and the closest existing target. The probability that cluster Z r k is generated by an existing target will be high ifd min is small and vice versa.

B. Probability of Birth
The probability of birth of a target depends on two different probabilities: 1) the probability of finding a cluster Z r k in the red region p(Z r k ) and 2) the probability that the cluster found in the red region is not due to an existing target. According to our assumption that at a given time step the cluster in the red region can only be either from an existing target or a new target, the probability that the cluster in the red region is due to a new target is 1 − p(Existing T arget|Z r k , x k ). The probability of birth can be calculated as p(birth) = p(Z r k )(1 − p(Existing T arget|Z r k , x k )). (67) Finally, (62) and (67) are applied in (60) and the number of targets is updated by using (66). This completes the full description of our algorithm. A summary of the overall proposed algorithm is described in Algorithm 2.

VI. EXPERIMENTAL RESULTS
To examine the robustness of the algorithm to close interactions, occlusions and varying number of targets, the algorithm is evaluated by tracking a variable number of people in three different publicly available video datasets: CAVIAR, PETS2006 and AV16.3. The test sequences are recorded at a resolution of 288 × 360 pixels at 25 frames/sec and in total there are 45 sequences. The performance of the proposed algorithm is compared with recently proposed techniques [4], [1] and [35]. All the parameters are discussed in the following subsections. The algorithm automatically detects and initializes the targets when they enter the monitored area. For evaluation the tracking results we convert the final tracked locations of targets from pixel coordinates to the ground coordinates by using four point homography.

A. Background Subtraction Results
The codebook background subtraction method [24] is implemented which is one of the best background subtraction methods since it is resistant to illumination changes and can capture the structure of background motion. Adaptation to the background changes and noise reduction is additionally achieved with the blob technique [36]. Figs. 1, 2 and 3 show the results obtained from the codebook background subtraction method [24] for a few selected video frames from the AV16.3, CAVIAR and PETS2006 datasets respectively. In the experiment, we set the shadow bound α = 0.5, highlight bound β = 2 and the color detection threshold ε = 20 (see [24] for further details about these parameters). These parameters are the same for all the sequences of all three datasets. The background subtraction results provide the coordinates of the foreground pixels (the moving objects) which represent data pointsỹ k . Frames 440 and 476 in Fig. 2 show that we get a few extra foreground pixels due to reflections on the floor and in the glass which can be eliminated at the clustering stage. Background subtraction results for certain frames of sequence "seq45-3p-1111 cam3 divx audio" of the AV16.3 dataset: codebook background subtraction is implemented to separate the foreground pixels from the background.

B. Clustering results
For the clustering process at each time step we assume that there are three types of clusters: 1) clusters for the existing Evaluate initial ρ j,q • by using initialized parameters α • , β • , m • W • and υ • in equations (81)  for q = 1, ..., κ do 15: Evaluate new N q k ,ȳ q k , α q k , β q , m q , W −1 q , υ q and S q k with equations (50), (55), (49), (53), (54), (56), (57) and (58) respectively. For first iteration use initial responsibilities r j,q • . 16: end for 17: Evaluate new responsibilities r j,q k for all j and q by using new N q k ,ȳ q k , α q k , β q , m q , W −1 q , υ q and S q k and repeating steps (4) through (12). 18: end while 19: Assign the l th cluster to measurementỹ j k , when r j,l k = max q=1:κ r j,q k and repeat it for all the measurements. 20: Delete the small clusters.

Data Association and Tracking
23: Evaluate the cost matrix by using equation (25). 24: Evaluate k-best hypotheses by using Murty's algorithm [29]. 25: for i = 1, ...., N (where N is the total number of targets) do 26: Draw N s samples by using the state transition equation (9) as explained in Algorithm 1.

27:
for q = 1, ..., κ do 28: Evaluate the set of hypotheses Ψ i,q k which assign the q th cluster to the i th target. 29: For every hypothesis ψ k ∈ Ψ i,q k evaluate ln p(ψ k |ỹ 1:k ) by using equations (23) and (26). 30: Evaluate A i,q k with equation (22). 31: end for 32: Weight all the particles using equation (24). 33: Update states of all the targets according to equation (21). 34: end for Background subtraction results for certain frames of sequence "S1-T1-C" of the PETS2006 dataset: codebook background subtraction is implemented to separate the foreground pixels from the background. targets, 2) clusters in the neighbors of the existing targets and 3) clusters near boundaries of the field of view. Due to the Dirichlet prior over the mixing coefficient and by setting α • = 0.6 the clustering algorithm converges automatically to the minimum possible number of clusters. The means m • of clusters are initialized with the tracked location of existing targets and hypothesized location of new targets. Other prior parameters are defined as: υ • = 3 and β • = 1. However, a full study about sensitivity of the system to the choice of different parameters is beyond the scope of the paper.
The prior parameter W • used in equation (56) determines the human shape which is modeled as an ellipse. It is defined with the help of the following equation where l1 and l2 are equatorial radii of the ellipse which models the human shape. The equatorial radii l1 and l2 are set to 500 and 300 respectively while U is defined as U = cos(π/2) − sin(π/2) − sin(π/2) cos(π/2) .
The clustering results for a few of the video frames are shown in Figs. 4, 5 and 6. Blue, red, green, magenta, cyan, yellow and black represent first, second, third, fourth, fifth, sixth and seventh clusters respectively. If there are more than 7 clusters we repeat the color scheme.  To eliminate the extra foreground pixels due to the reflections, the small clusters consisting of less than 100 pixels are deleted.

C. Data Association and Occlusion Handling
For data association the ten best hypotheses of joint association event ψ k are considered which are obtained with the help of Murty's algorithm [29]. The Bhattacharyya distance between the color histograms of cluster q and target i is used to calculate the distance where H i ref is the reference histogram which is created by using the cluster associated with target i at the time step when the target i first appears in the video. H q k is the histogram created for cluster q at the current time step k and ρ( where G represents the number of histogram bins and we have used 16 × 16 × 16 color histograms bins. A tracking algorithm with only a JPDAF, without the variational Bayesian clustering and without a social force model fails to identify an object during close interactions. In the proposed algorithm these tracking failures are overcome with the help of the proposed data association technique. Sometimes, due to the varying shape of the targets, the clustering stage may produce more than one cluster per target. Therefore, the proposed data association technique assigns multiple clusters to every target with some association probability.

D. Variable Number of Targets Results
The robustness of the algorithm for estimating the correct number of targets in a video sequence is compared with the framework developed in [4]. In [4]

E. Tracking Results
A minimum mean square error (MMSE) particle filter is used to estimate the states of the multiple targets. When a new target enters the room, the algorithm automatically initializes it with the help of the data association results by using the mean value of the cluster assigned to that target. Similarly, the algorithm removes the target which leaves the room. The particle size N s is chosen to be equal to 60 and the number S of social links are updated at every time step. Table I shows the values of the other parameters used in the social force model. The prediction step is performed by the social force model, described in Section II-A. The tracking results for a few selected frames from AV16.3, CAVIAR and PETS2006 datasets are given in Figs. 9, 10 and 11 respectively. Blue, red, green, magenta and cyan ellipses represent first, second, third, fourth and fifth targets, respectively. Fig. 9 shows that the algorithm has successfully initialized new targets in frames 225 and 278. In frame 278 it can be seen that the algorithm can cope with the occlusions between target 1 and target 2. Frame 288 shows that the tracker keeps tracking all targets even when they are very close to each other. Frames 320 and 326 show that the algorithm has successfully handled the occlusion between targets 2 and 3. In frame 375 we can see that target 2 has left the room and the algorithm has automatically removed its tracker, which is started again when the target has returned back in frame 420. Success in dealing with occlusions between targets 3 and 1 can be seen in frames 389 and 405.
Results in Fig. 10 demonstrate that new target appearing in frame 334 is initialized. In frame 476 it can be seen that the tracker has solved the occlusion between target three and four. Similarly, frames 524 and 572 show that the tracker copes with occlusion between target three and five. A similar behavior is also observed in Fig. 11.
1) CLEAR MOT matrices: The tracking performance of the proposed algorithm is measured based on two different performance measures: CLEAR multiple object tracking (MOT) matrices [37] and OSPAMT [38]. The performance is compared with the methods proposed in [4], [1] and [35].
The CLEAR MOT measure includes the multiple object tracking precision (MOTP) matrix, miss detection rate, false positive rate and mismatch rate.
Let  The threshold is set at 45cm, which is the distance between the estimated position of a target and its ground truth beyond which it is considered as a missed detection.
If we assume that m k , mme k , f p k and gt k are respectively the total number of missed detections, mismatches, false positives and ground truths at time k, the errors are calculated as [37] The precision is calculated as [37] where d i,k is the distance between the estimated location and the ground truth location of target i and c k is the total number of matches found at time step k. Table II presents the performance result of the proposed algorithm and the tracking algorithms proposed in [4], [1] and [35]. The performance results are obtained by using 45 sequences from AV16.3, CAVIAR and PETS2006 datasets. All three datasets have different environments and backgrounds. The results show that the proposed algorithm has significantly improved performance compared with other techniques [4], [1] and [35]. Table II shows that there is a significant reduction in missed detections, mismatches and false positives. The reduction in the missed detections is mainly due to the proposed clustering based approach along with the new social force model based estimation technique. A simple particle filter based algorithm without any interaction model proposed by [35] shows the worst results mainly due to the highest number of mismatches. The algorithm proposed in [1] shows better results because of a better interaction model. However, the number of missed detections, mismatches and false positives are higher than that of [4]. This is because [1] does not propose a solution when two targets partially occlude each other or appear again after full occlusion. The technique proposed in [4] gives improved results due to the clustering and the JPDAF however fails when targets reappear after occlusion. This is due to the fact that only location information is used for data association and a There is a significant 3.26% reduction in the wrong identifications which has improved the overall accuracy of the tracker. Average precision results for the video sequences from AV16.3, CAVIAR and PETS2006 data sets are shown in Table III. Precision plots for two video sequences against the video frame are shown in Fig. 12. Results show that the MOTP remains less than 4cm for most of the frames. It increases to 12cm for the frames where occlusions occur but it again drops when targets emerge from occlusion.  2) OSPAMT: A new performance metric, optimal subpattern assignment metric for multiple tracks (OSPAMT) [38] is also used to evaluate the performance of the proposed technique and approach of [4]. The OSPAMT metric calculates the localization distance D Loc k (ω, ω ) between a set of true tracks ω and a set of estimated tracks ω at time step k as follows where n t = max{n ω t , n ω t } and n ω t and n ω t are the number of targets in true and estimated set of tracks, respectively, at time step k. The distance d t (τ ω i , ω ) is between the i th true track and the set of estimated tracks ω , complete explanation of which can be found in Section IV of [38]. The distance D Card k (ω, ω ) at time index k is calculated as where S t is defined as follows wheren λ k,i is the number of targets at time k in ω assigned to target i in ω, λ represents an assignment between tracks in ω and the tracks in ω, ∆ is the assignment parameter, c is the cutoff parameter and p = 2. Details about all these parameters and the OSPAMT metric can be found in [38]. For results in Figures 13 and 14 we have used c = 80 and ∆ = 10. Note that D Card k (ω, ω ) is basically a distance between the objects. Figures 13 and 14 show the OSPAMT distances, (75) and (76), respectively, plotted against the video frame index. Figure  13 shows low localization errors compared to the technique proposed in [4]. Figure 14 indicates a peak dynamic error when a new target enters the field of view which is between frame numbers 63 and 65.  OSPAMT distance plot against the video frames of sequence "ThreePastShop2cor" of the CAVIAR dataset.

F. Computational complexity
A variable measurement size technique reduces significantly the computational complexity. The measurement size is increased or decreased on the basis of the distance between the targets. Downsampling of foreground pixels is performed when the distance between targets becomes 80cm or less. The decrease in the measurement size reduces the time and number of iterations for reaching clustering convergence which results in improved computational complexity. This is shown in Table  IV with the help of a few selected frames from sequence "seq45-3p-1111 cam3 divx audio" of the AV16.3 dataset.
Table IV also demonstrates that the proposed technique reduces the number of iterations needed for achieving convergence. The reduction in the number of iterations for convergence improves the run-time of the tracking algorithm. The average run-time (calculated using 45 video sequences from three datasets) of the proposed algorithm due to the reduction in number of iterations is 0.587 seconds per frame, as compared to the run-time without measurement reduction which is 2.611 seconds per frame. The run-time of the approach of The results in Tables II and III confirm that the overall performance of the proposed tracking algorithm results in an improvement as compared with the recently proposed techniques [4], [1] and [35].

VII. DISCUSSION
The quality of the clustering results influences the proposed data association technique. We initialize the clusters on the basis of the estimated locations of targets and by taking into account merging and splitting of targets. The motion of people is predicted with the social force model. This yields accurate clustering results in scenarios where one cluster contains regions of multiple targets. Therefore, the proposed data association technique performs well even during close interactions of targets. It associates multiple clusters to targets with the correct proportion of probability. These association results help in achieving highly accurate tracking of multiple targets even during their close interactions and partial occlusions. The extensive evaluation of the proposed technique on 45 sequences from three very different datasets validates the tracking accuracy of the proposed algorithm in such scenarios.
A further improvement in accuracy of cluster-to-target association can be achieved by calculating joint associations between multiple clusters. A several-to-one or several-toseveral correspondence strategy can be adopted. However, for these strategies the association problem is combinatorial and complexity may become intractable because it will grow geometrically. The cluster-to-target association can also be improved by considering the regions of multiple targets in a cluster as unique identities and associating those regions to respective targets [27].

VIII. CONCLUSIONS
A learned variational Bayesian clustering and social force model algorithm has been presented for multi-target tracking by using multiple measurements originating from a target. An improved data association technique is developed that exploits the clustering information to solve complex inter-target occlusions. The algorithm accurately and robustly handles a variable number of targets. Its is compared extensively with recent techniques [1], [4] and [35].

A. Derivation of equation (45)
The presented derivation is for time index k and for simplicity we omit the index k from the equations. According to (42) we can write ln Q * (B) = E C,µ,Σ [ln p(ỹ, B, C, µ, Σ)] + Const. (77) and using the decomposition of equation (37) The optimal distribution over B is calculated by considering different forms of q(B) and by keeping the remaining terms constant. Therefore, quantities not depending on B will merge into the constant and we can write equation (78) as By using equations (32) and (33) we can write Considering that (ρ j,q ) b j,q + Const.
where according to [32] whereα = κ q=1 α q andψ(·) is the digamma function. We can write Requiring that this distribution be normalized, r j,q = ρ j,q κ q=1 ρ j,q , B. Derivation of equation (48) Using again equation (41) We know that the right-hand side of this expression decomposes into a sum of terms involving only C and the terms only involving µ and Σ which implies that p(C, µ, Σ) factorizes to give p(C)p(µ, Σ).

C. Derivation of equation (52)
Consider again the factorization of p(C, µ, Σ) to give p(C), p(µ, Σ). By keeping only those terms of equation (91) which depend on µ and Σ we get ln Q * (µ, Σ) = whereas the Wishart distribution can be defined as also N (µ q |m • , β −1 • Σ q ) = ln Using these equations in (107), keeping only the terms which contain µ and Σ we get the final expression for (52) where β q , m q , W q and υ q are defined in (53)-(57).