简单多边形O(n)凸包算法：Melkman原理与C++实现

做计算几何的人应该都有这种经历拿到一个多边形轮廓第一反应就是用Graham扫描或者Andrew单调链去求凸包写起来不复杂但总觉得排序那一步有点浪费。直到我认真去翻Melkmans Algorithm才意识到在“顶点按顺序给出的简单多边形”这个特定前提下我们完全可以丢掉排序用O(n)的时间直接扫一遍出结果。这个算法不算冷门但在很多工程代码里反而见得少可能是因为教材里普遍只讲Graham和Andrew提到Melkman时也常常一笔带过。这篇文章我会从算法思路、几何原理到完整C实现讲透再把我实际调试中踩过的坑和排查经验一起整理出来。适合正在写几何算法、游戏碰撞体生成、图像轮廓分析或者地理围栏相关功能的人看完可以直接照着用。1. 为什么需要“多边形快速凸包算法”1.1 凸包到底在解决什么问题凸包在图形学、游戏、路径规划、图像处理里都太常见了。说人话就是给你一堆点或者一个多边形用一根无形的“橡皮筋”把所有点包住绷紧之后得到的那个凸多边形就是凸包。实际工程里我遇到最多的是两类需求。一类是碰撞检测引擎里如果用原始多边形做精确碰撞顶点多、计算量大把复杂凹多边形先转成凸包再套用凸多边形的分离轴定理速度会快很多。另一类是轮廓分析比如拿到图像分割得到的区域轮廓想获取它的外边界大致形状凸包就是最稳的第一步。还有一些场景比如障碍物膨胀、电子围栏生成、活动范围分析本质都是求一组有序顶点围成区域的凸包。这些场景有一个共同特征输入不是乱序散点而是多边形顶点按顺序给好的轮廓。1.2 已有凸包算法的主要瓶颈常见的凸包算法里Graham扫描和Andrew单调链都把排序作为前置步骤时间复杂度O(n log n)。排序对于“任意点集”来说无法避免因为点之间本来没有顺序关系。Jarvis步进法虽然不用排序但复杂度是O(n·h)h是凸包顶点个数最坏情况退化成O(n²)。这就在实际工程里造成一个奇怪的现象明明轮廓数据已经是按顺序给出来的多边形顶点却还要额外进行一次极角排序或坐标排序白白浪费了输入里已经包含的顺序信息。1.3 什么输入才适合MelkmanMelkman算法有一个明确的前提输入必须是简单多边形的顶点序列按顶点顺序顺时针或逆时针均可给出。简单多边形就是边与边不会自交的多边形。只要满足这个前提算法就能达到O(n)不需要排序也不依赖凸包顶点个数。如果拿到的输入是乱序散点那Melkman就不适用老老实实去排序做Andrew。这是选型时最重要的判断点。2. 算法为什么能O(n)从双端队列说起2.1 巧妙的“在线处理”思路Melkman处理顶点的方式是流式的从多边形的第一个顶点开始逐个读入每一步都维护“当前已读顶点集合”的凸包。你不需要等到所有点都读进来每处理一个点手里就已经有一个正确的凸包。这一点很像在线算法。你在一条折线上往前走手里的凸包是动态更新的。因为这个特性它甚至可以用在流式轮廓输入的场景里比如传感器一个点一个点往外吐数据的情况。它不需要排序的根本原因就是顶点顺序已经固定了。对于任意点集极角排序能帮我们确定谁先谁后但在多边形轮廓里先后关系天然存在算法只需要沿着这条顺序走下去即可。2.2 双端队列里到底存了什么实现的核心结构是双端队列deque。deque里维护的是当前凸包的一条闭合链front元素和back元素是同一个点相当于一个环形链首尾相接。举个例子初始化完成之后deque可能是这样d [p2, p0, p1, p2]这里的d.front()和d.back()都是p2从front到back遍历得到的顶点序列p2 → p0 → p1 → p2构成了一个闭合的凸多边形边界。每次来一个新顶点先把deque两端不合适的旧顶点弹出再把新顶点同时插入两端。插入之后front和back又变成同一个点都是新顶点闭合链继续成立。2.3 弹出规则背后的几何直觉处理每个新顶点时有两段循环判断。第一段看deque头部方向判断当前新顶点相对于边d[0] → d[1]的位置。如果新顶点在这条边的右侧或者与这条边共线就把d[0]从头部弹出一直弹到新顶点位于头部边的左侧为止。第二段看deque尾部方向判断新顶点相对于d[size-2] → d[size-1]的位置。同样如果新顶点在右侧或共线就从尾部弹出。理解这个规则最直观的类比是橡皮筋已有顶点被一根绷紧的橡皮筋包住新顶点从某个方向挤进来把皮筋顶出去原本支撑皮筋的旧顶点失去作用就被释放掉。deque两端的弹出就是在把失去支撑的顶点从边界上拿掉。弹出不是随意进行的。因为deque里维护的始终是一个凸链新顶点只可能影响与两端相连的那两条边。如果新顶点在最外层它会让某一段老边界失效我们把失效边界上的顶点从两头弹出直到新顶点与两端连成的边重新构成凸多边形此时插入就完成了。这个循环保证了每个顶点最多被弹出一次所以总复杂度是O(n)。3. 完整实现C版Melkman3.1 先搞定叉积和方向判断几何算法里方向判断是最基础的工具。对于三个点a、b、c叉积可以判断c在ab连线左侧、右侧还是线上long long 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); }如果cross大于0说明c在ab的左侧小于0说明c在右侧等于0说明三点共线。这里的“左侧”基于右手坐标系x轴向右、y轴向上。初始化时我们希望构成一个逆时针方向的三角形。判断前三个点的方向时如果cross小于0就交换第二个和第三个点保证最终三角形方向统一为逆时针。这个统一方向在后面输出结果时很有用。3.2 初始化前几个点共线怎么办很多教材里的Melkman实现默认前三个点不共线。但真实工程数据里多边形从某条直边中间开始的情况太常见了不处理前段共线点直接跑结果会丢掉凸包的极端端点。我的做法是找到第一个与起点不重合的点作为p1然后继续向后扫描凡是与p0、p1共线的点都跳过但在跳过过程中始终保留离p0最远的那个点作为p1。这样即使前几个点都在同一条边上也不会漏掉这条边的端点。如果整个多边形所有点共线直接返回p0和离p0最远的点作为结果——这就是一条线段的凸包。完整初始化代码如下int i 1; while (i n poly[i] poly[0]) i; if (i n) return {poly[0]}; Point p0 poly[0], p1 poly[i]; int j i 1; while (j n cross(p0, p1, poly[j]) 0) { if (dist2(p0, poly[j]) dist2(p0, p1)) p1 poly[j]; j; } if (j n) return {p0, p1}; Point p2 poly[j]; if (cross(p0, p1, p2) 0) swap(p1, p2);这里dist2计算的是两点距离的平方用整数运算避免精度问题。初始化deque时注意构建方式dequePoint d; d.push_back(p0); d.push_back(p1); d.push_back(p2); d.push_front(p2);得到的结果就是d [p2, p0, p1, p2]front和back重合中间按逆时针方向排列。3.3 主循环五个点的代码解决全部问题主循环非常短核心逻辑只有几行for (int k j 1; k n; k) { while (cross(d[0], d[1], poly[k]) 0) { d.pop_front(); } while (cross(d[d.size() - 2], d[d.size() - 1], poly[k]) 0) { d.pop_back(); } d.push_front(poly[k]); d.push_back(poly[k]); }第一个while处理头部弹出第二个while处理尾部弹出。两个循环都弹出之后再把新顶点同时插入到两端。注意这里用的是 0表示新顶点在边的右侧或共线时都会弹出旧顶点这样最终凸包只保留真正的角点共线的中间顶点会被去掉。如果你希望保留凸包边上的所有共线点即包含边上的中间点可以把 0改成 0。两种需求都有人遇到我建议默认用 0因为绝大多数碰撞检测和轮廓分析只需要角点。最后输出时由于deque的front和back是同一个点直接删除末尾元素vectorPoint hull(d.begin(), d.end()); hull.pop_back(); return hull;3.4 完整可运行代码示例把上面几段合在一起就是一个可以直接跑通的完整函数#include iostream #include vector #include deque #include algorithm using namespace std; struct Point { long long x, y; bool operator(const Point other) const { return x other.x y other.y; } }; long long 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); } long long dist2(const Point a, const Point b) { long long dx a.x - b.x; long long dy a.y - b.y; return dx * dx dy * dy; } vectorPoint melkman(const vectorPoint poly) { int n poly.size(); if (n 3) return poly; int i 1; while (i n poly[i] poly[0]) i; if (i n) return {poly[0]}; Point p0 poly[0], p1 poly[i]; int j i 1; while (j n cross(p0, p1, poly[j]) 0) { if (dist2(p0, poly[j]) dist2(p0, p1)) p1 poly[j]; j; } if (j n) return {p0, p1}; Point p2 poly[j]; if (cross(p0, p1, p2) 0) swap(p1, p2); dequePoint d; d.push_back(p0); d.push_back(p1); d.push_back(p2); d.push_front(p2); for (int k j 1; k n; k) { while (cross(d[0], d[1], poly[k]) 0) { d.pop_front(); } while (cross(d[d.size() - 2], d[d.size() - 1], poly[k]) 0) { d.pop_back(); } d.push_front(poly[k]); d.push_back(poly[k]); } vectorPoint hull(d.begin(), d.end()); hull.pop_back(); return hull; } int main() { vectorPoint poly {{0,0}, {4,0}, {2,2}, {4,4}, {0,4}}; vectorPoint hull melkman(poly); for (auto p : hull) { cout p.x p.y \n; } return 0; }用这个凹多边形测试顶点依次是(0,0)、(4,0)、(2,2)、(4,4)、(0,4)输出结果是(0,4)、(0,0)、(4,0)、(4,4)即整个多边形的最小凸包内部凹点(2,2)被正确剔除。4. 避坑经验这些细节不注意一定会翻车4.1 共线点处理的两难共线点是Melkman实现里最容易出问题的地方。第一段共线如果直接跳过不处理可能会丢凸包端点但如果过度处理又可能引入非必要的顶点。我的建议是只在初始化阶段做共线预处理主循环中不用额外关心共线点因为 0的判断会自然把凸包边中间的点去掉。需要特别注意的是 0会把所有共线点全部弹出也就是说如果多边形有一条很长的直边中间那些顶点不会出现在结果里。这也是正常的除非你的交付要求里明确写着“保留边上顶点”。4.2 循环弹出的顺序和“保底”判断有的教程代码会先处理尾部再处理头部理论上都能得到正确结果但我在调试时发现先处理头部再处理尾部的顺序更直观因为人可以顺着deque的链方向理解头部是链的起始尾部是链的终结。还有一个细节主循环中弹出可能让deque暂时失去“首尾重合”的状态。比如d从[p2, p0, p1, p2]开始头部弹出p2后变成[p0, p1, p2]front和back不重合。这没有问题因为紧接着的尾部判断和最后的双端插入会把闭合状态恢复。很多初学者在这里犹豫以为deque状态已经坏了其实这是算法正常工作的中间过程。4.3 整数溢出与坐标范围叉积计算里两个坐标差相乘如果坐标范围是1e9叉积结果可能到1e18。用int必炸用long long也要注意上限。坐标范围再大比如1e18级别就需要上高精度或改成浮点数但那已经是另一个话题了。实际项目中如果坐标来自经纬度或者屏幕像素long long完全够用。我习惯在结构体定义里直接写long long杜绝后续改类型带来的连锁修改。4.4 输入不是简单多边形时结果会怎样Melkman算法对输入多边形的自交非常敏感。如果顶点序列对应的多边形有自交边算法的凸性维护逻辑会被破坏输出结果可能是错误的而且不会报错。你只会得到一份看似正常的凸包但某些顶点明显说不过去。所以接入这个算法前务必确认你的输入来源。图像轮廓一般是简单多边形但某些GIS数据或者手工标注的围栏数据可能会出现自交排查方法可以单独写一篇文章我目前的做法是在调试期加一个自交检测开关跑出异常结果时提前暴露问题。4.5 输出方向和起始点由于初始化时强制改成逆时针最终凸包也是逆时针方向。但起始点不是固定的它取决于最后一个被处理的顶点。如果你要求凸包从某个指定顶点开始输出可以对结果做一次旋转。这种情况在GIS里比较常见比如要求首点等于原始轮廓的第一个点。5. 常见问题速查表现象原因解决办法结果里出现凹点初始化方向不对检查前三点是否统一为逆时针凸包端点丢失前段共线点被直接跳过按3.2节的逻辑保留最远点所有点共线时程序崩溃没有处理退化情况返回线段两端点凸包边上多个实心点 0改成 0确认需求再选择保留行为输出首点不一致起始顶点由最后处理点决定按需旋转结果数组自交多边形结果异常输入不是简单多边形预处理检测自交或换算法运行到一半deque只剩两个点一般是退化输入检查是否存在重复点或全共线区域坐标范围大导致WA叉积溢出改用long long或高精度这些坑我几乎都踩过一遍。第一次用的时候输出结果里多了一个内部点排查了很久最后发现是初始三角形方向没有统一头部边和尾部边的判断标准不一致导致的。后来我形成一个习惯所有几何模板里叉积函数和方向判断永远单独封装并且在初始化阶段显式统一逆时针方向这样无论输入是顺时针还是逆时针都不会出错。结语与个人体会这个小众算法最打动我的一点是它把“输入有序”这个信息利用到了极致。很多通用算法为了适配最差情况牺牲了输入中本来就存在的顺序结构而Melkman恰恰是那种“认清问题边界再把优势榨干”的典型。我个人的建议是如果你的轮廓来源稳定、确定是简单多边形直接优先考虑Melkman代码量不比其他算法多运行效率却实打实高一个档次。如果只是偶尔处理散点集那就继续用Andrew单调链。最怕的情况是输入类型不稳定混合场景下强行用Melkman反而要处理各种退化分支。最后补一个小技巧在主循环里临时打印deque的当前状态每处理一个点输出一次链表内容你会非常直观地看到“橡皮筋崩紧、旧点弹出、新点插入”的整个过程。我第一次调通之后把这个可视化循环跑了一遍很多疑问瞬间就解开了。