加速"find_simplex“从scipy.spatial.Delauna

Question

我使用Scypi的Delaunay三角剖分构建了一个应用程序。为了验证它，我想做一个保持一次的测试，这意味着下面提到的代码段被多次调用(~1e9)。因此，我想尽快做到这一点。

下面是我想要加快的最低工作示例：

from scipy.spatial import Delaunay as Delaunay
import numpy as np
import time

n_pts = 100      # around 1e9 in the real application
pts = np.random.random((n_pts, 2))

t = time.time()
for i in range(n_pts):
    delaunay = Delaunay(pts[np.arange(n_pts)!=i])
    simplex = delaunay.find_simplex(pts[i])
print(time.time()-t)

大部分时间都是由find_simplex方法使用的，在我的机器上大约有200 my的300 my。有办法加快速度吗？我已经看过了Delaunay构造函数中的qhull_options，但是没有成功。

请注意，我不能改变整体结构，因为“真正的”程序工作得很好，这个计算只是为了验证。非常感谢!

Answer 1

很难说find_simplex方法到底是什么，但我的猜测是，在第一个调用中，它构造了一些搜索结构，而且由于您只使用了一次构造初始化时间，所以没有对许多查询进行摊销。

一个简单的解决方案是在不调用find_simplex方法的情况下运行线性搜索。由于您的代码在每次迭代时构造一个Delaunay三角剖分，运行时将被三角剖分结构所支配，线性搜索时间可以忽略不计。

下面是一个矢量化的numpy函数。

 def is_in_triangle(pt, p0, p1, p2):
    """ Check in which of the triangles the point pt lies.
        p0, p1, p2 are arrays of shape (n, 2) representing vertices of triangles,
        pt is of shape (1, 2).
        Assumes p0, p1, p2 are oriented counter clockwise (as in scipy's Delaunay)
    """ 
    vp0 = pt - p0
    v10 = p1 - p0
    cross0 = vp0[:, 0] * v10[:, 1] - vp0[:, 1] * v10[:, 0]  # is pt to left/right of p0->p1

    vp1 = pt - p1
    v21 = p2 - p1
    cross1 = vp1[:, 0] * v21[:, 1] - vp1[:, 1] * v21[:, 0]  # is pt to left/right of p1->p2
    
    vp2 = pt - p2
    v02 = p0 - p2
    cross2 = vp2[:, 0] * v02[:, 1] - vp2[:, 1] * v02[:, 0]  # is pt to left/right of p2->p0

    return (cross0 < 0) & (cross1 < 0) & (cross2 < 0)  # pt should be to the left of all triangle edges

我修改了您的代码以使用以下函数运行：

t = time.time()
for i in range(n_pts):
    delaunay = Delaunay(pts[np.arange(n_pts)!=i])

    p0 = delaunay.points[delaunay.simplices[:, 0], :]
    p1 = delaunay.points[delaunay.simplices[:, 1], :]
    p2 = delaunay.points[delaunay.simplices[:, 2], :]
    pt = pts[i].reshape((1, 2))
    pt_in_simps_mask = is_in_triangle(pt, p0, p1, p2)
    simp_ind_lst = np.where(pt_in_simps_mask)[0]
    if len(simp_ind_lst) == 0:
        simp_ind = -1
    else:
        simp_ind = simp_ind_lst[0]

print("new time: {}".format(time.time()-t))

在我的机器上，当运行100分时，这段代码在大约0.036秒内运行，而原来代码的运行时间是0.13s (对于没有查询的代码，只有Delaunay构造)。