monotone cubic Hermite interpolation

   이전 글에서 설명한 cubic spline interpolation 은 원 데이터의 monotonicity를 보장해 주지 않는다. 말로 설명하는 것보다 그림으로 보면 쉽게 알 수 있다.

위의 그림에서 보면, 파란 열린 원이 원래의 데이터이다. 앞쪽 열에 있는 그림에서 보면 가장 뒤의 두 점을 보면 두 점은 증가하고 있음에도 cubic spline 으로 연결한 것은 아래로 갑자기 내려 갔다 올라가는 것을 볼 수 있다. 오른쪽 열에 있는 그림은, 앞쪽 두 데이터는 올라가고 있는데, cubic spline으로 연결한 곡선은 올라갔다 내려 오는 것을 볼 수 있다. 만약 이와 같이, 원래 데이터가 올라가는 구간이면 interpolation으로 한 것도 올라가고(최소한 내려오는 곳이 있지는 않고), 원래 데이터가 내려 가는 구간에서는 interpolation으로 한 것도 내려가고(최소한 올라가는 곳은 없고), 이렇게 monotonicity를 유지하게 하려면 어떻게 해야 할까?

위키피디아에 나와 있는 알고리즘(위키피디아의 monotone cubic interpolation) 을 이용하면 다음과 같이 구현할 수 있다.

double h00(double t)

{

    return 2*t*t*t - 3*t*t +1;

}

double h10(double t)

{

    return t*(1-t)*(1-t);

}

double h01(double t)

{

    return t*t*(3-2*t);

}

double h11(double t)

{

    return t*t*(t-1);

}

bool monotonic_cubic_Hermite_spline(int time_limit, const std::vector<double> x_src, const std::vector<double> y_src, std::vector<double>* destX, std::vector<double>* destY)

{

    // 0-based index 사용.

    int n = (int)x_src.size();

    int k = 0;

    double *m = new double[n];

    m[0] = (y_src[1] - y_src[0])/(x_src[1] - x_src[0]);

    m[n-1] = (y_src[n-1] - y_src[n-2])/(x_src[n-1]-x_src[n-2]);

 

    for(k = 1; k<n-1; k++){

        m[k] = (y_src[k] - y_src[k-1])/(2*(x_src[k] - x_src[k-1])) + (y_src[k+1] - y_src[k])/(2*(x_src[k+1]-x_src[k]));

    }

    for(k = 0; k<n-1; k++){

        double delta_k = (y_src[k+1]-y_src[k])/(x_src[k+1]-x_src[k]);

        if(fabs(delta_k) <= eps){

            m[k] = m[k+1] = 0;

        }

        else{

            double ak = m[k]/delta_k;

            double bk = m[k+1]/delta_k;

            if(ak*ak + bk*bk > 9){

                m[k] = 3/(sqrt(ak*ak+bk*bk))*ak*delta_k;

                m[k+1] = 3/(sqrt(ak*ak+bk*bk))*bk*delta_k;

            }

        }

    }

 

    for(k = 0; k<n-1; k++){

        double cur_x = (double)((int)(0.5 + x_src[k]));

        double next_x = (double)((int)(x_src[k+1]));

        double cur_y = y_src[k];

        double next_y = y_src[k+1];

        double h = next_x - cur_x;

        double x = 0;

        double inc = (next_x - cur_x)*0.1;

        for(x = cur_x; x<next_x; x+=inc){

            if(x > time_limit) break;

            double t = (x-cur_x)/h;

            if(destX != NULL){

                destX->push_back(x);

            }

            double y = cur_y*h00(t) + h*m[k]*h10(t) + next_y*h01(t) + h*m[k+1]*h11(t);

            destY->push_back(y);

        }

    }

 

    delete m;

 

    return true;

}

위 알고리즘을 구현하여 확인해 보면 다음과 같다.

위 그림에서 보면 가장 아랫 줄에 있는 것은 monotonicity가 유지되면서 부드럽게 연결된 것을 볼 수 있다.



ps. 원래의 논문은 Fritsch, F. N.; Carlson, R.E. SIAM Journal on Numerical Analysis 1980 Monotone piecewise cubic interpolation, 17 (2) 238‐246, http://dx.doi.org/10.1137/0717021, 으로 보이는데 내가 있는 곳에서 열리지 않아서 읽어 보지 못했다. 그래서 이 알고리즘은 원리를 잘 모르겠다.