# Computing bezier curves of degree n with recursive functions

I have been looking at homemade CNC and have been wondering how curves are drawn, so I looked into it and found this cool article. I then decided to try coming up with a bezier curve algorithm in C. It seems to work ok, and while I haven't tried plotting points with this exact implementation, it seems to match up with results from a previous implementation that I did plot points with.

#include <stdint.h>

static long double bbezier(long double t, long double p0, long double p1){
return ((p1 - p0) * t) + p0;
}

long double bezier(long double t, uint64_t *points, uint64_t n){
long double p0 = points[0], p1 = points[1];
if(n == 1) return bbezier(t, p0, p1);

long double q0 = bezier(t, points, n - 1),
q1 = bezier(t, points + 1, n - 1);
return bbezier(t, q0, q1);
}


I then quickly wrote this test program in C++.

#include <iostream>

extern "C" long double bezier(long double, uint64_t *, uint64_t);

uint64_t pointsx[] = {
0, 40, 100, 200
};
uint64_t pointsy[] = {
0, 150, 60, 100
};

int main(){
for(uint64_t i = 0; i <= 10000; ++i){
long double ii = i;
long double j = ii/10000;

long double x = bezier(j, pointsx, 3);
long double y = bezier(j, pointsy, 3);

std::cout << "X: " << x << ", Y: " << y << '\n';
}
return 0;
}


I wrote an implementation running in javascript from a lightly modified w3 schools canvas tutorial here to understand how bezier curves work but it only supports 3rd degree curves. It does plot the points though, and that's what I based the above implementation on.

It doesn't make any checks to ensure t is between 0 and 1 and n != 0 but I'm not too worried. The only thing I'm worried about is segfaults in cases where n is so high that you get a stack overflow but that will be a pretty crazy curve. Anyway, how does it look?

Recursion

In bezier(), the 2 recursive calls to bezier() is inefficient as it exponential grows with O(2n) and only O(n2) operations are needed. I suspect better efficiency (linear) can be had with a pre-computed weighing of the d[] terms below.

The concern about seg faulting due to excessive recursion would be mitigated with the above improvement.

I also change function bbezier() to code.

long double bezier_alt1(long double t, const uint64_t *points, size_t n) {
assert(n);
long double omt = 1.0 - t;
long double d[n];  // Save in between calculations.

for (size_t i = 0; i < n; i++) {
d[i] = omt * points[i] + t * points[i + 1];
}
while (n > 1) {
n--;
for (size_t i = 0; i < n; i++) {
d[i] = omt * d[i] + t * d[i + 1];
}
}
return d[0];
}


[Edit2]

A linear solution O(n) is possible with O(1) additional memory.

long double bezier_alt2(long double t, const uint64_t *points, size_t n) {
assert(n);
long double omt = 1.0 - t;
long double power_t = 1.0;
long double power_omt = powl(omt,n);
long double omt_div = omt != 0.0 ? 1.0/omt : 0.0;

long double sum = 0.0;
unsigned long term_n = 1;
unsigned long term_d = 1;
for (size_t i = 0; i < n; i++) {
long double y = power_omt*power_t*points[i]*term_n/term_d;
sum += y;
power_t *= t;
power_omt *= omt_div;
term_n *= (n-i);
term_d *= (i+1);
}
power_omt = 1.0;
long double y = power_omt*power_t*points[n]*term_n/term_d;
sum += y;
return sum;
}


Additional linear simplifications possible - perhaps another day.

Minor stuff

Use const

A const in the referenced data allows for some optimizations, wider application and better conveys code's intent.

// long double bezier(long double t, uint64_t *points, uint64_t n){
long double bezier(long double t, const uint64_t *points, uint64_t n){


Excessive wide type

With uint64_t n, there is no reasonable expectation that such an iteration will finish for large n.

Fortunately, n indicates the size of the array. For array indexing and sizing, using size_t. It is the right size - not too narrow, nor too wide a type.

// long double bezier(long double t, uint64_t *points, uint64_t n){
long double bezier(long double t, uint64_t *points, size_t n){


For this application, certainly unsigned would always suffice.

static

Good use of static in static long double bbezier() to keep that function local.

Missing "bezier.h"

I'd expect a bezier() declaration in a .h file and implementation in the .c file instead of extern "C" long double bezier(long double, uint64_t *, uint64_t); in main.c

// extern "C" long double bezier(long double, uint64_t *, uint64_t);
#include "bezier.h".


n range check

Perhaps in a debug build, test n.

long double bezier(long double t, const uint64_t *points, uint64_t n){
assert( n > 0);    // Low bound
assert( n < 1000); // Maybe a sane upper limit too

• Woah, those loops are smart. But isnt size_t the same as uint64_t on most machines? – user233009 Jan 10 at 18:46
• @user233009 "isn't size_t the same as uint64_t on most machines?" --> I very much doubt size_t is 64-bit on 32-bit or smaller machines. Most processors in 2019 are small embedded ones (billions per year). IAC, the assumption is not needed, serves scant benefit and incurs issues. – chux - Reinstate Monica Jan 10 at 19:11
• @user233009 I did not know better Bezier algorithms until researching due to this post. So we both LSNED. – chux - Reinstate Monica Jan 10 at 19:14
• I don't like bbezier. It collided too much with bezier, and is not very informative. It performs a linear interpolation, so why not call it interpolate?

• The p0 and p1 just add noise. Consider

    if (n == 1) {
return interpolate(t, points[0], points[1]);
}


I would seriously consider getting rid of q0 and q1:

    return interpolate(t,
bezier(t, points, n - 1),
bezier(t, points + 1, n - 1));


Don't take it as a recommendation.

• The recursion leads to the exponential time complexity. Way before you start having memory problems you'd face a performance problem. Consider computing the Bernstein form instead. It gives you linear time, and no memory problems.

• "It gives you linear time, and no memory problems." Evaluating them in Bernstein form through de Casteljau's algorithm is inefficient: $n(n+1)/2$ additions and $n(n+1)$ multiplications to calculate a point on a curve of degree $n$. Rather, I'd prefer Wang-Ball form. For reference: "Efficient algorithms for Bézier curves." – esote Jan 9 at 5:59
• @esote Did I ever mention de Casteljau? A not-so-naive implementation is linear. Two multiplications and two divisions per term. It will have some numerical issues with large n, sure. – vnp Jan 9 at 6:01