
简介scikit-geometry 是一套面向科学计算与几何算法开发的 Python 封装库适合具备一定 Python 基础、希望借助成熟计算几何能力解决实际问题的开发者与研究者。它主要封装 CGAL 的核心功能提供 Point2/Point3、Segment、Ray、Line、Vector、Circle、Plane、Triangle、Polyhedron 等基本几何类型并支持相交测试、凸包、Voronoi 图、布尔运算、骨架提取、主成分分析等操作可用于机器人路径规划、CAD 建模、地理信息处理等场景。资源包共 61 个文件约 273KB包含 18 个 cpp 源文件、12 个 py 脚本、10 个 ipynb 示例笔记以及 yml 配置、hpp 头文件、rst 文档等覆盖从底层实现到交互演示的完整结构。目前已有 2123 人学习下载读者可借助 Jupyter Notebook 示例快速上手结合测试脚本与文档理解各模块调用方式并参考 C 源码深入掌握几何算法的封装思路与工程组织。1. 当你的 Python 脚本开始算几何scikit-geometry 到底补了哪块短板用 Python 做几何计算的人大概都经历过这个阶段一开始用 shapely 处理点线面做做缓冲区、求求交集挺顺手等到需要做布尔运算的精确判定、多边形偏移、凸包、三角剖分甚至要处理带圆弧的曲线时shapely 就开始力不从心要么精度对不上要么干脆没有对应接口。这时候很多人会转向 C 的 CGAL但 CGAL 的学习曲线和编译成本又劝退一大片。scikit-geometry 就是在这个缝隙里出现的它把 CGAL 那套经过工业验证的计算几何内核用 Python 包装成一套符合 Python 习惯的 API让你在脚本里直接调用精确谓词、布尔运算、曲线偏移、三角剖分这些能力而不用去啃模板元编程。它适合谁做 CAD 二次开发、机器人路径规划、GIS 数据处理、图形学实验、甚至 3D 打印切片算法的从业者只要你的几何计算开始涉及“精确”和“退化情况”scikit-geometry 就值得放进工具箱。这一章不展开代码先把它的定位、和 shapely 的边界、以及为什么值得你花一个下午跑通它讲清楚。2. 装得上、跑得通scikit-geometry 的环境搭建与最小验证2.1 为什么推荐 conda 而不是 pip 直接装scikit-geometry 底层依赖 CGAL 的 C 库编译链涉及 Boost、GMP、MPFR 这些数学库。用 pip 从源码编译在 Windows 上基本是一场噩梦在 macOS 上也要处理头文件路径。常见做法是用 conda 安装预编译包把二进制依赖一次性解决。我一般会单独建一个环境避免和已有的 shapely、geopandas 环境互相污染。# 创建独立环境Python 版本建议 3.9 到 3.11 conda create -n skgeom python3.10 -y conda activate skgeom # 从 conda-forge 频道安装预编译包会带上 CGAL 依赖 conda install -c conda-forge scikit-geometry -y # 验证安装能打印版本号说明二进制链接没问题 python -c import skgeom; print(skgeom.__version__)这里的关键参数是-c conda-forgescikit-geometry 的官方维护频道就是 conda-forge不要混用 defaults 频道否则容易出现 GMP 库版本冲突。如果你在 Linux 上坚持用 pip需要先apt install libcgal-dev libgmp-dev libmpfr-dev再pip install scikit-geometry但编译时间可能超过十分钟且不保证一次过。2.2 用三个最小例子确认核心能力可用装完之后不要急着上业务代码先用三个小例子确认精确谓词、布尔运算、曲线偏移这三块能力都正常。这三个例子覆盖了 scikit-geometry 最常用的入口。import skgeom as sg from skgeom import Point2, Segment2, Polygon # 例子1精确的方向判定这是 CGAL 精确谓词的直接体现 p Point2(0, 0) q Point2(1, 1) r Point2(2, 2) # orientation 返回 COLLINEAR / LEFT_TURN / RIGHT_TURN print(三点方向:, sg.orientation(p, q, r)) # 例子2两个多边形的布尔交集 poly_a Polygon([Point2(0, 0), Point2(2, 0), Point2(2, 2), Point2(0, 2)]) poly_b Polygon([Point2(1, 1), Point2(3, 1), Point2(3, 3), Point2(1, 3)]) inter sg.boolean_set.intersect(poly_a, poly_b) print(交集面积:, sum(p.area() for p in inter)) # 例子3带圆弧的曲线偏移这是 shapely 做不了的 from skgeom import Curve, Point2 as P curve Curve([P(0, 0), P(1, 0), P(2, 1)]) offset curve.offset(0.1) # 向左偏移 0.1 print(偏移后曲线点数:, len(list(offset.vertices())))逻辑说明orientation是 CGAL 精确谓词的 Python 暴露它不会因为浮点误差把共线判成左转boolean_set.intersect返回的是多边形集合需要遍历求和Curve.offset支持直线段和圆弧段混合的曲线这是 scikit-geometry 相对 shapely 的差异化能力。参数上offset的正负号决定偏移方向具体方向取决于曲线走向建议先用简单直线段试出方向再上复杂曲线。2.3 和 shapely 的边界在哪里很多人会问既然有 shapely为什么还要多装一个库我的经验是把两者按“精度要求”和“几何类型”分工。shapely 适合处理大规模、精度要求不苛刻的矢量数据它的底层 GEOS 用浮点运算速度快、生态好。scikit-geometry 适合处理小规模、但要求精确判定和退化情况鲁棒的计算比如两个多边形共享一条边时的布尔运算、带圆弧的偏移、需要精确排序的扫描线算法。维度shapelyscikit-geometry底层内核GEOS浮点CGAL精确谓词圆弧支持无有布尔运算退化情况可能产生自交精确处理生态集成geopandas 无缝需手动转换大规模数据性能优一般这张表不是让你二选一而是让你在同一个项目里按模块选型。我一般会在数据清洗阶段用 shapely在核心算法阶段用 scikit-geometry中间用坐标列表做转换。3. 把 CGAL 的精确谓词用起来布尔运算与多边形处理的实操3.1 精确谓词为什么能解决“玄学”判定失败浮点几何最让人头疼的就是“玄学”失败明明两个点应该共线orientation却返回左转明明两个多边形应该刚好相切布尔运算却报错。根因是浮点表示无法精确表达大多数有理数累积误差让符号判定翻转。CGAL 的做法是用“过滤式精确谓词”先用浮点快速算一遍如果结果远离零就信任如果接近零就切换到精确有理数运算重新判定。scikit-geometry 把这套机制暴露成orientation、compare_distance这些函数。from skgeom import Point2, orientation, Orientation # 构造一个浮点下容易判错的共线场景 a Point2(0.1, 0.1) b Point2(0.2, 0.2) c Point2(0.3, 0.3) # 精确谓词会返回 COLLINEAR而不是受浮点误差影响 result orientation(a, b, c) print(判定结果:, result) assert result Orientation.COLLINEAR, 精确谓词应判定共线这段代码的价值在于如果你用(b.x-a.x)*(c.y-a.y) - (b.y-a.y)*(c.x-a.x)手算结果可能是 1e-17 这种非零值导致误判。scikit-geometry 的orientation内部做了精确化返回的是可靠的符号。参数上没有可调的它的行为就是精确的你只需要信任它。3.2 多边形布尔运算的三种调用方式和返回值处理scikit-geometry 的布尔运算在boolean_set模块下支持并、交、差、对称差。返回值是一个迭代器里面是若干多边形因为布尔运算的结果可能不连通。这一点和 shapely 返回单个 MultiPolygon 不同需要写循环处理。from skgeom import Point2, Polygon from skgeom import boolean_set # 构造两个有重叠的多边形 poly1 Polygon([Point2(0, 0), Point2(4, 0), Point2(4, 4), Point2(0, 4)]) poly2 Polygon([Point2(2, 2), Point2(6, 2), Point2(6, 6), Point2(2, 6)]) # 并集 union_result list(boolean_set.join(poly1, poly2)) print(并集多边形数:, len(union_result), 总面积:, sum(p.area() for p in union_result)) # 交集 inter_result list(boolean_set.intersect(poly1, poly2)) print(交集面积:, sum(p.area() for p in inter_result)) # 差集poly1 减去 poly2 diff_result list(boolean_set.difference(poly1, poly2)) print(差集面积:, sum(p.area() for p in diff_result)) # 对称差 sym_result list(boolean_set.symmetric_difference(poly1, poly2)) print(对称差面积:, sum(p.area() for p in sym_result))逻辑说明join、intersect、difference、symmetric_difference四个函数签名一致都是两个多边形进、迭代器出。必须用list()消费迭代器否则结果不会计算。面积用p.area()累加注意返回的多边形可能是带洞的area()已经扣除了洞的面积。参数上没有额外选项但要注意输入多边形的顶点顺序scikit-geometry 要求逆时针为外环、顺时针为洞顺序错了结果会异常。3.3 多边形偏移和三角剖分的组合用法多边形偏移在 CAM 路径生成、安全距离计算里很常见。scikit-geometry 提供offset方法但要注意它和 shapely 的buffer语义不同shapely 的 buffer 是向外扩张scikit-geometry 的 offset 是沿边界法向平移正负号决定内外。三角剖分用triangulation模块支持约束三角剖分适合做有限元网格或渲染预处理。from skgeom import Point2, Polygon from skgeom import triangulation # 构造一个凹多边形 poly Polygon([Point2(0, 0), Point2(4, 0), Point2(4, 2), Point2(2, 1), Point2(0, 2)]) # 向外偏移 0.5注意正负号需要根据顶点顺序试 offset_poly poly.offset(0.5) print(偏移后面积:, offset_poly.area()) # 约束三角剖分 tri triangulation.ConstrainedTriangulation() for i in range(len(poly.vertices)): tri.insert(poly.vertices[i]) # 插入约束边 for i in range(len(poly.vertices)): j (i 1) % len(poly.vertices) tri.insert_constraint(poly.vertices[i], poly.vertices[j]) faces list(tri.finite_faces()) print(三角形数量:, len(faces))逻辑说明offset的正负号取决于多边形顶点是逆时针还是顺时针逆时针时正数向外。ConstrainedTriangulation先插入所有顶点再插入约束边最后用finite_faces()取三角形。参数上insert_constraint要求两个顶点已经插入过否则会抛异常。三角剖分的输出是面对象需要自己提取顶点坐标。4. 避坑指南scikit-geometry 落地时最容易翻车的五个点4.1 坑一conda 环境混装导致 GMP 符号冲突现象import skgeom时报undefined symbol: __gmpz_init或类似链接错误。原因环境里同时装了多个来源的 GMP 库动态链接器找到了错误版本。解决新建干净环境只用 conda-forge 频道不要在里面 pip install 其他带 C 扩展的几何库。如果必须混用用conda list | grep gmp确认只有一个来源。4.2 坑二多边形顶点顺序错误导致布尔运算结果为空现象两个明显重叠的多边形做交集返回空迭代器。原因其中一个多边形的顶点是顺时针被解释成带洞的“负面积”多边形。解决统一用逆时针顺序构造外环可以用鞋带公式检查面积正负。scikit-geometry 不自动纠正顶点顺序这是和 shapely 的一个行为差异。4.3 坑三offset 方向搞反导致路径跑到多边形内部现象想做向外扩张的路径结果偏移后的多边形缩进去了。原因offset的方向由顶点顺序和正负号共同决定没有“向外”这个绝对语义。解决先用一个正方形测试逆时针顶点配正数看是向外还是向内记下规律再上业务数据。我一般会写一个assert offset_poly.area() poly.area()做方向校验。4.4 坑四大规模数据下性能断崖式下跌现象处理几千个多边形时scikit-geometry 比 shapely 慢一个数量级。原因精确谓词在远离零时虽然走浮点快路径但整体常数比 GEOS 大且 Python 层迭代器开销明显。解决把精确计算限制在真正需要的模块大规模预处理用 shapely 或 numpy 向量化只把退化情况多的核心步骤交给 scikit-geometry。4.5 坑五迭代器未消费导致结果“消失”现象调用boolean_set.intersect后直接打印显示一个生成器对象以为没结果。原因返回的是惰性迭代器不消费就不计算。解决始终用list()包一层或者写for p in result循环。这个坑在新手里出现频率极高因为 shapely 返回的是实体对象习惯很难改。5. 进阶技巧用 scikit-geometry 做精确的扫描线骨架提取5.1 骨架提取为什么需要精确谓词骨架提取skeletonization在形状分析、路径规划里很有用。常见做法是用中轴变换但中轴变换涉及大量“到两条边距离相等”的判定浮点下这些等距点会漂移导致骨架断裂或毛刺。scikit-geometry 的精确谓词能让等距判定稳定配合segment和voronoi模块可以搭一个简化版的骨架提取流程。from skgeom import Point2, Segment2 from skgeom import voronoi import math # 用多边形的边构造线段集合 def polygon_to_segments(vertices): segs [] n len(vertices) for i in range(n): j (i 1) % n segs.append(Segment2(vertices[i], vertices[j])) return segs # 构造一个 L 形多边形的边 verts [Point2(0, 0), Point2(4, 0), Point2(4, 2), Point2(2, 2), Point2(2, 4), Point2(0, 4)] segs polygon_to_segments(verts) # 用线段集构造 Voronoi 图取内部边作为骨架候选 vd voronoi.VoronoiDiagram() for s in segs: vd.insert(s) # 遍历 Voronoi 边筛选落在多边形内部的 skeleton_edges [] for edge in vd.edges(): src edge.source() tgt edge.target() if src is None or tgt is None: continue # 射线边跳过 mid Point2((src.x() tgt.x()) / 2, (src.y() tgt.y()) / 2) # 用射线法判断中点是否在多边形内 if point_in_polygon(mid, verts): skeleton_edges.append((src, tgt)) print(骨架边数量:, len(skeleton_edges))逻辑说明VoronoiDiagram支持插入线段作为站点生成的 Voronoi 边就是到相邻线段等距的轨迹落在多边形内部的部分构成中轴骨架的近似。point_in_polygon需要自己实现射线法这里省略。参数上vd.edges()返回的边可能有无界射线用source()或target()为None判断并跳过。5.2 验证骨架质量的三个指标骨架提取完不能直接信要用指标验证。我一般看三个连通性骨架应该是连通的除非形状本身断开、对称性对称形状的骨架应该对称、端点位置骨架端点应靠近形状的凸角。可以用下面的方式快速检查连通性。# 用并查集检查骨架连通性 parent {} def find(x): while parent[x] ! x: parent[x] parent[parent[x]] x parent[x] return x def union(a, b): ra, rb find(a), find(b) if ra ! rb: parent[ra] rb # 初始化所有端点 for src, tgt in skeleton_edges: key_src (round(src.x(), 6), round(src.y(), 6)) key_tgt (round(tgt.x(), 6), round(tgt.y(), 6)) parent.setdefault(key_src, key_src) parent.setdefault(key_tgt, key_tgt) union(key_src, key_tgt) components len({find(k) for k in parent}) print(连通分量数:, components) # 理想情况是 1逻辑说明把骨架边的端点哈希成整数坐标对用并查集合并最后数连通分量。理想情况下连通分量是 1如果大于 1 说明骨架断裂需要检查 Voronoi 边的筛选阈值。参数上round的精度要和坐标量级匹配坐标在 0 到 10 之间时保留 6 位足够。5.3 我踩过的那个坑最后说一个我自己的血泪经验。早期做骨架提取时我直接用浮点坐标做端点哈希结果因为 Voronoi 顶点有微小抖动同一个交点被哈希成两个键连通分量数虚高我以为是算法错了排查了一整天才发现是哈希精度问题。后来养成习惯凡是涉及精确几何的输出做任何基于坐标的索引或哈希之前先按坐标量级统一 round。这个习惯帮我省下了至少三次类似的排查时间。scikit-geometry 给了你精确的谓词但从精确谓词到你的业务数据结构之间还有一层需要你自己把控的精度纪律。希望帮到你。本文还有配套的精品资源点击获取