首页 / 资讯中心 / 文章详情

纯Python实现凸包算法:Graham扫描、Andrew单调链与Jarvis步进全解析

纯Python实现凸包算法:Graham扫描、Andrew单调链与Jarvis步进全解析 ★ FEATURED ARTICLE
简介这是一份用 Python 编写的凸包算法项目面向算法学习者、数据可视化爱好者以及需要在图形处理或地理信息系统中使用几何计算的开发者。资源包共 5 个文件包含 4 个 Python 脚本和 1 个 Markdown 说明文档压缩包仅 7KB轻量、易读可快速导入工程或直接阅读源码。脚本大致覆盖凸包核心实现、结果可视化、单元测试与基准对比几个模块说明文档用来梳理项目结构和调用方式。凸包在碰撞检测、数据聚类、区域边界提取与路径规划中都有典型应用这套代码提供了从几何计算到结果出图的完整链路。已有 261 人浏览学习。读者可以通过测试用例理解 Graham 扫描、Jarvis 步进等常见算法的实现差异并结合基准脚本观察不同数据规模下的耗时表现以便在实际项目中按需选用或扩展自己的几何工具。1. 重写凸包算法不满足于调 scipy.spatial 时的工程化选择凸包算法在计算几何里的地位有点像排序算法在数据结构里出镜率高、代码量不大但自己从头写一遍就会发现到处是边界条件。图像识别里对轮廓点集做简化、地图应用里提取行政区域的边界、游戏物理里把复杂多边形近似成碰撞体背后都是同一个数学问题——给一堆离散点圈出一个顺序正确的凸多边形。用 Python 重写凸包真正的门槛不在“背出算法流程”而在极角排序的符号处理、共线点的取舍、浮点精度带来的玄学报错。这个项目把三种常见实现Graham 扫描、Andrew 单调链、Jarvis 步进用纯 Python 各写了一遍不依赖 scipy.spatial 和 OpenCV自带随机点生成、可视化验证和性能对比脚本。想读懂计算几何底层逻辑的人、需要在离线环境里裁剪一段算法代码的人都能直接拿来跑。2. 项目架构与核心模块从几何原语到 Graham 扫描的完整链路2.1 工程结构与模块边界这套代码没有把三个算法揉在一个文件里而是按“几何原语、算法实现、可视化、脚本入口”拆成了四层。我拿到源码第一件事就是看目录结构清晰的项目往往后续改起来也顺手。convex_hull/ ├── hull/ │ ├── __init__.py # 包入口向上层暴露三个算法函数 │ ├── core.py # 叉积、距离等基础几何运算 │ ├── graham.py # Graham 扫描实现 │ ├── andrew.py # Andrew 单调链实现 │ └── jarvis.py # Jarvis 步进实现 ├── viz/ │ └── plot.py # matplotlib 可视化 ├── scripts/ │ ├── generate_points.py # 随机点集生成 │ ├── run_demo.py # 一键跑通完整流程 │ └── benchmark.py # 三种算法性能对比 ├── tests/ │ └── test_hull.py # 算法正确性回归测试 └── requirements.txtcore.py 独立放几何原语是我比较认可的取舍。叉积、欧氏距离这类函数会被三个算法反复用到放进公共模块能避免每个文件复制一遍而算法之间互不依赖之后想新增 QuickHull 或者 Melkman 算法只要在 hull 目录下加一个文件并在init.py 里登记即可。hull/core.py 的代码量不大但它是整个包的基石。2.2 Graham 扫描基准点选取与极角排序的实现Graham 扫描是所有凸包实现里最经典的一条路先找一个肯定在凸包上的点然后按极角排序最后用栈维护边界。第一步选基准点我一般会选 y 坐标最小的点y 相同时取 x 最小的直接用 min 函数搞定。# hull/core.py from math import hypot from typing import List, Tuple Point Tuple[float, float] def cross(o: Point, a: Point, b: Point) - float: 向量 oa 与 ob 的叉积正数表示 b 在 a 的逆时针方向 return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0]) def dist(a: Point, b: Point) - float: 两点欧氏距离用于极角排序时处理同角度点 return hypot(a[0] - b[0], a[1] - b[1])叉积是整套计算的“判向器”cross(o, a, b) 大于 0 说明从 oa 转向 ob 是逆时针小于 0 是顺时针等于 0 则三点共线。这里用浮点运算结果可能带微小误差所以后面避坑章节会提到加 epsilon 判断。dist 函数在极角相同时用来按距离排序保证排序结果稳定。排序这一步是 Graham 和 Andrew 最大的区别。Graham 用的是极角排序需要写一个比较器传给 functools.cmp_to_key比较器内部只做叉积和距离判断不调用三角函数原因很简单atan2 转角度会有额外浮点开销而且角度相等时还得再按距离排一次直接用叉积符号就能拿到严格的偏序关系。# hull/graham.py from functools import cmp_to_key from .core import cross, dist, Point def graham_scan(points: List[Point]) - List[Point]: pts list(set(points)) if len(pts) 3: return sorted(pts) pivot min(pts, keylambda p: (p[1], p[0])) def compare(a: Point, b: Point) - int: c cross(pivot, a, b) if c 0: return -1 if c 0: return 1 da, db dist(pivot, a), dist(pivot, b) return -1 if da db else (1 if da db else 0) pts.sort(keycmp_to_key(compare)) stack [pts[0], pts[1]] for p in pts[2:]: while len(stack) 2 and cross(stack[-2], stack[-1], p) 0: stack.pop() stack.append(p) return stack参数说明points 是任意长度的 Point 列表返回值是逆时针排列的凸包顶点列表。函数入口先用 set 去重既避免重复点影响比较器又把输入规模降下来len(pts) 3 时直接 sorted 返回因为 0 到 2 个点不需要算凸包。pivot 选最下点之后所有点按相对 pivot 的极角升序排列角度相同的按距离升序排完序后 pts[0] 就是 pivot 自己。排序完成后进入栈扫描栈里先放前两个点然后从第三个点开始迭代。每来一个新点 p检查栈顶两个点和 p 是否构成逆时针转向如果 cross 0说明栈顶点要么在 p 的顺时针方向要么与 p 共线这两种情况都不该留在凸包边界上于是弹栈。这个过程是 Graham 扫描最核心的收敛逻辑循环直到栈顶方向满足左转条件再把 p 压栈。2.3 栈维护与共线点取舍为什么弹栈条件必须带等号弹栈条件写成 cross 0 而不是 0这是有意的设计选择。带等号意味着凡是三点共线的情况全部忽略中间点只保留线段两端的拐点输出的是“严格凸包”——顶点序列里不会出现刚好落在边上的点。如果你做轮廓简化时希望保留共线边界点把 改成 就能把中间点留在结果里但那样顶点数量会增加后续如果用顶点数估算复杂度结果会偏大。从复杂度看Graham 扫描的瓶颈不在栈维护而在排序整体是 O(n log n)。栈扫描本身是线性时间每个点最多入栈一次、出栈一次。这个特性让它成为三种算法里综合表现最稳的也是我把 graham_scan 放在包导出函数第一位的原因——多数场景下默认选择它不会错。3. 算法谱系与可视化Andrew 单调链、Jarvis 步进和结果验证3.1 Andrew 单调链一次排序、两条链Andrew 单调链是 Graham 的“排序简化版”。它不需要选基准点也不需要极角比较器直接把点按 (x, y) 字典序排好然后从左到右构造下链从右到左构造上链最后拼起来。# hull/andrew.py from .core import cross, Point from typing import List def andrew_hull(points: List[Point]) - List[Point]: pts sorted(set(points)) if len(pts) 3: return pts def build(p_seq: List[Point]) - List[Point]: hull [] for p in p_seq: while len(hull) 2 and cross(hull[-2], hull[-1], p) 0: hull.pop() hull.append(p) return hull lower build(pts) upper build(pts[::-1]) return lower[:-1] upper[:-1]参数说明输入输出约定和 graham_scan 一致。set 去重之后排序保证了遍历顺序是从左往右、从下往上build 内部逻辑和 Graham 的栈扫描一模一样只是不需要极角排序每次迭代天然按 x 递增处理点。lower 返回下链upper 构造时把 pts 反转相当于从右往左再走一遍。拼接时各去掉最后一个点因为 lower 的终点和 upper 的终点是同一个点最右点upper 的终点和最左点也是同一个点最左点去掉重复才能得到完整的逆时针环。不少人在初学时会疑问既然两个算法复杂度都是 O(n log n)为什么还要区分区别隐藏在常数和实现风险上。Andrew 的排序用内置 sort天然稳定可靠Graham 的手写极角比较器一旦遇到重复点或浮点误差容易出现非传递性比较。Andrew 在绝大多数场景下代码更短、更不容易写错所以我个人在工程里其实更偏好 Andrew。3.2 Jarvis 步进礼品包装的直观与代价Jarvis 步进又叫做礼品包装法思路非常直观从最左下角的点出发每次在所有点里找一个“在当前点看来逆时针转角最小”的点作为下一个边界点一路包下去直到回到起点。# hull/jarvis.py from .core import cross, dist, Point from typing import List def jarvis_march(points: List[Point]) - List[Point]: pts list(set(points)) if len(pts) 3: return pts start min(pts, keylambda p: (p[1], p[0])) hull [start] current start while True: nxt pts[0] for p in pts: if p current: continue c cross(current, nxt, p) if c 0 or (c 0 and dist(current, p) dist(current, nxt)): nxt p if nxt start: break hull.append(nxt) current nxt return hull参数说明输出同样是逆时针顶点列表。每次迭代的 nxt 初始化为 pts[0]然后遍历全集找“更逆时针”的点cross 0 表示 p 比当前 nxt 更靠左c 0 且距离更远时用最远点这一步确保共线的情况下跳过中间点直接取端点。当前点 p 必须跳过否则会出现自己选自己的死循环。当 nxt 回到 start 说明已经包完一整圈。Jarvis 的复杂度是 O(nh)h 是凸包顶点数。均匀随机分布下 h 大约只有 O(n^1/3)所以小数据量表现不算差但如果输入是圆形分布且所有点都在边界上h 接近 n复杂度退化成 O(n²)。这也是我在工程里把 Jarvis 定位成“教学验证型”而不是生产首选的原因。3.3 可视化模块把算法输出叠到散点图上看算法写对了没有肉眼比断言更直观。项目里的 viz/plot.py 用 matplotlib 把原始点集和凸包边界画在同一张图上带填充区域一眼能看出顶点顺序是否正确。# viz/plot.py import matplotlib.pyplot as plt def plot_hull(points: list, hull: list, title: str Convex Hull): xs [p[0] for p in points] ys [p[1] for p in points] plt.figure(figsize(8, 6)) plt.scatter(xs, ys, s14, csteelblue, label原始点集) if len(hull) 2: hx [p[0] for p in hull] [hull[0][0]] hy [p[1] for p in hull] [hull[0][1]] plt.plot(hx, hy, ccrimson, lw2, label凸包边界) plt.fill(hx, hy, colororange, alpha0.12) plt.title(title) plt.legend(locbest) plt.gca().set_aspect(equal) plt.grid(alpha0.3) plt.show()逻辑说明scatter 画原始散点plot 把 hull 顶点依次连线注意末尾要追加 hull[0] 才能闭合成环fill 填充内部区域alpha 设小一点避免遮住原始点。参数说明points 是散点坐标列表hull 是算法返回的凸包顶点列表title 用于区分三个算法的输出图。set_aspect(equal) 保证横纵轴比例一致否则圆会被拉伸成椭圆凸包形状看起来是变形的。4. 从环境准备到基准测试让三种算法在本地跑起来4.1 依赖与启动三分钟跑通基础用例这个项目的外部依赖只有一个 matplotlib没有 numpy、没有 scipy这是刻意的——凸包核心逻辑只需要标准库就能完成依赖少意味着在离线环境或精简镜像里更容易部署。装依赖和跑 demo 的命令如下。pip install -r requirements.txt python scripts/run_demo.pyrun_demo.py 是项目的入口脚本它会生成 200 个随机点依次用三种算法计算凸包并弹出三张可视化图。脚本开头有一行路径处理值得留意# scripts/run_demo.py import sys from pathlib import Path sys.path.insert(0, str(Path(__file__).resolve().parents[1])) from hull import graham_scan, andrew_hull, jarvis_march from scripts.generate_points import random_points from viz.plot import plot_hull def main(): points random_points(200, seed7) for algo in (graham_scan, andrew_hull, jarvis_march): hull algo(points) print(f{algo.__name__}: {len(hull)} 个凸包顶点) plot_hull(points, hull, titlefConvex Hull via {algo.__name__}) if __name__ __main__: main()如果你直接用 python scripts/run_demo.py 运行Python 会把 scripts 目录当作 sys.path[0]此时 import hull 会报 ModuleNotFoundError。加一行 sys.path.insert 把项目根目录塞进搜索路径是常见做法。参数说明seed7 固定随机种子你可以改成任意整数复现同一点集三个算法函数签名完全一致循环里可以直接替换这也给后续做算法混用提供了便利。4.2 随机点集生成均匀分布与圆形退化场景随机点集的生成方式决定了测试覆盖范围。纯均匀分布能验证常规输入但凸包算法最容易翻车的是退化输入——全部共线、大量重复点、点分布在圆上。generate_points.py 里我建议至少支持两种分布模式。# scripts/generate_points.py import random import math def random_points(n: int 100, seed: int None, mode: str uniform): if seed is not None: random.seed(seed) if mode uniform: return [(random.uniform(0, 100), random.uniform(0, 100)) for _ in range(n)] if mode circle: return [ (100 50 * math.cos(2 * math.pi * i / n), 100 50 * math.sin(2 * math.pi * i / n)) for i in range(n) ] if mode line: return [(i * 1.0, 0.0) for i in range(n)] raise ValueError(funsupported mode: {mode})参数说明n 控制点数量seed 用于复现随机序列mode 支持 uniform、circle、line 三种。circle 模式下所有点都在圆周上凸包顶点数等于 n是测试 Jarvis 退化的理想输入line 模式让所有点落在一条直线上用来验证三个算法对共线输入是否都返回两端点。4.3 基准测试O(n log n) 与 O(nh) 的曲线长什么样benchmark.py 对三种算法做了简单的时间对比通过增大点集规模观察复杂度趋势。# scripts/benchmark.py import sys import time from pathlib import Path sys.path.insert(0, str(Path(__file__).resolve().parents[1])) from hull import graham_scan, andrew_hull, jarvis_march from scripts.generate_points import random_points def run_benchmark(): for n in (100, 1000, 10000): pts random_points(n, seed42, modeuniform) for algo in (graham_scan, andrew_hull, jarvis_march): t0 time.perf_counter() algo(pts) elapsed time.perf_counter() - t0 print(fn{n:6} {algo.__name__:14} {elapsed * 1000:8.2f} ms) if __name__ __main__: run_benchmark()在普通开发机上n 10000 时 Graham 和 Andrew 大概都在几十毫秒量级Jarvis 因为均匀分布的凸包顶点数远小于 n也不会太离谱但常数明显偏大。真正拉开差距的是 circle 模式把 benchmark 里的 mode 改成 circle 再跑n 10000 时 Jarvis 会慢到以秒计而 Graham 和 Andrew 依然稳定在几十毫秒。这个对比能直观验证 O(nh) 的退化风险。tests/test_hull.py 里还可以做一个一致性回归三个算法在同一随机点集上的输出经过归一化统一从最左下角的顶点开始按逆时针列举后应该完全相同。这个断言比单测某个固定用例更能捕捉实现里的隐藏 bug我建议保留下来后续每次改动算法都跑一遍。5. 避坑记录浮点精度、重复点与可视化环境的排障笔记5.1 浮点误差把共线点误判成左右转现象明明三个点理论上是共线的cross 算出来却是 1e-15 这种非零小值导致 Graham 扫描没有按预期弹栈凸包边界上多了一个“几乎共线”的顶点面积计算时误差被放大。 原因输入坐标经过旋转、缩放或从浮点文件读入后本来就带舍入误差叉积结果在零点附近抖动。 解决在 core.py 里加一个带 epsilon 的叉积包装函数绝对值小于阈值直接按 0 处理。EPS 1e-12 def cross_eps(o: Point, a: Point, b: Point) - float: c cross(o, a, b) return 0.0 if abs(c) EPS else c把 graham.py 和 andrew.py 里的 cross 调用换成 cross_eps弹栈条件改成 0 就能稳定处理共线退化。阈值 1e-12 不是绝对真理如果你的坐标量级是百万级别叉积量级会变大阈值要按坐标量级相应放大反之坐标在 0 到 1 之间1e-14 更合适。5.2 重复点让极角比较器失去传递性现象点集里有完全相同的坐标Graham 的 compare 函数返回结果不稳定sorted 结果里出现顺序错乱凸包顶点里混进重复点。 原因两个点与 pivot 的叉积和距离都相等时比较器返回 0 才是正确的等价类处理如果距离也相等却返回 1比较关系就会出现 a b 且 b a 的矛盾破坏传递性。 解决入口处先 set 去重并且 compare 里距离相等时必须返回 0。项目代码里已经做了这两件事如果你在自己的实现里抄漏了 set 去重小样本可能没事点集一上千就复现排序异常。5.3 全共线输入画出形态怪异的“多边形”现象所有点都在一条直线上时算法正确返回两个端点但可视化代码如果无条件调用 fill会画出面积为零却带两条边的奇怪图形或者报 ValueError。 原因fill 需要至少三个不共线的点才能构成合法多边形两个点无法填充。 解决plot_hull 里先判断 len(hull) 3 才执行 filllen(hull) 2 才执行 plot。这个防御性判断不算多余——真实业务里的输入点集永远可能退化。5.4 Jarvis 在退化输入上的死循环隐患现象jarvis_march 在某些共线点集上运行时间异常长甚至不退出。 原因如果遍历时不跳过 current 点nxt 有可能被设回当前点导致每轮都在局部打转或者共线时选了距离最近而不是最远的点导致一圈要经过所有中间点h 膨胀成 n复杂度退化成灾难级。 解决遍历前跳过 p current共线时取距离最远的点直接从一端跳到另一端。这样即使全部共线也能快速从起点跳到终点再跳回来保证终止。5.5 无显示环境下 matplotlib 直接抛 Invalid DISPLAY现象在只有终端的 Linux 服务器上跑 run_demo.pyplt.show() 抛 RuntimeError: Invalid DISPLAY。 原因matplotlib 默认后端依赖图形界面服务器没有 X Server。 解决给可视化脚本加上 Agg 后端切换并支持保存图片代替弹窗。import matplotlib matplotlib.use(Agg) import matplotlib.pyplot as plt plt.savefig(hull.png, dpi120)把 show 换成 savefig 后脚本在无头环境下也能批量出图方便做回归测试或写进 CI 流程。6. 二次开发思路把凸包算法改造成可插拔的几何工具6.1 策略模式注册表新增算法只需一行登记目前三个算法在 hull/init.py 里逐个导出调用方用哪个算法是写死的。更灵活的做法是维护一个算法注册表把函数名作为配置项传入业务代码。# hull/__init__.py from .graham import graham_scan from .andrew import andrew_hull from .jarvis import jarvis_march ALGORITHMS { graham: graham_scan, andrew: andrew_hull, jarvis: jarvis_march, }业务侧需要换算法时只改配置字符串hull ALGORITHMS algo_name 。要新增 QuickHull实现一个相同签名的函数在注册表里加一行即可所有调度逻辑不用动。6.2 业务落地GeoJSON 导出与碰撞体近似凸包顶点算出来后最常见的落地形态是导出成 GeoJSON 多边形或者作为物理引擎的碰撞包围体。导出代码很短但值得单独封装。import json def hull_to_geojson(hull): if len(hull) 3: return None coords [list(p) for p in hull] [list(hull[0])] return { type: Feature, properties: {}, geometry: {type: Polygon, coordinates: [coords]}, }在游戏项目里把复杂角色模型的顶点集先算凸包再拿凸包顶点去生成碰撞体物理检测的开销能小一个量级。我当时接手的某个模拟项目角色网格顶点超过一万直接用原始网格做碰撞检测每帧都在掉帧换成这个项目里的 Andrew 算法生成凸包后再做简化碰撞体帧率立刻回升。从那以后我每次往项目里加几何算法都会强制走一遍流程先跑一致性回归再切退化输入测边界最后才接业务逻辑。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助?
咨询建站