B样条曲线曲面和NURBS曲线曲1.doc

B样条曲线曲面和NURBS曲线曲面 (学习笔记和上机练习) 非均匀有理B样条，通常简称为NURBS(Non-Uniform Rational B-Splines)。 NURBS是非有理B样条、有理以及非有理Bzier曲线曲面的推广。一、 要对B样条曲线曲面和NURBS曲线曲面有所了解应先了解B样条基函数B样条基函数的定义和性质令是一个实数列，即，i=0,1,m-1。其中，称为节点，称为节点矢量，用表示第i个次（+1阶）B样条基函数，其定义为 若值为1，其他值为0 (2)由(2)式可知：（1）是一个阶梯函数，它在半开区间外都为零；（2）当0时，是两个-1次基函数的线性组合；（3）计算一组基函数时需要事先指定节点矢量和次数；（4）（2）式中可能出现0/0，我们规定0/0=0；（5）是定义在整个实数轴上的分段多项式函数，只在区间上有意义；（6）半开区间称为第i个节点区间，长度可以为零，因相邻节点可以相同；B样条基函数的一些重要性质：1 如果，则=0。 2 对于所有的和，有0.3 对于任意的节点区间，当时，。4 在节点区间内部，是无限次可微的。在节点处是次连续可微的，其中是节点的重复度。5 除的情况外，严格的达到一次最大值。B样条基函数的导数 基函数的求导公式为 （4）根据求导法则，对基函数 求导得到一般的求导公式： （5）另一个计算B样条基函数各阶导数的公式（参考Butt76）：， , (6)B样条基函数的一些重要计算a 确定节点()在节点矢量()的区间下标(据区间)/确定参数u所在的节点区间下标/n=m-p-1 /m为节点矢量U的最大下标/p为B样条函数次数int Buspan(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 计算节点参数u的所有非零B样条基函数值(据定义(2)式)/计算所有非零B样条基函数并返回其值,float N/i为参数u所在的节点区间下标void AllBNipu(int i,int p,float u,float U,float N) float tm1,tm2,TN100100; float left100100,right100100; for(int ti=0;ti=i+p;ti+) for(int j=0;j=Uti&u=i-p;t-) for(int k=1;k=i-p;r-) Nr=leftrp*TNrp-1+rightrp*TNr+1p-1; c 计算B样条基函数的导数(据(6)式)/计算基函数的1阶导数并保存在NP中/i为参数u所在的节点区间下标/p为B样条函数次数P2void BNipip(int i,int p,float u,float U,float NP)float tm1,tm2,TN100100; float left100100,right100100; for(int ti=0;ti=i+p+1;ti+) for(int j=0;j=Uti&u=i-p;t-) for(int k=1;k=i-p;l-) NPl=p*(u-Ul)*TNlp-1/(Ul+p-Ul)- (Ul+p+1-u)*TNl+1p-1/(Ul+p+1-Ul+1)/(p-1);二、B样条曲线曲面B样条曲线的定义次B样条曲线的定义为 ， aub (1)其中是控制点，是定义在非周期（并且非均匀）节点矢量 ，(a,a)和(b,b)为(+1)个节点，(包含m+1个节点)上的次B样条基函数。一般情况下，a=0,b=1。如下定义一条3次B样条曲线：U11=0.0,0.0,0.0,0.0,0.25,0.5,0.75,1.0,1.0,1.0,1.0;下图为定义在上节点矢量的一条3次B样条曲线(P0,P1,P6为控制多边形)。 P1 P2 P6 P3 P5 P4 图1、用三次Bezier曲线形式表示分段四次多项式3次B样条曲线B样条曲线的性质（1） 如果U=0.0,0.0,1.0,1.0;那么C(u)是Bezier曲线。 如下图2所示 P1 P2 P0 P3 图2、定义在U8=0.0,0.0,0.0,0.0,1.0,1.0,1.0,1.0;上的3次B样条曲线（2） C(u)是分段多项式曲线，因是分段多项式函数。次数p、控制点个数n+1和节点个数m+1应满足以下条件 m=n+p+1（3） 端点插值性：C(0)=P0 ，C(1)=Pn 。（4） 仿射不变性：对B样条曲线进行仿射变换，得到的仍为B样条曲线。（5） 强凸包性：曲线C(u)包含在它的控制多边形的凸包内（见图1和图2）。（6） 局部修改性：移动Pi只改变C(u)在区间上的形状。原因为，则=0。如下图所示定义在U10=0.0,0.0,0.0,0.0,0.34,0.7,1.0,1.0,1.0,1.0;上的三次曲线 P1 P2 P41 P5 P0 P3 P4 图3、移动P4到P41 （7） 控制多边形是对B样条曲线的一个分段性逼近。（8） 沿着曲线C(u)从u=0移动到u=1，基函数只与对应的u起作用。（9） C(u)的连续性和可微性（10） 变差减少性（11） 重控制顶点的利用B样条曲线的导矢 （7）其中表示的 阶导矢。B样条曲线在端点处的一阶导矢： ， （8）B样条曲线在端点处的二阶导矢： ， （9）对于Bezier曲线上面（8）、（9）式可简化为下面对应的表达式： B样条曲线的一些计算/计算B样条曲线上的点并返回其值/n=m-p-1 /m为节点矢量U的最大下标/p为B样条函数次数float BCurveP(int n,int p,float U,float P,float u)int span;float N100,Cp;span=Buspan(n,p,u,U);AllBNipu(span,p,u,U,N);Cp=0.0;for(int i=span-1;i=span-p;i-)Cp+=Ni*Pi;return Cp;/计算曲线上的1阶导矢并返回1阶的导矢/n为参数u所在的区间下标/p为次数float BCurveCN(int n,int p,float u,float U,float P) float CN,NP100; BNipip(n,p,u,U,NP);CN=0.0;for(int i=n-1;i=n-p;i-) CN+=NPi*Pi;return CN;B样条曲面的定义B样条曲面由两个方向的控制点网格、两个节点矢量和单变量B样条基函数的乘积来定义，方程如下 （10）节点矢量为 和 U中含有r+1个节点，V中有s+1个节点，据B样条曲线性质（2）有 r=n+p+1 , s=m+q+1B样条曲面的性质（1）曲面插值于控制网格的四个角点：S(0,0)=P0,0 ,S(1,0)=Pn,0 ,S(0,1)=P0,m ,S(1,1)=Pn,m 得到 （2） 仿射不变性（3） 局部修改性（4） 凸包性B样条曲面的偏导矢 （11）角点偏导矢公式（u,v）=(0,0): B样条曲面的一些计算/计算固定参数值(u,v),对应的B样条曲面上的点并返回其值/n,m分别为节点矢量U和V的最大下标/p,q分别为B样条曲面u,v两个方向上的次数float BSfacePoint(int n,int p,float U,int m,int q,float V, float P5,float u,float v)int uspan,vspan,vd,ud;float Nu100,Nv100,tmp100; float S;uspan=Buspan(n,p,u,U);AllBNipu(uspan,p,u,U,Nu);vspan=Buspan(m,q,v,V); 图4、双三次B样条曲面AllBNipu(vspan,q,v,V,Nv);for(vd=vspan;vd=vspan-q;vd-)tmpvd=0.0;for(ud=uspan;ud=uspan-p;ud-) tmpvd+=Nuud*Pudvd;S=0.0;for(vd=vspan;vd=vspan-q;vd-)S+=Nvvd*tmpvd;return S;下面以双二次（22次）B样条曲面为例说明曲面生成的算法#include windows.h#include glgl.h#include glglu.h#include glglut.h#include math.hfloat mx=0.0,32.0,32.0,-32.0,-32.0,0.0, 0.0,20.0,20.0,-20.0,-20.0,0.0, 0.0,40.0,40.0,-40.0,-40.0,0.0, 0.0,55.0,55.0,-55.0,-55.0,0.0, 0.0,50.0,50.0,-50.0,-50.0,0.0, 0.0,30.0,30.0,-30.0,-30.0,0.0;float my=70.0,70.0,70.0,70.0,70.0,70.0, 50.0,50.0,50.0,50.0,50.0,50.0, 40.0,40.0,40.0,40.0,40.0,40.0, 30.0,30.0,30.0,30.0,30.0,30.0, 10.0,10.0,10.0,10.0,10.0,10.0, -20.0,-20.0,-20.0,-20.0,-20.0,-20.0;float mz=32.0,32.0,-32.0,-32.0,32.0,32.0, 20.0,20.0,-20.0,-20.0,20.0,20.0, 40.0,40.0,-40.0,-40.0,40.0,40.0, 55.0,55.0,-55.0,-55.0,55.0,55.0, 50.0,50.0,-50.0,-50.0,50.0,50.0, 30.0,30.0,-30.0,-30.0,30.0,30.0;float MU9=0.0,0.0,0.0,0.25,0.5,0.75,1.0,1.0,1.0;float MV9=0.0,0.0,0.0,0.25,0.5,0.75,1.0,1.0,1.0; typedef struct BSPLINEPOINTSfloat x;float y;float z;float w;int u;int v;BSPLINEPOINTS;BSPLINEPOINTS suv161161;void lightm()GLfloat lamb4=0.0f,0.0f,0.0f,1.0f;GLfloat ldif4=0.1f,0.15f,0.13f,1.0f;GLfloat lspe4=0.0f,0.0f,0.0f,0.0f;GLfloat lpos4=0.0f,200.0f,200.0f,1.0f;GLfloat lamb14=0.2f,0.2f,0.2f,1.0f;GLfloat ldif14=0.15f,0.13f,0.2f,1.0f;GLfloat lspe14=0.0f,0.0f,0.0f,0.0f;GLfloat lpos14=200.0f,0.0f,0.0f,1.0f;GLfloat lamb24=0.2f,0.20f,0.2f,1.0f;GLfloat ldif24=0.00f,0.13f,0.15f,1.0f;GLfloat lspe24=0.0f,0.0f,0.0f,0.0f;GLfloat lpos24=0.0f,-200.0f,-200.0f,1.0f;GLfloat lamb34=0.0f,0.0f,0.0f,1.0f;GLfloat ldif34=0.1f,0.15f,0.1f,1.0f;GLfloat lspe34=0.0f,0.0f,0.0f,0.0f;GLfloat lpos34=-0.0f,80.0f,0.0f,1.0f;GLfloat mamb4=0.2f,0.2f,0.2f,1.0f;GLfloat mdif4=0.2f,0.2f,0.2f,1.0f;GLfloat mspe4=0.0f,0.0f,0.0f,0.0f;GLfloat memi4=0.1f,0.1f,0.1f,1.0f;GLfloat mshininess=0.0f;glLightfv(GL_LIGHT1,GL_AMBIENT,lamb);glLightfv(GL_LIGHT1,GL_DIFFUSE,ldif);glLightfv(GL_LIGHT1,GL_SPECULAR,lspe);glLightfv(GL_LIGHT1,GL_POSITION,lpos);glLightfv(GL_LIGHT2,GL_AMBIENT,lamb1);glLightfv(GL_LIGHT2,GL_DIFFUSE,ldif1);glLightfv(GL_LIGHT2,GL_SPECULAR,lspe1);glLightfv(GL_LIGHT2,GL_POSITION,lpos1);glLightfv(GL_LIGHT3,GL_AMBIENT,lamb2);glLightfv(GL_LIGHT3,GL_DIFFUSE,ldif2);glLightfv(GL_LIGHT3,GL_SPECULAR,lspe2);glLightfv(GL_LIGHT3,GL_POSITION,lpos2);glLightfv(GL_LIGHT4,GL_AMBIENT,lamb3);glLightfv(GL_LIGHT4,GL_DIFFUSE,ldif3);glLightfv(GL_LIGHT4,GL_SPECULAR,lspe3);glLightfv(GL_LIGHT4,GL_POSITION,lpos3);glMaterialfv(GL_FRONT_AND_BACK,GL_AMBIENT,mamb);glMaterialfv(GL_FRONT_AND_BACK,GL_DIFFUSE,mdif);glMaterialfv(GL_FRONT_AND_BACK,GL_SPECULAR,mspe);glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,memi);glMaterialf(GL_FRONT_AND_BACK,GL_SHININESS,mshininess);/计算多边形的外法线返回值tmNvoid getN(float x3,float y3,float z3,float tmN3)float p1,p2,p3,q1,q2,q3;float nx,ny,nz;p1=x1-x0;p2=y1-y0;p3=z1-z0;q1=x2-x1;q2=y2-y1;q3=z2-z1;nx=p2*q3-q2*p3;ny=q1*p3-p1*q3;nz=p1*q2-p2*q1;tmN0=nx;tmN1=ny;tmN2=nz;/确定参数u所在的节点区间下标/n=m-p-1 m为节点矢量U的最大下标/p为B样条函数次数int Buspan(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 AllBNipu(int i,int p,float u,float U,float N) float tm1,tm2,TN100100; float left100100,right100100; for(int ti=0;ti=i+p;ti+) for(int j=0;j=Uti&u=i-p;t-) for(int k=1;k=i-p;r-) Nr=leftrp*TNrp-1+rightrp*TNr+1p-1; /计算B样条曲线上的点并返回其值/n=m-p-1 m为节点矢量U的最大下标/p为B样条函数次数float BCurveP(int n,int p,float U,float P,float u)int span;float N100,Cp;span=Buspan(n,p,u,U);AllBNipu(span,p,u,U,N);Cp=0.0;for(int i=span-1;i=span-p;i-)Cp+=Ni*Pi;return Cp;/计算参数u的p次基函数值并存在BC中void BeBC(int p,double t,float BC)for(int i=0;i0&ip)BCi=p*(float)pow(t,i)*(float)pow(1-t,p-i);/计算二次Bezier曲面上所有的点并保存在Pt中/u和v分别为曲面(u,v)方向上的网格数void Bezier2Point(int u,int v,float px33,float py33,float pz33, BSPLINEPOINTS Pt161)float urx3161,ury3161,urz3161,BC10;int i,j,k,p=2;float dut,dvt,x,y,z;double t; if(u=0)u=1;if(v=0)v=1;dut=(float)u;dvt=(float)v;for(j=0;j3;j+) for(k=0;k=v;k+) t=k/dvt; BeBC(p,t,BC); x=0.0;y=0.0;z=0.0; for(i=0;ip+1;i+) x+=BCi*pxij; y+=BCi*pyij; z+=BCi*pzij; urxjk=x; uryjk=y; urzjk=z; for(i=0;i=v;i+) for(j=0;j=u;j+) t=j/dut; BeBC(p,t,BC); x=0.0;y=0.0;z=0.0; for(k=0;kp+1;k+) x+=BCk*urxki;y+=BCk*uryki;z+=BCk*urzki; Ptij.x=x; Ptij.y=y; Ptij.z=z; /计算2次B样条曲线所有控制多边形顶点保存在Pt/m为节点矢量U的最大下标void B2ControlPoints(int m,float U,float px,float py, float pz,BSPLINEPOINTS Pt)int i,j,n,sp,p=2;float tx20,ty20,tz20; n=m-p-1;sp=n-p;sp=n+sp;for(i=p+1;i=n;i+) txi-p=BCurveP(n,p,U,px,Ui); tyi-p=BCurveP(n,p,U,py,Ui); tzi-p=BCurveP(n,p,U,pz,Ui); for(i=0;i2;i+)Pti.x=pxi;Pti.y=pyi;Pti.z=pzi;Ptsp.x=pxn;Ptsp.y=pyn;Ptsp.z=pzn;i=2;for(j=2;jsp-1;j+=2)Ptj.x=txj/2;Ptj.y=tyj/2;Ptj.z=tzj/2;Ptj+1.x=pxi;Ptj+1.y=pyi;Ptj+1.z=pzi;i+;/计算双二次(2x2次)B样条曲面所有控制多边形顶点保存在Pt/nu和mv分别为节点矢量U和V的最大下标void B2FControlPoint(int nu,float U,int mv,float V,float px,float py,float pz,BSPLINEPOINTS Pt100) int i,j,p=2; float tpx30,tpy30,tpz30; BSPLINEPOINTS td50100; for(i=0;imv-p;i+) for(j=i*(nu-p);jnu-p+i*(nu-p);j+) tpxj-i*(nu-p)=pxj; tpyj-i*(nu-p)=pyj; tpzj-i*(nu-p)=pzj; B2ControlPoints(nu,U,tpx,tpy,tpz,tdi); for(i=0;i(nu-5)*p+3;i+) for(j=0;j=160)uk=160;if(vk=160)vk=160;bs00.u=uk;bs00.v=vk;du=nu-2*2;tm=uk%du;for(i=0;idu-1;i+) udki=(uk-tm)/du;udkdu-1=(uk-tm)/du+tm; dv=mv-2*2;tm=vk%dv;for(j=0;jdv-1;j+) vdkj=(vk-tm)/dv;vdkdv-1=(vk-tm)/dv+tm;B2FControlPoint(nu,U,mv,V,px,py,pz,tp);for(i=0;idv;i+) for(k=0;kdu;k+) for(j=i*2;j3+i*2;j+) for(l=k*2;l3+k*2;l+) tmxj-i*2l-k*2=tplj.x; tmyj-i*2l-k*2=tplj.y; tmzj-i*2l-k*2=tplj.z; Bezier2Point(udkk,vdki,tmx,tmy,tmz,td); for(sv=i*vdk0;sv=vdki+i*vdk0;sv+) for(hu=k*udk0;hu=udkk+k*udk0;hu+) bssvhu.x=tdsv-i*vdk0hu-k*udk0.x; bssvhu.y=tdsv-i*vdk0hu-k*udk0.y; bssvhu.z=tdsv-i*vdk0hu-k*udk0.z; /显示B样条曲面void ShowBSplineCurveFace(BSPLINEPOINTS bs161,int fill)int i,j;float x3,y3,z3,tmn3;for(i=0;i=bs00.v;i+) for(j=0;j=bs00.u;j+) if(fill!=0)x0=bsij.x;x1=bsij+1.x;x2=bsi+1j+1.x;y0=bsij.y;y1=bsij+1.y;y2=bsi+1j+1.y;z0=bsij.z;z1=bsij+1.z;z2=bsi+1j+1.z; getN(x,y,z,tmn); /计算多边形的外法线返回值tmNglEnable(GL_NORMALIZE);glBegin(GL_QUADS);glNormal3f(tmn0,tmn1,tmn2);if(jbs00.u)glVertex3f(bsij.x,bsij.y,bsij.z);glVertex3f(bsij+1.x,bsij+1.y,bsij+1.z); if(ibs00.v)glVertex3f(bsi+1j+1.x,bsi+1j+1.y,bsi+1j+1.z);glVertex3f(bsi+1j.x,bsi+1j.y,bsi+1j.z);glEnd();glDisable(GL_NORMALIZE); elseglBegin(GL_LINES);if(jbs00.u) glVertex3f(bsij.x,bsij.y,bsij.z); glVertex3f(bsij+1.x,bsij+1.y,bsij+1.z);if(ibs00.v)glVertex3f(bsij.x,bsij.y,bsij.z); glVertex3f(bsi+1j.x,bsi+1j.y,bsi+1j.z); glEnd(); void RenderScene() glTranslatef(-0.0,-20.0,0.0); glRotatef(-25.0,1.0,0.0,0.0); Bspline2Fac