Миллионы трехмерных точек: как найти 10 из них, ближайших к заданной точке?

Question

Точка в 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 точек из этой структуры данных.

Answer 1

Миллион очков - это небольшое число. Здесь работает самый простой подход (код, основанный на 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

Вот скрипт, который генерирует миллион 3D-точек:

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

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

Answer 2

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

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

Answer 3

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

Answer 4

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

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

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

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

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

Answer 5

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

Answer 6

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

Answer 7

Это не домашнее задание, не так ли? ; -)

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

Answer 8

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

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

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

Answer 9

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

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

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

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

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

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

Answer 10

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