Articulated Body Motion Capture by Annealed Particle Filtering
Jonathan Deutscher
University of Oxford
Dept. of Engineering Science
Oxford, OX13PJ
United Kingdom
jdeutsch@robots.ox.ac.uk
Andrew Blake
Microsoft Research
1 Guildhall St,
Cambridge, CB2 3NH
United Kingdom
ablake@microsoft.com
Ian Reid
University of Oxford
Dept. of Engineering Science
Oxford, OX13PJ
United Kingdom
ian@robots.ox.ac.uk
Abstract
The main challenge in articulated body motion track-
ing is the large number of degrees of freedom (around 30)
to be recovered. Search algorithms, either deterministic or
stochastic, that search such a space without constraint, fall
foul of exponentialcomputational complexity. One approach
is to introduce constraints — either labelling using markers
or colour coding, prior assumptions about motion trajec-
tories or view restrictions. Another is to relax constraints
arising from articulation, and track limbs as if their mo-
tions were independent. In contrast, here we aim for general
tracking without special preparation of subjects or restric-
tive assumptions.
The principal contribution of this paper is the develop-
ment of a modified particle filter for search in high dimen-
sional configuration spaces. It uses a continuation princi-
ple, based on annealing, to introduce the influence of narrow
peaks in the fitness function, gradually. The new algorithm,
termed annealed particle filtering, is shown to be capable of
recovering full articulated body motion efficiently.
1. Introduction
Marker-based human motion capture has been used commer-
cially [19] for a number of years with applications found
in special effects and biometrics. The use of markers how-
ever is intrusive, necessitates the use of expensivespecialised
hardware and can only be used on footage taken especially
for that purpose. A markerless system of human motion cap-
ture could be run usingconventionalcamerasand withoutthe
use of special apparel or other equipment. Combined with
today’s powerful off-the-shelf PC’s, cost-effective and real-
time markerless human motion capture has for the first time
become a possibility. Such a system would have a greater
number of applications than its marker based predecessor
ranging from intelligent surveillance to character animation
and computer interfacing. For this reason the field of human
motion capture has recently seen somewhat of a renaissance.
Research into human motion capture has so far failed to
produce a full-body tracker general enough to handle real-
istic real-world applications. This gives an insight into the
difficulty of the problem. Research has concentrated on the
articulated-model based approach. The reason this approach
is popular is the high level output it produces in the form of
a model configuration for each frame. This output can easily
be used by higher-order processes to perform tasks such as
character animation.
The problem with using articulated models is the high
dimensionality of the configuration space and the exponen-
tially increasing computational cost that results. A realistic
articulated model (see figure 4) of the human body usually
has at least 25 DOF. The model used in this paper for ex-
ample has 29 DOF, and models employed for commercial
character animation usually have over 40.
A number of effective 2D systems have been presented
[7] [10]. These are good for applications such as surveil-
lance, however they do not provide output in the form of 3D
model configurations that are needed for applications such
as 3D character animation.
There are several possible strategies for reducing the di-
mensionality of the configuration space. Firstly it is possible
to restrict the range of movement of the subject. This ap-
proach has been pursued by Hogg [8], Rohr [17] and Niyogi
[15]. All three assume the subject is walking. Rohr even
reduces the dimension of the problem to the phase of the
walking cycle. Goncalves [6] and Deutscher [3] assume a
constant angle of view of the subject as does Bregler [2] and
Rehg [16]. Such an approach greatly restricts the resulting
trackers generality.
Another way to constrain the configuration space is to
perform a hierarchical search. If one part of an articulated
model can be localised independently then it can be used as
a constraint for reducing the rest of the model. Gavrila [4]
does just this when he uses what he terms search space de-
composition. He is able to localise the torso using colour
cues and uses this information to constrain the search for
the limbs. Without the assistance of colour cues (or other
labelling cues) however it is very hard to independently
localise specific body parts in realistic scenarios. This is
mainly due to the problem of self occlusion and rules out
the use of a hierarchical search.
For a practical full body tracker to be developed it can-
not rely on assumptions about motion, angle of view or the
availability of labelling cues. The principal contribution of
this paper is the development of a modified particle filter for
searching high dimensional configuration spaces which does
not rely on such assumptions. It uses a continuation princi-
ple, based on annealing, to introduce the influence of narrow
peaks in the fitness function, gradually. The new algorithm,
termed annealed particle filtering, is shown to be capable of
recovering full articulated body motion efficiently.
2. Particle filters
Particle filtering (also known as the Condensation algorithm
[9]) provides a robust Bayesian framework for human mo-
tion capture. The Condensation algorithm was developed
for tracking objects in clutter, in which the posterior den-
sity
p
(
X
j
Z
k
)
and the observationprocess
p
(
Z
k
j
X
)
are often
non-Gaussian or even multi-modal (
X
denotes the model’s
configuration vector,
Z
k
=
f
Z
1
;:::;
Z
k
g
notates the his-
tory of observations at time
t
k
). The complicated nature of
the observation process during human motion capture causes
the posterior density to be non-Gaussian and multi-modal as
shown by Deutscher [3]. It is well knownthat a Kalman filter
will fail in this case. Deutscher et al were able to show that
the use of a particle filter will improvetracking performance.
The posterior density
p
(
X
j
Z
k
)
is represented by a set of
weighted particles
f
(
s
(0)
k
;
(0)
k
)
:::
(
s
(
N
)
k
;
(
N
)
k
)
g
where the
weights
(
n
)
k
/
p
(
Z
k
j
X
=
s
(
n
)
k
)
are normalised so that
P
N
(
n
)
k
= 1
. The state
X
k
at each time step
t
k
can be
estimated by
X
k
=
E
k
[
X
]=
N
X
n
=1
(
n
)
k
s
(
n
)
k
(1)
or the mode
X
k
=
M
k
[
X
]=
s
(
j
)
k
;
(
j
)
k
=max(
(
n
)
k
)
(2)
of the posterior density
p
(
X
j
Z
k
)
.
Particle filtering works well because it can model uncer-
tainty. Less likely model configurations will not be thrown
away immediately but given a chance to prove themselves
later on, resulting in more robust tracking. However a price
is paid for these attributes in computational cost. The most
expensive operation in the standard Condensation algorithm
is an evaluation of the likelihood function
p
(
Z
k
j
X
=
s
(
n
)
k
)
and this has to be done once at every time step for every
particle. To maintain a fair representation of
p
(
X
j
Z
k
)
a cer-
tain number of particles are required, and this number grows
with the size of the model’s configuration space. In fact it
has been shown by MacCormick and Blake [14] that
N
D
min
d
(3)
where
N
is the number of particles required,
d
is the number
of dimensions. The survival diagnostic
D
min
and the parti-
cle survival rate
are both constants with
<<
1
.Anex-
planation of both of these constants can be found in section
5. Clearly when
d
is large normal particle filtering becomes
infeasible.
Partitioned sampling was developed by MacCormick and
Blake[13] as a variationon Condensationto reduce the num-
ber of particles needed to track more than one object. Mac-
Cormick [14] has also now applied this technique to tracking
articulated objects. Using partitioned sampling reduces the
number of particles required to
N>
=
D
min
:
(4)
making the problem tractable. However, this assumes that
the configuration space can be sliced so that one can con-
struct an observation density
p
(
Z
k
j
x
i
k
)
for each dimension
x
i
k
of the model configuration vector
X
=
f
x
0
k
:::x
d
k
g
.This
assumption, that it is possible to independently localise sep-
arate parts of an articulated model, is similar to that make by
Gavrila to enable a hierarchical search. It has already been
argued that it is not possible to use this approach without the
use of labelling cues.
Another variation on the standard particle filter used to
reduce the number of particles needed to effectively repre-
sent a posterior density has been developed by Sullivan et
al [18]. Called layered sampling it is centered around the
concept of importance resampling. Experimental evidence
however suggests that this technique is not sufficient to solve
the problem of tracking with
d>
30
, reducing the number
of particles required by at best a factor of 5 to 10 before the
expected behaviour of the Condensation framework breaks
down.
The second reason why Bayesian particle filtering may
not be suitable for full body human motion capture is the
difficulties associated with constructing a valid observation
model
p
(
Z
k
j
X
k
)
as a normalised probability density distri-
bution. Another factor is the computational cost of calculat-
ing
p
(
Z
k
j
X
=
s
(
n
)
k
)
. Often an intuitive weighting function
w
(
Z
k
;
X
)
can be constructed that approximates the proba-
bilistic likelihood
p
(
Z
k
j
X
k
)
but which requires much less
computational effort to evaluate. Probabilistic observation
models also have a tendency to utilise only the information
that can be modelled well, discarding other available infor-
mation.
Given these factors it was decided to reduce the prob-
lem from propagating the conditional density
p
(
X
j
Z
k
)
us-
ing
p
(
Z
j
X
)
to finding the configuration
X
k
which returns
the maximum value from a simple and efficient weighting
function
w
(
Z
k
;
X
)
at each time
t
k
,given
X
k
,
1
. By doing
this gains will be made on two fronts. It should be possi-
ble to make do with fewer likelihood (or weighting function)
evaluations because the function
p
(
X
j
Z
k
)
no longer has to
be fully represented and an evaluation of a simple weighting
function
w
(
Z
k
;
X
)
should require minimal computational
effort when compared to an evaluation of the observation
model
p
(
Z
k
j
X
)
. The main disadvantage will be not being
able to work within a robust Bayesian framework.
It was decided to continue to use a particle based stochas-
tic framework because of its ability to handle multi-modal
likelihoods, or in the case of a weighting function, one with
many local maxima. The question is: What is an efficient
way to perform a particle based stochastic search for the
global maximum of a weighting function with many local
maxima? It was decided to use an approach which is sim-
ilar to that of simulated annealing.
3. Simulated annealing
The Markov chain based method of simulated annealing was
developed by Kirkpatrick et al [11] as a way of handling
multiple modes in an optimisation context. It employs a
series of distributions, with probability densities given by
p
0
(
x
)
to
p
M
(
x
)
, in which each
p
m
(
x
)
differs only slightly
from
p
m
+1
(
x
)
. Samples actually need to be drawn from
the distribution
p
0
(
x
)
. The distribution
p
M
is designed so
that the Markov chain used to sample from it allows move-
ment between all regionsof the state/search space. The usual
method is to set
p
m
(
x
)
/
p
0
(
x
)
m
,for
1 =
0
>
1
>
::: >
M
.
An annealing run is started in some initial state, from
which a Markov chain designed to converge to
p
M
is first
simulated. Some number of iterations of a Markov chain de-
signed to converge to
p
M
,
1
are simulated next, starting from
the final state of the previoussimulation. This process is con-
tinued in this fashion, using the final state of the simulation
for
p
m
as the initial state for the simulation for
p
m
,
1
, until
the chain designed to converge to
p
0
is finally simulated.
Note that if
p
0
contains isolated modes, simply simulat-
ing the Markov chain designed to converge to
p
0
starting
from some arbitrary point could give very poor results, as it
might become stuck in whatever mode is closest to the start-
ing point, even if that mode has little of the total probability
mass. The annealing process is a heuristic for avoiding this,
by taking advantage of the freer movement possible under
the other distributions. This is exactly the kind of behaviour
needed for the stochastic search. One wants to move towards
the global maximum of the weighting function
w
(
Z
k
;
X
)
,
using the overall trend of the matching function as a guide,
without becoming misguided by local maxima as seen in fig-
ure 1.
The idea of annealing for optimisation is now adapted to
perform a particle based stochastic search within the frame-
work of an annealed particle filter.
4. Annealed particle filter
A series of weighting functions
w
0
(
Z
;
X
)
to
w
M
(
Z
;
X
)
are
employed in which each
w
m
differs only slightly from
w
m
,
1
(see figure 2, where
M
=3
). The function
w
M
is designed
to be very broad, representing the overall trend of the search
space while
w
0
should be very peaked, emphasising local
features. This is achieved by setting
w
m
(
Z
;
X
)=
w
(
Z
;
X
)
m
;
(5)
for
0
>
1
> ::: >
M
,where
w
(
Z
;
X
)
is the original
weighting function. Because it is not the aim to sample from
w
(
Z
;
X
)
, but only to find its maximum it is not required that
0
=1
.
One annealing run is performedat each time
t
k
using im-
age observations
Z
k
. The state of the tracker after each layer
m
of an annealing run is represented by a set of
N
weighted
particles
S
k;m
=
f
(
s
(0)
k;m
;
(0)
k;m
)
:::
(
s
(
N
)
k;m
;
(
N
)
k;m
)
g
:
(6)
An unweighted set of particles will be denoted
S
k;m
=
f
(
s
(0)
k;m
)
:::
(
s
(
N
)
k;m
)
g
:
(7)
Each particle in the set
S
k;m
is considered as an
(
s
(
i
)
k;m
;
(
i
)
k;m
)
pair in which
s
(
i
)
k;m
is an instance of the multi-variate model
configuration
X
,and
(
i
)
k;m
is the corresponding particle
weighting. Each annealing run can be broken down as fol-
lows (the process is illustrated in figure 2).
1. For every time step
t
k
an annealing run is started at layer
M
, with
m
=
M
.
2. Each layer of an annealing run is initialised by a set of
un-weighted particles
S
k;m
.
3. Each of these particles is then assigned a weight
(
i
)
k;m
/
w
m
(
Z
k
;
s
(
i
)
k;m
)
(8)
which are normalised so that
P
N
(
i
)
k;m
=1
.Thesetof
weighted particles
S
k;m
has now been formed.
4.
N
particles are drawn randomly from
S
k;m
with replace-
ment and with a probability equal to their weighting
(
i
)
k;m
.Asthe
n
th
particle
s
(
n
)
k;m
is chosen it is used to
produce the particle
s
(
n
)
k;m
,
1
using
s
(
n
)
k;m
,
1
=
s
(
n
)
k;m
+
B
m
(9)
where
B
m
is a multi-variate gaussian random variable
with variance
P
m
and mean
0
.
5. The set
S
k;m
,
1
has now been produced which can be
used to initialise layer
m
,
1
. The process is repeated
until we arrive at the set
S
k;
0
.
6.
S
k;
0
is used to estimate the optimal model configuration
X
k
using
X
k
=
N
X
i
=1
s
(
i
)
k;
0
(
i
)
k;
0
:
(10)
7. The set
S
k
+1
;M
is then produced from
S
k;
0
using
s
(
n
)
k
+1
;M
=
s
(
n
)
k;
0
+
B
0
:
(11)
This set is then used to initialise layer
M
of the next
annealing run at
t
k
+1
.
w
(X)
0
X
S
(k,3)
Figure 1: Illustration of the annealed particle filter with M = 1.
Even though a large number of particles are used (so that an equiva-
lent number of weighting function evaluations are made as in figure
2), the search is misdirected by local maxima. From the resulting
weighted set it is very hard to tell where the global maximum of
w
0
lies.
5. Setting the tracking parameters
As stated previously the function
w
m
(
Z
k
;
X
)
, used in each
layer of the annealing process is determined by
w
m
(
Z
;
X
)=
w
(
Z
;
X
)
m
(12)
with
0
>
1
>:::>
M
.Thevalueof
m
will determine
the rate of annealing at each layer. A large
m
will produce
a peaked weighting function
w
m
resulting in a high rate of
annealing. Small values of
m
will have the opposite effect.
If the rate of annealing is too high the influence of local
maxima will distort the estimate of
X
k
as seen in figure 1. If
the rate is too low
X
k
will not be determined with enough
resolution (unless more layers are used wasting computa-
tional resources).
A good measure of the effective number of particles that
will be chosen for propagation to the next layer is the sur-
vival diagnostic
D
(taken from [14]) where
D
=
N
X
n
=1
(
(
n
)
)
2
!
,
1
(13)
and from this a good measure for the rate of annealing can
be derived, called the particle survival rate
[5] [12]
=
D
N
:
(14)
Now a measure for the rate of annealing has been derived
it is possible to set the values of
k
0
;:::;
k
M
at each time
step
t
k
. At layer
m
in an annealing run,
k
,
1
m
from
t
k
,
1
is used to calculate a preliminary set of particle weights for
S
k;m
. From this set an initial rate of annealing
init
can
S
S
π
S
S
2
w
(X)
w
(X)
w
(X)
w
(X)
3
1
0
X
X
X
X
k,3
k,3
k,2
π
S
k,2
π
S
k,1
π
S
k,0
k+1,3
S
k,1
S
k,0
Figure 2: Illustration of the annealed particle filter with M = 3.
With a multi-layered search the sparse particle set is able to gradu-
ally migrate towards the global maximum without being distracted
by local maxima. The final set
S
k;
0
provides a good indication of
the weighting function’s global maximum.
be calculated using equations 13 and 14. It can be shown
that
D
(
)
is monotonic decreasing in
so that, given
,the
equation
D
(
)=
N
(15)
has a unique solution for
. With this knowledge we can
minimise the error function
between the desired rate of
annealing
m
and the initial rate of annealing
init
(
)=(
m
,
init
(
))
;
(16)
using gradient descent to find the desired
k
m
. Note that
this does not mean the weights have to be completely re-
evaluated each time
k
m
is adjusted during gradient descent.
Since
w
m
(
Z
;
X
) =
w
(
Z
;
X
)
m
the values
w
(
Z
;
X
=
s
(
i
)
k;m
)
;i
: 1
:::N
can be stored for each set
S
k;m
and
k
m
applied to each individual weight as appropriate to produce
S
k;m
.
How then are the appropriate values for
0
:::
M
deter-
mined? There are also a number of other tracking parameters
1440 1540 1640
600
725
850
S
k−1,0
x
0
x
1
1440 1540 1640
600
725
850
x
0
x
1
S
k,3
1440 1540 1640
600
725
850
x
0
x
1
S
k,9
1440 1540 1640
600
725
850
x
0
x
1
S
k,2
1440 1540 1640
600
725
850
x
0
x
1
S
k,8
1440 1540 1640
600
725
850
x
0
x
1
S
k,1
1440 1540 1640
600
725
850
x
0
x
1
S
k,7
1440 1540 1640
600
725
850
x
0
x
1
S
k,0
Figure 3: Annealed particle filter in progress. The sets
S
k;m
are
plotted here, taken while tracking the walking person as seen in
figure 9. Only the horizontal translation components
x
0
and
x
1
of
the model configuration vector
X
are shown. Starting with
S
k
,
1
;
0
from the previous time step the particles are diffused to form
S
k;
9
which easily covers the expected range of translational movement
of the subject. The particles and are then slowly annealed over 10
layers (the sets
S
k;
6
to
S
k;
4
are omitted for brevity) to produce
S
k;
0
which is clustered around the maximum of the weighting function.
that need to be set before tracking can begin, including the
number of particles
N
, the number of annealing layers
M
and the diffusion variance vectors
P
0
:::
P
M
. A tentative
framework has been developed to allocate values to these pa-
rameters although it is acknowledged that more work needs
to be done in this area.
1. The first step is to decide on how many annealing layers
are needed. It was found that doubling the number of
annealing layers reduces the number of particles needed
for successful tracking by more than half. This will only
work up to a point however as there seems to be a min-
imum number (
N
) of particles needed for tracking no
matter how many layers are used. Using a 30 DOF model
it was found that setting
M
=10
with
N
200
worked
well.
2. Each element in the vector
P
0
is allocated a value equal
to half the maximum expected movement of the corre-
sponding model configuration parameter over one time
step. In this way the set
S
k
+1
;M
shouldcoverall possible
movements of the subject between time
t
k
and
t
k
+1
.The
amount of diffusion added to each successive annealing
layer should decrease at the same rate as the resolution
of the set
S
k;m
increases. It has been found that setting
P
m
=
P
0
(
M
M
,
1
:::
m
)
(17)
produces good results.
3. The appropriate rates of annealing
0
:::
M
are influ-
enced by the number of annealing layers used. With a
higher number of annealing layers a lower rate of an-
nealing can be used to obtain the desired resolution. It
was found that while using 10 annealing layers setting
0
=
1
=
:::
=
M
=0
:
5
provided sufficient resolu-
tion of
X
k
.
6. The model
The articulated model of the human body used in this pa-
per is built around the framework of a kinematic chain, as
seen in figure 4. Each limb is fleshed out using conic sec-
tions with elliptical cross-sections. It is believed that such a
model has a number of advantages including computational
simplicity, high-level interpretation of output and compact
representation.
(a) (b)
Figure 4: The model is based on a kinematic chain consisting of 17
segments (a). Six degrees of freedom are given to base translation
and rotation. The shoulder and hip joints are treated as sockets
with 3 degrees of freedom, the clavicle joints are given 2 degrees of
freedom (they are not allowed to rotate about their own axis) and the
remaining joints are modelled as hinges requiring only one. This
results in a model with 29 degrees of freedom and a configuration
vector
X
=
f
x
1
:::x
29
g
. The model is fleshed out by conical
sections (b).
7. The weighting function
When deciding which image features are to be used to con-
struct the weighting function a number of factors must be
taken into account.
Generality. The image features used should be invariant
under a wide range of conditions so that the same track-
ing framework will function well is a broad variety of
situations.
