LED curtain.

Barbara Bona --- http://www.barbarabona.com/

http://www.cnbeta.com/articles/110690.htm

Simulated windows.

每天早上醒来看到的窗外的景色都是一成不变的画面，会让人视觉疲劳的，如果窗外美景可以DIY的话，相信大家都愿意尝试。近日一位名叫Ryanhoagland的DIY高手就自己制作了一个仿真度超高的电子窗户.

http://www.cnbeta.com/articles/108797.htm



其实之前有一个建筑学院的老师跟我讨论了这种建筑的可能性，不过他的想法更大，就是一面墙整个都是虚拟场景，比如海滩等。看来，这个还是很有前途的。

A glance of Siggraph 2009 papers.

SIGGRAPH向来是Computer Graphics研究的风向标。今天粗略过了一遍SigGRAPH2009论文，乱弹一下心得：
1.hair modeling and simulation and rendering 主题真是耐做，有两篇都是关于此主题的：“Example-based hair geometry synthesis ”和“Detail preserving continuum simulation of straight hair ”。

2.”Fluid simulation“研究同样是常青树，专门有一个Session。
3.国内的强人“Kun Zhou“还是发了3篇论文。

4.最”牛气“的论文是”Semantic Deformation Transfer“，名字很短，正文只有6页，12篇参考文献；

5.有两个关于”Character animation”的session，8篇论文，占到总论文的10%；

6.有两个跟physics相关的session，“Physically based modeling: from contact to capture”和“Reduced physics for animation”，拓展一点，如果把“Creating natural variations”自然也看出跟物理相关的话，总共有三个session。物理是王道呀。

7. simulation方面，除了“Fluid simulation”session，还有“Directable, high-resolution simulation of fire on the GPU“，“Detail preserving continuum simulation of straight hair ”和“Interactive simulation of surgical needle insertion and steering ”。

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θ

Cg studying.......

1. GPU Pipeline.



2. The Transform


Modeling transfermation: Object Space --> World Space.It makes the model change its shape;

Viewing transfermation: World Space -->Eye Space . It puts eyes at the original.

Project transfermation: Eyes Space --> Clip(Project) Space. Changing View frustum to a CUBIC space, in which it is more convenient for clipping.


The Modelview Matrix
Most lighting and other shading computations involve quantities such as positions and surface normals. In general, these computations tend to be more efficient when performed in either eye space or object space. World space is useful in your application for establishing the overall spatial relationships between objects in a scene, but it is not particularly efficient for lighting and other shading computations.
For this reason, we typically combine the two matrices that represent the modeling and view transforms into a single matrix known as the modelview matrix. You can combine the two matrices by simply multiplying the view matrix by the modeling matrix

Computer Graphics Practical.

1. Point
CPoint(n)
2. Mesh
Quad or Triangle construction from Points
3. FillMesh
Using polygon filling algorithm to determine each point in polygon
4. Calculate each point's color with light & texture

For vertexs
Ambient light: Ambient Light Intensity * Ambient Reflection Coefficient * RGB
Diffuse light: Diffuse Light Intensity * Diffuse Reflection Coefficient *cos(θ)* RGB
-->L vector = Light Position - Point ; N vector = Point Normal;
--> cos(θ) = L · N/|L||N|

For other p0ints
Gourand: Each vertex color of the mesh -> BilinearInterpolation ->This Point Color
Phong: Each vertex normal of the mesh -> BilinearInterpolation -> This Point Normal -> Color

Problem of OnReceive of CSocket

以前用过CSocket这个类，主要是考虑可以用CSocketFile进行串行化数据，比较方便。这次写的时候竟然出了几个问题：

1. 网上下载了chatter例程，用VC6.0编译后运行，没有任何问题；切换到VS2008下，编译没有问题，运行时，发现OnReceive时而触发，时而不触发。到网上查了查，同样问题的人很多，看来是个老大难了。后来加上OnReceive结束之后，调用 AsyncSelect( FD_READFD_WRITE )，解决了此问题。几乎弄了一整天，跟http://www.ecjtu.org/forum/read.php?tid-19062.html这位仁兄情况几乎相同；

2. 其实，刚开始我就发现了这个解决方案，但是我一直觉得OnReceive函数由系统框架来调用，不应该需要施加“外力”影响。国外很多人说是OnReceive函数里面调用了多次Receive引发的。如果真是这样，另外一个问题就是CArchieve的load动作跟CSocket的Receive关系是什么？

3. CArchieve的Read和Write函数可以读写缓冲区，我原来一直用“>>”和“<<"符号一个一个地写DWORD，晕，效率很低。

结论：好像用CSocketFile来发送、接受数据有点得不偿失，问题很多，实际也没有带来很多好处。不如直接用CAsyncSocket类呢。

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

0x5f3759df, Magic?

昨天无意看到网上有一个开平方根的快速算法，其实是很老的一个帖子，大约在2000年出现的：
---------------------------------------------------
有人在Quake III的源代码里面发现这么一段用来求平方根的代码:

/*================SquareRootFloat================*/

float SquareRootFloat(float number) {
long i;
float x, y;
const float f = 1.5F;
x = number * 0.5F;
y = number;
i = * ( long * ) &y;
i = 0x5f3759df - ( i >> 1 ); //注意这一行
y = * ( float * ) &i;
y = y * ( f - ( x * y * y ) );
y = y * ( f - ( x * y * y ) );
return number * y;
}

0x5f3759df？ 这是个什么东西？ 学过数值分析就知道，算法里面求平方根一般采用
的是无限逼近的方法，比如牛顿迭代法，抱歉当年我数值分析学的太烂，也讲不清楚
。简单来说比如求5的平方根，选一个猜测值比如2，那么我们可以这么算

5/2 = 2.5; 2.5+2/2 = 2.25; 5/2.25 = xxx; 2.25+xxx/2 = xxxx ...
这样反复迭代下去，结果必定收敛于sqrt(5)，没错，一般的求平方根都是这么算的
。而卡马克的不同之处在于，他选择了一个神秘的猜测值0x5f3759df作为起始，使得
整个逼近过程收敛速度暴涨，对于Quake III所要求的精度10的负三次方，只需要一
次迭代就能够得到结果。

好吧，如果这还不算牛b，接着看。

普渡大学的数学家Chris Lomont看了以后觉得有趣，决定要研究一下卡马克弄出来的
这个猜测值有什么奥秘。Lomont也是个牛人，在精心研究之后从理论上也推导出一个
最佳猜测值，和卡马克的数字非常接近, 0x5f37642f。卡马克真牛，他是外星人吗？


传奇并没有在这里结束。Lomont计算出结果以后非常满意，于是拿自己计算出的起始
值和卡马克的神秘数字做比赛，看看谁的数字能够更快更精确的求得平方根。结果是
卡马克赢了... 谁也不知道卡马克是怎么找到这个数字的。

最后Lomont怒了，采用暴力方法一个数字一个数字试过来，终于找到一个比卡马克数
字要好上那么一丁点的数字，虽然实际上这两个数字所产生的结果非常近似，这个暴
力得出的数字是0x5f375a86。

Lomont为此写下一篇论文，"Fast Inverse Square Root"。
-------------------------------------------------

于是我写了一段程序来验证，程序参考了Lomont的程序，http://www.lomont.org/Math/Papers/2003/InvSqrtCode.zip。但是号称快四倍的神秘算法，我一直未能得到很多的结果，程序表明两个算法性能差不多。奇怪中！

---------------------------------------------------------
float InvSqrtExact(float x)
{
return (float)(1.0/sqrt(x));
} // InvSqrtExact

/// "Carmack" version of fast inverse sqrt,
/// based on Newton method, 1 iteration
float InvSqrt(float x)
{
float xhalf = 0.5f*x;
int i = *(int*)&x;
i = 0x5f3759df - (i>>1); // hidden initial guess, fast
x = *(float*)&i;
x = x*(1.5f-xhalf*x*x);
// x = x*(1.5f-xhalf*x*x); // add this in for added precision, or many more...
return x;
} // InvSqrt

/// Chris Lomont version of fast inverse sqrt,
/// based on Newton method, 1 iteration, more accurate
float InvSqrt_Lomont(float x)
{
float xhalf = 0.5f*x;
int i = *(int*)&x;
i = 0x5f375a86 - (i>>1); // hidden initial guess, fast - LOMONT
x = *(float*)&i;
x = x*(1.5f-xhalf*x*x);
// x = x*(1.5f-xhalf*x*x); // add this in for added precision, or many more...
return x;
} // InvSqrt_Lomont


/// given a function, a number of passes, and a value to test, do it
void TimeFunction(float (*func)(float), int passes, float v)
{
unsigned long used;
used = timeGetTime();
for (int pos = 0; pos < passes; pos++)
func(v);
used = timeGetTime() - used;

cout << used << " ms used for " << passes << " passes, avg ";
cout << (double)used/passes << " ms\n";
} // TimeFunction

int main(void)
{
// timing of versions
int passes = 10000000;
float v = 1.234f; // test value
cout << "Timing basic function\n";
TimeFunction(InvSqrtExact,passes,v);
cout << "Timing Carmack function\n";
TimeFunction(InvSqrt,passes,v);
cout << "Timing Lomont function\n";
TimeFunction(InvSqrt_Lomont,passes,v);

cout << "Done. By Chris Lomont 2003\n";
return 0;
} // main

---------------------------------------------------------

运算结果相差无几，难道我的程序有什么问题？

继续探索，发现同样的程序在VC6.0和VC20008下面竟然表现不同，怪事！

VC++.net 2008
-------------------------------------
Timing basic function
2187 ms used for 10000000 passes, avg 0.0002187 ms
Timing Carmack function
828 ms used for 10000000 passes, avg 8.28e-005 ms
Timing Lomont function
844 ms used for 10000000 passes, avg 8.44e-005 ms
Done. By Chris Lomont 2003

VC++ 6.0
--------------------------------------
Timing basic function
625 ms used for 10000000 passes, avg 6.25e-005 ms
Timing Carmack function
453 ms used for 10000000 passes, avg 4.53e-005 ms
Timing Lomont function
453 ms used for 10000000 passes, avg 4.53e-005 ms
Done. By Chris Lomont 2003
------------------------------------

貌似跟“using namespace std”有一定的关系，可能是sqrt的函数库版本不太一样，因为我是用VC++6.0，崩溃中......

看来0x5f3759df也不是那么神奇吧，另外，也可以看出.net程序的性能确实不怎样。

Taiwan problem.

今天看到新闻，“胡锦涛：必须更严厉警告台当局”和“胡锦涛会见布什 强调今明两年为台海局势高危期”，感觉到这好像是中央对台湾问题最严重的措辞了。或许是因为奥运会，台湾肆意妄为，估计是大陆在奥运前不会发动军事行动。胡主席这话将得很明白，台湾请自重。呵呵。下面是sohu网友的评论，觉得比较有道理：
“开战的日期是08奥运会后２００９年春天，目前解放军要开始做准备了，地方办奥运，军方搞战备。用奥运来强化外交，增进与世界各国友谊，赢得世界广泛支持，为战争做好外部环境准备。但现在还需一定的克制。光凭一时之勇是不行的。我认为台湾目前这样做是受了美国、日本等国的暗中指示，目的就是利用台湾让中国出手千方百计地破坏奥运会的顺利召开，如中国不够冷静，采取过激行动，势必会减弱中国的国际影响，损坏中国的国际形象，特别是让千幸万苦争取来的奥运会有可能导致流产，让中国失去一次提升国力的大好机会。所以目前要尽量采取一争和平方式，确保０８年奥运会顺利召开。”

New Chinese Phrases and Words of 2006 .

教育部今日在其官方网站发布了《中国语言生活状况报告(2006)》，列出了171条汉语新词语选目， 以下为今年公布的汉语新词语选目：
　　B
　　八荣八耻、白奴、白托、白银书、半糖夫妻、抱抱团、抱抱装、奔奔族、笔替、博斗、博客话剧、博文
　　C
　　草根网民、车奴、成考移民、城市依赖症、村证房
　　D
　　大肚子经济、倒扁、盗版党、等额配比基金、电话门、电子环保亭、吊瓶族、丁宠家庭、动能车、冻容、独二代、断背
　　E
　　★EMBA、2时歇业令、二奶专家
　　F
　　法商、返券黄牛、饭替、房魔、房奴、飞鱼族、废统、沸腾可乐、分手代理、福利腐败、祝寿螺患者、复古学堂
　　G
　　感恩红包、高薪跳蚤、搞笑、公司驻虫、谷歌、骨性、国际高考移民、国六条、国十条
　　H
　　海缆断网、海绵路、海啸音、寒促、汉芯造假事件、合吃族、红楼选秀、红衫军、换客、灰色技能、回购地、会议大使、婚嫁大年
　　J
　　急婚族、监控门、江选、奖骚扰、交强险、脚环鸡、节奴、解说门、禁电、掘客
　　K
　　卡神、考霸、科研包工头、啃椅族、★空调、骷髅门
　　L
　　赖校族、乐活族、梨花体、李娅空翻、两会博客、垄奴、绿色产房、裸考、裸替
　　M
　　M型社会、慢活族、美丽垃圾、迷卡、蜜月保姆、秒杀、明星枪手、墓产经济、墓奴
　　N
　　年后饭、暖巢管家、诺亚规则
　　P
　　跑酷、陪拼族、捧车族、拼卡、拼客
　　Q
　　7时代、擒人节、轻熟女、穷人跑、求学房、群租
　　R
　　让票区、人球、人户育婴师、润滑经济
　　S
　　三失、三手病、三限房、三支一扶、★晒、晒客、上海社保基金案、剩女、十五细则、世宗、试药族、手机幻听症、手机手、熟年、睡眠博客
　　T
　　她经济、痛快吧、图书漂移、土腐败、囤房捂盘、托业
　　W
　　网络晒衣族、威客、微笑北京、微笑圈、文化低保、文替、巫毒娃娃、捂盘惜售
　　X
　　新明星学者、新中间阶层、熊猫烧香、炫富、学术超男、学业预警
　　Y
　　压洲、洋腐败、洋漂族、医闹、医诉、移动商街、印客、游贿、游戏手
　　Z
　　砸票、择校税、证奴、职粉、终统、众包、住房痛苦指数、装嫩族、作弊克

看了一下，下出一身冷汗，因为只有红色的部分我知道一些，其他的很多都没有听说过呀。不会是我落伍了吧。你呢？如果我再去考语文，岂不是不及格？:-)。
还有，谷歌竟然列进去了，不知“百度”或者其他公司词汇为啥没有？Google强呀。

Chinese Cooperation.

联想收购Packard Bell成功。Packard Bell为欧洲第三大消费PC厂商，其总裁许立信是华人。前段时间，联想、宏基和许立信这三个Chinese之间的三角恋终于告一段落。当然，我也喜欢这个结果，毕竟宏基是台湾的，还是有些不同。呵呵。

想到之前的“联想收购IBM PC宣告成功”，联想是好样的！而国内的华为、中兴也在美国和欧洲不停的攻城略地，赞一个。

A Poem.

伟大的G-F-W，搞得无数Blogger彻夜难眠。Blogger昏睡至今，望眼欲穿。幸亏有CDN，不过硬是要在Blog后面加个尾部才可以访问（http://leepyzh.blogspot.com.nyud.net:8090 ）。我就是很奇怪，为什么别的Blog不封，如MS 的Live Space，就搞Google？Google的不作恶的信条，为什么找来了作恶的后果？世道呀......作诗一首，以解郁闷。

春晓

春眠没有晓，
夏天昏睡了。
长城挥手间，
博客死多少。