计算几何工程实践：叉积、相交判断与凸包的代码落地

1. 这不是数学课，是写代码时绕不开的“空间直觉”训练场

你有没有在写一个地图路径规划功能时，突然卡在“两条线段到底有没有交点”上？翻遍文档发现所有判断逻辑都依赖一个叫“叉积”的东西，但教材里全是向量公式和右手定则，根本不知道它和你正在敲的 if 语句有什么关系；或者调试一个轮廓提取模块，明明视觉上是个凸多边形，程序却总把某个内凹角当成顶点，最后发现是凸包算法没吃透——这些不是算法竞赛选手的专属烦恼，而是任何需要处理点、线、面关系的开发者日常：GIS系统、CAD插件、游戏碰撞检测、甚至UI中的不规则按钮热区判定，全都要和计算几何打交道。

这个标题里的“专题·计算几何入门”，说白了就是一套面向工程实现的坐标系生存手册。它不讲抽象代数，不推导高维空间，只聚焦四件事：怎么用两个浮点数（x, y）快速回答“点在直线哪一侧”（叉积）、“两线段是否真正在屏幕上相交”（相交判断）、“一堆散点围成的最外圈轮廓长什么样”（凸包）、以及“如何从轮廓里找出最宽的那条缝”（旋转卡壳）。POJ P1113 Wall 是个经典靶子——给一堆钉子坐标，要求围一圈墙，墙到钉子的最小距离必须≥L，求墙的最短长度。表面看是几何题，实则考的是：凸包怎么建、怎么处理共线点、怎么把凸包周长和圆弧长度拼起来。我带过几个刚转图形方向的新人，他们第一次独立写出 Wall 的 AC 代码时，眼神里的光不是解出题的兴奋，而是“原来叉积真的能当开关用”的顿悟。

关键词里藏着真实痛点：“叉积”是底层原子操作，“相交判断”是高频业务需求，“凸包”是轮廓抽象的核心工具，“旋转卡壳”则是性能敏感场景的加速器。它们共同指向一个事实：在屏幕坐标系里，点不是点，是状态；线不是线，是边界；面不是面，是决策域。这篇内容适合三类人：正在啃 OJ 题但卡在几何模块的算法学习者、需要快速实现空间关系判断的前端/客户端工程师、以及想补足图形学基础但被数学符号劝退的跨领域开发者。它不承诺让你秒变计算几何专家，但能确保你下次看到“判断点是否在多边形内”时，第一反应不是百度，而是打开编辑器，敲下 cross_product(p, a, b) > 0。

2. 核心设计思路：从“纸面公式”到“机器可执行”的三重降维

计算几何最容易踩的坑，是把教科书当 API 文档用。比如叉积公式：a × b = |a||b|sinθ，这在纸上画图很美，但扔进代码里就是灾难——浮点误差会让 sinθ 在临界点疯狂抖动，角度计算本身又慢又不必要。真正的工程化思路，是彻底抛弃角度和模长，只保留符号意义和坐标运算。我们重新定义叉积为一个二维向量的标量值：
cross(o, a, b) = (a.x - o.x) * (b.y - o.y) - (a.y - o.y) * (b.x - o.x)
这个表达式没有三角函数，只有四则运算；它的正负号直接对应“从 oa 到 ob 是逆时针还是顺时针”；它的零值严格对应三点共线。这才是机器真正能稳定执行的“叉积”。

相交判断同理。教科书喜欢用参数方程联立求解，但实际项目中，你永远要面对“端点重合算不算相交”“T型连接怎么处理”“浮点精度下如何避免误判”这些魔鬼细节。我们的方案是分层判断：先用叉积做快速排斥试验（如果两线段包围盒不重叠，直接返回 false），再用跨立试验（线段 ab 的两端点 c、d 必须在直线 ab 的两侧，且 cd 的两端点 a、b 也必须在直线 cd 的两侧），最后对共线情况做投影区间检查（把点坐标投影到线段方向轴上，看区间是否重叠）。三层过滤下来，既保证正确性，又兼顾效率——跨立试验失败率超 90%，大部分判断在第一步就结束了。

凸包算法选 Andrew 单调链而非 Graham 扫描，原因很实在：Graham 需要极角排序，涉及 atan2 计算，不仅慢还有精度陷阱；Andrew 只需按 x 坐标排序（x 相同时按 y），然后两次线性扫描，全程整数比较+叉积，稳定得像呼吸。而 POJ P1113 Wall 的关键，在于理解“墙的最短长度 = 凸包周长 + 一个完整圆周长”。为什么？因为墙必须离每个钉子至少 L 距离，等效于把每个钉子膨胀成半径为 L 的圆，求这些圆的并集的最小外接凸包——而这个外接凸包的直线段部分，恰好就是原钉子集的凸包；弯曲部分，则是所有圆在凸包顶点处的外切圆弧，所有圆弧拼起来正好是一个整圆。这个结论不是靠空间想象，而是用凸包顶点处的内外角差推导出来的：每个顶点处的外角 = π - 内角，所有外角和恒为 2π，所以圆弧总角度就是 2π。

旋转卡壳的本质，是利用凸包的单峰性（antipodal pairs 的距离函数是单峰的）来避免 O(n²) 暴力枚举。但直接实现“卡壳”容易陷入指针迷宫。我们的简化模型是：固定一条边作为底边，找离它最远的顶点（用叉积绝对值衡量高度），这个顶点就是当前边的对踵点；然后顺时针旋转底边到下一条边，对踵点也顺时针移动（不会回退），用一个 while 循环维护。整个过程就像用两把卡尺夹住凸包，一边固定，一边滑动寻找最大宽度。这种思路比“双指针同步旋转”更符合直觉，也更容易调试。

3. 核心细节解析与实操要点：那些教科书绝不会写的“手抖时刻”

3.1 叉积的符号陷阱与精度围栏

叉积的正负号是计算几何的命脉，但也是最易翻车的点。新手常犯的错误是：写cross(a, b, c) > 0判断 c 在 ab 左侧，却忘了 a、b、c 的顺序决定方向。叉积 cross(o, a, b) 的物理意义是：以 o 为原点，向量 oa 到 ob 的有向面积。所以cross(o, a, b) > 0表示从 oa 到 ob 是逆时针旋转，即 b 在 oa 的左侧；而cross(o, b, a) > 0就是相反方向。我在某次 GIS 项目中就栽在这儿——判断用户点击点是否在河流多边形左侧时，把点序写反了，结果所有判断全颠倒，测试时用鼠标划拉半天才意识到是符号逻辑错了。

更大的坑是浮点精度。当三点几乎共线时，叉积值可能是一个极小的正数或负数（如 1e-15），直接比较> 0会因舍入误差导致错误。解决方案是引入精度容差 eps，但 eps 不是随便设的。经验法则是：eps 应大于max(|x|,|y|) * 1e-12（相对误差）和1e-9（绝对误差）中的较大者。例如，若坐标范围在 [-1000, 1000]，则 eps 至少取1000 * 1e-12 = 1e-9，所以统一设const double eps = 1e-9是安全的。判断逻辑必须写成：

int dcmp(double x) { if (fabs(x) < eps) return 0; return x < 0 ? -1 : 1; } // 使用时：if (dcmp(cross(o, a, b)) > 0) // 正确！

dcmp函数返回 -1/0/1，彻底规避了浮点比较的不确定性。我见过太多人用abs(cross) < eps做共线判断，结果在边界 case 上反复失败——因为cross本身可能因精度丢失而产生虚假的非零值，dcmp才是唯一可靠的守门员。

3.2 相交判断的“四步法”与 T 型连接特例

线段相交判断必须拆解为四个明确步骤，缺一不可：

1. 快速排斥试验（Bounding Box Check）：
   若max(a.x,b.x) < min(c.x,d.x) || max(c.x,d.x) < min(a.x,b.x) || max(a.y,b.y) < min(c.y,d.y) || max(c.y,d.y) < min(a.y,b.y)，则无交。这是最廉价的剪枝，用四次比较干掉 80% 的无效 case。

2. 跨立试验（Straddle Test）：
   计算d1 = dcmp(cross(a,b,c)),d2 = dcmp(cross(a,b,d)),d3 = dcmp(cross(c,d,a)),d4 = dcmp(cross(c,d,b))。
   相交的充要条件是：(d1 * d2 < 0 && d3 * d4 < 0)—— 严格跨立。
   但业务中常需包含端点接触，此时条件变为：(d1 * d2 <= 0 && d3 * d4 <= 0)。

3. 共线点特判（Collinear Endpoint Handling）：
   当d1 == 0 || d2 == 0 || d3 == 0 || d4 == 0时，说明至少有一个端点在线段上。此时需做投影检查：例如d1 == 0表示 c 在 ab 上，需验证 c 的 x 坐标在[min(a.x,b.x), max(a.x,b.x)]且 y 坐标在[min(a.y,b.y), max(a.y,b.y)]内（注意用<=包含端点）。

4. T 型连接的业务取舍（The T-Junction Dilemma）：
   当d1 == 0 && d2 != 0 && d3 != 0 && d4 != 0时，是标准 T 型（c 是顶点，ab 是横杠）。是否算相交？没有标准答案，取决于业务。地图路网中，T 型路口是合法连接，应返回 true；而 CAD 中，T 型可能表示建模错误，需告警。我的做法是：封装为segment_intersect_strict()和segment_intersect_inclusive()两个函数，让调用者根据场景选择。

提示：所有坐标比较必须用dcmp，包括投影检查中的min/max边界判断。我曾在一个实时渲染项目中，因投影检查用了原始<导致 T 型连接在特定帧率下时有时无，排查三天才发现是浮点精度漏网。

3.3 凸包构建的 Andrew 算法：排序、去重、单调链

Andrew 算法的优雅在于其线性时间复杂度和极简逻辑。核心是构建下凸壳（从最左点到最右点，下半部分）和上凸壳（从最右点到最左点，上半部分），然后拼接。

排序阶段：按x升序，x相同时按y升序。这确保了点列从左到右排列，为单调链提供基础。注意：必须用stable_sort或自定义比较函数，避免x相同时y乱序。

去重阶段：排序后相邻点可能完全重合。必须在构建凸壳前剔除重复点，否则cross计算会得到 0，导致算法崩溃。简单循环即可：

vector<Point> unique_points; for (int i = 0; i < points.size(); i++) { if (i == 0 || dcmp(points[i].x - points[i-1].x) != 0 || dcmp(points[i].y - points[i-1].y) != 0) unique_points.push_back(points[i]); }

单调链构建：以下凸壳为例。初始化栈hull，压入前两点。对后续每个点p：

• 当栈中元素 ≥2 且cross(hull[second_top], hull[top], p) <= 0时（即新点使栈顶角内凹），弹出栈顶；
• 将p压入栈。
  关键细节：<= 0是为了剔除共线点（三点一线时，中间点不是凸包顶点）；若需保留共线点（如 Wall 题要求），则改为< 0。上凸壳同理，只是遍历顺序从右到左，且叉积条件符号相反（因方向反转）。

注意：Andrew 算法输出的凸包顶点是逆时针顺序，且首尾点不重复。这点至关重要——旋转卡壳和面积计算都依赖此顺序。我曾因忘记清空栈或搞错上下凸壳的叉积符号，导致凸包输出一团乱麻，调试时打印所有中间点坐标才定位到栈操作逻辑错误。

3.4 POJ P1113 Wall 的“膨胀-收缩”本质与圆弧处理

P1113 的题干看似是几何构造，实则是凸包应用的教科书级案例。其核心洞察是：满足“墙到任意钉子距离 ≥ L” 的最短闭合曲线，等价于将所有钉子视为半径为 L 的圆，求这些圆的并集的最小外接凸包的周长。

这个“外接凸包”的构成有两部分：

• 直线段部分：恰好是原始钉子集的凸包的每条边，向外平移 L 距离后的平行线段。由于平移距离相同，这些线段长度与原凸包边完全一致。
• 圆弧部分：在凸包的每个顶点处，由两个相邻平移边的延长线相交形成一个外角，该外角所对应的圆弧，正是半径为 L 的圆上的一段。

关键推导：设凸包有 n 个顶点，第 i 个顶点的内角为 α_i，则其外角为 π - α_i。所有外角之和恒为 2π（凸多边形性质）。因此，所有圆弧的圆心角之和 = Σ(π - α_i) = nπ - Σα_i = nπ - (n-2)π = 2π。这意味着，无论凸包形状如何，所有圆弧拼起来，恰好是一个完整的圆周，长度为 2πL。

所以最终墙长 = 凸包周长 + 2πL。
实操中，只需：

1. 对钉子坐标构建 Andrew 凸包；
2. 计算凸包周长（相邻顶点距离累加）；
3. 加上2 * PI * L。

但魔鬼在细节：Wall 题要求凸包必须是严格凸的，即不能有三点共线。Andrew 算法中若用< 0判断，则共线点会被剔除，符合要求；若用<= 0，则共线点保留，会导致周长计算偏大（因为共线点间的线段被计入，但实际墙在此处是直线，无需额外圆弧）。因此，Wall 必须用< 0版本的 Andrew。

实操心得：在提交 POJ 时，务必用double存储坐标和结果，printf("%.0f\n", ans)输出整数。我曾因用int存坐标导致中间计算溢出，WA 到怀疑人生；也因%.0f在 .5 时四舍五入规则不同（C++ 默认 round half to even），改用%.0f并确保输入数据无精度争议才 AC。

4. 实操过程与核心环节实现：从零开始跑通 POJ P1113 Wall

4.1 完整代码框架与数据结构定义

我们从最基础的Point结构体开始，它必须支持叉积、距离、输入输出：

#include <iostream> #include <vector> #include <algorithm> #include <cmath> #include <iomanip> using namespace std; const double eps = 1e-9; const double PI = acos(-1.0); struct Point { double x, y; Point(double x = 0, double y = 0) : x(x), y(y) {} bool operator < (const Point& p) const { return dcmp(x - p.x) < 0 || (dcmp(x - p.x) == 0 && dcmp(y - p.y) < 0); } bool operator == (const Point& p) const { return dcmp(x - p.x) == 0 && dcmp(y - p.y) == 0; } }; int dcmp(double x) { if (fabs(x) < eps) return 0; return x < 0 ? -1 : 1; } double cross(const Point& o, const Point& a, const Point& b) { return (a.x - o.x) * (b.y - o.y) - (a.y - o.y) * (b.x - o.x); } double dist(const Point& a, const Point& b) { return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y)); }

Point的operator<严格按dcmp比较，确保sort稳定；operator==同理。cross和dist是所有算法的基石。

4.2 Andrew 凸包算法的逐行实现与注释

凸包函数convex_hull返回逆时针顺序的顶点向量（首尾不重复）：

vector<Point> convex_hull(vector<Point>& points) { int n = points.size(); if (n <= 1) return points; // Step 1: Sort and remove duplicates sort(points.begin(), points.end()); vector<Point> unique_p; for (int i = 0; i < n; i++) { if (i == 0 || points[i] != points[i-1]) unique_p.push_back(points[i]); } n = unique_p.size(); if (n <= 2) return unique_p; // Step 2: Build lower hull vector<Point> hull; for (int i = 0; i < n; i++) { while (hull.size() >= 2 && dcmp(cross(hull[hull.size()-2], hull.back(), unique_p[i])) <= 0) { hull.pop_back(); } hull.push_back(unique_p[i]); } // Step 3: Build upper hull int k = hull.size(); // lower hull size for (int i = n-2; i >= 0; i--) { // from second last to first while (hull.size() > k && dcmp(cross(hull[hull.size()-2], hull.back(), unique_p[i])) <= 0) { hull.pop_back(); } hull.push_back(unique_p[i]); } // Step 4: Remove duplicate first/last point (hull[0] == hull.back()) if (hull.size() > 1) hull.pop_back(); return hull; }

关键注释：

• hull.size() >= 2是安全检查，防止访问越界；
• 下凸壳遍历i=0到n-1（从左到右），上凸壳遍历i=n-2到0（从右到左，跳过已加入的最右点）；
• k = hull.size()记录下凸壳结束位置，确保上凸壳只在新增点上操作；
• 最后pop_back()是因为上凸壳会把起点（最左点）再次加入，造成重复。

4.3 POJ P1113 Wall 的主逻辑与精度控制

主函数读入数据，调用凸包，计算周长，加上圆周：

int main() { int n, L; while (cin >> n >> L) { vector<Point> points(n); for (int i = 0; i < n; i++) { cin >> points[i].x >> points[i].y; } vector<Point> hull = convex_hull(points); int m = hull.size(); // Calculate convex hull perimeter double perimeter = 0.0; for (int i = 0; i < m; i++) { int j = (i + 1) % m; perimeter += dist(hull[i], hull[j]); } // Add the circle circumference: 2 * PI * L double ans = perimeter + 2.0 * PI * L; // Output as integer (round to nearest) cout << fixed << setprecision(0) << ans << endl; } return 0; }

这里setprecision(0)配合fixed，确保ans四舍五入到整数位。POJ 测试数据保证结果为整数，所以%.0f是安全的。

4.4 旋转卡壳求凸包直径的实战封装

虽然 Wall 不需要，但旋转卡壳是凸包的标配技能。我们封装一个rotating_calipers_diameter函数：

double rotating_calipers_diameter(const vector<Point>& hull) { int n = hull.size(); if (n == 1) return 0.0; if (n == 2) return dist(hull[0], hull[1]); double max_dist = 0.0; int j = 1; // antipodal point for i=0 // For each edge i->i+1 for (int i = 0; i < n; i++) { int next_i = (i + 1) % n; // Move j until distance from hull[j] to line hull[i]->hull[next_i] is maximized // Use cross product absolute value as "height" while (true) { int next_j = (j + 1) % n; double curr_h = fabs(cross(hull[i], hull[next_i], hull[j])); double next_h = fabs(cross(hull[i], hull[next_i], hull[next_j])); if (dcmp(next_h - curr_h) > 0) { j = next_j; } else { break; } } // Update max distance with current antipodal pair max_dist = max(max_dist, dist(hull[i], hull[j])); max_dist = max(max_dist, dist(hull[next_i], hull[j])); } return max_dist; }

逻辑说明：

• j是当前边i->next_i的对踵点；
• while循环中，比较hull[j]和hull[next_j]到当前边的“有向高度”（fabs(cross)），取更大的那个；
• 每次更新max_dist时，不仅要比较hull[i]和hull[j]，还要比较hull[next_i]和hull[j]（因为对踵点可能更靠近边的某一端点）；
• 时间复杂度 O(n)，因为j总共只循环一圈。

实测心得：在某次游戏开发中，用此函数实时计算角色碰撞体的“最大跨度”，替代了 O(n²) 的暴力枚举，CPU 占用从 12% 降到 1.5%。但要注意，此函数假设hull是逆时针且无重复点，必须用convex_hull的输出。

5. 常见问题与排查技巧实录：那些让我熬夜到凌晨三点的 Bug

5.1 “凸包输出为空”——排序与去重的隐形杀手

现象：输入有效点集，convex_hull返回空向量或只有 1 个点。
排查路径：

1. 检查输入：用cout打印原始points，确认坐标读入无误（常见错误：cin失败导致x,y为 0）；
2. 检查排序：打印排序后的points，看是否按x升序排列；若x全相同，是否按y升序；
3. 检查去重：打印unique_p大小，若为 0，说明所有点被判定为重复——大概率是operator==未用dcmp，导致1.0000000001 == 1.0返回 false，去重逻辑失效。

解决方案：在operator==中强制使用dcmp：

bool operator == (const Point& p) const { return dcmp(x - p.x) == 0 && dcmp(y - p.y) == 0; }

5.2 “相交判断总是 false”——跨立试验的符号迷宫

现象：明明两线段交叉，segment_intersect却返回 false。
最可能原因：跨立试验中d1 * d2 <= 0的<=写成了<，导致端点接触被忽略。
更隐蔽的原因：叉积参数顺序错误。例如判断ab和cd是否相交，错误地写了cross(a,b,c)和cross(a,b,d)，但cross(c,d,a)却写成cross(a,c,d)——参数o,a,b的顺序决定了参考系。

快速验证法：手动计算一个简单 case。取a(0,0), b(2,0), c(1,1), d(1,-1)，显然相交于(1,0)。计算：

• cross(a,b,c) = (1-0)*(1-0) - (1-0)*(2-0) = 1 - 2 = -1
• cross(a,b,d) = (1-0)*(-1-0) - (-1-0)*(2-0) = -1 + 2 = 1
• d1*d2 = -1*1 = -1 < 0✓
• cross(c,d,a) = (0-1)*(0-(-1)) - (0-(-1))*(1-1) = (-1)*(1) - (1)*(0) = -1
• cross(c,d,b) = (2-1)*(0-(-1)) - (0-(-1))*(1-1) = 1*1 - 1*0 = 1
• d3*d4 = -1*1 = -1 < 0✓
  全部通过。若任一计算不符，立刻检查cross函数实现。

5.3 “Wall 答案偏大”——共线点与圆弧的误解

现象：P1113 输出比标准答案大几单位。
根源几乎一定是凸包包含了共线点。例如三点(0,0), (1,0), (2,0)，若 Andrew 算法用<= 0，则(1,0)会被保留在凸包中，导致周长多算dist((0,0),(1,0)) + dist((1,0),(2,0)) = 2，而正确凸包只有(0,0)和(2,0)，边长为 2，但此时dist((0,0),(2,0)) = 2，周长相同？不！问题在圆弧：三点共线时，中间点(1,0)处的“外角”为 π，对应圆弧长度为 πL，而两端点处的外角各为 π/2，圆弧长各为 (π/2)L，总和仍是 2πL。但 Wall 的题意要求“墙必须离每个钉子 ≥ L”，三点共线时，中间钉子(1,0)的“安全区”已被两端点的墙覆盖，无需额外圆弧。因此，共线点必须剔除，只保留端点。

验证方法：对凸包输出hull，遍历每三个连续点hull[i-1], hull[i], hull[i+1]，计算dcmp(cross(hull[i-1], hull[i], hull[i+1]))，若为 0，则存在共线点，算法有误。

5.4 “旋转卡壳结果震荡”——对踵点移动的边界条件

现象：rotating_calipers_diameter在不同运行中返回微小差异的值（如10.0000001vs9.9999999）。
原因：while循环中，当next_h和curr_h极其接近时，dcmp(next_h - curr_h) > 0可能因精度误差返回 false，导致j停在次优位置。

解决方案：在比较前加入一个微小的“推进阈值”：

double diff = next_h - curr_h; if (dcmp(diff) > 0 || (dcmp(diff) == 0 && dcmp(dist(hull[i], hull[next_j]) - dist(hull[i], hull[j])) > 0)) { j = next_j; }

即：当高度相等时，选择距离更大的点作为对踵点，确保结果稳定。这个技巧在某次金融风控系统的空间聚类中救了急——客户要求“最远两点距离”必须幂等，否则触发告警。

5.5 综合排查速查表

我的终极建议：在项目根目录建一个geometry_test.cpp，里面塞满上述所有单元测试 case（三点共线、T 型、交叉、凸包退化等）。每次修改几何函数，先跑一遍测试。这比对着 POJ WA 列表猜原因快十倍。计算几何的确定性，就藏在这些看似琐碎的dcmp和<= 0里。