Миллионы 3D точек: Как найти 10 из них самыми близкими к данной точке?

Точка в 3-м определяется (x, y, z). Расстояние d между любыми двумя точками (X, Y, Z) и (x, y, z) является d = Sqrt [(X-x) ^2 + (Y-y)^2 + (Z-z)^2]. Теперь существует миллион записей в файле, каждая запись является некоторой точкой в пространстве ни в каком определенном порядке. Учитывая любую точку (a, b, c) находят ближайшие 10 точек к нему. Как Вы сохранили бы миллион точек и как Вы получите те 10 точек от той структуры данных.

67
задан MatrixFrog 21 March 2010 в 05:56
поделиться

12 ответов

Миллион очков - это небольшое число. Здесь работает наиболее простой подход (код на основе KDTree медленнее (для запроса только одной точки)).

Подход грубой силы (время ~ 1 секунда)

#!/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]]))

Выполните:

$ 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

Вот сценарий, который генерирует миллион трехмерных точек:

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

Вывод:

$ 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

Вы можете использовать этот код для проверить более сложные структуры данных и алгоритмы (например, действительно ли они потребляют меньше памяти или быстрее, чем описанный выше простейший подход). Стоит отметить, что на данный момент это единственный ответ, содержащий рабочий код.

Решение на основе KDTree (время ~ 1,4 секунды)

#!/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]

Запустить:

$ 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

Частичная сортировка в C ++ (время ~ 1,1 секунды)

// $ 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
}

Запустить:

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

Приоритетная очередь в C ++ (время ~ 1,2 секунды)

#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';  
  }
}

Запустить:

$ 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

Подход на основе линейного поиска (время ~ 1,15 секунды)

// $ 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));
  }
}

Измерения показывают, что большую часть времени тратится на чтение массива из файла, фактические вычисления выполняются по порядку величины меньше времени.

94
ответ дан 24 November 2019 в 14:32
поделиться

Это не вопрос домашнего задания, не так ли? ; -)

Моя мысль: перебрать все точки и поместить их в минимальную кучу или очередь с ограниченным приоритетом, привязанную к расстоянию от цели.

4
ответ дан 24 November 2019 в 14:32
поделиться

Если миллион записей уже находится в файле, нет необходимости загружать их все в структуру данных в памяти. Просто сохраните массив с десятью первыми точками, найденными на данный момент, и просматривайте миллион точек, обновляя свой список первой десятки по ходу дела.

Это O (n) по количеству очков.

20
ответ дан 24 November 2019 в 14:32
поделиться

Этот вопрос нуждается в дополнительном определении.

1) Решение относительно алгоритмов, которые предварительно индексируют данные, очень сильно зависит от того, можете ли вы хранить все данные в памяти или нет.

С kdtrees и octrees вам не придется хранить данные в памяти, и производительность выигрывает от этого факта, не только потому, что занимаемая память меньше, но и просто потому, что вам не придется читать весь файл.

При использовании bruteforce вам придется читать весь файл и заново вычислять расстояния для каждой новой точки, которую вы будете искать.

Тем не менее, для вас это может быть неважно.

2) Другой фактор - сколько раз вам придется искать точку.

Как утверждает J.F. Sebastian, иногда bruteforce быстрее даже на больших наборах данных, но учтите, что его бенчмарки измеряют чтение всего набора данных с диска (что не нужно, когда kd-tree или octree построено и записано куда-то) и что они измеряют только один поиск.

0
ответ дан 24 November 2019 в 14:32
поделиться

Рассчитайте расстояние для каждого из них и выполните Select (1..10, n) за O (n) раз. Думаю, это был бы наивный алгоритм.

0
ответ дан 24 November 2019 в 14:32
поделиться

Вы можете хранить точки в k-мерном дереве (kd-дерево). Kd-деревья оптимизированы для поиска ближайших соседей (нахождение точек n, ближайших к заданной точке).

14
ответ дан 24 November 2019 в 14:32
поделиться

Нет необходимости рассчитывать расстояние. Просто квадрат расстояния должен соответствовать вашим потребностям. Думаю, должно быть быстрее. Другими словами, вы можете пропустить бит sqrt .

5
ответ дан 24 November 2019 в 14:32
поделиться

в основном комбинация первых двух ответов выше меня. так как точки находятся в файле, нет необходимости хранить их в памяти. Вместо массива или минимальной кучи я бы использовал максимальную кучу, так как вы хотите проверить только расстояния меньше 10-й ближайшей точки. Для массива вам нужно будет сравнить каждое вновь рассчитанное расстояние со всеми 10 расстояниями, которые вы сохранили. За минимальную кучу вам нужно выполнить 3 сравнения с каждым вновь рассчитанным расстоянием. При максимальной куче выполняется только 1 сравнение, если вновь вычисляемое расстояние больше корневого узла.

1
ответ дан 24 November 2019 в 14:32
поделиться

Простой алгоритм:

Храните точки в виде списка кортежей и сканируйте точки, вычисляя расстояние и сохраняя «ближайший» список.

Более творческий:

Группируйте точки по областям (например, кубу, описанной от «0,0,0» до «50,50,50» или от «0,0,0» до «-20,-20,-20»), чтобы вы могли «индексировать» их из целевой точки. Проверьте, в каком кубе находится цель, и выполните поиск только по точкам в этом кубе. Если в этом кубе меньше 10 точек, проверьте «соседние» кубики и так далее.

При дальнейшем рассмотрении, это не очень хороший алгоритм: если ваша целевая точка находится ближе к стенке куба, чем 10 точек, то вам также придется искать в соседнем кубе.

Я бы пошел с подходом kd-tree и нашел ближайший, затем удалил (или помечал) этот ближайший узел и повторно искал новый ближайший узел. Смойте и повторите.

2
ответ дан 24 November 2019 в 14:32
поделиться

Для любых двух точек P1 (x1, y1, z1) и P2 (x2, y2, z2), если расстояние между ними равно d, то все следующее должно быть верно:

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

Удерживайте 10 ближайших по мере итерации по всему множеству, но также удерживайте расстояние до 10-й ближайшей. Сэкономите себе много сложностей, используя эти три условия перед вычислением расстояния для каждой рассматриваемой точки.

2
ответ дан 24 November 2019 в 14:32
поделиться

Этот вопрос, по сути, проверяет ваши знания и / или интуицию алгоритмов разделения пространства . Я бы сказал, что лучше всего хранить данные в октодереве . Он обычно используется в 3D-движках, которые решают именно такие проблемы (хранение миллионов вершин, трассировка лучей, поиск столкновений и т. Д.). Время поиска будет порядка log (n) в наихудшем сценарии (я полагаю).

4
ответ дан 24 November 2019 в 14:32
поделиться

Я думаю, что это сложный вопрос, который проверяет, не пытаетесь ли вы переусердствовать.

Рассмотрим простейший алгоритм, который люди уже привели выше: составьте таблицу из десяти лучших на данный момент кандидатов и просмотрите все пункты один за другим. Если вы найдете точку ближе, чем любая из десяти лучших на данный момент, замените ее. В чем сложность? Что ж, мы должны просмотреть каждую точку из файла один раз , вычислить ее расстояние (или фактически квадрат расстояния) и сравнить с 10-й ближайшей точкой. Если так лучше, вставьте его в соответствующее место в таблице 10 лучших на данный момент.

Так в чем же сложность? Мы смотрим на каждую точку один раз, так что это n вычислений расстояния и n сравнений. Если точка лучше, нам нужно вставить ее в правильную позицию, для этого потребуется еще несколько сравнений, но это постоянный коэффициент, поскольку таблица лучших кандидатов имеет постоянный размер 10.

В итоге мы получаем алгоритм, который выполняется за линейное время, O (n) по количеству точек.

Но теперь рассмотрим, какова нижняя граница такого алгоритма? Если во входных данных нет порядка, мы должны смотреть на каждую точку, чтобы увидеть, не является ли она одной из ближайших. Итак, насколько я понимаю, нижняя граница - это Omega (n), и, следовательно, приведенный выше алгоритм является оптимальным.

10
ответ дан 24 November 2019 в 14:32
поделиться
Другие вопросы по тегам:

Похожие вопросы: