'Natural Cubic Spline' 태그의 글 목록

Smoothing Spline

Computational Geometry 2024. 6. 29. 21:27

주어진 관측값 ${x i, y i} <math xmlns="http://www.w3.org/1998/Math/MathML"><mo fence="false" stretchy="false">{</mo><msub><mi>x</mi><mi>i</mi></msub><mo>,</mo><msub><mi>y</mi><mi>i</mi></msub><mo fence="false" stretchy="false">}</mo></math>$ 을 spline $y = S (x) <math xmlns="http://www.w3.org/1998/Math/MathML"><mi>y</mi><mo>=</mo><mi>S</mi><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo></math>$ 을 이용해서 모델링을 하자. 주어진 관측위치 $x i <math xmlns="http://www.w3.org/1998/Math/MathML"><msub><mi>x</mi><mi>i</mi></msub></math>$ 에서 spline은 관측값에 충분히 가까이 지나야 하고, 더블어 관측과 관측 사이에서는 너무 급격히 변하지 않는 곡선이어야 한다. 곡선의 변동은 곡률에 의해서 결정되므로 급격히 변하지 않기 위해서는 spline의 이차 미분값의 크기가 되도록이면 작아한다. 따라서 구하려고 하는 spline은 다음을 최소화시키는 조건에서 찾을 수 있을 것이다.

$L=α∑i(yi−S(xi)σi)2+β∫xnx0dx(S″<math xmlns="http://www.w3.org/1998/Math/MathML" display="block"><mi>L</mi><mo>=</mo><mi>α</mi><munder><mo data-mjx-texclass="OP">∑</mo><mi>i</mi></munder><msup><mrow data-mjx-texclass="INNER"><mo data-mjx-texclass="OPEN">(</mo><mfrac><mrow><msub><mi>y</mi><mi>i</mi></msub><mo>−</mo><mi>S</mi><mo stretchy="false">(</mo><msub><mi>x</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow><msub><mi>σ</mi><mi>i</mi></msub></mfrac><mo data-mjx-texclass="CLOSE">)</mo></mrow><mn>2</mn></msup><mo>+</mo><mi>β</mi><msubsup><mo data-mjx-texclass="OP">∫</mo><mrow data-mjx-texclass="ORD"><msub><mi>x</mi><mn>0</mn></msub></mrow><mrow data-mjx-texclass="ORD"><msub><mi>x</mi><mi>n</mi></msub></mrow></msubsup><mi>d</mi><mi>x</mi><mo stretchy="false">(</mo><msup><mi>S</mi><mo data-mjx-alternate="1">″</mo></msup><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo><msup><mo stretchy="false">)</mo><mn>2</mn></msup><mspace linebreak="newline"></mspace><mi>S</mi><mo>=</mo><munder><mtext>argmin</mtext><mrow></mrow></munder><mtext> </mtext><mi>L</mi></math>$

여기서 $σ_{i}$ 는 관측오차에 해당한다. $L$ 의 두 항은 서로 반대되는 경향으로 spline을 변화시키는데, 첫번째 항이 기여가 큰 경우는 spline이 주어진 관측값의 보간 곡선에 가까워지므로 곡선이 급격히 변할 수 있다. 두번째 항의 기여가 커지는 경우는 spline은 주어진 관측값을 통과하지는 못하지만 충분히 부드러워지게 만든다. 그리고 극단적으로 이 값만 기여하는 경우는 spline은 직선이 되게 된다. 두 항의 상대적인 가중치를 설정하기 위해서 $β = 1 - α$ 로 놓을 수 있지만, 첫번째 항과 두 번째 항의 물리적인 차원이 다르므로 단순 비교는 정확하지 않다. 첫번쨰 항은 차원이 없는 반면에 두 번째 항은 거리의 역수 차원을 가지므로 비교를 하기 위해서는 적당한 거리 척도를 곱해 $(1 - α) \times distance scale$ 로 하는 것이 맞다. 보통 $α = 1$ 로 잡고, $β$ 를 $0$ (interpolation에 가까워짐)에서부터 점점 큰 값으로 (직선에 가까워짐)으로 변화시키면 된다.

측정값 사이 구간에서 $n - 1$ 개의 삼차 spline

$S (x) = {S_{i} (x) = a_{i} x^{3} + b_{i} x^{2} + c_{i} x + d_{i}, x_{i} \leq x \leq x_{i + 1} | i = 0, 1, 2, . . ., n - 1}$

을 이용해서 smoothing spline을 구하자. 이 경우 한 구간에서 spline의 미분은 인접 spline 2차 계수의 선형보간으로 표현할 수 있다.

$S_{i}^{″} (x) = (1 - x / h_{i}) b_{i} + (x / h_{i}) b_{i + 1}, h_{i} = x_{i + 1} - x_{i}$

따라서 $L$ 은 다음과 같이 표현된다.

$L = α \sum_{i = 0}^{n} {(\frac{y_{i} - d_{i}}{σ_{i}})}^{2} + β \sum_{i = 0}^{n - 1} \frac{4 h_{i}}{3} (b_{i}^{2} + b_{i} b_{i + 1} + b_{i + 1}^{2})$

따라서 smoothing spline을 찾는 문제는 $L$ 을 최소화 시키는 $b_{i}, d_{i}$ 을 찾는 문제로 환원되었다. 나머지 계수 $c_{i}$ 는 spline의 연속성, $a_{i}$ 는 이차 도함수의 연속성을 요구하면 구할 수 있다.

$a_{i} = \frac{b_{i + 1} - b_{i}}{3 h_{i}} c_{i} = \frac{d_{i + 1} - d_{i}}{h_{i}} - \frac{1}{3} (b_{i + 1} - 2 b_{i}) h_{i}$

여기에 1차 도함수의 연속성을 요구하면 $b_{i}$ 와 $d_{i}$ 사이의 recursion relation을 구할 수 있으므로 최종적으로 $b_{i}$ 또는 $d_{i}$ 을 구하는 문제로 환원이 된다.

$b_{i - 1} h_{i - 1} + 2 b_{i} (h_{i - 1} + h_{i}) + b_{i + 1} h_{i} = \frac{3}{h_{i}} (d_{i + 1} - d_{i}) - \frac{3}{h_{i - 1}} (d_{i} - d_{i - 1})$

계속...

Ref: https://monoceros.physics.muni.cz/~jancely/NM/Texty/Numerika/CubicSmoothingSpline.pdf

wage.xlsx

0.06MB

나이에 따른 시급 분포: $σ_{i}$ 는 각 나이에서 시급의 표준편차(데이터가 1인 경우는 전체 편차의 평균을 사용함)을 때 평균값을 기준으로 추정한 곡선이다.

struct Spline {
    double x, y;
    double a, b, c, d;
    double eval(double t) {
        return d + (c + (b + a*t)*t)*t;
    }
    void toBezier(const double h, CPoint Q[4]) {
        double xp[4], yp[4];
        toBezier(h, xp, yp);
        for (int i = 0; i < 4; i++) 
            Q[i] = CPoint(int(xp[i]+.5), int(yp[i]+.5));
    }
    void toBezier(const double h, double xp[4], double yp[4]) {
        xp[0] = x; 
        xp[1] = x + h/3;
        xp[2] = x + 2*h/3;
        xp[3] = x + h;
        yp[0] = d;
        yp[1] = d + c * h / 3;
        yp[2] = d + (2 * c + b * h) * h / 3;
        yp[3] = d + (c + (b + a * h) * h) * h; 	
    }
};
void solveFiveDiag(std::vector<double> &u,
                   std::vector<double> &v,
                   std::vector<double> &w,
                   std::vector<double> &q) 
{
    const int n = u.size();
    u[0] = 0; //u[1]은 외부에서;
    v[1] = v[1]/u[1];
    w[1] = w[1]/u[1];
    for (int j = 2; j < n; j++) {
        u[j] = u[j] - u[j-2] * w[j-2] * w[j-2] - u[j-1] * v[j-1] * v[j-1];
        v[j] = (v[j] - u[j-1]*v[j-1]*w[j-1])/u[j];
        w[j] = w[j] / u[j];
    }
    //forward substitution
    q[1] = q[1] - v[0] * q[0];
    for (int j = 2; j < n; j++)
        q[j] = q[j] - v[j-1] * q[j-1] - w[j-2] * q[j-2];

    for (int j = 1; j < n; j++) 
        q[j] = q[j] / u[j];

    // back substitution
    q[n] = 0;
    // q[n-1]은 외부에서 정해짐;
    for (int j = n-2; j > 0; j--) 
        q[j] = q[j] - v[j] * q[j+1] - w[j]*q[j+2];
}
std::vector<Spline> SmoothingSpline1(
    std::vector<double> &x, //n+1
    std::vector<double> &y, //n+1
    std::vector<double> &sigma, //n+1
    double alpha, double beta)
{
    const int n = x.size()-1;
    std::vector<Spline> S(x.size());
    for (int i = x.size(); i-->0;)
        S[i].x = x[i], S[i].y = y[i];

    std::vector<double> h(n), r(x.size()), q(x.size()), u(n), v(n), w(n);
    h[0] = x[1] - x[0];
    r[0] = 3 / h[0];
    r[n] = 0;
    for (int i = 1; i < n; i++) {
        h[i] = x[i+1] - x[i];
        r[i] = 3/h[i];
        q[i] = (y[i+1] - y[i])*r[i] - (y[i] - y[i-1])*r[i-1];
    }
    //
    double mu = 2 * beta / (3 * alpha);
    for (int i = 1; i < n; i++) {
         double ui = SQR(r[i-1]) * sigma[i-1]
               + SQR(r[i-1] + r[i]) * sigma[i] + SQR(r[i]) * sigma[i+1];
         u[i] = mu * ui + 2 * (x[i+1] - x[i-1]);
         double vi = -(r[i-1] + r[i]) * r[i] * sigma[i] 
                - r[i] * (r[i] + r[i+1]) * sigma[i+1];
         v[i] = mu * vi + h[i];
         w[i] = mu * r[i] * r[i+1] * sigma[i+1];
    }
    solveFiveDiag(u, v, w, q);
    S[0].d = S[0].y - mu * r[0] * q[1] * sigma[0];
    double S1d = S[1].y - mu * ((-r[0] - r[1]) * q[1] + r[1] * q[2]) * sigma[1];
    S[0].a = q[1] / (3 * h[0]);
    S[0].b = 0;
    S[0].c = (S1d - S[0].d) / h[0] - q[1] * h[0] / 3;
    
    for (int j = 1; j < n; j++) {
        S[j].a = (q[j+1] - q[j]) / (3 * h[j]);
        S[j].b = q[j];
        S[j].c = (q[j] + q[j-1]) * h[j-1] + S[j-1].c;
        S[j].d = r[j-1] * q[j-1] - (r[j-1] + r[j]) * q[j] + r[j] * q[j+1];
        S[j].d = y[j] - mu * S[j].d * sigma[j];
    }
    return S;
}

저작자표시 비영리 변경금지

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

Alpha Shape (1)	2024.07.14
Approximate Convex Hull (0)	2024.06.29
Catmull-Rom Spline (2) (0)	2024.06.21
Minimum Volume Box (0)	2024.06.16
Natural Cubic Spline: revisited (0)	2024.06.14

Natural Cubic Spline: revisited

Computational Geometry 2024. 6. 14. 18:47

구간 $[a = t_{0}, . . . t_{n - 1} = b]$ 에서 샘플링된 데이터 ${f_{0}, . . . f_{n - 1}}$ 이 있다. $n - 1$ 개의 각 구간에서 데이터를 보간하는 삼차다항식 함수 모임 $S (t) = {S_{j} (t) | t_{j} \leq t < t_{j + 1}, j = 0, 1, 2. . ., n - 2}$ 을 찾아보자. 전체 구간에서 연속적이고 충분히 부드럽게 연결되기 위해서는 우선 각 node에서 연속이어야 하고, 또한 1차 도함수와 2차 도함수도 연속이 되어야 한다. 물론 각 node에서는 샘플링된 데이터값을 가져야 한다.

$\begin{aligned} (a) & S (t_{j}) = f_{j} \\ (b) & S_{j + 1} (t_{j + 1}) = S_{j} (t_{j + 1}) \\ (c) & S_{j + 1}^{'} (t_{j + 1}) = S_{j}^{'} (t_{j + 1}) \\ (d) & S_{j + 1}^{″} (t_{j + 1}) = S_{j}^{″} (t_{j + 1}) \end{aligned}$

$n - 1$ 개 구간에서 각각 정의된 3차 함수를 결정하려면 $4 \times (n - 1)$ 개의 조건이 필요하다. (a)에서 $n$ 개, (b), (c), (d)에서 $3 \times (n - 2)$ 개의 조건이 나오므로 총 $n + 3 \times (n - 2) = 4 n - 6$ 개의 조건만 생긴다. 삼차식을 완전히 결정하기 위해서는 2개의 추가적인 조건을 부여해야 하는데 보통 양끝에서 2차 도함수값을 0으로 하거나(natural boundary) 또는 양끝에서 도함수 값을 특정한 값으로 고정시킨다(clamped boundary). 여기서는 양끝의 2차 도함수를 0으로 한 natural cubic spline만 고려한다. 그리고 $S_{j} (t)$ 가 $n - 2$ 개의 구간에 대해서만 정의되어 있지만, 끝점을 포하는 구간에서도 정의하자. 이 경우 $S_{n - 1} (t), t \geq t_{n - 1}$ 에 부여된 조건은 $t_{n - 1}$ 에서 $S_{n - 2} (t)$ 와 연속, 미분연속 그리고 2차 도함수가 0인 조건만 부여된다.

$j -$ 번째 구간의 삼차함수를

$S_{j} (t) = a_{j} + b_{j} (t - t_{j}) + c_{j} (t - t_{j})^{2} + d_{j} (t - t_{j})^{3}$

로 쓰고, 구간의 간격을

$h_{j} = t_{j + 1} - t_{j}$ 로 놓으면 (a)에서

$S_{j} (t_{j}) = a_{j} = f_{i}, j = 0, 1, . . n - 2 (a)$

(b)에서

$a_{j + 1} = a_{j} + b_{j} h_{j} + c_{j} h_{h}^{2} + d_{j} h_{j}^{3}, j = 0, 1, . . ., n - 2 (b)$

(c)에서

$b_{j + 1} = b_{j} + 2 c_{j} h_{j} + 3 d_{j} h_{j}^{2}, j = 0, 1, . . ., n - 2 (c)$

(d)에서

$c_{j + 1} = c_{j} + 3 d_{j} h_{j}, j = 0, 1, . . ., n - 2 (d)$

이므로 $d_{j}$ 에 대해서 푼 후 (b)와 (c)에 대입하면

$a_{j + 1} = a_{j} + b_{j} h_{j} + \frac{h_{j}^{2}}{3} (2 c_{j} + c_{j + 1}) (e)$

$b_{j + 1} = b_{j} + h_{j} (c_{j} + c_{j + 1}) (f)$

을 얻는다. 그리고 (e)에서 $b_{j}$ 에 대해서 풀면

$b_{j} = \frac{1}{h_{j}} (a_{j + 1} - a_{j}) - \frac{h_{j}}{3} (2 c_{j} + c_{j + 1})$

$\to b_{j - 1} = \frac{1}{h_{j - 1}} (a_{j} - a_{j - 1}) - \frac{h_{j - 1}}{3} (2 c_{j - 1} + c_{j})$

이 식을 index가 1씩 줄어든 (f)에 대입하면

$h_{j - 1} c_{j - 1} + 2 (h_{j - 1} + h_{j}) c_{j} + h_{j} c_{j + 1} = \frac{3}{h_{j}} (a_{j + 1} - a_{j}) - \frac{3}{h_{j - 1}} (a_{j} - a_{j - 1})$ $j = 1, 2, . . ., n - 2$ 로 주어지는 $n$ 개의 계수 $c_{j} (j = 0, 1, . . ., n - 1)$ 에 대한 선형연립방정식을 얻는다. 물론 매개변수 구간 간격 $h_{j} (j = 0, 1, . . ., n - 2)$ 와 spline의 상수항 $a_{j} (j = 0, 1, 2, . . ., n - 2)$ 은 입력점으로 주어진다. 나머지 두 개의 조건은 양끝점에서 이차미분값을 0으로 설정하여 얻을 수 있다:

$S^{″} (t_{n - 1}) = 0 \to c_{n - 1} = 0$

$S^{″} (t_{0}) = 0 \to c_{0} = 0$

연립방정식을 $c = {c_{j}}$ 에 대한 행렬방정식 형태 $A \cdot x = b$ 로 정리하면 $n \times n$ 행렬 $A$ 는 tridiagonal 모양이어서 쉽게 답을 구할 수 있다.

그림은 매개변수 구간 간격을 일정하게 잡았을 때(black)와 입력점 사이의 직선거리를 사용할 때(blue)의 차이를 보여준다.

std::vector<Cubic> cubicspline1(const std::vector<double>& a, 
                                const std::vector<double>& h) {
    const int n = a.size();
    std::vector<double> w(n), l(n), m(n), z(n);

    for (int i = n-1; i-->1;) 
        w[i] = 3 * (a[i + 1] - a[i]) / h[i] - 3 * (a[i] - a[i - 1]) / h[i - 1];
    l[0]= 1; m[0] = 0; z[0] = 0;
    
    for (int i = 1; i < n-1; i++) {
        l[i] = 2 * (h[i] + h[i - 1]) - h[i - 1] * m[i - 1];
        m[i] = h[i] / l[i];
        z[i] = (w[i] - h[i - 1] * z[i - 1]) / l[i];
    }
    l[n-1] = 1;
    z[n-1] = 0;
    std::vector<double> c(n);
    c[n-1] = 0;
    std::vector<Cubic> spline(n-1);
    for (int i = n-1; i-->0;) {
        c[i] = z[i] - m[i] * c[i + 1];
        double b = (a[i + 1] - a[i]) / h[i] - h[i] * (c[i + 1] + 2 * c[i]) / 3;
        double d = (c[i + 1] - c[i]) / (3 * h[i]);
        spline[i] = Cubic(a[i], b, c[i], d);
    }
    return spline;
}
// 매개변수가 균일하지 않는 경우임; 간격 = 입력점 직선거리;
std::vector<CPoint> NatCubicSpline(const std::vector<CPoint>& pts) {
    if (pts.size() < 2) return std::vector<CPoint> (); //null_vector;
    std::vector<double> ax(pts.size());
    std::vector<double> ay(pts.size());
    std::vector<double>  h(pts.size()-1); // interval size;
    for (int i = pts.size(); i-->0;) {
        ax[i] = pts[i].x; 
        ay[i] = pts[i].y;
    };
    // parameter intervals;
    for (int i = pts.size()-1; i-->0;) {
        double dx = ax[i + 1] - ax[i];
        double dy = ay[i + 1] - ay[i];
        h[i] = hypot(dx, dy);
    }

    std::vector<Cubic> splinex = cubicspline1(ax, h);
    std::vector<Cubic> spliney = cubicspline1(ay, h);

    const int steps = 20;
    std::vector<CPoint> curve;
    curve.reserve(1 + splinex.size()*steps);
    curve.push_back(pts.front());
    for (int i = 0; i < pts.size()-1; i++) {
        const double ds = h[i] / steps;
        for (int k = 1; k <= steps; k++) {
            const double tk = k * ds;
            curve.push_back(CPoint(int(splinex[i].eval(tk) + 0.5), 
                                   int(spliney[i].eval(tk) + 0.5)));
        }
    }	
    return curve;
}

저작자표시 비영리 변경금지

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

Catmull-Rom Spline (2) (0)	2024.06.21
Minimum Volume Box (0)	2024.06.16
Polygon Decimation (0)	2024.06.08
Centripetal Catmull-Rom Spline (0)	2024.06.06
Bezier curves (1)	2024.04.19

Natural Cubic Spline

Computational Geometry 2024. 3. 11. 21:55

구간 $[a = t_{0}, . . . t_{n - 1} = b]$ 에서 일정한 간격(꼭 일정한 간격일 필요는 없지만 여기서는 1로 선택함)으로 샘플링된 데이터 ${f_{0}, . . . f_{n - 1}}$ 이 있다. $n - 1$ 개의 각 구간에서 데이터를 보간하는 삼차다항식 함수 모임 $S (t) = {S_{j} (t) | t_{j} \leq t < t_{j + 1}, j = 0, 1, 2. . ., n - 2}$ 을 찾아보자. 전체 구간에서 연속적이고 충분히 부드럽게 연결되기 위해서는 우선 각 node에서 연속이어야 하고, 또한 1차 도함수와 2차 도함수도 연속이 되어야 한다. 물론 각 node에서는 샘플링된 데이터값을 가져야 한다.

$j -$ 번째 구간의 삼차함수를

$S_{j} (t) = a_{j} + b_{j} (t - t_{j}) + c_{j} (t - t_{j})^{2} + d_{j} (t - t_{j})^{3}$

로 쓰면 (a)에서

$S_{j} (t_{j}) = a_{j} = f_{i}, j = 0, 1, . . n - 2$

(b)에서

$a_{j + 1} = a_{j} + b_{j} + c_{j} + d_{j}, j = 0, 1, . . ., n - 2$

(c)에서

$b_{j + 1} = b_{j} + 2 c_{j} + 3 d_{j}, j = 0, 1, . . ., n - 2$

(d)에서

$c_{j + 1} = c_{j} + 3 d_{j}, j = 0, 1, . . ., n - 2$

이므로 (b)와 (c)에서 $d_{j}$ 을 소거하면

$c_{j} = 3 (a_{j + 1} - a_{j}) - b_{j + 1} - 2 b_{j}$

그리고 (a)에서 $d_{j} = - 2 (a_{j + 1} - a_{j}) + b_{j} + b_{j + 1}$ 이므로 $b_{j}$ 에 대해 정리하면 다음과 같은 점화식을 얻을 수 있다.

$b_{j + 1} + 4 b_{j} + b_{j - 1} = 3 (a_{j + 1} - a_{j - 1}) = 3 (f_{j + 1} - f_{j}), j = 1, 2, . . ., n - 2$

물론 $j = 0$ 일 때는 (note, $S_{0}^{″} (t_{0}) = 0 \to c_{0} = 0$ ) $a_{1} = a_{0} + b_{0} + d_{0}, b_{1} = b_{0} + 3 d_{0}$ 이므로

$b_{1} + 2 b_{0} = 3 (f_{1} - f_{0})$

$j = n - 1$ 일 때 계수는 $S_{n - 1} (t_{n - 1}) = f_{n - 1}$ 이고 $S_{n - 1}^{″} (t_{n - 1}) = c_{n - 1} = 0$ 이므로 $2 b_{n - 1} + b_{n - 2} = 3 (f_{n - 1} - f_{n - 2})$

임을 알 수 있다. 따라서 계수를 구하는 과정은 $b_{j} (j = 0, . . ., n - 1)$ 을 구하는 것으로 결정된다. 이를 행렬로 표현하면

$[\begin{matrix} 2 & 1 \\ 1 & 4 & 1 \\ 1 & 4 & 1 \\ \dots \\ 1 & 4 & 1 \\ 1 & 2 \end{matrix}] [\begin{matrix} b_{0} \\ b_{1} \\ ⋮ \\ b_{n - 2} \\ b_{n - 1} \end{matrix}] = [\begin{matrix} 3 (f_{1} - f_{0}) \\ 3 (f_{2} - f_{0}) \\ ⋮ \\ 3 (f_{n - 1} - f_{n - 3}) \\ 3 (f_{n - 1} - f_{n - 2}) \end{matrix}]$

와 같다. band 행렬은 upper triangle로 변환한 후 역치환과정을 거치면 쉽게 해를 구할 수 있다.

평면에서 주어진 점들을 보간하는 곡선은 이들 점을 표현하는 곡선의 매개변수를 일정한 간격으로 나누어서 샘플링된 결과로 취급하면, $x, y$ 성분에 대해서 각각 natural cubic spline를 구하여 얻을 수 있다. 그러나 곡선이 급격히 변할 때는 매개변수를 일정하게 잡는 것보다는 입력점의 사이거리를 기준으로 선택하는 것이 더 유연한 보간곡선을 만들어 준다.

struct Cubic {
    double a,b,c,d;  /* a + b*t + c*t^2 +d*t^3 */
    Cubic() {}
    Cubic(double a_, double b_, double c_, double d_) 
        : a(a_), b(b_), c(c_), d(d_) { }
    /* evaluate;*/
    double eval(double t) {
        return (((d*t) + c)*t + b)*t + a;
    }
};
std::vector<Cubic> calcNaturalCubic(std::vector<double>& x) {
    std::vector<double> gamma(x.size());
    std::vector<double> delta(x.size());
    std::vector<double> D(x.size());
    /* solve the banded equation:
    [2 1       ] [  D[0]]   [3(x[1]   - x[0])  ]
    |1 4 1     | |  D[1]|   |3(x[2]   - x[0])  |
    |  1 4 1   | | .    | = |      .           |
    |    ..... | | .    |   |      .           |
    |     1 4 1| | .    |   |3(x[N-1] - x[N-3])|
    [       1 2] [D[N-1]]   [3(x[N-1] - x[N-2])]
    ** make the banded matrix to an upper triangle;
    ** and then back sustitution. D[i] are the derivatives at the nodes.
    */
    const int n = x.size() - 1;  // note n != x.size()=N;
    // gamma;
    gamma[0] = 0.5;
    for (int i = 1; i < n; i++)
        gamma[i] = 1/(4-gamma[i-1]);
    gamma[n] = 1/(2-gamma[n-1]);
    // delta;
    delta[0] = 3*(x[1]-x[0])*gamma[0];
    for (int i = 1; i < n; i++) 
        delta[i] = (3*(x[i+1]-x[i-1])-delta[i-1])*gamma[i];
    delta[n] = (3*(x[n]-x[n-1])-delta[n-1])*gamma[n];
    // D;
    D[n] = delta[n];
    for (int i = n; i-->0;)
        D[i] = delta[i] - gamma[i]*D[i+1];

    /* compute the coefficients;*/
    std::vector<Cubic> Spline(n);
    for (int i = 0; i < n; i++)
        Spline[i] = Cubic(x[i], D[i], 3*(x[i+1]-x[i])-2*D[i]-D[i+1],
                                      2*(x[i]-x[i+1]) + D[i] + D[i+1])) ;
    return Spline;
}
std::vector<CPoint> NaturalCubicSpline(const std::vector<CPoint>& points) {
    if (points.size() < 2) return std::vector<CPoint> (); // null vector;
    std::vector<double> xp(points.size()), yp(points.size());
    for (int i = points.size(); i-->0;)
        xp[i] = points[i].x, yp[i] = points[i].y;
        
    std::vector<Cubic> splineX = calcNaturalCubic(xp);
    std::vector<Cubic> splineY = calcNaturalCubic(yp);
#define STEPS 12
    std::vector<CPoint> curve;
    curve.reserve(splineX.size() * STEPS + 1);
    curve.push_back(points.front());
    for (int i = 0; i < splineX.size(); i++) {
        for (int j = 1; j <= STEPS; j++) {
            double t = double(j) / STEPS;
            curve.push_back(CPoint(int(splineX[i].eval(t)+0.5), int(splineY[i].eval(t)+0.5)));
        }
    }
    return curve;
}

저작자표시 비영리 변경금지

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

Least Squares Bezier Fit (0)	2024.04.05
Why Cubic Splines? (9)	2024.03.16
Approximate Distance Between Ellipse and Point (0)	2024.03.08
Distance from a Point to an Ellipse (0)	2024.03.06
Data Fitting with B-Spline Curves (0)	2021.04.30

이전 1 다음

내 블로그 - 관리자 홈 전환	`Q` `Q`
새 글 쓰기	`W` `W`

글 수정 (권한 있는 경우)	`E` `E`
댓글 영역으로 이동	`C` `C`

이 페이지의 URL 복사	`S` `S`
맨 위로 이동	`T` `T`
티스토리 홈 이동	`H` `H`
단축키 안내	`Shift` + `/` `⇧` + `/`

Geometry & Recognition

Smoothing Spline

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

Natural Cubic Spline: revisited

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

Natural Cubic Spline

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

카테고리

태그목록

최근에 올라온 글

최근에 달린 댓글

글 보관함

티스토리툴바

개인정보

단축키

내 블로그

블로그 게시글

모든 영역