我正在使用CGAL的(最新的)kd树实现来搜索点集中最近邻。维基百科和其他资源似乎也表明,kd树是可行的方法。但不知何故,他们太慢了,维基也建议他们的最坏情况时间为O(n),这远非理想。
[BEGIN-EDIT]我现在使用的是"nanoflann",对于k邻居搜索,它比等效的CGAL快100-1000倍。我使用"英特尔Embree"进行光线投射,它比CGAL的AABB树快100-200倍。(END-EDIT)
我的任务是这样的:
我有一个巨大的点集,比如说有几个亿。点! !它们的分布是在三角形几何的表面上(是的,光子示踪器)。所以可以说它们的分布在3D空间中是2D的,因为它在3D中是稀疏的,但在表面上是密集的……这可能就是问题所在,对吧?因为对我来说,这似乎触发了kd树的最坏性能,它可能会更好地处理3D密集点集…
CGAl做得很好,所以我有点怀疑他们的实现很糟糕。我用他们的AABB树进行光线追踪,在几分钟内在地面上直接燃烧了10亿条光线追踪。我想这很了不起。但另一方面,他们的kd树甚至不能处理一个mio。点和250k个样本(点查询)在任何合理的时间…
我想出了两个解决方案,把KD-trees踢出了屎:
1)使用纹理映射将光子存储在几何上的链表中。这总是一个0(1)的操作,因为你必须做光线投射…
2)使用视图相关的切片哈希集。也就是说你得到的越远,哈希集就越粗糙。所以基本上你可以用NDC坐标来考虑1024x1024x1024的栅格,但是用哈希集来节省稀疏区域的内存。这基本上具有0(1)复杂度,并且可以有效地并行化,无论是对于插入(微分片)还是查询(无锁)。
前者的缺点是几乎不可能对相邻的光子列表进行平均,这在较暗的区域对于避免噪声很重要。后一种解决方案没有这个问题,应该与kd树的功能相当,只是它有0(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);
}
}
答案不是提问的地方,但你的问题不是问题,而是CGAL的kd树糟糕的声明。
读取地质数据模型的180万个点,并为每个点计算50个最近的点,在我的Dell Precision, Windows7, 64位,VC10上的性能如下:
- 从文件中读取点:10秒
- 树的构建1.3秒
- 180万个查询报告最近的50个点:17秒
你们有类似的表演吗?你认为kd树会更快吗?
我也想知道你的查询点在哪里,即接近表面,或接近骨架。
几个月前我对快速kd树实现做了一些研究,我同意anonymous - mousse的观点,即库的质量(和"权重")差异很大。以下是我的一些发现:
kdtree2是一个鲜为人知的,非常直接的kd树实现,我发现对于3D问题来说非常快,特别是如果你允许它复制和重新排序你的数据。此外,它很小,很容易合并和适应。这是作者关于实现的一篇论文(不要因为在标题中提到Fortran而却步)。这是我最后使用的图书馆。我的同事将它的3D点速度与VLFeat的KD-trees和另一个我不记得的库(许多FLANN,见下文)进行了基准测试,结果它赢了。
FLANN以速度快而闻名,并且最近经常被使用和推荐。它针对高维情况,提供近似算法,但也用于处理3D问题的点云库。
我没有尝试过CGAL,因为我发现这个库太重量级了。我同意CGAL有良好的声誉,但我觉得它偶尔也有过于复杂的地方。
不幸的是,从我的经验来看,实现质量差异很大。但是,我从来没有看过CGAL的实现。
k-d树最坏的情况通常是由于增量更改,它变得太不平衡,应该重新加载。
然而,一般来说,当你不知道数据分布时,这样的树是最有效的。
在你的情况下,听起来好像一个简单的基于网格的方法可能是最好的选择。如果你愿意,你可以将纹理视为密集的2d网格。也许你可以找到一个二维的投影,在这个投影上用网格表示很好,然后与这个投影相交。
看看MPL许可下的ApproxMVBB库:
https://github.com/gabyx/ApproxMVBB:kdTree的实现应该与PCL(FLANN)相当,甚至可能更快。(使用PCL的测试似乎与我的实现更快!)
声明:我是这个库的所有者,绝不是这个库声称更快,严肃的性能测试还没有进行,但我在粒状刚体动力学中成功地使用了这个库,速度是王道!然而,这个库非常小,kdTree实现非常通用(参见示例),并允许您自定义拆分启发式和其他奇特的东西:-)。
实现了与nanoflann类似的改进和注意事项(直接数据访问等,通用数据,n维)…(见KdTree.hpp)头文件
关于时序的一些更新:
示例kdTreeFiltering
包含一些小基准:加载了35947个点的标准兔子(在现成的repo中有完整的工作示例):
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 model with 1400万points:
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)
注意解释!并查看示例文件中使用的设置。然而,与其他人的结果相比:~3100ms为14*10 26点是相当光滑的:-)
使用的处理器:Intel®Core™i7 CPU 970 @ 3.20GHz × 12, 12GB Ram
如果kdtree对于小集是快速的,但是对于大(>100000?)集"慢",您可能会刷新处理器缓存。如果最上面的几个节点与很少使用的叶节点交错,那么在处理器缓存中可以容纳更少的频繁使用的节点。这可以通过最小化节点的大小和在内存中仔细布局节点来改进。但最终还是会刷新相当数量的节点。最后,访问主存可能成为瓶颈。
Valgrind告诉我我的代码的一个版本是5%的指令少,但我相信秒表时,它告诉我它是大约10%慢相同的输入。我怀疑valgrind做一个完整的缓存模拟会告诉我正确的答案。
如果你是多线程的,你可能想鼓励线程在相似的区域做搜索来重用缓存…这假设一个多核处理器——多处理器可能需要相反的方法。
提示:内存中的32位索引比64位指针多。