Catmull-Rom Spline (2)

평면상에  주어진 점들을 보간하는 함수로써 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));
        }
    }
}