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

cp-algorithms 半平面交(Half-plane Intersection):从暴力到 Sort-and-Incremental 的 O(N log N) 完全指南

cp-algorithms 半平面交(Half-plane Intersection):从暴力到 Sort-and-Incremental 的 O(N log N) 完全指南 ★ FEATURED ARTICLE
文档教程知识库【免费下载链接】cp-algorithmsAlgorithm and data structure articles for https://cp-algorithms.com (based on http://e-maxx.ru)项目地址https://gitcode.com/GitHub_Trending/cp/cp-algorithms点击查看免费下载本篇文章基于 cp-algorithms 仓库的 半平面交文档系统讲解计算一组半平面交集的核心问题交集总是凸区域/凸多边形其中任意一点都属于所有半平面我们最终要构造的就是这个凸多边形。文章先给出问题的形式化定义与几何直觉再依次介绍 $O(N^3)$ 暴力法、$O(N^2)$ 增量法最后重点推导并给出Sort-and-Incremental 算法排序增量算法的完整 C 实现以及它在凸多边形求交、平面可见性、二分答案与二维线性规划中的典型应用。读完本篇你将能独立理解并写出该算法并将其应用到各类计算几何竞赛题中。问题定义与基础约定我们有一组数量为 $N$ 的半平面每个半平面指一条直线划分平面后得到的一侧区域含直线本身。所有半平面的交集形成一个凸区域该区域要么有界是一个凸多边形要么为空要么无界。文章的目标是构造这个交集多边形或者判定交集是否为空。在正式讨论算法之前需要统一以下约定除非特别说明$N$ 的定义$N$ 表示给定集合中半平面的数量。半平面的表示方式每条直线用一个“经过点 方向向量”表示即直线上任意一点与直线的方向向量。对于半平面我们约定其允许的区域是方向向量左侧的区域。同时定义半平面的“角度”为其方向向量的极角即与 $x$ 轴正方向的夹角。有界性假设我们假设交集结果总是有界或空的。若需要处理无界的情况只需额外加入 4 个构成一个足够大包围盒bounding box的半平面。无平行假设为了简化问题先假设给定集合中不存在相互平行的半平面。文章末尾会专门讨论如何处理平行半平面的情况。下图直观展示了第 2 条约定的表示方法例如半平面 $y \geq 2x - 2$ 可以表示为点 $P (1, 0)$ 与方向向量 $PQ Q - P (1, 2)$。建议读者预先熟悉基本几何原语点、向量、直线求交。另外仓库中关于 凸包构造 与 凸包技巧 / Li Chao 树 的知识有助于加深理解但并非本算法的必要前置条件。暴力方法$O(N^3)$最直白的思路是枚举所有半平面对共 $\binom{N}{2} O(N^2)$ 对的直线交点然后对每个交点逐一检查它是否落在其余所有半平面内部。每检查一个交点需要 $O(N)$ 次判定因此总时间复杂度为 $O(N^3)$。得到所有通过检查的交点后可以用凸包算法如 Graham 扫描或 Monotone chain重建出交集区域。为什么这样做是正确的因为交集凸多边形的所有顶点都必然是某两条半平面直线的交点而每个这样的顶点显然属于所有半平面。该方法的优点在于简单、容易记忆尤其适合在现场只需要判定交集是否为空的情况缺点是太慢无法应对大多数竞赛题因此需要更快的算法。增量方法$O(N^2)$另一种思路是把半平面逐一增量式地加入交集这等价于用一条直线把当前凸多边形“裁剪” $N$ 次并在每一步剔除冗余部分。具体做法是用线段列表表示当前的凸多边形用新半平面的直线与该多边形各条边求交若直线真正穿过多边形恰好有两个交点把这两交点之间的旧边替换为对应新半平面的新边。这个过程可以做到线性时间。因此从一个大包围盒出发用每个半平面依次裁剪总复杂度即为 $O(N^2)$。虽然比暴力法好很多但每一步都要遍历当前所有的 $O(N)$ 条边仍然浪费。下面的关键观察将把增量思路升级为 $O(N \log N)$ 算法。Sort-and-Incremental 算法$O(N \log N)$该算法最早有据可查的来源是Zeyuan Zhu朱泽园于 2006 年为中国国家队选拔赛撰写的论文New Algorithm for Half-plane Intersection and its Practical Value。本文描述的实现基于同一算法思想但将原论文中分别计算上半区与下半区两个交集的做法合并为一趟用双端队列deque构建全部交集。核心观察交集区域是凸的因此它由若干条半平面线段组成且这些线段按角度排序出现正是半平面在最终交集多边形中的排列顺序。由此得到关键结论如果按角度排序后的顺序增量加入半平面并将它们存于双端队列那么任何冗余的半平面只可能出现在队列的队首或队尾我们只需从两端弹出。为了直观理解假设半平面按角度从 $-\pi$ 到 $\pi$ 排序已构造出前 $k-1$ 个半平面的交集正准备处理第 $k$ 个。由于按角度排序第 $k$ 个半平面必然与第 $k-1$ 个半平面构成一个“凸转角”于是只会发生以下三种情况队尾的若干可能为零个半平面变得冗余需从队尾弹出队首的若干可能为零个半平面变得冗余需从队首弹出处理完以上 1、2 后交集变为空集此时直接报告交集为空并终止算法。“冗余”的定义一个半平面若对交集没有任何贡献即删去它后交集完全不变就称其为冗余。一个带图的示例设当前交集中的半平面集合为 $H {A, B, C, D, E}$相邻半平面直线的交点为 $P {p, q, r, s}$其中 $p A\cap B$、$q B\cap C$、$r C\cap D$、$s D\cap E$。现在想用新半平面 $F$ 切割这个交集可以看到半平面 $F$ 使得队首的 $A$ 与队尾的 $E$ 在交集中都变得冗余。于是我们从队首弹出 $A$、从队尾弹出 $E$并把 $F$ 加入队尾得到新的交集 $H {B, C, D, F}$新的相邻交点集合为 $P {q, r, t, u}$特殊情况的处理平行半平面有了包围盒无界的情况已经被自动解决剩下的棘手情形是平行半平面。平行的两条直线没有交点而我们正是依靠相邻半平面直线的交点来判断冗余与否因此需要特殊处理。平行可分为两种子情况方向相反的平行半平面由于所有半平面按角度排序且包围盒的 4 条边作为半平面参与其中排序后任何一对方向相反的相邻平行半平面之间必然隔着至少一条包围盒半平面因此不会直接相邻。不过在从队尾弹出若干半平面之后方向相反的两个平行半平面仍可能凑到一起——这种情况仅当这两个半平面构成空交集时才会发生因为新加入的半平面会清空整个队列。为避免问题必须手动检查平行性若方向相反直接终止算法并返回空交集。方向相同的平行半平面同一角度这是唯一真正需要常规处理的情形而且非常简单只保留最靠左的那个半平面其余全部删除因为它们完全冗余。算法总流程排序将半平面按角度排序耗时 $O(N \log N)$。增量扫描依次处理每个半平面必要时从双端队列队首、队尾弹出冗余项。由于每个半平面至多入队、出队一次整个扫描阶段总计线性时间。收尾重建扫描结束后计算队列中相邻半平面直线的交点即可得到交集凸多边形同样为线性时间。也可以在步骤 2 中就地保存这些交点以跳过本步但实现上即时计算交点通常更简单。综上总时间复杂度为 $O(N \log N)$。排序是唯一的瓶颈——若输入半平面已按角度有序例如直接由凸多边形各条边生成的半平面算法可以退化为线性时间 $O(N)$。直接实现完整 C 代码以下是可直接使用的实现。首先是基本的点/向量与半平面结构体// Redefine epsilon and infinity as necessary. Be mindful of precision errors. const long double eps 1e-9, inf 1e9; // Basic point/vector struct. struct Point { long double x, y; explicit Point(long double x 0, long double y 0) : x(x), y(y) {} // Addition, subtraction, multiply by constant, dot product, cross product. friend Point operator (const Point p, const Point q) { return Point(p.x q.x, p.y q.y); } friend Point operator - (const Point p, const Point q) { return Point(p.x - q.x, p.y - q.y); } friend Point operator * (const Point p, const long double k) { return Point(p.x * k, p.y * k); } friend long double dot(const Point p, const Point q) { return p.x * q.x p.y * q.y; } friend long double cross(const Point p, const Point q) { return p.x * q.y - p.y * q.x; } }; // Basic half-plane struct. struct Halfplane { // p is a passing point of the line and pq is the direction vector of the line. Point p, pq; long double angle; Halfplane() {} Halfplane(const Point a, const Point b) : p(a), pq(b - a) { angle atan2l(pq.y, pq.x); } // Check if point r is outside this half-plane. // Every half-plane allows the region to the LEFT of its line. bool out(const Point r) { return cross(pq, r - p) -eps; } // Comparator for sorting. bool operator (const Halfplane e) const { return angle e.angle; } // Intersection point of the lines of two half-planes. It is assumed theyre never parallel. friend Point inter(const Halfplane s, const Halfplane t) { long double alpha cross((t.p - s.p), t.pq) / cross(s.pq, t.pq); return s.p (s.pq * alpha); } };然后是算法的核心逻辑// Actual algorithm vectorPoint hp_intersect(vectorHalfplane H) { Point box[4] { // Bounding box in CCW order Point(inf, inf), Point(-inf, inf), Point(-inf, -inf), Point(inf, -inf) }; for(int i 0; i 4; i) { // Add bounding box half-planes. Halfplane aux(box[i], box[(i1) % 4]); H.push_back(aux); } // Sort by angle and start algorithm sort(H.begin(), H.end()); dequeHalfplane dq; int len 0; for(int i 0; i int(H.size()); i) { // Remove from the back of the deque while last half-plane is redundant while (len 1 H[i].out(inter(dq[len-1], dq[len-2]))) { dq.pop_back(); --len; } // Remove from the front of the deque while first half-plane is redundant while (len 1 H[i].out(inter(dq[0], dq[1]))) { dq.pop_front(); --len; } // Special case check: Parallel half-planes if (len 0 fabsl(cross(H[i].pq, dq[len-1].pq)) eps) { // Opposite parallel half-planes that ended up checked against each other. if (dot(H[i].pq, dq[len-1].pq) 0.0) return vectorPoint(); // Same direction half-plane: keep only the leftmost half-plane. if (H[i].out(dq[len-1].p)) { dq.pop_back(); --len; } else continue; } // Add new half-plane dq.push_back(H[i]); len; } // Final cleanup: Check half-planes at the front against the back and vice-versa while (len 2 dq[0].out(inter(dq[len-1], dq[len-2]))) { dq.pop_back(); --len; } while (len 2 dq[len-1].out(inter(dq[0], dq[1]))) { dq.pop_front(); --len; } // Report empty intersection if necessary if (len 3) return vectorPoint(); // Reconstruct the convex polygon from the remaining half-planes. vectorPoint ret(len); for(int i 0; i1 len; i) { ret[i] inter(dq[i], dq[i1]); } ret.back() inter(dq[len-1], dq[0]); return ret; }实现细节与注意事项重复顶点当多个半平面交于同一点时最终多边形可能出现相邻的重复点。这不会影响“交集是否为空”的判断也不影响多边形面积的计算。根据后续任务需求可以用std::unique轻松去重。但注意算法执行期间必须保留重复点这样才能正确处理面积为 0 的交集例如交成单点、线段的情形。建议读者自行构造一些交集退化为单点或线段的微小用例来验证。输入形式为线性约束 $ax by c \leq 0$ 怎么办此时有两种选择改造算法使其直接适配这种约束形式自行实现对应的半平面结构体若熟悉 凸包技巧 则并不困难或者取每条直线上任意 2 个点将其转换成文中使用的“点 方向向量”表示。一般而言推荐直接使用题目给定的表示形式以避免额外的精度损失。精度代码中使用long double与eps 1e-9。跨平台/极端数据下务必根据数值范围调整 epsilon 与无穷大inf的大小如本仓库实现中以inf 1e9构造包围盒并警惕交叉积、点积运算中的浮点误差。应用场景与例题思路许多本可用半平面交解决的问题通常也有替代解法但往往更复杂或更冷门。半平面交一般出现在与多边形主要是凸多边形、平面内可见性以及二维线性规划相关的问题中。以下是几类典型应用。凸多边形求交这是半平面交的经典应用给定 $N$ 个多边形求同时包含在所有多边形内部的区域。因为半平面交的结果是凸多边形反过来任何凸多边形都能表示为若干个半平面的交集每条边对应一个半平面。因此只需为每个多边形生成半平面再对全集求交即可。设 $S$ 为所有多边形边数总和总复杂度为 $O(S \log S)$。理论上也可以用堆合并 $N$ 组半平面后跳过排序步骤得到 $O(S \log N)$但常数因子远劣于直接排序仅对很小的 $N$ 有微弱提速。平面中的可见性Visibility形如“判断平面上某些线段能否从某些点看到”的问题通常可以归结为半平面交。例如给定一个简单多边形不一定凸判断多边形内是否存在一点使得从该点能观察到多边形的全部边界。这就是求多边形的核kernel即星形多边形的可见核。做法是把每条边当作半平面内测方向为允许区域然后求交集交集非空即存在这样的观察点。VasilyevArtem Vasilyev在巴西 ICPC 暑期学校的讲座中还给出过一个更有趣的问题给定平面上按编号排列的点集 $p_1, p_2, \dots, p_n$判断是否存在一个点 $q$站在该处可以从左到右按编号递增的顺序看到所有点。关键观察是能同时看到 $p_i$在 $p_j$ 左侧等价于能看到线段 $p_i p_j$ 的右侧或反向线段 $p_j p_i$ 的左侧。据此为每条线段 $p_i p_{i1}$或反向构造半平面检查整个集合的交集是否为空即可。半平面交 二分答案半平面交的另一个常用场景是充当二分验证器predicate。以 Vasilyev 同场讲座中的问题为例给定一个凸多边形 $P$求能被内接inscribe在其内部的最大半径圆。对半径 $r$ 二分。注意到半径为 $r$ 的圆能内接进 $P$当且仅当存在 $P$ 内一点它到 $P$ 边界所有点的距离都不小于 $r$。验证方式把 $P$ 各条边按逆时针顺序取半平面将它们沿允许区域方向即垂直于方向向量的方向向内平移距离 $r$再检查平移后半平面的交集是否非空多边形收缩后非退化哪怕缩成一个点或一条线段也行。显然若能内接半径 $r$ 的圆则能内接任意更小半径的圆因此二分单调性成立。这里还有个重要的优化点凸多边形各边生成的半平面天然按角度有序排序步骤可以直接跳过。设 $N$ 为多边形顶点数$K$ 为二分迭代次数总复杂度为 $O(NK)$$K$ 取决于答案范围与所需精度。二维线性规划半平面交还可用于求解两个变量的线性规划。所有二维线性约束都可写成 $Ax By C \leq 0$ 的形式不等号方向可变显然它们就是半平面。因此判断一组线性约束是否存在可行解直接用半平面交即可计算出可行域即半平面交后可以回答多个“在约束下最大化/最小化线性函数 $f(x, y)$”的查询每个查询用与 凸包技巧 类似的二分做到 $O(\log N)$。此外还存在一个相当简单的随机化算法可以判定线性约束是否有可行解并优化受约束的线性函数Vasilyev 在讲座中有详细讲解见下文参考资料。练习题以下题目可用于巩固对半平面交的理解难度分级参考原文档经典直接应用Codechef –Animesh decides to settle downCHN02POJ 3130 –How I Mathematician Wonder What You Are!POJ 3335 –Rotating ScoreboardPOJ 1474 –Video SurveillancePOJ 1279 –Art GalleryPOJ 2451 –Uyuws Concert进阶题目POJ 3525 –Most Distant Point from the Sea中等Baekjoon 3903 –Jejus Island同上但数据更强POJ 3384 –Feng Shui中等POJ 1755 –Triathlon中/难DMOJ ccoprep3p3 –Arrow中/难POJ 3968 –Jungle Outpost困难Codeforces Gym 101309 Problem J –Jungle Outpost困难Yandex Contest 2540 Problem F –Asymmetry Value很难需虚拟参赛查看补充第 40 届 Petrozavodsk 编程营2021 冬Day 1 的AlmostFair Cake-Cutting 问题 B。截至原文档撰写时该题处于私有状态仅参赛者可见。参考资料主要来源Zeyuan Zhu 的论文《New Algorithm for Half-plane Intersection and its Practical Value》——本算法的原始出处2006 年中国国家队选拔赛论文可在其 MIT CSAIL 主页 publications 处获取。Artem Vasilyev 的《Brazilian ICPC Summer School 2020》讲座——对半平面交及更多几何主题的精彩讲解。中文优质博客《计算几何基础——半平面交》《半平面交算法详解》《半平面交题目汇总》《半平面交的排序增量法》随机化算法《线性规划与半平面交》系列讲解第 4、5 部分Petr Mitrichev 的博客A half-plane week内含练习题列表中难度最高题目的解答注上述练习题与原文档均为外部在线资源实际提交请以对应评测网站为准。赞分享文档教程知识库【免费下载链接】cp-algorithmsAlgorithm and data structure articles for https://cp-algorithms.com (based on http://e-maxx.ru)项目地址https://gitcode.com/GitHub_Trending/cp/cp-algorithms点击查看免费下载相关推荐cp-algorithms 扫描线法查找相交线段对从 O(n²) 到 O(n log n) 的完整实现与原理剖析cp algorithms 扫描线法查找相交线段对从 O n² 到 O n log n 的完整实现与原理剖析 给定平面上的 $n$ 条线段需要判断其中是否存文档教程知识库cp-algorithms 凸包优化与李超线段树从 DP 加速到 O(n log n) 实战指南cp algorithms 凸包优化与李超线段树从 DP 加速到 O n log n 实战指南 导读 本文讲解 cp algorithms 仓库 src/g文档教程知识库cp-algorithms 分治 DP 优化Divide and Conquer DP从 O(mn²) 到 O(mn log n) 的动态规划递推加速cp algorithms 分治 DP 优化Divide and Conquer DP从 O mn² 到 O mn log n 的动态规划递推加速 分治文档教程知识库上一篇LAMDA 实战10 分钟跑通 Android 一键抓包与自动化控制下一篇DB-GPT 使用指南3 步让数据库听懂人话创作声明:本文部分内容由AI辅助生成(AIGC),仅供参考
阅读完成 · 觉得有帮助?
咨询建站