Перетин півплощин

У цій статті ми обговоримо задачу обчислення перетину набору півплощин. Такий перетин зручно подати у вигляді опуклої області/многокутника, де кожна точка всередині нього лежить також усередині всіх півплощин, — і саме цей многокутник ми намагаємося знайти чи побудувати. Ми дамо початкову інтуїцію для задачі, опишемо підхід зі складністю $O(N \log N)$, відомий як алгоритм «сортування та інкрементальної побудови» (Sort-and-Incremental), і наведемо кілька прикладів застосування цієї техніки.

Наполегливо рекомендуємо читачеві бути знайомим із базовими геометричними примітивами та операціями (точки, вектори, перетин прямих). Крім того, знання про опуклі оболонки або про трюк опуклої оболонки можуть допомогти краще зрозуміти поняття в цій статті, але вони аж ніяк не є обов'язковою передумовою.

Коли підходить цей алгоритм?

• Обмеження задачі — це набір лінійних нерівностей (півплощин), і потрібна їхня спільна допустима область?
• Шукаєш перетин півплощин, а не опуклу оболонку множини точок? (якщо ні → Побудова опуклої оболонки)
• Результуючу область можна вважати обмеженою (або готовий додати обмежувальний прямокутник із 4 півплощин)?

Початкові уточнення та означення

Упродовж усієї статті ми робитимемо деякі припущення (якщо не зазначено інше):

1. Позначимо через $N$ кількість півплощин у заданому наборі.
2. Прямі та півплощини ми подаватимемо однією точкою й одним вектором (будь-яка точка, що лежить на заданій прямій, та напрямний вектор прямої). У випадку півплощин ми вважаємо, що кожна півплощина дозволяє область ліворуч від свого напрямного вектора. Крім того, кутом півплощини ми називаємо полярний кут її напрямного вектора. Приклад дивіться на зображенні нижче.
3. Ми вважатимемо, що результуючий перетин завжди або обмежений, або порожній. Якщо нам потрібно обробити необмежений випадок, ми можемо просто додати 4 півплощини, що задають достатньо великий обмежувальний прямокутник.
4. Для простоти ми вважатимемо, що в заданому наборі немає паралельних півплощин. Ближче до кінця статті ми обговоримо, як працювати з такими випадками.

Півплощину $y \geq 2x - 2$ можна подати точкою $P = (1, 0)$ з напрямним вектором $PQ = Q - P = (1, 2)$

Підхід «в лоб» — $O(N^3)$

Одним із найпрямолінійніших і найочевидніших розв'язків було б обчислити точку перетину прямих для всіх пар півплощин і для кожної точки перевірити, чи лежить вона всередині всіх інших півплощин. Оскільки точок перетину є $O(N^2)$, а для кожної з них потрібно перевірити $O(N)$ півплощин, загальна часова складність становить $O(N^3)$. Саму область перетину потім можна відновити, застосувавши, наприклад, алгоритм побудови опуклої оболонки до набору точок перетину, які увійшли в усі півплощини.

Доволі легко бачити, чому це працює: вершини результуючого опуклого многокутника — це всі точки перетину прямих півплощин, і кожна з цих вершин, очевидно, є частиною всіх півплощин. Головна перевага цього методу в тому, що його легко зрозуміти, запам'ятати й написати «на льоту», якщо вам потрібно лише перевірити, чи порожній перетин. Проте він жахливо повільний і непридатний для більшості задач, тож нам потрібно щось швидше.

Інкрементальний підхід — $O(N^2)$

Інший доволі прямолінійний підхід — будувати перетин півплощин інкрементально, по одній за раз. Цей метод по суті еквівалентний відрізанню опуклого многокутника прямою $N$ разів із вилученням зайвих півплощин на кожному кроці. Щоб це зробити, ми можемо подати опуклий многокутник як список відрізків, а щоб відрізати його півплощиною, ми просто знаходимо точки перетину відрізків із прямою півплощини (точок перетину буде лише дві, якщо пряма коректно перетинає многокутник) і замінюємо всі відрізки між ними новим відрізком, що відповідає півплощині. Оскільки таку процедуру можна реалізувати за лінійний час, ми можемо просто почати з великого обмежувального прямокутника й відрізати його кожною з півплощин, отримавши загальну часову складність $O(N^2)$.

Цей метод — великий крок у правильному напрямку, але доводиться щоразу перебирати $O(N)$ півплощин, що видається марнотратним. Далі ми побачимо, що, зробивши кілька кмітливих спостережень, ідеї, які стоять за цим інкрементальним підходом, можна переробити так, щоб створити алгоритм складності $O(N \log N)$.

Алгоритм «сортування та інкрементальної побудови» — $O(N \log N)$

Найпершим належно задокументованим джерелом цього алгоритму, яке нам вдалося знайти, є дипломна робота Цеюаня Чжу (Zeyuan Zhu) для Chinese Team Selecting Contest під назвою New Algorithm for Half-plane Intersection and its Practical Value від 2006 року. Підхід, який ми опишемо далі, ґрунтується на тому самому алгоритмі, але замість обчислення двох окремих перетинів для нижньої та верхньої частин ми побудуємо все одразу за один прохід, використовуючи дек (двобічну чергу).

Сам алгоритм, як і натякає назва, користується тим фактом, що результуюча область перетину півплощин є опуклою, а отже вона складатиметься з деяких відрізків півплощин у порядку, відсортованому за їхніми кутами. Звідси випливає вирішальне спостереження: якщо ми інкрементально перетинаємо півплощини в порядку, відсортованому за кутом (так, як вони з'явилися б у фінальній, результуючій формі перетину), і зберігаємо їх у двобічній черзі, то нам коли-небудь доведеться видаляти півплощини лише з початку та кінця дека.

Щоб краще уявити цей факт, припустімо, що ми виконуємо описаний раніше інкрементальний підхід на наборі півплощин, відсортованому за кутом (у цьому випадку вважатимемо, що вони відсортовані від $-\pi$ до $\pi$), і припустімо, що ми ось-ось почнемо якийсь довільний $k$-й крок. Це означає, що ми вже побудували перетин перших $k-1$ півплощин. Тепер, оскільки півплощини відсортовані за кутом, то якою б не була $k$-та півплощина, ми можемо бути впевнені, що вона утворить опуклий поворот із $(K-1)$-ю півплощиною. З цієї причини може статися кілька речей:

1. Деякі (можливо, жодна) з півплощин у кінці перетину можуть стати зайвими. У цьому випадку нам потрібно виштовхнути ці тепер уже непотрібні півплощини з кінця дека.
2. Деякі (можливо, жодна) з півплощин на початку можуть стати зайвими. Аналогічно до випадку 1, ми просто виштовхуємо їх з початку дека.
3. Перетин може стати порожнім (після обробки випадків 1 та/або 2). У цьому випадку ми просто повідомляємо, що перетин порожній, і завершуємо алгоритм.

Ми кажемо, що півплощина є «зайвою», якщо вона нічого не додає до перетину. Таку півплощину можна було б видалити, і результуючий перетин зовсім не змінився б.

Ось невеликий приклад з ілюстрацією:

Нехай $H = \{ A, B, C, D, E \}$ — набір півплощин, наявних у перетині на цей момент. Крім того, нехай $P = \{ p, q, r, s \}$ — набір точок перетину сусідніх півплощин у $H$. Тепер припустімо, що ми хочемо перетнути його з півплощиною $F$, як показано на ілюстрації нижче:

Зауважте, що півплощина $F$ робить $A$ та $E$ зайвими в перетині. Тож ми вилучаємо як $A$, так і $E$ з початку та кінця перетину відповідно й додаємо $F$ у кінець. І нарешті отримуємо новий перетин $H = \{ B, C, D, F\}$ з $P = \{ q, r, t, u \}$.

З огляду на все це, у нас є майже все необхідне, щоб реально реалізувати алгоритм, але нам ще потрібно поговорити про деякі особливі випадки. На початку статті ми сказали, що додамо обмежувальний прямокутник, щоб подбати про випадки, коли перетин може бути необмеженим, тож єдиний хитрий випадок, який нам справді потрібно обробити, — це паралельні півплощини. Ми можемо мати два підвипадки: дві півплощини можуть бути паралельні з однаковим напрямком або з протилежним напрямком. Причина, через яку цей випадок потрібно обробляти окремо, полягає в тому, що нам доведеться обчислювати точки перетину прямих півплощин, щоб мати змогу перевірити, чи є півплощина зайвою, а дві паралельні прямі не мають точки перетину, тож нам потрібен особливий спосіб роботи з ними.

Для випадку паралельних півплощин протилежної орієнтації: зауважте, що оскільки ми додаємо обмежувальний прямокутник для роботи з необмеженим випадком, це також розв'язує випадок, коли після сортування ми маємо дві сусідні паралельні півплощини з протилежними напрямками, бо між цими двома обов'язково має бути щонайменше одна з півплощин обмежувального прямокутника (пам'ятаймо, що вони відсортовані за кутом).

• Однак можливо, що після вилучення деяких півплощин з кінця дека дві паралельні півплощини протилежного напрямку опиняться разом. Цей випадок трапляється, конкретно, лише тоді, коли ці дві півплощини утворюють порожній перетин, бо ця остання півплощина призведе до того, що з дека буде вилучено все. Щоб уникнути цієї проблеми, нам доводиться вручну перевіряти наявність паралельних півплощин, і якщо вони мають протилежний напрямок, ми просто негайно зупиняємо алгоритм і повертаємо порожній перетин.

Отже, єдиний випадок, який нам справді потрібно обробити, — це наявність кількох півплощин з однаковим кутом, і виявляється, що цей випадок доволі легко обробити: нам потрібно лише залишити найлівішу півплощину й стерти решту, оскільки вони все одно будуть повністю зайвими. Підсумовуючи, повний алгоритм виглядатиме приблизно так:

1. Ми починаємо з сортування набору півплощин за кутом, що займає $O(N \log N)$ часу.
2. Ми перебиратимемо набір півплощин, і для кожної виконуватимемо інкрементальну процедуру, виштовхуючи елементи з початку та кінця двобічної черги за потреби. Загалом це займе лінійний час, оскільки кожну півплощину можна додати або видалити лише один раз.
3. Наприкінці опуклий многокутник, що є результатом перетину, можна просто отримати, обчисливши точки перетину сусідніх півплощин у деку після завершення процедури. Це також займе лінійний час. Можна також зберігати такі точки під час кроку 2 й цілком пропустити цей крок, але ми вважаємо, що дещо легше (з погляду реалізації) обчислювати їх «на льоту».

Загалом ми досягли часової складності $O(N \log N)$. Оскільки сортування — це явне вузьке місце, алгоритм можна змусити працювати за лінійний час в особливому випадку, коли нам задано півплощини, заздалегідь відсортовані за їхніми кутами (прикладом такого випадку було б отримання півплощин, що задають опуклий многокутник).

Пряма реалізація

Ось приклад прямої реалізації алгоритму з коментарями, що пояснюють більшість частин:

Прості структури для точки/вектора та півплощини:

C++

// За потреби перевизначте epsilon та нескінченність. Зважайте на похибки точності.
const long double eps = 1e-9, inf = 1e9; 

// Базова структура точки/вектора.
struct Point { 

    long double x, y;
    explicit Point(long double x = 0, long double y = 0) : x(x), y(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 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; 
    }
};

// Базова структура півплощини.
struct Halfplane { 

    // 'p' — точка, через яку проходить пряма, а 'pq' — напрямний вектор прямої.
    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);    
    }

    // Перевірка, чи точка 'r' лежить поза цією півплощиною. 
    // Кожна півплощина дозволяє область ЛІВОРУЧ від своєї прямої.
    bool out(const Point& r) { 
        return cross(pq, r - p) < -eps; 
    }

    // Компаратор для сортування. 
    bool operator < (const Halfplane& e) const { 
        return angle < e.angle;
    } 

    // Точка перетину прямих двох півплощин. Вважається, що вони ніколи не паралельні.
    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);
    }
};

Python

import math


# За потреби перевизначте epsilon та нескінченність. Зважайте на похибки точності.
EPS = 1e-9
INF = 1e9


# Базова структура точки/вектора.
class Point:
    __slots__ = ("x", "y")

    def __init__(self, x: float = 0.0, y: float = 0.0):
        self.x = x
        self.y = y

    # Додавання, віднімання, множення на константу.
    def __add__(self, q: "Point") -> "Point":
        return Point(self.x + q.x, self.y + q.y)

    def __sub__(self, q: "Point") -> "Point":
        return Point(self.x - q.x, self.y - q.y)

    def __mul__(self, k: float) -> "Point":
        return Point(self.x * k, self.y * k)


# Скалярний добуток.
def dot(p: Point, q: Point) -> float:
    return p.x * q.x + p.y * q.y


# Векторний добуток.
def cross(p: Point, q: Point) -> float:
    return p.x * q.y - p.y * q.x


# Базова структура півплощини.
class Halfplane:
    __slots__ = ("p", "pq", "angle")

    # 'p' — точка, через яку проходить пряма, а 'pq' — напрямний вектор прямої.
    def __init__(self, a: Point = None, b: Point = None):
        if a is None:
            self.p = Point()
            self.pq = Point()
            self.angle = 0.0
        else:
            self.p = a
            self.pq = b - a
            self.angle = math.atan2(self.pq.y, self.pq.x)

    # Перевірка, чи точка 'r' лежить поза цією півплощиною.
    # Кожна півплощина дозволяє область ЛІВОРУЧ від своєї прямої.
    def out(self, r: Point) -> bool:
        return cross(self.pq, r - self.p) < -EPS

    # Компаратор для сортування (ключ за кутом).
    def __lt__(self, e: "Halfplane") -> bool:
        return self.angle < e.angle


# Точка перетину прямих двох півплощин. Вважається, що вони ніколи не паралельні.
def inter(s: Halfplane, t: Halfplane) -> Point:
    alpha = cross(t.p - s.p, t.pq) / cross(s.pq, t.pq)
    return s.p + s.pq * alpha

TypeScript

// За потреби перевизначте epsilon та нескінченність. Зважайте на похибки точності.
const EPS = 1e-9;
const INF = 1e9;

// Базова структура точки/вектора.
class Point {
    constructor(public x: number = 0, public y: number = 0) {}

    // Додавання, віднімання, множення на константу.
    add(q: Point): Point {
        return new Point(this.x + q.x, this.y + q.y);
    }

    sub(q: Point): Point {
        return new Point(this.x - q.x, this.y - q.y);
    }

    mul(k: number): Point {
        return new Point(this.x * k, this.y * k);
    }
}

// Скалярний добуток.
function dot(p: Point, q: Point): number {
    return p.x * q.x + p.y * q.y;
}

// Векторний добуток.
function cross(p: Point, q: Point): number {
    return p.x * q.y - p.y * q.x;
}

// Базова структура півплощини.
class Halfplane {
    // 'p' — точка, через яку проходить пряма, а 'pq' — напрямний вектор прямої.
    p: Point;
    pq: Point;
    angle: number;

    constructor(a?: Point, b?: Point) {
        if (a === undefined || b === undefined) {
            this.p = new Point();
            this.pq = new Point();
            this.angle = 0;
        } else {
            this.p = a;
            this.pq = b.sub(a);
            this.angle = Math.atan2(this.pq.y, this.pq.x);
        }
    }

    // Перевірка, чи точка 'r' лежить поза цією півплощиною.
    // Кожна півплощина дозволяє область ЛІВОРУЧ від своєї прямої.
    out(r: Point): boolean {
        return cross(this.pq, r.sub(this.p)) < -EPS;
    }
}

// Точка перетину прямих двох півплощин. Вважається, що вони ніколи не паралельні.
function inter(s: Halfplane, t: Halfplane): Point {
    const alpha = cross(t.p.sub(s.p), t.pq) / cross(s.pq, t.pq);
    return s.p.add(s.pq.mul(alpha));
}

Go

package main

import "math"

// За потреби перевизначте epsilon та нескінченність. Зважайте на похибки точності.
const eps = 1e-9
const inf = 1e9

// Базова структура точки/вектора.
type Point struct {
	x, y float64
}

// Додавання, віднімання, множення на константу.
func (p Point) Add(q Point) Point { return Point{p.x + q.x, p.y + q.y} }
func (p Point) Sub(q Point) Point { return Point{p.x - q.x, p.y - q.y} }
func (p Point) Mul(k float64) Point { return Point{p.x * k, p.y * k} }

// Скалярний добуток.
func dot(p, q Point) float64 { return p.x*q.x + p.y*q.y }

// Векторний добуток.
func cross(p, q Point) float64 { return p.x*q.y - p.y*q.x }

// Базова структура півплощини.
// 'p' — точка, через яку проходить пряма, а 'pq' — напрямний вектор прямої.
type Halfplane struct {
	p, pq Point
	angle float64
}

// Конструктор: будує півплощину через дві точки.
func NewHalfplane(a, b Point) Halfplane {
	pq := b.Sub(a)
	return Halfplane{p: a, pq: pq, angle: math.Atan2(pq.y, pq.x)}
}

// Перевірка, чи точка 'r' лежить поза цією півплощиною.
// Кожна півплощина дозволяє область ЛІВОРУЧ від своєї прямої.
func (h Halfplane) Out(r Point) bool {
	return cross(h.pq, r.Sub(h.p)) < -eps
}

// Точка перетину прямих двох півплощин. Вважається, що вони ніколи не паралельні.
func inter(s, t Halfplane) Point {
	alpha := cross(t.p.Sub(s.p), t.pq) / cross(s.pq, t.pq)
	return s.p.Add(s.pq.Mul(alpha))
}

Алгоритм:

C++

// Власне алгоритм
vector<Point> hp_intersect(vector<Halfplane>& H) { 

    Point box[4] = {  // Обмежувальний прямокутник у порядку проти годинникової стрілки
        Point(inf, inf), 
        Point(-inf, inf), 
        Point(-inf, -inf), 
        Point(inf, -inf) 
    };

    for(int i = 0; i<4; i++) { // Додаємо півплощини обмежувального прямокутника.
        Halfplane aux(box[i], box[(i+1) % 4]);
        H.push_back(aux);
    }

    // Сортуємо за кутом і запускаємо алгоритм
    sort(H.begin(), H.end());
    deque<Halfplane> dq;
    int len = 0;
    for(int i = 0; i < int(H.size()); i++) {

        // Видаляємо з кінця дека, поки остання півплощина зайва
        while (len > 1 && H[i].out(inter(dq[len-1], dq[len-2]))) {
            dq.pop_back();
            --len;
        }

        // Видаляємо з початку дека, поки перша півплощина зайва
        while (len > 1 && H[i].out(inter(dq[0], dq[1]))) {
            dq.pop_front();
            --len;
        }
        
        // Перевірка особливого випадку: паралельні півплощини
        if (len > 0 && fabsl(cross(H[i].pq, dq[len-1].pq)) < eps) {
        	// Протилежні паралельні півплощини, які опинилися перевіреними одна проти одної.
        	if (dot(H[i].pq, dq[len-1].pq) < 0.0)
        		return vector<Point>();
        	
        	// Півплощина того самого напрямку: залишаємо лише найлівішу півплощину.
        	if (H[i].out(dq[len-1].p)) {
        		dq.pop_back();
        		--len;
        	}
        	else continue;
        }
        
        // Додаємо нову півплощину
        dq.push_back(H[i]);
        ++len;
    }

    // Фінальне очищення: перевіряємо півплощини на початку проти кінця і навпаки
    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;
    }

    // За потреби повідомляємо про порожній перетин
    if (len < 3) return vector<Point>();

    // Відновлюємо опуклий многокутник із решти півплощин.
    vector<Point> ret(len);
    for(int i = 0; i+1 < len; i++) {
        ret[i] = inter(dq[i], dq[i+1]);
    }
    ret.back() = inter(dq[len-1], dq[0]);
    return ret;
}

Python

from collections import deque


# Власне алгоритм
def hp_intersect(H: list[Halfplane]) -> list[Point]:
    box = [  # Обмежувальний прямокутник у порядку проти годинникової стрілки
        Point(INF, INF),
        Point(-INF, INF),
        Point(-INF, -INF),
        Point(INF, -INF),
    ]

    for i in range(4):  # Додаємо півплощини обмежувального прямокутника.
        H.append(Halfplane(box[i], box[(i + 1) % 4]))

    # Сортуємо за кутом і запускаємо алгоритм
    H.sort()
    dq: deque[Halfplane] = deque()
    length = 0
    for i in range(len(H)):

        # Видаляємо з кінця дека, поки остання півплощина зайва
        while length > 1 and H[i].out(inter(dq[length - 1], dq[length - 2])):
            dq.pop()
            length -= 1

        # Видаляємо з початку дека, поки перша півплощина зайва
        while length > 1 and H[i].out(inter(dq[0], dq[1])):
            dq.popleft()
            length -= 1

        # Перевірка особливого випадку: паралельні півплощини
        if length > 0 and abs(cross(H[i].pq, dq[length - 1].pq)) < EPS:
            # Протилежні паралельні півплощини, які опинилися перевіреними одна проти одної.
            if dot(H[i].pq, dq[length - 1].pq) < 0.0:
                return []

            # Півплощина того самого напрямку: залишаємо лише найлівішу півплощину.
            if H[i].out(dq[length - 1].p):
                dq.pop()
                length -= 1
            else:
                continue

        # Додаємо нову півплощину
        dq.append(H[i])
        length += 1

    # Фінальне очищення: перевіряємо півплощини на початку проти кінця і навпаки
    while length > 2 and dq[0].out(inter(dq[length - 1], dq[length - 2])):
        dq.pop()
        length -= 1

    while length > 2 and dq[length - 1].out(inter(dq[0], dq[1])):
        dq.popleft()
        length -= 1

    # За потреби повідомляємо про порожній перетин
    if length < 3:
        return []

    # Відновлюємо опуклий многокутник із решти півплощин.
    ret = [Point() for _ in range(length)]
    for i in range(length - 1):
        ret[i] = inter(dq[i], dq[i + 1])
    ret[length - 1] = inter(dq[length - 1], dq[0])
    return ret

TypeScript

// Власне алгоритм
function hpIntersect(H: Halfplane[]): Point[] {
    const box: Point[] = [  // Обмежувальний прямокутник у порядку проти годинникової стрілки
        new Point(INF, INF),
        new Point(-INF, INF),
        new Point(-INF, -INF),
        new Point(INF, -INF),
    ];

    for (let i = 0; i < 4; i++) {  // Додаємо півплощини обмежувального прямокутника.
        H.push(new Halfplane(box[i], box[(i + 1) % 4]));
    }

    // Сортуємо за кутом і запускаємо алгоритм
    H.sort((a, b) => a.angle - b.angle);
    // Дек моделюємо масивом з індексами голови/хвоста.
    const dq: Halfplane[] = [];
    let len = 0;
    for (let i = 0; i < H.length; i++) {

        // Видаляємо з кінця дека, поки остання півплощина зайва
        while (len > 1 && H[i].out(inter(dq[len - 1], dq[len - 2]))) {
            dq.pop();
            --len;
        }

        // Видаляємо з початку дека, поки перша півплощина зайва
        while (len > 1 && H[i].out(inter(dq[0], dq[1]))) {
            dq.shift();
            --len;
        }

        // Перевірка особливого випадку: паралельні півплощини
        if (len > 0 && Math.abs(cross(H[i].pq, dq[len - 1].pq)) < EPS) {
            // Протилежні паралельні півплощини, які опинилися перевіреними одна проти одної.
            if (dot(H[i].pq, dq[len - 1].pq) < 0.0)
                return [];

            // Півплощина того самого напрямку: залишаємо лише найлівішу півплощину.
            if (H[i].out(dq[len - 1].p)) {
                dq.pop();
                --len;
            } else {
                continue;
            }
        }

        // Додаємо нову півплощину
        dq.push(H[i]);
        ++len;
    }

    // Фінальне очищення: перевіряємо півплощини на початку проти кінця і навпаки
    while (len > 2 && dq[0].out(inter(dq[len - 1], dq[len - 2]))) {
        dq.pop();
        --len;
    }

    while (len > 2 && dq[len - 1].out(inter(dq[0], dq[1]))) {
        dq.shift();
        --len;
    }

    // За потреби повідомляємо про порожній перетин
    if (len < 3) return [];

    // Відновлюємо опуклий многокутник із решти півплощин.
    const ret: Point[] = new Array(len);
    for (let i = 0; i + 1 < len; i++) {
        ret[i] = inter(dq[i], dq[i + 1]);
    }
    ret[len - 1] = inter(dq[len - 1], dq[0]);
    return ret;
}

Go

import "sort"

// Власне алгоритм
func hpIntersect(H []Halfplane) []Point {
	box := [4]Point{ // Обмежувальний прямокутник у порядку проти годинникової стрілки
		{inf, inf},
		{-inf, inf},
		{-inf, -inf},
		{inf, -inf},
	}

	for i := 0; i < 4; i++ { // Додаємо півплощини обмежувального прямокутника.
		H = append(H, NewHalfplane(box[i], box[(i+1)%4]))
	}

	// Сортуємо за кутом і запускаємо алгоритм
	sort.Slice(H, func(i, j int) bool { return H[i].angle < H[j].angle })
	// Дек моделюємо зрізом (pop_back/pop_front через зрізання).
	dq := make([]Halfplane, 0, len(H))
	length := 0
	for i := 0; i < len(H); i++ {

		// Видаляємо з кінця дека, поки остання півплощина зайва
		for length > 1 && H[i].Out(inter(dq[length-1], dq[length-2])) {
			dq = dq[:length-1]
			length--
		}

		// Видаляємо з початку дека, поки перша півплощина зайва
		for length > 1 && H[i].Out(inter(dq[0], dq[1])) {
			dq = dq[1:]
			length--
		}

		// Перевірка особливого випадку: паралельні півплощини
		if length > 0 && math.Abs(cross(H[i].pq, dq[length-1].pq)) < eps {
			// Протилежні паралельні півплощини, які опинилися перевіреними одна проти одної.
			if dot(H[i].pq, dq[length-1].pq) < 0.0 {
				return []Point{}
			}

			// Півплощина того самого напрямку: залишаємо лише найлівішу півплощину.
			if H[i].Out(dq[length-1].p) {
				dq = dq[:length-1]
				length--
			} else {
				continue
			}
		}

		// Додаємо нову півплощину
		dq = append(dq, H[i])
		length++
	}

	// Фінальне очищення: перевіряємо півплощини на початку проти кінця і навпаки
	for length > 2 && dq[0].Out(inter(dq[length-1], dq[length-2])) {
		dq = dq[:length-1]
		length--
	}

	for length > 2 && dq[length-1].Out(inter(dq[0], dq[1])) {
		dq = dq[1:]
		length--
	}

	// За потреби повідомляємо про порожній перетин
	if length < 3 {
		return []Point{}
	}

	// Відновлюємо опуклий многокутник із решти півплощин.
	ret := make([]Point, length)
	for i := 0; i+1 < length; i++ {
		ret[i] = inter(dq[i], dq[i+1])
	}
	ret[length-1] = inter(dq[length-1], dq[0])
	return ret
}

Обговорення реалізації

Окремо варто зауважити, що у випадку, коли є кілька півплощин, які перетинаються в одній точці, цей алгоритм може повертати повторювані сусідні точки в фінальному многокутнику. Однак це не повинно жодним чином впливати на правильність визначення того, порожній перетин чи ні, і також зовсім не впливає на площу многокутника. Можливо, ви захочете видалити ці дублікати залежно від того, які задачі вам потрібно виконати потім. Зробити це дуже легко за допомогою std::unique. Ми хочемо зберігати повторювані точки під час виконання алгоритму, щоб перетини з площею, що дорівнює нулю, можна було обчислити коректно (наприклад, перетини, що складаються з однієї точки, прямої чи відрізка). Я закликаю читача протестувати кілька невеликих створених вручну випадків, де перетин дає в результаті одну точку чи пряму.

Ще одну річ, про яку варто поговорити, — це що робити, якщо нам задано півплощини у формі лінійного обмеження (наприклад, $ax + by + c \leq 0$). У такому випадку є два варіанти. Ви можете або реалізувати алгоритм із відповідними модифікаціями, щоб працювати з таким поданням (по суті, створити власну структуру півплощини — це має бути доволі прямолінійно, якщо ви знайомі з трюком опуклої оболонки), або перетворити прямі в подання, яке ми використовували в цій статті, узявши будь-які 2 точки кожної прямої. Загалом рекомендується працювати з тим поданням, яке вам задано в задачі, щоб уникнути додаткових проблем із точністю.

Задачі, завдання та застосування

Багато задач, які можна розв'язати за допомогою перетину півплощин, також можна розв'язати без нього, але (зазвичай) складнішими чи менш поширеними підходами. Загалом перетин півплощин може з'являтися в задачах, пов'язаних з многокутниками (переважно опуклими), видимістю на площині та двовимірним лінійним програмуванням. Ось кілька прикладів завдань, які можна розв'язати цією технікою:

Перетин опуклих многокутників

Одне з класичних застосувань перетину півплощин: задано $N$ многокутників, обчислити область, яка міститься всередині всіх многокутників.

Оскільки перетин набору півплощин є опуклим многокутником, ми також можемо подати опуклий многокутник як набір півплощин (кожне ребро многокутника є відрізком півплощини). Згенеруємо ці півплощини для кожного многокутника й обчислимо перетин усього набору. Загальна часова складність становить $O(S \log S)$, де $S$ — загальна кількість сторін усіх многокутників. Теоретично задачу можна також розв'язати за $O(S \log N)$, об'єднавши $N$ наборів півплощин за допомогою купи й потім запустивши алгоритм без кроку сортування, але такий розв'язок має набагато гіршу сталу й дає лише незначний приріст швидкості для дуже малих $N$.

Видимість на площині

Задачі, які вимагають чогось на кшталт «визначити, чи видно деякі відрізки з якоїсь точки (точок) на площині», зазвичай можна сформулювати як задачі перетину півплощин. Візьмемо, наприклад, таке завдання: задано деякий простий многокутник (не обов'язково опуклий), визначити, чи є всередині многокутника якась точка, з якої видно всю межу многокутника. Це також відоме як знаходження ядра многокутника і може бути розв'язане простим перетином півплощин, узявши кожне ребро многокутника як півплощину й потім обчисливши їхній перетин.

Ось пов'язана, цікавіша задача, яку представив Артем Васильєв в одній зі своїх лекцій бразильської літньої школи ICPC: задано набір $p$ точок $p_1, p_2\ \dots \ p_n$ на площині, визначити, чи є якась точка $q$, у якій ви можете стати так, щоб бачити всі точки $p$ зліва направо у зростаючому порядку їхніх індексів.

Таку задачу можна розв'язати, помітивши, що можливість бачити деяку точку $p_i$ ліворуч від $p_j$ — це те саме, що можливість бачити правий бік відрізка від $p_i$ до $p_j$ (або, що еквівалентно, можливість бачити лівий бік відрізка від $p_j$ до $p_i$). З огляду на це, ми можемо просто створити півплощину для кожного відрізка $p_i p_{i+1}$ (або $p_{i+1} p_i$ залежно від обраної вами орієнтації) і перевірити, чи порожній перетин усього набору.

Перетин півплощин із бінарним пошуком

Інше поширене застосування — використання перетину півплощин як інструмента для перевірки предиката процедури бінарного пошуку. Ось приклад такої задачі, також представлений Артемом Васильєвим у тій самій лекції, що згадувалася раніше: задано опуклий многокутник $P$, знайти найбільше коло, яке можна вписати в нього.

Замість того щоб шукати якийсь розв'язок у замкненій формі, докучливі формули чи маловідомі алгоритмічні розв'язки, давайте натомість спробуємо застосувати бінарний пошук за відповіддю. Зауважте, що для деякого фіксованого $r$ коло радіуса $r$ можна вписати в $P$ лише тоді, коли всередині $P$ існує деяка точка, відстань від якої до всіх точок межі $P$ більша або дорівнює $r$. Цю умову можна перевірити, «стиснувши» многокутник усередину на відстань $r$ і перевіривши, що многокутник залишається невиродженим (або сам є точкою/відрізком). Таку процедуру можна змоделювати, узявши півплощини сторін многокутника в порядку проти годинникової стрілки, зсунувши кожну з них на відстань $r$ у напрямку області, яку вони дозволяють (тобто ортогонально до напрямного вектора півплощини), і перевіривши, чи перетин не порожній.

Очевидно, що якщо ми можемо вписати коло радіуса $r$, то можемо вписати й будь-яке інше коло меншого радіуса, ніж $r$. Тож ми можемо виконати бінарний пошук за радіусом $r$ і перевіряти кожен крок за допомогою перетину півплощин. Також зауважте, що півплощини опуклого многокутника вже відсортовані за кутом, тож крок сортування в алгоритмі можна пропустити. Таким чином ми отримуємо загальну часову складність $O(NK)$, де $N$ — кількість вершин многокутника, а $K$ — кількість ітерацій бінарного пошуку (фактичне значення залежатиме від діапазону можливих відповідей і бажаної точності).

Двовимірне лінійне програмування

Ще одне застосування перетину півплощин — лінійне програмування з двома змінними. Усі лінійні обмеження для двох змінних можна виразити у формі $Ax + By + C \leq 0$ (знак нерівності може відрізнятися). Очевидно, що це просто півплощини, тож перевірку того, чи існує допустимий розв'язок для набору лінійних обмежень, можна виконати за допомогою перетину півплощин. Крім того, для заданого набору лінійних обмежень можна обчислити область допустимих розв'язків (тобто перетин півплощин) і потім відповідати на кілька запитів максимізації/мінімізації деякої лінійної функції $f(x, y)$ за умови обмежень за $O(\log N)$ на запит, використовуючи бінарний пошук (дуже схоже на трюк опуклої оболонки).

Варто згадати, що існує також доволі простий рандомізований алгоритм, який може перевірити, чи має набір лінійних обмежень допустимий розв'язок, і максимізувати/мінімізувати деяку лінійну функцію за умови заданих обмежень. Цей рандомізований алгоритм також гарно пояснив Артем Васильєв у згаданій раніше лекції. Ось кілька додаткових ресурсів про нього, якщо читачеві цікаво: CG - Lecture 4, parts 4 and 5 та блог Петра Митричева (який містить розв'язок найскладнішої задачі зі списку задач для практики нижче).

Задачі для практики

Класичні задачі, пряме застосування

• Codechef - Animesh decides to settle down
• POJ - How I mathematician Wonder What You Are!
• POJ - Rotating Scoreboard
• POJ - Video Surveillance
• POJ - Art Gallery
• POJ - Uyuw's Concert

Складніші задачі

• POJ - Most Distant Point from the Sea - Medium
• Baekjoon - Jeju's Island - Same as above but seemingly stronger test cases
• POJ - Feng Shui - Medium
• POJ - Triathlon - Medium/hard
• DMOJ - Arrow - Medium/hard
• POJ - Jungle Outpost - Hard
• Codeforces - Jungle Outpost (alternative link, problem J) - Hard
• Yandex - Asymmetry Value (need virtual contest to see, problem F) - Very Hard

Додаткові задачі

• 40th Petrozavodsk Programming Camp, Winter 2021 - Day 1: Jagiellonian U Contest, Grand Prix of Krakow - Problem B: (Almost) Fair Cake-Cutting. На момент написання статті ця задача була приватною й доступною лише для учасників Programming Camp.

Джерела, бібліографія та інші ресурси

Основні джерела

• New Algorithm for Half-plane Intersection and its Practical Value. Оригінальна стаття про алгоритм.
• Artem Vasilyev's Brazilian ICPC Summer School 2020 lecture. Чудова лекція про перетин півплощин. Також охоплює інші теми геометрії.

Гарні блоги (китайською)

• Fundamentals of Computational Geometry - Intersection of Half-planes.
• Detailed introduction to the half-plane intersection algorithm.
• Summary of Half-plane intersection problems.
• Sorting incremental method of half-plane intersection.

Рандомізований алгоритм

• Linear Programming and Half-Plane intersection - Parts 4 and 5.
• Petr Mitrichev's Blog: A half-plane week.