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

比如下面的图示:

图一

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

对于凸多边形(Convex Polygon), 只需要像图中一样, 找一个起始点, 然后按顺序取点跟起始点组成三角形就行了, 非常简单. 代码都懒得写了.

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

图二

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

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

　　我是最近才有这个功能要做, 所以找了一下相关的实现方法, 发现很多人的实现基于很早以前的某篇论文 : 解放军测绘学院学报 1999年9月刊.

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

直接相信它吧.

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

　　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; 判断他们正负号都相同就行了( 判断同号不能这样算! 2019.06.25 ).

还是要老老实实判断 (cross1.y >= 0 && cross2.y >= 0 && cross3.y >= 0) || (cross1.y 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 /// /// Search a Convex point from vertices /// /// /// public static int SearchConvexIndex(List 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; } /// /// 检查是否相等, 解决计算机计算误差 /// /// /// /// public static bool IsTheSameValue(float a, float b) { const double Epsilon = 1e-5; var val = System.Math.Abs(a - b); return val 0) { point = Vector3.SqrMagnitude(p1 - ray.origin) < Vector3.SqrMagnitude(p2 - ray.origin) ? p1 : p2; return true; } } Debug.Log("it is Parallel Line, and has no intersection"); } return false; } point = new Vector3(variables[0], 0, variables[1]); bool intersect = Vector3.Dot(point - p3, rayDir) > 0; if(intersect) { intersect = IsPointInsideLine(point, p1, p2); } return intersect; } /// /// 判断点是否在多边形内 /// /// /// /// private static bool IsPointInsideGeometry(ref List 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(IsTheSamePoint(p1, intersectPoint) || IsTheSamePoint(p2, intersectPoint)) { var dir = IsTheSamePoint(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; } /// /// 判断点是否在三角形内 -- 快速 /// /// /// /// /// /// private static bool IsAnyPointInsideGeometryTriangle(List 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; } /// /// a fast way to check a point in triangle /// /// /// /// /// /// 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; bool inside = (cross1.y >= 0 && cross2.y >= 0 && cross3.y >= 0) || (cross1.y 3 && ((convexIndex = SearchConvexIndex(tempPoints)) != -1)) { int index0 = convexIndex == 0 ? tempPoints.Count - 1 : convexIndex - 1; int index1 = convexIndex; int index2 = (convexIndex == tempPoints.Count - 1) ? 0 : convexIndex + 1; GenericTriangulate_Mapping(index0, index1, index2, vertices, searchMap, normal_plane, normal_triangle); tempPoints.RemoveAt(convexIndex); for(int i = convexIndex; i < tempPoints.Count; i++) { searchMap[i] = searchMap[i] + 1; } Debug.Log("被剔除顶点 : " + convexIndex); } // last 3 point GenericTriangulate_Mapping(0, 1, 2, vertices, searchMap, normal_plane, normal_triangle); int baseCount = 0; List indexesCaculator = new List(); foreach(var planeData in normal_plane) { var dir = planeData.Key; var vertIndexes = planeData.Value; List plane_triangle = null; normal_triangle.TryGetValue(dir, out plane_triangle); if(plane_triangle != null && plane_triangle.Count >= 3) { int[] indexes = new int[vertIndexes.Count]; vertIndexes.CopyTo(indexes); indexesCaculator.Clear(); indexesCaculator.AddRange(indexes); for(int i = 0; i < indexesCaculator.Count; i++) { points.Add(vertices[indexesCaculator[i]]); } for(int i = 0; i < plane_triangle.Count; i++) { int planeIndex = indexesCaculator.IndexOf(plane_triangle[i]); triangles.Add(planeIndex + baseCount); } baseCount += vertIndexes.Count; } } }

　　它的实现逻辑其实就是在原有的三角剖分基础上, 添加一个顶点查找表, 因为在计算中会把顶点列表中的点一个个剔除, 所以查找表就表达当前列表的Index与原始列表的Index的对应关系. 然后每个平面上的顶点用一个HashSet存储顶点, 平面上的三角形用List 存储Index, 最后重新组合起来即可, 这样就实现了顶点共享.

(2019.06.26)

　　在前面的基础上, 添加一个顶点列表顺序设置, 按照顶点列表是顺时针或是逆时针来进行三角形三个顶点的顺序排布, 使最后生成的三角形是朝向Y轴正方向的. 这样在应用上就给三角剖分提供了便利性.

　　这样最终版本就到这里了, 经过优化使用了共享顶点, 并添加了顶点顺序标记来生成三角剖分.

　　PS : 顶点查找表错误, 后期进行了修改(2019.08.08)

　　PS : 这些过程在模型上是最优的, 可是程序上不是, 它中间计算过程new了很多Dictionaty / List / HashSet等, 非常不友好, 后期需要全部修改

扩展(2019.06.26)

　　在获得对顺序点进行的三角剖分之后, 能不能想办法把它扩展到更加一般的情况呢? 比如在三维空间中的点, 比如立体图形的三角剖分? 我觉得是可以做到的, 只不过对顶点列表的数据来源有严格要求.

先来看这张图 :

　　这张图左边是原始图形, 它由两个三角形构成( 绿色ABD和红色BCD ), 严格来说它是4个点构成的, 如果直接用4个点对它进行剖分, 可能得到的是右图的样子, 蓝色的是剖分后获得的三角形, 变成了ABC, ACD 肯定是不对的.

可是如果A点也是在X-Y平面上的话, 比如在y轴上, 那么按照ABC, ACD的剖分也是没有问题的, 看下图 :

　　在同一平面上的话, 只要能剖分出三角形就能在表现上正确了. 所以根据这个逻辑, 在三维空间内对于要剖分的点的集合, 需要根据它们所在的平面进行划分, 然后再进行剖分即可. 逻辑的主要部分在点集所在的平面上.

　　比如上图, 它是一个垂直于X-Z平面的面, 如果使用X-Z平面的话, 这个点集的所有点都会投影在X-Z平面的一条直线上, 就无法分割出三角形了, 所以计算应该基于当前面.

　　感觉这些意义不大, 始终是要它进行空间转换到同一个平面上的, 就用上面的逻辑就行了......

(2019.07.08)

　　这个博客没有时间线管理的功能, 每次发现问题然后往之前的错误里面更新的话看起来就乱七八糟了, 希望能有更好的管理工具像SVN一样可以把博客的时间线控制起来.

今天发现了几个问题:

　　1. 顶点列表顺序添加之后, 计算查找表的时候出错了, 查找表就是用当前的顶点列表顺序去对应原始顶点列表顺序的功能, 因为每次剖出三角形之后会去掉一个顶点, 然后就需要顶点查找表了.

　　2. Vector3的哈希值很迷, 在作为容器的Key的时候很不可靠, 可以看下面的测试:

　　按照逻辑相同的值的向量它的哈希值应该是相同的, 就像float, int这些值类型, Vector3代表它的唯一性的应该是它的三个值xyz, 所以不管你怎样 new Vector3() 只要值相等, 它都是一样的, 可是叉乘之后的就有问题. 看下图v1 / v2的各个xyz是相等的, 可是v1, v2的哈希值不等, 然后用这些值创建新的Vector3, d1 / d2 还是跟v1, v2的哈希值一样......

 　　看看它们的值相等, 可是哈希不等, 所以之前使用Vector3作为容器Key的地方都不对了, 所以做了修改把Vector3转成的String来获得唯一性