Scipy中的二维插值问题

Question

在我的应用程序中，数据是在一个扭曲的网格上采样的，我想把它重新采样到一个不扭曲的网格上。为了测试这个，我写了一个程序，里面有一些示例扭曲和一个简单的函数作为数据：

from __future__ import division

import numpy as np
import scipy.interpolate as intp
import pylab as plt

# Defining some variables:

quadratic = -3/128
linear = 1/16
pn = np.poly1d([quadratic, linear,0])

pixels_x = 50
pixels_y = 30
frame = np.zeros((pixels_x,pixels_y))

x_width= np.concatenate((np.linspace(8,7.8,57) , np.linspace(7.8,8,pixels_y-57)))

def data(x,y):
    z = y*(np.exp(-(x-5)**2/3) + np.exp(-(x)**2/5) + np.exp(-(x+5)**2))
    return(z)

# Generating grid coordinates

yt = np.arange(380,380+pixels_y*4,4)
xt = np.linspace(-7.8,7.8,pixels_x)
X, Y = np.meshgrid(xt,yt)
Y=Y.T
X=X.T

Y_m = np.zeros((pixels_x,pixels_y))
X_m = np.zeros((pixels_x,pixels_y))

# generating distorted grid coordinates:    

for i in range(pixels_y):
    Y_m[:,i] = Y[:,i] - pn(xt)
    X_m[:,i] = np.linspace(-x_width[i],x_width[i],pixels_x)


# Sample data:
for i in range(pixels_y):
    for j in range(pixels_x):
        frame[j,i] = data(X_m[j,i],Y_m[j,i])


Y_m = Y_m.flatten()
X_m = X_m.flatten()
frame = frame.flatten()
##
Y = Y.flatten()
X = X.flatten()
ipf = intp.interp2d(X_m,Y_m,frame)
interpolated_frame = ipf(xt,yt)

在这个时候，我有两个问题：

1. 代码是可以运行的，但我收到了以下警告：

   警告：不能再添加更多的节点，因为B样条系数的数量已经超过了数据点的数量m。可能的原因是：s或m太小了。(fp>s) kx,ky=1,1 nx,ny=54,31 m=1500 fp=0.000006 s=0.000000

另外，一些插值伪影出现了，我猜这些和警告有关——你们知道我哪里做错了吗？

1. 对于我的实际应用，框架需要大约500*100，但这样做时，我遇到了内存错误——除了把框架分成几个部分，还有什么办法可以解决这个问题吗？

谢谢！

Answer 1

你可能想看看basemap里的这个interp方法：

mpl_toolkits.basemap.interp http://matplotlib.sourceforge.net/basemap/doc/html/api/basemap_api.html

除非你真的需要基于样条的插值方法。

Answer 2

对于二维插值，griddata 是一个很不错的选择，它的特点是简单、快速，而且适合局部数据处理。你可以看看在StackOverflow上关于二维插值的问题。

Answer 3

这个问题很可能和在interp2d中使用的bisplrep和bisplev有关。文档提到它们使用了一个平滑因子s=0.0，如果你想对s有更多的控制，可以直接使用bisplrep和bisplev。相关的文档提到，s的值应该在(m-sqrt(2*m),m+sqrt(2*m))这个范围内，其中m是用来构建样条曲线的点的数量。我遇到过类似的问题，发现直接使用bisplrep和bisplev后问题解决了，这时s只是一个可选项。