给定两个 2d 点列表，如何为第一个列表中的每个点找到第二个列表中的最近点？
Given two lists of 2d points, how to find the closest point in the 2nd list for every point in the 1st list?
我有两个随机排序的 2d 点的大型 numpy 数组，假设它们是 A 和 B。我需要做的是找到两个数组之间 "matches" 的数量，其中匹配项是A 中的一个点（称为 A'）在某个给定的半径 R 内，B 中的一个点（称为 B'）。这意味着 A 中的每个点都必须与 B 中的 1 点或不匹配。 return 两个数组之间匹配的列表索引也很好，但这不是必需的。由于在这个半径 R 中可以有很多点，所以最好在 B 中找到最靠近 A' 的点，然后检查它是否在半径 R 内。这个简单地用距离公式 dx^2 + dy^2 进行测试.显然有循环遍历两个数组的蛮力 O(n^2) 解决方案，但我需要更快的东西，希望 O(n log n).
我看到的是 Voronoi 图可用于解决此类问题，但我不确定这将如何实现。我不熟悉 Voronoi 图，所以我用 scipy.spatial.Voronoi 生成它。是否有使用这些图解决此问题的快速算法，或者是否有其他算法？
您可能想要的是 KDTrees（它在高维度上很慢，但对于您的问题应该非常快。python 实现甚至实现了半径边界。
我认为有几种选择。我开始了一个小的比较测试来探索一些。这些中的第一对只是找出有多少点在彼此的半径范围内，以确保我在问题的主要部分得到一致的结果。它没有回答关于找到最近的问题的邮件，我认为这只是其中一些工作的更多工作 - 最后一个选项是否已完成，请参阅 post 的底部。问题的驱动因素是进行所有比较，我认为您可以通过一些排序（此处最后一个概念）来限制比较。
天真Python
使用暴力点对点比较。显然 O(n^2)。
Scipy 的 cdist 模块
对 "small" 数据的处理速度很快。对于大数据，由于内存中矩阵输出的大小，这开始爆炸。对于 1M x 1M 应用程序可能不可行。
Scipy 的 KDTree 模块
来自其他解决方案。快，但不如 cdist 或 "sectioning"（下图）快。也许有一种不同的方法可以使用 KDTree 来完成这项任务……我对此不是很有经验。这种方法（如下）似乎合乎逻辑。
对比较对象数组进行分段
这非常有效，因为您对 所有 的距离不感兴趣，您只需要半径内的距离。因此，通过对目标数组进行排序并仅在其周围的矩形 window 内查找 "contenders" 你可以获得非常快的性能 w/ native python 而没有 "memory explosion." 可能仍然是位 "left on the table" 可以通过在此实现中嵌入 cdist 或 (gulp) 尝试对其进行多线程处理来进行增强。
其他想法...
这是一个紧密的 "mathy" 循环，因此在 cython 中尝试一些东西或拆分其中一个数组和多线程会很新颖。并对结果进行酸洗，这样你就不必 运行 这通常看起来很谨慎。
我认为您可以很容易地使用数组中的索引扩充元组以获得匹配项列表。
我的旧 iMac 通过分段在 90 秒内完成 100K x 100K，因此这对于 1M x 1M 来说不是好兆头
比较：
# distance checker
from random import uniform
import time
import numpy as np
from scipy.spatial import distance, KDTree
from bisect import bisect
from operator import itemgetter
import sys
from matplotlib import pyplot as plt
sizes = [100, 500, 1000, 2000, 5000, 10000, 20000]
#sizes = [20_000, 30_000, 40_000, 50_000, 60_000] # for the playoffs. :)
naive_times = []
cdist_times = []
kdtree_times = []
sectioned_times = []
delta = 0.1
for size in sizes:
print(f'\n *** running test with vectors of size {size} ***')
r = 20 # radius to match
r_squared = r**2
A = [(uniform(-1000,1000), uniform(-1000,1000)) for t in range(size)]
B = [(uniform(-1000,1000), uniform(-1000,1000)) for t in range(size)]
# naive python
print('naive python')
tic = time.time()
matches = [(p1, p2) for p1 in A
for p2 in B
if (p1[0] - p2[0])**2 + (p1[1] - p2[1])**2 <= r_squared]
toc = time.time()
print(f'found: {len(matches)}')
naive_times.append(toc-tic)
print(toc-tic)
print()
# using cdist module
print('cdist')
tic = time.time()
dist_matrix = distance.cdist(A, B, 'euclidean')
result = np.count_nonzero(dist_matrix<=r)
toc = time.time()
print(f'found: {result}')
cdist_times.append(toc-tic)
print(toc-tic)
print()
# KDTree
print('KDTree')
tic = time.time()
my_tree = KDTree(A)
results = my_tree.query_ball_point(B, r=r)
# for count, r in enumerate(results):
# for t in r:
# print(count, A[t])
result = sum(len(lis) for lis in results)
toc = time.time()
print(f'found: {result}')
kdtree_times.append(toc-tic)
print(toc-tic)
print()
# python with sort and sectioning
print('with sort and sectioning')
result = 0
tic = time.time()
B.sort()
for point in A:
# gather the neighborhood in x-dimension within x-r <= x <= x+r+1
# if this has any merit, we could "do it again" for y-coord....
contenders = B[bisect(B,(point[0]-r-delta, 0)) : bisect(B,(point[0]+r+delta, 0))]
# further chop down to the y-neighborhood
# flip the coordinate to support bisection by y-value
contenders = list(map(lambda p: (p[1], p[0]), contenders))
contenders.sort()
contenders = contenders[bisect(contenders,(point[1]-r-delta, 0)) :
bisect(contenders,(point[1]+r+delta, 0))]
# note (x, y) in contenders is still inverted, so need to index properly
matches = [(point, p2) for p2 in contenders if (point[0] - p2[1])**2 + (point[1] - p2[0])**2 <= r_squared]
result += len(matches)
toc = time.time()
print(f'found: {result}')
sectioned_times.append(toc-tic)
print(toc-tic)
print('complete.')
plt.plot(sizes, naive_times, label = 'naive')
plt.plot(sizes, cdist_times, label = 'cdist')
plt.plot(sizes, kdtree_times, label = 'kdtree')
plt.plot(sizes, sectioned_times, label = 'sectioning')
plt.legend()
plt.show()
其中一种尺寸和地块的结果：
*** running test with vectors of size 20000 ***
naive python
found: 124425
101.40657806396484
cdist
found: 124425
2.9293079376220703
KDTree
found: 124425
18.166933059692383
with sort and sectioning
found: 124425
2.3414530754089355
complete.
注意：在第一个图中，cdist 覆盖了 sectioning。季后赛显示在第二个情节中。
"playoffs"
修改切片代码
此代码在半径内的点内找到最小值。运行时相当于上面的分段代码。
print('with sort and sectioning, and min finding')
result = 0
pairings = {}
tic = time.time()
B.sort()
def dist_squared(a, b):
# note (x, y) in point b will be inverted (below), so need to index properly
return (a[0] - b[1])**2 + (a[1] - b[0])**2
for idx, point in enumerate(A):
# gather the neighborhood in x-dimension within x-r <= x <= x+r+1
# if this has any merit, we could "do it again" for y-coord....
contenders = B[bisect(B,(point[0]-r-delta, 0)) : bisect(B,(point[0]+r+delta, 0))]
# further chop down to the y-neighborhood
# flip the coordinate to support bisection by y-value
contenders = list(map(lambda p: (p[1], p[0]), contenders))
contenders.sort()
contenders = contenders[bisect(contenders,(point[1]-r-delta, 0)) :
bisect(contenders,(point[1]+r+delta, 0))]
matches = [(dist_squared(point, p2), point, p2) for p2 in contenders
if dist_squared(point, p2) <= r_squared]
if matches:
pairings[idx] = min(matches)[1] # pair the closest point in B with the point in A
toc = time.time()
print(toc-tic)
我有两个随机排序的 2d 点的大型 numpy 数组，假设它们是 A 和 B。我需要做的是找到两个数组之间 "matches" 的数量，其中匹配项是A 中的一个点（称为 A'）在某个给定的半径 R 内，B 中的一个点（称为 B'）。这意味着 A 中的每个点都必须与 B 中的 1 点或不匹配。 return 两个数组之间匹配的列表索引也很好，但这不是必需的。由于在这个半径 R 中可以有很多点，所以最好在 B 中找到最靠近 A' 的点，然后检查它是否在半径 R 内。这个简单地用距离公式 dx^2 + dy^2 进行测试.显然有循环遍历两个数组的蛮力 O(n^2) 解决方案，但我需要更快的东西，希望 O(n log n).
我看到的是 Voronoi 图可用于解决此类问题，但我不确定这将如何实现。我不熟悉 Voronoi 图，所以我用 scipy.spatial.Voronoi 生成它。是否有使用这些图解决此问题的快速算法，或者是否有其他算法？
您可能想要的是 KDTrees（它在高维度上很慢，但对于您的问题应该非常快。python 实现甚至实现了半径边界。
我认为有几种选择。我开始了一个小的比较测试来探索一些。这些中的第一对只是找出有多少点在彼此的半径范围内，以确保我在问题的主要部分得到一致的结果。它没有回答关于找到最近的问题的邮件，我认为这只是其中一些工作的更多工作 - 最后一个选项是否已完成，请参阅 post 的底部。问题的驱动因素是进行所有比较，我认为您可以通过一些排序（此处最后一个概念）来限制比较。
天真Python
使用暴力点对点比较。显然 O(n^2)。
Scipy 的 cdist 模块
对 "small" 数据的处理速度很快。对于大数据，由于内存中矩阵输出的大小，这开始爆炸。对于 1M x 1M 应用程序可能不可行。
Scipy 的 KDTree 模块
来自其他解决方案。快，但不如 cdist 或 "sectioning"（下图）快。也许有一种不同的方法可以使用 KDTree 来完成这项任务……我对此不是很有经验。这种方法（如下）似乎合乎逻辑。
对比较对象数组进行分段
这非常有效，因为您对 所有 的距离不感兴趣，您只需要半径内的距离。因此，通过对目标数组进行排序并仅在其周围的矩形 window 内查找 "contenders" 你可以获得非常快的性能 w/ native python 而没有 "memory explosion." 可能仍然是位 "left on the table" 可以通过在此实现中嵌入 cdist 或 (gulp) 尝试对其进行多线程处理来进行增强。
其他想法...
这是一个紧密的 "mathy" 循环，因此在 cython 中尝试一些东西或拆分其中一个数组和多线程会很新颖。并对结果进行酸洗，这样你就不必 运行 这通常看起来很谨慎。
我认为您可以很容易地使用数组中的索引扩充元组以获得匹配项列表。
我的旧 iMac 通过分段在 90 秒内完成 100K x 100K，因此这对于 1M x 1M 来说不是好兆头
比较：
# distance checker
from random import uniform
import time
import numpy as np
from scipy.spatial import distance, KDTree
from bisect import bisect
from operator import itemgetter
import sys
from matplotlib import pyplot as plt
sizes = [100, 500, 1000, 2000, 5000, 10000, 20000]
#sizes = [20_000, 30_000, 40_000, 50_000, 60_000] # for the playoffs. :)
naive_times = []
cdist_times = []
kdtree_times = []
sectioned_times = []
delta = 0.1
for size in sizes:
print(f'\n *** running test with vectors of size {size} ***')
r = 20 # radius to match
r_squared = r**2
A = [(uniform(-1000,1000), uniform(-1000,1000)) for t in range(size)]
B = [(uniform(-1000,1000), uniform(-1000,1000)) for t in range(size)]
# naive python
print('naive python')
tic = time.time()
matches = [(p1, p2) for p1 in A
for p2 in B
if (p1[0] - p2[0])**2 + (p1[1] - p2[1])**2 <= r_squared]
toc = time.time()
print(f'found: {len(matches)}')
naive_times.append(toc-tic)
print(toc-tic)
print()
# using cdist module
print('cdist')
tic = time.time()
dist_matrix = distance.cdist(A, B, 'euclidean')
result = np.count_nonzero(dist_matrix<=r)
toc = time.time()
print(f'found: {result}')
cdist_times.append(toc-tic)
print(toc-tic)
print()
# KDTree
print('KDTree')
tic = time.time()
my_tree = KDTree(A)
results = my_tree.query_ball_point(B, r=r)
# for count, r in enumerate(results):
# for t in r:
# print(count, A[t])
result = sum(len(lis) for lis in results)
toc = time.time()
print(f'found: {result}')
kdtree_times.append(toc-tic)
print(toc-tic)
print()
# python with sort and sectioning
print('with sort and sectioning')
result = 0
tic = time.time()
B.sort()
for point in A:
# gather the neighborhood in x-dimension within x-r <= x <= x+r+1
# if this has any merit, we could "do it again" for y-coord....
contenders = B[bisect(B,(point[0]-r-delta, 0)) : bisect(B,(point[0]+r+delta, 0))]
# further chop down to the y-neighborhood
# flip the coordinate to support bisection by y-value
contenders = list(map(lambda p: (p[1], p[0]), contenders))
contenders.sort()
contenders = contenders[bisect(contenders,(point[1]-r-delta, 0)) :
bisect(contenders,(point[1]+r+delta, 0))]
# note (x, y) in contenders is still inverted, so need to index properly
matches = [(point, p2) for p2 in contenders if (point[0] - p2[1])**2 + (point[1] - p2[0])**2 <= r_squared]
result += len(matches)
toc = time.time()
print(f'found: {result}')
sectioned_times.append(toc-tic)
print(toc-tic)
print('complete.')
plt.plot(sizes, naive_times, label = 'naive')
plt.plot(sizes, cdist_times, label = 'cdist')
plt.plot(sizes, kdtree_times, label = 'kdtree')
plt.plot(sizes, sectioned_times, label = 'sectioning')
plt.legend()
plt.show()
其中一种尺寸和地块的结果：
*** running test with vectors of size 20000 ***
naive python
found: 124425
101.40657806396484
cdist
found: 124425
2.9293079376220703
KDTree
found: 124425
18.166933059692383
with sort and sectioning
found: 124425
2.3414530754089355
complete.
注意：在第一个图中，cdist 覆盖了 sectioning。季后赛显示在第二个情节中。
"playoffs"
修改切片代码
此代码在半径内的点内找到最小值。运行时相当于上面的分段代码。
print('with sort and sectioning, and min finding')
result = 0
pairings = {}
tic = time.time()
B.sort()
def dist_squared(a, b):
# note (x, y) in point b will be inverted (below), so need to index properly
return (a[0] - b[1])**2 + (a[1] - b[0])**2
for idx, point in enumerate(A):
# gather the neighborhood in x-dimension within x-r <= x <= x+r+1
# if this has any merit, we could "do it again" for y-coord....
contenders = B[bisect(B,(point[0]-r-delta, 0)) : bisect(B,(point[0]+r+delta, 0))]
# further chop down to the y-neighborhood
# flip the coordinate to support bisection by y-value
contenders = list(map(lambda p: (p[1], p[0]), contenders))
contenders.sort()
contenders = contenders[bisect(contenders,(point[1]-r-delta, 0)) :
bisect(contenders,(point[1]+r+delta, 0))]
matches = [(dist_squared(point, p2), point, p2) for p2 in contenders
if dist_squared(point, p2) <= r_squared]
if matches:
pairings[idx] = min(matches)[1] # pair the closest point in B with the point in A
toc = time.time()
print(toc-tic)