Bézier curve fitting with SciPy

Question

I have a set of points which approximate a 2D curve. I would like to use Python with numpy and scipy to find a cubic Bézier path which approximately fits the points, where I specify the exact coordinates of two endpoints, and it returns the coordinates of the other two control points.

I initially thought scipy.interpolate.splprep() might do what I want, but it seems to force the curve to pass through each one of the data points (as I suppose you would want for interpolation). I'll assume that I was on the wrong track with that.

My question is similar to this one: How can I fit a Bézier curve to a set of data?, except that they said they didn't want to use numpy. My preference would be to find what I need already implemented somewhere in scipy or numpy. Otherwise, I plan to implement the algorithm linked from one of the answers to that question, using numpy: An algorithm for automatically fitting digitized curves (pdf.page 622).

Thank you for any suggestions!

Edit: I understand that a cubic Bézier curve is not guaranteed to pass through all the points; I want one which passes through two given endpoints, and which is as close as possible to the specified interior points.

Answer 1

Here's a way to do Bezier curves with numpy:

import numpy as np from scipy.special 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() 

Answer 2

Here is a piece of python code for fitting points:

'''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

Generally on bezier curves check out Animated bezier and bezierinfo

Answer 3

@keynesiancross asked for "comments in [Roland's] code as to what the variables are" and others completely missed the stated problem. Roland started with a Bézier curve as input (to get a perfect match), which made it harder to understand both the problem and (at least for me) the solution. The difference from interpolation is easier to see for input that leaves residuals. Here is both paraphrased code and non-Bézier input -- and an unexpected outcome.

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()

This works well for fcn = np.cos but not for log. I kind of expected that the fit would use the t-component of the control points as additional degrees of freedom, as we would do by dragging the control points:

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()

The cause of failure, I guess, is that the norm measures the distance between points on the curves instead of the distance between a point on one curve to the nearest point on the other curve.