Millones de puntos 3D: ¿cómo encontrar los 10 más cercanos a un punto dado?

Un punto en 3-d se define por (x, y, z). La distancia d entre cualesquiera dos puntos (X, Y, Z) y (x, y, z) es d = Sqrt [(Xx) ^ 2 + (Yy) ^ 2 + (Zz) ^ 2]. Ahora hay un millón de entradas en un archivo, cada entrada es un punto en el espacio, sin un orden específico. Dado cualquier punto (a, b, c), encuentre los 10 puntos más cercanos a él. ¿Cómo almacenaría el millón de puntos y cómo recuperaría esos 10 puntos de esa estructura de datos?

Millones de puntos es un número pequeño. El enfoque más directo funciona aquí (el código basado en KDTree es más lento (para consultar solo un punto)).

Enfoque de fuerza bruta (tiempo ~ 1 segundo)

#!/usr/bin/env python import numpy NDIM = 3 # number of dimensions # read points into array a = numpy.fromfile('million_3D_points.txt', sep=' ') a.shape = a.size / NDIM, NDIM point = numpy.random.uniform(0, 100, NDIM) # choose random point print 'point:', point d = ((a-point)**2).sum(axis=1) # compute distances ndx = d.argsort() # indirect sort # print 10 nearest points to the chosen one import pprint pprint.pprint(zip(a[ndx[:10]], d[ndx[:10]])) 

Ejecutarlo:

 $ time python nearest.py point: [ 69.06310224 2.23409409 50.41979143] [(array([ 69., 2., 50.]), 0.23500677815852947), (array([ 69., 2., 51.]), 0.39542392750839772), (array([ 69., 3., 50.]), 0.76681859086988302), (array([ 69., 3., 50.]), 0.76681859086988302), (array([ 69., 3., 51.]), 0.9272357402197513), (array([ 70., 2., 50.]), 1.1088022980015722), (array([ 70., 2., 51.]), 1.2692194473514404), (array([ 70., 2., 51.]), 1.2692194473514404), (array([ 70., 3., 51.]), 1.801031260062794), (array([ 69., 1., 51.]), 1.8636121147970444)] real 0m1.122s user 0m1.010s sys 0m0.120s 

Aquí está el script que genera millones de puntos 3D:

 #!/usr/bin/env python import random for _ in xrange(10**6): print ' '.join(str(random.randrange(100)) for _ in range(3)) 

Salida:

 $ head million_3D_points.txt 18 56 26 19 35 74 47 43 71 82 63 28 43 82 0 34 40 16 75 85 69 88 58 3 0 63 90 81 78 98 

Podría usar ese código para probar estructuras de datos y algoritmos más complejos (por ejemplo, si realmente consumen menos memoria o más rápido que el enfoque más simple anterior). Vale la pena señalar que en este momento es la única respuesta que contiene código de trabajo.

Solución basada en KDTree (tiempo ~ 1.4 segundos)

 #!/usr/bin/env python import numpy NDIM = 3 # number of dimensions # read points into array a = numpy.fromfile('million_3D_points.txt', sep=' ') a.shape = a.size / NDIM, NDIM point = [ 69.06310224, 2.23409409, 50.41979143] # use the same point as above print 'point:', point from scipy.spatial import KDTree # find 10 nearest points tree = KDTree(a, leafsize=a.shape[0]+1) distances, ndx = tree.query([point], k=10) # print 10 nearest points to the chosen one print a[ndx] 

Ejecutarlo:

 $ time python nearest_kdtree.py point: [69.063102240000006, 2.2340940900000001, 50.419791429999997] [[[ 69. 2. 50.] [ 69. 2. 51.] [ 69. 3. 50.] [ 69. 3. 50.] [ 69. 3. 51.] [ 70. 2. 50.] [ 70. 2. 51.] [ 70. 2. 51.] [ 70. 3. 51.] [ 69. 1. 51.]]] real 0m1.359s user 0m1.280s sys 0m0.080s 

Clasificación parcial en C ++ (tiempo ~ 1.1 segundos)

 // $ g++ nearest.cc && (time ./a.out < million_3D_points.txt ) #include  #include  #include  #include  // _1 #include  // bind() #include  namespace { typedef double coord_t; typedef boost::tuple point_t; coord_t distance_sq(const point_t& a, const point_t& b) { // or boost::geometry::distance coord_t x = a.get<0>() - b.get<0>(); coord_t y = a.get<1>() - b.get<1>(); coord_t z = a.get<2>() - b.get<2>(); return x*x + y*y + z*z; } } int main() { using namespace std; using namespace boost::lambda; // _1, _2, bind() // read array from stdin vector points; cin.exceptions(ios::badbit); // throw exception on bad input while(cin) { coord_t x,y,z; cin >> x >> y >> z; points.push_back(boost::make_tuple(x,y,z)); } // use point value from previous examples point_t point(69.06310224, 2.23409409, 50.41979143); cout << "point: " << point << endl; // 1.14s // find 10 nearest points using partial_sort() // Complexity: O(N)*log(m) comparisons (O(N)*log(N) worst case for the implementation) const size_t m = 10; partial_sort(points.begin(), points.begin() + m, points.end(), bind(less(), // compare by distance to the point bind(distance_sq, _1, point), bind(distance_sq, _2, point))); for_each(points.begin(), points.begin() + m, cout << _1 << "\n"); // 1.16s } 

Ejecutarlo:

 g++ -O3 nearest.cc && (time ./a.out < million_3D_points.txt ) point: (69.0631 2.23409 50.4198) (69 2 50) (69 2 51) (69 3 50) (69 3 50) (69 3 51) (70 2 50) (70 2 51) (70 2 51) (70 3 51) (69 1 51) real 0m1.152s user 0m1.140s sys 0m0.010s 

Cola de prioridad en C ++ (tiempo ~ 1.2 segundos)

 #include  // make_heap #include  // binary_function<> #include  #include  // boost::begin(), boost::end() #include  // get<>, tuple<>, cout << namespace { typedef double coord_t; typedef std::tr1::tuple point_t; // calculate distance (squared) between points `a` & `b` coord_t distance_sq(const point_t& a, const point_t& b) { // boost::geometry::distance() squared using std::tr1::get; coord_t x = get<0>(a) - get<0>(b); coord_t y = get<1>(a) - get<1>(b); coord_t z = get<2>(a) - get<2>(b); return x*x + y*y + z*z; } // read from input stream `in` to the point `point_out` std::istream& getpoint(std::istream& in, point_t& point_out) { using std::tr1::get; return (in >> get<0>(point_out) >> get<1>(point_out) >> get<2>(point_out)); } // Adaptable binary predicate that defines whether the first // argument is nearer than the second one to given reference point template class less_distance : public std::binary_function { const T& point; public: less_distance(const T& reference_point) : point(reference_point) {} bool operator () (const T& a, const T& b) const { return distance_sq(a, point) < distance_sq(b, point); } }; } int main() { using namespace std; // use point value from previous examples point_t point(69.06310224, 2.23409409, 50.41979143); cout << "point: " << point << endl; const size_t nneighbours = 10; // number of nearest neighbours to find point_t points[nneighbours+1]; // populate `points` for (size_t i = 0; getpoint(cin, points[i]) && i < nneighbours; ++i) ; less_distance less_distance_point(point); make_heap (boost::begin(points), boost::end(points), less_distance_point); // Complexity: O(N*log(m)) while(getpoint(cin, points[nneighbours])) { // add points[-1] to the heap; O(log(m)) push_heap(boost::begin(points), boost::end(points), less_distance_point); // remove (move to last position) the most distant from the // `point` point; O(log(m)) pop_heap (boost::begin(points), boost::end(points), less_distance_point); } // print results push_heap (boost::begin(points), boost::end(points), less_distance_point); // O(m*log(m)) sort_heap (boost::begin(points), boost::end(points), less_distance_point); for (size_t i = 0; i < nneighbours; ++i) { cout << points[i] << ' ' << distance_sq(points[i], point) << '\n'; } } 

Ejecutarlo:

 $ g++ -O3 nearest.cc && (time ./a.out < million_3D_points.txt ) point: (69.0631 2.23409 50.4198) (69 2 50) 0.235007 (69 2 51) 0.395424 (69 3 50) 0.766819 (69 3 50) 0.766819 (69 3 51) 0.927236 (70 2 50) 1.1088 (70 2 51) 1.26922 (70 2 51) 1.26922 (70 3 51) 1.80103 (69 1 51) 1.86361 real 0m1.174s user 0m1.180s sys 0m0.000s 

Enfoque basado en búsqueda lineal (tiempo ~ 1.15 segundos)

 // $ g++ -O3 nearest.cc && (time ./a.out < million_3D_points.txt ) #include  // sort #include  // binary_function<> #include  #include  #include  // begin(), end() #include  // get<>, tuple<>, cout << #define foreach BOOST_FOREACH namespace { typedef double coord_t; typedef std::tr1::tuple point_t; // calculate distance (squared) between points `a` & `b` coord_t distance_sq(const point_t& a, const point_t& b); // read from input stream `in` to the point `point_out` std::istream& getpoint(std::istream& in, point_t& point_out); // Adaptable binary predicate that defines whether the first // argument is nearer than the second one to given reference point class less_distance : public std::binary_function { const point_t& point; public: explicit less_distance(const point_t& reference_point) : point(reference_point) {} bool operator () (const point_t& a, const point_t& b) const { return distance_sq(a, point) < distance_sq(b, point); } }; } int main() { using namespace std; // use point value from previous examples point_t point(69.06310224, 2.23409409, 50.41979143); cout << "point: " << point << endl; less_distance nearer(point); const size_t nneighbours = 10; // number of nearest neighbours to find point_t points[nneighbours]; // populate `points` foreach (point_t& p, points) if (! getpoint(cin, p)) break; // Complexity: O(N*m) point_t current_point; while(cin) { getpoint(cin, current_point); //NOTE: `cin` fails after the last //point, so one can't lift it up to //the while condition // move to the last position the most distant from the // `point` point; O(m) foreach (point_t& p, points) if (nearer(current_point, p)) // found point that is nearer to the `point` //NOTE: could use insert (on sorted sequence) & break instead //of swap but in that case it might be better to use //heap-based algorithm altogether std::swap(current_point, p); } // print results; O(m*log(m)) sort(boost::begin(points), boost::end(points), nearer); foreach (point_t p, points) cout << p << ' ' << distance_sq(p, point) << '\n'; } namespace { coord_t distance_sq(const point_t& a, const point_t& b) { // boost::geometry::distance() squared using std::tr1::get; coord_t x = get<0>(a) - get<0>(b); coord_t y = get<1>(a) - get<1>(b); coord_t z = get<2>(a) - get<2>(b); return x*x + y*y + z*z; } std::istream& getpoint(std::istream& in, point_t& point_out) { using std::tr1::get; return (in >> get<0>(point_out) >> get<1>(point_out) >> get<2>(point_out)); } } 

Las mediciones muestran que la mayor parte del tiempo se usa para leer una matriz del archivo, los cálculos reales toman un orden de magnitud menor en tiempo.

Si el millón de entradas ya están en un archivo, no hay necesidad de cargarlas todas en una estructura de datos en la memoria. Simplemente mantenga una matriz con los diez mejores puntos encontrados hasta el momento, y escanee más de un millón de puntos, actualizando su lista de los diez primeros sobre la marcha.

Este es O (n) en el número de puntos.

Puede almacenar los puntos en un árbol k-dimensional (kd-tree). Los árboles de Kd están optimizados para búsquedas de vecinos más cercanos (encontrando los n puntos más cercanos a un punto dado).

Creo que esta es una pregunta difícil que prueba si no intentas exagerar las cosas.

Considere el algoritmo más simple que las personas ya han dado anteriormente: mantenga una tabla con los diez candidatos más lejanos y revise todos los puntos uno por uno. Si encuentra un punto más cercano que cualquiera de los diez mejores hasta ahora, reemplácelo. ¿Cuál es la complejidad? Bueno, tenemos que mirar cada punto del archivo una vez , calcular su distancia (o el cuadrado de la distancia en realidad) y compararlo con el décimo punto más cercano. Si es mejor, insértelo en el lugar apropiado en la tabla de las 10 mejores hasta ahora.

Entonces, ¿cuál es la complejidad? Vemos cada punto una vez, por lo que n son los cálculos de la distancia y n las comparaciones. Si el punto es mejor, tenemos que insertarlo en la posición correcta, esto requiere algunas comparaciones más, pero es un factor constante ya que la tabla de mejores candidatos es de tamaño constante 10.

Terminamos con un algoritmo que se ejecuta en tiempo lineal, O (n) en el número de puntos.

Pero ahora considere cuál es el límite inferior de dicho algoritmo. Si no hay ningún orden en los datos de entrada, tenemos que mirar cada punto para ver si no es uno de los más cercanos. Por lo que puedo ver, el límite inferior es Omega (n) y, por lo tanto, el algoritmo anterior es óptimo.

No es necesario calcular la distancia. Solo el cuadrado de la distancia debe servir a sus necesidades. Debería ser más rápido, creo. En otras palabras, puede omitir el bit sqrt .

Esta no es una pregunta para la tarea, ¿verdad? 😉

Mi pensamiento: iterar sobre todos los puntos y ponerlos en un montón mínimo o cola de prioridad limitada, marcada por la distancia desde el objective.

Esta pregunta esencialmente prueba tu conocimiento y / o intuición de los algoritmos de partición del espacio . Diría que almacenar los datos en un octárbol es tu mejor opción. Se usa comúnmente en motores 3D que manejan solo este tipo de problema (almacenando millones de vértices, trazado de rayos, encontrando colisiones, etc.). El tiempo de búsqueda estará en el orden de log(n) en el peor de los casos (creo).

Algoritmo directo:

Almacene los puntos como una lista de tuplas, escanee los puntos, calcule la distancia y mantenga una lista “más cercana”.

Mas creativo:

Agrupe puntos en regiones (como el cubo descrito por “0,0,0” a “50,50,50”, o “0,0,0” a “-20, -20, -20”), por lo que puede “indexar” en ellos desde el punto de destino. Compruebe en qué cubo se encuentra el objective y solo busque en los puntos de ese cubo. Si hay menos de 10 puntos en ese cubo, verifique los cubos “vecinos”, y así sucesivamente.

Pensándolo bien, este no es un algoritmo muy bueno: si tu punto objective está más cerca de la pared de un cubo que 10 puntos, entonces también tendrás que buscar en el cubo vecino.

Me gustaría ir con el enfoque kd-tree y encontrar el más cercano, luego eliminar (o marcar) ese nodo más cercano, y volver a buscar el nuevo nodo más cercano. Enjuague y repita.

Para cualquier dos puntos P1 (x1, y1, z1) y P2 (x2, y2, z2), si la distancia entre los puntos es d, entonces todo lo siguiente debe ser verdadero:

 |x1 - x2| <= d |y1 - y2| <= d |z1 - z2| <= d 

Mantenga el 10 más cercano mientras itera sobre su conjunto completo, pero también mantenga la distancia hasta el décimo más cercano. Ahórrese mucha complejidad usando estas tres condiciones antes de calcular la distancia para cada punto que mira.

básicamente una combinación de las dos primeras respuestas sobre mí. dado que los puntos están en un archivo, no es necesario mantenerlos en la memoria. En lugar de una matriz, o un montón mínimo, usaría un montón máximo, ya que solo quiere verificar distancias menores que el décimo punto más cercano. Para una matriz, necesitaría comparar cada distancia recién calculada con las 10 distancias que mantuvo. Para un montón mínimo, debe realizar 3 comparaciones con cada distancia recién calculada. Con un montón máximo, solo realiza 1 comparación cuando la distancia recién calculada es mayor que el nodo raíz.

Esta pregunta necesita una mayor definición.

1) La decisión con respecto a los algoritmos que pre-indexan los datos varía mucho dependiendo de si puede contener toda la información en la memoria o no.

Con kdtrees y octrees no tendrá que mantener los datos en la memoria y el rendimiento se beneficia de ese hecho, no solo porque la huella de memoria es menor sino simplemente porque no tendrá que leer todo el archivo.

Con fuerza bruta tendrás que leer todo el archivo y volver a calcular las distancias para cada nuevo punto que estarás buscando.

Aún así, esto podría no ser importante para ti.

2) Otro factor es cuántas veces tendrá que buscar un punto.

Como afirma JF Sebastian, a veces bruteforce es más rápido incluso en grandes conjuntos de datos, pero ten en cuenta que sus puntos de referencia miden la lectura de todo el conjunto de datos del disco (que no es necesario una vez que kd-tree o octree se construyen y escriben) y que miden solo una búsqueda.

Calcule la distancia para cada uno de ellos, y seleccione Seleccionar (1..10, n) en O (n) tiempo. Eso sería el ingenuo algoritmo, supongo.