3차의 Hermite spline을 쓰는 interpolation에서 입력 데이터 점에서의 기울기 값을 수정하면 입력 데이터의 monotonic 한 성질을 그대로 유지하는 보간을 할 수 있다. 주어진 데이터가 $\{(x_1, y_1), (x_2, y_2), \dots, (x_n, y_n)\}$이고, $x/y$-좌표가 증가하는 순으로 정렬되었다고 하자. 구간 $[x_k, x_{k+1} ]$의 시작과 끝에서 값 ($p_0 = y_k $, $p_1=y_{k+1}$)과 기울기 ($m_0, m_1$) 가 주어지면 Cubic Hermite spline 보간 함수는 다음과 같이 매개변수 $t \in [0,1]$의 함수로 표현된다: 

$$f(t) = p_0 h_{00}(t) + m_1 h_{10}(t) + p_1 h_{01}(t) + m_1 h_{11}(t),            \quad \quad t= \frac {x-x_k}{x_{k+1}-x_k} $$

$$f(0) = p_0,\quad f'(0) = m_0,\quad f(1) = p_1,\quad f'(1)=m_1$$

Hermite 기저 함수의 정의는 아래와 같다:

$$ h_{00}(t) = 2t^3 - 3t^2+1, \quad h_{10}(t) = t^3 - 2t^2+1,\quad h_{01}(t)=-2t^3+3t^2,\quad h_{11}(t)=t^3-t^2$$

 

보간된 spline이 단조성(monotonicity)을 유지하도록 각 입력점에서 기울기 값($m_{k}$)을 조절해주어야 한다. 그럼 기울기 $m_{k}$는 어떻게 설정해야 인접 구간의 보간 함수와 smooth 하게 연결될 수 있을까? 입력점의 각 위치에서 기울기는 각각 왼쪽과 오른쪽 인접 입력점과의 평균 변화율($\delta_{k-1}, \delta_{k}$)을 구한 후 다시 이들의 평균값을 기울기로 사용하면 된다. (양 끝점에서는 오른쪽이나 왼쪽 구간의 평균 변화율로 대체)

$$ m_k = \frac {1}{2}(\delta_{k-1}+\delta_{k} )= \frac {1}{2}\left(\frac { y_k-y_{k-1}}{x_k - x_{k-1}}+\frac {y_{k+1}-y_k}{x_{k+1} + x_k}\right)$$

이것만으로는 단조성을 보증을 할 수 없으므로 다음 조건을 추가로 부여해야 됨이 관련 논문에서 증명이 되었다. 

$$1. \text{For}~~\delta_{k} = 0 \quad\Longrightarrow\quad m_k =m_{k+1} = 0$$

$$2. \text{For}~~\epsilon=\sqrt {\left(\frac {m_k }{\delta_k }\right)^2 + \left(\frac {m_{k+1}}{\delta_{k} }\right)^2 } > 3\\ \quad\Longrightarrow\quad m_{k} \leftarrow3\frac {m_{k}}{\epsilon}, ~m_{k+1} \leftarrow 3\frac {m_{k+1}}{\epsilon}$$

 

아래의 코드는 위키피디아에 제시된 알고리즘을 구현한 것이다

(입력 데이터가 한 방향으로 monotonic 할 때만 동작함)

참고: http://en.wikipedia.org/wiki/Monotone_cubic_interpolation

void estimateTangents(std::vector <double>& x, std::vector <double>& y, std::vector<double>& m);

void monotoneCubicSplineInterp(std::vector<CPoint>& points, std::vector<CPoint>& curve) {
    curve.clear();
    if (points.size() < 2) return ;
    //if (!isMonotonicPolyline(points)) return ;
    std::vector<double> xx(points.size()), yy(points.size());
    for (int i = points.size(); i-->0;) {
        xx[i] = points[i].x ;
        yy[i] = points[i].y;
    };
    //tangent table;
    std::vector<double> m(points.size(), 0);
    estimateTangents(xx, yy, m);
    //hermite spline;
#define STEPS 12
    curve.push_back(points[0]); //add first point;
    for (int i = 0; i < points.size() - 1; i++) {
        double h = xx[i + 1] - xx[i];
        for (int k = 1; k <= STEPS; k++) {
            double t = double(k) / STEPS;         // 0 < t <= 1;
            double t2 = t * t ;
            double t3 = t2 * t ;
            double h00 =  2. * t3 - 3. * t2     + 1.;
            double h10 =       t3 - 2. * t2 + t;
            double h01 = -2. * t3 + 3. * t2;
            double h11 =       t3      - t2;
            double yf = h00 * yy[i] + h10 * h * m[i] + h01 * yy[i + 1] + h11 * h * m[i + 1];
            curve.push_back(CPoint(int(xx[i] + t * h + 0.5), int(yf + 0.5)));
        }
    }
};

Natural cubic spline과의 비교:

 

 

 

 

 

 

 

 
 
728x90

'Computational Geometry' 카테고리의 다른 글

Chain Hull  (2) 2012.09.16
Quick Hull  (2) 2012.09.16
Brute Force Convex Hull Algorithm  (0) 2012.09.06
Curvature of a Curve  (0) 2012.08.07
Douglas-Peucker Algorithm  (0) 2012.05.03
Posted by helloktk
,