Анатомія KD-Дерев

image
Ця стаття повністю присвячена KD-Деревам: я описую тонкощі побудови KD-Дерев, тонкощі реалізації функцій пошуку 'близького' в KD-Дереві, а також можливі 'підводні камені", які виникають в процесі вирішення тих чи інших завдань алгоритму. Щоб не заплутувати читача термінологією(площина, гіпер-площина і т. п), та й взагалі для зручності, що належить основне дійство розгортається в тривимірному просторі. Однак же, де потрібно, я наголошую, що ми працюємо в просторі іншої розмірності. На мою думку стаття буде корисна як програмістам, так і всім тим, хто зацікавлений у вивченні алгоритмів: хтось знайде для себе щось нове, а хтось просто повторить матеріал і можливо, в коментарях доповнить статтю. У будь-якому випадку прошу всіх під кат.

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

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

Проілюструю даний процес на прикладі безлічі трикутників в 2D:

image
рис.1 зображено початкове безліч трикутників. Кожен трикутник поміщений в свій bounding box, а всі безліч трикутників вписано в один великий кореневої bounding box.

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

image
рис.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[7], отже, поділ буде проводитися саме цією площиною (хоча тут 2D простір-так що прямий)):

image

Значення SAH функції для поточної площині обчислюється наступним чином:

image

У своїй реалізації 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[iHigh] += 1, aLow[iLow] +=1.

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

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

image

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

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

Перебираючи константи cI cT, можна домогтися більш щільної структури дерева (або ж навпаки), жертвуючи часом роботи алгоритму. У багатьох статтях, присвячених в першу чергу побудови KD-Дерева для Ray Tracing рендера, автори використовують значення cI = 1., cT = [1; 2]: чим більше значення cT, тим швидше будується дерево, і тим гірші результати дослідження променя в такому дереві.

У своїй же реалізації я використовую дерево для пошуку 'близького', і, приділивши належну увагу отриманими результатами тестування і пошуку потрібних коефіцієнтів, можна помітити, що високі значення cT-дають нам абсолютно не щільно заповнені елементами вузли. Щоб уникнути подібної ситуації, мною було вирішено виставити значення cT 1., а значення cI протестувати на різних-великих-одиниці даних. В результаті-вдалося отримати досить щільну структуру дерева, розплатившись значним зростанням часу при побудові. На результатах ж пошуку 'близького' це дія відбилося позитивно-швидкість пошуку збільшилася.

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

image

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

Тепер, коли ми навчилися ділити вузол KD-Дерева, розповім про вихідні випадки, коли кількість елементів у вузлі виходить досить великим, а критерії зупинки за кількістю елементів відводять алгоритм у нескінченність. Власне, все увагу на картинку (на прикладі трикутників в 2D):

image

Я називаю такі ситуації 'віяловими' (мають спільну точку дотику, збігаються об'єкти, також відношу цю категорію). Видно, що як би ми не проводили січну площину, центральна точка так чи інакше потрапляє в один з вузлів, а разом з нею і трикутники, для яких вона є загальною.

Процес побудови дерева
Ми навчилися ділити вузол дерева, тепер залишилося застосувати отримані знання до процесу побудови всього дерева. Нижче я привожу опис своєї реалізації побудови 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. Запустимо від кореня дерева, пошук близького елемента в радіусі: спускаючись вниз по дереву, перевіряємо, перетинає чи правий або лівий вузли (крім того, вузол може лежати повністю всередині сфери або сфера всередині вузла...) дана сфера. Якщо який-небудь вузол був пересічений, то аналогічну перевірку виконуємо і для внутрішніх вузлів цього ж сайту. Якщо ми прийшли в листовий вузол, то виконаємо переборный пошук ближнього в даному вузлі і порівняємо отриманий результат з довжиною радіусу сфери. Якщо радіус сфери більше знайденого відстані, оновимо довжину радіуса сфери обчисленим значенням мінімуму. Подальші спуски по дереву відбуваються з оновленою довжиною радіуса сфери (якщо ми використовуємо рекурсивний алгоритм, то радіус просто передається у функцію за посиланням, а далі, де необхідно, оновлюється).

Зрозуміти суть описаного вище алгоритму допомагає наступний малюнок:

image

По малюнку: припустимо, ми знайшли близький листовий вузол (блакитний на малюнку) до даної точки (виділена червоним кольором), тоді, виконавши пошук ближнього у вузлі, отримаємо, що таким є трикутник з ичендексом 1, однак, як видно-це не так. Окружність з радіусом R перетинає суміжний вузол, отже, пошук потрібно виконати і в цьому вузлі, а потім порівняти знову знайдений мінімум з тим, що вже є. Як підсумок, стає очевидним, що ближнім є трикутник з індексом 2.

Тепер мені б хотілося розповісти про ефективної реалізації проміжних операцій, що використовуються в алгоритмі пошуку 'близького'.

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

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

Описаний алгоритм є не найефективнішим. Більш ефективний підхід спирається на пошук та аналіз (знаходження мінімуму градієнта тощо) функції, значеннями якої є всі можливі відстані від даної точки до будь-якої точки в трикутнику. Метод досить складний у сприйнятті і, як я вважаю, заслуговує окремої статті (поки що він реалізований у мене в коді, а посилання на код ви знайдете нижче). Ознайомитися з методом можна з [1]. Даний метод за результатами тестів опинився в 10 разів швидше того, що я описував раніше.

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

inline bool isSphereInBBox(const SBBox& bBox const D3& point, const double& radius)
{
return (bBox.m_minBB[0] < point[0] - radius && bBox.m_maxBB[0] > point[0] + radius &&
bBox.m_minBB[1] < point[1] - radius && bBox.m_maxBB[1] > point[1] + radius &&
bBox.m_minBB[2] < point[2] - radius && bBox.m_maxBB[2] > point[2] + radius);
}

З алгоритмом визначення перетину сфери з bounding box-го вузла, знаходження вузла всередині сфери або сфери усередині вузла, справи йдуть трохи інакше. Знову проілюструю (картинка взята з [2]) і наведу коректний код, що дозволяє виконати дану процедуру (в 2D, 3D-аналогічно):

image
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-Дерева:

image

Корисні посилання
[0]: Моя реалізація KD-Tree на GitHub
[1]: Пошук відстані від точки до трикутника
[2]: Визначення перетинів кола та прямокутники
[3]: Опис KD-Tree на Wikipedia
[4]: Цікава стаття присвячена SAH
[5]: Опис KD-Tree на ray-tracing.ru
Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.