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 " % format % H[0]) for p in H: f.write("%s lineto " % format % p) f.write("%s lineto " % format % H[0]) f.write("stroke ") # Save the whole list of points as individual dots for p in P: f.write("%s r circle " % format % p) f.write("stroke ") # Save footer f.write(" showpage ") # 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)