Вычислительная геометрия: Алгоритм Грэхема, Алгоритм Джарвиса, Проверка принадлежности точки многоугольнику

Алгоритм Грэхема

Выпуклая оболочка множества точек — наименьший выпуклый многоугольник, содержащий все эти точки. Алгоритм Грэхема строит её за $O(n\log n)$; основную стоимость составляет сортировка.

Опорная точка $p_0$ — самая нижняя из точек, а при равенстве $y$ — самая левая. Она заведомо лежит на оболочке. Остальные точки сортируем по полярному углу относительно $p_0$. Все углы лежат в одной полуплоскости, поэтому порядок задаёт знак векторного произведения: точка $a$ идёт раньше $b$, если $(a-p_0)\times(b-p_0)$ больше нуля. Точки на одном луче сортируем по расстоянию до $p_0$.

Затем проходим отсортированные точки и держим стек вершин. Перед добавлением новой точки $p$ снимаем с вершины стека точки, которые не дают строгого левого поворота: пока для двух последних вершин $a$ и $b$ выполнено $(b-a)\times(p-a)\le 0$, вершину $b$ удаляем. Так убираются вогнутости и точки на сторонах оболочки.

Все вычисления целочисленные: векторное произведение состоит только из сложений и умножений, без углов и корней. Повторы точек удаляются заранее.

Реализация на Python

# Алгоритм Грэхема, O(n log n)
from functools import cmp_to_key


def cross(o, a, b):
    """Векторное произведение (a - o) x (b - o): больше нуля, если поворот влево."""
    return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])


def dist2(a, b):
    return (a[0] - b[0]) ** 2 + (a[1] - b[1]) ** 2


def graham(points):
    """Выпуклая оболочка против часовой стрелки, начиная с нижней левой точки."""
    # Сортируем по (y, x) и удаляем повторы: первая точка — опорная
    pts = sorted(set(points), key=lambda p: (p[1], p[0]))
    if len(pts) <= 2:
        return pts
    pivot, rest = pts[0], pts[1:]

    def by_angle(a, b):
        # a раньше b, если поворот pivot -> a -> b влево; при равном угле ближняя раньше
        turn = cross(pivot, a, b)
        if turn != 0:
            return -turn
        return dist2(pivot, a) - dist2(pivot, b)

    rest.sort(key=cmp_to_key(by_angle))
    hull = [pivot]
    for p in rest:
        # Удаляем вершины, которые не дают строгого левого поворота
        while len(hull) >= 2 and cross(hull[-2], hull[-1], p) <= 0:
            hull.pop()
        hull.append(p)
    return hull


if __name__ == "__main__":
    print(graham([(0, 0), (4, 0), (4, 4), (0, 4), (2, 2), (2, 0), (1, 3), (0, 0)]))

Реализация на C++

// Алгоритм Грэхема, O(n log n)
#include <algorithm>
#include <cstdio>
#include <vector>

struct Point {
    long long x;
    long long y;
};

bool operator==(const Point& a, const Point& b) { return a.x == b.x && a.y == b.y; }

// Векторное произведение (a - o) x (b - o): больше нуля, если поворот влево
long long cross(const Point& o, const Point& a, const Point& b) {
    return (a.x - o.x) * (b.y - o.y) - (a.y - o.y) * (b.x - o.x);
}

long long dist2(const Point& a, const Point& b) {
    return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}

// Выпуклая оболочка против часовой стрелки, начиная с нижней левой точки
std::vector<Point> graham(std::vector<Point> pts) {
    // Сортируем по (y, x) и удаляем повторы: первая точка — опорная
    std::sort(pts.begin(), pts.end(), [](const Point& a, const Point& b) {
        return a.y != b.y ? a.y < b.y : a.x < b.x;
    });
    pts.erase(std::unique(pts.begin(), pts.end()), pts.end());
    if (pts.size() <= 2) return pts;

    // Остальные точки — по полярному углу вокруг опорной; при равном угле ближняя раньше
    const Point pivot = pts[0];
    std::sort(pts.begin() + 1, pts.end(), [&pivot](const Point& a, const Point& b) {
        long long turn = cross(pivot, a, b);
        return turn != 0 ? turn > 0 : dist2(pivot, a) < dist2(pivot, b);
    });

    std::vector<Point> hull;
    for (const Point& p : pts) {
        // Удаляем вершины, которые не дают строгого левого поворота
        while (hull.size() >= 2 && cross(hull[hull.size() - 2], hull.back(), p) <= 0) {
            hull.pop_back();
        }
        hull.push_back(p);
    }
    return hull;
}

int main() {
    const std::vector<Point> pts = {{0, 0}, {4, 0}, {4, 4}, {0, 4}, {2, 2}, {2, 0}, {1, 3}, {0, 0}};
    for (const Point& p : graham(pts)) {
        std::printf("(%lld, %lld)\n", p.x, p.y);
    }
    return 0;
}

Реализация на Rust

// Алгоритм Грэхема, O(n log n)
type Point = (i64, i64);

// Векторное произведение (a - o) x (b - o): больше нуля, если поворот влево
fn cross(o: Point, a: Point, b: Point) -> i64 {
    (a.0 - o.0) * (b.1 - o.1) - (a.1 - o.1) * (b.0 - o.0)
}

fn dist2(a: Point, b: Point) -> i64 {
    (a.0 - b.0).pow(2) + (a.1 - b.1).pow(2)
}

// Выпуклая оболочка против часовой стрелки, начиная с нижней левой точки
pub fn graham(mut pts: Vec<Point>) -> Vec<Point> {
    // Сортируем по (y, x) и удаляем повторы: первая точка — опорная
    pts.sort_by_key(|p| (p.1, p.0));
    pts.dedup();
    if pts.len() <= 2 {
        return pts;
    }
    // Остальные точки — по полярному углу вокруг опорной; при равном угле ближняя раньше
    let pivot = pts[0];
    pts[1..].sort_by(|&a, &b| {
        cross(pivot, a, b).cmp(&0).reverse().then(dist2(pivot, a).cmp(&dist2(pivot, b)))
    });

    let mut hull: Vec<Point> = Vec::new();
    for p in pts {
        // Удаляем вершины, которые не дают строгого левого поворота
        while hull.len() >= 2 && cross(hull[hull.len() - 2], hull[hull.len() - 1], p) <= 0 {
            hull.pop();
        }
        hull.push(p);
    }
    hull
}

fn main() {
    let pts = vec![(0, 0), (4, 0), (4, 4), (0, 4), (2, 2), (2, 0), (1, 3), (0, 0)];
    for p in graham(pts) {
        println!("({}, {})", p.0, p.1);
    }
}

Алгоритм Джарвиса («заворачивание подарка»)

Алгоритм Джарвиса строит ту же оболочку иначе: как будто обворачиваем точки верёвкой. Начинаем с самой левой точки $p_0$, она лежит на оболочке, и на каждом шаге находим следующую вершину.

Пусть $c$ — текущая вершина. Кандидат $n$ — такая точка, что все остальные точки лежат слева от направленного отрезка $c \to n$. Чтобы его найти, просматриваем все точки $p$: если $p$ лежит правее луча $c \to n$, она становится новым кандидатом. Если $p$ лежит на том же луче, берём более дальнюю точку, чтобы не класть в ответ точки, лежащие на стороне.

Каждый шаг занимает $O(n)$, а вершин в оболочке $h$, поэтому весь алгоритм работает за $O(n \cdot h)$. При малом $h$ он быстрее алгоритма Грэхема, а в худшем случае, когда все точки на оболочке, даёт $O(n^2)$.

Обход заканчивается, когда вернулись в $p_0$. Вычисления целочисленные, как и у алгоритма Грэхема: нужен только знак векторного произведения.

Реализация на Python

# Алгоритм Джарвиса («заворачивание подарка»), O(n·h)
def cross(o, a, b):
    """Векторное произведение (a - o) x (b - o): больше нуля, если поворот влево."""
    return (a[0] - o[0]) * (b[1] - o[1]) - (a[1] - o[1]) * (b[0] - o[0])


def dist2(a, b):
    return (a[0] - b[0]) ** 2 + (a[1] - b[1]) ** 2


def jarvis(points):
    """Выпуклая оболочка против часовой стрелки, начиная с самой левой точки."""
    pts = sorted(set(points))
    if len(pts) <= 2:
        return pts
    start = pts[0]
    hull = []
    current = start
    while True:
        hull.append(current)
        candidate = next(p for p in pts if p != current)
        for p in pts:
            turn = cross(current, candidate, p)
            # p правее луча current -> candidate, или на нём, но дальше: берём p
            if turn < 0 or (turn == 0 and dist2(current, p) > dist2(current, candidate)):
                candidate = p
        current = candidate
        if current == start:
            return hull


if __name__ == "__main__":
    print(jarvis([(0, 0), (4, 0), (4, 4), (0, 4), (2, 2), (2, 0), (1, 3), (0, 0)]))

Реализация на C++

// Алгоритм Джарвиса («заворачивание подарка»), O(n·h)
#include <algorithm>
#include <cstdio>
#include <vector>

struct Point {
    long long x;
    long long y;
};

bool operator==(const Point& a, const Point& b) { return a.x == b.x && a.y == b.y; }
bool operator<(const Point& a, const Point& b) { return a.x != b.x ? a.x < b.x : a.y < b.y; }

// Векторное произведение (a - o) x (b - o): больше нуля, если поворот влево
long long cross(const Point& o, const Point& a, const Point& b) {
    return (a.x - o.x) * (b.y - o.y) - (a.y - o.y) * (b.x - o.x);
}

long long dist2(const Point& a, const Point& b) {
    return (a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y);
}

// Выпуклая оболочка против часовой стрелки, начиная с самой левой точки
std::vector<Point> jarvis(std::vector<Point> pts) {
    std::sort(pts.begin(), pts.end());  // по (x, y)
    pts.erase(std::unique(pts.begin(), pts.end()), pts.end());
    if (pts.size() <= 2) return pts;

    const Point start = pts[0];
    std::vector<Point> hull;
    Point current = start;
    do {
        hull.push_back(current);
        Point candidate = pts[0] == current ? pts[1] : pts[0];
        for (const Point& p : pts) {
            long long turn = cross(current, candidate, p);
            // p правее луча current -> candidate или лежит на нём, но дальше: берём p
            if (turn < 0 || (turn == 0 && dist2(current, p) > dist2(current, candidate))) {
                candidate = p;
            }
        }
        current = candidate;
    } while (!(current == start));
    return hull;
}

int main() {
    const std::vector<Point> pts = {{0, 0}, {4, 0}, {4, 4}, {0, 4}, {2, 2}, {2, 0}, {1, 3}, {0, 0}};
    for (const Point& p : jarvis(pts)) {
        std::printf("(%lld, %lld)\n", p.x, p.y);
    }
    return 0;
}

Реализация на Rust

// Алгоритм Джарвиса («заворачивание подарка»), O(n·h)
type Point = (i64, i64);

// Векторное произведение (a - o) x (b - o): больше нуля, если поворот влево
fn cross(o: Point, a: Point, b: Point) -> i64 {
    (a.0 - o.0) * (b.1 - o.1) - (a.1 - o.1) * (b.0 - o.0)
}

fn dist2(a: Point, b: Point) -> i64 {
    (a.0 - b.0).pow(2) + (a.1 - b.1).pow(2)
}

// Выпуклая оболочка против часовой стрелки, начиная с самой левой точки
pub fn jarvis(mut pts: Vec<Point>) -> Vec<Point> {
    pts.sort();
    pts.dedup();
    if pts.len() <= 2 {
        return pts;
    }
    let start = pts[0];
    let mut hull = Vec::new();
    let mut current = start;
    loop {
        hull.push(current);
        let mut candidate = pts.iter().copied().find(|&p| p != current).unwrap();
        for &p in &pts {
            let turn = cross(current, candidate, p);
            // p правее луча current -> candidate или лежит на нём, но дальше: берём p
            if turn < 0 || (turn == 0 && dist2(current, p) > dist2(current, candidate)) {
                candidate = p;
            }
        }
        current = candidate;
        if current == start {
            return hull;
        }
    }
}

fn main() {
    let pts = vec![(0, 0), (4, 0), (4, 4), (0, 4), (2, 2), (2, 0), (1, 3), (0, 0)];
    for p in jarvis(pts) {
        println!("({}, {})", p.0, p.1);
    }
}

Принадлежность точки многоугольнику (правило луча)

Пусть простой многоугольник задан списком вершин в порядке обхода. Проведём из точки $p$ луч вправо, вдоль горизонтали, и посчитаем, сколько раз он пересекает рёбра. Нечётное число пересечений означает, что точка внутри, чётное — что снаружи.

Вершины, лежащие ровно на высоте $p_y$, нельзя считать дважды, поэтому используем полуинтервал: ребро $(a, b)$ учитываем, только если его концы лежат по разные стороны от горизонтали, то есть $(a_y > p_y) \neq (b_y > p_y)$. Горизонтальные рёбра при этом пропускаются автоматически.

Лежит ли точка пересечения правее $p$, проверяем без деления: это так тогда и только тогда, когда $(b-a)\times(p-a)$ и $b_y-a_y$ имеют один знак. Поэтому весь алгоритм целочисленный и работает за $O(n)$: один проход по рёбрам.

Точки на самой границе многоугольника этот алгоритм отдельно не выделяет, и результат для них не гарантирован. В демонстрации таких точек нет.

Реализация на Python

# Принадлежность точки многоугольнику (правило луча), O(n)
def cross(a, b, p):
    """Векторное произведение (b - a) x (p - a)."""
    return (b[0] - a[0]) * (p[1] - a[1]) - (b[1] - a[1]) * (p[0] - a[0])


def point_in_polygon(polygon, point):
    """Луч из точки вправо: нечётное число пересечений с рёбрами означает «внутри»."""
    inside = False
    for a, b in zip(polygon, polygon[1:] + polygon[:1]):
        # Ребро пересекает луч, только если его концы по разные стороны от горизонтали
        if (a[1] > point[1]) == (b[1] > point[1]):
            continue
        # Точка пересечения правее point, когда cross и dy одного знака
        if (cross(a, b, point) > 0) == (b[1] > a[1]):
            inside = not inside
    return inside


if __name__ == "__main__":
    poly = [(0, 0), (6, 0), (6, 6), (3, 2), (0, 6)]
    for point in [(3, 1), (3, 4), (1, 4), (5, 5)]:
        print(point, point_in_polygon(poly, point))

Реализация на C++

// Принадлежность точки многоугольнику (правило луча), O(n)
#include <cstdio>
#include <vector>

struct Point {
    long long x;
    long long y;
};

// Векторное произведение (b - a) x (p - a)
long long cross(const Point& a, const Point& b, const Point& p) {
    return (b.x - a.x) * (p.y - a.y) - (b.y - a.y) * (p.x - a.x);
}

// Луч из точки вправо: нечётное число пересечений с рёбрами означает «внутри»
bool point_in_polygon(const std::vector<Point>& polygon, const Point& point) {
    bool inside = false;
    const size_t n = polygon.size();
    for (size_t i = 0; i < n; ++i) {
        const Point& a = polygon[i];
        const Point& b = polygon[(i + 1) % n];
        // Ребро пересекает луч, только если его концы по разные стороны от горизонтали
        if ((a.y > point.y) == (b.y > point.y)) continue;
        // Точка пересечения правее point, когда cross и dy одного знака
        if ((cross(a, b, point) > 0) == (b.y > a.y)) inside = !inside;
    }
    return inside;
}

int main() {
    const std::vector<Point> poly = {{0, 0}, {6, 0}, {6, 6}, {3, 2}, {0, 6}};
    const std::vector<Point> tests = {{3, 1}, {3, 4}, {1, 4}, {5, 5}};
    for (const Point& p : tests) {
        std::printf("(%lld, %lld) %d\n", p.x, p.y, point_in_polygon(poly, p) ? 1 : 0);
    }
    return 0;
}

Реализация на Rust

// Принадлежность точки многоугольнику (правило луча), O(n)
type Point = (i64, i64);

// Векторное произведение (b - a) x (p - a)
fn cross(a: Point, b: Point, p: Point) -> i64 {
    (b.0 - a.0) * (p.1 - a.1) - (b.1 - a.1) * (p.0 - a.0)
}

// Луч из точки вправо: нечётное число пересечений с рёбрами означает «внутри»
pub fn point_in_polygon(polygon: &[Point], point: Point) -> bool {
    let mut inside = false;
    for (&a, &b) in polygon.iter().zip(polygon.iter().cycle().skip(1)) {
        // Ребро пересекает луч, только если его концы по разные стороны от горизонтали
        if (a.1 > point.1) == (b.1 > point.1) {
            continue;
        }
        // Точка пересечения правее point, когда cross и dy одного знака
        if (cross(a, b, point) > 0) == (b.1 > a.1) {
            inside = !inside;
        }
    }
    inside
}

fn main() {
    let poly = [(0, 0), (6, 0), (6, 6), (3, 2), (0, 6)];
    for point in [(3, 1), (3, 4), (1, 4), (5, 5)] {
        println!("{:?} {}", point, point_in_polygon(&poly, point));
    }
}