非单调数据的三次样条（不是一维函数）

Question

我有一条曲线如下图：

此图的 x 坐标和 y 坐标为：
path_x= (4.0, 5.638304088577984, 6.785456961280076, 5.638304088577984, 4.0)
path_y =(0.0, 1.147152872702092, 2.7854569612800755, 4.423761049858059, 3.2766081771559668)

我通过以下方式获得了上图：

x_min =min(path_x)-1
x_max =max(path_x)+1
y_min =min(path_y)-1
y_max =max(path_y)+1

num_pts = len(path_x)

fig = plt.figure(figsize=(8,8))
#fig = plt.figure()
plt.suptitle("Curve and the boundary")
ax = fig.add_subplot(1,1,1)

ax.set_xlim([min(x_min,y_min),max(x_max,y_max)])
ax.set_ylim([min(x_min,y_min),max(x_max,y_max)])
ax.plot(path_x,path_y)

现在我的意图是使用三次样条曲线绘制平滑曲线。但看起来对于三次样条曲线，您需要 x coordinates 以升序排列。而在这种情况下，x values 和 y values 都不是升序。
这也不是一个功能。也就是说，x value 映射了范围内的多个元素。

我还查看了this 的帖子。但是我想不出一个合适的方法来解决我的问题。

非常感谢您在这方面的帮助

Answer 1

按照 cmets 中的建议，您始终可以使用任意（和线性！）参数对任何曲线/曲面进行参数化。

例如，将t 定义为一个参数，这样您就可以得到x=x(t) 和y=y(t)。由于t 是任意的，因此您可以将其定义为：在t=0，您将获得第一个path_x[0],path_y[0]，在t=1，您将获得最后一对坐标path_x[-1],path_y[-1]。

这是使用scipy.interpolate的代码

import numpy
import scipy.interpolate
import matplotlib.pyplot as plt

path_x = numpy.asarray((4.0, 5.638304088577984, 6.785456961280076, 5.638304088577984, 4.0),dtype=float)
path_y = numpy.asarray((0.0, 1.147152872702092, 2.7854569612800755, 4.423761049858059, 3.2766081771559668),dtype=float)

# defining arbitrary parameter to parameterize the curve
path_t = numpy.linspace(0,1,path_x.size)

# this is the position vector with
# x coord (1st row) given by path_x, and
# y coord (2nd row) given by path_y
r = numpy.vstack((path_x.reshape((1,path_x.size)),path_y.reshape((1,path_y.size))))

# creating the spline object
spline = scipy.interpolate.interp1d(path_t,r,kind='cubic')

# defining values of the arbitrary parameter over which
# you want to interpolate x and y
# it MUST be within 0 and 1, since you defined
# the spline between path_t=0 and path_t=1
t = numpy.linspace(numpy.min(path_t),numpy.max(path_t),100)

# interpolating along t
# r[0,:] -> interpolated x coordinates
# r[1,:] -> interpolated y coordinates
r = spline(t)

plt.plot(path_x,path_y,'or')
plt.plot(r[0,:],r[1,:],'-k')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

有输出

Answer 2

对于非升序 x 样条线可以很容易地计算，如果你同时使用另一个参数 t 的 x 和 y 函数：x(t)，y(t)。

在您的情况下，您有 5 分，因此 t 应该只是这些点的枚举，即 t = 0, 1, 2, 3, 4 为 5 分。

所以如果x = [5, 2, 7, 3, 6] 那么x(t) = x(0) = 5、x(1) = 2、x(2) = 7、x(3) = 3、x(4) = 6。 y 也一样。

然后计算x(t) 和y(t) 的样条函数。然后计算所有许多中间t 点中的样条值。最后，只需将所有计算值 x(t) 和 y(t) 用作函数 y(x)。

在我使用 Numpy 从头开始​​实现三次样条计算之前，如果您不介意，我会在下面的示例中使用此代码（它可能对您了解样条数学很有用），替换为您的库函数。同样在我的代码中，您可以看到 numba 行被注释掉，如果您愿意，可以使用这些 Numba 注释来加快计算速度。

您必须查看代码底部的main()函数，它显示了如何计算和使用x(t)和y(t)。

Try it online!

import numpy as np, matplotlib.pyplot as plt

# Solves linear system given by Tridiagonal Matrix
# Helper for calculating cubic splines
#@numba.njit(cache = True, fastmath = True, inline = 'always')
def tri_diag_solve(A, B, C, F):
    n = B.size
    assert A.ndim == B.ndim == C.ndim == F.ndim == 1 and (
        A.size == B.size == C.size == F.size == n
    ) #, (A.shape, B.shape, C.shape, F.shape)
    Bs, Fs = np.zeros_like(B), np.zeros_like(F)
    Bs[0], Fs[0] = B[0], F[0]
    for i in range(1, n):
        Bs[i] = B[i] - A[i] / Bs[i - 1] * C[i - 1]
        Fs[i] = F[i] - A[i] / Bs[i - 1] * Fs[i - 1]
    x = np.zeros_like(B)
    x[-1] = Fs[-1] / Bs[-1]
    for i in range(n - 2, -1, -1):
        x[i] = (Fs[i] - C[i] * x[i + 1]) / Bs[i]
    return x
    
# Calculate cubic spline params
#@numba.njit(cache = True, fastmath = True, inline = 'always')
def calc_spline_params(x, y):
    a = y
    h = np.diff(x)
    c = np.concatenate((np.zeros((1,), dtype = y.dtype),
        np.append(tri_diag_solve(h[:-1], (h[:-1] + h[1:]) * 2, h[1:],
        ((a[2:] - a[1:-1]) / h[1:] - (a[1:-1] - a[:-2]) / h[:-1]) * 3), 0)))
    d = np.diff(c) / (3 * h)
    b = (a[1:] - a[:-1]) / h + (2 * c[1:] + c[:-1]) / 3 * h
    return a[1:], b, c[1:], d

# Spline value calculating function, given params and "x"
#@numba.njit(cache = True, fastmath = True, inline = 'always')
def func_spline(x, ix, x0, a, b, c, d):
    dx = x - x0[1:][ix]
    return a[ix] + (b[ix] + (c[ix] + d[ix] * dx) * dx) * dx

# Compute piece-wise spline function for "x" out of sorted "x0" points
#@numba.njit([f'f{ii}[:](f{ii}[:], f{ii}[:], f{ii}[:], f{ii}[:], f{ii}[:], f{ii}[:])' for ii in (4, 8)],
#    cache = True, fastmath = True, inline = 'always')
def piece_wise_spline(x, x0, a, b, c, d):
    xsh = x.shape
    x = x.ravel()
    ix = np.searchsorted(x0[1 : -1], x)
    y = func_spline(x, ix, x0, a, b, c, d)
    y = y.reshape(xsh)
    return y

def main():
    x0 = np.array([4.0, 5.638304088577984, 6.785456961280076, 5.638304088577984, 4.0])
    y0 = np.array([0.0, 1.147152872702092, 2.7854569612800755, 4.423761049858059, 3.2766081771559668])
    t0 = np.arange(len(x0)).astype(np.float64)
    plt.plot(x0, y0)
    vs = []
    for e in (x0, y0):
        a, b, c, d = calc_spline_params(t0, e)
        x = np.linspace(0, t0[-1], 100)
        vs.append(piece_wise_spline(x, t0, a, b, c, d))
    plt.plot(vs[0], vs[1])
    plt.show()

if __name__ == '__main__':
    main()

输出：