3次スプライン曲線の生成 #数学

だいぶ前に書きかけて限定公開になっていたのを公開してみる。（2014/07/01に書いたやつだった）

前に勉強しようと思って、あまり詳しく調べていなかったスプライン曲線について調べたのでそのメモと備忘録を残しておきます。
※ちなみにスプライン曲線はいくつかの方法があり、色々と内容が異なるようなので調べる際はその辺りに注意。

なお、今回の記事はこちらの記事を参考にさせてもらいました。というか、こちらの記事を自分なりに理解した範囲でメモを加えたものです。

スプライン曲線は、簡単に言えば与えられた任意の点（$n+1$個）の各点を滑らかに結ぶ曲線の一種です。
さらに3次スプライン曲線とは、その各点の間の補完された曲線を3次関数で表したものになります。

各点の間は「セグメント」と呼ばれます。
そのセグメント（S）の間を補完する曲線を表す3次関数は以下の式になります。

S_j(x) = a_j + b_j(x - x_j) + c_j(x - x_j)^2 + d_j(x - x_j)^3  \

ただし、x_j≦x≦x_{j+1}

前提条件

なお、5つの前提条件があります。

$S_j(x_j) = y_j$
$S_j(x_{j+1}) = S_{j+1}(x_{j+1}) = y_{j+1}$
$S_j'(x_{j+1}) = S_{j+1}'(x_{j+1})$
$S_j''(x_{j+1}) = S_{j+1}''(x_{j+1})$
$S_0''(0) = S_{n-1}''(x_n) = 0$

$S_j(x)$が点$(x_j, y_j)$と点$(x_{j+1}, y_{j+1})$を通り、各関数が滑らかに接続するように、上記関数の$a_j, b_j, c_j, d_j$を決定します。

※3, 4は隣合うセグメント同士が滑らかにつながるようにするための条件です。すなわち、ふたつの3次関数を微分したものが同一となる（$S_j'(x_{j+1}) = S_{j+1}'(x_{j+1})$、$S_j''(x_{j+1}) = S_{j+1}''(x_{j+1})$）ように条件を設定するものです。

上記前提を元に、係数$a_j, b_j, c_j, d_j$を計算すると以下のようになります。
（$h_j = x_{j+1} - x_j$）

$a_j = y_j$
$b_j = \frac{a_{j+1} - a_j}{h_j} - \frac{h_j(c_{j+1} + 2c_j)}{3}$
$d_j = \frac{c_{j+1} - c_j}{3h_j}$

参考記事のコードをJSで書き直したものが以下になります。

spline

function Spline() {
    this.a = [];
    this.b = [];
    this.c = [];
    this.d = [];

    this.num = 0;
}
Spline.prototype = {
    constructor: Spline,
    num: 0,
    a: null,
    b: null,
    c: null,
    d: null,
    init: function (sp) {
        var tmp;
        var w = [];
        var i;

        this.num = sp.length - 1;
        var num = this.num;

        var a = this.a;
        var b = this.b;
        var c = this.c;
        var d = this.d;

        // 3次多項式の0次係数(a)を設定
        for(i = 0; i <= num; i++) {
            a[i] = sp[i];
        }

        // 3次多項式の2次係数(c)を計算
        c[0] = c[num] = 0.0;
        for(i = 1; i < num; i++) {
            c[i] = 3.0 * (a[i - 1] - 2.0 * a[i] + a[i + 1]);
        }

        // 左下を消す
        w[0] = 0.0;
        for(i = 1; i < num; i++) {
            tmp = 4.0 - w[i - 1];
            c[i] = (c[i] - c[i - 1]) / tmp;
            w[i] = 1.0 / tmp;
        }

        // 右上を消す
        for(i = num - 1; i > 0; i--) {
            c[i] = c[i] - c[i + 1] * w[i];
        }

        // 3次多項式の1次係数(b)と3次係数(b)を計算
        b[num] = d[num] = 0.0;
        for (i = 0; i < num; i++) {
            d[i] = (c[i + 1] - c[i]) / 3.0;
            b[i] = a[i + 1] - a[i] - c[i] - d[i];
        }
    },

    // 媒介変数(0～num - 1の実数）に対する値を計算
    culc: function (t) {
        var j = floor(t);
        var dt;
        var num = this.num;

        var a = this.a;
        var b = this.b;
        var c = this.c;
        var d = this.d;

        if (j < 0) {
            j = 0;
        }
        else if (j >= num) {
            j = num - 1; // 丸め誤差を考慮
        }

        dt = t - j;

        return a[j] + (b[j] + (c[j] + d[j] * dt) * dt ) * dt;
    }
};

実際にこれを用いて簡単なデモを作成しました。

▼ jsdo.itで公開した、スプライン曲線を利用した彗星です

 http://jsdo.it/edo_m18/muEK

使い方

how-to-use

for (var i = 0; i < 20; i++) {
    arr.push(range(-200, 200));
}

var sp = new Spline();
sp.init(arr);

var t = /* any value as parameter */ ;
var val = sp.culc(t);

サンプルなのでランダムにポイントを生成しています。
そして Spline#init でスプラインのセットアップを行い、 culc メソッドによって該当の位置の値を取得する、という流れになります。

ここでの t はいわゆる「媒介変数」です。
なので、この値を様々に変えると曲線上の t の位置にある値を取得することができる、というわけです。
※ なお、実装にもコメントで書いていますが、この媒介変数 t の取りうる値は 0 〜 (arr.length - 1) となります。