Why are KD-trees so damn slow for nearest neighbor search in point sets?

Question

I am using CGAL's (the latest) KD-tree implementation for searching nearest neighbors in point sets. And also Wikipedia and other resources seem to suggest that KD-trees are the way to go. But somehow they are too slow and Wiki also suggests their worst-case time of O(n), which is far from ideal.

[BEGIN-EDIT] I am now using "nanoflann", which is about 100-1000 times faster than the equivalent in CGAL for K-neighbor search. And I am using "Intel Embree" for raycasting, which is about 100-200 times faster than CGAL's AABB trees. [END-EDIT]

My task looks like this:

I have a HUGE point set, say like up to a few 100 mio. points!! and their distribution is on surfaces of triangulated geometry (yes, a photon tracer). So one could say that their distribution is 2D in 3D space, because it is sparse in 3D but dense when looking at the surfaces... This could be the problem right? Because to me this seems to trigger the worst-case performance of a KD-tree which probably could deal much better with 3D dense point sets...

CGAl is quite good at what it does, so I have a bit doubt that their implementation just sucks. Their AABB tree I am using for raytracing burns a straight billion raytraces in a few mintues in the ground... That is remarkable I guess. But their KD-tree on the other hand can't even deal with a mio. points and 250k samples (point queries) in any reasonable time...

I came up with two solutions which kick the crap out of KD-trees:

1) Use texture maps to store the photons in a linked list on the geometry. This is always an O(1) operation, since you have to do the raycast anyway...

2) Use view dependent sliced hashsets. That is the farther away you get, the more coarse the hashsets get. So basically you can think of a 1024x1024x1024 raster in NDC coordinates, but with hashsets, to save memory in sparse areas. This basically has O(1) complexity and can be parallelized efficiently, both for inserts (micro-sharding) and queries (lock-free).

The former solution has the disadvantage that it is close to impossible to average over neighboring photon lists, which is important in darker regions to avoid noise. The latter solution doesn't have this problem and should be on par feature wise with KD-trees, just that it has O(1) worst case performance, lol.

So what do you think? Bad KD-tree implementation? If not, is there something "better" than a KD-tree for bounded nearest neighbor queries? I mean I have nothing against my second solution above, but a "proven" data-structure that delivers similar performance would be nicer!

Thanks!

Here is the code (not compilable though) that I used:

#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

From my experience, implementation quality varies widely, unfortunately. I have, however, never looked at the CGAL implementation.

The worst case for the k-d-tree usually is when due to incremental changes it becomes too unbalanced, and should be reloaded.

However, in general such trees are most efficient when you don't know the data distribution.

In your case it sounds as if a simple grid-based approach may be the best choice. If you want, you can consider a texture to be a dense 2d grid. So maybe you can find a 2d projection where a grid works good, and then intersect with this projection.