diff options
Diffstat (limited to 'add_advanced_objects_panels/DelaunayVoronoi.py')
| -rw-r--r-- | add_advanced_objects_panels/DelaunayVoronoi.py | 998 |
1 files changed, 998 insertions, 0 deletions
diff --git a/add_advanced_objects_panels/DelaunayVoronoi.py b/add_advanced_objects_panels/DelaunayVoronoi.py new file mode 100644 index 00000000..dcce7f68 --- /dev/null +++ b/add_advanced_objects_panels/DelaunayVoronoi.py @@ -0,0 +1,998 @@ +# -*- 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 + witdh = xmax - xmin + height = ymax - ymin + xmin = xmin - witdh * xpourcent / 100 + xmax = xmax + witdh * 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 corect 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 |