bezier曲线和bezier曲面

如果想理解贝塞尔曲面没有对其数学基本的认识是很难的，如果你不愿意读这一部分或者你已经知道了关于她的数学知识你可以跳过。首先我会描述贝塞尔曲线再介绍生成贝塞尔曲面。
奇怪的是，如果你用过一个图形程序，你就已经熟悉了贝塞尔曲线，也许你接触的是另外的名称。它们是画曲线的最基本的方法，而且通常被表示成一系列点，其中有两个点与两端点表示左右两端的切线。下图展示了一个例子。

这是最基础的贝塞尔曲线（长点的由很多点在一起（多到你都没发现））。这个曲线由4个点定义，有2个端点和2个中间控制点。对计算机而言这些点都是一样的，但是特意的我们通常把前后两对点分别连接，因为他们的连线与短点相切。曲线是一个参数化曲线，画的时候从曲线上平均找几点连接。这样你可以控制曲线曲面的精度（和计算量）。最通常的方法是远距离少细分近距离多细分，对视点，看上去总是很完好的曲面而对速度的影响总是最小。
贝塞尔曲面基于一个基本方程，其他复杂的都是基于此。方程为：

t + (1 - t) = 1

看起来很简单不是？的确是的，这是最基本的贝塞尔曲线，一个一维的曲线。你也许从术语中猜到，贝塞尔曲线是多项式形式的。从线性代数知，一个一维的多项式是一条直线，没多大意思。好，因为基本方程对所有t都成立，我们可以平方，立方两边，怎么都行，等式都是成立的，对吧？好，我们试试立方。

(t + (1-t))^3 = 1^3

t^3 + 3*t^2*(1-t) + 3*t*(1-t)^2 + (1-t)^3 = 1

这是我们最常用的计算贝塞尔曲面的方程，a)她是最低维的不需要在一个平面内的多项式（有4个控制点），而且b)两边的切线互相没有联系（对于2维的只有3个控制点）。那么你看到了贝塞尔曲线了吗？呵呵，我们都没有，因为我还要加一个东西。
好，因为方程左边等于1，可以肯定如果你把所有项加起来还是等于1。这是否意味着在计算曲线上一点时可以以此决定该用每个控制点的多少呢？（答案是肯定的）你对了！当我们要计算曲线上一点的值我们只需要用控制点（表示为向量）乘以每部分再加起来。基本上我们要用0<=t<=1，但不是必要的。不明白了把？这里有函数：

P1*t^3 + P2*3*t^2*(1-t) + P3*3*t*(1-t)^2 + P4*(1-t)^3 = Pnew

因为多项式是连续的，有一个很好的办法在4个点间插值。曲线仅经过P1，P4，分别当t=1,0。
好，一切都很好，但现在我怎么把这个用在3D里呢？其实很简单，为了做一个贝塞尔曲面，你需要16个控制点，（4*4），和2个变量t,v。你要做的是计算在分量v的沿4条平行曲线的点，再用这4个点计算在分量t的点。计算了足够的这些点，我们可以用三角带连接他们，画出贝塞尔曲面。

 

恩，我认为现在已经有足够的数学背景了，看代码把！

#include<math.h>       //数学库
#include<stdio.h>       //标准输入输出库
#include<stdlib.h>      //标准库

typedef struct point_3d{      //3D点的结构
 double x, y, z;
} POINT_3D;

typedef struct bpatch{      //贝塞尔面片结构
 POINT_3D anchors[4][4];     //由4x4网格组成
 GLuint  dlBPatch;     //绘制面片的显示列表名称
 GLuint  texture;     //面片的纹理
} BEZIER_PATCH;

BEZIER_PATCH  mybezier;     //创建一个贝塞尔曲面结构
BOOL   showCPoints=TRUE;    //是否显示控制点
int   divs =7;     //细分精度，控制曲面的显示精度

  
 以下是一些简单的向量数学的函数。如果你是C++爱好者你可以用一个顶点类（保证其为3D的）。

// 两个向量相加，p=p+q
POINT_3D pointAdd(POINT_3D p, POINT_3D q) {
 p.x +=q.x;  p.y +=q.y;  p.z += q.z;
 return p;
}

// 向量和标量相乘p=c*p
POINT_3D pointTimes(double c, POINT_3D p) {
 p.x *= c; p.y *=c; p.z *= c;
 return p;
}

// 创建一个3D向量
POINT_3D makePoint(double a, double b, double c) {
 POINT_3D p;
 p.x = a; p.y =b; p.z = c;
 return p;
}

  
 这基本上是用C写的3维的基本函数，她用变量u和4个顶点的数组计算曲线上点。每次给u加上一定值，从0到1，我们可得一个很好的近似曲线。
求值器基于Bernstein多项式定义曲线，定义p(u ')为：
p(u')=∑Bni(u')Ri

这里Ri为控制点
Bni(u')=[ni]u'i(1-u')n-i

且00=1,[n0]=1

u'=(u-u1)/(u2-u1)

当为贝塞尔曲线时，控制点为4，相应的4个Bernstein多项式为:
1、B30 =（1-u）3
2、B31 =3u(1-u)2
3、B32 =3u2(1-u)
4、B33 =u3

 

// 计算贝塞尔方程的值
// 变量u的范围在0-1之间
POINT_3D Bernstein(float u, POINT_3D *p) {
 POINT_3D a, b, c, d, r;

 a = pointTimes(pow(u,3), p[0]);
 b = pointTimes(3*pow(u,2)*(1-u), p[1]);
 c = pointTimes(3*u*pow((1-u),2), p[2]);
 d = pointTimes(pow((1-u),3), p[3]);

 r = pointAdd(pointAdd(a, b), pointAdd(c,d));

 return r;
}

  
 这个函数完成共享工作，生成所有三角带，保存在displaylist。我们这样就不需要每贞都重新计算曲面，除了当其改变时。另外，你可能想用一个很酷的效果，用MORPHING教程改变控制点位置。这可以做一个很光滑，有机的，morphing效果，只要一点点开销（你只要改变16个点，但要从新计算）。“最后”的数组元素用来保存前一行点，（因为三角带需要两行）。而且，纹理坐标由表示百分比的u,v来计算（平面映射）。
还有一个我们没做的是计算法向量做光照。到了这一步，你基本上有2种选择。第一是找每个三角形的中心计算X，Y轴的切线，再做叉积得到垂直与两向量的向量，再归一化，得到法向量。或者（恩，这是更好的方法）你可以直接用三角形的法矢（用你最喜欢的方法计算）得到一个近似值。我喜欢后者；我认为不值得为了一点点真实感影响速度。

// 生成贝塞尔曲面的显示列表
GLuint genBezier(BEZIER_PATCH patch, int divs) {
 int  u = 0,v;
 float  py, px,pyold;
 GLuint  drawlist=glGenLists(1);   //创建显示列表
 POINT_3D temp[4];
 POINT_3D *last =(POINT_3D*)malloc(sizeof(POINT_3D)*(divs+1)); //更具每一条曲线的细分数，分配相应的内存

 if (patch.dlBPatch !=NULL)     //如果显示列表存在则删除
  glDeleteLists(patch.dlBPatch,1);

 temp[0] =patch.anchors[0][3];    //获得u方向的四个控制点
 temp[1] = patch.anchors[1][3];
 temp[2] = patch.anchors[2][3];
 temp[3] = patch.anchors[3][3];

 for (v=0;v<=divs;v++){     //根据细分数，创建各个分割点额参数
  px =((float)v)/((float)divs);    
 // 使用Bernstein函数求的分割点的坐标
  last[v] = Bernstein(px,temp);
 }

 glNewList(drawlist,GL_COMPILE);    //创建一个新的显示列表
 glBindTexture(GL_TEXTURE_2D,patch.texture);   //邦定纹理

 for (u=1;u<=divs;u++) {
  py   =((float)u)/((float)divs);   //计算v方向上的细分点的参数
  pyold =((float)u-1.0f)/((float)divs);  //上一个v方向上的细分点的参数

  temp[0] = Bernstein(py,patch.anchors[0]);  //计算每个细分点v方向上贝塞尔曲面的控制点
  temp[1] = Bernstein(py,patch.anchors[1]);
  temp[2] = Bernstein(py,patch.anchors[2]);
  temp[3] = Bernstein(py,patch.anchors[3]);

  glBegin(GL_TRIANGLE_STRIP);    //开始绘制三角形带

  for(v=0;v<=divs;v++) {
   px =((float)v)/((float)divs);  //沿着u轴方向顺序绘制

   glTexCoord2f(pyold,px);   //设置纹理坐标
   glVertex3d(last[v].x,last[v].y, last[v].z); // 绘制一个顶点

   last[v] =Bernstein(px, temp);  //创建下一个顶点
   glTexCoord2f(py,px);   //设置纹理
   glVertex3d(last[v].x,last[v].y, last[v].z); // 绘制新的顶点
  }

  glEnd();      //结束三角形带的绘制
 }

 glEndList();      //显示列表绘制结束

 free(last);      //释放分配的内存
 returndrawlist;      //返回创建的显示列表
}

  
 这里我们调用一个我认为有一些很酷的值的矩阵。

void initBezier(void) {
 mybezier.anchors[0][0] =makePoint(-0.75, -0.75, -0.50); //设置贝塞尔曲面的控制点
 mybezier.anchors[0][1] =makePoint(-0.25, -0.75, 0.00);
 mybezier.anchors[0][2] = makePoint(0.25, -0.75,  0.00);
 mybezier.anchors[0][3] = makePoint(0.75, -0.75, -0.50);
 mybezier.anchors[1][0] =makePoint(-0.75, -0.25, -0.75);
 mybezier.anchors[1][1] =makePoint(-0.25, -0.25, 0.50);
 mybezier.anchors[1][2] = makePoint(0.25, -0.25,  0.50);
 mybezier.anchors[1][3] = makePoint(0.75, -0.25, -0.75);
 mybezier.anchors[2][0] =makePoint(-0.75,  0.25, 0.00);
 mybezier.anchors[2][1] =makePoint(-0.25, 0.25, -0.50);
 mybezier.anchors[2][2] = makePoint(0.25,  0.25, -0.50);
 mybezier.anchors[2][3] = makePoint(0.75,  0.25,  0.00);
 mybezier.anchors[3][0] =makePoint(-0.75, 0.75, -0.50);
 mybezier.anchors[3][1] =makePoint(-0.25, 0.75, -1.00);
 mybezier.anchors[3][2] = makePoint(0.25,  0.75, -1.00);
 mybezier.anchors[3][3] = makePoint(0.75,  0.75, -0.50);
 mybezier.dlBPatch =NULL;     //默认的显示列表为0
}

  
 这是一个优化的调位图的函数。可以很简单的把他们放进一个简单循环里调一组。

// 加载一个*.bmp文件，并转化为纹理

BOOL LoadGLTexture(GLuint *texPntr, char* name)
{
 BOOL success = FALSE;
 AUX_RGBImageRec *TextureImage = NULL;

 glGenTextures(1,texPntr);     //生成纹理1

 FILE* test=NULL;
 TextureImage = NULL;

 test = fopen(name,"r");     
 if (test != NULL){      
  fclose(test);      
  TextureImage =auxDIBImageLoad(name);   
 }

 if (TextureImage != NULL){     
  success = TRUE;

  // 邦定纹理
  glBindTexture(GL_TEXTURE_2D,*texPntr);
  glTexImage2D(GL_TEXTURE_2D, 0,3, TextureImage->sizeX,TextureImage->sizeY, 0, GL_RGB, GL_UNSIGNED_BYTE,TextureImage->data);
  glTexParameteri(GL_TEXTURE_2D,GL_TEXTURE_MIN_FILTER,GL_LINEAR);
  glTexParameteri(GL_TEXTURE_2D,GL_TEXTURE_MAG_FILTER,GL_LINEAR);
 }

 if (TextureImage->data)
  free(TextureImage->data);

 return success;
}

  
 仅仅加了曲面初始化在这。你每次建一个曲面时都会用这个。再一次，这里是一个用C++的好地方（贝塞尔曲面类？）。

intInitGL(GLvoid)       //初始化OpenGL
{
 glEnable(GL_TEXTURE_2D);     //使用2D纹理
 glShadeModel(GL_SMOOTH);     //使用平滑着色
 glClearColor(0.05f, 0.05f, 0.05f,0.5f);   //设置黑色背景
 glClearDepth(1.0f);     //设置深度缓存
 glEnable(GL_DEPTH_TEST);     //使用深度缓存
 glDepthFunc(GL_LEQUAL);     //设置深度方程
 glHint(GL_PERSPECTIVE_CORRECTION_HINT,GL_NICEST);  

 initBezier();      //初始化贝塞尔曲面
 LoadGLTexture(&(mybezier.texture),"./Data/NeHe.bmp");  //载入纹理
 mybezier.dlBPatch = genBezier(mybezier,divs);  // 创建显示列表

 returnTRUE;      //初始化成功
}

  
 首先调贝塞尔displaylist。再（如果边线要画）画连接控制点的线。你可以用SPACE键开关这个。

intDrawGLScene(GLvoid) {      //绘制场景
 int i, j;
 glClear(GL_COLOR_BUFFER_BIT |GL_DEPTH_BUFFER_BIT);  
 glLoadIdentity();      
 glTranslatef(0.0f,0.0f,-4.0f);     //移入屏幕4个单位
 glRotatef(-75.0f,1.0f,0.0f,0.0f);
 glRotatef(rotz,0.0f,0.0f,1.0f);     //旋转一定的角度

 glCallList(mybezier.dlBPatch);     //调用显示列表，绘制贝塞尔曲面

 if (showCPoints){       //是否绘制控制点
  glDisable(GL_TEXTURE_2D);
  glColor3f(1.0f,0.0f,0.0f);
  for(i=0;i<4;i++){      //绘制水平线
   glBegin(GL_LINE_STRIP);
   for(j=0;j<4;j++)
    glVertex3d(mybezier.anchors[i][j].x,mybezier.anchors[i][j].y, mybezier.anchors[i][j].z);
   glEnd();
  }
  for(i=0;i<4;i++){      //绘制垂直线
   glBegin(GL_LINE_STRIP);
   for(j=0;j<4;j++)
    glVertex3d(mybezier.anchors[j][i].x,mybezier.anchors[j][i].y, mybezier.anchors[j][i].z);
   glEnd();
  }
  glColor3f(1.0f,1.0f,1.0f);
  glEnable(GL_TEXTURE_2D);
 }

 returnTRUE;       //成功返回
}

  
 KillGLWindow()函数没有改动
 
  
  
 CreateGLWindow（）函数没有改动 
  
  
 我在这里加了旋转曲面的代码，增加/降低分辨率，显示与否控制点连线。

   if(keys[VK_LEFT]) rotz -=0.8f;  // 按左键，向左旋转
   if(keys[VK_RIGHT]) rotz += 0.8f; //按右键，向右旋转
   if(keys[VK_UP]){    //按上键，加大曲面的细分数目
    divs++;
    mybezier.dlBPatch= genBezier(mybezier, divs); // 更新贝塞尔曲面的显示列表
    keys[VK_UP]= FALSE;
   }
   if(keys[VK_DOWN] && divs> 1){    //按下键，减少曲面的细分数目
    divs--;
    mybezier.dlBPatch= genBezier(mybezier, divs); // 更新贝塞尔曲面的显示列表
    keys[VK_DOWN]= FALSE;
   }
   if(keys[VK_SPACE]){     //按空格切换控制点的可见性
    showCPoints= !showCPoints;
    keys[VK_SPACE]= FALSE;
   }