stepp0
stepp0

Reputation: 31

Least Squares Fit on Cubic Bezier Curve

I would like fit a cubic bezier curve on a set of 500 random points.

Here's the code I have for the bezier curve:

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

    nPoints = len(points)
    x = np.array([p[0] for p in points])
    y = 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(x, polynomial_array)
    yvals = np.dot(y, 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()

I'm aware that Numpy and Scipy have least squares methods: numpy.linalg.lstsq and scipy.optimize.least_squares

But I'm not sure how can I use them for fitting the curve on the 500 points. Can someone offer some assistance?

Thank you

Upvotes: 2

Views: 3412

Answers (2)

James Phillips
James Phillips

Reputation: 4647

The scipy documentation itself has a most excellent tutorial on using splines here:

https://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html

with lots of code, examples and cool graphs comparing different types of splines.

Upvotes: 0

Gaurav Dhama
Gaurav Dhama

Reputation: 1336

Use the function curve_fit in scipy, https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html

import numpy as np
from scipy.optimize import curve_fit
def func(x, a, b, c):
    return a * np.exp(-b * x) + c

xdata = np.linspace(0, 4, 50)
y = func(xdata, 2.5, 1.3, 0.5)
ydata = y + 0.2 * np.random.normal(size=len(xdata))

popt, pcov = curve_fit(func, xdata, ydata)

#Constrain the optimization to the region of 0 < a < 3, 0 < b < 2 and 0 < c < 1:


popt, pcov = curve_fit(func, xdata, ydata, bounds=(0, [3., 2., 1.]))

Upvotes: 1

Related Questions