curvefit.cpp

曲线拟合

// CurveFit.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include "math.h"
#include "stdio.h"
#define ABS(x) (x)>0?(x):-(x)
#define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}

void solve(double **a,double *b,double *x,int n)
{
	int i,j,k,ik;
	double mik,temp;
	for(k=0;k<n;k++)
	{
		mik=-1;
		for(i=k;i<n;i++)
			if(ABS(a[i][k])>mik)
			{
				mik=ABS(a[i][k]);
				ik=i;
			}
		for(j=k;j<n;j++)
			SWAP(a[ik][j],a[k][j]);
		SWAP(b[k],b[ik]);
		b[k]/=a[k][k];
		for(i=n-1;i>=k;i--)
			a[k][i]/=a[k][k];
		for(i=k+1;i<n;i++)
		{
			b[i]-=a[i][k]*b[k];
			for(j=n-1;j>=k;j--)
				a[i][j]-=a[i][k]*a[k][j];
		}
	}
	for(i=n-1;i>=0;i--)
	{
		x[i]=b[i];
		for(j=i+1;j<n;j++)
			x[i]-=a[i][j]*x[j];
	}
}

void linear(double **x,double *y,double *beta,int n,int p)
{
	double **a,*b;
	int i,j,k;
	a=new double*[p];
	for(i=0;i<p;i++)
		a[i]=new double[p];
	for(i=0;i<p;i++)
		for(j=0;j<p;j++)
		{
			a[i][j]=0;
			for(k=0;k<n;k++)
				a[i][j]+=x[k][i]*x[k][j];
		}
	b=new double[p];
	for(i=0;i<p;i++)
	{
		b[i]=0;
		for(j=0;j<n;j++)
			b[i]+=x[j][i]*y[j];
	}
	solve(a,b,beta,p);
	for(i=0;i<p;i++)
		delete a[i];
	delete a,b;
}

void parameter(double *x,double *y,double *u,int n)
{
	int i;
	u[0]=0;
	for(i=0;i<n-1;i++)
		u[i+1]=u[i]+
		sqrt((y[i+1]-y[i])*(y[i+1]-y[i])
		+(x[i+1]-x[i])*(x[i+1]-x[i]));
	double s;
	s=u[n-1]+sqrt((y[0]-y[n-1])*(y[0]-y[n-1])
		+(x[0]-x[n-1])*(x[0]-x[n-1]));
	for(i=0;i<n;i++)
		u[i]/=s;
}

void mintwo(double *p,double *u,double *beta,int n,int r)
{
	int i,j;
	double **x,*y;
	x=new double*[n];
	for(i=0;i<n;i++)
		x[i]=new double[r-1];
	for(i=0;i<n;i++)
		x[i][0]=1;
	for(j=1;j<r-1;j++)
		for(i=0;i<n;i++)
			x[i][j]=pow(u[i],j+2)
			-(j+2)*u[i]*u[i]*0.5+((j+2)*0.5-1)*u[i];
	y=new double[n];
	for(i=0;i<n;i++)
		y[i]=p[i];
	linear(x,y,beta,n,r-1);
	for(i=r;i>2;i--)
		beta[i]=beta[i-2];
	beta[1]=0;
	beta[2]=0;
	for(i=3;i<r+1;i++)
	{
		beta[1]+=(i*0.5-1)*beta[i];
		beta[2]-=0.5*i*beta[i];
	}
	for(i=0;i<n;i++)
		delete x[i];
	delete x,y;
}

double val(double *beta,double u,int r)
{
	int i;
	double rt=0;
	for(i=0;i<=r;i++)
		rt+=beta[i]*pow(u,i);
	return rt;
}

void main()
{
	int i;
	const int n=6;
	double x[n]={0 	,1.5 	,	2,  0 	,	-0.5 	,	-1};
	double y[n]={1 	,0.5 	,	-0.5,-1 	,-.8 	,0 	};
	double u[n];
	parameter(x,y,u,n);
	const int r=4;
	double betax[r+1],betay[r+1];

	mintwo(x,u,betax,n,r);
	mintwo(y,u,betay,n,r);
	FILE *fp=fopen("d:\\spline.txt","w");
	for(double t=0;t<=1;t+=0.001)
		fprintf(fp,"%f %f\n",val(betax,t,r),val(betay,t,r));
	fclose(fp);
}