Delaunay三角剖分
前言
最近的项目刚好用到了这个算法，但是我在网上找到的文章真的让我难以看懂。现在我好不容易理解了，那还是要写一篇文章记录一下，防止我忘了。
环境
windows11
Unity 2022.3.52f1c1
什么是Delaunay三角剖分
在几何中，三角剖分是指将平面对象细分为三角形，并且通过扩展将高维几何对象细分为三角形。对于一个给定的点集，有很多种三角剖分，其中“Delaunay三角剖分”就是其中一种。如果不理解的话，当做了解就好了。我认为重要的是它有什么样的性质和做用。我们知道了这些才能在开发中遇到对应问题时想起这个算法。下面的定义和性质我参考了三角剖分和德劳内三角形生成算法
Delaunay triangle generation algorithm。
定义
在数学和计算几何中，对于给定的平面中的离散点集 𝑃 ，其
Delaunay 三角剖分 DT(𝑃) 满足：
- 空圆性：DT(𝑃)是唯一的（任意四点不能共圆），在DT(𝑃)中，任意三角形的外接圆范围内不会有其它点存在。
- 最大化最小角：在点集 𝑃 可能形成的三角剖分中，DT(𝑃)
所形成的三角形的最小角最大，即更加接近正三角形。
性质
- 最接近：以最接近的三点形成三角形，且各线段（三角形的边）皆不相交。
- 唯一性：不论从区域何处开始构建，最终都将得到一致的结果（点集中任意四点不能共圆）。
- 最优性：任意两个相邻三角形构成的凸四边形的对角线如果可以互换的话，那么两个三角形六个内角中最小角度不会变化。
- 最规则：如果将三角剖分中的每个三角形的最小角进行升序排列，则
Delaunay 三角剖分的排列得到的数值最大。
- 区域性：新增、删除、移动某一个顶点只会影响邻近的三角形。
- 具有凸边形的外壳：三角剖分最外层的边界形成一个凸多边形的外壳。
实现
Delaunay三角的实现方法有很多，下面的实现是我从网络上看到的。这里我们认定所有的x轴正向都是从左往右，y轴正向是从下往上，且下面的算法只考虑二维平面。且默认在圆边界上的点属于圆内的范围。这个方法名为逐点法，我参照了三角剖分算法(delaunay)这篇文章的实现。
在一个给定的点集中，我们按照下面这样的步骤来实现。最终得到我们想要的三角形集合（下面称为tri）
将点集进行排序，排序后的点集是以x轴按从小到大进行排序的。
确定一个超级三角形。这个三角形将所有的点集包括其中。那么在一开始，我们就认定这超级三角形在我们要的三角形集合中。显然实际上这不一定是我们想要的，所以我们可以称这个三角形集合为临时三角形集合（在下面称为tmpTri，代码中也是一个意思）。
遍历点集中的点。对于每一个点，我们都将tmpTri和这些点做一个计算。对于我们tmpTri的每一个三角形而言。
- 3.1
如果现在的点在此三角形外接圆右侧（这里右侧的判断是仅用x轴来判断的），则这个三角形是我们要得到的三角形保存在
tri中并在tmpTri中删除。
- 3.2 如果不满足上述条件且点仍然处于三角形外接圆外，则不做处理。
- 3.3
如果这个点在三角形的外接圆内，我们要将这个三角形从
tmpTri中移除。我们将这三角形的边存储起来（我们设其被存储在edges中）。然后进行下一个三角形的判断。
- 3.4
对
tmpTri中所有的三角形都操作完后，我们对edges中所有的边进行判断，如果有多次出现的边，则在edges中移除掉这个边的所有数据。
- 3.5
最后，我们将当前判断的点和
edges中所有的边组成新的三角形存储在tmpTri中。
所有点位判断完毕后，我们将tmpTri整合到tri中后，删除tri中s所有和超级三角形有关的三角形，最终得到我们想要的三角形集合。
关于3.4，我这边举个例子，假设边A出现了4次，那么edges中所有的边A都要移除掉，也就是说最终处理过的edges中不包含边A。
关于3.5，我发现文章在说明的时候并没有考虑共线的情况。我并没有从算法的描述中想明白，它是如何能够避免在生成新三角形的时候不出现三点共线的情况。因此在后文的代码生成中，我加入了三点共线的判断。
由前面的定义我们可以知道，在Delaunay三角集合中对于任何的一个三角形，其外接圆中只能存在三个点。而在上面的文字说明中，我们却对点在三角形外接圆外的位置进行了额外的判断。这明显是定义中没有说明，按照定义只要三角形外接圆有且只有三个点被包括，则此三角形就是符合条件的三角形。而在所有点位都没有探明的情况下，上文说点在三角形外接圆外且在其右侧时，此三角形是符合条件的三角形。这是因为一开始我们对点集进行了预处理，处理后点集中的点是按x轴从小到大进行排序的。且我们仅使用x轴来判断是否在圆的右侧。如果这个点在外接圆外且在右侧，则后面的点不可能在这个三角形外接圆的内。
除了我们事先对点集进行排序外，我们还做了一个处理就是求得超级三角形。只要满足包含所有点集的三角形都是符合条件的超级三角形，理论上，这个超级三角形可以有无数个。但是我们只要取一个。在我参考的文章中，它的方法是先求得一个包含所有点位的最小矩形（矩形长宽分别平行于x轴和y轴）。然后以x轴为准对矩形进行分割，得到分割后矩形分别取对角线形成一个小三角形，其底是之前小矩形的底。最后对其三角形进行膨胀1倍多一点得到最终的超级三角形。其文章示意图如下：
光是理论上的描述或许并不能让你很好的理解，幸好在我参考的文章:三角剖分算法(delaunay)中，其也有对应的项目。下面我也会放出一部分自己也在Unity上实现的关键代码。下面代码中我虽然用的是Vector3，但是我也只用了x和y。
计算出超级三角形：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
|
protected Vector4 GetBoundary()
{
int maxPointNum = triPoints.Count;
Vector4 boundary = new Vector4(triPoints[0].x, triPoints[0].y, triPoints[0].x, triPoints[0].y);
for (int i = 1; i < maxPointNum; ++i)
{
var point = triPoints[i];
boundary.x = Mathf.Min(boundary.x, point.x);
boundary.y = Mathf.Min(boundary.y, point.y);
boundary.z = Mathf.Max(boundary.z, point.x);
boundary.w = Mathf.Max(boundary.w, point.y);
}
return boundary;
}
protected List<Vector3> CalSuperTri()
{
Vector4 boundary = GetBoundary();
var boundaryCenter = new Vector3((boundary.z + boundary.x) / 2, boundary.y);
var superPoint1 = boundaryCenter;
superPoint1.y += (boundary.w - boundary.y) * 2.1f;
var superPoint2 = boundaryCenter;
superPoint2.x += (boundary.z - boundary.x) / 2 * 2.1f;
superPoint2 = (superPoint2 - superPoint1) * 1.1f + superPoint1;
var superPoint3 = boundaryCenter;
superPoint3.x -= (boundary.z - boundary.x) / 2 * 2.1f;
superPoint3 = (superPoint3 - superPoint1) * 1.1f + superPoint1;
var ans = new List<Vector3>
{
superPoint1,
superPoint2,
superPoint3
};
return ans;
}
|
额外的操作
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
|
protected bool TriCirCheck(Vector3 triP1, Vector3 triP2, Vector3 triP3, Vector3 point)
{
var a1 = 2 * triP2.x - 2 * triP1.x;
var b1 = 2 * triP2.y - 2 * triP1.y;
var c1 = triP2.x * triP2.x + triP2.y * triP2.y - triP1.x * triP1.x - triP1.y * triP1.y;
var a2 = 2 * triP3.x - 2 * triP2.x;
var b2 = 2 * triP3.y - 2 * triP2.y;
var c2 = triP3.x * triP3.x + triP3.y * triP3.y - triP2.x * triP2.x - triP2.y * triP2.y;
var x = (c1 * b2 - c2 * b1) / (a1 * b2 - a2 * b1);
var y = (a1 * c2 - a2 * c1) / (a1 * b2 - a2 * b1);
var rad = (triP1.x - x) * (triP1.x - x) + (triP1.y - y) * (triP1.y - y);
var dis = (point.x - x) * (point.x - x) + (point.y - y) * (point.y - y);
return rad - dis > 1e-6;
}
protected bool InSameLine(Vector3 triP1, Vector3 triP2, Vector3 triP3)
{
var line1 = triP1 - triP2;
var line2 = triP3 - triP2;
return Mathf.Abs(line1.x/line1.y - line2.x / line2.y) <= 0.0001f;
}
|
三角形计算
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
| protected override void Delaunay()
{
if (points.Length > 2)
{
triPoints.Clear();
foreach (var point in points)
triPoints.Add(point.position);
triPoints.Sort((x, y) =>
{
int ans = x.x.CompareTo(y.x);
if (ans == 0)
return x.y.CompareTo(y.y);
return ans;
});
List<int> tmpTri = new List<int>();
tri.Clear();
int maxPointNum = triPoints.Count;
var superTri = CalSuperTri();
triPoints.AddRange(superTri);
tmpTri.Add(triPoints.Count - 3);
tmpTri.Add(triPoints.Count - 2);
tmpTri.Add(triPoints.Count - 1);
for (int i = 0;i < maxPointNum; ++i)
{
List<Vector2Int> edges = new List<Vector2Int>();
for(int j = tmpTri.Count - 1; j >= 0;j -= 3)
{
var point1 = triPoints[tmpTri[j]];
var point2 = triPoints[tmpTri[j - 1]];
var point3 = triPoints[tmpTri[j - 2]];
if (TriCirCheck(point1, point2, point3, triPoints[i]))
{
for(int k = 0;k < 3; ++k)
{
var point1Ind = tmpTri[j - k];
var point2Ind = tmpTri[j - ((k + 1) % 3)];
edges.Add(new Vector2Int(Mathf.Min(point1Indpoint2Ind), Mathf.Max(point1Ind, point2Ind)));
}
for (int k = 0; k < 3; ++k)
{
tmpTri.RemoveAt(j - 2);
}
}
else if(Mathf.Max(point1.x,point2.x,point3.x) triPoints[i].x)
{
tri.Add(tmpTri[j - 2]);
tri.Add(tmpTri[j - 1]);
tri.Add(tmpTri[j]);
tmpTri.RemoveAt(j);
tmpTri.RemoveAt(j - 1);
tmpTri.RemoveAt(j - 2);
}
}
edges.Sort((x, y) =>
{
return x.x.CompareTo(y.x);
});
for (int j = edges.Count - 1; j >= 0; j--)
{
bool need = true;
var curVal = edges[j];
for(int k = j - 1;k > -1; --k)
{
if (edges[k].x != curVal.x)
break;
else if (edges[k].y == curVal.y)
{
if(need)
edges.RemoveAt(j);
edges.RemoveAt(k);
need = false;
}
}
j = Mathf.Min(j, edges.Count);
}
for (int j = edges.Count - 1; j >= 0; j--)
{
if (!InSameLine(triPoints[edges[j].x], triPoints[edges[j].y], triPoints[i]))
{
tmpTri.Add(edges[j].x);
tmpTri.Add(edges[j].y);
tmpTri.Add(i);
}
}
}
tri.AddRange(tmpTri);
for(int i = tri.Count - 1;i >= 0; i -= 3)
{
bool needRemove = false;
for(int j = 0;j < 3; ++j)
{
if (tri[i - j] >= maxPointNum)
{
needRemove = true;
break;
}
}
if (needRemove)
{
for (int j = 0; j < 3; ++j)
{
tri.RemoveAt(i - j);
}
}
}
for(int i = 0;i < 3; ++i)
{
triPoints.RemoveAt(triPoints.Count - 1);
}
}
}
|
结语
我认为上面实现的方法有点复杂了。实际上这篇《三角剖分》文章中也展示了一个分治的方法。我没看懂，但是感觉上这可能比上面说明的方法更快一些。如果你有兴趣可以了解一下。
参考资料
闲言碎语
本来这篇文章应该是我努力去网络上找到各种各样的实现，然后我努力去理解，最终做一个方法的汇总。这样我必然要花很长的一段时间，要么做出一个年度文章（自评），要么就又多了一笔“烂账”。在加上最近受到AI的冲击有点大了。我发觉很多东西真的只要了解原理就可以了。于是这篇文章就这样收尾了。这样也好，我可以去研究新玩意了。
