Geometry & Recognition

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

$$ {\bf B}(t) = \sum _{k=0}^n b_{n,k}(t) {\bf Q}_k \\ b_{n,k}(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;
double Bezier(int n, double Q[], 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; 
}

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

void BezierSubdivision(int deg, CfPt *Q, 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 곡선이

$$ {\bf B}(t) = \sum_{k=0}^n b_{n,k}(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_{n-1,k}(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_{n-1,i}(t) {\bf Q}_{i+1} ,\qquad  {\bf B}_0(t) = \sum_{i=0}^{n-1} b_{n-1,i}(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;
}

$$ \{ {\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];
}
void basisMatrix(int degree, std::vector<double>& M) {
    const int n = degree + 1;
    double **C = new double* [n];
    C[0] = new double [n * n];
    for (int i = 1; i < n; i++) 
        C[i] = C[i-1] + n;
    // C[degree, k]; k=0,..,degree)
    binomialCoeffs(degree, C);
    // seting the diagonal;
    M.resize(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];
    
    delete [] C[0]; delete [] C;
}

주어진 데이터를 interpolation을 할 때 3차 spline을 많이 사용한다. 그럼 왜 3차 일까? 데이터를 연결하는 spline은 되도록이면  부드럽게 연결되어야 한다. 곡선이 부드럽게 그려지기 위해서는 급격한 꺾임이 없어야 하는데 얼마나 급격히 꺾이는가는 곡률에 비례하고, 곡률은 그 지점에서 함수의 두 번 미분값에 비례한다.(곡선의 곡률: https://kipl.tistory.com/105) 따라서 어떤 함수 $f(x)$가 구간 $[a,b]$에서 얼마나 부드럽게 연결되는가에 대한 척도는 곡선의 각부분에서 곡률의 크기의 제곱을 를 더한  다음 (음수가 아닌) 적분이 제공할 수 있다.

$$ \kappa[f] = \int_a^ b( f'' (x))^2 dx $$

$a \le x \le b$ 구간에서 균일하게 샘플링된 데이터가 있고 이들을 보간하는  3차 spline이 $S(x)=\{S_i(x) = a_i x^3 + b_i x^2 + c_i x +d_i\} $라고 하자. 또 두 번 이상 미분가능한 임의의 함수 $f(x)$도 주어진 샘플링 데이터를 지나가는 보간함수라고 하자. 이 경우 natural cubic spline이 주어진 sampling 데이터를 지나가면서 가장 부드럽게 이어지는 곡선임을 보일 수 있다.

$$ \kappa [f] \ge \kappa [S] \qquad \forall~ f$$

Spline $S$와 함수 $f$의 차이를 $h(x)= f(x) - S(x)$라 하면 각 node에서 $h(x_i)=0$이다. 이제 

$$ \kappa [f] = \int_a^b (h''(x) - S''(x))^2 = \kappa [h]  + \kappa [S]  - 2 \int_a^b h''(x) S''(s) dx $$

그런데, cubic spline의 경우 각 node에서 2차 미분이 연속이므로  

\begin{align} \int_a^b h''(x) S''(x) dx &= h'(x) S''(x)\Big|_a^b - \int _a^b h'(x) S'''(x)dx \\ &=h'(b)S''(b)-h'(a)S''(b) - \sum_i \int_{x_i}^{x_{i+1}} h'(x) S_i '''(x) dx \\  &=h'(b)S''(b)-h'(a)S''(a)- \sum_i \int_{x_i}^{x_{i+1} } h'(x) (6 a_i) dx \\ &= h'(b)S''(b)-h'(a)S''(a)- 6\sum_i a_i \left[  h(x_{i+1})- h(x_i) \right] \\&= h'(b)S''(b)-h'(a)S''(a) \end{align}  $$ \therefore  \quad \kappa [f] = \kappa[h]+\kappa[S]+ h'(b)S''(b)-h'(a) S''(a)$$

임의의 곡선 $h$에 대해서 $\kappa[h]\ge 0$이므로, 양끝에서 2차 미분이 $S''(a)=S''(b)=0$으로 고정된 natural cubic spline은  

$$ \kappa[f] \ge \kappa [S]$$

이므로 주어진 샘플링 데이터를 곡률척도가 가장 작게 즉, 가장 부드럽게 보간하는 곡선임을 알 수 있다.