Welcome, guest | Sign In | My Account | Store | Cart
```#!/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!

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

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

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"