每次您拨打电话时都会发生几件事scipy.interpolate.griddata
:
- 首先,致电
sp.spatial.qhull.Delaunay
用于对不规则网格坐标进行三角测量。
- 然后,对于新网格中的每个点,搜索三角剖分以查找它位于哪个三角形(实际上,在哪个单纯形中,在您的 3D 情况下将在哪个四面体中)。
- 计算每个新网格点相对于封闭单纯形的顶点的重心坐标。
- 使用重心坐标和封闭单纯形顶点处的函数值来计算该网格点的插值。
前三个步骤对于所有插值都是相同的,因此如果您可以为每个新网格点存储封闭单纯形的顶点索引和插值权重,则可以大大减少计算量。不幸的是,这并不容易直接使用可用的功能来完成,尽管它确实是可能的:
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
and 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