1
2
3 """convexhull.py
4
5 Calculate the convex hull of a set of n 2D-points in O(n log n) time.
6 Taken from Berg et al., Computational Geometry, Springer-Verlag, 1997.
7 Prints output as EPS file.
8
9 When run from the command line it generates a random set of points
10 inside a square of given length and finds the convex hull for those,
11 printing the result as an EPS file.
12
13 Usage: convexhull.py <numPoints> <squareLength> <outFile>
14
15 Dinu C. Gherman
16
17 Small Bug: Only works with a list of UNIQUE points, Evan Jones, 2005/05/18
18 If the list of points passed to this function is not unique, it will raise an assertion.
19 To fix this, remove these lines from the beginning of the convexHull function:
20
21 Taken from http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/66527
22 and modified to work with complex numbers.
23
24 """
25
26
27 import sys, random
28
32
33
34
35
36
37
39 """Calc. determinant of a special matrix with three 2D points.
40
41 The sign, "-" or "+", determines the side, right or left,
42 respectivly, on which the point r lies, when measured against
43 a directed vector from p to q.
44 """
45
46
47
48
49
50
51
52 sum = q.real*(r.imag-p.imag) + p.real*(q.imag-r.imag) + r.real*(p.imag-q.imag)
53
54
55 return sum
56
57
59 "Do the vectors pq:qr form a right turn, or not?"
60 if p==q or p==r or q==r:
61 if p==q:
62 dup = p
63 elif p==r:
64 dup = p
65 elif q==r:
66 dup = q
67 raise DuplicatePoints, 'Two points are identical at %s' % str(dup)
68 return _myDet(p, q, r) < 0
69
70
72 "Is point r inside a given polygon P?"
73
74
75 for i in xrange(len(P[:-1])):
76 p, q = P[i], P[i+1]
77 if not _isRightTurn((p, q, r)):
78 return False
79
80 return True
81
82
84 "Generate a list of random points within a square."
85
86
87 mn, mx = 0, sqrLength
88 P = []
89 for i in xrange(numPoints):
90 rand = random.randint
91 y = rand(mn+1, mx-1)
92 x = rand(mn+1, mx-1)
93 P.append(complex(x, y))
94
95
96 if addCornerPoints != 0:
97 P = P + [complex(min, min), complex(max, max), complex(min, max), complex(max, min)]
98
99 return P
100
101
102
103
104
105
106 epsHeader = """%%!PS-Adobe-2.0 EPSF-2.0
107 %%%%BoundingBox: %d %d %d %d
108
109 /r 2 def %% radius
110
111 /circle %% circle, x, y, r --> -
112 {
113 0 360 arc %% draw circle
114 } def
115
116 /cross %% cross, x, y --> -
117 {
118 0 360 arc %% draw cross hair
119 } def
120
121 1 setlinewidth %% thin line
122 newpath %% open page
123 0 setgray %% black color
124
125 """
126
128 "Save some points and their convex hull into an EPS file."
129
130
131 f = open(path, 'w')
132 f.write(epsHeader % (0, 0, boxSize, boxSize))
133
134 format = "%3d %3d"
135
136
137 if H:
138 f.write("%s moveto\n" % format % (H[0].real, H[0].imag))
139 for p in H:
140 f.write("%s lineto\n" % format % (p.real, p.imag))
141 f.write("%s lineto\n" % format % (H[0].real, H[0].imag))
142 f.write("stroke\n\n")
143
144
145 for p in P:
146 f.write("%s r circle\n" % format % (p.real, p.imag))
147 f.write("stroke\n")
148
149
150 f.write("\nshowpage\n")
151
152
153
154
155
156
158 """Calculate the convex hull of a set of complex points.
159 If the hull has a duplicate point, an exception will be raised.
160 It is up to the application not to provide duplicates.
161 """
162
163
164 points = list(set(P))
165 points.sort(lambda a, b: 2*cmp(a.real, b.real)+cmp(a.imag, b.imag))
166
167
168 upper = [points[0], points[1]]
169 for p in points[2:]:
170 upper.append(p)
171 while len(upper)>2 and not _isRightTurn(upper[-3:]):
172 del upper[-2]
173 while len(upper)>2 and not _isRightTurn(upper[-2:]+upper[:1]):
174 del upper[-1]
175
176 points.reverse()
177
178 lower = [points[0], points[1]]
179 for p in points[2:]:
180 lower.append(p)
181 while len(lower)>2 and not _isRightTurn(lower[-3:]):
182 del lower[-2]
183 while len(lower)>2 and not _isRightTurn(lower[-2:]+lower[:1]):
184 del lower[-1]
185 us = set(upper)
186 for q in lower:
187 if q not in us:
188 upper.append(q)
189
190 return tuple(upper)
191
192
193
194
195
196
202
203
204
205
206
207 if __name__ == '__main__':
208 try:
209 numPoints = int(sys.argv[1])
210 squareLength = int(sys.argv[2])
211 path = sys.argv[3]
212 except IndexError:
213 numPoints = 30
214 squareLength = 200
215 path = "sample.eps"
216
217 p = _makeRandomData(numPoints, squareLength, addCornerPoints=0)
218 c = convexHull(p)
219 saveAsEps(p, c, squareLength, path)
220