Least Squares Bezier Fit

$$\begin{gathered} \{{\mathbf{C}}_k\}=\text{argmin}(L) \\ L=\sum _{k=0}^{N-1}{|{\mathbf{P}}_k-\sum _{i=0}^n{b}_{i,n}(t_k){\mathbf{C}}_i|}^2 \end{gathered}$$

$$\begin{array}{c}{t}_k=d_k/d_{N-1},d_k=d_{k-1}+|{\mathbf{P}}_k-{\mathbf{P}}_{k-1}|\end{array}$$

$$\text{Bernstein basis polynomial: }{b}_{i,n}(t)=\left(\begin{array}{c}n\\ i\end{array}\right)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}\left(\begin{array}{c}j\\ i\end{array}\right)(j\ge i)$$

$$\text{note: }{M}_{jj}=\left(\begin{array}{c}n\\ j\end{array}\right),M_{ji}=0(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;
}