每日算法 — 使用java实现半平面交:计算几何中的区域求交算法

半平面交(Half-plane Intersection, HPI)是计算几何中的经典问题。给定一组由有向直线定义的半平面,求这些半平面的交集区域。该算法在可见性计算、机器人路径规划、线性规划等领域有广泛应用。本文用Java实现基于排序的O(n log n)半平面交算法,讲解极角排序、双端队列维护与交点计算。

问题定义与几何基础

半平面:由直线 ax + by + c = 0 将平面分为两部分,有向直线左侧(或右侧)的区域即为一个半平面。

半平面交:多个半平面的交集可能是一个凸多边形、凸无界区域、直线、线段、点或空集。

关键几何工具

/**
 * 二维点/向量
 */
class Point {
    double x, y;
    Point(double x, double y) { this.x = x; this.y = y; }
}

/**
 * 有向直线:点p指向点q,左侧为保留的半平面
 */
class Line {
    Point p, q;   // 直线上两点
    double angle; // 极角,用于排序

    Line(Point p, Point q) {
        this.p = p;
        this.q = q;
        this.angle = Math.atan2(q.y - p.y, q.x - p.x);
    }
}

算法核心思想

基于排序增量法(S&I算法):

  1. 极角排序:将所有半平面的边界直线按极角排序
  2. 双端队列维护:用双端队列维护当前半平面交的边界直线
  3. 冗余判断:新加入的直线可能导致队列首尾直线不再必要,及时弹出
  4. 交点计算:相邻直线的交点构成多边形的顶点

核心观察:半平面的交集是凸区域,其边界由若干直线段组成,这些直线在原直线集合中按极角有序排列。

Java完整实现

import java.util.*;

/**
 * 半平面交(Half-plane Intersection)
 * 基于排序增量法,时间复杂度 O(n log n)
 */
public class HalfplaneIntersection {

    static final double EPS = 1e-9;

    /**
     * 二维点/向量,支持基本向量运算
     */
    static class Point {
        double x, y;
        Point(double x, double y) { this.x = x; this.y = y; }
        Point() { this(0, 0); }

        Point add(Point p) { return new Point(x + p.x, y + p.y); }
        Point sub(Point p) { return new Point(x - p.x, y - p.y); }
        Point mul(double k) { return new Point(x * k, y * k); }

        double cross(Point p) { return x * p.y - y * p.x; }  // 叉积
        double dot(Point p) { return x * p.x + y * p.y; }     // 点积
        double len() { return Math.sqrt(x * x + y * y); }
    }

    /**
     * 有向直线:从点p到点q,保留的是直线的左侧半平面
     */
    static class Line implements Comparable<Line> {
        Point p, q;      // 直线上两点
        double angle;    // 极角

        Line(Point p, Point q) {
            this.p = p;
            this.q = q;
            this.angle = Math.atan2(q.y - p.y, q.x - p.x);
        }

        @Override
        public int compareTo(Line other) {
            // 按极角排序,极角相同的保留最靠左的
            if (Math.abs(angle - other.angle) < EPS) {
                // 极角相同,保留左侧更严格的(点在对方左侧)
                return onLeft(other.p) ? -1 : 1;
            }
            return Double.compare(angle, other.angle);
        }

        /**
         * 判断点是否在直线左侧(包含在直线上)
         */
        boolean onLeft(Point r) {
            return q.sub(p).cross(r.sub(p)) > -EPS;
        }
    }

    /**
     * 计算两条直线的交点
     * 假设两条直线不平行
     */
    static Point intersection(Line a, Line b) {
        Point u = a.p.sub(b.p);
        Point v = a.q.sub(a.p);
        Point w = b.q.sub(b.p);
        double t = w.cross(u) / v.cross(w);
        return a.p.add(v.mul(t));
    }

    /**
     * 平行判断
     */
    static boolean parallel(Line a, Line b) {
        return Math.abs(a.q.sub(a.p).cross(b.q.sub(b.p))) < EPS;
    }

    /**
     * 半平面交主算法
     * @param lines 输入的有向直线数组(左侧为保留区域)
     * @return 交集多边形的顶点列表(逆时针顺序),空集返回空列表
     */
    public static List<Point> halfplaneIntersection(Line[] lines) {
        int n = lines.length;

        // 1. 按极角排序,去除极角相同的直线(保留最严格的)
        Arrays.sort(lines);
        Line[] filtered = new Line[n];
        int m = 0;
        for (int i = 0; i < n; i++) {
            // 跳过与上一条极角相同的直线
            if (i > 0 && Math.abs(lines[i].angle - lines[i - 1].angle) < EPS) {
                continue;
            }
            filtered[m++] = lines[i];
        }

        // 2. 双端队列维护半平面交
        Line[] dq = new Line[m + 5];
        Point[] points = new Point[m + 5]; // 相邻直线交点
        int head = 0, tail = 0;

        for (int i = 0; i < m; i++) {
            Line cur = filtered[i];

            // 队尾冗余判断:新直线使队尾两条直线交点不在新半平面内
            while (tail - head > 1 && !cur.onLeft(points[tail - 1])) {
                tail--;
            }

            // 队首冗余判断
            while (tail - head > 1 && !cur.onLeft(points[head + 1])) {
                head++;
            }

            // 如果新直线与队尾直线平行且方向相同,则两者取更严格的
            if (tail > head && parallel(dq[tail - 1], cur)) {
                if (!cur.onLeft(dq[tail - 1].p)) {
                    // 当前直线更严格,替换队尾
                    tail--;
                } else {
                    // 队尾更严格,跳过当前
                    continue;
                }
            }

            dq[tail++] = cur;

            // 计算新加入直线与队尾前一条直线的交点
            if (tail - head > 1) {
                points[tail - 1] = intersection(dq[tail - 2], dq[tail - 1]);
            }
        }

        // 3. 最终处理:队首和队尾可能仍有冗余
        while (tail - head > 1 && !dq[head].onLeft(points[tail - 1])) {
            tail--;
        }
        while (tail - head > 1 && !dq[tail - 1].onLeft(points[head + 1])) {
            head++;
        }

        // 4. 收集多边形顶点
        if (tail - head <= 2) {
            return new ArrayList<>(); // 空集或面积为零
        }

        List<Point> polygon = new ArrayList<>();
        for (int i = head + 1; i < tail; i++) {
            polygon.add(points[i]);
        }
        // 队首和队尾的交点
        polygon.add(intersection(dq[tail - 1], dq[head]));

        return polygon;
    }

    /**
     * 计算多边形面积(鞋带公式)
     */
    public static double polygonArea(List<Point> polygon) {
        if (polygon.size() < 3) return 0;
        double area = 0;
        int n = polygon.size();
        for (int i = 0; i < n; i++) {
            int j = (i + 1) % n;
            area += polygon.get(i).x * polygon.get(j).y;
            area -= polygon.get(j).x * polygon.get(i).y;
        }
        return Math.abs(area) / 2.0;
    }

    /**
     * 将一般不等式 ax + by <= c 转换为有向直线
     * 保留的是不等式表示的半平面
     */
    static Line fromInequality(double a, double b, double c) {
        // ax + by <= c
        // 找直线上两点
        Point p, q;
        if (Math.abs(b) < EPS) {
            // 垂直线: x = c/a
            double x = c / a;
            p = new Point(x, 0);
            q = new Point(x, 1);
        } else {
            // y = (c - ax) / b
            p = new Point(0, c / b);
            q = new Point(1, (c - a) / b);
        }

        // 确保保留的是 ax + by <= c 的区域
        // 测试点(0,0):若 a*0 + b*0 <= c,即 c >= 0,则(0,0)在区域内
        // 需要确保有向直线的左侧是满足不等式的区域
        Line line = new Line(p, q);
        Point test = new Point(0, 0);
        boolean zeroInRegion = a * 0 + b * 0 <= c + EPS;
        boolean zeroOnLeft = line.onLeft(test);

        if (zeroInRegion != zeroOnLeft) {
            // 方向反了,交换p和q
            return new Line(q, p);
        }
        return line;
    }

    public static void main(String[] args) {
        // 示例1:求一个凸多边形(正方形)
        System.out.println("=== 示例1:正方形半平面交 ===");
        Line[] square = {
            new Line(new Point(0, 0), new Point(1, 0)),   // y >= 0(上方)
            new Line(new Point(1, 0), new Point(1, 1)),   // x <= 1(左侧)
            new Line(new Point(1, 1), new Point(0, 1)),   // y <= 1(下方)
            new Line(new Point(0, 1), new Point(0, 0))    // x >= 0(右侧)
        };

        // 注意:Line(p,q)保留的是p->q的左侧
        // 需要调整方向使内部在左侧
        Line[] squareFixed = {
            new Line(new Point(1, 0), new Point(0, 0)),   // y >= 0
            new Line(new Point(1, 1), new Point(1, 0)),   // x <= 1
            new Line(new Point(0, 1), new Point(1, 1)),   // y <= 1
            new Line(new Point(0, 0), new Point(0, 1))    // x >= 0
        };

        List<Point> result1 = halfplaneIntersection(squareFixed);
        System.out.println("交点数量: " + result1.size());
        for (Point p : result1) {
            System.out.printf("(%.2f, %.2f)%n", p.x, p.y);
        }
        System.out.printf("面积: %.4f%n", polygonArea(result1));

        // 示例2:使用不等式形式定义半平面
        System.out.println("\n=== 示例2:不等式形式 ===");
        // 定义三角形区域:
        // x >= 0, y >= 0, x + y <= 1
        Line[] triangle = {
            fromInequality(-1, 0, 0),   // -x <= 0 即 x >= 0
            fromInequality(0, -1, 0),   // -y <= 0 即 y >= 0
            fromInequality(1, 1, 1)     // x + y <= 1
        };

        List<Point> result2 = halfplaneIntersection(triangle);
        System.out.println("三角形顶点:");
        for (Point p : result2) {
            System.out.printf("(%.4f, %.4f)%n", p.x, p.y);
        }
        System.out.printf("三角形面积: %.4f%n", polygonArea(result2));

        // 示例3:无解情况(矛盾的半平面)
        System.out.println("\n=== 示例3:无解情况 ===");
        Line[] impossible = {
            new Line(new Point(0, 0), new Point(1, 0)),  // y >= 0
            new Line(new Point(0, 1), new Point(1, 1))   // y <= 0(矛盾,内部方向需调整)
        };
        // 实际上需要正确设置方向...
        Line[] impossibleFixed = {
            new Line(new Point(1, 0), new Point(0, 0)),  // y >= 0
            new Line(new Point(1, -1), new Point(0, -1)) // y <= -1(与y>=0矛盾)
        };
        List<Point> result3 = halfplaneIntersection(impossibleFixed);
        System.out.println("交点数量(应为0): " + result3.size());
    }
}

算法正确性分析

双端队列的维护基于以下关键观察:

  1. 极角单调性:按极角排序后,相邻直线在凸多边形上也是相邻的
  2. 冗余判定:如果队列中某条直线的两个相邻交点都在新半平面外,则该直线不再构成边界
  3. 凸性保证:半平面的交集是凸集,其边界由按极角有序的直线段组成

复杂度分析

步骤 时间复杂度 说明
极角排序 O(n log n) 主要开销
双端队列处理 O(n) 每条直线入队出队各一次
交点计算 O(n) 相邻直线求交
总复杂度 O(n log n) 由排序决定
空间复杂度 O(n) 存储直线和交点

应用场景

半平面交在多个领域有重要应用:

  • 线性规划的可行域计算:求约束条件的交集
  • 机器人可见性:计算光源或摄像头的可见区域
  • 运动规划:计算机器人可移动的安全区域
  • 计算机图形学:裁剪多边形的Weiler-Atherton算法的辅助工具
  • 核(Kernel)计算:计算点集的核,即可以看到所有点的区域

注意事项

  1. 精度问题:计算几何中浮点精度至关重要,使用EPS进行相等判断
  2. 边界情况:空集、无界区域、退化(直线/点)需要特殊处理
  3. 方向一致性:确保所有半平面的”保留侧”定义一致(本文统一使用左侧)
  4. 无界区域处理:本文实现返回有界多边形,处理无界区域需要添加一个大包围盒

掌握半平面交后,可以解决大量计算几何中的区域求交问题,是算法竞赛中的高级但实用的工具。