Выпуклые оболочки

Выпуклое множество — такое множество точек, что, для любых двух точек множества, все точки на отрезке между ними тоже принадлежат этому множеству.

Выпуклая оболочка множества точек — такое выпуклое множество точек, что все точки фигуры также лежат в нем.

Минимальная выпуклая оболочка множества точек — это минимальная по площади выпуклая оболочка.

Для экономии времени дальше минимальные выпуклые оболочки мы будем называть просто выпуклыми оболочками.

Для практических целей выпуклые оболочки полезны тем, что они компактно хранят всю необходимую информацию о множестве точек, что позволяет быстро отвечать на разнообразные запросы на этом множестве.

Выпуклые оболочки можно рассматривать в любом пространстве, но в этой статье мы ограничимся двумерным и научимся их эффективно строить по какому-то множеству из \(n\) точек на плоскости.

Алгоритм Джарвиса

Выберем какую-нибудь точку \(p_0\), которая гарантированно попадёт в выпуклую оболочку. Например, нижнюю, а если таких несколько, то самую левую из них.

Дальше будем действовать так: найдём самую «правую» точку от последней добавленной (то есть точку с минимальным полярным углом) и добавим её в оболочку. Будем так итеративно добавлять точки, пока не «замкнёмся», то есть пока самой правой точкой не станет \(p_0\).

Корректность алгоритма легко доказывается по индукции:

• На первом шаге мы выбрали точку, точно лежащую в МВО.

• На \(i\)-м шаге мы взяли такую точку, что все остальные лежат по «правильную» сторону отрезка \((p_i, p_{i+1})\).

Алгоритм Джарвиса также называют алгоритмом заворачивания подарка: мы каждый раз находим самый близкий «угол».

Для каждой точки выпуклой оболочки (обозначим их количество за \(h\)) мы из всех оставшихся \(O(n)\) точек будем искать оптимальную, что суммарно будет работать за \(O(n h)\).

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

vector<r> jarvis_convex_hull(vector<r> points) {
    r p0 = points[0];
    for (r p : points)
        if (p.x < p0.x || (p.x == p0.x && p.y < p0.y))
        p0 = p;
    vector<r> hull = {p0};
    while (true) {
        r t = p0; // кандидат на следующую точку
        for (r p : points)
            // лучше никакие полярные углы не считать
            if ((p - p0) ^ (t - p0) > 0)
                t = p;
        if (t == p0)
            continue;
        else {
            p0 = t;
            hull.push_back(t);
        }
    }
    return hull;
}

Важно помнить, что асимптотика именно \(O(nh)\), а не \(O(n^2)\): существуют задачи, где оболочка маленькая, и это существенно.

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

Алгоритм Грэхема — это оптимизация алгоритма Джарвиса, основанная на следующем наблюдении: если отсортировать все точки по полярному углу относительно точки \(p_0\), то выпуклая оболочка будет какой-то подпоследовательностью такого отсортированного массива точек.

Алгоритм последовательно строит выпуклые оболочки для каждого префикса этого отсортированного масива. При добавлении \(i\)-й точки в оболочку нужно удалить сколько-то последних добавленных точек, которые не будут входить в новую оболочку. Чтобы это делать эффективно, мы можем хранить выпуклую оболочку в стеке и в цикле while смотреть на три последние точки и проверять, образуют ли они правый поворот. Если это так, то среднюю следует удалить — мы нашли треугольник \((p_0, p_i, p_{i-2})\), который содержит \(p_{i-1}\), значит её можно удалить.

Каждая точка будет добавлена один раз удалена не более одного раза, что занимает константное количество операций. Соответственно, время работы будет упираться во время работы сортировки, то есть \(O(n \log n)\).

struct r {
    int x, y;
    r operator-(r a) {return {x - a.x; y - a.y};}
    // в векторном произведении возможно переполнение int
    int operator^(r a) {return x * a.y - y * a.x;}
};
vector<r> graham_convex_hull(vector<r> points) {
    // находим p0, как и раньше
    r p0 = points[0];
    for (r p : points)
        if (p.x < p0.x || (p.x == p0.x && p.y < p0.y))
            p0 = p;

    // TODO: удалить p0
    // сделаем p0 началом координат (для сокращения кода компаратора)
    for (r &x : points) {
        x.x -= p0.x;
        x.y -= p0.y;
    }
    // сортируем точки по полярному углу
    sort(points.begin(), points.end(), [&](r a, r b){
        return a ^ b > 0 || a ^ b == 0 && a.x * a.x + a.y * a.y < b.x * b.x + b.y * b.y;
    }); // если векторное произведение == 0 то сначала будем рассматривать точку, которая ближе к p0

    vector<r> hull;
    for (r p : points) {
        // удаляем последнюю точку МВО пока она образует невыпуклость (оператор "-" у точек должен быть перегружен, как и "^"
        /*
        Пусть т. A - последняя в оболочке, т. B - предпоследняя, т. C - точка, которую мы хотим добавить
        Тогда невыпуклость(или вырожденность) образуется, когда (C - A) ^ (B - A) <= 0
        */
        while (hull.size() >= 2 && ((p - hull.back()) ^ (hull[hull.size() - 2] - hull.back()) <= 0) {
             hull.pop_back();
        }
        hull.push_back(p);
    }
    // переместим начало координат обратно
    for (r &x : hull) {
        x.x += p0.x;
        x.y += p0.y;
    }
    return hull;
}

Верхние и нижние огибающие

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

Огибающие строить немного проще: можно отсортировать точки по \(x\) и пройтись по ним в таком порядке, поддерживая на стеке верхнюю огибающую для текущего префикса. При добавлении очередной точки нам нужно аналогичным образом проверить и удалить сколько-то верхних точек на стеке:

vector<r> upper_envelope(vector<r> points) {
    sort(points.begin(), points.end(), [](r a, r b){
        return a.x < b.x;
    });

    vector<r> s = {p0};
    for (r p : points) {
        while (...) {
            s[s.size()-2] = s[s.size()-1];
            s.pop_back();
        }
        s.push_back(p);
    }

    return s;
}

С огибающими работать легче, чем с целыми оболочками: их легко пересекать, объединять и делать бинпоиск (например, чтобы находить касательные). Кстати, объединение можно производить за линейное время (пройдясь двумя указателями по обеим оболочкам одновременно), что позволяет альтернативно строить огибающую методом «разделяй и властвуй», тоже за \(O(n \log n)\).

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

Алгоритм Эндрю

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

vector<r> andrew_convex_hull(vector<r> points) {
    vector<r> lower = // ...
    vector<r> upper = // ...
    vector<r> hull = // TODO     
}

Динамические выпуклые оболочки

Оперируя с оболочкой как с двумя огибающими, можно относительно просто обрабатывать запросы добавления точек, обернув огибающие, например, в std::set так, что они оказываются отсортированными по \(x\). При добавлении точки нужно просто найти её позицию в дереве и проверить и, возможно, удалить сколько-то её левых и правых соседей.

Обрабатывать удаление точки из множества сложнее. Если запросы известны заранее, то проще воспользоваться идеями корневой эвристики или dynamic connectivity problem и поддерживать только те точки, которые существуют на всем текущем блоке, и сводить удаление к добавлению.

Если же запросы удаления требуется обрабатывать онлайн (fully dynamic convex hull), то для этого есть очень неприятный алгоритм на 250-300 строк кода, заключающийся в поддержании огибающих для всех поддеревьев в этом дереве поиска и быстром объединении огибающих при перестроении дерева. Алгоритм работает за \(O(\log^2 n)\) на запрос: при слиянии огибающих при объединении вершин используется бинпоиск с разбором кучи случаев.

Алгоритм Чана

Как уже упоминалось ранее, иногда существенную роль играет то, что алгоритм Джарвиса работает за \(O(nh)\), а не \(O(n^2)\), то есть когда точек немного, он лучше алгоритмов Грэхэма и Эрдрю.

Алгоритм Чана пытается получить лучшее из двух и объединенить алгоритмы Джарвиса и Грэхэма (или Эндрю), чтобы получить асимптотику \(O(n \log h)\).

Алгоритм. Разделим все точки на группы по \(m\) точек. В каждой группе построим выпуклую оболочку за \(O(m \log m)\) алгоритмом Грэхэма. Точки никак не упорядочены, и эти оболочки могут пересекаться — это нормально. Суммарно для всех групп понадобится \(O(n \log m)\) операций.

Затем, начиная с самой левой нижей точки, мы будем строить общую выпуклую обочку аналогично алгоритму Джарвиса, но теперь «самую правую точку» можно находить каждый раз не за \(O(n)\), а за \(O(\frac{n}{m} \log m)\), если делать бинарный поиск в каждой из \(\frac{n}{m}\) оболочек.

Получается, что такое решение будет работать за \(O(h \frac{n}{m} \log m + n \log m)\). Если заранее приблизительно знать \(h\), то можно положить \(m = h\), и тогда асимптотика составит \(O(n \log h)\).

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

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

Пересечение полуплоскостей

Очень похожей задачей является пересечение полуплоскостей (англ. half-plane intersection). Требуется найти множество, удовлетворяющее набору неравенств:

\[ a_i x + b_i y + c \geq 0 \]

Например, через сведение к пересечению полуплоскостей можно (не самым оптимальным образом) строить подобные клёвые картинки:

Это диаграмма Вороного — для каждой точки выделено множество точек плоскости, которые являются к ней ближайшими. Чтобы построить это множество для данной точки \(p_i\), можно перебрать все осталньные точки \(p_j\), посчитать коэффициенты прямой, которая делит плоскость между ними пополам, и пересечь все ориентированные в нужную сторону полуплоскости, соответствующие этим прямым.

Результатом пересечения множества полуплоскостей тоже будет выпуклая оболочка — полуплоскость это выпуклое множество, а пересечение любых выпуклых множеств тоже выпуклое. Возможно, эта оболочка будет бесконечной по каким-то направлениям. Чтобы не обрабатывать это отдельно, часто добавляют bounding box — «коробку» из четырёх полуплоскостей на большом удалении.

Один из способов построения — построить отдельно верхнюю и нижнюю огибающую и пересечь. Для этого нужно разделить все полуплоскости на «смотрящие вверх» и «смотрящие вниз», для обеих групп отсортировать их по углу вектора нормали и пройтись по ним со стеком, удаляя «ненужные» полуплоскости — те, которые полностью покрываются новой и предпоследней полуплоскостью в стеке.

// TODO: я не хочу это кодить

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

Нахождение касательной

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

• Если это так, то она останется той же, и ничего дальше считать не надо.

• Если это не так, то новая точка должна лежать на пересечении новой полуплоскости и какой-то из старых. Важный факт здесь в том, что она лежит на границе новой полуплоскости. Чтобы её найти, пересечём все старые полуплоскости с этой границей и решим одномерную задачу пересечения множества интервалов, бесконечных в одну из сторон. Если получившийся интервал пуст, то и пересечение полуплоскостей пустое, а в противном случае один из его концов будет искомой точкой, если его положить обратно на границу новой полуплоскости.

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

Выясняется, что за счёт рандомизированности нам нужно выполнять это пересечение не так уж и часто. В случайном множестве из \(k\) полуплоскостей только две из них дают в пересечении оптимальную точку. Из этого следует, что вероятность того, что \(k\)-тая полуплоскость станет новой граничной (и понадобится всё пересчитывать) равна всего \(\frac{2}{k}\).

Таким образом, в среднем алгоритм совершит линейное количество операций:

\[ \frac{2}{2} \cdot 2 + \frac{2}{3} \cdot 3 + \frac{2}{4} \cdot 4 + \ldots + \frac{2}{n} \cdot n = \sum_k \frac{2}{k} \cdot k = O(n) \]

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