如何在高维numpy python数组中使用最近邻插值法进行插值

Question

我正在使用scipy和numpy编写python程序，我有一个数据查找表（LUT），我可以像这样访问它：

self.lut_data[n_iter][m_iter][l_iter][k_iter][j_iter][i_iter] 

我获取*_iter索引对应的值数组是通过一个字典实现的。例如，i_iter索引对应的是光波长，所以我有一个标签和值的字典可以这样获取：

labels['wavelength']

该函数将返回一个数组，其中每个i_iter对应一种波长。如果我想要直接查询500 nm处的lut_data，那么先找到labels['wavelength']中对应的索引，然后使用该索引进行索引。

lut_data[][][][][][wavelength_index]

我还要对其他维度进行同样的操作，包括观察角度等，它们对应于其他*_iters

我需要做的是在查找表中找到介于其值之间的数值，并且如果事先不知道查找表的维度，也需要让此方法可行。如果我已经知道了LUT的维度，则可以使用循环来解决每个维度的问题。但是，如果我不知道LUT有多少维度，那么我就不知道需要嵌套多少循环。

我认为我应该能够使用cKDTree来解决此问题，但我无法理解如何使其运作。希望您提供一个类似于我的结构的示例。

谢谢

- Caustic

---

我不是很理解这个数据结构，它是一个Numpy数组吗？但我认为你应该看一下scipy.interpolate.NearestNDInterpolator。它基于cKDTree... - user2379410

我看了一下，谢谢。但是我无法让它工作。是的，我的数据结构是一个numpy数组。文档不如其他地方好。Npoints是什么，是点还是点的数量还是点的维数？它需要知道每个维度的大小吗？Ndims也是同样的问题吗？我要将点作为元组传递吗？Ndim是一个nD数组data[x，y，z，k，j，l]，其中k，j，l比常见的3D xyz更高阶吗？关于值的问题也是同样的情况。我认为一个经过充分研究的（高于2D，最好是高于3D）会帮助我理解它。谢谢。 - Caustic

现在想一想，我不知道为什么要使用 [x][y][z] 的索引方式而不是 [x,y,z]。我正在使用继承来的代码并且一直沿用这种写法。这样做会有什么区别吗？LUT 是使用 numpy.zeros(<dimensions>) 初始化的，并且似乎使用任一种表示法都能得到相同的结果。抱歉问这些愚蠢的问题，我不是职业程序员，更像一个黑客。 - Caustic

好的，我明白了 ;) 但是仔细想想，我觉得使用 scipy.interpolate.RegularGridInterpolator 更好。使用 values 来作为您的 LUT，points 则是各自维度上的 np.arange 的元组（每个维度的大小）。或者...更好的方法是使用您已经有的“标签”。如果您愿意，我可以给出一个例子。 - user2379410

Answer 1

如果您有完整的信息数组进行插值，线性插值并不难。只是稍微耗费一些时间，但如果您可以将数组适配到RAM中，那么只需要几秒钟的时间即可完成。

诀窍在于线性插值可以逐个轴进行。因此，对于每个轴：

- 找到要进行插值的最近点 - 找到这些点之间的相对距离（d = 0..1），例如，如果您有540和550 nm，并且您想在548 nm处拥有数据，则 d = 0.8。 - 对所有轴重复此过程；每一轮都会将维度数量减少一个。

就像这样：

import numpy as np

def ndim_interp(A, ranges, p):
    # A: array with n dimensions
    # ranges: list of n lists or numpy arrays of values along each dimension
    # p: vector of values to find (n elements)

    # iterate through all dimensions
    for i in range(A.ndim):
        # check if we are overrange; if we are, use the edgemost values
        if p[i] <= ranges[i][0]:
            A = A[0]
            continue
        if p[i] >= ranges[i][-1]:
            A = A[-1]
            continue

        # find the nearest values
        right = np.searchsorted(ranges[i], p[i])
        left = right - 1

        # find the relative distance
        d = (p[i] - ranges[i][left]) / (ranges[i][right] - ranges[i][left])

        # calculate the interpolation
        A = (1 - d) * A[left] + d * A[right]            

    return A

作为一个例子：

# data axis points
arng = [1, 2, 3]
brng = [100, 200]
crng = [540, 550, 560]

# some data
A = np.array([
    [[1., 2., 3.], [2., 3., 4.]],
    [[0.5, 1.5, 2.], [1.5, 2.0, 3.0]],
    [[0., 0.5, 1.], [1., 1., 1.]]])

# lookup:
print ndim_interp(A, (arng, brng, crng), (2.3, 130., 542.))

---

如果您想要做更复杂的事情（立方样条等），则可以使用 scipy.ndimage.interpolation.map_coordinates。那么配方如下所示：

import numpy as np
import scipy.ndimage.interpolation

def ndim_interp(A, ranges, p):
    # A: array with n dimensions
    # ranges: list of n lists or numpy arrays of values along each dimension
    # p: vector of values to find (n elements)

    # calculate the coordinates into array positions in each direction
    p_arr = []
    # iterate through all dimensions
    for i in range(A.ndim):
        # check if we are overrange; if we are, use the edgemost values
        if p[i] <= ranges[i][0]:
            p_arr.append(0)
            continue
        if p[i] >= ranges[i][-1]:
            p_arr.append(A.shape[i] - 1)
            continue

        # find the nearest values to the left
        right = np.searchsorted(ranges[i], p[i])
        left = right - 1

        # find the relative distance
        d = (p[i] - ranges[i][left]) / (ranges[i][right] - ranges[i][left])

        # append the position
        p_arr.append(left + d)

    coords = np.array(p_arr).reshape(A.ndim, -1)
    return scipy.ndimage.interpolation.map_coordinates(A, coords, order=1, mode='nearest')[0]

当然，如果使用最简单的设置(order=1等于线性插值)，那么使用它是没有意义的，但即使是采用三次样条插值，编写自己的插值算法也并不简单。 当然，如果你的网格在所有方向上都是等距的，那么代码会更简单，因为不需要先进行正确位置的插值(一个简单的除法就可以了)。

- DrV

---

非常感谢，这非常有帮助。我不得不添加一个额外的步骤来压缩只有一个值的维度，但这似乎有效。好吧，它给我的值略微不同于下面的方法。我需要调查原因。 - Caustic

Answer 2

对于这个问题，scipy.interpolate.RegularGridInterpolator是非常好的选择。虽然它只在Scipy 0.14中可用（截至目前最新版本）。

如果你已经将*_iter存储在变量中，你可以执行以下操作：

from scipy.interpolate import RegularGridInterpolator

points = tuple([n_iter, m_iter, l_iter, k_iter, j_iter, i_iter])
interpolator = RegularGridInterpolator(points, lut_data, method='nearest')

或者你可以从你的字典中获取points：

keys = ['k1', 'k2', 'k3', 'k4', 'k5', 'wavelength']
points = tuple([labels[key] for key in keys])

如果你拥有插值器，那么你可以使用其__call__方法来进行插值。这基本上意味着你可以将你创建的类实例作为函数调用：

point_of interest = tuple([x1, x2, x3, x4, x5, some_wavelength])
interp_value = interpolator(point_of_interest)

插值器还允许同时插值多个值（即传递一个Numpy点数组），如果您的代码需要这样做，则可能会显着提高效率。

- user2379410

---

太棒了，这个看起来也有效，只是和上面的方法相比得到了稍微不同的答案。RegularGridInterpolator是否期望均匀间隔的值？我认为并不是所有的值都是这样的。我的叶绿素值是如此增加的：0.01、0.1、0.5、1... - Caustic

@Caustic，确实你的数据间距不均匀，但这没有问题。你尝试过使用 method='linear' 还是保持为 'nearest' 呢？另一个答案使用了线性插值，这可能解释了差异。 - user2379410

@Caustic，哦，现在我想我明白了你最初的想法。 你想找到最近的邻居，然后对它们进行一些插值，无论是线性的还是其他任何种类？ - user2379410

是的，你们两个回复都是正确的。谢谢你们的帮助。 - Caustic