在3D阵列上快速内插 [英] Fast interpolation over 3D array

问题描述

我有一个3D数组，需要在一个轴(最后一个维度)上进行插值.假设y.shape = (nx, ny, nz)，我想为每个(nx, ny)插入nz.但是，我想在每个[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属性和低级函数_bspleval()中使用_fitpack模块.这是代码:

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屋！