Alpha shape example¶

The example given in examples/alphashapes/alphashapes.py takes a filename containing points in 2D or 3D on the command line. It generates the alpha shape filtration of those points, and computes its persistence. It then outputs the persistence diagram in the format of a point (dimension, birth, death) per line.

# Computes the persistence diagram of the alpha shapes in both 2D and 3D 
# (decided dynamically based on the input file)

from    dionysus        import Filtration, StaticPersistence, data_dim_cmp, vertex_cmp, \
                               fill_alpha3D_complex, fill_alpha2D_complex, points_file
from    sys             import argv, exit
from    math            import sqrt


if len(argv) < 2:
    print "Usage: %s POINTS" % argv[0]
    exit()

points = [p for p in points_file(argv[1])]
f = Filtration()
if   len(points[0]) == 2:           # 2D
    fill_alpha2D_complex(points, f)
elif len(points[1]) == 3:           # 3D
    fill_alpha3D_complex(points, f)

print "Total number of simplices:", len(f)

f.sort(data_dim_cmp)
print "Filtration initialized"

p = StaticPersistence(f)
print "StaticPersistence initialized" 

p.pair_simplices()
print "Simplices paired"

print "Outputting persistence diagram"
smap = p.make_simplex_map(f)
for i in p:
    if i.sign():
        b = smap[i]
        if i.unpaired():
            print b.dimension(), sqrt(b.data[0]), "inf"
            continue

        d = smap[i.pair()]
        print b.dimension(), sqrt(b.data[0]), sqrt(d.data[0])

After the points are read into the list points, the functions fill_alpha*_complex fill the Filtration with the simplices of the Delaunay triangulation. Each one has its data attribute set to the tuple consisting of its alpha shape value (the minimum value of the squared distance function on its dual Voronoi cell) and whether the simplex is regular or critical.

The filtration then sorts the simplices with respect to their data and dimension (via data_dim_cmp()):

f = Filtration()
fill_alpha*_complex(points, f)
f.sort(data_dim_cmp)

We initialize StaticPersistence, and pair the simplices:

p = StaticPersistence(f)
p.pair_simplices()

Iterating over the StaticPersistence, we output the points of the persistence diagram (dimension, birth, death) in the last for loop. If the simplex is unpaired (i.unpaired()), the class it creates survives till infinity.