Главная » Заработок » Анатомия KD-Деревьев. Отрывок, характеризующий K-мерное дерево

Анатомия KD-Деревьев. Отрывок, характеризующий K-мерное дерево

– структура данных, представляющая собой древовидную структуру в виде набора связанных узлов.

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

Способ представления бинарного дерева:

A — корень дерева
В — корень левого поддерева
С — корень правого поддерева

Корень дерева расположен на уровне с минимальным значением.

Узел D , который находится непосредственно под узлом B , называется потомком B . Если D находится на уровне i , то B – на уровне i-1 . Узел B называется предком D .

Максимальный уровень какого-либо элемента дерева называется его глубиной или высотой .

Если элемент не имеет потомков, он называется листом или терминальным узлом дерева.

Остальные элементы – внутренние узлы (узлы ветвления).

Число потомков внутреннего узла называется его степенью . Максимальная степень всех узлов есть степень дерева.

Число ветвей, которое нужно пройти от корня к узлу x , называется длиной пути к x . Корень имеет длину пути равную 0 ; узел на уровне i имеет длину пути равную i .

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

Имеется много задач, которые можно выполнять на дереве.

Распространенная задача - выполнение заданной операции p с каждым элементом дерева. Здесь p рассматривается как параметр более общей задачи посещения всех узлов или задачи обхода дерева.

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

Способы обхода дерева

Пусть имеем дерево, где A - корень, B и C - левое и правое поддеревья.

Существует три способа обхода дерева:

Обход дерева сверху вниз (в прямом порядке): A, B, C — префиксная форма.
Обход дерева в симметричном порядке (слева направо): B, A, C — инфиксная форма.
Обход дерева в обратном порядке (снизу вверх): B, C, A — постфиксная форма.

Реализация дерева

Узел дерева можно описать как структуру:

struct tnode {

Int field; // поле данных

Struct tnode *right; // правый потомок
};

При этом обход дерева в префиксной форме будет иметь вид

if (tree!=NULL) {

cout << tree->field; //Отображаем корень дерева

Обход дерева в инфиксной форме будет иметь вид

void treeprint(tnode *tree) {

if (tree!=NULL) { //Пока не встретится пустой узел

cout << tree->field; //Отображаем корень дерева

Обход дерева в постфиксной форме будет иметь вид

void treeprint(tnode *tree) {

if (tree!=NULL) { //Пока не встретится пустой узел

cout << tree->field; //Отображаем корень дерева

Бинарное (двоичное) дерево поиска – это бинарное дерево, для которого выполняются следующие дополнительные условия (свойства дерева поиска):

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

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

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

Для составления бинарного дерева поиска рассмотрим функцию добавления узла в дерево.

Добавление узлов в дерево

struct tnode * addnode(int x, tnode *tree) {

If (tree == NULL) { // Если дерева нет, то формируем корень

tree =new tnode; // память под узел

tree->field = x; // поле данных

tree->left = NULL;

tree->right = NULL; // ветви инициализируем пустотой

}else if (x < tree->field) // условие добавление левого потомка

tree->left = addnode(x,tree->left);

Else // условие добавление правого потомка

tree->right = addnode(x,tree->right);

Return (tree);
}

Удаление поддерева

void freemem(tnode *tree) {

If (tree!=NULL) {

freemem(tree->left);

freemem(tree->right);

Delete tree;

Пример Написать программу, подсчитывающую частоту встречаемости слов входного потока.

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

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

В дереве каждый узел содержит:

указатель на текст слова;
счетчик числа встречаемости;
указатель на левого потомка;
указатель на правого потомка.

Рассмотрим выполнение программы на примере фразы

При этом дерево будет иметь следующий вид

Пример реализации

#include
#include
#include
#include
#define MAXWORD 100
struct tnode { // узел дерева

Char *word; // указатель на строку (слово)

Int count; // число вхождений

Struct tnode *left; // левый потомок

Struct tnode *right; // правый потомок
};
// Функция добавления узла к дереву
struct tnode *addtree(struct tnode *p, char *w) {

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

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)$ .