home *** CD-ROM | disk | FTP | other *** search
/ Chip 2011 November / CHIP_2011_11.iso / Programy / Narzedzia / Inkscape / Inkscape-0.48.2-1-win32.exe / share / extensions / voronoi.py < prev    next >
Text File  |  2011-07-08  |  27KB  |  790 lines
  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.  
  790.