Ytdyklly

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?

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

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 date 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 alway suffice.

static

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

Missing "bezier.h"

I'd expect a bezier() declarations 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 an sane upper limit too