概要

3次スプラインは、区間ごとにy=ax^3+bx^2+cx+dで表します。本プログラムは、係数a,b,c,dを算出します。コマンドライン上で動作します。
詳しい計算方法は、3次スプライン曲線の作成方法を参照してください。

ソースコード・プログラムのダウンロード spline1.zip

計算結果の見方

inputa dataが補間する前の既知点です。
coefficientのan,bn,cn,dnが各区間の補間係数です。
interpolationが補間された点の一覧です。

計算結果の出力例

input data
n,xn,yn
0,2041.68,1575.59
1,2830.84,2298.99
2,3685.07,1900.71
3,4191.15,1353.67
4,4889.14,1648.74
5,5403.39,2536.84

coefficient
y=an*(x-xn)^3+bn*(x-xn)^2+cn*(x-xn)+dn
n,an,bn,cn,dn
0,-4.54516e-007,0,1.19973,1575.59
1,1.4034e-007,-0.00107606,0.350549,2298.99
2,1.80486e-006,-0.000716409,-1.18063,1900.71
3,-9.6649e-007,0.0020238,-0.518986,1353.67
4,0,0,1.72698,1648.74

interpolation
No,xn,yn
0,2041.68,1575.59
1,2238.97,1808.79
2,2436.26,2021.06
3,2633.55,2191.44
4,2830.84,2298.99
5,3044.4,2326.14
6,3257.95,2263.35
7,3471.51,2118.8
8,3685.07,1900.71
9,3811.59,1743.52
10,3938.11,1585.33
11,4064.63,1448.07
12,4191.15,1353.67
13,4365.65,1319.6
14,4540.15,1377.96
15,4714.64,1497.94
16,4889.14,1648.74
17,5017.7,1870.77
18,5146.27,2092.79
19,5274.83,2314.82
20,5403.39,2536.84

ソースコードの説明

入力座標

座標は、構造体の配列sp2で定義しています。座標数は、マクロ NMAX で定義しています。
補完後の座標は、構造体の配列pに格納されます。

SPLINEクラス(spline1.cpp)

入力座標より係数a,b,c,d及び補完値を計算します。

coefficient

入力座標の設定及び係数a,b,c,dを計算します。
係数の計算には、MATRIXクラスを使います。

interpolation

coefficientによって計算された係数より補完値を計算します。

MATRIX_OWN_COLUMN,MAXRIXクラス(matrix.h)

ヘッダーファイルmatrix.hで定義されSPLINEクラス用に最低限の行列クラスを提供します。

alloc

行列の大きさを定義します。
MATRIX_OWN_COLUMNクラスは、1列固定で、行数が任意となります。
MATRIXクラスは、行数・列数を指定できますが、正方行列のみ正常に動作します。
行列のデーターは、newで確保されます。allocを2回以上実行するとサイズが不足する場合、メモリを再確保します。頻繁に行列サイズが確保する場合は、最初に大きめに確保しておくとよいでしょう。

set

double型の配列を渡すことにより行列に値を設定します。

rswap

行列の行同士を交換します。

Triangular

行列を上三角行列にガウスの前進消去法により変換します。引数によりピポット選択をしないようにも設定できます。

forward_sub

Triangularと同様、行列を上三角行列に変換できますが、このメンバー関数は、他の行列(1列の行列に限る)も同時にピポット選択時に行の交換及びガウスの消去法の対象にできます。引数違いで同名メンバー関数(行列のポインタ)がありますが、これは複数の行列を対象にできます。

backward_sub

後退代入法で連立方程式を解きます。あらかじめ行列は上三角行列に変換しておく必要があります。

det

行列の値を求めます。上三角行列でない場合は、行列のコピーを取ってから値を求めます。

zero

行列を0クリアします。

コンパイル方法

Win32コンソールアプリケーションで空のプロジェクトをつくり、spline1.cppとmatrix.hをプロジェクトに追加します。

ソースコード(spline1.cpp)

// スプライン曲線補間プログラムサンプル
#include <windows.h>
#include <stdio.h>
#include <math.h>
#include <tchar.h>
#include "matrix.h"

struct POINT1{
	double x,y;
};

// 2次座標上のスプライン曲線クラス

class SPLINE{
public:
	double* a;	//	区間ごとの係数
	double* b;	//	区間ごとの係数
	double* c;	//	区間ごとの係数
	double* d;
	int max;
	SPLINE(int n){
		a=new double[n-1];
		b=new double[n-1];
		c=new double[n-1];
		d=new double[n-1];
		max=n;
	}
	~SPLINE(){
		delete [] a;
		delete [] b;
		delete [] c;
		delete [] d;
	}
	void put(int n){
		_tprintf(TEXT("%i,%g,%g,%g,%g\n"),n,a[n],b[n],c[n],d[n]);
	}
	void put(void){
		for(int n=0;n<max-1;n++){
			put(n);
		}
	}
	void coefficient(POINT1* p){	// 座標一覧から係数の算定を行う。
		MATRIX g(max-3,max-3);
		MATRIX_OWN_COLUMN k(max-3);
		MATRIX_OWN_COLUMN s(max-3);
		double* h=new double[max];
		int i;
		for(i=0;i<max-1;i++){
			h[i]=p[i+1].x-p[i].x;
		}
		h[i]=0;
		for(i=1;i<max-2;i++){
			k[i-1]=6*( (p[i+1].y-p[i].y)/h[i] - (p[i].y-p[i-1].y)/h[i-1]); 
		}
		g.zero();
		for(i=0;i<=max-4;i++){
			if(0<i)
				g[i][i-1]=h[i];
			g[i][i]=2*(h[i+1]+h[i]);
			if(i<max-4)
				g[i][i+1]=h[i+1];
		}
		g.forward_sub(k);
		g.backward_sub(k,s);
		double s0,s1;
		for(i=0;i<=max-2;i++){
			if(i==0){
				s0=0;
				s1=s[i];
			}else if(i==max-3){
				s0=s[i-1];
				s1=0;
			}else if(i==max-2){
				s0=0;
				s1=0;
			}else{
				s0=s[i-1];
				s1=s[i];
			}
			a[i]=(s1-s0)/(6.0f*(p[i+1].x-p[i].x));
			b[i]=s0/2.0f;
			c[i]=(p[i+1].y-p[i].y)/(p[i+1].x-p[i].x)-(p[i+1].x-p[i].x)/6.0f*(2.0f*s0+s1);
			d[i]=p[i].y;
		}
		delete h;
	}
	double interpolation(POINT1* p,double x,int n){	//	指定した区間のdx=x-xnに対する補間値
		return a[n]*x*x*x + b[n]*x*x + c[n]*x + d[n];
	}
	void interpolation(POINT1* pb,POINT1* p,int d){		// 区間をd分割して補間値の座標一覧を作成しpに保存する
		int i;
		int u=0;
		for(i=0;i<max-1;i++){
			double dx,x;
			x=0;
			dx=pb[i+1].x-pb[i].x;
			for(int o=0;o<d;o++){
				p[u].x=pb[i].x+x;
				p[u].y=interpolation(pb,x,i);
				++u;
				x+=dx/double(d);
			}
			p[u].x=pb[i+1].x;
			p[u].y=pb[i+1].y;
		}
	}
};

#define NMAX 6	//	入力データー数
#define HMAX 4 //  区間の分割数
//	入力座標

POINT1 sp2[]={{	  2041.68,	1575.59 },
				{ 2830.84,	2298.99 },
				{ 3685.07,	1900.71 },
				{ 4191.15,	1353.67 },
				{ 4889.14,	1648.74 },
				{ 5403.39,	2536.84 }};

POINT1 p[NMAX*HMAX+1]; // 補間された座標が保存される

void main(void){
	int i;
	_tprintf(TEXT("input data\n"));
	_tprintf(TEXT("n,xn,yn\n"));
	POINT1 sp[NM
	for(i=0 ; i<NMAX ; i++){
		sp[i].y=sp2[i].y;
		sp[i].x=sp2[i].x;
		_tprintf(TEXT("%i,%g,%g\n"),i,sp[i].x,sp[i].y);
	}
	SPLINE s(NMAX);
	s.coefficient(sp);
	_tprintf(TEXT("\ncoefficient\ny=an*(x-xn)^3+bn*(x-xn)^2+cn*(x-xn)+dn\n"));
	_tprintf(TEXT("n,an,bn,cn,dn\n"));
	s.put();
	s.interpolation(sp,p,HMAX);
	_tprintf(TEXT("\ninterpolation\n"));
	_tprintf(TEXT("No,xn,yn\n"));
	for(i=0 ; i<(NMAX-1)*HMAX+1 ; i++){
		_tprintf(TEXT("%i,%g,%g\n"),i,p[i].x,p[i].y);
	}
}
 

ソースコード(matrix.h)

#ifndef MATRIX001_H


//	スプライン用の最小限の行列クラス


#define MATRIX001_H 1

// 行列の状態

enum MATRIX_TYPE{ 	undecided=0,upper_triangular,lowwer__triangular };

// 1列の行列クラス

class MATRIX_OWN_COLUMN{
	double* d;	//	行列データー
	double k;	//	行列の符号および乗数を示す 例えばスカラー倍した場合はこの係数のみ乗ずる。
	int rmax;	//	行列 行数
	int max;	//	確保されている配列のサイズ
	MATRIX_TYPE flag;	//	行列タイプ
public:
	MATRIX_OWN_COLUMN(){
		d=0;
		rmax=max=0;
	}
	MATRIX_OWN_COLUMN(int ar){
		d=0;
		rmax=max=0;
		alloc(ar);
	}
	~MATRIX_OWN_COLUMN(){
		delete [] d;
	}
	void alloc(int ar){
		if(max<ar){	//	配列サイズが不足している場合は再確保する
			max=rmax=ar;
			if(d)
				delete [] d;
			d=new double[rmax];
		}else{
			rmax=ar;
		}
		flag=undecided;
		k=1;
	}
	void rswap(int a,int b){	//	行の交換
		double t=d[a];
		d[a]=d[b];
		d[b]=t;
		k=-k;
	}
	void set(double* ad){	//	行列への定数の代入
		for(int n=0;n<rmax;n++)
			d[n]=ad[n];
	}
	void put(void){	//	行列を標準出力へ表示
		for(int rn=0;rn<rmax;rn++){
			_tprintf(TEXT("%i %g\n"),rn,d[rn]);
		}
	}
	double& operator[](int n){
		return d[n];
	}
	friend class MATRIX;
};

//	正方行列クラス


class MATRIX{
	double* d;	//	行列データー
	int cmax,rmax;	//	行列 列数 行数
	double k;	//	行列の符号および乗数を示す 例えばスカラー倍した場合はこの係数のみ乗ずる。det
	int max;	//	確保されている配列のサイズ
	MATRIX_TYPE flag;	//	行列タイプ
public:
	MATRIX(){
		d=0;
		cmax=rmax=0;
		max=0;
		k=1;
	}
	MATRIX(int ar,int ac){
		d=0;
		cmax=rmax=0;
		max=0;
		k=1;
		alloc(ar,ac);
	}
	~MATRIX(){
		delete [] d;
	}
	void alloc(int ar,int ac){
		cmax=ac;
		rmax=ar;
		if(max<ar*ac){	//	配列サイズが不足する場合は配列を再確保する
			if(d){
				delete [] d;
			}
			max=cmax*rmax;
			d=new double[max];
		}
		flag=undecided;
		k=1;
	}
	double* operator[](int n){
		return d+n*cmax;
	}
	void rswap(int a,int b){	//	行の交換
		for(int i=0;i<cmax;i++){
			double t=(*this)[a][i];
			(*this)[a][i]=(*this)[b][i];
			(*this)[b][i]=t;
		}
		k=-k;	//	符号を反転する
	}
	void set(double* ad){	//	行列への定数の代入
		for(int n=0;n<cmax*rmax;n++)
			d[n]=ad[n];
	}
	void put(void){	//	行列を標準出力へ表示
		for(int rn=0;rn<rmax;rn++){
			for(int cn=0;cn<cmax;cn++){
				_tprintf(TEXT("%g "),(*this)[rn][cn]);
			}
			_puttc(_T('\n'),stdout);
		}
	}
	void Triangular(bool pivoting=true){ // 上三角行列に変換(ガウスの前進消去法) 引数でピポット選択の有無を選択
		int i,j,k;
		double uk1;
		for(k=0;k<rmax-1;k++){
			double m=0;
			double u;
			int n=k;
			if(pivoting){	//	ピポット選択を行う
				for(i=k;i<rmax;i++){	//	絶対値が最大の行の検索
					u=fabs( (*this)[i][k] );
					if(m!=0 && m<u){
						m=u;
						n=i;
					}
				}
				if(n!=k)
					rswap(n,k);	//	絶対値が最大の行と入れ替え
			}
			for(i=k+1;i<rmax;i++){
				uk1 = (*this)[i][k] / (*this)[k][k];	//	  ui1=ai1/a11
				(*this)[i][k]=0;		//	上三角行列にすると参照されないので一般的には不要
				for(j=k+1;j<cmax;j++){
					 (*this)[i][j] = (*this)[i][j] - uk1 * (*this)[k][j];	 //aji = aji - ui1 * a1j
				}
			}
		}
		flag=upper_triangular;	//	行列タイプを設定
	}
	void forward_sub(MATRIX_OWN_COLUMN& b,bool pivoting=true){ // 上三角行列に変換(ガウスの前進消去法)  引数でピポット選択の有無を選択
		int i,j,k;
		double uk1;
		for(k=0;k<rmax-1;k++){
			double m=0;
			double u=0;
			int n=k;
			if(pivoting){
				for(i=k;i<rmax;i++){	//	絶対値が最大の行の検索
					u=fabs( (*this)[i][k] );
					if(m<u){
						m=u;
						n=i;
					}
				}
				if(m!=0 && n!=k){
					rswap(n,k);	//	絶対値が最大の行と入れ替えただしすべてが0の時は入れ替えない
					b.rswap(n,k);
				}
			}
			for(i=k+1;i<rmax;i++){
				uk1 = (*this)[i][k] / (*this)[k][k];	//	  ui1=ai1/a11
 				b[i] = b[i] - uk1 * b[k];
 				(*this)[i][k]=0;		//	上三角行列にすると参照されないので一般的には不要

				for(j=k+1;j<cmax;j++){
					 (*this)[i][j] = (*this)[i][j] - uk1 * (*this)[k][j];	 //aji = aji - ui1 * a1j
				}
			}
		}
		flag=upper_triangular;	//	行列タイプを設定
	}
	void forward_sub(MATRIX_OWN_COLUMN* b,int nmax, bool pivoting=true){ // 上三角行列に変換(ガウスの前進消去法) ピポット選択有り
		int i,j,k;
		double uk1;
		for(k=0;k<rmax-1;k++){
			double m=0;
			double u=0;
			int n=k;
			if(pivoting){
				for(i=k;i<rmax;i++){	//	絶対値が最大の行の検索
					u=fabs( (*this)[i][k] );
					if(m<u){
						m=u;
						n=i;
					}
				}
				if(m!=0 && n!=k){
					rswap(n,k);	//	絶対値が最大の行と入れ替えただしすべてが0の時は入れ替えない
					for(int o=0;o<nmax;o++){
						b[o].rswap(n,k);
					}
				}
			}

			for(i=k+1;i<rmax;i++){
				uk1 = (*this)[i][k] / (*this)[k][k];	//	  ui1=ai1/a11
				for(int o=0;o<nmax;o++){
					MATRIX_OWN_COLUMN& t=b[o];
					t[i] = t[i] - uk1 * t[k];
				}
 				(*this)[i][k]=0;		//	上三角行列にすると参照されないので一般的には不要

				for(j=k+1;j<cmax;j++){
					 (*this)[i][j] = (*this)[i][j] - uk1 * (*this)[k][j];	 //aji = aji - ui1 * a1j
				}
			}
		}
		flag=upper_triangular;	//	行列タイプを設定
	}
	void backward_sub(MATRIX_OWN_COLUMN& b,MATRIX_OWN_COLUMN& x){	//	後退代入法 あらかじめ上三角行列に変換しておく必要がある。
		int i,j;
		for(i=rmax-1;0<=i;i--){
			double c=b[i]*b.k;
			for(j=i+1;j<rmax;j++){
				c -= (*this)[i][j] * x[j];
			}
			x[i] = c / (*this)[i][i];
		}
	}
	void backward_sub(MATRIX_OWN_COLUMN* b,MATRIX_OWN_COLUMN* x,int nmax){	//	後退代入法 あらかじめ上三角行列に変換しておく必要がある。
		int i,j;
		for(i=rmax-1;0<=i;i--){
			for(int n=0;n<nmax;n++){
				double c=b[n][i]*b[n].k;
				for(j=i+1;j<rmax;j++){
					c -= (*this)[i][j] * x[n][j];
				}
				x[n][i] = c / (*this)[i][i];
			}
		}
	}
	double det(void){	//	行列の値を算出
		double a=k;
		if(flag==upper_triangular){	//	上三角行列
			for(int i=0;i<cmax;i++)
				a *= (*this)[i][i];
			return a;
		}else{	// 上三角行列に変換してから値を求める
			MATRIX x(rmax,cmax);
			x=(*this);
			x.Triangular();
			return x.det();
		}
	}
	void zero(void){	//	行列の初期化
		ZeroMemory((void*)d,sizeof(d[0])*rmax*cmax);
	}
};

#endif