不规则网格间多重插值的加速比网格数据
不规则网格间多重插值的加速比网格数据
提问于 2014-01-04 01:02:51
我有几个值是在相同的不规则网格(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进行此操作，并且运行良好。但是，由于我必须分别进行每一次插值，而且有很多点，所以速度很慢，计算中有大量的重复(即找出最接近的点，设置网格等)。
有没有办法加速计算，减少重复计算？也就是说，按照定义两个网格的方式，然后改变插值值？
回答 4
Stack Overflow用户
回答已采纳
发布于 2014-01-05 06:41:39
每次打电话到scipy.interpolate.griddata时，都会发生几件事
- 首先，调用
sp.spatial.qhull.Delaunay对不规则网格坐标进行三角剖分。 - 然后，对于新网格中的每个点，搜索三角剖分，以找出在哪个三角形(实际上，在哪个单纯形中，在您的三维情况下，该三角形将位于哪个四面体中)。
- 计算了每个新网格点相对于封闭单纯形顶点的重心坐标。
- 使用重心坐标和封闭单纯形顶点处的函数值，计算该网格点的插值值。
对于所有的插值，前三个步骤是相同的，所以如果您能够为每一个新的网格点存储包围单纯形的顶点的索引和插值的权重，那么计算量就会减少很多。不幸的是，使用可用的功能很难直接做到这一点，尽管这确实是可能的：
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 retStack Overflow用户
发布于 2015-08-11 12:55:35
非常感谢Jaime的解决方案(即使我并不真正理解重心计算是如何完成的.)
在这里，您将发现一个从2D中他的案例中改编的例子：
import scipy.interpolate as spint
import scipy.spatial.qhull as qhull
import numpy as np
def interp_weights(xy, uv,d=2):
tri = qhull.Delaunay(xy)
simplex = tri.find_simplex(uv)
vertices = np.take(tri.simplices, simplex, axis=0)
temp = np.take(tri.transform, simplex, axis=0)
delta = uv - 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)
m, n = 101,201
mi, ni = 1001,2001
[Y,X]=np.meshgrid(np.linspace(0,1,n),np.linspace(0,2,m))
[Yi,Xi]=np.meshgrid(np.linspace(0,1,ni),np.linspace(0,2,mi))
xy=np.zeros([X.shape[0]*X.shape[1],2])
xy[:,0]=Y.flatten()
xy[:,1]=X.flatten()
uv=np.zeros([Xi.shape[0]*Xi.shape[1],2])
uv[:,0]=Yi.flatten()
uv[:,1]=Xi.flatten()
values=np.cos(2*X)*np.cos(2*Y)
#Computed once and for all !
vtx, wts = interp_weights(xy, uv)
valuesi=interpolate(values.flatten(), vtx, wts)
valuesi=valuesi.reshape(Xi.shape[0],Xi.shape[1])
print "interpolation error: ",np.mean(valuesi-np.cos(2*Xi)*np.cos(2*Yi))
print "interpolation uncertainty: ",np.std(valuesi-np.cos(2*Xi)*np.cos(2*Yi)) 它可以应用于图像转换，例如图像映射，并加快速度。
您不能使用相同的函数定义，因为新的坐标在每次迭代时都会发生变化，但是您可以为所有人计算一次三角剖分。
import scipy.interpolate as spint
import scipy.spatial.qhull as qhull
import numpy as np
import time
# Definition of the fast interpolation process. May be the Tirangulation process can be removed !!
def interp_tri(xy):
tri = qhull.Delaunay(xy)
return tri
def interpolate(values, tri,uv,d=2):
simplex = tri.find_simplex(uv)
vertices = np.take(tri.simplices, simplex, axis=0)
temp = np.take(tri.transform, simplex, axis=0)
delta = uv- temp[:, d]
bary = np.einsum('njk,nk->nj', temp[:, :d, :], delta)
return np.einsum('nj,nj->n', np.take(values, vertices), np.hstack((bary, 1.0 - bary.sum(axis=1, keepdims=True))))
m, n = 101,201
mi, ni = 101,201
[Y,X]=np.meshgrid(np.linspace(0,1,n),np.linspace(0,2,m))
[Yi,Xi]=np.meshgrid(np.linspace(0,1,ni),np.linspace(0,2,mi))
xy=np.zeros([X.shape[0]*X.shape[1],2])
xy[:,1]=Y.flatten()
xy[:,0]=X.flatten()
uv=np.zeros([Xi.shape[0]*Xi.shape[1],2])
# creation of a displacement field
uv[:,1]=0.5*Yi.flatten()+0.4
uv[:,0]=1.5*Xi.flatten()-0.7
values=np.zeros_like(X)
values[50:70,90:150]=100.
#Computed once and for all !
tri = interp_tri(xy)
t0=time.time()
for i in range(0,100):
values_interp_Qhull=interpolate(values.flatten(),tri,uv,2).reshape(Xi.shape[0],Xi.shape[1])
t_q=(time.time()-t0)/100
t0=time.time()
values_interp_griddata=spint.griddata(xy,values.flatten(),uv,fill_value=0).reshape(values.shape[0],values.shape[1])
t_g=time.time()-t0
print "Speed-up:", t_g/t_q
print "Mean error: ",(values_interp_Qhull-values_interp_griddata).mean()
print "Standard deviation: ",(values_interp_Qhull-values_interp_griddata).std()在我的笔记本电脑上，速度在20到40倍之间！
希望能帮到一个人
Stack Overflow用户
发布于 2019-06-17 09:55:19
我也有同样的问题(网格数据非常慢，网格在许多插值中保持不变)，我最喜欢解决方案在此描述，主要是因为它非常容易理解和应用。
它使用的是LinearNDInterpolator，可以通过只需要计算一次的Delaunay三角剖分。从该帖子复制和粘贴(所有学分到xdze2)：
from scipy.spatial import Delaunay
from scipy.interpolate import LinearNDInterpolator
tri = Delaunay(mesh1) # Compute the triangulation
# Perform the interpolation with the given values:
interpolator = LinearNDInterpolator(tri, values_mesh1)
values_mesh2 = interpolator(mesh2)这使我的计算速度提高了大约2倍。
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/20915502
复制相关文章
相似问题
腾讯云开发者
