Delaunay三角剖分实现

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

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));
        }
    }
}