Приветствую, уважаемые читатели данной статьи! В статье я дам описание имплементации алгоритма Форчуна (англ. Fortune's algorithm) для построения диаграммы Вороного (англ. Voronoi diagram) с использованием нативных сбалансированных двоичных деревьев поиска (для уникальных элементов) (англ. BST, binary search tree), предусмотренных стандартом C++, — ассоциативных упорядоченных контейнеров std::map и std::set.

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

Репозиторий

https://github.com/tomilov/sweepline/

Пока нет релиза, но ветка develop содержит в данный момент стабильную версию. Потом появится и master, и релизы.

Определение диаграммы Вороного

Для некоторого конечного набора попарно различных точек на плоскости (здесь и далее N — количество точек) диаграмма Вороного представляет из себя разбиение плоскости на области — ячейки (англ. cell) диаграммы Вороного. Каждая ячейка содержит в себе лишь одну точку исходного набора точек, называемых сайтами диаграммы Вороного. Все точки ячейки находятся ближе к соответствующему сайту, нежели чем к любому другому. Для обычного определения расстояния на плоскости — метрики Евклида (корень квадратный из суммы квадратов разностей координат точек — из теоремы Пифагора) — форма ячеек в общем случае является выпуклым многоугольником. Конечно же существуют крайние случаи, как если в исходном множестве 1 (ячейка — вся плоскость), 2 точки (ячейка — полуплоскость), либо точки (N > 2) находящиеся на одной прямой (в этом случае внутренними ячейками будут полосы, а внешними — полуплоскости). Крайние из множества ячеек в общем случае являются частью выпуклых многоугольников с двумя сторонами, уходящими в бесконечность, то есть параллельными, либо расходящимися лучами. Стороны многоугольников ограничивающих ячейки (но при этом не являющимися их частью) — это рёбра (англ. edge) диаграммы Вороного. Они могут быть отрезками (англ. line segment), лучами (англ. ray), либо прямыми (англ. line). Каждое ребро — это множество точек, равноудалённых ровно от двух сайтов, то есть лежит на серединном перпендикуляре для двух сайтов. Вершины многоугольников, ограничивающих ячейки, так и называются — вершины диаграммы Вороного (англ. мн. ч. vertices, ед. ч. vertex). Вершины являются точками, равноудалёнными от трёх и более сайтов, то есть являются центрами (англ. circumcenter) описанных окружностей (англ. circumcircle) многоугольников, которые можно было бы построить на сайтах примыкающих ячеек, как на вершинах.

Неполное описание алгоритма Форчуна

Ранее, в опубликованных здесь другими авторами статьях ((1), (2)), уже было приведено описание алгоритма Форчуна и некоторых структур данных, которые используются для его реализации. Вкратце, суть алгоритма в следующем: по плоскости движется заметающая прямая (англ. sweepline). Движется скачками — от точки к точке. Первая часть этих точек — это точки со ввода, которые становятся сайтами. Вторая часть — это "виртуальные" точки, крайние по ходу движения заметающей прямой точки упомянутых описанных окружностей. При движении (параллельном переносе) заметающей прямой она касается любой такой описанной окружности дважды — второй раз эквивалентен событию, при котором диаграмма Вороного достраивается: к ней добавляется вершина, одно или более рёбер оканчиваются этой новой вершиной и одно или два новых ребра выходят из неё.

Не думаю, что приведённое выше описание алгоритма достаточно хорошо, тем более, в голове у меня своя картина того, как это всё происходит. Картина эта сформировалась под влиянием языка C++ и части его стандартной библиотеки — стандартной библиотеки шаблонов — STL (а именно, алгоритмов и контейнеров). И, здесь, мы плавно переходим к следующей части нашего повествования.

Алгоритм и структуры данных

"Натянуть" алгоритм Форчуна на STL — непростая задача. Точнее, совсем не очевидно как именно это сделать. Связано это с хитрой структурой используемой в алгоритме, которая называется "береговая линия" (англ. beachline). Она описывает геометрическое множество точек равноудалённых от заметающей прямой и от сайтов, ячейки которых ещё не были замкнуты. В нашем случае (плоскость с привычной Евклидовой метрикой) береговая линия состоит из дуг (англ. arc) — ограниченных с одной, либо с двух сторон частей парабол, фокусами (англ. ед. ч. focus, мн. ч. focii) которых являются сайты. В ходе перемещения заметающей прямой эти куски парабол растут, делятся (в случае если заметающая прямая пересекает очередной сайт и перпендикуляр к заметающей прямой из этого места попадает в кусок параболы где-то посередине, то этот кусок параболы делится на два и между ними "вклинивается" ещё один — с фокусом в новом сайте), и схлопываются в точку (в этом случае появляется новая вершина и одно или два новых ребра). Куски парабол формирующие береговую линию граничат между собой в endpoint-ах (рус. конец отрезка). Каждый данный endpoint движется по прямой (то есть по ребру) и является точкой, которая равноудалёна в каждый момент: от заметающей прямой и от каждого из двух сайтов, которые являются в свою очередь фокусами граничащей через этот endpoint двойки кусков парабол.

Так вот, эта структура, во-первых, должна быть упорядоченной, во-вторых, она обычно является у авторов существующих имплементаций гетерогенной, то есть содержит как куски парабол, так и endpoint-ы в качестве ключей/элементов. Для того, чтобы имплементация имела среднюю и максимальную асимптотическую сложность O(N * log(N)), с необходимостью придётся использовать именно сбалансированные сортирующие деревья. Для хранения гетерогенных данных можно использовать std::variant. Последний, впрочем, мне не понадобится.

std::map и неточный пользовательский "прозрачный" компаратор

Казалось бы (!) контейнеры STL std::set и std::map не очень подходят для хранения береговой линии, но, как оказывается, их можно использовать на полную катушку. Гетерогенности данных в узлах деревьев можно полностью избежать. Всё дело в том, что в качестве ключей хранить можно именно endpoint-ы — упорядоченные пары указателей (или итераторов, указывающих) на сайты — "левый" и "правый" (в некотором смысле, который станет ясным далее). В этом заключается неизбежная избыточность: каждая пара соседних endpoint-ов содержит одинаковую пару указателей (итераторов, указывающих) на сайт.

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

Современные реализации стандартных библиотек C++ (по крайней мере libc++ и libstdc++) имеют имплементацию интересующего меня контейнера std::map, которая позволяет производить гетерогенный поиск (англ. heterogeneous lookup), а именно, функции: count, find, equal_range, lower_bound, upper_bound, — имеют среди перегрузок (англ. overloadings) шаблон функции, который позволяет производить поиск элементов, эквивалентных значению, имеющему тип отличный от типа значения (или от типа ключа значения, в случае отображения) хранящегося в контейнере. Данная перегрузка каждой из функций может быть использована только если компаратором (функтором, используемым для сравнения элементов) является "прозрачный компаратор", то есть такой компаратор, который в своём пространстве имён определяет тип с именем is_transparent, будь то таг (имя класса (struct, class, union) или перечисления (enum или enum class)), typedef или псевдоним типа (англ. type alias), как это приведено в следующем коде:

struct less
{
    using is_transparent = void; // int, less, std::true_type, std::false_type...
    bool operator () (T const & lhs, T const & rhs) const
    {
        return /* вернуть истину, если lhs строго меньше rhs */;
    }
};

Как видно, тип имени (англ. qualified-id) less::is_transparent может быть незавершённым (англ. incomplete). Например, void или объявление имени может являться forward-declaration класса/перечисления.

Для того чтобы компаратор (less в примере) мог сравнивать значения типа T, находящиеся в контейнере, со значениями некоторого другого типа U, необходимо доопределить компаратор less ещё парой операций, эквивалентных less::operator () (T, U) const и less::operator () (U, T) const для каждого такого типа U.

Эквивалентность с допуском eps

Пользовательский (англ. custom) компаратор здесь ещё необходим в связи со следующими соображениями: в случае, если для промежуточных вычислении используются числа с плавающей запятой (англ. floating-point numbers) имеющие ограниченную точность, то алгоритм параметризуется положительной константой eps, которая задаёт точность сравнения. Числа, которые различаются на величину меньшую чем eps считаются эквивалентными:

double eps = 0.2;
double a = 2.0;
double b = 2.1;
assert(!(a + eps < b) && !(b + eps < a));

В assert-ионе стоит как раз утверждение (определение) эквивалентности a и b. В данном приложении — эквивалентности с точностью eps.

Вот пример вывода списка значений ключей (действительных (с нек. оговорками) чисел) всех элементов содержащихся в контейнере std::map или std::set и эквивалентных целому числу 1 (взято отсюда):

#include <set>
#include <map>
#include <iterator>
#include <algorithm>
#include <iostream>

double const eps = 0.2;

struct less
{
    bool operator () (double l, double r) const { return l < r; }
    using is_transparent = void;
    bool operator () (int l, double r) const { return l + eps < r; }
    bool operator () (double l, int r) const { return l + eps < r; }
};

int main()
{
    {
        std::set< double, less > s{0.0, 0.9, 1.0, 1.1, 2.0};
        auto p = s.equal_range(1);
        std::copy(p.first, p.second, std::ostream_iterator< double >(std::cout, " "));
        std::cout << std::endl;
    }
    {
        struct A {}; // irrelevant
        std::map< double, A, less > m{{0.0, {}}, {0.9, {}}, {1.0, {}}, {1.1, {}}, {2.0, {}}};
        auto p = m.equal_range(1);
        for (auto it = p.first; it != p.second; ++it) {
            std::cout << it->first << ' ';
        }
    }
}

Должно вывести 0.9 1 1.1 для каждого контейнера. Ясно, что даже в контейнере, содержащем уникальные значения какого-то типа T, может оказаться более одного элемента, эквивалентного некоторому значению некоторого типа U, отличного от T. При этом для работы алгоритмов std::map::{count,upper_bound,lower_bound,equal_range} должно выполняться условие — элементы контейнера должны быть хотя бы частично упорядочены (англ. partitioned) при сравнении со значением типа U, как это описано в этом ответе на SO.

Вставка с подсказкой

Полезная особенность реализации ассоциативных упорядоченных контейнеров STL — это возможность вставить элемент с подсказкой (англ. hint). Функции вставки std::map::insert передаётся параметром кроме вставляемого значения ещё и итератор на элемент hint, перед которым необходимо вставить элемент. В случае успешной вставки стандартом даётся гарантия, что будет произведено лишь ограниченное количество сравнений элементов контейнера (м/у собой (хотя зечем?) и со вставляемым). Здесь я, признаюсь, использую неопределённое поведение (англ. UB, undefined behaviour): дело в том, что есть некоторые требования к компаратору — Comparator concept — которые должны выполняться всегда, что бы мы ни делали с контейнером. Однако ж, логически можно заключить (и так на практике и получается с libc++), что при вставке уникального элемента в правильное место по подсказке необходимо максимум два сравнения вставляемого элемента с предыдущим *std::prev(hint) элементом и с элементом, на который указывает итератор hint, если таковые имеются. Если контейнер пустой, то сравнений не требуется, а если вставка в начало, либо в конец, то только лишь одно сравнение. Ни одна из двух реализаций стандартной библиотеки, протестированных мною, никогда не пытается сравнить ни один из элементов сам с собой или, в случае со вставкой с подсказкой, — выполнять сравнения элементов отличные от абсолютно логически необходимых для подтверждения правильности вставки. В общем везде и всюду торжествует здравый смысл. Но это всё же UB, если вы в своём коде полагаетесь на это. Я предупредил.

Лексикографическая упорядоченность

Для того, чтобы имплементация алгоритма Форчуна обладала внутренней красотой, необходимо, чтобы на вход алгоритму подавались точки упорядоченные лексикографически (англ. lexicographically ordered). Лексикографическая упорядоченность означает, что сначала две точки сравниваются по абсциссе (x, "икс"), а потом, в случае если их первые координаты эквивалентны, то по ординате (y, "игрек"). Причём в моём случае — с допуском eps.

Две точки плоскости с координатами (lx, ly) и (rx, ry) можно сравнить лексикографически с допуском eps так:

bool less(double l,
         double eps,
         double r)
{
    return l + eps < r;
}

bool less(double lx, double ly,
         double eps,
         double rx, double ry)
{
    if (less(lx, eps, rx)) {
        return true;
    } else if (less(rx, eps, lx)) {
        return false;
    } else {
        return less(ly, eps, ry);
    }
}

Такой компаратор задаёт квадратную (со стороной 2 * eps) окрестность.

Лексикографически без допуска можно сравнить две точки, используя готовую реализацию из STL так:

#include <tuple>

bool less(double lx, double ly,
         double rx, double ry)
{
    return std::tie(lx, ly) < std::tie(rx, ry);
}

Подробнее об этом здесь. Как std::tuple, так и std::pair имеют именно такие операции operator <.

Ввод и вывод, форматы данных

Весь ввод предполагается лексикографически упорядоченным и не содержащим эквивалентных с точностью eps точек. Так как сортировка ввода — это тоже O(N * log(N)), то эту (логически совершенно отдельную) часть я выношу за пределы имплементации. Это позволяет использовать на входе прокси-итераторы (есть в примере), то есть сортировать итераторы указывающие на тяжеловесные точки, а не сами точки, перед отправкой в алгоритм.

Алгоритм выглядит как алгоритмы STL — на вход принимает пару итераторов [first,last), задающих диапазон точек (сайтов) для обработки.

Формат точки:

struct point { value_type x, y; };

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

Результатом является граф: вершины содержатся в (автоматически упорядоченном в ходе работы алгоритма) списке vertices_, рёбра — в двусторонней очереди edges_. Контейнеры для вершин (std::list) и рёбер (std::deque) выбраны такими, чтобы не было проблем с недействительными итераторами (англ. invalid) при вставке и удалении.

Формат ребра:

struct edge { site l, r; pvertex b, e; };

Где site — указатель (или итератор, указывающий) на сайт, pvertex — итератор, указывающий на вершину. Сами вектора в паре (l, r) и (b, e) упорядочены по часовой стрелке, т.е. так, что [(l, r) x (b, e)] < 0 (здесь — векторное произведение). Последнее даёт возможность всё правильно нарисовать в конце и, вообще, хранить всю полноту информации о графе.

Вершина nov = std::end(vertices_) указывает на бесконечность. Копировать и перемещать данную пару контейнеров так просто нельзя. Так как ни при копировании (англ. copy), ни при перемещении (англ. move) действительные (англ. valid) итераторы, которые нельзя разыменовывать (англ. not dereferencable), в их числе и end-итератор nov, не сохраняются. Возможно только глубокое копирование (англ. deep-copying, cloning (есть в примере)), которое здесь имеет квадратичную сложность. Хотя есть один обходной манёвр. В самом начале сохранить фейковый элемент в vertices_ и использовать его итератор как указатель на бесконечность. Я не пробовал, но это принципиально осуществимо, но только для копирования/перемещения, не для клонирования.

Другие особенности имплементации

Обычно алгоритм Форчуна ведёт очередь для событий двух типов. Это так называемые события круга (англ. circle event) и события точки (англ. point event).

События точки возникают, когда заметающая прямая пересекает очередной сайт. Они неизбежны.

События круга возникают, когда схлопывается какой-то кусок параболы. В очередь они добавляются, когда три сайта a, b и c в порядке их перечисления в структуре береговой линии образуют двойку векторов, упорядоченную по часовой стрелке: [(a, b) x (b, c)] < 0 (здесь — векторное произведение). Но события круга могут быть удалены, если схлопывается какой-то из соседних кусков парабол, либо один из кусков парабол расчленяется куском параболы, вновь созданным для пересекаемого заметающей прямой сайта. Поэтому, вместо зачастую используемой обёртки контейнеров std::priority_queue, я использую всё тот же контейнер std::map или std::set. Сложности нужных операций для этих контейнеров те же самые, что и для std::priority_queue, но они позволяют удалить не только наибольший/наименьший элемент, но и любой другой (по итератору — за O(1) времени). Обычно события из очереди с приоритетами не удаляют, а помечают как недействительное (англ. disabled). Когда такое событие встречается — его пропускают.

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

struct vertex { struct point { value_type x, y; } c; value_type R; };

Где (x, y) — собственно сама вершина, а R — расстояние до сайтов соседних ячеек.

Для того, чтобы использовать лексикографическую упорядоченность всего и вся в моём алгоритме, заметающая прямая параллельна оси ординат и движется строго в направлении оси абсцисс. Таким образом все события круга возникают при координатах c.x + R для соответствующих вершин.

Особенные случаи

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

Особенные случаи, которые требуют особого внимания — это случаи низкой размерности ввода, либо локальные симметрии и одновременные события:

• 0, 1, 2 точки
• 2 и более первые точки со ввода лежат на одной вертикальной прямой (x = const) (при добавлении каждого следующего такого сайта создаётся только одно ребро и только один соответствующий endpoint добавляется в береговую линию).
• более, чем три точки расположены концентрически (при событии сходится более чем два endpoint-а, создаётся одно новое ребро и соответствующий endpoint добавляется в береговую линию)
• при движении заметающей прямой очередной сайт попадает в один endpoint (автоматически возникает и завершается событие круга порождается новая вершина и два новых ребра исходят из неё, два соответствующих endpoint-а добавляются в береговую линию)
• при движении заметающей прямой очередной сайт попадает в множество endpoint-ов, которые бы на следующем шагу алгоритма сошлись в новой вершине — в этом случае новый сайт эквивалентен уже лежащему в очереди событий круга самому ближайшему событию (это событие автоматически завершается, порождая новую вершину и два новых ребра, исходящих из неё, два соответственных endpoint-а добавляются в береговую линию)

Два последних особенных случая как раз выявляются гетерогенным сравнением сайта с endpoint-ами в структуре береговой линии. Посредством алгоритма std::map::equal_range за log(size(endpoints_)) времени можно найти весь диапазон endpoint-ов, сходящихся в вершине, эквивалентной текущему сайту. Если диапазон ([first, last)) не пустой, то мы попали в ребро (assert(std::next(first) == last);), либо в в несколько рёбер (assert(1 < std::distance(first, last));), а если пустой (assert(first == last)), то мы просто попали в кусок параболы где-то посередине. Сам этот кусок расчленяется на два куска, между которыми вставляется новый цельный кусок параболы (которая в этот момент вырождена в луч, перпендикулярный заметающей прямой) с фокусом в данном сайте. При этом создаётся только одно ребро и два соответствующих ему endpoint-а добавляется в береговую линию — это самый часто встречающийся в общем случае вариант. Только лишь в этом (и во втором особенном случае) создаётся ребро без хотя бы одного уже известного конца (т.е. ассоциированной вершины).

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

Вставка endpoint-ов

Как я уже упомянул выше, мой пользовательский компаратор, используемый в дереве представляющем береговую линию, является/приводит к UB, при сравнении endpoint-ов. Для гетерогенной части же — всё в строгом соответствии со Стандартом — language lawyer носа не подточит.

Алгоритм устроен так, что при вставке по подсказке hint, новые endpoint-ы вставляются справа налево в только правильные места (выявленные на предыдущих шагах), и, так как в случае вставки только одного endpoint-а указатели (итераторы) на хотя бы один разграничиваемый сайт присутствует в любом из соседних endpoint-ах (либо hint, либо std::prev(hint)), то всегда можно понять только сравнивая эти указатели (итераторы), находится ли сравниваемый endpoint левее (ниже) или правее (выше) в структуре береговой линии (прим. мысленно я всегда поворачиваю всю картину на +90°).

Однако же, есть и вариант, когда мы вставляем два endpoint-а, которые возникают в самом последнем особенном случае. Среди четырёх endpoint-ов, вовлечённых в операцию вставки по подсказке (два вставляемых плюс два граничных), определённо есть такие, которые не имеют общих разграничиваемых сайтов. Но и здесь есть простой выход. Когда мы вставляем endpoint-ы (как обычно: сверху вниз), то каждый из вставляемых endpoint-ов имеет одним из сайтов тот, который и вызвал это событие точки. Лексикографически сравнивая по одному сайту, лексикографически большему из двух для каждого endpoint-а, мы можем выяснить какой из endpoint-ов "правее", а какой — "левее" в структуре береговой линии.

Пучки

В случае концентрически расположенных сайтов, в списке, хранящем указатели на endpoint-ы, сходящиеся при соответствующем событии, сами эти указатели находятся в некотором беспорядке. Чтобы "задёшево" понять который из сходящихся endpoint-ов самый "левый", а который самый "правый", можно, пользуясь тем, что все соответствующие сайты находятся на одинаковом расстоянии от общей вершины, сравнивать углы векторов (либо углы серединных перпендикуляров пар сайтов, используя преобразование координат для поворота +pi / 2 <=> x <- y, y <- -x), для чего удобна C-функция std::atan2, правильно сопоставляющая знаки и квадранты.

Визуализация

Для визуализации я использую gnuplot. Он вполне быстро отрисовывает референтную для меня диаграмму Воронного для 100'000 сайтов.

Сам алгоритм обрезки рёбер по размерам ограничивающего прямоугольника (англ. bounding box) также есть в репозитории.

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

Сложный ввод

Сложными для обработки плохими реализациями являются следующие конфигурации точек на входе алгоритма:

• прямоугольная сетка (англ. rectangle grid) — точки расположены через равные шаги по обеим координатам без смещения
• диагональная сетка (англ. diagonal grid) — повёрнутая на pi / 4 прямоугольная сетка
• гексагональная сетка (англ. hexagonal grid) — точки расположены в вершинах правильных (равносторонних) треугольников, которыми замощена нужная часть плоскости
• треугольная сетка (англ. triangular grid) — точки расположены в вершинах правильных шестиугольников, которыми замощена нужная часть плоскости

Крепкий орешек (нужна помощь)

Каждая из перечисленных выше сеток выявила тонны багов на промежуточных стадиях моей имплементации. Последняя — это самая жесть, скажу я вам! Она очень быстро выявляет неточность в вычислениях с плавающей запятой. Уже на небольших разбросах диапазонов абсолютных значений координат разница в координатах, которые вычисляются для центра и радиуса (англ. circumradius) описанной окружности по этим формулам (вычисление вершины — ф-я make_vertex) и по выведенным вручную формулам для вычисления ординаты пересечения двух горизонтальных парабол (ф-я intersect из пользовательского компаратора для endpoint-ов), достигает зловеще больших значений! Бесконтрольное увеличение eps — это не выход, я уверен.

Не смотря на то, что в каждой из этих операций я добился использования только одного (SIMD :-) ) деления и одного корня квадратного, точность их разнится весьма сильно. Вроде бы и std::sqrt является точной операцией. Цитата отсюда:

std::sqrt is required by the IEEE standard be exact.

Однако существенное расхождение — установленный факт. Что делать? Я не знаю. Есть здесь специалисты по стабильным вычислениям с плавающей запятой? Помогите — в благодарность запишу в соавторы в README в репозитории. Цель — победить треугольную сетку на нативных типах чисел с плавающей запятой.

Тест производительности

Производительность я измеряю для 100'000 точек. Предвосхищая вопрос "Почему именно для ста тысяч?" я сразу сообщаю: потому что в этой задаче, которая является топ 5 в рейтинге Тимуса по сложности фигурирует эта цифра.

Так вот, производительность для 100'000 точек в общих позициях с аллокатором TCMalloc из Google Performance Tools на довольно древнем ноутбуке с Intel(R) Core(TM) i7-2670QM CPU @ 2.20GHz выражается в величине 205 миллисекунд. На серверном железе без этого хитрого аллокатора она будет примерно такой же (проверено).

Кстати, arena-аллокатор даёт такую же прибавку как и упомянутый TCMalloc.

ADDITIONAL

1. Поборол треугольную сетку этим коммитом. Но проблемы с точностью вычислений никуда не исчезли и всё так же присутствуют в алгоритме. Все соображения mcquay ничуть не потеряли актуальность. Это дало некоторую оптимизацию и теперь время для 100'000 точек — 197 мс.