如何使用最近邻居对高维 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。我先在标签['wavelength'] 中找到相应的索引，然后用它来索引

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

我对其他维度做同样的事情，包括视角等它们对应于其他 *_iters

我需要做的是找到查找表中的值之间的值，如果我事先不知道查找表的尺寸，我需要它来工作。如果我这样做了，那么我知道如何使用每个维度的循环来解决问题。但是如果我不知道 LUT 的维度是多少，那么我就不知道要嵌套多少个循环。

我认为我应该能够使用 cKDTree 来做到这一点，但我无法弄清楚如何让它工作。我非常感谢一个看起来与我的结构相似的示例

谢谢

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 等于线性插值）使用它是没有意义的，但即使是三次样条，编写自己的插值算法也不是那么简单。

当然，如果您的网格在所有方向上都是等距的，那么代码会更简单，因为无需先插入正确的位置（简单的除法即可）。

Answer 2

scipy.interpolate.RegularGridInterpolator 可以很好地解决这个问题。虽然它仅在 Scipy 0.14（目前最新版本）中可用。

如果您在变量中有 *_iters，您可以这样做：

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 点数组），如果您的代码需要这样做，这可能会非常有效。