Почему KD-деревья так чертовски медленны для поиска ближайшего соседа в наборах точек?

Я использую CGAL (новейшую) реализацию KD-дерева для поиска ближайших соседей в наборах точек. А также Википедия и другие ресурсы, похоже, предполагают, что KD-деревья - это путь. Но почему-то они слишком медленны, и Wiki также предлагает свое наихудшее время O (n), что далеко не идеально.

[НАЧАТЬ-EDIT] Теперь я использую "nanoflann", что примерно в 100-1000 раз быстрее, чем эквивалент в CGAL для поиска K-соседей. И я использую "Intel Embree" для raycasting, что примерно в 100-200 раз быстрее, чем деревья CGAL AABB. [END-EDIT]

Моя задача выглядит так:

У меня ОГРОМНЫЙ набор точек, скажем, до нескольких 100 миллионов. точки!! и их распределение на поверхностях триангулированной геометрии (да, фотонный трассировщик). Таким образом, можно сказать, что их распределение 2D в 3D-пространстве, потому что оно редкое в 3D, но плотное при взгляде на поверхности... Это может быть проблема правильно? Потому что для меня это, по-видимому, вызывает наихудшую производительность KD-дерева, которое, вероятно, может значительно улучшить работу с 3D плотными точками...

CGAl довольно хорош в том, что он делает, поэтому я немного сомневаюсь, что их реализация просто отстой. Их дерево AABB, которое я использую для raytracing, сжигает прямое миллиард raytraces в несколько минут в земле... Это замечательно, я думаю. Но их KD-дерево, с другой стороны, не может даже иметь дело с mio. точек и 250 тыс. выборок (точечных запросов) в любое разумное время...

Я придумал два решения, которые пинают дерьмо из KD-деревьев:

1) Используйте текстурные карты для хранения фотонов в связанном списке по геометрии. Это всегда операция O (1), так как вы все равно должны делать raycast...

2) Использовать зависящие от просмотра наборы хэш-наборов. Это тем дальше, что вы получаете, тем более грубый хэшсет. Таким образом, вы можете думать о растре 1024x1024x1024 в координатах NDC, но с хэш-настройками, чтобы сохранить память в разреженных областях. Это в основном имеет сложность O (1) и может эффективно распараллеливаться как для вставок (микросчет), так и для запросов (без блокировки).

Первое решение имеет тот недостаток, что почти невозможно усреднить по соседним спискам фотонов, что важно в более темных областях, чтобы избежать шума. Последнее решение не имеет этой проблемы и должно быть по парному признаку с KD-деревьями, так как оно имеет O (1) наихудшую производительность, lol.

И что ты думаешь? Реализация Bad KD-дерева? Если нет, есть ли что-то "лучше", чем KD-дерево для запросов ограниченного ближайшего соседа? Я имею в виду, что у меня нет ничего против моего второго решения выше, но "проверенная" структура данных, обеспечивающая схожую производительность, будет приятнее!

Спасибо!

Вот код (не компилируемый), который я использовал:

#include "stdafx.h"
#include "PhotonMap.h"

#pragma warning (push)
    #pragma warning (disable: 4512 4244 4061)
    #pragma warning (disable: 4706 4702 4512 4310 4267 4244 4917 4820 4710 4514 4365 4350 4640 4571 4127 4242 4350 4668 4626)
    #pragma warning (disable: 4625 4265 4371)

    #include <CGAL/Simple_cartesian.h>
    #include <CGAL/Orthogonal_incremental_neighbor_search.h>
    #include <CGAL/basic.h>
    #include <CGAL/Search_traits.h>
    #include <CGAL/point_generators_3.h>

#pragma warning (pop)

struct PhotonicPoint
{
    float vec[3];
    const Photon* photon;

    PhotonicPoint(const Photon& photon) : photon(&photon) 
    { 
        vec[0] = photon.hitPoint.getX();
        vec[1] = photon.hitPoint.getY();
        vec[2] = photon.hitPoint.getZ();
    }

    PhotonicPoint(Vector3 pos) : photon(nullptr) 
    { 
        vec[0] = pos.getX();
        vec[1] = pos.getY();
        vec[2] = pos.getZ();
    }

    PhotonicPoint() : photon(nullptr) { vec[0] = vec[1] = vec[2] = 0; }

    float x() const { return vec[0]; }
    float y() const { return vec[1]; }
    float z() const { return vec[2]; }
    float& x() { return vec[0]; }
    float& y() { return vec[1]; }
    float& z() { return vec[2]; }

    bool operator==(const PhotonicPoint& p) const
    {
        return (x() == p.x()) && (y() == p.y()) && (z() == p.z()) ;
    }

    bool operator!=(const PhotonicPoint& p) const 
    { 
        return ! (*this == p); 
    }
}; 

namespace CGAL 
{
    template <>
    struct Kernel_traits<PhotonicPoint> 
    {
        struct Kernel 
        {
            typedef float FT;
            typedef float RT;
        };
    };
}

struct Construct_coord_iterator
{
    typedef const float* result_type;

    const float* operator()(const PhotonicPoint& p) const
    { 
        return static_cast<const float*>(p.vec); 
    }

    const float* operator()(const PhotonicPoint& p, int) const
    { 
        return static_cast<const float*>(p.vec+3); 
    }
};

typedef CGAL::Search_traits<float, PhotonicPoint, const float*, Construct_coord_iterator> Traits;
typedef CGAL::Orthogonal_incremental_neighbor_search<Traits> NN_incremental_search;
typedef NN_incremental_search::iterator NN_iterator;
typedef NN_incremental_search::Tree Tree;


struct PhotonMap_Impl
{
    Tree tree;

    PhotonMap_Impl(const PhotonAllocator& allocator) : tree()
    {
        int counter = 0, maxCount = allocator.GetAllocationCounter();
        for(auto& list : allocator.GetPhotonLists())
        {
            int listLength = std::min((int)list.size(), maxCount - counter);
            counter += listLength; 
            tree.insert(std::begin(list), std::begin(list) + listLength);
        }

        tree.build();
    }
};

PhotonMap::PhotonMap(const PhotonAllocator& allocator)
{
    impl = std::make_shared<PhotonMap_Impl>(allocator);
}

void PhotonMap::Sample(Vector3 where, float radius, int minCount, std::vector<const Photon*>& outPhotons)
{
    NN_incremental_search search(impl->tree, PhotonicPoint(where));
    int count = 0;

    for(auto& p : search)
    {
        if((p.second > radius) && (count > minCount) || (count > 50))
            break;

        count++;
        outPhotons.push_back(p.first.photon);
    }

}

Ответы

Ответ 1

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

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

Однако, как правило, такие деревья наиболее эффективны, когда вы не знаете распределения данных.

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

Ответ 2

Ответы не могут задавать вопросы, но ваш вопрос не вопрос, а утверждение о том, что kd-дерево CGAL сосет.

Чтение 1,8 млн. точек геологической модели данных и вычисление 50 наиболее близких точек для каждой из этих точек имеет следующую производительность на моей Dell Precision, Windows7, 64bit, VC10:

чтение точек из файла: 10 секунд
Построение дерева 1.3 сек
1,8 млн. запросов, сообщающих о 50 ближайших точках: 17 секунд

Есть ли у вас подобные выступления. Ожидаете ли вы, что kd-дерево будет быстрее?

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

Ответ 3

Несколько месяцев назад я провел некоторое исследование быстрых реализаций KD-дерева, и я согласен с Anony-Mousse в том, что качество (и "вес" библиотек) сильно меняется. Вот некоторые из моих выводов:

kdtree2 - это немного известная и довольно простая реализация KD-дерева, которую я нашел довольно быстрым для 3D-проблем, особенно если вы разрешите это для копирования и повторной сортировки ваших данных. Кроме того, он очень мал и очень легко внедряется и адаптируется. Здесь - документ об исполнении автором (не отвлекайтесь на упоминание Fortran в названии). Это библиотека, которую я использовал. Мои коллеги сравнили скорость 3D-очков с VLFeat KD-деревья и другую библиотеку, которую я не помню (многие FLANN, см. Ниже), и она выиграла.

FLANN имеет репутацию быстрого и часто используется и рекомендуется совсем недавно. Он нацелен на более крупный размерный случай, где он предлагает приблизительные алгоритмы, но также используется в Cloud Cloud Library, которая занимается проблемами 3D.

Я не экспериментировал с CGAL, так как обнаружил, что библиотека слишком тяжелая. Я согласен с тем, что CGAL имеет хорошую репутацию, но я чувствую, что он иногда страдает от переосмысления.

Ответ 4

Взгляните на библиотеку ApproxMVBB под лицензией MPL:

https://github.com/gabyx/ApproxMVBB:

Реализация kdTree должна быть сопоставима с PCL (FLANN) и может быть еще быстрее. (тесты с PCL, по-видимому, были быстрее с моей реализацией!)

Diclaimer: Я являюсь владельцем этой библиотеки, и ни в коем случае эта библиотека не заявляет о каких-либо более быстрых и серьезных тестах производительности еще не проводилась, но я использую эту библиотеку успешно в гранулированной динамике твердого тела, где скорость короля! Тем не менее, эта библиотека очень мала, и реализация kdTree очень универсальна (см. Примеры) и позволяет создавать пользовательские разделительные эвристики и другие причудливые вещи: -).

Аналогичные улучшения и соображения, как в нанофланне (прямой доступ к данным и т.д., общие данные, n-мерные), реализованы... (см. заголовок KdTree.hpp).

Некоторые обновления по времени:
Пример kdTreeFiltering содержит некоторые небольшие контрольные показатели: Загружается лондонский кролик с 35947 очками (полностью рабочий пример в репо из коробки):

Результаты:

Bunny.txt

Loaded: 35947 points 
KDTree::  Exotic point traits , Vector3* +  id, start: =====
KdTree build took: 3.1685ms.
Tree Stats: 
     nodes      : 1199
     leafs      : 600
     tree level : 11
     avg. leaf data size : 29.9808
     min. leaf data size : 0
     max. leaf data size : 261
     min. leaf extent    : 0.00964587
     max. leaf extent    : 0.060337
SplitHeuristics Stats: 
     splits            : 599
     avg. split ratio  (0,0.5] : 0.5
     avg. point ratio  [0,0.5] : 0.22947
     avg. extent ratio (0,1]   : 0.616848
     tries / calls     : 599/716 = 0.836592
Neighbour Stats (if computed): 

     min. leaf neighbours    : 6
     max. leaf neighbours    : 69
     avg. leaf neighbours    : 18.7867
(Built with methods: midpoint, no split heuristic optimization loop)


Saving KdTree XML to: KdTreeResults.xml
KDTree:: Simple point traits , Vector3 only , start: =====
KdTree build took: 18.3371ms.
Tree Stats: 
     nodes      : 1199
     leafs      : 600
     tree level : 10
     avg. leaf data size : 29.9808
     min. leaf data size : 0
     max. leaf data size : 306
     min. leaf extent    : 0.01
     max. leaf extent    : 0.076794
SplitHeuristics Stats: 
     splits            : 599
     avg. split ratio  (0,0.5] : 0.448302
     avg. point ratio  [0,0.5] : 0.268614
     avg. extent ratio (0,1]   : 0.502048
     tries / calls     : 3312/816 = 4.05882
Neighbour Stats (if computed): 

     min. leaf neighbours    : 6
     max. leaf neighbours    : 43
     avg. leaf neighbours    : 21.11
(Built with methods: midpoint, median,geometric mean, full split heuristic optimization)

Модель Lucy.txt с 14 миллионами точек:

Loaded: 14027872 points 
KDTree::  Exotic point traits , Vector3* +  id, start: =====
KdTree build took: 3123.85ms.
Tree Stats: 
         nodes      : 999999
         leafs      : 500000
         tree level : 25
         avg. leaf data size : 14.0279
         min. leaf data size : 0
         max. leaf data size : 159
         min. leaf extent    : 2.08504
         max. leaf extent    : 399.26
SplitHeuristics Stats: 
         splits            : 499999
         avg. split ratio  (0,0.5] : 0.5
         avg. point ratio  [0,0.5] : 0.194764
         avg. extent ratio (0,1]   : 0.649163
         tries / calls     : 499999/636416 = 0.785648
(Built with methods: midpoint, no split heuristic optimization loop)

KDTree:: Simple point traits , Vector3 only , start: =====
KdTree build took: 7766.79ms.
Tree Stats: 
     nodes      : 1199
     leafs      : 600
     tree level : 10
     avg. leaf data size : 11699.6
     min. leaf data size : 0
     max. leaf data size : 35534
     min. leaf extent    : 9.87306
     max. leaf extent    : 413.195
SplitHeuristics Stats: 
     splits            : 599
     avg. split ratio  (0,0.5] : 0.297657
     avg. point ratio  [0,0.5] : 0.492414
     avg. extent ratio (0,1]   : 0.312965
     tries / calls     : 5391/600 = 8.985
Neighbour Stats (if computed): 

     min. leaf neighbours    : 4
     max. leaf neighbours    : 37
     avg. leaf neighbours    : 12.9233
(Built with methods: midpoint, median,geometric mean, full split heuristic optimization)

Позаботьтесь об интерпретации! и посмотрите настройки, используемые в примере Файл. Однако по сравнению с результатами других людей: ~ 3100 мс для 14 * 10⁶ точек довольно гладкий: -)

Используемый процессор: Intel® Core ™ i7 CPU 970 @3.20GHz × 12, 12GB Ram