Geometry & Recognition

Fixed-point version: 선분 길이가 4096 픽셀 정도까지는 1-pixel 이내의 오차로 그려진다.

void DDALine(CPoint A, CPoint B) {
    const int mulfac = 4096; // 2^12 = 4096;
    int dx = B.x - A.x;
    int dy = B.y - A.y;
    int adx = dx < 0 ? -dx : dx;
    int ady = dy < 0 ? -dy : dy;
    int steps = adx < ady ? ady : adx;
    if (steps == 0) SetPixel(A.x, A.y); //single point;
    else {
        int curX = A.x * mulfac;      
        int curY = A.y * mulfac;
        dx = (dx * mulfac) / steps;
        dy = (dy * mulfac) / steps;
        while (steps-- >= 0) {
            SetPixel(curX / mulfac, curY / mulfac);
            curX += dx;
            curY += dy;
        }
    }
}

void DDALine(int x0, int y0, int x1, int y1)  { 
    int dx = x1 - x0; 
    int dy = y1 - y0; 
    // calculate steps required for generating pixels 
    int steps = abs(dx) > abs(dy) ? abs(dx) : abs(dy); 
    float xinc = dx / (float)steps; 
    float yinc = dy / (float)steps; 
    float x = x0; 
    float y = x0; 
    for (int i = 0; i <= steps; i++) { 
        SetPixel(round(x), round(y));
        x += xinc;
        y += yinc;
    } 
}

$$\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));
}

입력점들을 부드럽게 연결하는 베지어 곡선을 찾아보자. 입력점 사이 구간을 3차 베지어 곡선으로 표현하려면  4개의 컨트롤점이 필요한데, 최소자승법을 이용하지 않고 순차적인 3 입력점 $P_0,P_1,P_2$이 주어질 때 이를 이용해서 컨트롤 점을 구성한다.(따라서 베지어 곡선은 시작과 끝이 아닌 입력점을 일반적으로 통과하지 않는다)  $$C_0=(P_0+P_1)/2$$

$$C_1=(P_0+5P_1)/6$$

$$C_2=(5P_1+P_2)/6$$

$$C_3=(P_1+P_2)/2$$

컨트롤점으로 하는 베지어 곡선은  $[P_0,P_1,P_2]$의 일부분을 표현한다. 이 과정은 다음 한 입력점을 포함시켜서 다시 반복을 하는데 이 떄  이전 베이어 곡선의 마지막 컨트롤점($C_2$)이 새로이 만들어지는 베지어 곡선의 시작점($C_0^{\text{new}}$)이 되므로 두 베지어 곡선은 부드럽게 연결이 된다. 입력점의 시작과 끝점은 베이어 곡선이 통과하도록 컨트롤점에 포함되게 설계한다.

$$\text{front}:C_0=P_0,C_1=(P_0+2P_1)/3$$

$$\text{back}:C_2=(2P_1+P_2)/3,C_3=P_3$$

그리고 주어진 베지어 곡선의 분할은 균일한 간격으로 하지 않고, 각각의 구간에서 베지어 곡선이 충분히 평탄하도록 분할해서 저장한다.

// test flatness of a bezier curve with control points {p1, c1, c2, p2};
bool FlatBezier(const double eps, 
                const CfPt& p1, const CfPt& c1, const CfPt& c2, const CfPt& p2) {
    CfPt U = c1 - (p1 * 2 + p2) / 3;
    CfPt V = c2 - (p1 + p2 * 2) / 3;
    double normU = hypot(U.x, U.y);
    double normV = hypot(V.x, V.y);
    return (3. * max(normU, normV) / 4.) <= eps;
}
// subdivision stage; De Casteljau's algorithm
void GetBezierPoints(const double eps, 
                     const CfPt &p1, const CfPt &c1, const CfPt &c2, const CfPt &p2,
                     std::vector<CfPt>& smoothed) {
    if (!FlatBezier(eps, p1, c1, c2, p2)) {
        // Subdivide the curve at t = 0.5
        // left;
        CfPt mid_segm = (p1  + c1 * 3 + c2 * 3 + p2) / 8.0; // = B(1/2);
        CfPt new_c1   = (p1 + c1) / 2.0;
        CfPt new_c2   = (p1 + c1 * 2 + c2) / 4.0;   // ((p1+c1)/2 + (c1+c2)/2)/2
        GetBezierPoints(eps, p1, new_c1, new_c2, mid_segm, smoothed);
        // right;
        new_c1  = (c1 + c2 * 2 + p2) / 4.0;
        new_c2  = (c2 + p2) / 2.0;
        GetBezierPoints(eps, mid_segm, new_c1, new_c2, p2, smoothed);
    } 
    else smoothed.push_back(p2); // flat enough;
        
}
// linear interpolation between a and b;
// a와 b 사이를  (1-t): t로 내분;
CfPt lerp(const double t, const CfPt &a, const CfPt &b) {
    return a * (1 - t) + b * t;
}
std::vector<CfPt> SmoothPathWithBezier(const std::vector<CfPt>& points) {
    if (points.size() < 3) return std::vector<CfPt> ();//null vector;
    CfPt cntl[4];
    CfPt *pts = &points[0];   // iterator;
    const int N = points.size();   
    std::vector<CfPt> smoothed;
    smoothed.push_back(pts[0]);    
    for (int i = 2; i < N; i++, pts++) {
        if (i == 2) {
            cntl[0] = pts[0];
            cntl[1] = lerp(2./ 3, pts[0], pts[1]);
        } else {
            //0번 control-pts= (0-1)의 중간점;
            //1번 control-pts= (0-1)의 5:1 내분점; 
            cntl[0] = lerp(1./ 2, pts[0], pts[1]);
            cntl[1] = lerp(5./ 6, pts[0], pts[1]);
        }
        if (i == N-1) {
            //1:2 내분점;
            cntl[2] = lerp(1./ 3, pts[1], pts[2]);
            cntl[3] = pts[2];
        } else {
            //2번 control-pts;(1-2)의 1:5 내분점;
            //3번 control-pts;(1-2)의 중간점;
            cntl[2] = lerp(1./ 6, pts[1], pts[2]);
            cntl[3] = lerp(1./ 2, pts[1], pts[2]);
        }
        //입력 중에 처음 두개가 같거나, 다음 두개가 같으면 마지막이 충분히 smooth하므로 
        // 마직막 control pt만 베지어 곡선에 삽입;
        if ((pts[0] == pts[1]) || (pts[1] == pts[2])) 
            smoothed.push_back(cntl[3]);
        else {
            double tolerance = 1.; // 1-pixel;
            GetBezierPoints(tolerance, cntl[0], cntl[1], cntl[2], cntl[3], smoothed);
        }
    }
    return smoothed;
}