【问题标题】:Why are KD-trees so damn slow for nearest neighbor search in point sets?为什么 KD 树对于点集中的最近邻搜索如此缓慢?
【发布时间】:2013-02-14 00:23:18
【问题描述】:

我正在使用 CGAL 的(最新)KD-tree 实现来搜索点集中的最近邻。而且 Wikipedia 和其他资源似乎表明 KD-trees 是要走的路。但不知何故,它们太慢了,而且 Wiki 还建议它们的最坏情况时间为 O(n),这远非理想。

[开始编辑] 我现在使用“nanoflann”,它比 CGAL 中的 K 邻居搜索快约 100-1000 倍。我使用“Intel Embree”进行光线投射,比 CGAL 的 AABB 树快 100-200 倍。 [结束编辑]

我的任务如下所示:

我有一个巨大的积分集,比如最多 100 米奥。积分!!并且它们的分布在三角几何的表面上(是的,光子示踪剂)。因此可以说它们在 3D 空间中的分布是 2D,因为它在 3D 中是稀疏的,但在查看表面时是密集的……这可能是问题所在吗?因为在我看来,这似乎触发了 KD-tree 的最坏情况性能,它可能可以更好地处理 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);
    }

}

【问题讨论】:

  • 你能具体说明你是如何使用kd-tree的吗?空间排序包描述了许多选项...
  • 如果您知道不希望超过 50 个邻居,为什么不使用 Orthogonal_k_nearest_neighbor?
  • 我们在 1M 点和 50 个最近邻 (eps=0) 上做了一个简单的基准测试,性能相当。因此,请提供一个显示这种差异的基准。

标签: c++ data-structures raytracing kdtree cgal


【解决方案1】:

答案不是提问的地方,但你的问题不是问题,而是 CGAL 的 kd-tree 很烂的陈述。

读取地质数据模型的 1.8mio 点,并计算每个点的 50 个最近点,在我的 Dell Precision、Windows7、64bit、VC10 上具有以下性能:

  • 从文件中读取点:10 秒
  • 树的构造 1.3 秒
  • 1.8 个 mio 查询报告 50 个最近点:17 秒

你有类似的表现吗?你认为 kd-tree 会更快吗?

我还想知道您的查询点在哪里,靠近地表还是靠近骨架。

【讨论】:

  • 是的,这些性能结果(尽管它们不匹配但属于同一个“联盟”)就是我得到的。对我来说,它至少要慢 1000 倍(假设它线性扩展到 100 mio 点,它不会大声笑)......
  • 所以我需要的是 100 mio。点(这是我已经计算出的射线交点),2 mio。全高清查询(甚至可能是超级采样/抗锯齿的倍数),具有无限的邻居,但半径为 apx。场景宽度的 1/1000,大概需要 20 秒。从字面上看,其他一切都是浪费时间。
  • 好的。有人向我解释一下 mio 到底是什么。这是表示度量 M 的某种奇怪方式,即 1e61,000,000?请使用那个。人们会认真对待你。但我知道什么。这可能是当地的会议或其他什么。
  • 我不住在美国或任何英语国家,我从来不知道 mio。是百万的缩写。哪个国家使用这个?
  • @thesaint "Mio" 肯定不是美国百万的缩写(至少在中西部不是)。根据Wikipedia,“mio”用于“德国、瑞士和荷兰市场”。
【解决方案2】:

几个月前我对快速 KD-tree 实现进行了一些研究,我同意 Anony-Mousse 的观点,即质量(和库的“重量”)差异很大。以下是我的一些发现:

kdtree2 是一个鲜为人知且非常简单的 KD-tree 实现,我发现它对于 3D 问题非常快,尤其是如果您允许它复制和重新排序您的数据。此外,它很小,很容易合并和适应。 Here 是作者关于实现的论文(不要因为标题中提到 Fortran 而被推迟)。这是我最终使用的库。我的同事将它的 3D 点速度与VLFeat's KD-trees 和另一个我不记得的库(许多 FLANN,见下文)进行了基准测试,结果它赢了。

FLANN 以速度快着称,最近经常被使用和推荐。它针对更高维的情况,提供近似算法,但也用于处理 3D 问题的Point Cloud Library

我没有尝试 CGAL,因为我发现这个库太重了。我同意 CGAL 有很好的声誉,但我觉得它偶尔会被过度复杂化。

【讨论】:

  • 感谢您的概述!但我认为问题似乎在于数据结构,而不是它的实现。所以 KD-tree 本身似乎太慢了,不管你实现它有多好(通常你只会获得一个常数因子)。
【解决方案3】:

不幸的是,根据我的经验,实施质量差异很大。然而,我从来没有看过 CGAL 的实现。

k-d-tree 最坏的情况通常是由于增量更改而变得过于不平衡,应该重新加载。

但是,一般而言,当您不知道数据分布时,此类树的效率最高。

在您的情况下,听起来简单的基于网格的方法可能是最佳选择。如果需要,您可以将纹理视为密集的 2d 网格。所以也许你可以找到一个网格工作良好的二维投影,然后与这个投影相交。

【讨论】:

  • 你知道,如果没有一本关于光子追踪的书中推荐使用 KD-trees 来完成这项任务,我永远不会考虑使用 KD-trees。所以我假设他们会很合适,并想知道我是否做错了什么!
  • 是实用书还是理论书?另外,实施真的很重要。我已经看到 kd-trees 的性能非常好,但不幸的是,我看到它们对性能的贡献很小。
  • 这本书肯定有实际的目的,因为它附带了一个光子追踪器,并使用源代码来说明所有的东西。这就是我买它的原因;)。但当然,它也非常详细地涵盖了该理论。那么你知道一个好的 KD-tree 实现(或点集搜索实现,因为我不关心数据结构)吗??
  • 我实际上不太相信 KD 树适合您的用例。我不确定你为什么需要调用那个慢速采样方法(但我不喜欢光子追踪)。您不应该模拟光子射线撞击某些几何体(最终是您的相机),而不是管理庞大的光子数据结构并在那里进行最近邻查询吗?
  • “最终撞到你的相机”......不,那是行不通的。你追踪光子直到“最后一英里”。之后,您将来自相机的样本光线拍摄到场景中,收集可见光子。由于各种原因,您都需要它们。我不会在这里写书,但简而言之,您需要能够跟踪返回所有光源的路径,并且您需要每个光子中的所有方向和颜色信息,以根据高斯分布对每个屏幕点进行采样每个光子通过其特定的材质着色器(取决于它所击中的对象)对图像做出贡献
【解决方案4】:

看看 MPL 许可下的 ApproxMVBB 库:

https://github.com/gabyx/ApproxMVBB

kdTree 实现应该与 PCL(FLANN) 相当,并且可能更快。 (使用 PCL 进行测试似乎更快!)

免责声明:我是这个库的所有者,并且绝不是这个库声称更快,还没有进行严格的性能测试,但我在速度为王的颗粒刚体动力学中成功地使用了这个库! 然而,这个库非常小,并且 kdTree 实现非常通用(参见示例),并且允许您自定义拆分启发式和其他奇特的东西:-)。

实现了与 nanoflann(直接数据访问等,通用数据,n 维)类似的改进和考虑...(参见 KdTree.hpp)标头。

一些时间更新:
示例kdTreeFiltering包含一些小基准: 加载了 35947 点的标准兔子(开箱即用的 repo 中的完整工作示例):

结果:

兔子.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 1400 万点的模型:

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)

注意解释!并查看示例File 中使用的设置。 然而,与其他人的结果相比:14*10⁶ 点约 3100 毫秒是相当光滑的 :-)

使用的处理器:Intel® Core™ i7 CPU 970 @ 3.20GHz × 12 , 12GB Ram

【讨论】:

  • 能否给个测试,需要多长时间ms,点数,你的cpu,系统等等。
  • 这在很大程度上取决于您的应用程序的用例,如果您想要简单的中点拆分,而不是疯狂的拆分启发式优化,那么它比评估每个潜在拆分的最佳拆分技术要快得多,您可以编译项目中的示例(我添加了一些计时功能,Lucy模型也包含在子模块中以测试1400万分)我稍微更新了答案
【解决方案5】:

如果 kdtree 对于小集合来说很快,但对于大集合(>100000?)来说“慢”,你可能会因为刷新处理器缓存而受苦。如果前几个节点与很少使用的叶节点交错,那么您将在处理器缓存中放置较少的频繁使用的节点。这可以通过最小化节点的大小和仔细布局内存中的节点来改善。但最终你还是会刷新相当数量的节点。访问主内存可能会成为瓶颈。

Valgrind 告诉我,我的代码的一个版本减少了 5% 的指令,但我相信秒表告诉我对于相同的输入它会慢 10% 左右。我怀疑 valgrind 进行完整的缓存模拟会告诉我正确的答案。

如果您是多线程的,您可能希望鼓励线程在类似区域中进行搜索以重用缓存...假设单个多核处理器 - 多个处理器可能需要相反的方法。

提示:与 64 位指针相比,您在内存中打包的 32 位索引更多。

【讨论】:

    猜你喜欢
    • 2015-03-17
    • 2011-05-24
    • 2012-06-11
    • 2012-06-27
    • 1970-01-01
    • 2018-03-12
    • 2015-02-03
    • 2011-07-05
    • 1970-01-01
    相关资源
    最近更新 更多