对于我的情况，最好的最近邻算法是什么？

Question

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

Answer 1

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 坐标到网格单元坐标计算，然后是一次内存访问，这意味着您不能真正希望更快的方法（每个查询可能需要 2-3 个 CPU 周期。）

我怀疑（有一些见解）这就是像 Google 或 Facebook 这样的大公司会如何解决这个问题（因为 #3 对他们来说甚至对整个世界来说都不是问题。）即使是较小的非营利组织也会使用像这样的计划这（就像 NASA。）尽管 NASA 使用的方案要复杂得多，具有多尺度的分辨率/精度。

澄清

从下面的评论来看，很明显最后一部分没有很好理解，所以我将包括更多细节。

假设您的点集由两个向量 x 和 y 给出，它们包含数据的 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 中的条目数决定了您尝试达到的精确度。显然，在我上面给出的示例中，精度很差。

Answer 2

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

我推荐OpenCV的implementation，正如this answer中提到的那样

在性能方面，您插入的点的顺序可以对查询时间产生影响，因为实现可能会选择是否重新平衡不平衡的树（例如，OpenCV 不会这样做）。一个简单的保障措施是以随机顺序插入点：首先打乱列表，然后以打乱顺序插入所有点。虽然不是最优的，但这可以确保以压倒性的概率产生的顺序不会是病态的。