Simplicial Complexes from Point Clouds

Vietoris-Rips Complex

Given a set of points S in (ℝ^n) and a parameter ϵ>0, the Vietoris-Rips complex (Subscript[VR, ϵ])(S) is the complex whose simplices consist of subsets of S consisting of points that are pairwise within ϵ. That is, {σ⊂S : d(a,b)<=ϵ ∀ a,b ∈σ}.

Here is a function that computes the Vietoris-Rips complex.
Input: A list of points S in (ℝ^n) and a parameter ϵ>0
Output: A list of all simplices in (Subscript[VR, ϵ])(S), specified as lists of vertex indexes

computeVRComplex[points_,ϵ_]:=Module[
{n=Length[points],complex,include},

(* start with all vertices *)
complex=Table[{k},{k,n}];

(* now consider higher-dimensional simplices *) 
Do[ (* loop over all sets of at least 2 vertices *) 
include=True;
Do[ (* check all pairwise distances among this set of vertices *) 
If[EuclideanDistance[points[[pair[[1]]]],points[[pair[[2]]]]]>ϵ,include=False;Break[]]
,{pair,Subsets[sub,{2}]}];

(* if all pairwise distances are less than ϵ, then include this simplex in the complex *) 
If[include,AppendTo[complex,sub]]

,{sub,Subsets[Range[n],{2,n}]}];

(* return the complex *) 
complex
]

Try it out. First generate some points:

numPoints=10;
points=RandomReal[1,{numPoints,2}];
pointPlot=ListPlot[points->Range[numPoints],PlotRange->{{0,1},{0,1}},AspectRatio->1]

complex=computeVRComplex[points,0.3]

Visualizing a 2D simplicial complex

Here is a version of our geometric realization algorithm from last week. This version is modified to label the points.

drawGeometricRealization[complex_,points_]:=Module[
{verts={}, edges={}, faces={}},

(* make lists of simplices of each dimension *) 
Do[
If[Length[simplex]==1,AppendTo[verts,points[[simplex[[1]]]]]];
If[Length[simplex]==2,AppendTo[edges,points[[simplex]]]];
If[Length[simplex]==3,AppendTo[faces,points[[simplex]]]]
,{simplex,complex}];

(* return the drawing *) 
Show[
ListPlot[points->Range[Length[points]],AspectRatio->1,PlotStyle->{PointSize[Medium],Black}],
Graphics[{Blue,Opacity[0.4],Triangle[faces],Black,Opacity[1],Line[edges]}] 
]
]

drawGeometricRealization[complex,points]

Čech Complex

Given a set of points S in (ℝ^n) and a parameter ϵ>0, the Čech complex (Subscript[C, ϵ])(S) is the complex whose simplices consist of subsets of S consisting of points that have mutually-intersecting ϵ-radius neighborhoods. That is, {σ⊂S : ∃y∈(ℝ^n) : d(a,y)<=ϵ ∀ a ∈σ}.

Smallest Enclosing Circle

To compute the Čech complex, we need to compute the smallest enclosing sphere for a set of points. The following code cell does this using a recursive algorithm known as Welzl’s algorithm.

(* function to return the center and radius of the circle passing through three points in the plane *) 
circ3pts[pts_]:=Module[{circ},
circ=Circumsphere[pts];
{circ[[1]],circ[[2]]}
]

(* function to return the center and radius of the circle defined by a diameter *) 
circ2pts[pts_]:=Module[{ctr,rad},
ctr=Midpoint[pts];
rad=EuclideanDistance[pts[[1]],pts[[2]]]/2;
{ctr,rad}
]

(* Welzl's algorithm returns the center and radius of the smallest circle enclosing a set of points in the plane *) 
smallestCircle[points_,boundary_:{}]:=Module[{p, center, radius},
(*Print[{points,boundary}];*)
(* if we have 3 points on the boundary, then return the circle they define *)
If[Length[boundary]==3,Return[circ3pts[boundary]]];
(* other end cases *) 
If[Length[boundary]==2&&Length[points]==0,Return[circ2pts[boundary]]];
If[Length[boundary]==1&&Length[points]==1,Return[circ2pts[Join[boundary,points]]]];
If[Length[boundary]==0&&Length[points]==1,Return[{points[[1]],0}]];

(* remove the last point and recurse *)
p=Last[points];
{center,radius}=smallestCircle[points[[;;-2]],boundary];

(* if p is in the circle, then return the circle *)
If[EuclideanDistance[p,center]<=radius,Return[{center,radius}]];

(* otherwise, p is not in the circle, so it must be on the boundary *)
Return[smallestCircle[points[[;;-2]],Append[boundary,p]]]
]

Try out Welzl’s algorithm:

numPoints=10;
points=RandomReal[1,{numPoints,2}];
circ=smallestCircle[points];
Show[ListPlot[points->Range[numPoints],PlotRange->{{0,1},{0,1}},AspectRatio->1],
Graphics[{Blue,Circle[circ[[1]],circ[[2]]]}]]

Now we can compute the Čech Complex

The following function computes the Čech complex up to 2-simplices (triangles).

computeCech[points_,ϵ_]:=Module[
{n=Length[points],complex},
(* vertices *)
complex=Table[{k},{k,n}];

(* edges *)
Do[
If[EuclideanDistance[points[[simp[[1]]]],points[[simp[[2]]]]]<2ϵ,AppendTo[complex,simp]]
,{simp,Subsets[Range[n],{2}]}];

(* triangles *)
Do[
If[        ]
,{simp,Subsets[Range[n],{3}]}];

(* return the complex *) 
complex
]

numPoints=16;
points=RandomReal[1,{numPoints,2}];
complex=computeCech[points,0.13]
drawGeometricRealization[complex,points]

Manipulating Points

Mathematica’s Canvas object gives you an interactive drawing environment:

Canvas[Graphics[Point[points]]]

You make your own point cloud by starting with one point and making copies of it:

input=Canvas[Graphics[Point[{0,0}]]]

input=(Canvas[(CanvasFrame)]);

Now extract the list of Point objects from the canvas:

pointList=input[[1,1]]

Extract the coordinates from the Point objects:

points=Part[#,1] &/@pointList

Discussion Question

What are the relationships between the Vietoris-Rips complex (Subscript[VR, ϵ])(S) and the Čech complex (Subscript[C, ϵ])(S)? As a function of the threshold parameter, when is one complex a subcomplex of the other?

A point cloud with multiscale features

input=(Canvas[(CanvasFrame)]);
pointsM=Part[#,1] &/@input[[1,1]]

complexM=computeVRComplex[pointsM,0.3];
drawGeometricRealization[complexM,pointsM]

Euler Characteristic

From the homework: The Euler characteristic of a simplicial complex K is integer χ(K)=(∑_(i=0)^dimK)(((-1))^i)(Subscript[K, i]), where (Subscript[K, i]) is the number of i-dimensional simplices in K. For example, if K={{1},{2},{3},{4},{1,2},{2,3},{2,4},{3,4},{2,3,4}}, then χ(K)=4-4+1=1 since K contains four 0-dimensional simplices, four 1-dimensional simplices, and one 2-dimensional simplex.

Exploration: Compute the Euler characteristic of some simplicial complexes of your choice. What does the Euler characteristic tell us about the simplicial complex?