3次スプライン補間

gistfile1.js

// 3次スプライン補間
// 参考：高橋大輔，数値計算，岩波書店，1996.
// 入力 N+1個の点列 points = [[x0,y0],[x1,x2],...,[xn,yn]]
// 出力 各データ点間の区間における多項式の係数を含む2次元配列
// x：x座標値で昇順ソートされた点列のx座標
// y：x座標値で昇順ソートされた点列のy座標
// a：区間[x_i+1 - x_i]における3次式の係数 配列サイズN-1
// b：区間[x_i+1 - x_i]における3次式の係数 配列サイズN-1
// c：区間[x_i+1 - x_i]における3次式の係数 配列サイズN-1
// d：区間[x_i+1 - x_i]における3次式の係数 配列サイズN-1
// ※三次式 y = a*x^3 + b*x^2 + c*x +d
function spline3Interpolation(points) {

	// 入力点数が3未満ではスプライン曲線を定義できないので終了
	if(points.length<3) {
		return null;
	}

	// x座標の値で昇順にソート
	var tmp=numeric.clone(points);
	var sortedPoints=tmp.sort(
		function (a, b) {
			if(a[0]<b[0]) {
				return -1;
			} else {
				return 1;
			}
		}
	);

	// データ点
	var x=new Array(sortedPoints.length);
	var y=new Array(sortedPoints.length);

	// 曲線が定義される区間の数
	var N=points.length-1;

	// 連立一次方程式を構成する行列
	var h=new Array(N);
	var H=numeric.rep([N-1, N-1], 0);
	var u=new Array(N+1);
	var v=new Array(N-1);

	// 各区間の係数 y = a*x^3 + b*x^2 + c*x +d
	var a=new Array(N);
	var b=new Array(N);
	var c=new Array(N);
	var d=new Array(N);

	for(var i=0; i<sortedPoints.length; ++i) {
		x[i]=sortedPoints[i][0];
		y[i]=sortedPoints[i][1];
	}

	for(var i=0; i<N; ++i) {
		h[i] = x[i+1] - x[i];
	}

	H[0][0]=2*(h[0]+h[1]);
	H[0][1]=h[1];
	for(var i=1; i<N-2; ++i) {
		H[i][i-1]=h[i];
		H[i][i]=2*(h[i]+h[i+1]);
		H[i][i+1]=h[i+1];
	}
	H[N-2][N-3]=h[N-2];
	H[N-2][N-2]=2*(h[N-2]+h[N-1]);

	for(var i=1; i<N; ++i) {
		v[i-1]=6*((y[i+1]-y[i])/h[i]-(y[i]-y[i-1])/h[i-1]);
	}

	var uBlock=numeric.solve(H, v);	// from u[1] to u[N-1]
	u[0]=0;
	u[N]=0;
	for(var i=0; i<N-1; ++i) {
		u[i+1] = uBlock[i];
	}

	for(var i=0; i<N; ++i) {
		a[i]=(u[i+1]-u[i])/(6*(x[i+1]-x[i]));
	}

	for(var i=0; i<N; ++i) {
		b[i]=0.5*u[i];
	}

	for(var i=0; i<N; ++i) {
		c[i]=(y[i+1]-y[i])/(x[i+1]-x[i])-(x[i+1]-x[i])*(2*u[i]+u[i+1])/6;
	}

	for(var i=0; i<N; ++i) {
		d[i]=y[i];
	}

	return [x,y,a,b,c,d];
}