问题
给出平面上 n=10000 个点和 m=1000 条线段。任务是确定距离每个点最近的线段。首选 Python 中的快速解决方案。
规格
线段由其端点定义,距离为欧几里德距离。线可以交叉。最小距离应计算到整个线段,而不仅仅是端点。到最近线段的距离是不必要的。
图解示例
这张图演示了如何找到少量点和线段的最近线段。这些点根据最近的线段着色。
可能的解决策略
实现此目的的一种可能方法是使用线段的 Voronoi 细分并确定点位于哪个单元格(link)。然而,计算这种镶嵌并不容易,并且需要外部包。跨线也可能是一个问题。此外,我并不关心沃罗诺伊单元的确切形状。虽然我不想排除这种方法,但通过向量化距离计算似乎可能有一个更简单的解决方案。
图片显示了上一个例子的Voronoi曲面细分。 Voronoi 单元根据线段颜色进行着色。
相关问题
这个问题已经在这篇post中提出过,但是这里的问题包含很多点和段。速度是一个问题。
解决方案提示
一个快速的解决方案已经被提出。对于给定的 CPU,n=10000、m=1000 需要 0.2 秒,n=100000、m=10000 需要 7.5 秒。
是否可以在 Python 中实现更快的算法,与现有解决方案相比,速度至少提高 30%? 更准确地说,如果现有的解决方案(函数
dist_fast
)在单个CPU上花费x时间,那么改进的解决方案最多应该花费0.7x时间。为了更好的比较,可以测试 n=10000,m=1000 和 n=100000,m=10000 的情况。
一次循环向量化可提升100倍速度
通过修改这篇post中的Python算法,我能够将速度提高100倍(下面代码中的函数
dist_fast
)。我对 m 进行向量化并提前计算重复步骤。
对于 n=10000 和 m=1000,当原始算法在 n 和 m 上循环时,给定的 CPU 需要 20 秒,但当我们仅在 m 上循环时,只需要 0.2 秒。
import numpy as np, time
def dist_fast(p,segment):
n=p.shape[1]
segment_closest=np.zeros(n,dtype=int)
dxs=segment[2,:]-segment[0,:]
dys=segment[3,:]-segment[1,:]
q=dxs**2+dys**2;dxsq=dxs/q;dysq=dys/q
for i in range(n):
dx=p[0,i]-segment[0];dy=p[1,i]-segment[1]
u=np.clip(dx*dxsq+dy*dysq,0,1)
sqdist = (dx-u*dxs)**2+(dy-u*dys)**2
segment_closest[i]=np.argmin(sqdist)
return segment_closest
def dist_slow(p,segment):
#use distance function from https://stackoverflow.com/a/2233538
def dist(x1, y1, x2, y2, x3, y3):
px = x2-x1
py = y2-y1
norm = px*px + py*py
u = ((x3 - x1) * px + (y3 - y1) * py) / float(norm)
if u > 1:
u = 1
elif u < 0:
u = 0
x = x1 + u * px
y = y1 + u * py
dx = x - x3
dy = y - y3
sqdist = dx*dx + dy*dy
return sqdist
n=p.shape[1];m=segment.shape[1]
segment_closest=np.zeros(n,dtype=int)
d=np.zeros(m)
for i in range(n):
x3,y3=p[0,i],p[1,i]
for j in range(m):
x1,y1,x2,y2=segment[0,j],segment[1,j],segment[2,j],segment[3,j]
d[j]=dist(x1,y1,x2,y2,x3,y3)
segment_closest[i]=np.argmin(d)
return segment_closest
n=10**4 #number of points
m=10**3 #number of line segments
#define repeatable random numbers
np.random.seed(0)
#create random points and line segments for testing
p=np.random.rand(2,n) #x,y of points
segment=np.random.rand(4,m) #x1,y1,x2,y2 start and end of line segments
#fast algorithm, vectorized, 1 loop
t1=time.time()
segment_closest1=dist_fast(p,segment)
time_fast=time.time()-t1
#slow algorithm, 2 loops
t1=time.time()
segment_closest2=dist_slow(p,segment)
time_slow=time.time()-t1
print('Time vectorized algorithm: ',time_fast)
print('Time non-vectorized algorithm: ',time_slow)
print('Results are equal?',np.array_equal(segment_closest1,segment_closest2))