Delaunay三角剖分实现

参考文章：https://www.cnblogs.com/zhiyishou/p/4430017.html

本文使用逐点插入法进行剖分，并使用Unity3D实现。

通过阅读文章《Triangulate》给出的伪代码进行具体编写，我加了些注释：

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 #从边缓冲区中删除所有双重指定的边(例如线段AB和线段BA被当成两个线段存入)
        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

通过Unity的Gizmo和协程，下面是我制作的逐点插入法动图流程：

灰色空心圆 - 输入点

红色十字 - 循环当前点

黄色阶段 - 新的一轮操作开始

红色阶段 - 进行外接圆的筛选

绿色阶段 - 点在外接圆内，拆除三角形，并将边加入边缓冲(EdgeBuffer)

蓝色阶段 - 通过边缓冲的数据和当前顶点，重新组合新的三角形

关于外接圆的计算，取三角形任意两条边的垂线的交点即可得到圆心，圆心和三角形任意一点的线段长度即为半径。

注意：原文链接github的js代码，外接圆半径部分忘了开方(也可能是优化操作，但若参考其实现需要补上)。

更多的具体流程就不写了，理解还是要看伪代码，并且参考原文也很详细。

最后给出实现：

using System;
using System.Collections;
using System.Collections.Generic;
using UnityEngine;

public class DelaunayTriangle : MonoBehaviour
{
    public struct EdgeBuffer
    {
        public Vector2 P0;
        public Vector2 P1;
    }

    public struct Triangle
    {
        public Vector2 A;
        public Vector2 B;
        public Vector2 C;
    }

    public Transform[] points;//通过层级面板放入测试点
    private Bounds mBounds;
    private List<Triangle> mTriangleList = new List<Triangle>();


    private void Update()
    {
        float kSuperTriangleScale = 5f;
        float kStepDuration = 0.5f;

        Vector3 min = Vector3.one * 10000f;
        Vector3 max = Vector3.one * -10000f;

        for (int i = 0; i < points.Length; i++)
        {
            Vector3 point = points[i].position;
            if (point.x < min.x) min = new Vector3(point.x, min.y, min.z);
            if (point.z < min.z) min = new Vector3(min.x, min.y, point.z);

            if (point.x > max.x) max = new Vector3(point.x, max.y, max.z);
            if (point.z > max.z) max = new Vector3(max.x, max.y, point.z);
        }

        mBounds = new Bounds() {min = min, max = max};
        mBounds.size *= kSuperTriangleScale;//此处做法比较粗暴
        mBounds.center += mBounds.extents * 0.5f;

        Vector2 superTriangleA = new Vector2(mBounds.min.x, mBounds.min.z);
        Vector2 superTriangleB = new Vector2(mBounds.min.x, mBounds.max.z);
        Vector2 superTriangleC = new Vector2(mBounds.max.x, mBounds.min.z);

        List<Vector2> vertexList = new List<Vector2>();
        List<EdgeBuffer> edgeBufferList = new List<EdgeBuffer>();
        mTriangleList.Clear();
        mTriangleList.Add(new Triangle() {A = superTriangleA, B = superTriangleB, C = superTriangleC});

        for (int i = 0; i < points.Length; i++)
        {
            Vector3 position = points[i].position;
            vertexList.Add(new Vector2(position.x, position.z));
        }

        vertexList.Add(superTriangleA);
        vertexList.Add(superTriangleB);
        vertexList.Add(superTriangleC);

        for (int i = 0; i < vertexList.Count; i++)//顶点遍历
        {
            Vector2 vertex = vertexList[i];

            edgeBufferList.Clear();

            for (int j = mTriangleList.Count - 1; j >= 0; j--)
            {
                Triangle triangle = mTriangleList[j];

                (Vector2 center, float radius) = Circumcircle(triangle.A, triangle.B, triangle.C);
                //外接圆计算

                if (Vector2.Distance(vertex, center) <= radius)
                {
                    edgeBufferList.Add(new EdgeBuffer() {P0 = triangle.A, P1 = triangle.B});
                    edgeBufferList.Add(new EdgeBuffer() {P0 = triangle.B, P1 = triangle.C});
                    edgeBufferList.Add(new EdgeBuffer() {P0 = triangle.C, P1 = triangle.A});

                    mTriangleList.RemoveAt(j);
                }//若点在外接圆内则移除三角形，并将三角形三个边加入边缓冲
            }

            Dedup(edgeBufferList);//边缓冲去重

            for (int j = 0; j < edgeBufferList.Count; j++)
            {
                EdgeBuffer edgeBuffer = edgeBufferList[j];
                Triangle triangle = new Triangle()
                {
                    A = edgeBuffer.P0,
                    B = edgeBuffer.P1,
                    C = vertex
                };

                mTriangleList.Add(triangle);
            }//重新组合三角形
        }

        for (int j = mTriangleList.Count - 1; j >= 0; j--)
        {
            Triangle triangle = mTriangleList[j];

            if (triangle.A == superTriangleA || triangle.B == superTriangleA || triangle.C == superTriangleA
                || triangle.A == superTriangleB || triangle.B == superTriangleB || triangle.C == superTriangleB
                || triangle.A == superTriangleC || triangle.B == superTriangleC || triangle.C == superTriangleC)
            {
                mTriangleList.RemoveAt(j);
            }
        }//移除连接超三角形的所有三角形
    }

    private void Dedup(List<EdgeBuffer> edgeBufferList)
    {
        for (int i = edgeBufferList.Count - 1; i >= 0; i--)
        {
            for (int j = i - 1; j >= 0; j--)
            {
                EdgeBuffer x = edgeBufferList[i];
                EdgeBuffer y = edgeBufferList[j];

                if ((x.P0 == y.P0 && x.P1 == y.P1) || (x.P0 == y.P1 && x.P1 == y.P0))
                {
                    edgeBufferList.RemoveAt(i);
                    edgeBufferList.RemoveAt(j);

                    i = edgeBufferList.Count - 1;
                    break;
                }
            }
        }
    }

    private (Vector2 center, float radius) Circumcircle(Vector2 a, Vector2 b, Vector3 c)
    {
        float kEps = 0.000001f;

        float x1 = a.x;
        float y1 = a.y;
        float x2 = b.x;
        float y2 = b.y;
        float x3 = c.x;
        float y3 = c.y;
        float fabsy1y2 = Mathf.Abs(y1 - y2);
        float fabsy2y3 = Mathf.Abs(y2 - y3);
        float xc = 0f;
        float yc = 0f;
        float m1 = 0f;
        float m2 = 0f;
        float mx1 = 0f;
        float mx2 = 0f;
        float my1 = 0f;
        float my2 = 0f;
        float dx = 0f;
        float dy = 0f;

        if (fabsy1y2 < kEps)
        {
            m2 = -((x3 - x2) / (y3 - y2));
            mx2 = (x2 + x3) / 2.0f;
            my2 = (y2 + y3) / 2.0f;
            xc = (x2 + x1) / 2.0f;
            yc = m2 * (xc - mx2) + my2;
        }

        else if (fabsy2y3 < kEps)
        {
            m1 = -((x2 - x1) / (y2 - y1));
            mx1 = (x1 + x2) / 2.0f;
            my1 = (y1 + y2) / 2.0f;
            xc = (x3 + x2) / 2.0f;
            yc = m1 * (xc - mx1) + my1;
        }

        else
        {
            m1 = -((x2 - x1) / (y2 - y1));
            m2 = -((x3 - x2) / (y3 - y2));
            mx1 = (x1 + x2) / 2.0f;
            mx2 = (x2 + x3) / 2.0f;
            my1 = (y1 + y2) / 2.0f;
            my2 = (y2 + y3) / 2.0f;
            xc = (m1 * mx1 - m2 * mx2 + my2 - my1) / (m1 - m2);
            yc = (fabsy1y2 > fabsy2y3) ? m1 * (xc - mx1) + my1 : m2 * (xc - mx2) + my2;
        }

        dx = x2 - xc;
        dy = y2 - yc;

        return (center: new Vector2(xc, yc), radius: Mathf.Sqrt(dx * dx + dy * dy));
    }

    private void OnDrawGizmos()
    {
        for (int i = 0; i < mTriangleList.Count; i++)
        {
            Triangle triangle = mTriangleList[i];

            Gizmos.DrawLine(new Vector3(triangle.A.x, 0f, triangle.A.y), new Vector3(triangle.B.x, 0f, triangle.B.y));
            Gizmos.DrawLine(new Vector3(triangle.B.x, 0f, triangle.B.y), new Vector3(triangle.C.x, 0f, triangle.C.y));
            Gizmos.DrawLine(new Vector3(triangle.C.x, 0f, triangle.C.y), new Vector3(triangle.A.x, 0f, triangle.A.y));
        }
    }
}