Geometry & Recognition

주어진 control point로 만들어지는 Bezier 곡선을 계산할 때 De Casteljas 알고리즘을 이용하여 계산한다. 이 경우 control points을 백업하는 과정이 필요하다. 여기서는 직접적으로 Bernstein 함수를 계산하여 Bezier 곡선을 구하자. 우선 control point가 ${\bf Q}_i$일 때 $n$차 Bezier 곡선을 Bernstein 다항식을 써서 표현하면

$$ {\bf B}(t) = \sum _{k=0}^n b_{k, n}(t) {\bf Q}_k \\ b_{k, n}(t) = \begin{pmatrix} n \\k \end{pmatrix} t^k (1-t)^{n-k}$$

이므로 이항계수의 계산이 요구된다. 

$$ \begin{pmatrix} n \\k \end{pmatrix} = \frac{ n!}{k! (n-k)!} $$

이항계수는 미리 계산된 lookup table을 이용하는 경우가 더 편리할 수 있다: https://kipl.tistory.com/575

// n = degree of a Bezier curve == Q.size()-1;
double Bezier(const int n, const std::vector<double>& Q, const double t) { 
    double b = 0;
    double tk = 1; 
    double onetk = pow(1-t, double(n)); 
    for (int k = 0; k <= n; k++) { 
        double blend = tk * onetk;   // bernstein polynomial; 
        // for the next stage;
        tk *= t;              
        onetk /= (1-t);     
        // multiply binomial coefficient;
        int nf = n;           // n! 계산;
        int kf = k;           // k! 계산; 
        int nkf = n-k;        // (n-k)! 계산;
        while (nf >= 1) {
            blend *= nf;      // multiply by n!
            nf--;
            if (kf > 1) {     // divide by k!
                blend /= kf; 
                kf--; 
            } 
            if (nkf > 1) {    // divide by (n-k)!
                blend /= nkf; 
                nkf--; 
            } 
        } 
        b += Q[k] * blend; 
    } 
    return b; 
}

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

$$ {\bf B}(t) = \sum_{k=0}^n b_{k, n}(t) {\bf Q}_k = \sum_{k=0}^n\frac{n!}{k! (n-k)!}t^k (1-t)^{n-k} {\bf Q}_k$$

로 표현되므로 이의 미분은 

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

$$~~~~~~~- \sum_{k=0}^{n-1} \frac{n!}{k! (n-1-k)!} t^k (1-t)^{n-1-k}{\bf Q}_{k}$$

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

$$= \sum_{k=0}^{n-1} b_{k, n-1}(t) (n {\bf Q}_{k+1} - n{\bf Q}_k)$$

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

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

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

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

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

$${\bf B}_0(t) = \sum_{i=0}^{n-1} b_{i, n-1}(t) {\bf Q}_i  $$처럼 분해를 하면 ${\bf B}_1(t)$는 컨트롤점이 $\{\bf P_1, P_2,...,P_n\}$으로 구성된 $(n-1)$차 Bezier 곡선이고, ${\bf B}_0(t)$는 컨트롤점이 $\{ \bf P_0, P_1,...,P_{n-1}\}$으로 만들어지는 $(n-1)$차 Bezier 곡선이다. 따라서 Casteljau 알고리즘을 이용하여 ${\bf B}_1(t)$와 ${\bf 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;
}

Bezier 곡선은 control points $\{ \mathbf {P}_i\}$의 선형 결합으로 주어진다:

$$\mathbf {B}(t) = \sum_{i=0}^{n} B_{ i, n} (t) \mathbf {P}_i , \quad B_{i,n}(t)=\left(\begin {array}{c} n \\ i \end {array} \right) t^i (1-t)^{n-i}.$$

Bernstein 다항식 $B_{i, n}(t)$이 control points 선형 결합의 가중치를 역할을 하는데 0과 1 사이의 양의 실수 값을 가진다. 그리고 이들 가중치의 합은 1이다:

$$ 0\le  B_{ i, n}(t) \le 1, \quad i=0,1,2,... n ,    \quad    0\le t\le 1 $$

$$\sum_{i=0}^{n}  B_{i, n}(t) = 1$$

이는 Bezier 곡선이 control points가 만드는 convex region 내부에 있음을 의미한다. Bezier 곡선의 convexity 성질은 여러 좋은 특성을 준다.  몇 가지만 나열하면, 첫째가 Bezier 곡선은 항상 컨트롤 포인트의 convex hull 내에 놓이게 되므로 곡선의 제어 가능성을 보장한다. 둘째는 교차 여부를 쉽게 확인할 수 있다. 또한 culling을 가능하게 한다.