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

1、學習小結(jié)：前面學習了Bezier曲線，B樣條基函數(shù)和B樣條曲線的一些基礎(chǔ)知識。掌握關(guān)鍵問題是一條B樣條曲線間的多段曲線的光滑連接。因為現(xiàn)在是用多段Bezier曲線來描繪一條B樣條曲線，所以問題變?yōu)閮啥蜝ezier曲線間光滑連接。兩段Bezier曲線段(3次)B1和B2光滑連接的條件：(1).要求B1和B2有共同的連接點，即G5連續(xù)。(2).要求B1和B2在連接點處有成比例的一階導(dǎo)數(shù)，即G連續(xù)。由端點處的一階導(dǎo)數(shù)1B1(1) =3(P3 P2), B2(0) =3(QQo),為實現(xiàn)G連續(xù)，則有：B2(0)=B1(1)即：Q1Qo=RP2這也表明，P2, P3(Qo), Q三點共線。如下圖表示了一

2、條3次B樣條曲線的所有控制多邊形:?；(P1) P -已P3P/4(P/ P火和P0P6P3P7圖5.3次B樣條曲線和所有控制多邊形圖5中，P0至P6為原始3次B樣條曲線控制多邊形頂點,圖6.雙二次(2x2) B樣條曲面6.B樣條曲線曲面和NURBJffi線曲面的C語言實現(xiàn)算法源程序 #ifndef _mynurbs_h#ifndef _MYNURBS_H#include glgl.h#include math.hP0至P12是計算后最終形成B樣條曲線控制多邊形頂點。11)俗如P12)七*-*-*-*-*-*-*-*-*-*-*-*/確定參數(shù)u所在的節(jié)點區(qū)間下標 /n=m-p-1/m為節(jié)點矢量

3、U的最大下標/p為B樣條函數(shù)次數(shù)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所在的節(jié)點區(qū)間的下標break; /退出二分搜索return mid; /返回區(qū)間下標/計算所有非零B樣條基函數(shù)并返回其值/i為參數(shù)u所在的節(jié)點區(qū)間下標void BasisFunction(int i,int p,float u,flo

4、at 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;elsetmpNdi0=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&

5、;uUdi+1)tmpNdi0=1;elsetmpNdi0=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; elseleft=0,dl=0;if(tur!=0)right=(Udi+j+1-u)/tur,dr=1/tur;elseright=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;/

6、*-*-*-*-*-*-*-*-*-*-*-*-*-*Bezier曲 線 曲 面 部 分*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*/計算參數(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個點的值void BezierPoint(int p,float px,float py,float pz,int lines,float tmp3)f

7、loat 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,

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

9、,pz,lines,pt);for(j=1;j=lines;j+)(gl Begin(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+) gl Begin(GL_LINES);g

10、lVertex3f(ptj-10,ptj-11,ptj-12);glVertex3f(ptj0,ptj1,ptj2);glEnd(); /-計算雙p次Bezier曲面上所有的點并保存在Pt皿中/u和v分別為曲面(u,v)方向上的網(wǎng)格數(shù)void BezierFacePoint(int p,int u,int v,float px4,float py4,floatpz4,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

11、+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分別為曲

12、面(u,v)方向上的網(wǎng)格數(shù)void NuBezierFacePoint(int p,int u,int v,float px4,float py4,floatpz4,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);fo

13、r(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,t z,tw,u,tmp);for(j=0;j=u;j+)(ptij0=tmpj0;ptij1=tmpj1;ptij2=tmpj2;/-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*B樣 條 曲 線 曲 面音K分-*-*-*-*-*-*-*-*-*-*-*-*-*-*-/計

14、算樣條曲線的1階導(dǎo)矢(u所對應(yīng)的所有點)保存在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;計算曲線上的點(u所對應(yīng)的所有點)保存在Poi中/n=m-p-1/p為曲線的次數(shù)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-

15、)tmp+=Nj*Pj;Poii-p=tmp;/計算3次樣條曲線上的所有控制多邊形保存在CP中/m為節(jié)點矢量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

16、=Derk/p;CPi=Poik;CPi-1=CPi-add;CPi+1=CPi+add;/計算2次樣條曲線上的所有控制多邊形保存在CP中/m為節(jié)點矢量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

17、+1=Ptm;tm+;/繪制3次B樣條曲線/m為節(jié)點矢量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+)(glB egi

18、n(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é)點矢量U的最大下標void BSpline2L(int m,float U,float px,float py,float p

19、z)(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樣條曲面所有控制多邊形頂點，

20、并保存在pt中/mu,mv分別為節(jié)點矢量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=px

21、j;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);B3SplineControlP

22、oint(mv,V,tmz,uvzi); for(k=0;k3*mv-17;k+) ( ptik0=uvxik; ptik1=uvyik;ptik2=uvzik; 計算雙二次(2x2) B樣條曲面所有控制多邊形頂點，并保存在pt中/mu,mv分別為節(jié)點矢量U,V的最大下標值void BS2FaceControlPoint(int mu,float U,int mv,float V,floatpx,floatpy,floatpz,floatpt1001003)(inti,j,k,dp;floattmx50,tmy50,tmz50;floattmpx50100,tmpy50100,tmpz5010

23、0; 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);B2Spli