基于网格的多元立方插值

Question

假设我对以下内容感兴趣：

y = f(n, m, o)

我为 n、m 和 o 创建了离散的网格：

import numpy as np
nGrid = np.arange(0,10)
mGrid = np.arange(100,110)
oGrid = np.arange(1000,1010)
N,M,O = np.meshgrid(nGrid, mGrid, oGrid, indexing='ij')

为了简单起见，我们假设 f(n,m,o) = n+m+p：

Y = N+M+O

所以，现在我在 Y 中有了 f 的网格近似值。从我查阅的 文档来看，

• 我需要关注多变量插值
• scipy.interpolate.Griddata 看起来不错，但它是为 非结构化 网格设计的。这意味着我得把所有结果转换成长数组——这样就没法利用我其实有结构化网格的优势了。
• 在 结构化 下列出的其他方法不支持同时处理三次和多变量（超过2维）的输入。

考虑到我已经创建的变量，怎样做插值比较好呢？

更新

这是第一个答案的实现：

import numpy as np
nGrid = np.arange(0,10)
mGrid = np.arange(100,110)
oGrid = np.arange(1000,1010)
N,M,O = np.meshgrid(nGrid, mGrid, oGrid, indexing='ij')
Y = N+M+O

n,m,o = 5, 103, 1007

m_i = np.argmin(np.abs(np.floor(m)-mGrid))
m_f = m - m_i

n_i = np.argmin(np.abs(np.floor(n)-nGrid))
n_f = n - n_i

o_i = np.argmin(np.abs(np.floor(o)-oGrid))
o_f = o - o_i

A = Y[m_i-1:m_i+3, n_i-1:n_i+3, o_i-1:o_i+3]

# cubic smoothing polynome in 1D
# Catmull-Rom style
def v(p0, p1, p2, p3, f):
    return 0.5 * (p0 * (-f**3 + 2*f**2 - f) +
                  p1 * (3*f**3 - 5*f**2 + 2) +
                  p2 * (-3*f**3 + 4*f**2 + f) +
                  p3 * (f**3 - f**2))



B = v(A[0], A[1], A[2], A[3], m_f)
C = v(B[0], B[1], B[2], B[3], n_f)
D = v(C[0], C[1], C[2], C[3], o_f)

# D is the final interpolated array

print m+n+o, D

不幸的是，f(n,m,o) = 1115，而 D = 2215。由于没有方法论的链接，我很难理解到底发生了什么，以及为什么近似值差得这么远。

Answer 1

我找到的一种方法是 ndimage.map_coordinates。虽然它有点隐藏，没有在插值模块中链接，但效果非常好，并且提供了多种近似的选择。

不过有两个主要缺点：

• 在网格的边缘处理得有点奇怪，使用 mode='nearest' 时，它会选择边界的值，而如果不这样做，它会假设网格外的点值为 0。
• 效率不是很好。

我仍然欢迎更好的建议，这只是我目前的情况汇报——也许这能在某个时候帮助到别人。

import numpy as np
nGrid = np.arange(0,10)
mGrid = np.arange(100,110)
oGrid = np.arange(1000,1010)
N,M,O = np.meshgrid(nGrid, mGrid, oGrid, indexing='ij')
Y = N+M+O
n,m,o = 0.01, 103, 1007
from scipy import ndimage
data = Y
# subtract lower end of grid
coords = np.array([[n-0,m-100,o-1000]])
ndimage.map_coordinates(data, coords.T, order=2, mode='nearest')

Answer 2

如果你找不到合适的方法，那就自己动手做吧！

不过要注意，三次立方插值有点复杂。下面会介绍一种简单的三次立方插值方法。其实还有很多更复杂、更快的方法（基于矩阵的），具体要看你的需求。另外，要知道有几种立方外推方法，我这里用的是最常见的那种（Catmull-Rom）。

在一个三维数组 M 中查找一个点 (m,n,o) 时，进行三次立方插值需要查看周围的一个 4x4x4 点（或者 3x3x3 单元）网格。每个轴上的点是通过以下方式定义的：floor(a)-1, floor(a), floor(a)+1, floor(a)+2 等等。（可以想象成一个魔方，点 (m,n,o) 在看不见的中心小方块里。）

三次立方插值的结果是这 64 个点的加权平均值，权重取决于该点与网格点的距离。

我们来定义一下：

m_i = np.floor(m).astype('int')    # integer part
m_f = m - m_i                      # fraction

对于 n 和 o 也有类似的定义。现在我们需要对数组进行操作：

A = M[m_i-1:m_i+3, n_i-1:n_i+3, o_i-1:o_i+3]

我们的点现在在 (1+m_f, 1+n_f, 1+o_f) 这个位置。（注意，这里假设矩阵是无限的。否则需要考虑边缘效应。）

这 64 个系数可以通过一些巧妙的矩阵运算来计算，但这里只需要知道插值是可结合的，也就是说我们可以逐轴进行计算：

# cubic smoothing polynome in 1D
# Catmull-Rom style
def v(p0, p1, p2, p3, f):
    return 0.5 * (p0 * (-f**3 + 2*f**2 - f) +
                  p1 * (3*f**3 - 5*f**2 + 2) +
                  p2 * (-3*f**3 + 4*f**2 + f) +
                  p3 * (f**3 - f**2))


# interpolate axis-by-axis
B = v(A[0], A[1], A[2], A[3], m_f)
C = v(B[0], B[1], B[2], B[3], n_f)
D = v(C[0], C[1], C[2], C[3], o_f)

# D is the final interpolated value

因此，插值是逐轴计算的。

---

现在让我们把这个变得更实用一些。首先，我们需要从自己的坐标空间计算到数组坐标的地址。如果 m 轴被定义为：

m0: position of the first grid plane on the axis
m1: position of the last grid plane on the axis
n_m: number of samples

我们可以从“真实世界”的位置 m 计算出等效的数组位置：

m_arr = (n_m - 1) * (m -  m0) / (m1 - m0)

这里的 -1 是因为“栅栏现象”（如果我们有 10 个样本，第一和最后之间的距离是 9）。

现在我们可以把所有东西放在一个函数里：

# cubic smoothing polynome in 1D
# Catmull-Rom style
def interp_catmull_rom(p0, p1, p2, p3, f):
    return 0.5 * (p0 * (-f**3 + 2*f**2 - f) +
                  p1 * (3*f**3 - 5*f**2 + 2) +
                  p2 * (-3*f**3 + 4*f**2 + f) +
                  p3 * (f**3 - f**2))

# linear interpolation
# only needs p1 and p2, the rest are there for compatibility
def interp_linear(p0, p1, p2, p3, f):
    return (1-f)*p1 + f*p2


# don't interpolate, use the nearest point
# only needs p1 and p2, the rest are there for compatibility
def interp_nearest(p0, p1, p2, p3, f):
    if f > .5:
        return p2
    else:
        return p1

# 3D interpolation done axis-by-axis
def tri_interp(M, f, m, n, o, m0, m1, n0, n1, o0, o1):
    # M: 3D array
    # f: interpolation function to use
    # m,n,o: coordinates where to interpolate
    # m0: real world minimum for m
    # m1: real world maximum for m
    # n0,n1,o0,o0: as with m0 and m1

    # calculate the array coordinates
    m_arr = (M.shape[0] - 1) * (m - m0) / (m1 - m0)
    n_arr = (M.shape[1] - 1) * (n - n0) / (n1 - n0)
    o_arr = (M.shape[2] - 1) * (o - o0) / (o1 - o0)

    # if we are out of our data box, return a nan
    if m_arr < 0 or m_arr > M.shape[0] or \
       n_arr < 0 or n_arr > M.shape[1] or \
       o_arr < 0 or o_arr > M.shape[2]:
        return np.nan

    # calculate the integer parts and fractions
    m_i = np.floor(m_arr).astype('int')
    n_i = np.floor(n_arr).astype('int')
    o_i = np.floor(o_arr).astype('int')
    m_f = m_arr - m_i
    n_f = n_arr - n_i
    o_f = o_arr - o_i

    # edge effects may be nasty, we may need elements outside of the array
    # there may be more efficient ways to avoid it, but we'll create lists of
    # coordinates:

    n_coords, m_coords, o_coords = np.mgrid[m_i-1:m_i+3, n_i-1:n_i+3, o_i-1:o_i+3]

    # these coordinate arrays are clipped so that we are in the available data
    # for example, point (-1,3,7) will use the point (0,3,7) instead
    m_coords = m_coords.clip(0, M.shape[0]-1)
    n_coords = n_coords.clip(0, M.shape[1]-1)
    o_coords = o_coords.clip(0, M.shape[2]-1)

    # get the 4x4x4 cube:
    A = M[m_coords, n_coords, o_coords]

    # interpolate along the first axis (3D to 2D)
    B = f(A[0], A[1], A[2], A[3], m_f)

    # interpolate along the second axis (2D to 1D)
    C = f(B[0], B[1], B[2], B[3], n_f)

    # interpolate along the third axis (1D to scalar)
    D = f(C[0], C[1], C[2], C[3], o_f)

    return D

这个函数可以让我们使用任何可以逐轴处理四个邻近点的插值方法。现在展示了三次插值、线性插值和“最近邻”插值。

---

那么这个方法有效吗？我们来通过生成一些随机数据和一条穿过这些数据的平面来测试一下。（后者出乎意料地困难。）

# random data
M = np.random.random((10, 12, 16))
m0,m1 = -10.,10.
n0,n1 = -7.,5.
o0,o1 = -4.,5.

# create a grid (grid from -15..15,-15..15 in n-m plane)
gr = mgrid[-15:15.01:.1, -15:15.01:.1]

# create two perpendicular unit vectors (forming a plane)
# vn: normal vector of the plane
# vp: some vector, whose projection determines one direction
# v0: unit vector on the plane (perpendicular to vn and vp)
# v1: unit vector on the plane (perpendicular to vn and v0)
vn = np.array([-.2, .3, 1])
vp = np.array([0, -1, 0])
v1 = np.cross(vn, vp)
v2 = np.cross(vn, v1)
v1 /= numpy.linalg.norm(v1)
v2 /= numpy.linalg.norm(v2)

# the grid and the unit vectors define the 3d points on the plane
gr3d = gr[0][:,:,None] * v1 + gr[1][:,:,None] * v2

# now we can fetch the points at grid points, just must flatten and reshape back
res = [ tri_interp(M, interp_catmull_rom, p[0], p[1], p[2], m0,m1,n0,n1,o0,o1) for p in gr3d.reshape(-1,3) ]
res = np.array(res).reshape((gr3d.shape[0], gr3d.shape[1]))

所以，我们有一条平面穿过我们的数据块。现在我们可以通过选择 interp_linear, interp_nearest 或 interp_catmull_rom 作为一维插值函数，实际看到三种不同插值方法之间的差异。下面的所有图片在尺寸和颜色上都是相同的比例。

1. 最近邻：取网格中最近的已知点（简单粗暴）

2. 线性：取一个 2x2x2 的立方体，使用线性插值计算结果；结果总是在已有数据点之间，不需要边缘处理。不过，第一次导数不光滑，这在图像中可以看到，有时在信号处理上会是个问题。

3. 三次：取一个 4x4x4 的立方体，使用三次插值计算结果。结果的第一次导数是连续的，也就是说结果在所有方向上都是“平滑”的。由于样本体积较大，可能会有边缘效应。（插值算法会重复最边缘的体素，但它们可能会被镜像处理、作为某个常数处理等。）

---

如果我们通过数据块取一条线段，可能会更清楚地看到这些插值方法的差异：

 # line segment from (-12,-12,-12) to (12,12,12)
 v = np.array([np.linspace(-12,12,1000)]*3).T

 res = [ tri_interp(M, interp_catmull_rom, p[0], p[1], p[2], m0,m1,n0,n1,o0,o1) for p in v ]

对于不同的插值函数，我们在这条线段上的值如下：

“最近邻”的块状和线性插值中的折线都很明显。可能会让人惊讶的是，“线性”插值并没有给出线性的结果。理解这一点可以通过想象一个 2D 插值，四个角的值分别是 1,0,1,0（即对角的两个角有相同的值）。没有办法插值使得所有部分都是线性的。

---

这里展示的算法是一个非常简单的算法，可以通过一些技巧加快速度。特别是如果在一个单元中需要多个点，基于矩阵的方法会快得多。不幸的是，没有非常快的方法，而基于矩阵的方法解释起来也更复杂。