Геометрия — 2

В прошлой лекции мы научились работать с точками, векторами, прямыми, отрезками и лучами. В этой лекции мы перейдём к более сложным объектам — многоугольникам. Мы научимся вычислять их площади, проверять принадлежность точки многоугольнику и познакомимся с формулой Пика.

Напомним основные инструменты из прошлой лекции, которыми мы будем пользоваться: скалярное произведение dot, векторное произведение cross, длина вектора length и расстояние между точками.

Note

В этой лекции мы изменим тип координат в структуре Point с int на double. Это связано с тем, что некоторые алгоритмы (например, метод суммы углов или метод лучей со случайным направлением) существенно используют вещественную арифметику. При сравнениях теперь придётся учитывать погрешность и использовать std::abs(x) < eps вместо x == 0.

struct Point {
    double x = 0.0;
    double y = 0.0;

    Point() = default;
    Point(double x, double y) : x(x), y(y) {}
};

Многоугольники

Многоугольник — часть плоскости, ограниченная замкнутой ломаной.

Примеры многоугольников

  • Многоугольник называется простым, если граничная ломаная не имеет самопересечений. На картинке все многоугольники, кроме последнего, простые.

  • Многоугольник называется выпуклым, если для любых его двух точек весь отрезок между ними целиком лежит внутри многоугольника. На картинке первые два многоугольника выпуклые.

  • Выпуклый многоугольник называется правильным, если у него все стороны и все углы равны. Все правильные многоугольники можно вписать в окружность; обратное неверно.

В трёхмерном пространстве обобщением многоугольника является многогранник.

В программировании многоугольник естественно хранить как последовательность вершин в порядке обхода:

std::vector<Point> polygon;

Площади

Площадь треугольника

Из прошлой лекции мы знаем, что модуль векторного произведения двух векторов равен площади параллелограмма, натянутого на них. Площадь треугольника — половина этого параллелограмма:

\[ S_{\triangle ABC} = \frac{1}{2} \cdot |(B - A) \times (C - A)| \]

Без модуля получается ориентированная площадь: положительная, если обход \(A \to B \to C\) идёт против часовой стрелки, и отрицательная — если по часовой стрелке.

double triangleArea(Point a, Point b, Point c) {
    return std::abs(cross(b - a, c - a)) / 2.0;
}

Площадь произвольного многоугольника

Для произвольного многоугольника, заданного последовательностью вершин \(A_1, A_2, \dots, A_n\) в порядке обхода, используется формула площади Гаусса (также называется формулой шнурования, англ. shoelace formula):

\[ S = \frac{1}{2} \left| \sum_{i=1}^{n} A_i \times A_{i+1} \right| \]

где \(A_{n+1} = A_1\) (обход замкнутый).

Идея вывода. Зафиксируем начало координат \(O = (0, 0)\) и рассмотрим серию треугольников \(OA_iA_{i+1}\) для каждого ребра многоугольника. Площадь такого треугольника с учётом знака равна \(\frac{1}{2} A_i \times A_{i+1}\) (радиус-векторы \(A_i\) и \(A_{i+1}\) играют роль двух сторон треугольника из начала координат).

Сумма ориентированных треугольников даёт площадь многоугольника

Когда мы суммируем эти ориентированные площади, области вне многоугольника учитываются с противоположными знаками и сокращаются — остаётся именно площадь многоугольника. Формула работает для любого простого многоугольника, в том числе невыпуклого.

Note

Знак результата зависит от направления обхода (по часовой или против). Поэтому в формуле мы берём модуль. Сам знак, впрочем, тоже полезен — по нему можно определить ориентацию многоугольника.

Как реализовать. Функция должна пройти по всем парам соседних вершин (включая пару последняя–первая), накапливая сумму векторных произведений:

double polygonArea(const std::vector<Point>& p) {
    int n = (int) p.size();
    double s = 0.0;
    for (int i = 0; i < n; ++i) {
        Point a = p[i];
        Point b = p[/* индекс следующей вершины циклически */];
        s += /* векторное произведение радиус-векторов a и b */;
    }
    return std::abs(s) / 2.0;
}

Приём для циклического перехода: (i + 1) % n — когда i == n - 1, результат равен \(0\), и мы возвращаемся к первой вершине. Этот приём пригодится практически во всех дальнейших алгоритмах.

Метод трапеций как альтернатива

Есть и более интуитивный (но более громоздкий) способ — метод трапеций. Мы так же проходимся по всем рёбрам, но теперь складываем ориентированные площади не треугольников с началом координат, а трапеций с основанием на оси \(Ox\):

Разложение площади на трапеции

«Нижние» трапеции при обходе по часовой стрелке имеют отрицательный знак и сокращают лишнюю площадь «верхних». Площадь одной трапеции на ребре \((A_i, A_{i+1})\):

\[ S_{\text{трап}} = \frac{(y_i + y_{i+1})}{2} \cdot (x_{i+1} - x_i) \]

Формула Гаусса и метод трапеций дают один и тот же ответ.

Целочисленные координаты

Из формулы для площади треугольника следует важное наблюдение: если все координаты вершин целые, то площадь любой фигуры будет либо целым, либо полуцелым числом (т.е. числом с \(2\) в знаменателе).

Если в задаче все входные координаты целые, часто имеет смысл заранее удвоить их (или работать с удвоенной площадью) — тогда все промежуточные значения останутся целыми, и можно не беспокоиться о погрешностях плавающей точки.

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

Точка в треугольнике

Начнём с простого случая — проверить, находится ли точка \(P\) внутри треугольника \(ABC\), заданного против часовой стрелки.

Идея: точка лежит внутри треугольника тогда и только тогда, когда она находится «слева» от всех трёх рёбер. «Слева» — потому что обход против часовой стрелки. Признак «слева» выражается знаком векторного произведения:

\[ P \in \triangle ABC \iff \begin{cases} (B - A) \times (P - A) \geq 0 \\ (C - B) \times (P - B) \geq 0 \\ (A - C) \times (P - C) \geq 0 \end{cases} \]

Если все три произведения строго положительные — точка строго внутри. Если какое-то из них равно нулю — точка лежит на соответствующем ребре.

Warning

Поскольку координаты теперь double, сравнения с нулём стоит делать с погрешностью: x > -eps вместо x >= 0. Типичное значение eps для таких задач — \(10^{-9}\).

bool pointInTriangle(Point p, Point a, Point b, Point c) {
    const double eps = 1e-9;
    double c1 = cross(b - a, p - a);    // знак относительно ребра AB
    double c2 = /* аналогично для ребра BC */;
    double c3 = /* аналогично для ребра CA */;
    return c1 > -eps && c2 > -eps && c3 > -eps;
}

Два оставшихся выражения пишутся по тому же образцу, что и первое. Обратите внимание: если обход треугольника задан по часовой стрелке, все знаки будут обратными — тогда нужно проверять c1 < eps && c2 < eps && c3 < eps.

Точка в выпуклом многоугольнике

Идея обобщается на произвольный выпуклый многоугольник. Если вершины обходятся против часовой стрелки, точка внутри тогда и только тогда, когда она «слева» от всех рёбер:

bool pointInConvex(Point p, const std::vector<Point>& poly) {
    const double eps = 1e-9;
    int n = (int) poly.size();
    for (int i = 0; i < n; ++i) {
        Point a = poly[i];
        Point b = poly[(i + 1) % n];
        double v = /* векторное произведение векторов AB и AP */;
        if (v < -eps) {
            return /* точка справа от ребра — ... */;
        }
    }
    return /* прошли все рёбра — ... */;
}

Проверка на выпуклость

Чтобы проверить, является ли многоугольник выпуклым, достаточно пройтись по всем тройкам соседних вершин \((a, b, c)\) и убедиться, что мы всегда «поворачиваем» в одну и ту же сторону. Это значит, что знак векторного произведения \((b - a) \times (c - b)\) не меняется вдоль всего обхода:

bool isConvex(const std::vector<Point>& poly) {
    int n = (int) poly.size();
    int sign = 0;  // знак первого ненулевого векторного произведения
    for (int i = 0; i < n; ++i) {
        Point a = poly[i];
        Point b = poly[(i + 1) % n];
        Point c = poly[/* индекс вершины через одну циклически */];
        double v = cross(b - a, c - b);
        if (v > 0) {
            if (sign < 0) return false;  // встретили разные знаки
            sign = 1;
        } else if (v < 0) {
            if (/* ... */) return false;
            sign = -1;
        }
        // случай v == 0 пропускаем: три точки на одной прямой
    }
    return true;
}

Основные места для заполнения: индекс третьей вершины (i + 2) % n и симметричная проверка смены знака в ветке v < 0.

Точка в произвольном многоугольнике

Для невыпуклых многоугольников проверка «слева от всех рёбер» не работает. Есть два стандартных подхода, оба за \(O(n)\) для одного запроса.

Метод суммы углов

Пройдёмся по всем вершинам многоугольника в порядке обхода и для каждой пары соседних вершин \(A_i\), \(A_{i+1}\) вычислим ориентированный угол, под которым отрезок \(A_i A_{i+1}\) виден из точки \(P\). Сумма этих углов по всем рёбрам:

\[ \theta = \sum_{i=1}^{n} \angle(A_i - P,\ A_{i+1} - P) \]

обладает замечательным свойством:

  • \(|\theta| \approx 2\pi\) — точка внутри многоугольника
  • \(\theta \approx 0\) — точка снаружи

Почему это работает: когда \(P\) внутри, луч \(P \to A_i\) делает полный оборот при обходе многоугольника, поэтому сумма углов равна \(\pm 2\pi\). Когда \(P\) снаружи, направление «отклоняется» и возвращается назад, итог — ноль.

Ориентированный угол между двумя векторами мы умеем считать из прошлой лекции:

double angleBetween(Point a, Point b) {
    return std::atan2(cross(a, b), dot(a, b));
}
bool pointInPolygonAngle(Point p, const std::vector<Point>& poly) {
    int n = (int) poly.size();
    double total = 0.0;
    for (int i = 0; i < n; ++i) {
        Point a = poly[i];
        Point b = poly[(i + 1) % n];
        total += /* угол, под которым отрезок ab виден из p, через angleBetween */;
    }
    return std::abs(total) > /* порог между 0 и 2*pi */;
}

В роли порога сравнения подойдёт любое значение строго между \(0\) и \(2\pi\) — например, \(\pi\): так мы отличаем сумму углов, близкую к \(\pm 2\pi\) (точка внутри), от суммы, близкой к \(0\) (снаружи).

Метод лучей

Второй метод: проведём из точки \(P\) произвольный луч и посчитаем, сколько раз он пересекает границу многоугольника.

Ключевое наблюдение: если точка внутри многоугольника, луч выйдет наружу нечётное число раз («внутри → снаружи → внутри → снаружи»); если снаружи — чётное число раз.

Метод лучей для проверки принадлежности точки
Tip

Чтобы избежать вырожденных случаев (луч проходит ровно через вершину или вдоль ребра), удобно взять луч в случайном направлении. Вероятность, что он случайно совпадёт с какой-то особой точкой, пренебрежимо мала.

bool pointInPolygonRay(Point p, const std::vector<Point>& poly) {
    Point d = {1.7, 2.3};  // "некрасивое" направление луча, чтобы избежать вырожденных случаев
    int n = (int) poly.size();
    int crossings = 0;
    for (int i = 0; i < n; ++i) {
        Point a = poly[i];
        Point b = poly[(i + 1) % n];
        if (/* луч p + t*d (t >= 0) пересекает отрезок [a, b] */) {
            ++crossings;
        }
    }
    return crossings % 2 == /* остаток, соответствующий "точка внутри" */;
}

Проверку пересечения луча с отрезком удобно вынести в отдельную функцию. Её можно свести к уже знакомой задаче пересечения двух прямых из прошлой лекции — надо только дополнительно проверить, что параметр \(t\) на луче неотрицательный, а параметр \(s\) на отрезке лежит в \([0, 1]\):

bool raySegmentIntersect(Point p, Point d, Point a, Point b) {
    // p + t * d = a + s * (b - a),  где t >= 0 и s in [0, 1]
    // это система из двух уравнений (по одному на каждую координату)
    // относительно двух неизвестных t и s — решите её, как в прошлой лекции
    // вернуть true, если t >= -eps и 0 - eps <= s <= 1 + eps
}

Формула Пика

Для многоугольника, все вершины которого имеют целые координаты, существует изящная формула, связывающая его площадь с числом целочисленных точек:

\[ S = V + \frac{E}{2} - 1 \]

где

  • \(V\) — количество целочисленных точек внутри многоугольника,
  • \(E\) — количество целочисленных точек на границе многоугольника (включая вершины).

Пример. Рассмотрим квадрат с вершинами \((0, 0)\), \((3, 0)\), \((3, 3)\), \((0, 3)\).

  • Площадь: \(S = 9\)
  • Внутренние целочисленные точки: \((1, 1), (1, 2), (2, 1), (2, 2)\), итого \(V = 4\)
  • Точки на границе: по \(4\) точки на каждой стороне (включая один конец), итого \(E = 12\)

Проверяем формулу Пика: \(V + \frac{E}{2} - 1 = 4 + 6 - 1 = 9\) — совпадает с прямым расчётом площади.

Формула полезна в задачах, где просят найти количество целочисленных точек внутри многоугольника: если мы можем вычислить \(S\) и \(E\) (что несложно — количество точек на отрезке \((x_1, y_1)\)—\((x_2, y_2)\) равно \(\gcd(|x_2 - x_1|, |y_2 - y_1|) + 1\)), то \(V\) находится из формулы Пика.

Выпуклая оболочка (опционально)

Выпуклой оболочкой множества точек называется наименьший выпуклый многоугольник, содержащий все эти точки. Представьте, что точки — это гвозди на доске; выпуклая оболочка — это то, как натянется резинка, наброшенная на них.

Выпуклая оболочка — резинка, натянутая на точки

Выпуклая оболочка — один из центральных объектов вычислительной геометрии. Существуют алгоритмы её построения за \(O(n \log n)\): обход Грэхема (Graham scan), алгоритм Эндрю (Andrew’s monotone chain), обход Джарвиса (Jarvis march, за \(O(nh)\), где \(h\) — число вершин оболочки).

Типичные применения:

  • Проверка принадлежности точки набору точек (достаточно проверить, что точка лежит внутри выпуклой оболочки — за \(O(\log n)\) с предпосчитанной оболочкой).

  • Нахождение диаметра множества точек (самой дальней пары).

  • Задачи о пересечении полуплоскостей.

Подробное изложение алгоритмов выходит за рамки нашего курса, но рекомендуем ознакомиться с ними самостоятельно — это интересная и полезная тема.