B样条曲线曲面和NURBS曲线曲面C语言算法源程序.doc

学习小结：前面学习了Bezier曲线，B样条基函数和B样条曲线的一些基础知识。掌握关键问题是一条B样条曲线间的多段曲线的光滑连接。因为现在是用多段Bezier曲线来描绘一条B样条曲线，所以问题变为两段Bezier曲线间光滑连接。两段Bezier曲线段（3次）B1和B2光滑连接的条件： （1）.要求B1和B2有共同的连接点，即G0连续。 （2）.要求B1和B2在连接点处有成比例的一阶导数，即G1连续。由端点处的一阶导数,为实现G1连续，则有： 即： 这也表明，三点共线。如下图表示了一条3次B样条曲线的所有控制多边形： (P1) P2 P3 P4 (P11) (P12) P5 P10 P0 P6 P9 P7 P8 图5.3次B样条曲线和所有控制多边形 图5中，P0至P6为原始3次B样条曲线控制多边形顶点，P0至P12是计算后最终形成B样条曲线控制多边形顶点。 图6.双二次（2x2）B样条曲面 6.B样条曲线曲面和NURBS曲线曲面的C语言实现算法源程序#ifndef _mynurbs_h#ifndef _MYNURBS_H#include glgl.h#include math.h/*-*-*-*-*-*-*-*-*-*-*-*-*-*-* B样条基函数计算部分 *-*-*-*-*-*-*-*-*-*-*-*-*-*/确定参数u所在的节点区间下标/n=m-p-1 /m为节点矢量U的最大下标/p为B样条函数次数int FindSource(int n,int p,float u,float U)int low,high,mid;if(u=Un+1)/特殊情况return n;/进行二分搜索low=p;high=n+1;mid=(int)(low+high)/2;while(uUmid)if(u=Umid&uUmid+1)/找到u所在的节点区间的下标 break; /退出二分搜索return mid; /返回区间下标/计算所有非零B样条基函数并返回其值/i为参数u所在的节点区间下标void BasisFunction(int i,int p,float u,float U,float N) int j,di,dp,k;float tul,tur,left,right;float tmpN5050;for(k=0;k=i-k;di-) if(u=Udi&uUdi+1) tmpNdi0=1; else tmpNdi0=0; dp+=1; for(j=1;j2void DerBasisFunc(int i,int p,float u,float U,float NP)int j,di,dp,k;float tul,tur,left,right,saved,dl,dr;float tmpN5050;for(k=0;k=i-k;di-) if(u=Udi&uUdi+1) tmpNdi0=1; else tmpNdi0=0; dp+=1; for(j=1;jdp;j+) tul=Udi+j-Udi; tur=Udi+j+1-Udi+1; if(tul!=0) left=(u-Udi)/tul,dl=1/tul; else left=0,dl=0; if(tur!=0) right=(Udi+j+1-u)/tur,dr=1/tur; else right=0,dr=0; tmpNdij=(left*tmpNdij-1+right*tmpNdi+1j-1); saved=p*(dl*tmpNdij-1-dr*tmpNdi+1j-1)/(p+p-1); NPi-k=saved;/*-*-*-*-*-*-*-*-*-*-*-*-*-* Bezier曲线曲面部分 *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*/计算参数u的p次基函数值并存在BC中void BernsteinFunc(int p,double t,float BC)for(int i=0;i0&ip)BCi=p*(float)pow(t,i)*(float)pow(1-t,p-i);/获取p次Bezier曲线上的lines个点的值void BezierPoint(int p,float px,float py,float pz,int lines,float tmp3)float BC20;int i,j; for(j=0;j=lines;j+) double t=j/(float)lines; BernsteinFunc(p,t,BC); tmpj0=tmpj1=tmpj2=0; for(i=0;ip+1;i+) tmpj0+=BCi*pxi; tmpj1+=BCi*pyi; tmpj2+=BCi*pzi; /获取p次有理Bezier曲线上的lines个点的值void NBezierPoint(int p,float px,float py,float pz,float pw,int lines,float tmp4)float x,y,z,w,BC20;int i,j; for(j=0;j=lines;j+) double t=j/(float)lines; BernsteinFunc(p,t,BC); x=y=z=w=0; for(i=0;ip+1;i+) x+=BCi*pxi*pwi; y+=BCi*pyi*pwi; z+=BCi*pzi*pwi; w+=BCi*pwi; tmpj0=x/w; tmpj1=y/w; tmpj2=z/w; tmpj3=w;/-/绘制p次的Bezier曲线void Bezier(int p,float px,float py,float pz,int lines)float pt1003;int j; BezierPoint(p,px,py,pz,lines,pt);for(j=1;j=lines;j+) glBegin(GL_LINES); glVertex3f(ptj-10,ptj-11,ptj-12); glVertex3f(ptj0,ptj1,ptj2); glEnd();/-/绘制p次的有理Bezier曲线void NBezier(int p,float px,float py,float pz,float w,int lines)float pt1004;int j; NBezierPoint(p,px,py,pz,w,lines,pt);for(j=1;j=lines;j+) glBegin(GL_LINES); glVertex3f(ptj-10,ptj-11,ptj-12); glVertex3f(ptj0,ptj1,ptj2); glEnd();/-/计算双p次Bezier曲面上所有的点并保存在Pt中/u和v分别为曲面(u,v)方向上的网格数void BezierFacePoint(int p,int u,int v,float px4,float py4,float pz4,float pt1611613)float urx11161,ury11161,urz11161;float tx11,ty11,tz11,tmp1613;int i,j,k; for(j=0;jp+1;j+) for(i=0;ip+1;i+) txi=pxij; tyi=pyij; tzi=pzij; BezierPoint(p,tx,ty,tz,v,tmp); for(k=0;k=v;k+) urxjk=tmpk0; uryjk=tmpk1; urzjk=tmpk2; for(i=0;i=v;i+) for(k=0;kp+1;k+) txk=urxki; tyk=uryki; tzk=urzki; BezierPoint(p,tx,ty,tz,u,tmp); for(j=0;j=u;j+) ptij0=tmpj0; ptij1=tmpj1; ptij2=tmpj2; /-/计算双p次有理Bezier曲面上所有的点并保存在Pt中/u和v分别为曲面(u,v)方向上的网格数void NuBezierFacePoint(int p,int u,int v,float px4,float py4,float pz4,float w4,float pt1611613)float urx11161,ury11161,urz11161,urw11161;float tx11,ty11,tz11,tw11,tmp1614;int i,j,k; for(j=0;jp+1;j+) for(i=0;ip+1;i+) txi=pxij; tyi=pyij; tzi=pzij; twi=wij; NBezierPoint(p,tx,ty,tz,tw,v,tmp); for(k=0;k=v;k+) urxjk=tmpk0; uryjk=tmpk1; urzjk=tmpk2; urwjk=tmpk3; for(i=0;i=v;i+) for(k=0;kp+1;k+) txk=urxki; tyk=uryki; tzk=urzki; twk=urwki; NBezierPoint(p,tx,ty,tz,tw,u,tmp); for(j=0;j=u;j+) ptij0=tmpj0; ptij1=tmpj1; ptij2=tmpj2; /-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* B样条曲线曲面部分 -*-*-*-*-*-*-*-*-*-*-*-*-*-*-/计算样条曲线的1阶导矢（u所对应的所有点）保存在Der中/n=m-p-1/p为曲线的次数void BSplineDer(int n,int p,float U,float P,float Der)float N100,tmp; int i,j; for(i=p+1;i=i-p;j-)tmp+=Nj*Pj;Deri-p=tmp;/计算曲线上的点（u所对应的所有点）保存在Poi中/n=m-p-1/p为曲线的次数void BSplinePoint(int n,int p,float U,float P,float Poi)float N100,tmp; int i,j; for(i=p+1;i=i-p;j-)tmp+=Nj*Pj;Poii-p=tmp;/计算3次样条曲线上的所有控制多边形保存在CP中/m为节点矢量U的最大下标void B3SplineControlPoint(int m,float U,float P,float CP)int n,k,i,cp,p;float Poi100,Der100,add;p=3;n=m-p-1;BSplinePoint(n,p,U,P,Poi);BSplineDer(n,p,U,P,Der);cp=(n-p)*3+p;for(i=0;i2;i+)CPi=Pi;CPcp-i=Pn-i;for(i=3;icp-1;i+=3) k=(int)i/3;add=Derk/p;CPi=Poik;CPi-1=CPi-add;CPi+1=CPi+add;/计算2次样条曲线上的所有控制多边形保存在CP中/m为节点矢量U的最大下标void B2SplineControlPoint(int m,float U,float P,float CP)int n,k,tm,i,cp,p;float Poi100;p=2;n=m-p-1;BSplinePoint(n,p,U,P,Poi);cp=(n-p)*2+p;for(i=0;i2;i+)CPi=Pi;CPcp=Pn;tm=2;for(i=2;icp-1;i+=2) k=(int)i/2;CPi=Poik;CPi+1=Ptm;tm+;/绘制3次B样条曲线/m为节点矢量U的最大下标void BSpline3L(int m,float U,float px,float py,float pz)float pcx100,pcy100,pcz100,drx4,dry4,drz4;int i,j,tmcp;B3SplineControlPoint(m,U,px,pcx);B3SplineControlPoint(m,U,py,pcy);B3SplineControlPoint(m,U,pz,pcz);/*glColor3f(0.0f,0.0f,0.0f);for(i=1;i3*m-17;i+)glBegin(GL_LINES);glVertex3f(pcxi-1,pcyi-1,pczi-1);glVertex3f(pcxi,pcyi,pczi);glEnd(); glColor3f(1.0f,0.0f,0.0f);*/tmcp=m-7;for(i=0;i=tmcp;i+)for(j=i*3;ji*3+4;j+)drxj-i*3=pcxj;dryj-i*3=pcyj;drzj-i*3=pczj;Bezier(3,drx,dry,drz,20);/绘制2次B样条曲线/m为节点矢量U的最大下标void BSpline2L(int m,float U,float px,float py,float pz)float pcx100,pcy100,pcz100,drx3,dry3,drz3;int i,j,tmcp;B2SplineControlPoint(m,U,px,pcx);B2SplineControlPoint(m,U,py,pcy);B2SplineControlPoint(m,U,pz,pcz);tmcp=m-5;for(i=0;i=tmcp;i+)for(j=i*2;ji*2+3;j+)drxj-i*2=pcxj;dryj-i*2=pcyj;drzj-i*2=pczj;Bezier(2,drx,dry,drz,20);/计算双三次（3x3）B样条曲面所有控制多边形顶点，并保存在pt中/mu,mv分别为节点矢量U,V的最大下标值void BS3FaceControlPoint(int mu,float U,int mv,float V,float px,float py,float pz,float pt1001003)int i,j,k,dp;float tmx50,tmy50,tmz50;float tmpx50100,tmpy50100,tmpz50100;float uvx100100,uvy100100,uvz100100;for(i=0;imv-3;i+)dp=i*(mu-3);for(j=dp;jmu-3+dp;j+)tmxj-dp=pxj;tmyj-dp=pyj;tmzj-dp=pzj;B3SplineControlPoint(mu,U,tmx,tmpxi);B3SplineControlPoint(mu,U,tmy,tmpyi);B3SplineControlPoint(mu,U,tmz,tmpzi);for(i=0;i3*mu-17;i+)for(j=0;jmv-3;j+)tmxj=tmpxji;tmyj=tmpyji;tmzj=tmpzji;B3SplineControlPoint(mv,V,tmx,uvxi);B3SplineControlPoint(mv,V,tmy,uvyi);B3SplineControlPoint(mv,V,tmz,uvzi);for(k=0;k3*mv-17;k+)ptik0=uvxik;ptik1=uvyik;ptik2=uvzik;/计算双二次（2x2）B样条曲面所有控制多边形顶点，并保存在pt中/mu,mv分别为节点矢量U,V的最大下标值void BS2FaceControlPoint(int mu,float U,int mv,float V,float px,float py,float pz,float pt1001003)int i,j,k,dp;float tmx50,tmy50,tmz50;float tmpx50100,tmpy50100,tmpz50100;float uvx100100,uvy100100,uvz100100;for(i=0;imv-2;i+)dp=i*(mu-2);for(j=dp;jmu-2+dp;j+)tmxj-dp=pxj;tmyj-dp=pyj;tmzj-dp=pzj;B2SplineControlPoint(mu,U,tmx,tmpxi);B2SplineControlPoint(mu,U,tmy,tmpyi);B2SplineControlPoint(mu,U,tmz,tmpzi);for(i=0;i2*mu-7;i+)for(j=0;jmv-2;j+)tmxj=tmpxji;tmyj=tmpyji;tmzj=tmpzji;B2SplineControlPoint(mv,V,tmx,uvxi);B2SplineControlPoint(mv,V,tmy,uvyi);B2SplineControlPoint(mv,V,tmz,uvzi);for(k=0;k2*mv-7;k+)ptik0=uvxik;ptik1=uvyik;ptik2=uvzik;/-/设置网格数void SetGridCount(int dt,int tu,int tmk)int i,tm;tm=tu%dt;for(i=0;idt-1;i+) tmki=(tu-tm)/dt;tmkdt-1=tmk0+tm;/-/计算双三次(3x3次)或双二次(2x2次)B样条曲面上所有的点并保存在bs中/nu,mv分别为节点矢量U,V的最大下标/uk,vk分别为B样条曲面(u,v)方向上的网格数/p为曲面的次数void BSplineFace(int p,int nu,float U,int uk,int mv,float V,int vk, float px,float py,float pz,float bs1611613)int udk20,vdk20,i,j,k,l,hu,sv,du,dv;float tp1001003,td1611613;float tmx44,tmy44,tmz44;du=nu-2*p;dv=mv-2*p; SetGridCount(du,uk,udk); SetGridCount(dv,vk,vdk);if(p=3) BS3FaceControlPoint(nu,U,mv,V,px,py,pz,tp);if(p=2) BS2FaceControlPoint(nu,U,mv,V,px,py,pz,tp);for(i=0;idv;i+) for(k=0;kdu;k+) for(j=i*p;jp+1+i*p;j+) for(l=k*p;lp+1+k*p;l+) tmxj-i*pl-k*p=tplj0; tmyj-i*pl-k*p=tplj1; tmzj-i*pl-k*p=tplj2; BezierFacePoint(p,udkk,vdki,tmx,tmy,tmz,td); for(sv=i*vdk0;sv=vdki+i*vdk0;sv+) for(hu=k*udk0;hu=udkk+k*udk0;hu+) bssvhu0=tdsv-i*vdk0hu-k*udk00; bssvhu1=tdsv-i*vdk0hu-k*udk01; bssvhu2=tdsv-i*vdk0hu-k*udk02; /-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* Nurbs 样条曲线曲面部分 -*-*-*-*-*-*-*-*-*-*-*-*-*-*-/计算Nurbs曲线上的点（u所对应的所有点）保存在Poi中/n=m-p-1/p为曲线的次数void NurbsPoint(int n,int p,float U,float P,float W,float Poi)float N100,tmp,tmw; int i,j; for(i=p+1;i=i-p;j-)tmp+=Nj*Pj*Wj;tmw+=Nj*Wj;Poii-p=tmp/tmw;/计算Nurbs曲线的1阶导矢（u所对应的所有点）保存在Der中/n=m-p-1/p为曲线的次数void NurbsDer(int n,int p,float U,float P,float W,float Der)float N100,CP100,NW100,tmp,tw; int i,j; NurbsPoint(n,p,U,P,W,CP);BSplinePoint(n,p,U,W,NW);for(i=p+1;i=i-p;j-)tmp+=Nj*Pj*Wj;tw+=Nj*Wj;Deri-p=(tmp-tw*CPi-p)/NWi-p;/计算3次Nurbs曲线上的所有控制多边形保存在CP中/m为节点矢量U的最大下标void Nurbs3ControlPoint(int m,float U,float P,float W,float CP)int n,k,i,cp,p;float Poi100,Der100,add;p=3;n=m-p-1;NurbsPoint(n,p,U,P,W,Poi);NurbsDer(n,p,U,P,W,Der);cp=(n-p)*3+p;for(i=0;i2;i+)CPi=Pi;CPcp-i=Pn-i;for(i=3;icp-1;i+=3) k=(int)i/3;add=Derk/p;CPi=Poik;CPi-1=CPi-add;CPi+1=CPi+add;/计算2次Nurbs曲线上的所有控制多边形保存在CP中/m为节点矢量U的最大下标void Nurbs2ControlPoint(int m,float U,float P,float W,float CP)int n,k,tm,i,cp,p;float Poi100;p=2;n=m-p-1;NurbsPoint(n,p,U,P,W,Poi);cp=(n-p)*2+p;for(i=0;i2;i+)CPi=Pi;CPcp=Pn;tm=2;for(i=2;icp-1;i+=2) k=(int)i/2;CPi=Poik;CPi+1=Ptm;tm+;/绘制3次Nurbs样条曲线/m为节点矢量U的最大下标void Nurbs3L(int m,float U,float px,float py,float pz,float W)float pcx100,pcy100,pcz100,drx4,dry4,drz4;float pcw100,drw4;int i,j,tmcp;Nurbs3ControlPoint(m,U,px,W,pcx);Nurbs3ControlPoint(m,U,py,W,pcy);Nurbs3ControlPoint(m,U,pz,W,pcz);B3SplineControlPoint(m,U,W,pcw);tmcp=m-7;for(i=0;i=tmcp;i+)for(j=i*3;ji*3+4;j+)drxj-i*3=pcxj;dryj-i*3=pcyj;drzj-i*3=pczj;drwj-i*3=pcwj;NBezier(3,drx,dry,drz,drw,20);/绘制2次Nurbs样条曲线/m为节点矢量U的最大下标void Nurbs2L(int m,float U,float px,float py,float pz,float W)float pcx100,pcy100,pcz100,drx3,dry3,drz3;float pcw100,drw3;int i,j,tmcp;Nurbs2ControlPoint(m,U,px,W,pcx);Nurbs2ControlPoint(m,U,py,W,pcy);Nurbs2ControlPoint(m,U,pz,W,pcz);B2SplineControlPoint(m,U,W,pcw);tmcp=m-5;for(i=0;i=tmcp;i+)for(j=i*2;ji*2+3;j+)drxj-i*2=pcxj;dryj-i*2=pcyj;drzj-i*2=pczj;drwj-i*2=pcwj;NBezier(2,drx,dry,drz,drw,20);/计算双三次（3x3）Nurbs样条曲面所有控制多边形顶点，并保存在pt中/mu,mv分别为节点矢量U,V的最大下标值void Nurbs3FControlPoint(int mu,float U,int mv,float V,float px,float py,float pz,float W,float pt1001004)int i,j,k,dp;float tmx50,tmy50,tmz50,tmw50;float tmpx50100,tmpy50100,tmpz50100,tmpw50100;float uvx100100,uvy100100,uvz100100,uvw100100;for(i=0;imv-3;i+)dp=i*(mu-3);for(j=dp;jmu-3+dp;j+)tmxj-dp=pxj;tmyj-dp=pyj;tmzj-dp=pzj;tmwj-dp=Wj;Nurbs3ControlPoint(mu,U,tmx,tmw,tmpxi);Nurbs3ControlPoint(mu,U,tmy,tmw,tmpyi);Nurbs3ControlPoint(mu,U,tmz,tmw,tmpzi);B3SplineControlPoint(mu,U,tmw,tmpwi);for(i=0;i3*mu-17;i+)for(j=0;jmv-3;j+)tmxj=tmpxji;tmyj=tmpyji;tmzj=tmpzji;tmwj=tmpwji;Nurbs3ControlPoint(mv,V,tmx,tmw,uvxi);Nurbs3ControlPoint(mv,V,tmy,tmw,uvyi);Nurbs3ControlPoint(mv,V,tmz,tmw,uvzi);B3SplineControlPoint(mv,V,tmw,uvwi);for(k=0;k3*mv-17;k+)ptik0=uvxik;ptik1=uvyik;ptik2=uvzik;ptik3=uvwik;/计算双二次（2x2）Nurbs样条曲面所有控制多边形顶点，并保存在pt中/mu,mv分别为节点矢量U,V的最大下标值void Nurbs2FControlPoint(int mu,float U,int mv,float V,float px,float py,float pz,float W,float pt1001004)int i,j,k,dp;float tmx50,tmy50,tmz50,tmw50;float tmpx50100,tmpy50100,tmpz50100,tmpw50100;float uvx100100,uvy100100,uvz100100,uvw100100;for(i=0;imv-2;i+)dp=i*(mu-2);for(j=dp;jmu-2+dp;j+)tmxj-dp=pxj;tmyj-dp=pyj;tmzj-dp=pzj;tmwj-dp=Wj;Nurbs2ControlPoint(mu,U,tmx,tmw,tmpxi);Nurbs2ControlPoint(mu,U,tmy,tmw,tmpyi);Nurbs2ControlPoint(mu,U,tmz,tmw,tmpzi);B2SplineControlPoint(mu,U,tmw,tmpwi);for(i=0;i2*mu-7;i+)for(j=0;jmv-2;j+)tmxj=tmpxji;tmyj=tmpyji;tmzj=tmpzji;tmwj=tmpwji;Nurbs2ControlPoint(mv,V,tmx,tmw,uvxi);Nurbs2ControlPoint(mv,V,tm