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(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));
        }
    }
}

일정한 간격(예를 들면 1)으로 sampling 된 물체의 윤곽정보가 있다고 하자. 이 윤곽정보를 좀 더 적은 데이터를 이용해서 표현하려고 할 때 가장 쉬운 방법 중 하나는 샘플링 간격을 2배로 늘리는 것이다(down sampling). 이는 더 큰 스케일에서 물체를 바라보는 것과 동일한 효과를 준다. 그러나 down sampling에서 일부의 원정보가 사라지게 되므로 이에 대한 보정이 필요하게 된다. Down sampling 된 데이터에서 버려진 중간지점(홀수)의 정보는 선형적으로 예측할 수 있지만 이는 충분하지 못하다. 따라서 down sampling 할 때 버리지는 부분이 지니고 있는 고주파 성분을 저장한 후 이를 down sampling 된 저주파 성분의 보정에 사용해야 할 필요가 생긴다.

그러면 다운 샘플링과정에서 잃어버리는 고주파 성분은 어떻게 알아낼 수 있을까? 이는 다운 샘플링된 데이터로 원 데이터를 구성하려면 보간 방법을 사용하는데, 원 데이터에서 보간 값을 빼주면 근사적으로 예측할 수 있다. 다운 샘플링된 저주파 성분에 예측된 고주파 성분을 일정 부분 추가함으로써 더 정교한 다운샘플링을 할 수 있다. 이는 wavelet 변환에서 lifting-update 과정에 해당한다. 여기서는 이 알고리즘을 복잡한 폴리곤을 단순화시키는 과정에 적용해 본다. 보간 방법으로 Catmull-Rom spline을 이용했다. Douglas-Peucker 알고리즘에서는 변화지 않는 두 개의 고정점이 있지만, 이 방법에서는 모든 점이 변할 수 있다. 알고리즘을 반복 적용하면 매회 1/2로 다운샘플링되므로 어느 단계를 넘어서면 과도하게 단순화되는 단점이 있다. 입력점의 개수가 2의 지수승이 아닌 경우는 마지막 점을 반복해서 넣어주면 된다. 

// at t=1/2 (p0, p1, p2, p3);
CfPt catmull_rom(const CfPt& p0, const CfPt& p1, const CfPt& p2, const CfPt& p3) {
    return ((p1 + p2) * 9 - (p0 + p3)) / 16;
}
// 홀수점에서 고주파 성분을 주변 4개의 짝수점으로 만든 
// catmull-rom 스프라인을 이용해서 예측;
void lift(const std::vector<CfPt>& points,
          std::vector<CfPt>& hf_comp) {
    ASSERT(!(points.size() & 1) &&  points.size() >= 2);
    hf_comp.resize(points.size() / 2);
    for (int i = 1; i < points.size(); i += 2) {
        const int im3 = max(i-3, 0);
        const int im1 = i-1;
        const int ip1 = min(i+1, points.size()-1);
        const int ip3 = min(i+3, points.size()-1);
        // high freq comp = observed - estimated values;
        hf_comp[i/2] = points[i] - catmull_rom(points[im3], points[im1], points[ip1], points[ip3]);
    }
}
// update
// 고주파 성분 예측 에러가 필터링된 저주파 성분(짝수항)을 계산함
// 짝수항(저주파 성분)에 근방에서 구한 고주파 성분 평균의 절반을 더해줌;
void update(const std::vector<CfPt>& points,
            const std::vector<CfPt>& hf_comp,
            std::vector<CfPt>& lf_comp) {
    ASSERT(!(points.size() & 1));
    lf_comp.resize(points.size() / 2);
    for (int i = 0; i < points.size(); i += 2)
        // high[i/2]은 points 기준으로 (i+1)위치임;
        lf_comp[i/2] = points[i] + ((i==0 ? CfPt(0,0): hf_comp[i/2-1]) + hf_comp[i/2]) / 4;
}
std::vector<CfPt> decimate_polygon(const std::vector<CfPt>& points, const int max_pass) {
    const int orig_sz = points.size();
    int two_power = 1; 
    while (two_power < orig_sz) two_power <<= 1;
    // last element padding;
    std::vector<CfPt> input = points;
    const CfPt q_padd = input.back();
    for (int k = input.size(); k < two_power; k++) 
        input.push_back(q_padd);
    // wavelet passes
    std::vector<CfPt> hf_comp;
    std::vector<CfPt> lf_comp;
    // max_pass = max(max_pass-2, 0); // points.size()<= 2^max_pass;
    int pass = 0;
    while (pass++ < max_pass) {
        if (pass > 1) input.swap(lf_comp);
        lift(input, hf_comp);
        update(input, hf_comp, lf_comp);
    }   
    lf_comp.resize(orig_sz>>pass);
    return lf_comp; //decimated polyline;
};

https://kipl.tistory.com/99

Spline은 일정한 매개변수 구간 $[a, b]$에서 분포된 값(다차원인 경우는 벡터)을 이용해서 만든 다항식을 의미한다. 매개변수 구간 $[a, b]$는 값이 주어진 부분구간들의 합으로 표현할 수 있는데 이 부분구간의 시작과 끝에 해당하는 매개변수 값을 knot이라고 한다. knot이 매개변수 구간 $[a, b]$에 균일하게 퍼져 있는 경우는 uniform spline, 그렇지 않은 경우는 non-uniform spline이라고 한다.

Catmull-Rom spline은 3차 다항식으로 써지는 uniform spline으로 각 knot 위치에서 주어진 값을 통과하므로 interpolation 함수이기도 한다. 그리고 1차 도함수가 연속성을 가지고 있고, knot의 위치 변화나 knot에서 주어진 값의 변화에 곡선이 민감하게 변하지 않는다. 그러나 2차원 이상에는 Catmull-Rom spline이 꼬이는 현상이 발생할 수 있고, 극단적으로는 cusp가 생기는 단점도 가지고 있다.

2차원 이상에서는 Catmull-Rom spline으로 주어진 점(벡터)들을 보간하는 곡선을 찾을 때 uniform이 아니면 knot을 정하는데 모호함이 존재한다. 보통은 순차적으로 주어진 입력점의 직선거리를 이용해서 knot을 정할 수 있는데 이 경우를 chordal Catmull-Rom spline이라고 한다.

주어진 입력점이 $\{{\bf P}_0, {\bf P}_1, {\bf P}_2, {\bf P}_3\}$일 때 knot은

$$ t_i = t_{i-1}+ ||{\bf P}_i – {\bf P}_{i-1}||$$

더 일반적으로 

$$ t_i = t_{i-1} + ||{\bf P}_i - {\bf P}_{i-1}||^\alpha, ~~~\alpha \in [0,1]$$

로 잡을 수 있는데 이 경우를 centripetal Catmull-Rom spline이라 하며 cusp를 피할 수 있다. 단, 입력점의 중복이 없도록 미리 제거되어야 한다. Spline ${\bf S}(t)$는 ${\bf P}_1$과 ${\bf P}_2$사이에서

• 3차 함수로 쓰여져야 하고,
• 입력점 통과: ${\bf S}(t_1) = {\bf P}_1$, ${\bf S}(t_2) ={\bf P}_2$ ,
• 기울기: ${\bf S}'(t_1)=\frac{1}{t_2-t_0} ({\bf P}_2 - {\bf P}_0)$, $ {\bf S}'(t_2) = \frac{1}{t_3-t_1}({\bf P}_3 - {\bf P}_1)$

인 조건을 사용하면 ${\bf S}(t)$에 대한 closed form을 얻을 수 있다(참고: 위키피디아). ${\bf S}(t)$은 $t_1\le t\le t_2$ 사이에서 정의되고, 곡선을 그릴 때 매개변수를 $t\in [0:1]$ 으로 잡으면 ${\bf S}(t)$를 호출할 때는 $t_1 +  (t_2-t_1)t$를 인자로 사용하면 된다.

ref: http://www.cemyuksel.com/research/catmullrom_param/catmullrom.pdf

https://kipl.tistory.com/226

double GetParam(double t, double alpha, const CfPt& p0, const CfPt& p1 ) {
    CfPt p = p1 - p0;
    double dsq = p.x * p.x + p.y * p.y;
	if (dsq == 0) return t;
    return pow(a, alpha / 2) + t;
}
//centripetal catmull-rom spline;
CfPt CRSplineA(double t /*[0:1]*/, 
               const CfPt& p0, 
               const CfPt& p1, 
               const CfPt& p2, 
               const CfPt& p3, double alpha=.5 /*[0:1]*/)
{
    double t0 = 0;
    double t1 = GetParam( t0, alpha, p0, p1 );
    double t2 = GetParam( t1, alpha, p1, p2 );
    double t3 = GetParam( t2, alpha, p2, p3 );
    t = t1 * (1 - t) + t2 * t;
    CfPt A1 = t1==t0 ? p0: p0*(t1-t)/(t1-t0) + p1*(t-t0)/(t1-t0);
    CfPt A2 = t2==t1 ? p1: p1*(t2-t)/(t2-t1) + p2*(t-t1)/(t2-t1);
    CfPt A3 = t3==t2 ? p2: p2*(t3-t)/(t3-t2) + p3*(t-t2)/(t3-t2);
    CfPt B1 = A1 * (t2-t)/(t2-t0) + A2 * (t-t0)/(t2-t0);
    CfPt B2 = A2 * (t3-t)/(t3-t1) + A3 * (t-t1)/(t3-t1);
    return B1 * (t2-t)/(t2-t1) + B2 * (t-t1)/(t2-t1);
}
void DrawCatmullRom(std::vector<CfPt>& Q, CDC *pDC) {  
#define STEPS (20)
    if (Q.size() < 4) return ;
    CPen red(PS_SOLID, 2, RGB(255, 0, 0));
    CPen *pOld = pDC->SelectObject(&red);
    const int n = Q.size();
    // centripetal Catmull-Rom spline;
    for (int i = 0; i < n - 1; i++) {
        pDC->MoveTo(Q[i]);
        // i = 0 인 경우에는 처음 점을 반복, i = n - 2인 경우에는 마지막 점을 반복..
        int ip = max(i - 1, 0);
        int inn = min(i + 2, n - 1);
        for (int it = 1; it <= STEPS; it++)
            pDC->LineTo(CRSplineA(double(it)/STEPS, Q[ip], Q[i], Q[i + 1], Q[inn]));
    };
    pDC->SelectObject(pOld);
    CPen blue(PS_SOLID, 2, RGB(0, 0, 255));
    pOld = pDC->SelectObject(&blue);
    // uniform Catmull-Rom spline;
    for (int i = 0; i < n - 1; i++) {
        pDC->MoveTo(Q[i]);
        // i = 0 인 경우에는 처음 점을 반복, i = n - 2인 경우에는 마지막 점을 반복..
        int ip = max(i - 1, 0);
        int inn = min(i + 2, n - 1);
        for (int it = 1; it <= STEPS; it++)
            pDC->LineTo(CRSpline(double(it)/STEPS, Q[ip], Q[i], Q[i + 1], Q[inn]));
    };
    pDC->SelectObject(pOld);
}