I l@ve RuBoard

### 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!

"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
P = P + [(min, min), (max, max), (min, max), (max, min)]
return P

# output

%%%%BoundingBox: %d %d %d %d

/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."

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"