Delaunay三角刨分算法理解及c#过程实现

0 引言

💻💻AI一下💻💻

三角剖分是什么？

  三角剖分是一种将平面或曲面划分成三角形集合的方法。在二维平面中，给定一个平面区域（可以是多边形等），通过连接区域内的一些点，使得整个区域被分割成若干个三角形，这些三角形彼此相邻，且它们的并集就是原来的平面区域。例如，对于一个简单的矩形，我们可以通过连接矩形的对角线，将其三角剖分成两个三角形。

三角剖分需满足的原则?

• 空外接圆原则（Empty - Circumcircle Criterion）
    这是 Delaunay 三角剖分最核心的原则。对于任意一个三角形，其外接圆内不包含点集中的其他点。也就是说，如果存在一个点在某个三角形的外接圆内部，那么这个三角形就不符合 Delaunay 三角剖分的要求。
    例如，假设有点 A、B、C 构成一个三角形，其外接圆为圆 O。如果在圆 O 内部存在点 D（这个点属于要进行三角剖分的点集），那么三角形 ABC 就需要重新调整，可能会通过连接 AD、BD 或 CD 等来改变三角剖分的结构，以满足外接圆内不包含其他点的条件。

• 最大最小角原则（Max - Min Angle Criterion）
    Delaunay 三角剖分在所有可能的三角剖分中，能够使三角形的最小角达到最大。这一原则使得剖分后的三角形形状相对比较 “规则”，尽量避免出现角度过小（如趋近于 0 度）的三角形，因为狭长的三角形在一些应用中（如有限元分析）可能会导致数值计算的不稳定等问题。
    例如，在给定一组点进行三角剖分的过程中，有多种连接方式可以形成不同的三角剖分结果。Delaunay 三角剖分算法会倾向于选择那种使得每个三角形的最小角尽可能大的连接方式。

• 局部优化原则（Local Optimization Procedure，LOP）
    在构建三角剖分的过程中，当插入一个新的点或者对已有的三角剖分进行调整时，通常会采用局部优化的方法。具体来说，当一个点被插入到现有的三角剖分中，会检查这个点周围的三角形，通过交换对角线等操作来优化局部的三角剖分结构，使其符合 Delaunay 三角剖分的原则。
    例如，在逐步构建三角剖分的过程中，新插入一个点 P 后，它周围可能会形成一些新的三角形。这时，会检查这些三角形及其相邻三角形的外接圆情况，如果发现有不符合空外接圆原则的三角形，就会在局部进行调整，比如交换某些三角形的边，直到满足 Delaunay 三角剖分的要求。

  本篇基于网络资料，对Delaunay三角剖分算法原理进行理解，并借鉴多篇优秀博文，对算法进行了实现，结果与Matlab内置函数进行比对，边缘存在划分差异，可以 满足实际需要。

1 关于三角剖分

  关于Delaunay三角剖分的算法网上有大量博文在介绍相关原理，这里不再展开。但在使用网上获得的各种Delaunay三角剖分算法的时候需要仔细斟酌，因为尝试了许多能下载到的代码，运行结果或多或少都有一些问题。下面是实现三角剖分的几个关键步骤。

  （1）构造一个超级三角形，包含所有散点，放入三角形链表；

  （2）将点集中的散点依次插入，在三角形链表中找出其外接圆包含插入点的三角形，删除影响三角形的公共边，将插入点同影响三角形的全部顶点连接起来；

  （3）根据优化准则对局部新形成的三角形进行优化，将形成的三角形放入 Delaunay 三角形链表；

  （4）循环执行上述步骤，直到所有散点插入完毕。

2 Delaunay三角剖分算法实现及对比

  本篇参考、使用的程序来自网站，网站主页如下图，浏览网站发现里边不仅有C#版本的triangulate源码，还有c、c++、Fortran、Java等版本的三角刨分源码，基于实际用途尝试运行了其中的一些程序，发现个别程序运行结果效果不好。本篇选择使用下图框选区域的程序进行三角刨分，该程序运行效率高、结果也相对准确，大多情况下与Matlab结果基本一致，但是受数据精度影响较大，可能会出现划分错误的情况，使用过程需注意。

  用相同的随机二维散点数据，分别代入Matlab程序和上述C#程序，比较c#程序计算结果是否准确。下面是调用Matlab自带算法进行三角抛分的简单过程：

point = rand(1000,2);
figure
T = delaunay(point);
triplot(T,point(:,1),point(:,2))

  以下是C#版本的三角刨分算法调用过程，C#源码中使用程序仅支持整数类型的坐标点，下面分享的方法是对源码进行的改造，可以支持实数类型的散点进行三角划分，下面仅为定义和关键过程。

## 测试代码，用于驱动计算过程
private void button1_Click(object sender, EventArgs e){Vertex.Clear();Triangle.Clear();// 计算前先释放Graphics g1 = pictureBox2.CreateGraphics();g1.Clear(Color.White);pictureBox2.BackColor = Color.White;string[] lines = File.ReadAllLines(@"basepoint.txt");int N = lines.Length;dVertex tnum;for (int i = 1; i <= N; i++){string line = lines[i - 1];// 拆分行string[] v = line.Split(new char[] { ' ' }, StringSplitOptions.RemoveEmptyEntries);tnum.x = double.Parse(v[0]);tnum.y = double.Parse(v[1]);tnum.z = 0.0;Vertex.Add(tnum);}int tPoints = N;HowMany = 0;HowMany = Triangulate(ref Vertex, ref Triangle);// 绘图Graphics g;g = pictureBox2.CreateGraphics();this.Show();g.CompositingMode = CompositingMode.SourceOver;g.SmoothingMode = SmoothingMode.None;Pen myPen = new Pen(Color.Blue, 1);for (int i = 0; i <= HowMany - 1; i++){double x1 = Vertex[Triangle[i].vv0].x * pictureBox2.Size.Width;double y1 = Vertex[Triangle[i].vv0].y * pictureBox2.Size.Height;double x2 = Vertex[Triangle[i].vv1].x * pictureBox2.Size.Width;double y2 = Vertex[Triangle[i].vv1].y * pictureBox2.Size.Height;double x3 = Vertex[Triangle[i].vv2].x * pictureBox2.Size.Width;double y3 = Vertex[Triangle[i].vv2].y * pictureBox2.Size.Height;g.DrawLine(myPen, (float)x1, (float)y1, (float)x2, (float)y2);g.DrawLine(myPen, (float)x2, (float)y2, (float)x3, (float)y3