将闭合曲线拟合到一组点
Fitting a closed curve to a set of points
我有一组点 pts 形成一个循环，它看起来像这样：
这有点类似于, but instead of putting points in between pairs of points, I would like to fit a smooth curve through the points (coordinates are given at the end of the question), so I tried something similar to scipy documentation on Interpolation：
values = pts
tck = interpolate.splrep(values[:,0], values[:,1], s=1)
xnew = np.arange(2,7,0.01)
ynew = interpolate.splev(xnew, tck, der=0)
但是我得到这个错误：
ValueError: Error on input data
有什么办法可以找到合适的吗？
点坐标：
pts = array([[ 6.55525 , 3.05472 ],
[ 6.17284 , 2.802609],
[ 5.53946 , 2.649209],
[ 4.93053 , 2.444444],
[ 4.32544 , 2.318749],
[ 3.90982 , 2.2875 ],
[ 3.51294 , 2.221875],
[ 3.09107 , 2.29375 ],
[ 2.64013 , 2.4375 ],
[ 2.275444, 2.653124],
[ 2.137945, 3.26562 ],
[ 2.15982 , 3.84375 ],
[ 2.20982 , 4.31562 ],
[ 2.334704, 4.87873 ],
[ 2.314264, 5.5047 ],
[ 2.311709, 5.9135 ],
[ 2.29638 , 6.42961 ],
[ 2.619374, 6.75021 ],
[ 3.32448 , 6.66353 ],
[ 3.31582 , 5.68866 ],
[ 3.35159 , 5.17255 ],
[ 3.48482 , 4.73125 ],
[ 3.70669 , 4.51875 ],
[ 4.23639 , 4.58968 ],
[ 4.39592 , 4.94615 ],
[ 4.33527 , 5.33862 ],
[ 3.95968 , 5.61967 ],
[ 3.56366 , 5.73976 ],
[ 3.78818 , 6.55292 ],
[ 4.27712 , 6.8283 ],
[ 4.89532 , 6.78615 ],
[ 5.35334 , 6.72433 ],
[ 5.71583 , 6.54449 ],
[ 6.13452 , 6.46019 ],
[ 6.54478 , 6.26068 ],
[ 6.7873 , 5.74615 ],
[ 6.64086 , 5.25269 ],
[ 6.45649 , 4.86206 ],
[ 6.41586 , 4.46519 ],
[ 5.44711 , 4.26519 ],
[ 5.04087 , 4.10581 ],
[ 4.70013 , 3.67405 ],
[ 4.83482 , 3.4375 ],
[ 5.34086 , 3.43394 ],
[ 5.76392 , 3.55156 ],
[ 6.37056 , 3.8778 ],
[ 6.53116 , 3.47228 ]])
要拟合通过 N 个点的平滑闭合曲线，您可以使用具有以下约束的线段：
- 每条线段都必须触及它的两个端点（每条线段2个条件）
- 对于每个点，左右线段必须具有相同的导数（每个点 2 个条件 == 每个线段 2 个条件）
为了能够对每条线段总共 4 个条件有足够的自由度，每条线段的方程应为 y = ax^3 + bx^2 + cx + d。 （所以导数是 y' = 3ax^2 + 2bx + c）
按照建议设置条件将为您提供 N * 4 个线性方程，其中 N * 4 个未知数 (a1..aN, b1..bN, c1..cN, d1..dN) 可通过矩阵求逆 (numpy ).
如果点在同一条垂直线上，则需要特殊（但简单）处理，因为导数将为 "infinite"。
您的问题是因为您试图直接使用 x 和 y。您调用的插值函数假定 x 值按排序顺序排列，并且每个 x 值都将具有唯一的 y 值。
相反，您需要创建一个参数化坐标系（例如顶点索引）并分别使用它对 x 和 y 进行插值。
首先，考虑以下几点：
import numpy as np
from scipy.interpolate import interp1d # Different interface to the same function
import matplotlib.pyplot as plt
#pts = np.array([...]) # Your points
x, y = pts.T
i = np.arange(len(pts))
# 5x the original number of points
interp_i = np.linspace(0, i.max(), 5 * i.max())
xi = interp1d(i, x, kind='cubic')(interp_i)
yi = interp1d(i, y, kind='cubic')(interp_i)
fig, ax = plt.subplots()
ax.plot(xi, yi)
ax.plot(x, y, 'ko')
plt.show()
我没有关闭多边形。如果您愿意，可以将第一个点添加到数组的末尾（例如 pts = np.vstack([pts, pts[0]])
如果这样做，您会注意到多边形闭合处存在不连续性。
这是因为我们的参数化没有考虑到 polgyon 的关闭。快速修复是用 "reflected" 点填充数组：
import numpy as np
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
#pts = np.array([...]) # Your points
pad = 3
pts = np.pad(pts, [(pad,pad), (0,0)], mode='wrap')
x, y = pts.T
i = np.arange(0, len(pts))
interp_i = np.linspace(pad, i.max() - pad + 1, 5 * (i.size - 2*pad))
xi = interp1d(i, x, kind='cubic')(interp_i)
yi = interp1d(i, y, kind='cubic')(interp_i)
fig, ax = plt.subplots()
ax.plot(xi, yi)
ax.plot(x, y, 'ko')
plt.show()
或者，您可以使用专门的曲线平滑算法，例如 PEAK 或切角算法。
使用 ROOT Framework 和 pyroot 接口我能够生成以下图像
使用以下代码（我将您的数据转换为名为 data.csv 的 CSV，因此将其读入 ROOT 会更容易，并为列标题提供 xp、yp）
from ROOT import TTree, TGraph, TCanvas, TH2F
c1 = TCanvas( 'c1', 'Drawing Example', 200, 10, 700, 500 )
t=TTree('TP','Data Points')
t.ReadFile('./data.csv')
t.SetMarkerStyle(8)
t.Draw("yp:xp","","ACP")
c1.Print('pydraw.png')
实际上，您离问题的解决方案不远了。
使用 scipy.interpolate.splprep 进行参数化 B 样条插值是最简单的方法。它还原生支持闭合曲线，如果您提供 per=1 参数，
import numpy as np
from scipy.interpolate import splprep, splev
import matplotlib.pyplot as plt
# define pts from the question
tck, u = splprep(pts.T, u=None, s=0.0, per=1)
u_new = np.linspace(u.min(), u.max(), 1000)
x_new, y_new = splev(u_new, tck, der=0)
plt.plot(pts[:,0], pts[:,1], 'ro')
plt.plot(x_new, y_new, 'b--')
plt.show()
从根本上说，这种方法与@Joe Kington 的回答中的方法没有太大区别。虽然，它可能会更健壮一些，因为默认情况下，选择 i 向量的等效项是基于点之间的距离，而不仅仅是它们的索引（请参阅 splprep documentation 了解 u参数）。
我有一组点 pts 形成一个循环，它看起来像这样：
这有点类似于scipy documentation on Interpolation：
values = pts
tck = interpolate.splrep(values[:,0], values[:,1], s=1)
xnew = np.arange(2,7,0.01)
ynew = interpolate.splev(xnew, tck, der=0)
但是我得到这个错误：
ValueError: Error on input data
有什么办法可以找到合适的吗？
点坐标：
pts = array([[ 6.55525 , 3.05472 ],
[ 6.17284 , 2.802609],
[ 5.53946 , 2.649209],
[ 4.93053 , 2.444444],
[ 4.32544 , 2.318749],
[ 3.90982 , 2.2875 ],
[ 3.51294 , 2.221875],
[ 3.09107 , 2.29375 ],
[ 2.64013 , 2.4375 ],
[ 2.275444, 2.653124],
[ 2.137945, 3.26562 ],
[ 2.15982 , 3.84375 ],
[ 2.20982 , 4.31562 ],
[ 2.334704, 4.87873 ],
[ 2.314264, 5.5047 ],
[ 2.311709, 5.9135 ],
[ 2.29638 , 6.42961 ],
[ 2.619374, 6.75021 ],
[ 3.32448 , 6.66353 ],
[ 3.31582 , 5.68866 ],
[ 3.35159 , 5.17255 ],
[ 3.48482 , 4.73125 ],
[ 3.70669 , 4.51875 ],
[ 4.23639 , 4.58968 ],
[ 4.39592 , 4.94615 ],
[ 4.33527 , 5.33862 ],
[ 3.95968 , 5.61967 ],
[ 3.56366 , 5.73976 ],
[ 3.78818 , 6.55292 ],
[ 4.27712 , 6.8283 ],
[ 4.89532 , 6.78615 ],
[ 5.35334 , 6.72433 ],
[ 5.71583 , 6.54449 ],
[ 6.13452 , 6.46019 ],
[ 6.54478 , 6.26068 ],
[ 6.7873 , 5.74615 ],
[ 6.64086 , 5.25269 ],
[ 6.45649 , 4.86206 ],
[ 6.41586 , 4.46519 ],
[ 5.44711 , 4.26519 ],
[ 5.04087 , 4.10581 ],
[ 4.70013 , 3.67405 ],
[ 4.83482 , 3.4375 ],
[ 5.34086 , 3.43394 ],
[ 5.76392 , 3.55156 ],
[ 6.37056 , 3.8778 ],
[ 6.53116 , 3.47228 ]])
要拟合通过 N 个点的平滑闭合曲线，您可以使用具有以下约束的线段：
- 每条线段都必须触及它的两个端点（每条线段2个条件）
- 对于每个点，左右线段必须具有相同的导数（每个点 2 个条件 == 每个线段 2 个条件）
为了能够对每条线段总共 4 个条件有足够的自由度，每条线段的方程应为 y = ax^3 + bx^2 + cx + d。 （所以导数是 y' = 3ax^2 + 2bx + c）
按照建议设置条件将为您提供 N * 4 个线性方程，其中 N * 4 个未知数 (a1..aN, b1..bN, c1..cN, d1..dN) 可通过矩阵求逆 (numpy ).
如果点在同一条垂直线上，则需要特殊（但简单）处理，因为导数将为 "infinite"。
您的问题是因为您试图直接使用 x 和 y。您调用的插值函数假定 x 值按排序顺序排列，并且每个 x 值都将具有唯一的 y 值。
相反，您需要创建一个参数化坐标系（例如顶点索引）并分别使用它对 x 和 y 进行插值。
首先，考虑以下几点：
import numpy as np
from scipy.interpolate import interp1d # Different interface to the same function
import matplotlib.pyplot as plt
#pts = np.array([...]) # Your points
x, y = pts.T
i = np.arange(len(pts))
# 5x the original number of points
interp_i = np.linspace(0, i.max(), 5 * i.max())
xi = interp1d(i, x, kind='cubic')(interp_i)
yi = interp1d(i, y, kind='cubic')(interp_i)
fig, ax = plt.subplots()
ax.plot(xi, yi)
ax.plot(x, y, 'ko')
plt.show()
我没有关闭多边形。如果您愿意，可以将第一个点添加到数组的末尾（例如 pts = np.vstack([pts, pts[0]])
如果这样做，您会注意到多边形闭合处存在不连续性。
这是因为我们的参数化没有考虑到 polgyon 的关闭。快速修复是用 "reflected" 点填充数组：
import numpy as np
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
#pts = np.array([...]) # Your points
pad = 3
pts = np.pad(pts, [(pad,pad), (0,0)], mode='wrap')
x, y = pts.T
i = np.arange(0, len(pts))
interp_i = np.linspace(pad, i.max() - pad + 1, 5 * (i.size - 2*pad))
xi = interp1d(i, x, kind='cubic')(interp_i)
yi = interp1d(i, y, kind='cubic')(interp_i)
fig, ax = plt.subplots()
ax.plot(xi, yi)
ax.plot(x, y, 'ko')
plt.show()
或者，您可以使用专门的曲线平滑算法，例如 PEAK 或切角算法。
使用 ROOT Framework 和 pyroot 接口我能够生成以下图像
使用以下代码（我将您的数据转换为名为 data.csv 的 CSV，因此将其读入 ROOT 会更容易，并为列标题提供 xp、yp）
from ROOT import TTree, TGraph, TCanvas, TH2F
c1 = TCanvas( 'c1', 'Drawing Example', 200, 10, 700, 500 )
t=TTree('TP','Data Points')
t.ReadFile('./data.csv')
t.SetMarkerStyle(8)
t.Draw("yp:xp","","ACP")
c1.Print('pydraw.png')
实际上，您离问题的解决方案不远了。
使用 scipy.interpolate.splprep 进行参数化 B 样条插值是最简单的方法。它还原生支持闭合曲线，如果您提供 per=1 参数，
import numpy as np
from scipy.interpolate import splprep, splev
import matplotlib.pyplot as plt
# define pts from the question
tck, u = splprep(pts.T, u=None, s=0.0, per=1)
u_new = np.linspace(u.min(), u.max(), 1000)
x_new, y_new = splev(u_new, tck, der=0)
plt.plot(pts[:,0], pts[:,1], 'ro')
plt.plot(x_new, y_new, 'b--')
plt.show()
从根本上说，这种方法与@Joe Kington 的回答中的方法没有太大区别。虽然，它可能会更健壮一些，因为默认情况下，选择 i 向量的等效项是基于点之间的距离，而不仅仅是它们的索引（请参阅 splprep documentation 了解 u参数）。