Geometry & Recognition

$n$차 Bezier 곡선을 기술하는 매개변수가 $t$일 때, 곡선을 $[0,t]$인 구간과 $[t,1]$인 두 구간으로 나누기 위해 필요한 control point을 구해보자.

// deg = Q.size()-1;
void BezierSubdivision(const int deg, const std::vector<CfPt>& Q, const double t, 
                       std::vector<CfPt>& Left,
                       std::vector<CfPt>& Right) {
    Left.resize(deg+1);	
    Right.resize(deg+1);
    std::vector<CfPt> Q1(deg+1);
    for (int i = 0; i <= deg; i++)/* copy of control points*/
        Q1[i] = Q[i];

    /* triangle computation	*/
    Left[0] = Q1[0];
    Right[deg] = Q1[deg];
    for (int i = 0; i < deg; i++) {
        for (int j = 0; j < (deg-i); j++)
            Q1[j] = Q1[j] * (1 - t) + Q1[j+1] * t;
        Left[i+1] = Q1[0];
        Right[deg-1-i] = Q1[deg-1-i]; 
    }
}

한 점에서 Bezier 곡선까지의 최단거리나, Bezier 곡선상의 한 지점에서 접선 또는 법선을 구하기 위해서는 도함수를 구해야 할 필요가 생긴다. 그런데 주어진 찻수에서 Bezier 곡선의 도함수는 한 찻수 낮은 Bezier 곡선으로 표현할 수 있어서 상대적으로 쉽게 계산할 수 있다. 찻수가 $n$인 Bezier 곡선이

$$\mathbf{B}(t)=\sum _{k=0}^n{b}_{k,n}(t){\mathbf{Q}}_k=\sum _{k=0}^n\frac{n!}{k!(n-k)!}{t}^k(1-t{)}^{n-k}{\mathbf{Q}}_k$$

로 표현되므로 이의 미분은 

$$\frac{d\mathbf{B}(t)}{dt}=\sum _{k=1}^n\frac{n!}{(k-1)!(n-k)!}{t}^{k-1}(1-t{)}^{n-k}{\mathbf{Q}}_k$$

$$-\sum _{k=0}^{n-1}\frac{n!}{k!(n-1-k)!}{t}^k(1-t{)}^{n-1-k}{\mathbf{Q}}_k$$

$$=\sum _{k=0}^{n-1}\frac{(n-1)!}{k!(n-1-k)!}{t}^k(1-t{)}^{n-1-k}(n{\mathbf{Q}}_{k+1}-n{\mathbf{Q}}_k)$$

$$=\sum _{k=0}^{n-1}{b}_{k,n-1}(t)(n{\mathbf{Q}}_{k+1}-n{\mathbf{Q}}_k)$$

따라서 $n$ Bezier 곡선의 미분은 control 점이

$${\tilde{\mathbf{Q}}}_k=n({\mathbf{Q}}_{k+1}-{\mathbf{Q}}_k)k=0,1,..,n-1$$

로 주어지는 $(n-1)$차 Bezier 곡선으로 표현된다. 미분을 

$$\frac{d\mathbf{B}(t)}{dt}=n[{\mathbf{B}}_1(t)-{\mathbf{B}}_0(t)]$$

$${\mathbf{B}}_1(t)=\sum _{i=0}^{n-1}{b}_{i,n-1}(t){\mathbf{Q}}_{i+1}$$

$${\mathbf{B}}_0(t)=\sum _{i=0}^{n-1}{b}_{i,n-1}(t){\mathbf{Q}}_i$$처럼 분해를 하면 ${\mathbf{B}}_1(t)$는 컨트롤점이 $\{{\mathbf{P}}_{\mathbf{1}}\mathbf{,}{\mathbf{P}}_{\mathbf{2}}\mathbf{,}\mathbf{.}\mathbf{.}\mathbf{.}\mathbf{,}{\mathbf{P}}_{\mathbf{n}}\}$으로 구성된 $(n-1)$차 Bezier 곡선이고, ${\mathbf{B}}_0(t)$는 컨트롤점이 $\{{\mathbf{P}}_{\mathbf{0}}\mathbf{,}{\mathbf{P}}_{\mathbf{1}}\mathbf{,}\mathbf{.}\mathbf{.}\mathbf{.}\mathbf{,}{\mathbf{P}}_{\mathbf{n}\mathbf{-}\mathbf{1}}\}$으로 만들어지는 $(n-1)$차 Bezier 곡선이다. 따라서 Casteljau 알고리즘을 이용하여 ${\mathbf{B}}_1(t)$와 ${\mathbf{B}}_0(t)$을 구하여 그 차이를 계산하면 미분값을 얻을 수 있다.

곡선의 curvature: https://kipl.tistory.com/105

// Bezier cureve evaluation at t;
// deg = the degree of the bezier curve;
double Bezier(int deg, double Q[], double t) {
    if (deg==0) return Q[0];
    else if (deg==1) return (1 - t) * Q[0] + t * Q[1];
    else if (deg==2) return (1-t)*((1-t)*Q[0] + t*Q[1]) + t*((1-t)*Q[1] + t*Q[2]);
    
    std::vector<double> Q1(deg + 1);
    for (int i = 0; i <= deg; i++)  Q1[i] = Q[i];
    // triangle computations;
    for (int k = 0; k < deg; k++)
        for (int j = 0; j < (deg - k); j++)
            Q1[j] = (1 - t) * Q1[j] + t * Q1[j + 1];
            
    return Q1[0];
 }
// derivative of a Bezier curve at t;
double BezierDerivative(int deg, double Q[], double t) {
    if (deg==0) return 0;
    else if (deg==1) return Q[1] - Q[0];
    else if (deg==2) return 2*((1-t)*(Q[1]-Q[0])+t*(Q[2]-Q[1]));
    
    std::vector<double> Q1(degree + 1);
    for (int i = 0; i <= deg; i++) Q1[i] = Q[i];
    // triangle computations;
    for (int k = 0; k < (deg - 1); k++)
        for (int j = 0; j < (deg - k); j++)
            Q1[j] = (1 - t) * Q1[j] + t * Q1[j + 1];
    
    return deg * (Q1[1] - Q1[0]);
};
// 2nd derivative of a Bezier curve at t;
double EvalBezier2ndDeriv(int deg, double Q[], double t) {
    if (2 > deg) return 0;
    std::vector<double> Q1(deg+1);
    for (int i = 0; i <= deg; i++) Q1[i] = Q[i];
    
    for (int i = 0; i < (deg-2); i++) 
        for (int j = 0; j < (deg-i); j++) 
            Q1[j] = (1-t)*Q1[j] + t*Q1[j+1];
    
    double v0 = 2*(Q1[1] - Q1[0]);
    double v1 = 2*(Q1[2] - Q1[1]);
    return v1 - v0;
}; 
// curvature of a Bezier curve at t;
double BezierCurvature(int deg, CfPt Q[], double t) {
    if (deg < 2) return 0;
    CfPt d1 = EvalBezierDeriv(deg, Q, t);
    CfPt d2 = EvalBezier2ndDeriv(deg, Q, t);
    double flen = hypot(d1.x, d1.y);
    return fabs(d1.x*d1.y - d2.x*d1.y)/ flen / flen / flen;
}

$$\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;
}