4
\$\begingroup\$

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?

\$\endgroup\$
4
\$\begingroup\$

Recursion

[Edit]
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
\$\endgroup\$
  • \$\begingroup\$ Woah, those loops are smart. But isnt size_t the same as uint64_t on most machines? \$\endgroup\$ – user233009 Jan 10 at 18:46
  • \$\begingroup\$ @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. \$\endgroup\$ – chux Jan 10 at 19:11
  • \$\begingroup\$ @user233009 I did not know better Bezier algorithms until researching due to this post. So we both LSNED. \$\endgroup\$ – chux Jan 10 at 19:14
4
\$\begingroup\$
  • 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.

\$\endgroup\$
  • \$\begingroup\$ "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." \$\endgroup\$ – esote Jan 9 at 5:59
  • 1
    \$\begingroup\$ @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. \$\endgroup\$ – vnp Jan 9 at 6:01
  • \$\begingroup\$ Fair enough, I should've read your answer more carefully. \$\endgroup\$ – esote Jan 9 at 6:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.