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

学习小结：前面学习了Bezier曲线，B样条基函数和B样条曲线的一些基础知识。掌握关键问题是一条B样条曲线间的多段曲线的光滑连接。因为现在是用多段Bezier曲线来描绘一条B样条曲线，所以问题变为两段Bezier曲线间光滑连接。两段Bezier曲线段（3次）B1和B2光滑连接的条件：（1）.要求B1和B2有共同的连接点，即G0连续。（2）.要求B1和B2在连接点处有成比例的一阶导数，即G1连续。由端点处的一阶导数)(3)0(2),(3)1(10123QQBPPB,为实现G1连续，则有：)1(1)0(2BB即：2301PPQQ这也表明，1032),(,QQPP三点共线。如下图表示了一条3次B样条曲线的所有控制多边形：(P1)P2P3P4(P11)(P12)P5P10P0P6P9P7P8图5.3次B样条曲线和所有控制多边形图5中，P0至P6为原始3次B样条曲线控制多边形顶点，P0至P12是计算后最终形成B样条曲线控制多边形顶点。图6.双二次（2x2）B样条曲面6.B样条曲线曲面和NURBS曲线曲面的C语言实现算法源程序#ifndef_mynurbs_h#ifndef_MYNURBS_H#includegl\gl.h#includemath.h//*-*-*-*-*-*-*-*-*-*-*-*-*-*-*B样条基函数计算部分*-*-*-*-*-*-*-*-*-*-*-*-*-*//确定参数u所在的节点区间下标//n=m-p-1//m为节点矢量U[]的最大下标//p为B样条函数次数intFindSource(intn,intp,floatu,floatU[]){intlow,high,mid;if(u==U[n+1])//特殊情况returnn;//进行二分搜索low=p;high=n+1;mid=(int)(low+high)/2;while(uU[mid]||uU[mid]){if(uU[mid])high=mid;elselow=mid;mid=(int)(low+high)/2;if(u=U[mid]&&uU[mid+1])//找到u所在的节点区间的下标break;//退出二分搜索}returnmid;//返回区间下标}//计算所有非零B样条基函数并返回其值//i为参数u所在的节点区间下标voidBasisFunction(inti,intp,floatu,floatU[],floatN[]){intj,di,dp,k;floattul,tur,left,right;floattmpN[50][50];for(k=0;k=p;k++){dp=0;for(di=i+p-k;di=i-k;di--){if(u=U[di]&&uU[di+1])tmpN[di][0]=1;elsetmpN[di][0]=0;dp+=1;for(j=1;jdp;j++){tul=U[di+j]-U[di];tur=U[di+j+1]-U[di+1];if(tul!=0)left=(u-U[di])/tul;elseleft=0;if(tur!=0)right=(U[di+j+1]-u)/tur;elseright=0;tmpN[di][j]=left*tmpN[di][j-1]+right*tmpN[di+1][j-1];}}N[i-k]=tmpN[i-k][p];}}//-----------------------------------------------------------------------//计算基函数的1阶导数并保存在NP[]中//i为参数u所在的节点区间下标//p为B样条函数次数P2voidDerBasisFunc(inti,intp,floatu,floatU[],floatNP[]){intj,di,dp,k;floattul,tur,left,right,saved,dl,dr;floattmpN[50][50];for(k=0;k=p;k++){dp=0;for(di=i+p-k;di=i-k;di--){if(u=U[di]&&uU[di+1])tmpN[di][0]=1;elsetmpN[di][0]=0;dp+=1;for(j=1;jdp;j++){tul=U[di+j]-U[di];tur=U[di+j+1]-U[di+1];if(tul!=0)left=(u-U[di])/tul,dl=1/tul;elseleft=0,dl=0;if(tur!=0)right=(U[di+j+1]-u)/tur,dr=1/tur;elseright=0,dr=0;tmpN[di][j]=(left*tmpN[di][j-1]+right*tmpN[di+1][j-1]);saved=p*(dl*tmpN[di][j-1]-dr*tmpN[di+1][j-1])/(p+p-1);}}NP[i-k]=saved;}}//*-*-*-*-*-*-*-*-*-*-*-*-*-*Bezier曲线曲面部分*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*//计算参数u的p次基函数值并存在BC[]中voidBernsteinFunc(intp,doublet,floatBC[]){for(inti=0;i=p;i++){if(i==0)BC[0]=(float)pow(1-t,p);if(i==p)BC[p]=(float)pow(t,p);if(i0&&ip)BC[i]=p*(float)pow(t,i)*(float)pow(1-t,p-i);}}//获取p次Bezier曲线上的lines个点的值voidBezierPoint(intp,floatpx[],floatpy[],floatpz[],intlines,floattmp[][3]){floatBC[20];inti,j;for(j=0;j=lines;j++){doublet=j/(float)lines;BernsteinFunc(p,t,BC);tmp[j][0]=tmp[j][1]=tmp[j][2]=0;for(i=0;ip+1;i++){tmp[j][0]+=BC[i]*px[i];tmp[j][1]+=BC[i]*py[i];tmp[j][2]+=BC[i]*pz[i];}}}//获取p次有理Bezier曲线上的lines个点的值voidNBezierPoint(intp,floatpx[],floatpy[],floatpz[],floatpw[],intlines,floattmp[][4]){floatx,y,z,w,BC[20];inti,j;for(j=0;j=lines;j++){doublet=j/(float)lines;BernsteinFunc(p,t,BC);x=y=z=w=0;for(i=0;ip+1;i++){x+=BC[i]*px[i]*pw[i];y+=BC[i]*py[i]*pw[i];z+=BC[i]*pz[i]*pw[i];w+=BC[i]*pw[i];}tmp[j][0]=x/w;tmp[j][1]=y/w;tmp[j][2]=z/w;tmp[j][3]=w;}}//-----------------------------------------------------------------------------------//绘制p次的Bezier曲线voidBezier(intp,floatpx[],floatpy[],floatpz[],intlines){floatpt[100][3];intj;BezierPoint(p,px,py,pz,lines,pt);for(j=1;j=lines;j++){glBegin(GL_LINES);glVertex3f(pt[j-1][0],pt[j-1][1],pt[j-1][2]);glVertex3f(pt[j][0],pt[j][1],pt[j][2]);glEnd();}}//------------------------------------------------------------------------------//绘制p次的有理Bezier曲线voidNBezier(intp,floatpx[],floatpy[],floatpz[],floatw[],intlines){floatpt[100][4];intj;NBezierPoint(p,px,py,pz,w,lines,pt);for(j=1;j=lines;j++){glBegin(GL_LINES);glVertex3f(pt[j-1][0],pt[j-1][1],pt[j-1][2]);glVertex3f(pt[j][0],pt[j][1],pt[j][2]);glEnd();}}//---------------------------------------------------------------------------------//计算双p次Bezier曲面上所有的点并保存在Pt[][][]中//u和v分别为曲面(u,v)方向上的网格数voidBezierFacePoint(intp,intu,intv,floatpx[][4],floatpy[][4],floatpz[][4],floatpt[161][161][3]){floaturx[11][161],ury[11][161],urz[11][161];floattx[11],ty[11],tz[11],tmp[161][3];inti,j,k;for(j=0;jp+1;j++){for(i=0;ip+1;i++){tx[i]=px[i][j];ty[i]=py[i][j];tz[i]=pz[i][j];}BezierPoint(p,tx,ty,tz,v,tmp);for(k=0;k=v;k++){urx[j][k]=tmp[k][0];ury[j][k]=tmp[k][1];urz[j][k]=tmp[k][2];}}for(i=0;i=v;i++){for(k=0;kp+1;k++){tx[k]=urx[k][i];ty[k]=ury[k][i];tz[k]=urz[k][i];}BezierPoint(p,tx,ty,tz,u,tmp);for(j=0;j=u;j++){pt[i][j][0]=tmp[j][0];pt[i][j][1]=tmp[j][1];pt[i][j][2]=tmp[j][2];}}}//--------------------------------------------------------------------------------//计算双p次有理Bezier曲面上所有的点并保存在Pt[][][]中//u和v分别为曲面(u,v)方向上的网格数voidNuBezierFacePoint(intp,intu,intv,floatpx[][4],floatpy[][4],floatpz[][4],floatw[][4],floatpt[161][161][3]){floaturx[11][161],ury[11][161],urz[11][161],urw[11][161];floattx[11],ty[11],tz[11],tw[11],tmp[161][4];inti,j,k;for(j=0;jp+1;j++){for(i=0;ip+1;i++){tx[i]=px[i][j];ty[i]=py[i][j];tz[i]=pz[i][j];tw[i]=w[i][j];}NBezierPoint(p,tx,ty,tz,tw,v,tmp);for(k=0;k=v;k++){urx[j][k]=tmp[k][0];ury[j][k]=tmp[k][