spine.cpp

样条类的C++,该类包括三次样条的插值,拟合和B样条的插值,拟和等功能 。自己编写的简单适用的函数

// Spine.cpp: implementation of the CSpine class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "Spine.h"

#include"memory.h"
#include"math.h"
#include"Matrix.h"

#define MAX_PAR_NUM 150//三次样条或B样条拟合中最大参数上限，避免计算时间过大

#define ST_CLAMPED  1 //表明三次样条插值采用的是指定一阶导数条件
#define ST_NATURAL  2 //表明三次样条插值采用的是自然参数条件 
#define ST_EXTRA    3 //表明三次样条插值采用的是外插条件 
#define ST_CONSTANT 4 //表明三次样条插值采用的二阶倒数是常数
#define ST_TWO      5 //表明三次样条插值采用的是指定二阶导数条件

#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////
//默认构造函数
//参数：无
//返回：无
//////////////////////////////////////////////////////////////////////
CSpine::CSpine()
{
	m_nNum=1;
	m_pTime=NULL;
	m_pData=NULL;
	m_p3SpineFit=NULL;
	m_p3SpineInsert=NULL;
	m_pBSpineFit=NULL;
	m_pBSpineInsert=NULL;
	BOOL bInit=Init(m_nNum);
	ASSERT(bInit);
}

//////////////////////////////////////////////////////////////////////
//构造函数,初附值
//参数：1 nNum------初始数据个数
//      2 Time------初始数据坐标点的值(注意是nNum+1维)
//      3 Data------初始数据坐标点的函数值(注意是nNum+1维)
//返回：无
//////////////////////////////////////////////////////////////////////
CSpine::CSpine(int nNum,double Time[],double Data[])
{
	m_nNum=nNum;

	m_pTime=NULL;
	m_pData=NULL;
	m_p3SpineFit=NULL;
	m_p3SpineInsert=NULL;
	m_pBSpineFit=NULL;
	m_pBSpineInsert=NULL;
	BOOL bInit=Init(m_nNum);
	ASSERT(bInit);
    SetTime(Time);
	SetData(Data);

}
///////////////////////////////////////////////////////////////////
// 拷贝构造函数
//////////////////////////////////////////////////////////////////
CSpine::CSpine(const CSpine& other)
{
	m_nNum=other.m_nNum;

	m_pTime=NULL;
	m_pData=NULL;
	m_p3SpineFit=NULL;
	m_p3SpineInsert=NULL;
	m_pBSpineFit=NULL;
	m_pBSpineInsert=NULL;	

	BOOL bInit=Init(m_nNum);
	ASSERT(bInit);
	memcpy(m_pTime,other.GetTime(),sizeof(double)*m_nNum);
	memcpy(m_pData,other.GetData(),sizeof(double)*m_nNum);
}
/////////////////////////////////////////////////////////////////////
//析构函数,释放内存
//参数:无
//返回:无
//////////////////////////////////////////////////////////////////////
CSpine::~CSpine()
{
	if(m_pTime)
		delete[] m_pTime;
	if(m_pData)
		delete[] m_pData;
	if(m_p3SpineFit)
		delete[] m_p3SpineFit;
	if(m_p3SpineInsert)
		delete[] m_p3SpineInsert;
	if(m_pBSpineFit)
		delete[] m_pBSpineFit;
	if(m_pBSpineInsert)
		delete[] m_pBSpineInsert;
}

////////////////////////////////////////////////////////////////
//初始函数,分配内存
//参数:1 nNum------初始数据个数
//返回:成功分配内存返回TRUE,否则FALSE
//////////////////////////////////////////////////////////////////
BOOL CSpine::Init(int nNum)
{
	
	m_nNum=nNum;
	//////////////////////////////////////////////////
	//释放内存
	if(m_pTime)
		delete[] m_pTime;
	if(m_pData)
		delete[] m_pData;
	if(m_p3SpineFit)
		delete[] m_p3SpineFit;
	if(m_p3SpineInsert)
		delete[] m_p3SpineInsert;
	if(m_pBSpineFit)
		delete[] m_pBSpineFit;
	if(m_pBSpineInsert)
		delete[] m_pBSpineInsert;
	//////////////////////////////////////////////////
	//分配内存
	m_pTime=new double[m_nNum];
	if(!m_pTime)
		return FALSE;
	m_pData=new double[m_nNum];
	if(!m_pData)
		return FALSE;
	

	/////////////////////////////////////////////////////////
	//初赋值
	memset(m_pTime,0,sizeof(double)*(m_nNum));
	memset(m_pData,0,sizeof(double)*(m_nNum));


	//////////////////////////////////////////////////////
	return TRUE;
}

////////////////////////////////////////////////////////////////
//设置初始坐标点的函数。注意时间是nNum+1维
//参数：1 Time------初始坐标点的数据集
//返回： 无
/////////////////////////////////////////////////////////////////
void CSpine::SetTime(double Time[])
{

	memset(m_pTime,0,sizeof(double)*(m_nNum));
	memcpy(m_pTime,Time,sizeof(double)*(m_nNum));
}

////////////////////////////////////////////////////////////////
//设置初始坐标点函数值的函数。注意是nNum+1维
//参数：1 Data------初始坐标点的函数值数据集
//返回： 无
/////////////////////////////////////////////////////////////////
void CSpine::SetData(double Data[])
{
	memset(m_pData,0,sizeof(double)*(m_nNum));
	memcpy(m_pData,Data,sizeof(double)*(m_nNum));
}

///////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////
double* CSpine::GetTime() const
{
	return m_pTime;
}

///////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////
double* CSpine::GetData() const
{
	return m_pData;
}


//////////////////////////////////////////////////////////////////
//B样条拟合函数
//参数：1 nInterval----用户要分的区间个数
//      2 nType--------用户指定的B样条阶数
//返回：残查
//////////////////////////////////////////////////////////////////
double CSpine::BSpineFit(int nParameter,int nType)
{
	if(nParameter<4)
	{
		AfxMessageBox("你输入的参数过小，请确认参数个数大于3！");
	}
	if(nParameter>MAX_PAR_NUM)
	{	AfxMessageBox("你输入的参数过大，将有可能造成数值上的困难！");
	}
	if(m_pBSpineFit)
		delete[] m_pBSpineFit;
	m_pBSpineFit=new double[nParameter];
	ASSERT(m_pBSpineFit);
	double *T=new double[nParameter];
	ASSERT(T);
	double a=m_pTime[0];
	double b=m_pTime[m_nNum-1];
	double h=(b-a)/(nParameter-3);
	for(int i=0;i<nParameter;i++)
		T[i]=a-h+i*h;
	/////////////////////////////////////////////////////////
	CMatrix mtxX(m_nNum,nParameter);
	CMatrix mtxY(m_nNum,1);
	for(i=0;i<m_nNum;i++)
		for(int j=0;j<nParameter;j++)
			mtxX.SetElement(i,j,ggBFuction((m_pTime[i]-T[j])/h,nType));
	mtxY.SetData(m_pData);
	CMatrix mtxL(nParameter);
	mtxL=mtxX.Transpose()*mtxX;

	mtxL.InvertGaussJordan();

	CMatrix mtxResult;
	mtxResult=mtxL*(mtxX.Transpose());
	mtxResult=mtxResult*mtxY;
	memcpy(m_pBSpineFit,mtxResult.GetData(),sizeof(double)*(nParameter));
	delete[] T;
	double er=0;
	for(i=0;i<m_nNum;i++)
		er+=(m_pData[i]-GetBSpineFitValue(m_pTime[i],nParameter,nType))*
		(m_pData[i]-GetBSpineFitValue(m_pTime[i],nParameter,nType));
	er=er/m_nNum;
	return sqrt(er);


}

///////////////////////////////////////////////////////////////////////////////
//求B样条拟合后给返回参数指针，注意：使用前务必先调用GetBSpineFit函数，切记
//参数：无
//返回：B样条拟合后给返回参数指针，个数为GetBSpineFit函数里面的nParameter参数 
//////////////////////////////////////////////////////////////////////////////
double* CSpine::GetBSpineFitPointer() const
{
	return m_pBSpineFit;
}

////////////////////////////////////////////////////////////////////////////////
//给定一个点，求B样条拟合后给定一个点求值，注意，求值前必须调用GetBSpineFit函数
//参数：1 x------------指定一个坐标
//      2 nInterval----用户要分的区间个数
//      3 nType--------用户指定的B样条阶数
//返回：该点拟合后的函数值
/////////////////////////////////////////////////////////////////////////////
double CSpine::GetBSpineFitValue(double x,int nParameter,int nType) 
{
	
	double *T=new double[nParameter];
	ASSERT(T);
	double a=m_pTime[0];
	double b=m_pTime[m_nNum-1];
	double h=(b-a)/(nParameter-3);
	for(int i=0;i<nParameter;i++)
		T[i]=a-h+i*h;
	double result=0;
	for(i=0;i<nParameter;i++)
		result+=m_pBSpineFit[i]*ggBFuction((x-T[i])/h,nType);
	delete[] T;
	return result;
	
}

//////////////////////////////////////////////////////////////////////////////////////////////////////
//内部函数，计算nType-1次B样条的值,暂时可取取 1 ，2 ，3, 4
//参数: 1 x-------指定一个坐标点
//      2 nType---指定B样条函数基的阶数
//返回：函数值
//////////////////////////////////////////////////////////////////////////////////////////////
double CSpine::ggBFuction(double x,int nType)
{
	switch(nType)
	{
	case 1:
		if(fabs(x)==0.5)
			return 0.5;
		else if(fabs(x)<0.5)
			return 1;
		else
			return 0;
		break;
	case 2:
		if(fabs(x)>1)
			return 0;
		else
			return 1-fabs(x);
		break;
	case 3:
		if(fabs(x)>=1.5)
			return 0;
		else if(fabs(x)<0.5)
			return 0.75-x*x;
		else
			return x*x/2.0-1.5*fabs(x)+9/8.0;
		break;
	case 4:
		if(fabs(x)>=2)
			return 0;
		else if(fabs(x)<1)
			return x*x*fabs(x)/2.0-x*x+2/3.0;
		else 
			return x*x-2*fabs(x)+4/3.0-x*x*fabs(x)/6.0;
		break;
	}
	return -1;

}

//////////////////////////////////////////////////////////////////////////
//三次样条插值函数，可以指定5种不同边界类型,返回参数指针（核心函数，必须调用)
double *CSpine::ThreeSpineInsert(UINT nType,double Border[2])
{
	if(m_p3SpineInsert)
		delete[] m_p3SpineInsert;
	m_p3SpineInsert=new double[4*m_nNum-4];
	ASSERT(m_p3SpineInsert);
	//////////////////////////////////////////////////////////////////////////
	//------------------------------------------------------------------------
	//以下调用了作者以前在做计算方法作业时编的代码
	double *H=new double[m_nNum-1];//步长，
	double *m=new double[m_nNum];//二阶导数
	double *d=new double[m_nNum-1];//差分d[0].....d[N-1]
	double *V=new double[m_nNum-1];//求解二阶导数线性方程组的右端向量V[1]...V[N-1]
	double *a=new double[m_nNum-1];//向量a,b,c,表示三对角方程组的系数
	double *b=new double[m_nNum-1];
	double *c=new double[m_nNum-1];
  //double S[MAX][4]; //存储每一个插值区间多项式的系数       
 
  d[0]=(m_pData[1]-m_pData[0])/H[0];
  for(int k=1;k<m_nNum-1;k++)
  {   H[k]=m_pTime[k+1]-m_pTime[k];
      d[k]=(m_pData[k+1]-m_pData[k])/H[k];//相关初赋值
	  a[k]=H[k];
	  b[k]=2*(H[k-1]+H[k]);
	  c[k]=H[k];
	  V[k]=6*(d[k]-d[k-1]);
  }
 double t;
 
  switch(nType)
  {  
      case ST_CLAMPED:
		  {
		   
		   b[1]=b[1]-H[0]/2;
		   t=d[0]-Border[0];
		   V[1]=V[1]-3*t;
	  	   b[m_nNum-2]=b[m_nNum-2]-H[m_nNum-2]/2;
	       t=Border[1]-d[m_nNum-2];
		   V[m_nNum-2]=V[m_nNum-2]-3*t;
		   break;
		  }
	  case 2:
		  {m[0]=0;
		   m[m_nNum-1]=0;
		   break;
		  }
      case 3:
		  {b[1]=b[1]+H[0]*H[0]/H[1];
		   c[1]=c[1]-H[0]*H[0]/H[1];
		   t=H[m_nNum-2]*H[m_nNum-2]/H[m_nNum-3];
		   b[m_nNum-2]=b[m_nNum-2]+H[m_nNum-2]+t;
		   a[m_nNum-3]=a[m_nNum-3]+t;
		   break;
		  }
      case 4:
		  {b[1]=b[1]+H[0];
		   b[m_nNum-2]=b[m_nNum-2]+H[m_nNum-2];
		   break;
		  }
      case 5:
		  {
           V[1]=V[1]-H[0]*Border[0];
		   V[m_nNum-2]=V[m_nNum-2]-H[m_nNum-2]*Border[1];
		   break;
		  }
  }
  for(k=2;k<=m_nNum-2;k++)//该循环是利用高斯主元法求解三对角方程
  {   t=a[k-1]/b[k-1];
      b[k]=b[k]-t*c[k-1];
	  V[k]=V[k]-t*V[k-1];
  }
  m[m_nNum-2]=V[m_nNum-2]/b[m_nNum-2];
  for(k=m_nNum-3;k>=1;k--)
	  m[k]=(V[k]-c[k]*m[k+1])/b[k];//求出每一点的二阶导数，不包括端点
  switch(nType)
  {       //这个switch语句用来求不同情况下端点的二阶导数
      case 1:
		  {t=(d[0]-Border[0])/H[0];
		   m[0]=3*t-m[1]/2;
		   t=(Border[1]-d[m_nNum-2])/H[m_nNum-2];
		   m[m_nNum-1]=3*t-m[m_nNum-3]/2;
		   break;
		  }
      case 2:
		  {m[0]=m[m_nNum-1]=0;
		  break;
		  }
      case 3:
		  {t=(m[2]-m[1])/H[1];
		   m[0]=m[1]-H[0]*t;
		   double ss=(m[m_nNum-1]-m[m_nNum-3])/H[m_nNum-3];
		   m[m_nNum-1]=m[m_nNum-2]+H[m_nNum-2]*ss;
		   break;
		  }
      case 4:
		  {m[0]=m[1];
		   m[m_nNum-1]=m[m_nNum-2];
		   break;
		  }
      case 5:
		  {m[0]=Border[0];
		   m[m_nNum-1]=Border[1];
		   break;
		  }
  }
  for(k=0;k<=m_nNum-2;k++)//回代求各个区间的三次多项式的系数
  {   m_p3SpineInsert[k*4+0]=m_pData[k];
      t=(2*m[k]+m[k+1])/6;
	 m_p3SpineInsert[k*4+1]=d[k]-H[k]*t;
	  m_p3SpineInsert[k*4+2]=m[k]/2;
	m_p3SpineInsert[k*4+3]=(m[k+1]-m[k])/(6*H[k]);
  }
  
  delete[] H;
  delete[] m;
  delete[] d;
  delete[] V;
  delete[] a;
  delete[] b;
  delete[] c;
  return m_p3SpineInsert;
}

double CSpine::GetThreeSpineInsertValue(double x) const
{
	ASSERT(x>=m_pTime[0]&&x<=m_pTime[m_nNum-1]);
	double z;
	int k;
	double t;
	for(int j=1;j<=m_nNum-1;j++)
	  if((m_pTime[j-1]<=x)&&(m_pTime[j]>=x))
	  {   k=j-1;j=m_nNum-1;
		  double w=x-m_pTime[k];
	      t=m_p3SpineInsert[k*4+3]*w+m_p3SpineInsert[k*4+2];
		  z=(t*w+m_p3SpineInsert[k*4+1])*w+m_p3SpineInsert[k*4+0];
	  }
  
       return z;                  
}