#!/usr/bin/env python """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. Prints 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, printing the result as an EPS file. Usage: convexhull.py Dinu C. Gherman """ import sys, string, random ###################################################################### # Helpers ###################################################################### def _myDet(p, q, r): """Calc. determinant of a special matrix with three 2D points. The sign, "-" or "+", determines the side, right or left, respectivly, 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 if _myDet(p, q, r) < 0: return 1 else: return 0 def _isPointInPolygon(r, P): "Is point r inside a given polygon P?" # We assume 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." # Fill a square with 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 != 0: 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 /cross %% cross, x, y --> - { 0 360 arc %% draw cross hair } 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 halfs 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)