kopia lustrzana https://github.com/vilemduha/blendercam
995 wiersze
36 KiB
Python
995 wiersze
36 KiB
Python
# -*- 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 math
|
|
import sys
|
|
import getopt
|
|
|
|
TOLERANCE = 1e-9
|
|
BIG_FLOAT = 1e38
|
|
|
|
if sys.version > '3':
|
|
PY3 = True
|
|
else:
|
|
PY3 = False
|
|
|
|
|
|
# ------------------------------------------------------------------
|
|
class Context(object):
|
|
def __init__(self):
|
|
self.doPrint = 0
|
|
self.debug = 0
|
|
self.extent = () # tuple (xmin, xmax, ymin, ymax)
|
|
self.triangulate = False
|
|
self.vertices = [] # list of vertex 2-tuples: (x,y)
|
|
self.lines = [] # equation of line 3-tuple (a b c), for the equation of the line a*x+b*y = c
|
|
self.edges = [] # edge 3-tuple: (line index, vertex 1 index, vertex 2 index) if either vertex index is -1, the edge extends to infinity
|
|
self.triangles = [] # 3-tuple of vertex indices
|
|
self.polygons = {} # a dict of site:[edges] pairs
|
|
|
|
########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 startPt[0] == endPt[0] or startPt[1] == endPt[
|
|
1]: # if start & end points are collinear then they are along an extent border
|
|
polyPts.append(polyPts[0]) # simple close
|
|
else: # close at extent corner
|
|
if (startPt[0] == xmin and endPt[1] == ymax) or (
|
|
endPt[0] == xmin and startPt[1] == ymax): # upper left
|
|
polyPts.append((xmin, ymax)) # corner point
|
|
polyPts.append(polyPts[0]) # close polygon
|
|
if (startPt[0] == xmax and endPt[1] == ymax) or (
|
|
endPt[0] == xmax and startPt[1] == ymax): # upper right
|
|
polyPts.append((xmax, ymax))
|
|
polyPts.append(polyPts[0])
|
|
if (startPt[0] == xmax and endPt[1] == ymin) or (
|
|
endPt[0] == xmax and startPt[1] == ymin): # bottom right
|
|
polyPts.append((xmax, ymin))
|
|
polyPts.append(polyPts[0])
|
|
if (startPt[0] == xmin and endPt[1] == ymin) or (
|
|
endPt[0] == xmin and startPt[1] == ymin): # bottom left
|
|
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:
|
|
if PY3:
|
|
return this.generator.__next__()
|
|
else:
|
|
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
|
|
|
|
# -----------------------------------------------------------------------------
|
|
# if __name__=="__main__":
|