C语言半平面交算法详解（从零开始掌握计算几何中的半平面交实现）

为什么需要半平面交？

半平面交是许多高级几何问题的基础。例如：

• 求多个不等式约束下的可行区域（线性规划）
• 构造多边形的核（Kernel of a Polygon）
• 光线投射与视野计算

掌握计算几何算法中的半平面交，能让你在算法竞赛或工程开发中游刃有余。

算法思路概述

最经典的半平面交算法基于双端队列（deque），其核心思想如下：

1. 将所有半平面按极角排序（通常使用向量的方向角）
2. 依次将半平面加入双端队列
3. 每次加入后，检查队首或队尾是否被新加入的半平面“切掉”，若是则弹出
4. 最后，检查队首是否被队尾“切掉”，进行闭环处理

整个过程时间复杂度为 O(n log n)，主要耗时在排序阶段。

C语言实现步骤

下面我们用 C 语言逐步实现一个简单的半平面交算法。为简化问题，我们假设输入的半平面不会退化（如平行、重合等边界情况暂不处理）。

1. 定义基本结构

#include <stdio.h>#include <stdlib.h>#include <math.h>#define EPS 1e-8#define MAXN 1005typedef struct {    double x, y;} Point;typedef struct {    Point p;      // 直线上一点    Point v;      // 方向向量（指向半平面内部）    double ang;   // 极角} HalfPlane;// 向量叉积double cross(Point a, Point b) {    return a.x * b.y - a.y * b.x;}// 向量减法Point sub(Point a, Point b) {    Point res = {a.x - b.x, a.y - b.y};    return res;}// 计算两直线交点Point intersect(HalfPlane a, HalfPlane b) {    double t = cross(sub(b.p, a.p), b.v) / cross(a.v, b.v);    Point res = {a.p.x + a.v.x * t, a.p.y + a.v.y * t};    return res;}// 比较函数：按极角排序int cmp_halfplane(const void *a, const void *b) {    HalfPlane *ha = (HalfPlane *)a;    HalfPlane *hb = (HalfPlane *)b;    if (fabs(ha->ang - hb->ang) < EPS)        return cross(ha->v, sub(hb->p, ha->p)) < 0 ? -1 : 1;    return ha->ang < hb->ang ? -1 : 1;}

2. 半平面交主函数

// 判断点 p 是否在半平面 h 的左侧（含边界）int on_left(Point p, HalfPlane h) {    return cross(h.v, sub(p, h.p)) > -EPS;}// 半平面交算法int half_plane_intersection(HalfPlane *hp, int n, Point *poly) {    // 按极角排序    qsort(hp, n, sizeof(HalfPlane), cmp_halfplane);    // 双端队列    HalfPlane dq[MAXN];    int head = 0, tail = 0;    dq[tail++] = hp[0];    for (int i = 1; i < n; i++) {        // 移除队尾无效半平面        while (head + 1 < tail && !on_left(intersect(dq[tail-2], dq[tail-1]), hp[i]))            tail--;        // 移除队首无效半平面        while (head + 1 < tail && !on_left(intersect(dq[head], dq[head+1]), hp[i]))            head++;        dq[tail++] = hp[i];    }    // 闭环检查    while (head + 1 < tail && !on_left(intersect(dq[tail-2], dq[tail-1]), dq[head]))        tail--;    while (head + 1 < tail && !on_left(intersect(dq[head], dq[head+1]), dq[tail-1]))        head++;    // 构建结果多边形    int m = 0;    for (int i = head; i < tail; i++) {        poly[m++] = intersect(dq[i], dq[(i+1) % (tail - head) + head]);    }    return m; // 返回交集多边形的顶点数}

3. 使用示例

int main() {    HalfPlane hp[4];    // 定义四个半平面：构成一个正方形 [0,2]x[0,2]    hp[0] = (HalfPlane){{0, 0}, {0, 1}, atan2(0, 1)};     // x >= 0    hp[1] = (HalfPlane){{2, 0}, {0, 1}, atan2(0, -1)};    // x <= 2    hp[2] = (HalfPlane){{0, 0}, {1, 0}, atan2(1, 0)};     // y >= 0    hp[3] = (HalfPlane){{0, 2}, {1, 0}, atan2(-1, 0)};    // y <= 2    // 计算极角    for (int i = 0; i < 4; i++) {        hp[i].ang = atan2(hp[i].v.y, hp[i].v.x);    }    Point result[MAXN];    int m = half_plane_intersection(hp, 4, result);    printf("交集多边形有 %d 个顶点:\n", m);    for (int i = 0; i < m; i++) {        printf("(%.2f, %.2f)\n", result[i].x, result[i].y);    }    return 0;}

注意事项与优化

上述代码是教学用途的简化版本。在实际项目中，你可能需要处理以下问题：

• 精度问题：使用 EPS 控制浮点误差
• 平行半平面：可能导致除零错误，需提前判断
• 空交集：算法应能检测并返回空结果
• 性能优化：对于大规模数据，可考虑增量式或分治法

通过不断练习和调试，你将能熟练掌握C语言几何编程中的这一核心技巧。

总结

本文详细介绍了C语言半平面交算法的原理与实现，从基础概念到完整代码，帮助初学者理解这一计算几何算法的核心思想。无论你是准备算法竞赛，还是从事图形学开发，掌握半平面交实现都将为你打下坚实基础。

希望这篇教程能助你在C语言几何编程的道路上更进一步！