300题目：求B样条曲线的控制点的算法！！

[nurbs]

这是我第三次开新贴，头两次都没有解决，这次如有能解决者，连同前两次一并给分！
http://expert.csdn.net/Expert/topic/1292/1292306.xml?temp=.3483698
http://expert.csdn.net/Expert/topic/1290/1290860.xml?temp=.2744562

问题如下：
我已有画B样条曲线的deboor算法，但需要知道控制点。但现在是反其道而行：求控制点。在鼠标绘图程序中，鼠标按中的点应为型值点，有什么算法可求出其控制点？
不要说哪本书上有，我没学过图形学，而且高数在学完以后都还给老师了！
最好提供代码！！

[nurbs]

up

[nurbs]

up

[FatGarfield]

实际上，Bazier曲线可以表示为如下形式，如：三次Bazier曲线
x = a0(t) * x0 + a1(t) * x1 + a2(t) * x2 + a3(t) * x3
y = a0(t) * y0 + a1(t) * y1 + a2(t) * y2 + a3(t) * y3
z = a0(t) * z0 + a1(t) * z1 + a2(t) * z2 + a3(t) * z3
其中， （x,y,z)为曲线上的点，(xi, yi, zi)为控制点,
a0(t)~a3(t)为Bernstein多项式，定义为：
ai(t) = [C(3,i)]* [t^i ] * [(1-t) ^( 3-i )] //C(3,i)表示组合
那么，实际上，上面三个参数方程可以化减为t的多项式方程
x = ( x0 - 3 * x1 + 3 * x2 - x3 )* t^3 +( 3 * x1 - 6 * x2 +
3 * x3 ) * t^2 + ( 3 * x2 -3 * x3 ) * t + x3 （1）
y = ( y0 - 3 * y1 + 3 * y2 - y3 )* t^3 +( 3 * y1 - 6 * y2 +
3 * y3 ) * t^2 + ( 3 * y2 -3 * y3 ) * t + y3 （2）
z = ( z0 - 3 * z1 + 3 * z2 - z3 )* t^3 +( 3 * z1 - 6 * z2 +
3 * z3 ) * t^2 + ( 3 * z2 -3 * z3 ) * t + z3 （3）
那么，现在的问题是给出一组曲线上的点(x,y,z)，需要求控制点的值
对于x,y或z,只要求出一组即可
方法是：
1.给出一组曲线上的点(x,y,z)，首先求出多项式插值曲线，比如，给出4个点，可以求出3次参数多项式插值曲线
x = b0 * u + b1 * u^2 + b2 * u^2 + b3 * u^3 （4）
将(1)与(4)比较，可得：

x0 - 3 * x1 + 3 * x2 - x3 = b0
3 * x1 - 6 * x2 + 3 * x3 = b1
3 * x2 -3 * x3 = b3
x3 = b4

4个方程，4个未知数，可以求得x0~x3
同理，可以求得y0~y3和z0~z3

[painache]

孙家广的《计算机图形学》上讲的很清楚。

高等数学的知识在绘制B-SPLINE的时候几乎用不到，如果你要计算曲线的长度或者要实现Snap功能的话，高数的知识必不可少。线形代数要懂，并且自己写程序去实践，总可以搞出来。

PS：我有VB的代码，完全符合你的要求，要的话我发给你。

[nurbs]

[e-mail]
谢谢

[fdz81]

我也要，也给我一份吧！谢谢！
[e-mail]

[painache]

过几天发，我现在不在家

[askcsdn]

我以前做过这样的VC程序，可以在n(n>=2)个点上画曲线，用鼠标绘图即可，将鼠标点中的点坐标代入程序就能画出曲线。
需说明一点：我画的是三次参数样条插值曲线，和B样条差不多，不只是不是你想要的，要得话给出邮箱

[cadinfo]

这是我本科时候辅修课程的作业，非常的简单，不过唯一的要求是必须给出起点和终点的法矢，否则方程组不定，这样你用0，0即可，因为B样条和三次样条不同，B样条有很好的性质，性质中的哪条就不记得了，也许叫做保形吧！

[dingdongfly]

给我也发一个好吗？[e-mail]

[jvsun]

http://www.csdn.net/cnshare/soft/15/15307.shtm

[vcshox]

反算过程要求给出边界条件，用的最多的是切矢边界条件。以非均匀的3次为例，只需要解方程组：(简化后的)
|b2 c2 ||d1 | |e1 |
|a3 b3 c3 ||d2 | |e2 |
| ...... ||.. | = |.. |
| a(n-1) b(n-1) c(n-1)||d(n) | |e(n-1)|
| a(n) b(n) ||d(n+1)| |e(n) |
其中：d0=p0 d(n+2)=p(n) d(i):控制点 p(i):B样条上的型值点
△(i)=u(i+1)-u(i)
a(i)=△(i+2)*△(i+2)/(△(i)+△(i+1)+△(i+1))
b(i)=△(i+2)*(△(i)+△(i+1))/(△(i)+△(i+1)+△(i+2))+
△(i+1)*(△(i+2)+△(i+3))/(△(i+1)+△(i+2)+△(i+3))
c(i)=△(i+1)*△(i+1)/(△(i+1)+△(i+2)+△(i+3))
e(i)=(△(i+1)+△(i+2))*p(i-1) i=2...n
e1=p0+△(3)*p0'/3 p0'=3(d1-d0)/△3
e(n+1)=p(n)-△(n+2)*p'(n)/3 p'(n)=3(d(n+2)-d(n+1))/△(n+2)
对于不同的边界条件，方程组的主体是一样的，只涉及边界上的方程的变化。
源程序可以到天网和GOOGLE上找。

[vcshox]

SORRY,上面的方程最后一列是从e1.....e(n),e(n+1)
另外请 painache() 把VB的源代码也给我发一份。[e-mail]
谢谢

[vcshox]

不好意思，我仔细看了一下方程，应该是：
d2.....d(n)和e(2)....e(n)

[painache]

已发，请查收

大家交个朋友如何？

QQ：23679411

[cadinfo]

我这里有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);
}