Quickhull (Python, arrays)

From LiteratePrograms

Other implementations: Javascript | Python, arrays

Compute the convex hull for a set of points in the plane using quickhull.

Contents 1 theory 2 practice 2.1 partitioning 2.2 pivoting 2.3 getting started 3 wrapping up

theory

Quickhull gets its name because of its strong resemblance to quicksort. Just as in quicksort, we pick pivots allowing us to partition the problem into smaller subsets, and just as in quicksort, the resulting call graph naturally forms a tree.

Instead of sorting a sublist, the basic recursive step of quickhull is building a dome on a base edge: given a base edge, the triangle formed by adding the furthest possible point outside of the edge produces two new candidate base edges; each of those may have a dome constructed on it, and the process continues until we have built a dome-like section of the convex hull in between the points on the original base edge.

Exercise: Create a picture demonstrating the process of pivot triangulation and dome construction for this article.

practice

Quickhull, in the abstract, can find the convex hull of data sets in any number of dimensions. For ease of exposition, we limit ourselves to planar (2D) data, and do not worry about dealing with non-generic geometries.

partitioning

To partition, we find the signed distances of each of the points in the sample from the given base edge, then discard all except those with positive distances — the outer points.

To calculate the distances, dists, we project each sample point (using the dot product) onto the vector normal to the edge.

<<dome>>=
def dome(sample,base): 
    h, t = base
    dists = dot(sample-h, dot(((0,-1),(1,0)),(t-h)))
    outer = repeat(sample, dists>0, 0)

Question: why subtract h from both sides of the dot product?

pivoting

To pivot, we find the point furthest outside of our base edge. We then form a cone on the base edge, with one new edge between the old head point and the pivot, and the other new edge between the pivot and old tail point, and return the result of linking together domes built on both of the new edges. If no pivot exists, we are already on the hull, and simply return the current edge.

<<dome>>=
    if len(outer):
        pivot = sample[argmax(dists)]
        return link(dome(outer, edge(h, pivot)),
                    dome(outer, edge(pivot, t)))
    else:
        return base

Question: In one dimension, adding a pivot is trivial: the new boundary is the pivot point itself. In two dimensions, adding a pivot is easy: the new boundary is built on the cone consisting of the two new edges of the triangle formed by the pivot and the old base edge. What extra work must be done to construct the facets of a cone from a pivot in three or more dimensions?

getting started

In order to get the recursion started, we must find an initial edge running between points that will eventually be on the convex hull. A simple way to do this is to pick an arbitrary direction (here, wlog, the x-axis) and to find the extremal points along that direction. Building a dome on this initial base edge will only find the hull above (outside) the edge, so we link it with a dome built on the same edge, running the opposite direction (with a stride of -1), thus adding the hull below (outside the reversed edge). By building domes on both sides we construct the complete hull.

<<qhull>>=
def qhull(sample):
    dome
    if len(sample) > 2:
    	axis = sample[:,0]
    	base = take(sample, [argmin(axis), argmax(axis)], 0)
    	return link(dome(sample, base),
                    dome(sample, base[::-1]))
    else:
	return sample

Question: Why is qhull of 2 points trivial but not qhull of 3?

Exercise: We could even start with the quadrilateral between the extremal points on both x and y axes. Rewrite qhull to link together the domes built starting from these four edges.

Exercise: What happens when the points all coincide on the x-axis? Fix the initialization of base to handle this degeneracy.

Question: In one dimension, starting with extremal points trivially produces the convex hull. In three or more dimensions, will this trick still work to pick an initial plane? Is it possible to start with a random initial flat?

wrapping up

Time to define the ancillary functions:

• link links two lists of edges. As the final point in the first list is always the same as the initial point of the second, we make sure to only include it once.
• edge constructs an edge between two points

<<ancillary functions>>=
link = lambda a,b: concatenate((a,b[1:]))
edge = lambda a,b: concatenate(([a],[b]))

Finally, if this module is run from the command line, it runs a quick smoke test by generating a Postscript visualization of the results of qhull applied to a random test vector.

<<quickhull2d.py>>=
from numpy import *
ancillary functions
qhull
if __name__ == "__main__":
    #sample = 10*array([(x,y) for x in arange(10) for y in arange(10)])
    sample = 100*random.random((32,2))
    hull = qhull(sample)
    print "%!"
    print "100 500 translate 2 2 scale 0 0 moveto"
    print "/tick {moveto 0 2 rlineto 0 -4 rlineto 0 2 rlineto"
    print "              2 0 rlineto -4 0 rlineto 2 0 rlineto} def"
    for (x,y) in sample:
	print x, y, "tick"
    print "stroke"
    print hull[0,0], hull[0,1], "moveto"
    for (x,y) in hull[1:]:
	print x, y, "lineto"
    print "closepath stroke showpage"

sample output:

In low dimensions, the convex hull can be a great simplification — for this sample, there are possible pairs of sample points, yet only 6 of these pairs suffice to describe the hull.

Question: how does the ratio of surface area to volume behave in higher dimensions?

Download code