Рассмотрим структуру бинарного пространственного разбиения, называемая kd-дерево. Эта структура представляет собой бинарное дерево ограничивающих параллелепипедов, вложенных друг в друга. Каждый параллелепипед в kd-дереве разбивается плоскостью, перпендикулярной одной из осей координат на два дочерних параллелепипеда.

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

Построение kd-дерева

Алгоритм построения kd-дерева можно представить следующим образом (будем называть прямоугольный параллелепипед англоязычным словом "бокс" (box)).

"Добавить" все примитивы в ограничивающий бокс. Т.е построить ограничивающий все примитивы бокс, который будет соответствовать корневому узлу дерева.

Если примитивов в узле мало или достигнут предел глубины дерева, завершить построение.

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

Добавить примитивы, пересекающиеся с боксом левого узла в левый узел, примитивы, пересекающиеся с боксом правого узла в правый.

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

1. Выбрать плоскость разбиения по центру.

Самый простой способ - выбирать плоскость по центру бокса. Сначала выбираем ось (x, y или z), в которой бокс имеет макимальный размер, затем разбиваем бокс по центру.

Рисунок 1 : Разбиение по центру

Этот способ имеет плохую адаптивность. Классическое октодерево имеет те же недостатки. Интуитивно можно описать причину плохой скорости на таких деревьях тем, что в каждом узле достаточно много пустого пространства, через которое луч траверсится (проходи через дерево во время поиска) в то время как пустое пространство должно сразу-же по-возможности отбрасываться. Более строгое объяснение будет дано чуть позже.

2. Выбрать плоскость по медиане.

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

Рисунок 2 : Разбиение по медиане

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

3. Surface Area Heuristic (SAH)

Каковы же критерии хорошо-построенного kd-дерева? На интуитивном уровне это на самом деле это не очень понятно, хотя выражение "как можно больше пустого пространство должно быть отброшено как можно быстрее" наверное подойдет.

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

SAH(x) = CostEmpty + SurfaceArea(Left)*N(Left) + SurfaceArea(Right)*N(Right)

Где CostEmpty - стоимость прослеживания пустого узла (некоторая константа), SurfaceArea - площадь поверхности соответствующего узла, а N - число примитивов в нем. Нужно выбирать плоскость разбиениея так, чтобы минимизировать эту функцию. Аргумент функции x является одномерной координатой плоскости разбиения.

Хорошими кандидатами на минимум SAH могут служить границы примитивов. Простой алгоритм посторения выглядит следующим образом. Каждый раз при выборе плоскости нужно перебрать всевозможные границы примитивов по трем измерениям, вычислить в них значение функции стоимости и найти минимум среди всех этих значений. Когда мы вычисляем SAH для каждой плоскости, то нам необходимо знать N(left) и N(right) - количества примитивов справа и слева от плоскости. Если вычислять N простым перебором, в итоге получится квадратичный по сложности алгоритм построения.

Рисунок 3 : Разбиние с учетом SAH

На рисунке 4 можно увидеть, что SAH сразу отбрасывает большие пустые пространства, плотно ограничивая геометрию.

Рисунок 4 : kd-дерево, построенное с учетом SAH

Surface Area Heuristic порождает более интеллектуальный критерий остановки, чем предел глубины дерева или малое количество примитивов. Если для выбранной плоскости разбиения суммарная стоимость прослеживания дочерних узлов больше, чем стоимость прослеживания текущего узла как листа (то есть по формуле SurfaceArea(Node)*N(Node)), то посторение дерева следует остановить.

4. Быстрое построение kd дерева

Вариант 0

Разбивайте по центру, выбирая ось, по которой размер бокса максимален. Этот способ самый быстрый, но на некоторых сценах трассировка лучей в таком дереве примерно в 2-3 раза медленнее.

Вариант 1

Самое простое что можно сделать чтобы ускорить построение - это снизить количество перебираемых плоскостей. Вместо того, чтобы перебирать N плоскойстей (где N - количество примитивов), можно вычислить SAH лишь в небольшом числе заранее фиксированных мест. Входной бокс равномерно разбивается по каждой оси на M частей. Обычно M лежит в пределах от нескольких десятков до пары сотен. N(left) и N(right) попрежнему вычисляются перебором, но теперь нам придется сделать полный перебор по массиву всего лишь M раз, а не N. Таким образом, алгоритм приобретает сложность N*M. Фактически, мы достигли ускорения за счет огрубления SAH, дискретизовав ее. Однако оказывается, что значение минимума полученной грубой функции как правило не сильно отличается от ее точного минимума.

Вариант 2

Что можно сделать, чтобы ускорить вариант 1? Нужно избавиться от полного перебора во время вычисления N(left) и N(right). Для этого построим дерево, в каждом узле которого записана некоторая разбивающая плоскость и количество примитивов справа и слева от плоскости. struct IntervalTree

Float split;

Int numPrimsOnLeftSide;

Int numPrimsOnRightSide;

IntervalTree* left;

IntervalTree* right;

На самом деле нам потребуется три таких дерева на каждом уровне - по одному на x,y,z. Входной интервал каждый раз будем делить пополам. Имея такое дерево, можно за log(N) операций достаточно точно определить количество примитивов справа и слева от плоскости. Алгоритм поиска в дереве выглядит достаточно просто.

vec2i IntervalTree::TraversePlane(

float a_split) const

// в нулевой компоненте будем хранить минимум, в первой - максимум

vec2i res; // вектор из двух целых чисел

if (a_split < split - EPSILON)

res = numPrimsOnRightSide;

if (left!=0)

res += Left()->TraversePlane(a_split);

// посчитаем хотя бы те примитивы, которые есть в данном узле, чтобы SAH не занулялся.

If (res.M == 0)

res.M = numPrimsOnLeftSide;

Else if (a_split > split + EPSILON)

res = numPrimsOnLeftSide;

If (right!=0)

res += right->TraversePlane(a_split);

// если по какой-то причине количество примитивов равно нулю, то

// посчитаем хотя бы те примитивы, которые есть в данном узле, чтобы SAH не занулялась.

If (res == 0)

res = numPrimsOnRightSide;

Else

// попали точно в плоскость узла

res = numPrimsOnLeftSide;

res = numPrimsOnRightSide;

return res;

Сложность постоения дерева плоскостей - N*log(N). Сложность оценки SAH - M*log(N), так как мы вычисляем SAH только в M фиксированных плоскостях. Таким образом, суммарная сложность алгоритма - N*log(N), т.к. M много меньше N.

Перед тем как выполнять построение kd дерева можно впринципе построит произвольную ускоряющую структуру только для того, чтобы ускорить впоследствие вычисление SAH.

Вариант 3 (Почти точно и за N*Log(N))

Используем сортировку. Для каждой оси (x,y,z) отсортируем примитивы сначала по максимумам а затем по минимумам их ограничивающих боксов. Назовем массив, отсортированный по максимумам sorted_max. Массив, отсортированный по минимумам - sorted_min.

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

for (int axis = 0; axis < 3; axis ++)

// sort primitives belong to current axis and min bounds

// CompareMin comparator (axis ); std :: sort (prims . begin (), prims . end (), comparator ); for (int i =1; i < prims . size (); i ++) int onLeftSide = i ; int onRightSide = prims . size ()- i ; float split = prims [ i ]. Box (). vmin [ axis ]; if (split == a_box . vmin [ axis ]) continue ;

// evaluate SAH

AABB3f box_left = a_box ;

AABB3f box_right = a_box ;

Box_left . vmax [ axis ] = split ;

Box_right . vmin [ axis ] = split ;

Float sah = EMPTY_COST + onLeftSide * SurfaceArea (box_left ) + onRightSide * SurfaceArea (box_right );

If (sah < min_sah )

Min_sah = sah ;

Min_split = split ;

Min_axis = axis ;

// sort primitives belong to current axis and max bounds

// CompareMax comparator2 (axis ); std :: sort (prims . begin (), prims . end (), comparator2 ); for (int i = prims . size ()-2; i >=0; i --) int onLeftSide = i +1; int onRightSide = prims . size ()- i -1; float split = prims [ i ]. Box (). vmax [ axis ]; if (split == a_box . vmax [ axis ]) continue ;

// evaluate SAH

AABB3f box_left = a_box ;

AABB3f box_right = a_box ;

Box_left .vmax [axis ] = split ;

Box_right .vmin [axis ] = split ;

Float sah = EMPTY_COST + onLeftSide *SurfaceArea (box_left ) + onRightSide *SurfaceArea (box_right );

If (sah < min_sah )

Min_sah = sah ;

Min_split = split ;

Min_axis = axis ;

В теории, у метода есть проблема. Рассмотрим массив sorted_min. Проходим по нему слева направо в цикле:

for (int i=0;isize ();i++)

onLeftSide = i ;

Int onRightSide = prims .size ()-i ;

// ...

Число onLeftSide мы знаем абсолютно точно, а вот число onRightSide - не совсем. Дело в том, что плоскость может разбивать некоторые примитивы и в этом случае один и тот же примитив лежит как справа от плоскости, так и слева, что данный алгоритм не учитывает. На практике эта проблема практически себя не проявляет.

Вариант 4

Можно ускорить поиск минимума ф-ии SAH, использовуя како-нибудь метод минимизации с несколькими начальными приближениями. Скажем, метод золотого сечения. В данном случае мы так же избавляемся от метода полного перебора. Нужно только учитывать, что SAH приобретает более-менее гладкую форму только при большом числе примитивов. Плюс этого подхода в том, что вычисление SAH брутфорсом легко перекладывается на GPU. Поэтому оценивая каждый раз SAH грубой силой и сведя количество оценок к небольшому числу (~10-20 max), можно строить kd дерево таким методом очень быстро.

Вариант 5 (Binning)

Часто используют именно этот вариант . Ссуть состоит в следующем:

Нужно поделить пространство на регулярные интервалы по x,y и z. Каждый такой интервал называется корзиной (bin). Обычно ограничиваются небольшым числом корзин ~32.

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

Прослеживание луча в kd-дереве на CPU

скачать алгоритм можно отсюда

Классический алгоритм бинарного поиска в kd деревьях (kd-tree traversal в англоязычной литературе), применяющийся в большинстве CPU реализаций, состоит примерно в следующем. На первом шаге алгоритма необходимо посчитать пересечение луча с ограничивающим сцену корневым параллелепипедом и запомнить информацию о пересечении в виде двух координат (в пространстве луча) – t_near и t_far, обозначающих пересечение с ближней и дальней плоскостями соответственно. На каждом следующем шаге необходима информация только о текущем узле (его адрес) и этих двух координатах.

При этом нет необходимости вычислять пересечение луча и дочернего параллелепипеда, достаточно лишь узнать пересечение с разбивающей параллелепипед плоскостью (обозначим соответствующую координату как t_split). Каждый нелистовой узел kd дерева имеет два дочерних узла. На рис. 5. изображены три варианта событий, которые могут возникать при прослеживании пути луча.

Рисунок 5 : Прослеживание луча в kd дереве

В случае A (t_split >= t_far) или B (t_split < t_near) луч пересекает только один дочерний узел, поэтому можно просто отбросить правый (соответственно левый) узел и продолжить поиск пересечения в оставшемся узле.

В случае C луч пересекает оба дочерних узла, поэтому необходимо сначала поискать пересечение в ближнем узле и если оно не найдено, искать его в дальнем. Так как в общем случае неизвестно, сколько раз произойдет последнее событие, необходим стек. Каждый раз, когда луч пересекает оба дочерних узла, адрес дальнего узла, t_near и t_far помещаются в стек и поиск продолжается в ближнем. Если в ближнем узле пересечение не найдено, из стека достаются адрес дальнего узла, t_near, t_far и поиск продолжается в дальнем узле.

Прослеживание луча в kd-дереве на GPU

Так как на GPU изначально отсутствует стек, появились бесстековые алгоритмы и алгоритмы, использующие искусственный стек небольшой длины. На GPU известны пять алгоритмов прослеживания пути луча - restart, backtrack, push- down, short stack и алгоритм прослеживания в kd дереве со связками.

kd-tree restart

Рисунок 2 :работа restart алгоритма .

Когда луч не находит пересечение в листе, его t_far приравнивается к t_near и поиск продолжается с корня kd-дерева (рис. 2). Грубо говоря, просто подвигается origin луча - то есть его точка испускания и поиск начинается сначала. Это приводит к тому, что луч проходит по одним и тем же узлам много раз, что неэффективно .

kd-tree push-down

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

То есть, если при спуске по ближайшим узлам kd дерева хотя бы один раз встретился узел, в котором луч пересекает и ближний и дальний дочерние узлы, то именно этот дальний дочерний узел должен выть выбран в качестве поддерева. Далее, если луч промахнулся, рестарт будет произведен с запомненного дальнего узла и можно опять попытаться найти новое поддерево . Фактически, это попытка завести стек длинной 1 элемент.

kd-tree short stack

Авторы также использовали короткий стек. Пока размера стека хватает, идет его заполнение аналогично классическому алгоритму. Когда стек переполняется, он начинает работать как кольцевой буффер. Если стек пуст, необходимо произвести рестарт. Например, если в стеке длинной 4 содержаться узлы с номерами (1), (4), (6), (8), то при добавлении нового элемента (12), стек будет иметь следующий вид: (12), (4), (6), (8). То есть будет затерт первый элемент. Удаляться элементы будут в том порядке в каком они добавлялись (то есть сначала 12, потом 8, 6 и наконец 4), но когда из стека будет удален элемент (4), нужно будет произвести рестарт, так как мы затерли элемент (1). Смысл короткого стека в том, что он сильно сокращает количество restart-ов, которое приходится делать лучу.

kd-tree backtrack

Если сохранять в узлах kd дерева ссылку на родительские узлы, в случае промаха по этой ссылке можно добраться до родителя. Кроме того, потребуется хранить в каждом узле его ограничивающий паралеллипипид, чтобы иметь возможность вычислить новый t_far в случае промаха.

Ct_near можно поступить так же, как и в случае restart алгоритма . Это потребует дополнительно 6*4 (для координат паралеллипипида) + 4 (для ссылки на родителя) = 28 байт. Так как память на GPU довольно ограничена, такое расточительство может создать проблемы. К тому же, каждый раз при подъеме придется считать пересечение луча и выровненного по осям паралеллипипида, что конечно не бесплатно по вычислительным ресурсам. Следует особо отметить, что kd дерево с сохраненными паралеллипипид будет занимать намного больше памяти, чем хорошо построенное BVH дерево в этой же сцене. Основная причина здесь в том, что в kd дереве паралеллипипиды будут иметь больше общих точек, которые в итоге будут дублироваться.

kd-tree со связками

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

Этот алгоритм, как и backtrack позволяет избежать многократных прохождений по одним и тем же узлам дерева. Однако, шесть ссылок требуют дополнительно 24 байта памяти, что в сумме с имеющимися 8 дает 32 байта.

Р исунок 3 : kd дерево со связками .

Преимущества kd деревьев

Очень простой и эффективный алгоритм траверса. Даже для GPU.
Занимают мало памяти (8 байт на узел).

Недостатки kd деревьев

Трудоемкое построение, а именно, поиск разбиения с минимальным SAH.

Имеет большую глубину, чем BVH. Больше шагов при построении.

Заключение

Резюмируя, можно сказать, что kd дерево идеально для трасссировки лучей. Это верно как для CPU так и для GPU.

Литература

Wald I. Realtime Ray Tracing and Interactive Global Illumination. PhD thesis, Saarland University, 2004.

Shevtsov M., Soupikov A., Kapustin A.
Highly Parallel Fast KD-tree Construction for Interactive Ray Tracing of Dynamic Scenes. In Proceedings of the EUROGRAPHICS conference, 2007.

Foley T., Sugerman J. KD-Tree Acceleration Structures for a GPU Raytracer. In Proceedings of the ACM SIGGRAPH/EUROGRAPHICS conference on Graphics hardware, p. 15-22, 2005.

Horn D., Sugerman J., Houston M., Hanrahan P. Interactive k-D Tree GPU Raytracing. Proceedings of the symposium on Interactive 3D graphics and games on Fast rendering, p. 167-174, 2007.

Popov S.,Günther J.,Seidel H.-P.,Slusallek P. Stackless KD-Tree Traversal for High Performance GPU Ray Tracing. In Proceedings of the EUROGRAPHICS conference, 2007.

Начало поста для разнообразия написано в духе "Серёжа, ты же учишь математику! Пиши строго.", не обращайте внимания, дальше всё нормально.
Да, и я не имел ввиду, что у меня куча опыта по этому вопросу, у меня совсем чуть-чуть опыта, но меня попросили рассказать.

Пусть у нас есть множество возможных элементов X, выбранное подмножество A, а задача заключается в том, чтобы по предъявленному x из X найти сколько-то "наиболее похожих" на него в A. Например, десять наиболее "похожих". Или все, чья "похожесть" не меньше некой заданной. Желательно не перебирая всё A, это долго.

Задача отлично и, главное, быстро решается, если элементы это числа, а "похожесть" тем больше, чем ближе друг к друг значения: просто сортированный массив и бинарный поиск, или дерево поиска (двоичное, n-ичное, B-дерево -- не важно).

Что в данном случае позволяет решать задачу быстро: наличие заданного на множестве порядка, согласованного с "похожестью". То есть, если искомый элемент x < a, и a < b, то x более похож на a, чем на b. Это позволяет не рассматривать сильно большие/меньшие элементы, так как они заведомо хуже подходят.

Если определить "похожесть" иначе, например через расстояние Левенштейна между десятичной записью, наличие порядка перестаёт помогать. То есть дело именно в согласованности одного с другим, а не в возможности задать порядок хоть как-нибудь, тем более что "хоть как-нибудь" можно всегда, есть такая теорема в теории множеств .

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

На случай n-мерного пространства и задача и решения обобщаются тривиально, ну и вроде бы можно даже ещё немного дальше, об этом ниже.

Идеи

Первая мысль, которая приходит в голову нормальному человеку, который не хочет выдумывать сложного, это "а давайте разобьём на клеточки". То есть, плоскость разбивается на квадраты, имеющиеся точки привязываются к соответствующим квадратам, по точке x определяется квадрат, в котором она лежит, поиск ведётся по нему и нескольким соседним. Работает не слишком хорошо. Проблема в том, что это способ разбиения "не зависящий от данных", за счёт этого он очень простой, за счёт этого же он может очень плохо соответствовать конкретному набору точек. Например, если выяснилось, что 99% точек сосредоточено в одной клетке, перебирать всё равно придётся всё, а мы так старались.

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

Начинаем с одной большой клетки.

Если в какой-то клетке слишком много точек для перебора, делим её на две части. Для того чтобы было просто приписывать точки к клеткам, делим вертикальными и горизонтальными прямыми, например по очереди то так, то этак (где именно их проводить -- отдельный вопрос).

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

Получилась очевидная иерархическая структура: клетки верхнего уровня содержат две клетки поменьше (правую и левую или верхнюю и нижнюю), клетки-листья содержат точки. Это двумерное KD-дерево (от K-dimensional, т.е. K-мерное), обобщение двоичного дерева для многомерного пространства. Если размерность больше чем два, вместо прямых будут плоскости (гиперплоскости), перпендикулярные осям координат. Основная идея, надеюсь, понятна.

Третья мысль : а зачем нам клетки там, где нет точек, давайте строить клетки только там где надо. Описать так же кратко и понятно, как KD-дерево, не удалось, но примерно так:

Нет точек, нет и клеточек.

С появлением первой точки вокруг неё точно по размеру строится прямоугольник (со сторонами 0 на 0).

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

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

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

Узлы продолжают распадаться. Когда переполняется очередной узел, он разваливается на части и все они становятся детьми корня.

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

И так, пока все точки не добавим

Это называется R-дерево (R значит rectangle), обобщение B-дерева на многомерный случай. Если размерностей больше чем две, получаются n-мерные параллелепипеды. Возможно вы заметили, что при построении может получиться так, что разные узлы будут пересекаться между собой; это не критично, хотя и неприятно. И опять я надеюсь, что какая-то общая идея ясна.

Как искать

О том из каких соображений строить деревья ещё пара слов будет. Но после того, как всё построено, оказывается, что KD-дерево это просто частный случай R-дерева: каждый узел содержит только два узла следующего уровня, и в пространстве они расположены весьма специфично, но всё это в рамках допустимого для R-деревьев, а других отличий нет. Так что алгоритм для поиска ближайших точек можно написать общий. Примерно такой:

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

Задача #1: найти все точки, лежащие от x на расстоянии не больше чем R
Для внутренних узлов, дети это узлы:
def getBall(self, x, R): return sum(, )
Для листьев, дети это точки:
def getBall(self, x, R): return sum(, )
При желании, конечно, можно определить функцию getBall() для точек и стереть разницу между листьями и корнями в этом месте.

Задача #2: найти n ближайших к x точек
Для начала нужно заиметь очередь с приоритетами, упорядоченную так, что наиболее удалённая от x точка находится наверху, первая на выход. Стандартный питоновский heapq не позволяет указать свою функцию сравнения (или я не в курсе и можно как-то просто и красиво налету заменять оператор "меньше" для конкретных экземпляров, или разработчики не очень подумали над интерфейсом), видимо нужно что-то другое, скорее всего кто-то уже написал. Пусть такая штука есть, она передаётся в параметре q. Тогда примерно так:

Для внутренних узлов:
def getN(self, x, q, n): self.children.sort(key = lambda a, x=x: a.dist(x)) # какая разница, в каком порядке идут дети, отсортируем. for c in self.children: if (len(q) < n) or q.top.dist(x) > c.dist(x): # если есть ещё совсем незаполненные места, q = c.getN(x, q, n) # или есть шанс встретить более близкую точку, чем текущая худшая else: break # т.к. только что отсортировали детей, дальше будет только хуже return q
Для листьев:
def getN(self, x, q, n): self.children.sort(lambda a,x=x: a.dist(x)) # какая разница, в каком порядке идут дети, отсортируем for c in self.children: if (len(q) < n) or q.top.dist(x) > c.dist(x): # если есть ещё совсем незаполненные места, q.push(c) # или точка ближе, чем текущая худшая else: break return q
Идея с сортировкой детей мало того что сомнительна идеологически, так ещё и эффективность под вопросом. Если n невелико, возможно, быстрее будет нужное количество раз выполнить min или что-нибудь типа того.

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

Зачем ещё нужны

Вопрос, конечно, кто тут "ещё". У меня нет статистики и я не знаю, какое применение основное. Тем не менее.

3D графика, алгоритм трассировки лучей, позволяющий получать изображения фотографического качества с игрой света, отражениями, преломлениями и прочей дифракцией. Работает не очень быстро. Основная задача, решаемая во время работы алгоритма -- найти пересечение заданного луча с объектами сцены. Если объекты сложной формы (например, не квадратные), это не очень просто, и чем сложнее форма объекта, тем сложнее. Сцена достаточно велика, объект -- какая-то поверхность, триангулированная (например) и заданная вершинами треугольников. Треугольников много, вершин тоже. Точки распределены по объёму очень неравномерно, по идее большинство лучей проходят мимо, осталось понять, какие именно.

Построение вокруг объекта KD-дерева или что-то типа R-дерева (насколько я понял более-менее аналог, они называют это Bounding Volume Hierarchy, BVH) позволяет достаточно быстро понять, какой луч куда попадает.

Картинка, взята отсюда

Как строить

Увы, это отдельная проблема, с которой я не разбирался. Может быть когда-нибудь, но не в этот раз. KD-дерево, которое тупо делит клеточки пополам (перпендикулярно разным осям, по кругу), не задумываясь о том, где бы лучше провести плоскость, тоже работает. При таком подходе его можно строить не по набору точек, а добавляя по одной, так же как я описывал R-дерево.

К описанию построения R-дерева нужно добавить только алгоритм разделения узла на части. Тоже в лоб: делим на две части, в качестве центров кристаллизации выбираем двух наиболее удалённые друг от друга детей (точки или узлы следующего уровня, это O(n^2), где n -- количество детей в узле), собираем оставшиеся по принципу "какому прямоугольнику расширяться меньше". Этот алгоритм квадратичен относительно максимального количества детей в узле, результат не оптимален, но оптимальный требует, кажется, экспоненту. В общем, и так тоже работает, а по соотношению цена/качество может и вообще лучший.

Советую посмотреть на то, как эту задачу решают в 3D графике: KD-дерево , BVH . Но, возможно, у них другие критерии оптимальности, тут надо думать. Upd: Про BVH неприменимо, увы. Там не просто множество точек, там известен исходный объект из которого это множество получено и связи между точками. Поэтому они могут себе позволить строить оболочки вокруг объекта. Если же есть только точки, для уточнения позиций прямоугольников можно разве что применить что-нибудть типа алгоритмов кластеризации. Но скорее всего это окажется довольно долго.

А, да. В вики про KD-деревья описано нечто совсем другое . Они предлагают не различать листья и внутренние узлы, а каждой точке сопоставить свою плоскость. Это получается совсем честное обобщение бинарного дерева на многомерное пространство. Мне кажется, так не надо делать. Но на самом деле я не пробовал.

Ещё пара мыслей

Ровно две:

Понятно, что прямоугольники в R-дереве обусловлены только скоростью подсчёта попадания точки внутрь. То есть, ничто в алгоритме не мешает заменить их на любые другие фигуры. Если мне не изменяет способность делать выводы, это позволяет обобщить задачу до ситуации, когда у нас нет осей и координат, а есть только элементы множества и годная функция расстояния между ними. Годная в математическом смысле: неотрицательная, симметричная, с правилом треугольника. Например, строки и упоминавшееся уже расстояние Левенштейна. KD-дерево тут очевидно неприменимо, так как тут нету понятия плоскости. И параллелепипедов тоже нет. А вот шарики есть. И построить дерево из шариков вроде бы можно. По идее оно должно довольно шустро предлагать варианты исправления опечаток, интересно, делают такое или нет? Upd: на этом пути есть

Совсем не обязательно хранить в деревьях именно точки. Можно, например, шарики, или любые другие объекты. Главное чтобы была функция вычисления расстояния и чтобы объект целиком помещался в ячейку, т.е. чтобы расстояние до него можно было грубо оценить как расстояние до ячейки. Этого будет достаточно для реализации алгоритма поиска.

Математическое описание

K-мерное дерево - это несбалансированное дерево поиска для хранения точек из . Оно предлагает похожую на R-дерево возможность поиска в заданном диапазоне ключей. В ущерб простоте запросов, требования к памяти вместо .

Существуют однородные и неоднородные k-d деревья. У однородных k-d деревьев каждый узел хранит запись. При неоднородном варианте внутренние узлы содержат только ключи, листья содержат ссылки на записи.

В неоднородном k-d дереве при параллельно оси -мерной гиперплоскости в точке . Для корня нужно разделить точки через гиперплоскость на два по возможности одинаково больших множества точек и записать в корень, слева от этого сохраняются все точки у которых , справа те у которых . Для левого поддерева нужно разделить точки опять на новую «разделенную плоскость» , а сохраняется во внутреннем узле. Слева от этого сохраняются все точки у которых . Это продолжается рекурсивно над всеми пространствами. Потом всё начинается снова с первого пространства, до того пока каждую точку можно будет ясно идентифицировать через гиперплоскость.

K-d дерево можно построить за . Поиск диапазона может быть выполнен за , при этом обозначает размер ответа. Требование к памяти для самого дерева ограничено .

Операции с k -d деревьями

Структура

Структура дерева, описанная на языке C++ :

Const N= 10 ; // количество пространств ключей Struct Item{ // структура элемента int key[ N] ; // массив ключей определяющих элемент char * info; // информация элемента } ; Struct Node{ // структура узла дерева Item i; // элемент Node * left; // левое поддерево Node * right; // правое поддерево }

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

Анализ поиска элемента

Очевидно, что минимальное количество просмотренных элементов равно , а максимальное количество просмотренных элементов - , где - это высота дерева. Остаётся посчитать среднее количество просмотренных элементов .

Заданный элемент.

Рассмотрим случай . Найденными элементами могут быть:

и так для каждого пространства ключей. При этом средняя длина поиска в одном пространстве составляет:

Средняя величина считается по формуле:

Остаётся найти вероятность . Она равна , где - число случаев, когда , и - общее число случаев. Не сложно догадаться, что

Подставляем это в формулу для средней величины:

то есть, , где - высота дерева.

Если перейти от высоты дерева к количеству элементов, то:

Где - количество элементов в узле.

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

Добавление элементов

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

Алгоритм продвижения по дереву:

For (int i = 0 ; tree; i++ ) // i - это номер пространства if (tree- > x[ i] < tree- > t) // t - медиана tree = tree- > left; // переходим в левое поддерево else tree = tree- > right; // переходим в правое поддерево

Добавление выполняется за , где - высота дерева.

Удаление элементов

При удалении элементов дерева может возникнуть несколько ситуаций.

Удаление листа дерева - довольно простое удаление, когда удаляется один узел, и указатель узла-предка просто обнуляется.

Удаление узла дерева (не листа) - очень сложная процедура, при которой приходится перестраивать всё поддерево для данного узла.

Иногда процесс удаления узла решается модификациями k-d дерева. К примеру, если у нас в узле содержится массив элементов, то при удалении всего массива узел дерева остаётся, но новые элементы туда больше не записываются.

Поиск диапазона элементов

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

Алгоритм

Z – узел дерева [ (x_0_min, x_1_min, x_2_min,..., x_n_min) ,(x_0_max, x_1_max, x_2_max,..., x_n_max) ] - заданный диапазон Функция Array(Node * & Z) { If ([ x_0_min, x_1_min, x_2_min,..., x_n_min] < Z) { Z= Z- > left; // левое поддерево } else If ([ x_0_max, x_1_max, x_2_max,..., x_n_max] > Z) { Z= Z- > right; // правое поддерево } Else{ // просмотреть оба поддерева Array(Z- > right) ; // запустить функцию для правого поддерева Z= Z- > left; // просмотреть левое поддерево } }

Анализ

Очевидно, что минимальное количество просмотренных элементов это , где - высота дерева. Так же очевидно, что максимальное количество просмотренных элементов это , то есть просмотр всех элементов дерева. Остаётся посчитать среднее количество просмотренных элементов .

Заданный диапозон.

Оригинальная статья про k-d деревья даёт такую характеристику: для фиксированного диапазона.

Если перейти высоты дерева к количеству элементов, то это будет:

Поиск ближайшего соседа

Поиск ближайшего элемента состоит из 2-х задачек. Определение возможного ближайшего элемента и поиск ближайших элементов в заданном диапазоне, то есть

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

При этом при нахождении более близкого элемента корректируется радиус поиска.

Алгоритм

Z – корень дерева| List – список ближайших элементов| [ x_0,x_1,x_2...,x_n] - элемент для которого ищется ближайшие Len – минимальная длина Функция Maybe_Near(Node * & Z) { // поиск ближайшего возможного элемента While(Z) { // проверка элементов в узле for (i= 0 ; i< N; i++ ) { len_cur= sqrt ((x_0- x[ i] _0) ^ 2 + (x_1- x[ i] _1) ^ 2 + ... + (x_n- x[ i] _n) ^ 2 ) ; // длина текущего элемента if (Len> // установление новой длины Delete(List) ; // очистка списка Add(List) ; // добавить новый элемент в список If((x_0= x[ i] _0) && (x_1= x[ i] _1) && ... && (x_n= x[ i] _n) ) Return 1 ; } If([ x_0,x_1,x_2...,x_n] < Z) Z= Z- > left; // левое поддерево If([ x_0,x_1,x_2...,x_n] > Z) Z= Z- > right; // правое поддерево } } Функция Near(Node * & Z) { // поиск наиближайшего элемента в заданном диапазоне While(Z) { // проверка элементов в узле for (i= 0 ; i< N; i++ ) { len_cur= sqrt ((x_0- x[ i] _0) ^ 2 + (x_1- x[ i] _1) ^ 2 + ... + (x_n- x[ i] _n) ^ 2 ) ; // длина текущего элемента if (Len> длины текущего элемента) { Len= len_cur; // установление новой длины Delete(List) ; // очистка списка Add(List) ; // добавить новый элемент в список } Else if (длины равны) Add(List) ; // добавить новый элемент в список } If([ x_0,x_1,x_2...,x_n] + len> Z) { // если диапазон больше медианы Near(Z- > right) ; // просмотреть оба дерева Z= Z- > left; } If([ x_0,x_1,x_2...,x_n] < Z) Z= Z- > left; // левое поддерево If([ x_0,x_1,x_2...,x_n] > Z) Z= Z- > right; // правое поддерево } }

Анализ

Очевидно, что минимальное количество просмотренных элементов это , где h - это высота дерева. Так же очевидно, что максимальное количество просмотренных элементов это , то есть просмотр всех узлов. Остаётся посчитать среднее количество просмотренных элементов.

Заданный элемент, относительно которого нужно найти ближайший. Эта задачка разбивается на 2. Нахождения ближайшего элемента в узле и нахождения ближайшего элемента в заданном диапазоне. Для первой части потребуется один спуск по дереву, то есть .

Для второй части, как мы уже вычислили, поиск элементов в заданном диапазоне равен . Что бы узнать среднее достаточно просто сложить эти 2 величины:

См. также

Примечания

Ссылки

libkdtree++ , an open-source STL-like implementation of k -d trees in C++.
FLANN and its fork nanoflann , efficient C++ implementations of k -d tree algorithms.
kdtree A simple C library for working with KD-Trees
libANN Approximate Nearest Neighbour Library includes a k -d tree implementation
Caltech Large Scale Image Search Toolbox : a Matlab toolbox implementing randomized k -d tree for fast approximate nearest neighbour search, in addition to LSH, Hierarchical K-Means, and Inverted File search algorithms.
Heuristic Ray Shooting Algorithms , pp. 11 and after
Into contains open source implementations of exact and approximate (k)NN search methods using k -d trees in C++.

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

Введение

KD-Tree (K-мерное дерево), специальная "геометрическая" структура данных, которая позволяет разбить K-мерное пространство на "меньшие части", посредством сечения этого самого пространства гиперплоскостями(K > 3 ), плоскостями (K = 3 ), прямыми (K = 2 ) ну, и в случае одномерного пространства-точкой (выполняя поиск в таком дереве, получаем что-то похожее на binary search ).
Логично, что такое разбиение обычно используют для сужения диапазона поиска в K-мерном пространстве. Например, поиск ближнего объекта (вершины, сферы, треугольника и т.д.) к точке, проецирование точек на 3D сетку, трассировка лучей (активно используется в Ray Tracing) и т.п. При этом объекты пространства помещаются в специальные параллелепипеды-bounding box-ы (bounding box-ом назовем такой параллелепипед, который описывает исходное множество объектов или сам объект, если мы строим bounding box лишь для него. У точек в качестве bounding box-а берется bounding box с нулевой площадью поверхности и нулевым объемом), стороны которых параллельны осям координат.

Процесс деления узла

Итак, прежде чем использовать KD-Дерево, нужно его построить. Все объекты помещаются в один большой bounding box, описывающий объекты исходного множества (при этом каждый объект ограничен своим bounding box-ом), который затем разбивается (делится) плоскостью, параллельной одной из его сторон, на два. В дерево добавляются два новых узла. Таким же образом каждый из полученных параллелепипедов разбивается на два и т.д. Процесс завершается либо по специальному критерию (см. SAH ), либо при достижении определенной глубины дерева, либо же при достижении определенного числа элементов внутри узла дерева. Некоторые элементы могут одновременно входить как в левый, так и в правый узлы (например, когда в качестве элементов дерева рассматриваются треугольники).

Проиллюстрирую данный процесс на примере множества треугольников в 2D:

На рис.1 изображено исходное множество треугольников. Каждый треугольник помещен в свой bounding box, а все множество треугольников вписано в один большой корневой bounding box.

На рис.2 делим корневой узел на два узла (по OX): в левый узел попадают треугольники 1, 2, 5, а в правый-3, 4, 5.

На рис.3 , полученные узлы снова делятся на два, а треугольник 5 входит в каждый из них. Когда процесс заканчивается, получаем 4 листовых узла.

Огромное значение имеет выбор плоскости для разделения узла дерева. Существует огромное количество способов сделать это, я привожу лишь некоторые, наиболее часто встречаемые на практике способы(предполагается, что начальные объекты помещены в один большой bounding box, а разделение происходит параллельно одной из осей координат):

⦁ Способ 1 : Выбрать наибольшую сторону bounding box. Тогда секущая плоскость будет проходить через середину выбранной стороны.

⦁ Способ 2 : Рассекать по медиане: отсортируем все примитивы по одной из координат, а медианой назовем элемент (или центр элемента), который находится на средней позиции в отсортированном списке. Секущая плоскость будет проходить через медиану так, что количество элементов слева и справа будет примерно равным.

⦁ Способ 3 : Использование чередования сторон при разбиении: на глубине 0 бьем через середину стороны параллельной OX, следующий уровень глубины-через середину стороны параллельной OY, затем по OZ и т.д. Если мы "прошлись по всем осям", то начинаем процесс заново. Критерии выхода описаны выше.

⦁ Способ 4 : Наиболее "умный" вариант-использовать SAH (Surface Area Heuristic) функцию оценки разделения bounding box. (Об этом подробно будет рассказано ниже). SAH также предоставляет универсальный критерий остановки деления узла.

Способы 1 и 3 хороши, когда речь идет именно о скорости построения дерева (так как найти середину стороны и провести через нее плоскость, отсеивая элементы исходного множества "налево" и "направо", тривиально). В то же время, они довольно часто дают неплотное представление разбиения пространства, что может негативно повлиять на основные операции в KD-Дереве (особенно, при прослеживании луча в дереве).

Способ 2 позволяет достигнуть хороших результатов, но требует немалого количества времени, которое тратится на сортировку элементов узла. Как и способы 1, 3, он довольно прост в реализации.

Наибольший же интерес представляет способ с использованием "умной" SAH эвристики (способ 4).

Bounding box узла дерева делится N (параллельными осям) плоскостями на N + 1 часть (части, как правило, равны) по каждой из сторон (на самом деле, количество плоскостей для каждой стороны можно быть любым, но для простоты и эффективности берут константу). Далее, в возможных местах пересечения плоскости с bounding box-ом вычисляют значение специальной функции: SAH. Деление производится плоскостью с наименьшим значением SAH функции (на рисунке ниже, я предположил, что минимум достигается в SAH, следовательно, деление будет производиться именно этой плоскостью (хотя здесь 2D пространство-так что прямой)):

Значение SAH функции для текущей плоскости вычисляется следующим образом:

В своей реализации KD-Дерева я делю каждую сторону на 33 равные части, используя 32 плоскости. Таким образом, по результатам тестов мне удалось найти "золотую" середину-скорость работы основных функций дерева/скорость построения дерева.

В качестве SAH эвристики я использую функцию, представленную на рисунке выше.

Стоит отметить, что для принятия решения необходим лишь минимум данной функции по всем секущим плоскостям. Следовательно, используя простейшие математические свойства неравенств, а также отбросив умножение на 2 при расчете площади поверхности узла (в 3D)(SAR, SAL, SA ), данную формулу можно существенно упростить. В полной же мере расчеты производятся лишь один раз на узел: при оценке критерия выхода из функции деления. Столь простая оптимизация дает существенный прирост скорости построения дерева (x2.5 ).

Для эффективного вычисления значения SAH функции необходимо уметь быстро определить, сколько узловых элементов находятся справа от данной плоскости, а сколько-слева. Результаты будут неудовлетворительны, если в качестве алгоритма использовать грубый, так называемый brute force подход с квадратичной асимптотикой. Однако же при использовании "корзиночного" (binned) метода ситуация значительно меняется в лучшую сторону. Описание данного метода приведено ниже:

Предположим, что мы разбиваем сторону bounding box-а на N равных частей (количество плоскостей-(N-1)), bounding box храним парой координат(pointMin, pointMax-см. рис. 1 ), тогда за один проход по всем элементам узла мы можем для каждой плоскости точно определить, сколько элементов лежит справа, а сколько-слева от нее. Создадим два массива по N элементов в каждом, и проинициализируем нулями. Пусть это будут массивы с именами aHigh и aLow . Последовательно пробегаем по всем элементам узла. Для текущего элемента проверяем, в какую из частей попадают pointMin и pointMax его bounding box-а. Соответственно, получаем два индекса на элемент множества. Пусть это будут индексы с именами iHigh (для pointMax) и iLow (для pointMin). После этого выполним следующее: aHigh += 1, aLow +=1.

После прохода по всем элементам, получаем заполненные массивы aHigh и aLow. Для массива aHigh посчитаем частичные-постфикс (suffix) суммы, а для aLow посчитаем частичные-префикс (prefix) суммы (см. рисунок). Получается, что количество элементов справа (и только справа! ) от плоскости с индексом i будет равно aLow, количество элементов, лежащих слева (и только слева! ): aHigh[i], количество же элементов, которые входят как в левый, так и в правый узлы: AllElements-aLow-aHigh[i].

Поставленная задача решена, а иллюстрация этого нехитрого процесса приведена ниже:

Хотелось бы отметить, что получение заранее известного количества элементов слева и справа от "бьющей" плоскости позволяет нам заранее выделить нужное количество памяти под них (ведь далее, после получения минимума SAH, нам нужно еще раз пройтись по всем элементам и каждый разместить в нужном массиве, (а использование банального push_back(если reserve-не вызывался), приводит к постоянному аллоцированию памяти-весьма затратная операция), что положительно влияет на скорость работы алгоритма построения дерева (x3.3).

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

Перебирая константы cI и cT , можно добиться более плотной структуры дерева (или же наоборот), жертвуя временем работы алгоритма. Во многих статьях, посвященных в первую очередь построению KD-Дерева для Ray Tracing рендера, авторы используют значения cI = 1., cT = : чем больше значение cT, тем быстрее строится дерево, и тем хуже результаты прослеживания луча в таком дереве.

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

Критерий остановки деления узла, довольно прост:

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

Теперь, когда мы научились делить узел KD-Дерева, расскажу про исходные случаи, когда количество элементов в узле получается довольно большим, а критерии остановки по количеству элементов уводят алгоритм в бесконечность. Собственно, все внимание на картинку (на примере треугольников в 2D):

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

Процесс построения дерева

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

Дерево строится от корня. В каждом узле дерева я храню указатели на левое и правое поддеревья, если таковых у узла нет, то он считается листовым (листом т.е). Каждый узел хранит bounding box, который описывает объекты данного узла. В листовых(и только в листовых! ) узлах я храню индексы тех объектов, которые входят в данный узел. Однако, в процессе построения память под узлы выделяется порциями (т.е. сразу для K узлов: во-первых, так эффективнее работать с менеджером памяти, во-вторых,-подряд идущие элементы идеальны для кэширования). Такой подход запрещает хранить узлы дерева в векторе, ибо добавление новых элементов в вектор может приводить к реаллоцированию памяти под все существующие элементы в другое место.

Соответственно, указатели на поддеревья теряют всякий смысл. Я использую контейнер типа список (std::list), каждый элемент которого-вектор (std::vector), с заранее заданным размером (константа). Построение дерева выполняю многопоточно (использую Open MP), то есть каждое поддерево строю в отдельном потоке(x4 к скорости). Для копирования индексов в листовой узел идеально подходит использование move семантики (C++11)(+10% к скорости).

Поиск "ближнего" к точке в KD-Дереве

Итак, дерево построено, перейдем к описанию реализации операции поиска в KD-Дереве.

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

Решать поставленную задачу с применением bruteforce невыгодно: для множества из N точек и M треугольников ассимптотика составит O(N * M). Кроме того, хотелось бы, чтобы алгоритм вычислял расстояние от точки до треугольника "по-честному" и делал это быстро.

Воспользуемся KD-Деревом и выполним следующее:

⦁ Шаг 1 . Найдем самый ближний листовой узел к данной точке (в каждом узле, как известно, мы храним bounding box, а расстояние можем смело вычислять до центра((pointMax + pointMin)*0.5) bounding box-а узла) и обозначим его nearestNode.

⦁ Шаг 2 . Методом перебора среди всех элементов найденного узла (nearestNode). Полученное расстояние обозначим minDist.

⦁ Шаг 3 . Построим сферу с центром в исходной точке и радиусом длины minDist. Проверим, лежит ли эта сфера полностью внутри (то есть без каких либо пересечений сторон bounding box-узла) nearestNode. Если лежит, то ближний элемент найден, если нет, то перейдем к пункту 4.

⦁ Шаг 4 . Запустим от корня дерева, поиск ближнего элемента в радиусе: спускаясь вниз по дереву, проверяем, пересекает ли правый или левый узлы (кроме того, узел может лежать полностью внутри сферы или сфера внутри узла...) данная сфера. Если какой-либо узел был пересечен, то аналогичную проверку выполняем и для внутренних узлов этого же узла. Если мы пришли в листовой узел, то выполним переборный поиск ближнего в данном узле и сравним полученный результат с длиной радиуса сферы. Если радиус сферы больше найденного расстояния, обновим длину радиуса сферы вычисленным значением минимума. Дальнейшие спуски по дереву происходят с обновленной длиной радиуса сферы (если мы используем рекурсивный алгоритм, то радиус просто передается в функцию по ссылке, а далее, где необходимо, обновляется).

Понять суть описанного выше алгоритма помогает следующий рисунок:

По рисунку : предположим, мы нашли ближний листовой узел (голубой на рисунке) к данной точке (выделена красным цветом), тогда, выполнив поиск ближнего в узле, получим, что таковым является треугольник с ичёндексом 1, однако, как видно-это не так. Окружность с радиусом R пересекает смежный узел, следовательно, поиск нужно выполнить и в этом узле, а затем сравнить вновь найденный минимум с тем, что уже есть. Как итог, становится очевидным, что ближним является треугольник с индексом 2.

Теперь мне бы хотелось рассказать об эффективной реализации промежуточных операций, используемых в алгоритме поиска "ближнего".

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

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

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

Определить находится ли сфера с центром в точке O и радиусом R внутри конкретного узла, представленного bounding box, просто(3D):

Inline bool isSphereInBBox(const SBBox& bBox, const D3& point, const double& radius) { return (bBox.m_minBB < point - radius && bBox.m_maxBB > < point - radius && bBox.m_maxBB > point + radius && bBox.m_minBB < point - radius && bBox.m_maxBB > point + radius); }
С алгоритмом же определения пересечения сферы с bounding box-ом узла, нахождения узла внутри сферы или сферы внутри узла, дела обстоят несколько иначе. Снова проиллюстрирую (картинка взята с ) и приведу корректный код, позволяющий выполнить данную процедуру (в 2D, 3D-аналогично):

bool intersects(CircleType circle, RectType rect) { circleDistance.x = abs(circle.x - rect.x); circleDistance.y = abs(circle.y - rect.y); if (circleDistance.x > (rect.width/2 + circle.r)) { return false; } if (circleDistance.y > (rect.height/2 + circle.r)) { return false; } if (circleDistance.x <= (rect.width/2)) { return true; } if (circleDistance.y <= (rect.height/2)) { return true; } cornerDistance_sq = (circleDistance.x - rect.width/2)^2 + (circleDistance.y - rect.height/2)^2; return (cornerDistance_sq <= (circle.r^2)); }
Сперва (первая пара строк) мы сводим вычисления 4-ёх квадрантов к одному. В следующей паре строк проверяем, лежит ли окружность в зеленой области. Если лежит, то пересечения нет. Следующая пара строк проверит, входит ли окружность в оранжевую или серую области рисунка. Если входит, то пересечение обнаружено.

Собственно, этот расчет возвращает "false" для всех окружностей, центр которых находится в пределах красной области, и "true" для всех окружностей, центр которых находится в белой области.

В целом, данный код обеспечивает то, что нам и нужно (я привел здесь реализацию кода для 2D, однако в моем коде KD-Дерева я использую 3D вариант).

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

Как было сказано выше, "веерные" ситуации порождают большое число элементов внутри узла, чем их больше-тем медленнее поиск. Более того, если все элементы равноудалены от данной точки, то поиск будет осуществляться за O(N) (множество точек, которые лежат на поверхности сферы, а ближний ищем к центру сферы). Однако, если данные ситуации убрать, то поиск в среднем будет равносилен спуску по дереву с перебором элементов в нескольких узлах т.е. за O(log(N)) . Понятно, что скорость работы поиска зависит от числа посещенных листовых узлов дерева.

Рассмотрим следующие два рисунка:

Суть этих рисунков заключается в том, что если точка к которой мы ищем ближний элемент расположена очень-очень далеко от исходного bounding box-а множества, то сфера с радисом длины minDist(расстояние до ближнего), будет пересекать много больше узлов, чем если бы мы рассматривали эту же сферу, но с центром в точке гораздо ближней к исходному bounding box-у множества (естественно, minDist изменится). В общем, поиск ближних к сильно удаленной точке, осуществляется медленне, чем поиск к точке, расположенной близко к исходному множесту. Мои тесты эту информацию подтвердили.

Результаты и подведение итогов

В качестве итогов хотел бы добавить еще пару слов о своей реализации KD-Дерева и привести полученные результаты. Собственно, дизайн кода разрабатывался так, чтобы его можно было легко адаптировать под любые объекты исходного множества(треугольники, сферы, точки и т.п). Все что нужно-создать класс наследник, с переопределенными виртуальными функциями. Более того, моя реализация также предусматривает передачу специального класса Splitter-а , определяемого пользователем. Этот класс, а точнее его виртуальный метод split, определяют, где конкретно будет проходить секущая плоскость. В своей реализации я предоставляю класс, принимающий решение на основе SAH. Здесь, отмечу, что во многих статьях посвященных ускорению построения KD-Дерева на основе SAH, многие авторы на начальных значениях глубины дерева(вообще, когда количество элементов в узле дерева велико) используют более простые техники поиска секущей плоскости(вроде деления по центру или медиане), а SAH эвристику применяют лишь в тот момент, когда число элементов в узле небольшое.

Моя реализация таких ухищрений не содержит, но позволяет быстро добавить их(нужно лишь расширить новым параметром конструктор KD-Дерева и вызвать функцию-член построения с нужным сплитером, контролируя нужные ограничения). Поиск в дереве выполняю многопоточно. Все вычисления произвожу в числах с двойной точностью(double ). Максимальная глубина дерева задается константой(по дефолту-32). Определены некоторые #defines , позволяющие, например, выполнять обход дерева при поиске без рекурсии(с рекурсией, правда быстрее выходт-все узлы это элементы некоторого вектора(то есть расположены рядом по памяти), а значит, хорошо кэшируются). Вместе с кодом я предоставляю тестовые наборы данных(3D модели "измененного формата OFF" с разным количеством треугольников внутри(от 2 до 3 000 000)). Пользователь может скинуть дамп построенного дерева(формат DXF) и просмотреть его в соответствующей графической программе. Также программа ведет лог(можно включить/выключить) качества построения дерева: сбрасывается глубина дерева, максимальное количество элементов в листовом узле, среднее количесвто элементов в листовых узлах, время работы . Ни в коем случае не утверждаю что полученная реализация идеальна, да что уж там, я и сам знаю места где я промахнулся(например, не передаю аллокатор в параметр шаблона, частое присутствие С-кода(чтение и запись файлов выполняю не с помощью потоков), возможны незамеченные баги и т.п-it"s time to fix it). Ну и конечно, дерево сделано и оптимизировано строго для работы в 3D пространстве.

Тестирование кода я проводил на ноутбуке, со следующими характеристиками: Intel Core I7-4750HQ, 4 core(8 threads) 2 GHz, RAM-8gb, Win x64 app на Windows 10 . Коэффициенты для вычисления SAH я брал следующие: cT = 1., cI = 1.5 . И, если говорть о результатах, то получилось, что на 1, 5млн. треугольников дерево строится менее чем за 1,5 секунды. На 3-ех миллионах за 2,4 секунды. Для 1,5млн. точек и 1.5млн треугольников(точки расположены не очень далеко от исходной модели), поиск выполнился за 0.23 секунды, а если точки удалять от модели-время растет, вплоть до 3 секунд. Для 3-ех миллионов точек(снова, расположены близко к модели) и 3-ех миллионов треугольников поиск занимает около 0,7 секунд. Надеюсь, ничего не перепутал. Напоследок, иллюстрация визуализации построенного KD-Дерева:

). k -d деревья - особый вид двоичных деревьев поиска .

Математическое описание

K-мерное дерево - это несбалансированное дерево поиска для хранения точек из $\mathbb{R}^k$ . Оно предлагает похожую на R-дерево возможность поиска в заданном диапазоне ключей. В ущерб простоте запросов, требования к памяти $O(kn)$ вместо $O((log(n))^{k-1})$ .

Существуют однородные и неоднородные k-d деревья. У однородных k-d деревьев каждый узел хранит запись . При неоднородном варианте внутренние узлы содержат только ключи, листья содержат ссылки на записи.

В неоднородном k-d дереве $H_i(t) = (x_1, x_2, \ldots , x_{i-1}, t , x_{i+1}, \ldots , x_k)$ при $1 \leq i \leq k$ параллельно оси $(k-1)$ -мерной гиперплоскости в точке $t$ . Для корня нужно разделить точки через гиперплоскость $H_1(t)$ на два по возможности одинаково больших множества точек и записать $t$ в корень, слева от этого сохраняются все точки у которых $x_1, справа те у которых x_1>t . Для левого поддерева нужно разделить точки опять на новую «разделенную плоскость» H_2(t), а t сохраняется во внутреннем узле. Слева от этого сохраняются все точки у которых x_2 . Это продолжается рекурсивно над всеми пространствами. Потом всё начинается снова с первого пространства до того, пока каждую точку можно будет ясно идентифицировать через гиперплоскость.$

K-d дерево можно построить за $O(n(k+log(n)))$ . Поиск диапазона может быть выполнен за $O(n^{1-\frac{1}{k}}+a)$ , при этом $a$ обозначает размер ответа. Требование к памяти для самого дерева ограничено $O(kn)$ .