最近计算几何课程要求交一个结课论文,借这个机会,我就认真地学习了下Delaunay三角剖分。
Delaunay三角剖分其实并不是一种算法,它只是给出了一个“好的”三角网格的定义,它的优秀特性是空圆特性和最大化最小角特性,这两个特性避免了狭长三角形的产生,也使得Delaunay三角剖分应用广泛。
空圆特性其实就是对于两个共边的三角形,任意一个三角形的外接圆中都不能包含有另一个三角形的顶点,这种形式的剖分产生的最小角最大。
假如以AC来剖分四边形ABCD,这里B点在三角形ACD的外接圆内,D点在三角形ABC的外接圆内,所以不具有空圆特性
以BD来剖分四边形ABCD,C点不在三角形ABD的外接圆内,A点也不在三角形BCD的外接圆内,具有空圆特性,而且它的最小角要比不满足空圆特性的最小角大。
问题描述
现在给定了平面上的一个点集V,求出它的Delaunay三角剖分
Bowyer逐点插入法
Bowyer逐点插入法是一个很经典的实现方法,网上的资料大多数都是采用的这种算法。这里使用的算法是按照 Paul Bourke在论文中提供的伪代码实现的,这个伪代码也写的很好,非常容易懂
以下摘自他的论文
subroutine triangulate
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 把边数组中所有重复的边都删除,注意这里是把所有的重复边都删除,比如有边e1,e2,e2,e3,删除后应该只剩e1和e3
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
在使用过程中,发现“把超级三角形的顶点加入到顶点链表中”似乎没什么必要,另外,顶点用数组存好一点,因为不需要添加和删除。
以四个点A、B、C、D举例,首先建立一个超级三角形PQR,这个三角形要把所有点都包含进去。
接着分析A点,因为A点在三角形PQR的外接圆内部,所以利用A点将三角形PQR分拆成三个子三角形。
再考虑B点,它只在三角形AQR的内部,再将三角形AQR分拆成三个子三角形。
接着是C点,这时我们有5个三角形,对这5个三角形每一个检查C点在不在它的外接圆内。经过检测,发现它在三角形APR和三角形ABR的外接圆内。
删除三角形APR和三角形ABR的公共边AR,得到四边形ABPR,并将C点与每一个边组成一个三角形。
最后对D点进行分析,它在三角形ABC和三角形BCR的外接圆内,所以应删除公共边BC,再用D点与这两个三角形的其他边形成子三角形。
最后把含有超级三角形的顶点的三角形全部删除,就得到这四个点的三角剖分
我编写程序的时候使用三个点进行测试,有时候结果出不来,是因为在最后一步删除超级三角形的时候是这样的
另外关于这个方法,我在Github上看到了一个实现,作者是jeanbroid,用openGL做显示,里面用了很多C++标准库的东西,让我大开眼界,非常感谢他,不过他用到了双缓存,需要把demo.cpp中的GLUT_SINGLE改成GLUT_DOUBLE才能运行。
分治法
分治法是我在查看维基百科的时候看到的,据维基百科说分治法是效率最高的一种方法。在维基百科的引用中发现了一个介绍delaunay三角剖分的分治法的一个非常好的网站,网站作者是Samuel Peterson,这个网站主要介绍了在给定条件下的Delaunay三角剖分。
分治法进行三角剖分需要对点集中所有点按照x坐标的升序进行排序,x坐标相同时,按照y坐标进行排序,接着将整个点集递归地划分数量相近的两个子集,直到子集中点的数目小于等于3。
先给出一个定义方便后面叙述。对于一条边,如果它的两个端点都在分割线的左边,称它为LL边,一端在左边一端在右边称为LR边,两个端点都在右边称为RR边。
分治法的重点在于如何递归地合并三角网格。
首先是确定一条LR基线,这条线位于最下方,且不与其他边相交,如27边
接着,确定下一条LR边。以右侧三角网为例,按照顺时针方向,计算出∠672,∠872,∠1072,∠972的值,并根据角度大小对顶点进行排序,角度最小的点为第一候选点,以此类推,这里第一候选点为6,第二候选点为8。
作三角形276的外接圆,发现第二候选点8不在该外接圆内,且∠672小于180°,故右侧的候选点为6。候选点的要求是以候选点与该侧基线端点形成的角小于180°,且不包含下一候选点。若大于180°,则该侧无候选点,若包含下一候选点,则考虑下一候选点是否满足该条件。按照这种方法对左侧三角网格进行分析,可得到左侧候选点为点3。
当左右两侧都存在候选点时,如当前情况,这时需要考虑下一条LR边是37还是26。作三角形237的外接圆,可发现右侧候选点6在此外接圆外,而左侧候选点在三角形276外接圆内,故下一条LR边为边37。
当只有一侧存在候选点时,将此候选点与基线另一侧的端点连线作为下一条LR边。
将边37作为下一条基线重复该过程,直到两边都不存在候选点,程序结束。
附一张Bowyer算法最后的显示效果