使用scipy.interpolate.LSQBivariateSplines拟合带噪声和缺失的2D样条

Question

我有一个numpy数组，里面是一个网格上的矩形数据，我想给它拟合一个二维样条曲线，这样可以保留大范围的变化，同时去掉大部分噪声。数据中还有一些区域被标记为无效，值是NaN。

我试着用scipy.interpolate.RectBivariateSpline这个函数，但因为有空缺的地方，结果就不太好。所以我又尝试了同一个包里的LSQBivariateSpline函数，希望把所有NaN像素的权重设置为0，这样它就能忽略这些无效值。然而，这时我遇到了一个我不知道怎么解决的错误：

我的代码是：

# some preparation, loading data and stuff
# all my data is stored in 'data'

# Create the knots (10 knots in each direction, making 100 total
xcoord = numpy.linspace(5, data.shape[0]-5, 10)
ycoord = numpy.linspace(5, data.shape[1]-5, 10)

# Create all weights, and set them to 0 when the data is NaN
weights = numpy.ones(data.shape)
weights[numpy.isnan(data)] = 1e-15  # weights must be >0

# LSQBivariateSpline needs x and y coordinates as 1-D arrays
x, y = numpy.indices(data.shape)
spline_fit = scipy.interpolate.LSQBivariateSpline(x.ravel(), y.ravel(), data.ravel(), 
                                               xcoord, ycoord, 
                                               w=weights.ravel(),
                                               bbox=[None, None, None, None], 
                                               kx=2, ky=2)

运行代码后出现的错误信息是：

The coefficients of the spline returned have been computed as the
minimal norm least-squares solution of a (numerically) rank deficient
system (deficiency=25). If deficiency is large, the results may be
inaccurate. Deficiency may strongly depend on the value of eps.
done!
Traceback (most recent call last):
  File "./fitpg.py", line 513, in <module>
    fit_pupilghost(prebinned, (center_y, center_x), (r_inner, r_outer), dr)
  File "./fitpg.py", line 298, in fit_pupilghost
    fitout = pupil2d(radius[:,y], angle[:,y])
  File "/usr/local/lib64/python2.7/site-packages/scipy/interpolate/fitpack2.py", line 545, in __call__
    raise ValueError("Error code returned by bispev: %s" % ier)
ValueError: Error code returned by bispev: 10

我输入的矩阵（'data'）大约是1000 x 1000像素，这应该足够用来约束100个节点。把每个方向的节点数增加到100会让代码运行得慢很多，但除了缺陷数量外，其他都没有变化。我还尝试过把eps值调高调低，范围在1-e30到0.9之间（默认是1e-16）。

我也试着在网上搜索这个错误代码，但没有找到有用的信息，所以这也没什么帮助。

你觉得这里可能出什么问题吗？或者有没有什么解决方法或更好的办法来处理这个问题？

任何帮助都非常感谢。

谢谢

Answer 1

这个样条拟合的代码对NaN（不是一个数字）没有特别的处理。因为任何与NaN接触的数字都会变成NaN，这就意味着NaN的存在会影响整个计算，导致你得不到任何结果。

你可以做的是，把NaN值替换成一些（随便的）有限值，同时设置一个零权重，比如说：

weights = numpy.ones(data.shape)
mask = numpy.isnan(data)
weights[mask] = 0
data[mask] = 0   # arbitrary

由于权重很小，选择的值其实没什么关系。如果代码对零权重有问题，你也可以尝试把对应的权重设置成一个小值，比如你的 1e-15。