(整理)优化了的过关键点的光滑曲线拟合算法.

优化了的过关键点的光滑曲线拟合算法

作者：中国论坛网收集来源： 加入时间：2004-8-25

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

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;