share/extensions/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
 114     def circle(self,x,y,rad):
 115         pass
 116
 117     def clip_line(self,edge):
 118         pass
 119
 120     def line(self,x0,y0,x1,y1):
 121         pass
 122
 123     def outSite(self,s):
 124         if(self.debug):
 125             print "site (%d) at %f %f" % (s.sitenum, s.x, s.y)
 126         elif(self.triangulate):
 127             pass
 128         elif(self.plot):
 129             self.circle (s.x, s.y, cradius)
 130         elif(self.doPrint):
 131             print "s %f %f" % (s.x, s.y)
 132
 133     def outVertex(self,s):
 134         self.vertices.append((s.x,s.y))
 135         if(self.debug):
 136             print  "vertex(%d) at %f %f" % (s.sitenum, s.x, s.y)
 137         elif(self.triangulate):
 138             pass
 139         elif(self.doPrint and not self.plot):
 140             print "v %f %f" % (s.x,s.y)
 141
 142     def outTriple(self,s1,s2,s3):
 143         self.triangles.append((s1.sitenum, s2.sitenum, s3.sitenum))
 144         if(self.debug):
 145             print "circle through left=%d right=%d bottom=%d" % (s1.sitenum, s2.sitenum, s3.sitenum)
 146         elif(self.triangulate and self.doPrint and not self.plot):
 147             print "%d %d %d" % (s1.sitenum, s2.sitenum, s3.sitenum)
 148
 149     def outBisector(self,edge):
 150         self.lines.append((edge.a, edge.b, edge.c))
 151         if(self.debug):
 152             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)
 153         elif(self.triangulate):
 154             if(self.plot):
 155                 self.line(edge.reg[0].x, edge.reg[0].y, edge.reg[1].x, edge.reg[1].y)
 156         elif(self.doPrint and not self.plot):
 157             print "l %f %f %f" % (edge.a, edge.b, edge.c)
 158
 159     def outEdge(self,edge):
 160         sitenumL = -1
 161         if edge.ep[Edge.LE] is not None:
 162             sitenumL = edge.ep[Edge.LE].sitenum
 163         sitenumR = -1
 164         if edge.ep[Edge.RE] is not None:
 165             sitenumR = edge.ep[Edge.RE].sitenum
 166         self.edges.append((edge.edgenum,sitenumL,sitenumR))
 167         if(not self.triangulate):
 168             if self.plot:
 169                 self.clip_line(edge)
 170             elif(self.doPrint):
 171                 print "e %d" % edge.edgenum,
 172                 print " %d " % sitenumL,
 173                 print "%d" % sitenumR
 174
 175 #------------------------------------------------------------------
 176 def voronoi(siteList,context):
 177     edgeList  = EdgeList(siteList.xmin,siteList.xmax,len(siteList))
 178     priorityQ = PriorityQueue(siteList.ymin,siteList.ymax,len(siteList))
 179     siteIter = siteList.iterator()
 180
 181     bottomsite = siteIter.next()
 182     context.outSite(bottomsite)
 183     newsite = siteIter.next()
 184     minpt = Site(-BIG_FLOAT,-BIG_FLOAT)
 185     while True:
 186         if not priorityQ.isEmpty():
 187             minpt = priorityQ.getMinPt()
 188
 189         if (newsite and (priorityQ.isEmpty() or cmp(newsite,minpt) < 0)):
 190             # newsite is smallest -  this is a site event
 191             context.outSite(newsite)
 192
 193             # get first Halfedge to the LEFT and RIGHT of the new site
 194             lbnd = edgeList.leftbnd(newsite)
 195             rbnd = lbnd.right
 196
 197             # if this halfedge has no edge, bot = bottom site (whatever that is)
 198             # create a new edge that bisects
 199             bot  = lbnd.rightreg(bottomsite)
 200             edge = Edge.bisect(bot,newsite)
 201             context.outBisector(edge)
 202
 203             # create a new Halfedge, setting its pm field to 0 and insert
 204             # this new bisector edge between the left and right vectors in
 205             # a linked list
 206             bisector = Halfedge(edge,Edge.LE)
 207             edgeList.insert(lbnd,bisector)
 208
 209             # if the new bisector intersects with the left edge, remove
 210             # the left edge's vertex, and put in the new one
 211             p = lbnd.intersect(bisector)
 212             if p is not None:
 213                 priorityQ.delete(lbnd)
 214                 priorityQ.insert(lbnd,p,newsite.distance(p))
 215
 216             # create a new Halfedge, setting its pm field to 1
 217             # insert the new Halfedge to the right of the original bisector
 218             lbnd = bisector
 219             bisector = Halfedge(edge,Edge.RE)
 220             edgeList.insert(lbnd,bisector)
 221
 222             # if this new bisector intersects with the right Halfedge
 223             p = bisector.intersect(rbnd)
 224             if p is not None:
 225                 # push the Halfedge into the ordered linked list of vertices
 226                 priorityQ.insert(bisector,p,newsite.distance(p))
 227
 228             newsite = siteIter.next()
 229
 230         elif not priorityQ.isEmpty():
 231             # intersection is smallest - this is a vector (circle) event
 232
 233             # pop the Halfedge with the lowest vector off the ordered list of
 234             # vectors.  Get the Halfedge to the left and right of the above HE
 235             # and also the Halfedge to the right of the right HE
 236             lbnd  = priorityQ.popMinHalfedge()
 237             llbnd = lbnd.left
 238             rbnd  = lbnd.right
 239             rrbnd = rbnd.right
 240
 241             # get the Site to the left of the left HE and to the right of
 242             # the right HE which it bisects
 243             bot = lbnd.leftreg(bottomsite)
 244             top = rbnd.rightreg(bottomsite)
 245
 246             # output the triple of sites, stating that a circle goes through them
 247             mid = lbnd.rightreg(bottomsite)
 248             context.outTriple(bot,top,mid)
 249
 250             # get the vertex that caused this event and set the vertex number
 251             # couldn't do this earlier since we didn't know when it would be processed
 252             v = lbnd.vertex
 253             siteList.setSiteNumber(v)
 254             context.outVertex(v)
 255
 256             # set the endpoint of the left and right Halfedge to be this vector
 257             if lbnd.edge.setEndpoint(lbnd.pm,v):
 258                 context.outEdge(lbnd.edge)
 259
 260             if rbnd.edge.setEndpoint(rbnd.pm,v):
 261                 context.outEdge(rbnd.edge)
 262
 263
 264             # delete the lowest HE, remove all vertex events to do with the
 265             # right HE and delete the right HE
 266             edgeList.delete(lbnd)
 267             priorityQ.delete(rbnd)
 268             edgeList.delete(rbnd)
 269
 270
 271             # if the site to the left of the event is higher than the Site
 272             # to the right of it, then swap them and set 'pm' to RIGHT
 273             pm = Edge.LE
 274             if bot.y > top.y:
 275                 bot,top = top,bot
 276                 pm = Edge.RE
 277
 278             # Create an Edge (or line) that is between the two Sites.  This
 279             # creates the formula of the line, and assigns a line number to it
 280             edge = Edge.bisect(bot, top)
 281             context.outBisector(edge)
 282
 283             # create a HE from the edge
 284             bisector = Halfedge(edge, pm)
 285
 286             # insert the new bisector to the right of the left HE
 287             # set one endpoint to the new edge to be the vector point 'v'
 288             # If the site to the left of this bisector is higher than the right
 289             # Site, then this endpoint is put in position 0; otherwise in pos 1
 290             edgeList.insert(llbnd, bisector)
 291             if edge.setEndpoint(Edge.RE - pm, v):
 292                 context.outEdge(edge)
 293
 294             # if left HE and the new bisector don't intersect, then delete
 295             # the left HE, and reinsert it
 296             p = llbnd.intersect(bisector)
 297             if p is not None:
 298                 priorityQ.delete(llbnd);
 299                 priorityQ.insert(llbnd, p, bot.distance(p))
 300
 301             # if right HE and the new bisector don't intersect, then reinsert it
 302             p = bisector.intersect(rrbnd)
 303             if p is not None:
 304                 priorityQ.insert(bisector, p, bot.distance(p))
 305         else:
 306             break
 307
 308     he = edgeList.leftend.right
 309     while he is not edgeList.rightend:
 310         context.outEdge(he.edge)
 311         he = he.right
 312
 313 #------------------------------------------------------------------
 314 def isEqual(a,b,relativeError=TOLERANCE):
 315     # is nearly equal to within the allowed relative error
 316     norm = max(abs(a),abs(b))
 317     return (norm < relativeError) or (abs(a - b) < (relativeError * norm))
 318
 319 #------------------------------------------------------------------
 320 class Site(object):
 321     def __init__(self,x=0.0,y=0.0,sitenum=0):
 322         self.x = x
 323         self.y = y
 324         self.sitenum = sitenum
 325
 326     def dump(self):
 327         print "Site #%d (%g, %g)" % (self.sitenum,self.x,self.y)
 328
 329     def __cmp__(self,other):
 330         if self.y < other.y:
 331             return -1
 332         elif self.y > other.y:
 333             return 1
 334         elif self.x < other.x:
 335             return -1
 336         elif self.x > other.x:
 337             return 1
 338         else:
 339             return 0
 340
 341     def distance(self,other):
 342         dx = self.x - other.x
 343         dy = self.y - other.y
 344         return math.sqrt(dx*dx + dy*dy)
 345
 346 #------------------------------------------------------------------
 347 class Edge(object):
 348     LE = 0
 349     RE = 1
 350     EDGE_NUM = 0
 351     DELETED = {}   # marker value
 352
 353     def __init__(self):
 354         self.a = 0.0
 355         self.b = 0.0
 356         self.c = 0.0
 357         self.ep  = [None,None]
 358         self.reg = [None,None]
 359         self.edgenum = 0
 360
 361     def dump(self):
 362         print "(#%d a=%g, b=%g, c=%g)" % (self.edgenum,self.a,self.b,self.c)
 363         print "ep",self.ep
 364         print "reg",self.reg
 365
 366     def setEndpoint(self, lrFlag, site):
 367         self.ep[lrFlag] = site
 368         if self.ep[Edge.RE - lrFlag] is None:
 369             return False
 370         return True
 371
 372     @staticmethod
 373     def bisect(s1,s2):
 374         newedge = Edge()
 375         newedge.reg[0] = s1 # store the sites that this edge is bisecting
 376         newedge.reg[1] = s2
 377
 378         # to begin with, there are no endpoints on the bisector - it goes to infinity
 379         # ep[0] and ep[1] are None
 380
 381         # get the difference in x dist between the sites
 382         dx = float(s2.x - s1.x)
 383         dy = float(s2.y - s1.y)
 384         adx = abs(dx)  # make sure that the difference in positive
 385         ady = abs(dy)
 386
 387         # get the slope of the line
 388         newedge.c = float(s1.x * dx + s1.y * dy + (dx*dx + dy*dy)*0.5)
 389         if adx > ady :
 390             # set formula of line, with x fixed to 1
 391             newedge.a = 1.0
 392             newedge.b = dy/dx
 393             newedge.c /= dx
 394         else:
 395             # set formula of line, with y fixed to 1
 396             newedge.b = 1.0
 397             newedge.a = dx/dy
 398             newedge.c /= dy
 399
 400         newedge.edgenum = Edge.EDGE_NUM
 401         Edge.EDGE_NUM += 1
 402         return newedge
 403
 404
 405 #------------------------------------------------------------------
 406 class Halfedge(object):
 407     def __init__(self,edge=None,pm=Edge.LE):
 408         self.left  = None   # left Halfedge in the edge list
 409         self.right = None   # right Halfedge in the edge list
 410         self.qnext = None   # priority queue linked list pointer
 411         self.edge  = edge   # edge list Edge
 412         self.pm     = pm
 413         self.vertex = None  # Site()
 414         self.ystar  = BIG_FLOAT
 415
 416     def dump(self):
 417         print "Halfedge--------------------------"
 418         print "left: ",    self.left
 419         print "right: ",   self.right
 420         print "edge: ",    self.edge
 421         print "pm: ",      self.pm
 422         print "vertex: ",
 423         if self.vertex: self.vertex.dump()
 424         else: print "None"
 425         print "ystar: ",   self.ystar
 426
 427
 428     def __cmp__(self,other):
 429         if self.ystar > other.ystar:
 430             return 1
 431         elif self.ystar < other.ystar:
 432             return -1
 433         elif self.vertex.x > other.vertex.x:
 434             return 1
 435         elif self.vertex.x < other.vertex.x:
 436             return -1
 437         else:
 438             return 0
 439
 440     def leftreg(self,default):
 441         if not self.edge:
 442             return default
 443         elif self.pm == Edge.LE:
 444             return self.edge.reg[Edge.LE]
 445         else:
 446             return self.edge.reg[Edge.RE]
 447
 448     def rightreg(self,default):
 449         if not self.edge:
 450             return default
 451         elif self.pm == Edge.LE:
 452             return self.edge.reg[Edge.RE]
 453         else:
 454             return self.edge.reg[Edge.LE]
 455
 456
 457     # returns True if p is to right of halfedge self
 458     def isPointRightOf(self,pt):
 459         e = self.edge
 460         topsite = e.reg[1]
 461         right_of_site = pt.x > topsite.x
 462
 463         if(right_of_site and self.pm == Edge.LE):
 464             return True
 465
 466         if(not right_of_site and self.pm == Edge.RE):
 467             return False
 468
 469         if(e.a == 1.0):
 470             dyp = pt.y - topsite.y
 471             dxp = pt.x - topsite.x
 472             fast = 0;
 473             if ((not right_of_site and e.b < 0.0) or (right_of_site and e.b >= 0.0)):
 474                 above = dyp >= e.b * dxp
 475                 fast = above
 476             else:
 477                 above = pt.x + pt.y * e.b > e.c
 478                 if(e.b < 0.0):
 479                     above = not above
 480                 if (not above):
 481                     fast = 1
 482             if (not fast):
 483                 dxs = topsite.x - (e.reg[0]).x
 484                 above = e.b * (dxp*dxp - dyp*dyp) < dxs*dyp*(1.0+2.0*dxp/dxs + e.b*e.b)
 485                 if(e.b < 0.0):
 486                     above = not above
 487         else:  # e.b == 1.0
 488             yl = e.c - e.a * pt.x
 489             t1 = pt.y - yl
 490             t2 = pt.x - topsite.x
 491             t3 = yl - topsite.y
 492             above = t1*t1 > t2*t2 + t3*t3
 493
 494         if(self.pm==Edge.LE):
 495             return above
 496         else:
 497             return not above
 498
 499     #--------------------------
 500     # create a new site where the Halfedges el1 and el2 intersect
 501     def intersect(self,other):
 502         e1 = self.edge
 503         e2 = other.edge
 504         if (e1 is None) or (e2 is None):
 505             return None
 506
 507         # if the two edges bisect the same parent return None
 508         if e1.reg[1] is e2.reg[1]:
 509             return None
 510
 511         d = e1.a * e2.b - e1.b * e2.a
 512         if isEqual(d,0.0):
 513             return None
 514
 515         xint = (e1.c*e2.b - e2.c*e1.b) / d
 516         yint = (e2.c*e1.a - e1.c*e2.a) / d
 517         if(cmp(e1.reg[1],e2.reg[1]) < 0):
 518             he = self
 519             e = e1
 520         else:
 521             he = other
 522             e = e2
 523
 524         rightOfSite = xint >= e.reg[1].x
 525         if((rightOfSite     and he.pm == Edge.LE) or
 526            (not rightOfSite and he.pm == Edge.RE)):
 527             return None
 528
 529         # create a new site at the point of intersection - this is a new
 530         # vector event waiting to happen
 531         return Site(xint,yint)
 532
 533
 534
 535 #------------------------------------------------------------------
 536 class EdgeList(object):
 537     def __init__(self,xmin,xmax,nsites):
 538         if xmin > xmax: xmin,xmax = xmax,xmin
 539         self.hashsize = int(2*math.sqrt(nsites+4))
 540
 541         self.xmin   = xmin
 542         self.deltax = float(xmax - xmin)
 543         self.hash   = [None]*self.hashsize
 544
 545         self.leftend  = Halfedge()
 546         self.rightend = Halfedge()
 547         self.leftend.right = self.rightend
 548         self.rightend.left = self.leftend
 549         self.hash[0]  = self.leftend
 550         self.hash[-1] = self.rightend
 551
 552     def insert(self,left,he):
 553         he.left  = left
 554         he.right = left.right
 555         left.right.left = he
 556         left.right = he
 557
 558     def delete(self,he):
 559         he.left.right = he.right
 560         he.right.left = he.left
 561         he.edge = Edge.DELETED
 562
 563     # Get entry from hash table, pruning any deleted nodes
 564     def gethash(self,b):
 565         if(b < 0 or b >= self.hashsize):
 566             return None
 567         he = self.hash[b]
 568         if he is None or he.edge is not Edge.DELETED:
 569             return he
 570
 571         #  Hash table points to deleted half edge.  Patch as necessary.
 572         self.hash[b] = None
 573         return None
 574
 575     def leftbnd(self,pt):
 576         # Use hash table to get close to desired halfedge
 577         bucket = int(((pt.x - self.xmin)/self.deltax * self.hashsize))
 578
 579         if(bucket < 0):
 580             bucket =0;
 581
 582         if(bucket >=self.hashsize):
 583             bucket = self.hashsize-1
 584
 585         he = self.gethash(bucket)
 586         if(he is None):
 587             i = 1
 588             while True:
 589                 he = self.gethash(bucket-i)
 590                 if (he is not None): break;
 591                 he = self.gethash(bucket+i)
 592                 if (he is not None): break;
 593                 i += 1
 594
 595         # Now search linear list of halfedges for the corect one
 596         if (he is self.leftend) or (he is not self.rightend and he.isPointRightOf(pt)):
 597             he = he.right
 598             while he is not self.rightend and he.isPointRightOf(pt):
 599                 he = he.right
 600             he = he.left;
 601         else:
 602             he = he.left
 603             while (he is not self.leftend and not he.isPointRightOf(pt)):
 604                 he = he.left
 605
 606         # Update hash table and reference counts
 607         if(bucket > 0 and bucket < self.hashsize-1):
 608             self.hash[bucket] = he
 609         return he
 610
 611
 612 #------------------------------------------------------------------
 613 class PriorityQueue(object):
 614     def __init__(self,ymin,ymax,nsites):
 615         self.ymin = ymin
 616         self.deltay = ymax - ymin
 617         self.hashsize = int(4 * math.sqrt(nsites))
 618         self.count = 0
 619         self.minidx = 0
 620         self.hash = []
 621         for i in range(self.hashsize):
 622             self.hash.append(Halfedge())
 623
 624     def __len__(self):
 625         return self.count
 626
 627     def isEmpty(self):
 628         return self.count == 0
 629
 630     def insert(self,he,site,offset):
 631         he.vertex = site
 632         he.ystar  = site.y + offset
 633         last = self.hash[self.getBucket(he)]
 634         next = last.qnext
 635         while((next is not None) and cmp(he,next) > 0):
 636             last = next
 637             next = last.qnext
 638         he.qnext = last.qnext
 639         last.qnext = he
 640         self.count += 1
 641
 642     def delete(self,he):
 643         if (he.vertex is not None):
 644             last = self.hash[self.getBucket(he)]
 645             while last.qnext is not he:
 646                 last = last.qnext
 647             last.qnext = he.qnext
 648             self.count -= 1
 649             he.vertex = None
 650
 651     def getBucket(self,he):
 652         bucket = int(((he.ystar - self.ymin) / self.deltay) * self.hashsize)
 653         if bucket < 0: bucket = 0
 654         if bucket >= self.hashsize: bucket = self.hashsize-1
 655         if bucket < self.minidx:  self.minidx = bucket
 656         return bucket
 657
 658     def getMinPt(self):
 659         while(self.hash[self.minidx].qnext is None):
 660             self.minidx += 1
 661         he = self.hash[self.minidx].qnext
 662         x = he.vertex.x
 663         y = he.ystar
 664         return Site(x,y)
 665
 666     def popMinHalfedge(self):
 667         curr = self.hash[self.minidx].qnext
 668         self.hash[self.minidx].qnext = curr.qnext
 669         self.count -= 1
 670         return curr
 671
 672
 673 #------------------------------------------------------------------
 674 class SiteList(object):
 675     def __init__(self,pointList):
 676         self.__sites = []
 677         self.__sitenum = 0
 678
 679         self.__xmin = pointList[0].x
 680         self.__ymin = pointList[0].y
 681         self.__xmax = pointList[0].x
 682         self.__ymax = pointList[0].y
 683         for i,pt in enumerate(pointList):
 684             self.__sites.append(Site(pt.x,pt.y,i))
 685             if pt.x < self.__xmin: self.__xmin = pt.x
 686             if pt.y < self.__ymin: self.__ymin = pt.y
 687             if pt.x > self.__xmax: self.__xmax = pt.x
 688             if pt.y > self.__ymax: self.__ymax = pt.y
 689         self.__sites.sort()
 690
 691     def setSiteNumber(self,site):
 692         site.sitenum = self.__sitenum
 693         self.__sitenum += 1
 694
 695     class Iterator(object):
 696         def __init__(this,lst):  this.generator = (s for s in lst)
 697         def __iter__(this):      return this
 698         def next(this):
 699             try:
 700                 return this.generator.next()
 701             except StopIteration:
 702                 return None
 703
 704     def iterator(self):
 705         return SiteList.Iterator(self.__sites)
 706
 707     def __iter__(self):
 708         return SiteList.Iterator(self.__sites)
 709
 710     def __len__(self):
 711         return len(self.__sites)
 712
 713     def _getxmin(self): return self.__xmin
 714     def _getymin(self): return self.__ymin
 715     def _getxmax(self): return self.__xmax
 716     def _getymax(self): return self.__ymax
 717     xmin = property(_getxmin)
 718     ymin = property(_getymin)
 719     xmax = property(_getxmax)
 720     ymax = property(_getymax)
 721
 722
 723 #------------------------------------------------------------------
 724 def computeVoronoiDiagram(points):
 725     """ Takes a list of point objects (which must have x and y fields).
 726         Returns a 3-tuple of:
 727
 728            (1) a list of 2-tuples, which are the x,y coordinates of the
 729                Voronoi diagram vertices
 730            (2) a list of 3-tuples (a,b,c) which are the equations of the
 731                lines in the Voronoi diagram: a*x + b*y = c
 732            (3) a list of 3-tuples, (l, v1, v2) representing edges of the
 733                Voronoi diagram.  l is the index of the line, v1 and v2 are
 734                the indices of the vetices at the end of the edge.  If
 735                v1 or v2 is -1, the line extends to infinity.
 736     """
 737     siteList = SiteList(points)
 738     context  = Context()
 739     voronoi(siteList,context)
 740     return (context.vertices,context.lines,context.edges)
 741
 742 #------------------------------------------------------------------
 743 def computeDelaunayTriangulation(points):
 744     """ Takes a list of point objects (which must have x and y fields).
 745         Returns a list of 3-tuples: the indices of the points that form a
 746         Delaunay triangle.
 747     """
 748     siteList = SiteList(points)
 749     context  = Context()
 750     context.triangulate = true
 751     voronoi(siteList,context)
 752     return context.triangles
 753
 754 #-----------------------------------------------------------------------------
 755 if __name__=="__main__":
 756     try:
 757         optlist,args = getopt.getopt(sys.argv[1:],"thdp")
 758     except getopt.GetoptError:
 759         usage()
 760         sys.exit(2)
 761
 762     doHelp = 0
 763     c = Context()
 764     c.doPrint = 1
 765     for opt in optlist:
 766         if opt[0] == "-d":  c.debug = 1
 767         if opt[0] == "-p":  c.plot  = 1
 768         if opt[0] == "-t":  c.triangulate = 1
 769         if opt[0] == "-h":  doHelp = 1
 770
 771     if not doHelp:
 772         pts = []
 773         fp = sys.stdin
 774         if len(args) > 0:
 775             fp = open(args[0],'r')
 776         for line in fp:
 777             fld = line.split()
 778             x = float(fld[0])
 779             y = float(fld[1])
 780             pts.append(Site(x,y))
 781         if len(args) > 0: fp.close()
 782
 783     if doHelp or len(pts) == 0:
 784         usage()
 785         sys.exit(2)
 786
 787     sl = SiteList(pts)
 788     voronoi(sl,c)
 789