HOME/Articles/

スプライン補間をしたい

Article Outline

スプライン補間をしたい

とにかく補間をしなければいけないことがあります.特に方式にこだわりもない,そんな時にとりあえず選ばれるオーソドックスな方式がスプライン補間だと思います.スプライン補間で描かれた曲線は,必ず既知のデータ上を通ります.
今回もとりあえずコピペして簡単に使えるC++のコードを用意しました.もちろん動作保証はないです.数式とか細かいことは他のサイトを当たってください.

基本的な使い方については,まず既知の等間隔のy軸方向のデータをvector型で用意します.次に,そのデータを引数にしてInitParameter関数を呼びます. 計算結果は,0~データ数の間を引数としてCalc関数を呼ぶことで取得することができます.
以下のサンプルプログラムは0.01間隔でxが0~6の間の結果を取得することで,元データに対して100倍の細かさで補間データを取得する例になります.

CSpline calc;
vector<double> sy;//spline用
float skip = 0.01f;//skip幅
sy.clear();

//等間隔のyデータをセット
sy.push_back(1);//データを入力
sy.push_back(6);//データを入力
sy.push_back(4);//データを入力
sy.push_back(7);//データを入力
sy.push_back(0);//データを入力
sy.push_back(3);//データを入力

calc.InitParameter(sy); //データをセット
for (float i = 0; i < sy.size(); i += skip)
{
    cout << calc.Calc(i) << endl;
}

使い方は以上で,以下が実際のコードになります.間違いもあるかと思いますが適当に使ってください.

header file

#include <vector>

using namespace std;
class CSpline
{
public:
    CSpline();
    ~CSpline();

    //データセット関数
    void InitParameter(const vector<double> &y);

    //演算関数
    double Calc(float t);

private:
    vector<double> a_;
    vector<double> b_;
    vector<double> c_;
    vector<double> d_;
    vector<double> w_;
};

cpp file

#include "CSpline.h"

CSpline::CSpline()
{
}

CSpline::~CSpline()
{
}

//CSplineデータセット関数
void CSpline::InitParameter(const vector<double> &y) {
    int ndata = (int)y.size() - 1;

    for (int i = 0; i <= ndata; i++)
    {
        a_.push_back(y[i]);
    }

    for (int i = 0; i <= ndata; i++)
    {
        if (i == 0)
        {
            c_.push_back(0.0);
        }
        else if (i == ndata)
        {
            c_.push_back(0.0);
        }
        else
        {
            c_.push_back(3.0*(a_[i - 1] - 2.0*a_[i] + a_[i + 1]));
        }
    }

    for (int i = 0; i < ndata; i++)
    {
        if (i == 0)
        {
            w_.push_back(0.0);
        }
        else
        {
            double tmp = 4.0 - w_[i - 1];
            c_[i] = (c_[i] - c_[i - 1]) / tmp;
            w_.push_back(1.0 / tmp);
        }
    }

    for (int i = (ndata - 1); i > 0; i--)
    {
        c_[i] = c_[i] - c_[i + 1] * w_[i];
    }

    for (int i = 0; i <= ndata; i++) {
        if (i == ndata)
        {
            d_.push_back(0.0);
            b_.push_back(0.0);
        }
        else
        {
            d_.push_back((c_[i + 1] - c_[i]) / 3.0);
            b_.push_back(a_[i + 1] - a_[i] - c_[i] - d_[i]);
        }
    }
}

//CSpline演算関数
double CSpline::Calc(float t)
{
    int j = int(floor(t));
    if (j < 0) {
        j = 0;
    }
    else if (j >= (int)a_.size()) {
        j = ((int)a_.size() - 1);
    }

    double dt = t - j;
    double result = a_[j] + (b_[j] + (c_[j] + d_[j] * dt)*dt)*dt;
    return result;
}