C++ open source library for curve fitting

I'm searching for the most minimalist C++ open source library that allows to obtain a curve parameters (e.g. Bezier) given a set of points.

In my application the user gets some strokes by drawing on the screen. I have a track of all the points of each stroke, and would like to smooth out each stroke after it is drawn. I realized it can be achieved by curve fitting so that to turn polygon strokes into curvy strokes.

I found some papers on how to implement curve fitting, e.g. this or this, but they are not C++. While I could implement a simple curve fitter, I thought to see if there is already a good C++-based library for it, so that I could use it right away. It needs to be open source.

It would be nice if the library has no dependencies on other larger libraries.

Four points are required to uniquely describe a cubic curve (the first article you've linked covers that case). You have more than four points so are unlikely to get a perfect fit - some kind of compromise or trade-off will be required. Welcome to the black art of numerical optimisation!

The second article you've linked introduces this topic using C# and Math.NET.

If you're happy to treat the general purpose optimiser as a black box, then I prefer the Excel Solver (and a good tutorial) for learning about what an optimiser can do. Basically we define a function, then the optimiser calls it repeatedly with different parameters attempting to get the lowest return value it can.

I recommend dlib. Here I use dlib to fit a few points to a Bezier curve. I got it to compile fairly easily in Visual Studio 2013, but it will probably need tweaking to get it to work with a different compiler. The parts of dlib used here are all in headers, there's no .lib file you need to link to.

#include <stdio.h>
#include "../dlib/optimization.h"

typedef dlib::matrix<double,0,1> column_vector;

struct Point
{
double x;
double y;
};

static Point points[] =
{
{ 254 , 102 },
{ 226 , 63  },
{ 185 , 49  },
{ 146 , 74  },
{ 142 , 119 },
{ 117 , 169 },
{ 86  , 214 },
{ 40  , 200 },
};

//
// cubicBezier
//
// p1 - start point
// c1 - first control point
// c2 - second control point
// p2 - end point
//
double cubicBezier(double p1, double c1, double c2, double p2, double t)
{
double s = (1 - t);

double v = 0;
v = v + (1 * p1 * s * s * s);
v = v + (3 * c1 * s * s * t);
v = v + (3 * c2 * s * t * t);
v = v + (1 * p2 * t * t * t);

return v;
};

Point cubicBezier(double p1x, double p1y, double c1x, double c1y,
double c2x, double c2y, double p2x, double p2y,
double t)
{
Point pt;
pt.x = cubicBezier(p1x, c1x, c2x, p2x, t);
pt.y = cubicBezier(p1y, c1y, c2y, p2y, t);
return pt;
}

// Any distance function can be used for optimisation.  This one, where we want
// to find the least-squares is most common.
double dist(double x1, double y1, double x2, double y2)
{
double x = x2 - x1;
double y = y2 - y1;

return (x * x) + (y * y);
}

// This is the function that the optimiser calls repeatedly with different
// parameters, attempting to get the lowest return value it can.
double rateCurve(const column_vector& params)
{
double p1x = points.x;
double p1y = points.y;
double c1x = params(0,0);
double c1y = params(1,0);
double c2x = params(2,0);
double c2y = params(3,0);
double p2x = points.x;
double p2y = points.y;

double distances = 0;

for (Point target : points)
{
double distance = _DMAX;

for (double t = 0; t <= 1; t += 0.02)
{
Point pt = cubicBezier( p1x, p1y, c1x, c1y, c2x, c2y, p2x, p2y, t);

double dCandidate = dist(pt.x, pt.y, target.x, target.y);

distance = std::min(distance, dCandidate);
}

distances += distance;
}

// Thats the curve-fitting done.  Now incentivise slightly smoother curves.
double p1c1 = dist(p1x, p1y, c1x, c1y);
double p2c2 = dist(p2x, p2y, c2x, c2y);

return distances + pow(p1c1, 0.6) + pow(p2c2, 0.6);
}

int main(int argc, char* argv[])
{
column_vector params(4);
params = points.x, points.y, points.x, points.y;

dlib::find_min_using_approximate_derivatives(
dlib::cg_search_strategy(),
dlib::objective_delta_stop_strategy(1).be_verbose(),
rateCurve,
params,
-1);

printf("p1x = %f;\n", points.x);
printf("p1y = %f;\n", points.y);
printf("c1x = %f;\n", params(0,0));
printf("c1y = %f;\n", params(1,0));
printf("c2x = %f;\n", params(2,0));
printf("c2y = %f;\n", params(3,0));
printf("p2x = %f;\n", points.x);
printf("p2y = %f;\n", points.y);

return 0;
}

Here's the output:

iteration: 0   objective: 9740.42
iteration: 1   objective: 4880.37
iteration: 2   objective: 2872.77
iteration: 3   objective: 2523.82
iteration: 4   objective: 2048.95
iteration: 5   objective: 1680.86
iteration: 6   objective: 1519.74
iteration: 7   objective: 1366.39
iteration: 8   objective: 1330.56
iteration: 9   objective: 1285.79
iteration: 10   objective: 1275.27
iteration: 11   objective: 1274.82
p1x = 254.000000;
p1y = 102.000000;
c1x = 127.524342;
c1y = -86.849427;
c2x = 146.034795;
c2y = 283.099363;
p2x = 40.000000;
p2y = 200.000000;

And here's a visualisation showing my points and the attempt to fit them: • Thank you for pointing out to dlib. Although I ended up implementing the curve fitting myself, I guess it can be useful for other people who are in search. Sep 5 '16 at 21:51
• @vicrucann Well done - that can't have been an easy implementation. StackExchange is an indispensable resource for Q&A, but your experience of having to find a timely solution yourself seems relatively common :-/ Sep 6 '16 at 8:25

Graphics Gems has a simple C code example of Bezier curve fitting with no other library dependencies: https://github.com/erich666/GraphicsGems/blob/master/gems/FitCurves.c

(The code is public domain; see the readme in the repo.)