huangyhg

求B样条曲线的控制点的算法！！
>

我这里有vc的代码，先反求控制顶点，然后用不同的方法拟合，比如deboor
我自己写的样本程序参考如下：（参考书籍就是那本计算机辅助设计B样条）
-----------------------------------------------------
void CSpline::RealizeD(CDC*pDC)
{
int i; //存放循环变量
double w; //存放临时分母
double temp[VER_MAX+2]={0}; //存放解三对角方程时用到的临时数组

CPoint P[VER_MAX+2]={CPoint(0,0)}; //存放三对角方程系数D，反求控制顶点

/////////////////////////////////////////////////////////////////////////////////
// 非均匀3次B样条曲线节点划分 //
double x[VER_MAX+6]={0}; //存放节点矢量
double l[VER_MAX]={0}; //存放型值点间距
double L=0;
for(i=1;i<=n;i++)
{
l[i]=sqrt((v[i].x-v[i-1].x)*(v[i].x-v[i-1].x)
+(v[i].y-v[i-1].y)*(v[i].y-v[i-1].y));
L+=l[i]; //计算控制顶点总距
}
for(i=1;i<n;i++)
x[i+3]=x[i+3-1]+l[i]/L;
x[10]=x[11]=x[12]=x[13]=1;
///////////////////////////////////////////////////

double m11[VER_MAX-1],m12[VER_MAX-1],m13[VER_MAX-1],m14[VER_MAX-1],
m21[VER_MAX-1],m22[VER_MAX-1],m23[VER_MAX-1],m24[VER_MAX-1],
m31[VER_MAX-1],m32[VER_MAX-1],m33[VER_MAX-1],m34[VER_MAX-1],
m41[VER_MAX-1],m42[VER_MAX-1],m43[VER_MAX-1],m44[VER_MAX-1];
for(i=0;i<n;i++)
{
m11[i]=((x[i+3]-x[i+4])*(x[i+3]-x[i+4]))/((x[i+4]-x[i+1])*(x[i+4]-x[i+2]));
m13[i]=((x[i+2]-x[i+3])*(x[i+2]-x[i+3]))/((x[i+2]-x[i+5])*(x[i+2]-x[i+4]));
m12[i]=1-(m11[i]+m13[i]);
m14[i]=0;

m21[i]=-3.0*m11[i];
m23[i]=3.0*((x[i+4]-x[i+3])*(x[i+3]-x[i+2]))/((x[i+2]-x[i+5])*(x[i+2]-x[i+4]));
m22[i]=-(m21[i]+m23[i]);
m24[i]=0;

m31[i]= 3.0*m11[i];
m33[i]=3.0*((x[i+4]-x[i+3])*(x[i+4]-x[i+3]))/((x[i+2]-x[i+5])*(x[i+2]-x[i+4]));
m32[i]=-(m31[i]+m33[i]);
m34[i]=0;

m41[i]=-m11[i];
m43[i]=-((x[i+4]-x[i+3])*(x[i+4]-x[i+3]))*(1/((x[i+2]-x[i+4])*(x[i+2]-x[i+5]))+1/((x[i+3]-x[i+5])*(x[i+3]-x[i+6]))+1/((x[i+5]-x[i+2])*(x[i+5]-x[i+3])));
m44[i]=((x[i+4]-x[i+3])*(x[i+4]-x[i+3]))/((x[i+3]-x[i+6])*(x[i+3]-x[i+5]));
m42[i]=-(m41[i]+m43[i]+m44[i]);
}

double matrix[VER_MAX+2][VER_MAX+2]={0.0};
for(i=1;i<=n;i++)
{
matrix[i][i-1]=m11[i-1];
matrix[i][i]=m12[i-1];
matrix[i][i+1]=m13[i-1];
}
matrix[0][0]=m21[0]/(x[4]-x[3]);
matrix[0][1]=m22[0]/(x[4]-x[3]);
matrix[0][2]=m23[0]/(x[4]-x[3]);

matrix[n+1][n-1]=m11[n-1]+m21[n-1]+m31[n-1]+m41[n-1];
matrix[n+1][n]=m12[n-1]+m22[n-1]+m32[n-1]+m42[n-1];
matrix[n+1][n+1]=m13[n-1]+m23[n-1]+m33[n-1]+m43[n-1];
matrix[n+1][n+2]=m14[n-1]+m24[n-1]+m34[n-1]+m44[n-1];

matrix[n+2][n-1]=(m21[n-1]+2*m31[n-1]+3*m41[n-1])/(x[n+3]-x[n+2]);
matrix[n+2][n]=(m22[n-1]+2*m32[n-1]+3*m42[n-1])/(x[n+3]-x[n+2]);
matrix[n+2][n+1]=(m23[n-1]+2*m33[n-1]+3*m43[n-1])/(x[n+3]-x[n+2]);
matrix[n+2][n+2]=(m24[n-1]+2*m34[n-1]+3*m44[n-1])/(x[n+3]-x[n+2]);
//////////////////以上代码赋值，以下代码整理成系数矩阵
matrix[n+2][n]=matrix[n+2][n]-matrix[n+1][n]*(matrix[n+2][n]/matrix[n+1][n]);
matrix[n+2][n+1]=matrix[n+2][n+1]-matrix[n+1][n+1]*(matrix[n+2][n]/matrix[n+1][n]);
matrix[n+2][n+2]=matrix[n+2][n+2]-matrix[n+1][n+2]*(matrix[n+2][n]/matrix[n+1][n]);

//整理系数，化成n+3阶方程
for(i=1;i<=n+1;i++)
{
P[i].x=v[i-1].x;
P[i].y=v[i-1].y;
}
P[0].x=long(1000/sqrt(1+m[0]*m[0]));
P[0].y=long(1000*m[0]/sqrt(1+m[0]*m[0]));
P[n+2].x=long(1000/sqrt(1+m[n]*m[n]));
P[n+2].y=long(1000*m[n]/sqrt(1+m[n]*m[n]));

///////////////// 根据端点条件，可得出CtrlVertex[0]和CtrlVertex[n+2]
///////////////// 并计算 CtrlVertex[1]CtrlVertex[n+1]
CtrlVertex[0]=P[1];
CtrlVertex[1].x=long((P[0].x-matrix[0][0]*CtrlVertex[0].x)/matrix[0][1]);
CtrlVertex[1].y=long((P[0].y-matrix[0][0]*CtrlVertex[0].y)/matrix[0][1]);

CtrlVertex[n+2]=P[n+1];
CtrlVertex[n+1].x=long((P[n+2].x-matrix[n+2][n+2]*CtrlVertex[n+2].x)/matrix[n+2][n+1]);
CtrlVertex[n+1].y=long((P[n+2].y-matrix[n+2][n+2]*CtrlVertex[n+2].y)/matrix[n+2][n+1]);
///化三对角方程
P[2].x=long(P[2].x-matrix[2][1]*CtrlVertex[1].x);
P[2].y=long(P[2].y-matrix[2][1]*CtrlVertex[1].y);

P[n].x=long(P[n].x-matrix[n][n+1]*CtrlVertex[n+1].x);
P[n].y=long(P[n].y-matrix[n][n+1]*CtrlVertex[n+1].y);

/////////////////////////////////////////////////////////////////
//////// 解三对角方程，反求控制顶点CtrlVertex. 0～n+2 /////
w=matrix[2][2];
CtrlVertex[2].x=long(P[2].x/w);
CtrlVertex[2].y=long(P[2].y/w);
//消去过程
for(i=3;i<=n;i++)
{
temp[i-1]=matrix[i-1][i]/w;
w=matrix[i][i]-matrix[i][i-1]*temp[i-1];
CtrlVertex[i].x=long((P[i].x-matrix[i][i-1]*CtrlVertex[i-1].x)/w);
CtrlVertex[i].y=long((P[i].y-matrix[i][i-1]*CtrlVertex[i-1].y)/w);
}
//回代过程
for(i=n-1;i>=2;i--)
{
CtrlVertex[i].x=long(CtrlVertex[i].x-temp[i]*CtrlVertex[i+1].x);
CtrlVertex[i].y=long(CtrlVertex[i].y-temp[i]*CtrlVertex[i+1].y);
}
///////////////////////////////////////////////////////////////
/////////在图上标出控制顶点CtrlVertex. 0～n+2 /////
///////////////////////////////////////////////////////////
CPen *oldpp,dotpen(PS_DASH,1,RGB(100,100,100));
oldpp=pDC->SelectObject(&dotpen);
for(i=0;i<=n+2;i++)
{
pDC->MoveTo(CtrlVertex[i]);
CRect rect(CtrlVertex[i].x-20,CtrlVertex[i].y-20,
CtrlVertex[i].x+20,CtrlVertex[i].y+20);
pDC->FillRect(rect,&CBrush(RGB(255,0,0)));
CString str;
str.Format("V[%d]",i);
pDC->SetTextColor(RGB(0,100,0));
pDC->TextOut(CtrlVertex[i].x,CtrlVertex[i].y-30,str);
}
pDC->SelectObject(oldpp);
//////////////以下代码连接控制顶点///////////////////////
// pDC->MoveTo(CtrlVertex[0]);
// for(i=1;i<=n+2;i++)
// pDC->LineTo(CtrlVertex[i]);
// pDC->SelectObject(oldpp);


///////////////////////////////////////////////////////////////////////////////////////////////
//////////////////运用DeBoor算法计算插值点
double t; //累加变量
double alpha[VER_MAX+2][4]; //存放临时变量
CPoint dp[VER_MAX+2][4]={CPoint(0,0)}; //存放计算插值点
//////// 绘制n段曲线 ////////////////////////////////////
CPen *oldPen,newPen;
newPen.CreatePen(PS_SOLID,1,RGB(0,0,255));
oldPen=pDC->SelectObject(&newPen);

pDC->MoveTo(v[0]);
for(i=3;i<n+3;i++)
for(t=x[i];t<=x[i+1];t+=(x[i+1]-x[i])/10)
for(int m=0;m<=3;m++)
for(int j=i-3+m;j<=i;j++)
{
if(m==0) dp[j][m]=CtrlVertex[j];
else
{

alpha[j][m]=(t-x[j])/(x[j+4-m]-x[j]);
dp[j][m].x=long((1-alpha[j][m])*dp[j-1][m-1].x+alpha[j][m]*dp[j][m-1].x);
dp[j][m].y=long((1-alpha[j][m])*dp[j-1][m-1].y+alpha[j][m]*dp[j][m-1].y);
if(j==i&&m==3)
pDC->LineTo(dp[j][3]);
}
}

pDC->SetTextColor(RGB(255,0,255));
pDC->TextOut(8500,1700,"非均匀B-样条曲线：");
pDC->SelectObject(oldPen);
}

__________________
借用达朗贝尔的名言：前进吧，你会得到信心!
[url="http://www.dimcax.com"]几何尺寸与公差标准[/url]