在3D阵列上快速插值

Question

我有一个3D数组,我需要插入一个轴(最后一个维度).让我们说y.shape = (nx, ny, nz),我想插入nz每一个(nx, ny).但是,我想在每个中插入不同的值[i, j].

这里有一些代码可以举例说明.如果我想,内插一个值,说new_z,我会使用scipy.interpolate.interp1d这样的

# 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)

但是,对于这个问题,我真正想要的是new_z为每个插入不同的内容y[i, j].所以我这样做:

# 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中手动完成.

Answer 1

提速高阶插值,可以调用interp1d()一次,然后使用_spline属性和低级别功能_bspleval()中_fitpack的模块.这是代码:

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