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

Question

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.7f\n" % 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.7f\t%0.7f,%0.7f\n" % (point_i[0], point_i[1], point_j[0], point_j[1]))
    sys.exit(0)

Australian

A mere list of points cannot describe an alpha shape, as it is not even necessarily connected. – Alteration

I just want a nice boundary of a point cloud. But if this description is not exact enough see my update in the question. – Australian

a connected alpha shape may still be not very nice, i.e. contain holes. – Alteration

Ok, well I would prefer a "nice" shape, but I'm not versed in the terminology of computational geometry so I don't know how to describe what a "nice shape" is. But I feel it has intuitive properties, such as being connected, not containing holes, etc. I am also willing to just try different values of alpha until I get one that is "nice." – Australian

Answer 1

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]

Quarterphase

That prints out a long list of tuples. Most of the tuples have one member, some have two or three. There are around five thousand of these tuples, even though my original input consisted of around one thousand points. Instead, I was expecting a list of two dimensional points that had fewer members than the input. – Australian

After some explanation from the author of dionysus, I've concluded that this answer is right. The integers in each of the tuples are the indices of the points. The tuples with two members are edges, those with three are triangles. The ".5" is the value of alpha. There is a way to determine the smallest value of alpha such that the resulting shape is connected. If I have time, I'll place this code in the question. – Australian

ok, I posted it. in the end I went for "hull" over dionysus for reasons I outline in the solution (which I've added to the question) above. – Australian

Excellent question & research. Thanks for that. I got as far as getting the list of vertex indices from the dionysus route. Now that I have them, how do I find the alpha such that I can form a concave polygon from the points? You said there is a way to figure out the alpha -- could you elaborate? – Scourge

I tried it but I got an error: NameError: name 'fill_alpha2D_complex' is not defined – Gingivitis

Answer 2

import alphashape
import numpy as np
from descartes import PolygonPatch
import matplotlib.pyplot as plt

n = 1000
points = np.random.rand(n, 2)

# Define alpha parameter
alpha= 20

# Generate the alpha shape
alpha_shape = alphashape.alphashape(points, alpha)

# Initialize plot
fig, ax = plt.subplots()

# Plot input points
ax.scatter(*zip(*points))

# Plot alpha shape
ax.add_patch(PolygonPatch(alpha_shape, alpha=.5, fc='blue', ec='red'))


pnts = [x for x in alpha_shape.boundary.coords]
plt.scatter(*zip(*pnts),marker='s')



plt.show()

Evidently