Delaunay Triangulation 互动算法

import math EPSILON = math.pow(2,-52) EDGE_STACK =[None] * 512 class Delaunator: def __init__(self,points): n = len(points) if (len(points) < 3): raise ValueError("Need at least 3 points") coords = [None] * n * 2 for i in range(0,n): p = points[i] coords[2 * i] = (p[0]) coords[2 * i+1] = (p[1]) triangles = self.constructor(coords) def constructor(self, coords): n = len(coords) >> 1 self.coords = coords # arrays that will store the triangulation graph maxTriangles = max(2 * n - 5, 0) self._triangles = [None] * maxTriangles * 3 self._halfedges = [None] * maxTriangles * 3 # temporary arrays for tracking the edges of the advancing convex hull self.hashSize = math.ceil(math.sqrt(n)) self.hullPrev = [None] * n # edge to prev edge self.hullNext = [None] * n # edge to next edge self.hullTri = [None] * n # edge to adjacent triangle self.hullHash = [-1] * self.hashSize # angular edge hash # temporary arrays for sorting points self._ids = [None] * n self._dists = [None] * n triangles = self.update(coords) return triangles def update(self,coords): n = len(coords) >> 1 # populate an array of point indices; calculate input data bbox minX = math.inf minY = math.inf maxX = -math.inf maxY = -math.inf for i in range(0,n): x = coords[2 * i] y = coords[2 * i + 1] if (x < minX): minX = x if (y < minY): minY = y if (x > maxX): maxX = x if (y > maxY): maxY = y self._ids[i] = i cx = (minX + maxX) / 2 cy = (minY + maxY) / 2 minDist = math.inf i0 = 0 i1 = 0 i2 = 0 # pick a seed point close to the center for i in range(0,n): d = dist(cx, cy, coords[2 * i], coords[2 * i + 1]) if (d < minDist): i0 = i minDist = d i0x = coords[2 * i0] i0y = coords[2 * i0 + 1] minDist = math.inf # find the point closest to the seed for i in range(0,n): if (i == i0): continue d = dist(i0x, i0y, coords[2 * i], coords[2 * i + 1]) if (d < minDist and d > 0): i1 = i minDist = d i1x = coords[2 * i1] i1y = coords[2 * i1 + 1] minRadius = math.inf # find the third point which forms the smallest circumcircle with the first two for i in range(0,n): if (i == i0 or i == i1): continue r = circumradius(i0x, i0y, i1x, i1y, coords[2 * i], coords[2 * i + 1]) if (r < minRadius): i2 = i minRadius = r i2x = coords[2 * i2] i2y = coords[2 * i2 + 1] if (minRadius == math.inf): # order collinear points by dx (or dy if all x are identical) # and return the list as a hull for i in range(0,n): self._dists[i] = (coords[2 * i] - coords[0]) or (coords[2 * i + 1] - coords[1]) quicksort(self._ids, self._dists, 0, n - 1) hull = [None] * n j = 0 d0 = -math.inf for i in range(0,n): id = self._ids[i] if (self._dists[id] > d0): hull[j] = id j+=1 d0 = self._dists[id] self.hull = hull[0:j] self.triangles = [] self.halfedges = [] # swap the order of the seed points for counter-clockwise orientation if (orient(i0x, i0y, i1x, i1y, i2x, i2y)): i = i1 x = i1x y = i1y i1 = i2 i1x = i2x i1y = i2y i2 = i i2x = x i2y = y center = circumcenter(i0x, i0y, i1x, i1y, i2x, i2y) self._cx = center[0] self._cy = center[1] for i in range(0,n): self._dists[i] = dist(coords[2 * i], coords[2 * i + 1], center[0], center[1]) # sort the points by distance from the seed triangle circumcenter quicksort(self._ids, self._dists, 0, n - 1) # set up the seed triangle as the starting hull self._hullStart = i0 hullSize = 3 self.hullNext[i0] = self.hullPrev[i2] = i1 self.hullNext[i1] = self.hullPrev[i0] = i2 self.hullNext[i2] = self.hullPrev[i1] = i0 self.hullTri[i0] = 0 self.hullTri[i1] = 1 self.hullTri[i2] = 2 self.hullHash[self._hashKey(i0x, i0y)] = i0 self.hullHash[self._hashKey(i1x, i1y)] = i1 self.hullHash[self._hashKey(i2x, i2y)] = i2 self.trianglesLen = 0 self._addTriangle(i0, i1, i2, -1, -1, -1) xp=0 yp=0 for k in range(0,len(self._ids)): i = self._ids[k] x = coords[2 * i] y = coords[2 * i + 1] # skip near-duplicate points if (k > 0 and abs(x - xp) <= EPSILON and abs(y - yp) <= EPSILON): continue xp = x yp = y # skip seed triangle points if (i == i0 or i == i1 or i == i2): continue # find a visible edge on the convex hull using edge hash start = 0 key = self._hashKey(x, y) for j in range(0,self.hashSize): start = self.hullHash[(key + j) % self.hashSize] if (start != -1 and start != self.hullNext[start]): break start = self.hullPrev[start] e = start while True: q = self.hullNext[e] if orient(x, y, coords[2 * e], coords[2 * e + 1], coords[2 * q], coords[2 * q + 1]): break e = q if (e == start): e = -1 break if (e == -1): continue # likely a near-duplicate point; skip it # add the first triangle from the point t = self._addTriangle(e, i, self.hullNext[e], -1, -1, self.hullTri[e]) # recursively flip triangles from the point until they satisfy the Delaunay condition self.hullTri[i] = self._legalize(t + 2,coords) self.hullTri[e] = t # keep track of boundary triangles on the hull hullSize+=1 # walk forward through the hull, adding more triangles and flipping recursively n = self.hullNext[e] while True: q = self.hullNext[n] if not (orient(x, y, coords[2 * n], coords[2 * n + 1], coords[2 * q], coords[2 * q + 1])): break t = self._addTriangle(n, i, q, self.hullTri[i], -1, self.hullTri[n]) self.hullTri[i] = self._legalize(t + 2,coords) self.hullNext[n] = n # mark as removed hullSize-=1 n = q # walk backward from the other side, adding more triangles and flipping if (e == start): while True: q = self.hullPrev[e] if not (orient(x, y, coords[2 * q], coords[2 * q + 1], coords[2 * e], coords[2 * e + 1])): break t = self._addTriangle(q, i, e, -1, self.hullTri[e], self.hullTri[q]) self._legalize(t + 2,coords) self.hullTri[q] = t self.hullNext[e] = e # mark as removed hullSize-=1 e = q # update the hull indices