python - 为两个不规则网格之间的多个插值加速 scipy griddata

Question

我在同一个不规则网格上定义了几个值 (x, y, z)我想插入一个新的网格(x1, y1, z1) .即，我有 f(x, y, z), g(x, y, z), h(x, y, z)我想计算f(x1, y1, z1), g(x1, y1, z1), h(x1, y1, z1) .

目前我正在使用 scipy.interpolate.griddata它运作良好。但是，因为我必须单独执行每个插值并且有很多点，所以速度很慢，计算中有大量重复(即找到最接近的点，设置网格等......)。

有没有办法加快计算速度，减少重复计算？即类似于定义两个网格，然后更改插值的值？

Answer 1

每次调用 scipy.interpolate.griddata 时都会发生几件事情:

1. 首先，调用 sp.spatial.qhull.Delaunay 对不规则网格坐标进行三角测量。
2. 然后，对于新网格中的每个点，搜索三角剖分以找到它位于哪个三角形(实际上是在哪个单纯形中，在您的 3D 情况下是哪个四面体)。
3. 计算每个新网格点相对于封闭单纯形顶点的重心坐标。
4. 使用重心坐标和封闭单纯形顶点处的函数值计算该网格点的插值。

前三个步骤对于所有插值都是相同的，因此如果您可以为每个新网格点存储封闭单纯形的顶点索引和插值的权重，则可以通过以下方式最小化计算量很多。不幸的是，这并不容易直接使用可用的功能来完成，尽管这确实是可能的:

import scipy.interpolate as spint
import scipy.spatial.qhull as qhull
import itertools

def interp_weights(xyz, uvw):
    tri = qhull.Delaunay(xyz)
    simplex = tri.find_simplex(uvw)
    vertices = np.take(tri.simplices, simplex, axis=0)
    temp = np.take(tri.transform, simplex, axis=0)
    delta = uvw - temp[:, d]
    bary = np.einsum('njk,nk->nj', temp[:, :d, :], delta)
    return vertices, np.hstack((bary, 1 - bary.sum(axis=1, keepdims=True)))

def interpolate(values, vtx, wts):
    return np.einsum('nj,nj->n', np.take(values, vtx), wts)

函数 interp_weights 为我上面列出的前三个步骤进行计算。然后函数 interpolate 使用这些计算值非常快地执行第 4 步:

m, n, d = 3.5e4, 3e3, 3
# make sure no new grid point is extrapolated
bounding_cube = np.array(list(itertools.product([0, 1], repeat=d)))
xyz = np.vstack((bounding_cube,
                 np.random.rand(m - len(bounding_cube), d)))
f = np.random.rand(m)
g = np.random.rand(m)
uvw = np.random.rand(n, d)

In [2]: vtx, wts = interp_weights(xyz, uvw)

In [3]: np.allclose(interpolate(f, vtx, wts), spint.griddata(xyz, f, uvw))
Out[3]: True

In [4]: %timeit spint.griddata(xyz, f, uvw)
1 loops, best of 3: 2.81 s per loop

In [5]: %timeit interp_weights(xyz, uvw)
1 loops, best of 3: 2.79 s per loop

In [6]: %timeit interpolate(f, vtx, wts)
10000 loops, best of 3: 66.4 us per loop

In [7]: %timeit interpolate(g, vtx, wts)
10000 loops, best of 3: 67 us per loop

所以首先，它和 griddata 一样，这很好。其次，设置插值，即计算 vtx 和 wts 与调用 griddata 大致相同。但第三，您现在几乎可以立即在同一个网格上插入不同的值。

griddata 唯一没有在这里考虑的是将 fill_value 分配给必须外推的点。您可以通过检查至少有一个权重为负的点来做到这一点，例如:

def interpolate(values, vtx, wts, fill_value=np.nan):
    ret = np.einsum('nj,nj->n', np.take(values, vtx), wts)
    ret[np.any(wts < 0, axis=1)] = fill_value
    return ret