为什么kd树在点集的最近邻搜索中这么慢



我正在使用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位指针多。

相关内容

  • 没有找到相关文章

最新更新