17.19 Module: Finding the Convex Hull of a Set of 2D Points

Credit: Dinu C. Gherman

Convex hulls of point sets are an important building block in many computational-geometry applications. Example 17-1 calculates the convex hull of a set of 2D points and generates an Encapsulated PostScript (EPS) file to visualize it. Finding convex hulls is a fundamental problem in computational geometry and is a basic building block for solving many problems. The algorithm used here can be found in any good textbook on computational geometry, such as Computational Geometry: Algorithms and Applications, 2nd edition (Springer-Verlag). Note that the given implementation is not guaranteed to be numerically stable. It might benefit from using the Numeric package for gaining more performance for very large sets of points.

Example 17-1. Finding the convex hull of a set of 2D points

""" convexhull.py

Calculate the convex hull of a set of n 2D points in O(n log n) time.
Taken from Berg et al., Computational Geometry, Springer-Verlag, 1997.
Emits output as EPS file.

When run from the command line, it generates a random set of points
inside a square of given length and finds the convex hull for those,
emitting the result as an EPS file.
Usage:
    convexhull.py <numPoints> <squareLength> <outFile>

Dinu C. Gherman
"""
import sys, string, random

# helpers

def _myDet(p, q, r):
    """ Calculate determinant of a special matrix with three 2D points.

    The sign, - or +, determines the side (right or left, respectively) on which
    the point r lies when measured against a directed vector from p to q.
    """
    # We use Sarrus' Rule to calculate the determinant
    # (could also use the Numeric package...)
    sum1 = q[0]*r[1] + p[0]*q[1] + r[0]*p[1]
    sum2 = q[0]*p[1] + r[0]*q[1] + p[0]*r[1]
    return sum1 - sum2

def _isRightTurn((p, q, r)):
    "Do the vectors pq:qr form a right turn, or not?"
    assert p != q and q != r and p != r
    return _myDet(p, q, r) < 0

def _isPointInPolygon(r, P):
    "Is point r inside a given polygon P?"
    # We assume that the polygon is a list of points, listed clockwise
    for i in xrange(len(P)-1):
        p, q = P[i], P[i+1]
        if not _isRightTurn((p, q, r)):
            return 0 # Out!
    return 1 # It's within!

def _makeRandomData(numPoints=10, sqrLength=100, addCornerPoints=0):
    "Generate a list of random points within a square (for test/demo only)"
    # Fill a square with N random points
    min, max = 0, sqrLength
    P = []
    for i in xrange(numPoints):
        rand = random.randint
        x = rand(min+1, max-1)
        y = rand(min+1, max-1)
        P.append((x, y))
    # Add some "outmost" corner points
    if addCornerPoints:
        P = P + [(min, min), (max, max), (min, max), (max, min)]
    return P

# output

epsHeader = """%%!PS-Adobe-2.0 EPSF-2.0
%%%%BoundingBox: %d %d %d %d

/r 2 def                %% radius

/circle                 %% circle, x, y, r --> -
{
    0 360 arc           %% draw circle
} def

1 setlinewidth          %% thin line
newpath                 %% open page
0 setgray               %% black color
"""

def saveAsEps(P, H, boxSize, path):
    "Save some points and their convex hull into an EPS file."

    # Save header
    f = open(path, 'w')
    f.write(epsHeader % (0, 0, boxSize, boxSize))

    format = "%3d %3d"

    # Save the convex hull as a connected path
    if H:
        f.write("%s moveto\n" % format % H[0])
        for p in H:
            f.write("%s lineto\n" % format % p)
        f.write("%s lineto\n" % format % H[0])
        f.write("stroke\n\n")

    # Save the whole list of points as individual dots
    for p in P:
        f.write("%s r circle\n" % format % p)
        f.write("stroke\n")

    # Save footer
    f.write("\nshowpage\n")

# public interface

def convexHull(P):
    "Calculate the convex hull of a set of points."

    # Get a local list copy of the points and sort them lexically
    points = map(None, P)
    points.sort(  )

    # Build upper half of the hull
    upper = [points[0], points[1]]
    for p in points[2:]:
        upper.append(p)
        while len(upper) > 2 and not _isRightTurn(upper[-3:]):
            del upper[-2]

    # Build lower half of the hull
    points.reverse(  )
    lower = [points[0], points[1]]
    for p in points[2:]:
        lower.append(p)
        while len(lower) > 2 and not _isRightTurn(lower[-3:]):
            del lower[-2]

    # Remove duplicates
    del lower[0]
    del lower[-1]

    # Concatenate both halves and return
    return tuple(upper + lower)

# Test

def test(  ):
    a = 200
    p = _makeRandomData(30, a, 0)
    c = convexHull(p)
    saveAsEps(p, c, a, file)

if _ _name_ _ == '_ _main_ _':
    try:
        numPoints = string.atoi(sys.argv[1])
        squareLength = string.atoi(sys.argv[2])
        path = sys.argv[3]
    except IndexError:
        numPoints = 30
        squareLength = 200
        path = "sample.eps"

    p = _makeRandomData(numPoints, squareLength, addCornerPoints=0)
    c = convexHull(p)
    saveAsEps(p, c, squareLength, path)

17.19.1 See Also

Computational Geometry: Algorithms and Applications, 2nd edition, by M. de Berg, M. van Kreveld, M. Overmars, and O. Schwarzkopf (Springer-Verlag).