我正在尝试从一系列点中找到三次贝塞尔曲线的控制点和手柄。我当前的代码如下(归功于 Python Discord 上的零零)。 Cubic Spline 正在创建所需的拟合,但手柄(橙色)不正确。我怎样才能找到这条曲线的句柄?
谢谢!
import numpy as np
import scipy as sp
def fit_curve(points):
# Fit a cubic bezier curve to the points
curve = sp.interpolate.CubicSpline(points[:, 0], points[:, 1], bc_type=((1, 0.0), (1, 0.0)))
# Get 4 control points for the curve
p = np.zeros((4, 2))
p[0, :] = points[0, :]
p[3, :] = points[-1, :]
p[1, :] = points[0, :] + 0.3 * (points[-1, :] - points[0, :])
p[2, :] = points[-1, :] - 0.3 * (points[-1, :] - points[0, :])
return p, curve
ypoints = [0.0, 0.03771681353260319, 0.20421680080883106, 0.49896111463402026, 0.7183501026981503, 0.8481517096346528, 0.9256128196832564, 0.9705404287079152, 0.9933297674379904, 1.0]
xpoints = [x for x in range(len(ypoints))]
points = np.array([xpoints, ypoints]).T
from scipy.interpolate import splprep, splev
tck, u = splprep([xpoints, ypoints], s=0)
#print(tck, u)
xnew, ynew = splev(np.linspace(0, 1, 100), tck)
# Plot the original points and the Bézier curve
import matplotlib.pyplot as plt
#plt.plot(xpoints, ypoints, 'x', xnew, ynew, xpoints, ypoints, 'b')
plt.axis([0, 10, -0.05, 1.05])
plt.legend(['Points', 'Bézier curve', 'True curve'])
plt.title('Bézier curve fitting')
# Get the curve
p, curve = fit_curve(points)
# Plot the points and the curve
plt.plot(points[:, 0], points[:, 1], 'o')
plt.plot(p[:, 0], p[:, 1], 'o')
plt.plot(np.linspace(0, 9, 100), curve(np.linspace(0, 9, 100)))
plt.show()
我的案例的答案是 Bezier
best fit
函数,它接受点值的输入,将这些点拟合到三次样条曲线,并通过找到它们的系数来输出曲线的 Bézier 句柄。
这是一个这样的脚本,fitCurves,可以像这样使用:
import numpy as np
from fitCurve import fitCurve
import matplotlib.pyplot as plt
y = [0.0,
0.03771681353260319,
0.20421680080883106,
0.49896111463402026,
0.7183501026981503,
0.8481517096346528,
0.9256128196832564,
0.9705404287079152,
0.9933297674379904,
1.0]
x = np.linspace(0, 1, len(y))
pts = np.array([x,y]).T
bezier_handles = fitCurve(points=pts , maxError=20)
x_bez = []
y_bez = []
for bez in bezier_handles:
for pt in bez:
x_bez.append(pt[0])
y_bez.append(pt[1])
plt.plot(pts[:,0], pts[:,1], 'bo-', label='Points')
plt.plot(x_bez[:2], y_bez[:2], 'ro--', label='Handle') # handle 1
plt.plot(x_bez[2:4], y_bez[2:4], 'ro--') # handle 2
plt.legend()
plt.show()
fitCurve.py
from numpy import *
""" Python implementation of
Algorithm for Automatically Fitting Digitized Curves
by Philip J. Schneider
"Graphics Gems", Academic Press, 1990
"""
# evaluates cubic bezier at t, return point
def q(ctrlPoly, t):
return (1.0-t)**3 * ctrlPoly[0] + 3*(1.0-t)**2 * t * ctrlPoly[1] + 3*(1.0-t)* t**2 * ctrlPoly[2] + t**3 * ctrlPoly[3]
# evaluates cubic bezier first derivative at t, return point
def qprime(ctrlPoly, t):
return 3*(1.0-t)**2 * (ctrlPoly[1]-ctrlPoly[0]) + 6*(1.0-t) * t * (ctrlPoly[2]-ctrlPoly[1]) + 3*t**2 * (ctrlPoly[3]-ctrlPoly[2])
# evaluates cubic bezier second derivative at t, return point
def qprimeprime(ctrlPoly, t):
return 6*(1.0-t) * (ctrlPoly[2]-2*ctrlPoly[1]+ctrlPoly[0]) + 6*(t) * (ctrlPoly[3]-2*ctrlPoly[2]+ctrlPoly[1])
# Fit one (ore more) Bezier curves to a set of points
def fitCurve(points, maxError):
leftTangent = normalize(points[1] - points[0])
rightTangent = normalize(points[-2] - points[-1])
return fitCubic(points, leftTangent, rightTangent, maxError)
def fitCubic(points, leftTangent, rightTangent, error):
# Use heuristic if region only has two points in it
if (len(points) == 2):
dist = linalg.norm(points[0] - points[1]) / 3.0
bezCurve = [points[0], points[0] + leftTangent * dist, points[1] + rightTangent * dist, points[1]]
return [bezCurve]
# Parameterize points, and attempt to fit curve
u = chordLengthParameterize(points)
bezCurve = generateBezier(points, u, leftTangent, rightTangent)
# Find max deviation of points to fitted curve
maxError, splitPoint = computeMaxError(points, bezCurve, u)
if maxError < error:
return [bezCurve]
# If error not too large, try some reparameterization and iteration
if maxError < error**2:
for i in range(20):
uPrime = reparameterize(bezCurve, points, u)
bezCurve = generateBezier(points, uPrime, leftTangent, rightTangent)
maxError, splitPoint = computeMaxError(points, bezCurve, uPrime)
if maxError < error:
return [bezCurve]
u = uPrime
# Fitting failed -- split at max error point and fit recursively
beziers = []
centerTangent = normalize(points[splitPoint-1] - points[splitPoint+1])
beziers += fitCubic(points[:splitPoint+1], leftTangent, centerTangent, error)
beziers += fitCubic(points[splitPoint:], -centerTangent, rightTangent, error)
return beziers
def generateBezier(points, parameters, leftTangent, rightTangent):
bezCurve = [points[0], None, None, points[-1]]
# compute the A's
A = zeros((len(parameters), 2, 2))
for i, u in enumerate(parameters):
A[i][0] = leftTangent * 3*(1-u)**2 * u
A[i][1] = rightTangent * 3*(1-u) * u**2
# Create the C and X matrices
C = zeros((2, 2))
X = zeros(2)
for i, (point, u) in enumerate(zip(points, parameters)):
C[0][0] += dot(A[i][0], A[i][0])
C[0][1] += dot(A[i][0], A[i][1])
C[1][0] += dot(A[i][0], A[i][1])
C[1][1] += dot(A[i][1], A[i][1])
tmp = point - q([points[0], points[0], points[-1], points[-1]], u)
X[0] += dot(A[i][0], tmp)
X[1] += dot(A[i][1], tmp)
# Compute the determinants of C and X
det_C0_C1 = C[0][0] * C[1][1] - C[1][0] * C[0][1]
det_C0_X = C[0][0] * X[1] - C[1][0] * X[0]
det_X_C1 = X[0] * C[1][1] - X[1] * C[0][1]
# Finally, derive alpha values
alpha_l = 0.0 if det_C0_C1 == 0 else det_X_C1 / det_C0_C1
alpha_r = 0.0 if det_C0_C1 == 0 else det_C0_X / det_C0_C1
# If alpha negative, use the Wu/Barsky heuristic (see text) */
# (if alpha is 0, you get coincident control points that lead to
# divide by zero in any subsequent NewtonRaphsonRootFind() call. */
segLength = linalg.norm(points[0] - points[-1])
epsilon = 1.0e-6 * segLength
if alpha_l < epsilon or alpha_r < epsilon:
# fall back on standard (probably inaccurate) formula, and subdivide further if needed.
bezCurve[1] = bezCurve[0] + leftTangent * (segLength / 3.0)
bezCurve[2] = bezCurve[3] + rightTangent * (segLength / 3.0)
else:
# First and last control points of the Bezier curve are
# positioned exactly at the first and last data points
# Control points 1 and 2 are positioned an alpha distance out
# on the tangent vectors, left and right, respectively
bezCurve[1] = bezCurve[0] + leftTangent * alpha_l
bezCurve[2] = bezCurve[3] + rightTangent * alpha_r
return bezCurve
def reparameterize(bezier, points, parameters):
return [newtonRaphsonRootFind(bezier, point, u) for point, u in zip(points, parameters)]
def newtonRaphsonRootFind(bez, point, u):
"""
Newton's root finding algorithm calculates f(x)=0 by reiterating
x_n+1 = x_n - f(x_n)/f'(x_n)
We are trying to find curve parameter u for some point p that minimizes
the distance from that point to the curve. Distance point to curve is d=q(u)-p.
At minimum distance the point is perpendicular to the curve.
We are solving
f = q(u)-p * q'(u) = 0
with
f' = q'(u) * q'(u) + q(u)-p * q''(u)
gives
u_n+1 = u_n - |q(u_n)-p * q'(u_n)| / |q'(u_n)**2 + q(u_n)-p * q''(u_n)|
"""
d = q(bez, u)-point
numerator = (d * qprime(bez, u)).sum()
denominator = (qprime(bez, u)**2 + d * qprimeprime(bez, u)).sum()
if denominator == 0.0:
return u
else:
return u - numerator/denominator
def chordLengthParameterize(points):
u = [0.0]
for i in range(1, len(points)):
u.append(u[i-1] + linalg.norm(points[i] - points[i-1]))
for i, _ in enumerate(u):
u[i] = u[i] / u[-1]
return u
def computeMaxError(points, bez, parameters):
maxDist = 0.0
splitPoint = len(points)/2
for i, (point, u) in enumerate(zip(points, parameters)):
dist = linalg.norm(q(bez, u)-point)**2
if dist > maxDist:
maxDist = dist
splitPoint = i
return maxDist, splitPoint
def normalize(v):
return v / linalg.norm(v)