曲线拟合

1.最小二乘法

using System;

namespace eking.generic
{
 /// <summary>
 /// CurveSimulate 的摘要说明。
 /// 基于最小二乘法拟合曲线
 /// eking   2005-07-08  my birthday
 /// </summary>
 public abstract class CurveSimulate
 {
  public CurveSimulate()
  {
  }

  /// <summary>
  /// 拟合曲线－基于最小二乘法
  /// </summary>
  /// <param name="X">X两轴的坐标</param>
  /// <param name="Y">Y两轴的坐标</param>
  /// <param name="Result">结果参数数组</param>
  /// <returns></returns>
  public static bool CalculateCurveParameter(double[] X,double[] Y, ref double[] Result)
  {
   long i,j,k;
   double Z,D1,D2,C,P,G,Q=0;

   long N=X.Length;
   long M=Result.Length;
   
   double[] B=new double[N];
   double[] T=new double[N];
   double[] S=new double[N];

   if(M>N)
    M=N;
   for(i=0;i<M;i++)
    Result[i]=0;
   Z=0;
   B[0]=1;
   D1=N;
   P=0;
   C=0;
   for(i=0;i<N;i++)
   {
    P=P+X[i]-Z;
    C=C+Y[i];
   }
   C=C/D1;
   P=P/D1;
   Result[0]=C*B[0];
   if(M>1)
   {
    T[1]=1;
    T[0]=-P;
    D2=0;
    C=0;
    G=0;
    for(i=0;i<N;i++)
    {
     Q=X[i]-Z-P;
     D2=D2+Q*Q;
     C=Y[i]*Q+C;
     G=(X[i]-Z)*Q*Q+G;
    }
    C=C/D2;
    P=G/D2;
    Q=D2/D1;
    D1=D2;
    Result[1]=C*T[1];
    Result[0]=C*T[0]+Result[0];
   }
   for(j=2;j<M;j++)
   {
    S[j]=T[j-1];
    S[j-1]=-P*T[j-1]+T[j-2];
    if(j>=3)
    {
     for(k=j-2;k>=1;k--)
      S[k]=-P*T[k]+T[k-1]-Q*B[k];
    }
    S[0]=-P*T[0]-Q*B[0];
    D2=0;
    C=0;
    G=0;
    for(i=0;i<N;i++)
    {
     Q=S[j];
     for(k=j-1;k>=0;k--)
      Q=Q*(X[i]-Z)+S[k];
     D2=D2+Q*Q;
     C=Y[i]*Q+C;
     G=(X[i]-Z)*Q*Q+G;
    }
    C=C/D2;
    P=G/D2;
    Q=D2/D1;
    D1=D2;
    Result[j]=C*S[j];
    T[j]=S[j];
    for(k=j-1;k>=0;k--)
    {
     Result[k]=C*S[k]+Result[k];
     B[k]=T[k];
     T[k]=S[k];
    }
   }
   return true;
  }
 }
}

2. 光滑曲线拟合

这个是我一个数学老师（教授，数学高手，经常自己做算法）给我的例子，用于多个离散点拟合光滑曲线的，他优化了追赶法，这个例子适用于闭合和不闭合两种情况。当时由于工程情况，写的急，代码不好看，但是很好用。为了方便传递参数，我做了一个链表，用时候根据自己情况可以修改，核心算法不动即可。

class CFoldPoint

{public:

    double X;    double Y;

};

typedef CTypedPtrList CFoldPointList;

typedef CArray CDoubleArray;

三个函数，SPLine 调用另外两个。用时候直接调用SPLine函数，入口pList是已知离散点链表，pDestList是生成的点的链表。SM是在两个点中间插入点的数目，continue=0是采样点无规律，要求生成闭合曲线。1是采样点x坐标连续 2是y连续

void ZG(CDoubleArray *A,CDoubleArray *B,CDoubleArray *C,CDoubleArray *G,int &LOGI)
{
      //追赶法
 register long I;
 int N;
 N=A->GetSize();
 if(LOGI==0)
 {
  (*C)[0]=(*C)[0]/(*B)[0];
   for(I=1;I   {
    (*B)[I]=(*B)[I]-(*A)[I]*(*C)[I-1];
   (*C)[I]=(*C)[I]/(*B)[I];
  }
  (*A)[0]=0.;
   (*C)[N-1]=0.;
   LOGI=1;
  }
  (*G)[0]=(*G)[0]/(*B)[0];
  for(I=1;I  {
   (*G)[I]=((*G)[I]-(*A)[I]*(*G)[I-1])/(*B)[I];
  }
  for(I=N-2;I>-1;I--)//DO 30 I=N-1,1,-1
  {
   (*G)[I]=(*G)[I]-(*C)[I]*(*G)[I+1];
  }
  return;
}

void SPLine4(CDoubleArray *X,CDoubleArray *Y,double &XI,double&YI,CDoubleArray *A,CDoubleArray *B,CDoubleArray *C,CDoubleArray *G,int &LOGI,int MD)
{
 
 register long I;
 double W1,W2,H;
 int N=X->GetSize();
 
 if(LOGI==0)
 {
  for(I=1;I  {
   (*B)[I]=(*X)[I]-(*X)[I-1];
   (*C)[I]=((*Y)[I]-(*Y)[I-1])/(*B)[I];
  }
  for(I=1;I  {
   (*A)[I]=(*B)[I]+(*B)[I+1];
   (*G)[I]=6.*((*C)[I+1]-(*C)[I])/(*A)[I];
   (*A)[I]=(*B)[I]/(*A)[I];
  }
  for(I=1;I  {
   (*C)[I]=1.-(*A)[I];
   (*B)[I]=2.;
  }
  (*B)[0]=2.;
  (*B)[N-1]=2.;
  if(MD==3)
  {
   (*C)[0]=-1.;
   (*A)[N-1]=-1.;
   (*A)[0]=0.;
   (*C)[N-1]=0.;
  }
  ZG(A,B,C,G,LOGI);
 }
 for(I=1;I {
  if(XI>=(*X)[I-1] && XI<=(*X)[I])//GE LE
  {
   H=(*X)[I]-(*X)[I-1];
   W1=(*X)[I]-XI;
   W2=XI-(*X)[I-1];
   YI=W1*W1*W1*(*G)[I-1]/6./H;
   YI=YI+W2*W2*W2*(*G)[I]/6./H;
   YI=YI+W1*((*Y)[I-1]-(*G)[I-1]*H*H/6.)/H;
   YI=YI+W2*((*Y)[I]-(*G)[I]*H*H/6.)/H;
  }
 }
}
void SPLine(CFoldPointList *pList,CFoldPointList *pDestList,int SM,int Continue=0)
{
 CFoldPoint *pFoldHead,*pFoldTail;
 POSITION pos;
 CDoubleArray A,B,C,G,X,Y,T;
 double XI,YI,XX,YY;
 register long i;
 long N;
 int LOGI;
 long RealSM;
 long Bei,Yu;
 CFoldPoint *pFold;
 //初值
 N=pList->GetCount();
 A.SetSize(N);
 B.SetSize(N);
 C.SetSize(N);
 G.SetSize(N);
 X.SetSize(N);
 Y.SetSize(N);
 T.SetSize(N);
 RealSM=(N-1)*SM+N;
 pos=pList->GetHeadPosition();
 for(i=0;i {
  pFold=pList->GetNext(pos);
  X[i]=pFold->X;
  Y[i]=pFold->Y;
 }
 
 pFoldHead=pList->GetHead();
 pFoldTail=pList->GetTail();
 if(Continue==0)//pFoldHead->X==pFoldTail->X && pFoldHead->Y==pFoldTail->Y)
 { //闭合
  T[0]=0;
  for(i=0;i  {
   T[i+1]=T[i]+CalculateDistance(X[i],Y[i],X[i+1],Y[i+1])+0.000000001;
  }
  LOGI=0;
  YI=0;
  for(i=0;i  {
   Bei=i/(SM+1);
   Yu=i%(SM+1);
   if(Yu!=0)
   {
    XI=T[Bei]+(T[Bei+1]-T[Bei])/(SM+1)*Yu;
    SPLine4(&T,&Y,XI,YI,&A,&B,&C,&G,LOGI,3);
    YY=YI;//+Y[Bei];
   }
   else
   {
    YY=Y[Bei];
   }
   pFold=new CFoldPoint;
   pFold->Y=YY;
   pDestList->AddTail(pFold);
  }
  LOGI=0;
  YI=0;
  pos=pDestList->GetHeadPosition();
  for(i=0;i  {
   Bei=i/(SM+1);
   Yu=i%(SM+1);
   if(Yu!=0)
   {
    XI=T[Bei]+(T[Bei+1]-T[Bei])/(SM+1)*Yu;
    SPLine4(&T,&X,XI,YI,&A,&B,&C,&G,LOGI,3);
    YY=YI;//+X[Bei];
   }
   else
   {
    YY=X[Bei];
   }
   pFold=pDestList->GetNext(pos);
   pFold->X=YY;
  }
 }
 else if(Continue==1)
 {
//连续
  LOGI=0;
  YI=0;
  for(i=0;i  {
   Bei=i/(SM+1);
   Yu=i%(SM+1);
   if(Yu!=0)
   {
    XI=X[Bei]+(X[Bei+1]-X[Bei])/(SM+1)*Yu;
    SPLine4(&X,&Y,XI,YI,&A,&B,&C,&G,LOGI,3);
    XX=XI;
    YY=YI;
   }
   else
   {
    XX=X[Bei];
    YY=Y[Bei];
   }
   pFold=new CFoldPoint;
   pFold->X=XX;
   pFold->Y=YY;
   pDestList->AddTail(pFold);
  }
 }
 else
 {
//连续
  LOGI=0;
  YI=0;
  for(i=0;i  {
   Bei=i/(SM+1);
   Yu=i%(SM+1);
   if(Yu!=0)
   {
    XI=Y[Bei]+(Y[Bei+1]-Y[Bei])/(SM+1)*Yu;
    SPLine4(&Y,&X,XI,YI,&A,&B,&C,&G,LOGI,3);
    XX=YI;
    YY=XI;
   }
   else
   {
    XX=X[Bei];
    YY=Y[Bei];
   }
   pFold=new CFoldPoint;
   pFold->X=XX;
   pFold->Y=YY;
   pDestList->AddTail(pFold);
  }
 }
 return;
}

那个例子中用了CFoldPointList和CDoubleArray

class CFoldPoint

{

public:

  double x;double y;

}

typedef CTypedPtrList<CPtrList,CFoldPoint*> CFoldPointList;
typedef CDoubleArray<double,double> CDoubleArray;