Анатомия KD-Деревьев. KD-деревья и R-деревья Зачем ещё нужны

Эта статья полностью посвящена 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-Дерева:

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

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

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

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++.

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

Пусть у нас есть множество возможных элементов 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: на этом пути есть

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