我正在使用 CGAL 的（最新的）KD-tree 实现来搜索点集中的最近邻。而且 Wikipedia 和其他资源似乎表明 KD-trees 是要走的路。但不知何故，它们太慢了，Wiki 还建议它们的最坏情况时间为 O(n)，这远非理想。
[BEGIN-EDIT] 我现在使用“nanoflann”，它比 CGAL 中的 K-neighbor 搜索快 100-1000 倍。我正在使用“Intel Embree”进行光线投射，它比 CGAL 的 AABB 树快大约 100-200 倍。 [结束编辑]
我的任务如下所示：
我有一个巨大的积分集，比如最多 100 mio。积分！！并且它们的分布在三角几何的表面上（是的，光子示踪剂）。所以可以说它们在 3D 空间中的分布是 2D，因为它在 3D 中是稀疏的，但在查看表面时是密集的……这可能是问题所在吗？因为对我来说，这似乎触发了 KD 树的最坏情况性能，它可能可以更好地处理 3D 密集点集......
CGAl 非常擅长它的工作，所以我有点怀疑他们的实现很糟糕。我用于光线追踪的他们的 AABB 树在几分钟内在地面上燃烧了十亿条光线......我想这很了不起。但另一方面，他们的 KD-tree 甚至无法处理 mio。在任何合理的时间点和 250k 样本（点查询）......
我想出了两个解决方案，它们可以摆脱 KD-trees：
1) 使用纹理贴图将光子存储在几何体的链表中。这始终是 O(1) 操作，因为无论如何您都必须进行光线投射......
2) 使用视图相关的切片哈希集。那就是你离得越远，哈希集就越粗糙。所以基本上你可以在 NDC 坐标中考虑一个 1024x1024x1024 的栅格，但是使用哈希集来节省稀疏区域的内存。这基本上具有 O(1) 复杂性并且可以有效地并行化，对于插入（微分片）和查询（无锁）都是如此。
前一种解决方案的缺点是几乎不可能对相邻的光子列表进行平均，这在较暗的区域避免噪声很重要。后一种解决方案没有这个问题，应该在功能方面与 KD-trees 相提并论，只是它具有 O(1) 最坏情况下的性能，哈哈。
所以你怎么看？糟糕的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);
}
}