Алгоритм Quickhull для нахождения выпуклой оболочки

Как гласит определение, выпуклая оболочка некоторого множества — это наименьшее выпуклое множество , содержащее в себе множество . Выпуклой оболочкой конечного множества попарно различных точек является многогранник.
Для реализации одномерного случая алгоритма Quickhull годится функция std::minmax_element. В сети можно найти множество реализаций алгоритма Quickhull для плоского случая. Однако, для случая произвольной размерности сходу находится лишь одна тяжёловесная реализация с сайта qhull.org.

В скобках рядом с терминами я буду давать перевод на английский (в стиле Туґбо, уж извините). Это может быть полезно для тех, кто ещё не знаком с терминами и будет читать англоязычные тексты по ссылкам, либо решит разобраться с текстами исходного кода.
Упомянутая выше «каноническая» реализация написана на C (есть C++ биндинги), используется давно и широко (в составе пакетов GNU Octave и MATLAB) и, следовательно, хорошо проверена. Но если вы решите изолировать код, относящийся только лишь к одному алгоритму Quickhull, то столкнётесь со сложностями. Для простейших применений хотелось бы иметь возможность обходится чем-то более компактным. Немного постаравшись, я написал свою реализацию (https://github.com/tomilov/quickhull): это единственный класс, без зависимостей. Всего порядка 750 строк.

Геометрические понятия.

Начиная работать с многомерными сущностями, понимаешь, что справедливы какие-то аналогии с привычными двухмерными и трёхмерными объектами. Необходимо несколько упорядочить знания: ввести некоторые новые термины и прояснить отношения между ними.
Для программиста геометрические объекты — это прежде всего структуры данных. Для начала я дам геометрические определения (не очень строгие). Дальше по тексту будут приведены описания соответственных структур данных, удобных для того, чтобы хранить данные и оперировать с ними достаточно эффективно.
Обозначим размерность пространства как .

• Симплекс (англ. simplex) задаётся аффинно независимой (англ. affinely independent) точкой. Об аффинной независимости поясню далее. Эти точки называются вершинами (англ. ед. vertex, мн. vertices).
• Многогранник (син. политоп, полиэдр, англ. polytope, polyhedron) определяется как минимум точкой (вершиной), из которых должны являться аффинно независимыми; симплекс (англ. simplicial polytope) — простейший случай -мерного тела: многогранники с меньшим числом вершин обязательно будут иметь нулевой -мерный объём.
• Параллелотоп (англ. parallelotope) — обобщение плоского параллелограмма и объёмного параллелепипеда. Для симплекса можно построить соответственных параллелотопов (все они будут иметь одинаковый объём, равный объёмам симплекса), взяв в качестве образующих (векторов) параллелотопа вектора, выходящие из одной фиксированной вершины в каждую из остальных вершин.
• Понятие выпуклого многогранника (англ. convex polytope, convex polyhedron) я здесь приводить не буду, сгодится ваше интуитивное представление о нём. Симплекс, рассматриваемый как частный случай многогранника, всегда является выпуклым.
• Понятие плоскости является исходным и неопределяемым. Замечу, что она имеет на единицу меньшую размерность, чем пространство, в которое вложена.
• Симплициальная грань (англ. simplicial facet) определяется аффинно независимыми точками (вершинами). Далее речь будет идти только о симплициальных объектах (кроме выпуклого многогранника), так что определение «симплициальный» я буду опускать. Многие утверждения в дальнейшем будут справедливы только для невырожденных (все точки попарно различны) случаев симплициальных геометрических объектов.
• Ребро (англ. ridge) определяется как пересечение двух граней и имеет вершину. Две грани могут иметь только одно общее ребро. Одна грань, таким образом, может иметь соседей через рёбра.
• Пик (англ. peak) определяется точками. Здесь интуитивная аналогия с трёхмерным пространством может поплыть, так как понятие пика и вершины не совпадают, но такими объектами мы не будем оперировать. Грань может иметь любое количество соседних граней через пик, отсюда можно сделать вывод, что хранить граф смежности граней через пики накладно и по памяти и по затратам на поддержку актуальности структур данных при каких-либо преобразованиях.

Грань многогранника, её границы, границы её границ и т.д. на английском именуются face. Дам ссылку на хороший ресурс, где логически обосновывается, что в -мерном пространстве видимыми или наблюдаемыми объектами могут являться только лишь -мерные объекты, то есть границы (англ. boundaries) -мерных объектов. Английский термин facet — это именно наблюдаемая грань. Прямая (англ. straight line) — это всегда одномерная сущность. Точку (вершину) можно в этой градации считать нуль-мерной.
Имея дело с некоторым набором из точек, мы можем перейти к рассмотрению соответственного набора из векторов, вычтя одну из точек из всех других. Таким образом мы как бы назначим эту точку началом координат (англ. origin). Если все точки лежали в одной плоскости, то плоскость смещается таким образом, что начинает проходить через начало координат. Аффинная независимость точек выполняется, когда выполняется линейная независимость -го соответственного вектора. Определение линейной независимости векторов вам скорей всего уже знакомо. Таким образом в -мерном пространстве нельзя выбрать больше точки так, чтобы все они были аффинно независимы. Поясню это: в трёхмерном пространстве треугольник (грань трёхмерного симплекса — тетраэдра) задан тремя точками (не лежащими на одной прямой, конечно); через эти три точки проходит единственная плоскость; любая точка этой плоскости, добавленная ко множеству из трёх вершин треугольника превратит можнество трёх аффинно независимых точек в множество 4 аффинно зависимых. Аналогия справедлива для других размерностей.
Кроме всего прочего необходимо ввести понятие гиперобъёма (англ. hypervolume).

Алгоритм.

Сам алгоритм Quickhull для случая произвольной размерности был предложен в труде Barber, C. Bradford; Dobkin, David P.; Huhdanpaa, Hannu (1 December 1996). «The quickhull algorithm for convex hulls». ACM Transactions on Mathematical Software 22 (4): 469–483. «Каноническая» реализация на C от авторов располагается на уже упомянутом сайте http://qhull.org/, репозиторий с C++ интерфейсом https://gitorious.org/qhull/qhull. Процитирую алгоритм из первоисточника:

create a simplex of d + 1 points
for each facet F
	for each unassigned point p
		if p is above F
			assign p to F`s outside set
for each facet F with a non-empty outside set
	select the furthest point p of F`s outside set
	initialize the visible set V to F
	for all unvisited neighbours N of facets in V
		if p is above N
			add N to V
	the set of horizon ridges H is the boundary of V
	for each ridge R in H
		create a new facet from R and P
		link the new facet to its neighbours
	for each new facet F'
		for each unassigned point q in an outside set of a facet in V
			if q is above F'
				assign q to F'`s outside set
	delete the facets in V

Выпуклая оболочка представляет собой список граней.

Стартовый симплекс.

Как видно, мы должны начать с некоторого стартового симплекса. Для этого можно выбрать любые аффинно независимые точки, но лучше использовать какую-нибудь дешёвую эвристику. Основное требование на этом шаге — чтобы наш симплекс по возможности не получился «узким» (англ. narrow). Это важно в случае, если используется арифметика с плавающей точкой (англ. floating-point arithmetic). В оригинальной реализации используются точки с максимальными и минимальными значениями координат (очевидно, они же войдут во множество вершин выпуклой оболочки). Начав со стартового симплекса, на каждом последующем шаге мы будем получать всё больший и больший многогранник, неизменным свойством которого всегда будет его выпуклость. Я буду называть его временный многогранник.


Далее мы должны распределить оставшиеся точки в так называемые списки внешних точек (англ. outside sets) граней — это такие списки точек, которые ведутся для каждой грани и в которые попадают ещё не рассмотренные (далее — неназначенные, англ. unassigned) точки, находящиеся относительно плоскости грани по другую сторону, нежели внутренность симплекса (или временного многогранника в дальнейшем). Здесь мы пользуемся тем свойством выпуклого многогранника (и симплекса в частности), что весь он лежит целиком только в одном из двух полупространств, получаемых при разбиении всего пространства плоскостью любой грани. Из этого свойства следует тот факт, что если точка не попала ни в один из списков внешних точек граней выпуклой оболочки, то она достоверно находится внутри этой выпуклой оболочки. В дальнейшем в алгоритме понадобится информация о том, какая из точек, содержащихся в списке внешних точек, наиболее удалена от плоскости грани, поэтому, добавляя точки в список внешних точек, необходимо сохранять информацию о том, которая из них самая дальняя. Любая точка попадает только в один список внешних точек. На асимптотическую сложность не влияет то как именно распределяются точки по спискам.
На данном этапе мне необходимо сделать отступление и рассказать о том, как задать грань, как вычислять расстояние до неё, как определить ориентацию точки относительно грани.


Теперь описывать плоскость грани мы умеем, но осталась одна неясность: как же выбрать порядок перечисления точек в списках вершин граней стартового симплекса? Авторы не хранят список вершин упорядоченным. Они поступают иначе. Задав уравнение плоскости для произвольного, но фиксированного порядка вершин, они выясняют ориентацию какой-либо внутренней точки (англ. inner point) относительно этой плоскости. То есть вычисляют ориентированное расстояние. В случае, если расстояние отрицательное они меняют знак у координат нормали к гиперплоскости и у её смещения от начала координат. Либо задают значение флага, сигнализирующее о том, что грань перевёрнута (англ. flipped). Внутренней точкой для выпуклого многогранника будет среднее арифметическое хотя бы его вершин (т.н. центр масс, если взять все вершины).


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

Главный цикл.

Теперь есть временный многоугольник с которым можно работать в главном цикле. В главном цикле этот временный многоугольник достраивается, «захватывая» всё больше и больше пространства, а вместе с ним и точки исходного множества, пока не настанет такой момент, когда вне многогранника не останется ни одной точки. Далее опишу этот процесс более подробно. А именно то, как он устроен в моей реализации.
Главный цикл — это цикл по граням. На каждой итерации среди граней с непустыми списками внешних точек выбирается самая лучшая грань (англ. best facet); лучшая в том смысле, что самая дальняя точка из списка внешних точек находится дальше от грани, чем у других граней. Цикл завершается, когда граней с непустыми списками внешних точек не осталось.
В теле главного цикла выбранная грань удаляется как заведомо не являющаяся гранью выпуклой оболочки. Но предварительно она «разбирается на составляющие». Из списка внешних точек извлекается самая дальняя точка (англ. furthest point) — она наверняка является вершиной выпуклой оболочки. Остальные точки перемещаются в отдельный временный список неназначенных точек . Далее по спискам соседних граней, начиная со списка грани , обходится граф смежности граней временного многогранника и каждая грань добавляется в список уже посещённых (англ. visited) граней и тестируется на видимость (англ. visibility) из точки . Если грань невидимая из точки , то её соседние грани далее не обходятся, иначе обходятся. Если у очередной видимой грани среди соседних граней есть невидимые из точки , то она попадает в список граничных видимых граней (англ. before-the-horizon как противоположность over-the-horizon — «загоризонтный»). Но если все соседние грани являются видимыми, то такая грань добавляется в список граней на удаление. После обхода графа смежности граней список граничных видимых граней и список граней на удаление в сумме содержат все видимые из точки грани. Теперь удаляются все грани из списка на удаление. Причём каждая удаляемая грань не удаляется из списков соседних граней у своих соседей. Так делать безопасно, потому что даже в случае если соседняя грань является граничной видимой гранью, из списка её соседних граней в дальнейшем используется только информация о соседстве с невидимыми из точки гранями. Следом за списком граней на удаление просматривается список граничных видимых граней. Каждая граничная видимая грань «разбирается» (далее понадобится её список соседних граней и список вершин) и удаляется. Для каждой граничной видимой грани создаются новые грани (от одной до штук), которые добавляются в список новых граней. При этом просматривается список её соседних граней и для каждой грани из этого списка, являющейся невидимой, ищется общее ребро с данной граничной видимой гранью. В список вершин новой грани выписываются вершины общие для невидимой и граничной видимой граней (т.е. вершины общего ребра), причём именно в том порядке, в котором они записаны в списке вершин граничной видимой грани. Ребро (англ. horizon ridge), кстати, является частью горизонта той части поверхности временного многогранника, которая видна из точки . Та единственная вершина граничной видимой грани, которая не пренадлежит этому ребру, заменяется на точку в списке вершин новой грани, остальные же вершины — те же что и в (упорядоченном) списке вершин граничной видимой грани и стоят на тех же местах. Таким образом построенный список вершин новой грани гарантирует её правильную ориентацию, то есть при расчёте её нормированного уравнения плоскости нормаль будет направлена вовне временного многогранника и значение смещения получит верный знак. После того, как все граничные видимые грани обработаны мы имеем список новых граней. Уравнения гиперплоскостей для этих граней возможно рассчитать сразу после того, как мы имеем списки вершин для них. Так же, при создании новых граней, в их списки сосдених граней можно добавить по одной грани, а именно: соответствующую невидимую грань, которая была соседней для граничной видимой грани при формировании списка вершин новой грани. Очевидно, что остальные соседние грани необходимо искать среди остальных новых граней.
Как оказалось (при семплировании профайлером из Google Performance Tools), больше всего времени в случае больших размерностей и/или плохих входных данных (т.е. когда большая часть точек является вершинами выпуклой оболочки) алгоритм проводит в поиске соседних граней для новых граней. В моей реализации изначально алгоритм был устроен предельно просто: перебирались все возможные пары (два вложенных цикла) граней из списка новых граней и производились сравнения упорядоченных (по номеру из входного списка точек) списков вершин с одним несовпадением (модифицированный алгоритм std::set_difference или std::set_symmetric_difference). То есть сравнений, где — количество новых граней. Но в дальнейшем мне удалось выиграть пару (двоичных =) порядков в скорости поиска соседних использовав упорядоченные ассоциативные массивы (std::set) в этом узком месте. Вообще, есть целая область знаний, которая связана с этой задачей (определения соседей) — это locality-sensitive hashing. В оригинальной реализации qhull используются как раз хеши (не LSH). В ней для каждого ребра (двух граней) хранится отдельная структура данных (список вершин ребра и информация о гранях «над» и «под» ребром), каждая грань содержит список своих рёбер. Для определения соседства граней (среди граней из списка новых граней) для всех рёбер граней создаётся ассоциативный неупорядоченный массив хешей массивов их вершин (с одним пропуском (англ. skip) и одной отброшеной вершиной ). То есть савнений хешей. Если происходит коллизия (и последующее полное совпадение списков) хешей, то (очевидно) такие две грани имеют общее ребро (хеш, кстати, при этом можно удалить для уменьшения хеш-таблицы) и каждую из них можно добавить в список соседних граней другой. Этот способ я не использую, так как ни в стандартной библиотеке, ни в STL нет такого средства, как hash combiner (upd. поступила информация, что лучший выбор — это просто XOR). Самому специализировать без теоретико-информационного обоснования корректности std::hash для std::vector< std::size_t > мне не хотелось. Именно поэтому я использовал ассоциативный упорядоченный контейнер рёбер для поиска соседних граней. В результате (судя по всему) моя реализация алгоритма проигрывает qhull в скорости асимптотически в от до раз (для разного рода входных данных).

Сравнительное тестирование.

Образцом служит «каноническая» реализация алгоритма Quickhull qhull, установленная через менеджер пакетов (в Ubuntu 14.04 LTS). Для генерации входных данных из этого же пакета qhull-bin использовалась утилита rbox. Команда rbox t n D4 30000 s > /tmp/qhuckhull.txt создаёт файл с координатами 30000 точек на четырёхмерной сфере. Команда rbox t n D10 30 s > /tmp/quickhull.txt создаёт файл с координатами 30 точек на десятимерной сфере. Количество памяти, которое затрачивает программа можно посмотреть в выводе утилиты /usr/bin/time с ключом -v. Таким образом в выводе /usr/bin/time -v bin/quickhull /tmp/quickhull.txt | head -7 можно узнать потребление и памяти и процессорного времени (очищенного от времени на чтение файлов и вывод на терминал) для моей реализации, а в выводе /usr/bin/time -v qconvex s Qt Tv TI /tmp/quickhull.txt — для «канонической» реализации qhull.
Корректность своей реализации я определяю по совпадению количества граней выпуклых оболочек. Но для отладки (в режимах двух- и трёхмерном) на некотором этапе я реализовывал пошагово (через команду pause) анимированный вывод в формате gnuplot. Где-то в коммитах он есть. Вывод программы — это выпуклая оболочка и пронумерованное множество входных точек представленные в формате gnuplot.
Кроме того, как-то было начал писать (но не закончил) утилиту randombox — аналог rbox. randombox задумывалась как утилита, которая генерирует равномерно распределённые в пространстве (англ. uniform spatial distribution) точки, в отличие от rbox, которая генерирует точки распределённые неравномерно. randombox может генерировать наборы точек, ограниченные симплексом (любой размерности, которая может быть вложена в пространство с точками), на единичной сфере (англ. unit sphere), в шаре (англ. ball), на поверхности единичного «ромбовидного» многогранника (англ. unit diamond), внутри единичного куба (англ. unit cube), внутри параллелотопа, в пределах стандартного единичного симплекса (англ. unit simplex), а также проецирует множество точек в объём конуса (англ. cone), либо цилиндра (англ. cylinder). Для генерации координат, равномерно распределённых в пределах произвольного сиплекса (вложимого в пространство), я нашёл численно стабильный алгоритм в интернете (через распределение Дирихле), а такеже придумал свой алгоритм, который не отличается таким свойством. В итоге выбрал первый.
Для того, чтобы «пощупать» выпуклые оболочки взглядом, используя gnuplot, из корня проекта введите команду rbox n D3 40 t | bin/quickhull > /tmp/quickhull.plt && gnuplot -p -e «load '/tmp/quickhull.plt'». rbox n D3 40 t сгенерирует 40 точек внутри ограничивающего куба (англ. bounding box). Ключ t задаёт использование текущего времени (в секундах) в качестве начального значения ГПСЧ. Чтобы получить точки на сфере, необходимо добавить ключ s. Также интересно посмотреть как разбиваются несимплициальные грани: rbox n D3 c — куб, rbox n D3 729 M1,0,1 — точек, расположенных в целых узлах. Красиво выглядит выпуклая оболочка спирали (англ. helix): rbox n D3 50 l. В случае одновременного расположения большого количества точек, компланарных граням, большое значение для корректности результатов начинает играть константа eps, которая задаётся при конструировании объекта алгоритма. Например, если использовать в качестве малого значения std::numeric_limits< value_type >::epsilon(), то в случае «повёрнутого» куба из 64 целых точек rbox n D3 64 M3,4 в результате работы алгоритма может получиться некорректный геометрический объект, не являющийся выпуклым многоугольником.

Итог или запоздалый TL;DR.

Итогом моих трудов явилась релизация алгоритма Quickhull на C++, довольно компактная и не сильно более медленная, чем реализация с qhull.org. Алгоритм получает на входе значение размерности пространства, значение малой константы eps и набор точек, представляемый как диапазон, задаваемый парой итераторов так, как это принято в STL. На первом этапе create_simplex строится стартовый симплекс и возвращается точечный базис аффинного (под)пространства, вмещающего входные точки. Если количество точек в базисе больше, чем размерность Евклидова пространства, вмещающего точки, то далее запускается алгоритм достройки выпуклой оболочки. На выходе алгоритм даёт массив структур данных, описывающих грани выпуклой оболочки входного множества, который и является ответом. В случае неудачи получается точечный базис некоторого подпространства меньшей размерности, которое вмещает все точки. Используя алгоритм Хаусхолдера можно некоторым образом повернуть входное множество, занулив при этом старшие координаты точек. Такие координаты можно отбросить и применить алгоритм Quickhull уже в пространстве меньшей размерности. Последнее не имплементировано.

Применения.

У данного алгоритма есть множество приложений. Кроме нахождения выпуклой оболочки как таковой, благодаря тому, что существует связь между выпуклой оболочкой, триангуляцией Делоне и диаграммой Вороного, несложно найти применения алгоритму «опосредованно».
Данная реализация уже кое-кому пригодилась, и была оплачена благодарностью (в Acknowledgements), что очень приятно.

Обновления.

24.06.2015 г.: Добавлена проверка выпуклости получаемой геометрической структуры по алгоритму Kurt Mehlhorn, Stefan Näher, Thomas Schilz, Stefan Schirra, Michael Seel, Raimund Seidel, and Christian Uhrig. Checking geometric programs or verification of geometric structures. In Proc. 12th Annu. ACM Sympos. Comput. Geom., pages 159–165, 1996.
24.06.2015 г.: Грани получаемого многогранника теперь содержат наборы вершин и индексов соседних граней упорядоченные таким образом, что грань является противолежащей (англ. opposite) соответственной вершине. Это делает структуру намного более удобной при использовании в некоторых случаях.
24.06.2015 г.: Теперь для определения соседства вновь получаемых граней используются их хэши.