1Long-term Prediction of Motion Trajectories Using Path Homology
Clusters
J. Frederico Carvalho1 , Mikael Vejdemo-Johansson2,3 , Florian T. Pokorny1 and Danica Kragic1

Abstract— In order for robots to share their workspace withpeople, they need to reason about human motion efficiently. In
this work we leverage large datasets of paths in order to inferlocal models that are able to perform long-term predictions ofhuman motion. Further, since our method is based on simple
dynamics, it is conceptually simple to understand and allows
one to interpret the predictions produced, as well as to extract a
cost function that can be used for planning. The main difference
between our method and similar systems, is that we employ
a map of the space and translate the motion of groups ofpaths into vector fields on that map. We test our method on
synthetic data and show its performance on the Edinburgh
forum pedestrian long-term tracking dataset [1] where we wereable to outperform a Gaussian Mixture Model tasked withextracting dynamics from the paths.

I. I NTRODUCTION
In order for robots to share their environment with people,
it is necessary for them to reason about human motion,and use this capability to not inconvenience those theyshare their workspace with. We propose to address thisproblem through a motion prediction algorithm that leverages
previously observed paths to predict the future motion ofpartial paths in an unsupervised manner.
Intuitively, we can formulate the problem as a model
learning problem, where we hypothesize that human paths
follow a certain model M (θ) and we want to find θ so thatthe error between the path predicted by the model, and theoriginal path is minimized.
Current methods for this type of motion prediction fall
into two categories, i) those which rely on the dynamics ofmotion, assuming that a subject moves in efficient paths, i.e.
straight lines, and any deviations from such motions are due
to obstacle avoidance which may include other agents [2]–[5].
ii) those that rely on training parametric models such as neuralnetworks or Hidden Markov Models (HMM) that produce a
prediction given past observations [6]–[9].
In the first instance, by assuming that human motion only
deviates from a straight line to avoid obstacles, the methods
ignore the possibility that subjects may have preferredpathways, and place the onus of predicting deviations fromefficiency on the exactness of the tracker, or the environment
model. Whereas in the second case, the obtained predictions
may be represented in the form of a probability density
as in [9] or a static model of occupancy as in [6] which
1 CAS/RPL, KTH, Royal Institute of Technology, Stocholm, Sweden.

{jfpbdc,dani,ftpokorny}@kth.se
2 Mathematics Department, CUNY College of Staten Island, New York,USA. mvj@math.csi.cuny.edu
3 Computer Science, CUNY Graduate Center, New York, USA.

Fig. 1. Illustration of our method, a pedestrian path in an environment,
together with two vector fields generated by pedestrian paths (red and blue),and the resulting different predicted paths (dashed).

necessitate one to be aware of the threshold of occupancy
probabilities when using such methods for planning.
In this paper we expand our previous work in pathclustering [10], [11] and compute vector fields which represent
the average motion of paths in each cluster. We exploit
the nature of our distance metric, which ensures a cluster
should cover a small area of the environment to ensure thatdifferent vector fields encode motions in different parts of theworkspace. We employ the resulting fields to predict motion
by integrating the field that most closely resembles the pathobserved so far.
To the best of our knowledge, our method is the first to
employ insights obtained from the topology of the spaceof paths to extract higher order dynamics and use them for
motion prediction.
II. R ELATED WORK
Motion prediction is a fairly well-studied task and is closely
related to object tracking, therefore rather similar methods
can be used to address it. In [5] the authors present a Kalmanfilter based method for predicting the motion of pedestrians
in 2D scene. In [8] The authors use HMMs to predict the
motion of construction workers in a cluttered environment
in order to improve building site safety. Similarly, in [9]the authors employ an HMM to the motion of people in an
office environment, using both the immediate and far pastobservations to obtain a distribution over future observations.
Similarly recurrent neural network models, have been usedin [7] to infer driver intention in a roundabout. In [12] theauthors employ Support Vector Machines (SVM) and Kmeans clustering with respect to the edit distance on a datasetof paths as a training step, subsequently future motions are
predicted from a partial path by matching it to one fo theclusters.

2

One important application of motion prediction is the taskof compliant motion planning, since by predicting the motionof a pedestrian in a scene, for example, a robot is able toplan so as to not interfere with the pedestrian’s future path.
In [13] the authors cluster paths to extract motion patterns
using an expectation-maximization (EM) algorithm, thesemodels are then used for prediction using a HMM which areused to synthesize compliant motion. An alternative approach
was presented in [6] where the authors use a CNN to infer
which parts of an environment are most traversed by humans,
and use this in order to generate a robot motion that is nonobstructive. In [14] the authors use trajectory clustering inorder to predict the motion of an agent, so that they are able
to plan an intersecting path.
An area related to motion prediction is learning bydemonstration. Here instead of using a model to infer wherean object is moving toward, the model is used to instruct a
robot how to move in order to perform some task. Popular
approaches to this problem employ Gaussian Processes (GP)such as in [15] where the authors use local GP regression
to perform a given task from a handful of demonstrations.
Another common choice of model is a Gaussian Mixture
Model (GMM) such as in [16] where the authors train a
GMM from a handful demonstrations in order to obtain acompact model that generalizes a given task.
One way to extract motion primitives is through pathclustering, where one intends to group a dataset of paths
into subsets of paths that are pairwise similar. Generally,
this is achieved through distance-based clustering methods,
and a recent survey of such distance based methods canbe found in [17]. Another possibility is to concentrate onpath clusters that are fundamentally different with respect totheir relationship to obstacles in the environments as in [10],and such methods can even be combined with a compatible
distance function [18] to obtain more efficient clustering, as
we showed in [11].
In this paper we present an algorithm for extracting motion
primitives from datasets of paths in an efficient way. We do
this by employing our path clustering pipeline from [11] toproduce groups of paths that are deemed to be similar, andaggregating the observed velocities into a set of vector fields.
To predict motion of a path up to a point, we choose thevector field which is most similar to the path observed so far,and extrapolate the motion by integrating that vector field.
The remainder of the paper is structured as follows: In
Section III we give a brief explanation of simplicial complexes
and their homology, followed by a short explanation of ournotation conventions on vector fields. These will be used toextend our path clustering pipeline from [11], and to establish
the establish the main contributions of the paper, respectively.
In Section IV, we provide the aforementioned extensions
to our path clustering pipeline, followed by an algorithm tocalculate vector fields based on those path clusters as well
as describe how we employ these vector fields for motion
prediction. In Section V we describe our construction ofsynthetic datasets, and present the results from running ourpipeline on such datasets, followed by the results of running

the aforementioned pipeline on the dataset from [1].
III. BACKGROUND
In this section we give a short overview of the backgroundand terminology that will be used throughout the rest of thepaper. Namely, we provide a short introduction to simplicial
complexes and their homology, which will be necessary inorder to prove that our alteration to the path distance metricfrom [11] produces the same clustering results. We also
provide a short overview of our notation for vector fields
which form the basis of our method for motion prediction.
A. Simplicial Complexes
Here we give a very short introduction to the notion of
simplicial complexes. More detailed summaries can be foundin [11], [19] as well as in more classical sources such as
chapter 2 of [20].
Intuitively, a finite simplicial complex is a generalizationof a triangulation. Formally it can be specified by a set X
of subsets of {1, . . . , N } such that for any σ ∈ X, and any
non-empty τ ⊂ σ, τ ∈ X. Each element τ ∈ X is called aface of X.
If we consider a set of N points in R2 , X can be seen
as a triangulation where {xi | i ∈ σ} is the set of verticesof a triangle or edge when |σ| = 3, 2. However, we haveto make a further restriction that convex hulls of [σ] :=Conv({xi | i ∈ σ}) satisfy [σ] ∩ [τ ] = Conv({xi | i ∈(σ ∩ τ )}) only intersect at shared faces. The dimension of asimplex σ ∈ X is defined given by dim(σ) = |σ| − 1 andthe dimension of the complex X is given by dim(X) =maxσ∈X dim(σ). We denote X (i) as the i-th level of X,
containing
all the i-dimensional simplices of X and X i =Si(j), which is a simplicial complex called the i-th
j=0 X
skeleton.
A simplicial complex X gives rise to a chain complex
X{Ci (X)}dim
i=1 . This corresponds to a sequence of vector
spaces, where Ci (X) has X (i) as a base together with afamily of boundary maps ∂i : Ci (X) → Ci−1 (X) which senda basis element to a formal sum of its boundary elementswith signs shifted according to the induced orientation, amore extensive description can be found in chapter 2 of [20].
An element c ∈ Ci (X) is called an i-chain of the complex,
furthermore if ∂i c = 0 then c is called a i-cycle and if there
exists some b ∈ Ci+1 (X) such that ∂i+1 b = c, then c iscalled a boundary.
It can be shown that in this setup every boundary is a cycle,
and that cycles and boundaries form subspaces of the Ci (X).
The i-th homology group of X consists of equivalence classesof cycles in Ci (X) which are not boundaries.
When dealing with paths, we view these as an ordered
list of (oriented) edges on the simplicial complex which areplaced end-to-end. These we regard as 1-chains of X, withinteger weights. For example, a weight of 2 in a particularedge indicates that the edge is traversed twice, and a weightof −1 indicates that the edge is traversed once but in thenegative direction. We call a 1-chain obtained from sucha path a path chain. Note that our definition of simplicial

3complex implies that any two vertices share at most one edge,
therefore a path chain p is uniquely determined by an ordered
list of vertices p0 , p1 , . . . , pn .
We say that two path chains p, q are homologous if p − q isa boundary, i.e. if there exists a 2-chain cp,q so that ∂cp,q =p − q. In this case cp,q is said to be the bounding chain ofp − q. However, the requirement that p − q forms a cycle
is in general too restrictive, therefore it will be useful toconsider quotients of X. The quotient X/A is defined fromthe simplicial complex X by identifying every simplex in thesubcomplex A ⊂ X with a single point in the quotient.
B. Vector Fields
Given a subset of X ⊆ R2 we define a vector field V on X
as a continuous function V : X → R2 . Such a vector field canmodel the motion of a particle by setting ẋ = V (x). We definethe (positive) flow of a vector field starting at x0 as a function
ΦV,x0 (t) : J → X where J is an interval in R, containingd0 and ΦV,x0 (t) satisfies dt
ΦV,x0 (t) = V (ΦV,x0 (t)) andΦV,x0 (0) = x0 . The existence of J and ΦV,x0 in some region
around x0 is guaranteed by the Picard-Lindelöf theorem [21]provided that the vector field V (x) is Lipschitz continuousin x in some neighborhood.
We employ a discrete approximation of vector fields inthe plane that consists of an assignment P 7→ R2 where P
is a finite subset of R2 . This can then be extrapolated to avector field in each simplex of a triangulation of P by linear
interpolation. Due to this approximation we guarantee thatthe obtained vector field is Lipschitz continuous.
IV. M ETHOD
We begin by describing a small alteration to our pathclustering pipeline. In [11] we used a model of the quotient of
the base space X by a subspace A ⊂ X which encompasses
the endpoints of the paths in the dataset. As noted Section IIIA, this allows us to treat path chains as cycles and use thatto find if they are homologous. We then defined the distance
between the paths to be either the area of the bounding chain
if tey are homologous, and infinity otherwise.
In the following we show that under mild assumptions a pairof paths in X/A can only be homologous if they begin in thesame component of A and end in the same component of A.
Therefore, under this condition the results from [11] are valid
when instead of considering X/A we use the decomposition of
A into its connected components A1 ∪, . . . , ∪Ak and consider
instead the iterated quotient (. . . (X/A1 )/ . . .)/Ak which wedenote by X/(A1 , . . . , Ak ). All proofs except for Lemma 1
are contained in appendix.
The following lemma establishes that the results we obtain
are independent of the ordering of A1 , . . . , Ak .
Lemma 1: Let A1 , . . . , Ak be disjoint subsets of X and
let Sk be the set of permutations of k elements, then for any
σ ∈ Sk we have X/(A1 , . . . , Ak ) ∼= X/(Aσ1 , . . . , Aσk ).
For convenience, we will further restrict the subcomplexes
A1 , . . . , Ak to satisfy the following:

Definition 1: We say that a sequence of subcomplexes
A1 , . . . , Ak ⊆ X are well disconnected if for every simplex
σ ∈ X, σ ∩ Ai 6= ∅ implies σ ∩ Aj = ∅ for all j 6= i.
The following proposition is a reformulation of Proposition 2 from [11] which establishes the uniqueness of bounding
chains in our setting and guarantees that the distance function
is well defined.
Proposition 1: Let X be a simplicial complex which canbe embedded in R2 and A1 , . . . , Ak ⊆ X disconnectedsubcomplexes such that H1 (Ai ) = 0, then given any 1-cycle
c of X/(A1 , . . . , Ak ), there exists at most one 2-chain s suchthat ∂s = c.
The following lemmas characterize the notion of cycle andboundary arising from path chains in an iterated quotient
X/(A1 , . . . , Ak ) when compared to the quotient X/A.
Lemma 2: Given two path chains p = p0 , . . . , pn and q =q0 , . . . , qm in a simplicial complex, then p − q is a cycle if
and only if both p0 = q0 and pn = qm or p and q are bothcycles.
Lemma 3: Let X be a simplicial complex embedded in R2
and A1 , . . . , Ak be well disconnected subcomplexes so that
H1 (Ai ) = 0 . Let c be a path chain in X, then the associatedchain [c] is a boundary in C1 (X/(A1 ∪ . . . ∪ Ak )) if and only
if it is a boundary in C1 (X/(A1 , . . . , Ak )).
These results summarize the situation for path chainsin the iterated quotient X/(A1 , . . . , Ak ) for X ⊂ R2 andA1 , . . . , Ak well disconnected. Namely, two path chains p, q
are homologous if and only if they were homologous in X/A,
and if they are homologous, there is at most one bounding
chain. This entails that our notion of area homology distancecarries over to this setting without change.
A. Calculating Vector Fields
We define a vector field starting from the first orderdifferences of the path and extrapolating to the rest of
the space by using an underlying simplicial complex. In
Algorithm 1 we provide the pseudocode for generating sucha vector field from a path. The remainder of this subsection
is dedicated to explaining how this algorithm works, and wedefer to the next subsection to explain how we use the vector
fields for motion prediction.
Definition 2: Let γ = (t0 , p0 ), (t1 , p1 ), . . . , (tn , pn ) be apath, then we define the differential of γ at a given point pi
to be given by:
pi+1 − pidγ(ti ) =.
ti+1 − ti
The differential defines a vector field associated to γ but
only at the points pi . To extrapolate this field to other points in
R2 we first extrapolate it to neighboring points by averaging
the contributions from neighbors (as illustrated in Fig. 2).
The main idea of Algorithm 1 is to establish the vector
field along the path, and then use the edges of the mesh S toexpand it to the rest of the space. This allows one to define
the vector field as a simple N × 2 matrix where N is thenumber of vertices S. For points which are not vertices of S
the value of the vector field can be established through linear
interpolation.

4

Algorithm 1: Calculate the vector field associated to a path
input : S = (V, E)
input : P = (p1 , t1 ), ..., (pn , tn )1 Lvi ← ∞where i ← 1, . . . , size(V );2 N ptsi ← 0
where i ← 1, . . . , size(V );pi+1 −pi
3 dPi ← (pi , t)where i ← 1, . . . , n;i+1 −ti
4 vPi ← (0, 0)
where i ← 1, . . . , size(V );5 for (p, dp) ← dP0 , dP1 , . . . do
6q ← nearest(V, p);
7i ← index(q, V );8d ← dist(q, p);
9vPi ← vPi + DECAY(d)dp;
10
N ptsi ← N ptsi + 1;

Fig. 2. Illustration of the process by which the vector field is expanded toneighboring vertices. The vertices along the line in red are in the path, andthe black arrows on them indicate the value of the vector field along that
path. The lines in light gray are the edges of the mesh, and the arrows inblue indicate which elements are averaged together before into the arrow indark gray that represents the vector field at the opposite (w.r.t. the red path)end of the edge.

N extLv ← {};for i ← 0, . . . , size(V ) do
13n ← N ptsi ;14
if n > 0 then15
vPi ← vPi /n;
16
Lvi ← 0;17
append(N extLv, neighborIdx(S, Vi ));11
12

CurrLv ← 1;
while notEmpty(N extLv) do
20
N ewLv ← {};21
for i ← N extLv0 , N extLv1 , . . . do
22
if Lvi < CurrLv then
23
skip

18
19

24
25
26
27
28
29
30

N bh ← neighborIdx(S, Vi );N ← 0;
for j ← N bh0 , N bh1 , . . . do
if Lvj < CurrLv thend ← dist(Vi , Vj );vPi ← vPi + DECAY(d)vPj ;N ← N + 1;
else

31

append(N ewLv, j);

32
33
34
35
36

i
vPi ← vP;N

N extLv ← eliminateRepeated(N ewLv);
CurrLv ← CurrLv + 1;

Fig. 3. Example of a path and the associated vector fields. Depicted are asimple curve γ in red with its associated vector field Vγ in blue computed
by Algorithm 1. As we can see, the field forms a “corridor” around the path.

the vector field vP at those neighbors that are at a level lower
than the current level CurrLv. The neighbors that have notbeen assigned a level yet are then added to the new level set
N ewLv. At the end of this loop the new level set replaces
the level set and the current level is incremented.
Finally, in order to obtain the vector field associatedto a cluster of paths Γ = {γ1 , . . . , γn } we will have toaverage the different vector fields. To this end, we need tomake the different fields comparable, which we achieve byrenormalizing them so that the highest velocity is 1:
Wγ (x) =return vP ;Algorithm 1 begins by initializing the vector field vP inline 3 to be zero at every point, and the variable dP inline 2 which stores the differential of the path P . Two further
auxiliary variables are set: N ptsi , which is the number ofpoints that are averaged when producing the estimation at thei-th vertex of the mesh, and Lvi which indicates at which“level” the i-th point has been calculated. Here i is assumed
to range over the index of all points in V .
In lines 5 to 10, the vector field is calculated at the nearest
neighbors of the points in the path, using to this effect thevalues of dP . The expansion of the vector field to the mesh
begins in earnest with the loop in lines 12 to 17. In it, thepoints which were previously assigned a vector field value
through the nearest neighbor query, are assigned level zero,
whereas their neighbors are inserted into the variable N extLv
which contains the next “level” of points to be analyzed.
Finally, the main loop in lines 19 to 35 proceeds bychecking the neighbors of each element of N extLv which
has not been previously analyzed, and averaging the values ofVγ (x)
supy∈X Vγ (y)Similarly, we obtain the vector field VΓ associated to the cluster Γ by combining the vector fields Wγ and renormalizingthem the same way:WΓ (x) =X

Wγ (x)VΓ (x) =γ∈ΓWΓ (x)
supy∈X WΓ (y)B. Cluster assignment and motion prediction
With the vector fields calculated we can finally assign
a hypothesis cluster to a path that has only been partiallyobserved. Namely, for each vector field we measure “howparallel” it is to the path that is being examined, and thecorresponding cluster is then the cluster that is assigned i.e. fora given path η = (t1 , η1 ), . . . , (tk , ηl ) we assign the cluster:Ξ(η) = arg maxlX

hdη(tj ), VΓi (η(tj ))ii=1,...,k j=1

and we predict the path to be the integral path of the vector
field starting at the last point in the path while maintaining

5

the current velocity by (numerically) solving the integral
equation:
Z tVΞ(η) (η̃(s))η̃(t) = η(tl ) + kdη(tl )kds
tl kVΞ(η) (η̃(s))kV. E XPERIMENTS
In this section, we present our experimental results. We startwith two synthetic examples, where the data is drawn from aknown distribution. The first set exemplifies an ideal scenario
where the data is well separated, and our score function is
ideally suited to distinguish between the different sets ofdata. The second dataset contains more complex dynamics,
and showcases situations in which our score function canunderperform and be led to produce misclassification errors.
Finally, in the real world data example, we show how ourmethod performs on the dataset from [1].
We compare our method with simple dynamics based
motion prediction methods, the first one simply integrates
the current speed of the path into the future, i.e. given thepath γ up to time t we predict that at time t + s the pathis at γ(t) + sγ̇(t). The second can only be used when there
is a known goal region, which is the case for the synthetic
datasets, and works by setting an attractor on a target in the
S in the goal region but mantains the norm of the velocity,
so in essence we solve the initial value problem:
d
dt





kγ̇(t)kγ(t + s)γ̇(t + s)=γ̇(t + s)kγ̇(t + s)k β(S − γ(t + s))

where the parameter β is chosen to minimize the prediction error. In the real world dataset we use an alternative vector fieldmodel, that consists of modeling the vector field associatedto each path cluster using a GMM as in [22]. All integrations
are calculated using an explicit fourth order Runge-Kutta
method with a step size of 0.01 on the synthetic datasets,and (1/90)s−1 on the real-world dataset, this corresponds toten steps per frame, as the paths were drawn from a tracking
process running at 9 frames per second.
A. Simple Synthetic Data
In this test we employed a simple model where we draw
paths in a 100 × 100 square in R2 (the distance units can betaken to be meters) and demark four regions:• left near point (0, 50) where all paths begin.
• right near point (100, 50) where all paths end.
• top near point (50, 100) which contains the midpointsof all paths in cluster 1.
• bottom near point (50, 0) which contains the midpointsof all paths in cluster 2.
Each path is interpolated from three points (begining, middle
and end with timestamps 0, 0.5 and 1 respectively) using aGP, the resulting paths can be seen in Figure 4. To draw apath from the distribution we first randomly choose either the
top or bottom region, then we draw the prescribed three way
points that define the path, and interpolate them as specified.
To test in this scenario we sampled randomly 1000 paths
from the distribution and clustered them into two clusters,
and calculated the associated vector fields. We tested theFig. 4. Simple synthetic data set: The first cluster of paths is depicted inred, and the second one is depicted in yellow. To the right we show thestream plot with the vectorfields associated to each cluster in the respective
color.

Fig. 5. Top row: Both plots contain a path in red drawn from the same
distribution as the dataset used to calculate the vector fields. Each of theblue paths is predicted with our method starting from the large blue point at
its extremity. Bottom row: Error for the predicted paths in the top row. The
horizontal axis represents time, and the vertical axis represents the distancefrom the ground truth. Each of the blue,red and green lines represents themeasured error when motion is predicted using our method (blue), constantvelocity prediction (red) and an attractor at the goal region (green).

performance of the resulting model on 100 more paths also
drawn from the distribution.
In Figure 5 we show an example of predictions from our
method, together with the error incurred by both our method,
and the constant-velocity and the attractor-towards goalmodels. The error is represented in time (of the path/predicted
path which starts at the point where the prediction starts)
versus distance between the ground-truth and the predicted
path. As we can see our method is able to track the overall
shape of a simple path, which on these instances allows it tooutperform both other predictors.
In Figure 6 we present the statistics of the errors for thethree different models, namely we plot the prediction horizon(that is, amount of time after the beginning of the prediction)versus error (again, distance between the predicted path andthe ground truth at the given prediction horizon). For eachmodel we plot the average (solid line) and standard deviation(shaded region) of the error at any given time. As we cansee our observations from Figure 5 remain valid when we
look at statistics of the error.
Note that in this case we used partial information from how
the dataset was generated to calculate the vector fields, namely
we used two path clusters, which generated two vector fields,
each of which closely tracks the associated cluster, therefore

6Fig. 6. Error statistics for the prediction models for the dataset in Figure 4.
Here we plot the statistics for the three different models, our model is inblue, the constant velocity prediction is in red, and the attractor model is ingreen. The statistics shown are the mean error (solid line) and the standard
deviation region (shaded region).

Fig. 8. Top row: Two paths from the dataset Figure 7 together withpredictions performed using our method, as in the top row of Figure 5.
Bottom row: Error for the predicted paths in the top row, the legend is thesame as in the bottom row of Figure 5.

as long as the correct vector field is identified, the prediction
error should remain low. Furthermore, since the vector fields
almost do not overlap in space, and when they do, they aremostly perpendicular to each other, identifying the correct
vector field becomes a simple task. Next we show anothersynthetic example where this separation between fields is
not well demarcated, and therefore incurs in higher error
predictions for our model.

compare it with the two simple dynamics methods used in
the previous example in the bottom row. From loooking at
the predictions we notice that one of them follows the wrong
vector field leading to a gigh prediction error. This happens
because at the point of prediction all the fields are similar,
this implies that the score functions will also have a similarvalue and leads to the possibility of missclassifications which
may result in higher error than the simple models.

B. More Complex Synthetic Data
With this example we aim to show that our model is able toencode more complex dynamics than a simple single curved
motion towards a target, and to show where prediction errors
may arise from in our model. To generate the data we usean almost identical procedure to the one before, except wehave an extra “foreword” region near the point (25, 0) where
all paths go through. Once again, to draw a path from one
cluster we draw four points (one from the beginning region,
one from the foreword region, one from the top or bottom
region, and one from the end region). The points are thengiven timestamps 0, 0.2 , 0.7 and 1 respectively, and areinterpolated with a Gaussian process.
Fig. 9. Error statistics for the prediction models for the dataset in Figure 7.
Here we plot the statistics for the three different models, our model is inblue, the constant velocity prediction is in red, and the attractor model is ingreen. The statistics shown are the mean error (solid line) and the standard
deviation region (shaded region).

Fig. 7. A more complicated dataset, comprising four path clusters (red,blue, black and yellow). To the right we see the associated vector fields.

From the description of the generated dataset it would
be expected that there would be only two clusters as in theprevious example, however, our clustering method proposed amodel with four clusters as seen in Figure 7 when maximizing
the ratio of inter-cluster distance to in-cluster distance. In
Figure 8 both demonstrate our method in the top row, andHowever, when we once again analyze the error statistics
for the predictions in this dataset, we see that our model once
again outperforms the simple dynamics models. However, themean and standard deviations of the error at all horizons are
now larger than they were in the simple dataset which comesfrom the fact that the possibility of misclassification leads tolarger errors at all prediction horizons.
C. Real-world Data
We now repeat the analysis performed for the synthetic
datasets, using instead the Edinburgh dataset [1]. This dataset

7comprises tracks of pedestrians in Edinburgh University, andpart of it is pictured in Figure 10 (left).

Fig. 10. Subset of the path clusters (left) and corresponding vector fields(right) associated to paths from the Edinburgh pedestrian tracking dataset [1]which share origin and destination regions.

However, due to high-frequency noise from measurement
errors. For example, the positions output by the tracker aremeasured in pixels, so even when an object is stationary
the position is expected to change due to naturally shiftinglighting conditions. We thus need to preprocess the data sothat velocities vary smoothly as expected. To this end weassume the noise is white, so that it can be filtered by using asimple 5-point moving average (x̂t = (xt−2 + · · · + xt+2 )/5).
In Figure 10 (right) we show a picture of the vector
fields generated from the paths in the left, as we cansee the vectorfields have a large overlap, which results inmisclassification errors, as we previously pointed out. As
we also pointed out before, we compare our method tousing the simple constant velocity prediction and a GMM
vector field approach [22]. In this model we calculate allpairs (γ(t), γ̇(t)) from paths γ in some cluster Γ, which arethen assumed to be samples drawn from a random variable(XΓ , VΓ ) with domain R2 ×R2 and with according to a GMM.
To choose the vector field to integrate in order to predicta new pathPγ we select the cluster Γ that maximizes theTlikelihood t=0 log(p(XΓ ,VΓ ) (γ(t), γ̇(t))) and we calculate
the value of eth vector field at point x as being the expectedvalue E(VΓ |XΓ = x).
However, when analyzing the predictions in Figure 11, we
see that once again, in the long term our method is able
to attain a smaller error than simply following the current
velocity or using the GMM vector field model. Upon a closer
inspection and comparing with Figure 10, we see that thepaths on the top row are not always attributed to the same
vector field. Namely the first two predictions of the path on
the right follow the red cluster, the third prediction follows
the blue cluster as the path curves upward more abruptly ina region where the red cluster does not.
Finally, the statistics in Figure 12, second our observations
from the analysis of Figure 11 and we we can see that our
method outperforms on average constant-velocity extrapolation and integrating the GMM vector field model, while also
attaining lower variance.
VI. C ONCLUSIONS AND F UTURE W ORK
In this paper we presented a method to extract motion
primitives in a scalable way from large datasets of paths
using piece-wise linear vector fields. We compared ourFig. 11. Top row: two paths from the Edinburgh pedestrian tracking dataset
and associated predictions using our method, starting at different points
along the path. Bottom row: comparison of the prediction error using our
method (blue) and constant velocity prediction (red), and GMM vector field
integration (green).

Fig. 12. Error statistics for the prediction models for the Edinburgh
pedestrian tracking dataset. Depicted are the error statistics for our model(blue), constant velocity prediction (red) and the GMM vector field (green).
The statistics shown are the mean error (solid line) and the standard deviation
region (shaded region).

method to other dynamics based methods, and observed that
our approach is able to effectively capture more complex
dynamics than conceptually simpler primitives, and in thetested datasets lead to more faithful predictions than a GMM
vector field approach at all time horizons.
In future work we hope to improve our framework inseveral ways, namely by using sparse vector field models,instead of the dense model currently used, and using insightsfrom sheaf theory [23] in order to provide more resilient
prediction dynamics. We also intend to integrate our method
into a tracking framework, for example a bayesian filter totest the system performance with our system dynamics.
A PPENDIXP ROOFS OF RESULTS FROM S ECTION IV
A. Proof of Proposition 1:
We follow an argument similar to the proof of Proposition 2
in [11]. Namely when A ⊆ X is a subcomplex, we can build
the exact sequence of the pair [20].
· · · → Hn (A) → Hn (X) → Hn (X/A) → Hn−1 (A) → · · ·8Since this sequence is in reduced homology, and each
Ai is of dimension 2 and embedded in R2 , therefore has
H2 (Ai ) = 0, furthermore it satisfies H1 (Ai ) = 0 and since
it is connected, it also satisfies H0 (Ai ) = 0 (since it isreduced homology). Therefore Hn (Ai ) = 0 for every n
and there is an isomorphism in homology between Hn (X)and Hn (X/Ai ) furthermore since the Ai are disjoint, thesubcomplex Ai /Aj ⊆ X/Aj is isomorphic to Ai and so stillsatisfies Hn (Ai ) = 0.
Now, this implies that in the conditions of the proposition
Hn (X/(A1 , . . . , Ak )) ∼= Hn (X) for every n, and therefore
H2 (X/(A1 , . . . , Ak )) = 0 which implies that for any 1-chainc there exists at most one 2-chain s so that ∂s = c.
B. Proof of Lemma 2:
Since p is a path chain, ∂p = (p1 − p0 ) + (p2 − p1 ) +. . . + (pn − pn−1 ) = pn − p0 . In a similar fashion we have∂q = qm − q0 and so ∂(p − q) = 0 if and only if pn − p0 =qm − q0 , i.e. if and only if pn = qm and p0 = q0 or pn = p0
and qm = q0 .

C. Proof of Lemma 3:
Assume, without loss of generality that k = 2, and therefore
the relative region has only two connected components A1 ,and A2 . This decomposition induces a factorization of thequotient X → X/(A1 ∪ A2 ) through the iterated quotient
X/(A1 , A2 ) (which corresponds simply to the identificationof the points representing each of the connected components).
Call these maps φ and ψ:
φ

X −→

X
X
ψ
−→
A1 , A2A1 ∪ A2

Since φ and ψ are simplicial maps, they induce chain
homomorphisms φ∗ ([c]) = [φ(c)], and ψ∗ ([d]) = [ψ(d)]. By
definition of chain homomorphism, they commute with thedifferential, meaning ∂ ◦ ψ∗ ([d]) = ψ∗ (∂[d]). Therefore ifφ∗ ([c]) is a boundary, then it can be written as ∂[b] for some2-chain [b] which means.

R EFERENCES
[1]
[2]
[3]
[4]

[5]

[6]
[7]
[8]
[9]

[10]

[11]
[12]
[13]
[14]

ψ∗ (φ∗ ([c])) = ψ∗ (∂[b]) = ∂(ψ∗ ([b]))[15]i.e. ψ∗ ◦ φ∗ ([c]) = [ψ ◦ φ(c)] is also a boundary. To prove thereciprocal assume now that c is a path chain, and [ψ◦φ(c)] is aboundary. Note once again, that all that ψ does is identify the
points representing A1 and A2 , leaving every other simplex in
X/(A1 , A2 ) intact. Furthermore if we let the points in A1 , A2
be the first two points in the ordering of (X/(A1 , A2 ))(0) ,we note that since A1 , A2 are well disjoint, so they do notshare an edge, and hence ψ preserves simplex orientations.
With that in mind, since any [b] such that ∂[b] = ψ∗ ◦ φ∗ ([c])
is a 2-chain, we can calculate the preimage of [b] by ψ∗ ,which is itself a well defined 2-chain [b∗ ] satisfying:∗∗ψ∗ (φ∗ ([c])) = ∂ψ∗ ([b ]) = ψ∗ (∂[b ])And again since ψ∗ identifies the 0-simplices corresponding
to A1 and A2 it is an isomorphism on 1-chains, and thereforeφ∗ ([c]) = ∂[b∗ ] meaning, φ∗ ([c]) is a boundary.

ACKNOWLEDGEMENTS
The authors would like to thank Patric Jensfelt for his help
in the writing of this paper. This work was supported by
the Knut and Alice Wallenberg Foundation and the Swedish
Research Council.

[16]

[17]
[18]
[19]
[20]
[21]
[22]
[23]

B. Majecka, “Statistical models of pedestrian behaviour in the forum,”
Master’s thesis, School of Informatics, University of Edinburgh, 2009.
S. Kim, A. Bera, A. Best, R. Chabra, and D. Manocha, “Interactive
and adaptive data-driven crowd simulation,” in 2016 IEEE Virtual
Reality (VR), Mar. 2016.
D. Helbing and P. Molnár, “Social force model for pedestrian
dynamics,” Physical Review E, 1995.
S. Kim, S. J. Guy, W. Liu, D. Wilkie, R. W. Lau, M. C. Lin, and
D. Manocha, “Brvo: Predicting pedestrian trajectories using velocityspace reasoning,” The International Journal of Robotics Research,
vol. 34, no. 2, 2015.
A. Bera, S. Kim, T. Randhavane, S. Pratapa, and D. Manocha, “Glmp- realtime pedestrian path prediction using global and local movement
patterns,” in 2016 IEEE International Conference on Robotics andAutomation (ICRA), May 2016.
J. Doellinger, M. Spies, and W. Burgard, “Predicting occupancy
distributions of walking humans with convolutional neural networks,”
IEEE Robotics and Automation Letters, vol. 3, no. 3, Jul. 2018.
A. Zyner, S. Worrall, and E. Nebot, “A recurrent neural network
solution for predicting driver intentions at unsignalized intersections,”
Robotics and Automation Letters, vol. 3, no. 3, Jul. 2018.
K. M. Rashid and A. H. Behzadan, “Enhancing motion trajectory
prediction for site safety by incorporating attitude toward risk,” inComputing in Civil Engineering 2017, 2017.
Z. Wang, P. Jensfelt, and J. Folkesson, “Multi-scale conditionaltransition map: Modeling spatial-temporal dynamics of human
movements with local and long-term correlations,” in 2015 IEEE/RSJ,
IROS 2015, Hamburg, Germany, September 28 - October 2, 2015,
2015.
F. T. Pokorny, K. Goldberg, and D. Kragic, “Topological trajectory
clustering with relative persistent homology,” in 2016 IEEE International Conference on Robotics and Automation (ICRA), IEEE, May
2016.
J. F. Carvalho, D. Kragic, M. Vejdemo-Johansson, and F. T. Pokorny,“Path clustering with homology area,” in 2018 IEEE International
Conference on Robotics and Automation (ICRA), IEEE, May 2018.
S. Xiao, Z. Wang, and J. Folkesson, “Unsupervised robot learning topredict person motion,” in 2015 IEEE International Conference on
Robotics and Automation (ICRA), May 2015.
M. Bennewitz, W. Burgard, G. Cielniak, and S. Thrun, “Learning motion patterns of people for compliant robot motion,” The International
Journal of Robotics Research, vol. 24, no. 1, 2005.
C. Sung, D. Feldman, and D. Rus, “Trajectory clustering for motion
prediction,” in 2012 IEEE/RSJ International Conference on Intelligent
Robots and Systems, Oct. 2012.
M. Schneider and W. Ertel, “Robot learning by demonstration withlocal gaussian process regression,” in 2010 IEEE/RSJ International
Conference on Intelligent Robots and Systems, Oct. 2010.
S. Calinon, T. Alizadeh, and D. G. Caldwell, “On improving theextrapolation capability of task-parameterized movement models,” in2013 IEEE/RSJ International Conference on Intelligent Robots andSystems, Nov. 2013.
S. Aghabozorgi, A. Seyed Shirkhorshidi, and T. Ying Wah, “Timeseries clustering – A decade review,” Information Systems, vol. 53,
Oct. 2015.
E. W. Chambers and M. Vejdemo-Johansson, “Computing minimum
area homologies,” Comput. Graph. Forum, vol. 34, no. 6, 2015.
J. F. Carvalho, M. Vejdemo-Johansson, D. Kragic, and F. T. Pokorny,“An algorithm for calculating top-dimensional bounding chains,” PeerJ
Computer Science, vol. 4, May 2018.
A. Hatcher, Algebraic Topology. Cambridge University Press, 2002.
S. Lang, Fundamentals of Differential Geometry. New York, NY:
Springer US, 1999.
S. Calinon, F. Guenter, and A. Billard, “On learning, representing
and generalizing a task in a humanoid robot,” IEEE Transactions on
Systems, Man and Cybernetics, Part B, vol. 37, no. 2, 2007.
M. Robinson, “Sheaves are the canonical data structure for sensor
integration,” Information Fusion, vol. 36, 2017.