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

주어진 데이터를 잘 피팅하는 직선을 찾기 위해서는 데이터를 이루는 각 점의 $y$ 값과 같은 위치에서 구하려는 직선과의 차이 (residual)의 제곱을 최소화시키는 조건을 사용한다. 그러나 직선의 기울기가 매우 커지는 경우에는  데이터에 들어있는 outlier에 대해서는 그 차이가 매우 커지는 경우가 발생할 수도 있어 올바른 피팅을 방해할 수 있다. 이를 수정하는 간단한 방법 중 하나는 $y$값의 차이를 이용하지 않고 데이터와 직선 간의 최단거리를 이용하는 것이다. 

수평에서 기울어진 각도가 $\theta$이고 원점에서 거리가 $s$인 직선의 방정식은 

$$ x \sin \theta - y \cos \theta + s=0$$

이고, 한 점 $(x_i, y_i)$에서 이 직선까지 거리는

$$ d_i = | \sin \theta y_i - \cos \theta x_i + s|$$

이다. 따라서 주어진 데이터에서 떨어진 거리 제곱의 합이 최소인 직선을 구하기 위해서는 다음을 최소화시키는 $\theta, s$을 구해야 한다. $$ L = \sum_i \big( \sin \theta x_i - \cos \theta y_i +s \big)^2 \\ (\theta, s)=\text{argmin}(L)$$

$\theta$와 $s$에 대한 극값 조건에서 

$$\frac{1}{2} \frac{\partial L}{\partial \theta} = \frac{1}{2} \sin 2 \theta \sum_i (x_i^2 - y_i^2) - \cos 2 \theta \sum x_i y_i + s \sin \theta \sum_i x_i + s \cos \theta \sum_i y_i = 0$$

$$ \frac{1}{2}\frac{\partial L}{\partial s}=\cos \theta \sum y_i  - \sin \theta \sum_i x_i - N s=0$$

주어진 데이터의 질량중심계에서 계산을 수행하면 $\sum_i x_i = \sum_i y_i =0$ 이므로 데이터의 2차 모멘트를 $$ A= \sum_i (x_i^2 - y_i^2), \qquad B = \sum_i x_i y_i $$로 놓으면 직선의 파라미터를 결정하는 식은

$$ \frac{1}{2} A \sin 2 \theta   - B \cos 2 \theta = 0  \qquad \to \quad  \tan 2\theta = \frac{2B}{A} \\ s = 0$$

두 번째 식은 직선이 질량중심(질량중심계에서 원점)을 통과함을 의미한다. 첫번째 식을 풀면

$$ \tan \theta = \frac{- A \pm \sqrt{A^2 + (2B)^2 }}{2B}$$

두 해 중에서 극소값 조건을 만족시키는 해가 직선을 결정한다. 그런데

$$ \frac{1}{2}\frac{\partial^2 L}{\partial \theta^2}=  A \cos 2 \theta + 2B \sin 2 \theta = \pm \sqrt{A^2 + (2B)^2}  >0$$

이므로 위쪽 부호로 직선($x\sin \theta = y\cos \theta$)이 정해진다. 질량중심계에서는 원점을 지나지만 원좌표계로 돌아오면 데이터의 질량중심을 통과하도록 평행이동시키면 된다.

$$  \left(-A+ \sqrt{A^2+ (2B)^2} \right)  (x-\bar{x}) = 2B (y - \bar{y})  $$

여기서 주어진 데이터의 질량중심은 원좌표계에서

$$ \bar{x} = \frac{1}{N} \sum_i x_i, \quad \bar{y} = \frac{1}{N} \sum_i y_i$$

이다. 또한 원좌표계에서 $A$와 $B$의 계산은 

$$ A = \sum_i [ (x_i - \bar{x})^2 - (y_i - \bar{y})^2], \qquad B = \sum (x_i  - \bar{x})(y_i - \bar{y})$$ 

이 결과는 데이터 분포에 PCA를 적용해서 얻은 결과와 동일하다. PCA에서는 공분산 행렬의 고유값이 큰 쪽에 해당하는 고유벡터가 직선의 방향을 결정했다. https://kipl.tistory.com/211.  또한 통계에서는 Deming regression이라고 불린다.