山本ワールド
Windowsプログラミング
アルゴリズム Vitual C++ 2008/2013によるWin32/Win64 APIレベルのプログラム 基礎 Vitual C++ 2008/2013によるAPIレベルのプログラム(32/64bit) Wix3でインストーラーを作る Visual C++ 2008 Standard Editonによるフォームアプリケーションのプログラム(32/64bit) Vitual C++ 2008 Standard EditonによるAPIレベルのプログラム(32/64bit) Windows 7対応 Visual C++ 2008 ExpressによるAPIレベルのプログラム Visual C++ 2005 ExpressによるAPIレベルのプログラム Visual C++ Versiosn 5 BORLAND C++ Windowsプログラム全般 Excel VBA その他3次スプライン曲線の係数算定プログラム(32/64bit)
概要
3次スプラインは、区間ごとにy=ax3+bx2+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と同様、行列を上三角行列に変換できますが、このメンバー関数は、他の行列(１列の行列に限る)も同時にピポット選択時に行の交換及びガウスの消去法の対象にできます。引数違いで同名メンバー関数(行列のポインタ)がありますが、これは複数の行列を対象にできます。
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
Copyright (C) 2012 山本ワールド All Rights Reserved.