Bowyer-Watson快速入门：
function BowyerWatson (pointList)
// pointList is a set of coordinates defining the points to be triangulated
triangulation := empty triangle mesh data structure
add super-triangle to triangulation // must be large enough to completely contain all the points in pointList
for each point in pointList do // add all the points one at a time to the triangulation
badTriangles := empty set
for each triangle in triangulation do // first find all the triangles that are no longer valid due to the insertion
if point is inside circumcircle of triangle
add triangle to badTriangles
polygon := empty set
for each triangle in badTriangles do // find the boundary of the polygonal hole
for each edge in triangle do
if edge is not shared by any other triangles in badTriangles
add edge to polygon
for each triangle in badTriangles do // remove them from the data structure
remove triangle from triangulation
for each edge in polygon do // re-triangulate the polygonal hole
newTri := form a triangle from edge to point
add newTri to triangulation
for each triangle in triangulation // done inserting points, now clean up
if triangle contains a vertex from original super-triangle
remove triangle from triangulation
return triangulation
这相对容易地实现到python 3.7.2和pygame 1.7〜中。
在此之后，我在另一篇文章中发现this comment，指出大多数Bowyer-Watson算法都缺乏对无穷远点的顶点的计算以及克服这一问题的实现方法：
检查是否有任何三角形顶点位于无穷大处。换一种说法：
检查三角形是否与边界三角形共享某些顶点。
如果要共享所有三个顶点：琐碎的。
如果共享零个顶点：经典方法-检查是否
点到外接心的距离短于外接半径。
如果共享一个顶点：检查点是否位于顶点的左侧/右侧
由其他两个顶点定义的线。 one vertex in infinity
如果共享两个顶点：检查点是否位于左侧/右侧
这两个顶点定义的线的线，但移至第三个
点。换句话说：您仅从直线获取斜率矢量
在这些共享顶点之间移动它，以便线通过
通过第三点。 two vertices in infinity
因此，我创建了一个新的
Point()类，而不是内置的pygame.Vector2()类，并添加了isInfinite布尔元素。对于
Triangle()类，它接受3个Point()，然后分为2个列表：finiteVertices和infiniteVertices基于每个isInfinite中的Point元素，并根据有多少个来执行上述过程。元素在每个列表中。如果不把我的三角剖分变成意大利面条，那就太好了。
我的代码是：
import pygame
import pygame.gfxdraw
import math
import random
pygame.init()
def Sign(value):
if value < 0:
return -1
if value > 0:
return 1
if value == 0:
return 0
def SideOfLineOfPoint(x1,y1,x2,y2,posX,posY):
d = (posX-x1)*(y2-y1) - (posY-y1)*(x2-x1)
return Sign(d)
def LineIsEqual(line1,line2): # Detect congruence of line, no matter which end is defined first
if (line1[0] == line2[0] and line1[1] == line2[1]) or (line1[0] == line2[1] and line1[1] == line2[0]):
return True
return False
class Point:
def __init__(self,x,y,isInfinite):
self.x = x
self.y = y
self.isInfinite = isInfinite
def distanceTo(self,other):
return math.sqrt( (self.x-other.x)**2 + (self.y-other.y)**2 )
class Triangle:
def __init__(self,a,b,c):
self.vertices = [a,b,c] # a,b,c are vertices defining the triangle
self.edges = [[a,b],
[b,c],
[c,a]] # Manually defining all edges of triangle ([])
self.CalculateCircumcenter()
self.infiniteVertices = []
self.finiteVertices = []
for vertex in self.vertices:
if vertex.isInfinite:
self.infiniteVertices.append(vertex)
else:
self.finiteVertices.append(vertex)
def CalculateCircumcenter(self): # Copied from Delaunator page
a = [self.vertices[0].x , self.vertices[0].y]
b = [self.vertices[1].x , self.vertices[1].y]
c = [self.vertices[2].x , self.vertices[2].y]
ad = a[0] * a[0] + a[1] * a[1]
bd = b[0] * b[0] + b[1] * b[1]
cd = c[0] * c[0] + c[1] * c[1]
D = 2 * (a[0] * (b[1] - c[1]) + b[0] * (c[1] - a[1]) + c[0] * (a[1] - b[1]))
self.circumcenter = Point(1 / D * (ad * (b[1] - c[1]) + bd * (c[1] - a[1]) + cd * (a[1] - b[1])),
1 / D * (ad * (c[0] - b[0]) + bd * (a[0] - c[0]) + cd * (b[0] - a[0])),
False)
def IsPointInCircumcircle(self,point):
if len(self.infiniteVertices) == 3:
return True # Any point is within the circumcircle if all therr vertices are infinite
elif len(self.infiniteVertices) == 2: # If two infinite vertices: check if point lies to the left/right of line defined by these two vertices but shifted to the third point.
x1 = self.finiteVertices[0].x
y1 = self.finiteVertices[0].y
x2 = self.infiniteVertices[0].x - self.infiniteVertices[1].x + x1
y2 = self.infiniteVertices[0].y - self.infiniteVertices[1].y + y1
sideOfLineOfVertex = SideOfLineOfPoint(x1,y1,x2,y2,point.x,point.y)
sideOfLineOfPoint = SideOfLineOfPoint(x1,y1,x2,y2,self.infiniteVertices[0].x,self.infiniteVertices[0].x)
if sideOfLineOfVertex == sideOfLineOfPoint:
return False
else:
return True
elif len(self.infiniteVertices) == 1: # If one infinite vertex: check if point lies to the left/right of line defined by the other two vertices.
x1 = self.finiteVertices[0].x
y1 = self.finiteVertices[0].y
x2 = self.finiteVertices[1].x
y2 = self.finiteVertices[1].y
sideOfLineOfVertex = SideOfLineOfPoint(x1,y1,x2,y2,point.x,point.y)
sideOfLineOfPoint = SideOfLineOfPoint(x1,y1,x2,y2,self.infiniteVertices[0].x,self.infiniteVertices[0].x)
if sideOfLineOfVertex == sideOfLineOfPoint:
return False
else:
return True
elif len(self.infiniteVertices) == 0: # For triangle with finite vertices
if self.vertices[0].distanceTo(self.circumcenter) > point.distanceTo(self.circumcenter):
return True # If point is closer to circumcenter than any vertices, point is in circumcircle
else:
return False
def HasVertex(self,point):
if point in self.vertices:
return True
return False
def Show(self,screen,colour):
for edge in self.edges:
pygame.draw.aaline(screen,colour,(edge[0].x,edge[0].y),(edge[1].x,edge[1].y))
class DelaunayTriangulation:
def __init__(self,points,width,height):
self.triangulation = [] # Create empty list
self.superTriangleA = Point(-100,-100,True)
self.superTriangleB = Point(2*width+100,-100,True)
self.superTriangleC = Point(-100,2*height+100,True)
superTriangle = Triangle(self.superTriangleA,self.superTriangleB,self.superTriangleC)
self.triangulation.append(superTriangle) # Create super-triangle
for point in points: # For every single point to be triangulated
self.addPoint(point)
def addPoint(self,point):
invalidTriangles = [] # Invalid triangle list
for triangle in self.triangulation: # For every existing triangle
if triangle.IsPointInCircumcircle(point): # If new point is in the circumcenter of triangle
invalidTriangles.append(triangle) # Triangle is invalid and added to invalid triangle list
polygon = [] # List for remaining edges after removal of invalid triangles
for triangle in invalidTriangles: # For every invalid triangle
for triangleEdge in triangle.edges: # For each invalid triangle's edges
isShared = False # Assume no edges are shared
for other in invalidTriangles: # For every other invalid triangle
if triangle == other: # If both are the same triangle
continue
for otherEdge in other.edges: # For every edge in other triangle
if LineIsEqual(triangleEdge,otherEdge):
isShared = True # Congruent edges are shared
if isShared == False: # Only append edges not shared by invalid triangles to polygon hole
polygon.append(triangleEdge)
for triangle in invalidTriangles: # Remove invalid triangles
self.triangulation.remove(triangle)
for edge in polygon:
newTriangle = Triangle(edge[0],edge[1],point) # Retriangulate based on edges of polygon hole and point
self.triangulation.append(newTriangle)
def Show(self,screen,colour):
self.shownTriangulation = self.triangulation
superTriangles = [] # List for triangles that are part of super-triangle
for triangle in self.triangulation:
if (triangle.HasVertex(self.superTriangleA) or triangle.HasVertex(self.superTriangleB) or triangle.HasVertex(self.superTriangleC)) and (triangle in self.triangulation):
superTriangles.append(triangle) # Add triangles that have any super-triangle vertex
for triangle in superTriangles:
self.triangulation.remove(triangle) # Remove super-triangles
for triangle in self.shownTriangulation:
triangle.Show(screen,colour)
background = 20,40,100
red = 255,0,0
white = 255,255,255
width = int(500)
height = int(500)
amount = int(5)
screen = pygame.display.set_mode((width,height))
screen.fill(background)
points = []
for i in range(amount):
x = random.randint(1,width-1)
y = random.randint(1,height-1)
points.append(Point(x,y,False))
delaunay = DelaunayTriangulation(points,width,height)
delaunay.Show(screen,white)
pygame.display.update()
在我看来，可能导致此问题的函数是
Triangle.IsPointInCircumcircle()和SideOfLineOfPoint()，尽管同样有可能原始算法并非一开始就支持无限顶点计算。如果取消了整个无限顶点计算并使用了正常的外接圆检测，该代码将起作用，尽管这将比我的目标落后一步。
我希望有人能在我的代码中发现任何错误来解决此问题，甚至只是指向正确的方向开始调试此问题。
提前致谢。
最佳答案
通过避免昂贵的math.sqrt操作，可以最大程度地提高性能。
而不是比较点之间的Euclidean distances，而是比较距离的平方：
class Point:
# [...]
def distanceToSquare(self,other):
dx, dy = self.x-other.x, self.y-other.y
return dx*dx + dy*dy
class Triangle:
# [...]
def IsPointInCircumcircle(self,point):
return (self.vertices[0].distanceToSquare(self.circumcenter) >
point.distanceToSquare(self.circumcenter))
除此之外，您的代码中还有一些简单的错字。
x的infiniteVertices[0]组件两次传递给方法SideOfLineOfPoint，但是y组件丢失（在两种情况下）：SideOfLineOfPoint(...,self.infiniteVertices[0].x,self.infiniteVertices[0].x)
此外，如果点在同一侧，方法
IsPointInCircumcircle必须返回True。您执行相反的操作：if sideOfLineOfVertex == sideOfLineOfPoint:
return False
else:
return True
我建议在方法
IsPointInCircumcircle中颠倒案例的顺序。当添加第一个点时，len(self.infiniteVertices) == 3情况仅发生一次。比较而言，len(self.infinite Vertices) == 0是最常见的情况，尤其是当点数增加时。查看更正方法的结果（20个随机点）：
class Triangle:
# [...]
def IsPointInCircumcircle(self,point):
if len(self.infiniteVertices) == 0: # For triangle with finite vertices
if self.vertices[0].distanceToSquare(self.circumcenter) > point.distanceToSquare(self.circumcenter):
return True # If point is closer to circumcenter than any vertices, point is in circumcircle
else:
return False
elif len(self.infiniteVertices) == 1: # If one infinite vertex: check if point lies to the left/right of line defined by the other two vertices.
x1 = self.finiteVertices[0].x
y1 = self.finiteVertices[0].y
x2 = self.finiteVertices[1].x
y2 = self.finiteVertices[1].y
sideOfLineOfVertex = SideOfLineOfPoint(x1,y1,x2,y2,point.x,point.y)
sideOfLineOfPoint = SideOfLineOfPoint(x1,y1,x2,y2,self.infiniteVertices[0].x,self.infiniteVertices[0].y)
if sideOfLineOfVertex == sideOfLineOfPoint:
return True
else:
return False
elif len(self.infiniteVertices) == 2: # If two infinite vertices: check if point lies to the left/right of line defined by these two vertices but shifted to the third point.
x1 = self.finiteVertices[0].x
y1 = self.finiteVertices[0].y
x2 = self.infiniteVertices[0].x - self.infiniteVertices[1].x + x1
y2 = self.infiniteVertices[0].y - self.infiniteVertices[1].y + y1
sideOfLineOfVertex = SideOfLineOfPoint(x1,y1,x2,y2,point.x,point.y)
sideOfLineOfPoint = SideOfLineOfPoint(x1,y1,x2,y2,self.infiniteVertices[0].x,self.infiniteVertices[0].y)
if sideOfLineOfVertex == sideOfLineOfPoint:
return True
else:
return False
return True # Any point is within the circumcircle if all there vertices are infinite
关于python - 尝试对无限远处的顶点执行外接圆计算时，Bowyer-Watson的三角剖分不正确，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/58203812/