Catmull-Rom Spline (2)

Catmull-Rom Spline은 Edwin Catmull과 Raphael Rom이 제안한 interpolating spline으로, 주어진 제어점을 모두 통과하면서 매끄러운 곡선을 생성하는 방법이다. 컴퓨터 그래픽스, 애니메이션, 로봇 경로 계획, CAD 등에서 널리 사용된다.

평면상에  주어진 점들을 보간하는 함수로써 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(x_i) = \frac{y_{i+1}-y_{i-1}}{x_{i+1} - x_{i-1}},~ S'_{i}(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));
        }
    }
}