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).