★ FEATURED ARTICLE 简介这是一套基于C的增量德劳内三角剖分算法实现面向计算几何、计算机图形学及网格生成方向的开发者与学生。算法通过最大化三角剖分中所有三角形的最小角来生成较高质量的网格项目采用石溪大学顾贤峰教授提出的半边数据结构组织网格拓扑。主要源码集中在配套的示例工程中并包含MeshLib基础库方便读者理解算法细节并扩展到曲率计算等应用。资源压缩包共89个文件既有核心的源文件和头文件也有Visual Studio工程配置、可执行程序、调试日志等辅助文件整体大小约4.17MB目录结构清晰便于快速导入和重新编译。目前已有1001人学习下载适合正在学习计算几何或需要实现网格剖分的读者参考。通过运行演示和阅读代码读者可以掌握增量插入点的完整流程、半边结构的构建方法以及在实际工程中调试几何算法的技巧。1. 增量 Delaunay 三角剖分为什么值得用 C 自己实现Delaunay 三角剖分是计算几何里的老问题给定平面上一组点构造三角形网格使得任意三角形的外接圆内不包含其他点。增量算法逐点插入是其中最直观、最适合工程落地的一条路线——每插入一个新点只影响局部若干三角形通过边翻转恢复 Delaunay 性质复杂度在随机点集下可达到 O(n log n) 的期望水平。对做网格生成、路径规划、地形建模、插值渲染的从业者来说自己用 C 实现一遍增量 Delaunay能彻底搞懂剖分过程的每一步而不是把第三方库当黑匣子。这篇文章从几何判定写到完整代码把结构设计、参数取舍和踩过的坑一并讲透。2. 剖分原理与数据结构从几何判据到 C 选型2.1 外接圆测试Delaunay 性质的唯一判据要写增量 Delaunay第一个绕不开的操作是外接圆测试业内通常叫 inCircle 测试。给定三个点 a、b、c 构成一个三角形再给一个待判断的点 d如果 d 落在三角形 abc 的外接圆内部那么 d 与 abc 组成的四边形就不满足 Delaunay 性质需要翻转对角线。外接圆测试在实现上有两种常见做法一种是用几何坐标直接求解圆心和半径然后比较距离另一种是构造行列式做符号判断。第一种做法看起来直观但隐藏着数值隐患计算圆心需要解线性方程组涉及浮点除法分子分母都很接近时误差会明显放大。第二种做法只需要做四次乘法和若干次减法最后判断一个行列式的符号中间没有除法数值稳定性好得多。实际工程里我几乎只用行列式版本配合 double 精度在常规坐标范围内完全够用。行列式版本长这样struct Point { double x, y; }; // 返回 det 的符号det 0 表示 d 在 abc 的外接圆内 // 前提a, b, c 按逆时针顺序排列 double inCircleDet(const Point a, const Point b, const Point c, const Point d) { double adx a.x - d.x, ady a.y - d.y; double bdx b.x - d.x, bdy b.y - d.y; double cdx c.x - d.x, cdy c.y - d.y; double ad2 adx * adx ady * ady; double bd2 bdx * bdx bdy * bdy; double cd2 cdx * cdx cdy * cdy; return adx * (bdy * cd2 - cdy * bd2) - bdx * (ady * cd2 - cdy * ad2) cdx * (ady * bd2 - bdy * ad2); }逻辑说明把四个点都平移到以 d 为原点然后构造一个 3×3 行列式前三列是 a、b、c 相对 d 的向量第三列是向量长度的平方。行列式的符号直接给出内外关系不需要计算圆心坐标。参数和边界条件当 det 恰好等于 0 时表示四点共圆。这种情况在严格意义上是退化情形工程上我一般把 0 判为不在圆内避免在共圆时做无意义的翻转。另外逆时针顺序这个前提非常关键——如果三角形 abc 是顺时针行列式符号会反转原本在圆内的点会判成在圆外。所以我在调用 inCircleDet 之前都用 orient2d 检查一遍方向必要时交换顶点。这也是后续所有写三角形、翻转边的代码都要遵循的统一约定所有三角形顶点一律逆时针存储。2.2 三角形表与邻接索引十万点级别的最优选型增量 Delaunay 的核心操作是找包含新点的三角形、拆三角形、然后沿边翻转。这个过程中需要频繁查询三角形的邻接关系、遍历边、修改拓扑。数据结构直接决定实现的复杂度和运行速度。常见的选择有三种。第一种是直接用三角形表每个三角形存三个顶点索引再附加三个邻接三角形索引。这种结构简单、内存紧凑处理十万级别点集时性能完全足够缺点是删除和插入三角形时维护邻接关系要仔细。第二种是半边结构Half-Edge在有限元网格和三维网格处理里更通用但实现量大做二维增量剖分有点像杀鸡用牛刀。第三种是只维护边到三角形的哈希映射适合需要频繁按边查询的场景但多一层哈希开销。对于这里的目标我推荐三角形表加邻接索引。原因有三第一增量算法的主要操作是局部替换三角形表做局部替换最直观第二内存连续缓存友好第三最后输出网格时往往就要三角形列表省去一次格式转换。下面给出简洁定义struct Triangle { int v[3]; // 三个顶点索引逆时针排列 int t[3]; // 邻接三角形索引-1 表示边界边 bool removed; // 延迟删除标记 }; std::vectorPoint points; std::vectorTriangle tris;邻接索引的约定是t[i] 记录与顶点 v[i] 相对的邻接三角形也就是共享边 (v[(i1)%3], v[(i2)%3]) 的那个三角形。这个约定让翻转操作可以直接通过 t[edgeIdx] 拿到对面的三角形不需要遍历查找。用 -1 表示边界初始化超三角形时把三个邻接都设为 -1。这个数据结构的直接后果是删除三角形时不真正从 vector 里 erase而是打上 removed 标记。这样能避免频繁搬移内存也避免索引失效。代价是 vector 里会堆积一些死三角形当死三角形数量超过总量一半时做一次紧凑化压缩即可。这个策略在绝大多数网格生成库里都能看到属于工程上的基本操作。3. 增量插入实现超三角形、定位、翻转与清理3.1 超三角形初始化与步行定位增量算法需要一个初始包围结构。常见做法是构造一个足够大的超三角形把输入的所有点都包含在内。超三角形的三个顶点是虚拟的不属于最终网格算法结束时删除所有与超三角形顶点相连的三角形。超三角形怎么选我见过有人取点集包围盒的三个角也有人用重心加半径构造一个等边三角形。我习惯用后者先算出点集的包围盒取中心 c取半径 r max(xmax - xmin, ymax - ymin) 的 20 倍然后构造等边三角形的三个顶点。为什么要放大这么多我踩过的坑是如果超三角形不够大某些插入点会落在超三角形某条边的小邻域外步行法定位时会直接走出网格边界更隐蔽的是边界三角形的外接圆仍然可能包含超三角形之外的区域在翻转时把虚拟边搞乱。放大 20 倍后这些问题几乎不会出现代价只是几个浮点乘法。void initSuperTriangle(const std::vectorPoint pts, double margin 20.0) { double xmin pts[0].x, xmax pts[0].x; double ymin pts[0].y, ymax pts[0].y; for (const auto p : pts) { xmin std::min(xmin, p.x); xmax std::max(xmax, p.x); ymin std::min(ymin, p.y); ymax std::max(ymax, p.y); } double cx (xmin xmax) * 0.5; double cy (ymin ymax) * 0.5; double r std::max(xmax - xmin, ymax - ymin) * margin; // 等边三角形三个顶点逆时针排列 Point p0 {cx, cy 2.0 * r}; Point p1 {cx - std::sqrt(3.0) * r, cy - r}; Point p2 {cx std::sqrt(3.0) * r, cy - r}; int i0 (int)points.size(); points.push_back(p0); points.push_back(p1); points.push_back(p2); tris.clear(); tris.push_back(Triangle{{i0, i1, i2}, {-1, -1, -1}, false}); }参数说明margin 是放大系数默认 20 倍。对于坐标范围在 [-1e6, 1e6] 以内的数据20 倍足够让超三角形外接圆覆盖所有可能的插入点如果你处理的是经纬度坐标或者数值范围极大/极小的数据要适当调大这个系数。定位三角形我用步行法。思路是从当前任意一个三角形出发沿着目标点所在方向逐个三角形跳过去。对于随机顺序插入的点步行法的期望步数是常数整体复杂度依然是 O(n log n) 级别。// 返回包含 p 的三角形索引找不到返回 -1 int locateTriangle(const Point p, int startTri) { int cur startTri; int guard 0; int maxGuard (int)tris.size() * 4 1024; while (cur 0 guard maxGuard) { const Triangle tri tris[cur]; if (tri.removed) { cur -1; break; } // 逆时针三角形内部在三条边的左侧 if (orient2d(points[tri.v[0]], points[tri.v[1]], p) 0) { cur tri.t[0]; } else if (orient2d(points[tri.v[1]], points[tri.v[2]], p) 0) { cur tri.t[1]; } else if (orient2d(points[tri.v[2]], points[tri.v[0]], p) 0) { cur tri.t[2]; } else { return cur; // 点在三角形内部或边上 } } return -1; }逻辑说明orient2d(a, b, p) 0 表示点 p 在向量 a→b 的右侧。因为三角形内部在每条边的左侧一旦 p 在某边的右侧说明 p 越过了这条边应该跳到这条边对面的邻接三角形。三次判断都不越过任何边说明 p 就落在当前三角形内部。这里 guard 保护上限是必要的。三角形数量是动态增长的所以用当前 tris.size() * 4 再加一个常量的形式。返回 -1 时不要硬着头皮继续打印出当前坐标和 cur 三角形的顶点坐标基本一眼能看出是超三角形太小还是数据结构坏了。我把 -1 视为可恢复错误在工程里我会重新扩大超三角形并重建整个剖分而不是尝试修复。3.2 拆分三角形与邻接修复定位到目标三角形后把新点插进去。插入操作的本质是把旧三角形拆成三个新三角形每个新三角形都以新点为公共顶点。关键点在于邻接关系的修复——旧三角形原来的三个邻接三角形它们的邻接索引还指向旧三角形必须全部改成指向对应位置的新三角形。int insertPoint(int pIdx) { int loc locateTriangle(points[pIdx], 0); if (loc 0) return -1; Triangle old tris[loc]; if (old.removed) return -1; int a old.v[0], b old.v[1], c old.v[2]; int ta old.t[0], tb old.t[1], tc old.t[2]; int i0 (int)tris.size(); int i1 i0 1; int i2 i0 2; // 新三角形都保持逆时针顶点顺序 tris.push_back(Triangle{{pIdx, a, b}, {i2, i1, ta}, false}); tris.push_back(Triangle{{pIdx, b, c}, {i0, i2, tb}, false}); tris.push_back(Triangle{{pIdx, c, a}, {i1, i0, tc}, false}); old.removed true; // 修复三个外部邻接把指向 loc 的索引改成对应新三角形 if (ta 0) setNeighbor(ta, loc, i0); if (tb 0) setNeighbor(tb, loc, i1); if (tc 0) setNeighbor(tc, loc, i2); // 对新产生的三条边做翻转检查 flipCheck(i0, 0); flipCheck(i1, 0); flipCheck(i2, 0); return i0; }setNeighbor 的实现是在旧邻接三角形的三个 t 索引里找到原来指向 loc 的那一个替换成新三角形索引。这一步看起来很机械但出了 bug 往往很隐蔽因为邻接关系是双向的——你改了这头那头也必须同步改。void setNeighbor(int triIdx, int oldNeighbor, int newNeighbor) { if (triIdx 0) return; Triangle t tris[triIdx]; for (int k 0; k 3; k) { if (t.t[k] oldNeighbor) { t.t[k] newNeighbor; return; } } }拆完三角形后网格局部出现了三条新边和三个新三角形。但这些新三角形不一定满足 Delaunay 性质需要用边翻转把网格修回合法状态。3.3 边翻转恢复 Delaunay 性质边翻转是增量 Delaunay 的灵魂。具体来说如果两个相邻三角形共享一条边把其中一个三角形的第三个顶点放进另一个三角形的外接圆里测试发现它在圆内那么这条公共边必须翻转——删掉公共边换成连接两个第三顶点的新边。我倾向于把翻转逻辑拆成两个函数flipCheck 负责判断是否需要翻转flipEdge 负责执行翻转。判断的精度直接影响网格质量执行的正确性直接影响拓扑完整性。flipCheck 的实现如下// 检查三角形 idx 的第 edgeIdx 条边是否需要翻转 void flipCheck(int idx, int edgeIdx) { Triangle tri tris[idx]; if (tri.removed) return; int oppIdx tri.t[edgeIdx]; if (oppIdx 0) return; Triangle opp tris[oppIdx]; if (opp.removed) return; // 找到 opp 中不属于公共边的顶点 int oppVertex -1; for (int i 0; i 3; i) { if (opp.t[i] idx) { oppVertex opp.v[i]; break; } } if (oppVertex 0) return; // 邻接关系不一致数据结构已损坏 const Point p1 points[tri.v[(edgeIdx 1) % 3]]; const Point p2 points[tri.v[(edgeIdx 2) % 3]]; const Point pThis points[tri.v[edgeIdx]]; const Point pOpp points[oppVertex]; // tri 的顶点是逆时针的inCircleDet 正号表示 pOpp 在圆内 if (inCircleDet(p1, p2, pThis, pOpp) 0) { flipEdge(idx, oppIdx); } }flipEdge 的核心是重建两个三角形的顶点和邻接关系。这一步最容易写乱我踩过的坑包括翻完边后忘了修外部邻接、新三角形的顶点顺序变成顺时针、递归检查时没有防止死循环。为了不让自己在顶点顺序上反复推导我在工程里用了一个小工具函数每次构造完三角形后检查 orient2d方向不对就交换 v[1] 和 v[2]。void makeCCW(Triangle t) { if (orient2d(points[t.v[0]], points[t.v[1]], points[t.v[2]]) 0) { std::swap(t.v[1], t.v[2]); } }有了 makeCCW 兜底翻转函数可以专注于拓扑关系。整体流程是找公共边、找各自非公共顶点、构造两个新三角形、修复外部邻接、递归检查新形成的边。这个函数完整写下来大约需要六七十行我给出骨架和关键顺序其余按语义补齐即可void flipEdge(int idx, int oppIdx) { Triangle tri tris[idx]; Triangle opp tris[oppIdx]; // 第一步收集四个顶点确认公共边 int common[2], cnt 0; for (int i 0; i 3; i) { for (int j 0; j 3; j) { if (tri.v[i] opp.v[j]) common[cnt] tri.v[i]; } } if (cnt ! 2) return; // 理论上不会发生 // 第二步找各自非公共顶点 int p -1, q -1; for (int i 0; i 3; i) { if (tri.v[i] ! common[0] tri.v[i] ! common[1]) p tri.v[i]; if (opp.v[i] ! common[0] opp.v[i] ! common[1]) q opp.v[i]; } // 第三步保存两边各自的外部邻接 int edgeIdx -1; for (int i 0; i 3; i) { int v1 tri.v[i], v2 tri.v[(i 1) % 3]; if ((v1 common[0] v2 common[1]) || (v1 common[1] v2 common[0])) { edgeIdx i; break; } } int extA tri.t[(edgeIdx 1) % 3]; int extB tri.t[(edgeIdx 2) % 3]; int oppEdgeIdx -1; for (int i 0; i 3; i) { if (opp.t[i] idx) { oppEdgeIdx i; break; } } int extC opp.t[(oppEdgeIdx 1) % 3]; int extD opp.t[(oppEdgeIdx 2) % 3]; // 第四步构造新三角形暂时不关心邻接用 makeCCW 保证方向 tris[idx] {{common[0], p, q}, {-1, -1, -1}, false}; tris[oppIdx] {{common[1], q, p}, {-1, -1, -1}, false}; makeCCW(tris[idx]); makeCCW(tris[oppIdx]); // 第五步重新设置邻接 setPairNeighbor(idx, oppIdx); setPairNeighbor(idx, extA); setPairNeighbor(idx, extB); setPairNeighbor(oppIdx, extC); setPairNeighbor(oppIdx, extD); // 第六步对 idx 和 oppIdx 各自的两条新边继续检查 flipCheck(idx, 0); flipCheck(idx, 1); flipCheck(idx, 2); }setPairNeighbor 做的事情是在 triA 的三个邻接槽位里找到指向 triB 的那一个设为新的值。这个函数和 setNeighbor 类似只是它同时保证双向关系最终由下一次调用建立。这里我简化了细节实际工程中我建议把邻接关系修复写成一个统一的函数每次修改都同步两边。递归检查这里有一个风险两个相邻三角形可能在翻转后又被翻回来形成死循环。我的处理方式是在递归函数里维护一个翻转计数器超过三角形容量的 10 倍就报错并终止。在已经正确的实现里这个计数器永远不会触发但它能在数据结构被意外破坏时给出快速反馈。3.4 移除超三角形边界与输出所有原始点插入完成后网格中还残留与超三角形三个虚拟顶点相连的三角形。这些三角形需要剔除否则输出网格会包含三个极远处的虚拟点。清理逻辑很简单遍历所有三角形只要三个顶点索引里有任何一个大于等于原始点数 n即属于超三角形顶点就打上 removed 标记。最后再做一次紧凑化把活着的三角形搬到 vector 前部顺便重映射邻接索引。void removeSuperTriangle(int originalPointCount) { for (auto tri : tris) { if (tri.removed) continue; if (tri.v[0] originalPointCount || tri.v[1] originalPointCount || tri.v[2] originalPointCount) { tri.removed true; } } // 紧凑化把未删除的三角形前移并重映射邻接索引 int writeIdx 0; for (int i 0; i (int)tris.size(); i) { if (tris[i].removed) continue; tris[writeIdx] tris[i]; } tris.resize(writeIdx); // 重映射邻接索引 std::vectorint remap((int)tris.size(), -1); // 具体实现需要记录旧索引到新索引的映射 // 这里省略细节但要强调邻接索引必须同步更新否则输出结果一查邻接就崩 }清理完成后三角形表里存的就是最终剖分。输出格式可以是简单的三个顶点索引也可以带上邻接关系供后续有限元求解或渲染使用。我一般会额外生成一份边列表方便可视化验证。4. 避坑与排查浮点、退化、边界与调试4.1 浮点误差导致的外接圆误判现象剖分结果整体看起来正常但个别三角形特别狭长或者明明两点距离很近却生成了跨越整个区域的边。用可视化工具放大后能看到一些三角形的顶点几乎共线。原因外接圆测试的行列式结果接近零时double 精度不足以区分在圆内和在圆上。尤其在坐标值非常大比如经纬度的 1e7 量级或非常小1e-6 量级时平方项会放大或缩小数值差异导致符号翻转。还有另一种情况是输入点中有三点接近共线此时的三角形外接圆半径极大测试结果对微小扰动极其敏感。解决第一输入点先做归一化或平移把坐标中心移到原点附近范围缩放到大致 [-1, 1] 区间这个操作能显著改善数值稳定性。第二在 inCircleDet 返回值的绝对值小于某个阈值比如 1e-12时按不在圆内处理不做翻转。阈值怎么定我的经验是先跑一遍剖分打印所有 det 的绝对值分布取一个比最小非零值小几个数量级的数作为阈值。// 带 epsilon 的版本调用时传入相对容差 bool inCircleEps(const Point a, const Point b, const Point c, const Point d, double eps) { double det inCircleDet(a, b, c, d); return det eps; // det 绝对值低于 eps 视为共圆 }4.2 共线与重复点最容易被忽略的输入陷阱现象程序在插入某个点后陷入死循环或者翻转检查无限递归最终栈溢出。去掉几个点后又恢复正常。原因输入点集中存在重复坐标或三点共线。重复点会导致定位三角形的外接圆测试永远返回 0翻转逻辑在两个三角形之间来回切换。三点共线更隐蔽共线的三个点不能构成有效的 Delaunay 三角形但算法在插入过程中会产生这种退化三角形如果不特殊处理它会一直参与翻转检查。解决在算法开始前先对输入做预处理。去重的标准做法是按坐标排序相邻点距离小于阈值则合并。共线的处理有两种思路一种是在三角形构造函数里加检查发现三点共线就跳过另一种是微扰法给每个点加上一个极小的随机偏移破坏共线性。我建议在工程里用第一种因为微扰会改变原始数据对需要精确保留坐标的应用不可接受。// 简单的分组去重按格子哈希格子大小为 eps std::vectorPoint deduplicate(const std::vectorPoint pts, double eps) { std::vectorPoint result; std::unordered_maplong long, bool seen; for (const auto p : pts) { long long key ((long long)(p.x / eps) 32) | (long long)(p.y / eps); if (!seen[key]) { seen[key] true; result.push_back(p); } } return result; }这个去重方法的思路是把点按 eps 大小的格子分组同一个格子里的点只保留第一个。注意格子哈希有边界效应——两个点恰好在格子边界两侧时距离很小却不会去重。更严谨的做法是检查相邻格子的点但工程上这个误差通常可接受。4.3 超三角形尺寸与边界拓扑泄漏现象剖分结果里出现了连接原始点与超三角形顶点的长边而且这些长边不在点集凸包上。偶尔还会出现三角形重叠或者某些原始点没有出现在任何三角形里。原因超三角形选小了。当插入点靠近超三角形边界时定位可能失败即使定位成功该点的外接圆可能超出超三角形范围翻转操作会把虚拟边当成真实边处理造成边界区域拓扑混乱。解决把超三角形半径放大。我之前提到的 20 倍系数在实际使用中很少出问题但如果点集分布极度不均匀——比如大部分点集中在一个小区域少量点散布在极远处——就需要把系数上调到 100 甚至更大。另一个补救办法是插入结束后检查所有三角形的外接圆如果某个外接圆覆盖了凸包边界之外的大片区域说明超三角形不够大需要重新计算。还有一个细节构造超三角形后三个虚拟顶点的坐标不要用整数类型存直接用 double。某些实现为了省内存用 float放大 20 倍后 float 的精度在大坐标下严重不足反而引入新的数值问题。4.4 可视化调试与不变量断言现象代码写完了但剖分结果时好时坏肉眼观察网格看不出规律。原因缺少系统的验证手段。Delaunay 三角剖分有一个很好的自检性质任意三角形外接圆内不能包含其他顶点。这个性质可以写成 O(n) 或 O(n log n) 的校验函数每次跑完数据都验证一遍。另一条是拓扑不变量欧拉公式 V - E T 1 C其中 V 是顶点数E 是边数T 是三角形数C 是连通分量数。对于平面单连通区域C 1。我的调试习惯是先跑随机点集数量从 10 到 1e4 递增每跑完一个规模都用校验函数验证不通过就打印出具体出错的三角形和邻接关系。单独看索引很难定位问题我会把出错三角形的三个顶点坐标打印出来再用 MATLAB 脚本或者直接输出为 SVG 文件可视化基本一眼就能看出是翻转逻辑还是邻接修复的问题。// 校验函数检查所有三角形邻接关系是否一致 bool validateTopology() { std::vectorint incidentCount(tris.size(), 0); for (int i 0; i (int)tris.size(); i) { if (tris[i].removed) continue; for (int k 0; k 3; k) { int n tris[i].t[k]; if (n 0) continue; // 邻接三角形的邻接索引里必须能找到 i bool found false; for (int j 0; j 3; j) { if (tris[n].t[j] i) { found true; break; } } if (!found) return false; } } return true; }这个校验函数应该在每次插入后都调用直到确认逻辑稳定后再关掉。保留它在 release 版本里成本也不高十万三角形大约几十毫秒对于非实时应用完全可以接受。5. 收尾验证正确性与性能优化的具体技巧5.1 用随机点集和欧拉公式验证网格拓扑做完一个版本后我习惯写一个独立的验证程序输入是随机生成的点集输出是检测报告。验证分三层第一层检查拓扑完整性每个三角形的邻接索引双向一致边界数为凸包边数第二层检查 Delaunay 性质任意三角形的外接圆内不包含其他点第三层用欧拉公式做全局校验V - E T 1。欧拉公式的代码很直接V 是原始点数T 是剖分后的三角形数量E 需要统计所有非边界边的数量。边可以从三角形表里提出来去重也可以用邻接关系统计每条内部边被两个三角形共享边界边只被一个三角形共享所以 E (3T 边界边数) / 2。代进欧拉公式就能验证网格有没有在拓扑层面出错。这个验证看起来简单但在开发阶段帮我抓到了不止一次的翻转 bug。有一次我对一个 5000 点随机集跑剖分Delaunay 性质校验全部通过但欧拉公式差了一个数。排查后发现是 removeSuperTriangle 之后邻接索引没有重映射边界边的计数被污染了。如果不做这个全局校验这个 bug 几乎不可能用肉眼发现。5.2 两点性能优化内存布局与减少无效翻转十万点级别的剖分在默认实现下大约需要几百毫秒对大多数场景够用。但如果数据量上到百万级就需要关注性能了。我做过两个回报很高的优化第一是预分配内存在读取完所有点之后按照估计的三角形数量 reserve 三倍的内存避免 vector 反复扩容。第二是减少翻转检查的次数。减少翻转的核心观察是新插入点只影响它周围的一个小区域离它远的三角形完全不需要重新检查。我在实现里维护一个待检查栈每次插入后只压入受影响的三角形索引而不是对整个网格做全量检查。对于随机点集这个策略能把翻转相关的开销降低一个数量级。还有一种优化是把步行法的起点设为上一个插入点所在的位置。这个优化基于一个经验观察如果输入点本身具有空间局部性比如按扫描顺序采样上一个点和下一个点往往距离很近步行从上次位置出发只需几步就能到达目标三角形。配合这个优化整体性能还能再提升 20% 到 30%。// 插入主循环利用空间局部性减少定位开销 int hint 0; for (int i 0; i n; i) { int h hint; hint insertPointWithHint(points, i, h); if (hint 0) { // 记录错误信息必要时重建 } }我踩过的最后一个坑和编译器有关在 Release 开启 -O2 后inCircleDet 里的浮点操作顺序可能被编译器改变导致同一份代码在 Debug 下正确、Release 下出现个别翻转遗漏。这个问题很玄学解决方法是把 inCircleDet 标记为不优化或者用 volatile 关键字把中间变量固定下来。虽然现在的新编译器和标准浮点模式很少再触发这个问题但我还是习惯在数值敏感的代码文件顶部写上统一的编译选项说明防止团队里的其他人踩到。这个项目就像很多计算几何算法一样——理论和实现之间隔着无数细节。把每一步都想清楚再落代码比写好再回头调试要快得多。我在做第一个版本时吃了不少亏现在再去写类似的东西一定会先花半小时把数据结构约定和翻转语义写清楚再动手写第一行 C。希望这篇笔记能帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
