input : vertex list output : triangle list initialize the triangle list determine the supertriangle add supertriangle vertices to the end of the vertex list add the supertriangle to the triangle list for each sample point in the vertex list initialize the edge buffer for each triangle currently in the triangle list calculate the triangle circumcircle center and radius if the point lies in the triangle circumcircle then add the three triangle edges to the edge buffer remove the triangle from the triangle list endif endfor delete all doubly specified edges from the edge buffer this leaves the edges of the enclosing polygon only add to the triangle list all triangles formed between the point and the edges of the enclosing polygon endfor remove any triangles from the triangle list that use the supertriangle vertices remove the supertriangle vertices from the vertex list end
def__init__(self, points: List[Point]): self.points: List[Point] = points self.super_triangle = self.getSuperTriangle() self.triangles: List[Triangle] = [self.super_triangle] # add super-triangle to triangulation for point in points: # add all the initial points one at a time to the triangulation self.addPoint(point) self.removeSuperTriangle()
defaddPoint(self, point: Point) -> None: bad_triangles: List[Triangle] = [] # first find all the triangles that are no longer valid due to the insertion for triangle in self.triangles: if point.isInCircumcircleOf(triangle): bad_triangles.append(triangle)
polygon: List[Edge] = [] # Find the boundary of the polygonal hole for triangle1 in bad_triangles: for edge in triangle1.edges: edge_shared = False for triangle2 in bad_triangles: if triangle1 == triangle2: continue if edge in triangle2.edges: edge_shared = True ifnot edge_shared: polygon.append(edge)
# Remove broken triangles from the triangulation list for triangle in bad_triangles: self.triangles.remove(triangle)
# Create triangles with the newly created edges for edge in polygon: new_triangle = Triangle(edge.begin, edge.end, point) self.triangles.append(new_triangle)
defremoveSuperTriangle(self) -> None: # Remove the triangles that has connection to the super-triangle super_verticies = self.super_triangle.verticies for i inrange(len(self.triangles) - 1, -1, -1): triangle = self.triangles[i] has_common = False for vertex1 in triangle.verticies: for vertex2 in super_verticies: if vertex1 == vertex2: has_common = True if has_common: self.triangles.remove(triangle)
defexportTriangles(self): ps = [p for t in self.triangles for p in t.verticies] xs = [p.x for p in ps] ys = [p.y for p in ps] ts = [(ps.index(t.verticies[0]), ps.index(t.verticies[1]), ps.index(t.verticies[2])) for t in self.triangles] return xs, ys, ts
from random import randint, seed from Delaunay_Triangulation import Delaunay_Triangulation, Point import time
if __name__ == '__main__': seed(5) n = 10 xs = [randint(1, 98) for x inrange(n)] ys = [randint(1, 98) for x inrange(n)] seted_points = set(zip(xs, ys)) print("The actual number of input points: ", len(seted_points)) points = [Point(x, y) for x, y in seted_points] start_time = time.time() dt = Delaunay_Triangulation(points) print(f"Triangulating {len(seted_points)} points takes {time.time() - start_time} s") # number of DT triangles print(len(dt.triangles), "Delaunay triangles")
import matplotlib.pyplot as plt import matplotlib.tri as tri
classVertex: def__init__(self, x, y, _id = None): self.id = _id self.x = x self.y = y
self.name = f'v_{self.id}'# for debugging
classQuad_Edge: """A directed edge: org -> dest. When traversing edge ring: Next is CCW, Prev is CW.""" def__init__(self, org, dest): self.org = org # Origin self.dest = dest # Destination self.onext = None# next edge around origin,with same origin self.oprev = None# prev edge around origin,with same origin self.sym = None# edge pointing opposite dest this edge self.deleted = False# Deleted flag
self.name = f'e_{self.org.id}_{self.dest.id}'# for debugging
defcreate_edge(org, dest, edges) -> Quad_Edge: """ Creates an edge, add it dest edges, and return it. """ e = Quad_Edge(org, dest) es = Quad_Edge(dest, org) e.sym, es.sym = es, e # make edges mutually symmetrical e.onext, e.oprev = e, e es.onext, es.oprev = es, es edges.append(e) return e
defupdate_next_prev(e1: Quad_Edge, e2: Quad_Edge): """ Either combines e1 and e2 into a single edge, or seperates them. Which one is determined by the orientation of e1 and e2. """ if e1 == e2: return e1.onext.oprev = e2 e2.onext.oprev = e1 # Swap a.onext and b.onext e1.onext, e2.onext = e2.onext, e1.onext
defconnect(e1: Quad_Edge, e2: Quad_Edge, edges) -> Quad_Edge: """ Connecting destination of e1 with the origin of e2 with an edge O(1) time and O(1) space """ e = create_edge(e1.dest, e2.org, edges) # Maintain the onext and oprev values update_next_prev(e, e1.sym.oprev) update_next_prev(e.sym, e2) return e
defmark_edge_deleted(e: Quad_Edge): """ Delete edge from the edge list O(1) time and O(1) space """ # Update the e.onext' and e.oprev's values update_next_prev(e, e.oprev) update_next_prev(e.sym, e.sym.oprev) # Mark the edge dest be deleted e.deleted = True e.sym.deleted = True
点$d$在三角形$(a,b,c)$外接圆内满足如下不等式： $$ ret =\begin{vmatrix}a_x&a_y&a_x^2+a_y^2&1\[0.3em]b_x&b_y&b_x^2+b_y^2&1\[0.3em]c_x&c_y&c_x^2+c_y^2&1\[0.3em]d_x&d_y&d_x^2+d_y^2&1\end{vmatrix}>0\tag{4-1} $$ 行列式每行减去第四行，得： $$ ret =\begin{vmatrix}a_x-d_x&a_y-d_y&a_x^2+a_y^2-(d_x^2+d_y^2)&0\[0.3em]b_x-d_x&b_y-d_y&b_x^2+b_y^2-(d_x^2+d_y^2)&0\[0.3em]c_x-d_x&c_y-d_y&c_x^2+c_y^2-(d_x^2+d_y^2)&0\[0.3em]d_x&d_y&d_x^2+d_y^2&1\end{vmatrix}>0\tag{4-2} $$ 简化得到下式： $$ ret =\begin{vmatrix}a_x-d_x&a_y-d_y&a_x^2+a_y^2-(d_x^2+d_y^2)\[0.3em]b_x-d_x&b_y-d_y&b_x^2+b_y^2-(d_x^2+d_y^2)\[0.3em]c_x-d_x&c_y-d_y&c_x^2+c_y^2-(d_x^2+d_y^2)\end{vmatrix}>0\tag{4-3} $$ 接下来，将式（4-3）中行列式第三列加上第一列和第二列的$-2d_x$倍，然后将第三列元素整理化成平方差形式，得到下式： $$ ret =\begin{vmatrix}a_x-d_x&a_y-d_y&(a_x-d_x)^2+(a_y-d_y)^2\[0.3em]b_x-d_x&b_y-d_y&(b_x-d_x)^2+(b_y-d_y)^2\[0.3em]c_x-d_x&c_y-d_y&(c_x-d_x)^2+(c_y-d_y)^2\end{vmatrix}>0\tag{4-4} $$ 形如$ad_x=a_x-d_x$替换相同元，得到下式： $$ ret =\begin{vmatrix}ad_x&ad_y&ad_x^2+ad_y^2\[0.3em]bd_x&bd_y&bd_x^2+bd_y^2\[0.3em]cd_x&cd_y&cd_x^2+cd_y^2\end{vmatrix}>0\tag{4-5} $$
defleft_test(p, e): """ Left test for point p relative to the line of edge e. """ a, b = e.org, e.sym.org det1 = (a.x - p.x) * (b.y - p.y) det2 = (a.y - p.y) * (b.x - p.x) return det1 - det2
deftoRight(p, e) -> bool: """Does point p lie to the right of the line of edge e?""" return left_test(p, e) < 0
deftoLeft(p, e) -> bool: """Does point p lie to the left of the line of edge e?""" return left_test(p, e) > 0
# Remove edges that are not part of the triangulation self.edges = [e for e in self.edges if e.deleted isFalse]
definit_points(self): # Validate the input size iflen(self.points) < 2: return
# Sort points by x coordinate, y is a tiebreaker self.points.sort(key = lambda point: (point[0], point[1]))
# Remove duplicates i = 0 while i < len(self.points) - 1: if self.points[i] == self.points[i + 1]: del self.points[i] else: i += 1
# Vertex naming for i, point inenumerate(self.points): self.verticies.append(Vertex(point[0], point[1], i))
defdiv_and_conq_triangulate(self, points): """ Computes the Delaunay triangulation of self.points and returns two edges, le and re, which are the counterclockwise convex hull edge out of the leftmost vertex and the clockwise convex hull edge out of the rightmost vertex, respectively. """ n = len(points) # Base case: 2 points if n == 2: edge = create_edge(points[0], points[1], self.edges) return edge, edge.sym
# Base case: 3 points elif n == 3: # Create edge S[0]-S[1] and edge S[1]-S[2] edge1 = create_edge(points[0], points[1], self.edges) edge2 = create_edge(points[1], points[2], self.edges) update_next_prev(edge1.sym, edge2)
# Create edge S[2]-S[0] if toRight(points[2], edge1): # Right connect(edge2, edge1, self.edges) return edge1, edge2.sym elif toLeft(points[2], edge1): # Left edge3 = connect(edge2, edge1, self.edges) return edge3.sym, edge3 else: # Points are linear return edge1, edge2.sym
# Recurively triangulate the left and right halves else: m = n // 2 ldo, ldi = self.div_and_conq_triangulate(points[:m]) rdi, rdo = self.div_and_conq_triangulate(points[m:]) ldo_r, rdo_r = self.merge(ldo, ldi, rdi, rdo)
return ldo_r, rdo_r
defmerge(self, ldo, ldi, rdi, rdo): """ Takes 2 halves of the triangulation and merges them into a single triangulation. While doing so it uses previosly calculated values of these halves. Reference: https://github.com/alexbaryzhikov/triangulation """ # Compute the upper common tangent of L and R. whileTrue: if toRight(rdi.org, ldi): # Advance dest the next edge on the convex hull of L. ldi = ldi.sym.onext elif toLeft(ldi.org, rdi): # Advance dest the next edge on the convex hull of R. rdi = rdi.sym.oprev else: break
# Create a first cross edge base. base = connect(ldi.sym, rdi, self.edges)
# Adjust ldo and rdo if ldi.org.x == ldo.org.x and ldi.org.y == ldo.org.y: ldo = base if rdi.org.x == rdo.org.x and rdi.org.y == rdo.org.y: rdo = base.sym
# Merge two halves whileTrue: # Locate the first R and L points dest be encountered by the diving bubble. rcand, lcand = base.sym.onext, base.oprev
# If both lcand and rcand are invalid, then base is the lower common tangent. v_rcand, v_lcand = toRight(rcand.dest, base), toRight(lcand.dest, base) ifnot (v_rcand or v_lcand): break
# Delete R edges out of base.dest that fail the circle test. if v_rcand: while toRight(rcand.onext.dest, base) and inCircle(base.dest, base.org, rcand.dest, rcand.onext.dest): t = rcand.onext mark_edge_deleted(rcand) rcand = t
# Symmetrically, delete L edges. if v_lcand: while toRight(lcand.oprev.dest, base) and inCircle(base.dest, base.org, lcand.dest, lcand.oprev.dest): t = lcand.oprev mark_edge_deleted(lcand) lcand = t
# The next cross edge is dest be connected dest either lcand.dest or rcand.dest. # If both are valid, then choose the appropriate one using the in_circle test. ifnot v_rcand or (v_lcand and inCircle(rcand.dest, rcand.org, lcand.org, lcand.dest)): # Add cross edge base from rcand.dest dest base.dest. base = connect(lcand, base.sym, self.edges) else: # Add cross edge base from base.org dest lcand.dest base = connect(base.sym, rcand.sym, self.edges)
