c++ - ¿Por qué los árboles KD son tan lentos para la búsqueda del vecino más cercano en conjuntos de puntos?

data-structures raytracing

Question

Estoy usando la implementación del árbol KD de CGAL (la última) para buscar vecinos más cercanos en conjuntos de puntos. Y también Wikipedia y otros recursos parecen sugerir que los árboles KD son el camino a seguir. Pero de alguna manera son demasiado lentos y Wiki también sugiere su peor momento de O (n), lo cual está lejos de ser ideal.

[BEGIN-EDIT] Ahora estoy usando "nanoflann", que es aproximadamente 100-1000 veces más rápido que el equivalente en CGAL para la búsqueda de K-neighbor. Y estoy usando "Intel Embree" para raycasting, que es aproximadamente 100-200 veces más rápido que los árboles AABB de CGAL. [FIN-EDICION]

Mi tarea se ve así:

Tengo un conjunto de puntos ENORME, digamos como unos 100 mio. ¡¡puntos!! y su distribución está en superficies de geometría triangulada (sí, un trazador de fotones). Entonces, uno podría decir que su distribución es 2D en el espacio 3D, porque es escasa en 3D pero densa al mirar las superficies ... Este podría ser el problema, ¿verdad? Porque para mí esto parece desencadenar el peor de los casos de un árbol KD que probablemente podría funcionar mucho mejor con conjuntos de puntos densos en 3D ...

CGAl es bastante bueno en lo que hace, así que tengo un poco de duda de que su implementación simplemente apesta. Su árbol AABB que estoy usando para trazado de rayos quema mil millones de rayas lineales en unos pocos minutos en el suelo ... Eso es notable, supongo. Pero su árbol KD, por otro lado, ni siquiera puede tratar con un mio. puntos y 250k muestras (consultas puntuales) en cualquier tiempo razonable ...

Se me ocurrieron dos soluciones que sacan a patadas a los KD-trees:

1) Use mapas de textura para almacenar los fotones en una lista vinculada en la geometría. Esta es siempre una operación O (1), ya que tienes que hacer el raycast de todos modos ...

2) Use hashsets cortados dependientes de la vista. Cuanto más te alejas, más gruesas se vuelven las hashsets. Así que, básicamente, puede pensar en un ráster de 1024x1024x1024 en coordenadas NDC, pero con hashsets, para ahorrar memoria en áreas dispersas. Básicamente tiene una complejidad O (1) y se puede paralelizar de manera eficiente, tanto para inserciones (micro-sharding) como para consultas (sin bloqueos).

La primera solución tiene la desventaja de que es casi imposible promediar sobre las listas de fotones vecinas, lo cual es importante en las regiones más oscuras para evitar el ruido. La última solución no tiene este problema y debería estar a la par de las características de los árboles KD, solo que tiene O (1) peor rendimiento de caso, lol.

¿Entonces, qué piensas? ¿Mala implementación del árbol KD? Si no es así, ¿hay algo "mejor" que un árbol KD para las consultas delimitadas más cercanas? Quiero decir que no tengo nada en contra de mi segunda solución anterior, ¡pero una estructura de datos "comprobada" que ofrezca un rendimiento similar sería más agradable!

¡Gracias!

Aquí está el código (no compilable) que utilicé:

#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

Eche un vistazo a la biblioteca ApproxMVBB bajo la licencia MPL:

https://github.com/gabyx/ApproxMVBB :

La implementación de kdTree debería ser comparable a PCL (FLANN) y podría ser incluso más rápida. (¡Las pruebas con PCL parecían ser más rápidas con mi implementación!)

Diclaimer: soy el dueño de esta biblioteca y de ninguna manera esta biblioteca afirma ser más rápida y aún no se han realizado pruebas de rendimiento serio, pero estoy usando esta biblioteca con éxito en una dinámica granular rígida donde la velocidad es el rey. Sin embargo, esta biblioteca es muy pequeña, y la implementación de kdTree es muy genérica (vea los ejemplos) y le permite tener heurísticos de división personalizados y otras cosas sofisticadas :-).

Se implementan mejoras y consideraciones similares a las de nanoflann (acceso directo a datos, datos genéricos, n-dimensional) ... (ver el encabezado KdTree.hpp).

Algunas actualizaciones sobre el tiempo:
El ejemplo kdTreeFiltering contiene algunos puntos de referencia pequeños: se carga el conejito standord con 35947 puntos (ejemplo completamente funcional en el repositorio desde el primer momento):

Los resultados:

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)

Modelo Lucy.txt con 14 millones de puntos:

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)

¡Cuida la interpretación! y mira la configuración utilizada en el archivo de ejemplo. Sin embargo, comparando con los resultados de otras personas: ~ 3100ms para 14 * 10⁶ puntos es bastante hábil :-)

Procesador utilizado: CPU Intel® Core ™ i7 970 a 3.20 GHz × 12, 12 GB de RAM

Answer 2

Hace unos meses investigué un poco sobre las implementaciones rápidas del árbol KD, y estoy de acuerdo con Anony-Mousse en que la calidad (y el "peso" de las bibliotecas) varía mucho. Estos son algunos de mis hallazgos:

kdtree2 es una implementación de KD-tree poco conocida y bastante sencilla. Descubrí que era bastante rápido para problemas 3D, especialmente si permitía copiar y reordenar sus datos. Además, es pequeño y muy fácil de incorporar y adaptar. Aquí hay un documento sobre la implementación por parte del autor (no se deje intimidar por la mención de Fortran en el título). Esta es la biblioteca que terminé usando. Mis colegas compararon la velocidad de los puntos 3D con los árboles KD de VLFeat y otra biblioteca que no recuerdo (muchos FLANN, ver abajo) y ganó.

FLANN tiene la reputación de ser rápido, y se utiliza y recomienda con bastante frecuencia recientemente. Tiene como objetivo el caso dimensional más alto, donde ofrece algoritmos aproximados, pero también se utiliza en la Biblioteca de nubes de puntos que se ocupa de problemas en 3D.

No experimenté con CGAL porque encontré que la biblioteca era demasiado pesada. Estoy de acuerdo en que CGAL tiene una buena reputación, pero creo que ocasionalmente sufre demasiada sofisticación.

Answer 3

Las respuestas no son el lugar para hacer preguntas, pero su pregunta no es una pregunta, sino una afirmación que el árbol kd de CGAL apesta.

La lectura de 1,8 millones de puntos de un modelo de datos geológicos, y el cálculo de los 50 puntos más cerrados para cada uno de estos puntos tiene el siguiente rendimiento en mi Dell Precision, Windows 7, 64 bits, VC10:

• leyendo los puntos de un archivo: 10 segundos
• Construcción del árbol 1.3 sec
• 1.8 millones de consultas que informan los 50 puntos más cercanos: 17 segundos

¿Tienes actuaciones similares? ¿Esperarías que un árbol kd sea más rápido?

También me pregunto dónde están tus puntos de consulta, que está cerca de la superficie o cerca del esqueleto.

Answer 4

Según mi experiencia, la calidad de implementación varía ampliamente, lamentablemente. Sin embargo, nunca miré la implementación de CGAL.

El peor caso para el árbol kd generalmente es cuando debido a cambios incrementales se vuelve demasiado desequilibrado y debe ser recargado.

Sin embargo, en general, dichos árboles son más eficientes cuando no se conoce la distribución de datos.

En su caso, parece que un enfoque simple basado en la cuadrícula puede ser la mejor opción. Si lo desea, puede considerar que una textura es una cuadrícula 2d densa. Así que tal vez puedas encontrar una proyección 2d donde una cuadrícula funcione bien, y luego te cruces con esta proyección.