The really stupid way to evaluate bezier curves is with recursion, which is \$O(2^n)\$, which can be lowered to \$O(n^2)\$ with memoization. De Casteljau's algorithm improves on this, and is still \$O(n^2)\$ but faster.
I misinterpreted De Casteljau's and thought it was the Bernstein form:
$$\sum_{i=0}^{n-1}{\binom{n}{i}p_i (1-t)^{n-i-1} t^i}$$
Where there are \$n\$ control points including the start and end named \$p_0, p_1,p_2, \dots, p_{n-1}\$.
Calculating the combination is normally \$O(n)\$ which would make this algorithm as a whole \$O(n^2)\$ since there are \$n\$ terms to sum, however, there is a cheat:
$$\binom{n}{k}=\prod_{i=1}^{k}{\frac{n+1-i}{k}}$$
We can store the result after each iteration, resulting in the entire sequence calculated in \$O(n)\$ time, making the algorithm as a whole \$O(n)\$.
\$2\leq n\leq 5\$ is the general use case (2 and 4 probably actually), so it's hardcoded. For \$n\geq 6\$, which uses this algorithm, I've done a ridiculous number of optimizations.
The smallest double increment (one mantissa bit) is \$2.2E{-16}\$, and based on some quick benchmarks with random data it never was more inaccurate than \$1E{-14}\$, so the inaccuracy is acceptable.
For very large \$n\$ (hundreds, 1200 was where I first saw it) the combination function exceeded the max value for doubles, causing the algorithm to return NaN, whereas De Casteljau's, though much slower, returned what was more or less the correct value.
The final code:
public static double bezier2(double a,double b,double t){
//Total of 3 floating point operations
// 1 2 3
return a+t*(b-a);
}
public static double bezier(double t,double... ds){
return bezier(ds,t);
}
public static double bezier(double[] ds,double t){
int count = ds.length;
switch(count){
case 0:throw new IllegalArgumentException("Must have at least two items to interpolate between");
case 1:return ds[0];
case 2:return bezier2(ds[0],ds[1],t);
case 3:{
double a = ds[0];
double b = ds[1];
double c = ds[2];
double t1 = 1d - t;
/*
* Hardcoded for n=3
* Total of 8+1=9 floating point operations
* 1 2 3 4 5 6 7 8
*/
return (a*t1+2d*b*t)*t1+c*t*t;
}
case 4:{
double a = ds[0];
double b = ds[1];
double c = ds[2];
double d = ds[3];
double t1 = 1d - t;
/*
* Hardcoded for n=4
* Total of 13+1=9 floating point operations
* 1 2 3 4 5 6 7 8 9 10 11 12 13
*/
return (a*t1+3d*b*t)*t1*t1+(3d*c*t1+d*t)*t*t;
}
case 5:{
double a = ds[0];
double b = ds[1];
double c = ds[2];
double d = ds[3];
double e = ds[4];
double t1 = 1d - t;
double i = t1*t1;
double j = t*t;
double k = t1*t;
double l = 4d*k;
/*
* Hardcoded for n=5
* Total of 12+5=17 floating point operations
* 1 2 3 4 5 6 7 8 9 10 11 12
*/
return (a*i + b*l + 6d*c*j) * i + (d*l + e*j) * j;
}
case 6:
case 7:
case 8:
case 9:
case 10:
case 11:return decasteljauBezier(ds,t);
default:{
double t1 = 1d - t;
int n1 = count - 1;
int halfn = (n1>>1)+1;
int halfn1 = halfn+1;
double[] choose;
if(count<29){
choose = new double[halfn1];
int[] chooseInt = chooseIntRange(n1,halfn);
for(int i=0;i<halfn1;i++){
choose[i]=chooseInt[i];
}
}else if(count<60){
choose = new double[halfn1];
long[] chooseLong = chooseLongRange(n1,halfn);
for(int i=0;i<halfn1;i++){
choose[i]=chooseLong[i];
}
}else{
choose = chooseDoubleRange(n1,halfn);
}
double[] terms = new double[count];
double power = 1d;
for(int i=0;i<halfn;i++){
terms[i] = ds[i] * choose[i] * power;
power *= t;
}
for(int i=halfn;i<count;i++){
terms[i] = ds[i] * choose[n1-i] * power;
power *= t;
}
power = t1;
for(int i=1;i<count;i++){
terms[n1-i] *= power;
power *= t1;
}
double sum = 0d;
for(double v:terms){
sum += v;
}
return sum;
}
}
}
public static double decasteljauBezier(double[] ds,double t){
int n = ds.length;
double[] result = Arrays.copyOf(ds, n);
for(int i=n-1;i>0;i--){
for(int j=0;j<i;j++){
result[j]+=(result[j+1]-result[j])*t;
}
}
return result[0];
}
public static int[] chooseIntRange(int n,int k){
int n1 = n+1;
int product = 1;
int[] result = new int[k+1];
result[0] = 1;
for(int i=1;i<=k;i++){
product = product*(n1-i)/i;
result[i] = product;
}
return result;
}
public static long[] chooseLongRange(int n,int k){
int n1 = n+1;
int product = 1;
long[] result = new long[k+1];
result[0] = 1;
for(int i=1;i<=k;i++){
product = product*(n1-i)/i;
result[i] = product;
}
return result;
}
public static double[] chooseDoubleRange(int n,int k){
int n1 = n+1;
double product = 1d;
double[] result = new double[k+1];
result[0] = 1;
for(int i=1;i<=k;i++){
product = (product * (n1-i)) / i;
result[i] = product;
}
return result;
}
Notes:
This was optimized by hand by me, and likely again by the compiler. It would explain how \$2\leq n\leq 5\$ was all roughly the same speed but \$n=6\$ then suddenly was \$\frac{1}{3}\$ of the speed.
These numbers aren't randomly chosen. They're based on the max values for
ints andlongs. It might be possible to increase these a bit (and therefore improve the algorithm ever so slightly) but I like safety.Semi-in-place De Casteljau's with all the optimizations I could think of. It's pretty fast but still \$O(n^2)\$.
De Casteljau's was faster for \$6\leq n\leq 11\$, so I thought of redirecting. But as it turns out, the extra cost of that single if statement and a redirect call meant that it was only worth redirecting for \$6\leq n\leq 9\$. Even that redirect might be gone in the future with more optimizing. It's not like the time isn't comparable for such a small \$n\$.
Beating \$O(n)\$ is definitely impossible, but I'm certain there is a better algorithm out there, or if not, a better way to implement this one. If not even that, surely there are optimizations I've missed.
