Computing bezier curves of degree n with recursive functions

Question

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?

Answer 1

Recursion

[Edit]
In bezier(), the 2 recursive calls to bezier() is inefficient as it exponential grows with O(2ⁿ) and only O(n²) 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

Answer 2

• 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.