B样条曲线曲面和NURBS曲线曲面C语言算法源程序.docx
《B样条曲线曲面和NURBS曲线曲面C语言算法源程序.docx》由会员分享，可在线阅读，更多相关《B样条曲线曲面和NURBS曲线曲面C语言算法源程序.docx（38页珍藏版）》请在冰豆网上搜索。
B样条曲线曲面和NURBS曲线曲面C语言算法源程序
学习小结：
前面学习了Bezier曲线，B样条基函数和B样条曲线的一些基础知识。
掌握关键问题是一条B样条曲线间的多段曲线的光滑连接。
因为现在是用多段Bezier曲线来描绘
一条B样条曲线，所以问题变为两段Bezier曲线间光滑连接。
两段Bezier曲线段（3次）
B1和B2光滑连接的条件：
（1）.要求B1和B2有共同的连接点，即&连续。
（2）.要求B1和B2在连接点处有成比例的一阶导数，即G连续。
由端点处的一阶导数
1
B1
（1）3（F3P2）,B2（0）3（QiQ。
），为实现G连续，则有:
B2（0）B1
（1）即:
Q1QoRP2
这也表明，P2,B（Qo）,Q三点共线。
如下图表示了一条
3次B样条曲线的所有控制多边形:
E（P1）
P2
（P
/P
用P0
10
P4
图5.3次B样条曲线和所有控制多边形
图5中，P0至P6为原始3次B样条曲线控制多边形顶点，P0至P12是计算后最终形成
样条曲线控制多边形顶点。
图6.双二次（2x2）B样条曲面
6.B样条曲线曲面和NURBS曲线曲面的C语言实现算法源程序#ifndef_mynurbs_h
#ifndefMYNURBSH
#include"gl\gl.h"
#include"math.h"
if（u==U[n+1]）//特殊情况returnn;
//进行二分搜索
low=p;
high=n+1;
mid=（int）（low+high）/2;while（uU[mid]）
{if（umid=（int）（low+high）/2;
if（u>=U[mid]&&u}
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]&&utmpN[di][0]=1;
else
tmpN[di][0]=0;
dp+=1;
for（j=1;j{
tul=U[di+j]-U[di];
tur=U[di+j+1]-U[di+1];
if（tul!
=0）
left=（u-U[di]）/tul;
else
left=0;
if（tur!
=0）
right=（U[di+j+1]-u）/tur;
else
right=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样条函数次数P>2
voidDerBasisFunc（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]&&utmpN[di][0]=1;
else
tmpN[di][0]=0;
dp+=1;
for（j=1;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;
else
left=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（i>0&&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;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个点的值
pw[],int
voidNBezierPoint（intp,floatpx[],floatpy[],floatpz[],floatlines,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;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）方向上的网格数
py[][4],float
voidBezierFacePoint（intp,intu,intv,floatpx[][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;j{
for（i=0;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;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）方向上的网格数
py[][4],float
voidNuBezierFacePoint（intp,intu,intv,floatpx[][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;j{
for（i=0;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][1];
urz[j][k]=tmp[k][2];
urw[j][k]=tmp[k][3];
}
}
for（i=0;i<=v;i++）
{
for（k=0;k{
tx[k]=urx[k][i];
ty[k]=ury[k][i];
tz[k]=urz[k][i];tw[k]=urw[k][i];
}
NBezierPoint（p,tx,ty,tz,tw,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];
}
}
}
//
********
*******
B样条曲线曲面部
**************
//计算样条曲线的1阶导矢（u所对应的所有点）保存在Der[]中//n=m-p-1
//p为曲线的次数
voidBSplineDer（intn,intp,floatU[],floatP[],floatDer[]）
{
floatN[100],tmp;
inti,j;
for（i=p+1;i<=n;i++）
{
DerBasisFunc（i,p,U[i],U,N）;
tmp=0;
for（j=i;j>=i-p;j--）tmp+=N[j]*P[j];
Der[i-p]=tmp;
//计算曲线上的点（u所对应的所有点）保存在Poi[]中//n=m-p-1
//p为曲线的次数
voidBSplinePoint（intn,intp,floatU[],floatP[],floatPoi[]）
{
floatN[100],tmp;
inti,j;
for（i=p+1;i<=n;i++）
{
BasisFunction（i,p,U[i],U,N）;
tmp=0;for（j=i;j>=i-p;j--）tmp+=N[j]*P[j];
Poi[i-p]=tmp;
}
}
//计算3次样条曲线上的所有控制多边形保存在CP[]中//m为节点矢量U[]的最大下标
voidB3SplineControlPoint（intm,floatU[],floatP[],floatCP[]）
{
intn,k,i,cp,p;
floatPoi[100],Der[100],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;i<2;i++）
{
CP[i]=P[i];CP[cp-i]=P[n-i];
}
for（i=3;i{
k=（int）i/3;add=Der[k]/p;CP[i]=Poi[k];
CP[i-1]=CP[i]-add;CP[i+1]=CP[i]+add;
}
}
//计算2次样条曲线上的所有控制多边形保存在CP[]中
//m为节点矢量U[]的最大下标
voidB2SplineControlPoint（intm,floatU[],floatP[],floatCP[]）
{
intn,k,tm,i,cp,p;
floatPoi[100];
p=2;
n=m-p-1;
BSplinePoint（n,p,U,P,Poi）;
cp=（n-p）*2+p;
for（i=0;i<2;i++）
CP[i]=P[i];
CP[cp]=P[n];
tm=2;
for（i=2;i{
k=（int）i/2;
CP[i]=Poi[k];
CP[i+1]=P[tm];
tm++;
}
}
//绘制3次B样条曲线
//m为节点矢量U[]的最大下标
voidBSpline3L（intm,floatU[],floatpx[],floatpy[],floatpz[]）
{
floatpcx[100],pcy[100],pcz[100],drx[4],dry[4],drz[4];
inti,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;i<3*m-17;i++）
{
glBegin（GL_LINES）;glVertex3f（pcx[i-1],pcy[i-1],pcz[i-1]）;glVertex3f（pcx[i],pcy[i],pcz[i]）;
glEnd（）;
}
glColor3f（1.0f,0.0f,0.0f）;*/
tmcp=m-7;
for（i=0;i<=tmcp;i++）
for（j=i*3;j{
drx[j-i*3]=pcx[j];
dry[j-i*3]=pcy[j];
drz[j-i*3]=pcz[j];
}
Bezier（3,drx,dry,drz,20）;
}
}
//绘制2次B样条曲线
//m为节点矢量U[]的最大下标
voidBSpline2L（intm,floatU[],floatpx[],floatpy[],floatpz[]）
{
floatpcx[100],pcy[100],pcz[100],drx[3],dry[3],drz[3];
inti,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;j{
drx[j-i*2]=pcx[j];
dry[j-i*2]=pcy[j];
drz[j-i*2]=pcz[j];
}
Bezier（2,drx,dry,drz,20）;
}
}
//计算双三次（3x3）B样条曲面所有控制多边形顶点，并保存在pt[][][]中
//mu,mv分别为节点矢量U[],V[]的最大下标值
px[],float
voidBS3FaceControlPoint（intmu,floatU[],intmv,floatV[],floatpy[],floatpz[],floatpt[100][100][3]）
{
inti,j,k,dp;
floattmx[50],tmy[50],tmz[50];
floattmpx[50][100],tmpy[50][100],tmpz[50][100];
floatuvx[100][100],uvy[100][100],uvz[100][100];
for（i=0;i{
dp=i*（mu-3）;
for（j=dp;j{
tmx[j-dp]=px[j];
tmy[j-dp]=py[j];
tmz[j-dp]=pz[j];
}
B3SplineControlPoint（mu,U,tmx,tmpx[i]）;
B3SplineControlPoint（mu,U,tmy,tmpy[i]）;
B3SplineControlPoint（mu,U,tmz,tmpz[i]）;
}
for（i=0;i<3*mu-17;i++）
{
for（j=0;j{
tmx[j]=tmpx[j][i];
tmy[j]=tmpy[j][i];
tmz[j]=tmpz[j][i];
}
B3SplineControlPoint（mv,V,tmx,uvx[i]）;
B3SplineControlPoint（mv,V,tmy,uvy[i]）;
B3SplineControlPoint（mv,V,tmz,uvz[i]）;for（k=0;k<3*mv-17;k++）
{
pt[i][k][0]=uvx[i][k];
pt[i][k][1]=uvy[i][k];
pt[i][k][2]=uvz[i][k];
}
}
}
//计算双二次（2x2）B样条曲面所有控制多边形顶点，并保存在pt[][][]中
//mu,mv分别为节点矢量U[],V[]的最大下标值
px[],float
voidBS2FaceControlPoint（intmu,floatU[],intmv,floatV[],floatpy[],floatpz[],floatpt[100][100][3]）
{
inti,j,k,dp;
floattmx[50],tmy[50],tmz[50];
floattmpx[50][100],tmpy[50][100],tmpz[50][100];
floatuvx[100][100],uvy[100][100],uvz[100][100];
for（i=0;i{
dp=i*（mu-2）;
for（j=dp;j{
tmx[j-dp]=px[j];
tmy[j-dp]=py[j];
tmz[j-dp]=pz[j];
}
B2SplineControlPoint（mu,U,tmx,tmpx[i]）;
B2SplineControlPoint（mu,U,tmy,tmpy[i]）;
B2SplineControlPoint（mu,U,tmz,tmpz[i]）;
}
for（i=0;i<2*mu-7;i++）
{
for（j=0;j{
tmx[j]=tmpx[j][i];
tmy[j]=tmpy[j][i];
tmz[j]=tmpz[j][i];
}
B2SplineControlPoint（mv,V,tmx,uvx[i]）;
B2SplineControlPoint（mv,V,tmy,uvy[i]）;
B2SplineControlPoint（mv,V,tmz,uvz[i]）;for（k=0;k<2*mv-7;k++）
{
pt[i][k][0]=uvx[i][k];
pt[i][k][1]=uvy[i][k];
pt[i][k][2]=uvz[i][k];
}
}
}
////设置网格数
voidSetGridCount（intdt,inttu,inttmk[]）
{
inti,tm;
tm=tu%dt;
for（i=0;itmk[i]=（tu-tm）/dt;
tmk[dt-1]=tmk[0]+tm;
}
//
//计算双三次（3x3次）或双二次（2x2次）B样条曲面上所有的点并保存在bs[][][]中
//nu,mv分别为节点矢量U[],V[]的最大下标
//uk,vk分别为B样条
升级会员 