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