找到一组平面点的最近线段的快速方法

Question

问题

给出平面上 n=10000 个点和 m=1000 条线段。任务是确定距离每个点最近的线段。首选 Python 中的快速解决方案。

规格

线段由其端点定义，距离为欧几里德距离。线可以交叉。最小距离应计算到整个线段，而不仅仅是端点。到最近线段的距离是不必要的。

图解示例

这张图演示了如何找到少量点和线段的最近线段。这些点根据最近的线段着色。

可能的解决策略

实现此目的的一种可能方法是使用线段的 Voronoi 细分并确定点位于哪个单元格（link）。然而，计算这种镶嵌并不容易，并且需要外部包。跨线也可能是一个问题。此外，我并不关心沃罗诺伊单元的确切形状。虽然我不想排除这种方法，但通过向量化距离计算似乎可能有一个更简单的解决方案。

图片显示了上一个例子的Voronoi曲面细分。 Voronoi 单元根据线段颜色进行着色。

相关问题

这个问题已经在这篇post中提出过，但是这里的问题包含很多点和段。速度是一个问题。

Answer 1

一次循环向量化可提升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))