Algorithm for deriving control points of a bezier curve from points along that curve?

Question

I've been looking for, but obviously not finding, an algorithm that will allow me to plug in a list of x,y coordinates that are known to be along a curve so as to get the 4 control points for a cubic bezier curve spit out.

To be more precise, I'm looking for an algorithm that will give me the two control points required to shape the curve while inputting a series of discrete points including the two control points which determine the start and end of the curve.

Thanks!

Edit: Okay, due to math, an old foe, I need to ask for the bezier curve of best fit to a polynomial function.

Answer 1

So I assume that the endpoints are fixed, and then you have a number of (x,y) sample points that you want to fit with a cubic Bezier.

The number of sample points that you have will determine what approach to take. Let's look through a few cases:

2 points

2 sample points is the simplest case. That gives you a total of 4 points, if you count the end points. This is the number of CVs in a cubic Bezier. To solve this, you need a parameter (t) value for both of the sample points. Then you have a system of 2 equations and 2 points that you need to solve, where the equation is the parametric equation of a Bezier curve at the t values you've chosen.

The t values can be whatever you like, but you will get better results by using either 1/3 and 2/3, or looking at relative distances, or relative distances along a baseline, depending on your data.

1 point

This is similar to 2 points, except that you have insufficient information to uniquely determine all your degrees of freedom. What I would suggest is to fit a quadratic Bezier, and then degree elevate. I wrote up a detailed example of quadratic fitting in this question.

More than 2 points

In this case, there isn't a unique solution. I have used least-squares approximation with good results. The steps are:

• Pick t values for each sample
• Build your system of equations as a matrix
• Optionally add fairing or some other smoothing function
• Solve the matrix with a least-squares solver

There is a good description of these steps in this free cagd textbook, chapter 11. It talks about fitting b-splines, but a cubic bezier is a type of b-spline (knot vector is 0,0,0,1,1,1 and has 4 points).

Answer 2

Let's say you have a curve y = f(x)

To define a bezier curve you need 4 points, like: P1x, P1y, P2x, P2y, P3x, P3y and P4x and P4y

P1 and P4 you are the begin/end points of the curve. P2 and P3 are control points. You already know where the beginning and end of the curve is. You have to calculate P2 and P3. The x coordinate P2x and P3x are easy, because you just pick them by selecting the curve's t to be eg 1/3 and 2/3. So you have P2x and P3x Then, you end up with a system of two equations and two unknowns (the P2y and P3y). After crunching some math you end up with something like this:

(My f(x) was a cubic polynomial, which also guaranteed that I would be able to fit one cubic Bezier curve to it exactly.)

/**
    @params {Object} firstPoint = {x:...,y...}
    @params {Object} lastPoint = {x:...,y...}
    @params {Object} cubicPoly Definition of a cubic polynomial in the form  y=ax^3+bx^2+c. 
  Has a method EvaluateAt, which calculates y for a particular x

*/
var CalcBezierControlPoints = function(firstPoint, lastPoint, cubicPoly) {
    var xDiff = lastPoint.X - firstPoint.X;
    var x1 = firstPoint.X + xDiff / 3.0;
    var x2 = firstPoint.X + 2.0 * xDiff / 3.0;

    var y1 = cubicPoly.EvaluateAt(x1);
    var y2 = cubicPoly.EvaluateAt(x2);

    var f1 = 0.296296296296296296296; // (1-1/3)^3
    var f2 = 0.037037037037037037037; // (1-2/3)^3
    var f3 = 0.296296296296296296296; // (2/3)^3

    var b1 = y1 - firstPoint.Y * f1 - lastPoint.Y / 27.0;
    var b2 = y2 - firstPoint.Y * f2 - f3 * lastPoint.Y;

    var c1 = (-2 * b1 + b2) / -0.666666666666666666;
    var c2 = (b2 - 0.2222222222222 * c1) / 0.44444444444444444;

    var p2 = {};
    var p3 = {};
    p2.X = x1;
    p2.Y = c1;

    p3.X = x2;
    p3.Y = c2;

    return ([p2, p3]);
}