二維空間內的三角剖分 -- (給出邊緣頂點的例子)

　　三角剖分的種類不少, 根據不一樣需求有各類各樣的算法, 這裏講的是按順序給出邊緣頂點以後, 如何對這個頂點集合進行三角剖分.

好比下面的圖示:

圖一

　　給出了邊緣點並按順序排列, 將它剖分紅三角形, 雖然有多種方法, 不過咱們只須要得到正確的三角形就好了, 不用關心它剖成怎麼樣.

對於凸多邊形(Convex Polygon), 只須要像圖中同樣, 找一個起始點, 而後按順序取點跟起始點組成三角形就好了, 很是簡單. 代碼都懶得寫了.

　　而對於凹多邊形(Concave Polygon)就不能簡單劃分三角形了, 由於會超出原有多邊形的範圍 :

圖二

　　HAB組成的三角形明顯超出了邊界, 不能簡單剖分.

　　PS: 不要把這個和 delaunay triangulation 搞混了, delaunay作的是離散點的三角剖分, 得到的是這些離散點組成的凸多邊形.

　　我是最近纔有這個功能要作, 因此找了一下相關的實現方法, 發現不少人的實現基於很早之前的某篇論文<<三維模型重建中的凹多邊形三角剖分>> : 解放軍測繪學院學報 1999年9月刊.

並且看了下基本都是抄來抄去, 代碼很多, 錯誤很多, 基本都沒有實現這個剖分邏輯的...... 論文也是頗有問題, 給出了幾個定義, 幾個性質, 而後就給出了算法, 至於算法的正確性證實也沒有, 我就

直接相信它吧.

　　把論文的算法精簡出來就有那麼幾步 ( 好比頂點列表爲 List<Vector3> points, 排列順序就是在List中的順序 ) :

　　一. 尋找出一個凸頂點

　　　　1.1 定義 : 臨時從列表中去掉某個頂點, 而後列表剩下的頂點構成新的的多邊形, 這個被去掉的頂點不在新多邊形內部, 就是凸頂點 ( 好比圖二中的B點是凸頂點. 點D被CE組成的多邊形包在裏面, 因此不是凸頂點 )

　　　　1.2 功能 : 若是不是凸頂點, 把這個頂點插回列表原來位置. 找下一個頂點, 若是是凸頂點, 測試是否可劃分頂點.

　　二. 肯定該凸頂點是可劃分頂點

　　　　2.1 定義 : 凸頂點和在列表中相鄰的兩點組成的三角形, 若是不包含列表中的其它頂點, 就是可劃分點 ( 好比圖二中的ABC組成的三角形裏面沒有其它點, B點就是可劃分點. HAB 中有D點, A點就不可劃分 )

　　　　2.2 功能 : 若是不是可劃分頂點, 把這個頂點插回列表原來位置. 找下一個頂點

　　三. 可劃分就把頂點從列表中刪除

　　　　3.1 功能 : 把以前該點組成的三角形記錄到返回值列表中, 它就是一個剖分出來的三角形. 不插回去就等於刪除了.

　　四. 一直重複直到只剩最後3個頂點, 最後3個頂點就是最後一個三角形.

　　　　4.1 功能 : 把最後三個頂點做爲一個三角形加到返回值列表中, 剖分就完成了.

　　

　　看似簡單的邏輯, 拆分出來的功能也不簡單, 先把各個大項中的算法邏輯整理一下:

　　一. 尋找出一個凸頂點

　　最重要的是測試一個點(P)是否在另外一個多邊形內, 這裏按照論文邏輯, 在P點處按某個方向引一條射線, 穿過這個多邊形, 而後看射線與多邊形的邊的交點個數,

若是是奇數就說明點在多邊形內部, 是偶數就說明在多邊形外部. 看看下面一些情形:

　　PS : 射線與線段的交點計算方法在 : 功能實現的數學模型選擇 -- (射線與線段交點的例子)

　　　　　　　　　　　　圖三

　　　　　　　　　　　　　　　　　　圖四

 

 　　上圖的狀況在射線通過線段端點的時候(圖四-1, 圖四-2), 會有歧義, 由於無論怎樣計算, 射線都是跟兩條線段相交, 偶數相交點就判斷在多邊形外部的話是不對的.

還有在射線跟線段平行的狀況(圖四-3), 理論上來講不該該有這樣的狀況, 由於在建模時三個點在一條直線上還不如去掉中間的點, 節省資源, 然而這種情形確實存在, 也會影響判斷.

因此按照論文要給這些歧義的地方添加特殊判斷 :

　　1. 在交點爲線段端點的時候, 按照射線方向計算出線段位於射線的左邊仍是右邊, 只對在左邊或右邊的線段進行計數. 論文中沒有給出證實來講明這個方法的正確性, 若是有時間我會進行一下數學證實, 在附錄補上.

　　　　計算端點與射線相交時, 線段在射線左右邊的方法 : 只要求出相交端點到另外一個端點的向量, 而後跟射線方向叉乘一下, 左右判斷就看大於0或者小於0來判斷便可. 大於0在左邊, 小於0在右邊.

　　2. 在射線與某條線段重合的時候, 不進行計數, 也就是說該點的狀況不受這條線段影響. 這個結論也沒有通過證實, 不過根據個人直覺它是對的, 由於像上面說的建模時原本就應該去掉的頂點, 不參與計算也對.(2019.06.19)

　　圖五 -- 反對上面的直覺, 射線跟線段重疊並不必定是多餘的點!(2019.06.20)  若是有時間我會進行一下數學證實爲什麼不進行計數, 在附錄補上.

　　PS : 這些理論在數學上即便可以證實是正確的, 但是用在計算機系統裏面就不必定正確, 由於小數點的精度問題, 在這些計算中確定會出現錯誤, 好比射線跟線段端點相交, 若是有一點誤差就是不相交或是相交在線段內了.

 　　3. 怎樣肯定點P的射線方向, 讓它穿過多邊形? 最簡單的方法就在多邊形上隨便選一個點, 而後點P到該點的方向便可, 爲了減小計算偏差, 我會選擇隨機某條邊的中間點來做爲射線方向.

     

　　　　　　　　　　圖六　　　　　　　　　　　　　　　　　　　　　　　　圖七

　　對圖2的多邊形頂點進行凸頂點測試, A點與B點都在新多邊形以外, 都是凸頂點.

 

　　二. 肯定該凸頂點是可劃分頂點

　　1. 當一個頂點是凸頂點的時候, 還須要確認它是可劃分頂點, 可劃分頂點就是剖出的三角形不包含剩下的各個頂點的狀況. 看下圖:

　　　　　　　　　　　　　　　　　　　　　　　　圖八

　　這是對圖二的多邊形進行的一次三角剖分, 這裏A, B頂點都是是凸頂點, 用HAB組成了一個新三角形, 而後剩下的頂點是C,D,E,F,G顯然D點在HAB裏面, 因此如今A不是一個可劃分點, 而後下一個點B, ABC沒有包含其它頂點, B點就是可劃分點.

核心算法就是怎樣測試其他點是否在要剖分的三角形以內. 

　　2. 怎樣計算一個點是否在三角形以內, 這裏使用的是一個相似幾何的方法, 看下圖:

　　　　　　　　　　　　圖九

　　若是點P在一個三角形以內, 能夠得出 : 角A, 角B, 角C 任意角都小於180度.

　　數學證實 : 角A + 角B + 角C = 360度, P與各個端點構成的三角形好比BPA, 最大的角度就是P點在AB的線段上, 那麼角BPA最大就是180度, 最小就是角BCA, 因此 角BPA  >= 角BCA && 角BPA <= 180度, 那麼其他兩角的和就大於等於180度.

　　計算邏輯 : 使用叉乘來計算 ( 用⊙ 表明叉乘符號 ) =>

　　cross1 = (向量PA ⊙ 向量PB)

　　cross2 = (向量PB ⊙ 向量PC)

　　cross3 = (向量PC ⊙ 向量PA)

　　這裏用的是X-Z平面做爲二維平面, 那麼叉積獲得的就是y軸的值, 用Vector3向量的話很容易計算, 那麼若是點在三角形內, 能夠得出它們的叉積在y軸的值必定同時都是負的或者同時都是正的. 只有點在三角形以外的時候纔會有超過180度的角,

纔會有同時有正負的狀況出現. 因此斷定就能夠很簡單: (cross1.y * cross2.y) >= 0 && (cross2.y * cross3.y) >= 0; 判斷他們正負號都相同就好了.

 

　　三. 可劃分就把頂點從列表中刪除

　　當可劃分頂點找到後, 就能夠把該頂點從原有列表中刪除了. 按照圖七的劃分, 記錄下剖分三角形ABC, 頂點列表刪除B點, 剩下的列表若是大於三個頂點, 繼續重複這個過程便可. 直到最後三個點, 把三個點也加入返回列表就完成了.

 

　　理論過程就到這裏, 下來就是代碼實現了.

　　說明一下我用的是Unity3D, 原始頂點列表做爲輸入, 返回值輸出的是剖分後的三角形頂點, 而不是組成三角形的Index, 由於Mesh須要normal的數量跟頂點數量一致, 若是Mesh使用原始頂點列表做爲頂點的話(好比正方體只有8個頂點的話),

它自動計算出來的normal也是隻有8個, 也就是說它把各個共享頂點的向量進行了插值, 光照會出錯的.

 

    public static class Triangulation
    {
        #region Main Funcs
        /// <summary>
        /// 泛用的三角剖分
        /// </summary>
        /// <param name="vertices"> 按照順序排列的邊界點 </param>
        /// <returns></returns>
        public static List<Vector3> GenericTriangulate(List<Vector3> vertices)
        {
            if(vertices.Count < 3)
            {
                throw new System.Exception("No Enough Points!!!");
            }
            List<Vector3> tempPoints = new List<Vector3>(vertices);
            List<Vector3> points = new List<Vector3>();
            int convexIndex = -1;
            while(tempPoints.Count > 3 && ((convexIndex = SearchConvexIndex(tempPoints)) != -1))
            {
                var p0 = convexIndex == 0 ? tempPoints[tempPoints.Count - 1] : tempPoints[convexIndex - 1];
                var p1 = tempPoints[convexIndex];
                var p2 = (convexIndex == tempPoints.Count - 1) ? tempPoints[0] : tempPoints[convexIndex + 1];
                points.Add(p0);
                points.Add(p1);
                points.Add(p2);

                tempPoints.RemoveAt(convexIndex);
                Debug.Log("被剔除頂點 : " + convexIndex);
            }
            points.AddRange(tempPoints);
            return points;
        }
        #endregion

        #region Help Funcs
        /// <summary>
        /// Search a Convex point from vertices
        /// </summary>
        /// <param name="vertices"></param>
        /// <returns></returns>
        public static int SearchConvexIndex(List<Vector3> vertices)
        {
            if(vertices.Count > 3)
            {
                for(int i = 0; i < vertices.Count; i++)
                {
                    if(IsPointInsideGeometry(ref vertices, i) == false)
                    {
                        int index0 = i == 0 ? vertices.Count - 1 : i - 1;
                        int index1 = i;
                        int index2 = i == vertices.Count - 1 ? 0 : i + 1;
                        if(IsAnyPointInsideGeometryTriangle(vertices, index0, index1, index2))
                        {
                            Debug.LogWarning("有其它頂點在剔除三角形內, 沒法剔除:" + i);
                            continue;
                        }
                        return i;
                    }
                }
            }
            return -1;
        }

        /// <summary>
        /// 檢查是否相等, 解決計算機計算偏差
        /// </summary>
        /// <param name="a"></param>
        /// <param name="b"></param>
        /// <returns></returns>
        public static bool IsTheSameValue(float a, float b)
        {
            const double Epsilon = 1e-5;
            var val = System.Math.Abs(a - b);
            return val <= Epsilon;
        }

        /// <summary>
        /// 射線與線段相交性判斷
        /// </summary>
        /// <param name="ray"></param>
        /// <param name="p1"></param>
        /// <param name="p2"></param>
        /// <param name="point"></param>
        /// <param name="isParallel"></param>
        /// <returns></returns>
        public static bool RayAndLineIntersection(Ray ray, Vector3 p1, Vector3 p2, out Vector3 point, out bool isParallel)
        {
            point = Vector3.zero;
            isParallel = false;
            Vector3 p3 = new Vector3(ray.origin.x, 0, ray.origin.z);
            Vector3 p4 = ray.GetPoint(1.0f);
            p4.y = 0;

            float a1, b1, c1;
            float[] variables1;
            GeometryMath.MultivariateLinearEquation(new float[2] { p1.x, p1.z }, new float[2] { p2.x, p2.z }, out variables1);
            a1 = variables1[0];
            b1 = variables1[1];
            c1 = variables1[2];

            float a2, b2, c2;
            float[] variables2;
            GeometryMath.MultivariateLinearEquation(new float[2] { p3.x, p3.z }, new float[2] { p4.x, p4.z }, out variables2);
            a2 = variables2[0];
            b2 = variables2[1];
            c2 = variables2[2];

            DeterminantEquation determinantEquation = new DeterminantEquation(2);
            determinantEquation.determinant.SetRow(a1, b1);
            determinantEquation.determinant.SetRow(a2, b2);
            var variables = determinantEquation.CaculateVariables(-c1, -c2);
            if(variables == null)
            {
                // no result -- check is Parallel line with same
                float ea, eb;
                ea = a1 * c2 - a2 * c1;
                eb = b1 * c2 - b2 * c1;
                if((ea == 0) && (eb == 0))
                {
                    point = p1;
                    isParallel = true;
                    Debug.Log("Determinant No Result, it is Parallel Line.");
                    return true;
                }
                return false;
            }
            point = new Vector3(variables[0], 0, variables[1]);

            var rayDir = ray.direction;
            rayDir.y = 0;

            bool intersect = Vector3.Dot(point - p3, rayDir) > 0;
            if(intersect)
            {
                bool xClamp = IsTheSameValue(point.x, p1.x) || IsTheSameValue(point.x, p2.x) || (point.x >= Mathf.Min(p1.x, p2.x) && point.x <= Mathf.Max(p1.x, p2.x));
                bool zClamp = IsTheSameValue(point.z, p1.z) || IsTheSameValue(point.z, p2.z) || (point.z >= Mathf.Min(p1.z, p2.z) && point.z <= Mathf.Max(p1.z, p2.z));
                intersect = xClamp && zClamp;
            }
            return intersect;
        }

        /// <summary>
        /// 判斷點是否在多邊形內
        /// </summary>
        /// <param name="points"></param>
        /// <param name="checkPointIndex"></param>
        /// <returns></returns>
        private static bool IsPointInsideGeometry(ref List<Vector3> points, int checkPointIndex)
        {
            var p = points[checkPointIndex];
            var point = new Vector3(p.x, 0, p.z);

            points.RemoveAt(checkPointIndex);

            int interCount = 0;
            int randIndex = UnityEngine.Random.Range(0, points.Count - 1);
            var randPoint = (points[randIndex] + points[randIndex + 1]) * 0.5f;
            randPoint.y = 0.0f;

            Vector3 randDir = (randPoint - point).normalized;
            Ray ray = new Ray(point, randDir);          // random direction

            for(int i = 0; i < points.Count; i++)
            {
                var p1 = points[i];
                var p2 = (i + 1) >= points.Count ? points[0] : points[i + 1];
                Vector3 intersectPoint;
                bool isParallel;
                if(RayAndLineIntersection(ray, p1, p2, out intersectPoint, out isParallel))
                {
                    if(isParallel == false)
                    {
                        if(p1 == intersectPoint || p2 == intersectPoint)
                        {
                            var dir = (p1 == intersectPoint) ? (p2 - p1) : (p1 - p2);
                            var cross = Vector3.Cross(dir, randDir);
                            if(cross.y < 0)
                            {
                                continue;
                            }
                        }
                        interCount++;
                    }
                }
            }
            points.Insert(checkPointIndex, p);
            bool inside = ((interCount % 2) == 1);
            if(inside == false)
            {
                Debug.Log("Index 不在多邊形內, 嘗試剔除此頂點 : " + checkPointIndex + " 測試頂點:" + randPoint);
            }
            return inside;
        }

        /// <summary>
        /// 判斷點是否在三角形內 -- 快速
        /// </summary>
        /// <param name="points"></param>
        /// <param name="triangleIndex0"></param>
        /// <param name="triangleIndex1"></param>
        /// <param name="triangleIndex2"></param>
        /// <returns></returns>
        private static bool IsAnyPointInsideGeometryTriangle(List<Vector3> points, int triangleIndex0, int triangleIndex1, int triangleIndex2)
        {
            var p0 = points[triangleIndex0];
            var p1 = points[triangleIndex1];
            var p2 = points[triangleIndex2];
            for(int i = 0; i < points.Count; i++)
            {
                if(i != triangleIndex0 && i != triangleIndex1 && i != triangleIndex2)
                {
                    var point = points[i];
                    if(IsPointInsideTriangle(p0, p1, p2, point))
                    {
                        return true;
                    }
                }
            }
            return false;
        }

        /// <summary>
        /// a fast way to check a point in triangle
        /// </summary>
        /// <param name="triangleP0"></param>
        /// <param name="triangleP1"></param>
        /// <param name="triangleP2"></param>
        /// <param name="point"></param>
        /// <returns></returns>
        public static bool IsPointInsideTriangle(Vector3 triangleP0, Vector3 triangleP1, Vector3 triangleP2, Vector3 point)
        {
            point.y = 0.0f;
            triangleP0.y = 0.0f;
            triangleP1.y = 0.0f;
            triangleP2.y = 0.0f;

            var dir1 = triangleP0 - point;
            var dir2 = triangleP1 - point;
            var dir3 = triangleP2 - point;

            var cross1 = Vector3.Cross(dir1, dir2);
            var cross2 = Vector3.Cross(dir2, dir3);
            var cross3 = Vector3.Cross(dir3, dir1);

            bool inside = (cross1.y * cross2.y) >= 0 && (cross2.y * cross3.y) >= 0;
            return inside;
        }
        #endregion
    }

 

　　最後直接測試一下生成一個Mesh

        List<Vector3> points = new List<Vector3>();
        points.Add(new Vector3(0.5f, 0, 0.5f));
        points.Add(new Vector3(0.5f, 0, 1.5f));
        points.Add(new Vector3(1.5f, 0, 1.5f));
        points.Add(new Vector3(1.5f, 0, -1f));
        points.Add(new Vector3(-1.5f, 0, -1f));
        points.Add(new Vector3(-1.5f, 0, 1.5f));
        points.Add(new Vector3(-0.5f, 0, 1.5f));
        points.Add(new Vector3(-0.5f, 0, 0.5f));