二维空间内的三角剖分 -- (给出边缘顶点的例子)

时间 2019-12-09

标签二维空间三角给出边缘顶点例子繁體版

原文原文链接

　　三角剖分的种类不少, 根据不一样需求有各类各样的算法, 这里讲的是按顺序给出边缘顶点以后, 如何对这个顶点集合进行三角剖分.html

好比下面的图示:算法

图一dom

　　给出了边缘点并按顺序排列, 将它剖分红三角形, 虽然有多种方法, 不过咱们只须要得到正确的三角形就好了, 不用关心它剖成怎么样.ide

对于凸多边形(Convex Polygon), 只须要像图中同样, 找一个起始点, 而后按顺序取点跟起始点组成三角形就好了, 很是简单. 代码都懒得写了.测试

　　而对于凹多边形(Concave Polygon)就不能简单划分三角形了, 由于会超出原有多边形的范围 :spa

图二3d

　　HAB组成的三角形明显超出了边界, 不能简单剖分.code

　　PS: 不要把这个和 delaunay triangulation 搞混了, delaunay作的是离散点的三角剖分, 得到的是这些离散点组成的凸多边形.orm

　　我是最近才有这个功能要作, 因此找了一下相关的实现方法, 发现不少人的实现基于很早之前的某篇论文<<三维模型重建中的凹多边形三角剖分>> : 解放军测绘学院学报 1999年9月刊.htm

并且看了下基本都是抄来抄去, 代码很多, 错误很多, 基本都没有实现这个剖分逻辑的...... 论文也是颇有问题, 给出了几个定义, 几个性质, 而后就给出了算法, 至于算法的正确性证实也没有, 我就

直接相信它吧.

　　把论文的算法精简出来就有那么几步 ( 好比顶点列表为 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));