Хабрахабр

[Перевод] Алгоритм Форчуна, подробности реализации

Последние несколько недель я работал над реализацией алгоритма Форчуна на C++. Этот алгоритм берёт множество 2D-точек и строит из них диаграмму Вороного. Если вы не знаете, что такое диаграмма Вороного, то взгляните на рисунок:

Для каждой входной точки, которая называется «местом» (site), нам нужно найти множество точек, которые ближе к этому месту, чем ко всем остальным. Такие множества точек образуют ячейки, которые показаны на изображении выше.

В алгоритме Форчуна примечательно то, что он строит такие диаграммы за время $O(n\log n)$ (что оптимально для использующего сравнения алгоритма), где $n$ — это количество мест.

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

Как обычно, код выложен на github, а все использованные мной справочные материалы перечислены в конце статьи.

Я не буду объяснять, как работает алгоритм, потому что с этим уже хорошо справились другие люди. Могу порекомендовать изучить эти две статьи: здесь и здесь. Вторая очень любопытна — автор написал интерактивное демо на Javascript, которое полезно для понимания работы алгоритма. Если вам нужен более формальный подход со всеми доказательствами, то рекомендую прочитать главу 7 книги Computational Geometry, 3rd edition.

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

Я только написал псевдокод алгоритма, чтобы вы получили понятие о его глобальной структуре:

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

Первая проблема, с которой я столкнулся — выбор способа хранения диаграммы Вороного.

Я решил использовать структуру данных, широко применяемую в вычислительной геометрии, которая называется двусвязным списком рёбер (doubly connected edge list, DCEL).

Мой класс VoronoiDiagram использует в качестве полей четыре контейнера:

class VoronoiDiagram
{
public: // ... private: std::vector<Site> mSites; std::vector<Face> mFaces; std::list<Vertex> mVertices; std::list<HalfEdge> mHalfEdges;
}

Я подробно расскажу о каждом из них.

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

struct Site
{ std::size_t index; Vector2 point; Face* face;
};

Вершины ячейки представлены классом Vertex, у них есть только поле координат:

struct Vertex
{ Vector2 point;
};

Вот реализация полурёбер:

struct HalfEdge
{ Vertex* origin; Vertex* destination; HalfEdge* twin; Face* incidentFace; HalfEdge* prev; HalfEdge* next;
};

Вы можете задаться вопросом — что такое полуребро? Ребро в диаграмме Вороного является общим для двух соседних ячеек. В структуре данных DCEL мы разделяем эти рёбра на два полуребра, по одному для каждой ячейки, и они связываются указателем twin. Более того, у полуребра есть исходная и конечная точки. Поле incidentFace указывает на грань, которой принадлежит полуребро. Ячейки в DCEL реализуются как циклический двусвязный список полурёбер, в котором соседние полурёбра соединены вместе. Поэтому поля prev и next указывают на предыдующие и следующие полурёбра в ячейке.

На рисунке ниже показаны все эти поля для красного полуребра:

Наконец, класс Face задаёт ячейку. Он просто содержит указатель на своё место и ещё один на одно из его полурёбер. Не важно, какое из полурёбер выбрано, потому что ячейка — это замкнутый полигон. Таким образом, мы получаем доступ ко всем полурёбрам при обходе циклического связанного списка.

struct Face
{ Site* site; HalfEdge* outerComponent;
};

Стандартный способ реализации очереди событий — очередь с приоритетами. В процессе обработки событий места и окружностей нам может понадобиться удалить события окружностей из очереди, потому что они больше не являются действительными. Но большинство стандартных реализаций очереди с приоритетами не позволяют удалять элемент, не находящийся на вершине. В частности это относится и к std::priority_queue.

Первый, более простой — добавить к событиям флаг valid. Эту проблему можно решить двумя способами. Тогда вместо удаления события окружности из очереди, мы можем просто присвоить его флагу значение false. Изначально valid имеет значение true. Наконец, при обработке всех событий в основном цикле мы проверяем, равно ли значение флага valid события false, и если это так, то просто пропускаем его и обрабатываем следующее.

Вместо неё я реализовал собственную очередь с приоритетами, которая поддерживает удаление любого содержащегося в ней элемента. Второй способ, который применил я — не использовать std::priority_queue. Я выбрал этот способ, потому что он делает код алгоритма чище. Реализация такой очереди довольно проста.

Структура данных береговой линии — сложная часть алгоритма. В случае неверной реализации нет никаких гарантий, что алгоритм будет выполняться за $O(n\log n)$. Ключ к получению такой временной сложности — использование самобалансирующегося дерева. Но это проще сказать, чем сделать!

Но в них ничего не говорится о том, как сбалансировать дерево. В большинстве ресурсов, которые я изучал (две вышеупомянутые статьи и книга Computational Geometry), советуют реализовывать береговую линию как дерево, в котором внутренние узлы обозначают точки излома, а листья обозначают дуги. Думаю, что такая модель не самая лучшая, и вот почему:

  • в ней есть избыточная информация: мы знаем, что между двумя соседними дугами есть точка излома, поэтому необязательно представлять эти точки в виде узлов
  • она неадекватна для самобалансировки: можно сбалансировать только поддерево, образованное точками излома. Мы и в самом деле не можем сбалансировать всё дерево целиком, потому что в противном случае дуги могут стать внутренними узлами и листьями точек излома. Написание алгоритма для балансировки только поддерева, образованного внутренними узлами, мне кажется кошмарной задачей.

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

Вот определение дуги Arc в моей реализации:

struct Arc
; // Hierarchy Arc* parent; Arc* left; Arc* right; // Diagram VoronoiDiagram::Site* site; VoronoiDiagram::HalfEdge* leftHalfEdge; VoronoiDiagram::HalfEdge* rightHalfEdge; Event* event; // Optimizations Arc* prev; Arc* next; // Only for balancing Color color;
};

Первые три поля используются для структуры дерева. Поле leftHalfEdge указывает на полуребро, отрисованное крайней левой точкой дуги. А rightHalfEdge — на полуребро, отрисованное крайней правой точкой. Два указателя, prev и next используются для получения прямого доступа к предыдущей и следующей дуге береговой линии. Кроме того, они также позволяют обойти береговую линию как двусвязанный список. И наконец, у каждой дуги есть цвет, который используется для балансировки береговой линии.

При написании кода я вдохновлялся книгой Introduction to Algorithms. Для балансировки береговой линии я решил воспользоваться красно-чёрной схемой. В главе 13 описаны два интересных алгоритма, insertFixup и deleteFixup, которые балансируют дерево после вставки или удаления.

В алгоритме Форчуна нет ключей, мы только знаем, что нужно вставлять дугу до или после другой в береговой линии. Однако я не могу использовать приведённый в книге метод insert, потому что для нахождения места вставки узла в нём используются ключи. Чтобы реализовать это, я создал методы insertBefore и insertAfter:

void Beachline::insertBefore(Arc* x, Arc* y)
{ // Find the right place if (isNil(x->left)) { x->left = y; y->parent = x; } else { x->prev->right = y; y->parent = x->prev; } // Set the pointers y->prev = x->prev; if (!isNil(y->prev)) y->prev->next = y; y->next = x; x->prev = y; // Balance the tree insertFixup(y); }

Вставка y перед x выполняется в три этапа:

  1. Находим место для вставки нового узла. Для этого я воспользовался следующим наблюдением: левый дочерний элемент x или правый дочерний элемент x->prev равен Nil, и тот из них, который равен Nil находится перед x и после x->prev.
  2. Внутри береговой линии мы сохраняем структуру двусвязного списка, поэтому мы должны соответствующим образом обновлять указатели prev и next элементов x->prev, y и x.
  3. Наконец, мы просто вызываем для балансировки дерева метод insertFixup, описанный в книге.

insertAfter реализуется аналогично.

Взятый из книги метод удаления можно реализовать без изменений.

Вот выходные данные описанного выше алгоритма Форчуна:

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

Например, если мы возьмём три точки, находящиеся на одной прямой, то у средней точки будет два бесконечных полуребра, не связанных вместе. Хуже того — ячейка может и не быть одним фрагментом. Это не очень нас устраивает, потому что мы не сможем получить доступ к одному из полурёбер, ведь ячейка является связанным списком рёбер.

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

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

Мой алгоритм ограничения получает в качестве входных данных прямоугольник (box) и состоит из трёх этапов:

  1. Он обеспечивает размещение каждой вершины диаграммы внутри прямоугольника.
  2. Отсекает каждое бесконечное ребро.
  3. Замыкает ячейки.

Этап 1 тривиален — нам просто нужно расширить прямоугольник, если он не содержит вершину.

Этап 2 тоже довольно прост — он состоит из вычисления пересечений между лучами и прямоугольником.

Я выполняю его в два этапа. Этап 3 также не очень сложен, только требует внимательности. Затем я делаю так, чтобы все вершины ячейки были соединены полурёбрами. Сначала я добавляю точки углов прямоугольника к ячейкам, в вершинах которых они должны быть.

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

Вот выходная диаграмма ограничивающего алгоритма:

Теперь мы видим, что все рёбра отрисовываются. И если уменьшить масштаб, то можно убедиться, что все ячейки замкнуты:

Отлично! Но первое изображение из начала статьи лучше, правда?

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

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

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

  1. Полуребро полностью находится внутри прямоугольника: такое полуребро мы сохраняем
  2. Полуребро полностью находится снаружи прямоугольника: такое полуребро мы отбрасываем
  3. Полуребро выходит наружу прямоугольника: мы усекаем полуребро и сохраняем его как последнее полуребро, выходившее наружу.
  4. Полуребро входит внутрь прямоугольника: мы усекаем полуребро, чтобы связать его с последним полуребром, выходившим наружу (мы сохраняем его в случае 3 или 5)
  5. Полуребро пересекает прямоугольник дважды: мы усекаем полуребро и добавляем полурёбра, чтобы связать ег с последним полуребром, выходившим наружу, а затем сохраняем его как новое последнее полуребро, выходившее наружу.

Да, получилось много случаев. Я создал картинку, чтобы все их показать:

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

Применив этот алгоритм к ограниченной диаграмме, мы получаем ожидаемый результат:

Статья оказалась довольно длинной. И я уверен, что многие аспекты до сих пор вам непонятны. Тем не менее, надеюсь, она будет вам полезной. Изучите код, чтобы разобраться в подробностях.

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

  • $n = 1000$: 2 мс
  • $n = 10000$: 33 мс
  • $n = 100000$: 450 мс
  • $n = 1000000$: 6600 мс

Мне не с чем сравнить эти показатели, но кажется, что это невероятно быстро!

Теги
Показать больше

Похожие статьи

Добавить комментарий

Ваш e-mail не будет опубликован. Обязательные поля помечены *

Кнопка «Наверх»
Закрыть