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