B-Spline

$$\mathbf{C}(u)=\sum _{i=0}^n{N}_{i,p}(u){\mathbf{P}}_i$$

where $N_{i,p}(u)$'s are B-spline basis functions of degree $p$, $n+1$ control points $\{{\mathbf{P}}_i\}$, and a knot vector $\mathbf{U}=(u_0,u_1,...,u_m)$. Note, $m=n+p+1$.

B-spline의 basis 함수 $N_{i,p}(t)$는 유한한 구간에서만 0이 아닌 값을 갖는(local support) 함수로 다음과 같은  recursion relation을 가진다:

$$\begin{array}{rl}{N}_{i,0}(t)& =\{\begin{array}{cc}1& \mathrm{if}{u}_i\le x<u_{i+1}\\ 0& \mathrm{otherwise}\end{array}(0\le i\le n+p),\\ N_{i,p}(t)& =\frac{t-u_i}{u_{i+p}-u_i}{N}_{i,p-1}(t)+\frac{u_{i+p+1}-t}{u_{i+p+1}-u_{i+1}}{N}_{i+1,p-1}(t).\end{array}$$

$N_{i,p}(t)$의 local support는 $[u_i,u_{i+p+1}]$로 주어진다.

일반적으로 knot vector가 특별한 구조를 가지지 않으면 B-spline곡선은 Bezier 곡선의 경우와 달리 처음과 끝 control 점을 통과하지 않는다(open B-spline). 처음과 끝점을 통과하도록 만들기 위해서는 knot 벡터의 시작 $(p+1)$개 성분이 같은 값을 갖고, 마지막 $(p+1)$개 성분도 같은 값을 갖도록 조정하면 된다(clamped B-spline):

$$clamped:u[0]=u[1]=...=u[p],...,u[n+1]=...=u[n+p]=u[n+p+1]$$

knot는 균일할 필요는 없지만, 균일한 경우는 보통

$$u_i=\{\begin{array}{rcl}0& ,& 0\le i\le p,\\ \frac{i-p}{n+1-p}& ,& p+1\le i\le n\\ 1& ,& n+1\le i\le n+p+1\end{array}$$

로 선택한다. 

ref: ae.sharif.edu/~aerocad/curve%20modelling-B_Spline.pdf

clamped(black), open(blue)

// Note, p = degree;
void BSpline(int p, std::vector<CPoint>& control, int resolution, CDC *pDC) {
    const int n = control.size() - 1;  ASSERT(n > 0);
    // knot vector; 
    std::vector<double> u = calculateKnots(n, p);             // uniform knot;
    double delta = (u[n + p + 1] - u[0]) / (resolution - 1);  // parameter increment;
    for (double t = u[0] ; t <= u[n+p+1]; t += delta) {
        CPoint Q = calculatePoint(&u[0], n, p, t, &control[0]);
        if (t == u[0]) pDC->MoveTo(Q);
        pDC->LineTo(Q);
    }
    pDC->LineTo(control[n]);
}
// de Boor Cox's algorithm;
double Blend(int i, int p, double *u, double t) {  
    if (p == 0) { // termination condition of recursion;
        if ((u[i] <= t) && (t < u[i+1]))        return 1;
        else                                    return 0;
    }
    else {
        double coef1, coef2;
        if (u[i+p] == u[i]) {
            if (t == u[i])  coef1 = 1;
            else            coef1 = 0;
        } 
        else                coef1 = (t - u[i]) / (u[i+p] - u[i]); 
        
        if (u[i+p+1] == u[i+1]) {
           if (t == u[i+1]) coef2 = 1;
           else             coef2 = 0;
        } 
        else                coef2 = (u[i+p+1] - t) / (u[i+p+1] - u[i+1]);
        return coef1 * Blend(i, p-1, u, t) + coef2 * Blend(i+1, p-1, u, t);
    }
}
// clamped b-spline knot vector; uniform example;
// u[0]=u[1]= ... =u[p],.....,u=[n+p-2]=u[n+p-1]=u[n+p]=u[n+p+1]
std::vector<double> calculateKnots(int n, int p) {
    std::vector<double> u(n+p+2);
    int j = 0; 
    while (j <= p)     u[j++] = 0;                 // multiplicity = p+1;
    while (j <= n)     u[j++] = j - p;             // m = n+p+1;
    while (j <= n+p+1) u[j++] = n - p + 1;         // multiplicity = p+1;
    return u;
}
CPoint calculatePoint(double *u, int n, int p, double t, CPoint *control) {
    double x = 0, y = 0;
    for (int i = 0; i <= n; i++) {
        double b = Blend(i, p, u, t);
        x += control[i].x * b;
        y += control[i].y * b;
    }
    return CPoint(int(x + 0.5), int(y + 0.5));
}