C言語関係掲示板

過去ログ

No.971 スプライイン補間

[戻る] [ホームページ]
No.12522

スプライン補間
投稿者---困ってます(2004/02/06 13:30:43)


時間と共に変化する2次元データを1秒毎にX,Yの値が観測される。
(X0,Y0,X1,Y1,X3,Y3,・…)=(0.000000,1.000000,0.295520,0.764842,0.564642,0.169967,0.783327,-0.504846,0.932039,-0.942222,0.997495,-0.936457,0.973848,-0.490261,0.863209,0.186512,0.675463,0.775566,0.427380,0.999859,0.141120,0.753902,-0.157746,0.153374,-0.442520,-0.519289,-0.687766,-0.947722,-0.871576,-0.930426,-0.977530,-0.475537,-0.996165,0.203005,-0.925815,0.786070,-0.772764,0.999435,-0.550686,0.742749,-0.279415,0.136737,0.016814,-0.533584,0.311541,-0.952953,0.578440,-0.924133,0.793668,-0.460679,0.938000,0.219440,0.99854
3,0.796352,0.969890,0.998728,0.854599,0.731386,0.662969,0.120062)
データの数値はどのように使っても良いです。スプライイン補間でときたいです。どのように組み立てればいいんですか?
ヒント下さい。

No.12524

Re:スプライン補間
投稿者---瓜図(2004/02/06 15:14:06)


http://www-cc.ee.tokushima-u.ac.jp/cgi-bin/cboard.cgi?log=&v=676&e=msg&lp=676&st=0

No.12542

Re:スプライン補間
投稿者---困ってます(2004/02/07 00:58:44)


/* スプライン補間 */
#include <math.h>
#include <cgi.h>
#define MM 256
#define NN 256

double Xmin,Xmax,Ymin,Ymax;

Point point(ix, iy)
int ix,iy;
{
Point p;
p.x=ix;
p.y=iy;
return(p);
}

/* メインプログラムの例 */
main()
{
int i,np,mm=100;
double xp[NN],yp[NN],x[MM],y[MM],s[NN];
double alpha[NN],beta[NN],gamma[NN];
printf("xmin,xmax=\n");
scanf("%lf %lf",&Xmin,&Xmax);
printf("ymin,ymax=\n");
scanf("%lf %lf",&Ymin,&Ymax);
printf("np=\n");
scanf("%d",&np);
for ( i=0 ; i<=np-1 ; ++i )
{
printf("x[%d],y[%d]=\n",i,i);
scanf ("%lf %lf",&xp[i],&yp[i]);
}
spline(xp,yp,np,s,alpha,beta,gamma);
hokan(xp,yp,np,alpha,beta,gamma,x,y,mm);
graph(x,y,mm);
}

/* 連続補間 */
hokan(xp,yp,np,alpha,beta,gamma,x,y,mm)
int np,mm;
double xp[NN],yp[NN],x[MM],y[MM];
double alpha[NN],beta[NN],gamma[NN];
{
int i,k;
double xmin,xmax,dx,xi,t;
xmin=xp[0];
xmax=xp[np-1];
dx=(xmax-xmin)/mm;
k=0;
printf(" x y \n");
for ( i=0 ; i<=mm ; ++i )
{
xi=xmin+dx*(double)i;
x[i]=xi;
if (xi>=xp[k+1]) k=k+1;
t=xi-xp[k];
y[i]=((gamma[k]*t+beta[k])*t+alpha[k])*t+yp[k];
/* printf("%15.8lf %15.8lf\n",x[i],y[i]); */
}
}

/* 係数の計算 */
spline(xp,yp,np,s,alpha,beta,gamma)
int np;
double xp[NN],yp[NN],s[NN];
double alpha[NN],beta[NN],gamma[NN];
{
int i,k,n;
double a[NN],b[NN],c[NN],d[NN];
double h,hl,hr,c1,c2;
n=np-2;
/* 連立1次方程式の係数を作る */
for ( k=1 ; k<=n ; ++k )
{
hl=xp[k]-xp[k-1];
hr=xp[k+1]-xp[k];
a[k]=hl/(hl+hr);
b[k]=1-a[k];
c1=(yp[k+1]-yp[k])/hr;
c2=(yp[k]-yp[k-1])/hl;
c[k]=6.0*((c1-c2)/(hl+hr));
d[k]=2.0;
}
/* 前進消去 */
for ( k=1 ; k<n ; ++k )
{
b[k]=b[k]/d[k];
c[k]=c[k]/d[k];
d[k+1]=d[k+1]-a[k+1]*b[k];
c[k+1]=c[k+1]-a[k+1]*c[k];
}
s[n]=c[n]/d[n];
/* 後退代入 */
for ( k=n-1 ; k>=1 ; --k )
s[k]=c[k]-b[k]*s[k+1];
/* 微係数の計算 */
s[0]=s[n+1]=0.0;
for ( k=0 ; k<=n ; ++k )
{
h=xp[k+1]-xp[k];
beta[k]=s[k]/2.0;
gamma[k]=(s[k+1]-s[k])/(6.0*h);
alpha[k]=(yp[k+1]-yp[k])/h
-beta[k]*h-gamma[k]*h*h;
}
}

graph(xx,yy,mm)
int mm;
double xx[NN],yy[NN];
{
cg_init();
cg_mesh();
cg_line(xx,yy,mm);
cg_term();
}

cg_line(xx,yy,nn)
int nn;
double xx[NN],yy[NN];
{
int i,ix,iy;
double cx,cy;
Point pointlist[1280];

cx=1000.0/(Xmax-Xmin);
cy=1000.0/(Ymax-Ymin);
for ( i=0 ; i<=nn ; ++i )
{
ix=(int)((xx[i]-Xmin)*cx)+140;
iy=1000-(int)((yy[i]-Ymin)*cy);
pointlist[i]=point(ix,iy);
}
CGI_Polyline(nn+1,pointlist);
}

cg_mesh()
{
int ix,iy;
Point pointlist[1280];

for ( ix=140 ; ix<=1140 ; ix+=100 )
{
pointlist[0]=point(ix,0);
pointlist[1]=point(ix,1000);
CGI_Polyline(2,pointlist);
}
for ( iy=0 ; iy<=1000 ; iy+=100 )
{
pointlist[0]=point( 140,iy);
pointlist[1]=point(1140,iy);
CGI_Polyline(2,pointlist);
}
}

cg_init()
{
CGI_Initialise();
CGI_VDC_Extent(point(0, 0), point(1279, 1023));
CGI_Clip_Rectangle(point(0, 0), point(1279, 1023));
}

cg_term()
{
CGI_Reset_To_Defaults();
CGI_Terminate();
}

vout(c,v,nn)
char c;
double v[NN];
int nn;
{
int i;
printf("%c ",c);
for ( i=0 ; i<=nn ; ++i )
printf("%9.5lf ",v[i]);
printf("\n");
}

mout(c,a,m,nn)
char c;
double a[MM][NN];
int m,nn;
{
int i,j;
printf("%c\n",c);
for ( j=0 ; j<=nn ; ++j )
{
for ( i=0 ; i<=m ; ++i )
printf("%9.5lf ",a[i][j]);
printf("\n");
}
}


スプライン補間をさがしてきました。
どう書き換えればよいのでしょうか?

No.12549

Re:スプライン補間
投稿者---たか(2004/02/07 17:27:38)


開曲線のスプライン補完です。

/*
 * opened curve spline interpolation
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define MAXN 100
#define MAXP 1000 /* 補完するデータ数 */

struct parameter {
  int n;
  double x[MAXN], y[MAXN];
  double tx[MAXN], ty[MAXN];
  double para[MAXN];
};

void maketable(int n, double x[], double y[], double z[]);
double spline(double t, int n, double x[], double y[], double z[]);
void maketable2(struct parameter *pp);
void spline2(double t, double *px, double *py, struct parameter *pp);

int main(void)
{
  FILE *fi, *fo;
  int i;
  double u, v;
  struct parameter pm;
  
  if ((fi = fopen("spline.txt", "r")) == NULL)
    exit(1);
  
  fscanf(fi, "%d", &pm.n);
  if (pm.n > MAXN) {
    printf("データ数が多すぎます\n");
    exit(1);
  }
  for (i = 0; i < pm.n; i++)
    fscanf(fi, "%lf %lf", &pm.x[i], &pm.y[i]);
  fclose(fi);
  
  maketable2(&pm);
  
  if ((fo = fopen("spline_res.txt", "w")) == NULL)
    exit(1);

  fprintf(fo, "%d\n", MAXP);
  for (i = 0; i <= MAXP; i++) {
    spline2(i / (double)MAXP, &u, &v, &pm);
    fprintf(fo, "%f %f\n", u, v);
  }
  
  fclose(fo);
  return 0;
}

void maketable(int n, double x[], double y[], double z[])
{
  int i;
  double t;
  double h[MAXN], d[MAXN];

  z[0] = 0;  z[n - 1] = 0;
    for (i = 0; i < n - 1; i++) {
    h[i] =  x[i + 1] - x[i];
    d[i + 1] = (y[i + 1] - y[i]) / h[i];
  }
  z[1] = d[2] - d[1] - h[0] * z[0];
  d[1] = 2 * (x[2] - x[0]);
  for (i = 1; i < n - 2; i++) {
    t = h[i] / d[i];
    z[i + 1] = d[i + 2] - d[i + 1] - z[i] * t;
    d[i + 1] = 2 * (x[i + 2] - x[i]) - h[i] * t;
  }
  z[n - 2] -= h[n - 2] * z[n - 1];
  for (i = n - 2; i > 0; i--)
    z[i] = (z[i] - h[i] * z[i + 1]) / d[i];
}

double spline(double t, int n, double x[], double y[], double z[])
{
  int i, j, k;
  double d, h;

  i = 0;  j = n - 1;
  while (i < j) {
    k = (i + j) / 2;
    if (x[k] < t) i = k + 1;  else j = k;
  }
  if (i > 0) i--;
  h = x[i + 1] - x[i];  d = t - x[i];
  return (((z[i + 1] - z[i]) * d / h + z[i] * 3) * d + ((y[i + 1] - y[i]) / h - (z[i] * 2 + z[i + 1]) * h)) * d + y[i];
}

void maketable2(struct parameter *pp)
{
  int i;
  double t1, t2;

  pp->para[0] = 0;
  for (i = 1; i < pp->n; i++) {
    t1 = pp->x[i] - pp->x[i - 1];
    t2 = pp->y[i] - pp->y[i - 1];
    pp->para[i] = pp->para[i - 1] + sqrt(t1 * t1 + t2 * t2);
  }
  for (i = 1; i < pp->n; i++) pp->para[i] /= pp->para[pp->n - 1];
  maketable(pp->n, pp->para, pp->x, pp->tx);
  maketable(pp->n, pp->para, pp->y, pp->ty);
}

void spline2(double t, double *px, double *py, struct parameter *pp)
{
  *px = spline(t, pp->n, pp->para, pp->x, pp->tx);
  *py = spline(t, pp->n, pp->para, pp->y, pp->ty);
}


No.12550

spline.txt(入力データ例)
投稿者---たか(2004/02/07 17:28:29)


30
0.000000 1.000000
0.295520 0.764842
0.564642 0.169967
0.783327 -0.504846
0.932039 -0.942222
0.997495 -0.936457
0.973848 -0.490261
0.863209 0.186512
0.675463 0.775566
0.427380 0.999859
0.141120 0.753902
-0.157746 0.153374
-0.442520 -0.519289
-0.687766 -0.947722
-0.871576 -0.930426
-0.977530 -0.475537
-0.996165 0.203005
-0.925815 0.786070
-0.772764 0.999435
-0.550686 0.742749
-0.279415 0.136737
0.016814 -0.533584
0.311541 -0.952953
0.578440 -0.924133
0.793668 -0.460679
0.938000 0.219440
0.998543 0.796352
0.969890 0.998728
0.854599 0.731386
0.662969,0.120062


No.12552

Re:作ってみたんですがコンパイラが
投稿者---困ってます。(2004/02/07 20:11:34)


見てください。
#include <stdio.h>
#include <math.h>
#include "chinougl.h"
#define NMAX 33
#define PAI 3.1415926
#define OFF (PAI/4)
#define SC 0.5
#define NN 2
void circle(double x0,double y0,double r);
int main(){
int i,j,k;
double x[NMAX];
double y[NMAX][NN];
double sig[NMAX][NN],h[NMAX],del[NMAX][NN],a[NMAX][NN],bet[NMAX][NN];
double b[NMAX][NN],c[NMAX],d[NMAX][NN],
double xx,yy[NN],xx0,yy0[NN],xr,hh,s=0;

FILE*fp;
fp=fopen("spl3.dat","r");

for (i=0;i<NMAX; i++ ){
x[i]=i;
for (i=0;j<NMAX; j++ ){
y[i][j]=sin(x[i]/10.0*(4*j+3.0)+j*PAI/2.0);

fscanf(fp,"%lf,%lf,",&y[i][0],&y[i][1]);
}


for (i=0;i<NMAX-1;++i){
x[i]=i;
h[i]=x[i+1]-x[i];
for (i=0;i<NN;++j){
del[i][j]=(y[i+1][j]-y[i][j])/h[i];
}

//alha & beta


for(j=0;j<NN-1;j++){
al[1][j]=2.0*(h[0]+h[1]);
bet[1][j]=del[1][j]-del[0][j];
}


for(i=2;i<NMAX-1;i++){
al[i][j]=2.0*(h[i-1]+h[i])-h[i-1]*h[i-1]/al[i-1][j];
bet[i][j]=del[i][j]-del[i-1][j]-bet[i-1][j]/al[i-1][j];
}
}
for(j=0;i<NN;j++){
sig[0][j]=0.0;
sig[NMAX-1][j]=0.0;
sig[NMAX-2][j]=bet[NMAX-2][j]/al[NAX-2][j];


for (i=0;i<NMAX-3>0;i--)
for (j=0;j<NN;j++)
sig[i][j]=(bet[i][j]-h[i]*sig[i+1][j])/al[i][j];

for (i=0;i<NMAX-1>0;i++){
for (j=0;j<NN;j++){
b[i][j]=del[i][j]-h[i]*(sig[i+1][j]+2.0*sig[i][j]);
c[i][j]=3.0*sig[i][j];
d[i][j]=(sig[i+1][j]-sig[i][j])/h[i];
}
}




k=0;
for (xx=0.0;xx<=NMAX-1;xx+=PAI/100){
for (i=0;i<NMAX;i++){
if (xx<=x[i+1])
break;
}
xr=xx-x[i];
for (j=0;j<NN;j++)
yy[j]=y[i][j]+(c[i][j]+d[i][j]*xr)*xr)*xr;
if (k!=0){
line(yy[0]*SC,yy[1]*SC,yy0[0]*SC,yy0[1]*SC);
color(0,1,1);
}
for(j=0;j<NN;j++)
yy0[j]=yy[j];
k++;
}
for (i=0;i<NMAX;i++){
circle(y[i][0]*SC,0.02);
fprintf(fp,"%f,%f,",y[i][0],y[i][1]);
}
fclose(fp);
gotogl();
return 0;
}

void circle(double x0,double y0,double r){
int i,n=10;
double x,y,xx,yy,pi=3.1415926;
for(i=0;i<360;i+-n){
x=x0+r*cos((double)i/180.0*pi);
y=y0+r*sin((double)i/180.0*pi);
xx=x0+r*cos((double)(i+n)/180.0*pi);
yy=y0+r*sin((double)(i+n)/180.0*pi);
line(x,y,xx,yy);
colol(1.0,1.0,0.0);
}
}

No.12553

追加
投稿者---困ってます。(2004/02/07 21:00:23)


データ(x0,y0,z0,x1,y1,z1,x2,y2,z2,…)自然スプライン補間して
3次元データを表示する場合、Y=z-0.2y ,X=x-0.4yのように変換してX-Y面表示するために、上のプログラムを書き換えるとすればどうすればいいですか?




No.12554

Re:追加2
投稿者---困ってます。(2004/02/07 21:04:19)


ちなみにデータは32
 x     y    z
0.000000,1.000000,-1.000000,
0.909297,-0.416147,-0.935484,
-0.756802,-0.653644,-0.870968,
-0.279415,0.960170,-0.806452,
0.989358,-0.145500,-0.741935,
-0.544021,-0.839072,-0.677419,
-0.536573,0.843854,-0.612903,
0.990607,0.136737,-0.548387,
-0.287903,-0.957659,-0.483871,
-0.750987,0.660317,-0.419355,
0.912945,0.408082,-0.354839,
-0.008851,-0.999961,-0.290323,
-0.905578,0.424179,-0.225806,
0.762558,0.646919,-0.161290,
0.270906,-0.962606,-0.096774,
-0.988032,0.154251,-0.032258,
0.551427,0.834223,0.032258,
0.529083,-0.848570,0.096774,
-0.991779,-0.127964,0.161290,
0.296369,0.955074,0.225806,
0.745113,-0.666938,0.290323,
-0.916522,-0.399985,0.354839,
0.017702,0.999843,0.419355,
0.901788,-0.432178,0.483871,
-0.768255,-0.640144,0.548387,
-0.262375,0.964966,0.612903,
0.986628,-0.162991,0.677419,
-0.558789,-0.829310,0.741935,
-0.521551,0.853220,0.806452,
0.992873,0.119180,0.870968,
-0.304811,-0.952413,0.935484,
-0.739181,0.673507,1.000000,



No.12556

Re:追加2
投稿者---とおり(2004/02/07 21:19:23)


いい加減【掲示板ご利用上の注意】を読んだらどうですか?
↓とか。

>ソースを添付する際には「HTML変換ツール」で字下げしてください。


No.12569

Re:助けてください
投稿者---困ってます。(2004/02/09 04:10:07)


まだ、出来上がりません。助けてください。



No.12571

出来ましたがあってるかふあんで…
投稿者---困ってます。(2004/02/09 10:50:22)


プロラム見てください。

#include <stdio.h>
#include <math.h>
#include "chinougl.h"
#define nmax 32
#define pai 3.1415926
#define off(pai/4)
#define sc 0.4
#define nn 3

void circle(double x0, double y0, double r);

int main(){
    int i,j,k;
    double x[nmax];
    double y[nmax][nn];
    double sig[nmax][nn],h[nmax],del[nmax][nn],al[nmax][nn],bet[nmax][nn];
    double b[nmax][nn],c[nmax][nn],d[nmax][nn];
    double xx,yy[nn],xx0,xr,hh;
    double xt[nmax],yt[nmax],xt0,yt0;
    //h&delta
    FILE *fp;
    fp=fopen("spl3.dat","r");

    for(i=0;i<nmax;i++){
        x[i]=i;
        for(j=0;j<nn;j++)
            y[i][j]=sin(x[i]/10.0*(4*j+3.0)+j*pai/2.0);

        fscanf(fp,"%lf,%lf,%lf,",&y[i][0],&y[i][1],&y[i][2]);
    }
    //printf("x=%lf,y=%lf,z=%lf",&y[0][0],&y[0][1],&y[0][2]);

    for(i=0;i<nmax-1;i++){
        h[i]=x[i+1]-x[i];
        for(j=0;j<nn;j++)
            del[i][j]=(y[i+1][j]-y[i][j]/h[i]);
    }

    //alpha&beta
    for(j=0;j<nn;j++){
        al[1][j]=2.0*(h[0]+h[1]);
        bet[1][j]=del[1][j]-del[0][j];
    }
    for(i=2;i<nmax-1;i++){
        for(j=0;j<nn;j++){
            al[i][j]=2.0*(h[i-1]+h[i])-h[i-1]*h[i-1]/al[i-1][j];
            bet[i][j]=del[i][j]-del[i-1][j]-bet[i-1][j]*h[i-1]/al[i-1][j];
        }
    }

    for(j=0;j<nn;j++){
        sig[0][j]=0.0;
        sig[nmax-1][j]=0.0;
        sig[nmax-2][j]=bet[nmax-2][j]/al[nmax-2][j];
    }

    for(i=nmax-3;i>0;i--)
        for(j=0;j<nn;j++){
            sig[i][j]=(bet[i][j]-h[i]*sig[i+1][j])/al[i][j];
        }

    for(i=0;i<nmax-1;i++){
        for(j=0;j<nn;j++){
            b[i][j]=del[i][j]-h[i]*(sig[i+1][j]+2.0*sig[i][j]);
            c[i][j]=3.0*sig[i][j];
            d[i][j]=(sig[i+1][j]-sig[i][j])/h[i];
        }

    }

    //gl

    k=0;
    for(xx=0.0;xx<=nmax-1;xx+=pai/100){
        for(i=0;i<nmax;i++){
            if(xx<=x[i+1])
                break;
        }

        xr=xx-x[i];
        for(j=0;j<nn;j++){
            yy[j]=y[i][j]+(b[i][j]+(c[i][j]+d[i][j]*xr)*xr)*xr;

        }
        xt[i]=yy[0]-0.4*yy[1];
        yt[i]=yy[2]-0.2*yy[1];


        if(k!=0){
            line(xt[i]*sc,yt[i]*sc,xt0*sc,yt0*sc);
            color(0,1,1);
        }

        //for(j=0;j<nn;j++)
        //  yy0[j]=yy[j];

        xt0=xt[i];
        yt0=yt[i];

        k++;
    }

    for(i=0;i<nmax;i++){
        circle(xt[i]*sc,yt[i]*sc,0.02);
    }

    fclose(fp);
    gotogl();
    return 0;
}

void circle(double x0,double y0, double r){
    int i,n=10;
    double x,y,xx,yy;

    for(i=0;i<360;i+=n){
        x=x0+r*cos((double)i/180.0*pai);
        y=y0+r*sin((double)i/180.0*pai);
        xx=x0+r*cos((double)(i+n)/180.0*pai);
        yy=y0+r*sin((double)(i+n)/180.0*pai);
        line(x,y,xx,yy);
        color(1.0,1.0,0.0);
    }
}