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 return ( 392 self.y < other.y 393 or self.y <= other.y and self.x < other.x 394 ) 395 396 def distance(self, other): 397 dx = self.x - other.x 398 dy = self.y - other.y 399 return math.sqrt(dx * dx + dy * dy) 400 401 402# ------------------------------------------------------------------ 403 404 405class Edge(object): 406 LE = 0 407 RE = 1 408 EDGE_NUM = 0 409 DELETED = {} # marker value 410 411 def __init__(self): 412 self.a = 0.0 413 self.b = 0.0 414 self.c = 0.0 415 self.ep = [None, None] 416 self.reg = [None, None] 417 self.edgenum = 0 418 419 def dump(self): 420 # fix_print_with_import 421 print("(#%d a=%g, b=%g, c=%g)" % (self.edgenum, self.a, self.b, self.c)) 422 # fix_print_with_import 423 print("ep", self.ep) 424 # fix_print_with_import 425 print("reg", self.reg) 426 427 def setEndpoint(self, lrFlag, site): 428 self.ep[lrFlag] = site 429 return self.ep[Edge.RE - lrFlag] is not None 430 431 @staticmethod 432 def bisect(s1, s2): 433 newedge = Edge() 434 newedge.reg[0] = s1 # store the sites that this edge is bisecting 435 newedge.reg[1] = s2 436 437 # to begin with, there are no endpoints on the bisector - it goes to infinity 438 # ep[0] and ep[1] are None 439 440 # get the difference in x dist between the sites 441 dx = float(s2.x - s1.x) 442 dy = float(s2.y - s1.y) 443 adx = abs(dx) # make sure that the difference in positive 444 ady = abs(dy) 445 446 # get the slope of the line 447 newedge.c = float(s1.x * dx + s1.y * dy + (dx**2 + dy**2) * 0.5) 448 if adx > ady: 449 # set formula of line, with x fixed to 1 450 newedge.a = 1.0 451 newedge.b = dy / dx 452 newedge.c /= dx 453 else: 454 # set formula of line, with y fixed to 1 455 newedge.b = 1.0 456 newedge.a = dx / dy 457 newedge.c /= dy 458 459 newedge.edgenum = Edge.EDGE_NUM 460 Edge.EDGE_NUM += 1 461 return newedge 462 463 464# ------------------------------------------------------------------ 465class Halfedge(object): 466 467 def __init__(self, edge=None, pm=Edge.LE): 468 self.left = None # left Halfedge in the edge list 469 self.right = None # right Halfedge in the edge list 470 self.qnext = None # priority queue linked list pointer 471 self.edge = edge # edge list Edge 472 self.pm = pm 473 self.vertex = None # Site() 474 self.ystar = BIG_FLOAT 475 476 def dump(self): 477 # fix_print_with_import 478 print("Halfedge--------------------------") 479 # fix_print_with_import 480 print("left: ", self.left) 481 # fix_print_with_import 482 print("right: ", self.right) 483 # fix_print_with_import 484 print("edge: ", self.edge) 485 # fix_print_with_import 486 print("pm: ", self.pm) 487 # fix_print_with_import 488 print("vertex:") 489 if self.vertex: 490 self.vertex.dump() 491 else: 492 # fix_print_with_import 493 print("None") 494 # fix_print_with_import 495 print("ystar: ", self.ystar) 496 497 def __eq__(self, other): 498 return (self.vertex.x == other.vertex.x) and (self.ystar == other.ystar) 499 500 def __lt__(self, other): 501 return ( 502 self.ystar < other.ystar 503 or self.ystar <= other.ystar 504 and self.vertex.x < other.vertex.x 505 ) 506 507 def leftreg(self, default): 508 if not self.edge: 509 return default 510 elif self.pm == Edge.LE: 511 return self.edge.reg[Edge.LE] 512 else: 513 return self.edge.reg[Edge.RE] 514 515 def rightreg(self, default): 516 if not self.edge: 517 return default 518 elif self.pm == Edge.LE: 519 return self.edge.reg[Edge.RE] 520 else: 521 return self.edge.reg[Edge.LE] 522 523 # returns True if p is to right of halfedge self 524 def isPointRightOf(self, pt): 525 e = self.edge 526 topsite = e.reg[1] 527 right_of_site = pt.x > topsite.x 528 529 if (right_of_site and self.pm == Edge.LE): 530 return True 531 532 if (not right_of_site and self.pm == Edge.RE): 533 return False 534 535 if (e.a == 1.0): 536 dyp = pt.y - topsite.y 537 dxp = pt.x - topsite.x 538 fast = 0 539 if ((not right_of_site and e.b < 0.0) or (right_of_site and e.b >= 0.0)): 540 above = dyp >= e.b * dxp 541 fast = above 542 else: 543 above = pt.x + pt.y * e.b > e.c 544 if (e.b < 0.0): 545 above = not above 546 if (not above): 547 fast = 1 548 if (not fast): 549 dxs = topsite.x - (e.reg[0]).x 550 above = e.b * (dxp * dxp - dyp * dyp) < dxs * dyp * (1.0 + 2.0 * dxp / dxs + e.b * e.b) 551 if (e.b < 0.0): 552 above = not above 553 else: # e.b == 1.0 554 yl = e.c - e.a * pt.x 555 t1 = pt.y - yl 556 t2 = pt.x - topsite.x 557 t3 = yl - topsite.y 558 above = t1 * t1 > t2 * t2 + t3 * t3 559 560 if (self.pm == Edge.LE): 561 return above 562 else: 563 return not above 564 565 # -------------------------- 566 # create a new site where the Halfedges el1 and el2 intersect 567 def intersect(self, other): 568 e1 = self.edge 569 e2 = other.edge 570 if (e1 is None) or (e2 is None): 571 return None 572 573 # if the two edges bisect the same parent return None 574 if e1.reg[1] is e2.reg[1]: 575 return None 576 577 d = e1.a * e2.b - e1.b * e2.a 578 if isEqual(d, 0.0): 579 return None 580 581 xint = (e1.c * e2.b - e2.c * e1.b) / d 582 yint = (e2.c * e1.a - e1.c * e2.a) / d 583 if (cmp(e1.reg[1], e2.reg[1]) < 0): 584 he = self 585 e = e1 586 else: 587 he = other 588 e = e2 589 590 rightOfSite = xint >= e.reg[1].x 591 if ((rightOfSite and he.pm == Edge.LE) 592 or (not rightOfSite and he.pm == Edge.RE)): 593 return None 594 595 # create a new site at the point of intersection - this is a new 596 # vector event waiting to happen 597 return Site(xint, yint) 598 599 600# ------------------------------------------------------------------ 601class EdgeList(object): 602 603 def __init__(self, xmin, xmax, nsites): 604 if xmin > xmax: 605 xmin, xmax = xmax, xmin 606 self.hashsize = int(2 * math.sqrt(nsites + 4)) 607 608 self.xmin = xmin 609 self.deltax = float(xmax - xmin) 610 self.hash = [None] * self.hashsize 611 612 self.leftend = Halfedge() 613 self.rightend = Halfedge() 614 self.leftend.right = self.rightend 615 self.rightend.left = self.leftend 616 self.hash[0] = self.leftend 617 self.hash[-1] = self.rightend 618 619 def insert(self, left, he): 620 he.left = left 621 he.right = left.right 622 left.right.left = he 623 left.right = he 624 625 def delete(self, he): 626 he.left.right = he.right 627 he.right.left = he.left 628 he.edge = Edge.DELETED 629 630 # Get entry from hash table, pruning any deleted nodes 631 def gethash(self, b): 632 if (b < 0 or b >= self.hashsize): 633 return None 634 he = self.hash[b] 635 if he is None or he.edge is not Edge.DELETED: 636 return he 637 638 # Hash table points to deleted half edge. Patch as necessary. 639 self.hash[b] = None 640 return None 641 642 def leftbnd(self, pt): 643 # Use hash table to get close to desired halfedge 644 bucket = int(((pt.x - self.xmin) / self.deltax * self.hashsize)) 645 646 bucket = max(bucket, 0) 647 648 if (bucket >= self.hashsize): 649 bucket = self.hashsize - 1 650 651 he = self.gethash(bucket) 652 if (he is None): 653 i = 1 654 while True: 655 he = self.gethash(bucket - i) 656 if (he is not None): 657 break 658 he = self.gethash(bucket + i) 659 if (he is not None): 660 break 661 i += 1 662 663 # Now search linear list of halfedges for the correct one 664 if (he is self.leftend) or (he is not self.rightend and he.isPointRightOf(pt)): 665 he = he.right 666 while he is not self.rightend and he.isPointRightOf(pt): 667 he = he.right 668 he = he.left 669 else: 670 he = he.left 671 while (he is not self.leftend and not he.isPointRightOf(pt)): 672 he = he.left 673 674 # Update hash table and reference counts 675 if (bucket > 0 and bucket < self.hashsize - 1): 676 self.hash[bucket] = he 677 return he 678 679 680# ------------------------------------------------------------------ 681class PriorityQueue(object): 682 683 def __init__(self, ymin, ymax, nsites): 684 self.ymin = ymin 685 self.deltay = ymax - ymin 686 self.hashsize = int(4 * math.sqrt(nsites)) 687 self.count = 0 688 self.minidx = 0 689 690 self.hash = [ 691 Halfedge() for _ in range(self.hashsize) 692 ] 693 694 def __len__(self): 695 return self.count 696 697 def isEmpty(self): 698 return self.count == 0 699 700 def insert(self, he, site, offset): 701 he.vertex = site 702 he.ystar = site.y + offset 703 last = self.hash[self.getBucket(he)] 704 next = last.qnext 705 while ((next is not None) and cmp(he, next) > 0): 706 last = next 707 next = last.qnext 708 he.qnext = last.qnext 709 last.qnext = he 710 self.count += 1 711 712 def delete(self, he): 713 if (he.vertex is not None): 714 last = self.hash[self.getBucket(he)] 715 while last.qnext is not he: 716 last = last.qnext 717 last.qnext = he.qnext 718 self.count -= 1 719 he.vertex = None 720 721 def getBucket(self, he): 722 bucket = int(((he.ystar - self.ymin) / self.deltay) * self.hashsize) 723 bucket = max(bucket, 0) 724 if bucket >= self.hashsize: 725 bucket = self.hashsize - 1 726 if bucket < self.minidx: 727 self.minidx = bucket 728 return bucket 729 730 def getMinPt(self): 731 while (self.hash[self.minidx].qnext is None): 732 self.minidx += 1 733 he = self.hash[self.minidx].qnext 734 x = he.vertex.x 735 y = he.ystar 736 return Site(x, y) 737 738 def popMinHalfedge(self): 739 curr = self.hash[self.minidx].qnext 740 self.hash[self.minidx].qnext = curr.qnext 741 self.count -= 1 742 return curr 743 744 745# ------------------------------------------------------------------ 746class SiteList(object): 747 748 def __init__(self, pointList): 749 self.__sites = [] 750 self.__sitenum = 0 751 752 self.__xmin = pointList[0].x 753 self.__ymin = pointList[0].y 754 self.__xmax = pointList[0].x 755 self.__ymax = pointList[0].y 756 for i, pt in enumerate(pointList): 757 self.__sites.append(Site(pt.x, pt.y, i)) 758 if pt.x < self.__xmin: 759 self.__xmin = pt.x 760 if pt.y < self.__ymin: 761 self.__ymin = pt.y 762 if pt.x > self.__xmax: 763 self.__xmax = pt.x 764 if pt.y > self.__ymax: 765 self.__ymax = pt.y 766 self.__sites.sort() 767 768 def setSiteNumber(self, site): 769 site.sitenum = self.__sitenum 770 self.__sitenum += 1 771 772 class Iterator(object): 773 774 def __init__(this, lst): 775 this.generator = iter(lst) 776 777 def __iter__(this): 778 return this 779 780 def __next__(this): 781 try: 782 return next(this.generator) 783 except StopIteration: 784 return None 785 786 def iterator(self): 787 return SiteList.Iterator(self.__sites) 788 789 def __iter__(self): 790 return SiteList.Iterator(self.__sites) 791 792 def __len__(self): 793 return len(self.__sites) 794 795 def _getxmin(self): 796 return self.__xmin 797 798 def _getymin(self): 799 return self.__ymin 800 801 def _getxmax(self): 802 return self.__xmax 803 804 def _getymax(self): 805 return self.__ymax 806 807 xmin = property(_getxmin) 808 ymin = property(_getymin) 809 xmax = property(_getxmax) 810 ymax = property(_getymax) 811 812 813# ------------------------------------------------------------------ 814def computeVoronoiDiagram(points): 815 """ Takes a list of point objects (which must have x and y fields). 816 Returns a 3-tuple of: 817 818 (1) a list of 2-tuples, which are the x,y coordinates of the 819 Voronoi diagram vertices 820 (2) a list of 3-tuples (a,b,c) which are the equations of the 821 lines in the Voronoi diagram: a*x + b*y = c 822 (3) a list of 3-tuples, (l, v1, v2) representing edges of the 823 Voronoi diagram. l is the index of the line, v1 and v2 are 824 the indices of the vetices at the end of the edge. If 825 v1 or v2 is -1, the line extends to infinity. 826 """ 827 siteList = SiteList(points) 828 context = Context() 829 voronoi(siteList, context) 830 return (context.vertices, context.lines, context.edges) 831 832 833# ------------------------------------------------------------------ 834 835 836def computeDelaunayTriangulation(points): 837 """ Takes a list of point objects (which must have x and y fields). 838 Returns a list of 3-tuples: the indices of the points that form a 839 Delaunay triangle. 840 """ 841 siteList = SiteList(points) 842 context = Context() 843 context.triangulate = True 844 voronoi(siteList, context) 845 return context.triangles 846 847 848# ----------------------------------------------------------------------------- 849if __name__ == "__main__": 850 try: 851 optlist, args = getopt.getopt(sys.argv[1:], "thdp") 852 except getopt.GetoptError: 853 usage() 854 sys.exit(2) 855 856 doHelp = 0 857 c = Context() 858 c.doPrint = 1 859 for opt in optlist: 860 if opt[0] == "-d": 861 c.debug = 1 862 if opt[0] == "-p": 863 c.plot = 1 864 if opt[0] == "-t": 865 c.triangulate = 1 866 if opt[0] == "-h": 867 doHelp = 1 868 869 if not doHelp: 870 pts = [] 871 fp = sys.stdin 872 if len(args) > 0: 873 fp = open(args[0], 'r') 874 for line in fp: 875 fld = line.split() 876 x = float(fld[0]) 877 y = float(fld[1]) 878 pts.append(Site(x, y)) 879 if len(args) > 0: 880 fp.close() 881 882 if doHelp or len(pts) == 0: 883 usage() 884 sys.exit(2) 885 886 sl = SiteList(pts) 887 voronoi(sl, c) 888 889 890def cmp(a, b): 891 """Compare the two objects x and y and return an integer according to the 892 outcome. The return value is negative if x < y, zero if x == y and strictly 893 positive if x > y. 894 895 In python 2 cmp() was a built in function but in python 3 is gone. 896 """ 897 return (b < a) - (a < b) 898