对于我的情况,最佳的最近邻居算法是什么? 澄清

问题描述

我有一个预定义的gps位置列表,基本上可以定义一个预定义的汽车轨迹。列表中大约有 15000 个点。整个列表是事先已知的,之后不需要插入任何点。然后,我绕了 100万个额外的采样gps位置,因此需要在预定义列表中查找最近的邻居。我需要在一次迭代中处理所有100万个项目,并且需要尽快完成。在这种情况下,最佳的最近邻居算法是什么? 我可以根据需要尽可能多地预处理预定义列表,但是处理100万个项目应尽可能快。
我已经测试了KDTree c#实现,但是性能似乎很差,也许对于我的2D数据存在一个更合适的算法。 (我的情况忽略了GPS高度) 谢谢您的任何建议!

My gps track

enter image description here

解决方法

K-D树确实非常适合该问题。您应该首先使用已知的良好实现再试一次,如果性能不够好,则可以轻松地并行化查询-由于每个查询完全独立于其他查询,因此可以通过并行处理N个查询来实现N个加速,如果您有足够的硬件。

我推荐OpenCV的implementation,如this answer

所述

在性能方面,插入 的点的顺序可能会影响查询时间,因为实现可能会选择是否重新平衡不平衡的树(例如,OpenCV不会这样做)。一种简单的保护措施是按随机顺序插入点:首先对列表进行随机排序,然后按随机顺序插入所有点。尽管这不是最佳方法,但可以确保以绝对的可能性使所得到的顺序不致于病态。

,

CGAL具有2d point library,用于基于Delaunay三角剖分数据结构的最近邻和范围搜索。

以下是您的用例的库基准:

// file: cgal_benchmark_2dnn.cpp
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Point_set_2.h>
#include <chrono>
#include <list>
#include <random>

typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Point_set_2<K>::Vertex_handle Vertex_handle;
typedef K::Point_2 Point_2;

/**
 * @brief Time a lambda function.
 *
 * @param lambda - the function to execute and time
 *
 * @return the number of microseconds elapsed while executing lambda
 */
template <typename Lambda>
std::chrono::microseconds time_lambda(Lambda lambda) {
  auto start_time = std::chrono::high_resolution_clock::now();
  lambda();
  auto end_time = std::chrono::high_resolution_clock::now();
  return std::chrono::duration_cast<std::chrono::microseconds>(end_time -
                                                               start_time);
}

int main() {
  const int num_index_points = 15000;
  const int num_trials = 1000000;

  std::random_device
      rd; // Will be used to obtain a seed for the random number engine
  std::mt19937 gen(rd()); // Standard mersenne_twister_engine seeded with rd()
  std::uniform_real_distribution<> dis(-1,1.);
  std::list<Point_2> index_point_list;

  {
    auto elapsed_microseconds = time_lambda([&] {
      for (int i = 0; i < num_index_points; ++i) {
        index_point_list.emplace_back(dis(gen),dis(gen));
      }
    });
    std::cout << " Generating " << num_index_points << " random points took "
              << elapsed_microseconds.count() << " microseconds.\n";
  }

  CGAL::Point_set_2<K> point_set;
  {
    auto elapsed_microseconds = time_lambda([&] {
      point_set.insert(index_point_list.begin(),index_point_list.end());
    });
    std::cout << " Building point set took " << elapsed_microseconds.count()
              << " microseconds.\n";
  }

  {
    auto elapsed_microseconds = time_lambda([&] {
      for (int j = 0; j < num_trials; ++j) {
        Point_2 query_point(dis(gen),dis(gen));
        Vertex_handle v = point_set.nearest_neighbor(query_point);
      }
    });
    auto rate = elapsed_microseconds.count() / static_cast<double>(num_trials);
    std::cout << " Querying " << num_trials << " random points took "
              << elapsed_microseconds.count()
              << " microseconds.\n >> Microseconds / query :" << rate << "\n";
  }
}

在我的系统(Ubuntu 18.04)上,可以使用以下命令进行编译

g++ cgal_benchmark_2dnn.cpp -lCGAL -lgmp -O3

并在运行时产生性能:

 Generating 15000 random points took 1131 microseconds.
 Building point set took 11469 microseconds.
 Querying 1000000 random points took 2971201 microseconds.
 >> Microseconds / query :2.9712

这是相当快的。请注意,使用N个处理器,您可以将其加速大约N倍。

最快的实施方式

如果以下两个或多个是正确的:

  1. 您有一个小边框,可容纳150000个索引点
  2. 您只关心一个小数点以下的精度(请注意,对于经纬度坐标,多于6个小数点会产生厘米/毫米的比例精度)
  3. 您的系统上有大量的内存

然后缓存所有内容!。您可以在索引点的边界框上预先计算所需精度的网格。将每个网格单元映射到一个唯一的地址,该地址可以在知道查询点的2D坐标的情况下进行索引。

然后简单地使用任何最近的邻居算法(例如我提供的算法)将每个网格单元映射到最近的索引点。请注意,只需初始化一次该步骤即可初始化网格中的网格单元。

要运行查询,这将需要一个2D坐标到网格单元坐标的计算,再进行一次内存访问,这意味着您不能真正希望有一个更快的方法(每个查询可能需要2-3个CPU周期。)

我怀疑(有一些洞察力)这是像Google或Facebook这样的大公司如何解决该问题的方法(因为即使对于整个世界,#3对他们来说也不是问题。)甚至较小的非营利组织也使用类似的方案尽管与NASA一样,但NASA使用的方案要复杂得多,具有多种分辨率/精度。

澄清

从下面的评论中,最后一节的内容很不清楚,因此我将提供更多详细信息。

假设您的点集由两个向量xy给出,它们包含数据的x和y坐标(或lat&long或您正在使用的任何东西)。

然后,数据的边界框定义为维度width = max(x)-min(x)height=max(y)-min(y)。 现在,使用一组测试点(x_t,y_t)的映射,使用NxM点创建一个精细的网格网格来表示整个边界框

u(x_t) = round((x_t - min(x)) / double(width) * N)
v(y_t) = round((y_t - min(y)) / double(height) * M)

然后只需使用indices = grid[u(x_t),v(y_t)],其中indices是最接近[x_t,y_t]的索引点的索引,而grid是预先计算的查找表,它映射网格中的每个项目到最近的索引点[x,y]

例如,假设您的索引点是[0,0][2,2](按此顺序。)您可以将网格创建为

grid[0,0] = 0
grid[0,1] = 0
grid[0,2] = 0 // this is a tie
grid[1,0] = 0
grid[1,1] = 0 // this is a tie
grid[1,2] = 1 
grid[2,0] = 1 // this is a tie
grid[2,1] = 1
grid[2,2] = 1

其中上面的右手边是索引0(映射到点[0,0])或1(映射到点[2,2])。注意:由于这种方法具有离散性,因此您将拥有到一个点的距离等于到另一个索引点的距离的联系,因此您必须想出一些方法来确定如何打破这些联系。请注意,grid中的条目数决定了您要达到的精确度。显然,在我上面给出的示例中,精度非常糟糕。