5

I have a dataset with about 50 points (x,y) and I would like to draw a smooth curve that can pass as closer as possible on those points.

I have heard about Casteljau's algorithm for splines but after hours searching on google I was not able to find a single piece of code I can use.

As far as I understood, to use this algorithm, I have to divide my dataset in groups of 4 points, right? 1234 5678 etc.. and as far as I noticed, my only problem is to find points in the middle. I mean, if I am calculating a curve for points 1234, I already have points 1 and 4 and I need to calculate 2 and 3, right? But it is a mystery to me how to do that.

I would like to ask you guys if you know some code in C, C++ or Objective-C that computes the curves based on datasets with any amount of number.

What I need is: I send the code an array with the dataset and I receive back an array with the points to draw.

My math is crap. So, please give me practical examples. No need to send me to pages with math theory. Looking at these pages makes my brain hurt...

Answer as you would ask a 10 year old child... :D

thanks.

Duck
  • 218
  • Since it's a programming example you're after, this is a good candidate to migrate to StackOverflow. – hardmath Jun 07 '11 at 20:54
  • there they tell me to post here, here you tell me to post there... – Duck Jun 07 '11 at 21:23
  • 4
    As I understand it, Casteljau's algorithm does not tell you how to compute a Bezier curve approximating a set of points. Rather, once you know the equation for the Bezier curve, Casteljau's algorithm gives you a numerically stable way of computing the points on the curve. Also, it seems quite unlikely that someone is going to write code from scratch for you... – Chris Taylor Jun 07 '11 at 21:35
  • You need to specify the problem a bit more. When you say that you want a smooth curve passing as closely as possible to the points, what do you mean by "smooth", and what is the maximal polynomial order permitted for the curve? Overhauser's cubic will give an interpolation through the points which IIRC has a continuous second derivative. – Peter Taylor Jun 07 '11 at 22:29
  • I want it to be as much as closer to the real curve (=points) as possible, just smoother. It can be any degree, but it has to be fast and work from any number of points (x,y). – Duck Jun 07 '11 at 22:36
  • To add to my comment above - if you take your $n$ points $(x_i, y_i)$ and let them be the control points of a Bezier curve, then Casteljau's algorithm will evaluate the curve at a specific point $t$ in $[0,1]$ along the curve. I have no idea if the resulting curve will be close enough to your points for you. Having said that it's unlikely someone will write code for you, I've now gone and done exactly that. Call it boredom. – Chris Taylor Jun 07 '11 at 23:05

6 Answers6

5

Here's an implementation of the Casteljau algorithm that I just wrote (in Java, though you should be able to convert to C/C++/Ob-C with little effort - I didn't use any high-level language features).

I haven't tried it so I don't know if it's correct - the best I can say is that it compiles. Here your list of $x_i$ is given in the array x, your $y_i$ are given in y and $n$ is the number of points. It takes $O(n^2)$ time to compute the location of one point along the Bezier curve fitting your points, and if you want to plot the curve with $k$ points, it will take $O(kn^2)$ time to run.

public class Casteljau {

    private double[] x;
    private double[] y;
    private int n;

    private double[][] b;

    public Casteljau(double[] x, double[] y, int n) {
        //require x.length = y.length = n
        this.x = x;
        this.y = y;
        this.n = n;
        this.b = new double[n][n];
    }

    private void init(double[] initialValues) {
        for (int i = 0; i < n; i++) {
            b[0][i] = initialValues[i];
        }
    }

    private double evaluate(double t, double[] initialValues) {
        init(initialValues);
        for (int j = 1; j < n; j++) {
            for (int i = 0; i < n - j; i++) {
                b[j][i] = b[j-1][i] * (1-t) + b[j-1][i+1] * t;
            }
        }
        return(b[n-1][0]);
    }

    public double[] getXYvalues(double t) {
        double xVal = evaluate(t, x);
        double yVal = evaluate(t, y);
        return new double[] {xVal, yVal};
    }

}
Chris Taylor
  • 28,955
  • Don't classes count as "high-level"? Or at least, they should be "middle-level". – Mateen Ulhaq Jun 08 '11 at 01:06
  • The algorithm as written here can be streamlined a bit. In particular, a triangular recursion such as de Casteljau's can be implemented to use a one-dimensional array instead of the triangular array that is normally used to display the algorithm. This requires careful overwriting, though. – J. M. ain't a mathematician Nov 24 '11 at 08:24
  • If I understand correctly, that would reduce the memory requirements to $O(n)$ but keep the time requirements at $O(n^2)$? – Chris Taylor Nov 24 '11 at 08:53
  • Yep, the routine is still $O(n^2)$ in effort. At least the storage is reduced. – J. M. ain't a mathematician Dec 05 '11 at 15:28
  • It's java, a class is the only level.

    The time requirement is really O((n-1)(n)/2) the space requirement is still O((n)(n+1)/2)You can certainly put it in a 1 d array, and index the points by calculating their location directly. http://stackoverflow.com/a/38290902/631911 But, the reduction of the storage isn't to O(n) It's still O((n)*(n+1)/2). Which is still n^2.

     – Tatarize Jul 10 '16 at 22:13
3

@Chris Taylor's implementation in C++:

struct XY
{
    double x;
    double y;
};

class Casteljau
{
private:
    double x[];
    double y[];
    int n;

    double b[][];

    void init(double initialValues[])
    {
        for(int i = 0; i < n; i++)
        {
            b[0][i] = initialValues[i];
        }
    }

    double evaluate(double t, double initialValues[])
    {
        init(initialValues);
        for(int j = 1; j < n; j++)
        {
            for(int i = 0; i < n - j; i++)
            {
                b[j][i] = b[j-1][i] * (1-t) + b[j-1][i+1] * t;
            }
        }
        return(b[n-1][0]);
    }

public:
    Casteljau(double in_x[], double in_y[], int in_n)
    {
        x = in_x;
        y = in_y;
        n = in_n;
        b = new double[n][n];
    }

    XY getXYvalues(XY &xy, double t)
    {
        double xVal = evaluate(t, x);
        double yVal = evaluate(t, y);
        xy.X = xVal;
        xy.Y = yVal;
    }
};
Mateen Ulhaq
  • 1,211
3

To interpolate a set of points smoothly I wouldn't use Bezier curves in the way you suggest, I would use something like a Catmull-Rom spline. See this link for details and source code, albeit using DirectX. I'm sure you could find other public libraries available for constructing Catmull-Rom splines if that doesn't suit your requirements.

Tom
  • 131
2

The premise of this question is wrong. Always has been. If you have 50 points and you want to smoothly curve to all of them. Do a smooth curve. It's even built into things like SVG. Basically you set one control point and interpolate control points for all the other points as the reflection of the point vs the previous point. Then the curve goes through every given point and has a nice curvy bit, that can be set based on the position of the one original control point. Bezier curves at n(50) do not actually go that close to all 50 of those points.

It's how one implements the s and t commands for smooth cubic and smooth quad bezier curves as expressed in the SVG specifications. What's more basically everything lets you use SVG so you could just use an SVG rather than implement anything. Just path d = q0,0,?,?,1,1,t x y x y x y x y x y (for all 50 points). w3.org/TR/SVG/paths.html#PathDataQuadraticBezierCommands

In fact, here's something that just randomly makes up 50 points and draws it as an svg. jsfiddle.net/Tatarize/3n0stp46 -- Keep in mind though that most applications want to smooth curve only a little. If we reflect the control point across the anchor point completely we get a control point that is a bit far away and can cause it to smooth the entire curve. More common to programs like Illustrator and Inkscape are smooth point which vary in degree. Which is to say you can draw a line between the previous control point of the last curve and the new control point of the next curve but you don't go as far. So you reflect across the anchor point by scale it back significantly so that it's smooth but only a bit so it smooths the connection parts but not the entire length of the curve. Just a bit smooth to make sure the transitions are reasonable.

Addendum: For a practical examples and understanding of splines see Freya Holmer's fantastic video on the subject. The goal of the project is to smoothly curve to a lot of different points, there's a few very reasonable methods to do that.

Tatarize
  • 177
  • Nice, I didn't know of that reflection trick. –  Jul 10 '16 at 08:01
  • It's how one implements the s and t commands for smooth cubic and smooth quad bezier curves as expressed in the SVG specifications. What's more basically everything lets you use SVG so you could just use an SVG rather than implement anything. Just path d = q0,0,?,?,1,1,t x y x y x y x y x y (for all 50 points) ...

    https://www.w3.org/TR/SVG/paths.html#PathDataQuadraticBezierCommands

     – Tatarize Jul 10 '16 at 08:31
  • In fact, here's something that just randomly makes up 50 points and draws it as an svg. https://jsfiddle.net/Tatarize/3n0stp46/ – Tatarize Jul 10 '16 at 08:58
2

Splines and Bezier curves are two different things. The former are usually interpolating, i.e. they pass through the given control points, while the latter pass "close to" them, in a certain sense.

As said by others, a Bezier of degree 50 is certainly not the right approach, it will strongly filter out the variations and won't be that close.

For interpolation, on prefers "piecewise curves", i.e. curves decomposed in several arcs with different equations, ensuring a smooth connection at the control points.

You can refere here: http://www.codeproject.com/Articles/560163/Csharp-Cubic-Spline-Interpolation. Beware to avoid the confusion between function interpolation (you want to reconstruct a function $Y = f(X)$) and curve interpolation (parametric equations $X = f(t)$, $Y = g(t)$). You are in the second case.

2

Check this out.

Quixotic
  • 22,431
user11869
  • 370