B樣條曲線曲面和NURBS曲線曲面C語言算法源程序(總27頁

1、 學(xué)習(xí)小結(jié)：前面學(xué)習(xí)了Bezier曲線，B樣條基函數(shù)和B樣條曲線的一些基礎(chǔ)知識。掌握關(guān)鍵問題是一條B樣條曲線間的多段曲線的光滑連接。因?yàn)楝F(xiàn)在是用多段Bezier曲線來描繪一條B樣條曲線，所以問題變?yōu)閮啥蜝ezier曲線間光滑連接。兩段Bezier曲線段（3次）B1和B2光滑連接的條件： （1）.要求B1和B2有共同的連接點(diǎn)，即G0連續(xù)。 （2）.要求B1和B2在連接點(diǎn)處有成比例的一階導(dǎo)數(shù)，即G1連續(xù)。由端點(diǎn)處的一階導(dǎo)數(shù),為實(shí)現(xiàn)G1連續(xù)，則有： 即： 這也表明，三點(diǎn)共線。如下圖表示了一條3次B樣條曲線的所有控制多邊形： (P1) P2 P3 P4 (P11) (P12) P5 P10 P0 P6

2、 P9 P7 P8 圖5.3次B樣條曲線和所有控制多邊形 圖5中，P0至P6為原始3次B樣條曲線控制多邊形頂點(diǎn)，P0至P12是計(jì)算后最終形成B樣條曲線控制多邊形頂點(diǎn)。 圖6.雙二次（2x2）B樣條曲面 6.B樣條曲線曲面和NURBS曲線曲面的C語言實(shí)現(xiàn)算法源程序#ifndef _mynurbs_h#ifndef _MYNURBS_H#include glgl.h#include math.h/*-*-*-*-*-*-*-*-*-*-*-*-*-*-* B樣條基函數(shù)計(jì)算部分 *-*-*-*-*-*-*-*-*-*-*-*-*-*/確定參數(shù)u所在的節(jié)點(diǎn)區(qū)間下標(biāo)/n=m-p-1 /m為節(jié)點(diǎn)矢量U的最

3、大下標(biāo)/p為B樣條函數(shù)次數(shù)int FindSource(int n,int p,float u,float U)int low,high,mid;if(u=Un+1)/特殊情況return n;/進(jìn)行二分搜索low=p;high=n+1;mid=(int)(low+high)/2;while(uUmid)if(u=Umid&uUmid+1)/找到u所在的節(jié)點(diǎn)區(qū)間的下標(biāo) break; /退出二分搜索return mid; /返回區(qū)間下標(biāo)/計(jì)算所有非零B樣條基函數(shù)并返回其值/i為參數(shù)u所在的節(jié)點(diǎn)區(qū)間下標(biāo)void BasisFunction(int i,int p,float u,float U,f

4、loat 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) t

5、mpNdi0=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=sa

6、ved;/*-*-*-*-*-*-*-*-*-*-*-*-*-* Bezier曲線曲面部分 *-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*/計(jì)算參數(shù)u的p次基函數(shù)值并存在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個點(diǎn)的值void BezierPoint(int p,float px,float py,float pz,int lines,float tmp3)floa

7、t 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個點(diǎn)的值void NBezierPoint(int p,float px,float py,float pz,float pw,int lines,float tmp4)float x,y,z,w,BC20;int

8、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; BezierP

9、oint(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

10、(GL_LINES); glVertex3f(ptj-10,ptj-11,ptj-12); glVertex3f(ptj0,ptj1,ptj2); glEnd();/-/計(jì)算雙p次Bezier曲面上所有的點(diǎn)并保存在Pt中/u和v分別為曲面(u,v)方向上的網(wǎng)格數(shù)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+

11、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; /-/計(jì)算雙p次有

12、理Bezier曲面上所有的點(diǎn)并保存在Pt中/u和v分別為曲面(u,v)方向上的網(wǎng)格數(shù)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; NB

13、ezierPoint(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; /-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-

14、* B樣條曲線曲面部分 -*-*-*-*-*-*-*-*-*-*-*-*-*-*-/計(jì)算樣條曲線的1階導(dǎo)矢（u所對應(yīng)的所有點(diǎn)）保存在Der中/n=m-p-1/p為曲線的次數(shù)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;/計(jì)算曲線上的點(diǎn)（u所對應(yīng)的所有點(diǎn)）保存在Poi中/n=m-p-1/p為曲線的次數(shù)void BSplinePoint(int n,int p,float U,float P,float P

15、oi)float N100,tmp; int i,j; for(i=p+1;i=i-p;j-)tmp+=Nj*Pj;Poii-p=tmp;/計(jì)算3次樣條曲線上的所有控制多邊形保存在CP中/m為節(jié)點(diǎn)矢量U的最大下標(biāo)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;CP

16、cp-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;/計(jì)算2次樣條曲線上的所有控制多邊形保存在CP中/m為節(jié)點(diǎn)矢量U的最大下標(biāo)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;f

17、or(i=2;icp-1;i+=2) k=(int)i/2;CPi=Poik;CPi+1=Ptm;tm+;/繪制3次B樣條曲線/m為節(jié)點(diǎn)矢量U的最大下標(biāo)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.

18、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為節(jié)點(diǎn)矢量U的最大下標(biāo)void BSpline2L(

19、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

20、,drz,20);/計(jì)算雙三次（3x3）B樣條曲面所有控制多邊形頂點(diǎn)，并保存在pt中/mu,mv分別為節(jié)點(diǎn)矢量U,V的最大下標(biāo)值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);

21、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

22、,V,tmy,uvyi);B3SplineControlPoint(mv,V,tmz,uvzi);for(k=0;k3*mv-17;k+)ptik0=uvxik;ptik1=uvyik;ptik2=uvzik;/計(jì)算雙二次（2x2）B樣條曲面所有控制多邊形頂點(diǎn)，并保存在pt中/mu,mv分別為節(jié)點(diǎn)矢量U,V的最大下標(biāo)值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

23、 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+)

24、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;/-/設(shè)置網(wǎng)格數(shù)void SetGridCount(int dt,int tu,int tmk)int i,tm;tm=tu%dt;for(i=0;idt-1;i+) tmki=(tu-tm)/dt;t

25、mkdt-1=tmk0+tm;/-/計(jì)算雙三次(3x3次)或雙二次(2x2次)B樣條曲面上所有的點(diǎn)并保存在bs中/nu,mv分別為節(jié)點(diǎn)矢量U,V的最大下標(biāo)/uk,vk分別為B樣條曲面(u,v)方向上的網(wǎng)格數(shù)/p為曲面的次數(shù)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,tmy4

26、4,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;

27、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 樣條曲線曲面部分 -*-*-*-*-*-*-*-*-*-*-*-*

28、-*-*-/計(jì)算Nurbs曲線上的點(diǎn)（u所對應(yīng)的所有點(diǎn)）保存在Poi中/n=m-p-1/p為曲線的次數(shù)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;/計(jì)算Nurbs曲線的1階導(dǎo)矢（u所對應(yīng)的所有點(diǎn)）保存在Der中/n=m-p-1/p為曲線的次數(shù)void NurbsDer(int n,int p,float U,float P,float W,f

29、loat 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;/計(jì)算3次Nurbs曲線上的所有控制多邊形保存在CP中/m為節(jié)點(diǎn)矢量U的最大下標(biāo)void Nurbs3ControlPoint(int m,float U,float P,float W,float CP)int n,k,i,cp,p;float Poi100

30、,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;/計(jì)算2次Nurbs曲線上的所有控制多邊形保存在CP中/m為節(jié)點(diǎn)矢量U的最大下標(biāo)void Nurbs2ControlPoint(int m,float U,float P,float W,float C

31、P)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為節(jié)點(diǎn)矢量U的最大下標(biāo)void Nurbs3L(int m,float U,float px,float py,float pz,float W)float pcx100,pcy100,pcz100,drx4,dry4,dr

32、z4;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);/繪

33、制2次Nurbs樣條曲線/m為節(jié)點(diǎn)矢量U的最大下標(biāo)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;

34、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);/計(jì)算雙三次（3x3）Nurbs樣條曲面所有控制多邊形頂點(diǎn)，并保存在pt中/mu,mv分別為節(jié)點(diǎn)矢量U,V的最大下標(biāo)值void Nurbs3FControlPoint(int mu,float U,int mv,float V,float px,float py,float pz,float W,float pt1001004)int i,j,k,d

35、p;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);N

36、urbs3ControlPoint(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,

37、tmw,uvwi);for(k=0;k3*mv-17;k+)ptik0=uvxik;ptik1=uvyik;ptik2=uvzik;ptik3=uvwik;/計(jì)算雙二次（2x2）Nurbs樣條曲面所有控制多邊形頂點(diǎn)，并保存在pt中/mu,mv分別為節(jié)點(diǎn)矢量U,V的最大下標(biāo)值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 tmpx5010

38、0,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,tmy,