Python中的规则网格插值

Question

我一直在努力插值我的2D矩阵中的“空”像素数据。基本上，我了解（但不深入）诸如反距离加权、克里金、双三次等插值技术。我不确定起点确切位置（无论是在问题陈述中还是在Python案例中）。

问题定义： 我有一个MxN矩阵（常规网格），其中每个像素表示某个测量值（下图和此图中使用的数据在这里）。我想要对“问号空间”（也由相同大小但为空的像素组成的白色空间）区域的数据进行插值，使用我已有的蓝色像素作为现有数据。

我的问题：

1）如何插值这些数据？谁能给我一个简单的例子（例如3x3矩阵），以便我可以清楚地理解它？

2）谁能指导我如何在Python环境中执行解决方案步骤？

3）如何使用Python比较插值技术的精度？

4）您认为根据数据密度使用不同的插值方法是一个好主意吗？

我会感激您的答案和建议。

- Spider

---

这里有一种方法可以做到：http://stackoverflow.com/questions/17115030/want-to-smooth-a-contour-from-a-masked-array/17125125#17125125 尽管我认为你的空洞太大了。 - theta

对于中间的小空洞，你可能能够做出一些明智的事情，但你已有的数据变化范围小于边缘大面积丢失的数据。如果你只是为了视觉相似性，我建议随机填充并进行平滑处理。 - Gabriel

这个问题可能在https://stats.stackexchange.com/上得到更好的回答。 - jb.

对于一种计算所有缺失值的原则性方法，请参见 Numerical Recipes 中的 Laplace 插值，第 150-153 页。 - denis

@Spider，你还能访问“measurements.txt”吗？Dropbox链接已经失效了，这使得复制这些示例变得困难。 - Robbi Bishop-Taylor

Answer 1

对于插值像素，一个合理的解决方案主要取决于你想用它来回答什么问题——请注意：在缺失数据上进行外推可能会导致非常误导人的答案！

径向基函数插值/核平滑

就可在Python中实现的实际解决方案而言，一种填充这些像素的方法是使用Scipy的径向基函数插值实现（请参见这里），该方法旨在对分散数据进行平滑/插值。

鉴于您的矩阵M和基础1D坐标数组r和c（使得M.shape == (r.size, c.size)），其中M的缺失条目设置为nan，这似乎可以使用线性RBF内核运行得相当好：

import numpy as np
import scipy.interpolate as interpolate

with open('measurement.txt') as fh:
    M = np.vstack(map(float, r.split(' ')) for r in fh.read().splitlines())
r = np.linspace(0, 1, M.shape[0]) 
c = np.linspace(0, 1, M.shape[1])

rr, cc = np.meshgrid(r, c)
vals = ~np.isnan(M)
f = interpolate.Rbf(rr[vals], cc[vals], M[vals], function='linear')
interpolated = f(rr, cc)

以下是您链接到的数据进行插值后的结果。虽然看起来合理，但它确实突出了缺失样本与真实数据比例不利的问题:

高斯过程回归/Kriging

Kriging插值可以通过scikit-learn库中的高斯过程回归实现来提供（它本身基于Matlab的DACE Kriging工具箱）。可以按如下方式调用：

from sklearn.gaussian_process import GaussianProcess

gp = GaussianProcess(theta0=0.1, thetaL=.001, thetaU=1., nugget=0.01)
gp.fit(X=np.column_stack([rr[vals],cc[vals]]), y=M[vals])
rr_cc_as_cols = np.column_stack([rr.flatten(), cc.flatten()])
interpolated = gp.predict(rr_cc_as_cols).reshape(M.shape)

这会产生与上面径向基函数示例非常相似的插值。在这两种情况下，有许多参数需要探索 - 这些选择在很大程度上取决于您可以对数据做出什么样的假设。（在上面的RBF示例中使用的线性内核的一个优点是它没有自由参数）

Inpainting

最后，一个完全视觉动机的解决方案是使用OpenCV的inpainting功能，尽管这假定为8位数组（0-255），并且没有简单的数学解释。

- rroowwllaanndd