Least Squares Bezier Fit

$$ \{ {\bf C}_k \} =\text{argmin}(L) \\ L = \sum_{k=0}^{N-1} \left|  {\bf P}_k - \sum_{i=0}^n  b_{i,n} (t_k ) {\bf C} _i   \right|^2 $$

\begin{gather} t_k = d_k / d_{N-1}, \quad d_k = d_{k-1} + | {\bf P}_k -{\bf P}_{k-1}|\end{gather}

$$\text{Bernstein basis polynomial: } b_{i, n} (t)= \begin{pmatrix} n \\i \end{pmatrix} t^{i} (1-t)^{n-i} = \sum_{j=0}^n t^j  M_{ji} $$

$$ M_{ji} = \frac{n!}{(n-j)!} \frac{ (-1)^{j+i}}{ i! (j-i)!} =(-1)^{j+i}  M_{jj} \begin{pmatrix} j\\i  \end{pmatrix} \quad (j\ge i)$$

$$ \text{note: }~M_{jj} = \begin{pmatrix} n\\j \end{pmatrix},\qquad M_{ji}=0\quad (j<i)$$

// calculate binomial coefficients using Pascal's triangle
void binomialCoeffs(int n, double** C) {
    // binomial coefficients
    for (int i = 0; i <= n; i++)
        for (int j = 0; j <= i; j++)
            if (j == 0 || j == i) C[i][j] = 1;
            else                  C[i][j] = C[i - 1][j - 1] + C[i - 1][j];
}
std::vector<double> basisMatrix(int degree) {
    const int n = degree + 1;
    std::vector<double* > C(n);
    std::vector<double> Cbuff(n*n);
    for (int i = n; i-->0;) 
        C[i] = &Cbuff[i * n];
    // C[degree, k]; k=0,..,degree)
    binomialCoeffs(degree, &C[0]);
    // seting the diagonal;
    std::vector<double> M(n * n, 0); // lower triangle; 
    for (int j = 0; j <= degree; j++)
        M[j * n + j] = C[degree][j];
    // compute the remainings;
    for (int i = 0; i <= degree; i++) 
        for (int j = i + 1; j <= degree; j++)
            M[j * n + i] = ((j + i) & 1 ? -1 : 1) * C[j][i] * M[j * n + j];
    return M;
}