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

B 样条曲线曲面和NURBS 曲线曲面
(学习笔记和上机练习)
非均匀有理B 样条，通常简称为NURBS(Non-Uniform Rational B-Splines)。

NURBS 是非有理B 样条、有理以及非有理B ézier 曲线曲面的推广。

一、要对B 样条曲线曲面和NURBS 曲线曲面有所了解应先了解B 样条基函数 B 样条基函数的定义和性质
令{}m u u u U ,...,,10=是一个实数列，即i u ≤1+i u ，i=0,1,…,m-1。

其中，i u 称为节点，
U 称为节点矢量，用)(,u N p i 表示第i 个p 次（p +1阶）B 样条基函数，其定义为
⎩
⎨
⎧=,0,
1)(0,u N i 若i u ≤u ＜1+i u 值为1，其他值为0 )()()(1,11
111,,u N u u u u u N u u u u u N p i i p i p i p i i p i i
p i -++++++-+--+--= (2)
由(2)式可知：
（1）)(0,u N i 是一个阶梯函数，它在半开区间),[1+∈i i u u u 外都为零； （2）当p ＞0时，)(,u N p i 是两个p -1次基函数的线性组合；
（3）计算一组基函数时需要事先指定节点矢量U 和次数p ； （4）（2）式中可能出现0/0，我们规定0/0=0；
（5）)(,u N p i 是定义在整个实数轴上的分段多项式函数，只在区间][,0m u u 上有意义； （6）半开区间),[1+i i u u 称为第i 个节点区间，长度可以为零，因相邻节点可以相同； B 样条基函数的一些重要性质：
1 如果),[1++∉p i i u u u ，则)(,u N p i =0。

2 对于所有的p i ,和u ，有)(,u N p i ≥0.
3 对于任意的节点区间),[1+i i u u ，当),[1+∈i i u u u 时，
∑-==i
p
i j p
j u N
1)(,。

4 在节点区间内部，)(,u N p i 是无限次可微的。

在节点处)(,u N p i 是k p -次连续可微
的，其中k 是节点的重复度。

5 除0=p 的情况外，)(,u N p i 严格的达到一次最大值。

B 样条基函数的导数
基函数的求导公式为
)()(1,11
11,,u N u u p
u N u u p N p i i p i p i i p i p i -++++-+-+-=
' （4）
根据求导法则g f g f fg '+'=')(，对基函数 )()()(1,11
111,,u N u u u u u N u u u u u N p i i p i p i p i i p i i
p i -++++++-+--+--=
求导得到一般的求导公式： )(
1
1)
1(1
,1)1(1
,)
(,+++--++----
-=i p i k p i i
p i k p i k p
i u u N u u N p N
（5）
另一个计算B 样条基函数各阶导数的公式（参考[Butt76]）：
)()
(1,11
11)(1,)(,k p i i p i p i k p i i p i i k p i N u u u u N u u u u k p p N -++++++-+------=
， ,1,0=k …,1-p (6) B 样条基函数的一些重要计算
a – 确定节点(u )在节点矢量(U )的区间下标(据),[1+∈i i u u u 区间) //确定参数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==U[n+1])//特殊情况 return n; //进行二分搜索 low=p; high=n+1;
mid=(int)(low+high)/2; while(u<U[mid]||u>U[mid]) {
if(u<U[mid]) high=mid; else
low=mid;
mid=(int)(low+high)/2;
if(u>=U[mid]&&u<U[mid+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,TN[100][100];
float left[100][100],right[100][100];
for(int ti=0;ti<=i+p;ti++)
for(int j=0;j<=p;j++)
{
tm1=U[ti+j]-U[ti];
tm2=U[ti+j+1]-U[ti+1];
if(tm1!=0) left[ti][j]=(u-U[ti])/tm1;
else
left[ti][j]=0.0;
if(tm2!=0)
right[ti][j]=(U[ti+j+1]-u)/tm2;
else
right[ti][j]=0.0;
if(u>=U[ti]&&u<U[ti+1])
TN[ti][0]=1.0;
else
TN[ti][0]=0.0;
}
N[0]=TN[i][0];
for(int t=i+p-1;t>=i-p;t--)
for(int k=1;k<p;k++)
TN[t][k]=left[t][k]*TN[t][k-1]+right[t][k]*TN[t+1][k-1];
for(int r=i;r>=i-p;r--)
N[r]=left[r][p]*TN[r][p-1]+right[r][p]*TN[r+1][p-1];
}
c –计算B样条基函数的导数(据(6)式)
//计算基函数的1阶导数并保存在NP[]中
//i为参数u所在的节点区间下标
//p为B样条函数次数P>2
void BNipip(int i,int p,float u,float U[],float NP[])
{
float tm1,tm2,TN[100][100];
float left[100][100],right[100][100];
for(int ti=0;ti<=i+p+1;ti++)
for(int j=0;j<=p;j++) {
tm1=U[ti+j]-U[ti];
tm2=U[ti+j+1]-U[ti+1];
if(tm1!=0) left[ti][j]=(u-U[ti])/tm1; else
left[ti][j]=0.0; if(tm2!=0)
right[ti][j]=(U[ti+j+1]-u)/tm2; else
right[ti][j]=0.0; if(u>=U[ti]&&u<U[ti+1]) TN[ti][0]=1.0; else
TN[ti][0]=0.0; }
for(int t=i+p-1;t>=i-p;t--) for(int k=1;k<p;k++)
TN[t][k]=left[t][k]*TN[t][k-1]+right[t][k]*TN[t+1][k-1]; for(int l=i;l>=i-p;l--)
NP[l]=p*((u-U[l])*TN[l][p-1]/(U[l+p]-U[l])-
(U[l+p+1]-u)*TN[l+1][p-1]/(U[l+p+1]-U[l+1]))/(p-1); }
二、B 样条曲线曲面 B 样条曲线的定义
p 次B 样条曲线的定义为
∑==
n
i i p
i P u N
u C 0
,)()( ， a ≤u ≤b (1)
其中{}i P 是控制点，{}
)(,u N p i 是定义在非周期（并且非均匀）节点矢量
{}
b b u u a a U p m p ,...,,...,,...11--+=，(a,…a)和(b,…b)为(p +1)个节点，(包含m+1个节点)上的p 次B 样条基函数。

一般情况下，a=0,b=1。

如下定义一条3次B 样条曲线： U[11]={0.0,0.0,0.0,0.0,0.25,0.5,0.75,1.0,1.0,1.0,1.0};
下图为定义在上节点矢量的一条3次B 样条曲线(P 0,P 1,…,P 6为控制多边形)。

P 1 P 2
P 6
0P P 3 P 5 P 4
图1、用三次Bezier 曲线形式表示分段四次多项式3次B 样条曲线 B 样条曲线的性质
（1） 如果U[]={0.0,…,1.0,…,1.0};那么C(u)是Bezier 曲线。

如下图2所示
P 1 P 2
P 0 P 3
图2、定义在U[8]={0.0,0.0,0.0,0.0,1.0,1.0,1.0,1.0};上的3次B 样条曲线 （2） C(u)是分段多项式曲线，因)(,u N p i 是分段多项式函数。

次数p 、控制点个数n+1和节点个数m+1应满足以下条件
m=n+p+1 （3） 端点插值性：C(0)=P 0 ，C(1)=P n 。

（4） 仿射不变性：对B 样条曲线进行仿射变换，得到的仍为B 样条曲线。

（5） 强凸包性：曲线C(u)包含在它的控制多边形的凸包内（见图1和图2）。

（6） 局部修改性：移动P i 只改变C(u)在区间),[1++p i i u u 上的形状。

原因为),[1++∉p i i u u u ，则)(,u N p i =0。

如下图所示定义在U[10]={0.0,0.0,0.0,0.0,0.34,0.7,1.0,1.0,1.0,1.0};上的三次曲线
P 1 P 2 P 41
P 5 P 0 P 3 P 4
图3、移动P 4到P 41
（7） 控制多边形是对B 样条曲线的一个分段性逼近。

（8） 沿着曲线C(u)从u=0移动到u=1，基函数)(,u N p i 只与对应的u 起作用。

（9） C(u)的连续性和可微性 （10） 变差减少性
（11） 重控制顶点的利用 B 样条曲线的导矢 ∑==n
i i k p i k P u N u C
)(,)
()()( （7）
其中)()
(u C k 表示)(u C 的k 阶导矢。

B 样条曲线在端点处的一阶导矢： )()0(011
P P u p C p -=
'+ ，)(1)1(11
-----=
'n n p m P P u p
C （8）
B 样条曲线在端点处的二阶导矢：
⎥⎥⎦
⎤⎢⎢⎣⎡++--=
+++++++222112110
1)
2()()1()0(p p p p p p p u P u u P u u u P u p p C
，
⎥⎥⎦
⎤
⎢⎢⎣⎡-+--------=
----------------222112111
)
2(1)1)(1()2(11)
1()1(p m n p m p m n p m p m p m n p m u P u u P u u u P u p p C
（9）
对于Bezier 曲线上面（8）、（9）式可简化为下面对应的表达式：
)()0(01P P n C -=' )()1(1--='n n P P n C
)2)(1()0(012P P P n n C +--='' )2)(1()1(21--+--=''n n n P P P n n C
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 N[100],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+=N[i]*P[i]; return Cp; }
//计算曲线上的1阶导矢并返回1阶的导矢 //n 为参数u 所在的区间下标 //p 为次数
float BCurveCN(int n,int p,float u,float U[],float P[]) {
float CN,NP[100];
BNipip(n,p,u,U,NP); CN=0.0;
for(int i=n-1;i>=n-p;i--) CN+=NP[i]*P[i]; return CN; }
B 样条曲面的定义
B 样条曲面由两个方向的控制点网格、两个节点矢量和单变量B 样条基函数的乘积来定义，方程如下
∑∑===n i m
j j i q j p
i P v N u N
v u S 00
,,,)()(),( （10）
节点矢量为
}1,...,1,,...,0,...,0{1--=p r u U 和 }1,...,1,,...,0,...,0{1--=q s v V
U 中含有r+1个节点，V 中有s+1个节点，据B 样条曲线性质（2）有 r=n+p+1 , s=m+q+1 B 样条曲面的性质
（1）曲面插值于控制网格的四个角点：S(0,0)=P 0,0 ,S(1,0)=P n,0 ,S(0,1)=P 0,m ,S(1,1)=P n,m 得到
1)1()1()1()0()0()1()0()0(,,,,0,0,,0,0====q m p n q m p q p n q p N N N N N N N N （2） 仿射不变性 （3） 局部修改性 （4） 凸包性
B 样条曲面的偏导矢
j i l q j n i m
j k p i l k l k P v N u N v u S v u ,)
(,00
)(,)()(),(∑∑==+=∂∂∂ （11）
角点偏导矢公式（u,v ）=(0,0): )()0,0(0,00,11
P P u p S p u -=
+ )()0,0(0,01,01
P P v q S q v -=
+
)()0,0(0,00,11,01,11
1P P P P v u pq
S q p uv +--=
++
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 P[][5],float u,float v) {
int uspan,vspan,vd,ud;
float Nu[100],Nv[100],tmp[100]; 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--)
{
tmp[vd]=0.0;
for(ud=uspan;ud>=uspan-p;ud--)
tmp[vd]+=Nu[ud]*P[ud][vd];
}
S=0.0;
for(vd=vspan;vd>=vspan-q;vd--)
S+=Nv[vd]*tmp[vd];
return S;
}
下面以双二次（2×2次）B样条曲面为例说明曲面生成的算法
#include "windows.h"
#include "gl\gl.h"
#include "gl\glu.h"
#include "gl\glut.h"
#include "math.h"
float 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,
,-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 MU[9]={0.0,0.0,0.0,0.25,0.5,0.75,1.0,1.0,1.0};
float MV[9]={0.0,0.0,0.0,0.25,0.5,0.75,1.0,1.0,1.0}; typedef struct BSPLINEPOINTS
{
float x;
float y;
float z;
float w;
int u;
int v;
}BSPLINEPOINTS;
BSPLINEPOINTS suv[161][161];
void lightm()
{
GLfloat lamb[4]={0.0f,0.0f,0.0f,1.0f};
GLfloat ldif[4]={0.1f,0.15f,0.13f,1.0f};
GLfloat lspe[4]={0.0f,0.0f,0.0f,0.0f};
GLfloat lpos[4]={0.0f,200.0f,200.0f,1.0f};
GLfloat lamb1[4]={0.2f,0.2f,0.2f,1.0f};
GLfloat ldif1[4]={0.15f,0.13f,0.2f,1.0f};
GLfloat lspe1[4]={0.0f,0.0f,0.0f,0.0f};
GLfloat lpos1[4]={200.0f,0.0f,0.0f,1.0f};
GLfloat lamb2[4]={0.2f,0.20f,0.2f,1.0f};
GLfloat ldif2[4]={0.00f,0.13f,0.15f,1.0f};
GLfloat lspe2[4]={0.0f,0.0f,0.0f,0.0f};
GLfloat lpos2[4]={0.0f,-200.0f,-200.0f,1.0f};
GLfloat lamb3[4]={0.0f,0.0f,0.0f,1.0f};
GLfloat ldif3[4]={0.1f,0.15f,0.1f,1.0f};
GLfloat lspe3[4]={0.0f,0.0f,0.0f,0.0f};
GLfloat lpos3[4]={-0.0f,80.0f,0.0f,1.0f};
GLfloat mamb[4]={0.2f,0.2f,0.2f,1.0f};
GLfloat mdif[4]={0.2f,0.2f,0.2f,1.0f};
GLfloat mspe[4]={0.0f,0.0f,0.0f,0.0f};
GLfloat memi[4]={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); }
//计算多边形的外法线返回值tmN[]
void getN(float x[3],float y[3],float z[3],float tmN[3]) {
float p1,p2,p3,q1,q2,q3;
float nx,ny,nz;
p1=x[1]-x[0];
p2=y[1]-y[0];
p3=z[1]-z[0];
q1=x[2]-x[1];
q2=y[2]-y[1];
q3=z[2]-z[1];
nx=p2*q3-q2*p3;
ny=q1*p3-p1*q3;
nz=p1*q2-p2*q1;
tmN[0]=nx;
tmN[1]=ny;
tmN[2]=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==U[n+1])//特殊情况
return n;
//进行二分搜索
low=p;
high=n+1;
mid=(int)(low+high)/2;
while(u<U[mid]||u>U[mid])
{
if(u<U[mid])
high=mid;
else
low=mid;
mid=(int)(low+high)/2;
if(u>=U[mid]&&u<U[mid+1])//找到u所在的节点区间的下标
break; //退出二分搜索
}
return mid; //返回区间下标
}
//计算所有非零B样条基函数并返回其值
//i为参数u所在的节点区间下标
void AllBNipu(int i,int p,float u,float U[],float N[])
{
float tm1,tm2,TN[100][100];
float left[100][100],right[100][100];
for(int ti=0;ti<=i+p;ti++)
for(int j=0;j<=p;j++)
{
tm1=U[ti+j]-U[ti];
tm2=U[ti+j+1]-U[ti+1];
if(tm1!=0) left[ti][j]=(u-U[ti])/tm1;
else
left[ti][j]=0.0;
if(tm2!=0)
right[ti][j]=(U[ti+j+1]-u)/tm2;
else
right[ti][j]=0.0;
if(u>=U[ti]&&u<U[ti+1])
TN[ti][0]=1.0;
else
TN[ti][0]=0.0;
}
N[0]=TN[i][0];
for(int t=i+p-1;t>=i-p;t--)
for(int k=1;k<p;k++)
TN[t][k]=left[t][k]*TN[t][k-1]+right[t][k]*TN[t+1][k-1];
for(int r=i;r>=i-p;r--)
N[r]=left[r][p]*TN[r][p-1]+right[r][p]*TN[r+1][p-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 N[100],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+=N[i]*P[i];
return Cp;
}
//计算参数u的p次基函数值并存在BC[]中
void BeBC(int p,double t,float BC[])
{
for(int i=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)
BC[i]=p*(float)pow(t,i)*(float)pow(1-t,p-i);
}
}
//计算二次Bezier曲面上所有的点并保存在Pt[][]中
//u和v分别为曲面(u,v)方向上的网格数
void Bezier2Point(int u,int v,float px[3][3],float py[3][3],float pz[3][3], BSPLINEPOINTS Pt[][161])
{
float urx[3][161],ury[3][161],urz[3][161],BC[10];
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;j<3;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;i<p+1;i++)
{
x+=BC[i]*px[i][j];
y+=BC[i]*py[i][j];
z+=BC[i]*pz[i][j];
}
urx[j][k]=x;
ury[j][k]=y;
urz[j][k]=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;k<p+1;k++)
{
x+=BC[k]*urx[k][i];
y+=BC[k]*ury[k][i];
z+=BC[k]*urz[k][i];
}
Pt[i][j].x=x;
Pt[i][j].y=y;
Pt[i][j].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 tx[20],ty[20],tz[20];
n=m-p-1;
sp=n-p;
sp=n+sp;
for(i=p+1;i<=n;i++)
{
tx[i-p]=BCurveP(n,p,U,px,U[i]);
ty[i-p]=BCurveP(n,p,U,py,U[i]);
tz[i-p]=BCurveP(n,p,U,pz,U[i]);
}
for(i=0;i<2;i++)
{
Pt[i].x=px[i];
Pt[i].y=py[i];
Pt[i].z=pz[i];
}
Pt[sp].x=px[n];
Pt[sp].y=py[n];
Pt[sp].z=pz[n];
i=2;
for(j=2;j<sp-1;j+=2)
{
Pt[j].x=tx[j/2];
Pt[j].y=ty[j/2];
Pt[j].z=tz[j/2];
Pt[j+1].x=px[i];
Pt[j+1].y=py[i];
Pt[j+1].z=pz[i];
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 Pt[][100])
{
int i,j,p=2;
float tpx[30],tpy[30],tpz[30];
BSPLINEPOINTS td[50][100];
for(i=0;i<mv-p;i++)
{
for(j=i*(nu-p);j<nu-p+i*(nu-p);j++)
{
tpx[j-i*(nu-p)]=px[j];
tpy[j-i*(nu-p)]=py[j];
tpz[j-i*(nu-p)]=pz[j];
}
B2ControlPoints(nu,U,tpx,tpy,tpz,td[i]);
}
for(i=0;i<(nu-5)*p+3;i++)
{
for(j=0;j<mv-p;j++)
{
tpx[j]=td[j][i].x;
tpy[j]=td[j][i].y;
tpz[j]=td[j][i].z;
}
B2ControlPoints(mv,V,tpx,tpy,tpz,Pt[i]);
}
}
//计算双二次(2x2次)B样条曲面上所有的点并保存在bs[][]中
//nu,mv分别为节点矢量U[],V[]的最大下标
//uk,vk分别为B样条曲面(u,v)方向上的网格数
void Bspline2Face(int nu,float U[],int uk,int mv,float V[],int vk,
float px[],float py[],float pz[],
BSPLINEPOINTS bs[][161])
{
int udk[20],vdk[20],i,j,k,l,hu,sv,du,dv,tm;
BSPLINEPOINTS tp[50][100],td[161][161];
float tmx[3][3],tmy[3][3],tmz[3][3];
if(uk>=160)
uk=160;
if(vk>=160)
vk=160;
bs[0][0].u=uk;
bs[0][0].v=vk;
du=nu-2*2;
tm=uk%du;
for(i=0;i<du-1;i++)
udk[i]=(uk-tm)/du;
udk[du-1]=(uk-tm)/du+tm;
dv=mv-2*2;
tm=vk%dv;
for(j=0;j<dv-1;j++)
vdk[j]=(vk-tm)/dv;
vdk[dv-1]=(vk-tm)/dv+tm;
B2FControlPoint(nu,U,mv,V,px,py,pz,tp);
for(i=0;i<dv;i++)
{
for(k=0;k<du;k++)
{
for(j=i*2;j<3+i*2;j++)
for(l=k*2;l<3+k*2;l++)
{
tmx[j-i*2][l-k*2]=tp[l][j].x;
tmy[j-i*2][l-k*2]=tp[l][j].y;
tmz[j-i*2][l-k*2]=tp[l][j].z;
}
Bezier2Point(udk[k],vdk[i],tmx,tmy,tmz,td);
for(sv=i*vdk[0];sv<=vdk[i]+i*vdk[0];sv++)
for(hu=k*udk[0];hu<=udk[k]+k*udk[0];hu++)
{
bs[sv][hu].x=td[sv-i*vdk[0]][hu-k*udk[0]].x;
bs[sv][hu].y=td[sv-i*vdk[0]][hu-k*udk[0]].y;
bs[sv][hu].z=td[sv-i*vdk[0]][hu-k*udk[0]].z;
}
}
}
}
//显示B样条曲面
void ShowBSplineCurveFace(BSPLINEPOINTS bs[][161],int fill) {
int i,j;
float x[3],y[3],z[3],tmn[3];
for(i=0;i<=bs[0][0].v;i++)
for(j=0;j<=bs[0][0].u;j++)
{
if(fill!=0)
{
x[0]=bs[i][j].x;
x[1]=bs[i][j+1].x;
x[2]=bs[i+1][j+1].x;
y[0]=bs[i][j].y;
y[1]=bs[i][j+1].y;
y[2]=bs[i+1][j+1].y;
z[0]=bs[i][j].z;
z[1]=bs[i][j+1].z;
z[2]=bs[i+1][j+1].z;
getN(x,y,z,tmn);//计算多边形的外法线返回值tmN[] glEnable(GL_NORMALIZE);
glBegin(GL_QUADS);
glNormal3f(tmn[0],tmn[1],tmn[2]);
if(j<bs[0][0].u)
{
glVertex3f(bs[i][j].x,bs[i][j].y,bs[i][j].z);
glVertex3f(bs[i][j+1].x,bs[i][j+1].y,bs[i][j+1].z);
}
if(i<bs[0][0].v)
{
glVertex3f(bs[i+1][j+1].x,bs[i+1][j+1].y,bs[i+1][j+1].z);
glVertex3f(bs[i+1][j].x,bs[i+1][j].y,bs[i+1][j].z);
}
glEnd();
glDisable(GL_NORMALIZE);
}
else{
glBegin(GL_LINES);
if(j<bs[0][0].u)
{
glVertex3f(bs[i][j].x,bs[i][j].y,bs[i][j].z);
glVertex3f(bs[i][j+1].x,bs[i][j+1].y,bs[i][j+1].z);
}
if(i<bs[0][0].v)
{
glVertex3f(bs[i][j].x,bs[i][j].y,bs[i][j].z);
glVertex3f(bs[i+1][j].x,bs[i+1][j].y,bs[i+1][j].z);
}
glEnd();
}
}
}
void RenderScene()
{
glTranslatef(-0.0,-20.0,0.0);
glRotatef(-25.0,1.0,0.0,0.0);
Bspline2Face(8,MU,160,8,MV,160,mx,my,mz,suv);
glEnable(GL_BLEND);
glBlendFunc(GL_ONE,GL_SRC_ALPHA_SATURATE);
glColor3f(1.0f,0.8f,0.72f);
ShowBSplineCurveFace(suv,1);
glutSwapBuffers();
}
void SetupRC()
{
glClearColor(1.0f,1.0f,1.0f,0.0f);
lightm();
glEnable(GL_DEPTH_TEST);
glEnable(GL_LIGHTING);
glLightModelf(GL_LIGHT_MODEL_LOCAL_VIEWER,1.0);
glEnable(GL_LIGHT1);
glEnable(GL_LIGHT2);
glEnable(GL_LIGHT3);
glEnable(GL_LIGHT4);
glEnable(GL_COLOR_MATERIAL);
}
void WindowSize(GLsizei w,GLsizei h)
{
GLfloat aspectRatio;
GLfloat tmb=100.0;
if(h==0)
h=1;
glViewport(0,0,w,h);
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
aspectRatio=(GLfloat)w/(GLfloat)h;
if(w<=h)
glOrtho(-tmb,tmb,-tmb/aspectRatio,tmb/aspectRatio,tmb,-tmb);
else
glOrtho(-tmb*aspectRatio,tmb*aspectRatio,-tmb,tmb,tmb,-tmb);
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
}
int main(int argc,char *argv[])
{
glutInit(&argc,argv);
glutInitDisplayMode(GLUT_DOUBLE|GLUT_RGB);
glutInitWindowSize(1024,768);
glutCreateWindow("bcurved");
glutDisplayFunc(RenderScene);
glutReshapeFunc(WindowSize);
SetupRC();
glutMainLoop();
return 0;
}
运行程序结果如下
图5、双二次B 样条曲面渲染图 图6、双二次B 样条曲面网格图
三、NURBS 曲线曲面
NURBS 曲线的定义
一条p 次NURBS 曲线定义为 ∑∑===n
i i
p i n i i i p i w u N P w u N u C 0,0
,)()()( ， a ≤u ≤b （3.1） 其中i P 是控制点，i w 是权因子，{})(,u N p i 是定义在非周期（并且非均匀）节点矢量 {}
b b u u a a U p m p ,...,,...,,...11--+=，(a,…a)和(b,…b)为(p +1)个节点，(包含m+1个节点)上的p 次B 样条基函数。

一般情况下，a=0,b=1，并且对所有的i w i ,＞0 。

NURBS 曲线的性质
（1） 端点 C(0)=P 0 ,C(1)=P n
（2） 仿射不变性
（3） 强凸包性
（4） 变差减少性
（5） 局部修改性：如下图所示定义在权因子w[7]={1.0,9.0,7.0,2.0,5.0,10.0,1.0}; 节点矢量U[11]={0.0,0.0,0.0,0.0,0.25,0.5,0.75,1.0,1.0,1.0,1.0};的一条3次
NURBS 曲线(P 0,…,P 6为控制多边形)。

P 1 P 2
P 6
P 0 P 3 P 5
P 4
图7、三次NURBS 曲线
NURBS 曲线的导矢 )()
()()()(0,u w u C u w P w u N u C n i i i p i ∑='-'=' （3.2）
NURBS 曲线在端点处的一阶导矢表达式 )()0(01011P P w u pw C p -='+ ， )()1()1(111------='n n n
p m n P P w u pw C （3.3） 有理Bezier 曲线在端点处的一阶导矢表达式 )()0(0101P P w w C -=' ， )()1(11---='n n n
n P P w w C
NURBS曲线的一些计算
//计算有理B样条曲线上的点返回其值
//n=m-p-1
//m为节点矢量U[]的最大下标
//p为有理B样条函数次数
float NurbsP(int n,int p,float U[],float P[],float W[],float u) {
int span;
float N[100],Np,nb,wb;
span=Buspan(n,p,u,U);
AllBNipu(span,p,u,U,N);
nb=0.0;wb=0.0;
for(int i=span-1;i>=span-p;i--)
{
nb+=N[i]*P[i]*W[i];
wb+=N[i]*W[i];
}
Np=nb/wb;
return Np;
}
//计算有理样条曲线上的1阶导矢返回其值
//n=m-p-1;
//m为节点矢量U[]的最大下标
//p为次数
float NurbsCN(int n,int p,float u,float U[],float P[],float W[]) {
int span;
float CN,NP[100],nb,wb,Np,Wp;
span=Buspan(n,p,u,U);
BNipip(span,p,u,U,NP);
Np=NurbsP(n,p,U,P,W,u);
Wp=BCurveP(n,p,U,W,u);
nb=0.0;wb=0.0;
for(int i=span-1;i>=span-p;i--)
{
nb+=NP[i]*P[i]*W[i];
wb+=NP[i]*W[i];
}
CN=(nb-wb*Np)/Wp;
return CN;
}
NURBS曲面的定义
∑∑∑∑=====n i m
j j
i q j p i n i m j j i j i q j p i w v N u N P w v N u N
v u S 00,,,00
,,,,)()()()(),( ， 0≤v u ,≤1 （3.4） 其中，j i P ,两个方向上的控制网格，j i w ,是权因子，)(),(,,v N u N q j p i 是分别定义在节点矢量V U ,上的非有理B 样条基函数，
节点矢量为
}1,...,1,,...,0,...,0{1--=p r u U 和 }1,...,1,,...,0,...,0{1--=q s v V
图8、双三次NURBS 曲面 图9、双三次NURBS 曲面。