|
For a double the max is in fact 170!
So, make the array bigger and top it at 170.
I did it fo c++ with a vector:
std::vector<double> factorialLookup;
double factorial(size_t n)
{
while (n >= factorialLookup.size())
{
auto s = factorialLookup.size();
if (s == 0)
factorialLookup.push_back(1);
else if (s > 170)
return factorialLookup[170];
else
{
factorialLookup.push_back( factorialLookup[s - 1] * s );
}
}
return factorialLookup[n];
}
|
|
|
|
|
Here in C# with a List
List<double> factorialLookup = new List<double>();
double factorial_N(int n)
{
while (n >= factorialLookup.Count)
{
var s = factorialLookup.Count;
if (s == 0)
factorialLookup.Add(1);
else if (s > 170) // max n! for a double
return factorialLookup[170];
else
{
factorialLookup.Add(factorialLookup[s - 1] * s);
}
}
return factorialLookup[n];
}
|
|
|
|
|
Hi Tolga,
I've been rewriting the code in C in order to get to know about Bezier curves. It's very helpful, but I have a question for you.
In the function public void Bezier2D(double[] b, int cpts, double[] p)", there is a line:
double basis = Bernstein(npts - 1, i, t); The Bernstein() function is being passed npts - 1, not npts. Are the implications fully recognized? I.e. 32 points (the maximum allowable number of points) will become a value of 31 as passed to the function. This value passes through Berstein() onwards to the Ni() and factorial() functions, and the FactorialLookup.
Is this OK..? Please let me know.
Many thanks,
Anthony
|
|
|
|
|
I incorporated your neat little class in a VB.Net app. I'm now looking at the limitation of 32 points (as mentioned above in the Q&A's) and also a way to influence the tension of the Bezier towards the control points. Any idea on that?
Here's my VB.net code for those of you interested. I've made some changes, major one is moving from an array of interlaced x and y's to a list of points (PointD is similar to the pointF structure except it uses doubles instead of singles.
Public Class BezierApprox
Public Sub New()
End Sub
Public Function BezierApproximate(points As List(Of PointD), controlPoints As Integer) As List(Of PointD)
If points Is Nothing OrElse points.Count < 2 Then Return Nothing
If controlPoints < 2 Then controlPoints = 2
Try
Dim newpoints As New List(Of PointD)
Dim x, y, basis As Double
Dim t As Double = 0
Dim stepsz As Double = 1/(controlPoints-1)
For i As Integer = 0 To controlPoints-1
If (1-t) < 5E-06 Then t = 1
x = 0
y = 0
For p As Integer = 0 To points.Count-1
basis = calcBernstein(points.Count-1,p,t)
x += basis * points(p).X
y += basis * points(p).Y
Next
newpoints.Add(New PointD(x,y))
t += stepsz
Next
Return newpoints
Catch exc As Exception
Debug.Print(exc.ToString)
Return Nothing
End Try
End Function
Public Function GetFactorial(number As Integer) As Integer
If number < 2 Then
Return 1
Else
Try
Dim fact As Integer = 1
For i As Integer = 1 To number
fact = fact * i
Next
Return fact
Catch exc As Exception
Debug.Print(exc.ToString)
Return 0
End Try
End If
End Function
Private Function calcNi(n As Integer, i As Integer) As Double
Dim ni As Double
Dim a1 As Double = Veet.Geometry.GetFactorial(n)
Dim a2 As Double = Veet.Geometry.GetFactorial(i)
Dim a3 As Double = Veet.Geometry.GetFactorial(n-i)
ni = a1/(a2*a3)
Return ni
End Function
Private Function calcBernstein(n As Integer, i As Integer, t As Double) As Double
Dim basis As Double
Dim ti As Double
Dim tni As Double
If t = 0.0 AndAlso i = 0 Then
ti = 1.0
Else
ti = Math.Pow(t,i)
End If
If n = i AndAlso t = 1.0 Then
tni = 1.0
Else
tni = Math.Pow((1-t),(n-i))
End If
basis = calcNi(n,i) * ti * tni
Return basis
End Function
End Class
|
|
|
|
|
After having implemented you class with my above VB.NET code, I ended up in the end doing it differently. See: 2D Polyline Vertex Smoothing. Anyways, your article certainly put me on the right path and the nice thing of coding is that you never know where you will end up... 
|
|
|
|
|
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 re-sampling a signal.
Polynomial Interpolators for High-Quality 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 control-points. 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 control-points.
The next most obvious approach is to use the bezeir function and treat the curve in a piece-wise 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 piece-wise one. If instead of the bezier function, we choose a function that requires 4 control-points 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 control-points 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 control-points and can do so at an arbitrary number of points along it's length.
I really can't stress what a gem and how well-written 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 re-sampled to provide the full-range 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 very-well considering. This is how Mod-Trackers 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
|
|
|
|
|
gnuplot?
Regards
Sascha
"Plexar was here"
|
|
|
|
|
Hello it was very easy to make it work in 3D, you already know how to get the X and Y and for Z it will be exactly the same. You have to extend the datastructure so you have it like XYZXYZXYZ...and then you just calc the Z by multiply the controlpoint.Z with the basic value (a little simplified explanation).
I defined a class to hold the data (Vertex) Here it is...
Class Vertex;
{
float X, Y, Z;
public static Vertex operator +(Vertex v, Vertex u)
{
return new Vertex(v.X + u.X, v.Y + u.Y, v.Z + u.Z);
}
public static Vertex operator *(float r, Vertex v)
{
return new Vertex(v.X * r, v.Y * r, v.Z * r);
}
}
public Vertex[] Bezier3D(Vertex[] b, Vertex[] p)
{
double step, t;
// Calculate points on curve
t = 0;
step = (double)1.0 / (p.Length - 1);
for (int i = 0; i < p.Length; i++)
{
if ((1.0 - t) < 5e-6)
t = 1.0;
p[i].X = p[i].Y = p[i].Z = 0f;
for (int j = 0; j < b.Length; j++)
{
double basis = Bernstein(b.Length - 1, j, t);
p[i] += (float)basis * b[j];
}
t += step;
}
return p;
}
modified 6-Jun-20 1:55am.
|
|
|
|
|
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 in-between 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) < 5e-6) {
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 8-bit 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*(n-i)/[(1-t)*(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 look-up 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
|
|
|
|
|
General 
News 
Suggestion 
Question 
Bug 
Answer 
Joke 
Praise 
Rant 
Admin 
Submit an article or tip
