Voronoi 图的两条构造路径：69270 个单元逐一对拍

一个几何结构只写一遍，很难知道自己错在哪。页面能画出来、看起来没有洞，不等于数字是对的。Voronoi 图尤其如此：单元边界差几个像素，肉眼几乎看不出来，但它会一路影响下游的邻接关系、面积统计和插值结果。

所以这个项目里 Voronoi 写了两遍，最后逐单元比面积。需要先说清楚独立到什么程度：两条路径的输入不同、依赖不同，一条直接吃点集，另一条必须先把三角剖分算对；但底层的半平面裁剪原语是同一个函数。这意味着对拍能覆盖「选哪些站点来切、从哪里起切」，覆盖不到这个原语本身——那部分只能靠外部成熟实现对拍补。下面是对拍的设计、抓到的错，以及几条可以迁移的验证思路。

先把对偶关系摆清楚

给定一组站点，Delaunay 三角剖分和 Voronoi 图是一对对偶结构：

• Delaunay 三角形的外接圆心，是 Voronoi 图的一个顶点；
• 两个站点在 Delaunay 里相邻（共享一条边），它们的垂直平分线上的一段就是 Voronoi 边；凸包边对应射线，共圆时还会退化成多条边交于一点；
• 凸包边界上的点（包括共线边上的中间点）对应无界单元，它的边界会一直延伸到无穷远。

对偶关系给了两条依赖不同的构造路径：一条只看站点集合本身，一条必须先把三角剖分算对。它们的输入不同、中间结构不同，正好拿来互相检查。

路径一：半平面裁剪

最直接的构造方式：每个站点对应一个单元，这个单元等于「所有到该站点比到其他站点更近的点」的集合。每个「比其他站点更近」的条件，就是一个半平面，边界是两站点的垂直平分线。

于是算法可以写成：

1. 从画布矩形出发（它同时承担了「把无穷远截断」的职责）；
2. 对每个其他站点，用垂直平分线切一刀，保留靠近当前站点的一侧；
3. 全部切完，剩下的多边形就是这个站点的 Voronoi 单元。

多边形裁剪用 Sutherland–Hodgman：逐条边遍历，按顶点在直线哪一侧决定保留、丢弃还是插入交点。实现只有几十行，代价却不低：n 个站点，每个站点要切 n−1 刀，裁剪次数是 O(n²)；每刀还要扫描当前多边形的所有边，最坏情况下总代价到 O(n³)。实际点集里单元边数远小于 n，两百个点跑起来没问题，上千点就会明显变慢。

这条路径的优点是不依赖任何别的结构。它不关心 Delaunay 算得对不对，甚至可以先于三角剖分存在。

路径二：只沿 Delaunay 邻居切

第二条路径利用对偶关系：一个站点的 Voronoi 单元，只由它在 Delaunay 里的邻居决定。理由是单元的每一条边都对应一条 Delaunay 边——非邻居站点对应的平分线不会穿过这个单元的内部，可以完全跳过。这个结论依赖三角剖分本身正确：漏一条边、错一条边，单元就会少切一刀，面积立刻不一致。

于是：

1. 从 Delaunay 三角剖分取出每个站点的邻居表；
2. 从一个放大若干倍的包围盒出发（放大是为了容纳无界单元，不必真的算到无穷远）；
3. 只对邻居站点做垂直平分线裁剪；
4. 最后裁回画布矩形，这样才和路径一可比。

这条路径的输入是三角剖分的结果，所以它顺带检验了一件事：剖分有没有漏掉邻居。

两条路径的复杂度不同、依赖不同、连起始多边形都不同。它们算出的单元面积必须一致，这就是对拍的基础。

教科书式的对偶构造，在近共线点集上塌了

第二条路径最初不是这么写的。更「像教科书」的写法是：从 Delaunay 三角形取外接圆心，作为 Voronoi 顶点；对无界单元，沿垂直平分线方向从外接圆心向外延长一条射线，再想办法闭合多边形。

一般位置的点集上，这个写法完全没问题。但近共线的点集会把它击穿：三个近乎共线的点，外接圆心可以飞到十万像素以外（分母是接近零的叉积）；靠这些远点连弦闭合多边形，弦的位置对输入极其敏感，稍微偏一点就从站点内侧穿过去，把单元切掉一大半。

这个错误在页面上看不出来。页面渲染的是半平面裁剪路径算出的单元，对偶路径只参与数值比对、不单独画出来，所以图上不会显示任何异常。它是靠逐单元比面积暴露的：压力对照跑了 5 种点集、11 种点数、每种 30 个种子，共 69270 个单元；修复前的记录里，794 个单元与裁剪路径不一致，其中 787 个集中在近共线预设。

改成「只用邻居的平分线去裁一个放大的包围盒」之后，同一批算例的不一致数量归零。原因很朴素：半平面裁剪从头到尾只做「直线切多边形」，没有接近零的叉积做分母，也没有远点闭合。这次之后，看到需要靠很远的点来闭合的几何构造，我会先打个问号。

对拍脚本怎么定判据

上面反复提到的「不一致」，判据是这么定的。对拍最容易犯的错，是拿平均误差或者「看起来一样」当标准。这里用的是逐单元、相对误差、固定阈值：

const rel = Math.abs(clip - dual) / Math.max(clip, 1e-9);
if (rel > 1e-6) mismatches.push({ site, rel, clip, dual });

几个细节值得说明：

• 逐单元比，不做汇总。总面积相等可能只是正负误差互相抵消；逐个单元比，才能定位到具体是哪个站点、哪种点集出的问题。
• 用相对误差，不用绝对误差。坐标跨度变化时，绝对误差没有可比性。分母加一个小常数，是为了避免面积趋近于零时除零。
• 阈值要写明。1e-6 是这套实现的选择，换一套浮点策略可能需要调整；写清楚阈值，别人才知道「通过」是在什么标准下通过。
• 固定种子。同一批种子跑出来的点集完全一样，问题才能复现。

压力对照的规模是 5 种点集（随机、网格、圆环、三簇、近共线）× 11 种点数（3 到 200）× 30 个种子。网格是刻意加进去的：大量四点共圆，Delaunay 剖分本身不唯一；圆环和三簇压不同的密度分布，近共线专门压退化路径。

外部对拍与页内恒等式

独立实现互相比对，只能证明两份实现不矛盾，不能证明它们都对。真正的裁判在项目外面：同一批算例导出成 JSON，用 scipy 和 shapely 逐项核对。

• 凸包顶点集与面积：对 scipy.spatial.ConvexHull；
• 三角形集合（一般位置点集）：对 scipy.spatial.Delaunay，连三角形下标都能对上；
• 三角形面积和等于凸包面积：对 ConvexHull.volume，相对差最大 4.74e-16；
• 每个 Voronoi 单元的面积：对 shapely.voronoi_polygons（底层是 GEOS），最大相对差 4.83e-14；
• Voronoi 顶点等于外接圆心：对 scipy.spatial.Voronoi。

11 个算例、91 条判据（含两条内部路径的交叉核对）全部通过。但这里要说清楚它证明了什么、没证明什么：它证明这几套算法在一般位置和常见退化点集上，输出与成熟实现逐项吻合；它不证明对所有病态输入都稳。遇到共圆点集，三角形集合这一项只比面积和与三角形数，不比具体下标；坐标跨度极端的点集没有测。

除了外部对拍，还有几条不依赖外部库的恒等式，随页面实时计算：

• 三角形面积和 = 凸包面积；
• 欧拉关系 3F = 2E − h，h 是只属于一个三角形的边数（凸包边界边数）；
• Voronoi 单元面积和 = 画布面积，单元两两不重叠，站点落在自己的单元内；
• 对凸包顶点再求一次凸包，结果不变（幂等）。

这些恒等式成本低、覆盖面广，适合放在每次重算之后跑一遍。

空圆性检查：绝对容差会骗人

Delaunay 的定义是「每个三角形的外接圆内没有其他点」。按定义直接写检查，很容易给一个固定的绝对容差：EPS 是个很小的常数，跟外接圆半径的平方没有关系。一般点集上没问题，网格点集上立刻翻车：5×5 网格的外接圆半径平方在 10⁴ 以上，一旦 EPS 小于这个尺度下的浮点舍入误差，圆周上的点就会被算成「进了圆内」，一次报出 32 处违规。这不是剖分错了，是判据错了。

改成按半径平方取相对容差：

const r2 = cc.r * cc.r;
const tol = r2 * 1e-6;
const dx = p.x - cc.cx, dy = p.y - cc.cy;
if (dx * dx + dy * dy < r2 - tol) bad.push({ triangle: ti, point: i });

同一批网格点集的违规数归零。几何判据里的容差最好跟着尺度走，固定绝对容差等于给坐标范围设了一个隐含上限。

面积恒等式抓住过百万分之一量级的缺口

最初的三角剖分用的是「超三角形 + Bowyer–Watson 坏三角形集合」。当时的记录里，1920 个随机算例有 4 个（约 0.2%），三角形面积和比凸包面积小一点点。缺口是贴着凸包边的一条细条，肉眼几乎看不见。

根因是贴近凸包边的点会跟超三角形顶点连成三角形，最后删除超三角形时，连那一小块一起删掉了。换成「凸包扇形起步 + 内部点增量插入 + Lawson 翻边」之后，同一批算例零失败。

这件事的教训不是「Bowyer–Watson 不好」，而是：面积和等于凸包面积这条恒等式，比「看起来没洞」可靠得多。它把百万分之一量级的缺口也抓出来了，而且不需要任何外部依赖。

几条可以迁移的验证思路

1. 同一个结构写两遍，先想清楚独立到什么程度。共用底层原语可以接受，但要承认对拍覆盖不到它；共享得越多，同时错的概率越高。
2. 对拍逐项、定量、可复现，并专门构造退化输入。固定种子、固定阈值、逐单元比较，失败时能定位到具体站点；共线、共圆、重复点、网格这些才是浮点代码容易翻车的地方。
3. 恒等式和判据本身都要被检验。面积和、欧拉关系、两两不重叠很便宜也很敏感；空圆性用绝对容差误报 32 处，问题不在算法，在判据。

适用范围与局限

这套实现有明确的适用范围：不支持带权 Voronoi（每个站点权重不同），也不做约束 Delaunay（强制某些边必须出现）；输入是点集而不是自交多边形；重复点会被去重，全共线的点集退化成两个端点，并显式给出退化标记，不会硬算出一堆零面积三角形。浮点退化靠容差挡，不保证任意病态输入都正确；工业级实现要上精确谓词（比如 Shewchuk 那套自适应精度算术），这里没有做。

对拍能做的，是在自己声明的范围内，把「两份实现不一致」和「与成熟实现不一致」这两类问题尽量提前暴露出来。它不能替代精确谓词，也不能证明代码对所有输入都正确——但比起只看页面上画得像不像，它至少把错误变成了一个可以定位的数字。