计算几何第九讲：凸包与最近点对的核心算法与实战

如果你是一个学算法的人，走到“计算几何”这一讲，大概率已经有了一种隐约的心理准备：前面那些排序、图论、动态规划，好歹还能靠“背模板+理解状态转移”硬啃下来。但计算几何不一样，它把坐标、向量、浮点数全部堆到你面前，很多代码在本地跑得好好的，交上去就莫名其妙 WA。第九讲选了“凸包”和“最近点对”这两个经典问题，其实就是在给你补两样最核心的能力：一是用向量判断几何位置关系，二是把“看起来只能暴力”的问题用分治或扫描拆到可接受的复杂度。这篇文章没有华丽的背景介绍，直接讲原理、代码、边界条件和我在实际练习里踩过的坑，适合正在学算法设计、准备机试或者刷竞赛题的人。

1. 别急着写凸包，先花半小时解决“精度”和“方向判断”

很多人学计算几何第一个动作就是打开编辑器抄一份凸包模板，结果怎么调都调不通。我个人的建议是：先别管凸包是什么，先把“判断一个点在一个向量的左侧还是右侧”这件事搞到滚瓜烂熟。因为凸包、最近点对、线段相交、多边形面积，几乎所有计算几何问题最终都会落到这类基本谓词上。

1.1 叉积是一个“带符号的面积”，不是单纯的长度

平面里有三个点 A、B、C，我们想判断 C 在 AB 的顺时针方向还是逆时针方向。最稳的做法不是算角度，而是算叉积：

struct Point { double x, y; }; // 计算向量 AB 与 AC 的叉积 // 返回 (B.x-A.x)*(C.y-A.y) - (B.y-A.y)*(C.x-A.x) double cross(const Point& A, const Point& B, const Point& C) { return (B.x - A.x) * (C.y - A.y) - (B.y - A.y) * (C.x - A.x); }

这个值的几何意义是“向量 AB 到向量 AC 的有向面积的两倍”。如果返回值大于 0，说明 C 在 AB 的逆时针方向；小于 0，说明 C 在顺时针方向；等于 0，说明三点共线。

你可能会问：为什么不用 atan2 算出角度再比大小？因为角度计算涉及反三角函数，速度慢，而且浮点误差更大。叉积只做几次加减乘，在同样的 double 精度下误差累积要小得多，而且不需要处理“360° 回绕”这种麻烦问题。这是计算几何所有模板的基石。

1.2 EPS：所有浮点比较都要留一条“容错带”

直接比较两个 double 是否相等，是计算几何最常见的翻车现场。比如两个点理论上应该重合，但在连续运算后可能分别是 (0.999999999999, 1.000000000001) 和 (1.0, 1.0)。如果你用a.x == b.x判断，结果就是 false，后面整个逻辑全乱。

所以标准做法是引入一个极小量 eps，通常取 1e-8 到 1e-10：

const double EPS = 1e-8; bool same(const Point& a, const Point& b) { return fabs(a.x - b.x) < EPS && fabs(a.y - b.y) < EPS; } int sgn(double x) { if (fabs(x) < EPS) return 0; return x > 0 ? 1 : -1; }

sgn函数会把“接近 0”的数统一当作 0 处理。平时判断叉积时，不要直接写if (cross > 0)，而是写if (sgn(cross) > 0)。我在初学阶段经常省掉这一步，结果构造凸包时遇到三点近乎共线的数据，栈顶元素会被错误地弹出，最后出来的多边形缺一个角。加了 eps 之后，这类“玄学错误”立刻减少了一大半。

1.3 排序比较器：严格弱序里不能有“相等”

无论是凸包算法还是最近点对，第一步都是对点集排序。排序的比较器必须满足“严格弱序”，也就是说，对任意两个点，比较器不能同时返回 a < b 和 b < a。

一个典型的反例是直接按欧几里得距离排序，或者用fabs(a.x - b.x) < EPS当作“相等”后再比 y。这种做法在极少数数据下会导致sort行为未定义，甚至死循环。最安全的写法是：先按 x 排序，如果 x 的差值大于 EPS，直接返回 x 较小；否则再按 y 排序。

bool cmp(const Point& a, const Point& b) { if (fabs(a.x - b.x) > EPS) return a.x < b.x; return a.y < b.y; }

这种比较器看起来朴素，但它保证了“唯一确定的全序关系”，是后续所有算法稳定的前提。我推荐你在写完比较器之后，随手构造一些点，跑一遍sort，再unique去重，确保没有重复点。重复点不处理干净，凸包会遇到除零或死循环，最近点对会得到距离 0 但位置诡异的答案。

2. 凸包：用“弹栈”模拟橡皮筋收紧的过程

凸包的定义很直观：平面上有 n 个点，用一个最小的凸多边形把所有点包住，这个多边形的顶点就是凸包。可以想象成在钉了很多钉子的木板上拉一根橡皮筋，松手后橡皮筋的形状就是凸包。

凸包的用途比你想的更广：求点集的最小包围范围、判断点是否在多边形内、计算凸多边形的直径（旋转卡壳）、做简单的碰撞检测，甚至在地理信息系统里勾勒一片区域的轮廓。面试和竞赛里最常见的考题是“给定点集求凸包顶点”和“求凸包周长 / 面积”。

2.1 三种主流凸包算法怎么选

表格里对比的是三种最常见的实现：

我自己写题和比赛时基本只用 Andrew 单调链，原因很实际：Graham 扫描要先求极角再排序，遇到多个点在同一极角方向上的时候，排序规则要写得很小心，否则顺序错了整个凸包就错。而 Andrew 只需要按 x-y 排序，然后从左到右扫一遍构造下凸壳，从右到左扫一遍构造上凸壳，代码短，逻辑直接。

2.2 Andrew 单调链完整实现

算法分三步：

1. 按 x 升序（x 相同按 y 升序）排序，并去除重复点。
2. 从左到右扫描所有点，维护一个栈，用“新点是否使得栈顶转弯方向不合法”来判断是否弹栈。
3. 从右到左再扫一遍，构建上半部分，最后把两部分拼起来。

所谓“转弯方向不合法”，就是检查：当前栈顶是 p，下一个点是 q，新点是 r。如果向量 p->q 与 q->r 的叉积小于等于 0（或小于 0，取决于是否保留共线点），说明 r 在 p->q 的顺时针方向或者同向，这时候不应该继续留着 q 作为凸包顶点，而要弹栈。

这里有一个常见分歧：保留共线点还是不保留。如果题目要求输出所有在凸包边上的点，比较时用<= 0弹栈；如果只要凸包顶点，不关心共线点，就用< 0弹栈。下面的代码采用“不保留共线点”的常用版：

#include <bits/stdc++.h> using namespace std; const double EPS = 1e-8; struct Point { double x, y; bool operator<(const Point& other) const { if (fabs(x - other.x) > EPS) return x < other.x; return y < other.y; } bool operator==(const Point& other) const { return fabs(x - other.x) < EPS && fabs(y - other.y) < EPS; } }; double cross(const Point& A, const Point& B, const Point& C) { return (B.x - A.x) * (C.y - A.y) - (B.y - A.y) * (C.x - A.x); } vector<Point> convexHull(vector<Point> pts) { int n = pts.size(); if (n <= 1) return pts; sort(pts.begin(), pts.end()); pts.erase(unique(pts.begin(), pts.end()), pts.end()); n = pts.size(); if (n <= 2) return pts; vector<Point> hull(2 * n); int k = 0; // 下凸壳 for (int i = 0; i < n; ++i) { while (k >= 2 && cross(hull[k-2], hull[k-1], pts[i]) <= EPS) { --k; } hull[k++] = pts[i]; } // 上凸壳 for (int i = n - 2, t = k + 1; i >= 0; --i) { while (k >= t && cross(hull[k-2], hull[k-1], pts[i]) <= EPS) { --k; } hull[k++] = pts[i]; } hull.resize(k - 1); return hull; }

注意看第二个循环，我从n-2开始，因为最右边的点已经在第一个循环里放进去了。t = k + 1保留了下凸壳的右端点，避免在上凸壳构建时把它弹掉。最后resize(k - 1)是因为起点被重复加了一次，去掉重复的起点正好把首尾相连的多边形顶点序列凑齐。

2.3 一个小例子，手动跑一遍

假设点集是：(0, 0), (1, 1), (2, 2), (1, 0), (0, 2)。按 x-y 排序后是 (0,0), (0,2), (1,0), (1,1), (2,2)。

第一个循环：

• 压入 (0,0)，压入 (0,2)。
• 遇到 (1,0)，检查 (0,0)->(0,2)->(1,0) 的叉积，这个向量是顺时针方向，所以弹掉 (0,2)，压入 (1,0)。
• 压入 (1,1) 和 (2,2)。

第二个循环从右往左，会把中间凹陷的点弹掉，最终得到 (0,0), (1,0), (2,2), (0,2)。这就是正确凸包。

这个例子虽然小，但足以说明“弹栈”的本质：栈里保存的是目前已经确定可能在凸包上的点，当新点让最近两个点的转弯方向变“错”，那中间那个点就不可能是凸包顶点。整个过程和橡皮筋收拢的思路一致。

2.4 边界情况：不是所有点集都能生成“多边形”

很多新手栽在边界上。点数只有 1 或 2，凸包退化成点或线段，这时候hull.size()返回 1 或 2，代码里不能去访问hull[2]。所有点共线时，第一个循环会把所有点顺序压进去，第二个循环可能因为k >= t的条件而完全不弹栈，最后返回的是一个包含两个端点、中间还夹着一些共线点的大集合。如果你只需要两端点，建议在凸包构建后单独判断：如果所有点都满足sgn(cross(pts[0], pts[1], pts[i])) == 0，直接返回{pts.front(), pts.back()}。

另外，题目如果只要求周长，别忘了计算的周长是凸包上相邻点的欧几里得距离之和，而hull最后一个点和第一个点也是相邻的。我记得第一次写凸包周长时候，末尾点距丢掉了，硬是差了一个边的长度，最后对拍才查出来。

3. 最近点对：分治的真正威力在于“剪掉”没用的距离

如果说凸包练的是“向量和栈”，那最近点对练的就是“分治 + 合并时如何减少无效比较”。问题很简单：给 n 个点，找一对欧氏距离最小的点。暴力做法是 O(n²)，n 到 1e5 就完全跑不动了。分治做法的复杂度是 O(n log n)，代码量并不大，但合并时的边界条件是很多人卡住的核心。

3.1 分治流程：三分治，七分合并

先把点按 x 坐标排序，然后递归处理。

假设现在处理的是区间 [l, r] 内的点：

• 如果区间里只剩一个点，返回无穷大。
• 如果区间里剩两个点，直接返回它们距离。
• 取中点 mid，递归求左半部分最小距离 d1、右半部分最小距离 d2，当前已知最小距离 d = min(d1, d2)。
• 合并时，只有横坐标与中线距离小于 d 的点才可能刷新答案。
• 把这些点按 y 坐标排序，然后对每个点只和它后面“y 差距小于 d”的点比较。

很多人会问：合并时为什么只比较中线附近 d 距离以内的点？因为如果某个点对的距离小于 d，但它们的横坐标都距离中线超过 d，那么两者横坐标差就超过 d，距离必然超过 d，不可能成为新的最优解。所以我们可以大胆地把距离中线超过 d 的点全部排除。

3.2 为什么内层比较只需要看常数个点

这是整个算法最精巧的地方。假设当前最小距离是 d，我们选取中线附近的带状区域，把区域内所有点按 y 排序。对于其中一个点 p，只需要检查它后面的那些“y 坐标差小于 d”的点。

为什么不会出现一个点后面有几十个点都要检查？因为一旦某两个点的距离小于 d，它们不会同时出现在一个 2d × d 的矩形区域内。更具体地说，把 d × 2d 的区域划分成若干个 d/2 × d/2 的小方格，每个小方格内最多只能有一个点；在 2d × d 的范围内，点的数量被限制在常数级别。实际编码里，很多人直接写“每个点最多检查后面的 7 个点”或者“当 y 坐标差 >= d 就 break”，两种写法都对。我推荐用 y 差 break 的写法，不容易出错，而且不依赖具体的点数量上界。

这部分如果不理解也不影响你写出能 AC 的代码，但如果面试官问你“为什么复杂度是 O(n log n) 而不是 O(n²)”，你需要能讲清楚“常数点”这个核心。

3.3 最近点对完整实现

#include <bits/stdc++.h> using namespace std; const double INF = 1e20; const double EPS = 1e-8; struct Point { double x, y; }; double dist(const Point& a, const Point& b) { double dx = a.x - b.x; double dy = a.y - b.y; return sqrt(dx * dx + dy * dy); } bool cmpX(const Point& a, const Point& b) { if (fabs(a.x - b.x) > EPS) return a.x < b.x; return a.y < b.y; } bool cmpY(const Point& a, const Point& b) { if (fabs(a.y - b.y) > EPS) return a.y < b.y; return a.x < b.x; } double closestPair(vector<Point>& pts, int l, int r) { if (r - l <= 1) { return INF; } int mid = (l + r) / 2; double midX = pts[mid].x; double d = min(closestPair(pts, l, mid), closestPair(pts, mid, r)); vector<Point> strip; for (int i = l; i < r; ++i) { if (fabs(pts[i].x - midX) < d) { strip.push_back(pts[i]); } } sort(strip.begin(), strip.end(), cmpY); int m = strip.size(); for (int i = 0; i < m; ++i) { for (int j = i + 1; j < m; ++j) { if (strip[j].y - strip[i].y >= d) break; d = min(d, dist(strip[i], strip[j])); } } return d; } double solve(vector<Point>& pts) { sort(pts.begin(), pts.end(), cmpX); return closestPair(pts, 0, pts.size()); }

这段代码的关键细节有两个。第一，递归写的是closestPair(pts, l, mid)和closestPair(pts, mid, r)，半个闭半个开，避免mid点被重复计算。第二，合并时strip用的是midX，这是“当前区间中点的 x 坐标”，不是整个点集的中线。如果你误用了全局中线，会在递归子区间时把一些实际不在中线附近的点拖进来，导致复杂度退化甚至答案错误。

3.4 一个容易忽略的性能优化

上面的代码每次合并都对新产生的strip做一次sort，排序复杂度是 O(k log k)。在分治每一层，所有 strip 加起来的长度是 O(n)，所以总复杂度是 O(n log² n)。理论上仍然能过绝大多数题，但如果是 1e5 以上的点，可能有点悬。

经典的优化方法是：在递归过程中维护一个按 y 排序的辅助数组，通过归并的方式让每层合并维持在 O(n)，总复杂度严格降到 O(n log n)。不过说实话，我大部分时候用上面这版已经能过题了。面试时可以先讲 O(n log² n) 版本，再提一句“可以用归并优化到严格 O(n log n)”，这反而显得你考虑周到。

3.5 手推一个随机小数据

点集是： (1, 2), (4, 7), (3, 1), (9, 8), (5, 3), (8, 2)

按 x 排序后： (1,2), (3,1), (4,7), (5,3), (8,2), (9,8)

mid 取 3，左半是 (1,2),(3,1),(4,7)，递归后左半最小距离应该是 (1,2) 和 (3,1)，距离约 2.236；右半是 (5,3),(8,2),(9,8)，递归后最小距离约 3.162，所以当前 d = 2.236。

合并时，只有横坐标在 4±2.236 范围内的点进 strip。结果里 (4,7) 和 (5,3) 的距离约 4.123，(1,2) 和 (3,1) 已经在左半算过，不会重复计算。最终答案是 2.236。

这个例子虽然简单，但能看到合并的本质：横向距离过大的点被早早淘汰，不需要全部两两比较。

4. 从模板到实战：凸包和最近点对怎么写才不白练

很多人刷题有一个毛病：看完模板，觉得懂了，复制到编辑器里，找一道模板题 AC，然后就不管了。计算几何要是这样练，下次遇到变体题必挂。

4.1 一定要学会造随机数据对拍

我强烈建议你写一个naive版本（暴力 O(n²) 做最近点对，或直接暴力枚举所有点对求最小距离），然后用脚本生成随机点，不断比较暴力结果和分治结果是否一致。

对拍是让你“发现边界 bug”的最快方式。比如我曾经在最近点对代码里只对 x 排序、忘了重排 y，导致少数随机数据下结果偏大，这种错误靠肉眼根本看不出来，只有对拍能逮住。

下面是一个简单的 Python 对拍脚本思路：

import random for test in range(1000): n = random.randint(1, 20) points = [(random.randint(-100, 100), random.randint(-100, 100)) for _ in range(n)] # 用 C++ 程序跑一份输出，和暴力结果对照

别追求测试数据大，20 个点以内就能暴露全部逻辑问题。等你对拍通过 1000 组，再上 1e5 的大数据测性能，这样心里才有底。

4.2 共线点、重复点、浮点误差是三座大山

• 重复点：最近点对的最短距离会变成 0，如果没有去重，可能导致分治代码把结果算成“两个重复点之间”的距离，虽然数值正确，但在某些输出要求“输出点对”的问题里会输出一堆相同的点。
• 共线点：凸包如果不处理共线，会出现凸包顶点上多出一些中间点。如果是求周长或面积，问题不大；如果是要求输出顶点序列，就会 WA。
• 浮点误差：别忘了用sgn判断叉积符号。我在一次比赛里，因为一个cross > 0少写了sgn，导致凸包在极大数据下偶尔会漏掉一个点，那个 case 恰好卡在 0.0000001 的误差上，调了整整一个晚上。

4.3 两道题帮你串联这讲的两个算法

一道是“Wall”（POJ 1113）：给 n 个点，求一个多边形，多边形上的点到所有点的距离不小于 L，求这个多边形的最小周长。答案其实就是凸包周长加上一个半径为 L 的圆的周长。你看，凸包算完，还要想到“外扩距离”用圆弧补上，而不是简单地周长 += L。

另一道是“平面最近点对”（洛谷 P1429 或 POJ 3714）。前者是纯模板，后者要求在一个点集和另一个点集之间求最近距离，做合并时要注意区分两个集合里的点，这是对分治的进一步理解：递归时不能把同集合的点对算进去。

建议先 AC 这两道，再回头读一遍本文的边界情况部分，很多当时没感觉的知识点就串起来了。

5. 进一步的下一步：旋转卡壳、扫描线，以及我的学习心得

如果第九讲的内容你已经完全拿下了，下一步通常会接旋转卡壳（Rotating Calipers），用它求凸包直径、最大三角形面积、最小宽度等等。旋转卡壳的思想并不难：在凸包上拿一条边作为基准，让另一个点沿着凸包移动，直到距离达到极值，然后旋转到下一条边。它和凸包一样，也是“维护极值、利用单调性避免重复扫描”的典型。

还有一种常见扩展是扫描线类问题，比如判断一堆线段是否相交、求矩形面积并。它们的共同点都是“先排序，再用事件点/单调结构维护当前扫描线状态”，和最近点对里“按 y 排序后只查有限范围”的思路一脉相承。所以我常说，计算几何的第一课不是背模板，而是理解“怎么减少没必要的比较”。

我这里再分享一个个人体会：计算几何是最容易靠“画图”来 debug 的算法领域。用 matplotlib 把点、凸包节点、中线和最近点对用不同颜色画出来，再标上数字，很多 bug 一眼就能看出来。我写的每个凸包模板都会先用随机点画一次图确认形状，再上对拍。这比盯着代码空想要高效得多。

如果你刚学到第九讲，我的建议是：先不用追求学完所有算法，把叉积、eps、凸包、最近点对这四样彻底吃透，配合对拍脚本，刷 5 道左右的相关题目，计算几何的地基就稳了。之后再去碰旋转卡壳、半平面交，你会发现很多“新概念”实际上都是旧思维的变体。