这个是我一个数学老师（教授，数学高手，经常自己做算法）给我的例子，用于多个离散点拟合光滑曲线的，他优化了追赶法，这个例子适用于闭合和不闭合两种情况。当时由于工程情况，写的急，代码不好看，但是很好用。为了方便传递参数，我做了一个链表，用时候根据自己情况可以修改，核心算法不动即可。
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; file://赋初值 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) { file://闭合 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) { file://x连续 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 { file://y连续 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; }
