从 scipy.spatial.Delauna 加速“find_simplex”

Posted 2023-03-16

【中文标题】从 scipy.spatial.Delauna 加速“find_simplex”

【英文标题】：Speed up "find_simplex" from scipy.spatial.Delauna

【发布时间】：2021-12-23 03:30:24

【问题描述】：

我使用 Scypi 的 Delaunay 三角测量构建了一个应用程序。为了验证它，我想做一个 Hold-One-Out 测试，这意味着下面提到的代码 sn-p 被调用了很多次（~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 毫秒到 300 毫秒。有什么办法可以加快速度吗？我已经查看了 Delaunay 构造函数中的 qhull_options，但是我没有成功。

请注意，我无法更改整体结构，因为“真实”程序运行良好，并且此计算仅用于验证。非常感谢！

【参考方案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.13 秒（完全没有查询的代码，只有 Delaunay 结构，运行时间为 0.030 秒）。

【讨论】：

多么好的答案，非常感谢！实现线性搜索的想法正是我正在寻找的。虽然速度没有你机器上的那么快（0.39s -> 0.13s），但它仍然是一件很棒的事情。也许我会尝试用 Numba 或类似的工具编译新函数，看看是否有更多帮助。谢谢！