A lot of the examples of cubic bezier fitting code that I can find online are based on that "Algorithm for Automatically Fitting Digitized Curves" Graphics Gems article by Philip J. Schneider.
The problem is that the algorithm only fits "aligned" knots, knots where the tangents have geometric continuity.

This has three input points, the blue curve was fit to them (it's from this demo by the way).

Is there an algorithm that can fit a piecewise cubic bezier curve that can have sharp ("cusp") type knots? On some cases (like those 3 input points) it would make the curve fit the polyline formed by the 3 points perfectly.
Something like this (a mockup):

You could drop the number of points down to 2 for those lines that are to be straight. https://en.wikipedia.org/wiki/Bézier_curve#Linear_Bézier_curves

For reference, below is the C++ code that I use to calculate Bezier curves of arbitrary degree. I found the C version at https://stackoverflow.com/questions/785097/how-do-i-implement-a-bézier-curve-in-c

vector_4 getBezierPoint(vector<vector_4> points, float t)
{
    int i = points.size() - 1;
	
    while (i > 0)
    {
        for (int k = 0; k < i; k++)
        {
            points[k].x += t * (points[k + 1].x - points[k].x);
            points[k].y += t * (points[k + 1].y - points[k].y);
            points[k].z += t * (points[k + 1].z - points[k].z);
            points[k].w += t * (points[k + 1].w - points[k].w);
        }
	        
        i--;
    }
	    
    return points[0];
}
	

@cowcow Thanks for the post and code.

I also found this answer which suggests making all data points as knots in the curve, and then iteratively removing the knots that don't contribute too much to the curve:

https://stackoverflow.com/a/39781094