Add your own alternative version
Stats
192.9K views 9.5K downloads 90 bookmarked
Posted
14 Apr 2008

Comments and Discussions



Great Article!, Love your work!
I would like to understand if the limit of 32 points will impact me for what I'm hoping to use your example for.
I have heart rate data from people ranging from 10 to 70 yo and its pretty noisy. I'm hoping to use Bezier Curves to smoothen the curve for further calculations. However, I think I may have 60 points, which is well above the 32. Is the limit only on 32 bit machines? Or am I misunderstanding the limit?
I think I do have a workaround bby doing several parses, starting off with 5 year increments, which means I only need 60Y / 2Y increments = 30 points. Then, to get more detail, I iterate over the increments using 5 coordinates before and 5 after to calculate another 30 points within that "Zoomed" range.
What do you reckon?
Thanks!





It may be of no use to you. I've found it to be very helpful for both sound and graphics applications.
Here's a paper that discusses different methods of interpolation when resampling a signal.
Polynomial Interpolators for HighQuality Resampling of Oversampled Audio[^]
Do be aware however, that when using splines, if the number of output points is fewer than the number of input points, things get pretty messy pretty quickly  messy as in, results are not worth the electrons needed to store the bits.
The bezier method computes points in the interval [0..1], where t=0 is the first point and t=1 is the last point. Naively, we may think that we can simply use an arbitrary number N as the number of controlpoints. Unfortunately, the equation relies on the computation of factorial  N! in fact. As we can see if we calculate that, it's huge  32! is 263130836933693530167218012160000000.0  not practical at all for any but the smallest sets of controlpoints.
The next most obvious approach is to use the bezeir function and treat the curve in a piecewise manner. With this method, we break the set of points up into groups of 4 and then use them to control our curves. This is great, but ignores one minor detail  the curve produced with this method is not continuous at all  we have abrupt changes to the slope of the line(the derivative) as each small spline finishes and the next begins. For audio applications  this is utterly horrible. Utterly.
The other approach is to again use a piecewise one. If instead of the bezier function, we choose a function that requires 4 controlpoints and only computes the interval [0..1] between points 1 and 2, we've got points 0 and 3 to control the direction that the curve leaves/enters points 1 and 2.  Basically, our overlapping use of controlpoints means that we can now ensure that the derivative of the curve remains unchanged as it transitions from one piece of the overall curve to the next. We can now approximate curves with an arbitrary number of controlpoints and can do so at an arbitrary number of points along it's length.
I really can't stress what a gem and how wellwritten the paper is!
I can toss you the interpolation code I arrived at if it's any use. I commonly throw it several seconds of 44.1khz audio  100s of 1000s of control points. A single Middle C note played on a guitar can be resampled to provide the fullrange of available notes on the instrument. It's not as good sounding as 100 high quality samples, but a single 6kb sample can be coerced into performing verywell considering. This is how ModTrackers work  a type of music that works quite similarly to .MIDI music in fact.





hi !
Given the same number of bezier control points, is there a way to get a list of curve points between two consecutive control points ?
so :
public void Bezier2D(double[] b, int cpts, double[] p)
would become :
public void Bezier2D(double[] b, int controlPointIndex, int cpts,double[] p)
where p[] would be filled with cpts points between the b[controlPointIndex] and b[controlPointIndex+1]
Thanks !
Cyril





This shouldn't be hard. It is only a matter of indexing between the points. Since your input is already ordered, you could simply take the points in between.
Another workaround is the following: Generate the full spline and sample some points in between the two. Then use the desired control points as begin and end. Then, generate another curve.





Hello
I'm searching for plotting 3D Bezier surfaces. Would you have any guess on how to do so?
(I have my own ideas, but maybe you have an exemple or working code...)
Thanks





I don't know if anyone is interested, but I translated Tolga's excellent code over to simple C. I needed it to dump out some smooth inbetween values for a throttle control system. I used it to save a microcontroller from doing expensive math zillions of times a second by precomputing a lookup table of 1024 points from just 32. Nice smooth curve. Here's the code:
#include <stdio.h>
#include <math.h>
#define ARRAY_SIZE( array ) sizeof( array ) / sizeof( array[0] )
double factorialLookup[33] = {
1.0,
1.0,
2.0,
6.0,
24.0,
120.0,
720.0,
5040.0,
40320.0,
362880.0,
3628800.0,
39916800.0,
479001600.0,
6227020800.0,
87178291200.0,
1307674368000.0,
20922789888000.0,
355687428096000.0,
6402373705728000.0,
121645100408832000.0,
2432902008176640000.0,
51090942171709440000.0,
1124000727777607680000.0,
25852016738884976640000.0,
620448401733239439360000.0,
15511210043330985984000000.0,
403291461126605635584000000.0,
10888869450418352160768000000.0,
304888344611713860501504000000.0,
8841761993739701954543616000000.0,
265252859812191058636308480000000.0,
8222838654177922817725562880000000.0,
263130836933693530167218012160000000.0
};
double factorial( int n ) {
if ( n < 0 ) { printf( "ERROR: n is less than 0\n" ); }
if ( n > 32 ) { printf( "ERROR: n is greater than 32\n" ); }
return factorialLookup[n];
}
double Ni( int n, int i ) {
double ni;
double a1 = factorial( n );
double a2 = factorial( i );
double a3 = factorial( n  i );
ni = a1 / ( a2 * a3 );
return ni;
}
double Bernstein( int n, int i, double t ) {
double basis;
double ti;
double tni;
if ( t == 0.0 && i == 0 ) {
ti = 1.0;
} else {
ti = pow( t, i );
}
if ( n == i && t == 1.0 ) {
tni = 1.0;
} else {
tni = pow( ( 1  t ), ( n  i ) );
}
basis = Ni( n, i ) * ti * tni;
return basis;
}
void Bezier2D( double b[], int bCount, int cpts, double p[] ) {
int npts = bCount / 2;
int icount, jcount;
double step, t;
icount = 0;
t = 0;
step = (double)1.0 / ( cpts  1 );
for ( int i1 = 0; i1 != cpts; i1++ ) {
if ((1.0  t) < 5e6) {
t = 1.0;
}
jcount = 0;
p[icount] = 0.0;
p[icount + 1] = 0.0;
for ( int i = 0; i != npts; i++ ) {
double basis = Bernstein(npts  1, i, t);
p[icount] += basis * b[jcount];
p[icount + 1] += basis * b[jcount + 1];
jcount = jcount +2;
}
icount += 2;
t += step;
}
}
int main( int argc, const char * argv[] ) {
double outputs[2048];
double inputs[66] = {
0, 0.0,
32, 9.0,
64, 16.0,
96, 22.0,
128, 28.0,
160, 0.0,
192, 41.0,
224, 46.0,
256, 50.0,
288, 56.0,
320, 59.0,
352, 60.0,
384, 61.0,
416, 62.0,
448, 63.0,
480, 63.0,
512, 64.0,
544, 64.0,
576, 65.0,
608, 65.0,
640, 66.0,
672, 67.0,
704, 69.0,
736, 71.0,
768, 75.0,
800, 82.0,
832, 87.0,
864, 92.0,
896, 98.0,
928, 105.0,
960, 111.0,
992, 119.0,
1023, 128.0
};
Bezier2D( inputs, ARRAY_SIZE( inputs ), 1024, outputs );
for ( int x = 0; x < 2048; x += 2 ) {
printf( "%li\t%li\n", (long)outputs[x], (long)outputs[x+1] );
}
return 0;
}

Andy
http://StuffAndyMakes.com
@StuffAndyMakes






Of course, after a little more time thinking through what I was doing AND with the realization that the 8bit AVR microcontroller I'm using only has 8K of SRAM to work with, I've tuned the curve down to just 4 points and I get amazing and perfect data:
float inputs[8] = {
1, minPower,
(1024/4), maxPower,
(1024/4)*3, minPower,
1023, maxPower
};
My curve has X values from 0 to 1023 (the analog joystick axes values coming into the microcontroller) and Y values from 0 to 127 (the possible command byte values to the motors). There are two motors, each with its own set of bytes (motor A gets 1 to 127 and motor B gets 128 to 255). Again, for optimization, only one curve is generated (the one for motor A) and the other motor's command is the motor A value with 127 added. So, as the stick is moved, the firmware maps the stick's analog value (0 to 1023) to data in the curve (0 to 127) and depending on the axis, the motor command value is passed to the motor untouched (motor A) or has 127 added to and then passed to the motor (motor B).
Here is the curve deisgn:
Office Chairiot Mark II Throttle Response Curve (Advanced)
The smooth area through the middle of the curve is to allow wiggle and sloppy stick movement while the vehicle is "parked" or idle. All of the three starter curves have that same smooth center:
Office Chairiot Mark II Throttle Response Curve (Regular)
Office Chairiot Mark II Throttle Response Curve (Kiddie/Safe)
Your code has been a TREMENDOUS help in making my code and memory usage more efficient. Thank you so much!





I intend to use your BezierCurve class to adjust the volume levels coming out of my electronic piano keyboard. I recently bought a nice keyboard but some notes are louder or softer than the rest. I'm going to write an application that applies custom velocity curves to the offending notes.
Thanks
Todd Piltingsrud





Todd,
This is an exciting application, let me know if you got any results.
Bests,





many thanks for the detailed description and the useful advice





why only 32 points? i don't understand. I have a list of point and I have an error in factorial can You help me?





Because:
I have a very simple factorial computation routine and 32 is already as much as I can compute as an int value. Since this code is only a demonstration, I have not bothered to go beyond that. Normally, there are smarter ways.
Cheers,





First of all this is great work!
The limitation of 32 can be easily removed. The issue is in the calculation of Bernstein basis, you used the formula with all the big numbers. If you canculate the basis succsseively, say
B(i+1) = B(i)*t*(ni)/[(1t)*(i+1)]
Then you can bypass all the factorial parts (remove the table, Ni).
Chengyou





Thank you Chengyou,
This is a design choice to speed up the calculations. As you mentioned there are many ways to get rid of such limitation. I think the best way would be to compute those coefficients by a table lookup to 32 and then compute the rest with a recursive formula, like the on you put.






Good work! I just wrote an article on codeproject and used your work. Of course I referenced to you (in code and in the article). You can find it here: Bezier curve for the Microsoft Web Chart control[^]
____________________________________________________________________
Never underestimate the power of stupid people in large groups.





Thanks Rainer, it's great to see such good use of it. Nice idea.
Tolga





Hi,
I found an interesting easy solution to solve Beziercurves: http://www.cubic.org/docs/bezier.htm[^]
Regards, Rainer
____________________________________________________________________
Never underestimate the power of stupid people in large groups.





Makalenizi çok beğendim.
Saygılarımla,
Erol Esen
Eroldynamic





Tesekkur ederim.
Saygılarımla,
Tolga





i'm trying to use the algorithm to generate curvy paths that i will move characters on these paths, but the problem is that this algorithm don't give a smooth curve, how can i fix that ?





Well, to my opinion, it clearly generates smooth curves... It doesn't even pass through the control points, like catmull rom splines do.
How smooth do you want?
Also, one quick idea that I can think of is to reduce the amount of control points. The less control points you have, the smoother the curve will look like.





This is driving me nuts  I can't figure out how to draw filled antialiased Béziers in C#/GDI+.
GDI+ offers nonfilled Béziers, and also offers other kinds of filled curves, but not filled Béziers. And I think I can forget about antialiasing.
So I'll probably have to roll my own Bézier drawing component. Google isn't helping me here. Anyone have any suggestions to offer?
Nice article btw
c





I am faced with the same problem.
The solution ought to be to create the boundary with nonfilled Beziers, choose any point inside the boundary and use FloodFill. I find this sometimes fails because the Bezier functions in my (admittedly old) Microsoft GDI library sometimes leak. That is, they sometimes create lines with a gap through which FloodFill leaks to the exterior.
To solve this I have written my own Bezier function which I am sure doesn't leak. I have experimented with different versions. My fastest version uses fixed point arithmetic and a buffer to hold all the points and then Polyline to draw them. It is only slightly (less than 20%) slower than the library version. Using LineTo for each point to avoid the buffer is many times slower. Using floating point instead of fixed point doubles the time but avoids a possible arithmetic overflow problem with too many segments. This limits my simplest version with a 1024 bit screen to 64 segments and a slightly more complex version to 1024 segments. Many suggest that 32 segments is enough anyway. Incidentally I set a maximum number of segments but unlike other code I have seen I use a recursive binary chop to save time by avoiding calculating and drawing zero length segments in small curves and avoiding calculating and drawing successive segments in a straight line in large curves.





Hello Burisch,
We are now some 6 month further after your cry for help.
If you still need some hints:
I recently needed a Bézier that should exactly fit into a gap in an existing multipoint graph.
Birdal's code, although very good and helpfull, didn't do that job and always left a hole at the end of the Bézier where it should, but didn't, connect to the other end of the gap in the existing graph.
The solution was found in a small change in the code:
change:
// step = (double)1.0 / (cpts  1);
into:
step = (double)1.0 / cpts;
Just the problem of trees and inbetween's.
Till sofar the simple approach.
A far more sophisticated solution to your antialiasing problem you might find in:
"Interpolation with Bezier Curves. A very simple method of smoothing polygons"
which is to be found in http://www.antigrain.com/agg2.5.zip
Good luck,
jvd





Very useful; well explained; nicely simplified.






Thanks for a great article!
Sorry but could you possibly supply the "Bernstein polynomials" function you are using as well? I've looked around for a way to do this but my math is quite rusty and I can't get it to work without it. If you don't have it, maybe a simple description of what I should code to get the same effect? (without the formula as that is all greek to me)
Thanks!





He can't because he's not really using them. He's faked the article and ripped off some code from gamedev.





Please look into the comments below for my answers to fresi.
Isn't Bernstein basis is in the code? It is in a separate function. Please download and check.





Big thanks for your article. Good roundup and not too long, considering the depth of the subject. I'd recommend making the FactorialLookupTable static so that every instance of BezierCurve shares the same table.
Regards,
André





Thank you André,
You are right about static look up table. I should update it. I will post it with the next update.
Tolga





I modified the code for my needs, to make it more easy to use when handling with Point structs:
public void Bezier2D(Point[] bezierPoints, int numOfPointsToCalculate, out List<Point> calculatedPoints) {
calculatedPoints = new List<Point>();
double t = 0;
double step = 1.0 / (numOfPointsToCalculate  1);
int numOfKnots = bezierPoints.Length  1;
for (int i1 = 0; i1 != numOfPointsToCalculate; i1++) {
if ((1.0  t) < 5e6) {
t = 1.0;
}
Point newPoint = new Point(0.0, 0.0);
for (int i = 0; i < bezierPoints.Length; i++) {
double basis = CalculateBernstein(numOfKnots, i, t);
newPoint.X += basis * bezierPoints[i].X;
newPoint.Y += basis * bezierPoints[i].Y;
}
calculatedPoints.Add(newPoint);
t += step;
}
}
Regards,
André





Thanks Andre,
I get happy when people work on such codes, modify it and post it back.
One key note is: I try not to allocate pointers inside functions. You either allocate the whole structure inside or allocate it outside and pass it as an argument. It is good in terms of convention (Saying as a c++ programmer).
Tolga





Hi Tolga,
thanks for your reply. I think you got a good point there and in fact I wasn't really happy with the initialization inside the Beziermethod. Thanks.
Regards,
André





thanks to the author!
here is my version, adopted for my own needs:
using System;
using System.Collections.Generic;
using System.Drawing;
namespace Curves
{
public static class BezierCurve
{
private static double[] FactorialLookup = new double[]
{
1.0,
1.0,
2.0,
6.0,
24.0,
120.0,
720.0,
5040.0,
40320.0,
362880.0,
3628800.0,
39916800.0,
479001600.0,
6227020800.0,
87178291200.0,
1307674368000.0,
20922789888000.0,
355687428096000.0,
6402373705728000.0,
121645100408832000.0,
2432902008176640000.0,
51090942171709440000.0,
1124000727777607680000.0,
25852016738884976640000.0,
620448401733239439360000.0,
15511210043330985984000000.0,
403291461126605635584000000.0,
10888869450418352160768000000.0,
304888344611713860501504000000.0,
8841761993739701954543616000000.0,
265252859812191058636308480000000.0,
8222838654177922817725562880000000.0,
263130836933693530167218012160000000.0
};
// just check if n is appropriate, then return the result
//private static double factorial(int n)
//{
// if (n < 0) { throw new Exception("n is less than 0"); }
// if (n > 32) { throw new Exception("n is greater than 32"); }
// return FactorialLookup[n]; /* returns the value n! as a SUMORealing point number */
//}
private static double Ni(int n, int i)
{
double ni;
double a1 = FactorialLookup[n];
double a2 = FactorialLookup[i];
double a3 = FactorialLookup[n  i];
//double a1 = factorial(n);
//double a2 = factorial(i);
//double a3 = factorial(n  i);
ni = a1 / (a2 * a3);
return ni;
}
// Calculate Bernstein basis
private static double Bernstein(int n, int i, double t)
{
double basis;
double ti; /* t^i */
double tni; /* (1  t)^i */
/* Prevent problems with pow */
if (t == 0.0 && i == 0)
ti = 1.0;
else
ti = Math.Pow(t, i);
if (n == i && t == 1.0)
tni = 1.0;
else
tni = Math.Pow((1  t), (n  i));
//Bernstein basis
basis = Ni(n, i) * ti * tni;
return basis;
}
public static List<PointF> Bezier2D(List<PointF> bezierPoints, int numOfPointsToCalculate)
{
List<PointF> calculatedPoints = new List<PointF>();
double t = 0;
double step = 1.0 / (numOfPointsToCalculate  1);
int numOfKnots = bezierPoints.Count  1;
for (int i1 = 0; i1 != numOfPointsToCalculate; i1++)
{
if ((1.0  t) < 5e6)
{
t = 1.0;
}
double x = 0.0;
double y = 0.0;
int i = 0;
foreach (PointF bezPoint in bezierPoints)
{
double basis = Bernstein(numOfKnots, i, t);
x += basis * (double)bezPoint.X;
y += basis * (double)bezPoint.Y;
i++;
}
calculatedPoints.Add(new PointF((float)x, (float)y));
t += step;
}
return calculatedPoints;
}
}
}





You clearly have a large gap in your understanding of splines. Proof of point, you describe the use of Bernstein polynomials, yet your digram is clearly explaining the De casteljau algorithm.





From Wikipedia:
In the mathematical subfield of numerical analysis the de Casteljau's algorithm, named after its inventor Paul de Casteljau, is a recursive method to evaluate polynomials in Bernstein form or Bézier curves. The de Casteljau's algorithm can also be used to split a single Bézier curve into two Bézier curves at an arbitrary parameter value.
So obviously the algorithm is used to construct Bezier Curves. Yet, if you read carefully, in my article I have two quotes:
"What this equation tells us is nothing but the formulation of the above algorithm (the midpoint iterations). It is very important in the sense that a whole algorithm could be summarized into a formula and a straightforward implementation would yield correct results."
"The objective here is to find points in the middle of 2 nearby points and iterate this until we have no more iteration. The new values of points will give us the curve. The famous Bezier equation is the exact formulation of this idea."
The algorithm may be called De Casteljau's algorithm but it is also a way of describing Bezier Curves.
Tolga





ok then seeing as though you think you understand it, using your method in your article code and all interpolate t between 01 determine the points then draw lines between said point as you would if you were to draw a spline.
If you really think your description and the image are the same thing then you should see a totally uniform distance between each line segment if your delta (or step as you call it) between t_i and t_i+1 is constant, lets say it is 0.001 and the control points will be (100,0)(50,100)(500,100)(100,0), come on show me using your code that the distance between consecutive points is equal, come on mr wikiman show me I am wrong and that you and wikipedia are correct...





Hi,
As I understand you are certainly not convinced by my statement: Berstein Polynomials is the mathematical formulation of Casteljau algorithm. I won't put the detailed proof here, but you can find any proof you want on De Casteljau algorithm here:
http://timotheegroleau.com/Flash/articles/cubic_bezier/casteljau_proof.pdf[^]
Additionally, here you will see some geometric interpretations:
http://cec.wustl.edu/~cse452/lectures/lect18_Approximation.pdf[^]
If you are not satisfied with these, please read the book:
Geometric Modeling with Splines  An Introduction
Gershon Elber, Elaine Cohen, and Richard F. Riesenfeld
For even more info you can also check this:
http://www.springerlink.com/content/b765t53217738103/fulltext.pdf[^]
As you see, it is not only me and Wikipedia telling this. These resources would give you a better proof than I do.
PS: Even though De Casteljau's algorithm is exactly stating what Bernstein Polynomials do, in many cases it is not implemented directly. The reason for this is that because of its recursive nature it is slower than implementing the formula directly. Additionally, in the Bernstein formulation you can do many precomputations which can improve performance.
Tolga





Instead of babbling on why didn't you do as I asked? was it really that hard? so you couldn't do it or you did it and found you were wrong all along. now in my my opinion the best thing you can do is to either way cleanup the article or delete it.





First of all, I won't respond to any further question looking comments of you, because they are nonsense. This is the last one I am answering.
Fresric Bouemont wrote: Instead of babbling on why didn't you do as I asked? was it really that hard? so you couldn't do it or you did it and found you were wrong all along
1) Why didn't you read the proofs? Or you read and didn't understand? I think they are straightforward and totally mathematical. They answer your question exactly. If you want to prove it is wrong, do it and show me.
2) I don't have time to write code for every dummy question that is asked around.
3) How can you comment so unconsciously on a proven mathematical fact? Don't you know any math?
4) You don't have any articles written. Take your time and spend it to write a few articles if you would like to prove me wrong that much.
Fresric Bouemont wrote: now in my my opinion the best thing you can do is to either way cleanup the article or delete it.
5) There are people who will use this article and code for learning and other purposes. Why would I remove it? Just because a timewaster wanted me to do it?
Tolga





So what you're saying is that you can't reproduce the fancy graphics you have plagiarized in your article.





One reason, I am also not responding is that I have read your other comments. Please don't do what you are trying to do and respect people's work. You are using codes from articles written here...
Tolga





Wow what an a**hole fresric bouement is! sleeze bag get off the boards please! Nice artice Tolga







General News Suggestion Question Bug Answer Joke Praise Rant Admin Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.

