UVa 137 Polygons这道题是我刷计算几何专题时印象很深的一道题。题面一句话——给你两个凸多边形求它们重叠区域的面积——看起来人畜无害实际上手全是细节顶点顺序乱给、叉积方向搞反、浮点精度翻车、输出格式被卡任何一个坑踩进去都够你调试一个晚上。这篇文章我会把完整推导、可提交的C代码和我在踩坑过程中的实际经验一次性写清楚适合正在刷UVa/ACM专题的竞赛选手也适合GIS叠加分析、图形裁剪、2D碰撞检测这类工作中需要处理多边形相交的开发同学参考。1. 两分钟看懂题目两个凸多边形的交集面积到底怎么算1.1 题意还原与解题突破口UVa 137的输入结构很直接每行先给一个点数nn0表示输入结束接着给这个多边形的n个顶点坐标然后给第二个多边形的点数m和m个顶点坐标。两个多边形都保证是凸的顶点按顺序给出顺时针或逆时针都有可能。要求输出的核心是交集面积也就是两个多边形共同覆盖的那部分有多大。突破口其实就一句话既然两个多边形都是凸的交集一定还是凸多边形。这意味着你不需要处理带洞区域也不需要管凹多边形裁剪的复杂情况。拿到这个结论后解题路径就清晰了——找到交集多边形的顶点集合然后用鞋带公式算面积。而找到交集多边形的顶点这件事最经典的做法就是Sutherland-Hodgman多边形裁剪算法。1.2 为什么选择Sutherland-Hodgman而不是半平面交模板很多同学看到求凸多边形交集第一反应是上通用的半平面交算法。半平面交确实能解但那是把问题复杂化了。半平面交需要维护双端队列、判平行、判交点顺序代码量接近一百行调试成本高。而Sutherland-Hodgman的思路简单得多用第二个多边形的每条边当作一把刀按顺序去切第一个多边形切剩下的部分就是交集。这个算法原本是图形学里做视口裁剪的经典方案后来被广泛用到多边形求交问题里。它要求裁剪多边形是凸的UVa 137恰好满足所以是完美匹配。算法复杂度是O(n*m)n和m是两个多边形的顶点数本题量级下完全无所谓。如果两个多边形顶点数都到了上万级别再考虑优化也不迟后面我会单独讲。1.3 面积关系交集、并集、非重叠面积一眼分清UVa 137的输出里通常涉及两个面积数字一个是交集面积另一个是某种差集或并集相关面积。这里先把面积关系理清楚后面编码不会懵交集面积 I 两个多边形共同部分并集面积 U areaA areaB - I非重叠面积 D areaA areaB - 2 * I两个多边形各自独有部分之和我在代码里会把这三个量全部算出来输出时按自己手上的题面要求挑选排列。至于具体输出哪个在前哪个在后不同版本OJ的题面措辞可能不同以题面PDF为准。后面代码部分的写法可以很方便地调换。2. 动手前的三块基石叉积、方向、点在凸多边形内2.1 叉积是计算几何的瑞士军刀计算几何里叉积是一个绕不开的概念。二维向量a(ax, ay)b(bx, by)它们的叉积定义为cross(a, b) ax * by - ay * bx这个数值的绝对值等于两个向量围成的平行四边形面积符号则表示b相对a的旋转方向cross 0说明b在a的逆时针方向cross 0说明b在a的顺时针方向等于0说明共线。多边形裁剪中叉积用在了三个地方算多边形面积、判断点在多边形的哪一侧、求两条线段交点。理解了叉积的几何意义后面所有代码都是顺手的事。我习惯把它理解成方向的裁判——它告诉你在某个参照向量看来另一个向量是往左偏还是往右偏。2.2 鞋带公式算面积顺带统一多边形方向已知多边形所有顶点按顺序给出面积不需要三角剖分直接用鞋带公式S 0.5 * | Σ (x_i * y_(i1) - x_(i1) * y_i) |求和里的每一项其实就是一个叉积cross(p_i, p_(i1))。代码实现double directed_area(const Polygon poly) { double s 0; int n poly.size(); for (int i 0; i n; i) { int j (i 1) % n; s cross(poly[i], poly[j]); } return s / 2; } double poly_area(const Polygon poly) { return fabs(directed_area(poly)); }模板里我留了个有向面积的版本它返回带符号的面积正数说明顶点是逆时针顺序负数说明是顺时针顺序。这个符号在Sutherland-Hodgman算法里至关重要因为点在多边形内侧的判断依赖方向。代码里直接用它来统一方向void make_ccw(Polygon poly) { if (directed_area(poly) 0) { reverse(poly.begin(), poly.end()); } }2.3 点在凸多边形内部判定利用叉积符号一致凸多边形有个非常优美的性质如果顶点按逆时针顺序排列那么多边形内部始终在每条有向边的左侧。因此判断点P是否在多边形内部只需要依次检查每条边a-b看看cross(b-a, P-a)是否都大于等于0。bool inside(const Point p, const Point a, const Point b) { return cross(b - a, p - a) -eps; }这里用 -eps而不是 0是为了容忍浮点误差。点恰好在边上的情况既算内也算外其实不影响最终面积但会把裁剪结果搞出一堆重复点所以统一认为在边上是内侧。我用一个简单例子验证正方形(0,0)-(1,0)-(1,1)-(0,1)逆时针。点(0.5, 0.2)依次检查四条边cross都是正的点(0.5, -0.1)第一条边cross(edge, p-a) 1 * (-0.1) - 0 * 0.5 -0.1负号一出直接判为外部。简单又高效。3. Sutherland-Hodgman裁剪算法拆解用第二个多边形的每条边去切第一个3.1 算法思想半平面求交Sutherland-Hodgman的核心思路可以这样类比假设你有一块黏土第一个多边形要把它修成某个模子第二个多边形的形状。模子的每条边都是一把无限长的刀拿这把刀沿着边所在直线切一刀把外侧的黏土切掉。四条边依次切下去剩下的黏土就是交集。更准确地说多边形的一条有向边a-b定义了一个半平面所有位于这条边左侧的点都在半平面内右侧的点在半平面外。拿这个半平面去裁剪当前多边形保留内侧部分。第二个多边形的每条边都做一次这样的裁剪最终留下的多边形就是两个多边形交集。注意这里用的是边所在的无限直线不是线段本身。因为裁剪多边形是凸的用无限半平面裁剪和用线段裁剪效果完全一致而且计算更简单。如果裁剪多边形是凹的这条就不成立了这是算法适用范围的边界。3.2 四种情况的边处理规则每一轮裁剪都要遍历当前多边形的每条边p-q判断p和q相对于裁剪线的内外状态。只有四种组合起点状态终点状态处理动作内侧内侧把终点q加入结果内侧外侧把p-q与裁剪线的交点加入结果外侧内侧先加入交点再加入终点q外侧外侧什么都不加为什么起点在内、终点在外时只加交点不加别的因为当前多边形的边从内侧穿到外侧中间必然与裁剪线相交交点就是多边形被切掉的边界外侧部分的边已经不属于保留区域不能加终点。反过来起点在外、终点在内时先加交点再加终点相当于从裁剪线上重新进入多边形内部。两个点都在内侧时整条边都在保留区域内加终点即可——因为起点在前一轮已经处理过了不需要重复加。3.3 交点计算公式线段p1-p2和裁剪线a-b所在直线求交设交点P p1 t * (p2 - p1)因为P在直线a-b上所以cross(b-a, P-a) 0代入解tcross(b-a, p1 t*(p2-p1) - a) 0 t cross(b-a, a-p1) / cross(b-a, p2-p1)代码实现Point intersect(const Point p1, const Point p2, const Point a, const Point b) { double t cross(b - a, a - p1) / cross(b - a, p2 - p1); return p1 (p2 - p1) * t; }分母是cross(b-a, p2-p1)如果两条线段平行分母为0。但在Sutherland-Hodgman的正常流程里出现跨越边界的分支就意味着p1和p2在裁剪线的两侧这两条线段不可能平行所以分母不会为0。不过为了绝对安全还是建议判断一下fabs(denominator) eps的情况我实际测试中没触发过但防御性写代码没有坏处。3.4 完整流程与状态变化外层循环遍历裁剪多边形的每一条边内层循环遍历当前多边形的每一条边。每轮裁剪都生成一个新的顶点序列作为下一轮裁剪的输入。Polygon clip_polygon(const Polygon subject, const Polygon clipper) { Polygon cur subject; int m clipper.size(); for (int i 0; i m; i) { const Point a clipper[i]; const Point b clipper[(i 1) % m]; Polygon nxt; int n cur.size(); if (n 0) break; for (int j 0; j n; j) { const Point p cur[j]; const Point q cur[(j 1) % n]; bool pin inside(p, a, b); bool qin inside(q, a, b); if (pin qin) { nxt.push_back(q); } else if (pin !qin) { nxt.push_back(intersect(p, q, a, b)); } else if (!pin qin) { nxt.push_back(intersect(p, q, a, b)); nxt.push_back(q); } } cur nxt; } return cur; }注意每轮裁剪开始时我会检查cur是否已经为空。如果当前多边形被某条裁剪线完全切掉了说明两个多边形没有任何重叠可以直接返回空没必要继续循环。4. 可直接提交的C实现从读入到输出的完整代码4.1 数据结构设计多边形就用vector 表示Point封装x和y坐标重载减法、加法和数乘运算符这样求交点、算叉积的表达式会干净很多。UVa的输入坐标可能是整数也可能是浮点数统一用double读。4.2 核心函数列表整个代码需要五个核心部分点的加法、减法、数乘运算叉积函数两向量叉积、三点叉积略过不写也行有向面积 统一方向inside判断clip_polygon裁剪函数主函数读入两个多边形后先make_ccw统一方向再调用裁剪函数最后算面积输出。4.3 完整代码#include cstdio #include cmath #include vector #include algorithm using namespace std; const double eps 1e-9; struct Point { double x, y; Point() {} Point(double x, double y) : x(x), y(y) {} Point operator (const Point p) const { return Point(x p.x, y p.y); } Point operator - (const Point p) const { return Point(x - p.x, y - p.y); } Point operator * (double t) const { return Point(x * t, y * t); } }; typedef vectorPoint Polygon; double cross(Point a, Point b) { return a.x * b.y - a.y * b.x; } double directed_area(const Polygon poly) { double s 0; int n poly.size(); for (int i 0; i n; i) { int j (i 1) % n; s cross(poly[i], poly[j]); } return s / 2; } double poly_area(const Polygon poly) { return fabs(directed_area(poly)); } void make_ccw(Polygon poly) { if (directed_area(poly) 0) { reverse(poly.begin(), poly.end()); } } bool inside(const Point p, const Point a, const Point b) { return cross(b - a, p - a) -eps; } Point intersect(const Point p1, const Point p2, const Point a, const Point b) { double t cross(b - a, a - p1) / cross(b - a, p2 - p1); return p1 (p2 - p1) * t; } Polygon clip_polygon(const Polygon subject, const Polygon clipper) { Polygon cur subject; int m clipper.size(); for (int i 0; i m; i) { const Point a clipper[i]; const Point b clipper[(i 1) % m]; Polygon nxt; int n cur.size(); if (n 0) break; for (int j 0; j n; j) { const Point p cur[j]; const Point q cur[(j 1) % n]; bool pin inside(p, a, b); bool qin inside(q, a, b); if (pin qin) { nxt.push_back(q); } else if (pin !qin) { nxt.push_back(intersect(p, q, a, b)); } else if (!pin qin) { nxt.push_back(intersect(p, q, a, b)); nxt.push_back(q); } } cur nxt; } return cur; } int main() { int n; while (scanf(%d, n) 1 n) { Polygon A(n); for (int i 0; i n; i) { scanf(%lf%lf, A[i].x, A[i].y); } int m; scanf(%d, m); Polygon B(m); for (int i 0; i m; i) { scanf(%lf%lf, B[i].x, B[i].y); } make_ccw(A); make_ccw(B); Polygon inter clip_polygon(A, B); double areaA poly_area(A); double areaB poly_area(B); double I poly_area(inter); double non_overlap areaA areaB - 2 * I; if (fabs(non_overlap) 1e-6) non_overlap 0; if (fabs(I) 1e-6) I 0; // 这里输出顺序按你手上题面要求调换即可 printf(%7.2f%7.2f\n, non_overlap, I); } return 0; }4.4 用两个正方形验证正确性写题解最怕代码看着对但样例过不了。我用一个手工可算的例子来验证。正方形A(0,0)-(1,0)-(1,1)-(0,1)正方形B(0.5,0)-(1.5,0)-(1.5,1)-(0.5,1)。两个正方形的交集是(0.5,0)-(1,0)-(1,1)-(0.5,1)面积0.5。用一个简化思路来人工模拟一遍以B的最后一条边(0.5,1)-(0.5,0)裁剪A时这条边是竖直向上的内侧在左边也就是x0.5。A的顶点(0,0)和(0,1)都在x0.5一侧会被切掉。第一次交点出现在(0,0)-(1,0)这条底边与x0.5的交点即(0.5,0)第二次交点出现在(0,1)-(1,1)这条顶边与x0.5的交点即(0.5,1)。裁剪后得到(0.5,0)-(1,0)-(1,1)-(0.5,1)面积用鞋带公式一算就是0.5。这个例子放进代码输出就是预期的两个数字。我建议拿到模板后先构造几个自己口算得出来的数据跑一遍确认AC之后再提交。常见测试包括完全包含、完全分离、仅触碰顶点、两个完全相同的多边形。全部通过后这个问题就稳了。5. 提交前必须绕开的三个坑方向、精度、输出格式5.1 坑一多边形顶点方向不统一裁剪结果直接反掉这是新手最容易踩的坑。inside函数里的cross(b-a, p-a) -eps依赖裁剪多边形是逆时针方向。如果输入是顺时针原本在左侧的点其实在多边形外部裁剪出来要么是空要么是一堆乱七八糟的顶点。统一方向的方法前面已经说了用有向面积判断符号负则reverse。但要注意两个多边形要分别判断、分别调整不能默认它们方向一致。UVa的输入数据有时候两个都是逆时针有时候两个都是顺时针也可能一个顺一个逆。我在第一次写这道题时默认了输入都是逆时针结果第二个多边形是顺时针样例直接挂了。调试时如果发现裁剪结果面积全为0或者面积比实际小先检查方向。5.2 坑二浮点误差和-0.00问题坐标是浮点数叉积计算会累积浮点误差。eps取1e-9是我在UVa上的经验值太小起不到容错作用太大会把明显不重合的点判成重合。面积计算完以后fabs小于1e-6的可以直接置0避免输出-0.00这种让人抓狂的结果。另外交点计算里的除法也要注意。虽然跨越边界时分母理论上不为0但如果两个多边形存在共线边某些极端case里分母可能非常接近0导致交点坐标变成天文数字。我自己遇到过一次两个多边形的边刚好平行且几乎重叠交点计算溢出最终面积变成了一亿多。后来我在intersect函数里加了分母判断小于eps直接返回p1问题就消失了。稳妥起见这个保护值得加。5.3 坑三输出格式精确到像素级UVa老题的输出格式非常严格。常见的格式是每个数字占7位右对齐保留两位小数两个数字之间没有空格。我建议先直接用printf(%7.2f%7.2f\n, non_overlap, I)去交如果Presentation Error再根据题面提示调整空格。很多人一上来自作聪明加了一个空格结果PE到怀疑人生。这算是我用血泪换来的教训。还需要注意输出顺序。我个人代码里先输出non_overlap再输出I如果你的题面要求交集面积在前把printf里的两个参数换一下就行。因为输出的两个数在大部分版本里就是非重叠区域面积和重叠区域面积这两个语义具体顺序保持一致即可。5.4 复杂度与进一步优化Sutherland-Hodgman的复杂度是O(n*m)每次裁剪都要遍历当前多边形的所有顶点而当前多边形顶点数最坏情况会每一轮增加。UVa 137的数据规模很小这个复杂度毫无压力。但如果换到数据规模大的场景比如两个多边形各有十万个顶点这个方法就扛不住了。优化方向有两种一是改成标准半平面交把每个凸多边形看作若干半平面的交集一次性求交集复杂度O((nm) log(nm))但代码复杂度也上来了二是用平面扫描法求多边形交处理凹多边形和含洞多边形时更通用。竞赛里遇到更大的数据我通常直接选择半平面交模板。不过要是只是实战项目里做一次性的多边形叠加分析Sutherland-Hodgman这种简单粗暴的方案反而是最不容易出错的工程上的可维护性有时候比理论复杂度更重要。我个人的习惯是上模板之前先写一个基于随机生成凸多边形的对拍脚本用暴力方法验证裁剪面积是否正确。方向问题和精度问题在这种暴力对拍下几乎当场就会暴露。UVa 137的代码虽然不长但把这里的每个细节都抠透了后面再碰半平面交、多边形核、Voronoi图这些专题时底子就算打扎实了。