Я использую 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);
}
}
