为什么KD树在点集中最近邻搜索的速度如此之慢？

Question

我正在使用CGAL(最新的)KD树实现来搜索点集中的最近邻居.维基百科和其他资源似乎也暗示KD树是可行的方式.但不知怎的,它们太慢了,而且Wiki也暗示了O(n)的最坏情况,这远非理想.

[BEGIN-EDIT] 我现在正在使用"nanoflann",这比用于K邻居搜索的CGAL中的等效物快约100-1000倍.我正在使用"英特尔Embree"进行光线投射,这比CGAL的AABB树快约100-200倍. [END-EDIT]

我的任务看起来像这样:

我有一个巨大的点集,比如高达几百万.点!它们的分布在三角形几何体的表面上(是的,光子示踪剂).所以可以说它们的分布在3D空间中是2D的,因为它在3D中是稀疏的,但在观察表面时是密集的......这可能是问题对吗？因为对我而言,这似乎触发了KD树的最坏情况性能,这可能会更好地处理3D密集点...

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

我提出了两个解决方案,从KD树中剔除废话:

1)使用纹理贴图将光子存储在几何体上的链接列表中.这总是一个O(1)操作,因为你无论如何都必须进行光线投射......

2)使用视图相关的切片哈希集.这是你得到的越远,哈希集越粗.所以基本上你可以想到NDC坐标中的1024x1024x1024光栅,但是使用哈希集来节省稀疏区域的内存.这基本上具有O(1)复杂度,并且可以有效地并行化,用于插入(微分片)和查询(无锁).

前一解决方案的缺点在于,几乎不可能对相邻光子列表进行平均,这在较暗区域中是重要的,以避免噪声.后一种解决方案没有这个问题,并且应该与KD树特征相同,只是它具有O(1)最坏情况性能,lol.

所以你怎么看？糟糕的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树很糟糕的说法.

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

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

你有类似的表现吗？你期望kd树更快吗？

此外,我想知道您的查询点在哪里,接近表面,或接近骨架.

Answer 2

几个月前,我对快速KD树实现做了一些研究,我同意Anony-Mousse认为质量(和图书馆的"权重")变化很大.以下是我的一些发现:

kdtree2是一个鲜为人知且非常简单的KD树实现,我发现3D问题非常快,特别是如果你允许它复制和重新排序你的数据.此外,它很小,很容易合并和适应.这是一篇关于作者实施的论文(不要因为在标题中提到Fortran而推迟).这是我最终使用的库.我的同事们根据VLFeat的 KD树以及另一个我不记得的图书馆(许多FLANN,见下文)对3D点的速度进行了基准测试,并获胜.

FLANN的声誉很快,最近经常使用和推荐.它针对更高维度的情况,提供近似算法,但也用于处理3D问题的点云库.

我没有尝试CGAL,因为我发现图书馆太重了.我同意CGAL有良好的声誉,但我觉得它偶尔会受到过度批评的影响.

Answer 3

不幸的是，根据我的经验，实施质量差异很大。然而，我从未研究过 CGAL 的实现。

kd 树最糟糕的情况通常是由于增量更改而变得过于不平衡，并且应该重新加载。

但是，一般来说，当您不知道数据分布时，此类树是最有效的。

就您而言，听起来简单的基于网格的方法可能是最佳选择。如果需要，您可以将纹理视为密集的二维网格。因此，也许您可​​以找到一个网格效果良好的二维投影，然后与该投影相交。