voronoi.py
1 |
#############################################################################
|
---|---|
2 |
#
|
3 |
# Voronoi diagram calculator/ Delaunay triangulator
|
4 |
# Translated to Python by Bill Simons
|
5 |
# September, 2005
|
6 |
#
|
7 |
# Calculate Delaunay triangulation or the Voronoi polygons for a set of
|
8 |
# 2D input points.
|
9 |
#
|
10 |
# Derived from code bearing the following notice:
|
11 |
#
|
12 |
# The author of this software is Steven Fortune. Copyright (c) 1994 by AT&T
|
13 |
# Bell Laboratories.
|
14 |
# Permission to use, copy, modify, and distribute this software for any
|
15 |
# purpose without fee is hereby granted, provided that this entire notice
|
16 |
# is included in all copies of any software which is or includes a copy
|
17 |
# or modification of this software and in all copies of the supporting
|
18 |
# documentation for such software.
|
19 |
# THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
|
20 |
# WARRANTY. IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
|
21 |
# REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
|
22 |
# OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
|
23 |
#
|
24 |
# Comments were incorporated from Shane O'Sullivan's translation of the
|
25 |
# original code into C++ (http://mapviewer.skynet.ie/voronoi.html)
|
26 |
#
|
27 |
# Steve Fortune's homepage: http://netlib.bell-labs.com/cm/cs/who/sjf/index.html
|
28 |
#
|
29 |
#############################################################################
|
30 |
|
31 |
def usage(): |
32 |
print """ |
33 |
voronoi - compute Voronoi diagram or Delaunay triangulation
|
34 |
|
35 |
voronoi [-t -p -d] [filename]
|
36 |
|
37 |
Voronoi reads from filename (or standard input if no filename given) for a set
|
38 |
of points in the plane and writes either the Voronoi diagram or the Delaunay
|
39 |
triangulation to the standard output. Each input line should consist of two
|
40 |
real numbers, separated by white space.
|
41 |
|
42 |
If option -t is present, the Delaunay triangulation is produced.
|
43 |
Each output line is a triple i j k, which are the indices of the three points
|
44 |
in a Delaunay triangle. Points are numbered starting at 0.
|
45 |
|
46 |
If option -t is not present, the Voronoi diagram is produced.
|
47 |
There are four output record types.
|
48 |
|
49 |
s a b indicates that an input point at coordinates a b was seen.
|
50 |
l a b c indicates a line with equation ax + by = c.
|
51 |
v a b indicates a vertex at a b.
|
52 |
e l v1 v2 indicates a Voronoi segment which is a subsegment of line number l
|
53 |
with endpoints numbered v1 and v2. If v1 or v2 is -1, the line
|
54 |
extends to infinity.
|
55 |
|
56 |
Other options include:
|
57 |
|
58 |
d Print debugging info
|
59 |
|
60 |
p Produce output suitable for input to plot (1), rather than the forms
|
61 |
described above.
|
62 |
|
63 |
On unsorted data uniformly distributed in the unit square, voronoi uses about
|
64 |
20n+140 bytes of storage.
|
65 |
|
66 |
AUTHOR
|
67 |
Steve J. Fortune (1987) A Sweepline Algorithm for Voronoi Diagrams,
|
68 |
Algorithmica 2, 153-174.
|
69 |
"""
|
70 |
|
71 |
#############################################################################
|
72 |
#
|
73 |
# For programmatic use two functions are available:
|
74 |
#
|
75 |
# computeVoronoiDiagram(points)
|
76 |
#
|
77 |
# Takes a list of point objects (which must have x and y fields).
|
78 |
# Returns a 3-tuple of:
|
79 |
#
|
80 |
# (1) a list of 2-tuples, which are the x,y coordinates of the
|
81 |
# Voronoi diagram vertices
|
82 |
# (2) a list of 3-tuples (a,b,c) which are the equations of the
|
83 |
# lines in the Voronoi diagram: a*x + b*y = c
|
84 |
# (3) a list of 3-tuples, (l, v1, v2) representing edges of the
|
85 |
# Voronoi diagram. l is the index of the line, v1 and v2 are
|
86 |
# the indices of the vetices at the end of the edge. If
|
87 |
# v1 or v2 is -1, the line extends to infinity.
|
88 |
#
|
89 |
# computeDelaunayTriangulation(points):
|
90 |
#
|
91 |
# Takes a list of point objects (which must have x and y fields).
|
92 |
# Returns a list of 3-tuples: the indices of the points that form a
|
93 |
# Delaunay triangle.
|
94 |
#
|
95 |
#############################################################################
|
96 |
import math |
97 |
import sys |
98 |
import getopt |
99 |
TOLERANCE = 1e-9
|
100 |
BIG_FLOAT = 1e38
|
101 |
|
102 |
#------------------------------------------------------------------
|
103 |
class Context( object ): |
104 |
def __init__(self): |
105 |
self.doPrint = 0 |
106 |
self.debug = 0 |
107 |
self.plot = 0 |
108 |
self.triangulate = False |
109 |
self.vertices = [] # list of vertex 2-tuples: (x,y) |
110 |
self.lines = [] # equation of line 3-tuple (a b c), for the equation of the line a*x+b*y = c |
111 |
self.edges = [] # edge 3-tuple: (line index, vertex 1 index, vertex 2 index) if either vertex index is -1, the edge extends to infiinity |
112 |
self.triangles = [] # 3-tuple of vertex indices |
113 |
self.extra_edges = [] # list of additional vertex 2-tubles (x,y) based on bounded voronoi tesselation |
114 |
self.set_bounds(None) |
115 |
self.use_bound = False |
116 |
self.xmin = self.ymin = self.xmax = self.ymax = None |
117 |
self.poly_vertices = [] # AJS |
118 |
|
119 |
def circle(self,x,y,rad): |
120 |
pass
|
121 |
|
122 |
def clip_line(self,edge,lid,rid): |
123 |
# here is where I will create false verticies if
|
124 |
# the voronoi line extends to infinity...
|
125 |
# the extra verticies will be added to the
|
126 |
# extra edges list as 2-tuples
|
127 |
a,b,c = edge.a,edge.b,edge.c |
128 |
if lid == -1: |
129 |
x = self.xmin
|
130 |
y = (c-a*x) / b |
131 |
if y < self.ymin or y > self.ymax: |
132 |
if y < self.ymin: y = self.ymin |
133 |
elif y > self.ymax: y = self.ymax |
134 |
x = (c-b*y) / a |
135 |
self.extra_edges.append((x,y))
|
136 |
ang0 = math.atan2((y - edge.reg[0].y),(x - edge.reg[0].x)) # AJS Add false vertices to list of polygon vertices |
137 |
ang1 = math.atan2((y - edge.reg[1].y),(x - edge.reg[1].x)) # AJS |
138 |
self.poly_vertices.append((edge.reg[0].sitenum,ang0,x,y)) # AJS |
139 |
self.poly_vertices.append((edge.reg[1].sitenum,ang1,x,y)) # AJS |
140 |
lid = -(len(self.extra_edges)-1) |
141 |
if rid == -1: |
142 |
x = self.xmax
|
143 |
y = (c-a*x) / b |
144 |
if y < self.ymin or y > self.ymax: |
145 |
if y < self.ymin: y = self.ymin |
146 |
elif y > self.ymax: y = self.ymax |
147 |
x = (c-b*y) / a |
148 |
self.extra_edges.append((x,y))
|
149 |
ang0 = math.atan2((y - edge.reg[0].y),(x - edge.reg[0].x)) # AJS As above |
150 |
ang1 = math.atan2((y - edge.reg[1].y),(x - edge.reg[1].x)) # AJS |
151 |
self.poly_vertices.append((edge.reg[0].sitenum,ang0,x,y)) # AJS |
152 |
self.poly_vertices.append((edge.reg[1].sitenum,ang1,x,y)) # AJS |
153 |
rid = -(len(self.extra_edges)-1) |
154 |
print lid,rid
|
155 |
return (lid,rid)
|
156 |
|
157 |
def line(self,x0,y0,x1,y1): |
158 |
pass
|
159 |
|
160 |
def outSite(self,s): |
161 |
if(self.debug): |
162 |
print "site (%d) at %f %f" % (s.sitenum, s.x, s.y) |
163 |
elif(self.triangulate): |
164 |
pass
|
165 |
elif(self.plot): |
166 |
self.circle (s.x, s.y, 3) #cradius) |
167 |
elif(self.doPrint): |
168 |
print "s %f %f" % (s.x, s.y) |
169 |
|
170 |
def outVertex(self,s,bot,top,mid,botang,topang,midang): # AJS |
171 |
self.poly_vertices.append((bot,botang,s.x,s.y)) # AJS |
172 |
self.poly_vertices.append((top,topang,s.x,s.y)) # AJS |
173 |
self.poly_vertices.append((mid,midang,s.x,s.y)) # AJS |
174 |
self.vertices.append((s.x,s.y))
|
175 |
if s.x < self.xmin: self.xmin = s.x |
176 |
elif s.x > self.xmax: self.xmax = s.x |
177 |
if s.y < self.ymin: self.ymin = s.y |
178 |
elif s.y > self.ymax: self.ymax = s.y |
179 |
|
180 |
if(self.debug): |
181 |
print "vertex(%d) at %f %f" % (s.sitenum, s.x, s.y) |
182 |
elif(self.triangulate): |
183 |
pass
|
184 |
elif(self.doPrint and not self.plot): |
185 |
print "v %f %f" % (s.x,s.y) |
186 |
|
187 |
def outTriple(self,s1,s2,s3): |
188 |
self.triangles.append((s1.sitenum, s2.sitenum, s3.sitenum))
|
189 |
if(self.debug): |
190 |
print "circle through left=%d right=%d bottom=%d" % (s1.sitenum, s2.sitenum, s3.sitenum) |
191 |
elif(self.triangulate and self.doPrint and not self.plot): |
192 |
print "%d %d %d" % (s1.sitenum, s2.sitenum, s3.sitenum) |
193 |
|
194 |
def outBisector(self,edge): |
195 |
self.lines.append((edge.a, edge.b, edge.c))
|
196 |
if(self.debug): |
197 |
print "line(%d) %gx+%gy=%g, bisecting %d %d" % (edge.edgenum, edge.a, edge.b, edge.c, edge.reg[0].sitenum, edge.reg[1].sitenum) |
198 |
elif(self.triangulate): |
199 |
if(self.plot): |
200 |
self.line(edge.reg[0].x, edge.reg[0].y, edge.reg[1].x, edge.reg[1].y) |
201 |
elif(self.doPrint and not self.plot): |
202 |
print "l %f %f %f" % (edge.a, edge.b, edge.c) |
203 |
|
204 |
def outEdge(self,edge): |
205 |
sitenumL = -1
|
206 |
if edge.ep[Edge.LE] is not None: |
207 |
sitenumL = edge.ep[Edge.LE].sitenum |
208 |
sitenumR = -1
|
209 |
if edge.ep[Edge.RE] is not None: |
210 |
sitenumR = edge.ep[Edge.RE].sitenum |
211 |
|
212 |
if (sitenumL == -1 or sitenumR == -1) and self.use_bound: # AJS (Logic error - OR condition need braces) |
213 |
sitenumL,sitenumR = self.clip_line(edge,sitenumL,sitenumR) # AJS |
214 |
|
215 |
self.edges.append((edge.edgenum,sitenumL,sitenumR))
|
216 |
if(not self.triangulate): |
217 |
if self.plot: |
218 |
self.clip_line(edge)
|
219 |
elif(self.doPrint): |
220 |
print "e %d" % edge.edgenum, |
221 |
print " %d " % sitenumL, |
222 |
print "%d" % sitenumR |
223 |
|
224 |
def set_bounds(self,bounds): |
225 |
if not bounds == None: |
226 |
self.xmin = bounds.xmin
|
227 |
self.ymin = bounds.ymin
|
228 |
self.xmax = bounds.xmax
|
229 |
self.ymax = bounds.ymax
|
230 |
else:
|
231 |
self.xmin = self.ymin = self.xmax = self.ymax = None |
232 |
|
233 |
|
234 |
#------------------------------------------------------------------
|
235 |
def voronoi(siteList,context): |
236 |
edgeList = EdgeList(siteList.xmin,siteList.xmax,len(siteList))
|
237 |
priorityQ = PriorityQueue(siteList.ymin,siteList.ymax,len(siteList))
|
238 |
siteIter = siteList.iterator() |
239 |
|
240 |
bottomsite = siteIter.next() |
241 |
context.outSite(bottomsite) |
242 |
newsite = siteIter.next() |
243 |
minpt = Site(-BIG_FLOAT,-BIG_FLOAT) |
244 |
while True: |
245 |
if not priorityQ.isEmpty(): |
246 |
minpt = priorityQ.getMinPt() |
247 |
|
248 |
if (newsite and (priorityQ.isEmpty() or cmp(newsite,minpt) < 0)): |
249 |
# newsite is smallest - this is a site event
|
250 |
context.outSite(newsite) |
251 |
|
252 |
# get first Halfedge to the LEFT and RIGHT of the new site
|
253 |
lbnd = edgeList.leftbnd(newsite) |
254 |
rbnd = lbnd.right |
255 |
|
256 |
# if this halfedge has no edge, bot = bottom site (whatever that is)
|
257 |
# create a new edge that bisects
|
258 |
bot = lbnd.rightreg(bottomsite) |
259 |
edge = Edge.bisect(bot,newsite) |
260 |
context.outBisector(edge) |
261 |
|
262 |
# create a new Halfedge, setting its pm field to 0 and insert
|
263 |
# this new bisector edge between the left and right vectors in
|
264 |
# a linked list
|
265 |
bisector = Halfedge(edge,Edge.LE) |
266 |
edgeList.insert(lbnd,bisector) |
267 |
|
268 |
# if the new bisector intersects with the left edge, remove
|
269 |
# the left edge's vertex, and put in the new one
|
270 |
p = lbnd.intersect(bisector) |
271 |
if p is not None: |
272 |
priorityQ.delete(lbnd) |
273 |
priorityQ.insert(lbnd,p,newsite.distance(p)) |
274 |
|
275 |
# create a new Halfedge, setting its pm field to 1
|
276 |
# insert the new Halfedge to the right of the original bisector
|
277 |
lbnd = bisector |
278 |
bisector = Halfedge(edge,Edge.RE) |
279 |
edgeList.insert(lbnd,bisector) |
280 |
|
281 |
# if this new bisector intersects with the right Halfedge
|
282 |
p = bisector.intersect(rbnd) |
283 |
if p is not None: |
284 |
# push the Halfedge into the ordered linked list of vertices
|
285 |
priorityQ.insert(bisector,p,newsite.distance(p)) |
286 |
|
287 |
newsite = siteIter.next() |
288 |
|
289 |
elif not priorityQ.isEmpty(): |
290 |
# intersection is smallest - this is a vector (circle) event
|
291 |
|
292 |
# pop the Halfedge with the lowest vector off the ordered list of
|
293 |
# vectors. Get the Halfedge to the left and right of the above HE
|
294 |
# and also the Halfedge to the right of the right HE
|
295 |
lbnd = priorityQ.popMinHalfedge() |
296 |
llbnd = lbnd.left |
297 |
rbnd = lbnd.right |
298 |
rrbnd = rbnd.right |
299 |
|
300 |
# get the Site to the left of the left HE and to the right of
|
301 |
# the right HE which it bisects
|
302 |
bot = lbnd.leftreg(bottomsite) |
303 |
top = rbnd.rightreg(bottomsite) |
304 |
|
305 |
# output the triple of sites, stating that a circle goes through them
|
306 |
mid = lbnd.rightreg(bottomsite) |
307 |
context.outTriple(bot,top,mid) |
308 |
|
309 |
# get the vertex that caused this event and set the vertex number
|
310 |
# couldn't do this earlier since we didn't know when it would be processed
|
311 |
v = lbnd.vertex |
312 |
siteList.setSiteNumber(v) |
313 |
botang = math.atan2((v.y - bot.y),(v.x - bot.x)) # AJS Calculate angles from vertex to source layer points
|
314 |
topang = math.atan2((v.y - top.y),(v.x - top.x)) # AJS
|
315 |
midang = math.atan2((v.y - mid.y),(v.x - mid.x)) # AJS
|
316 |
context.outVertex(v,bot.sitenum,top.sitenum,mid.sitenum,botang,topang,midang) # AJS
|
317 |
|
318 |
# set the endpoint of the left and right Halfedge to be this vector
|
319 |
if lbnd.edge.setEndpoint(lbnd.pm,v):
|
320 |
context.outEdge(lbnd.edge) |
321 |
|
322 |
if rbnd.edge.setEndpoint(rbnd.pm,v):
|
323 |
context.outEdge(rbnd.edge) |
324 |
|
325 |
|
326 |
# delete the lowest HE, remove all vertex events to do with the
|
327 |
# right HE and delete the right HE
|
328 |
edgeList.delete(lbnd) |
329 |
priorityQ.delete(rbnd) |
330 |
edgeList.delete(rbnd) |
331 |
|
332 |
|
333 |
# if the site to the left of the event is higher than the Site
|
334 |
# to the right of it, then swap them and set 'pm' to RIGHT
|
335 |
pm = Edge.LE |
336 |
if bot.y > top.y:
|
337 |
bot,top = top,bot |
338 |
pm = Edge.RE |
339 |
|
340 |
# Create an Edge (or line) that is between the two Sites. This
|
341 |
# creates the formula of the line, and assigns a line number to it
|
342 |
edge = Edge.bisect(bot, top) |
343 |
context.outBisector(edge) |
344 |
|
345 |
# create a HE from the edge
|
346 |
bisector = Halfedge(edge, pm) |
347 |
|
348 |
# insert the new bisector to the right of the left HE
|
349 |
# set one endpoint to the new edge to be the vector point 'v'
|
350 |
# If the site to the left of this bisector is higher than the right
|
351 |
# Site, then this endpoint is put in position 0; otherwise in pos 1
|
352 |
edgeList.insert(llbnd, bisector) |
353 |
if edge.setEndpoint(Edge.RE - pm, v):
|
354 |
context.outEdge(edge) |
355 |
|
356 |
# if left HE and the new bisector don't intersect, then delete
|
357 |
# the left HE, and reinsert it
|
358 |
p = llbnd.intersect(bisector) |
359 |
if p is not None: |
360 |
priorityQ.delete(llbnd); |
361 |
priorityQ.insert(llbnd, p, bot.distance(p)) |
362 |
|
363 |
# if right HE and the new bisector don't intersect, then reinsert it
|
364 |
p = bisector.intersect(rrbnd) |
365 |
if p is not None: |
366 |
priorityQ.insert(bisector, p, bot.distance(p)) |
367 |
else:
|
368 |
break
|
369 |
|
370 |
he = edgeList.leftend.right |
371 |
while he is not edgeList.rightend: |
372 |
context.outEdge(he.edge) |
373 |
he = he.right |
374 |
|
375 |
#------------------------------------------------------------------
|
376 |
def isEqual(a,b,relativeError=TOLERANCE): |
377 |
# is nearly equal to within the allowed relative error
|
378 |
norm = max(abs(a),abs(b)) |
379 |
return (norm < relativeError) or (abs(a - b) < (relativeError * norm)) |
380 |
|
381 |
#------------------------------------------------------------------
|
382 |
class Site(object): |
383 |
def __init__(self,x=0.0,y=0.0,sitenum=0): |
384 |
self.x = x
|
385 |
self.y = y
|
386 |
self.sitenum = sitenum
|
387 |
|
388 |
def dump(self): |
389 |
print "Site #%d (%g, %g)" % (self.sitenum,self.x,self.y) |
390 |
|
391 |
def __cmp__(self,other): |
392 |
if self.y < other.y: |
393 |
return -1 |
394 |
elif self.y > other.y: |
395 |
return 1 |
396 |
elif self.x < other.x: |
397 |
return -1 |
398 |
elif self.x > other.x: |
399 |
return 1 |
400 |
else:
|
401 |
return 0 |
402 |
|
403 |
def distance(self,other): |
404 |
dx = self.x - other.x
|
405 |
dy = self.y - other.y
|
406 |
return math.sqrt(dx*dx + dy*dy)
|
407 |
|
408 |
#------------------------------------------------------------------
|
409 |
class Edge(object): |
410 |
LE = 0
|
411 |
RE = 1
|
412 |
EDGE_NUM = 0
|
413 |
DELETED = {} # marker value
|
414 |
|
415 |
def __init__(self): |
416 |
self.a = 0.0 |
417 |
self.b = 0.0 |
418 |
self.c = 0.0 |
419 |
self.ep = [None,None] |
420 |
self.reg = [None,None] |
421 |
self.edgenum = 0 |
422 |
|
423 |
def dump(self): |
424 |
print "(#%d a=%g, b=%g, c=%g)" % (self.edgenum,self.a,self.b,self.c) |
425 |
print "ep",self.ep |
426 |
print "reg",self.reg |
427 |
|
428 |
def setEndpoint(self, lrFlag, site): |
429 |
self.ep[lrFlag] = site
|
430 |
if self.ep[Edge.RE - lrFlag] is None: |
431 |
return False |
432 |
return True |
433 |
|
434 |
@staticmethod
|
435 |
def bisect(s1,s2): |
436 |
newedge = Edge() |
437 |
newedge.reg[0] = s1 # store the sites that this edge is bisecting |
438 |
newedge.reg[1] = s2
|
439 |
|
440 |
# to begin with, there are no endpoints on the bisector - it goes to infinity
|
441 |
# ep[0] and ep[1] are None
|
442 |
|
443 |
# get the difference in x dist between the sites
|
444 |
dx = float(s2.x - s1.x)
|
445 |
dy = float(s2.y - s1.y)
|
446 |
adx = abs(dx) # make sure that the difference in positive |
447 |
ady = abs(dy)
|
448 |
|
449 |
# get the slope of the line
|
450 |
newedge.c = float(s1.x * dx + s1.y * dy + (dx*dx + dy*dy)*0.5) |
451 |
if adx > ady :
|
452 |
# set formula of line, with x fixed to 1
|
453 |
newedge.a = 1.0
|
454 |
newedge.b = dy/dx |
455 |
newedge.c /= dx |
456 |
else:
|
457 |
# set formula of line, with y fixed to 1
|
458 |
newedge.b = 1.0
|
459 |
newedge.a = dx/dy |
460 |
newedge.c /= dy |
461 |
|
462 |
newedge.edgenum = Edge.EDGE_NUM |
463 |
Edge.EDGE_NUM += 1
|
464 |
return newedge
|
465 |
|
466 |
|
467 |
#------------------------------------------------------------------
|
468 |
class Halfedge(object): |
469 |
def __init__(self,edge=None,pm=Edge.LE): |
470 |
self.left = None # left Halfedge in the edge list |
471 |
self.right = None # right Halfedge in the edge list |
472 |
self.qnext = None # priority queue linked list pointer |
473 |
self.edge = edge # edge list Edge |
474 |
self.pm = pm
|
475 |
self.vertex = None # Site() |
476 |
self.ystar = BIG_FLOAT
|
477 |
|
478 |
def dump(self): |
479 |
print "Halfedge--------------------------" |
480 |
print "left: ", self.left |
481 |
print "right: ", self.right |
482 |
print "edge: ", self.edge |
483 |
print "pm: ", self.pm |
484 |
print "vertex: ", |
485 |
if self.vertex: self.vertex.dump() |
486 |
else: print "None" |
487 |
print "ystar: ", self.ystar |
488 |
|
489 |
|
490 |
def __cmp__(self,other): |
491 |
if self.ystar > other.ystar: |
492 |
return 1 |
493 |
elif self.ystar < other.ystar: |
494 |
return -1 |
495 |
elif self.vertex.x > other.vertex.x: |
496 |
return 1 |
497 |
elif self.vertex.x < other.vertex.x: |
498 |
return -1 |
499 |
else:
|
500 |
return 0 |
501 |
|
502 |
def leftreg(self,default): |
503 |
if not self.edge: |
504 |
return default
|
505 |
elif self.pm == Edge.LE: |
506 |
return self.edge.reg[Edge.LE] |
507 |
else:
|
508 |
return self.edge.reg[Edge.RE] |
509 |
|
510 |
def rightreg(self,default): |
511 |
if not self.edge: |
512 |
return default
|
513 |
elif self.pm == Edge.LE: |
514 |
return self.edge.reg[Edge.RE] |
515 |
else:
|
516 |
return self.edge.reg[Edge.LE] |
517 |
|
518 |
|
519 |
# returns True if p is to right of halfedge self
|
520 |
def isPointRightOf(self,pt): |
521 |
e = self.edge
|
522 |
topsite = e.reg[1]
|
523 |
right_of_site = pt.x > topsite.x |
524 |
|
525 |
if(right_of_site and self.pm == Edge.LE): |
526 |
return True |
527 |
|
528 |
if(not right_of_site and self.pm == Edge.RE): |
529 |
return False |
530 |
|
531 |
if(e.a == 1.0): |
532 |
dyp = pt.y - topsite.y |
533 |
dxp = pt.x - topsite.x |
534 |
fast = 0;
|
535 |
if ((not right_of_site and e.b < 0.0) or (right_of_site and e.b >= 0.0)): |
536 |
above = dyp >= e.b * dxp |
537 |
fast = above |
538 |
else:
|
539 |
above = pt.x + pt.y * e.b > e.c |
540 |
if(e.b < 0.0): |
541 |
above = not above
|
542 |
if (not above): |
543 |
fast = 1
|
544 |
if (not fast): |
545 |
dxs = topsite.x - (e.reg[0]).x
|
546 |
above = e.b * (dxp*dxp - dyp*dyp) < dxs*dyp*(1.0+2.0*dxp/dxs + e.b*e.b) |
547 |
if(e.b < 0.0): |
548 |
above = not above
|
549 |
else: # e.b == 1.0 |
550 |
yl = e.c - e.a * pt.x |
551 |
t1 = pt.y - yl |
552 |
t2 = pt.x - topsite.x |
553 |
t3 = yl - topsite.y |
554 |
above = t1*t1 > t2*t2 + t3*t3 |
555 |
|
556 |
if(self.pm==Edge.LE): |
557 |
return above
|
558 |
else:
|
559 |
return not above |
560 |
|
561 |
#--------------------------
|
562 |
# create a new site where the Halfedges el1 and el2 intersect
|
563 |
def intersect(self,other): |
564 |
e1 = self.edge
|
565 |
e2 = other.edge |
566 |
if (e1 is None) or (e2 is None): |
567 |
return None |
568 |
|
569 |
# if the two edges bisect the same parent return None
|
570 |
if e1.reg[1] is e2.reg[1]: |
571 |
return None |
572 |
|
573 |
d = e1.a * e2.b - e1.b * e2.a |
574 |
if isEqual(d,0.0): |
575 |
return None |
576 |
|
577 |
xint = (e1.c*e2.b - e2.c*e1.b) / d |
578 |
yint = (e2.c*e1.a - e1.c*e2.a) / d |
579 |
if(cmp(e1.reg[1],e2.reg[1]) < 0): |
580 |
he = self
|
581 |
e = e1 |
582 |
else:
|
583 |
he = other |
584 |
e = e2 |
585 |
|
586 |
rightOfSite = xint >= e.reg[1].x
|
587 |
if((rightOfSite and he.pm == Edge.LE) or |
588 |
(not rightOfSite and he.pm == Edge.RE)): |
589 |
return None |
590 |
|
591 |
# create a new site at the point of intersection - this is a new
|
592 |
# vector event waiting to happen
|
593 |
return Site(xint,yint)
|
594 |
|
595 |
|
596 |
|
597 |
#------------------------------------------------------------------
|
598 |
class EdgeList(object): |
599 |
def __init__(self,xmin,xmax,nsites): |
600 |
if xmin > xmax: xmin,xmax = xmax,xmin
|
601 |
self.hashsize = int(2*math.sqrt(nsites+4)) |
602 |
|
603 |
self.xmin = xmin
|
604 |
self.deltax = float(xmax - xmin) |
605 |
self.hash = [None]*self.hashsize |
606 |
|
607 |
self.leftend = Halfedge()
|
608 |
self.rightend = Halfedge()
|
609 |
self.leftend.right = self.rightend |
610 |
self.rightend.left = self.leftend |
611 |
self.hash[0] = self.leftend |
612 |
self.hash[-1] = self.rightend |
613 |
|
614 |
def insert(self,left,he): |
615 |
he.left = left |
616 |
he.right = left.right |
617 |
left.right.left = he |
618 |
left.right = he |
619 |
|
620 |
def delete(self,he): |
621 |
he.left.right = he.right |
622 |
he.right.left = he.left |
623 |
he.edge = Edge.DELETED |
624 |
|
625 |
# Get entry from hash table, pruning any deleted nodes
|
626 |
def gethash(self,b): |
627 |
if(b < 0 or b >= self.hashsize): |
628 |
return None |
629 |
he = self.hash[b]
|
630 |
if he is None or he.edge is not Edge.DELETED: |
631 |
return he
|
632 |
|
633 |
# Hash table points to deleted half edge. Patch as necessary.
|
634 |
self.hash[b] = None |
635 |
return None |
636 |
|
637 |
def leftbnd(self,pt): |
638 |
# Use hash table to get close to desired halfedge
|
639 |
bucket = int(((pt.x - self.xmin)/self.deltax * self.hashsize)) |
640 |
|
641 |
if(bucket < 0): |
642 |
bucket =0;
|
643 |
|
644 |
if(bucket >=self.hashsize): |
645 |
bucket = self.hashsize-1 |
646 |
|
647 |
he = self.gethash(bucket)
|
648 |
if(he is None): |
649 |
i = 1
|
650 |
while True: |
651 |
he = self.gethash(bucket-i)
|
652 |
if (he is not None): break; |
653 |
he = self.gethash(bucket+i)
|
654 |
if (he is not None): break; |
655 |
i += 1
|
656 |
|
657 |
# Now search linear list of halfedges for the corect one
|
658 |
if (he is self.leftend) or (he is not self.rightend and he.isPointRightOf(pt)): |
659 |
he = he.right |
660 |
while he is not self.rightend and he.isPointRightOf(pt): |
661 |
he = he.right |
662 |
he = he.left; |
663 |
else:
|
664 |
he = he.left |
665 |
while (he is not self.leftend and not he.isPointRightOf(pt)): |
666 |
he = he.left |
667 |
|
668 |
# Update hash table and reference counts
|
669 |
if(bucket > 0 and bucket < self.hashsize-1): |
670 |
self.hash[bucket] = he
|
671 |
return he
|
672 |
|
673 |
|
674 |
#------------------------------------------------------------------
|
675 |
class PriorityQueue(object): |
676 |
def __init__(self,ymin,ymax,nsites): |
677 |
self.ymin = ymin
|
678 |
self.deltay = ymax - ymin
|
679 |
self.hashsize = int(4 * math.sqrt(nsites)) |
680 |
self.count = 0 |
681 |
self.minidx = 0 |
682 |
self.hash = []
|
683 |
for i in range(self.hashsize): |
684 |
self.hash.append(Halfedge())
|
685 |
|
686 |
def __len__(self): |
687 |
return self.count |
688 |
|
689 |
def isEmpty(self): |
690 |
return self.count == 0 |
691 |
|
692 |
def insert(self,he,site,offset): |
693 |
he.vertex = site |
694 |
he.ystar = site.y + offset |
695 |
last = self.hash[self.getBucket(he)] |
696 |
next = last.qnext |
697 |
while((next is not None) and cmp(he,next) > 0): |
698 |
last = next
|
699 |
next = last.qnext |
700 |
he.qnext = last.qnext |
701 |
last.qnext = he |
702 |
self.count += 1 |
703 |
|
704 |
def delete(self,he): |
705 |
if (he.vertex is not None): |
706 |
last = self.hash[self.getBucket(he)] |
707 |
while last.qnext is not he: |
708 |
last = last.qnext |
709 |
last.qnext = he.qnext |
710 |
self.count -= 1 |
711 |
he.vertex = None
|
712 |
|
713 |
def getBucket(self,he): |
714 |
bucket = int(((he.ystar - self.ymin) / self.deltay) * self.hashsize) |
715 |
if bucket < 0: bucket = 0 |
716 |
if bucket >= self.hashsize: bucket = self.hashsize-1 |
717 |
if bucket < self.minidx: self.minidx = bucket |
718 |
return bucket
|
719 |
|
720 |
def getMinPt(self): |
721 |
while(self.hash[self.minidx].qnext is None): |
722 |
self.minidx += 1 |
723 |
he = self.hash[self.minidx].qnext |
724 |
x = he.vertex.x |
725 |
y = he.ystar |
726 |
return Site(x,y)
|
727 |
|
728 |
def popMinHalfedge(self): |
729 |
curr = self.hash[self.minidx].qnext |
730 |
self.hash[self.minidx].qnext = curr.qnext |
731 |
self.count -= 1 |
732 |
return curr
|
733 |
|
734 |
|
735 |
#------------------------------------------------------------------
|
736 |
class SiteList(object): |
737 |
def __init__(self,pointList): |
738 |
self.__sites = []
|
739 |
self.__sitenum = 0 |
740 |
|
741 |
self.__xmin = pointList[0].x() |
742 |
self.__ymin = pointList[0].y() |
743 |
self.__xmax = pointList[0].x() |
744 |
self.__ymax = pointList[0].y() |
745 |
for i,pt in enumerate(pointList): |
746 |
self.__sites.append(Site(pt.x(),pt.y(),i))
|
747 |
if pt.x() < self.__xmin: self.__xmin = pt.x() |
748 |
if pt.y() < self.__ymin: self.__ymin = pt.y() |
749 |
if pt.x() > self.__xmax: self.__xmax = pt.x() |
750 |
if pt.y() > self.__ymax: self.__ymax = pt.y() |
751 |
self.__sites.sort()
|
752 |
|
753 |
def setSiteNumber(self,site): |
754 |
site.sitenum = self.__sitenum
|
755 |
self.__sitenum += 1 |
756 |
|
757 |
class Iterator(object): |
758 |
def __init__(this,lst): this.generator = (s for s in lst) |
759 |
def __iter__(this): return this |
760 |
def next(this): |
761 |
try:
|
762 |
return this.generator.next()
|
763 |
except StopIteration: |
764 |
return None |
765 |
|
766 |
def iterator(self): |
767 |
return SiteList.Iterator(self.__sites) |
768 |
|
769 |
def __iter__(self): |
770 |
return SiteList.Iterator(self.__sites) |
771 |
|
772 |
def __len__(self): |
773 |
return len(self.__sites) |
774 |
|
775 |
def _getxmin(self): return self.__xmin |
776 |
def _getymin(self): return self.__ymin |
777 |
def _getxmax(self): return self.__xmax |
778 |
def _getymax(self): return self.__ymax |
779 |
xmin = property(_getxmin)
|
780 |
ymin = property(_getymin)
|
781 |
xmax = property(_getxmax)
|
782 |
ymax = property(_getymax)
|
783 |
|
784 |
|
785 |
#------------------------------------------------------------------
|
786 |
def computeVoronoiDiagram( points ): |
787 |
""" Takes a list of point objects (which must have x and y fields).
|
788 |
Returns a 3-tuple of:
|
789 |
|
790 |
(1) a list of 2-tuples, which are the x,y coordinates of the
|
791 |
Voronoi diagram vertices
|
792 |
(2) a list of 3-tuples (a,b,c) which are the equations of the
|
793 |
lines in the Voronoi diagram: a*x + b*y = c
|
794 |
(3) a list of 3-tuples, (l, v1, v2) representing edges of the
|
795 |
Voronoi diagram. l is the index of the line, v1 and v2 are
|
796 |
the indices of the vetices at the end of the edge. If
|
797 |
v1 or v2 is -1, the line extends to infinity.
|
798 |
"""
|
799 |
siteList = SiteList( points ) |
800 |
context = Context() |
801 |
context.use_bound = True # AJS |
802 |
context.set_bounds( siteList ) |
803 |
voronoi( siteList, context ) |
804 |
return ( context.vertices, context.lines, context.edges, (context.xmin,context.ymin,context.xmax,context.ymax),context.poly_vertices) # AJS |
805 |
|
806 |
#------------------------------------------------------------------
|
807 |
def computeDelaunayTriangulation( points ): |
808 |
""" Takes a list of point objects (which must have x and y fields).
|
809 |
Returns a list of 3-tuples: the indices of the points that form a
|
810 |
Delaunay triangle.
|
811 |
"""
|
812 |
siteList = SiteList( points ) |
813 |
context = Context() |
814 |
context.triangulate = True
|
815 |
voronoi( siteList, context ) |
816 |
return context.triangles
|
817 |
|
818 |
#-----------------------------------------------------------------------------
|
819 |
if __name__=="__main__": |
820 |
try:
|
821 |
optlist,args = getopt.getopt(sys.argv[1:],"thdp") |
822 |
except getopt.GetoptError:
|
823 |
usage() |
824 |
sys.exit(2)
|
825 |
|
826 |
doHelp = 0
|
827 |
c = Context() |
828 |
c.doPrint = 1
|
829 |
for opt in optlist: |
830 |
if opt[0] == "-d": c.debug = 1 |
831 |
if opt[0] == "-p": c.plot = 1 |
832 |
if opt[0] == "-t": c.triangulate = 1 |
833 |
if opt[0] == "-h": doHelp = 1 |
834 |
|
835 |
if not doHelp: |
836 |
pts = [] |
837 |
fp = sys.stdin |
838 |
if len(args) > 0: |
839 |
fp = open(args[0],'r') |
840 |
for line in fp: |
841 |
fld = line.split() |
842 |
x = float(fld[0]) |
843 |
y = float(fld[1]) |
844 |
pts.append(Site(x,y)) |
845 |
if len(args) > 0: fp.close() |
846 |
|
847 |
if doHelp or len(pts) == 0: |
848 |
usage() |
849 |
sys.exit(2)
|
850 |
|
851 |
sl = SiteList(pts) |
852 |
voronoi(sl,c) |
853 |
|