Кривая Безье с SciPy
У меня есть набор точек, которые приближаются к 2D-кривой. Я хотел бы использовать Python с numpy и scipy, чтобы найти кубический путь Безье, который приблизительно соответствует точкам, где я указываю точные координаты двух конечных точек, и он возвращает координаты двух других контрольных точек.
Я изначально думал scipy.interpolate.splprep() может делать то, что я хочу, но кажется, что кривая должна пройти через каждую из точек данных (как я полагаю, вы захотите для интерполяции). Я предполагаю, что я был на неправильном пути с этим.
Мой вопрос похож на этот: как я могу подогнать кривую Безье к набору данных? за исключением того, что они сказали, что не хотят использовать NumPy. Я предпочел бы найти то, что мне нужно, уже реализованное где-нибудь в scipy или numpy. В противном случае я планирую реализовать алгоритм, связанный с одним из ответов на этот вопрос, используя numpy: алгоритм для автоматической подгонки оцифрованных кривых (pdf.page 622).
Спасибо за любые предложения!
Редактировать: Я понимаю, что кубическая кривая Безье не обязательно пройдет через все точки; Мне нужен тот, который проходит через две заданные конечные точки и как можно ближе к указанным внутренним точкам.
8 ответов
Вот фрагмент кода 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
Вообще на кривых Безье посмотрите Анимированные Безье и Безьеинфо
Вот способ сделать кривые Безье с NumPy:
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()
Основываясь на ответах @reptilicus и @Guillaume P., вот полный код:
- Получите параметры Безье, то есть контрольные точки из списка точек.
- Создайте кривую Безье из параметров Безье, то есть контрольных точек.
- Постройте исходные точки, контрольные точки и полученную кривую Безье.
Получение параметров Безье, то есть контрольных точек из набора точек X,Y или координат. Другой необходимый параметр - это степень аппроксимации, и результирующие контрольные точки будут (степень + 1).
import numpy as np
from scipy.special import comb
def get_bezier_parameters(X, Y, degree=3):
""" Least square qbezier fit using penrose pseudoinverse.
Parameters:
X: array of x data.
Y: array of y data. Y[0] is the y point for X[0].
degree: degree of the Bézier curve. 2 for quadratic, 3 for cubic.
Based on https://stackoverflow.com/questions/12643079/b%C3%A9zier-curve-fitting-with-scipy
and probably on the 1998 thesis by Tim Andrew Pastva, "Bézier Curve Fitting".
"""
if degree < 1:
raise ValueError('degree must be 1 or greater.')
if len(X) != len(Y):
raise ValueError('X and Y must be of the same length.')
if len(X) < degree + 1:
raise ValueError(f'There must be at least {degree + 1} points to '
f'determine the parameters of a degree {degree} curve. '
f'Got only {len(X)} points.')
def bpoly(n, t, k):
""" Bernstein polynomial when a = 0 and b = 1. """
return t ** k * (1 - t) ** (n - k) * comb(n, k)
#return comb(n, i) * ( t**(n-i) ) * (1 - t)**i
def bmatrix(T):
""" Bernstein matrix for Bézier curves. """
return np.matrix([[bpoly(degree, t, k) for k in range(degree + 1)] for t in T])
def least_square_fit(points, M):
M_ = np.linalg.pinv(M)
return M_ * points
T = np.linspace(0, 1, len(X))
M = bmatrix(T)
points = np.array(list(zip(X, Y)))
final = least_square_fit(points, M).tolist()
final[0] = [X[0], Y[0]]
final[len(final)-1] = [X[len(X)-1], Y[len(Y)-1]]
return final
Создайте кривую Безье с учетом параметров Безье, то есть контрольных точек.
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=50):
"""
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
Использованы образцы данных (можно заменить любыми данными, это данные GPS).
points = []
xpoints = [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]
ypoints = [-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]
for i in range(len(xpoints)):
points.append([xpoints[i],ypoints[i]])
Постройте исходные точки, контрольные точки и полученную кривую Безье.
import matplotlib.pyplot as plt
# Plot the original points
plt.plot(xpoints, ypoints, "ro",label='Original Points')
# Get the Bezier parameters based on a degree.
data = get_bezier_parameters(xpoints, ypoints, degree=4)
x_val = [x[0] for x in data]
y_val = [x[1] for x in data]
print(data)
# Plot the control points
plt.plot(x_val,y_val,'k--o', label='Control Points')
# Plot the resulting Bezier curve
xvals, yvals = bezier_curve(data, nTimes=1000)
plt.plot(xvals, yvals, 'b-', label='B Curve')
plt.legend()
plt.show()
@keynesiancross попросил "прокомментировать в коде [Роланда], что это за переменные", а другие полностью пропустили указанную проблему. Роланд начал с кривой Безье в качестве входных данных (чтобы получить идеальное соответствие), что усложнило понимание как проблемы, так и (по крайней мере для меня) решения. Отличие от интерполяции легче увидеть для входных данных, которые оставляют остатки. Здесь представлен как перефразированный код, так и не входные данные Безье - и неожиданный результат.
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()
Я полагаю, что причиной отказа является то, что норма измеряет расстояние между точками на кривых, а не расстояние между точкой на одной кривой до ближайшей точки на другой кривой.
Короткий ответ: нет, потому что кривые Безье работают не так. Более длинный ответ: взгляните на сплайны Катмулла-Рома. Их довольно легко сформировать (касательный вектор в любой точке P, за исключением начала и конца, параллелен линиям {P-1,P+1}, поэтому их тоже легко программировать) и всегда проходит через точки, которые их определяют, в отличие от кривых Безье, которые интерполируют "где-то" внутри выпуклой оболочки, установленной всеми контрольными точками.
Кривая Безье не обязательно пройдет через каждую точку, которой вы ее снабжаете; контрольные точки являются произвольными (в том смысле, что не существует специального алгоритма для их нахождения, вы просто выбираете их сами) и только вытягиваете кривую в направлении.
Если вам нужна кривая, которая будет проходить через каждую точку, которой вы ее снабжаете, вам понадобится что-то вроде естественного кубического сплайна и из-за его ограничений (вы должны снабжать их растущими координатами x, или они стремятся к бесконечности) Вы, вероятно, захотите параметрический естественный кубический сплайн.
Здесь есть хорошие уроки:
У меня была та же проблема, что и в вопросе. Я взял код, предоставленный Роландом Пунтайером, и смог заставить его работать. Вот:
def get_bezier_parameters(X, Y, degree=2):
""" Least square qbezier fit using penrose pseudoinverse.
Parameters:
X: array of x data.
Y: array of y data. Y[0] is the y point for X[0].
degree: degree of the Bézier curve. 2 for quadratic, 3 for cubic.
Based on https://stackru.com/questions/12643079/b%C3%A9zier-curve-fitting-with-scipy
and probably on the 1998 thesis by Tim Andrew Pastva, "Bézier Curve Fitting".
"""
if degree < 1:
raise ValueError('degree must be 1 or greater.')
if len(X) != len(Y):
raise ValueError('X and Y must be of the same length.')
if len(X) < degree + 1:
raise ValueError(f'There must be at least {degree + 1} points to '
f'determine the parameters of a degree {degree} curve. '
f'Got only {len(X)} points.')
def bpoly(n, t, k):
""" Bernstein polynomial when a = 0 and b = 1. """
return t ** k * (1 - t) ** (n - k) * comb(n, k)
def bmatrix(T):
""" Bernstein matrix for Bézier curves. """
return np.matrix([[bpoly(degree, t, k) for k in range(degree + 1)] for t in T])
def least_square_fit(points, M):
M_ = np.linalg.pinv(M)
return M_ * points
T = np.linspace(0, 1, len(X))
M = bmatrix(T)
points = np.array(list(zip(X, Y)))
return least_square_fit(points, M).tolist()
Чтобы зафиксировать конечные точки кривой, игнорируйте первый и последний параметры, возвращаемые функцией, и используйте свои собственные точки.
То, что сказал Майк Камерманс, верно, но я также хотел бы отметить, что, насколько мне известно, сплайны катмулла можно определить в терминах кубических Безье. Итак, если у вас есть только библиотека, которая работает с кубиками, вы все равно сможете выполнять сплайны catmull-rom: