Memo for myself about the implementation of Bezier curves/surfaces.

They are defined by control points, which are vectors of same dimensions as the output dimension. There are \(order+1\) control points per input dimension, meaning there are \((order+1)^{dimIn}\) control points in total. \(order\) is the 3rd attribute (with the numbers of input and output dimension) defining a Bezier curve/surface.

The evalution of a Bezier curve/surface for a given input \(\overrightarrow{x}\) is equal to the sum of the weighted control points, where the weight of a control point is equal to the product of the weight in each dimension of this control point at the given input. The weight of the \(j\)-th control point in a given dimension \(i\) is equal to \(B(j)*(1-x_i)^{order-j}*x_i^j\), where \(B(j)\) is the \(j\)-th binomial value for \(order\). Binomial values are the Pascal triangle values (\(order\equiv row\)):

1 1 1 1 2 1 1 3 3 1 1 4 6 4 1 ...

Using an array of vectors for the control point, and look up tables for the binomial values and indices of each control point in each dimension allows an efficient implementation of the evaluation of a Bezier curve/surface.

Given that the Bezier curve/surface structure is implemented as follow:

// Bezier object typedef struct CapyBezier { // Order uint8_t const order; // Input dimension uint8_t const dimIn; // Output dimension uint8_t const dimOut; // Number of control points uint64_t const nbCtrls; // Look up table for the binomial values double const* const binomialLUT; // Look up table for the index in a given dimension of a given control // points uint8_t const** const idxCtrlLUT; // Control points. Ordered as follow: // (iDimIn_0, iDimIn_1, ..., iDimIn_{dimIn-1}) // 0 <= iDim_j <= order, 0 <= iDimOut < dimOut CapyArrVec ctrls; // Vectors to memorise the control points' weight during evaluation. // Ordered as that->ctrls CapyVec* ctrlWeights; // Vectors to memorise the result during evaluation. Ordered as // follow: // (iCtrl, iDimOut) CapyVec* out; ... } CapyBezier;

The look up table for the binomial can be prepared as follow:

static double* CreateBinomialLUT(uint8_t const order) { // Static look up table static double binomial[255][255]; // If the LUT hasn't been calculated yet static bool flag = false; if (flag == false) { // Initialise the look up table loop (n, 255){ binomial[n][0] = 1.0; for (int j = 1; j < n; ++j) binomial[n][j] = binomial[n - 1][j - 1] + binomial[n - 1][j]; binomial[n][n] = 1.0; } // Mark the LUT as calculated flag = true; } // Return the releavant row in the look up table return binomial[order]; }

The loop up table for the control point indices can be prepared as follow:

static void CreateIdxCtrlLUT(CapyBezier* const that) { // Allocate memory for the look up table uint8_t** lut = NULL; safeMalloc(lut, sizeof(uint8_t*) * (that->nbCtrls)); loop (iCtrl, that->nbCtrls) safeMalloc(lut[iCtrl], sizeof(uint8_t) * that->dimIn); // Initialise the look up table uint64_t steps = 1; loop (iDim, that->dimIn) { uint64_t idxDim = 0; uint64_t iCtrl = 0; while (iCtrl < that->nbCtrls) { loop (iStep, steps) lut[iCtrl + iStep][iDim] = idxDim; iCtrl += steps; idxDim = (idxDim + 1) % (that->order + 1); } steps *= (that->order + 1); } // Set the look up table memcpy((uint8_t**)&(that->idxCtrlLUT), &lut, sizeof(uint8_t**)); }

Then, the CapyBezier object can be instanciated as follow:

CapyBezier CapyBezierCreate( uint8_t const order, uint8_t const dimIn, uint8_t const dimOut) { // Ensure there will be no overflow during calculation due to uint8_t assert(order < 254 && "CapyBezier.order must be less than 254"); assert(dimIn < 255 && "CapyBezier.dimIn must be less than 255"); assert(dimOut < 255 && "CapyBezier.dimOut must be less than 255"); // Create the Bezier CapyBezier that = { .order=order, .dimIn=dimIn, .dimOut=dimOut, .eval = Eval, .getCtrl = GetCtrl, .setCtrl = SetCtrl, }; uint64_t nbCtrls = GetNbCtrls(&that); memcpy((uint64_t*)&(that.nbCtrls), &nbCtrls, sizeof(uint64_t)); CapyArrVec arr = CapyVecAllocArr(dimOut, that.nbCtrls); memcpy(&(that.ctrls), &arr, sizeof(CapyArrVec)); that.ctrlWeights = CapyVecAlloc(that.nbCtrls); that.out = CapyVecAlloc(dimOut); double* binomial = CreateBinomialLUT(that.order); memcpy((double*)&(that.binomialLUT), &binomial, sizeof(double*)); CreateIdxCtrlLUT(&that); return that; }

and the evaluation method becomes:

static void Eval(CapyVec const* const x) { CapyThatBezier; // Reset the result values loop (iOut, that->dimOut) that->out->vals[iOut] = 0.0; // Loop on the control points #pragma omp parallel for loop (iCtrl, that->nbCtrls) { // Update the control point's weight according to the input value // in each dimension that->ctrlWeights->vals[iCtrl] = 1.0; loop (iIn, that->dimIn) { double t = x->vals[iIn]; uint8_t jCtrl = that->idxCtrlLUT[iCtrl][iIn]; double weight = that->binomialLUT[jCtrl] * pow(1.0 - t, that->order - jCtrl) * pow(t, jCtrl); that->ctrlWeights->vals[iCtrl] *= weight; } // Update the result values with the weighted contribution of the // control point CapyBezierCtrl* ctrl = $(&(that->ctrls), getPtr)(iCtrl); #pragma omp critical loop (iOut, that->dimOut) that->out->vals[iOut] += that->ctrlWeights->vals[iCtrl] * ctrl->vals[iOut]; } }