在有限集中寻找到另一个点的最近点的有效算法

Posted 2023-02-22

【中文标题】在有限集中寻找到另一个点的最近点的有效算法

【英文标题】：Efficient algorithm for finding closest point in a finite set to another point

【发布时间】：2017-04-12 14:42:44

【问题描述】：

我有一个包含约 30k 个位置的列表 L（写为经度/纬度对）和一个包含约 1m 个事件的列表 E（位置写为经度/纬度对），每个事件都发生在 L 中的一个点。我想用它在 L 中的相应位置标记 E 中的每个事件。但是 L 和 E 中的坐标四舍五入不同——E 到小数点后五位，L 到十三位——所以表面上相同的坐标实际上可以相差 ~10^-5 度，或约 1 米。 （L 中的点至少相隔约 10 m。）

因此我需要 L 中离 E 的每个点最近的点；明显的 O(|L||E|) 蛮力算法太慢了。 L 与 E 相比足够小，因此预处理 L 并将预处理时间摊销到 E 上的算法很好。这是一个经过充分研究的问题吗？我能找到的链接是针对相关但不同的问题，例如找到一组中一对点之间的最小距离。

可能相关：Voronoi diagrams，虽然我看不出将 L 预处理成 Voronoi 图如何节省我的计算时间。

【问题讨论】：

好吧..我已经跳过了你的大部分帖子，因为有太多无用的信息。基本上，您有一个包含 30k XY 点的列表“L”（如果需要，可以称为纬度/经度），并且您有一个包含一百万个 XY 点的列表“E”（同上），并且您想知道每个点E 最接近“L”中的哪个。是这样吗？请确认。

等等 - 如果舍入误差导致最大偏差为 1m，但点之间的最小距离为 10m，那么您可以将两者舍入到相同的精度并比较相等，对吧？

@AlexG：可以确认（我已经编辑了一些无用的信息）。

@ConnorHarris 谢谢。看起来四叉树应该可以完成这项工作。 fr.wikipedia.org/wiki/Quadtree 可能有很多开源库可以实现它，但我认为您可以很容易地自己完成。按 x 对所有参考点进行排序。将前半部分和后半部分按 y 排序。您刚刚将您的点分为 4 个区域。只需在构建树结构时递归执行此操作，直到每个区域仅包含 1 个点。之后，您可以通过比较范围和遍历树来搜索最近的点。

【参考方案1】：

是的，你是对的。首先，您可以使用Furtune's Sweep Line 方法在 O(|L| log |L|) 时间内构建位置 L 的点集的 Voronoi 图。有多种实现可供使用，Triangle 是最常见的实现之一。

现在你有一个 O(|L|) 大小的平面分区。要允许 O(log |L|) 最近邻查询，您需要在 Voronoi 图之上的搜索结构。一种常见的方法是使用 Dobkin-Kirkpatrick Hierarchy，详细信息可以在各种lecture notes 中找到。此方法支持 O(log |L|) 查询，并且只需要 O(|L|) 大小。 （在this post中也提到过。）

然后是 |E|查询可以在 O(|E| log |L|) 时间内完成。

另一种方法是使用k-d trees。从实现的角度来看，它们的工作量可能更少并且提供相同的复杂性（据我所知）。 快速搜索发现这两个可能值得测试的实现：C++、Java。

【讨论】：

Voronoi + 分区理论上是实点的最佳方法。在实践中，它是一个庞大而复杂的机器，特殊情况很难处理。在实践中，我还会推荐一个二维树，它更容易实现。

是的，我同意你的观点，这就是我最后提到树木的原因。不过，很高兴知道如何在这方面使用 Voronoi 图。

【参考方案2】：

这适用于space-partitioning 的相当直接的应用。大多数 GIS 库都应提供执行此操作的工具。

我自己的经验仅限于使用R-Trees。我创建了 L 个项目的索引。您可以直接使用这些点，也可以使用表示该点周围不确定性的边界框。然后，R-tree 支持对 n 个最近的相邻点/边界框进行有效的（log(L) 时间）查询。

我成功使用的实现是将libspatialindex包装到一个名为Rtree的python库中。

但是，请注意，此特定工具仅在使用欧几里得 x,y 坐标时才准确。使用远离赤道的 lat long 很容易出错，特别是如果你想覆盖一个地理上很大的区域。就我而言，我被限制在一个使用东/北的国家。我不知道哪些库支持使用大圆距离处理问题，但它们当然应该可用。

【参考方案3】：

根据您的描述：

E 中的每个点都与 L 中的点相同，四舍五入到小数点后五位。因此，E 中的每个点与 L 中的匹配点最多偏离约 1 米。 L 中的点至少相隔约 10 米。

解决方案：将 L 中的坐标四舍五入到与 E 相同的精度并匹配相等的对。

解释：四舍五入有效地将每个点映射到正方形网格上的最近邻居。只要网格分辨率（舍入到小数点后约 1 米）低于 L 中两点之间最小距离的一半（约 10/2 米），您就可以开始了。

为了获得最佳性能和利用 |L|

【讨论】：

有趣的方法。尽管在某些情况下，使用不同策略四舍五入的值会相差一个单位，但存在危险。万一错过，也必须尝试相邻的垃圾箱。

【参考方案4】：

我的库 GeographicLib 包含一个类 NearestNeighbor 它实现了vantage-point trees。这适用于任何数据 有真实距离度量的地方；测地距离明显 满足这个标准。 基本上你用你的L 位置初始化这个类 然后对于每个E 事件，你 询问最近的位置。计算上的细分 费用如下

set-up cost: L log L
cost per query: log L
cost of all queries: E log L
total cost: (L + E) log L

（这些是 base-2 日志。）这里的代码将进行查找。在 30k 位置和 1M 事件上运行此程序大约需要 40 秒，并涉及总共 16M 测地线距离计算。 （蛮力方式大约需要 21 小时。）

// Read lat/lon locations from locations.txt and lat/lon events from
// events.txt.  For each event print to closest.txt: the index for the
// closest location and the distance to it.

// Needs GeographicLib 1.47 or later

// compile and link with
// g++ -o closest closest.cpp -lGeographic \
//   -L/usr/local/lib -Wl,-rpath=/usr/local/lib

#include <iostream>
#include <vector>
#include <fstream>
#include <GeographicLib/NearestNeighbor.hpp>
#include <GeographicLib/Geodesic.hpp>

using namespace std;
using namespace GeographicLib;

// A structure to hold a geographic coordinate.
struct pos 
  double lat, lon;
  pos(double lat = 0, double lon = 0)
    : lat(lat), lon(lon) 
;

// A class to compute the distance between 2 positions.
class DistanceCalculator 
private:
  Geodesic _geod;
public:
  explicit DistanceCalculator(const Geodesic& geod)
    : _geod(geod) 
  double operator() (const pos& a, const pos& b) const 
    double s12;
    _geod.Inverse(a.lat, a.lon, b.lat, b.lon, s12);
    return s12;
  
;

int main() 
  // Define a distance function object
  DistanceCalculator distance(Geodesic::WGS84());

  // Read in locations
  vector<pos> locs;
  
    double lat, lon;
    ifstream is("locations.txt");
    while (is >> lat >> lon)
      locs.push_back(pos(lat, lon));
  

  // Construct NearestNeighbor object
  NearestNeighbor<double, pos, DistanceCalculator>
    locationset(locs, distance);

  ifstream is("events.txt");
  ofstream os("closest.txt");

  double lat, lon, d;
  vector<int> k;
  while (is >> lat >> lon) 
    pos event(lat, lon);
    d = locationset.Search(locs, distance, event, k);
    os << k[0] << " " << d << "\n";