I l@ve RuBoard |

## 17.19 Module: Finding the Convex Hull of a Set of 2D PointsCredit: 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
## 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 AlsoComputational Geometry: Algorithms and Applications, 2nd edition, by M. de Berg, M. van Kreveld, M. Overmars, and O. Schwarzkopf (Springer-Verlag). |

I l@ve RuBoard |