2D 轨迹上点的稳健样条拟合

Question

我有一些数据可以绘制二维轨迹。我想为这些数据拟合一条插值曲线，我认为每组三个连续点之间的三次样条曲线会很好用。但是，scipy 中的插值函数让我有点困惑。这是我目前拥有的：

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

xs = [398.01543948 400.99034244 402.36995272 401.05813953 398.65277778
 395.97260274 393.08474576 390.325      387.42105263 384.17073171
 380.42028986 377.20754717 373.80769231 370.04545455 366.796875
 363.33823529 359.63636364 356.22033898 352.95555556 349.41176471
 345.87878788 341.89189189 337.91666667 334.84482759 331.60273973
 328.         296.51515152 293.91176471 290.31111111 287.16666667
 283.97222222 281.56       278.79591837 276.32631579 273.65195849
 271.53191489 270.25503356 279.75497404 276.09359445 270.42035064
 298.7761194  298.74285714]
ys = [172.76204179 176.63910967 179.49095377 180.34710744 180.82075472
 181.04255319 181.07368421 181.06382979 181.1875     181.21875
 181.15909091 181.06410256 181.         181.01923077 180.84615385
 180.8125     180.69135802 180.69565217 180.68       180.68571429
 180.69117647 180.62264151 180.87323944 180.71666667 180.63333333
 180.64788732 181.4137931  180.94871795 181.31372549 181.25
 181.02       180.53030303 180.63541667 180.73529412 180.72631579
 180.00884956 179.03915758 172.17751105 171.69548635 173.73376084
 156.93617021 156.52671756]

tck, u = interpolate.splprep([xs, ys])
x_i, y_i = interpolate.splev(np.linspace(0, 1, 100), tck)

plt.plot(x_i, y_i)
plt.scatter(xs, ys, c='r')
plt.show()

结果是这样的：

现在我知道一切都像宣传的那样工作，但我想知道是否有某种方法可以对导数或生成的样条曲线的平滑度设置限制，以消除像那个丑陋的尖峰这样的伪影。

非常感谢！

Answer 1

尝试单调插值，例如Pchip 或 Akima1DInterpolator。第一个镜头可能是制作曲线的参数化（例如，通过笛卡尔距离的线长度 --- 但你需要使用各种参数化，因为它们确实会产生不同的曲线）并使用单调插值。它们支持 multidim y 值，因此您可以对坐标与弧长进行插值。

或者，您可以使用控制平滑量的splrep 的s 参数。

Answer 2

问题是您的最后五个点（下面标记为绿色）的顺序“不正确”，样条曲线需要按该顺序扔掉它们。如果您按ys 值排序最后五个点并稍微调整您的s 值，您会得到以下结果。

idx = np.arange(ys.shape[0])
idx[-5:] = np.flip(idx[-5:][np.argsort(ys[-5:])])

tck, u = interpolate.splprep([xs[idx], ys[idx]], s=40)
x_i, y_i = interpolate.splev(np.linspace(0, 1, 100), tck)

plt.plot(x_i, y_i)
plt.scatter(xs, ys, c='r')
plt.scatter(xs[-5:], ys[-5:], c='g')

所以你是说我必须手动删除/排序异常值？

没有。问题是样条插值尊重顺序。所以它不是最接近你的点的线，而是它首先接近你的第一点，第二点到你的第二点等等。一个解决方案是像我上面那样订购你的点。或者，您可以使用一种众所周知的启发式方法来解决旅行推销员问题。因为我们从一个非常合理的顺序开始，所以局部优化应该给我们一个更好的可能最优顺序。通过局部优化，我的意思是如果这会缩短总距离，我只需交换两个点，直到不再有两个点可以这样做。在你的情况下，这正是我上面提出的顺序

def combinations(arr):
    n = arr.shape[0]
    upper = np.tri(n,n,-1,dtype='bool').T
    a,b = np.meshgrid(arr,arr)
    return b[upper].reshape(-1), a[upper].reshape(-1)

idx = np.arange(xs.shape[0])
grid = np.meshgrid(idx,idx)
distances = ((xs[grid[0]]-xs[grid[1]])**2+(ys[grid[0]]-ys[grid[1]])**2)**(1/2)

idx = np.arange(distances.shape[0])

def local_opt(idx, dist):
    current_distance = calc_dist(idx, dist)
    for i,j in itertools.combinations(idx, 2):
            idx[i], idx[j] = idx[j], idx[i]
            if calc_dist(idx, dist) < current_distance:
                return local_opt(idx,dist)
            idx[i], idx[j] = idx[j], idx[i]
    return idx

def calc_dist(idx, dist):
    return dist[idx[:-1],np.roll(idx,-1)[:-1]].sum()

idx = local_opt(idx,distances)