Delaunay三角剖分實作

參考文章：javascript:void(0)

本文使用逐點插入法進行剖分，并使用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));
        }
    }
}