手写增量式 Delaunay：凸包扇形起步、边上点拆分与 Lawson 翻边

Delaunay 三角剖分有好几种构造方式：分治、扫描线、增量插入、先建任意三角网再翻边。它的核心要求是空圆性：每个三角形的外接圆内不含其他输入点。演示页里用的是增量插入，因为它比较贴近「一步一步算出来」的直觉：先有一个合法的三角网，再把点一个个放进去，每次只修受影响的一小片。

这篇按代码的执行顺序拆开：初始三角网怎么来、点怎么定位、落在边上怎么拆、翻边栈怎么工作、最后怎么收尾。关键步骤都贴了实现节选，不是伪代码。

起点：为什么从凸包扇形开始

增量插入需要一个初始三角网。常见的做法是造一个超三角形，把全部点包住，插完再删掉与超三角形顶点相连的三角形。这个做法很常见，但它引入了三个临时顶点，也就引入了「哪些三角形是临时的、哪些是真实的」这个问题；删除阶段的判断一旦有偏差，缺口往往就出现在凸包边附近，而随机测试不容易复现。

凸包扇形起步避开了这个问题：凸包是一个凸多边形，从任意一个顶点向其余顶点连对角线，得到的扇形三角化一定是合法的。初始三角网只由真实输入点组成，没有任何临时顶点，也就不需要「删除临时结构」这个阶段。

const hullRes = convexHull(pts);
const hullIdx = hullRes.hullIndices;   // 去重后下标，逆时针
for (let k = 1; k + 1 < hullIdx.length; k++) {
  addTri(hullIdx[0], hullIdx[k], hullIdx[k + 1]);
}

代价是初始三角化依赖凸包的正确性。凸包漏一个拐点，扇形就少一块，后面的插入补不回来。所以凸包和剖分必须一起测，不能假设「凸包反正很简单，不会错」。

数据结构：三角形表 + 边索引

代码里有两个核心结构：

const T = [];                 // { v: [a, b, c], alive }
const emap = new Map();       // "a,b" -> [三角形下标]

T 存所有三角形，包括已经被删除的；alive 标记它是否还在。保留死三角形是为了让下标稳定，翻边时不需要搬数组。

emap 是边到三角形的索引。键用两个端点下标排序后拼成字符串，保证 (a,b) 和 (b,a) 落到同一个键：

const ekey = (a, b) => (a < b ? a + ',' + b : b + ',' + a);

addTri 在插入前统一方向：如果三个点算出来的叉积为负，就交换后两个顶点，保证每个三角形都是逆时针。这个约定后面很重要：外接圆行列式要求前三个点逆时针，方向反了，判断结果会整体反过来。

killTri 把三角形标记为死，并从三条边的索引里移除它；某条边的列表空了，就把这个键删掉。这样 emap 里不会留下空数组，遍历时少一层判断。

第一步：插入顺序

凸包上的拐点已经在初始扇形里了，不需要再插。剩下的点按输入顺序逐个插入：

const onHull = new Set(hullIdx);
const insertOrder = [];
for (let i = 0; i < n; i++) if (!onHull.has(i)) insertOrder.push(i);

注意 onHull 只包含凸包拐点。凸包边上的共线点不在这个集合里，它们会走插入流程，并在「点落在三角形边上」的分支里被拆开。这一点如果漏了，这些点就不会被任何三角形覆盖，面积和、覆盖面检查都会失败。

第二步：定位

对每个待插入点 p，遍历当前还活着的三角形，用三个有向面积判断它是否在内部：

const d1 = cross(all[a], all[b], p);
const d2 = cross(all[b], all[c], p);
const d3 = cross(all[c], all[a], p);
if (d1 >= -onEdgeEps && d2 >= -onEdgeEps && d3 >= -onEdgeEps) {
  host = ti;
  // 记录哪条边的 d 接近 0
}

onEdgeEps = 1e-9 * dmax，dmax 是包围盒长边。它随点集尺度放大，不是一个固定的绝对小数；固定阈值在坐标整体放大后会失效。

严格按量纲，有向面积是长度平方，容差更自然的写法是 dmax² 的量级。这里用 dmax 是工程上的启发式，量级上仍然随尺度走，够这个演示用；如果坐标跨度变化很大，应该把它改成和叉积同量纲的阈值。

定位是线性扫描：每次插入都从第一个三角形扫到最后一个。200 个点、几百个三角形时没问题；上千点就会明显变慢。工业实现会用 walk（从上一个三角形出发，沿邻接关系走过去）或者空间索引，避免每次都扫全表。

第三步：拆三角形

找到宿主之后，分两种情况。

点在三角形内部：把宿主一分为三，新点分别和三条边组成三角形。

const [a, b, c] = T[host].v;
killTri(host);
addTri(a, b, i); addTri(b, c, i); addTri(c, a, i);
stack.push([a, b], [b, c], [c, a]);

点落在边上：这条边两侧可能各有一个三角形。内部边两侧各拆成两个，一共四个；凸包边只有一侧，拆成两个。

const [a, b] = onEdge;
const t1 = host, c1 = apex(t1, a, b);
const t2 = neighborOf(a, b, t1);
killTri(t1);
if (t2 !== null) {
  const c2 = apex(t2, a, b);
  killTri(t2);
  addTri(a, i, c1); addTri(i, b, c1);
  addTri(b, i, c2); addTri(i, a, c2);
  stack.push([a, c1], [b, c1], [a, c2], [b, c2]);
} else {
  addTri(a, i, c1); addTri(i, b, c1);
  stack.push([a, c1], [b, c1]);
}

apex(t, a, b) 返回三角形 t 里既不是 a 也不是 b 的那个顶点；neighborOf(a, b, t1) 从边索引里找另一侧的三角形。

这个分支是凸包边共线点的入口。如果只按「点在三角形内部」处理，这些点会被当成内部点拆成三块：零面积的那块随后被清掉，边上点与对顶点连了线，凸包边界被折成两段，后续翻边还可能造出重叠。对拍时表现为面积和大于凸包面积、空圆性出现违规——不是「少一块」，而是「多出来一块」。

第四步：Lawson 翻边

翻边就是交换四边形的一条对角线：把不合法的边拿掉，换成另一条。新三角形加进去之后，局部可能不满足空圆条件。把新产生的边压进栈，逐条检查：

const [a, b] = stack.pop();
const arr = emap.get(ekey(a, b)) || [];
let t1 = -1, t2 = -1;
for (const ti of arr) {
  if (!T[ti].alive) continue;
  if (T[ti].v.includes(i)) t1 = ti; else t2 = ti;
}
if (t1 < 0 || t2 < 0) continue;
const c = apex(t1, a, b);
const q = apex(t2, a, b);

t1 是包含新点 i 的一侧，t2 是另一侧。如果 q 落在三角形 (a,b,c) 的外接圆内，这条边就不满足局部 Delaunay 条件，把它翻成 (c,q)。

判断用外接圆行列式：

const ax = a.x - p.x, ay = a.y - p.y;
const bx = b.x - p.x, by = b.y - p.y;
const cx = c.x - p.x, cy = c.y - p.y;
const det =
  (ax * ax + ay * ay) * (bx * cy - cx * by) -
  (bx * bx + by * by) * (ax * cy - cx * ay) +
  (cx * cx + cy * cy) * (ax * by - bx * ay);
return det > tol;

这里有两个容易写错的地方。

第一，前三个点必须逆时针。建三角形时统一过方向，但栈里存的 (a,b) 来自边索引的键，顺序未必和三角形边序一致，所以不能假设 (a,b,c) 已经是逆时针。代码先归一化：

const e0 = cross(all[a], all[b], all[c]) < 0 ? [b, a] : [a, b];
if (!inCircumcircle(all[e0[0]], all[e0[1]], all[c], all[q], flipTol)) continue;

第二，容差要跟着行列式的次数走。这个行列式对坐标是四次的，所以 flipTol = 1e-12 * dmax^4；用 dmax 或 dmax² 都不对。容差太小，共圆点会在边界上来回翻；太大，真正该翻的边翻不动。

翻边替换了两个三角形，四条外边界上的邻接关系变了，继续压栈复查：

killTri(t1); killTri(t2);
addTri(a, q, c); addTri(q, b, c);
stack.push([a, q], [q, b], [b, c], [c, a]);

新对角线 (c,q) 本身不压栈：它是按「对面点在外接圆内」这个条件翻出来的，翻完在这对三角形之间已经合法；需要复查的是它周围那四条边。

栈空时，局部合法化结束。

第五步：全局合法化与收尾

局部翻边只处理新点周围的边，而且代码里也只对「包含新点的一侧」做判断。凸包扇形作为初始网是合法的三角化——三角形互不重叠、刚好覆盖凸包——但一般不是 Delaunay 的；插入过程中，远处的边如果本来就不合法，局部过程也不会碰它。所以代码在全部点插完之后再扫一遍所有内部边：

for (let round = 0; round < 60; round++) {
  let flipped = 0;
  for (const [k, arr] of Array.from(emap.entries())) {
    if (arr.length !== 2) continue;
    // 两侧顶点 c、q；若 q 在 (a,b,c) 外接圆内则翻边
  }
  if (!flipped) break;
}

这一步不是「为了稳妥」的保险，而是正确性所必需。我做过对照：把这段循环整个关掉，凸五边形会出现 6 处空圆违规，一条 40 点近共线算例出现 19 处；恢复后都归零。

上限 60 轮是防御性的：共圆点在容差边界上可能来回翻，没有上限就会卡死。正常情况下，一轮或两轮就收敛了。

最后，把面积接近零的三角形剔掉：

if (Math.abs(cross(all[a], all[b], all[c])) <= 1e-12 * dmax * dmax) continue;

这些零面积三角形通常来自凸包边上的共线点：三个点共线，几何上不构成面。留在结果里会让三角形数、最小角、欧拉关系全部失真。

复杂度与边界

• 定位：每次插入都从表头扫到表尾，而且扫的是 T.length，包含已经死掉的三角形；翻边越多，表越长。最坏 O(n²)，实际常数也比「只扫活三角形」大；
• 翻边：代码按输入顺序插入，没有随机化，所以只能谈常见输入下的经验表现，没有「期望接近常数」的理论保证；全局合法化还要多扫几轮；
• 空间：T 保留死三角形让下标稳定，规模与累计翻边次数同阶，不只是 O(n)；emap 同理。

这套实现适合演示和几百个点的交互，不适合上万点的离线剖分。要上规模，得换定位结构、减少全局扫描，或者直接上分治。

已知边界：

• 共圆点集的剖分本身不唯一，对拍这类算例（网格、正圆、正方形）只比三角形数量和面积；其余算例比三角形集合；
• 不支持带权 Delaunay，也不做约束 Delaunay。

验证

验证分三层。

一是内核断言：40 项，覆盖退化输入、外接圆心、面积和、欧拉关系等。不依赖浏览器，改完代码几秒钟能跑完。

二是外部对拍——拿一套独立实现跑同一批输入，逐项比对。固定 11 个算例、91 条判据：凸包 10/10 一致；三角形集合 7/7 完全相同；面积和的相对误差在 1e-16 量级。固定算例的价值是可复现：同一个输入，谁跑都是同一组判据。

三是页面实测：四个页签、拖点、200 点压力下没有控制台错误，页内 14 项自检全部通过。

三层各管一段：内核断言管边界条件，外部对拍管「两套实现不会一起错」，页面实测管交互路径。单靠任何一层都不够。

几条实现上的经验

1. 从一个合法但不必是 Delaunay 的结构开始。凸包扇形保证互不重叠、刚好覆盖凸包，插入和翻边再把它逐步修成 Delaunay；这样不需要临时顶点，但收尾的全局合法化不能省。
2. 方向统一放在建三角形的时候。每个三角形建的时候就保证逆时针，后面的外接圆判断、面积计算、翻边都不用再考虑方向。
3. 容差按量纲分组。外接圆行列式是四次量，零面积判断是二次；不要用一个固定 epsilon 套所有地方。点在三角形内的判定这里用了 dmax 的启发式，和严格量纲的差别前面已经说明过。
4. 退化输入要能返回状态，不要抛异常。全共线、点数不足都是合法输入，返回 degenerate 标记，调用方才好决定是画一条线还是跳过剖分；重复点则在入口直接去重并计数，不进入退化分支。
5. 外部对拍要固定算例。同一批点、同一组判据，每次改完都能复现；只靠随机点和肉眼，低概率错误很难被抓到。

（完）