C++ open source library for curve fitting

Question

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.

Answer 1

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[0].x;
    double p1y = points[0].y;
    double c1x = params(0,0);
    double c1y = params(1,0);
    double c2x = params(2,0);
    double c2y = params(3,0);
    double p2x = points[7].x;
    double p2y = points[7].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[7].x, points[7].y, points[0].x, points[0].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[0].x);
    printf("p1y = %f;\n", points[0].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[7].x);
    printf("p2y = %f;\n", points[7].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:

Answer 2

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