3D 数组上的快速插值 [英] Fast interpolation over 3D array
问题描述
我有一个 3D 数组,需要在一个轴(最后一个维度)上进行插值.假设 y.shape = (nx, ny, nz)
,我想在 nz
中为每个 (nx, ny)
进行插值.但是,我想在每个 [i, j]
中插入不同的值.
I have a 3D array that I need to interpolate over one axis (the last dimension). Let's say y.shape = (nx, ny, nz)
, I want to interpolate in nz
for every (nx, ny)
. However, I want to interpolate for a different value in each [i, j]
.
这是一些示例代码.如果我想插入单个值,比如 new_z
,我会像这样使用 scipy.interpolate.interp1d
Here's some code to exemplify. If I wanted to interpolate to a single value, say new_z
, I'd use scipy.interpolate.interp1d
like this
# y is a 3D ndarray
# x is a 1D ndarray with the abcissa values
# new_z is a number
f = scipy.interpolate.interp1d(x, y, axis=-1, kind='linear')
result = f(new_z)
然而,对于这个问题,我真正想要的是为每个 y[i, j]
插入一个不同的 new_z
.所以我这样做:
However, for this problem what I actually want is to interpolate to a different new_z
for each y[i, j]
. So I do this:
# y is a 3D ndarray
# x is a 1D ndarray with the abcissa values
# new_z is a 2D array
result = numpy.empty(y.shape[:-1])
for i in range(nx):
for j in range(ny):
f = scipy.interpolate.interp1d(x, y[i, j], axis=-1, kind='linear')
result[i, j] = f(new_z[i, j])
不幸的是,如果有多个循环,这会变得低效且缓慢.有没有更好的方法来做这种插值?线性插值就足够了.一种可能性是在 Cython 中实现这一点,但我试图避免这种情况,因为我希望能够灵活地更改为三次插值,并且不想在 Cython 中手动完成.
Unfortunately, with multiple loops this becomes inefficient and slow. Is there a better way to do this kind of interpolation? Linear interpolation is sufficient. A possibility is to implement this in Cython, but I was trying to avoid that because I want to have the flexibility of changing to cubic interpolation and don't want to do it by hand in Cython.
推荐答案
为了加速高阶插值,你可以只调用一次 interp1d()
,然后使用 _spline
> 属性和 _fitpack
模块中的低级函数 _bspleval()
.代码如下:
To speedup high order interpolate, you can call interp1d()
only once, and then use the _spline
attribute and the low level function _bspleval()
in the _fitpack
module. Here is the code:
from scipy.interpolate import interp1d
import numpy as np
nx, ny, nz = 30, 40, 50
x = np.arange(0, nz, 1.0)
y = np.random.randn(nx, ny, nz)
new_x = np.random.random_integers(1, (nz-1)*10, size=(nx, ny))/10.0
def original_interpolation(x, y, new_x):
result = np.empty(y.shape[:-1])
for i in xrange(nx):
for j in xrange(ny):
f = interp1d(x, y[i, j], axis=-1, kind=3)
result[i, j] = f(new_x[i, j])
return result
def fast_interpolation(x, y, new_x):
from scipy.interpolate._fitpack import _bspleval
f = interp1d(x, y, axis=-1, kind=3)
xj,cvals,k = f._spline
result = np.empty_like(new_x)
for (i, j), value in np.ndenumerate(new_x):
result[i, j] = _bspleval(value, x, cvals[:, i, j], k, 0)
return result
r1 = original_interpolation(x, y, new_x)
r2 = fast_interpolation(x, y, new_x)
>>> np.allclose(r1, r2)
True
%timeit original_interpolation(x, y, new_x)
%timeit fast_interpolation(x, y, new_x)
1 loops, best of 3: 3.78 s per loop
100 loops, best of 3: 15.4 ms per loop
这篇关于3D 数组上的快速插值的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!