We can use standard polynomial curve fitting to find a single cubic Bezier curve through any set of four points, but you're always going to be left with a free parameter problem: two of the points are going to be the start and end of the curve, so we know their time values are 0 and 1, but the two points are entirely free in terms of time value, and so you're going to have to come up with a reasonable assumption on what they should be before you can evaluate your curve fit.
See https://pomax.github.io/bezierinfo/#curvefitting for the maths if you want to implement this yourself, or find a library that does polynomial fitting (which means any half-decent statistics package), and then for timing values you have a few options:
- Assume the points line up with t=0, t=1/3, t=2/3, and t=1. This is almost always a bad idea because points are rarely evenly spaced with respects to the cubic distribution.
- Guess the t values based on linear distance, with the full length of the polyline p0-p1-p2-p3 set to 1, and the value at each point simply equal to the distance along the polyline, so t0=0, t1=dist(p0,p1), t2=dist(p0,p1)+dist(p1,2), t3=1.
- Get fancier and pick an initial set of t values, then analyze the resulting Bezier and iteratively generate new curves, optimizing on some quality like curvature homogeneity, bounding box/hull area, aligning p1/p2 on minimum radius, etc.
Of these, obviously 2 is going to give "reasonable" results with the least amount of effort, but if you're writing a graphics application, "reasonable" depends on what your users need, not what is easy for you.