Journal name

Topological Trajectory Classification with
Filtrations of Simplicial Complexes andPersistent Homology

000(00):1–13©The Author(s) 2010
Reprints and permission:
sagepub.co.uk/journalsPermissions.nav
DOI:doi number
http://mms.sagepub.com

Florian T. Pokorny∗, Majd Hawasly† and Subramanian Ramamoorthy†
Abstract

1. Introduction

In this work, we present a sampling-based approach to tra-For robots to autonomously operate in a wide variety ofjectory classification which enables automated high-levelreasoning about topological classes of trajectories. Our
approach is applicable to general configuration spaces andrelies only on the availability of collision free samples.
Unlike previous sampling-based approaches in robotics
which use graphs to capture information about the pathconnectedness of a configuration space, we construct a multiscale approximation of neighborhoods of the collision free
configurations based on filtrations of simplicial complexes.
Our approach thereby extracts additional homological information which is essential for a topological trajectory classification. We propose a multiscale classification algorithm for
trajectories in configuration spaces of arbitrary dimension
and for sets of trajectories starting and ending in two fixed
points. Using a cone construction, we then generalize thisapproach to classify sets of trajectories even when trajectory
start and end points are allowed to vary in path-connected
subsets. We furthermore show how an augmented filtrationof simplicial complexes based on an arbitrary function on
the configuration space, such as a costmap, can be defined toincorporate additional constraints. We present an evaluationof our approach in 2, 3, 4 and 6 dimensional configuration
spaces in simulation and in real-world experiments using a
Baxter robot and motion capture data1 .

environments, we require algorithms and data structures thatenable such systems to reason rigorously about the collisionfree subset Cf of their configuration space. Furthermore, it isdesirable for such systems to be able to classify continuous
trajectories in their configuration space in order to reason
about the space of possible motions.
In the case where Cf is explicitly describable in a noise-freemanner in terms of semi-algebraic functions, analytic methods (Canny, 1988; Latombe, 1991; Schwartz and Sharir,
1983) provide an elegant classical avenue to studying Cf .
With the increase in computational resources and the availability of large data sets, the data-driven paradigm has
more recently provided an alternative to analytical methods; providing a more powerful and usable solution in many
scenarios where analytic models are either infeasible or computationally too expensive to obtain. An example of thisapproach is the use of Gaussian Mixture Models in the
Programming by Demonstration framework (Billard et al.,2008), where trajectories are modelled as integral curvesof vector fields obtained from a trained Gaussian Mixture
Model. Similarly, in motion planning, an explicit descriptionof Cf ⊂ Rd is often not available and popular algorithms
such as Rapidly-exploring Random Trees (RRT) and Probabilistic Roadmaps (PRM), (Kavraki et al., 1996; LaValle,
2006; LaValle and Kuffner, 2001), are based on the idea of∗ Centre for Autonomous Systems, CAS/CVAP,

KTH Royal Institute of Technology, Sweden,
fpokorny@kth.se† IPAB, School of Informatics, University of Edinburgh, UK
1 Parts of this work were presented in (Pokorny et al., 2014) at Robotics:
Science and Systems 2014. This work contains an extended discussion,
theory section and experimental evaluation.

utilizing a set of random samples X ⊂ Cf to construct agraph Γ with vertices in Cf and where edges correspond to
local paths which can be determined by a local path planner. The graph Γ can then be used to efficiently carry out

2

Journal name 000(00)motion planning. When Cf is a tame space, Γ, for sufficiently

The main benefit of working with homology is that homol-

large X, provides an approximation of Cf which allows us

ogy groups are purely linear algebraic objects that can beto answer basic questions about the path-connectivity of Cf .

computed efficiently while homotopy groups are generally

We observe that both, a large class of current probabilistic

complicated noncommutative groups. Our filtration F isapproaches such as Gaussian Mixture Models and Gaussian

based on Delaunay-Čech complexes which depend on a

Processes (Deisenroth et al., 2013; Rasmussen and Williams,

scale parameter r and which have recently been proven2006), as well as graph-based approaches in motion planningare currently not able to make use of homological informa-by Bauer and Edelsbrunner (2014) to provide a homotopyS
equivalent reconstruction of the space Xr = x∈X Br (x),tion about Cf . In particular, these methods are not able towhere Br (x) = {y ∈ Rd : kx − yk 6 r}. Our work utilizes

detect whether trajectories are homotopy inequivalent – i.e.

persistent homology as described in Carlsson (2009); Edels-

pairs of trajectories that cannot be continuously deformed

brunner and Harer (2008, 2010) which generalizes classicalinto each other. The reason for this is that, besides approx-

homology groups to a multiscale setting – meaning that we

imating the path-connectedness of Cf , unaugmented RRT

are able to compute topological information about the ana-

and PRM graphs do not store information about higher order

logue Fr of Γ for all scales r > 0 simultaneously without

topological information contained in higher homology andhaving to choose a particular scale upfront.

homotopy groups and, similarly, probabilistic methods haveso far not incorporated such information.
In this paper, we use filtrations F = {Fr : r > 0} of

Additionally, the 1-skeleton Fr1 ⊆ Fr is a graph which
can be used for path-planning. To summarize, our main
contributions are the following:

simplicial complexes defined in terms of random samples
X ⊂ Cf as a key technique for representing configuration spaces in practical robotics applications. Such filtrations•

cial complexes as well as persistent homology for mod-are common in topological data analysis (Carlsson, 2009),

elling robot configuration spaces and sampled trajectory

but have to the best of our knowledge not been utilized tomodel multi-joint robot configuration spaces prior to ourwork (Pokorny et al., 2014). We in particular propose a novel

data arising in robotics.
•

topologically. We then show how, using a cone construc-

simplicial complexes defined in terms of random samples

tion, our approach can classify trajectories whose start

X ⊂ Cf . From such filtrations, we then extract higher-orderand classifying equivalence classes of trajectories in Cf .

and end-points lie in connected sets S, T respectively.
•

jectories with the property that no trajectory belonging to
one equivalence class can be continuously deformed to any
trajectory in any of the other equivalence classes. Notethat there is a subtle difference to an exact classification
by homotopy classes in the statement above. Our classification is coarser in the following precise sense: we use
the first homology group for classification, which yields a
weaker classification than a classification by the first homotopy group. In particular, there exist cases where two trajectories which are considered equivalent in homology are
not homotopy equivalent, but two trajectories in different
homology classes will always be homotopy inequivalent.

We show how to use our approach when a general filtration function is considered, and we provide an example

Given a sufficiently good approximation of Cf by Fr , ourapproach yields a finite set of equivalence classes of tra-We propose an algorithm based on persistent homology
to classify trajectories with fixed start and end points

approach to trajectory classification based on filtrations of

topological information for the purpose of understanding

We introduce and evaluate the use of filtrations of simpli-of filtrations arising from a cost-function.
•

We provide an experimental evaluation of our approachin configuration spaces of dimension 2, 3, 4 and 6
using Delaunay-Čech complexes both in simulation, withmotion capture data of a human-robot interaction, andwith trajectories recorded using a Baxter robot.

The paper is structured as follows: In Sec. 2 we discuss the
motivation of our approach and the differences and similarities with prior works. Sec. 3 reviews the theoreticalbackground on persistent homology and filtrations. We discuss our approach to trajectory classification in Sec. 4 andpresent an experimental evaluation in Sec. 5. Finally, we
conclude with Sec. 6 and discuss future work.

3

Running head right side

Fig. 1. We display a rectangular configuration space C of side-length 500 and with two obstacles (in black). Several trajectories are
depicted in red and blue. The gray area displays an approximation of Cf by a Delaunay-Čech complex DCr (X) from 10000 samples
X ⊂ Cf and with filtration value r = 11.07 (left and middle figure) and r = 73.76 (right figure). All of our figures are best viewed incolor.

2. Motivation and Related WorkFor a robot to reason efficiently about trajectories within its
own free configuration space Cf , or about the motions ofother human or robot agents in its environment, a suitablepartitioning of continuously varying families of trajectories
into a discrete set of equivalence classes is desirable.
Clustering trajectories is difficult in general since trajectories can have varying length and are not immediately
representable as vectors in a vector space of fixed dimension
as required by commonly used algorithms such as support
vector machines (Cristianini and Shawe-Taylor, 2000).
Several approaches to the classification of trajectories, as
reviewed in (Zheng and Zhou, 2011), are based on various
approaches to measuring the dissimilarity between trajectories, such as the Hausdorff distance, edit distances anddynamic time warping. For the purpose of activity analysis,
the work of Morris and Trivedi (2009) reviews trajectory
clustering approaches based on various clustering algorithms and distance measures. These methods for clustering
are related to other approaches for simplifying trajectories
by taking into account topological attributes. For instance,
Katsikouli et al. (2014) present an approach to eliminate spurious features in trajectories, such as GPS traces of humanmobility data, in order to reduce the burden on further
analysis.
In robotics, the knowledge of classes of trajectories is beneficial for example in the programming by demonstration
framework (Billard et al., 2008) where movement primitivesof a robot’s behavior are constructed from initial trajectorydemonstrations provided by a human teacher. Equivalenceclasses of robot trajectories can furthermore be useful inorder to reason about alternative trajectories when a subsetof trajectories becomes invalid due to changing environment
conditions. The recent work of Knepper et al. (2012) has
demonstrated the usefulness of trajectory classes of local
paths to improve the efficiency of a motion planning algorithm in particular. Another reason to adopt an approachsuch as ours is the need to reason about task level attributesof robot motion, abstracting over, e.g., the variability implicitin the many ways a robot could traverse between two regionsin the configuration space. In Sisbot and Alami (2012), complex human-robot manipulation tasks are made tractable byreasoning at the level of object hand-over points and task
segments, for example. Our experiments, in section 5, illustrate how our classification procedure can be used to extractsuch qualitative movement primitives.
Purely topological approaches to the analysis of trajectories in Cf focus on notions of equivalence which do notdepend on a metric. There, two paths α, β : [0, 1] → Cf with
α(0) = β(0), α(1) = β(1) are called homotopy equivalentif α can continuously be deformed to β in Cf while keeping
the end-points fixed. We use Fig. 1 as a running example forillustration. Note that all figures in this paper are best viewed
in color. For now, consider Cf as being approximated by thegray region, while the black regions correspond to obstacles. Several trajectories with identical start and end-points
are depicted in red and blue. The blue trajectory in the leftmost figure is homotopy equivalent to the red trajectory,
while the trajectories in the middle figure are not homotopy
equivalent. Note that, while trajectories can have identical

4

Journal name 000(00)distance in R2 under e.g. the Hausdorff distance, they may or

us to study the homotopy equivalence of paths within the

may not be homotopy equivalent. When the lower obstacleneighborhood Xr of a set of samples X ⊂ Cf more gen-is removed in the right figure, the red and blue paths from

erally, e.g. when Xr does not yield a reconstruction of thethe middle figure become homotopy equivalent, for example.

full space Cf . In a more recent related work, Zhang et al.

While homotopies between selected trajectories can in some(2007) construct an approximate cell decomposition using acases be constructed explicitly, for example using a curve

recursively refined decomposition of Cf into hypercubes toshortening algorithm (Chou and Zhu, 2001), this approachensure a sufficiently fine reconstruction of Cf . However, only

can fail in several cases. Firstly, an explicit homotopy con-the path-connectivity of this decomposition is then used forstruction, can only establish the existence of a homotopy inmotion planning and homological properties are not further

some cases, but may fail, for example when there exist just ainvestigated.

single but multiple geodesics between two points. Further-

Our work is also related to sampling-based algorithms con-more, numerical stability and resolution pose challenging

structing a graph Γ from X to answer questions about the

problems for such methods, and the failure to construct apath-connectivity of Cf . RRTs and PRMs (Kavraki et al.,homotopy between two curves does not imply homotopy

1996; LaValle, 2006; LaValle and Kuffner, 2001), in partic-

inequivalence. Our approach instead does allow us to proveular, are examples of these which have attracted unabatedthe homotopy inequivalence of trajectories when the underlying configuration space is sufficiently well approximatedinterest since their invention (Karaman and Frazzoli, 2011;by a simplicial complex.

2007). The graph Γ can be thought of as an approxima-

Lindemann and LaValle, 2005; Masehian and Sedighizadeh,

Topological concepts such as retractions and cell decom-tion of Cf from X. The filtrations of simplicial complexes

positions have played a key role in classical approaches toused in our work extend the concept of a graph to a multi-

motion planning (Latombe, 1991). There, Cf is typicallyscale approach which can recover more detailed information

assumed to have a known algebraic or semi-algebraic struc-about Cf . Such filtrations depending on a scale parame-

ture. The visibility graph in 2D and retraction-based methods
rely on constructing a graph using which motion planning is

ter r have been used in topological data analysis (TDA)(Carlsson, 2009; Edelsbrunner and Harer, 2008, 2010) to

performed. The general roadmap method of Canny (1988)

study the persistent homology groups which we use here

uses ideas closely related to Morse theory and projections

and which capture information about the topology of data

to lower-dimensional spaces to obtain a complete motion

at all scales simultaneously. The origins of practical multi-

planner for semi-algebraic sets. Similarly, the seminal work

scale simplicial complex filtrations trace back to early workof Schwartz and Sharir (1983) proceeds by constructing an

by Edelsbrunner (1992); Edelsbrunner and Mücke (1994) on

exact cell decomposition by means of a cylindrical alge-

alpha-shapes and the introduction of persistence in (Edels-braic decomposition of Cf . This is related to our approachbrunner et al., 2002). In the field of sensor networks, persis-

since our simplicial complexes form a particular type oftence has also been utilized to solve multi-sensor coverageapproximate cell decomposition. In Schwartz and Sharir

problems (de Silva and Ghrist, 2006, 2007). In machine(1983), homology groups of Cf are computed by an exact

learning, persistent homology has been investigated as a toolcell decomposition Z and the general path planning prob-

for topology aware kernel density estimation, topologicallem is solved using Z. These classical works have to thefeature extraction and inference (Chazal et al., 2013; Fasy

best of our knowledge however not considered the use of

et al., 2013; Pokorny et al., 2012).

the first homology group of Cf for trajectory classification,

One of the advantages of the knowledge of homotopy

and the focus has been on motion planning and not classifi-

classes is that a motion planning algorithm can utilize effi-

cation. An important difference to our work is the fact thatcient replanning within each such class (Brock and Khatib,

we only assume the knowledge of potentially noisy point-

2000). Since local variational or gradient based methods cansamples from Cf using which we build a simplicial complex

continuously deform trajectories towards local optimality

filtration rather than assuming a known description of Cf

only within the homotopy class of an initial trajectory, it is

as a (semi-)algebraic set. Furthermore, our approach allows

advantageous to maintain a set of homotopy inequivalent

5

Running head right side
trajectories each of which can then be optimized using thesecan take into account many different aspects such as safety

methods. Topological information about path classes hence

distances and visibility (Mainprice et al., 2011) in particu-

allows us to incorporate non-trivial global information withlar. We show how our trajectory classification method canthese local methods.

be applied within this general costmap-based problem for-

Recent approaches which attempt to obtain equivalence

mulation, so as to obtain solution classes dependent both onclasses of paths include the works of Jaillet and Siméon

the topology of the configuration space and the cost function(2008a,b) on path deformation roadmaps where a graph-defined over it. A related line of work involves characterisingbased representation to plan in the space of paths up toenvironmental uncertainty probabilistically and then pos-a class of continuous deformations is proposed. Recently,

ing sampling-based motion planning in terms of paths that

researchers have in particular investigated homotopy classes

exceed a specified probability threshold in terms of robust-of trajectories in explicitly described spaces. Using theness, e.g. Aoude et al. (2013) and Luders et al. (2010). Our

residue theorem of complex analysis (Bhattacharya et al.,formulation of trajectory classification is also related to this2010) studied an application of homology classes to motion

approach since the sub or super-levelsets of a probability

planning in 2D in the case where the obstacles in C can bedensity function can be used with our approach.

contracted into representative and explicitly defined points.
In Bhattacharya et al. (2011), this was extended using electromagnetism theory and Ampère’s law to the case of 3D.
There, obstacles were assumed to be contracted into skeletons and then modeled as current-carrying wires. Similarly
to our work, the authors argue that homological information is useful and computationally favorable to moregeneral homotopy invariants in robotics. In Bhattacharya
et al. (2012), a generalization to arbitrary dimension is proposed and an integration of differential 1-forms over cycles
is shown to be sufficient to determine topological classes
using the language of de Rham cohomology theory. In Kim
et al. (2012), motion planning in 2D with homology constraints is formulated as a mixed-integer quadratic program
by endowing path segments with binary labels that identify
their relation to the domain obstacles. A problem the above
recent approaches suffer from is that they require an explicit
description of the obstacles in the configuration space, e.g.
in 2D as a union of shapes each of which is contractible
to a geometrically specified point p ∈ C − Cf . In manycases, such information is however not easily available forreal robotic systems or too expensive to compute. We instead
propose a data-driven, sampling-based approach to building
a representation of Cf from which topological information
about trajectories can be extracted.
Prior work such as Jaillet et al. (2010) addresses path planning under the consideration of a cost function defined overthe configuration space. In such work, a user-given cost function is defined over the configuration space as an additional
input to the standard path-planning problem. Cost functions

3. Theoretical BackgroundOur work builds on techniques from the field of algebraic
topology (Hatcher, 2002; Munkres, 1984) which studies theproperties of topological spaces by constructing algebraic
objects, such as the homology groups which we will discussshortly. Special cases of topological spaces are, for example,
smooth manifolds, arbitrary subsets of a vector space as well
as arbitrary metric spaces. Crucially, topological properties
are independent of metric information. In mathematics, algebraic topology has played a key role in the classificationsmooth manifolds. Homology groups in particular are homotopy invariants, meaning that homology groups do not vary
under continuous deformations, called homotopies, of the
topological space under consideration. A popular exampleof two spaces that are homotopy equivalent are a torus anda tea-cup which, when imagined as consisting of a rubberlike material can be stretched and bent (without tearing) todeform into one another. The homology groups of these two
spaces are identical and measure in particular the fact that
there exists a single ’hole/tunnel’ in both spaces and thateach space is connected.
In this work, we will not provide a full introduction toalgebraic topology, since this is not within the scope of thecurrent paper, but we will try to define and discuss some of
the key concepts – focussing mainly on providing those readers with a robotics background that might be unfamiliar withthe subject with an intuition for these techniques, as well as
providing sufficient detail for the readers to implement the6

Journal name 000(00)methods proposed here. For general persistent homology,
the works of Carlsson (2009) and Edelsbrunner and Harer(2008, 2010) provide an excellent introduction, while a comprehensive discussion of algebraic topology is provided in(Hatcher, 2002; Munkres, 1984).

3.1. Simplicial Complexes
A key notion in computational algebraic topology is thatof a simplicial complex which is used to model a large
class of topological spaces of interest. We recall several
basic definitions (see Edelsbrunner and Harer (2010)): k + 1
points v0 , . . . , vk ⊂ Rd are called affinely independent if
the vectors v1 − v0 , . . . , vk − v0 are linearly independent. A

Fig. 2. We depict an example simplicial complex K in the top left
figure. In the top right and bottom left figure, a 1-cycle c and c0 isshown in blue respectively. The two 1-cycles are equivalent in thefirst homology group H1 (K) since their difference is given by theshaded 2-boundary displayed in the bottom right figure.

k-simplex σ ⊂ Rd is a convex hull of k + 1 affinely independent points v0 , . . . , vk ∈ Rd , σ = Conv(v0 , . . . , vk )and we also denote σ by [v0 , . . . , vk ]. A face τ of a simplex
σ = [v0 , . . . , vk ] is a convex hull of a non-empty subset of{v0 , . . . , vk } and we denote this relationship by τ 6 σ. Notethat 2-simplices hence correspond to triangles, 1-simplices
to finite line segments, and 0-simplices to points in Rd .
Using unions of such simplices, we can construct interesting topological spaces by means of a simplicial complex: A
simplicial complex is a finite collection of simplices K suchthat σ ∈ K and τ 6 σ implies τ ∈ K, and σ, σ 0 ∈ K impliesσ ∩ σ is either empty of a face of both σ and σ . A subcom00plex L ⊂ K of a simplicial complex K is a subset L of thesimplices of K such that L is itself a simplicial complex.
One of the aims of this work, besides our trajectory classification algorithm, is the development and use of simplicial
complexes as a key data-structure in robotics, in particular as a powerful mechanism to accurately represent robot
configuration spaces.
In the top left part of Fig. 2, we illustrate an example of

3.2. Simplicial Homology
We consider simplicial homology over a field F for a simplicial complex K. In general, non-field coefficients, such
as Z can be considered, but for the discussion of persistenthomology later, we shall require field coefficients. With field
coefficients, the homology groups which we shall define
shortly are simply F-vector spaces that carry topologicalinformation. In our application, we will consider the binary
field F = Z2 = {0, 1}, in particular since Z2 coefficientscan be most efficiently implemented on a computer.
Pk
A p-chain c is a formal sum c = i=1 λi σi of p-simplices{σ1 , . . . , σk } in K with λi ∈ F and Cp (K) denotes the
F-vector space of all p-chains. In the case of F = Z2 , 1chains can simply be considered as finite sets of 1-simplices(edges) and 2-chains are finite sets of 2-simplices (triangles),in particular. For every p-simplex σ = [v0 , . . . , vp ] let ∂p σ
be the p − 1-chain formed by the formal sum of all p − 1
dimensional faces of σ as follows:a very simple simplicial complex in R2 consisting of six 2simplices (shaded triangles), fifteen 1-simplices (edges) andnine 0-simplices (the vertices). Note that simplicial com-∂p σ =pX(−1)i [v0 , . . . , vi−1 , v̂i , vi+1 , . . . , vp ],i=0

plexes are natural generalizations of graphs (which only

where v̂i indicates that vi is omitted. In the case of F = Z2 ,contain 0 and 1-simplices) as well as triangulations, both

−1 = 1 and the sign in the above sum can be ignored andof which are popular geometric data structures in robotics.

the boundary of a 1-simplex consists of the two end-points,

The union |K| of all simplices in a simplicial complex K is awhile the boundary of a 2-simplex consists of the threesubset of R which yields a topological space whose prop-

boundary edges. ∂p extends to a linear map ∂p : Cp (K) →derties can be studied using the concept of homology which

Cp−1 (K). A chain c ∈ Cp (K) such that c = ∂p+1 ω forwe discuss now.

some ω ∈ Cp+1 (K) is called a p-boundary, and we call

7

Running head right sidec a p-cycle if ∂p c = 0. The set of p-boundaries and p-

topological information which remains invariant under con-cycles is denoted by Bp (K) and Zp (K) respectively andtinuous deformations of |K| called homotopies (Hatcher,

Bp (K) ⊆ Zp (K) since ∂p ∂p+1 = 0. The quotient vector

2002).

space Hp (K) = Zp (K)/Bp (K) is called the pth homology
group of K. We denote the equivalence class of a p-cycle cin the quotient vector space Hp (K) by [c].
To illustrate these concepts, consider Fig. 2: The top rightfigure displays a 1-cycle c over the field Z2 in blue. Closed
loops correspond to 1-cycles because the each vertex isincluded an even number of times in a set of edges in thecycle. The bottom left figure corresponds to another 1-cycle
c0 . Observe that, when we add c + c0 in the Z2 -vector spaceof 1-chains, we obtain the blue closed curve γ surrounding
the shaded area in the bottom right figure. This curve γ is infact a 1-boundary, since it arises as the boundary γ = ∂2 ω ofthe shaded collection of triangles ω which form a 2-chain.
Note that it is easy to see why the interior edges of thesetriangles are not part of the boundary: they occur precisely
twice each and are hence zero modulo Z2 . The boundary
operator hence exactly assigns the geometric boundary to
our 2-chain ω. The same also holds more generally for 2chains over a general field F, when the correct orientation istaken into account (Hatcher, 2002).
With this intuition at hand, we finally observe that c −c = γ = ∂ω, so that c and c0 differ only by a 1-boundary.
0

Hence, in homology [c] = [c0 ], while [γ] = 0 in homology
since it arises as the boundary of the shaded triangles. We
furthermore have [c] 6= 0 because there does not exist any
2-cycle ν such that c = ∂2 ν. To summarize: two p-cycles
represent the same homology class if they differ by a pboundary, and a homology class [c] ∈ Hp (K) is non-trivial
precisely if c does not arise as the boundary of some p +1-chain.
We denote bp (K) = dim(Hp (K)), which is called thepth Betti number of K. In the example just considered, thefact that there exists only a single enclosed tunnel/void in|K| is measured by b1 (K) = 1. Similarly, b0 (K) is equal
to the number of connected components of K, so in our
example b0 (K) = 1 and higher Betti numbers are zero.
These Betti numbers capture crucial topological information
about the space |K| which do not depend on the particularsimplicial representation of |K| that was chosen. Furthermore it is a remarkable fact and one of the key theoremsof algebraic topology that these homology groups capture

3.3. Delaunay-Čech Complexes
In the preceding discussion, we assumed that a simplicial
complex K was provided to us. In this section we discussa particular family of simplicial complexes which can beconstructed from point-cloud data in Rd .
Consider a set of uniformly sampled points X={x1 , . . . , xn } ⊂ Y from a subset Y ⊆ R . The rSnneighborhood Xr = i=1 Br (xi ), where Br (xi ) = {x ∈d

Rd : kx − xi k 6 r} for r > 0, forms an interesting topological space. Unfortunately, Xr is not itself a simplicial
complex, but we can compute the homology of Xr by representing Xr by any simplicial complex Kr which is homotopy
equivalent to Xr . A well-known simplicial complex which
has this property is the Alpha complex Ar (X) of Edelsbrunner (1992) which can be defined in arbitrary dimension
and is a sub-complex of the Delaunay triangulation of X. In
this work, we shall instead use another recently studied simplicial complex, called the Delaunay-Čech complex (Bauerand Edelsbrunner, 2014) which is also a sub-complex of theDelaunay triangulation.
To proceed, let us first generalize the notion of a simplex
and simplicial complex to that of an abstract simplex andabstract simplicial complex: an abstract k-simplex σ is justa set of k + 1 elements, and we call k the dimension of σ.
An abstract simplicial complex K is a finite non-empty setof abstract simplices such that if σ ∈ K and ∅ =6 τ ⊆σ∈K, then τ ∈ K. If τ ⊆ σ ∈ K, τ is called a face of σ
(Edelsbrunner and Harer, 2008).
The Čech complex Cr with parameter r > 0 is an example of an abstract simplicial complex, defined by Cr (X) ={σ ⊆ X : ∩x∈σ Br (x) 6= ∅}. Note that Cr (X) has nodirect representation as a subset of Rd and hence does notdirectly form a geometric simplicial complex as previously
considered.
Given X, we can instead consider the complex D(X) ={σ ⊆ X : ∩x∈σ Vx 6= ∅} (the Delaunay triangulation ofX) where Vx denotes the Voronoi cell containing x. We
consider D(X) for points in X which are in general position, which occurs with probability one and which can also

8

Journal name 000(00)Fig. 3. We illustrate a point-cloud X ⊂ [0, 1]2 ⊂ R2 (left figure) as well as the resulting Delaunay-Čech complex DCr (X) at filtration
parameter r = 0, r = 0.05, r = 0.076, r = 0.2 (from left to right). For each filtration parameter value r > 0, DCr (X) is topologicallyequivalent (homotopy equivalent) to the union of balls Xr of radius r around the samples, which is depicted in blue. Using DCr (X)we can hence compute the homology groups of Xr at all scales.

be enforced by a small perturbation of X. Note here thatDCr2 (X) and consisting only of simplices of DCr (X) up

D(X), as defined above, is a simplicial complex containing

to dimension 2, directly without having to first compute the0-simplices up to d-simplices in R .

filtration values for the higher skeleta.

d

The Delaunay-Čech complex DCr (X), for r > 0, is
the subcomplex of D(X) defined by DCr (X) = {σ ∈D(X) : ∩x∈σ Br (x) 6= ∅}. The recent work Bauer and3.4. Filtrations

Edelsbrunner (2014) establishes that DCr (X) is homo-

Consider a simplicial complex K and a function f : K → R

topy equivalent to Xr , so that topological information about

defined on the simplices of K. If f satisfies f (τ ) 6 f (σ)

Xr can be extracted from DCr (X) directly. We define

whenever τ 6 σ, then we observe that Kr = f −1 ((−∞, r])f : D(X) → R by f (σ) = min{r : ∩x∈σ Br (x) 6= ∅},is a simplicial complex and Kr ⊆ Kr0 whenever r 6 r0 ,so that DCr (X) = f −1 ((−∞, r]) and DCr (X) changes

yielding a filtration of simplicial complexes with filtration

only at finitely many r1 < . . . < rm which can be com-parameter r. In the previous section, DCr (X), for r > 0,puted at all scales by determining f (σ) for each simplex
σ ∈ D(X). Every k-simplex σ = {v0 , . . . , vk } ∈ DCr (X)provides an example of such a filtration. Considering Fig. 3,
we observe that as r is increased, the topology of the union ofcorresponds to the geometric simplex given by the convex

balls space Xr changes. Initially there exists one connectedhull Conv(σ), so that 0-simplices are points, 1-simplices arecomponent per data-point. These components then merge

edges, and 2-simplices are triangles.
In Fig. 3, we illustrate a point-cloud X ⊂ R2 as well

as r increases and we can recover the fact that the point-

as the associated DCr (X) for several parameters r > 0.

third and fourth figure, where the first homology group of

The previously discussed Fig. 1 also displayed a Delaunay-

DCr (X) is one-dimensional. Persistent homology, whichČech complex constructed from sampled points. Note that,

we shall discuss now allows us to rigorously make statements

instead of DCr (X), we could also have considered theabout which homological features exist at various filtration

Alpha complexes Ar (X) of Edelsbrunner (1995) sinceparameters.

cloud ‘looks like there exists a hole in the middle’ in the

Ar (X) is also homotopy equivalent to Xr and furthermore
Ar (X) ⊆ DCr (X), and for sufficiently large parameter
R the complex AR (X) and DCR (X) reach their maximal

3.5. Persistent Homology

size, and AR (X) = DCR (X) = D(X). Since DCr (X)For a filtration of simplicial complexes, where f : K → R,

does not change for r > R, we shall use the notation

K is a finite simplicial complex, and Kr = f −1 ((−∞, r]),

DCR (X) = DC∞ (X). An advantage of DCr (X) over

we denote the finitely many filtration values at which Kr

Ar (X) in our application is that we can compute the fil-changes by r1 < . . . < rm . The inclusion αij : Kri →tration values for the 2-skeleton of DCr (X), denoted by

Krj , for i 6 j, induces a linear map hji : Hp (Kri ) →Hp (Krj ). We say that a homology class α ∈ Hp (Kri ) is

9Running head right side

0.07

200
0.403
150

0.05

0.2

100

0.02
50

0

0.02

0.05

0

0.076

0.2

0.403

Fig. 4. The left figure displays the inessential intervals of the zeroth
persistence diagram for the point-cloud displayed in the left part of
Fig. 3. Each point is of the form (0, ri ) and corresponds to a connected component that merges with another component at filtration
value ri . Besides the displayed points, there exists an additional
essential point (0, ∞) corresponding to the final connected component persisting for all filtration values. Note in particular that, forr > 0.07, only this essential component survives when the red persistence interval (0, 0.07), marked in red, has died and DCr (X),for r > 0.07 has only a single connected component corresponding
to (0, ∞). The right figure displays the first persistence diagram.
Note how the fact that a single hole exists in the data corresponds toa single persistence interval (0.076, 0.403), with large persistence
marked in red. Observe also how, in Fig. 3, at r = 0.076, the holeis first enclosed by Xr and by simplices. Similarly, at 0.403 thishole is filled for the first time.

born at ri if α ∈/ im(hii−1 ). A class α ∈ Hp (Kri ) born atj
ri is said to die at rj if hj−1(α) ∈/ im(hj−1ii−1 ), but hi (α) ∈im(hji−1 ). The difference rj − ri is called the persistence

0

0

50

100

150

200

Fig. 5. A reconstruction of Cf with two obstacles in black and from1000 collision free samples (red points) on a square of side-length500. DC25 (X) is displayed in the top-left yielding a good approximation to Cf . The first persistence diagram for DCr (X) is shownin the top-right. The two marked red points p1 = (10.58, 74.0),
p2 = (12.97, 90.38) with large persistence correspond to the birthand death filtration of the two holes in Cf . The bottom row displays
DC10.58 (X), DC12.97 (X) and DC74.0 (X) which correspond tothe birth of the smaller and larger hole (the first time they are
enclosed by edges), and finally to the death filtration value of thesmaller hole (the hole is covered at r = 74.0).

In Fig. 4, we display the zeroth and first persistence diagram corresponding to the point-cloud shown in Fig. 3,
allowing us to understand how the number of connectedof α: it measures how long a homological feature survives

components and holes in the point-cloud evolve over all fil-in the filtration. Classes born at ri which do not die are

tration values. The fact that a large hole exists in the middle

associated to (ri , ∞) and are called essential. The remainingof this point-cloud is recovered by the existence of the single

classes are called inessential. Similarly, if a cycle represents

large persistence interval in the first persistence diagram. In

an essential (inessential) class, we call the cycle essential

Fig. 5, we display the first persistence diagram and DCr (X)(inessential). For i 6 j, the p-th persistent homology groupfor DCr (X) arising from collision free samples in a sim-

is defined as Hpi,j = Zp (Kri )/(Bp (Krj ) ∩ Zp (Kri )). Non-ple configuration space with two obstacles shown in black.

trivial elements of Hpi,j correspond to equivalence classes

Observe that the two obstacles correspond to the two redof p-cycles born at or before ri and which persist, i.e. dopoints in the diagram which are far from the diagonal. The

not die in the filtration for r ∈ [ri , rj ). For i = j, thisremaining points correspond to holes which are due to noise

recovers the usual notion of homology Hpi,i = Hp (Kri ) =and which do not persist for a large filtration interval.

Zp (Kri )/Bp (Kri ). A graphical representation is obtainedby the p-th persistence diagram which associates (ri , rj ) toclasses born at ri and dying at rj and (ri , ∞) to essential

3.6. Computation via Matrix Reduction

classes born at ri (with multiplicity). The number of points

To compute the persistence diagrams of a filtration Kr1 ⊂in (−∞, ri ] × (rj , ∞] equals dim(Hpi,j ) and the verticalKr2 ⊂ . . . ⊂ Krm , it is convenient to refine the filtration

distance of a point to the diagonal indicates how long theas follows: we pick an ordering σ1 , . . . , σn of the simplices

feature persists (see Edelsbrunner and Harer (2008)).

of Krm such that, for all i ∈ {1, . . . , n}, Ki = ∪il=1 σl is a

10

Journal name 000(00)simplicial complex and there exist indices 1 6 i1 < i2 <3.7. H1 (Y ) and Homotopy Classes of Trajectories. . . < im = n such that Kij = Krj . Such a simplexwise

The final piece of background work we require is the connec-

filtration can be obtained by inserting simplices in Kri before
simplices in Krj if i < j and by inserting the faces τ ⊂ σof any simplex σ before inserting σ itself (Edelsbrunner andHarer, 2008).
Sn
Let K = i=1 σi be such a simplexwise filtration. The
boundary operator ∂ : ⊕dp=0 Cp (K) → ⊕dp=0 Cp (K) is a linear map which we express in the ordered basis σ1 , . . . , σnyielding an n × n matrix D with Z2 entries. For a matrix M ,we denote by Mj the j th column and by Mij the (i, j)-entry.
Note that D is upper triangular and Dij = 1 if σi is a codimension 1 face of σj . We let low(Mj ) = max{i : Mij 6= 0}
if Mj 6= 0 and low(Mj ) is undefined otherwise. A left-toright column addition Mj ← Mj + Mi , i < j is called
reducing if it decreases low(Mj ) and M is called reduced if
no reducing left-to-right column addition can be performedon any of its columns. The standard persistence algorithm
Edelsbrunner and Harer (2010) applies left-to-right column
additions to D until D is reduced, yielding a reduced matrix
R. Note that this reduction approach is very similar tothe standard Gaussian elimination algorithm, although performed over Z2 and not including all reduction steps. The
algorithm can further be generalized to F coefficients for an
arbitrary field F.
We can keep track of the column additions by initializing
the algorithm with R = D, V = In , so that R = DV .
For each left-to-right column addition Rj ← Rj + Ri for
i < j, we perform the column addition Vj ← Vj + Vi .
This algorithm terminates when R is reduced and we have
R = DV , where V is the matrix relating R to its unreduced version D. One defines (Chen and Kerber, 2011)
P = {(i, j) : Rj 6= 0 and i = low(Rj )}, E = {i : Ri =0 and low(Rj ) 6= i for all j ∈ {1, . . . , n}}. Returning to
Kr = f −1 ((−∞, r]), each (i, j) ∈ P with dim(σi ) = p
corresponds to (f (σi ), f (σj )) in the p-th persistence diagram and is generated by the p-cycle Rj which dies with
the introduction of the simplex σj . Similarly, each i ∈ Ewith dim(σi ) = p corresponds to (f (σi ), ∞) and the pcycle Vi which is still alive in the final filtration Kn = Krm .
Note that the cycles Vi , Rj do not correspond to canonicalchoices, but the persistence diagrams determine the ranks ofall persistent homology groups.

tion between the first homology group and homotopy classesof paths in a topological space Y . The obvious case to keep
in mind is Y = Cf ⊂ Rd . Recall that the first fundamentalgroup π1 (Y, x0 ) (Edelsbrunner and Harer, 2010) is a wellknown group whose elements consist of equivalence classesof closed continuous curves through x0 ∈ Y and lying
entirely in Y . Two closed paths α, β : [0, 1] → Y throughx0 lie in the same equivalence class if there exists a homotopy (i.e. a continuous deformation) between them which isconstant at the base-point x0 . When Y is path-connected,
π1 (Y, x0 ) is independent of the chosen base-point x0 andhence often denoted simply by π1 (Y ). Furthermore, if thespaces Y, Y 0 are homotopy equivalent spaces, π1 (Y ) andπ1 (Y 0 ) are isomorphic as groups. Two paths γ1 , γ2 in Y withthe same start point x and end point y can be deformed into
each other via a homotopy if the closed curve γ following γ1
from x to y and then γ2 from y to x is trivial in π1 (Y ). Hence,
π1 (Y ) is a natural group to consider for the purpose of trajectory classification. Unfortunately, to the best of our knowledge, currently no sufficiently efficient method for generalconfiguration spaces exists to compute the group structureof π1 (Y ) (e.g. for large scale simplicial complexes), since
π1 (Y ) can be complicated and non-commutative. To extract
topological information about homotopy classes, we canturn to the first singular homology group H1 (Y ) with binary
Z2 = {0, 1} coefficients, yielding a vector space which canbe explicitly computed via simplicial homology when Y is
homotopy equivalent to a simplicial complex K. The closed
curve γ can be represented explicitly as a 1-cycle in a sufficiently fine subdivision of K when a deformation retraction
from Y to K is computable, and γ then corresponds to avector [γ] in H1 (Y ) ∼= H1 (K). Finally, [γ] 6= 0 implies
that γ1 and γ2 are not homotopy equivalent, allowing us todiscern homotopy classes of continuous paths. Note however that homology is a weaker concept than homotopy,
so [γ] = 0 ∈ H1 (Y ) does not imply that γ1 and γ2 are
homotopy equivalent. To gain somewhat more granularity,
one can further replace Z2 coefficients for example with Zp
coefficients for a large prime p. In this work we choose Z2
coefficients due to their computational advantages for large
simplicial complexes.

11

Running head right side

4. Methodology
We consider a configuration space C ⊂ Rd and the set Cf ⊆C of collision-free configurations. We do not assume that
we have an explicit description of Cf or C available, and wewould like to study homotopy classes of a set of trajectories
T = {γ1 , . . . , γk } ⊂ Cf with a fixed starting point x ∈ Cf
and end point y ∈ Cf . In order to classify the trajectories,
we shall exploit the connection between homotopy classes

lesser degree from the ‘curse of dimensionality’, especiallywhen approximations of Xr , such as Vietoris-Rips or Witness complexes, which we shall discuss in the following, areconsidered.
For a sample set X, let R be the minimal r > 0 such thatγi ⊂ Xr for all i ∈ {1, . . . , k}. Our approach in both cases
above will now be to study the homotopy classes of thesepaths in the topological spaces Xr ' DCr (X), for r > R.

and the first homology group which we just discussed. We
now consider two multiscale settings:

4.1. Trajectory Discretization

1) X is a sufficiently dense sample We assume that X ={x1 , . . . , xn } ⊂ Cf yields a sufficiently dense sample, forexample sampled via rejection sampling from the uniform
distribution on C, or via a randomized exploration of theconfiguration space. We can then ask about a likely approximation of Cf from X. Our working hypothesis is that theSfamily of spaces {Xr = x∈X Br (x) : r > 0} contain
good such estimates. If X was sampled uniformly and Cf

Fig. 6. We display two continuous trajectories in red and blue (solidlines) and their homotopy equivalent discretizations as sequencesof 1-simplices (as dashed lines).

is a smooth compact submanifold M ⊂ Rd , this intuition
is in fact well-founded due to the reconstruction theorem ofIn order to compute properties of a trajectory γ : [0, 1] →Niyogi et al. (2008) which guarantees that, for a sufficiently

Cf , we first need to represent γ by a homotopy equivalentdense sample set, Xr deformation retracts to the manifold M

path of edges (i.e. 1-simplices) in the simplicial complex

for appropriately chosen r. Using the previously introduced

used to model Xr . Figure 6 illustrates this in an exam-

Delaunay-Čech complex and the fact that DCr (X) is homo-ple. Starting with a collection of samples X formed by the

topy equivalent to Xr (Bauer and Edelsbrunner, 2014), we

centers of the indicated shaded gray balls, we consider thewill then compute homological information about Xr fromspace XR , where R is the radius of the shaded balls. For

DCr (X). This approach is currently applicable mostly inan arbitrary continuous trajectory in XR such as the smooth

lower dimensions, since the ‘curse of dimensionality’ phe-trajectories in red and blue in the figure, we require a homo-nomenon implies that an increasingly large sample set would

topy equivalent trajectory of 1-simplices in the simplicial

be required in higher dimensions to faithfully reconstruct Cf .

complex DCR (X) homotopy equivalent to XR . In the figure these are indicated by red and blue dashed trajectories

2) X ⊂ T Here, we assume only the availability of therespectively.

trajectories T . We then discretize each trajectory γi as a

A fast heuristic to obtain a discretization of a trajectory γ

piecewise linear curve and use the vertex positions of all theis to consider vi = γ(i/N ), for some large N ∈ N, to map

piecewise linear segments in T as our sample set X. We

vi to a closest 0-simplex vi0 ∈ DCR (X) and to then replace

study the homotopy classes of these trajectories within thethe path segment between vi , vi+1 by a shortest edge-path

topological spaces Xr which constitute an approximation0
between vi0 and vi+1
in the 1-skeleton of the simplicial com-of the r-neighborhoods around T . This then allows us toclassify trajectories within Xr . In this framework, holes canplex. Alternatively, one can attempt to construct an explicit

arise either due to obstacles in the configuration space (as infirst to a path contained in DCR (X) and then approximat-

the dense case), or due to the distribution of the trajectories ining γ by a homotopy equivalent sequence of 1-simplices on

Cf . We consider applications of this case in our experimentsa sufficiently fine subdivision of DCR (X). The Alpha complexes Ar (X) of Edelsbrunner (1995) are subcomplexes ofwith a Baxter robot. Observe that this approach suffers to a

deformation retraction from XR to DCR (X) mapping γ

12Journal name 000(00)DCr (X) for all r > 0 which are also homotopy equivalentlarge r any subset of X will be included in Vr (X), resulting

to Xr and onto which an explicit such deformation retrac-

in a complex that is infeasible to compute even for moder-tion from Xr has been described in Edelsbrunner (1995),

ately sized X. Furthermore, the Vietoris-Rips complex is notfor example. While the study of efficient and theoreticallyalways homotopy equivalent to Xr , but provides instead ansound homotopy equivalent trajectory discretizations could

approximation (Carlsson, 2009). If computational efficiency

be explored further, we will instead focus on the classifica-

is required at a cost of the accuracy of reconstruction, sparsetion problem here, assuming that each trajectory has been

approximations to Vr (X) can be considered, such as thediscretized as a path of edges in DCR (X).

Witness complexes discussed by Carlsson (2009).
To compute 2-skeleta of DCr (X), we first compute a

4.2. Computation of DCr (X) and alternative Delaunay triangulation using (Barber et al., 1996; CGAL,
2013). This results in d-simplices [v0 , . . . , vd ] from which
complexes
As previously outlined, our key motivation for using
Delaunay-Čech complexes in this work is the fact that
DCr (X) is homotopy equivalent to the natural union of rballs space Xr . A drawback of this approach is however thefact that we rely on the computation of a Delaunay triangulation of X ⊂ Rd , which is computationally challenging as
the dimension d increases. The worst case arises for n points
distributed along the curve γ(t) = (t, t2 , . . . , td ), called themoment curve, in which case the complexity of the Delaunay triangulation as measured by the number of simplices isO(ndd/2e ) (Amenta et al., 2007; McMullen, 1970), while, on
the other hand the complexity is O(n) with a constant factor
exponential in dimension for points that are uniformly distributed in a unit ball (Dwyer, 1991). Given these results andthe state of the art implementations provided in (Barber et al.,
1996; CGAL, 2013), Delaunay-based simplicial complexes
currently face a serious ‘curse of dimensionality’ challenge
as d is increased. In this work, we will demonstrate the feasibility of Delaunay-Čech complexes in dimensions up to 6,

we extract all 0-simplices vi , 1-simplices [vi , vj ], 0 6 i <j 6 d and 2-simplices [vi , vj , vk ], for 0 6 i < j < k 6 dwithout repetitions. The 0-simplices are given by vertices,and for the large simplicial complexes in our experiments,
we found that extracting these k-simplices (k = 1, 2) was
sufficiently efficiently achieved as follows: first, we extract
the k-faces (k = 1, 2) for each d-simplex, which resultsin multiple occurances of simplices that are faces of morethan one simplex. Next, we lexicographically sorted theseresulting simplices by vertex indices via quicksort and finally
removed multiple occurances of simplices. Alternatively, if
sufficient memory is available, a set data-structure could be
used for this purpose. In the final step, the Delaunay-Čech filtration parameter for each simplex has to be computed. Here,0-simplices are assigned filtration value 0, 1-simplices [v, v 0 ]0)are assigned d(v,vand for general k-simplices we can apply2

the MiniBall algorithm (see Edelsbrunner and Harer (2010))
to determine f (σ), which is particularly simple in the caseof 2-simplices.

but for very high dimensional data, other approximations for
Xr are required.
The Vietoris-Rips complex, in particular, is suitable forsuch data. The Vietoris-Rips complex Vr (X) for a set ofsamples X ⊂ Rd and with parameter r > 0 is an abstract
simplicial complex defined as follows:

4.3. Homological Trajectory Classification
Consider a set of edge-paths {α0 , . . . , αm } in DCR (X)starting and ending at 0-simplices s, t ∈ DCR (X) respecdef

tively. We consider the 1-cycle cα0 (αu ) = α0 + αu ∈Z1 (DCR (X)), where the addition is performed over Z2

Vr (X) = {σ ⊆ X : diam(σ) 6 2r},

and we have reduced all edge-paths modulo Z2 , so that

where diam(σ) denotes the supremum over the pairwise[cα0 (αu )] 6= [cα0 (αw )] ∈ H1i,j = H1i,j (DC(X)) implies

Euclidean distances between points in σ. The Vietoris-Rips[αu +αw ] 6= 0, so that αu , αw are not homotopy equivalent incomplex can be efficiently constructed in high dimensions

DCr (X), R 6 ri 6 r < rj , where r1 < . . . < rm denote

if the maximal r parameter is chosen sufficiently small

the critical filtration values at which DCr (X) changes. We(Zomorodian, 2010). Note however that, for sufficiently

hence have trajectory classes {[cα0 (α0 )], . . . , [cα0 (αm )]} ∈edges occurring an even number of time are removed. Now

13Running head right side
H1i,j and the class membership can be computed once wewhile V has full rank. Since low is injective on the set S n ,have determined a basis for H1i,j . Note that α0 corresponds

we order elements of S n (for p = 1) by their low value andto the zero vector 0 = [cα0 (α0 )] and, for Z2 coefficients,

we store low−1 = l as a map such that l(k) is the element

there can be up to 2 trajectory classes for fixed s, t ands ∈ S n with low(s) = k. For any cycle c ∈ Z1 (Ki ), wei, j when dim(H1i,j ) = k. We can now compute a basis for

can then trivially solve for the coefficients in the basis S n by

H1i,j :

iterating c ← c + l(low(c)). Each iteration reduces low(c)kLemma 4.1. Let K1 ⊂ . . . ⊂ Kn be a simplexwise filtrationof simplicial complexes, let R = DV denote the reduced
boundary matrix after applying the left-to-right reduction
algorithm, and let Ep ⊆ E, Pp ⊆ P denote those elements
corresponding to p-cycles only. For 1 6 i 6 n, a basis ofZp (Ki ) is given by S i = {Rt : (s, t) ∈ Pp , s 6 i} ∪ {Vs :s ∈ Ep , s 6 i}, and, for 1 6 i 6 j 6 n, the image of thesetuntil we arrive at the zero vector. In the ordered basis S n , cnthen has non-zero coefficients F (c) ∈ Z#S
exactly at those2basis elements s ∈ S n for which low(s) = low(c) during
the execution of the above loop. Again, n can be very large(millions), but the vector F (c) is in our experiments very
sparse so that the algorithm does not exhibit its worst castnO(n2 ) computation time. We call F (cα0 (αu )) ∈ Z#S
the2
persistent cycle coordinates of αu with respect to α0 .
If we want to determine a trajectory class at scales corre-

T i,j ={Rt : (s, t) ∈ Pp , s 6 i, t > j}
∪ {Vs : s ∈ Ep , s 6 i}sponding to filtration values ri < rj , we select the coordinates F i,j (cα0 (αu )) of F (cα0 (αu )) corresponding to thebasis T i,j . Two trajectories αu , αw are then not homotopy equivalent if F i,j (cα0 (αu )) 6= F i,j (cα0 (αw ). Eachunder the quotient map

non-zero coordinate of F (cα0 (αu )) corresponds to a col-q : Zp (Ki )  Hpi,j = Zp (Ki )/(Bp (Kj ) ∩ Zp (Ki ))umn Rt of R which has a death filtration value f (σt ).

forms a basis of Hpi,j . Finally #Ep = dim(Hp (Kn )).

have been born and have not died yet contribute to the

Proof. This follows from the reduction algorithm (Edels-classification of cycles. We hence obtain an agglomerative
clustering of trajectories lying in a common DCR (X) as

brunner and Harer, 2010). Please see the appendix for further
details.

At filtration value r, only those non-zero coordinates that

we increase the filtration value r > R. Finally, at rm ,DCrm (X) = DC∞ (X) = Conv(X) and all trajectories

In order to classify {α0 , . . . , αm }, we first select athen lie in the same class.

simplex-wise refinement {Ki }ni=1 of the filtration given by
DCr (X), r > 0. Next, we compute the Z2 coordinates ofIllustrationcα0 (αu ) for 0 6 u 6 m in the basis S once. To classify tra-Consider Fig. 1. The red trajectory corresponds to α0 , thenjectories at a scale given by the filtration value ri = f (σ), we
simply look up the binary coordinates of cα0 (αu ) restricted
to the basis elements T i,i ⊆ S n . Similarly, we can checkif two trajectories αu , αw are homotopy inequivalent for all
ri 6 r < rj by looking up whether the coordinates ofcα0 (αu ) and cα0 (αw ) differ in the basis T i,j ⊆ S n .
Note now that DCr (X) = D(X) for sufficiently large r,
where D(X) denotes the full Delaunay triangulation, andtwo blue trajectories in the left and middle figure represent
α1 , α2 respectively, and all trajectories lie in DCri (X), ri =11.07. We have [cα0 (α0 )] = [cα0 (α1 )] = 0 ∈ H1a,b for alli 6 a 6 b, but [cα0 (α2 )] 6= 0 ∈ H1a,b , for i 6 a 6 b 6 j,where rj = 73.76 is the critical filtration value at which thehole surrounded by α0 , α2 gets filled in.

H1 (D(X)) = {0} since D(X) is contractible. Hence E1

4.4. Arbitrary filtration functions

is empty implying that we do not need to keep track of the

Suppose now that we have sampled Cf sufficiently densely

to determine a basis of H1i,j . This is importantmatrix V
since, in our experiments, these matrices have up to millions

and that DCR (X), for some fixed R, provides a goodof columns and R is typically very sparse and of low rank,

Our aim now is not only to classify trajectories in theapproximation of Cf . Consider a function c : Cf → R.

14

Journal name 000(00)space Cf , but to take into account a threshold for the function c. Naturally, such c : Cf → R can arise as a costor risk function associating a certain cost or risk to eachx ∈ X. Similarly, one might have defined a probabilistic model associating a certain likelihood or density value
to each x ∈ X. Given c, we define the value ĉ of a ksimplex σ = {v0 , . . . , vk } ∈ DCR (X) to be ĉ(σ) =max(c(v0 ), . . . , c(vk )). Then ĉ satisfies ĉ(τ ) 6 ĉ(σ) whenever τ ⊆ σ. In order to apply our algorithm without having
to keep track of the potentially non-sparse matrix V , wefurthermore let ĉ(σ) = max(c(v0 ), . . . , c(vk )) + C for any
k-simplex σ ∈ D(X) − DCR (X) and for C larger thanthe cost of any v ∈ X. Then LR,λ = ĉ−1 ((−∞, λ]) yields
a filtration as λ varies and LR,λ ⊆ DCR (X) for λ < C
and LR,∞ = D(X), ensuring H1 (LR,∞ ) = 0. We think ofLR,λ as an approximation to MR,λ = XR ∩ c−1 ((−∞, λ])for λ < C. Note that this approach yields a filtration by

Fig. 7. Examples of trajectory classes obtained using the cone construction for samples in [0, 1]2 ⊂ R2 . All the trajectories depicted
here start in the red part of the complex and terminate in the blue
part of the complex. The trajectories in the left part of the figure
belong to a single topological class and the trajectories to the rightsimilarly belong to a single topological class at the depicted filtration value of R = 0.1. The subcomplexes corresponding to S
and T are shaded in red and blue respectively. The square has aside-length one and we used 3000 sample points X. The depicted[R (X) whichsimplicial complex shows only the simplices of DC
do not involve the special points s, t.

sub-levelsets. In the case of a probability density c, superlevelsets c−1 ([λ, ∞)) might be more appropriate. In thatcase, we simply replace c by −c in the above framework toyield the desired filtration.

now propose an alternative approach that can be applied in
all cases and which relies only on a cone construction.
Intuitively, we would like to work with a quotient space
where the regions S, T are identified with points s, t respectively. We hence define a simplicial complex filtration [DC

4.5. General start and goal regionswith two additional 0-simplices at filtration value zero,

Real-world trajectories are often noisy and varied and might

called s, t. We then reassign filtration value zero to any sim-

not naturally start and end in a fixed start point s andend-point t. In this section, we discuss an extension of
our approach which enables us to work with connected
regions S, T generalizing the points s, t respectively. We
hence assume that there exist connected regions S, T ⊂ Cf ,S ∩ T = ∅ such that, each of our trajectories starts in a pointin S and terminates in a point in T .

plex in the filtration whose vertices all lie either in S or in
Sd+1
T . We furthermore define simplices p=1 (S p ∪ T p ), whereS p = {{v0 , . . . , vp−1 , s} : σ = [v0 , . . . , vp−1 ] ∈ DC∞ (X),vj ∈ S for all j ∈ {0, . . . , p − 1}}
T p = {{v0 , . . . , vp−1 , t} : σ = [v0 , . . . , vp−1 ] ∈ DC∞ (X),vj ∈ T for all j ∈ {0, . . . , p − 1}}.

Note that, in the case where S, T correspond to subcomplexes with trivial first homology, the following pro-and also assign filtration value zero to these and we cre-cedure could be employed: we select a representative 0-

ate a simplex-wise filtration where we also insert faces ofsimplex s ∈ S and t ∈ T and extend each trajectory

simplices before their containing simplices whenever these{α0 , . . . , αm } to a trajectory from s to t by first determin-have identical filtration value. This yields a new filtration

ing an edge-path from s to αj (0) inside S, then following

where the two subcomplexes of DC∞ (X) corresponding to

αj and finally concatenating the resulting path by another

S, T and the newly added cones over these subcomplexesedge-path from αj (1) to t. The resulting augmented trajec-

appear at filtration value zero and additional simplices fromtories could then be classified with our previous approach.
In the case where H1 (S) 6= 0 or H1 (T ) 6= 0, this procedure

DCR (X) are only added in afterwards. Any simplex in[
DC R (X) which involves s, t does not correspond to a spe-will however depend on how the trajectories are augmented,

cific geometric simplex in Cf , but is abstract in nature. Noterequiring motion planning under homology constraints. We

that s corresponds to the apex of the cone over S and t to theRunning head right side15

apex of the cone over t. Since cones are contractible topological spaces, this yields an explicit identification of the
sub-complexes corresponding to S and T with the points s,
t respectively.
Each trajectory α can now first be discretized as a sequenceof edges in DCR (X) as before. Assume that α starts in pointa ∈ S ∩ X and terminates in b ∈ T ∩ X. By appending theedges (s, a) and (b, t) to each such trajectory, we obtain aset of lifted trajectories α̂0 , . . . α̂m in [DC R (X) which now
all start and end in the same (non-geometric) points s and t
respectively.
Since we will only be interested in the first homology
group, we can in fact ignore all simplices of dimension
larger than 2 in the above construction since these do not
influence the first homology group. We can then apply ourprevious algorithm on these augmented trajectories to obtaina topological classification. If two trajectories α̂i , α̂j areinequivalent at filtration R > 0 using this approach, thenthe corresponding trajectories αi , αj cannot be continuously
deformed into one another, even when the end-points areallowed to move within the sub-complexes corresponding to
S, T . Fig. 7 illustrates two examples of classes of trajectories
which are considered inequivalent using this approach. The
computational cost incurred by adding the cone construction
is dependent only the number of additional simplices in the
2-skeleton of [DC. In particular, we add two 0-simplices s, t,
and one edge for each point in S ∪ T as well as one triangle
for each edge lying either entirely in S or entirely in T .

Fig. 8. We display example worlds and examples of paths whichwere determined to lie in a single class (in blue) at filtration value r2 .
DCr2 (X) was constructed from 100000 samples and the classes
are computed using the indicated red trajectories corresponding to
α0 .

5. Experiments
Our experiments were performed on an Intel Core i7 lap-

for all but our Baxter experiments where we used QHull

top with 8GB of RAM. We present only the computation

Barber et al. (1996) which was faster in higher dimensions.

times of core algorithms and disregard the time required toload data into memory. We used the matrix and binary tree
column vector data structure of the PHAT library (Bauer

Trajectory classification in 2D We generated the setof 2D worlds W1 , . . . , W10 displayed in Fig 8 and of

et al., 2013) to efficiently manipulate large boundary matri-size 512 × 512 by sampling Gaussian Random Fields

ces. Instead of working with the full simplicial complex

and defining those regions above a threshold to be obsta-

DCr (X), we extracted the 2-skeleton DCr2 (X) from thecles. From the resulting free space Cf , we sampled

Delaunay triangulation D(X). The 2-skeleton is sufficientfor our purposes since H1i,j (DCr (X)) does not depend on

N ∈ {1000, 10000, 100000, 1000000} uniform samples.
We computed the Delaunay-Čech filtration for all examples

higher dimensional simplices. We reduced only the sub-and recorded the computation times of the Delaunay triangu-

matrix of the boundary matrix corresponding to the firstlation, for the construction of the filtration, as well as the time

homology group. D(X) was computed with CGAL (2013)

required to reduce the boundary matrix D to its reduced form

16

Journal name 000(00)

0.5
0.382
0.301

0

Fig. 9. We display the example world W1 with DCr2 (X) for 1000,
10000 and 100000 sample points per row. In each column, we plot
paths α1 , . . . , αs (in blue) which belong to a fixed trajectory class at
filtration value r2 . The fixed reference path α0 is plotted in red. As
expected, we can clearly see that two paths in different classes also
lie in different homotopy classes. In our experiments, paths withina class are furthermore homotopy equivalent in DCr2 (X), but thequality of the approximation DCr2 (X) ' Cf is only sufficientfor 10000 or more sample points as can be seen in the right figurein the first row. There, some 2-simplices (triangles) cover the thinobstacle region to the right.

00.5

Fig. 10. The top figure displays the robot arm in start configuration(blue) on the right and in goal configuration (red) on the left. The
bottom right figure displays the first persistence diagram for ourreconstruction with one red point far above the diagonal. A projection of the samples onto θ1 , θ2 is shown in the bottom left and an
illustration of the difference between the two trajectory classes forr ∈ [0.301, 0.382] is shown in the middle figure. In the first trajectory class (in red), the arm is extended to the left when passing
under the narrow passage while in the second class (in blue), thearm is extended to the right.

from v1 to wi and then to v2 utilizing Dijkstra’s algorithm on
the 1-skeleton graph of DCr1 (X). The computation times
for the persistent cycle coordinates for these paths were

R. The Delaunay triangulation took 1ms, 2ms, 76ms, 810ms,

1.8ms, 10ms, 115ms and 1.75s for a batch of 100 query paths

the construction of the filtration took 11ms, 31ms, 278ms and3.27s and the reduction of the boundary matrix took 14ms,

and for the respective sample sizes on average. These encouraging timings suggest that our framework could be used as

13ms, 76ms, 981ms on average as the sample size increased.

a classification ‘black box’ e.g. for continuous trajectory

We investigated the filtration DCr (X) at various thresholds.
p
At a filtration value of r1 = 25 1000/N , we found that Cfp
was conservatively covered, while at r2 = 35 1000/N , theoptimization engines.

space was well covered with a minimum number of holes inrobot arm displayed in the top part of Fig. 10 which iscollision free areas. In order to investigate interesting pathclasses, we generated a set of 1000 paths per world andattached to the central black disk and with 4 joints θ1 , . . . , θ4 .
We constrain θ1 ∈ [− π2 , π2 ], θ2 , θ3 , θ4 ∈ [−0.9π, 0.9π] andsample setting as follows: In 10 trials, we selected two sam-furthermore disallow self-collisions and collisions with theple points v1 , v2 at random and, for each such setting, we

environment (the black rectangle and the floor), yielding

selected another 100 random waypoints w1 , . . . , w100 fromCf ⊂ R4 . The robot now has the task of moving from thethe sampled point-cloud. We determined shortest edge-pathstart configuration displayed in blue to the red goal joint

Trajectory classification in 4D We consider the planar

17

Running head right side
configuration as shown in the top left figure. We sampled

at λ = 70 both hills are truncated in the remaining figures.

100000 poses uniformly in Cf using OpenRave (Diankov

We sampled 100 random paths in this configuration space2
and Kuffner, 2008) and applied our framework. DC∞(X)

by fixing the initial and terminal vertex at the start and end-

had about 6.2 million triangles and 1.8 million edges. The

point of the drawn red reference trajectory and by sampling

bottom right part of Fig. 10 displays the resulting first persis-

random waypoints as before. The figure displays example

tence diagram which clearly shows that a single homologicaltrajectory classes for differing cost filtration values. Notefeature has large persistence in Cf . The projection of thehow, at a cost threshold of λ = 90 in the right plot, the two

joint configurations onto the first two angles, as shown in theclasses depicted in the two leftmost parts of the figure merge.

middle figure, confirms the existence of a single hole. We
computed 1000 edge-paths in DC0.25 between the start andend-configuration using 1000 random waypoints as before.
For filtration values r ∈ [0.301, 0.382] only two trajectory
classes existed. The reduction of the boundary matrix took
0.46s, while the persistent cycle coordinates for all 1000
paths were calculated in 0.55s. The Delaunay triangulationin R4 took 251s, partially due to the increased dimension.
Note however that these results are not directly comparable
to the 2D case since methods for 2D Delaunay triangulations
in CGAL (2013) are especially optimized. We inspected the
trajectories in each homology class and found that they wereclassified according to whether the second link was positioned to the left or to the right of the base link of the arm
when θ1 = 0 as the arm passed the narrow passage (seeFig. 10). Our framework hence allows the robot to discover
the fact that two fundamentally different solution trajectory
classes exist. Note that it is also intuitively clear that no
trajectory in the red class can be continuously deformed to
any trajectory in the blue class since the second link would
have to collide with the marked black obstacle above thearm at some point in the course of any such deformation.
Our method hence allows us to automatically recover thisinformation.

Baxter robot, 3D and 6D We now investigate a kinesthetic
demonstration scenario where the Baxter robot in Fig. 12 istaught a set of trajectories which we then classify topologically. In the 1st experiment (E1, Fig. 12, column 1-2), therobot is shown two ways to reach from one point above its
head to a point in front of its torso. Only one arm is moved
in each demonstration while the other arm remains still. We
used the end-effector positions of the moving arm to represent trajectories in R3 . In the 2nd (column 3-4) and 3rd(column 5-6) experiment E2 and E3, we record the positions of both end-effectors during dual arm manipulations
resulting in a 6D configuration space. In E2, the robot istaught to pick up a cylindrical object with both hands from atable and to move it to one of two positions, one higher
and one lower than the table. We also vary the distance
between the hands during grasping between demonstrations.
The trajectories in E2 are periodic. The motions start withthe arms in a rest position on the sides, the object is thengrasped and moved, and the rest position is visited again.
In experiment E3, the robot moves the same object froma horizontal to a vertical configuration, but a metal bar
is located between the robot and the object. Two intuitive
motion classes are based on whether the left arm crosses
in front of the obstacle, or behind it. Note that, in exper-

Filtrations with cost functions We consider the free con-iment E1 and E2, no obvious obstacles lie directly in Cf ,figuration space Cf ⊂ R of size 250 by 500 with two

but due to the type of demonstrations and the robot’s joint

obstacles (in white) displayed in Fig. 11. We would now like

limits, the space Xr ' DCr (X) ⊂ Cf exhibits interest-

to distinguish not only between homotopy classes depending

ing voids which we can exploit for classification. We foundon the obstacles in the configuration space, but also discernthat trajectories were well-approximated using the described

how trajectories behave with respect to the two peaks of theheuristic mapping to nearby edge-paths in DCR (X) (Sec. 4)

cost function. The simplicial complex L10,λ (X) is displayedfor R = 0.08m (experiment E1) and R = 0.15m (experi-

for 10000 samples X and height values are determined byment E2 and E3) respectively. For smaller r, DCr (X) was

the cost function. At cost threshold λ = 90, the top of one

either not path-connected, or the edge path approximationof the hills defined by the cost-function is removed from thecomplex in the rightmost figure (indicated in blue), while

deteriorated significantly. We hence investigated classifying

2

paths in Xr for r > R. Fig. 13 displays how the number of18

Journal name 000(00)Fig. 11. We display a cost function and classes of trajectories (in blue) depending on a cost threshold and a path α0 (in red). For thehigher threshold in the rightmost plot, the two classes in the two leftmost figures merge.

Fig. 12. Baxter robot experiments (E1: column 1-2, E2: column 3-4, E3: column 5-6). Trajectory classes are illustrated in the top
row and details are provided in the bottom row. We recorded 97, 18 and 40 trajectories resulting in point clouds with 9594, 3650,
4326 points in dimension 3, 6, and 6, and in filtrations with 0.19, 2.72, 3.05 million edges and triangles for E1, E2, E3 respectively. The computation times for (Delaunay triangulation, 2-skeleton and filtration computation, boundary matrix reduction) were(0.22, 0.26, 0.12), (41.84, 46.95, 44.17) and (44.92, 51.99, 52.44) seconds for E1, E2, E3 respectively, while the classification of all
trajectories in each experiment took no more than 0.05s. Bottom row: The first two images show the end-effector trajectories in E1.
While we obtain 3 classes for r = 0.08m (first image) the green class merges with the blue one at r ≈ 0.084m (second image, seeFig 13). The 3rd , 4th and 5th , 6th image show the right and left hand trajectories for E2 and E3 respectively. The trajectory classes at
r = 0.19m are indicated in red and blue.

topological trajectory classes changes with varying r, and theCone construction and human-robot collaboration In

second row of Fig. 12 displays trajectory classes at various

this experiment, we demonstrate the cone construction withfiltration values. In all experiments, there exists a large filtra-

real world data and provide initial experiments explor-tion interval with just two trajectory classes corresponding

ing how our approach could be incorporated in a human-to the two intuitive classes we just described for each exper-

robot collaboration scenario. In particular, we illustrateiment. The first image in the second row of Fig. 12 also

how, by topologically clustering trajectories, a robot canillustrates the three classes one obtains in E1 for a choice ofautomatically extract high-level motion behaviors.

r = 0.08m. More details on the computation times and the

We consider the top part of Fig. 14. A PR2 robot is posi-

sizes of the simplicial complexes is provided in Fig. 12.

tioned at the depicted position in the photo and does not
move, while a human stands at the green disk marked on

3
2
1

0.08

0.25

3
2
1

0.15

0.3

3
2
1

the floor. Using a Nest of Birds magnetic motion tracker,
0.15

0.19

0.25

Fig. 13. Number of path classes (vertical axis) vs. filtration value(horizontal axis) for experiments E1, E2, E3 from left to right.
All paths exist starting at filtration value R = 0.08m for E1 and
R = 0.15m for E2 and E3.

we recorded the position of an object that is initially placedin the PR2’s right gripper and which is then handed over
to a human worker who then places it on the blue table byturning his upper body while keeping his feet on the shaded
green disk on the floor. In the first experiment, we recorded

19

Running head right sidePR2 handhuman

PR2 handred table

1

red t

-0.5
-1

0.5

able

Fig. 15. Top and 3D-view of our second set of recorded trajectories
between the robot’s endeffector and the work surface marked in
red in Fig. 14. The depicted classification was obtained at filtration
R ≈ 15cm and uses the cone construction.
blue table

PR2 handpart of Fig. 14 (i.e. to perform an additional work step), fol-

red table

1

-1

d

-0.5

blue table

re

PR2
hand

ta

bl

e

lowed by placing the object on the table marked in blue.
Our method is able to identify these outliers as belonging toa seperate class (red) from the motions that move directly
towards the blue table region (blue trajectories). Note also0Fig. 14. The top figure displays a PR2 robot in an environment with
two tables. The marked red and blue regions simulate two target
work surfaces. A human is positioned to stand on the green marked
disk on the floor. We recorded human-robot hand-over trajectories
simulating a human worker collaborating with the robot. The 3D
position of an object placed in the PR2’s right endeffector is tracked
and the human worker takes the object and places it on the blue work
surface, but in some cases first moves the object to the red work
surface for an intermediate simulated ‘assembly step’. The bottom
part of the figure illustrates an example topological classification ofsuch trajectories at filtration R ≈ 16cm. The found red trajectories
correspond to sequences where the object is first placed on the
red work surface, while blue trajectories correspond to sequences
without this intermediate step. The bottom left figure is displayed
from a bird’s eye viewpoint, while the bottom right figure displaysa viewpoints from behind the blue work surface.

how the cone construction allows us to identify the wide
variations of terminal positions of the blue trajectories.
In a second experiment, we recorded 18 trajectories resulting in 4228 points in R3 (and 45895 edges and 83334
triangles in [DC). This time, the object is picked up from thePR2’s hand and placed on the red table. The human turns
either left or right by 180 degrees during this motion anddoes not leave the marked green spot in Fig. 14. The resulting trajectories are depicted from two perspectives in Fig. 15.
Again, we apply the cone construction with points marked
in green and cyan corresponding to S and T respectively.
At filtration value R ≈ 15cm, we obtain the 4 trajectory
classes indicated in color. At R ≈ 20cm the red and blueclasses merge and at R ≈ 23cm the green and black classes
merge and only two classes exist. One corresponding to the

7 trajectories of the object being handed over and placed

black and green trajectories and one corresponding to theon the blue table, resulting in 2486 points in R . The bot-blue and red trajectories. These finally merge into one class

tom part of Fig. 14 indicates the trajectories. Green shaded

at R ≈ 35cm. The two classes found in the significant filtra-

points correspond to points which we classified as lying intion interval (23cm, 35cm) exactly correspond to whether

the start region S in the vicinity of the robot’s right gripper.

the worker turned to the left or to the right after obtaining the

The shaded cyan points correspond to goal positions T in

object. The green trajectories correspond to a simulated casethe vicinity of the blue table-region.

where the object was dropped or placed on the ground before
being put on the table marked in red. Again, our method3

We applied our cone construction and classified the trajectories. For a large filtration interval between R ≈ 16cm

is able to identify such interesting trajectories as well as

and R ≈ 21cm, we obtained two classes as indicated in

the high-level information about the turning behavior in thiscolor the figure. Here, the red trajectories imitate a worker

trajectory dataset. We believe that such classification infor-

first placing the object on the table marked in red in the top

mation could serve as input to further algorithms to allow the

20

Journal name 000(00)robot to reason about the behaviors of the observed humanrobotics applications. In future work, we plan to further

worker, for example in order to optimally assist the humaninvestigate the potential applications of these data-structures

worker or to automatically extract clusters of trajectories

and techniques in robotics.

based on which dynamic motion primitives can be learned.

Acknowledgments
6. Conclusion and Future Work

This work was supported by the EU projects TOMSY (IST-In this work, we have proposed a novel sampling-based

FP7-270436) and TOPOSYS (ICT-FP7-318493).

approach to studying topological classes of trajectories
in general configuration spaces which utilizes persistenthomology and filtrations of simplicial complexes as key

References

ingredients. We have shown the applicability of our approach N. Amenta, D. Attali, and O. Devillers. Complexity of delaunay
in real world scenarios and believe that our approach could
furthermore be incorporated with many existing algorithms.

triangulation for points on lower-dimensional polyhedra. In Pro-For example, the integration of local trajectory optimiza-

Discrete algorithms, pages 1106–1113. Society for Industrial

tion based algorithms with our approach which extracts

and Applied Mathematics, 2007.

ceedings of the eighteenth annual ACM-SIAM symposium on

global information about trajectories could be of interest. G. S. Aoude, B. D. Luders, J. M. Joseph, N. Roy, and J. P. How.
Another promising future application of our method could be

Probabilistically safe motion planning to avoid dynamic obsta-

a class-dependent generation of dynamic motion primitives.

cles with uncertain motion patterns. Autonomous Robots, 35:

While our method scaled well to large sample sizes in 2 to 4

51–76, 2013.

dimensions and was applicable also for a smaller set of sam- C. B. Barber, D. P. Dobkin, and H. Huhdanpaa. The quickhull
ples in 6D, Delaunay triangulations and hence DCr (X) havealgorithm for convex hulls. ACM Trans. Math. Software, 22(4),a worst-case complexity of O(n

1996.

dd/2e

) (McMullen, 1970)

in dimension d and sample size n making triangulations U. Bauer and H. Edelsbrunner. The morse theory of Čech and Delauin higher dimensions computationally very challenging. It

nay filtrations. In Proc. of the Thirtieth Annual Symp. on Comp.

remains to be seen where the boundaries of feasibility ofGeometry, SOCG’14, pages 484:484–484:490, New York, NY,

simplicial complex based methods lie exactly in this respect

USA, 2014. ACM.

in robotics. The generalization of computer graphics based U. Bauer, M. Kerber, and J. Reininghaus. PHAT (Persistent Homolworks on mesh-compression such as Blandford et al. (2005)

ogy Algorithm Toolbox), 2013. http://code.google.com/p/phat/.

to higher dimensions might enable us to scale the pre- S. Bhattacharya, V. Kumar, and M. Likhachev. Search-based pathsented techniques to very large complexes by reducing theplanning with homotopy class constraints. In Proc. of The

required memory footprint, for example. Another direction

Twenty-Fourth AAAI Conf. on Artificial Intelligence, Atlanta,of future research will be an investigation of alternative

Georgia, 11-15 July 2010.

simplicial complexes, such as Vietoris-Rips and (weak) wit- S. Bhattacharya, M. Likhachev, and V. Kumar. Identification andness complexes (Carlsson, 2009) which approximate Xr

representation of homotopy classes of trajectories for search-

and which could be used to mitigate the ‘curse of dimen-

based path planning in 3D. In Proc. of Robotics: Science andsionality’ in higher dimensions to some extent. Another

Systems, 27-30 June 2011.

potential direction could be the use of projections or non- S. Bhattacharya, D. Lipsky, R. Ghrist, and V. Kumar. Invariants
linear dimensionality reduction techniques to accommodate

for homology classes with application to optimal search andhigher dimensional data.
Our work has established that filtrations of simplicial

planning problem in robotics. Electronic pre-print, Aug 2012.
arXiv:1208.0573.

complexes can provide a rich alternative data-structure for A. Billard, S. Calinon, R. Dillmann, and S. Schaal. Robot programming by demonstration. In Springer handbook of robotics, pages
robotics which is capable to generalize the notion of a graphto extract novel topological information for the purpose of1371–1394. Springer, 2008.

21

Running head right side
D. K. Blandford, G. E. Blelloch, D. E. Cardoze, and C. Kadow.
Compact representations of simplicial meshes in two and three
dimensions. International journal of computational geometry &and Comp. Geometry, 13(1):415–440, 1995.
H. Edelsbrunner and J. Harer.

Persistent homology-a survey.

Contemporary mathematics, 453:257–282, 2008.
H. Edelsbrunner and J. L. Harer. Computational topology: anapplications, 15(01):3–24, 2005.
O. Brock and O. Khatib. Real-time re-planning in high-dimensional
configuration spaces using sets of homotopic paths. In Proc.
of the IEEE Int. Conf. on Robotics and Automation (ICRA’00),
2000., volume 1, pages 550–555. IEEE, 2000.
J. Canny. The complexity of robot motion planning. MIT press,

introduction. AMS Bookstore, 2010.
H. Edelsbrunner and E. P. Mücke. Three-dimensional alpha shapes.
ACM Transactions on Graphics (TOG), 13(1):43–72, 1994.
H. Edelsbrunner, D. Letscher, and A. Zomorodian. Topological persistence and simplification. Discrete and Computational
Geometry, 28(4):511–533, 2002.

1988.
G. Carlsson. Topology and data. Bull. Amer. Math. Soc. (N.S.), 46

B. T. Fasy, F. Lecci, A. Rinaldo, L. Wasserman, S. Balakrishnan, and A. Singh. Statistical inference for persistent homol-(2):255–308, 2009.
CGAL. Cgal, Computational Geometry Algorithms Library, 2013.

ogy: Confidence sets for persistence diagrams. arXiv preprint
arXiv:1303.7117, 2013.

http://www.cgal.org.
F. Chazal, B. T. Fasy, F. Lecci, A. Rinaldo, A. Singh, and L. Wasserman. On the bootstrap for persistence diagrams and landscapes.
arXiv preprint arXiv:1311.0376, 2013.
C. Chen and M. Kerber. Persistent homology computation with atwist. In Proc. of the 27th European Workshop on Computational
Geometry, 2011.

A. Hatcher. Algebraic topology. Cambridge University Press,
Cambridge, 2002.
L. Jaillet and T. Siméon. Path deformation roadmaps: Compact
graphs with useful cycles for motion planning. Int. Journal of
Robotics Research, 27(11-12):1175–1188, 2008a.
L. Jaillet and T. Siméon. Path deformation roadmaps. In S. Akella,

K.-S. Chou and X.-P. Zhu. The curve shortening problem. CRC

N. M. Amato, W. H. Huang, and B. Mishra, editors, Algorithmic
Foundation of Robotics VII, volume 47 of Springer Tracts in

Press, 2001.
N. Cristianini and J. Shawe-Taylor.

An introduction to sup-port vector machines and other kernel-based learning methods.
Cambridge university press, 2000.
V. de Silva and R. Ghrist.

Coordinate-free coverage in sen-

Advanced Robotics, pages 19–34. Springer, 2008b.
L. Jaillet, J. Cortés, and T. Siméon. Sampling-Based Path Planning on Configuration-Space Costmaps. IEEE Transactions on
Robotics, 26:635–646, 2010. .

sor networks with controlled boundaries via homology. The

S. Karaman and E. Frazzoli. Sampling-based algorithms for optimal

International Journal of Robotics Research, 25(12):1205–1222,

motion planning. Int. Journal of Robotics Research, 30(7):846–
894, 2011.

2006.
V. de Silva and R. Ghrist. Coverage in sensor networks via persis-

P. Katsikouli, R. Sarkar, and J. Gao. Persistence based online signal

tent homology. Algebraic & Geometric Topology, 7(1):339–358,

and trajectory simplification for mobile devices. In International

2007.

Conference on Advances in Geographic Information Systems

M. P. Deisenroth, D. Fox, and C. E. Rasmussen. Gaussian pro-

(ACM SIGSPATIAL). 2014.

cesses for data-efficient learning in robotics and control. IEEE

L. E. Kavraki, P. Svestka, J.-C. Latombe, and M. H. Overmars.

Transactions on Pattern Analysis and Machine Intelligence, 99

Probabilistic roadmaps for path planning in high-dimensional(PrePrints):1, 2013. ISSN 0162-8828. .

configuration spaces. IEEE Trans. on Robotics and Automation,

R. Diankov and J. Kuffner. OpenRAVE: A Planning Architecture forAutonomous Robotics. Technical Report CMU-RI-TR-08-34,
Robotics Institute, Pittsburgh, PA, July 2008.
R. A. Dwyer. Higher-dimensional voronoi diagrams in linear
expected time.

Discrete & Computational Geometry, 6(1):343–367, 1991.

12(4):566–580, 1996.
S. Kim, K. Sreenath, S. Bhattacharya, and V. Kumar. Optimal trajectory generation under homology class constraints. In 51stIEEE Conf. on Decision and Control, 10-13 Dec 2012.
R. A. Knepper, S. S. Srinivasa, and M. T. Mason. Toward a deeper
understanding of motion alternatives via an equivalence relation

H. Edelsbrunner. Weighted alpha shapes. 1992.

on local paths. Int. Journal of Robotics Research, 31(2):167–186,

H. Edelsbrunner. The union of balls and its dual shape. Discrete

2012.

22

Journal name 000(00)Press, 2006.

J.-C. Latombe. Robot motion planning. Springer, 1991.
S. M. LaValle. Planning algorithms. Cambridge University Press,

J. T. Schwartz and M. Sharir. On the “piano movers” problem. II.
General techniques for computing topological properties of real

2006.
S. M. LaValle and J. J. Kuffner. Rapidly-Exploring Random Trees:
Progress and Prospects. In B. R. Donald, K. M. Lynch, andalgebraic manifolds. Advances in Applied Mathematics, 4(3):298–351, 1983.

D. Rus, editors, Algorithmic and Computational Robotics: New

E. Sisbot and R. Alami. A human-aware manipulation planner.

Directions, pages 293–308, Wellesley, MA, 2001. A K Peters.

Robotics, IEEE Transactions on, 28(5):1045–1057, Oct 2012.

S. R. Lindemann and S. M. LaValle. Current issues in samplingbased motion planning. In Robotics Research, pages 36–54.

L. Zhang, Y. J. Kim, and D. Manocha. A hybrid approach forcomplete motion planning. In Proc. of the IEEE/RSJ Int. Conf.

Springer, 2005.
B. Luders, M. Kothari, and J. How. Chance constrained rrt for probabilistic robustness to environmental uncertainty. In AIAA Guidance, Navigation, and Control Conference. American Instituteon Intelligent Robots and Systems, (IROS’07), pages 7–14. IEEE,2007.
Y. Zheng and X. Zhou.

Computing with spatial trajectories.

Springer, 2011.

of Aeronautics and Astronautics, 2014/11/07 2010. .
J. Mainprice, E. A. Sisbot, L. Jaillet, J. Cortés, R. Alami, andT. Siméon. Planning human-aware motions using a samplingbased costmap planner. In Robotics and Automation (ICRA),
2011 IEEE International Conference on, pages 5012–5017.
IEEE, 2011.
E. Masehian and D. Sedighizadeh. Classic and heuristic approaches
in robot motion planning - a chronological review.

ISSN 1552-3098. .

World

Academy of Science, Engineering and Technology, 23:101–106,
2007.
P. McMullen. The maximum numbers of faces of a convex polytope.
Mathematika, 17(02):179–184, 1970.
B. Morris and M. Trivedi. Learning trajectory patterns by clustering: Experimental studies and comparative evaluation. In IEEE
Int. Conf. on Comp. Vision and Pattern Recognition (CVPR’09),
pages 312–319. IEEE, 2009.
J. R. Munkres. Elements of algebraic topology. volume 2. AddisonWesley Reading, 1984.
P. Niyogi, S. Smale, and S. Weinberger. Finding the homologyof submanifolds with high confidence from random samples.
Discrete and Comp. Geometry, 39(1-3):419–441, 2008.
F. T. Pokorny, C. H. Ek, H. Kjellström, and D. Kragic. Topologicalconstraints and kernel-based density estimation. In Advancesin Neural Information Processing Systems 25, Workshop on
Algebraic Topology and Machine Learning. 2012.
F. T. Pokorny, M. Hawasly, and S. Ramamoorthy. Multiscale topological trajectory classification with persistent homology. In
Proceedings of Robotics: Science and Systems, Berkeley, USA,
July 2014.
C. E. Rasmussen and C. Williams. Gaussian processes for machinelearning. In Gaussian Processes for Machine Learning. MIT

A. Zomorodian. Fast construction of the vietoris-rips complex.
Computers & Graphics, 34(3):263–271, 2010.

23

Running head right side

Lemma. Let K1 ⊂ . . . ⊂ Kn be a simplexwise filtration

AppendixHere, we provide further details of the proof of Lemma 4.1.
While we could not find the exact statement of this lemmain the literature, this result is a direct consequence of theabstract decomposition theorems of persistent homology(Edelsbrunner and Harer, 2010). We sketch the details of
direct proofs which are obtained simply by reasoning about
reduced matrices and which do not require further theory toprovide a ‘hands on’ practical approach to these results. Let
us recall that we denote by R the reduced boundary matrix,
by D the original boundary matrix, and R = DV for thetransformation matrix V . When M is a matrix Mi denotes
the ith column and low(Mi ) is defined to be the index of
the lowest non-zero element of Mi when Mi 6= 0 and low
is undefined otherwise. We begin with two claims:
Claim 1. Let (i, j) ∈ P = {(i, j) : Rj 6= 0 and i =low(Rj )}. Then Ri = 0.
Proof. Suppose Ri 6= 0 instead. Then Ri 6∈ B(Ki−1 ) =span(R1 , . . . , Ri−1 ) since low is injective on the non-zerocolumns of R since R is reduced. We have ∂Rj = ∂(σi +w) = ∂σi + η = ∂∂Vj = 0, for some w ∈ C(Ki−1 ),η = ∂w ∈ B(Ki−1 ) and σi denotes the ith simplex in
the filtration. Hence ∂σi ∈ B(Ki−1 ), but Ri = ∂Vi =∂σi +w for some w ∈ B(Ki−1 ). But then Ri ∈ B(Ki−1 ).
00

A contradiction.

of simplicial complexes, let R = DV denote the reduced
boundary matrix after applying the left-to-right reduction
algorithm, and let Ep ⊆ E, Pp ⊆ P denote those elements
corresponding to p-cycles only. For 1 6 i 6 n, a basis ofZp (Ki ) is given by Spi = {Rt : (s, t) ∈ Pp , s 6 i} ∪ {Vs :s ∈ Ep , s 6 i}, and, for 1 6 i 6 j 6 n, the image of theset
Tpi,j ={Rt : (s, t) ∈ Pp , s 6 i, t > j}
∪ {Vs : s ∈ Ep , s 6 i}under the quotient map
q : Zp (Ki )  Hpi,j = Zp (Ki )/(Bp (Kj ) ∩ Zp (Ki ))forms a basis of Hpi,j . Finally #Ep = dim(Hp (Kn )).
Proof. Let us first show that Spi forms a basis for Zp (Ki ).
Note that each column of of R only contains non-zerocoefficients corresponding to simplices of a single fixed
dimension each since column additions are only performed
among columns of the same dimension. It is hence clearthat Spi ⊆ Cp (Ki ). Note furthermore that Spi ⊆ Zp (Ki )since, for each Rt corresponding to (s, t) ∈ Pp , s 6 i, we
have ∂Rt = ∂∂Vt = 0. Similarly, for Vs , for s ∈ Ep ands 6 i, we have ∂Vs = Rs = 0 by the definition of Ep .
Next, we show that the elements of Spi are linearly inde-

Claim 2. With the same notation as in the main lemma, so

pendent. Define Aip = {Rt : (s, t) ∈ Pp , s 6 i} andK = ∪ni=1 σi , we have

Bpi = {Vs : s ∈ Ep , s 6 i}. We observe that low is injectivea) #P = #non-zero columns of R, where (i, j) ∈ P
corresponds to column Rj .
b) #P + #E = #zero columns of R, where (i, j) ∈ P
corresponds to column Ri and s ∈ E to column Rs .
c) 2#P + #E = n

on Aip . To see this, assume that (i, j), (i, j 0 ) ∈ Pp for j 6= j 0 .
Without loss of generality, assume j < j 0 , then we obtain acontradiction since R would not be reduced as we could add
Rj to Rj0 to further reduce R. Similarly, low(Vi ) = i for all
i 6 n, as can be seen by observing that V is the result ofa sequence of left-to-right column additions performed on
the identity matrix, so low is also injective on Bpi . Finally

Proof. Part a): If Ri 6= 0, low(Ri ) = s for some s, solow(Aip )∩low(Bpi ) = ∅ by the definition of E and P . Hence(s, i) ∈ P and s is unique since R is reduced. Also, for each

low is injective on Spi = Aip ∪ Bpi and the set is hence lin-(s, i) ∈ P , Ri 6= 0. Part b): If Ri = 0, either there exists

early independent. To show that Spi is a basis for Zp (Ki ),a j such that i = low(Rj ) in which case (i, j) ∈ P and

we now just need to show dim(Zp (Ki ) = #Spi .

this j is unique since R is reduced, or there is no such j, in

Since the upper-left submatrices of the reduced matrix R

which case i ∈ E by definition of E. Part c) follows fromwhich correspond to the subcomplexes Ki for i 6 n are

the above since each column of R is either zero or non-zerothemselves reduced, and due to the grading by p, it is not

and R has n columns.

difficult to extend the reasoning in Claim 2 to show that

24

Journal name 000(00)dim Bp (Ki ) = #Pp6i and dim Cp (Ki ) = 2#Pp6i +#Ep6i ,for all i 6 n, where we define Pp6i = {(a, j) : Rj 6= 0, a =low(Rj ), dim(σa ) = p and j 6 i} and Ep6i = {t : Rt =0, low(Rj ) 6= t for all j ∈ {1, . . . , i} and dim(σt ) = p}.
By the rank-nullity theorem, it follows that dim(Zp (Ki )) =#Pp6i + #Ep6i for all i 6 n. Since it is clear that this isalso the cardinality of the set Spi , the result follows.
Next, we prove the statement about Tpi,j . Note that Tpi,jis a linearly independent set since low is injective on Tpi,j .
We need to show that Bp (Kj ) ∩ Zp (Ki ), i 6 j, is spannedby Spi − Tpi,j = {Rt : (s, t) ∈ Pp , s 6 i and t 6 j}. First,
observe that
Bp (Kj ) = span{Rt : (s, t) ∈ Pp , t 6 j}.
This follows since span(V1 , . . . , Vj ) = C(Kj ) and if ∂Vi 6=0, then Ri = ∂Vi is in the above spanning set which is alsolinearly independent.
Next, we know that Zp (Ki ) has basis Spi . From this,
we observe, Spi − Tpi,j ⊂ Bp (Kj ) ∩ Zp (Ki ). Finally,
Bp (Kj ) + Zp (Ki ) has basis {Rt : (s, t) ∈ Pp , t 6j} ∪ Spi = {Vs : s ∈ Ep , s 6 i} ∪ {Rt : (s, t) ∈Pp , t 6 j or s 6 i} which can be written as a disjoint unionof Spi and {Rt : (s, t) ∈ Pp , i < s 6 t 6 j}. Counting dimensions, we find that dim(Bp (Kj ) ∩ Zp (Ki )) =dim(Bp (Kj ))+dim(Zp (Ki ))−dim(Bp (Kj )+Zp (Ki )) =#{Rt : (s, t) ∈ Pp , t 6 j} − #{Rt : (s, t) ∈ Pp : i <s 6 t 6 j} = #{Rt : (s, t) ∈ Pp : s 6 i and t 6 j} =#(Spi − Tpi,j ). To see that #Ep = dim(Hp (Kn )), observethat Hp (Kn ) = Hpn,n (Kn ) and Tpn,n = Ep .