ARTICLE DETAIL

资讯详情

深耕网站建设、视觉设计与SEO优化的一线实战洞察。

计算几何四模块:点积叉积、线段相交、凸包与多边形面积实战

计算几何四模块:点积叉积、线段相交、凸包与多边形面积实战

计算几何里的点积叉积、线段相交、凸包算法、多边形面积,这四个东西几乎覆盖了二维几何编程八成的日常需求。我最早接触它们是在做一个图形编辑器的碰撞检测模块,当时以为把公式抄下来就能跑,结果上线之后遇到一堆"肉眼看着明明相交却判成不相交"的诡异问题,排查了两天才发现是浮点误差和退化情形没处理。后来做地图标注、做物料排版、做游戏里的小地图裁剪,又反复用到这套东西,慢慢攒了一些真正能落地的经验。网上大部分教程只讲公式不讲边界,代码跑得过样例,跑不过真实数据。这篇就按我自己的实践顺序把四块内容串一遍,从最基础的向量运算讲到能直接抄去用的完整实现,顺带把踩过的坑和排查思路摊开说清楚。适合刚接触计算几何的同学当入门到进阶的过渡材料,也适合已经写过一些几何代码但总在边界情况上翻车的朋友对照检查。

1. 先把这几个模块的关系捋清楚

1.1 为什么它们总是被放在一起讲

这四个知识点看起来是并列的,实际上是有依赖关系的。点积和叉积是原子运算,线段相交判定是叉积的第一个正经应用,凸包是线段拐向判断的批量版本,多边形面积则是叉积的另一个经典用法。你把叉积吃透了,后面三个基本就是"换个场景调用同一个函数"。

我在实际项目里搭几何库的时候,第一件事永远是写一个稳健的cross和dot,然后在它们外层包一个带EPS比较的版本。这个习惯救过我好几次。因为几何问题的 bug 有个特点:它不会在单点测试里暴露,而是会在某个特定输入组合下突然炸掉,而且往往是在线上环境,数据量大、坐标分布杂的时候。

举个我遇到过的真实场景:某个排版工具需要判断两个矩形贴纸有没有重叠、如果有重叠需要沿边缘切开。最初的实现用的是简单的坐标区间比较,矩形轴对齐时完全没问题,后来产品要求支持旋转一定角度的贴纸,代码立刻失效。改成线段相交判定加凸包相交判定之后才稳下来。这个转变让我彻底明白一件事:几何问题不要试图用特例去糊,老老实实回到向量运算才是长期方案。

1.2 选型时的三个决策点

真正动手前,我一般会先确认三件事,这三件事的选择会直接影响后面所有代码的写法。

第一是坐标的量级。你的坐标是做屏幕像素级的(0 到几千),还是做地理经纬度的(小数点后六七位),还是做 CAD 那种高精度的。量级决定了EPS该取多少。用 1e-9 去判定像素级坐标,等于没判;用 1e-2 去判定经纬度,精度又全丢了。

第二是要不要保留共线点。做凸包的时候,多边形的边上如果有三点共线,你是把中间那个点踢掉,还是保留?大部分情况下踢掉更省事,但如果这个凸包后面要用来做碰撞的接触点,保留反而更好,因为接触面可能是一整条边。这个决策要在写第一行代码前定下来,不然改起来牵连一大片。

第三是用整数还是浮点。如果输入本身就是整数坐标(比如网格地图、棋盘、像素图片),我强烈建议全程用整数运算,只在最后算距离或者面积的时候转浮点。这样叉积的符号判断是绝对精确的,一个EPS都不需要。这个技巧能消灭掉至少一半的精度问题。

1.3 一个被低估的前置工作:统一数据结构

我见过太多项目里每个函数自己定义一套点结构,Point、Vec2、Vector2D、pt到处都是,字段名有的叫x/y,有的叫px/py,有的叫first/second。写几何代码本来就容易出错,再加上类型转换,出错的概率指数级上升。

我的做法是:整个项目只保留一个二维点类型,同时提供向量和点的双重语义(数学上它们本来就是同构的),然后所有几何函数都围绕它展开。C++ 里就是一个struct加一堆运算符重载,Python 里可以用dataclass或者干脆用元组加纯函数。下面是我常用的 Python 版本骨架,后面所有代码都基于它:

from dataclasses import dataclass import math EPS = 1e-9 @dataclass(frozen=True) class Pt: x: float y: float def __add__(self, o): return Pt(self.x + o.x, self.y + o.y) def __sub__(self, o): return Pt(self.x - o.x, self.y - o.y) def __mul__(self, k): return Pt(self.x * k, self.y * k) def sgn(v: float) -> int: """带容差的三态符号函数,整个库的核心守卫""" if v > EPS: return 1 if v < -EPS: return -1 return 0

sgn这个函数看着不起眼,但它是我所有几何代码的地基。所有需要"判断正负"的地方一律调它,绝不直接写> 0。这样做的好处是精度策略集中在一个地方,哪天发现EPS不合适,改一个常量就够了,不用去翻几百行代码。

2. 点积与叉积:所有判断的地基

2.1 点积在回答什么问题

点积的定义很朴素:a · b = ax*bx + ay*by,几何意义是|a||b|cosθ。我习惯把它理解成"a 在 b 方向上有多少分量",也就是投影。

这个理解比背公式有用得多。比如你要判断两个向量是不是锐角,取点积看正负就够了,因为|a||b|恒为正,符号完全由cosθ决定。同理钝角是负,直角是零。

我在项目里点积用得最多的地方是点到直线的最短距离。给定直线上两点 A、B 和任意点 P,先算t = dot(P-A, B-A) / dot(B-A, B-A)。这个t就是垂足在直线上的参数位置:t在 0 到 1 之间说明垂足落在线段内部,小于 0 说明离 A 端最近,大于 1 说明离 B 端最近。算出t之后把 P 投影过去,两点距离就是答案。

def point_seg_dist(p: Pt, a: Pt, b: Pt) -> float: ab = b - a if ab.x == 0 and ab.y == 0: # 退化成点 return math.hypot(p.x - a.x, p.y - a.y) t = ((p.x - a.x) * ab.x + (p.y - a.y) * ab.y) / (ab.x*ab.x + ab.y*ab.y) t = max(0.0, min(1.0, t)) # 关键:钳制到线段范围内 proj = Pt(a.x + ab.x * t, a.y + ab.y * t) return math.hypot(p.x - proj.x, p.y - proj.y)

这里那个max/min钳制是整个函数最有价值的一行。很多人写点到直线距离很顺,一到点到线段距离就翻车,原因就是忘了钳制。不钳制的话,垂足跑到线段延长线上去了,算出来的是点到直线的距离,不是到线段的距离。这个 bug 在碰撞检测里尤其致命——玩家明明离墙还有一段距离,系统判定你已经贴上了。

注意:钳制之后t的精度损失会累积,如果你要做的是"点是否在线段上"这种判定,不要用距离阈值,用叉积判共线加坐标区间判,精度高得多。

2.2 叉积的符号比大小更重要

叉积是二维里的"伪标量":a × b = ax*by - ay*bx,几何意义是|a||b|sinθ,同时它的绝对值等于以 a、b 为邻边的平行四边形面积。

我在代码里 90% 的叉积调用只关心符号,不关心数值。符号告诉你的信息量非常大:cross(b-a, c-a)为正说明 c 在向量 ab 的左侧,为负说明在右侧,为零说明三点共线。就这么一句话,撑起了后面线段相交和凸包的全部逻辑。

为什么是左侧?想象你站在 a 点面朝 b 点,c 在你的左手边,那就是正。这个直觉在调试的时候特别好用,你可以直接在纸上画一画验证代码算出来的符号对不对,不用去推导。

叉积的另一个高频用法是算三角形面积:S = |cross(b-a, c-a)| / 2。由此可以扩展到任意多边形,这就是第 5 节要讲的鞋带公式。顺带提一句,叉积还能判断多边形的顶点是顺时针还是逆时针排列,把相邻边两两叉积加起来看总符号,为正就是逆时针,为负就是顺时针。这个判断在渲染和填充的时候必须先做,不然填充区域会跑到外面去。

2.3 EPS 到底该取多少

这个问题我被问过很多次,也踩过不少坑。先说结论:没有一个放之四海皆准的值,得看坐标量级和运算次数。

一个比较可靠的估算方式是:如果你的坐标量级是M(比如最大坐标是 1e6),双精度浮点的相对误差约 1e-16,那么单次叉积运算的绝对误差上限大概是M * M * 1e-16。当M = 1e6时,误差上限就是 1e-4,这时候你取EPS = 1e-9完全没意义,因为误差本身就比它大三个数量级。

所以我一般在库里提供一个"按量级自适应"的 EPS:

def eps_for(pts) -> float: scale = max(max(abs(p.x), abs(p.y)) for p in pts) return max(1e-12, scale * scale * 1e-12)

这个策略偏保守,但比我见过的直接写死 1e-9 要稳得多。当然如果你的坐标是整数或者量级很小(比如像素级),完全可以直接取 1e-9,甚至用整数运算彻底绕开。

还有一个技巧是比较前先做归一化。判断两个叉积结果是否相等时,不要比绝对值,而是把两个叉积都除以参与的向量长度,比较"正弦值"是否相等。这样结果是无量纲的,不受坐标量级影响。这个做法在判断角度关系时特别好用。

3. 线段相交判定:跨立实验的完整推导

3.1 跨立实验的核心思想

线段 AB 和线段 CD 相交,本质上要求两件事同时成立:C 和 D 分居直线 AB 的两侧,并且 A 和 B 分居直线 CD 的两侧。这就是所谓的"跨立"。

用叉积表达极其简洁。d1 = cross(B-A, C-A),d2 = cross(B-A, D-A),如果d1 * d2 < 0,说明符号相反,C、D 在 AB 两侧。同理算d3 = cross(D-C, A-C)、d4 = cross(D-C, B-C),如果d3 * d4 < 0,A、B 在 CD 两侧。两个条件都满足,就是严格相交。

但这里有个陷阱:直接乘d1 * d2可能会溢出。如果坐标量级大,两个叉积都是 1e200 这个级别,乘起来就 inf 了,符号判断直接失效。稳妥的写法是比较符号:

def straddle(a: Pt, b: Pt, c: Pt, d: Pt) -> bool: d1 = cross(a, b, c) d2 = cross(a, b, d) d3 = cross(c, d, a) d4 = cross(c, d, b) return (sgn(d1) * sgn(d2) < 0) and (sgn(d3) * sgn(d4) < 0)

sgn返回值只有 -1、0、1,乘起来绝对安全。

3.2 四种退化情形必须分开处理

严格相交只覆盖了"两条线段在内部某点交叉"这一种情况。真实数据里至少还有四种变体:

情形几何特征判定条件
端点接触一个端点落在另一条线段内部叉积为 0 且点在段内
端点重合两条线段共用端点坐标相等
共线分离四点共线但区间不重叠坐标区间不交
共线重叠四点共线且有公共部分坐标区间相交

这四种情况如果不处理,你的相交判定会在"两条线段首尾相接"时返回 false,然后整条路径的连通性判断就断了。我做地图路径规划的时候就被这个坑过一次,两条道路在路口共用角点,结果判定不连通,绕了一大圈才排查到。

完整的实现我一般写成这样:

def cross(a: Pt, b: Pt, c: Pt) -> float: return (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x) def on_seg(a: Pt, b: Pt, p: Pt) -> bool: """前提:p 与 a、b 共线,判断 p 是否落在线段 ab 上""" return (min(a.x, b.x) - EPS <= p.x <= max(a.x, b.x) + EPS and min(a.y, b.y) - EPS <= p.y <= max(a.y, b.y) + EPS) def seg_intersect(a: Pt, b: Pt, c: Pt, d: Pt) -> bool: d1 = sgn(cross(a, b, c)) d2 = sgn(cross(a, b, d)) d3 = sgn(cross(c, d, a)) d4 = sgn(cross(c, d, b)) # 严格跨立 if d1 * d2 < 0 and d3 * d4 < 0: return True # 退化:端点落在对方线段上 if d1 == 0 and on_seg(a, b, c): return True if d2 == 0 and on_seg(a, b, d): return True if d3 == 0 and on_seg(c, d, a): return True if d4 == 0 and on_seg(c, d, b): return True return False

这个版本我用了好几年,只要EPS选得合理,基本没出过问题。

3.3 快速排斥实验:要不要加

很多教材会先做一步"快速排斥":如果两条线段的外接矩形都不重叠,那肯定不相交,直接返回 false,省掉后面的叉积计算。

要不要加取决于你的数据。如果大部分线段对都是远离的(比如空间索引的初步筛选),这个剪枝能带来 3 到 5 倍的提速。但如果你的数据本来就密集,大部分线段对都靠得很近,那这步剪枝的命中率很低,反而多了一堆比较开销。

我实测下来的经验是:在暴力两两判定的场景(n 小于 2000)加上它,在已经有空间索引(网格、四叉树)的场景不要加。因为索引已经把大部分无关的过滤掉了,剩下的候选本来就大概率相交,排斥实验几乎不生效。

注意:快速排斥必须用带容差的比较,max(a.x, b.x) < min(c.x, d.x) - EPS这样写,不能直接用<。否则刚好贴边的情况会被误判成不相交,又是一类难查的 bug。

3.4 一组必要的测试用例

写完相交判定,我强烈建议先跑一遍下面这些用例。它们覆盖了我这些年遇到过的全部退化情形,能过这组测试基本就不会有大问题了。

tests = [ ((0,0),(2,2),(0,2),(2,0), True), # 标准交叉 ((0,0),(2,0),(3,0),(5,0), False), # 共线分离 ((0,0),(4,0),(2,0),(6,0), True), # 共线重叠 ((0,0),(2,0),(2,0),(4,0), True), # 端点重合 ((0,0),(4,0),(2,0),(2,2), True), # T 型接触 ((0,0),(1,0),(2,0),(3,0), False), # 首尾相接但之间有缝 ((0,0),(1,1),(2,2),(3,3), True), # 完全共线重叠 ((0,0),(1,0),(0,1),(1,1), False), # 平行不相交 ]

这组用例里我最想强调后三个。平行不相交看着简单,但因为d1 == d2 == 0,如果没有处理好会走到"端点落在线段上"的分支里误判成相交。测试跑通了再谈优化。

4. 凸包:从朴素解到单调链

4.1 为什么最终都收敛到单调链

凸包的朴素算法是 O(n³):枚举每一对点组成的直线,判断所有其他点是不是都在同一侧。n 到了一千就卡得跑不动了。

中间层的方案是 Graham 扫描,O(n log n),核心是按极角排序之后用栈维护拐向。这个算法我实现过三次,每次都卡在同一个地方:极角排序的基准点怎么选、共线点怎么处理、精度不够时排序不稳定怎么办。特别是极角排序要用atan2的时候,浮点误差会让两个本该相邻的点顺序颠倒,然后整个栈就崩了。

现在我一律用Andrew 单调链。它只按坐标排序,(x, y)字典序,完全不涉及角度和三角函数,精度天然稳健。思路也直观:把点排序之后,从左到右扫一遍构建下凸包,再从右到左扫一遍构建上凸包,拼起来就是完整的凸包。

def convex_hull(points): pts = sorted(set(points)) # 排序 + 去重 if len(pts) <= 2: return pts def build(seq): stack = [] for p in seq: while len(stack) >= 2 and sgn(cross(stack[-2], stack[-1], p)) <= 0: stack.pop() stack.append(p) return stack lower = build(pts) upper = build(reversed(pts)) return lower[:-1] + upper[:-1]

这十几行代码能处理绝大多数情况,我把它抄在笔记里用了很多年。

4.2 那个<=和<的差别决定成败

上面代码里最关键的是那个sgn(...) <= 0。这个符号决定了共线点怎么处理:

  • 用<= 0:遇到共线就弹出,最终凸包只保留"极端"顶点,边上的中间点被删掉。
  • 用< 0:只弹出右转的点,共线点保留在栈里。

选哪个完全取决于你的下流需求。如果凸包后面要算面积、做包含判定,用<= 0更干净,顶点数最少。如果凸包要做碰撞接触面、要做贴边吸附,必须用< 0保留共线点,否则一条直线边会被压缩成两个端点,中间的接触点全丢了。

我第一次做贴边吸附的时候就是用了<= 0,结果物体贴着斜坡滑动时抖得厉害,因为斜面上的接触点一直在跳。改成< 0之后立刻就稳了。这个教训挺深的:凸包的算法细节和业务语义是绑定的,不能闭着眼睛抄。

提示:lower[:-1] + upper[:-1]这个写法是为了避免首尾两个点重复出现。如果你用的是< 0版本并且输入里有共线点,这句需要额外检查一下重复,用dict.fromkeys去个重更保险。

4.3 单调链的手工推演

算法光看代码容易懵,我拿一个具体例子走一遍。假设点集是(0,0), (1,1), (2,0), (3,1), (4,0), (2,-1),排序后是(0,0), (1,1), (2,-1), (2,0), (3,1), (4,0)。

构建下凸包时依次处理:

  1. 压入(0,0)、(1,1)
  2. 来(2,-1),cross((0,0),(1,1),(2,-1)) = 1*(-1) - 1*2 = -3,负,右转,保留。栈变(0,0),(1,1),(2,-1)
  3. 来(2,0),cross((1,1),(2,-1),(2,0)) = 1*0 - (-2)*1 = 2,正,左转,弹出(2,-1)。再检查cross((0,0),(1,1),(2,0)) = 0,用<= 0会弹出(1,1)。栈变(0,0),(2,0)
  4. 来(3,1),cross((0,0),(2,0),(3,1)) = 2*1 - 0 = 2,正,保留。栈(0,0),(2,0),(3,1)
  5. 来(4,0),cross((2,0),(3,1),(4,0)) = 1*0 - 1*2 = -2,负,保留。栈(0,0),(2,0),(3,1),(4,0)

下凸包完成。上凸包反向扫一遍,最后拼起来。手工推一遍之后,这个算法在我脑子里就再也没忘过,比看十遍代码都管用。建议每个人至少手推一次。

4.4 凸包的几个实用扩展

凸包本身不难,难的是它衍生出来的一堆需求。我按使用频率排一下。

点到凸包的距离。如果是判断点在凸包外有多远,用二分找最近边。如果是点在内,距离为 0。我一般偷懒直接用凸包的每条边跑一遍point_seg_dist取最小值,O(h) 的复杂度,h 是凸包顶点数,通常很小,实际够用。

凸包包含判定。除了暴力叉积,还有 O(log h) 的二分法:以凸包第一个顶点为原点,按极角二分找到点落在哪个三角形扇区里,再判一次叉积。点数上千的时候这个优化才有意义,否则直接暴力更快。

动态凸包。如果点集一直在增删,单调链就不好用了。这时候要么用垂直分解的树结构,要么用可以维护上下凸包的平衡树。我做过一版,代码量比静态版本大五倍,除非真有性能瓶颈,否则不建议上。

5. 多边形面积:鞋带公式及其边界

5.1 鞋带公式的推导与符号含义

鞋带公式的名字很形象,因为它算的时候像是在系鞋带:交叉相乘再相减。

Area = 0.5 * |Σ(x_i * y_{i+1} - x_{i+1} * y_i)|

注意下标是绕圈的,最后一个点和第一个点也要算一次。推导思路很简单:把多边形从一个原点出发切成若干个三角形,每个三角形的有向面积用叉积算,加起来。由于内部有向面积会互相抵消,最终剩下的正好是多边形面积。

def polygon_area(poly) -> float: n = len(poly) s = 0.0 for i in range(n): x1, y1 = poly[i].x, poly[i].y x2, y2 = poly[(i + 1) % n].x, poly[(i + 1) % n].y s += x1 * y2 - x2 * y1 return abs(s) * 0.5

这里我特意保留了abs之外的符号信息,因为符号本身有语义:正值表示逆时针,负值表示顺时针。这个符号在判断多边形朝向时非常有用,比单独写一个遍历叉积的函数要快。

5.2 别急着取绝对值

我见过不少代码上来就abs,把符号丢了。然后下游要判断顶点顺序的时候,又得重新扫一遍所有边算叉积和。其实鞋带公式已经算出来了,白扔。

我的习惯是写一个signed_area返回带符号的值,再包一个polygon_area返回绝对值。朝向判断直接看signed_area > 0。代码里少一次遍历,也就少一次出错机会。

def signed_area(poly) -> float: n = len(poly) return sum(poly[i].x * poly[(i+1) % n].y - poly[(i+1) % n].x * poly[i].y for i in range(n)) * 0.5

另外坐标原点离多边形很远的时候,鞋带公式会有精度损失,因为每个叉积都很大,互相抵消之后剩一个小数,有效位就丢了。稳妥做法是先把所有点平移到以第一个点为原点的坐标系,再算:

def signed_area_stable(poly) -> float: ox, oy = poly[0].x, poly[0].y n = len(poly) s = 0.0 for i in range(1, n - 1): ax, ay = poly[i].x - ox, poly[i].y - oy bx, by = poly[i+1].x - ox, poly[i+1].y - oy s += ax * by - bx * ay return s * 0.5

这个版本实际上是在做三角形扇形求和,数值稳定性明显好于直接套鞋带公式,尤其是在坐标量级 1e6 以上、多边形面积又不大的情况下差别很明显。我现在默认用这个版本。

5.3 自交多边形与凹多边形怎么处理

鞋带公式对凹多边形完全正确,这一点很多人有误解。凹多边形的某些三角形有向面积是负的,正好抵消掉多出来的部分,最终结果依然准确。

真正有问题的是自交多边形。比如画一个"8"字,鞋带公式算出来的结果是两瓣面积的差,不是和。这时候的"面积"怎么定义要看业务:你是想做布尔运算求并集,还是想求绝对覆盖面积。这两种需求对应的算法完全不同,前者需要先做多边形求交,后者需要把多边形拆成若干简单多边形。

我做排版工具的时候就遇到过这个,用户能随手乱画,画出来的图形经常自交。最后的处理方案是先用扫描线把多边形拆成若干个单调多边形,再分别算面积求和。这一步的实现量不小,如果不是必须支持任意手绘输入,我建议在交互层就限制住,比如只允许拖拽出凸多边形。

5.4 面积相关的几个延伸计算

鞋带公式还能顺手带出几个常用量。

形心(重心)。多边形的形心不是顶点的平均值,那个只对三角形成立。正确做法是用带符号面积加权:

def centroid(poly): a = signed_area(poly) if abs(a) < EPS: return Pt(sum(p.x for p in poly)/len(poly), sum(p.y for p in poly)/len(poly)) cx = cy = 0.0 n = len(poly) for i in range(n): p, q = poly[i], poly[(i+1) % n] f = p.x * q.y - q.x * p.y cx += (p.x + q.x) * f cy += (p.y + q.y) * f return Pt(cx / (6 * a), cy / (6 * a))

这个公式我用得很多,做标签居中、做拖拽锚点都要它。注意面积接近 0 的时候要兜底,否则除零。

多边形周长。简单,遍历相邻边求距离累加。但我提一个容易忽略的点:如果首尾两点重复(很多数据格式里多边形首尾是闭合的),累加时会多算一段长度为 0 的边,不影响结果但会多一次 sqrt。养成在入口处检查并去掉重复末点的习惯。

6. 常见问题与排查实录

6.1 精度问题速查表

这类问题占了我在几何代码上花的时间的一大半。整理成表方便对照。

现象可能原因排查手段处理方案
交点算成 NaN两直线平行导致除零打印叉积分母先判平行再求交
明明相交判 falseEPS 太小,误差淹没打印四个叉积值按量级放大 EPS
明明不相交判 true共线分支误命中检查sgn返回值修正on_seg区间
凸包顶点顺序乱排序不稳定打印排序后的点列用元组字典序排序
面积偶尔负顶点顺序不一致检查signed_area符号统一成逆时针输入
大面积算偏小抵消误差平移原点重算用扇形求和版本

我最想强调第二行和第三行。这两个现象看起来矛盾,实际上经常同时存在于一个项目里:一部分数据 EPS 太小漏判,另一部分数据 EPS 太大误判。这时候不要盲目调 EPS,正确做法是把四种退化情形分别写测试用例,逐个确认,找到真正出问题的那一条分支。

6.2 性能排查的两个方向

几何代码慢,绝大多数情况不是算法复杂度问题,而是两个细节没做好。

第一个是开方次数。距离比较里频繁用hypot或sqrt,其实比较距离大小的时候完全可以比平方距离,省掉开方。我在一个 10 万点的最近点对检测里把开方去掉之后,耗时从 1.2 秒降到 0.3 秒。这个优化几乎零成本,只要注意别把平方距离当真距离传给下游就行。

第二个是重复计算。向量减法在循环里被反复算,比如b - a在四条边里各算一次。把它提到循环外面存变量,性能能提升 20% 以上。这个在凸包构建里尤其明显,因为栈操作会反复访问相邻两个点。

至于空间索引,我的经验是点数到 5000 以上再考虑上网格或四叉树。低于这个量级,建索引的开销比重都超过收益。我做过一次对比测试,2000 个点的暴力凸包是 3 毫秒,加索引反而要 5 毫秒。数据量小的时候,简单方案永远是对的。

6.3 我在实际项目里踩过的坑

挑几个印象深的说说。

坑一:整数坐标混用浮点 EPS。有一次做像素级的地图标注,坐标全是整数,但我图省事用了浮点版本加 1e-9 的 EPS。结果两条像素级线段共线判定频繁出错。原因是整数转浮点本身是精确的,但叉积算完之后有些值恰好是 1 附近,sgn判断没问题,问题出在后面的向量缩放上引入了误差。后来改成全程整数运算,一条 EPS 都不需要,问题消失。

坑二:顺时针多边形输入传给期望逆时针的函数。这个错误我犯过不止一次。鞋带公式对两种顺序都能算出正确的面积绝对值,但形心公式会算出一个完全相反的坐标,凸包的单调链会构建出错误结果。现在我在所有接受多边形的函数入口都加了一句朝向检查,不一致就自动翻转顶点顺序。这个防御性编程省了我很多调试时间。

坑三:abs用在了该保留符号的地方。上面提过,但值得再强调。判断点在多边形内还是外、判断多边形是否自交,都需要原始符号信息。养成先 signed 后 abs的习惯。

坑四:重合点没有去重。输入数据里偶尔有连续重复的顶点(通常来自数据转换),凸包算法会因此产生零长度的边,然后叉积算出来全是 0,判断全落在共线分支里,输出结果看着对但顶点全是脏数据。在convex_hull入口加一个set去重就解决了。

心得:几何代码的稳健性不来自聪明的算法,来自对输入数据的怀疑。我现在写任何一个几何函数,第一件事是把所有退化输入列一遍:空集、单点、两点、全共线、全重合、坐标超大量级。把这些都覆盖了,剩下的常规情况基本不会出问题。

7. 把这几个模块拼成一个可用的几何库

单独的四个函数不难写,难的是它们一起工作的时候不出岔子。我现在的做法是配一套统一的约定,整个库内部严格遵守。

第一条约定:所有函数输入的多边形统一为逆时针。入口处用signed_area检查,逆时针直接过,顺时针就reversed一下。这一条消灭了我在朝向问题上 90% 的 bug。

第二条约定:所有返回多个点的函数都保证顶点不重复,包括首尾不重复。凸包函数返回的第一个点和最后一个点不相等,多边形函数不做自动闭合。这条约定让下游代码少一层判断。

第三条约定:所有容差比较走统一的sgn。禁止在任何几何判断里直接写裸的比较运算符。这是纪律问题,违反一次就可能埋一个三个月后爆发的雷。

按这三条约定写出来的库,代码量其实不大,一个文件两三百行,但用起来非常省心。我现在新开项目,直接把这套东西复制过去,改改类型声明就能用。

最后分享一个小技巧,关于几何 bug 的定位。这类 bug 最难的地方是难以可视化,你看到一堆数字,不知道哪里错了。我的办法是写一个异步的小工具,把出问题的输入点列画成 SVG。不用装任何图形库,直接拼字符串在浏览器里打开就行:

def dump_svg(polys, w=800, h=600, scale=40, pad=20): def tx(p): return (p.x * scale + pad, h - (p.y * scale + pad)) out = [f'<svg xmlns="http://www.w3.org/2000/svg" width="{w}" height="{h}" style="background:#fff">'] colors = ["#e63946", "#457b9d", "#2a9d8f", "#e9c46a"] for i, poly in enumerate(polys): pts = " ".join(f"{tx(p)[0]:.2f},{tx(p)[1]:.2f}" for p in poly) out.append(f'<polyline points="{pts}" fill="none" stroke="{colors[i%4]}" stroke-width="1.5"/>') for p in poly: x, y = tx(p) out.append(f'<circle cx="{x:.2f}" cy="{y:.2f}" r="2.5" fill="#333"/>') out.append("</svg>") return "\n".join(out)

把所有参与计算的原始点和结果多边形丢进去,一眼就能看出是哪一步算歪了。这个方法帮我在半小时内定位过好几个此前查了大半天的问题,比打日志猜数字高效太多。你可以在里面加标注文字、加点序号,反正就是拼字符串,想怎么改都行。

实际动手的时候,我还是建议按这个顺序推进:先把点积叉积和sgn写稳,用几十组手算的用例验证;再上线段相交,把四种退化情形全测一遍;然后凸包,记得手推一遍单调链;最后多边形面积,顺手把形心和朝向判断带出来。每一步都验证过再往下走,比一口气写完再统一调试要快得多。

返回列表