如何使用scipy获得非平滑2D样条插值 [英] How to get a non-smoothing 2D spline interpolation with scipy

问题描述

我希望2D三次样条曲线适合一些不规则间隔的数据-即一个函数，该函数完全适合给定点的数据-但也可以返回介于两者之间的值.

对于不规则间隔的数据，我只能找到 scipy.interpolate.SmoothBivariateSpline .我不知道如何关闭平滑"功能(无论我在 s 参数中输入什么值.

但是，我确实发现，使用 scipy.interpolate.griddata 我基本上可以得到我想要的东西-尽管每次都必须重新计算(即不只是生成一个函数).从根本上说，这两者之间是否有任何区别-即 griddata 在做的事情与样条曲线"不同吗?无论如何，是否可以在 SmoothBivariateSpline 中关闭平滑功能或等效的不平滑功能?

以下是我用来测试样条曲线与多项式拟合的脚本

 将numpy导入为np从mpl_toolkits.mplot3d导入Axes3D导入scipy.optimize导入scipy.interpolate导入matplotlib.pyplot作为plt将numpy.polynomial.polynomial导入为poly#网格和测试功能N = 9；x，y = np.linspace(-1,1，N)，np.linspace(-1,1，N)X，Y = np.网格(x，y)F =λX，Y:X + Y-1 * X * Y-(X * Y)** 2 -2 * X * Y ** 2 + X ** 2 * Y + 3 * np.exp(-((X + 1)** 2+(Y + 1)** 2)* 5)Z = F(X，Y)噪音= 0.4Z * = 1+(np.random.random(Z.shape)* 2-1)* noise#noise#更细的网格和测试功能N2 = 19；x2，y2 = np.linspace(-1,1，N2)，np.linspace(-1,1，N2)X2，Y2 = np.meshgrid(x2，y2)Z2 = F(X2，Y2)#将数据放入列表Xl = X.reshape(X.size)Yl = Y.reshape(Y.size)Zl = Z.reshape(Z.size)#多项式拟合#polyval(x，y，p)= p [0,0] + p [0,1] y + p [1,0] x + p [1,1] xy + p [1,2] xy ^ 2...， 等等#我为p使用了一个平面(1D)数组，因此需要先将其重塑为2D数组#传递给polyval顺序= 3p0 = np.zeros(order ** 2)#猜测参数(目前全部为0)f_poly = lambda x，y，p:poly.polyval2d(x，y，p.reshape((order，order)))#包装我们的多项式errf = lambda p:np.mean((f_poly(Xl，Yl，p.reshape((order，order)))-Zl)** 2)#误差函数以寻找最小平方误差溶胶= scipy.optimize.minimize(errf，p0)psol = sol ['x']#样条插值#双变量(2D)，平滑(不完全适合*点*)三次(3阶-即kx = ky = 3)样条spl = scipy.interpolate.SmoothBivariateSpline(Xl，Yl，Zl，kx = 3，ky = 3)f_spline = spl.ev#常规插值f_interp = lambda x，y:scipy.interpolate.griddata((Xl，Yl)，Zl，(x，y)，method ='cubic')# 阴谋无花果= plt.figure(1，figsize =(7,8))plt.clf()#聚适合ax = fig.add_subplot(311，投影='3d')ax.scatter3D(X2，Y2，Z2，s = 3，color ='red'，label ='actual data')适合= f_poly(X2，Y2，psol)l ='order {} poly fit'.format(order)ax.plot_wireframe(X2，Y2，fit，color ='black'，label = l)ax.scatter3D(X，Y，Z，color ='blue'，label ='noisy data')plt.legend()print(平均{}错误:{}".format(l，np.sqrt(np.mean((fit-Z2)** 2))))#样条拟合ax = fig.add_subplot(312，投影='3d')ax.scatter3D(X2，Y2，Z2，s = 3，color ='red'，label ='actual data')l =平滑花键"适合= f_spline(X2，Y2)ax.plot_wireframe(X2，Y2，fit，color ='black'，label = l)ax.scatter3D(X，Y，Z，color ='blue'，label ='noisy data')plt.legend()print(平均{}错误:{}".format(l，np.sqrt(np.mean((fit-Z2)** 2))))#适合ax = fig.add_subplot(313，projection ='3d')ax.scatter3D(X2，Y2，Z2，s = 3，color ='red'，label ='actual data')l ='三阶插值'fit = f_interp(X2，Y2)ax.plot_wireframe(X2，Y2，fit，color ='black'，label = l)ax.scatter3D(X，Y，Z，color ='blue'，label ='noisy data')plt.legend()print(平均{}错误:{}".format(l，np.sqrt(np.mean((fit-Z2)** 2))))plt.show(假)请暂停(1)raw_input('按键继续')#如果使用python3则更改为input() 

解决方案

对于非结构化网格， griddata 是正确的插值工具.但是，每次执行三角剖分(Delaunay)和插值.一种解决方法是使用

样条曲线法和Clough-Tocher插值法都是基于在网格元素上构造分段多项式函数.区别在于，对于样条曲线，网格是规则的并由算法指定(请参见

解决方案

For unstructured mesh, griddata is the right interpolation tool. However, the triangulation (Delaunay) and the interpolation is performed each time. One workaround is to use either CloughTocher2DInterpolator for a C1 smooth interpolation or LinearNDInterpolator for a linear interpolation. These are the functions actually used by griddata. The difference is that it is possible to use as input a Delaunay object and it returns an interpolation function.

Here is an example based on your code:

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

import numpy as np
from scipy.interpolate import CloughTocher2DInterpolator
from scipy.spatial import Delaunay

# Example unstructured mesh:
nodes = np.array([[-1.        , -1.        ],
       [ 1.        , -1.        ],
       [ 1.        ,  1.        ],
       [-1.        ,  1.        ],
       [ 0.        ,  0.        ],
       [-1.        ,  0.        ],
       [ 0.        , -1.        ],
       [-0.5       ,  0.        ],
       [ 0.        ,  1.        ],
       [-0.75      ,  0.4       ],
       [-0.5       ,  1.        ],
       [-1.        , -0.6       ],
       [-0.25      , -0.5       ],
       [-0.5       , -1.        ],
       [-0.20833333,  0.5       ],
       [ 1.        ,  0.        ],
       [ 0.5       ,  1.        ],
       [ 0.36174242,  0.44412879],
       [ 0.5       , -0.03786566],
       [ 0.2927264 , -0.5411368 ],
       [ 0.5       , -1.        ],
       [ 1.        ,  0.5       ],
       [ 1.        , -0.5       ]])

# Theoretical function:
def F(x, y):
    return x + y -  x*y - (x*y)**2 - 2*x*y**2 + x**2*y + 3*np.exp( -((x+1)**2 + (y+1)**2)*5 )

z = F(nodes[:, 0], nodes[:, 1])

# Finer regular grid:
N2 = 19
x2, y2 = np.linspace(-1, 1, N2), np.linspace(-1, 1, N2)
X2, Y2 = np.meshgrid(x2, y2)

# Interpolation:
tri = Delaunay(nodes)
CT_interpolator = CloughTocher2DInterpolator(tri, z)
z_interpolated = CT_interpolator(X2, Y2)

# Plot
fig = plt.figure(1, figsize=(8,14))

ax = fig.add_subplot(311, projection='3d')
ax.scatter3D(nodes[:, 0], nodes[:, 1], z, s=15, color='red', label='points')

ax.plot_wireframe(X2, Y2, z_interpolated, color='black', label='interpolated')
plt.legend();

The obtained graph is:

Both the spline method and Clough-Tocher interpolation are based on constructing a piecewise polynomial function on mesh elements. The difference is that for the spline the mesh is regular and given by the algorithm (see .get_knots()). And the coefficients are fitted so that the function is close as possible to the points and smooth (fit). For the Clough-Tocher interpolation, the mesh elements are the ones given as input. The resulting function is therefore guaranteed to go through the points.

这篇关于如何使用scipy获得非平滑2D样条插值的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！