Geometry & Recognition

평면상에  주어진 점들을 보간하는 함수로써 Catmull-Rom spline을 사용할 때는 매개변수를 써서 표현하는 것보다는 $x$-좌표를 사용해서 표현하는 것이 더 편리할 수 있다. 이 경우 control point의 $x$-좌표 간격이 균일하지 않으므로 이를 고려해야 한다. 

• 주어진 구간에서 spline은 $x$의 삼차함수로  써진다:
  $$S_i(x)=a+b(x-x_i)+c(x-x_i{)}^2+d(x-x_i{)}^3,x_i\le x\le x_{i+1}$$
• 곡선은 모든 control point를 통과해야 한다:
  $$S_i(x_i)=y_i,S_i(x_{i+1})=y_{i+1}$$
• 1차 미분값은 양 이웃점의 기울기로 결정한다: 
  $$S_i^{\prime }(x_i)=\frac{y_{i+1}-y_{i-1}}{x_{i+1}-x_{i-1}},S_i^{\prime }(x_{i+1})=\frac{y_{i+2}-y_i}{x_{i+2}-x_i}$$ 

// Catmull-Rom spline interpolating {q[i]};
void CatmullRom(const std::vector<CfPt> &q, CDC *pDC) {
    if (q.size() < 2) return;
    for (int k = 0; k < q.size()-1; k++) {
        int km1 = max(k-1, 0);
        int kp2 = min(k+2, q.size()-1);

        /* compute spline coefficients */
        double dx  = q[k+1].x - q[k].x; 
        double dy  = (q[k+1].y - q[k].y) / dx;
        double dy1 = (q[k+1].y - [km1].y) / (q[k+1].x - q[km1].x);
        double dy2 = (q[kp2].y - q[k].y) / (q[kp2].x - q[k].x);
        double a = q[k].y;
        double b = dy1;
        double c = (3 * dy - 2 * dy1 - dy2) / dx;
        double d = (-2 * dy +   dy1 + dy2) / dx / dx;
#define STEPS 20
        if (k == 0) pDC->MoveTo(q[0]);
        const double deltax = q[k+1].x - q[k].x;
        for (double j = 1; j <= STEPS; j++)  { 
            double x = j * deltax / STEPS;
            /* use Horner's rule to calculate cubic polynomial */
            double y = ((d * x + c)*x + b)*x + a;
            pDC->LineTo(CfPt(q[k].x + x, y));
        }
    }
}

주어진 점들을 감싸는 최소  면적의 직사각형은 먼저 convex hull을 구한 후 rotating caliper을 써서 구할 수 있다. 다음 코드는 반시계방향으로 정렬된 convex hull을 이용해서 minimum volume box을 구한다. 처음 좌표축 방향으로 정렬된 bounding box을 이용해서 가장 작은 각을 형성하는 convex hull의 한 변을 찾고, 그 변으로 만들어지는 bounding box를 구한다. 반시계방향으로 bounding box가 convex hull의 변을 순차적으로 접하도록 회전시키면서 최소 면적인 경우를 찾으면 된다. 이는 convex hull을 바닥에서 굴리면서 bounding box을 찾는 경우로 생각해 보면 더 쉽다. bounding box가 convex hull의 모든 변을 다 순회하면 알고리즘이 종료되므로  time complexity는 $O(N)$이다.

void MinVolumeBox(const std::vector<CfPt>& chull, CfPt mvbox[4]) {
    if (chull.size() < 3) return;
    std::vector<CfPt> edges(chull.size());
    std::vector<int>  visited(chull.size(), 0);  // initialized = 0;
    enum {BB_NONE, BB_LEFT, BB_RIGHT, BB_BOTTOM, BB_TOP};
    for (int i = 0; i < chull.size(); i++) {
        const CfPt e = chull[(i+1) % chull.size()] - chull[i];
        edges[i] = e / e.Norm();// make a unit vector;
    }
    // 먼저 좌표의 최소-최대 값을 구한다;
    // index을 순환할 때, 항상 최대-최소가 (index+chull.size)%chull.size가 되게
    double xmin = chull[0].x, xmax = xmin;
    double ymin = chull[0].y, ymax = curr;
    int left = 0, bot = 0, right = 0, top = 0 ;
    for (int i = 1; i < chull.size(); i++) {
        if (chull[i].x <= xmin)      { xmin = chull[i].x; left = i;}
        else if (chull[i].x >= xmax) { xmax = chull[i].x; right = i;}
        if (chull[i].y <= ymin)      { ymin = chull[i].y; bot = i;}
        else if (chull[i].y >= ymax) { ymax = chull[i].y; top = i;}
    }
    if (chull[0].x <= xmin)       { xmin = chull[0].x; left = 0;}
    else if ( chull[0].x >= xmax) { xmax = chull[0].x; right = 0;}
    if (chull[0].y <= ymin)       { ymin = chull[0].y; bot = 0;}
    else if (chull[0].y >= ymax)  { ymax = chull[0].y; top = 0;}

    /*____Initial mvbox Setting____*/
    CfPt center = CfPt(xmin + xmax, ymin + ymax)/2;
    CfPt ax0 = CfPt(1,0), ax1 = CfPt(0,1);	// mvbox의 unit vectors;
    CfPt eh = ax0, ev = ax1;				// working unit vectors;
    double hscale  = (xmax - xmin) / 2; 
    double vscale  = (ymax - ymin) / 2;
    double minarea = hscale * vscale;		// area / 4;
    int done = 0;
    while (!done) {
        // edges[bot]는 chull[bot+1]-chull[bot];
        // edges[bot],edges[right],edges[top],edges[left] 중에서
        // 현재의 mvbox의 변(eh, ev,-eh,-ev)과 가장 작은 각(max_cosine)을 이루는 
        // edge를 찾음;
        int iflag = BB_NONE;
        double	max_cosine = 0;
        double	cosine = eh * edges[bot];
        if (cosine > max_cosine) { max_cosine = cosine; iflag = BB_BOTTOM;}

        cosine = ev * edges[right]; ;
        if (cosine > max_cosine) { max_cosine = cosine; iflag = BB_RIGHT;}
        
        cosine = -eh * edges[top]; 
        if (cosine > max_cosine) { max_cosine = cosine; iflag = BB_TOP;}
        
        cosine = -ev * edges[left]; 
        if (cosine > max_cosine) { max_cosine = cosine; iflag = BB_LEFT;}
        
        switch (iflag) {
        case BB_BOTTOM:
            if (!visited[bot]) {// edges[bot]가 mvbox의 한변에 포함됨;
                eh = edges[bot]; ev = eh.Rot90(); 
                // 다음 변으로 회전;
                visited[bot] = 1;
                if (++bot == chull.size()) bot = 0;
            } else done = 1;
            break;
        case BB_RIGHT:
            if (!visited[right]) {// edges[right]가 mvbox의 한변에 포함됨;
                ev = edges[right]; eh = -ev.Rot90(); 
                // 다음 변으로 회전;
                visited[right] = 1;
                if (++right == chull.size()) right = 0;
            } else done = 1;
            break;
        case BB_TOP:
            if (!visited[top]) { // edges[top]이 mvbox의 한변에 포함됨;
                eh = -edges[top]; ev = eh.Rot90(); 
                // 다음 변으로 회전;
                visited[top] = 1;
                if (++top == chull.size()) top = 0;
            } else done = 1;
            break;
        case BB_LEFT:
            if (!visited[left]) {// edges[left]가 mvbox의 한변에 포함됨;
                ev = -edges[left]; eh = -ev.Rot90(); 
                // 다음 변으로 회전;
                visited[left] = 1;
                if (++left == chull.size()) left = 0;
            } else done = 1;
            break;
        case BB_NONE: // polygon이 직사각형임;
            done = 1;
            break;
        }
        // check area of mvbox;
        double len0 = eh * (chull[right] - chull[left]) / 2;
        double len1 = ev * (chull[top] - chull[bot]) / 2;
        double area = len0 * len1; //면적의 1/4 임;
        if (area < minarea) {
            minarea	= area;
            ax0	= eh; ax1 = ev;
            hscale = len0; vscale = len1;
            // box center: tt = left 기준 top-bot을 연결하는 vector/2;
            CfPt tt = (chull[top] + chull[bot]) / 2 - chull[left];
            double zz = ev * tt; 
            center = eh * len0 + ev * zz + chull[left];
        }
    }
    eh = ax0 * hscale; // recover original scale;
    ev = ax1 * vscale; // recover original scale;
    mvbox[0] = center + eh + ev;
    mvbox[1] = center - eh + ev;
    mvbox[2] = center - eh - ev;
    mvbox[3] = center + eh - ev;
}

구간 $[a=t_0,...t_{n-1}=b]$에서  샘플링된 데이터 $\{f_0,...f_{n-1}\}$이 있다. $n-1$ 개의 각 구간에서 데이터를 보간하는 삼차다항식 함수 모임 $$S(t)=\{S_j(t)|t_j\le t<t_{j+1},j=0,1,2...,n-2\}$$을 찾아보자. 전체 구간에서 연속적이고 충분히 부드럽게 연결되기 위해서는 우선 각 node에서 연속이어야 하고, 또한 1차 도함수와 2차 도함수도 연속이 되어야 한다.  물론 각 node에서는 샘플링된 데이터값을 가져야 한다. 

$$\begin{array}{rl}(a)& S(t_j)=f_j\\ (b)& S_{j+1}(t_{j+1})=S_j(t_{j+1})\\ (c)& S_{j+1}^{\prime }(t_{j+1})=S_j^{\prime }(t_{j+1})\\ (d)& S_{j+1}^{″}(t_{j+1})=S_j^{″}(t_{j+1})\end{array}$$

$n-1$개 구간에서 각각 정의된  3차 함수를 결정하려면 $4\times (n-1)$개의 조건이 필요하다. (a)에서 $n$개, (b), (c), (d)에서 $3\times (n-2)$개의 조건이 나오므로 총 $n+3\times (n-2)=4n-6$개의 조건만 생긴다. 삼차식을 완전히 결정하기 위해서는 2개의 추가적인 조건을 부여해야 하는데 보통 양끝에서 2차 도함수값을 0으로 하거나(natural boundary) 또는 양끝에서 도함수 값을 특정한 값으로 고정시킨다(clamped boundary). 여기서는 양끝의 2차 도함수를 0으로 한 natural cubic spline만 고려한다. 그리고 $S_j(t)$가 $n-2$개의 구간에 대해서만 정의되어 있지만, 끝점을 포하는 구간에서도 정의하자. 이 경우 $S_{n-1}(t),t\ge t_{n-1}$에 부여된 조건은 $t_{n-1}$에서 $S_{n-2}(t)$와 연속, 미분연속 그리고 2차 도함수가 0인 조건만 부여된다.

$j-$번째 구간의 삼차함수를 

$$S_j(t)=a_j+b_j(t-t_j)+c_j(t-t_j{)}^2+d_j(t-t_j{)}^3$$

로 쓰고, 구간의 간격을 

$$h_j=t_{j+1}-t_j$$로 놓으면 (a)에서 

$$S_j(t_j)=a_j=f_i,j=0,1,..n-2(a)$$

(b)에서 

$$a_{j+1}=a_j+b_j{h}_j+c_j{h}_h^2+d_j{h}_j^3,j=0,1,...,n-2(b)$$

(c)에서 

$$b_{j+1}=b_j+2c_j{h}_j+3d_j{h}_j^2,j=0,1,...,n-2(c)$$

(d)에서 

$$c_{j+1}=c_j+3d_j{h}_j,j=0,1,...,n-2(d)$$

이므로  $d_j$에 대해서 푼 후 (b)와 (c)에 대입하면

$$a_{j+1}=a_j+b_j{h}_j+\frac{h_j^2}{3}(2c_j+c_{j+1})(e)$$

$$b_{j+1}=b_j+h_j(c_j+c_{j+1})(f)$$

을 얻는다. 그리고 (e)에서 $b_j$에 대해서 풀면

$$b_j=\frac{1}{h_j}(a_{j+1}-a_j)-\frac{h_j}{3}(2c_j+c_{j+1})$$

$$\to b_{j-1}=\frac{1}{h_{j-1}}(a_j-a_{j-1})-\frac{h_{j-1}}{3}(2c_{j-1}+c_j)$$

이 식을 index가 1씩 줄어든 (f)에 대입하면

$$h_{j-1}{c}_{j-1}+2(h_{j-1}+h_j)c_j+h_j{c}_{j+1}=\frac{3}{h_j}(a_{j+1}-a_j)-\frac{3}{h_{j-1}}(a_j-a_{j-1})$$$$j=1,2,...,n-2$$로 주어지는 $n$개의 계수 $c_j(j=0,1,...,n-1)$에 대한 선형연립방정식을 얻는다. 물론 매개변수 구간 간격 $h_j(j=0,1,...,n-2)$와  spline의 상수항 $a_j(j=0,1,2,...,n-2)$은 입력점으로 주어진다. 나머지 두 개의 조건은 양끝점에서 이차미분값을 0으로 설정하여 얻을 수 있다:

$$S^{″}(t_{n-1})=0\to c_{n-1}=0$$

$$S^{″}(t_0)=0\to c_0=0$$

연립방정식을 $\mathbf{c}=\{c_j\}$에 대한 행렬방정식 형태 $\mathbf{A}\cdot \mathbf{x}\mathbf{=}\mathbf{b}$로 정리하면 $n\times n$ 행렬 $\mathbf{A}$는 tridiagonal 모양이어서 쉽게 답을 구할 수 있다.

그림은 매개변수 구간 간격을 일정하게 잡았을 때(black)와 입력점 사이의 직선거리를 사용할 때(blue)의 차이를 보여준다.

std::vector<Cubic> cubicspline1(const std::vector<double>& a, 
                                const std::vector<double>& h) {
    const int n = a.size();
    std::vector<double> w(n), l(n), m(n), z(n);

    for (int i = n-1; i-->1;) 
        w[i] = 3 * (a[i + 1] - a[i]) / h[i] - 3 * (a[i] - a[i - 1]) / h[i - 1];
    l[0]= 1; m[0] = 0; z[0] = 0;
    
    for (int i = 1; i < n-1; i++) {
        l[i] = 2 * (h[i] + h[i - 1]) - h[i - 1] * m[i - 1];
        m[i] = h[i] / l[i];
        z[i] = (w[i] - h[i - 1] * z[i - 1]) / l[i];
    }
    l[n-1] = 1;
    z[n-1] = 0;
    std::vector<double> c(n);
    c[n-1] = 0;
    std::vector<Cubic> spline(n-1);
    for (int i = n-1; i-->0;) {
        c[i] = z[i] - m[i] * c[i + 1];
        double b = (a[i + 1] - a[i]) / h[i] - h[i] * (c[i + 1] + 2 * c[i]) / 3;
        double d = (c[i + 1] - c[i]) / (3 * h[i]);
        spline[i] = Cubic(a[i], b, c[i], d);
    }
    return spline;
}
// 매개변수가 균일하지 않는 경우임; 간격 = 입력점 직선거리;
std::vector<CPoint> NatCubicSpline(const std::vector<CPoint>& pts) {
    if (pts.size() < 2) return std::vector<CPoint> (); //null_vector;
    std::vector<double> ax(pts.size());
    std::vector<double> ay(pts.size());
    std::vector<double>  h(pts.size()-1); // interval size;
    for (int i = pts.size(); i-->0;) {
        ax[i] = pts[i].x; 
        ay[i] = pts[i].y;
    };
    // parameter intervals;
    for (int i = pts.size()-1; i-->0;) {
        double dx = ax[i + 1] - ax[i];
        double dy = ay[i + 1] - ay[i];
        h[i] = hypot(dx, dy);
    }

    std::vector<Cubic> splinex = cubicspline1(ax, h);
    std::vector<Cubic> spliney = cubicspline1(ay, h);

    const int steps = 20;
    std::vector<CPoint> curve;
    curve.reserve(1 + splinex.size()*steps);
    curve.push_back(pts.front());
    for (int i = 0; i < pts.size()-1; i++) {
        const double ds = h[i] / steps;
        for (int k = 1; k <= steps; k++) {
            const double tk = k * ds;
            curve.push_back(CPoint(int(splinex[i].eval(tk) + 0.5), 
                                   int(spliney[i].eval(tk) + 0.5)));
        }
    }	
    return curve;
}