ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

凸包算法全解析:从Andrew单调链到点云处理实战

凸包算法全解析:从Andrew单调链到点云处理实战 1. 凸包到底是什么——先用直觉建立一个图像如果你在一张木板上钉下一堆钉子拿一根橡皮筋把所有钉子都圈住然后松手橡皮筋会收缩成一个紧贴最外层钉子的多边形。这个多边形围出来的区域就是这一堆点的凸包。稍微严谨一点说给定平面上的点集 S凸包是包含 S 中所有点的最小凸集。所谓凸集就是集合内任意两点连线上的所有点仍然在这个集合内。你可以在脑子里做个实验——取凸包内任意两个点拉一条直线段这条线段必然完全落在凸包内部不会跑出去。满足这个性质的“包裹住所有点”的形状里面积最小的那个就是凸包。从计算几何的角度来看凸包是一系列有方向、有序的顶点序列。注意“有序”这两个字后面所有算法都在解决同一个问题找出这些顶点并且按顺序把它们排好。这个顺序通常指逆时针方向。这个问题的古典味道很重但直到今天它依然是计算机图形学、路径规划、点云处理、碰撞检测等领域的底层地基。你后面会看到很多看起来复杂的任务第一步都是把散乱的点集外包络求出来。2. 为什么凸包值得你花一下午认真研究很多人学凸包是被算法题逼的刷完就忘了。但实际工程里凸包出现的频率远超你的想象。2.1 碰撞检测与物理引擎游戏引擎里复杂的物体形状会被近似成一个或多个凸多面体。判断两个凹物体是否碰撞很麻烦但判断两个凸体是否碰撞就快得多——分离轴定理SAT就是建立在凸体前提上的。求物体点集的凸包相当于给物理引擎喂一个干净、可靠的输入。2.2 点云处理与三维重建这是我工作中接触最多的一块。激光雷达扫描得到的点云是散乱的、带噪声的。在对点云做配准、分割、特征提取之前经常需要先求出点云的外边界也就是三维点云的凸包。它能帮你快速估算物体占据的空间范围、计算体积或者剔除离群点。后面我会专门讲三维点云的情况。2.3 规划与控制移动机器人路径规划时需要把障碍物膨胀成凸区域来简化碰撞检测无人机的飞行走廊也常用凸多边形来近似可用空间。这些场景里凸包算法是标准前置步骤。2.4 数据分析与可视化做聚类分析时用凸包圈出每个簇的可视化边界评估一组数据的分布范围时凸包比“四个角的最大最小值”更准确。比如 GPS 轨迹的外边界提取本质上就是个带噪声的凸包问题。2.5 一个容易被忽略的价值算法思维的训练凸包算法是你遇到的第一个“几何 算法”深度结合的问题。它要求你把几何直觉左转还是右转转化成程序逻辑叉积符号还要考虑数值误差、退化情况共线、重复点。这些能力是刷多少动态规划题都练不出来的。3. 基于常见实践的补充说明需要坦白一点凸包的算法变体很多没有哪个实现是万能药。不同场景下时间复杂度和实现复杂度要做一个权衡。接下来的几章我会把我实际用过的、验证过的几种主流算法逐一拆开讲提供可以“抄作业”的代码并指出每个算法的致命弱点在哪里。4. Jarvis 步进法——最符合直觉的算法也叫 Gift Wrapping 算法中文通常叫卷包裹法。它的思路跟“给礼物盒子缠一圈胶带”一模一样。4.1 核心思想先找到最左下角的点 P0y 坐标最小y 相同时取 x 最小。这个点必然在凸包上。从 P0 出发找下一个点要求是在所有点中这个点相对于当前点的极角最小。一直重复这个过程直到回到 P0。你想象自己在绕着点集边缘走每到一个点环顾四周选一个“最靠左”的朋友继续走。走完一圈你留下的足迹就是凸包。4.2 时间复杂度最好情况是 O(nh)n 是总点数h 是凸包顶点数。最坏情况下所有点都在凸包上比如所有点均匀分布在一个圆周上h n复杂度退化为 O(n²)。4.3 Python 实现def cross(o, a, b): # 向量 oa 和 ob 的叉积 return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0]) def jarvis_march(points): n len(points) if n 3: return points # 找最左下角的点 start min(points, keylambda p: (p[1], p[0])) hull [] current start while True: hull.append(current) endpoint points[0] for p in points[1:]: # 如果 p 在 current-endpoint 的右侧更新 endpoint if endpoint current or cross(current, endpoint, p) 0: endpoint p current endpoint if current start: break return hull注意上面这段代码有个隐藏 bug——如果endpoint初始值points[0]恰好等于current或者叉积为 0共线处理起来会出问题。工程上更稳妥的写法是用一个on_right函数并显式处理共线情况。4.4 适合什么场景Jarvis 的 O(nh) 在 h 远小于 n 时非常快。比如一个 10 万个点的圆形点集你要求凸包h 可能只有几十这比所有排序类算法都快。反过来如果点均匀撒在圆周上它就是灾难。5. Graham 扫描法——排序是算法的灵魂Graham 扫描法的历史地位很高也是很多教材的标配。它的想法很妙如果我们能按极角把点排好序那么凸包的构造就可以在线性时间内完成。5.1 算法步骤拆解选基准点同样是取最左下角的点 P0。极角排序把其余所有点按相对于 P0 的极角从小到大排序。如果极角相同保留距离最远的那个因为更远的点会“罩住”近的点。扫描用一个栈维护当前凸包候选。依次访问排序后的点每次判断“栈顶第二个点 → 栈顶点 → 当前点”构成的折线是左转还是右转。如果是右转或直行说明栈顶点不是凸包顶点弹出继续检查直到满足左转条件再把当前点入栈。5.2 极角排序的关键细节排序是 Graham 的灵魂也是最容易出浮点误差的地方。计算极角最常用的方式是atan2(dy, dx)但浮点误差和边界问题很烦人。实测建议如果坐标是整数用叉积比较做排序键能规避大量精度坑。但叉积比较不满足严格弱序需要小心处理共线情况。如果是工程项目而非竞赛直接用atan2排序通常就够了。5.3 Python 实现import math def cross(o, a, b): return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0]) def graham_scan(points): points sorted(set(points)) # 去重并排序 if len(points) 3: return points # 以最左下点为基准按极角排序 p0 points[0] rest sorted(points[1:], keylambda p: math.atan2(p[1] - p0[1], p[0] - p0[0])) stack [p0, rest[0]] for p in rest[1:]: while len(stack) 2 and cross(stack[-2], stack[-1], p) 0: stack.pop() stack.append(p) return stack注意第 11 行的 0这个细节决定了你是否保留凸包边上的共线点。用 0会保留用 0会去掉。根据你要不要“最外层边上的中间点”灵活切换。5.4 Graham 的软肋从工程角度Graham 的软肋在于极角排序中的atan2调用。10 万个点就要调用 10 万次atan2性能开销不小。另外浮点误差在排序阶段可能把两个本来应该顺序互换的点排反导致后续栈操作出现微小偏差。6. Andrew 单调链算法——我实际项目里最常用的方案如果你只打算记住一种凸包算法我推荐 Andrew 单调链Monotone Chain也叫水平序算法。6.1 为什么它比 Graham 更实用不需要按极角排序直接按 x 坐标排序。整个排序过程只涉及比较操作误差积累更少。不需要显式算atan2只用叉积判断转向。思路极其对称代码不容易写错。上下链分别构造边界情况好处理。6.2 算法流程把所有点按 x 坐标排序x 相同按 y 排序。构造下凸壳lower hull从左往右扫维护栈保证栈中相邻三点的转向始终是逆时针叉积 0即左转。构造上凸壳upper hull从右往左扫同理。合并下凸壳和上凸壳去掉首尾重复的点得到完整凸包。6.3 Python 完整实现def cross(o, a, b): return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0]) def convex_hull(points): points sorted(set(points)) if len(points) 3: return points lower [] for p in points: while len(lower) 2 and cross(lower[-2], lower[-1], p) 0: lower.pop() lower.append(p) upper [] for p in reversed(points): while len(upper) 2 and cross(upper[-2], upper[-1], p) 0: upper.pop() upper.append(p) # 去掉 lower 的最后一个点和 upper 的最后一个点它们重复 return lower[:-1] upper[:-1]这段代码很紧凑但信息量很大。逐行说明sorted(set(points))同时做了三件事去重、按 x 排序、x 相同时按 y 排序。cross(lower[-2], lower[-1], p) 0意味着如果当前点是右转或直行就把栈顶弹出去。因为我们从左往右扫要求的是下凸壳保持“逆时针左转”。上链从右往左扫逻辑完全对称。6.4 C 参考实现如果你在写高性能代码C 版本如下using Point pairlong long, long long; long long cross(const Point o, const Point a, const Point b) { return (a.first - o.first) * (b.second - o.second) - (a.second - o.second) * (b.first - o.first); } vectorPoint convexHull(vectorPoint pts) { sort(pts.begin(), pts.end()); pts.erase(unique(pts.begin(), pts.end()), pts.end()); if (pts.size() 3) return pts; vectorPoint lower, upper; for (const auto p : pts) { while (lower.size() 2 cross(lower[lower.size()-2], lower.back(), p) 0) lower.pop_back(); lower.push_back(p); } for (int i pts.size() - 1; i 0; --i) { const auto p pts[i]; while (upper.size() 2 cross(upper[upper.size()-2], upper.back(), p) 0) upper.pop_back(); upper.push_back(p); } lower.pop_back(); upper.pop_back(); lower.insert(lower.end(), upper.begin(), upper.end()); return lower; }long long 的用法要特别注意叉积的中间结果可能达到坐标差的平方量级用 int 会溢出。坐标范围在 1e9 时差值平方是 1e18正好卡在 long long 的边界附近如果坐标到 1e18 那还得用 __int128这是很多人在竞赛和工程里踩过的坑。7. 常见算法的横向对比与选型建议算法时间复杂度核心思想优点缺点适合场景Jarvis 步进法O(nh)每次找极角最小点实现最简单、h小时极快h大时退化O(n²)凸包顶点少的场合Graham 扫描法O(n log n)极角排序栈维护思路经典、代码简洁依赖atan2、浮点排序误差竞赛题、对性能不敏感Andrew 单调链O(n log n)x排序上下链构造无三角函数、数值稳健需要理解对称构造工程首选、点云处理QuickHull平均O(n log n)分治划分左右点高维推广容易最坏O(n²)三维凸包的基础分治算法O(n log n)左右递归再合并理论基础扎实实现复杂教学、并行化场景我的选型逻辑99% 的情况下直接上 Andrew 单调链。原因不是别的就是它最不容易出 bug而且浮点误差可控。如果你是在竞赛环境下刷题Graham 模板背熟也行但 Andrew 同样简单为什么不选数值上更稳的那个8. 数值稳定性——凸包算法里最隐蔽的杀手很多初学者把凸包算法跑通一次就认为自己会了。但实际工程里的点集不会像教科书那么规整坐标可能极大、极小可能布满共线点也可能有重复点。这里面的坑远比想象的多。8.1 叉积的浮点误差def cross(o, a, b): return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])这个函数说起来简单但如果坐标是浮点数比如(0.1, 0.2)这样的值差乘结果可能在真正的 0 附近震荡。你拿 0判断左转可能一会儿为1e-16一会儿为-1e-16凸包直接穿洞或者多出莫名其妙的顶点。做法是引入一个小量 epsEPS 1e-10 def is_left_turn(o, a, b): return cross(o, a, b) EPS def is_right_turn(o, a, b): return cross(o, a, b) -EPS def is_collinear(o, a, b): return abs(cross(o, a, b)) EPSEPS 的取值很讲究。太大会把真正的左转误判成共线丢掉凸包关键顶点太小又起不到滤波作用。我常用的经验公式EPS 1e-9 * max(1, max_coord)^2其中 max_coord 是坐标绝对值的最大值。这样 EPS 量级随数据范围缩放适应性更好。8.2 用整数坐标规避浮点误差如果你的输入坐标本身是整数很多点云体素化后就是整数坐标完全可以别用浮点直接用 long long 做叉积彻底告别精度问题。这是竞赛选手的常规操作也是工程里常见的取巧手段。多数情况下浮点数坐标也能通过缩放取整来近似处理但要评估误差可接受范围。8.3 重复点与共线点的处理输入里出现重复点最直接的办法是排序后unique去重——Andrew 实现里我已经内置了这个处理。共线点则看你的需求如果你只需要凸包的“骨架顶点”用 0把所有共线中间点弹掉。如果你需要凸包边上所有点比如做栅格地图边界提取改用 0让共线点保留在凸包边上。这个选择必须在算法入口处作为参数暴露出来不要写死在函数里否则后面调用方会很痛苦。9. 从二维到三维点云凸包到底是怎么回事接下来这一节面向的是做点云处理、三维视觉、机器人感知的读者。很多时候二维凸包不够用你需要求三维点云的凸包——比如给一堆积木重建出它的外表面。9.1 三维凸包的难度骤增二维凸包可以在 O(n log n) 内解决三维凸包的增量法复杂度最优也就是 O(n²)而且实现难度比二维高一个量级。三维修空间里“转向”的概念变成了“朝向”叉积变成了面法向量的点积判断。你不再维护一个栈而是维护一个由三角形面片组成的网格结构每插入一个新点要删除所有“能看到”这个点的面。9.2 工程上的高阶替代方案PCL 库实际工作中不要自己造三维凸包的轮子。PCLPoint Cloud Library中pcl::ConvexHull已经帮你封装好了。pcl::PointCloudpcl::PointXYZ::Ptr cloud(new pcl::PointCloudpcl::PointXYZ); pcl::ConvexHullpcl::PointXYZ hull; hull.setInputCloud(cloud); hull.setDimension(3); std::vectorpcl::Vertices polygons; pcl::PointCloudpcl::PointXYZ::Ptr surface_hull(new pcl::PointCloudpcl::PointXYZ); hull.Reconstruct2D(surface_hull, polygons); // 注意这个函数名是历史的坑这段代码返回的是凸包表面三角形网格可以配合pcl::visualization画出来也可以直接用于体积估计或碰撞检测。9.3 三维点云凸包的常见应用体素滤波后的外边界规整体素化后的点云边缘是锯齿状的凸包可以给出平滑的外包面。地面滤波在自动驾驶中用最低点集的凸包拟合地面边界。抓取规划机械臂抓取前用凸包近似物体的包围体积用于避碰计算。3D IoU 计算两个三维凸包的重叠体积比是物体检测中常用的评价指标。10. 工程中真正会用到的几个心得这一节是我的实战备忘写出来主要是方便你在自己的项目里少走弯路。10.1 输入数据先做归一化如果坐标范围差异极大比如 x 在 0~1e6y 在 0~1e-4叉积计算很容易出现严重的量级失衡。先用平移缩放把点云归一化到原点附近统一尺度计算完再映射回去。这样数值误差可控不出奇怪问题。10.2 测试集至少准备这五类随机圆上的点凸包顶点数 n随机正方形内的点h 远小于 n大量共线点 少量随机点考验共线处理完全共线的点集退化到只剩一条线只有一个点和两个点的极端数据每次改算法实现先把这五类全跑一遍再谈优化。10.3 排序稳定性对凸包有影响吗值得说明的是在 Andrew 单调链里排序后相同 x 坐标的点会按 y 排序如果坐标是浮点数且 x 坐标完全相同的点不少要看清楚栈的弹出逻辑是否会把这些点处理干净。实测下来只要去重做得好排序稳定性对最终凸包影响不大。但如果共线点比较多建议用稳定排序Python 的sorted保证稳定C 的sort不是。10.4 一个容易被忽略的性能细节如果你需要循环求很多个点集的凸包比如每帧点云都求一次提前把点集按 x 排序并压缩到连续内存能让排序阶段快很多。Python 里尽量用 numpy 数组存储点集排序和叉积都可以向量化比纯 Python 循环快几个量级。import numpy as np def convex_hull_np(points): pts np.unique(points, axis0) if len(pts) 3: return pts def cross2d(o, a, b): return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])但这种向量化写法牺牲了可读性我的经验是先写清楚纯 Python 版本确认逻辑正确再按需优化成 numpy 版本不要一开始就陷入性能优化。11. 一个完整的实操案例从散乱点云到凸包可视化为了让你对整体流程有个完整的感知我提供一个可直接复现的小案例生成一个随机点云计算凸包然后画出来。import matplotlib.pyplot as plt import numpy as np import random def convex_hull(points): points sorted(set(points)) if len(points) 3: return points lower [] for p in points: while len(lower) 2 and cross(lower[-2], lower[-1], p) 0: lower.pop() lower.append(p) upper [] for p in reversed(points): while len(upper) 2 and cross(upper[-2], upper[-1], p) 0: upper.pop() upper.append(p) return lower[:-1] upper[:-1] def cross(o, a, b): return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0]) # 生成测试数据一个圆形轮廓 内部随机点 random.seed(42) points [] for _ in range(200): angle random.uniform(0, 2 * np.pi) radius random.uniform(0.8, 1.2) points.append((radius * np.cos(angle), radius * np.sin(angle))) for _ in range(100): points.append((random.uniform(-1, 1), random.uniform(-1, 1))) hull convex_hull(points) hull.append(hull[0]) # 闭合多边形 # 可视化 xs, ys zip(*points) hx, hy zip(*hull) plt.figure(figsize(6, 6)) plt.scatter(xs, ys, s5, cblue, labelpoints) plt.plot(hx, hy, cred, linewidth2, labelconvex hull) plt.legend() plt.axis(equal) plt.show()这个案例里我把 200 个点撒在圆环上再混入 100 个内部随机点。输出结果你应该会看到红色凸包完美贴合外部圆环上的点内部随机点全部被包裹。这个测试场景能有效检验算法对“凸包顶点完全由外围点构成”的处理能力。12. 关于快速排序的一个联想——顺带聊一句看到热搜词里有“冒泡排序算法c”和“排序算法”这些词我多说一句凸包算法的复杂度瓶颈在排序所以对排序稳定性和性能的追求最终会传导到凸包性能上。你在学凸包时如果顺便把归并排序、快速排序的底层差异搞明白对写凸包工程代码有不小帮助。这也是为什么很多算法教材把凸包安排在排序章节之后讲——排序真的是它的“前置技能”。13. 下一步你可以做什么陆陆续续讲了这么多核心的代码也都在上面了。你在自己的项目里去用就好了。真要我说什么总结出来的经验我觉得最值钱的是三条第一默认用 Andrew 单调链。别在基础选型上花过多时间反复横跳。Tom 我见过太多人纠结半天最后选了 Jarvis结果数据一多就崩又回来换算法。第二把 EPS 处理当成算法的一部分写进去。哪怕是整数坐标的输入你也要写清楚边界条件。这不是过度设计而是让自己以后被人问起“这里为什么这么判断”时能答得清楚。第三多准备几组特殊测试数据。凸包这种几何算法最常见的 bug 不是逻辑错而是数学退化情况没处理干净。全共线、只有一个点、两个点重合这些情况你的函数都必须能优雅地返回结果而不是抛异常或者死循环。把这些边界测试写进自动化测试里以后重构算法实现时你会庆幸自己这么做过。另外提一句如果你学了凸包之后觉得计算几何有意思值得往深挖的两个方向是 Voronoi 图 和 Delaunay 三角剖分。它们跟凸包有着千丝万缕的联系——Delaunay 三角剖分的外边界其实就是凸包。理解凸包之后这两个概念的上手成本会大大降低。
返回列表