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

puntos mas cercanos divide y venceras

Question

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?

Answer 1

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 objetivo 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 objetivo 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.

Answer 2

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.

Answer 3

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.

Answer 4

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 objetivo.

Answer 5

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 octree 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).

Answer 6

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.

Answer 7

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 <algorithm>
#include <iostream>
#include <vector>

#include <boost/lambda/lambda.hpp>  // _1
#include <boost/lambda/bind.hpp>    // bind()
#include <boost/tuple/tuple_io.hpp>

namespace {
  typedef double coord_t;
  typedef boost::tuple<coord_t,coord_t,coord_t> 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<point_t> 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<coord_t>(), // 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 <algorithm>           // make_heap
#include <functional>          // binary_function<>
#include <iostream>

#include <boost/range.hpp>     // boost::begin(), boost::end()
#include <boost/tr1/tuple.hpp> // get<>, tuple<>, cout <<

namespace {
  typedef double coord_t;
  typedef std::tr1::tuple<coord_t,coord_t,coord_t> 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 T>
  class less_distance : public std::binary_function<T, T, bool> {
    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<point_t> 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 <algorithm>           // sort
#include <functional>          // binary_function<>
#include <iostream>

#include <boost/foreach.hpp>
#include <boost/range.hpp>     // begin(), end()
#include <boost/tr1/tuple.hpp> // get<>, tuple<>, cout <<

#define foreach BOOST_FOREACH

namespace {
  typedef double coord_t;
  typedef std::tr1::tuple<coord_t,coord_t,coord_t> 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<point_t, point_t, bool> {
    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.

Answer 8

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 .

Answer 9

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.

Answer 10

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).

Answer 11

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.

Answer 12

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.