下载文档原格式
四、B 样条曲线与曲面Bezier 曲线具有很多优越性，但有二点不足：1）特征多边形顶点数决定了它的阶次数，当n 较大时，不仅计算量增大，稳定性降低，且控制顶点对曲线的形状控制减弱；2）不具有局部性，即修改一控制点对曲线产生全局性影响。
1972年Gordon 等用B 样条基代替Bernstein 基函数，从而改进上述缺点。
Ｂ样条曲线的数学表达式为：∑=+⋅=nk n k ki n i u N Pu P 0,,)()(在上式中，0 ≤ u ≤ 1； i= 0, 1, 2, …, m 所以可以看出：Ｂ样条曲线是分段定义的。
如果给定 m+n+1 个顶点 Pi ( i=0, 1, 2,…, m+n)，则可定义 m+1 段 n 次的参数曲线。
在以上表达式中：N k,n (u) 为 n 次B 样条基函数，也称Ｂ样条分段混合函数。
其表达式为：∑-=+--+⋅⋅-=kn j nj n j n k j k n u C n u N 01,)()1(!1)(式中：0 ≤ u ≤1k = 0, 1, 2, …, n1．均匀B 样条曲线1一次均匀B 样条曲线的矩阵表示空间n+1个顶点i P （i = 0，1，…，n ）定义n 段一次（k ＝0,1，n=1）均匀B 样条曲线，即每相邻两个点可构造一曲线段P i （u ），其定义表达为：[]10 ;,...,1 0111 1)(1≤≤=⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡-=-u n i u u P i i i P P＝（1－u ）P i －1 ＋ u P i＝ N 0，1（u ）P i －1 ＋ N 1，1（u ）P i第i 段曲线端点位置矢量：i i i i P P P P ==-)1(,)0(1，且一次均匀B 样条曲线就是控制多边形。
2 二次均匀B 样条曲线的空间n+1个顶点的位置矢量i P （i=0，1，…，n ）定义n －1段二次（k ＝0,1,2， n=2）均匀B 样条曲线，每相邻三个点可构造一曲线段P i （u ）（i=1，…，n －1），其定义表达为：[]10 ;1,...,1 011022121 121)(112≤≤-=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--=+-u n i u u u P i i i i P P P＝ !21（1 － 2 u ＋ u 2）P i －1 ＋ !21（1 ＋ 2 u － 2u 2）P i ＋ !21u 2 P i ＋1＝ N 0，2（u ）P i －1 ＋ N 1，2（u ）P i ＋ N 2，2（u ）P i ＋1端点位置矢量：)(5.0)0(1i i i P P P +=-，)(5.0)1(1++=i i i P P P ，即曲线的起点和终点分别位于控制多边形P i-1P i 和P i P i+1的中点。
Bezier曲线、B样条曲线和NURBS曲线0.概述1. 贝塞尔曲线（Bezier Curve）：贝塞尔曲线由一组控制点和控制点上的权重组成。
贝塞尔曲线的阶数由控制点的数量决定，阶数为n的贝塞尔曲线需要n+1个控制点。
贝塞尔曲线具有局部控制的特性，即曲线上的一段由相邻的几个控制点决定，不受其他控制点的影响。
贝塞尔曲线的计算相对简单，但在变形过程中可能会出现形状扭曲的问题。
2. B样条（B-Spline）： B样条曲线是一种基于分段多项式的曲线表示方法。
与贝塞尔曲线不同，B样条曲线的每个控制点都有一个关联的基函数。
这些基函数决定了曲线上每一点的形状。
B样条曲线的阶数可以是任意的，较高阶的B样条曲线能够更灵活地描述复杂的曲线形状。
B样条曲线具有良好的局部控制性和平滑性，可以很好地避免贝塞尔曲线的形状扭曲问题。
3. NURBS曲线（Non-Uniform Rational B-Spline Curve）：NURBS曲线是对B样条曲线的扩展，它引入了有理权重的概念。
NURBS曲线的每个控制点都有一个关联的权重，这些权重可以调节曲线上各个点的影响程度。
NURBS曲线能够表示更复杂的曲线形状，如圆弧和椭圆等。
总的来说Bezier曲线中的每个控制点都会影响整个曲线的形状，而B样条中的控制点只会影响整个曲线的一部分，显然B样条提供了更多的灵活性；Bezier和B样条都是多项式参数曲线，不能表示一些基本的曲线，比如圆，所以引入了NURBS，即非均匀有理B样条来解决这个问题；贝塞尔曲线适用于简单的曲线形状设计，B样条曲线具有更好的局部控制和平滑性，适用于复杂曲线的建模而NURBS曲线在B样条的基础上引入了有理权重，可以更准确地描述各种曲线形状Bezier曲线是B样条的一个特例，而B样条又是NURBS的一个特例1.Bezier曲线1.1 贝塞尔曲线的历史:贝塞尔曲线于 1962 年，由法国工程师皮埃尔·贝济埃（PierreBézier）所广泛发表，他运用贝塞尔曲线来为汽车的主体进行设计,贝塞尔曲线最初由保尔·德·卡斯特里奥于1959年运用德卡斯特里奥算法开发，以稳定数值的方法求出贝塞尔曲线。
八、NURBS曲线（面）1、NURBS曲线（面）的提出Bezier曲线（面）B样条曲线（面）B样条方法在表示与设计自由型曲线、曲面形状时显示了强大的威力，然而在表示与设计初等曲线、曲面时时却遇到了麻烦B样条曲线（面）和Bezier曲线（面）不能精确表示出抛物线以外的二次曲线（面）。
如圆弧、圆而在工程应用特别是在工业设计当中，都用到图纸。
而表示图纸的元素就是直线、圆弧和一些标注信息NURBS（Non-Uniform Rational B-Splines）——非均匀有理B样条方法，主要是为了找到与描述自由型曲线曲面的B样条方法既相统一，又能精确表示二次曲线曲面（圆弧、抛物线、椭圆等）的数学方法NURBS一种最一般最有代表性的表示方法。
既为标准解析形状(即前面提到的初等曲线曲面)，又为自由型曲线曲面的精确表示与设计提供了一个公共的数学形式NURBS广泛应用在各种3D造型系统中，如3DMAX，Maya等非均匀——Non-Uniform：指节点向量的值与间距可以为任意值。
这样我们可以在不同区间上得到不同的混合函数形状，为自由控制曲线形状提供了更大自由有理——Rational：是指每个NURBS物体都可以用有理多项式形式表达式来定义Some years ago a few researchers joked about NURBS, saying that the acronym really stands for NOBODY Understands Rational B-Splines2、NURBS曲线的定义NURBS（Non-Uniform Rational B-Splines）——非均匀有理B样条的缩写有理（Rational）函数是两个多项式之比。
因此，有理B 样条可描述为：∑∑∑=====n i k i i n i k i i n i k i i i u R P u B u B P u P 0,0,0,)()()()(ωω∑==n j k j j k i i k i u B u B u R 0,,,)()()(ωω:i ω是控制顶点的权因子。
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,11111,,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 ip 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 时，∑-==ipi j pj u N1)(,。
第二章三维形态基本建模方法第一节形体的空间定位及表示方法一、空间、物体和结构我们每天的生活发生在三维环境中，而且充满着三维物体，我们总是看到、感到三维。
当设计实体模型时，我们通常认为许多事情理所当然。
但在用计算机对三维场景模型化时，那么我们不得不熟悉大量的计算机软件工具，这些工具可用于模型化物体和环境。
在描述三维场景的三维模型化软件中使用的许多基本约定是基于各种行业中使用的传统约定。
例如，建筑师为了用一个简明的方法表达他们设计的空间，使用各种涉及测量、构图和定序的约定。
即使简单的矩形房间设计也要测量多次，以便于房间的所有构件放在被设计放置的地方。
此外，为了准确地按照设计师的图纸来建造，泥瓦工需要进行测量。
多年来泥瓦工和建筑师已形成约定，如何测量空间、建造物体、在结构中安装，它们的约定是精确、简洁的。
我们用类似的约定来描述用一个计算机程序模拟的三维空间中物体的尺寸、位置和次序。
让我们从定义空间或场景的边界开始三维空间的定义，最简单的方法是想象我们是在一个大立方体内工作。
可以将这个立方体当作我们的空间或环境。
在这个立方体中的物体是可见的，在其外部的物体是不可见的。
在这个空间中的主参考点称为主空间原点。
这个原点通常位于这个空间的中心。
根据模型需要和方案，该点也可放在或重新放在其他点上。
所有三维空间都有3个基本的维：宽度、高度和深度。
表达三维空间中这些维的普遍方法是使用箭头或轴。
通常用字母X表示标记三维空间宽度的轴；用Y表示标记三维空间高度的轴；用Z表示标记三维空间深度的轴。
这三个轴交叉的空间点就是主坐标原点。
直角坐标系可以用来定义三维空间中特定的位置，精确定位三维空间中物体的点。
René Descartes是一位18世纪法国的哲学家和数学家，他正式使用标记为X、Y、Z的3个轴表示三维空间中维的思想。
他推导出的坐标系称为笛卡尔坐标系，在该系统中每个轴被分成许多测量单位。
原理上，这些单位是抽象的值，它可表示不同的测量单位和维刻度。
图6.双二次（2x2）B 样条曲面6.B 样条曲线曲面和 NURBS 曲线曲面的C 语言实现算法源程序 #ifndef _mynu rbs_h#ifndef MYNURBS H学习小结：前面学习了 Bezier 曲线，B 样条基函数和 B 样条曲线的一些基础知识。
掌握关 键问题是一条B 样条曲线间的多段曲线的光滑连接。
因为现在是用多段 Bezier 曲线来描绘 一条B 样条曲线，所以问题变为两段 Bezier 曲线间光滑连接。
两段 Bezier 曲线段（3次） B1和B2光滑连接的条件： （1） .要求B1和B2有共同的连接点，即 &连续。
（2） .要求B1和B2在连接点处有成比例的一阶导数，即 G 连续。
由端点处的一阶导数 1 B1(1) 3(F 3 P 2), B 2(0) 3(Q i Q 。
)，为实现 G 连续，则有: B 2(0) B1(1) 即: Q 1 Q o R P 2 这也表明，P 2, B （Q o ）,Q 三点共线。
如下图表示了一条 3次B 样条曲线的所有控制多边形: E(P 1) P2 (P/ P 用P 0 10 P4 图5.3次B 样条曲线和所有控制多边形图5中，P0至P6为原始3次B 样条曲线控制多边形顶点， P 0至P 12是计算后最终形成样条曲线控制多边形顶点。
#include "gl\gl.h"#include "math.h"//* ************* B样条基函数计算部分************** // 确定参数u 所在的节点区间下标//n=m-p-1//m 为节点矢量U[] 的最大下标//p 为B 样条函数次数int FindSource(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 BasisFunction(int i,int p,float u,float U[],float N[]){int j,di,dp,k;float tul,tur,left,right;float tmpN[50][50];for(k=0;k<=p;k++){dp=0;for(di=i+p-k;di>=i-k;di--){if(u>=U[di]&&u<U[di+1])tmpN[di][0]=1;elsetmpN[di][0]=0;dp+=1;for(j=1;j<dp;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 样条函数次数P>2void 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 tmpN[50][50];for(k=0;k<=p;k++){dp=0;for(di=i+p-k;di>=i-k;di--){if(u>=U[di]&&u<U[di+1])tmpN[di][0]=1;elsetmpN[di][0]=0;dp+=1;for(j=1;j<dp;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;else right=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[] 中void BernsteinFunc(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);}}// 获取p 次Bezier 曲线上的lines 个点的值void BezierPoint(int p,float px[],float py[],float pz[],int lines,float tmp[][3]){float BC[20];int i,j;for(j=0;j<=lines;j++){double t=j/(float)lines;BernsteinFunc(p,t,BC); tmp[j][0]=tmp[j][1]=tmp[j][2]=0;for(i=0;i<p+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 个点的值void NBezierPoint(int p,float px[],float py[],float pz[],float lines,float tmp[][4])pw[],int {float x,y,z,w,BC[20];int i,j;for(j=0;j<=lines;j++){double t=j/(float)lines;BernsteinFunc(p,t,BC);x=y=z=w=0;for(i=0;i<p+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 曲线void Bezier(int p,float px[],float py[],float pz[],int lines){float pt[100][3];int j;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 曲线void NBezier(int p,float px[],float py[],float pz[],float w[],int lines){float pt[100][4];int j;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) 方向上的网格数void BezierFacePoint(int p,int u,int v,float px[][4],float pz[][4],float pt[161][161][3])py[][4],float {float urx[11][161],ury[11][161],urz[11][161];float tx[11],ty[11],tz[11],tmp[161][3];int i,j,k;for(j=0;j<p+1;j++){for(i=0;i<p+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;k<p+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) 方向上的网格数void NuBezierFacePoint(int p,int u,int v,float px[][4],float pz[][4],float w[][4],floatpy[][4],float pt[161][161][3]){float urx[11][161],ury[11][161],urz[11][161],urw[11][161];float tx[11],ty[11],tz[11],tw[11],tmp[161][4];int i,j,k;for(j=0;j<p+1;j++){for(i=0;i<p+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][1];urz[j][k]=tmp[k][2];urw[j][k]=tmp[k][3];}}for(i=0;i<=v;i++){for(k=0;k<p+1;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 为曲线的次数void BSplineDer(int n,int p,float U[],float P[],float Der[]){float N[100],tmp;int i,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 为曲线的次数void BSplinePoint(int n,int p,float U[],float P[],float Poi[]){float N[100],tmp;int i,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[] 的最大下标void B3SplineControlPoint(int m,float U[],float P[],float CP[]){int n,k,i,cp,p;float Poi[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<cp-1;i+=3){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[] 的最大下标void B2SplineControlPoint(int m,float U[],float P[],float CP[]){int n,k,tm,i,cp,p;float Poi[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<cp-1;i+=2){k=(int)i/2;CP[i]=Poi[k];CP[i+1]=P[tm];tm++;}}// 绘制3 次B 样条曲线//m 为节点矢量U[] 的最大下标void BSpline3L(int m,float U[],float px[],float py[],float pz[]){float pcx[100],pcy[100],pcz[100],drx[4],dry[4],drz[4];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;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<i*3+4;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[] 的最大下标void BSpline2L(int m,float U[],float px[],float py[],float pz[]){float pcx[100],pcy[100],pcz[100],drx[3],dry[3],drz[3];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;j<i*2+3;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[] 的最大下标值void BS3FaceControlPoint(int mu,float U[],int mv,float V[],float py[],float pz[],floatpx[],float pt[100][100][3]){int i,j,k,dp;float tmx[50],tmy[50],tmz[50];float tmpx[50][100],tmpy[50][100],tmpz[50][100];float uvx[100][100],uvy[100][100],uvz[100][100];for(i=0;i<mv-3;i++){dp=i*(mu-3);for(j=dp;j<mu-3+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<mv-3;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[] 的最大下标值void BS2FaceControlPoint(int mu,float U[],int mv,float V[],float py[],float pz[],floatpx[],float pt[100][100][3]){int i,j,k,dp;float tmx[50],tmy[50],tmz[50];float tmpx[50][100],tmpy[50][100],tmpz[50][100];float uvx[100][100],uvy[100][100],uvz[100][100];for(i=0;i<mv-2;i++){dp=i*(mu-2);for(j=dp;j<mu-2+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<mv-2;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];}}}//---------------------------------------------------------------------------- // 设置网格数void SetGridCount(int dt,int tu,int tmk[]){int i,tm;tm=tu%dt;for(i=0;i<dt-1;i++)tmk[i]=(tu-tm)/dt;tmk[dt-1]=tmk[0]+tm;}//----------------------------------------------------------------------------//计算双三次（3x3次）或双二次（2x2次）B样条曲面上所有的点并保存在bs[][][] 中//nu,mv 分别为节点矢量U[],V[] 的最大下标//uk,vk 分别为B 样条曲面(u,v) 方向上的网格数//p 为曲面的次数void BSplineFace(int p,int nu,float U[],int uk,int mv,float V[],int vk,float px[],float py[],float pz[],float bs[161][161][3]){int udk[20],vdk[20],i,j,k,l,hu,sv,du,dv;float tp[100][100][3],td[161][161][3];float tmx[4][4],tmy[4][4],tmz[4][4];du=nu-2*p;dv=mv-2*p;SetGridCount(du,uk,udk);SetGridCount(dv,vk,vdk);if(p==3)BS3FaceControlPoint(nu,U,mv,V,px,py,pz,tp);if(p==2)BS2FaceControlPoint(nu,U,mv,V,px,py,pz,tp);for(i=0;i<dv;i++){for(k=0;k<du;k++){for(j=i*p;j<p+1+i*p;j++) for(l=k*p;l<p+1+k*p;l++) {tmx[j-i*p][l-k*p]=tp[l][j][0];tmy[j-i*p][l-k*p]=tp[l][j][1];tmz[j-i*p][l-k*p]=tp[l][j][2];}BezierFacePoint(p,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][0]=td[sv-i*vdk[0]][hu-k*udk[0]][0];bs[sv][hu][1]=td[sv-i*vdk[0]][hu-k*udk[0]][1];bs[sv][hu][2]=td[sv-i*vdk[0]][hu-k*udk[0]][2];}}}}//-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-* Nurbs 样条曲线曲面部************** // 计算Nurbs 曲线上的点（u 所对应的所有点）保存在Poi[] 中//n=m-p-1//p 为曲线的次数void NurbsPoint(int n,int p,float U[],float P[],float W[],float Poi[]){float N[100],tmp,tmw;int i,j;for(i=p+1;i<=n;i++){BasisFunction(i,p,U[i],U,N); tmp=0;tmw=0;for(j=i;j>=i-p;j--){tmp+=N[j]*P[j]*W[j]; tmw+=N[j]*W[j];}Poi[i-p]=tmp/tmw;}}// 计算Nurbs 曲线的1 阶导矢( u 所对应的所有点)保存在Der[] 中//n=m-p-1 //p 为曲线的次数void NurbsDer(int n,int p,float U[],float P[],float W[],float Der[]){float N[100],CP[100],NW[100],tmp,tw;int i,j;NurbsPoint(n,p,U,P,W,CP);BSplinePoint(n,p,U,W,NW); for(i=p+1;i<=n;i++){DerBasisFunc(i,p,U[i],U,N); tmp=0;tw=0;for(j=i;j>=i-p;j--){tmp+=N[j]*P[j]*W[j]; tw+=N[j]*W[j];}Der[i-p]=(tmp-tw*CP[i-p])/NW[i-p];}}// 计算3 次Nurbs 曲线上的所有控制多边形保存在CP[] 中//m 为节点矢量U[] 的最大下标void Nurbs3ControlPoint(int m,float U[],float P[],float W[],float CP[]){int n,k,i,cp,p;float Poi[100],Der[100],add;p=3;n=m-p-1;NurbsPoint(n,p,U,P,W,Poi);{NurbsDer(n,p,U,P,W,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<cp-1;i+=3){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 次Nurbs 曲线上的所有控制多边形保存在CP[] 中//m 为节点矢量U[] 的最大下标void Nurbs2ControlPoint(int m,float U[],float P[],float W[],float CP[]){int n,k,tm,i,cp,p;float Poi[100];p=2;n=m-p-1;NurbsPoint(n,p,U,P,W,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<cp-1;i+=2){k=(int)i/2;CP[i]=Poi[k];CP[i+1]=P[tm];tm++;// 绘制3 次Nurbs 样条曲线//m 为节点矢量U[] 的最大下标void Nurbs3L(int m,float U[],float px[],float py[],float pz[],float W[]){float pcx[100],pcy[100],pcz[100],drx[4],dry[4],drz[4]; float pcw[100],drw[4];int i,j,tmcp;Nurbs3ControlPoint(m,U,px,W,pcx);Nurbs3ControlPoint(m,U,py,W,pcy);Nurbs3ControlPoint(m,U,pz,W,pcz);B3SplineControlPoint(m,U,W,pcw); tmcp=m-7;for(i=0;i<=tmcp;i++){{for(j=i*3;j<i*3+4;j++)drx[j-i*3]=pcx[j];dry[j-i*3]=pcy[j];drz[j-i*3]=pcz[j]; drw[j-i*3]=pcw[j];}NBezier(3,drx,dry,drz,drw,20);}}// 绘制2 次Nurbs 样条曲线//m 为节点矢量U[] 的最大下标void Nurbs2L(int m,float U[],float px[],float py[],float pz[],float W[]){float pcx[100],pcy[100],pcz[100],drx[3],dry[3],drz[3]; float pcw[100],drw[3];int i,j,tmcp;Nurbs2ControlPoint(m,U,px,W,pcx);Nurbs2ControlPoint(m,U,py,W,pcy);Nurbs2ControlPoint(m,U,pz,W,pcz);B2SplineControlPoint(m,U,W,pcw); tmcp=m-5;for(i=0;i<=tmcp;i++){for(j=i*2;j<i*2+3;j++){drx[j-i*2]=pcx[j];dry[j-i*2]=pcy[j];drz[j-i*2]=pcz[j];drw[j-i*2]=pcw[j];}NBezier(2,drx,dry,drz,drw,20);}}// 计算双三次( 3x3) Nurbs 样条曲面所有控制多边形顶点，并保存在pt[][][] //mu,mv 分别为节点矢量U[],V[] 的最大下标值void Nurbs3FControlPoint(int mu,float U[],int mv,float V[],float py[],float pz[],float W[],float pt[100][100][4]){int i,j,k,dp;float tmx[50],tmy[50],tmz[50],tmw[50];float tmpx[50][100],tmpy[50][100],tmpz[50][100],tmpw[50][100];float uvx[100][100],uvy[100][100],uvz[100][100],uvw[100][100];for(i=0;i<mv-3;i++){dp=i*(mu-3);for(j=dp;j<mu-3+dp;j++){tmx[j-dp]=px[j];tmy[j-dp]=py[j];tmz[j-dp]=pz[j];tmw[j-dp]=W[j];}Nurbs3ControlPoint(mu,U,tmx,tmw,tmpx[i]);Nurbs3ControlPoint(mu,U,tmy,tmw,tmpy[i]);Nurbs3ControlPoint(mu,U,tmz,tmw,tmpz[i]);B3SplineControlPoint(mu,U,tmw,tmpw[i]);}for(i=0;i<3*mu-17;i++){for(j=0;j<mv-3;j++){tmx[j]=tmpx[j][i];tmy[j]=tmpy[j][i];tmz[j]=tmpz[j][i];tmw[j]=tmpw[j][i]; 中px[],float}Nurbs3ControlPoint(mv,V,tmx,tmw,uvx[i]); Nurbs3ControlPoint(mv,V,tmy,tmw,uvy[i]); Nurbs3ControlPoint(mv,V,tmz,tmw,uvz[i]); B3SplineControlPoint(mv,V,tmw,uvw[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];pt[i][k][3]=uvw[i][k];}}}// 计算双二次( 2x2) Nurbs 样条曲面所有控制多边形顶点，并保存在pt[][][] //mu,mv 分别为节点矢量U[],V[] 的最大下标值void Nurbs2FControlPoint(int mu,float U[],int mv,float V[],float py[],float pz[],float W[],float pt[100][100][4]){int i,j,k,dp;float tmx[50],tmy[50],tmz[50],tmw[50];float tmpx[50][100],tmpy[50][100],tmpz[50][100],tmpw[50][100];float uvx[100][100],uvy[100][100],uvz[100][100],uvw[100][100];for(i=0;i<mv-2;i++){dp=i*(mu-2);for(j=dp;j<mu-2+dp;j++){tmx[j-dp]=px[j];tmy[j-dp]=py[j];tmz[j-dp]=pz[j];tmw[j-dp]=W[j];}Nurbs2ControlPoint(mu,U,tmx,tmw,tmpx[i]);Nurbs2ControlPoint(mu,U,tmy,tmw,tmpy[i]);Nurbs2ControlPoint(mu,U,tmz,tmw,tmpz[i]);B2SplineControlPoint(mu,U,tmw,tmpw[i]);}for(i=0;i<2*mu-7;i++){for(j=0;j<mv-2;j++){tmx[j]=tmpx[j][i];tmy[j]=tmpy[j][i];tmz[j]=tmpz[j][i]; 中px[],floattmw[j]=tmpw[j][i];}Nurbs2ControlPoint(mv,V,tmx,tmw,uvx[i]);Nurbs2ControlPoint(mv,V,tmy,tmw,uvy[i]);Nurbs2ControlPoint(mv,V,tmz,tmw,uvz[i]);B2SplineControlPoint(mv,V,tmw,uvw[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];pt[i][k][3]=uvw[i][k];}}}//----------------------------------------------------------------------------//计算双三次(3x3次)或双二次(2x2次)Nurbs样条曲面上所有的点并保存在//nu,mv 分别为节点矢量U[],V[] 的最大下标〃uk,vk 分别为B样条曲面(u,v)方向上的网格数//p 为曲面的次数void NurbsFace(int p,int nu,float U[],int uk,int mv,float V[],int vk,float px[],float py[],float pz[],float bs[161][161][3]){int udk[20],vdk[20],i,j,k,l,hu,sv,du,dv;float tp[100][100][4],td[161][161][3];float tmx[4][4],tmy[4][4],tmz[4][4],tmw[4][4];du=nu-2*p;dv=mv-2*p;SetGridCount(du,uk,udk);SetGridCount(dv,vk,vdk);if(p==3)Nurbs3FControlPoint(nu,U,mv,V,px,py,pz,w,tp);if(p==2)Nurbs2FControlPoint(nu,U,mv,V,px,py,pz,w,tp);for(i=0;i<dv;i++){for(k=0;k<du;k++){for(j=i*p;j<p+1+i*p;j++)for(l=k*p;l<p+1+k*p;l++){tmx[j-i*p][l-k*p]=tp[l][j][0];tmy[j-i*p][l-k*p]=tp[l][j][1];tmz[j-i*p][l-k*p]=tp[l][j][2]; bs[][][] 中w[],floattmw[j-i*p][l-k*p]=tp[l][j][3];}NuBezierFacePoint(p,udk[k],vdk[i],tmx,tmy,tmz,tmw,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][0]=td[sv-i*vdk[0]][hu-k*udk[0]][0];bs[sv][hu][1]=td[sv-i*vdk[0]][hu-k*udk[0]][1];bs[sv][hu][2]=td[sv-i*vdk[0]][hu-k*udk[0]][2];}}}} for(j=0;j<=u;j++) //* ************* 绘制曲面部*******************// 计算多边形的外法线返回值 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;} //----------------------------------------------------------------------------// 显示 B 样条曲面//fill 取值为 0 或 1void ShowSurface(int u,int v,float bs[161][161][3],int fill) {int i,j;float x[3],y[3],z[3],tmn[3];for(i=0;i<=v;i++){if(fill!=0){x[0]=bs[i][j][0]; x[1]=bs[i+1][j][0]; x[2]=bs[i+1][j+1][0];y[0]=bs[i][j][1]; y[1]=bs[i+1][j][1]; y[2]=bs[i+1][j+1][1];z[0]=bs[i][j][2];z[1]=bs[i+1][j][2]; z[2]=bs[i+1][j+1][2];getN(x,y,z,tmn); glEnable(GL_NORMALIZE); glBegin(GL_QUADS);glNormal3f(tmn[0],tmn[1],tmn[2]); if(j<u){glVertex3f(bs[i][j][0],bs[i][j][1],bs[i][j][2]);glVertex3f(bs[i][j+1][0],bs[i][j+1][1],bs[i][j+1][2]);}if(i<v){glVertex3f(bs[i+1][j+1][0],bs[i+1][j+1][1],bs[i+1][j+1][2]);glVertex3f(bs[i+1][j][0],bs[i+1][j][1],bs[i+1][j][2]);}glEnd(); glDisable(GL_NORMALIZE);}else{ glBegin(GL_LINES); if(j<u) {glVertex3f(bs[i][j][0],bs[i][j][1],bs[i][j][2]);glVertex3f(bs[i][j+1][0],bs[i][j+1][1],bs[i][j+1][2]); }if(i<v){glVertex3f(bs[i][j][0],bs[i][j][1],bs[i][j][2]);glVertex3f(bs[i+1][j][0],bs[i+1][j][1],bs[i+1][j][2]);} glEnd();}#endif#endif#include 建立库文件“ myNurbs.h ”保存在include 目录下，在文件开始处写上"myNurbs.h" 即可用各函数的功能。
