Voronoi diagram......

　　Voronoi图，又叫泰森多边形或Dirichlet图，它是由一组由连接两邻点直线的垂直平分线组成的连续多边形组成。N个在平面上有区别的点，按照最邻近原则划分平面；每个点与它的最近邻区域相关联。Delaunay三角形是由与相邻Voronoi多边形共享一条边的相关点连接而成的三角形。Delaunay三角形的外接圆圆心是与三角形相关的Voronoi多边形的一个顶点。Voronoi三角形是Delaunay图的偶图；

　　对于给定的初始点集P，有多种三角网剖分方式，其中Delaunay三角网具有以下特征：

　　1、Delaunay三角网是唯一的；
　　2、三角网的外边界构成了点集P的凸多边形“外壳”；
　　3、没有任何点在三角形的外接圆内部，反之，如果一个三角网满足此条件，那么它就是Delaunay三角网。
　　4、如果将三角网中的每个三角形的最小角进行升序排列，则Delaunay三角网的排列得到的数值最大，从这个意义上讲，Delaunay三角网是“最接近于规则化的“的三角网。

　　Delaunay三角形网的特征又可以表达为以下特性：

　　1、在Delaunay三角形网中任一三角形的外接圆范围内不会有其它点存在并与其通视，即空圆特性；
　　2、在构网时，总是选择最邻近的点形成三角形并且不与约束线段相交；
　　3、形成的三角形网总是具有最优的形状特征，任意两个相邻三角形形成的凸四边形的对角线如果可以互换的话，那么两个三角形6个内角中最小的角度不会变大；
　　4、不论从区域何处开始构网，最终都将得到一致的结果，即构网具有唯一性。

Reference
http://www.docin.com/p-22007734.html
http://dragonlhw.blog.163.com/blog/static/20806092006916579332/

Gradient, divergenz,rotation

梯度: 运算的对像是纯量,运算出来的结果会是向量在一个纯量场中,梯度的计算结果会是"在每个位置都算出一个向量,而这个向量的方向会是在任何一点上从其周围纯量值最小处指向周围纯量值最大处.
例子：举例子来讲会比较简单,如果现在的纯量场用一座山来表示,纯量值越大的地方越高,反之则越低.经过梯度这个运操作数的运算以后,会在这座山的每一个点上都算出一个向量,这个向量会指向每个点最陡的那个方向,而向量的大小则代表了这个最陡的方向到底有多陡.

散度: 运算的对像是向量,运算出来的结果会是纯量。散度的作用对像是向量场,如果现在我们考虑任何一个点(或者说这个点的周围极小的一块区域),在这个点上,向量场的发散程度,如果是正的,代表这些向量场是往外散出的；如果是负的,代表这些向量场是往内集中的.
例子: 因为散度的作用对像是向量场,所以就不能用上面所讲的山来想象,这次要想象一个大广场里挤了很多人,如果每个人都在到处走动,是不是可以把每个人的行动都看成是一个向量,假如现在某人放了一个屁,周围的人(可能包含他自己)都想要赶快闪远一点,就会发现,在这块区域的人都往这小块区域以外的方向移动.对啦..这就是散度(你也可以想说是闪远一点的闪度....冷....),
大家如果散得越快,散得人越多,这个散度算出来就就越大.

旋度: 运算的对像是向量,运算出来的结果会是向量。旋度的作用对象也是向量场。
例子:如果现在散开的众人都是直直的往那个屁的反方向散开,这时候你看到这些人的动线是不是就是一个标准的幅射状??不过事实上,每个人在闻到屁的时候是不会确切的知道屁到底是来自哪个方向的.而可能会走错方向,试过之后才发现不对劲,越找越臭.这时候你看到众人的走向不见得就是一个幅射状(大家都径向移动),而可能有一些切向移动的成份在(以屁发点为中心来看)旋度对应的就是这些切向移动的情况,相对来讲,散度对应的其实就是径向移动的情况.而一个屁,虽然可能会像上述的造成一些切向的移动,但理论上来讲,并不会使散开的众人较趋向于顺时钟转,或逆时钟转.在这种情况,顺时钟转的情况可以看作与逆时钟转的情况抵消,因此,在这情况下,旋度仍然是零.也就是说,一个屁能造成散度,而不会造成旋度....
而甚么时候是有旋度的呢??如果这时候音乐一放,大家开始围着中间的营火手拉手跳起土风舞(当然是要绕着营火转的那种啦)这时候就会有旋度没有散度啦.(刚刚一直放屁的那位跑出去找厕所的除外)

以上这三个,有一点一定要记得的.
不论是梯度,散度,旋度,都是一种local的量(纯量,向量),所考虑的都是任何一点(其周围极接近,极小的小范围)的情况.以上举的例子因为要容易了解,所以都是针对二度空间向量为例,而且都是很大的东西,但广场是一个点,营火晚会也是一个点,纳须弥于芥子,这就请自行想象吧.

————————————————————————

简单来说 在三维空间里 对一个面做三维的梯度计算得到得向量就是垂直那个面的法线向量 因为沿著法线的「值」变化最大.

散度的感觉像是点的通量(还是通量密度?) 比方在三维空间中,圈出一个封闭曲面; 通量即该向量(的垂直平面分量)穿过平面的大小 如果缩小封闭曲面的大小到一个点; 该点的通量(通量密度)就是散度.一般点的散度为0 当散度不为0的点表示该点有提供source.
从水流的角度来看 就是那里有水源不断冒出(或流入)

旋度(原PO的卷度?)像是没有流出的量 比方在三维空间中圈出一个很小很小的曲面,旋度即该向量(的平行平面分量)延平面的大小密度(即大小/面积),旋度不为0表示有量在该平面「逗留.
从水流的角度来看 就是有涡流的现象.

__________________________________________________

http://blog.sina.com.cn/s/blog_5512f7650100al2p.html
http://uwb.blog.hexun.com/1573987_d.html

Oops, I got it.

前天出了一个问题（http://leepyzh.blogspot.com/2009/09/cg-lighting.html），是在看Cg例子“09_vertex_lighting”的时候，为了计算光照，他有一个对modelview矩阵求逆的过程，这是由于GPU顶点程序的顶点位于model坐标系，将灯光和眼睛的位置逆变换到物体的位置空间，这样在同一空间里面可以求解。

看了这个例子后，总觉得求逆是一件费时的事，干吗不把顶点坐标和向量变换到世界坐标系？这样也可以在同一空间求解呀。但是修改程序后，感觉灯光的位置有些偏差，百思不得其解。Cg程序如下：

#if 1
//here is my code to change the vertex coord & normal to world space

float4 worldcoordinate = mul(modelToWorld, position);
float3 P = worldcoordinate.xyz;
float4 Normal4;
Normal4[0] = normal[0];
Normal4[1] = normal[1];
Normal4[2] = normal[2];
Normal4[3] = 1;
float4 N4 = mul(modelToWorld, Normal4);
float3 N = N4.xyz;
N = normalize(N);

#else
 //here is the formal code

float3 P = position.xyz;
float3 N = normal;
#endif

今天，看到“18_cube_map_reflection”的时候发现他用到这个想法。呵呵，至少心里放心了。看看例子的代码，比我的简洁多了：
// Compute position and normal in world space

float3 P = mul(modelToWorld, position).xyz;
float3 N = mul((float3x3)modelToWorld, normal);
N = normalize(N);

经过盘查，我终于找到了问题所在，原来是向量相乘的时候出了问题：

Normal4[3] = 1;

其实对于向量来说，第四维是没有意义的。如果赋值为1，那么需要将原点看成起点，在两个变换完成后，然后相减得到新的向量。最简单的方法，将第四维赋值为0。从几何上可以这么理解：平移不改变向量。

到此，似乎问题已经解决。当接着往下看Cg tutorial时，又让我吃了一大惊，他说：

If the modeling transform scales positions nonuniformly, you must multiply normal by the inverse transpose of the modeling matrix ( modelToWorldInvTrans ), rather than simply by modelToWorld.

是呀。旋转和平移时，点的相对位置不会改变，所以直接变换向量没有问题；然后，不均匀的缩放时，平面发生变形，此时，法向不能简单地乘以变换矩阵而得到。

网上一位老兄给出了总结，我想应该是对的：

1.如果变换模型位置的仅仅是rotation 或 translation, 可以用该Matrix直接作用于法线,获得最终的法线.

2.如果传输模型的Matrix 是同比例scale, 用该matrix作用法线后,获得的新法线需要重新normalize

3.如果非同比例scale的matrix，必须用该matrix的inverse后,再transpose来作用于该法线。

4. 具体内容可参考: Real-Time Shader Programming 82-83页。

Finally, 看来“09_vertex_lighting”的求逆阵还是个最好的、不会出错的解决方法。

Cg Vertex prgram & fragment program.

从上图可以看到，顶点程序和片段程序在整个GPU处理过程和流程中所处的位置。顶点程序每个顶点都要处理一次；片段程序是每个片段（顶点插值得到的下一级）要处理一次。

顶点程序输入：
1.未经过任何变换的顶点，坐标和法向为均为物体空间，不是世界坐标系；
2. 纹理坐标；
3. 顶点被（OpenGL程序）设置的颜色。
顶点程序输出：
1. 变换到投影坐标系的顶点坐标【必须】；
2. 纹理坐标【可选，或者是顶点程序要传给片段程序的参数，封装在纹理坐标里】；
3. 顶点光照后的颜色【可选，也可在fragment程序中完成】；

片段程序的输入：
1.顶点经过二次线性插值后的片段的颜色【可选】；
2.纹理坐标【可选】；


片段程序的输出：
1. 片段颜色【必须，纹理映射后的颜色或者按片段进行“精确”光照后的颜色或者直接pass through】。

Cg Lighting......

今天看到Cg中OpenGL的例子“09_vertex_lighting”，例子简洁明了，用生动地展示了基本光照过程。真是爱死Cg了，“The Cg Tutorial"这本书的作者真是太牛了。

light.......
{
float3 P = position.xyz; //current point
float3 N = normal;

// Compute emissive term
float3 emissive = Ke;

// Compute ambient term
float3 ambient = Ka * globalAmbient;

// Compute the diffuse term
float3 L = normalize(lightPosition - P);
float diffuseLight = max(dot(N, L), 0);
float3 diffuse = Kd * lightColor * diffuseLight;

// Compute the specular term
float3 V = normalize(eyePosition - P);
float3 H = normalize(L + V);
float specularLight = pow(max(dot(N, H), 0), shininess);
if (diffuseLight <= 0) specularLight = 0;
float3 specular = Ks * lightColor * specularLight;
color.xyz = emissive + ambient + diffuse + specular;
color.w = 1;
}

几点理解：
1）emissive并不是场景中的光源，它不能照亮其他物体，而自身将会呈现一种单一的颜色；
2）镜面反射取出眼睛和灯光中间向量，与法向求夹角。shininess = 0，光晕最大；shininess增加，呈指数衰减，最后shininess->∽时，只有dot(N, H)=1也就是N和H重合时，采用光斑;
3)在GPU的顶点程序中，传进来的顶点position是未经过任何变换的坐标值，它处于模型坐标系中；而眼睛和光源的位置是在世界坐标系中定义的。所以，GPU计算光照时，必须将眼睛和光源的坐标进行反推到模型坐标系，具体是对modelMatrix求逆，在乘上眼睛和光源的坐标向量；
4)我尝试将modelMatrix传入顶点程序，以将position和normal转入世界坐标系求光照，但是得到错误的结果，思考中......
5)本例子提供了求矩阵逆的函数，弓虽！

/* Invert a row-major (C-style) 4x4 matrix. */
static void invertMatrix(float *out, const float *m)
{
/* Assumes matrices are ROW major. */
#define SWAP_ROWS(a, b) { GLdouble *_tmp = a; (a)=(b); (b)=_tmp; }
#define MAT(m,r,c) (m)[(r)*4+(c)]

double wtmp[4][8];
double m0, m1, m2, m3, s;
double *r0, *r1, *r2, *r3;

r0 = wtmp[0], r1 = wtmp[1], r2 = wtmp[2], r3 = wtmp[3];

r0[0] = MAT(m,0,0), r0[1] = MAT(m,0,1),
r0[2] = MAT(m,0,2), r0[3] = MAT(m,0,3),
r0[4] = 1.0, r0[5] = r0[6] = r0[7] = 0.0,

r1[0] = MAT(m,1,0), r1[1] = MAT(m,1,1),
r1[2] = MAT(m,1,2), r1[3] = MAT(m,1,3),
r1[5] = 1.0, r1[4] = r1[6] = r1[7] = 0.0,

r2[0] = MAT(m,2,0), r2[1] = MAT(m,2,1),
r2[2] = MAT(m,2,2), r2[3] = MAT(m,2,3),
r2[6] = 1.0, r2[4] = r2[5] = r2[7] = 0.0,

r3[0] = MAT(m,3,0), r3[1] = MAT(m,3,1),
r3[2] = MAT(m,3,2), r3[3] = MAT(m,3,3),
r3[7] = 1.0, r3[4] = r3[5] = r3[6] = 0.0;

/* Choose myPivot, or die. */
if (fabs(r3[0])>fabs(r2[0])) SWAP_ROWS(r3, r2);
if (fabs(r2[0])>fabs(r1[0])) SWAP_ROWS(r2, r1);
if (fabs(r1[0])>fabs(r0[0])) SWAP_ROWS(r1, r0);
if (0.0 == r0[0]) {
assert(!"could not invert matrix");
}

/* Eliminate first variable. */
m1 = r1[0]/r0[0]; m2 = r2[0]/r0[0]; m3 = r3[0]/r0[0];
s = r0[1]; r1[1] -= m1 * s; r2[1] -= m2 * s; r3[1] -= m3 * s;
s = r0[2]; r1[2] -= m1 * s; r2[2] -= m2 * s; r3[2] -= m3 * s;
s = r0[3]; r1[3] -= m1 * s; r2[3] -= m2 * s; r3[3] -= m3 * s;
s = r0[4];
if (s != 0.0) { r1[4] -= m1 * s; r2[4] -= m2 * s; r3[4] -= m3 * s; }
s = r0[5];
if (s != 0.0) { r1[5] -= m1 * s; r2[5] -= m2 * s; r3[5] -= m3 * s; }
s = r0[6];
if (s != 0.0) { r1[6] -= m1 * s; r2[6] -= m2 * s; r3[6] -= m3 * s; }
s = r0[7];
if (s != 0.0) { r1[7] -= m1 * s; r2[7] -= m2 * s; r3[7] -= m3 * s; }

/* Choose myPivot, or die. */
if (fabs(r3[1])>fabs(r2[1])) SWAP_ROWS(r3, r2);
if (fabs(r2[1])>fabs(r1[1])) SWAP_ROWS(r2, r1);
if (0.0 == r1[1]) {
assert(!"could not invert matrix");
}

/* Eliminate second variable. */
m2 = r2[1]/r1[1]; m3 = r3[1]/r1[1];
r2[2] -= m2 * r1[2]; r3[2] -= m3 * r1[2];
r2[3] -= m2 * r1[3]; r3[3] -= m3 * r1[3];
s = r1[4]; if (0.0 != s) { r2[4] -= m2 * s; r3[4] -= m3 * s; }
s = r1[5]; if (0.0 != s) { r2[5] -= m2 * s; r3[5] -= m3 * s; }
s = r1[6]; if (0.0 != s) { r2[6] -= m2 * s; r3[6] -= m3 * s; }
s = r1[7]; if (0.0 != s) { r2[7] -= m2 * s; r3[7] -= m3 * s; }

/* Choose myPivot, or die. */
if (fabs(r3[2])>fabs(r2[2])) SWAP_ROWS(r3, r2);
if (0.0 == r2[2]) {
assert(!"could not invert matrix");
}

/* Eliminate third variable. */
m3 = r3[2]/r2[2];
r3[3] -= m3 * r2[3], r3[4] -= m3 * r2[4],
r3[5] -= m3 * r2[5], r3[6] -= m3 * r2[6],
r3[7] -= m3 * r2[7];

/* Last check. */
if (0.0 == r3[3]) {
assert(!"could not invert matrix");
}

s = 1.0/r3[3]; /* Now back substitute row 3. */
r3[4] *= s; r3[5] *= s; r3[6] *= s; r3[7] *= s;

m2 = r2[3]; /* Now back substitute row 2. */
s = 1.0/r2[2];
r2[4] = s * (r2[4] - r3[4] * m2), r2[5] = s * (r2[5] - r3[5] * m2),
r2[6] = s * (r2[6] - r3[6] * m2), r2[7] = s * (r2[7] - r3[7] * m2);
m1 = r1[3];
r1[4] -= r3[4] * m1, r1[5] -= r3[5] * m1,
r1[6] -= r3[6] * m1, r1[7] -= r3[7] * m1;
m0 = r0[3];
r0[4] -= r3[4] * m0, r0[5] -= r3[5] * m0,
r0[6] -= r3[6] * m0, r0[7] -= r3[7] * m0;

m1 = r1[2]; /* Now back substitute row 1. */
s = 1.0/r1[1];
r1[4] = s * (r1[4] - r2[4] * m1), r1[5] = s * (r1[5] - r2[5] * m1),
r1[6] = s * (r1[6] - r2[6] * m1), r1[7] = s * (r1[7] - r2[7] * m1);
m0 = r0[2];
r0[4] -= r2[4] * m0, r0[5] -= r2[5] * m0,
r0[6] -= r2[6] * m0, r0[7] -= r2[7] * m0;

m0 = r0[1]; /* Now back substitute row 0. */
s = 1.0/r0[0];
r0[4] = s * (r0[4] - r1[4] * m0), r0[5] = s * (r0[5] - r1[5] * m0),
r0[6] = s * (r0[6] - r1[6] * m0), r0[7] = s * (r0[7] - r1[7] * m0);

MAT(out,0,0) = r0[4]; MAT(out,0,1) = r0[5],
MAT(out,0,2) = r0[6]; MAT(out,0,3) = r0[7],
MAT(out,1,0) = r1[4]; MAT(out,1,1) = r1[5],
MAT(out,1,2) = r1[6]; MAT(out,1,3) = r1[7],
MAT(out,2,0) = r2[4]; MAT(out,2,1) = r2[5],
MAT(out,2,2) = r2[6]; MAT(out,2,3) = r2[7],
MAT(out,3,0) = r3[4]; MAT(out,3,1) = r3[5],
MAT(out,3,2) = r3[6]; MAT(out,3,3) = r3[7];

#undef MAT
#undef SWAP_ROWS
}

Z Buffer or Depth Buffer.

在投影空间中，z的范围是[-1,+1]。在OpenGL中，默认情况下，最接近眼睛的片断（在近截面上）被映射到0.0，离眼睛最远的片断（在远截面上）映射到1.0。这个映射调整应该是线性的（没有查到正式文献）。当然，还可用glDepthRange调节深度值范围。

在正交投影中距离和Z值的关系是线性的，但在透视投影中却不是的。在透视投影中，情况发生了变化，还是看看映射公式：

z'= f(z)/w' = （Zfar+Znear） / ( Zfar – Znear ) + [2* Zfar*Znear / ( Zfar – Znear )]/z

这个映射是个双曲线：

下面来理解（http://www.cppblog.com/zmj/archive/2008/03/14/14122.html）给出的这几句话：

1. 在透视投影中这种关系是非线性的，而且非线性的程度与Frustum函数中的far/near（或gluPerspective函数中的zFar/zNear）成比例——也不是严格成比例，因为有一个Zfar*Znear 在。

2. 这种非线性在靠近近截面时增加了Z值的精度，而且增加了深度缓存的效率；但是在视见体的其它部分则降低了精度，也就减少了深度缓存的精确性——在上图，将[-Zfar,-Znear]之间等分，呵呵，靠近-Znear得到的Z’的范围显然大得多。

3. 根据经验，far/near比值大于1000会有这种不好的效果。所以一般far/near比值应小于1000。要想解决这个问题，最简单的方法是通过将近截面远离眼睛来降低far/near比值，其唯一的副作用是离眼睛很近的物体可能会被裁减掉，但在特定的应用程序中这很少是个问题，近截面的位置对X、Y坐标的投影没有影响，因此这对图像的影响很小——这个要牢记，呵呵。

Inside CG transformation (III).

最后看一下投影矩阵，这是最复杂的一个变换。我们只看透视投影矩阵，平行投影会简单很多。

在前面Cg Studying（http://leepyzh.blogspot.com/2009/09/cg-studying.html）一文中，我们已经看到了透视投影变换的过程：将平头视锥体变换了边长为2，中心在原点的立方体，以便于后面的裁剪和消隐工作。过程如下：

下面将变换设定在常用的情况上，就是投影平面的中心和xy平面的中心重合，即在上面右图的Z轴上。OpenGL中构建投影矩阵的函数是gluPerspective，其原型是：

void gluPerspective(
　　GLdouble fovy, //角度
　　GLdouble aspect,//视景体的宽高比
　　GLdouble zNear,//沿z轴方向的两裁面之间的距离的近处
　　GLdouble zFar //沿z轴方向的两裁面之间的距离的远处
　　)

向x轴负方向看过去，这两个变换空间如下图所示。

考虑上边界情况，对于眼坐标系来说：y=-ztan(fovy/2), 而对于投影坐标系来说:y'=1.

很显然,y的映射为:y'=-y/(ztan(fovy/2)=-ycot(fovy/2)/z。 考虑z相同的情况，此时y‘是y的单调增函数（考虑y>0情况)。所以这个映射也符合视锥内部的点。

对于y的负边界来说，情况类似。对于x来说，只需额外除以aspect这个值即可，因为：

aspect = width/height = x宽度/y高度

所以：x' = -xcot(fovy/2)/z/aspect.

由于这是非线性变换，要解决除以z，可以借助齐次坐标w来进行。令w'=-z，于是：

x' = xcot(fovy/2)/aspect.

y' = ycot(fovy/2)

w'=-z

那么z呢？其实在屏幕上显示后，z值对于颜色没有任何贡献。但是，消隐少不了z值，也就是所谓的z缓冲。对于z变换来说，一定要保留两个点的z值大小顺序不变，另外，为了归一化处理，将其映射到[-1,1]之间。而最后我们还可以在OpenGL中用glDepthRange调节深度值，这一点说明从这里算出来的Z值到最终的深度值之间还会有一个映射调整。好，话说远了，回来看看z的变换。

那么所需的Z变换有以下约束：

1)z=-znear--->z'=-1;

2)z=-zfar---->z'=+1;

3)映射具有单调减性质（因为有负号，所以单调减）

4)最终z'还要除以w'=-z

根据以上条件，构造映射：

f(z) = -z*（Zfar+Znear） / ( Zfar – Znear ) – 2* Zfar*Znear / ( Zfar – Znear )

z'= f(z)/w' = （Zfar+Znear） / ( Zfar – Znear ) + [2* Zfar*Znear / ( Zfar – Znear )]/z

验证一下，上述几个条件都满足。你要是想正向推导这个公式，也是可以的。可以参考后面的文献。

还是直接给出代码，一看就清楚：

static const double myPi = 3.14159265358979323846;

static void buildPerspectiveMatrix(double fieldOfView,
double aspectRatio,
double zNear, double zFar,
float m[16])
{
double sine, cotangent, deltaZ;
double radians = fieldOfView / 2.0 * myPi / 180.0;

deltaZ = zFar - zNear;
sine = sin(radians);
/* Should be non-zero to avoid division by zero. */
assert(deltaZ);
assert(sine);
assert(aspectRatio);
cotangent = cos(radians) / sine;

m[0*4+0] = cotangent / aspectRatio;
m[0*4+1] = 0.0;
m[0*4+2] = 0.0;
m[0*4+3] = 0.0;

m[1*4+0] = 0.0;
m[1*4+1] = cotangent;
m[1*4+2] = 0.0;
m[1*4+3] = 0.0;

m[2*4+0] = 0.0;
m[2*4+1] = 0.0;
m[2*4+2] = -(zFar + zNear) / deltaZ;
m[2*4+3] = -2 * zNear * zFar / deltaZ;

m[3*4+0] = 0.0;
m[3*4+1] = 0.0;
m[3*4+2] = -1;
m[3*4+3] = 0;
}

至此，几个矩阵都推导完成了。

题外话：其实网上有很介绍矩阵的推导过程，包括很多教材上也有。但要真正理解，只有一种方法：实践出真知，自己动手。

Reference：
1. http://game.chinaitlab.com/arithmetic/28004.html；
2.http://blog.csdn.net/popy007/archive/2007/09/23/1797121.aspx；
3.http://blog.csdn.net/popy007/archive/2009/04/19/4091967.aspx。

Inside CG transformation (II).

上一篇先解决了CG中的第一步变换，接下来来看viewMatrix。

视图变换实际上是坐标系的变换，将世界坐标系变换到眼坐标系（view ），这样做可以方便后面的投影(Projection)操作，因为投影是在眼坐标系中进行的。

建立眼坐标系用gluLookAt函数来实现，我们看看函数原型：
gluLookAt(double eyex, double eyey, double eyez, double centerx, double centery, double centerz, double upx, double upy, double upz）

这些参数中包括三种信息：眼睛的位置eye，视点中心位置center和向上的向量up。注意up不是y轴，只是指明了y轴的大致方向。

设眼睛的位置为原点，向量（eye-center）所得的向量为Z轴，y叉乘z得到x轴，再由x轴和z轴叉乘，反算出y轴。这就是眼坐标系的三个方向，设其单位向量为u、v、n。

然后就是坐标系之间的变换，这一点可以通过两种思路来完成：
1. 将眼坐标系变换到同世界坐标系重合。
2.代数的方法。

对于1而言，又有两种思路：
1. 代数上已经证明的公式：先移动eye的位置到世界坐标系的原点，再进行正交变换；
2. 由基本的变换合成，如下：
1）移动眼睛的位置到世界坐标系原点；
2）将眼坐标系绕世界坐标系x轴旋转，使得眼坐标系的z轴位于世界坐标系xoz平面；
3）将眼坐标系绕世界坐标系y轴旋转，使得眼坐标系的z轴和世界坐标系z轴重合；
4）将眼坐标系绕世界坐标系z轴旋转，使得眼坐标系的x、y轴和世界坐标系x、y轴重合；
这个步骤也需要很多的草稿纸，O(∩_∩)O

对于2，代数的方法，显得更简洁一些，看下图：


图中，对于任意的P点，向量减OP-Oeye= eyeP。再考虑eyeP向量，将其对眼坐标系的x/y/z轴进行投影，呵呵，就能得到P点在新坐标系中的坐标。投影很简单，用点乘即可。
于是，对于P点的新坐标P'（x',y',z'):
x' =(OP-Oeye).u = OP.u - Oeye.u =(x,y,z).(u[0],u[1],u[2]) -Oeye.(u[0],u[1],u[2]);
其中，u为眼坐标x轴方向的单位向量, x为P点的x坐标，u[0]为u向量的x值；其他类似。

y、z类似，我想你应该能写出矩阵来了吧，此处矩阵略过，直接看代码吧。

/* Build a row-major (C-style) 4x4 matrix transform based on the parameters for gluLookAt. */
static void buildLookAtMatrix(double eyex, double eyey, double eyez,
double centerx, double centery, double centerz,
double upx, double upy, double upz,float m[16])
{
double x[3], y[3], z[3], mag;

/* Difference eye and center vectors to make Z vector. */
z[0] = eyex - centerx;
z[1] = eyey - centery;
z[2] = eyez - centerz;
/* Normalize Z. */
mag = sqrt(z[0]*z[0] + z[1]*z[1] + z[2]*z[2]);
if (mag) {
z[0] /= mag;
z[1] /= mag;
z[2] /= mag;
}

/* Up vector makes Y vector. */
y[0] = upx;
y[1] = upy;
y[2] = upz;

/* X vector = Y cross Z. */
x[0] = y[1]*z[2] - y[2]*z[1];
x[1] = -y[0]*z[2] + y[2]*z[0];
x[2] = y[0]*z[1] - y[1]*z[0];

/* Recompute Y = Z cross X. */
y[0] = z[1]*x[2] - z[2]*x[1];
y[1] = -z[0]*x[2] + z[2]*x[0];
y[2] = z[0]*x[1] - z[1]*x[0];

/* Normalize X. */
mag = sqrt(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]);
if (mag) {
x[0] /= mag;
x[1] /= mag;
x[2] /= mag;
}

/* Normalize Y. */
mag = sqrt(y[0]*y[0] + y[1]*y[1] + y[2]*y[2]);
if (mag) {
y[0] /= mag;
y[1] /= mag;
y[2] /= mag;
}

/* Build resulting view matrix. */
m[0*4+0] = x[0]; m[0*4+1] = x[1];
m[0*4+2] = x[2]; m[0*4+3] = -x[0]*eyex + -x[1]*eyey + -x[2]*eyez;

m[1*4+0] = y[0]; m[1*4+1] = y[1];
m[1*4+2] = y[2]; m[1*4+3] = -y[0]*eyex + -y[1]*eyey + -y[2]*eyez;

m[2*4+0] = z[0]; m[2*4+1] = z[1];
m[2*4+2] = z[2]; m[2*4+3] = -z[0]*eyex + -z[1]*eyey + -z[2]*eyez;

m[3*4+0] = 0.0; m[3*4+1] = 0.0; m[3*4+2] = 0.0; m[3*4+3] = 1.0;
}

到此，view矩阵已经OK了。

Inside CG transformation (I).

对于学习计算机图形学的人来说，对其中的几个变换开始时往往有些难以理解。就算理解了，也是一知半解。不过，一般的图形学编程来说，不用深入这些变换照样能够编写程序。有些人学了4、5年图形学，但是还没有完整的推导或者实践一下这几个关键的变换，因为这耗时，而且得沉住气，但往往收获也是巨大的。

前几天，看到Cg的OpenGL examples中有一个“08_vertex_transform”例子，所以的这几个变换都有程序实现，我就决心来深入一下这几个变换。

先要明白两点：
首先，OpenGL采用列向量矩阵，所以几个相应矩阵按照左乘的方式进行，也就是：
最终的投影矩阵 modelViewProj = projectionMatrix * viewMatrix * modelMatrix。modelMatrix中平移、旋转缩放同样符合这种顺序。

第二，对于3D程序来说，一般情况下viewMaxtrix和projectMatrix只需设置一次即可, 产生动画通过变换modelMatrix达到。当然，如果你变换viewMaxtrix，也能达到动画的效果。在gl中采用两个函数gluLookAt和gluPerspective来设置这两个矩阵。

下面来看看这几个变换。

模型变换时最容易理解的，涉及旋转、平移和缩放。在OpenGL中的glRotate，glTranslate和glScale.目前，Cg的这个例子没有完成Scale的变换，我自己写的，很简单。先给出平移和缩放的矩阵构建代码。
static void makeTranslateMatrix(float x, float y, float z, float m[16])
{
m[0] = 1; m[1] = 0; m[2] = 0; m[3] = x;
m[4] = 0; m[5] = 1; m[6] = 0; m[7] = y;
m[8] = 0; m[9] = 0; m[10] = 1; m[11] = z;
m[12] = 0; m[13] = 0; m[14] = 0; m[15] = 1;
}

static void makeScaleMatrix(float ax, float ay, float az, float m[16])
{
m[0] = ax; m[1] = 0; m[2] = 0; m[3] = 0;
m[4] = 0; m[5] = ay; m[6] = 0; m[7] = 0;
m[8] = 0; m[9] = 0; m[10] = az; m[11] = 0;
m[12] = 0; m[13] = 0; m[14] = 0; m[15] = 1;
}

而旋转要复杂一点，因为要考虑任意向量V，当然这个向量是从原点出发的（要是任意轴旋转，轴的端点可以是任意两点）。5个步骤如下：
1. 绕x轴将向量V旋转a角度到xoz平面，记为Tx(a)；
2. 绕y轴将向量V旋转b角度到与x轴重合，记为Ty(b);
3. 将物体绕x轴（向量V）旋转θ角度，记为Tx(θ);
4. 2的逆过程，记为Ty(-b);
5. 1的逆过程，记为Tx(-a)；

为方便起见，将向量V归一化, 得到Vn（v0,v1,v2）。由此 可知
cosa = v2/sqrt(v1^2+v2^2);
cosb = v0;

再加这个Tx(a)变换，其他类似：
x′＝x
y′＝ycosθ－zsinθ
z′＝ysinθ＋zcosθ

这些信息都清楚了，何不手工算一下这个旋转矩阵呢？多准备一点草稿纸，肯定能得出最后结果。旋转变换矩阵有些麻烦，直接给代码：

static void makeRotateMatrix(float angle, float ax, float ay, float az,float m[16])
{
float radians, sine, cosine, ab, bc, ca, tx, ty, tz;
float axis[3];
float mag;

axis[0] = ax;
axis[1] = ay;
axis[2] = az;
mag = sqrt(axis[0]*axis[0] + axis[1]*axis[1] + axis[2]*axis[2]);
if (mag)
{
axis[0] /= mag;
axis[1] /= mag;
axis[2] /= mag;
}
//normalizing above

radians = angle * myPi / 180.0;
sine = sin(radians);
cosine = cos(radians);
ab = axis[0] * axis[1] * (1 - cosine);
bc = axis[1] * axis[2] * (1 - cosine);
ca = axis[2] * axis[0] * (1 - cosine);
tx = axis[0] * axis[0];
ty = axis[1] * axis[1];
tz = axis[2] * axis[2];

m[0] = tx + cosine * (1 - tx);
m[1] = ab + axis[2] * sine;
m[2] = ca - axis[1] * sine;
m[3] = 0.0f;

m[4] = ab - axis[2] * sine;
m[5] = ty + cosine * (1 - ty);
m[6] = bc + axis[0] * sine;
m[7] = 0.0f;

m[8] = ca + axis[1] * sine;
m[9] = bc - axis[0] * sine;
m[10] = tz + cosine * (1 - tz);
m[11] = 0;

m[12] = 0;
m[13] = 0;
m[14] = 0;
m[15] = 1;
}

至此，模型变换已完成。

Math Basis

点乘 dot product

点乘，也叫向量的内积、数量积。顾名思义，求下来的结果是一个数。
向量a·向量b=abcosθ
在物理学中，已知力与位移求功，实际上就是求向量F与向量s的内积，即要用点乘。
将向量用坐标表示（三维向量），
若向量a=(a1,b1,c1)，向量b=(a2,b2,c2)，
则向量a·向量b=a1a2+b1b2+c1c2

叉乘 cross product

叉乘，也叫向量的外积、向量积。顾名思义，求下来的结果是一个向量，记这个向量为c。
|c| = |axb|=|a||b|sinθ
向量c的方向与a,b所在的平面垂直，且方向要用“右手法则”判断（用右手的四指先表示向量a的方向，然后手指朝着手心的方向摆动到向量b的方向，大拇指所指的方向就是向量c的方向）。因此向量的外积不遵守乘法交换率，因为向量a×向量b= -向量b×向量a。在物理学中，已知力与力臂求力矩，就是向量的外积，即叉乘。
将向量用坐标表示（三维向量），
若向量a=(a1,b1,c1)，向量b=(a2,b2,c2)，
则 向量a×向量b=
i j k
a1 b1 c1
a2 b2 c2
=(b1c2-b2c1,c1a2-a1c2,a1b2-a2b1)

三角函数公式

1.诱导公式
sin(-a)=-sin(a)
cos(-a)=cos(a)

2.两角和与差的三角函数
sin(a+b)=sin(a)cos(b)+cos(α)sin(b)
cos(a+b)=cos(a)cos(b)-sin(a)sin(b)
sin(a-b)=sin(a)cos(b)-cos(a)sin(b)
cos(a-b)=cos(a)cos(b)+sin(a)sin(b)
tan(a+b)=[tan(a)+tan(b)]/[1-tan(a)tan(b)]
tan(a-b)=[tan(a)-tan(b)]/[1+tan(a)tan(b)]

3.和差化积公式
sin(a)+sin(b)=2sin((a+b)/2)cos((a-b)/2)
sin(a)−sin(b)=2cos((a+b)/2)sin((a-b)/2)
cos(a)+cos(b)=2cos((a+b)/2)cos((a-b)/2)
cos(a)-cos(b)=-2sin((a+b)/2)sin((a-b)/2)

4.二倍角公式
sin(2a)=2sin(a)cos(b)
cos(2a)=cos^2(a)-sin^2(a)=2cos^2(a)-1=1-2sin^2(a)

5. 三角恒等式
sin^2θ+cos^2θ=1；
1+tan^2θ=sec^2θ；
1+cot^2θ=csc^2θ

Study of photography

扫盲一下照相相关的几个术语，这个解释很到位。

光圈

光圈的功能就如同我们人类眼睛的虹蟆，主要用来调整数码相机的进光量，一般以f/2、F2、1：2来表示，举例来说：f/2表示光圈的大小为镜头直径的1/2 ，而f/8 则表示光圈为镜头直径的1/8而已，所以较小的f值表示较大的光圈，一般镜头上的标示会以数码相机的最大光圈值表示，变焦镜头若有显示2个f值，则表示此相机最大及最小的光圈。光圈除了可控制光线的明暗外，对于图像的景深也是会有影响。景深是影响拍摄主体与背景之间清晰程度的关键，也可以说是图像的锐利度，大光圈时拍摄出来的图像锐利度较小，小光圈拍摄则较大。从数码上而言，光圈数值越小表示光圈越大，也代表透光的孔径大、透光量大，不论要拍摄快速移动的物体或在昏暗的空间拍摄，都很方便。而且光圈也决定了画面的景深(锐利度)，如果是设定为大光圈，那么画面中除了主题清晰，其它景物都会呈现模糊、柔美的感觉。时下多数数码相机的光圈值最大都在2.8左右。

景深

在进行拍摄时，调节相机镜头，使距离相机一定距离的景物清晰成像的过程，叫做对焦，那个景物所在的点，称为对焦点，因为"清晰"并不是一种绝对的概念，所以，对焦点前（靠近相机）、后一定距离内的景物的成像都可以是清晰的，这个前后范围的总和，就叫做景深，意思是只要在这个范围之内的景物，都能清楚地拍摄到。景深的大小，首先与镜头焦距有关，焦距长的镜头，景深小，焦距短的镜头景深大。其次，景深与光圈有关，光圈越小（数值越大，例如f16的光圈比f11的光圈小），景深就越大；光圈越大（数值越小，例如f2.8的光圈大于f5.6）景深就越小。其次，前景深小于后后景深，也就是说，精确对焦之后，对焦点前面只有很短一点距离内的景物能清晰成像，而对焦点后面很长一段距离内的景物，都是清晰的。

F值

就是光圈值
光圈是一个用来控制光线透过镜头，进入机身内感光面的光量的装置，它通常是在镜头内。表达光圈大小我们是用f值。

光圈f值=镜头的焦距/镜头口径的直径

单反相机

单反就是指单镜头反光,即SLR（Single Lens Reflex）。在这种系统中，反光镜和棱镜的独到设计使得摄影者可以从取景器中直接观察到通过镜头的影像。单镜头反光照相机的构造图中可以看到，光线透过镜头到达反光镜后，折射到上面的对焦屏并结成影像，透过接目镜和五棱镜，我们可以在观景窗中看到外面的景物。拍摄时，当按下快门钮，反光镜便会往上弹起，软片前面的快门幕帘便同时打开，通过镜头的光线（影像）便投影到软片上使胶片感光，尔后反光镜便立即恢复原状，观景窗中再次可以看到影像。单镜头反光相机的这种构造，确定了它是完全透过镜头对焦拍摄的，它能使观景窗中所看到的影像和胶片上永远一样，它的取景范围和实际拍摄范围基本上一致，消除了旁轴平视取景照相机的视差现象，从学习摄影的角度来看，十分有利于直观地取景构图。 单镜头反光相机还有一个很大的特点就是可以交换不同规格的镜头。

Radiosity Rendering Algorithm.

辐射度算法是一种全局光照算法。基本思路如下：
1. 首先，根据外部光源，得到一幅直接光照生成的光照图。
2. 然后，以此光照图+外部光源为光源，对场景进行辐射。
3. 迭代2过程。

Hugo Elias 的解释非常形象，以场景中的每个面片看到的“场景”为光源，计算此面片的亮度。

算法伪代码：

load scene

divide each surface into roughly equal sized patches

initialise_patches:
for each Patch in the scene
if this patch is a light then
patch.emmision = some amount of light
else
patch.emmision = black
end if
patch.excident = patch.emmision
end Patch loop

Passes_Loop:each patch collects light from the scene
for each Patch in the scene
render the scene from the point of view of this patch
patch.incident = sum of incident light in rendering
end Patch loop

calculate excident light from each patch:
for each Patch in the scene
I = patch.incident
R = patch.reflectance
E = patch.emmision
patch.excident = (I*R) + E
end Patch loop

Have we done enough passes?
if not then goto Passes_Loop

代码解释

initialize patches:(初始化面片)
一开始，所有的面片都是黑的，除了那些能自身辐射出光线的面片。因此，那些能辐射光线的面片的出射光强的初始值应被初始化伪它的辐射光强。其他面片的辐射光强都应为0。

Passes Loop(遍历循环):
代码多次重复这个循环直到场景有了可接受的光照效果。每次循环之后，也就多模拟了一次光在场景中的反射。

each patch collects light from the scene(每个面片从场景中收集光强)
如果我之前在文章中解释的那样，每个面片都被它能够看见的其他面片照亮。要达到这个目的，可以简单地把把场景渲染到面片的视角，然后把它所看的光强相加。我将在下一小节更详细地解释这一步。

calculate excident light from each patch(为每个面片计算出射光强):
计算出有多少光强抵达面片之后，我们现在可以计算出有多少光强离开面片(被反射)。
这个过程必须被循环多次以达到一个好的效果。如果渲染器还需要一个循环，我们就调转到标记"Passes Loop"。

Reference:
1. English:http://freespace.virgin.net/hugo.elias/radiosity/radiosity.htm
2. Chinese:http://dev.gameres.com/Program/Visual/3D/Radiosity_Translation.htm

Illumination Tech of CG.

Direct Illumination Problems and Advantages

Direct Illumination is a term that covers the principal lighting methods used by old school rendering engines such as 3D Studio and POV. A scene consists of two types of entity: Objects and Lights. Lights cast light onto Objects, unless there is another Object in the way, in which case a shadow is left behind.

Ray Tracing:
- Can render both mathematically described objects and polygons
- Allows you to do some cool volumetric effects
- Slow
- Very sharp shadows and reflections

Shadow Volumes:
- Can be modified to render soft shadows (very tricky)
- Tricky to implement
- Very sharp shadows
- Polygons only

Z-Buffer:
- Easy to implement
- Fast (real-time)
- Sharp shadows with aliasing problems


Global Illumination Problems and Advantages

Images produced by global illumination methods can look very convincing indeed; in a league of their own, leaving old skool renderers to churn out sad cartoons. But, and it's a big 'but': 'BUT!' they are slower. Just as once you may have left your ray tracer all day, and come back to be thrilled by the image it produced, you will be doing the same here.

Radiosity:
- Very realistic lighting for diffuse surfaces
- Conceptually simple and easy to implement
- Easy to optimise with 3D hardware
- Slow
- Does not handle point sources well - nor shiny surfaces
- Always over complicated and poorly explained in books

Monte Carlo Method:
- Very, very good results.
- Can simulate pretty well any optical phenomenon
- Slow
- Slightly difficult
- Requires some cleverness to optimise
- Always over complicated and poorly explained in books

Reference:http://freespace.virgin.net/hugo.elias/radiosity/radiosity.htm

Study of Computer Graphic.

Reference: 3D Games - Real time rendering and software technology Volume I.

1. 仿射变换(affine transformations):旋转、缩放、平移、错切变换；

2. 法线：向量叉积；
向量夹角：点积(V.W>0,θ<90;v.w=0,θ=90...)；
投影：点积（v单位向量，w对v进行投影：|x|=v.w）。
3. 双线性插值(bi-linear interpolation):顶点-->边；边-->多边形内部；
4. 模型表示方法：

• polygonal--ploygon or triangle(Problem:continuous LOD is hard, avoiding popping).
• bi-cubic parametric patches(双三次曲面）= curved quadrilaterals(曲面四边形）(Problem:smoothness between quadrilaterals,inapposite for complex object. BUT easy for LOD).
• constructive solid geometry,CSG(intersection,union,subtraction).
• spatial subdivision techniques(空间细分技术）.
• implicit representation(隐式表示）.

 5. Economy of mesh represatation

• Triangle chains
• LOD

6.Terrain modeling

• fractal geometry(分形几何学)
• terrain LOD: triangle bintree(三角二叉树）

7. Curve and curved face

• Bezier curve:
• 平行六面体
• 四个control point Pi(i=0,1,2,3);
• 基函数Bi为(1-u)^3,3u(1-u)^2,3u^2(1-u),u^3; 曲线Q(u)= ∑PiBi.
• Disadvantage: global effection when manipulating one control point; smoothness of multi Bezier connection;
• B样条：任意个控制点，任意4个一组；内在连续性；局部性（任意个控制点改变，影响4个曲线段）。
• 非均匀有理B样条（NURBS）

8. Rendering

• 视见体裁剪；
• 局部反射模型：计算顶点的光强；明暗算法：由定点得到面片的每一点的光强；
• 局部反射模型：Phong—反射光=环境光+漫反射+镜面反射，这个是Phong光照模型；
• 明暗算法：
• Gouraud—Ip=interpolation(Ia,Ib),对顶点a，b的亮度进行插值，不会有高光出现;
• Phong—对法向插值，得到插值点的Np，再求亮度。这个是Phong着色模型。
• 一般Gouraud球漫反射分量，Phong求镜面发射分量。
• 隐面剔除：Z缓冲—在搜索的同时，将Z值最小的像素写入Frame Buffer；需要x*y*n大小。

9. texture mapping

• function
• 1)common color of pixel: color of texture与局部反射模型计算得到的漫反射系数相乘；
• 2)specular color: 进行环境映射贴图，避免完全光线跟踪；沿着反射后的视见向量，在场景中寻找纹理，可倒映出环境中有光泽的物体；
• 3)凹凸纹理（bump mapping）: normal vector perturbation;
• 4)transparent:控制透明物体的不透明程度；
• bitmap-->planar, cylinder, sphere-->object

Study about Fourier Transform

Reference：http://blog.csdn.net/dznlong/archive/2008/04/08/2261150%20.aspx。

1.傅里叶变换是一种逼近表示，而非精确表示；
2.选用正弦、余弦波形表示是否唯一？不是，只是为了方便。
3. 傅里叶变换的四种分类：
非周期性连续信号 --- 傅立叶变换（Fourier Transform）
周期性连续信号 --- 傅立叶级数(Fourier Series)
非周期性离散信号 --- 离散时域傅立叶变换（Discrete Time Fourier Transform）
周期性离散信号 --- 离散傅立叶变换(Discrete Fourier Transform)
4. 傅立叶变换是针对正无穷大和负无穷大的信号，即信号的的长度是无穷大的. 非无穷怎么办？扩展即可。
5. 傅里叶变换实质：将信号分解为若干个离散的、不同频率的正弦、余弦信号分量。