NumPy：如何计算多轴上的分段线性插值

Question

给定以下 ndarray t -

In [26]: t.shape                                                                                     
Out[26]: (3, 3, 2)

In [27]: t                                                                                           
Out[27]: 
array([[[ 0,  1],
        [ 2,  3],
        [ 4,  5]],

       [[ 6,  7],
        [ 8,  9],
        [10, 11]],

       [[12, 13],
        [14, 15],
        [16, 17]]])

点 t[:, 0, 0] 的分段线性插值可以使用 numpy.interp -

对 [0 , 0.66666667, 1.33333333, 2.] 进行如下计算

In [38]: x = np.linspace(0, t.shape[0]-1, 4)                                                         

In [39]: x                                                                                           
Out[39]: array([0.        , 0.66666667, 1.33333333, 2.        ])

In [30]: xp = np.arange(t.shape[0])                                                                  

In [31]: xp                                                                                          
Out[31]: array([0, 1, 2])

In [32]: fp = t[:,0,0]                                                                               

In [33]: fp                                                                                          
Out[33]: array([ 0,  6, 12])

In [40]: np.interp(x, xp, fp)                                                                        
Out[40]: array([ 0.,  4.,  8., 12.])

对于 fp -

的所有值，如何有效计算并一起返回所有插值

array([[[ 0,  1],
        [ 2,  3],
        [ 4,  5]],

       [[ 4,  5],
        [ 6,  7],
        [ 8,  9]],

       [[ 8,  9],
        [10, 11],
        [12, 13]],

       [[12, 13],
        [14, 15],
        [16, 17]]])

Answer 1

你可以试试 scipy.interpolate.interp1d:

from scipy.interpolate import interp1d
import numpy as np

 t = np.array([[[ 0,  1],
                [ 2,  3],
                [ 4,  5]],

               [[ 6,  7],
                [ 8,  9],
                [10, 11]],

               [[12, 13],
                [14, 15],
                [16, 17]]])

# for the first slice
f = interp1d(np.arange(t.shape[0]), t[..., 0], axis=0)
# returns a function which you call with values within range np.arange(t.shape[0])

# data used for interpolation
t[..., 0]
>>> array([[ 0,  2,  4],
           [ 6,  8, 10],
           [12, 14, 16]])

f(1)
>>> array([ 6.,  8., 10.])

f(1.5)
>>> array([ 9., 11., 13.])

Answer 2

由于插值是 1d 并改变 y 值，因此对于 t 的每个 1d 切片必须是运行。显式循环可能更快，但使用 np.apply_along_axis

循环更整洁

import numpy as np

t = np.arange( 18 ).reshape(3,3,2)

x = np.linspace( 0, t.shape[0]-1, 4)
xp = np.arange(t.shape[0])


def interfunc( arr ):
    """ Function interpolates a 1d array. """
    return np.interp( x, xp, arr )

np.apply_along_axis( interfunc, 0, t ) # apply function along axis 0

"""  Result
array([[[ 0.,  1.],
        [ 2.,  3.],
        [ 4.,  5.]],

       [[ 4.,  5.],
        [ 6.,  7.],
        [ 8.,  9.]],

       [[ 8.,  9.],
        [10., 11.],
        [12., 13.]],

       [[12., 13.],
        [14., 15.],
        [16., 17.]]]) """

有显式循环

result = np.zeros((4,3,2))
for c in range(t.shape[1]):
    for p in range(t.shape[2]):
       result[:,c,p] = np.interp( x, xp, t[:,c,p])

在我的机器上，第二个选项运行只用了一半的时间。

编辑使用np.nditer

由于结果和参数具有不同的形状，我似乎必须创建两个 np.nditer 对象，一个用于参数，一个用于结果。这是我第一次尝试将 nditer 用于任何事情，所以它可能过于复杂。

def test( t ):                                                              
    ts = t.shape              
    result = np.zeros((ts[0]+1,ts[1],ts[2]))
    param = np.nditer( [t], ['external_loop'], ['readonly'],  order = 'F')
    with np.nditer( [result], ['external_loop'], ['writeonly'],  order = 'F') as res:       
        for p, r in zip( param, res ):
            r[:] = interfunc(p)
    return result

它比显式循环稍慢，并且比其他任何一种解决方案都更难理解。

Answer 3

根据@Tis Chris 的要求，这是一个使用 np.nditer 和 multi_index 标志的解决方案，但我更喜欢上面的显式嵌套 for 循环方法，因为它快 10%

In [29]: t = np.arange( 18 ).reshape(3,3,2)                                                          

In [30]: ax0old = np.arange(t.shape[0])                                                              

In [31]: ax0new = np.linspace(0, t.shape[0]-1, 4)                                                    

In [32]: tnew = np.zeros((len(ax0new), t.shape[1], t.shape[2]))                                      

In [33]: it = np.nditer(t[0], flags=['multi_index'])                                                 

In [34]: for _ in it: 
    ...:     tnew[:, it.multi_index[0], it.multi_index[1]] = np.interp(ax0new, ax0old, t[:, it.multi_
    ...: index[0], it.multi_index[1]]) 
    ...:                                                                                             

In [35]: tnew                                                                                        
Out[35]: 
array([[[ 0.,  1.],
        [ 2.,  3.],
        [ 4.,  5.]],

       [[ 4.,  5.],
        [ 6.,  7.],
        [ 8.,  9.]],

       [[ 8.,  9.],
        [10., 11.],
        [12., 13.]],

       [[12., 13.],
        [14., 15.],
        [16., 17.]]])

NumPy：如何计算多轴上的分段线性插值

NumPy: How to calulate piecewise linear interpolant on multiple axes

numpy

linear-interpolation