#!/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 <numPoints> <squareLength> <outFile>
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)
