c++ - 为什么 KD 树对于点集中的最近邻搜索如此缓慢？

Question

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

[开始编辑]我现在使用“nanoflann”，它比 CGAL 中的等效项快约 100-1000 倍用于 K-neighbor 搜索。我使用“Intel Embree”进行光线转换，比 CGAL 的 AABB 树快 100-200 倍。[结束编辑]

我的任务如下所示:

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

CGAl 非常擅长它的工作，所以我有点怀疑他们的实现很糟糕。我用于光线追踪的他们的 AABB 树在几分钟内在地面上燃烧了十亿条光线......我想这很了不起。但另一方面，他们的 KD-tree 甚至无法处理 mio。在任何合理的时间点和 250k 样本(点查询)...

我想出了两个解决方案，可以摆脱 KD-trees:

1) 使用纹理贴图将光子存储在几何体的链表中。这始终是一个 O(1) 操作，因为无论如何您都必须进行光线转换...

2) 使用 View 相关的切片哈希集。那就是你离得越远，哈希集就越粗糙。所以基本上你可以在 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);
    }

}

Answer 1

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

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

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

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

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