Data Fitting with B-Spline Curves

$m$개의 control point $\{\mathbf{Q}_i\}$가 주어지면  $p$차의 B-Spline curve는 basis 함수를 $N_{i, p}(t)$을 써서

$$\mathbf{B}(t) = \sum_{i=0}^{m-1} N_{i, p}(t) \mathbf{Q}_i$$

로 표현할 수 있다. 이를 이용해서 일정한 순서로 샘플링된 평면상의 $N$개의 입력점 $\{ \mathbf{P}_i \}$을 찾아보자. B-spline 곡선을 이용하면 이 문제는 control 점을 찾는 문제가 된다. 곡선이 입력점을 잘 표현하기 위해서는 곡선과 입력점과의 차이를 최소로 하는 control 점을 찾아야 한다:

$$\mathbf{ Q}^* = \text{argmin}(L), \quad\quad L:= \sum_{k = 0}^{N-1} | \mathbf{B}(t_k) - \mathbf{P}_k|^2$$

여기서 $\{ t_k| k=0,1,...,N-1\}$는 입력점이 얻어지는 sample time으로 $0= t_0\le t_1\le...\le t_{N-1}= 1$로 rescale 할 수 있다. 

행렬을 이용해서 식을 좀 더 간결하게 표현할 수 있다. 

$$L = \sum_{k = 0}^{N-1} | \mathbf{B}(t_k) - \mathbf{P}_k|^2 = \sum_{k = 0}^{N-1} \left| \sum_{i=0}^{m-1} N_{i, p}(t_k) \mathbf{Q}_i -  \mathbf{P}_k \right|^2 = \sum_{k = 0}^{N-1} \left| \sum_{i=0}^{m-1} A_{ki} \mathbf{Q}_i - \mathbf{P}_k \right|^2, \\ \quad A_{ki} = N_{i, p}(t_k)$$

로 쓸 수 있으므로, $\hat {Q}= (\mathbf{Q}_0, \mathbf{Q}_1,..., \mathbf{Q}_{m-1})^t$, $\hat {P} = (\mathbf{P}_0, \mathbf{P}_1,..., \mathbf{P}_{N-1})^t$인 벡터, 그리고 $\hat {A} = (A_{ki})$인 행렬로 보면

$$L= \left| \hat{A} \cdot \hat {Q} - \hat {P} \right|^2 =  (\hat{Q}^t \cdot \hat {A}^t -\hat{P}^t ) \cdot (\hat{A}\cdot \hat{Q} - \hat{P}).$$

위 식의 값을 최소로 하는 최소자승해는 $\hat {Q}^t$에 대한 미분 값이 0인 벡터를 찾으면 된다; $$\frac {\partial L}{\partial \hat {Q}^t } = \hat {A}^t \cdot \hat {A} \cdot \hat {Q} - \hat {A}^{t} \cdot \hat {P} = 0.$$ 이 행렬 방정식의 해는

$$\hat {Q}^* = ( \hat {A} ^t \cdot \hat {A})^{-1} \cdot ( \hat {A}^t \cdot \hat {P})$$ 로 표현된다. $\hat {A} ^t \cdot \hat {A}$가 (banded) real symmetric ($m\times m$) 행렬이므로 Cholesky decomposion을 사용하면 쉽게 해를 구할 수 있다. $\hat{A}$가 banded matrix 형태를 가지는 이유는 basis가 local support에서만 0이 아닌 값을 취하기 때문이다. 

open b-spline(cubic)

std::vector<CfPt> BSplineFit_LS(const std::vector<CfPt>& data,
                                const int degree,    // cubic(3); 
                                const int nc)        // num of control points;
{
    // open b-spline;
    std::vector<double> knot((nc - 1) + degree + 2);
    for (int i = 0; i <= nc + degree; i++) knot[i] = i;
    
    std::vector<double> t(data.size());       // parameter;
    double scale = (knot[nc] - knot[degree]) / (data.size()-1);
    for (int i = data.size(); i-->0;) 
        t[i] = knot[degree] + scale * i;

    // design matrix;
    std::vector<double> A(ndata * nc); 
    for (int i = data.size(); i-->0;)
        for (int j = 0; j < nc; j++)
            A[i*nc + j] = Basis(j, degree, &knot[0], t[i]); //A(i,j)=N_j(t_i)

    // scattering matrix: S = A^t * A; real-symmetric matrix;
    std::vector<double> S(nc * nc); 
    for (int i = 0; i < nc; ++i)
        for (int j = i; j < nc; ++j) {
            double s = 0;
            for (int k = data.size(); k-->0;)
                s += A[k*nc + i] * A[k*nc + j];
            S[i*nc + j] = s;
        }
    // fill lower triangle;
    for (int i = 0; i < nc; ++i) 
        for (int j = 0; j < i; j++) 
            S[i*nc + j] = S[j*nc + i];
    
    std::vector<CfPt> X(nc);
    // X = A^t * P;
    for (int i = 0; i < nc; i++) {
        double sx = 0, sy = 0;
        for (int k = data.size(); k-->0;) {
            sx += A[k*nc + i] * data[k].x;
            sy += A[k*nc + i] * data[k].y;
        };
        X[i] =CfPt(sx, sy);
    };
    // psinv(A, n) = inverse of real symmetric matrix A;
    // ccmath-2.2.1 version;
    if (0 == psinv(&S[0], nc)) {
        std::vector<CfPt> Q(nc);
        for (int i = 0; i < nc; ++i) {
            souble sx = 0, sy = 0;
            for (int k = 0; k < nc; ++k) {
            	sx += S[i*nc + k] * X[k].x;
                sy += S[i*nc + k] * X[k].y;
            }
            Q[i] = CfPt(sx, sy);
        }
        return Q;  // estimated control points;
    }
    return std::vector<CfPt> ();
};

**네이버 블로그 이전;