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과의 비교:
'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 |