Milhões de pontos 3D: como encontrar os 10 deles mais próximos de um determinado ponto?

Um ponto em 3-d é definido por (x, y, z). Distância d entre quaisquer dois pontos (X, Y, Z) e (x, y, z) é d = Sqrt [(Xx) ^ 2 + (Yy) ^ 2 + (Zz) ^ 2]. Agora há um milhão de inputs em um arquivo, cada input é algum ponto no espaço, em nenhuma ordem específica. Dado qualquer ponto (a, b, c), encontre os 10 pontos mais próximos. Como você armazenaria o milhão de pontos e como você recuperaria esses 10 pontos dessa estrutura de dados?

Milhões de pontos é um pequeno número. A abordagem mais direta funciona aqui (o código baseado no KDTree é mais lento (para consultar apenas um ponto)).

Abordagem de força bruta (tempo ~ 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]])) 

Executá-lo:

 $ 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 

Aqui está o script que gera milhões de pontos 3D:

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

Saída:

 $ 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 

Você poderia usar esse código para testar estruturas de dados e algoritmos mais complexos (por exemplo, se eles realmente consomem menos memory ou mais rapidamente do que a abordagem mais simples acima). Vale a pena notar que no momento é a única resposta que contém código de trabalho.

Solução baseada no KDTree (tempo ~ 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] 

Executá-lo:

 $ 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 

Classificação parcial em C ++ (tempo ~ 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 } 

Executá-lo:

 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 

Fila de prioridade em C ++ (tempo ~ 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'; } } 

Executá-lo:

 $ 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 

Abordagem baseada em pesquisa linear (tempo ~ 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)); } } 

As medições mostram que a maior parte do tempo é gasto lendo array a partir do arquivo, computações reais levam em ordem de grandeza menos tempo.

Se o milhão de inputs já estiver em um arquivo, não há necessidade de carregá-las em uma estrutura de dados na memory. Apenas mantenha uma matriz com os dez pontos mais importantes encontrados até o momento e examine mais de um milhão de pontos, atualizando sua lista dos dez mais.

Este é O (n) no número de pontos.

Você pode armazenar os pontos em uma tree k-dimensional (kd-tree). Kd-trees são otimizadas para pesquisas de vizinhos mais próximos (encontrando os n pontos mais próximos de um determinado ponto).

Eu acho que essa é uma pergunta complicada que testa se você não tentar exagerar.

Considere o algoritmo mais simples que as pessoas já deram acima: mantenha uma tabela com dez candidatos mais bem-sucedidos e passe por todos os pontos, um por um. Se você encontrar um ponto mais próximo do que qualquer um dos dez melhores, substitua-o. Qual a complexidade? Bem, temos que olhar para cada ponto do arquivo uma vez , calcular sua distância (ou quadrado da distância, na verdade) e comparar com o décimo ponto mais próximo. Se for melhor, insira-o no local apropriado na tabela 10-melhor-até-agora.

Então, qual é a complexidade? Nós olhamos para cada ponto uma vez, então são n computações da distância e n comparações. Se o ponto é melhor, precisamos inseri-lo na posição correta, isso requer mais algumas comparações, mas é um fator constante, já que a tabela de melhores candidatos é de tamanho constante 10.

Acabamos com um algoritmo que é executado em tempo linear, O (n) no número de pontos.

Mas agora considere qual é o limite inferior de tal algoritmo? Se não há ordem nos dados de input, temos que olhar para cada ponto para ver se não é um dos mais próximos. Então, tanto quanto eu posso ver, o limite inferior é Omega (n) e, portanto, o algoritmo acima é ótimo.

Não há necessidade de calcular a distância. Apenas o quadrado da distância deve atender às suas necessidades. Deve ser mais rápido eu acho. Em outras palavras, você pode pular o bit sqrt .

Esta não é uma questão de lição de casa, é? 😉

Meu pensamento: iterar sobre todos os pontos e colocá-los em um heap mínimo ou fila de prioridade limitada, codificada pela distância do destino.

Esta questão é essencialmente testar seu conhecimento e / ou intuição de algoritmos de particionamento de espaço . Eu diria que armazenar os dados em uma octree é sua melhor aposta. É comumente usado em motores 3D que lidam com esse tipo de problema (armazenando milhões de vértices, traçando raios, encontrando colisões, etc.). O tempo de pesquisa será na ordem de log(n) no pior cenário (acredito).

Algoritmo direto:

Armazene os pontos como uma lista de tuplas e escaneie os pontos, calculando a distância e mantendo uma lista “mais próxima”.

Mais criativo:

Agrupe pontos em regiões (como o cubo descrito por “0,0,0” a “50,50,50”, ou “0,0,0” a “-20, -20, -20”), para que você pode “indexar” para eles a partir do ponto alvo. Verifique em qual cubo o alvo está e procure apenas pelos pontos nesse cubo. Se houver menos de 10 pontos nesse cubo, verifique os cubos “vizinhos” e assim por diante.

Pensando melhor, esse não é um algoritmo muito bom: se o seu ponto alvo estiver mais próximo da parede de um cubo do que 10 pontos, você também terá que pesquisar no cubo vizinho.

Eu iria com a abordagem kd-tree e encontraria o mais próximo, então removesse (ou marcasse) aquele nó mais próximo, e procuraria novamente pelo novo nó mais próximo. Enxague e repita.

Para quaisquer dois pontos P1 (x1, y1, z1) e P2 (x2, y2, z2), se a distância entre os pontos for d, então todos os itens a seguir devem ser verdadeiros:

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

Segure os 10 mais próximos à medida que você percorrer todo o seu conjunto, mas também mantenha a distância até o 10º mais próximo. Poupe muita complexidade usando estas três condições antes de calcular a distância para cada ponto que você olha.

basicamente uma combinação das duas primeiras respostas acima de mim. Como os pontos estão em um arquivo, não há necessidade de mantê-los na memory. Em vez de uma matriz, ou um heap mínimo, eu usaria um heap máximo, já que você só deseja verificar as distâncias menores que o 10º ponto mais próximo. Para um array, você precisaria comparar cada distância recém-calculada com todas as 10 das distâncias que você manteve. Para um heap mínimo, você precisa realizar 3 comparações com cada nova distância calculada. Com um heap máximo, você só executa uma comparação quando a distância recém-calculada é maior que o nó raiz.

Esta questão precisa de mais definições.

1) A decisão sobre os algoritmos que pré-indexam os dados varia muito, dependendo se você pode manter os dados inteiros na memory ou não.

Com kdtrees e octrees, você não precisará manter os dados na memory e o desempenho se beneficiará desse fato, não apenas porque o tamanho da memory é menor, mas simplesmente porque você não precisará ler o arquivo inteiro.

Com o bruteforce, você terá que ler o arquivo inteiro e recompilar as distâncias para cada novo ponto que você estará procurando.

Ainda assim, isso pode não ser importante para você.

2) Outro fator é quantas vezes você terá que procurar por um ponto.

Como afirmado por JF Sebastian, às vezes, o bruteforce é mais rápido até mesmo em grandes conjuntos de dados, mas tenha em mente que seus benchmarks medem a leitura de todo o dataset do disco (o que não é necessário uma vez que o kd-tree ou octree é construído e escrito em algum lugar) e que eles medem apenas uma pesquisa.

Calcule a distância para cada um deles e selecione (1..10, n) no tempo O (n). Isso seria o algoritmo ingênuo, eu acho.