C#，计算几何，随机点集之三角剖分的德劳内（Delaunay）算法的源代码

一、三角剖分Delaunay算法简介

���集的三角剖分（Triangulation），对数值分析（比如有限元分析）以及图形学来说，都是极为重要的一项预处理技术。尤其是Delaunay三角剖分，由于其独特性，关于点集的很多种几何图都和Delaunay三角剖分相关，如Voronoi图，EMST树，Gabriel图等。Delaunay三角剖分有最大化最小角，“最接近于规则化的“的三角网和唯一性（任意四点不能共圆）两个特点。

 EMST(Euclidean minimum spanning tree)

Delaunay 三角剖分广泛应用于许多不同应用程序中的科学计算。虽然有大量的计算三角剖分的算法，但 Delaunay 三角剖分以其实用的几何属性广受欢迎。

 Gabriel Graph

基本属性是 Delaunay 规则。如果是二维三角剖分，通常将其称为空外接圆规则。对于一组二维点而言，这些点的 Delaunay 三角剖分可确保与每个三角形相关的外接圆的内部都不包含其他点。这种三角剖分便是 Delaunay 三角剖分。

Delaunay 三角剖分堪称“外形整齐”，原因在于为满足空外接圆属性，优先选择带有较大内角的三角形，而不是带有较小内角的三角形。非 Delaunay 三角剖分中的三角形在顶点 V2 和 V4 处呈锐角。如果将 {V2, V4} 边替换为连接 V1 和 V3 的边，会实现最小角的最大化并且使得该三角剖分变为 Delaunay 三角剖分。另外，Delaunay 三角剖分将最近邻点的点连接在一起。这两个特征（外形整齐和最近邻点关系）在实践中具有重要的作用，有助于促进在散点数据插值中使用 Delaunay 三角剖分。

虽然 Delaunay 属性定义明确，但存在退化点集时三角剖分的拓扑并不唯一。在二维中，4 个或更多特征点位于同一圆中时会引发退化。例如，正方形的顶点不具有唯一的 Delaunay 三角剖分。

二、三角剖分Delaunay算法的源代码

namespace Legalsoft.Truffer.Algorithm{    public struct Vertex    {        public int x;        public int y;        public int z;    }    public struct Triangle    {        public int vv0;        public int vv1;        public int vv2;    }    public class Delaunay    {        public const int MaxVertices = 500;        public const int MaxTriangles = 1000;        public Vertex[] Vertex = new Vertex[MaxVertices];        public Triangle[] Triangle = new Triangle[MaxTriangles];        private bool InCircle(int xp, int yp, int x1, int y1, int x2, int y2, int x3, int y3, double xc, double yc, double r)        {            double eps;            double m1;            double m2;            double mx1;            double mx2;            double my1;            double my2;            double dx;            double dy;            double rsqr;            double drsqr;            eps = 0.000000001;            if (Math.Abs(y1 - y2)  ymax)                {                    ymax = Vertex[i].y;                }            }            dx = Convert.ToDouble(xmax) - Convert.ToDouble(xmin);            dy = Convert.ToDouble(ymax) - Convert.ToDouble(ymin);            if (dx > dy)            {                dmax = dx;            }            else            {                dmax = dy;            }            xmid = (xmax + xmin) / 2;            ymid = (ymax + ymin) / 2;            Vertex[nvert + 1].x = Convert.ToInt64(xmid - 2 * dmax);            Vertex[nvert + 1].y = Convert.ToInt64(ymid - dmax);            Vertex[nvert + 2].x = xmid;            Vertex[nvert + 2].y = Convert.ToInt64(ymid + 2 * dmax);            Vertex[nvert + 3].x = Convert.ToInt64(xmid + 2 * dmax);            Vertex[nvert + 3].y = Convert.ToInt64(ymid - dmax);            Triangle[1].vv0 = nvert + 1;            Triangle[1].vv1 = nvert + 2;            Triangle[1].vv2 = nvert + 3;            Complete[1] = false;            ntri = 1;            for (i = 1; i