How can I find the alpha shape (concave hull) of a 2d point cloud?

I am looking for an implementation that calculates alpha shapes in two dimensions. I am running ubuntu. I would prefer a command line utility for this task, but will also be fine with a python library.

In Google I have found many implementations that calculate alpha shapes. But none of them output what I want. As input I have a list of two dimensional points (e.g., one pair of floats per line in a text file). As output I want another list of two dimensional points with the same scale.

I have tried installing the latest python bindings of cgal, but these have not been supported in a while and no longer compile on Ubuntu 11.04 (I also tried on Ubuntu 10.04 and had no luck). Clustr, a project developed at flickr by Aaron Straup Cope will also not compile on Ubuntu 11.04 (probably because it is also tied to older CGAL libraries).

I also tried this implementation from a Ken Clarkson at bell labs. It outputs nearly what I want, the output seems to be in another scale and it turns floats into ints.

I also tried the python bindings of dionysus. These compiled, but when I fed the function fill_alpha2D_complex(points, f) with my list of points, the output was not what I expected. It was not a list of two dimensional points, but rather seems to be a “persistence diagram” I don’t know what that means.

Anyone know of a simple solution to this problem?

UPDATE: I want to print out the points associated with the alpha shape where it is on the verge of not being connected anymore. I think this means “give me the points associated with the smallest alpha value such that the shape is connected.”

UPDATE I now found out how to get what I want from Ken Clarkson’s implementation and (more or less what I want) from the dionysus implementation. Clarkson’s implementation was doing the right thing, it just output the indices of the points rather than the points themselves (same story with dionysus), and I needed to get some optional flags right. The wrapper I wrote is below. This solution is ideal because it produces an alpha shape that is both connected and does not contain holes. Alpha is set automatically. Dionysus, on the other hand, does not automatically discover this values of alpha. Plus Clarkson’s implementation can be set to output a ps image of the shape (with the -afps flag). To get Clarkson’s code to compile with non-ancient version of GCC, you need to follow the step outlined here. The following code can be used as a library or as a stand alone wrapper:

#!/usr/bin/python -O

import sys, os
import subprocess
import tempfile

hull_path = "./hull.exe"

def get_alpha_shape(points):
    # Write points to tempfile
    tmpfile = tempfile.NamedTemporaryFile(delete=False)
    for point in points:
        tmpfile.write("%0.7f %0.7fn" % point)
    tmpfile.close()

    # Run hull
    command = "%s -A -m1000000 -oN < %s" % (hull_path, tmpfile.name)
    print >> sys.stderr, "Running command: %s" % command
    retcode = subprocess.call(command, shell=True)
    if retcode != 0:
        print >> sys.stderr, "Warning: bad retcode returned by hull.  Retcode value:" % retcode
    os.remove(tmpfile.name)

    # Parse results
    results_file = open("hout-alf")
    results_file.next() # skip header
    results_indices = [[int(i) for i in line.rstrip().split()] for line in results_file]
#    print "results length = %d" % len(results_indices)
    results_file.close()
    os.remove(results_file.name)

    return [(points[i], points[j]) for i,j in results_indices]

if __name__ == "__main__":
    points = [tuple([float(i) for i in line.rstrip().split()]) for line in sys.stdin]
    for point_i, point_j in get_alpha_shape(points):
        sys.stdout.write("%0.7f,%0.7ft%0.7f,%0.7fn" % (point_i[0], point_i[1], point_j[0], point_j[1]))
    sys.exit(0)

Answer

I found this in the dionysus docs which might give you the alpha shape:

complex = Filtration()
fill_alpha2D_complex(points, complex)
alphashape = [s for s in complex if s.data[0] <= .5]

Then I believe you need to do something like:

for simplex in alphashape:
    print [v for v in simplex.vertices]

Advertisement

Answer