This simple code calculates the convex hull of a set of 2D points and generates EPS files to visualise them. The algorithm was taken from a textbook on Computional Geometry.
Python, 196 lines
Download
Copy to clipboard1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 | #!/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)
|
Finding convex hulls is a fundamental problem in Computational Geometry and a basic building block for solving many problems. The given implementation is not guaranteed to be numerically very stable. It might benefit from using the Numeric package for gaining more performance for very large sets of points.
Tags: algorithms
Small Bug: Only works with a list of UNIQUE points. If the list of points passed to this function is not unique, it will raise an assertion. To fix this, remove these lines from the beginning of the convexHull function:
Replace this with these lines:
string.atoi()is deprecated. you should useint()instead