pythonnumpymatplotlibscipybezier
Deriving Cubic Bezier Curve control points & handles from series of points
I am trying to find the control points and handles of a Cubic Bezier curve from a series of points. My current code is below (credit to Zero Zero on the Python Discord). The Cubic Spline is creating the desired fit, but the handles (in orange) are incorrect. How may I find the handles of this curve?
Thank you!
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()
Solution
The answer for my case was a Bezier best fit function that accepts an input of point values, fits the points to a Cubic Spline, and outputs the Bézier handles of the curve by finding their coefficients.
Here is one such script, fitCurves, which can be used like so:
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)
