Python中的曲线插值算法

在实际应用中需要对路径或者曲线进行重采样，重采样的过程就是"曲线拟合->重采样曲线点"的过程。

1.待解决问题

如下一系列点组成的曲线，我们需要对曲线进行拟合重采样。

代码语言：javascript复制

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

from scipy.interpolate import UnivariateSpline

from scipy.interpolate import interp1d

# Define some points:
theta = np.linspace(-3, 2, 40)
points = np.vstack( (np.cos(theta), np.sin(theta)) ).T

# add some noise:
points = points   0.05 * np.random.randn(*points.shape)

plt.plot(*points.T, 'ok')
plt.plot(*points.T, '-r')

plt.show()

2. 拟合重采样过程遇到的问题

下面的方法都不Work!!

2.1 Cubic Interpolation

代码语言：javascript复制

def cubic_interpolation(fig, axs):
    xnew = np.linspace(x_arr.min(), x_arr.max(), len(x_arr))
    f = interpolate.interp1d(xnew , y_arr, kind='cubic')

    axs.plot(xnew, f(xnew))
    axs.set_title('cubic')

2.2 Cubic Spline Interpolation

代码语言：javascript复制

def cubic_spline_interpolation(fig, axs):
    xnew = np.linspace(x_arr.min(), x_arr.max(), len(x_arr))
    tck = interpolate.splrep(x_arr, y_arr, s=0)
    ynew = interpolate.splev(xnew, tck, der=0)

    axs.plot(xnew, ynew)
    axs.set_title('cubic spline')

Always Fail (ValueError: Error on input data)

原因: 函数要求x严格从小到大有序。

2.3 Parametric Spline Interpolation

代码语言：javascript复制

def parametric_spline_interpolation(fig, axs):
    xnew = np.linspace(x_arr.min(), x_arr.max(), len(x_arr))
    tck, u = interpolate.splprep([x_arr, y_arr], s=0)
    out = interpolate.splev(xnew, tck)

    axs.plot(out[0], out[1])
    axs.set_title('parametric spline')

2.4 Univariate Spline Interpolation

代码语言：javascript复制

def univariate_spline_interpolated(fig, axs):   
    s = interpolate.InterpolatedUnivariateSpline(x_arr, y_arr)# ValueError: x must be strictly increasing
    xnew = np.linspace(x_arr.min(), x_arr.max(), len(x_arr))
    ynew = s(xnew)

    axs.plot(xnew, ynew)
    axs.set_title('univariate spline')

错误：ValueError: x must be strictly increasing

2.5 Rbf Interpolation

代码语言：javascript复制

def rbf(fig, axs):
    xnew = np.linspace(x_arr.min(), x_arr.max(), len(x_arr))
    rbf = interpolate.Rbf(x_arr, y_arr) # numpy.linalg.linalg.LinAlgError: Matrix is singular.
    fi = rbf(xnew)

    axs.plot(xnew, fi)
    axs.set_title('rbf')

2.6 Linear Interpolation

代码语言：javascript复制

def linear_interpolation(fig, axs):
    xnew = np.linspace(x_arr.min(), x_arr.max(), len(x_arr))
    f = interpolate.interp1d(xnew , y_arr)

    axs.plot(xnew, f(xnew))
    axs.set_title('linear')

3.UnivariateSpline曲线拟合采样

将x和y作为曲线offset的函数分别拟合，解决了拟合函数对自变量必须严格从小到大有序的要求。

代码语言：javascript复制

# Linear length along the line:
distance = np.cumsum(np.sqrt(np.sum(np.diff(points, axis = 0)**2, axis = 1 )) )

distance = np.insert(distance, 0, 0)/distance[-1]

# Build a list of the spline function, one for each dimension:
splines = [UnivariateSpline(distance, coords, k=3, s=0) for coords in points.T]

# Computed the spline for the asked distances:
alpha = np.linspace(np.min(distance), np.max(distance), 75)
points_fitted = np.vstack(spl(alpha) for spl in splines).T

# Graph:
plt.plot(*points_fitted.T, 'ok', label='original points')
plt.plot(*points_fitted.T, '-r', label='fitted spline k=3, s=0')
plt.axis('equal')
plt.legend();
plt.xlabel('x')
plt.ylabel('y')

plt.show()

data error interpolation spline valueerror

0 人点赞