我有一组近似于2D曲线的点.我想使用带有numpy和scipy的Python来找到一个大致适合点的三次Bézier路径,其中我指定了两个端点的精确坐标,并返回其他两个控制点的坐标.
我最初认为scipy.interpolate.splprep()可能会做我想要的,但它似乎迫使曲线通过每个数据点(因为我想你想要插值).我会假设自己走错了路.
我的问题与此类似:如何将Bézier曲线拟合为一组数据？,除了他们说他们不想使用numpy.我的偏好是找到我需要已经在scipy或numpy中实现的东西.否则,我计划使用numpy实现从该问题的一个答案链接的算法:一种自动拟合数字化曲线的算法(pdf.page 622).
谢谢你的任何建议!
编辑:据我所知,立方Bézier曲线无法保证通过所有点; 我想要一个通过两个给定端点,并且尽可能接近指定内部点的端点.
rep*_*cus 18
这是一种用numpy做Bezier曲线的方法:
import numpy as np
from scipy.misc import comb
def bernstein_poly(i, n, t):
"""
The Bernstein polynomial of n, i as a function of t
"""
return comb(n, i) * ( t**(n-i) ) * (1 - t)**i
def bezier_curve(points, nTimes=1000):
"""
Given a set of control points, return the
bezier curve defined by the control points.
points should be a list of lists, or list of tuples
such as [ [1,1],
[2,3],
[4,5], ..[Xn, Yn] ]
nTimes is the number of time steps, defaults to 1000
See http://processingjs.nihongoresources.com/bezierinfo/
"""
nPoints = len(points)
xPoints = np.array([p[0] for p in points])
yPoints = np.array([p[1] for p in points])
t = np.linspace(0.0, 1.0, nTimes)
polynomial_array = np.array([ bernstein_poly(i, nPoints-1, t) for i in range(0, nPoints) ])
xvals = np.dot(xPoints, polynomial_array)
yvals = np.dot(yPoints, polynomial_array)
return xvals, yvals
if __name__ == "__main__":
from matplotlib import pyplot as plt
nPoints = 4
points = np.random.rand(nPoints,2)*200
xpoints = [p[0] for p in points]
ypoints = [p[1] for p in points]
xvals, yvals = bezier_curve(points, nTimes=1000)
plt.plot(xvals, yvals)
plt.plot(xpoints, ypoints, "ro")
for nr in range(len(points)):
plt.text(points[nr][0], points[nr][1], nr)
plt.show()
Akl*_*lka 16
\n基于 @reptilicus 和 @Guillaume P. 的答案，以下是完整的代码：
\n获取贝塞尔曲线参数，即来自一组 X、Y 点或坐标的控制点。所需的另一个参数是近似的度数，所得控制点将为 (度 + 1)
\n
import numpy as np\nfrom scipy.special import comb\n\ndef get_bezier_parameters(X, Y, degree=3):\n """ Least square qbezier fit using penrose pseudoinverse.\n\n Parameters:\n\n X: array of x data.\n Y: array of y data. Y[0] is the y point for X[0].\n degree: degree of the B\xc3\xa9zier curve. 2 for quadratic, 3 for cubic.\n\n Based on /sf/ask/885015561/%C3%A9zier-curve-fitting-with-scipy\n and probably on the 1998 thesis by Tim Andrew Pastva, "B\xc3\xa9zier Curve Fitting".\n """\n if degree < 1:\n raise ValueError(\'degree must be 1 or greater.\')\n\n if len(X) != len(Y):\n raise ValueError(\'X and Y must be of the same length.\')\n\n if len(X) < degree + 1:\n raise ValueError(f\'There must be at least {degree + 1} points to \'\n f\'determine the parameters of a degree {degree} curve. \'\n f\'Got only {len(X)} points.\')\n\n def bpoly(n, t, k):\n """ Bernstein polynomial when a = 0 and b = 1. """\n return t ** k * (1 - t) ** (n - k) * comb(n, k)\n #return comb(n, i) * ( t**(n-i) ) * (1 - t)**i\n\n def bmatrix(T):\n """ Bernstein matrix for B\xc3\xa9zier curves. """\n return np.matrix([[bpoly(degree, t, k) for k in range(degree + 1)] for t in T])\n\n def least_square_fit(points, M):\n M_ = np.linalg.pinv(M)\n return M_ * points\n\n T = np.linspace(0, 1, len(X))\n M = bmatrix(T)\n points = np.array(list(zip(X, Y)))\n \n final = least_square_fit(points, M).tolist()\n final[0] = [X[0], Y[0]]\n final[len(final)-1] = [X[len(X)-1], Y[len(Y)-1]]\n return final\n给定贝塞尔参数（即控制点）创建贝塞尔曲线。
\ndef bernstein_poly(i, n, t):\n """\n The Bernstein polynomial of n, i as a function of t\n """\n return comb(n, i) * ( t**(n-i) ) * (1 - t)**i\n\n\ndef bezier_curve(points, nTimes=50):\n """\n Given a set of control points, return the\n bezier curve defined by the control points.\n\n points should be a list of lists, or list of tuples\n such as [ [1,1], \n [2,3], \n [4,5], ..[Xn, Yn] ]\n nTimes is the number of time steps, defaults to 1000\n\n See http://processingjs.nihongoresources.com/bezierinfo/\n """\n\n nPoints = len(points)\n xPoints = np.array([p[0] for p in points])\n yPoints = np.array([p[1] for p in points])\n\n t = np.linspace(0.0, 1.0, nTimes)\n\n polynomial_array = np.array([ bernstein_poly(i, nPoints-1, t) for i in range(0, nPoints) ])\n\n xvals = np.dot(xPoints, polynomial_array)\n yvals = np.dot(yPoints, polynomial_array)\n\n return xvals, yvals\n使用的示例数据（可以用任何数据替换，这是 GPS 数据）。
\npoints = []\nxpoints = [19.21270, 19.21269, 19.21268, 19.21266, 19.21264, 19.21263, 19.21261, 19.21261, 19.21264, 19.21268,19.21274, 19.21282, 19.21290, 19.21299, 19.21307, 19.21316, 19.21324, 19.21333, 19.21342]\nypoints = [-100.14895, -100.14885, -100.14875, -100.14865, -100.14855, -100.14847, -100.14840, -100.14832, -100.14827, -100.14823, -100.14818, -100.14818, -100.14818, -100.14818, -100.14819, -100.14819, -100.14819, -100.14820, -100.14820]\nfor i in range(len(xpoints)):\n points.append([xpoints[i],ypoints[i]])\n绘制原始点、控制点和生成的贝塞尔曲线。
\nimport matplotlib.pyplot as plt\n# Plot the original points\nplt.plot(xpoints, ypoints, "ro",label=\'Original Points\')\n# Get the Bezier parameters based on a degree.\ndata = get_bezier_parameters(xpoints, ypoints, degree=4)\nx_val = [x[0] for x in data]\ny_val = [x[1] for x in data]\nprint(data)\n# Plot the control points\nplt.plot(x_val,y_val,\'k--o\', label=\'Control Points\')\n# Plot the resulting Bezier curve\nxvals, yvals = bezier_curve(data, nTimes=1000)\nplt.plot(xvals, yvals, \'b-\', label=\'B Curve\')\nplt.legend()\nplt.show()\nRol*_*ier 11
这是一段用于拟合点的python代码:
'''least square qbezier fit using penrose pseudoinverse
>>> V=array
>>> E, W, N, S = V((1,0)), V((-1,0)), V((0,1)), V((0,-1))
>>> cw = 100
>>> ch = 300
>>> cpb = V((0, 0))
>>> cpe = V((cw, 0))
>>> xys=[cpb,cpb+ch*N+E*cw/8,cpe+ch*N+E*cw/8, cpe]
>>>
>>> ts = V(range(11), dtype='float')/10
>>> M = bezierM (ts)
>>> points = M*xys #produces the points on the bezier curve at t in ts
>>>
>>> control_points=lsqfit(points, M)
>>> linalg.norm(control_points-xys)<10e-5
True
>>> control_points.tolist()[1]
[12.500000000000037, 300.00000000000017]
'''
from numpy import array, linalg, matrix
from scipy.misc import comb as nOk
Mtk = lambda n, t, k: t**(k)*(1-t)**(n-k)*nOk(n,k)
bezierM = lambda ts: matrix([[Mtk(3,t,k) for k in range(4)] for t in ts])
def lsqfit(points,M):
M_ = linalg.pinv(M)
return M_ * points
通常在贝塞尔曲线上查看 Animated bezier和 bezierinfo
@keynesiancross要求“在[Roland]代码中就变量是什么发表评论”，其他人则完全忽略了上述问题。Roland从贝塞尔曲线开始输入（以获得完美匹配），这使得很难理解问题和（至少对我而言）解决方案。对于留下残差的输入，更容易看到与插值的区别。这既是释义代码，又是非贝塞尔（Bézier）输入-以及意外的结果。
import matplotlib.pyplot as plt
import numpy as np
from scipy.special import comb as n_over_k
Mtk = lambda n, t, k: t**k * (1-t)**(n-k) * n_over_k(n,k)
BézierCoeff = lambda ts: [[Mtk(3,t,k) for k in range(4)] for t in ts]
fcn = np.log
tPlot = np.linspace(0. ,1. , 81)
xPlot = np.linspace(0.1,2.5, 81)
tData = tPlot[0:81:10]
xData = xPlot[0:81:10]
data = np.column_stack((xData, fcn(xData))) # shapes (9,2)
Pseudoinverse = np.linalg.pinv(BézierCoeff(tData)) # (9,4) -> (4,9)
control_points = Pseudoinverse.dot(data) # (4,9)*(9,2) -> (4,2)
Bézier = np.array(BézierCoeff(tPlot)).dot(control_points)
residuum = fcn(Bézier[:,0]) - Bézier[:,1]
fig, ax = plt.subplots()
ax.plot(xPlot, fcn(xPlot), 'r-')
ax.plot(xData, data[:,1], 'ro', label='input')
ax.plot(Bézier[:,0],
Bézier[:,1], 'k-', label='fit')
ax.plot(xPlot, 10.*residuum, 'b-', label='10*residuum')
ax.plot(control_points[:,0],
control_points[:,1], 'ko:', fillstyle='none')
ax.legend()
fig.show()
这适用于，fcn = np.cos但不适用于log。我有点期望拟合将控制点的t分量用作附加的自由度，就像我们拖动控制点所做的那样：
manual_points = np.array([[0.1,np.log(.1)],[.27,-.6],[.82,.23],[2.5,np.log(2.5)]])
Bézier = np.array(BézierCoeff(tPlot)).dot(manual_points)
residuum = fcn(Bézier[:,0]) - Bézier[:,1]
fig, ax = plt.subplots()
ax.plot(xPlot, fcn(xPlot), 'r-')
ax.plot(xData, data[:,1], 'ro', label='input')
ax.plot(Bézier[:,0],
Bézier[:,1], 'k-', label='fit')
ax.plot(xPlot, 10.*residuum, 'b-', label='10*residuum')
ax.plot(manual_points[:,0],
manual_points[:,1], 'ko:', fillstyle='none')
ax.legend()
fig.show()
我认为失败的原因是，该规范衡量的是曲线上各点之间的距离，而不是一条曲线上的点到另一条曲线上最近的点之间的距离。
| 归档时间： |
|
| 查看次数： |
26100 次 |
| 最近记录： |