text: # -*- coding: utf-8 -*-
  
  # Voronoi diagram calculator/ Delaunay triangulator
  #
  # - Voronoi Diagram Sweepline algorithm and C code by Steven Fortune,
  #   1987, http://ect.bell-labs.com/who/sjf/
  # - Python translation to file voronoi.py by Bill Simons, 2005, http://www.oxfish.com/
  # - Additional changes for QGIS by Carson Farmer added November 2010
  # - 2012 Ported to Python 3 and additional clip functions by domlysz at gmail.com
  #
  # Calculate Delaunay triangulation or the Voronoi polygons for a set of
  # 2D input points.
  #
  # Derived from code bearing the following notice:
  #
  #  The author of this software is Steven Fortune.  Copyright (c) 1994 by AT&T
  #  Bell Laboratories.
  #  Permission to use, copy, modify, and distribute this software for any
  #  purpose without fee is hereby granted, provided that this entire notice
  #  is included in all copies of any software which is or includes a copy
  #  or modification of this software and in all copies of the supporting
  #  documentation for such software.
  #  THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
  #  WARRANTY.  IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
  #  REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
  #  OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
  #
  # Comments were incorporated from Shane O'Sullivan's translation of the
  # original code into C++ (http://mapviewer.skynet.ie/voronoi.html)
  #
  # Steve Fortune's homepage: http://netlib.bell-labs.com/cm/cs/who/sjf/index.html
  #
  # For programmatic use, two functions are available:
  #
  #   computeVoronoiDiagram(points, xBuff, yBuff, polygonsOutput=False, formatOutput=False):
  #   Takes :
  #       - a list of point objects (which must have x and y fields).
  #       - x and y buffer values which are the expansion percentages of the
  #         bounding box rectangle including all input points.
  #       Returns :
  #       - With default options :
  #         A list of 2-tuples, representing the two points of each Voronoi diagram edge.
  #         Each point contains 2-tuples which are the x,y coordinates of point.
  #         if formatOutput is True, returns :
  #               - a list of 2-tuples, which are the x,y coordinates of the Voronoi diagram vertices.
  #               - and a list of 2-tuples (v1, v2) representing edges of the Voronoi diagram.
  #                 v1 and v2 are the indices of the vertices at the end of the edge.
  #       - If polygonsOutput option is True, returns :
  #         A dictionary of polygons, keys are the indices of the input points,
  #         values contains n-tuples representing the n points of each Voronoi diagram polygon.
  #         Each point contains 2-tuples which are the x,y coordinates of point.
  #         if formatOutput is True, returns :
  #               - A list of 2-tuples, which are the x,y coordinates of the Voronoi diagram vertices.
  #               - and a dictionary of input points indices. Values contains n-tuples representing
  #                 the n points of each Voronoi diagram polygon.
  #                 Each tuple contains the vertex indices of the polygon vertices.
  #
  #   computeDelaunayTriangulation(points):
  #       Takes a list of point objects (which must have x and y fields).
  #       Returns a list of 3-tuples: the indices of the points that form a Delaunay triangle.
  
  import bpy
  import math
  
  # Globals
  TOLERANCE = 1e-9
  BIG_FLOAT = 1e38
  
  
  class Context(object):
  
      def __init__(self):
          self.doPrint = 0
          self.debug = 0
  
          # tuple (xmin, xmax, ymin, ymax)
          self.extent = ()
          self.triangulate = False
          # list of vertex 2-tuples: (x,y)
          self.vertices = []
          # equation of line 3-tuple (a b c), for the equation of the line a*x+b*y = c
          self.lines = []
  
          # edge 3-tuple: (line index, vertex 1 index, vertex 2 index)
          # if either vertex index is -1, the edge extends to infinity
          self.edges = []
          # 3-tuple of vertex indices
          self.triangles = []
          # a dict of site:[edges] pairs
          self.polygons = {}
  
  
  # Clip functions #
      def getClipEdges(self):
          xmin, xmax, ymin, ymax = self.extent
          clipEdges = []
          for edge in self.edges:
              equation = self.lines[edge[0]]       # line equation
              if edge[1] != -1 and edge[2] != -1:  # finite line
                  x1, y1 = self.vertices[edge[1]][0], self.vertices[edge[1]][1]
                  x2, y2 = self.vertices[edge[2]][0], self.vertices[edge[2]][1]
                  pt1, pt2 = (x1, y1), (x2, y2)
                  inExtentP1, inExtentP2 = self.inExtent(x1, y1), self.inExtent(x2, y2)
                  if inExtentP1 and inExtentP2:
                      clipEdges.append((pt1, pt2))
                  elif inExtentP1 and not inExtentP2:
                      pt2 = self.clipLine(x1, y1, equation, leftDir=False)
                      clipEdges.append((pt1, pt2))
                  elif not inExtentP1 and inExtentP2:
                      pt1 = self.clipLine(x2, y2, equation, leftDir=True)
                      clipEdges.append((pt1, pt2))
              else:  # infinite line
                  if edge[1] != -1:
                      x1, y1 = self.vertices[edge[1]][0], self.vertices[edge[1]][1]
                      leftDir = False
                  else:
                      x1, y1 = self.vertices[edge[2]][0], self.vertices[edge[2]][1]
                      leftDir = True
                  if self.inExtent(x1, y1):
                      pt1 = (x1, y1)
                      pt2 = self.clipLine(x1, y1, equation, leftDir)
                      clipEdges.append((pt1, pt2))
          return clipEdges
  
      def getClipPolygons(self, closePoly):
          xmin, xmax, ymin, ymax = self.extent
          poly = {}
          for inPtsIdx, edges in self.polygons.items():
              clipEdges = []
              for edge in edges:
                  equation = self.lines[edge[0]]       # line equation
                  if edge[1] != -1 and edge[2] != -1:  # finite line
                      x1, y1 = self.vertices[edge[1]][0], self.vertices[edge[1]][1]
                      x2, y2 = self.vertices[edge[2]][0], self.vertices[edge[2]][1]
                      pt1, pt2 = (x1, y1), (x2, y2)
                      inExtentP1, inExtentP2 = self.inExtent(x1, y1), self.inExtent(x2, y2)
                      if inExtentP1 and inExtentP2:
                          clipEdges.append((pt1, pt2))
                      elif inExtentP1 and not inExtentP2:
                          pt2 = self.clipLine(x1, y1, equation, leftDir=False)
                          clipEdges.append((pt1, pt2))
                      elif not inExtentP1 and inExtentP2:
                          pt1 = self.clipLine(x2, y2, equation, leftDir=True)
                          clipEdges.append((pt1, pt2))
                  else:  # infinite line
                      if edge[1] != -1:
                          x1, y1 = self.vertices[edge[1]][0], self.vertices[edge[1]][1]
                          leftDir = False
                      else:
                          x1, y1 = self.vertices[edge[2]][0], self.vertices[edge[2]][1]
                          leftDir = True
                      if self.inExtent(x1, y1):
                          pt1 = (x1, y1)
                          pt2 = self.clipLine(x1, y1, equation, leftDir)
                          clipEdges.append((pt1, pt2))
              # create polygon definition from edges and check if polygon is completely closed
              polyPts, complete = self.orderPts(clipEdges)
              if not complete:
                  startPt = polyPts[0]
                  endPt = polyPts[-1]
                  # if start & end points are collinear then they are along an extent border
                  if startPt[0] == endPt[0] or startPt[1] == endPt[1]:
                      polyPts.append(polyPts[0])  # simple close
                  else:  # close at extent corner
                      # upper left
                      if (startPt[0] == xmin and endPt[1] == ymax) or (endPt[0] == xmin and startPt[1] == ymax):
                          polyPts.append((xmin, ymax))  # corner point
                          polyPts.append(polyPts[0])    # close polygon
                      # upper right
                      if (startPt[0] == xmax and endPt[1] == ymax) or (endPt[0] == xmax and startPt[1] == ymax):
                          polyPts.append((xmax, ymax))
                          polyPts.append(polyPts[0])
                      # bottom right
                      if (startPt[0] == xmax and endPt[1] == ymin) or (endPt[0] == xmax and startPt[1] == ymin):
                          polyPts.append((xmax, ymin))
                          polyPts.append(polyPts[0])
                      # bottom left
                      if (startPt[0] == xmin and endPt[1] == ymin) or (endPt[0] == xmin and startPt[1] == ymin):
                          polyPts.append((xmin, ymin))
                          polyPts.append(polyPts[0])
              if not closePoly:  # unclose polygon
                  polyPts = polyPts[:-1]
              poly[inPtsIdx] = polyPts
          return poly
  
      def clipLine(self, x1, y1, equation, leftDir):
          xmin, xmax, ymin, ymax = self.extent
          a, b, c = equation
          if b == 0:       # vertical line
              if leftDir:  # left is bottom of vertical line
                  return (x1, ymax)
              else:
                  return (x1, ymin)
          elif a == 0:     # horizontal line
              if leftDir:
                  return (xmin, y1)
              else:
                  return (xmax, y1)
          else:
              y2_at_xmin = (c - a * xmin) / b
              y2_at_xmax = (c - a * xmax) / b
              x2_at_ymin = (c - b * ymin) / a
              x2_at_ymax = (c - b * ymax) / a
              intersectPts = []
              if ymin <= y2_at_xmin <= ymax:  # valid intersect point
                  intersectPts.append((xmin, y2_at_xmin))
              if ymin <= y2_at_xmax <= ymax:
                  intersectPts.append((xmax, y2_at_xmax))
              if xmin <= x2_at_ymin <= xmax:
                  intersectPts.append((x2_at_ymin, ymin))
              if xmin <= x2_at_ymax <= xmax:
                  intersectPts.append((x2_at_ymax, ymax))
              # delete duplicate (happens if intersect point is at extent corner)
              intersectPts = set(intersectPts)
              # choose target intersect point
              if leftDir:
                  pt = min(intersectPts)  # smaller x value
              else:
                  pt = max(intersectPts)
              return pt
  
      def inExtent(self, x, y):
          xmin, xmax, ymin, ymax = self.extent
          return x >= xmin and x <= xmax and y >= ymin and y <= ymax
  
      def orderPts(self, edges):
          poly = []  # returned polygon points list [pt1, pt2, pt3, pt4 ....]
          pts = []
          # get points list
          for edge in edges:
              pts.extend([pt for pt in edge])
          # try to get start & end point
          try:
              startPt, endPt = [pt for pt in pts if pts.count(pt) < 2]  # start and end point aren't duplicate
          except:  # all points are duplicate --> polygon is complete --> append some or other edge points
              complete = True
              firstIdx = 0
              poly.append(edges[0][0])
              poly.append(edges[0][1])
          else:  # incomplete --> append the first edge points
              complete = False
              # search first edge
              for i, edge in enumerate(edges):
                  if startPt in edge:  # find
                      firstIdx = i
                      break
              poly.append(edges[firstIdx][0])
              poly.append(edges[firstIdx][1])
              if poly[0] != startPt:
                  poly.reverse()
          # append next points in list
          del edges[firstIdx]
          while edges:  # all points will be treated when edges list will be empty
              currentPt = poly[-1]  # last item
              for i, edge in enumerate(edges):
                  if currentPt == edge[0]:
                      poly.append(edge[1])
                      break
                  elif currentPt == edge[1]:
                      poly.append(edge[0])
                      break
              del edges[i]
          return poly, complete
  
      def setClipBuffer(self, xpourcent, ypourcent):
          xmin, xmax, ymin, ymax = self.extent
          width = xmax - xmin
          height = ymax - ymin
          xmin = xmin - width * xpourcent / 100
          xmax = xmax + width * xpourcent / 100
          ymin = ymin - height * ypourcent / 100
          ymax = ymax + height * ypourcent / 100
          self.extent = xmin, xmax, ymin, ymax
  
      # End clip functions #
  
      def outSite(self, s):
          if(self.debug):
              print("site (%d) at %f %f" % (s.sitenum, s.x, s.y))
          elif(self.triangulate):
              pass
          elif(self.doPrint):
              print("s %f %f" % (s.x, s.y))
  
      def outVertex(self, s):
          self.vertices.append((s.x, s.y))
          if(self.debug):
              print("vertex(%d) at %f %f" % (s.sitenum, s.x, s.y))
          elif(self.triangulate):
              pass
          elif(self.doPrint):
              print("v %f %f" % (s.x, s.y))
  
      def outTriple(self, s1, s2, s3):
          self.triangles.append((s1.sitenum, s2.sitenum, s3.sitenum))
          if (self.debug):
              print("circle through left=%d right=%d bottom=%d" % (s1.sitenum, s2.sitenum, s3.sitenum))
          elif (self.triangulate and self.doPrint):
              print("%d %d %d" % (s1.sitenum, s2.sitenum, s3.sitenum))
  
      def outBisector(self, edge):
          self.lines.append((edge.a, edge.b, edge.c))
          if (self.debug):
              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)
                  )
          elif(self.doPrint):
              print("l %f %f %f" % (edge.a, edge.b, edge.c))
  
      def outEdge(self, edge):
          sitenumL = -1
          if edge.ep[Edge.LE] is not None:
              sitenumL = edge.ep[Edge.LE].sitenum
          sitenumR = -1
          if edge.ep[Edge.RE] is not None:
              sitenumR = edge.ep[Edge.RE].sitenum
  
          # polygons dict add by CF
          if edge.reg[0].sitenum not in self.polygons:
              self.polygons[edge.reg[0].sitenum] = []
          if edge.reg[1].sitenum not in self.polygons:
              self.polygons[edge.reg[1].sitenum] = []
          self.polygons[edge.reg[0].sitenum].append((edge.edgenum, sitenumL, sitenumR))
          self.polygons[edge.reg[1].sitenum].append((edge.edgenum, sitenumL, sitenumR))
  
          self.edges.append((edge.edgenum, sitenumL, sitenumR))
  
          if (not self.triangulate):
              if (self.doPrint):
                  print("e %d" % edge.edgenum)
                  print(" %d " % sitenumL)
                  print("%d" % sitenumR)
  
  
  def voronoi(siteList, context):
      context.extent = siteList.extent
      edgeList = EdgeList(siteList.xmin, siteList.xmax, len(siteList))
      priorityQ = PriorityQueue(siteList.ymin, siteList.ymax, len(siteList))
      siteIter = siteList.iterator()
  
      bottomsite = siteIter.next()
      context.outSite(bottomsite)
      newsite = siteIter.next()
      minpt = Site(-BIG_FLOAT, -BIG_FLOAT)
      while True:
          if not priorityQ.isEmpty():
              minpt = priorityQ.getMinPt()
  
          if (newsite and (priorityQ.isEmpty() or newsite < minpt)):
              # newsite is smallest -  this is a site event
              context.outSite(newsite)
  
              # get first Halfedge to the LEFT and RIGHT of the new site
              lbnd = edgeList.leftbnd(newsite)
              rbnd = lbnd.right
  
              # if this halfedge has no edge, bot = bottom site (whatever that is)
              # create a new edge that bisects
              bot = lbnd.rightreg(bottomsite)
              edge = Edge.bisect(bot, newsite)
              context.outBisector(edge)
  
              # create a new Halfedge, setting its pm field to 0 and insert
              # this new bisector edge between the left and right vectors in
              # a linked list
              bisector = Halfedge(edge, Edge.LE)
              edgeList.insert(lbnd, bisector)
  
              # if the new bisector intersects with the left edge, remove
              # the left edge's vertex, and put in the new one
              p = lbnd.intersect(bisector)
              if p is not None:
                  priorityQ.delete(lbnd)
                  priorityQ.insert(lbnd, p, newsite.distance(p))
  
              # create a new Halfedge, setting its pm field to 1
              # insert the new Halfedge to the right of the original bisector
              lbnd = bisector
              bisector = Halfedge(edge, Edge.RE)
              edgeList.insert(lbnd, bisector)
  
              # if this new bisector intersects with the right Halfedge
              p = bisector.intersect(rbnd)
              if p is not None:
                  # push the Halfedge into the ordered linked list of vertices
                  priorityQ.insert(bisector, p, newsite.distance(p))
  
              newsite = siteIter.next()
  
          elif not priorityQ.isEmpty():
              # intersection is smallest - this is a vector (circle) event
              # pop the Halfedge with the lowest vector off the ordered list of
              # vectors.  Get the Halfedge to the left and right of the above HE
              # and also the Halfedge to the right of the right HE
              lbnd = priorityQ.popMinHalfedge()
              llbnd = lbnd.left
              rbnd = lbnd.right
              rrbnd = rbnd.right
  
              # get the Site to the left of the left HE and to the right of
              # the right HE which it bisects
              bot = lbnd.leftreg(bottomsite)
              top = rbnd.rightreg(bottomsite)
  
              # output the triple of sites, stating that a circle goes through them
              mid = lbnd.rightreg(bottomsite)
              context.outTriple(bot, top, mid)
  
              # get the vertex that caused this event and set the vertex number
              # couldn't do this earlier since we didn't know when it would be processed
              v = lbnd.vertex
              siteList.setSiteNumber(v)
              context.outVertex(v)
  
              # set the endpoint of the left and right Halfedge to be this vector
              if lbnd.edge.setEndpoint(lbnd.pm, v):
                  context.outEdge(lbnd.edge)
  
              if rbnd.edge.setEndpoint(rbnd.pm, v):
                  context.outEdge(rbnd.edge)
  
              # delete the lowest HE, remove all vertex events to do with the
              # right HE and delete the right HE
              edgeList.delete(lbnd)
              priorityQ.delete(rbnd)
              edgeList.delete(rbnd)
  
              # if the site to the left of the event is higher than the Site
              # to the right of it, then swap them and set 'pm' to RIGHT
              pm = Edge.LE
              if bot.y > top.y:
                  bot, top = top, bot
                  pm = Edge.RE
  
              # Create an Edge (or line) that is between the two Sites.  This
              # creates the formula of the line, and assigns a line number to it
              edge = Edge.bisect(bot, top)
              context.outBisector(edge)
  
              # create a HE from the edge
              bisector = Halfedge(edge, pm)
  
              # insert the new bisector to the right of the left HE
              # set one endpoint to the new edge to be the vector point 'v'
              # If the site to the left of this bisector is higher than the right
              # Site, then this endpoint is put in position 0; otherwise in pos 1
              edgeList.insert(llbnd, bisector)
              if edge.setEndpoint(Edge.RE - pm, v):
                  context.outEdge(edge)
  
              # if left HE and the new bisector don't intersect, then delete
              # the left HE, and reinsert it
              p = llbnd.intersect(bisector)
              if p is not None:
                  priorityQ.delete(llbnd)
                  priorityQ.insert(llbnd, p, bot.distance(p))
  
              # if right HE and the new bisector don't intersect, then reinsert it
              p = bisector.intersect(rrbnd)
              if p is not None:
                  priorityQ.insert(bisector, p, bot.distance(p))
          else:
              break
  
      he = edgeList.leftend.right
      while he is not edgeList.rightend:
          context.outEdge(he.edge)
          he = he.right
      Edge.EDGE_NUM = 0  # CF
  
  
  def isEqual(a, b, relativeError=TOLERANCE):
      # is nearly equal to within the allowed relative error
      norm = max(abs(a), abs(b))
      return (norm < relativeError) or (abs(a - b) < (relativeError * norm))
  
  
  class Site(object):
  
      def __init__(self, x=0.0, y=0.0, sitenum=0):
          self.x = x
          self.y = y
          self.sitenum = sitenum
  
      def dump(self):
          print("Site #%d (%g, %g)" % (self.sitenum, self.x, self.y))
  
      def __lt__(self, other):
          if self.y < other.y:
              return True
          elif self.y > other.y:
              return False
          elif self.x < other.x:
              return True
          elif self.x > other.x:
              return False
          else:
              return False
  
      def __eq__(self, other):
          if self.y == other.y and self.x == other.x:
              return True
  
      def distance(self, other):
          dx = self.x - other.x
          dy = self.y - other.y
          return math.sqrt(dx * dx + dy * dy)
  
  
  class Edge(object):
      LE = 0  # left end indice --> edge.ep[Edge.LE]
      RE = 1  # right end indice
      EDGE_NUM = 0
      DELETED = {}  # marker value
  
      def __init__(self):
          self.a = 0.0  # equation of the line a*x+b*y = c
          self.b = 0.0
          self.c = 0.0
          self.ep = [None, None]  # end point (2 tuples of site)
          self.reg = [None, None]
          self.edgenum = 0
  
      def dump(self):
          print("(#%d a=%g, b=%g, c=%g)" % (self.edgenum, self.a, self.b, self.c))
          print("ep", self.ep)
          print("reg", self.reg)
  
      def setEndpoint(self, lrFlag, site):
          self.ep[lrFlag] = site
          if self.ep[Edge.RE - lrFlag] is None:
              return False
          return True
  
      @staticmethod
      def bisect(s1, s2):
          newedge = Edge()
          newedge.reg[0] = s1  # store the sites that this edge is bisecting
          newedge.reg[1] = s2
  
          # to begin with, there are no endpoints on the bisector - it goes to infinity
          # ep[0] and ep[1] are None
  
          # get the difference in x dist between the sites
          dx = float(s2.x - s1.x)
          dy = float(s2.y - s1.y)
          adx = abs(dx)  # make sure that the difference in positive
          ady = abs(dy)
  
          # get the slope of the line
          newedge.c = float(s1.x * dx + s1.y * dy + (dx * dx + dy * dy) * 0.5)
          if adx > ady:
              # set formula of line, with x fixed to 1
              newedge.a = 1.0
              newedge.b = dy / dx
              newedge.c /= dx
          else:
              # set formula of line, with y fixed to 1
              newedge.b = 1.0
              newedge.a = dx / dy
              newedge.c /= dy
  
          newedge.edgenum = Edge.EDGE_NUM
          Edge.EDGE_NUM += 1
          return newedge
  
  
  class Halfedge(object):
  
      def __init__(self, edge=None, pm=Edge.LE):
          self.left = None    # left Halfedge in the edge list
          self.right = None   # right Halfedge in the edge list
          self.qnext = None   # priority queue linked list pointer
          self.edge = edge    # edge list Edge
          self.pm = pm
          self.vertex = None  # Site()
          self.ystar = BIG_FLOAT
  
      def dump(self):
          print("Halfedge--------------------------")
          print("left: ", self.left)
          print("right: ", self.right)
          print("edge: ", self.edge)
          print("pm: ", self.pm)
          print("vertex: "),
          if self.vertex:
              self.vertex.dump()
          else:
              print("None")
          print("ystar: ", self.ystar)
  
      def __lt__(self, other):
          if self.ystar < other.ystar:
              return True
          elif self.ystar > other.ystar:
              return False
          elif self.vertex.x < other.vertex.x:
              return True
          elif self.vertex.x > other.vertex.x:
              return False
          else:
              return False
  
      def __eq__(self, other):
          if self.ystar == other.ystar and self.vertex.x == other.vertex.x:
              return True
  
      def leftreg(self, default):
          if not self.edge:
              return default
          elif self.pm == Edge.LE:
              return self.edge.reg[Edge.LE]
          else:
              return self.edge.reg[Edge.RE]
  
      def rightreg(self, default):
          if not self.edge:
              return default
          elif self.pm == Edge.LE:
              return self.edge.reg[Edge.RE]
          else:
              return self.edge.reg[Edge.LE]
  
      # returns True if p is to right of halfedge self
      def isPointRightOf(self, pt):
          e = self.edge
          topsite = e.reg[1]
          right_of_site = pt.x > topsite.x
  
          if(right_of_site and self.pm == Edge.LE):
              return True
  
          if(not right_of_site and self.pm == Edge.RE):
              return False
  
          if(e.a == 1.0):
              dyp = pt.y - topsite.y
              dxp = pt.x - topsite.x
              fast = 0
              if ((not right_of_site and e.b < 0.0) or (right_of_site and e.b >= 0.0)):
                  above = dyp >= e.b * dxp
                  fast = above
              else:
                  above = pt.x + pt.y * e.b > e.c
                  if(e.b < 0.0):
                      above = not above
                  if (not above):
                      fast = 1
              if (not fast):
                  dxs = topsite.x - (e.reg[0]).x
                  above = e.b * (dxp * dxp - dyp * dyp) < dxs * dyp * (1.0 + 2.0 * dxp / dxs + e.b * e.b)
                  if(e.b < 0.0):
                      above = not above
          else:  # e.b == 1.0
              yl = e.c - e.a * pt.x
              t1 = pt.y - yl
              t2 = pt.x - topsite.x
              t3 = yl - topsite.y
              above = t1 * t1 > t2 * t2 + t3 * t3
  
          if(self.pm == Edge.LE):
              return above
          else:
              return not above
  
      # create a new site where the Halfedges el1 and el2 intersect
      def intersect(self, other):
          e1 = self.edge
          e2 = other.edge
          if (e1 is None) or (e2 is None):
              return None
  
          # if the two edges bisect the same parent return None
          if e1.reg[1] is e2.reg[1]:
              return None
  
          d = e1.a * e2.b - e1.b * e2.a
          if isEqual(d, 0.0):
              return None
  
          xint = (e1.c * e2.b - e2.c * e1.b) / d
          yint = (e2.c * e1.a - e1.c * e2.a) / d
          if e1.reg[1] < e2.reg[1]:
              he = self
              e = e1
          else:
              he = other
              e = e2
  
          rightOfSite = xint >= e.reg[1].x
          if((rightOfSite and he.pm == Edge.LE) or
                  (not rightOfSite and he.pm == Edge.RE)):
              return None
  
          # create a new site at the point of intersection - this is a new
          # vector event waiting to happen
          return Site(xint, yint)
  
  
  class EdgeList(object):
  
      def __init__(self, xmin, xmax, nsites):
          if xmin > xmax:
              xmin, xmax = xmax, xmin
          self.hashsize = int(2 * math.sqrt(nsites + 4))
  
          self.xmin = xmin
          self.deltax = float(xmax - xmin)
          self.hash = [None] * self.hashsize
  
          self.leftend = Halfedge()
          self.rightend = Halfedge()
          self.leftend.right = self.rightend
          self.rightend.left = self.leftend
          self.hash[0] = self.leftend
          self.hash[-1] = self.rightend
  
      def insert(self, left, he):
          he.left = left
          he.right = left.right
          left.right.left = he
          left.right = he
  
      def delete(self, he):
          he.left.right = he.right
          he.right.left = he.left
          he.edge = Edge.DELETED
  
      # Get entry from hash table, pruning any deleted nodes
      def gethash(self, b):
          if(b < 0 or b >= self.hashsize):
              return None
          he = self.hash[b]
          if he is None or he.edge is not Edge.DELETED:
              return he
  
          #  Hash table points to deleted half edge.  Patch as necessary.
          self.hash[b] = None
          return None
  
      def leftbnd(self, pt):
          # Use hash table to get close to desired halfedge
          bucket = int(((pt.x - self.xmin) / self.deltax * self.hashsize))
  
          if(bucket < 0):
              bucket = 0
  
          if(bucket >= self.hashsize):
              bucket = self.hashsize - 1
  
          he = self.gethash(bucket)
          if(he is None):
              i = 1
              while True:
                  he = self.gethash(bucket - i)
                  if (he is not None):
                      break
                  he = self.gethash(bucket + i)
                  if (he is not None):
                      break
                  i += 1
  
          # Now search linear list of halfedges for the correct one
          if (he is self.leftend) or (he is not self.rightend and he.isPointRightOf(pt)):
              he = he.right
              while he is not self.rightend and he.isPointRightOf(pt):
                  he = he.right
              he = he.left
          else:
              he = he.left
              while (he is not self.leftend and not he.isPointRightOf(pt)):
                  he = he.left
  
          # Update hash table and reference counts
          if(bucket > 0 and bucket < self.hashsize - 1):
              self.hash[bucket] = he
          return he
  
  
  class PriorityQueue(object):
  
      def __init__(self, ymin, ymax, nsites):
          self.ymin = ymin
          self.deltay = ymax - ymin
          self.hashsize = int(4 * math.sqrt(nsites))
          self.count = 0
          self.minidx = 0
          self.hash = []
          for i in range(self.hashsize):
              self.hash.append(Halfedge())
  
      def __len__(self):
          return self.count
  
      def isEmpty(self):
          return self.count == 0
  
      def insert(self, he, site, offset):
          he.vertex = site
          he.ystar = site.y + offset
          last = self.hash[self.getBucket(he)]
          next = last.qnext
          while((next is not None) and he > next):
              last = next
              next = last.qnext
          he.qnext = last.qnext
          last.qnext = he
          self.count += 1
  
      def delete(self, he):
          if (he.vertex is not None):
              last = self.hash[self.getBucket(he)]
              while last.qnext is not he:
                  last = last.qnext
              last.qnext = he.qnext
              self.count -= 1
              he.vertex = None
  
      def getBucket(self, he):
          bucket = int(((he.ystar - self.ymin) / self.deltay) * self.hashsize)
          if bucket < 0:
              bucket = 0
          if bucket >= self.hashsize:
              bucket = self.hashsize - 1
          if bucket < self.minidx:
              self.minidx = bucket
          return bucket
  
      def getMinPt(self):
          while(self.hash[self.minidx].qnext is None):
              self.minidx += 1
          he = self.hash[self.minidx].qnext
          x = he.vertex.x
          y = he.ystar
          return Site(x, y)
  
      def popMinHalfedge(self):
          curr = self.hash[self.minidx].qnext
          self.hash[self.minidx].qnext = curr.qnext
          self.count -= 1
          return curr
  
  
  class SiteList(object):
  
      def __init__(self, pointList):
          self.__sites = []
          self.__sitenum = 0
  
          self.__xmin = min([pt.x for pt in pointList])
          self.__ymin = min([pt.y for pt in pointList])
          self.__xmax = max([pt.x for pt in pointList])
          self.__ymax = max([pt.y for pt in pointList])
          self.__extent = (self.__xmin, self.__xmax, self.__ymin, self.__ymax)
  
          for i, pt in enumerate(pointList):
              self.__sites.append(Site(pt.x, pt.y, i))
          self.__sites.sort()
  
      def setSiteNumber(self, site):
          site.sitenum = self.__sitenum
          self.__sitenum += 1
  
      class Iterator(object):
  
          def __init__(this, lst):
              this.generator = (s for s in lst)
  
          def __iter__(this):
              return this
  
          def next(this):
              try:
                  # Note: Blender is Python 3.x so no need for 2.x checks
                  return this.generator.__next__()
              except StopIteration:
                  return None
  
      def iterator(self):
          return SiteList.Iterator(self.__sites)
  
      def __iter__(self):
          return SiteList.Iterator(self.__sites)
  
      def __len__(self):
          return len(self.__sites)
  
      def _getxmin(self):
          return self.__xmin
  
      def _getymin(self):
          return self.__ymin
  
      def _getxmax(self):
          return self.__xmax
  
      def _getymax(self):
          return self.__ymax
  
      def _getextent(self):
          return self.__extent
  
      xmin = property(_getxmin)
      ymin = property(_getymin)
      xmax = property(_getxmax)
      ymax = property(_getymax)
      extent = property(_getextent)
  
  
  def computeVoronoiDiagram(points, xBuff=0, yBuff=0, polygonsOutput=False,
                            formatOutput=False, closePoly=True):
      """
      Takes :
      - a list of point objects (which must have x and y fields).
      - x and y buffer values which are the expansion percentages of the bounding box
          rectangle including all input points.
      Returns :
      - With default options :
        A list of 2-tuples, representing the two points of each Voronoi diagram edge.
        Each point contains 2-tuples which are the x,y coordinates of point.
        if formatOutput is True, returns :
                      - a list of 2-tuples, which are the x,y coordinates of the Voronoi diagram vertices.
                      - and a list of 2-tuples (v1, v2) representing edges of the Voronoi diagram.
                        v1 and v2 are the indices of the vertices at the end of the edge.
      - If polygonsOutput option is True, returns :
        A dictionary of polygons, keys are the indices of the input points,
        values contains n-tuples representing the n points of each Voronoi diagram polygon.
        Each point contains 2-tuples which are the x,y coordinates of point.
        if formatOutput is True, returns :
                      - A list of 2-tuples, which are the x,y coordinates of the Voronoi diagram vertices.
                      - and a dictionary of input points indices. Values contains n-tuples representing
                        the n points of each Voronoi diagram polygon.
                        Each tuple contains the vertex indices of the polygon vertices.
      - if closePoly is True then, in the list of points of a polygon, last point will be the same of first point
      """
      siteList = SiteList(points)
      context = Context()
      voronoi(siteList, context)
      context.setClipBuffer(xBuff, yBuff)
      if not polygonsOutput:
          clipEdges = context.getClipEdges()
          if formatOutput:
              vertices, edgesIdx = formatEdgesOutput(clipEdges)
              return vertices, edgesIdx
          else:
              return clipEdges
      else:
          clipPolygons = context.getClipPolygons(closePoly)
          if formatOutput:
              vertices, polyIdx = formatPolygonsOutput(clipPolygons)
              return vertices, polyIdx
          else:
              return clipPolygons
  
  
  def formatEdgesOutput(edges):
      # get list of points
      pts = []
      for edge in edges:
          pts.extend(edge)
      # get unique values
      pts = set(pts)  # unique values (tuples are hashable)
      # get dict {values:index}
      valuesIdxDict = dict(zip(pts, range(len(pts))))
      # get edges index reference
      edgesIdx = []
      for edge in edges:
          edgesIdx.append([valuesIdxDict[pt] for pt in edge])
      return list(pts), edgesIdx
  
  
  def formatPolygonsOutput(polygons):
      # get list of points
      pts = []
      for poly in polygons.values():
          pts.extend(poly)
      # get unique values
      pts = set(pts)  # unique values (tuples are hashable)
      # get dict {values:index}
      valuesIdxDict = dict(zip(pts, range(len(pts))))
      # get polygons index reference
      polygonsIdx = {}
      for inPtsIdx, poly in polygons.items():
          polygonsIdx[inPtsIdx] = [valuesIdxDict[pt] for pt in poly]
      return list(pts), polygonsIdx
  
  
  def computeDelaunayTriangulation(points):
      """ Takes a list of point objects (which must have x and y fields).
              Returns a list of 3-tuples: the indices of the points that form a
              Delaunay triangle.
      """
      siteList = SiteList(points)
      context = Context()
      context.triangulate = True
      voronoi(siteList, context)
      return context.triangles
  