正文


至少,目前,工程论文想要做一个Material Editor。

当然,对于别人商用Engine的Editor只有膜拜的份,也不指望咱这小身板子能搞出一个那么厉害的东西,不过至少模仿一下吧。

大概看了一下Robert L. Cook的那篇《Shade Tree》,结合曾经对RenderMonkey为数不多的使用以及“只见过猪跑”的Editor,大致有了一些认识,这里先记下来。

暂且不说各种Shader效果,这些东西准备最终定下来之后逐步学习,大致浏览了一下相关资料,内容还是很多的,特大一壶,以至于我有点担心能否喝的完……

这一篇不会很长,只说明一下Shade Tree到底是个什么东西。

Shade Tree

还是先上一张图,看着图好说话,该图片展示的是U3D的Strumpy Shader Editor:

U3D

先来说明一下图表示的意义:

每一个方框框代表一个“过程”节点,框的右边是输入,左边是输出,例如图中的 TexCUBE 节点,这个过程接收一个名为 Sampler 与 Normal 的参数,输出两个颜色 RGBA 和 A。(这里的输入输出原来左右写反了,原谅一下左撇子)

Master部分是最终结果,基本每个Shader都会有。

其实我个人更倾向于 TexCUBE 接受了两个“过程”作为参数,而并非只是两个运算结果,因为Sampler可以有2D,3D,Cube等多种类型;而Normal,则根据情况可以有Reflection,Refraction等类型;而对于一个Surafec来讲,法线作为其自身属性,单独分离出去给其他的“过程”并非一个好的方案。

我对上面的整个过程的理解如下(当然有可能是错误的,因为图片是随意找的……):

使用一个 SimpleWorldReflection 来计算反射(个人感觉应该是折射),纹理采样使用Cube采样来对Surface进行处理,得到Cube纹理的运算结果,然后这个运算结果呗传入到Master中,与Albedo和Alpha量来得到最终的结果。

追加:我找到了与这张图片相关的描述:第966楼。不想点链接的可以看下面,我始终以为Reflection有点不太合适,不过CG上要达到一个看起来一样或者差别很小的效果可以有多种手段,所以不必在此钻牛角,当然,我自己Shader玩的不好也有很大原因。

What you want to do (to start with) is sample the cube from the reflection vector of the
surface using the ‘Vertex Reflection’ node. In this example I have also made the alpha be
based on an alpha from a texture and set the blend mode to be srcalpha /
oneminussourcealpha (so you can configure the alpha from the texture), and also given the
windows a ‘blue’ hint in the albedo. It’s nothing amazing, but it’s a decent start for
some glass. You will need to set the queue to transparent if you want to use this shader
also.

如果把上面图中的Flow换个方式:

U3D Simple

再来看看《Shade Tree》中的图:

Shade tree for copper

是不是很相似?

Cook文中的图看起来更像是编译原理中的语法树,但其实是无妨的:一个的N元operation本身就可以抽象为一个接受N个参数的processor。

万幸自己最终还是选修了编译,虽然学的不算好,理论的东西终究没剩下多少,都还回去了。本科专业没这门课导致我很长时间对很多东西理解困难,包括这里的Shade Tree。

求解final color的过程实际上就是一个bottom-up的过程,这就不用多做解释了。

基本上Shade Tree的本质就是这些了,当然对于具体的效果,比如光照,大气,雾,阴影等效果需要具体的工程实现,这些会随着自己的逐渐深入慢慢写出来。

实现?

实现倒是有那么一点点的小想法,但总感觉不太对,后面要找点相关资料才行,先说说想法吧:

首先都要将一些基本的过程给拆开,成为比如前面的Reflection等一个个单独的单元,对于这些单独的单元有两种方式组合起来:

  • shade编译之后有链接的步骤,那么可以首先将需要的各个单元编译好,然后在编辑结束后将对应的单元连接起来。但是分的太细可能会导致片元过多,不太好。

  • 直接在编辑结束的时候生成一个大的Shader进行编译链接。不过这种方式可能会出现指令条数的超出限制。新的硬件不太清楚,但旧硬件是一定会出问题的。

或者将两者折中一下?比如所有的Sampler可以放到一个文件中编译,最后作为一个库连接过去。

结束


话说最近的一次写Shader代码都快一年了,况且本来就一般般……不过还是先笔记想法慢慢记:不怕慢差,就怕不干嘛。

正文


过年把以前写的那个自娱自乐的破烂Software Renderer又写了一遍,加了一些新的Features。代码依旧稀烂,以至于拿出去给人看都嫌丢脸——虽然差不多算是明了设计模式到底是个什么玩意儿,但这跟能写出好的代码完全是两码事情——不过好在是代码虽然是一坨,但还是能跑起来的。

这里对比一下,姑且称之为版本一和版本二吧。

数学

两次的数学运算基本一致,版本二由于添加了近裁剪面裁剪,因此增加了对三角形进行切割的方法。

数学的内容基本都在前面图形相关的文里面谈到了,当然碰撞检测以及物理相关都没有涉及(较早的时候有参考过《Real Time Collision Detection》写过一些碰撞检测的代码,但时间太久,换电脑把代码倒腾丢了,具体内容也忘得差不多了,基本上就是OBB,AABB,Capsule,Sphere,Plane,Triangle,Ray,Segment之间的各种关系判定,求解特定交点等等)。

两个版本都不包含四元数相关运算。

总览

我们先按下光栅化部分关于Z-Buffer的话题,把它稍微放后一些。

版本一是真的写死了,只能贴一张纹理,也不支持光照,说白了整个Renderer就是靠一个光栅化函数撑起来的(drawTriangleWithTexture),当然不使用纹理也可以。不包含三角形裁剪与剔除,没有对每次提交的三角形列表进行任何修改,将三角形的顶点信息放进绘制列表之后直接就进行World->View->Perspective->ViewPort变换,然后交给光栅化函数进行绘制处理。由于不裁剪,因此穿过nearZ的坐标会被错误投影,必须保证要绘制的物体的顶点坐标都在nearZ之前。

版本二比版本稍微强了一丁点,但是还不够好。第二个版本加了近裁剪平面裁剪的功能,穿过近裁剪平面的三角形会被裁剪,位于裁剪面前面的三角形会传递到后续的三角形列表中,这里使用了两个vector来完成这个工作,但不好,一个Queue就行。裁剪是在View-Space中进行的,意味着必须先将World-Sapce中的顶点进行一次View-Space变换;光照则是直接在World-Space中进行,计算出来光对顶点颜色产生的影响,添加到顶点颜色中(只是用了Diffuse,也就是Normal与LightDir的夹角的cos值,整个过程属于Vertex光照而非Phong光照);裁剪后的三角形经过投影和视口变换,最后被光栅化。

版本一的顶点缓冲处理只有一种方式,钉死的:

1
2
3
4
5
6
struct Vertex
{
float x,y,z;
DWORD color;
float u,v;
};

版本二则稍微灵活了一点,模仿了GL,用了VertexArray、NormalArray、ColorArray、TexCoordArray*2来进行处理,至于过程,咱这写得自然比不上GL。但说老实话,现在觉得这不算是个好方案:至少我写的代码里面,每个绘制批次提交的数据最少两块(即便只有VertexArray的提交,但最终后台还是给加了个ColorArray,这个确实不好),多则四块,数据的局部性不能保证,而且不同属性的Array的Stride也不一样,写代码的时候一不注意就给了Bug……

下次准备用VertexDeclaration,如果还写的话。

从前面的描述中也能看出来版本二增加了一个Texture纹理,这样就可以对一个Surface使用两个2D纹理了,纹理之间的颜色运算是写死的Add,没有使用其他的。

纹理采样没有什么好说的,版本一采用的是POINT和LINEAR。版本二是采用POINT 和 0.2 * (LEFT + RIGHT + UP + BOT) + 0.8 * POINT,这么写的理由……想到了拉普拉斯算子,算个理由么?

1
2
3
4
5
6
7
8
9
10
11
12
13
14
关于POINT 与 LINEAR:

例如算出来实际的u = 78.4, v = 30.8

那么最近的点POINT就是mc = [78, 31],

取包含POINT的四个点mc0 = [78, 30], mc1 = [79, 31],mc2 = [79, 30],

rc1 = InterpColor(mc, mc1, 0.4);
rc2 = InterpColor(mc0, mc2, 0.4);

双线性插值的结果:final = InterpColor(rc1, rc2, 0.2);

追加:其实双线性过滤的百度百科讲得蛮好 :)

上几张图:

这是两个纹理叠加的表面,可以看到纹理投影是正确的。
SR-T2

上面的结果是下面两张纹理叠加的:

SR-T2.1
SR-T2.2

被近裁剪面裁剪之后的表面线框:

ClippedFrame

方向光照,设置的不好,结果看起来像是FLAT模式的光照……

DirLighting

Alpha混合没做,实际上知道怎么混合不同的纹理,也就知道怎么处理Alpha混合了。

1/Z-Buffer深度缓冲,基于Pixel的画家算法:

Z-Buffer

还有后续的小结

不满意,瑕疵很多,是不能投入使用的一个Renderer:比如图中三角形相接部分的黑线,实际上是没有填充入颜色,原因应该在于颜色区块的填充没有设置好,导致Float->Int的时候有的大一点,有的小一点,不连贯;比如当表面充满整个视野的时候,视口右下角会有一部分无法得到光栅化结果…两个版本多多少少都有一点这个问题;比如3D坐标齐次化并没有得到很有效的处理;再比如光照效果,从上面就能看出来没有能注意到Ambient的影响;然后就是效率有点惨……

哎,这个轮子造的甚是……

后面把重心放到Shader上面吧,终归是要跟毕业课题挂钩的……

深度缓冲、透视纹理与1/Z

其实基于Pixel的画家算法就是对每一个Pixel都保存一个深度值,每次会执行的像素的时候是对比新像素的深度值与旧像素的深度值,如果通过比较,则可以绘制新像素,否则不会绘制新像素。对于Z缓存来讲,深度值小的在前面,需要绘制,对于1/Z缓存来讲,则是相反的。

主要的问题是,光栅化插值是在2D空间中进行的线性插值,得到的深度数据并不符合3D中的实际情况,观察下图:

1/Z

稍微有点乱……

在投影平面上的线性插值是以Δx和Δy为基准的,图中可以看到,虽然沿着X进行了线性插值,但是反应到3D空间中的Z坐标上却并非是线性的结果(红色的点),实际的正确结果在投影平面的的插值如绿色点和线条所示,在X方向上并非线性的。

前面的文章中我们有提到过投影矩阵的生成过程,

1
2
X’ = (d*X0)/Z0
Y’ = (d*Y0)/Z0

并且也知道线性插值是根据X’和Y’的Δx和Δy进行的,也就是说,屏幕插值是受到1/Z的影响的,只有在包含1/Z的时候线性插值才是正确的。

事实上根本不用重新计算Z值,直接使用1/Z就可以了。

1/Z的影响对于越接近于平行投影平面的表面影响越小,反之越大,这个情况只需要调整一下上图中的红色线段的位置就知道了(如果是垂直于Z轴的话,确实是一点影响都没有的)。

理解了1/Z的作用之后,透视纹理也就不是一个大问题了:透视纹理的纹理坐标其实可以视作纹理在3D空间中的坐标位置——Z坐标一样,在屏幕坐标系中的插值计算同样受到1/Z的影响,因此需要对u/z和v/z进行线性插值,最后再将结果换算为正确的uv即可,当然在换算的过程中免不了使用到同步插值得到的1/Z。

版本一中是在投影后,由[0~1.0]的Z重新计算了实际的Z,进而得到1/Z的,版本而则是在投影之前保存了实际的Z,省去了一些步骤。

突然想到由于自己对于Shader并不是很熟练,以至于忽略的一个信息:在使用Shader来绘制深度信息到RTT的时候,大部分时间是在Pixel Shader中把插值得到的Z写入RTT的,这个时候需要注意一下写入的是实际的Z还是位于[0,1.0]的Z还是1/Z,这对重新计算实际的Z值有很大影响。

正文


本次问题主要集中在HiHo一下第16~21周,从Tarjan的Sparse-Table(ST)算法开始到线段树的使用主要讲述RMQ问题的解决方案以及离散化思想。RMQ问题在《算法竞赛入门经典——训练指南》一书中第三章也有专门讨论,并包含了这里不会讨论的树状数组。

HiHo一下第16周

与前面的一样,题目不再详述。解决这个问题使用了Tarjan的Sparse-Table方法,并且在题目的提示中已经讲解了,但是并没有给出图示,只给出了递推关系式,这里补充给出题目给出的case的解决过程来具体说明ST方法的工作原理。

首先需要明确的是ST方法实际上依旧是分治法的一种,并且需要一个表来保存计算结果以备查询(看起来与动态规划很像哈)。

来看看使用ST方法对已知case的解决过程吧:

先上主要工作,也就是构建Table的过程图:

RMQ_ST

图中的第i行第j列的数据记录以j开始,长度为2^i的区间内的最小值。很显然,我们的第0行就是原数据了。

第1行则是两个元素中的最小值,注意这个结果可以保存下来用来计算第2行的数据。

第2行则是四个元素中的最小值,同样可以用来计算第3行的数据。

第三行则是八个元素中的最小值。

由于case中的数据量为10,所以我们没有第4行了(注意我们是从0开始计数的,HiHo中的数据则是从1开始计数的)。

计算每一行的最大耗费为O(N),总共需要计算log(N)行,因此ST算法总的时间复杂度上界为O(N*logN)

既然已经计算出了查询表,那么接下来就应该使用个这表来处理查询了,处理查询的思想是这样的:

假定想知道从下标为1的元素开始,长度为7的区间内数据的最小值,那么我们可以将这个区间表示为这样的连个区间的并:

  • 下标为1的元素开始,长度为4的区间;
  • 下标为4的元素开始,长度为4的区间;

如图示:

Query1

由于区间长度为4 = 2^2,我们去找第2行中下表为14的数据:1556和981,两者的最小值981就是想要的结果。

再来一个例子:

假定想知道从下标为1的元素开始,长度为9的区间内数据的最小值,那么我们可以将这个区间表示为这样的连个区间的并:

  • 下标为1的元素开始,长度为8的区间;
  • 下标为2的元素开始,长度为8的区间;

如图示:

Query2

由于区间长度为8 = 2^3,我们去找第3行中下表为12的数据:981和981,两者的最小值981就是想要的结果。

case中的查询按照上面的思想就可以求解了,很显然经过预处理之后,查询操作的时间复杂度是O(1)

接下来给出我的TLE代码……

真是够了,跟《训练指南》中的代码比较过,除了内循环那部分,别的没出入,依旧TLE,然而改成书中的格式之后反而更慢(sad face)。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
#include <iostream>
#include <stdlib.h>
#include <string.h>
using namespace std;

#define MAX_N 1000000
#define MAX_WEIGHT 10001
#define MIN(a,b) ((a)>(b))?(b):(a)

int d[20][MAX_N];

void RMQ_ST_INIT()
{

int n = 0;
cin >> n;
for (int p = 0; p < n; ++p)
{
cin >> d[0][p];
}

for(int j = 1; (1 << j) <= n; ++j)
{
int stride = 1 << (j - 1);
for(int k = 0 ; k < n; ++k)
{
if(k + stride >= n)
{
break;
}
d[j][k] = MIN(d[j-1][k],d[j-1][ k + stride]);
}
}
}

int query(int lb, int rb)
{

lb--;
rb--;
int len = rb - lb + 1;
int k = 1;
while((1 << k) <= len){
++k;
}
k--;
return MIN(d[k][lb],d[k][lb + len - (1 << k)]);
}

int main(int argc, char** argv)
{

RMQ_ST_INIT();
int qn = 0;
cin >> qn;
while(qn > 0)
{
int l = 0, r = 0;
cin >> l >> r;
cout << query(l,r) << endl;
--qn;
}
}

HiHO一下第17周

前面讨论了RMQ问题的ST算法,这次就是它在LCA问题上的具体应用了,其中的一个关键思想就是:将树转化为数组

具体方法是这样的:我们知道DFS在工作的时候将会进入一个节点一次,退出一个节点一次(实际上就是该节点的Environment入栈和出栈啦),在这个过程中我们可以得到节点的深度信息,那么我们可以根据DFS的路径将每次遇到的深度信息放到一个数组中:

1
2
3
4
5
6
7
8
9
void dfs(Node* n, vector<int>& v)
{

foreach( node in n->children)
{
v.push_back(n->depth);
dfs(node,v);
v.push_back(n->depth);
}
}

于是对于树:

tree

上图说明了DFS的路径。

有下面的数组:

1
2
3
4
DFS路径:
[A, B, D, D, B, E, I, I, E, B, A, C, F, J, J, F, C, G, G, C, H, H, C, A]
得到的对应数组:
[0, 1, 2, 2, 1, 2, 3, 3, 2, 1, 0, 1, 2, 3, 3, 2, 1, 2, 2, 1, 2, 2, 1, 0]

回顾LCA问题,实际上就是求两个节点在树中深度最低的那个公共祖先,结社我们要求D(深度2)I(深度3)的LCA,已知结果为B,的对应深度为1

记下DI 最后一次 出现的位置下标:37,在这段区间中的最小值为1对于节点为B,正是所求的结果。
d
同理,对于节点IJ的LCA,下标714之间的区间内,最小值为0,对应节点为A,同样为所求结果。

再有,对于节点CJ的LCA,下标2214之间的区间内,最小值为1,对应节点为C,同样为所求结果。

事实上,上面的DFS路径其实是将树的无向边变为双向边得到的。

线段树

很显然,ST算法在原数组中的数据不会发生改动的时候工作的很好,而且代码也不难写,然而当原数组中的数据发生改变的时候,即便是只改变一个数据,就会导致表中数据大量的更新(看一下上面构建表过程图中的关系表示就知道了),这是不可接受的。

本想用《算法导论》上面的线段树例子,不过对比之后感觉不太合适进行应用方面的说明。还是按照《训练指南》上的说法来吧。

事实上,线段树跟前面图形方面的文中的Quad Tree很像的:子树的并集“覆盖”(之所以用引号是因为多少还是不同的,不过这种类比可以更好的理解线段树,或者反之理解Quad Tree)了父节点所代表的空间;这种类似甚至包含了对于树的查询过程:满足要求的节点直接得到结果,而完全不满足要求的则拒绝向下求解,部分满足要求的则继续向下细分。

直接借用HiHo上面的图来说明这一点:

线段树

一图胜千言,说白了就是Quad Tree的一维形式,然后标准的Quad Tree是二维的,三位的则是Oct Tree

下面我们想要查询区间[3, 9]:

查询

橘黄色的节点就是我们想要询问的节点了,这里的每个节点都保存着与该节点相关的区间的“目标数据”比如:最大值,最小值,和等等。查询最坏情况也只对两个边界来一个深度查询,因此也是O(logN)的耗费。

回想一下我们讨论的Quad/Oct Tree的工作方式,是不是几乎一样的?将空间中的物体集合换成其他信息了而已,思路都是一样的。

既然知道了查询,那么修改怎么办呢?区间修改是一个down-top的过程,我们对每个节点维护一个指向父节点的指针,自下而上修改相关区间的“目标数据”就可以了,非相关区间很明显是不受影响的,于是修改的时间耗费就成了树的高度,线段树很明显是一个二叉树,那么修改操作的时间耗费就是O(logN)

关于修改,还有一个 比较重要 的方法就是“懒惰标识”,这个思想跟“懒惰求值”是一样的,也就是不到万不得已不去修改,只做一个标记说明此处在下一次“问询”的时候需要重新计算。这种情况一般多见于对于未知问询的保守策略,不问的东西永远不会去求解,从而能提高效率(不需要每次修改把所有的点都计算一遍,只计算参与问询的点)。

线段树基本上介绍完毕了,接下来是与线段树相关的一个思想,在解决“区间相关”问题上非常有帮助(再比如2/3D空间的包含、重叠之类的问题)。

离散化

离散化思想的应用在HiHo一下第21周。下面根据题目跟出的case来说明一下离散化思想的应用。

元数据

首先是最基本的“元数据”离散化解法:

meta

将长度为N的区间分成N个单位元,这个样我们可以得到一个如下所示的线段树:

tree

绿色的叶子结点完全覆盖了整个区间,接下来我们按照顺序将给出的区间覆盖上去,这里用颜色来标注每一条线段覆盖上去之后的变化:

首先是区间[4, 10]的线段。注意,严格来讲,代表区间[5, 10]的节点以及它下属的所有节点都应该被标注为红色,但是我们这里不这样做,只是加上一个懒惰标识,在接下来的过程中会看到懒惰标识是如何发挥作用的。

4,10

接下来是区间[0, 2]

0, 2

区间[1, 6]。这个过程需要注意的是:一旦出现包含懒惰标识但依旧需要细分的节点,那么它必须先向下给自己的子节点染色,例如代表区间[5, 10]的节点。

1, 6

区间[5, 9]

5, 9

区间[3, 4]

3, 4

最后,我们深度遍历整棵树,遇到有颜色的节点就终止深度递归,记录这个过程中遇到的不同的颜色数,就可以得到结果了。结果是5,正应了case的正确结果。

然而上面这个过程需要首先对整个区间建树,其耗费与区间长度N正相关,是一个O(NlogN)的过程,如果区间很大而实际出现的线段数量很少,那么就得不偿失了。事实上,我们只需要根据线段来建树,这样可以使得问题的求解只与线段数目相关,而这只要对线段做一个预处理就可以了。

超元数据

“超元数据”离散化解法:

先上图,注意与原来的划分方法比较,并注意缺失的3条垂直虚线与结果的联系:

hyper-meta

首先根据线段的左右端点来对整个区间进行划分:

1
[0, 1] [1, 2] [2, 3] [3, 4] [4, 5] [5, 6] [6, 9] [9, 10]

按照[0, 1, 2, 3, 4, 5, 6, 9, 10]划分,对应的线段树:

tree

是不是比前面的少了一些节点呢?注意虽然少了一些节点,但是绿色的叶子节点一样覆盖了整个区间。

然后就是跟前面一样按顺序覆盖线段了:

区间[4, 10]

4 10

区间[0, 2]

0 2

区间[1, 6]

1 6

区间[5, 9]

5 9

区间[3, 4]

3 4

完毕,比前面按照“元数据”划分要来的省心多了。

最后统计的时候依旧按照前面的方法就可以了。

超元数据元数据的概念借鉴自:陈宏《数据结构的选择与算法效率——从IOI98试题PICTURE谈起》,1999.

至于代码……我的代码一直WA,讨论区的测试用例都过了(sad…),程序用的是链表,没用数组,必定TLE了(实际上通过的那些case都已经超时了)……

给出其他人的AC代码:点击这里

离散策略

在上面的题目对离散区间采用按照端点连续划分的,也就是相邻区间以[m,n] [n, k] [k, l]...这样的区间划分(我自己的代码就是这样的)

还有一种划分方法是按区间编号划分:

离散

那么对于线段[5, 9],它所对应的区间编号应为{6, 7, 8, 9},除了划分区间不同之外,其他方法一致,上面链接中的AC代码就是按照这种方式进行区间划分的。


题目没有AC感觉真是不好……自己功夫不到家,拿不到HiHo的测试数据,代码TLE了还好,WA了真是找不到错在哪里了,唉……

正文


[HiHo一下]第13~15周主要通过例子来阐述LCA问题的朴素解法以及Tarjan离线算法,同时兼顾介绍Tarjan算法使用到的一个数据结构——并查集。

HiHo一下第13周

问题描述参见标题链接,不再重复叙述。使用朴素解法即可AC,下面图示说明一下朴素解法的过程并给出相应AC代码:

LCA朴素解法

上图中的继承关系再明显不过了,求两个子节点的最近公共祖先,比如H 和 I

我们一眼就可以看出结果是A。事实上,注意从HA的路径与从IA的路径,LCA问题的朴素解法实际上是求解两条交叉链表的第一个交点。

下面是AC代码,请不要吐槽为什么这么长以及代码复用率低下…… (T- T)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
#include <iostream>
#include <string>
#include <stdlib.h>
#include <map>
using namespace std;

struct Node
{
Node* p;
string name;
};

map<string, Node*> people;

int main(int argc, char** argv)
{

int dataNum = 0;
int reqNum = 0;
cin >> dataNum;
while (dataNum > 0)
{
string father, child;
cin >> father >> child;
Node* fn = NULL;
Node* cn = NULL;
auto fit = people.find(father);
auto cit = people.find(child);
if (fit == people.end())
{
fn = new Node();
fn->name = father;
fn->p = NULL;
people.insert(pair<string, Node*>(father, fn));
}
else
{
fn = fit->second;
}

if (cit == people.end())
{
cn = new Node();
cn->name = child;
cn->p = NULL;
people.insert(pair<string, Node*>(child, cn));
}
else
{
cn = cit->second;
}
cn->p = fn;
--dataNum;
}

cin >> reqNum;
for (; reqNum > 0; --reqNum)
{
string _1, _2;
cin >> _1 >> _2;
if (_1 == _2)
{
cout << _1 << endl;
continue;
}
unsigned int step1 = 0, step2 = 0;
auto _1it = people.find(_1);
auto _2it = people.find(_2);
if (_1it == people.end() || _2it == people.end())
{
cout << "-1" << endl;
continue;
}
Node* _1n = _1it->second;
Node* _2n = _2it->second;
Node* _1temp = _1n;
Node* _2temp = _2n;
while (_1temp->p) { _1temp = _1temp->p; ++step1; }
while (_2temp->p) { _2temp = _2temp->p; ++step2; }
if (_1temp != _2temp)
{
cout << "-1" << endl;
continue;
}
_1temp = _1n;
_2temp = _2n;
bool find = false;
if (step1 > step2)
{
int s = step1 - step2;
while (s > 0)
{
_1temp = _1temp->p;
--s;
}
while (_2temp)
{
if (_1temp == _2temp)
{
cout << _1temp->name << endl;
find = true;
break;
}
_2temp = _2temp->p;
_1temp = _1temp->p;
}
}
else
{
int s = step2 - step1;
while (s > 0)
{
_2temp = _2temp->p;
--s;
}
while (_1temp)
{
if (_1temp == _2temp)
{
cout << _1temp->name << endl;
find = true;
break;
}
_2temp = _2temp->p;
_1temp = _1temp->p;
}
}
if (false == find)
{
cout << "-1" << endl;
}
}
return 0;
}

HiHo一下第14周

第14周主要是说明并查集这一数据结构的应用,题目详见标题链接。

在《算法导论》中有关于并查集的论述,参见【不相交数据结构】一章。这里只对其做基本的说明并给出题目的AC代码。

首要一点:森林与集合之间是一个等价的关系,借用《算法导论》中的一幅图来说明:

并查集与森林

上图包含了三棵树,对于每一棵树,我们可以将其视为一个集合——于是我们有{c, h, e, b},{f, d, g},{f, c, h, e, b, d, g}三个集合,并且很明显第三个集合是前两个集合的并集。每一棵树的树根可以视为是对应集合的代表元素。

合并两棵树的操作只需要将一棵树的根节点作为另一棵树的根节点的孩子插入即可。集合之间的并运算也可以这样实现。

于是在进行元素归属操作的判定问题中,我们可以对给定元素,查找到它所在的树的根节点即可。

然而树的合并带来了一个这样的问题:在极端情况下,树的深度将会被拉长为一个单链表的结构,前面归属判定的操作将会退化为对单链表的遍历操作。

事实上,考虑到纯粹的集合元素之间并没有父——子这样的关系,我们可以通过对查找(或者称为上溯)路径的压缩来加快归属判定的求解:

路径压缩

压缩过程实际上就是将判定链上的所有元素直接接到代表元素上的过程,这一过程可以通过递归实现:

1
2
3
4
5
6
7
8
9
10
11
12
ElemReference compress(ElemReference elem)
{
if (elem is set's represent element)
{
return elem;
}
else
{
elem.represent = compress(elem.represent);
return elem.represent;
}
}

了解了上面的并查集的基本内容之后,就可以对题目进行编码了,下面给出AC代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
#include <map>
#include <string>
#include <iostream>
using namespace std;

struct Node
{
Node* p;
string name;
};

map<string, Node*> mhash;

void compress(Node* n)
{

// 最后一个节点不用压缩
if (n == n->p)
{
return;
}

compress(n->p);

if (n->p && n->p->p)
{
n->p = n->p->p;
}
}

Node* getRoot(Node* n)
{

while (n != n->p)
{
n = n->p;
}
return n;
}

int main(int argc, char** argv)
{

unsigned int lineNum = 0;
cin >> lineNum;

for (;lineNum > 0;--lineNum)
{
int op = -1;
string n1, n2;
cin >> op >> n1 >> n2;
auto i1 = mhash.find(n1);
auto i2 = mhash.find(n2);
Node* p1 = NULL;
Node* p2 = NULL;
if (op == 0)
{
//known
if (i1 == mhash.end())
{
p1 = new Node();
p1->p = p1;
p1->name = n1;
mhash.insert(pair<string, Node*>(n1, p1));
}
else
{
p1 = i1->second;
}

if (i2 == mhash.end())
{
p2 = new Node();
p2->p = p2;
p2->name = n2;
mhash.insert(pair<string, Node*>(n2, p2));
}
else
{
p2 = i2->second;
}

getRoot(p1)->p = getRoot(p2);

compress(p1);
compress(p2);

}
else if (op == 1)
{
//unknown
if (i1 == mhash.end() || i2 == mhash.end())
{
cout << "no" << endl;
continue;
}
p1 = i1->second;
p2 = i2->second;


if (getRoot(p1) == getRoot(p2))
{
cout << "yes" << endl;
}
else
{
cout << "no" << endl;
}
}
}

return 0;
}

注意压缩操作的位置:可以在构建集合的过程中进行,也可以在查询的时候进行压缩,当然也可以不进行压缩——下面是三种策略对于同样的测试数据耗费的时间:

Comparative

事实上,上述代码依旧可以有优化的空间:压缩过程需要走一个查找链,同样查找链越长,递归压缩也就越耗费时间。这就要求在合并的时候尽量将比较“矮”的树合并到比较“高”的树上——即启发性合并。

HiHo一下第15周

第15周的问题主要是利用Tarjan算法来解决LCA问题,Tarjan算法用到了在14周中的并查集这一数据结构,同时该算法对递归的不同阶段进行了充分的利用,十分具有学习价值——在以后对动态语言使用的Garbage Collector的介绍中(三色标记法)将再次看到这一手段。下面对Tarjan算法进行分析,并给出相应的AC代码。

我们知道递归分为三个阶段:进入->处理->返回,Tarjan算法就是根据这三个阶段进行的:根据当前DFS的状态,给树中的节点赋予不同的颜色:未访问过的为白色,尚未递归返回的为红色,访问过的为绿色。以下图为例,我们来看看对于查询[D, E], [D, B], [I, J], [J, H]的求解过程。

Tarjan

初始情况下,所有的节点都是白色,我们DFS递归查找到节点D:

Find 'D'

此时与ID相关的查询[D, E], [D, B]中,B节点为红色,E节点为白色。通过观察,查询[D, B]的结果应该是B(__红色__)。由于查询无关元素的顺序,因此我们将查询[D, E]更改为查询[E, D],并将其延后。

DFS继续进行,到达节点E

Find 'E'

查询[E, D]的结果为B(红色),注意此时D(__绿色__)

DFS继续进行,到达节点I

Find 'I'

此时与I相关的查询[I, J]中,J为白色,将其转为查询[J, I]并延后。

DFS进行到节点J

Find 'J'

相关查询[J, I], [J, H]中,I为绿色,H为白色。与以前一样,查询[J, H],转为[H ,J]并延后。而对于查询[J, I]我们知道结果是A,而A则是由I向上遇到的__第一个红色节点__(可以给A添加一个父节点来验证一下)。

最后DFS进行到H

Find 'H'

相关查询[H, J]中,J是绿色,顺着J向上找第一个红色节点,为C,正是我们想要的结果。

至此,Tarjan算法的主要过程已经完成了,按照这个过程已经可以编写能够正确工作的代码了。但是我们应当注意到整个过程依旧有需要改进的地方:如何比较高效地找到 第一个红色节点 。至此,可以引出我们前面讨论过的并查集了。

注意前面几张图中的绿色子树与他的父跟节点,前面提到过集合可以通过树来呈现出来,我们现在将绿色的子树(单独一个节点也算)视为一个集合,将他的红色的父根节点视为该集合的代表元素(注意这里的代表元素是集合外的一个元素了),于是我们就可以通过一个O(1)的查询来得到我们想要的那个红色节点了。

前面提到过,想要得到一个O(1)的查询,需要一个递归的压缩过程。回头再看压缩操作,它实际上也是一个DFS的操作,而这恰好与我们遍历树的方式一致,因此可以在遍历的过程中完成这个操作。

下面是我的代码,没有使用数组去解这个问题,在OJ上超时了…

另附一个AC的代码链接:点击这里

再附加另外一个讲解该算法的文:点击这里

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
#include <iostream>
#include <string>
#include <vector>
#include <list>
#include <map>
using namespace std;

#define RED 2
#define GREEN 1
#define WHITE 0

struct Node
{
string name; // 名字
int color; // 颜色
vector<Node*> children; // 孩子们
Node* parent; // 父节点
Node* represent; // 所处集合的代表节点
Node(string& s)
{
name = s;
color = WHITE;
parent = NULL;
represent = this;
}
};

struct Query
{
Node* n1;
Node* n2;
string result; // 查询结果
Query()
{
n1 = NULL;
n2 = NULL;
}

Query(Node* n_1, Node* n_2)
{
n1 = n_1;
n2 = n_2;
}

void swapOrd() // 当n2是白色的时候交换n1和n2位置
{

Node* t = n1;
n1 = n2;
n2 = t;
}
};

map<string, Node*> myhash; // for find
Node* root= NULL;
vector<Query*> querys; // 原始查询列表
list<Query*> queReq; // 查询队列

Node* compress(Node* n)
{

if (n->color == RED)
{
return n;
}
n->represent = compress(n->represent);
return n-> represent;
}

void printResult()
{

for (unsigned int i = 0; i < querys.size(); ++ i)
{
cout << querys[i]->result << endl;
}
}

void dfs(Node* n)
{

if (n == 0)
{
return;
}
// 标记为红色
n->color = RED;

// 处理请求
auto it = queReq.begin();
for (; it != queReq.end();)
{
auto temp = it;
it++;

Query* q = * temp;
if (q->n1 == q->n2)
{
q->result = q->n1->name;
queReq.erase(temp);
continue;
}
if (q->n1 == n)
{
if (q->n2->color == WHITE)
{
// 对于白色的n2,交换n1与n2
q->swapOrd();
}
else if (q->n2->color == RED)
{
// 对于红色的n2,直接得到结果
q->result = q->n2->name;
// 移除当前队列元素
queReq.erase(temp);
}
else
{
// 对于绿色的n2,其结果是它的represent
compress(q->n2);
q->result = q->n2->represent->name;
// 移除当前队列元素
queReq.erase(temp);
}
}
}

// dfs
for (unsigned int i = 0; i < n->children.size(); ++i)
{
Node* child = n->children[i];

dfs(child);
}

n->color = GREEN;
};

int main(int args, char** argv)
{

int data = 0, req = 0;
cin >> data;
string s1, s2;
// 建树
Node* father = NULL;
Node* son = NULL;
while (data > 0)
{
cin >> s1 >> s2;
auto it1 = myhash.find(s1);
auto it2 = myhash.find(s2);

if (it1 == myhash.end())
{
father = new Node(s1);
myhash.insert(pair<string, Node*>(s1, father));
}
else
{
father = it1->second;
}

if (it2 == myhash.end())
{
son = new Node(s2);
myhash.insert(pair<string, Node*>(s2, son));
}
else
{
son = it2->second;
}

son->parent = father;
son->represent = father;
// a son has only one father. Not check here.
father->children.push_back(son);
data--;
}
// 确定所有人的根(根据题目说明,所有人都有一个公共祖先)
while (son->parent)
son = son->parent;
root = son;
// Query collection
cin >> req;
while (req > 0)
{
cin >> s1 >> s2;

auto it1 = myhash.find(s1);
auto it2 = myhash.find(s2);

Node* n1 = it1->second;
Node* n2 = it2->second;

Query* p = new Query(n1, n2);

querys.push_back(p);
queReq.push_back(p);

req--;
}
dfs(root);
printResult();
return 0;
}

开坑啦~


想“嫁个好家户”还是得刷题,然而不搞比赛的弱鸡刷题其实是一件蛮痛苦的事情,各种挫败遭不住……

之前有刷了一部分的LeetCode,但是感觉不够系统,对于新题不能有效地找到切入点(不包括朴素解法,然而OJ么,大家都懂的),不参考别人的解法渐渐就搞不动了。这当然也跟自己练得少,总结得少有关。因此开个坑,打算把HiHoCoder的“HiHo一下”系列搞掉,顺带做一些总结。

嗯,从试做情况来看,在HiHo上做题感觉比LeetCode要更“扎手”一些。

正文


前面有提到过之前做的Software Renderer是不包含裁剪的,因为裁剪的过程比较繁琐,无论是2D的还是3D的都有可能涉及。这次整理后的Software Renderer是要包含裁剪的过程的,因此这里就管线中的裁剪做一次说明。

3D空间下的三角形裁剪与剔除

背面剔除(背面消隐)

这个方法是最为简单有效的三角形剔除方法之一了:只需要知道三角形的法线朝向N以及相机的朝向L即可——只有当三角形是朝向观察者的时候,观察者才能看到它,否则三角形是不可见的,也就不会被绘制。

已知向量点乘结果的符号是向量投影的方向,那么点乘结果为负的两个向量,它们相互的投影必然是方向相反的。这个性质可以用来判定三角形可见与否:当Dot(N, L)的结果为负的时候,三角形朝向观察者,可以绘制;否则不予绘制。这样就完成了背面剔除。

至于三角形的法线方向是如何得到的,一般而言,我们以坐标系的标准环绕方向来计算。在本Software Renderer中是按照左手系,也就是顺时针方向为准。

视锥体裁剪/剔除

在前面介绍图形管线中的矩阵的时候有提到透视变换,透视变换实际上就是将一个3D空间中的视锥体规整为一个方体的过程。在这个过程中我们可以看到实际上我们最终在投影平面上看到场景的就是视锥体中包含的场景,那么很自然地,我们希望能够通过这个视锥体对3D空间中的物体进行剔除以及裁剪,以减小对管线的压力。

除上述目的外,视锥体裁剪的一个很重要的作用就是对穿过近裁剪平面的三角形进行裁剪,否则,这下在裁剪平面后端的点将会在投影之后产生一个小于0的深度值,从而影响深度排序的结果,造成错误的绘制(绘制一些不应该看到的东西);远裁剪平面则可以选择不去裁剪,这并不影响投影结果,对于超过远裁剪平面的点,将会在投影后得到一个大于1.0的深度值。

这里简单提一下基于画家算法的深度排序:光栅化的过程中我们需要对投影屏幕上的每一个像素进行深度值的排序,以便距离视平面近的点能够覆盖掉视平面较远的的点。这可以采用两种方法来解决问题:第一:直接使用像素的Z值,很显然,Z值越大,表示距离越远;第二,使用像素的1/Z值,那么1/Z越小,表示越远。当前我们不去深究这个问题,只需要保有这个概念即可,后面的文章会谈到。

平面与三角形

平面在本Software Renderer中是采用的法线——距离的形式定义的,法线自然指平面的单位法向量,而距离则是指平面距离远点的距离,从而我们可以通过一个点乘计算来判断空间中一个点与平面之间的关系:

1
2
3
// 3D Plane ==> { N = [Px, Py, Pz] , d }
// 3D vertex ==> { [ V = (Vx, Vy, Vz) ] }
float result = Dot(N, V) + d

通过判断上面result与0之间的关系可以得到:result = 0,则说明点在平面上,result > 0,则说明点在平面前面,否则点在平面后面。

TriangleClip

上面的的图片可以用来说明三角形裁剪的过程,这张图片的视角是贴着3D平面看过去的。让我们从点A开始,沿着顺时针方向来说明裁剪过程。

三角形的边AB与平面相交,这可以得到一个交点E;边BC不与平面相交,无交点;边CA与平面相交,得到另一个交点D;回到点A完成交点的计算过程。

几乎没什么难度是不是么?稍微麻烦一点的也就是要计算交点的纹理坐标、法向量、颜色之类的属性,但这都可以通过线性插值来完成——由于这是在3D空间中线性插值,因此不会产生任何不良影响(2D空间中插值产生的影响将会在后面透视纹理以及深度缓存部分说明)。

上面的过程推广到普通多边形就是Weiler-Atherton算法。经过上面的过程,我们可以讲一个三角形分为两个部分:ADE和CDEB。注意图中的虚线CE,如此将三角形分为三个新的三角形:DAE,CDE和EBC。

BSP(Binary Space Partition, 二元空间划分)

有了上面的平面与三角形的剪切过程,接下来我们就可以进入到另外一种绘制方式的介绍了——虽然这个方式只限于在静态场景之中使用,然而它的效益则是十分巨大的——它可以让我们抛弃Z缓存,但是同样能够得到正确的完美的绘制结果。这便是大名鼎鼎的BSP,二元空间划分。虽然在本Software Renderer中不准备去实现BSP,但是为了以后考虑,还是介绍一下吧。

至于Z缓存的具体情形将会在下一篇文章中进行介绍,届时将会看到Z缓存带来的好处以及相应的代价,这个代价说大不大说小不小,足以让DX与OpenGL的管线将其与Alpha混合、光照以及纹理采样平级,并设置一个状态开关来控制是否需要启用它。

借用《Michael Abrash’s Graphics Programming Black Book trees》, Chapter 59: The idea of BSP中的几张图,这些图是在2D情况下的BSP模拟,3D下的过程与之一致,以此为基础可以较为容易地理解3D情况。

创建

FinalBSP

图中的实线段代表实际的多边形,虚线代表多边形所在的平面,箭头指向代表平面的法向量方向。途中左边是一个场景实例,右边就是该场景对应的BSP树。

可以看到BSP树实际上就是一个二叉树。一个节点的左右子树分别代表了它前后空间的多边形集合。回忆画家算法,它实际上要求我们在绘制场景的时候先绘制较远的多边形,然后绘制近处的多边形,使用了深度缓冲的画家算法实际上是将绘制的粒度细化到一个像素——而在多边形没有相交的情形下,以多边形为粒度的画家算法同样能正常工作,BSP就是利用了这一事实。

BSP树的创建过程实际上就是利用静态场景中的多边形来对场景进行自我划分的。例如图中的A和E是被B所在的平面分割后的结果。

在创建BSP的时候,一个要点就是尽量使得每次分割的节点将自己所在的“子空间”中的多边形平分两堆,使得最后得到的的BSP树尽量保持平衡(Balanced Binary Tree),以提高效率。由于BSP树适用于静态场景,那么也就意味着这个步骤可以离线进行,从而不至影响实时渲染的速率。常用的方式是对当前划分步骤中的所有多变形进行遍历并进行评估,最后选取评估结果最好的多边形所在的平面做划分平面进行场景分割。

遍历

BSPViewTraverseOrder
假定相机位置位于图中的眼睛位置。

通过观察我们可以知道正确的多边形绘制顺序应为:A, E-> B -> C -> D,回顾上面的BSP树的形态,我们会发现这这个顺序正好与二叉树的“先右子树再左子树最后父节点”的后序遍历顺序相同。

当然这并不能完全避免像素重绘,事实上只要是以多边形为单位的绘制过程都不能完全避免像素重绘——绘制粒度太大了。BSP+扫描线渲染则可以完全避免重绘,这里不予介绍,详细情形可参见《3D游戏编程大师技巧》第13章。

由于BSP属于空间划分算法,而空间划分算法的一个重要作用就是将当前无用大部分多边形信息从流水线中剔除,因此BSP树结合视锥体剔除可以极大改善对静态场景的绘制效率。
一个简单的策略就是对于视点之后的空间直接完全剔除:注意这里的“视点之后”意指:划分平面与视向量垂直,该策略完全不考虑空间包围体,如图:

BSPReject

上面的策略不够完整,完整的剔除策略应当引入包围体:当一个子树的包围体完全被视锥体检测拒绝的时候,该子树可被完全剔除。

BSP的介绍就先到这里,不过既然开了剔除这个话题,那么下面再介绍几种常用方法来说明可见性集合的选择策略。BSP策略是针对静态场景,那么接下来将要介绍的策略不但适用于静态场景,也适用于动态场景。

四叉树(Quad-Tree)与八叉树(Oct-Tree)

四叉树与八叉树可以放在一起进行说明,他们只是维度数目上不不同而已——八叉树比四叉树多了一个维度。这里只对四叉树进行说明,八叉树可以由四叉树扩展得到。

QuadTree

上图中的黑色边框为第一层的子节点,红色划分得到的是第二层子节点,紫色划分的为第三层子节点;蓝色的两条线是视景范围,点A是视点位置。图中没有标识出远近裁剪面,不过这不影响我们的理解。

我们在创建四叉树的时候,每创建一个层级的节点,就要对应创建该节点对应的包围体(AABB或者OBB),每一个层级都将上一个层级分为四个子部分,这四个子部分的并集应当完全“覆盖”它对应的上层节点(图中是均分的,实际应用可用采用其他策略);递归进行这一过程,直到达到一个预定的阈值结束递归。

在进行剔除的时候,我们需要遍历整个树,对每个节点的包围体进行视锥体检测。完全在视景体外部的节点及其对应的子树将会被剔除,例如图中的白色部分;完全在视景体内部的节点则可以直接将其包含的渲染对象组织给渲染列表;相交的节点进一步递归进行上面的过程(图中棕色的部分)。

八叉树的剔除过程与四叉树极为类似,不再赘述。

K-D树

K-D树,也就是多维树(K-dimension tree)。一般作为数据库数据的一种组织方式,实际上它是二叉树的一种形式,只不过在树的层次上附加了根据属性剔除的意义。在某种程度上,BSP树也可以视作K-D树的一种形式(非轴对齐,而是平面对齐)。我们这里说明的是轴对齐的K-D树。

K-D树

图中所示就是轴对齐K-D树的结构——每一层都是上面一层在另一个轴上的划分,于是每次剔除都可以完全剔除当前的半个空间。

入口(Portal)

Portal

上图说明了Portal的工作过程:最初的视景体F1经过入口P1得到新的视景体F2F2经过P2和P3得到F3和F4……这是一个递归的过程,该过程中最重要的部分是新的视景体的生成。

入口需要每一帧都需要重新进行计算,最终使用生成的一系列的视景体来对场景进行剔除。

遮挡剔除

遮挡剔除主要是为了消除一部分的重复绘制——在标准的3D绘制流程中,我们是由远及近对多边形进行绘制。设想两个物体,A距离投影平面的距离大于B,也就是A距离视距较远;并且投影之后A会被B完全覆盖掉,也就是视野中只能看到B。

在上面的情况下,采取由远及近的绘制方式会导致A一定被绘制一遍,而实际上我们只需要绘制B即可,此时便可以采用遮挡剔除来避免对A的绘制了。

遮挡剔除的一个关键是要处理好物体相应的遮挡关系——一般是以一个内接体和一个外接体来实现,通过求出遮挡体在投影之后的相互遮挡关系就可以了,或者对此类情况直接采用从前到后的绘制方式绘制(这种方式只是在光栅化部分消除了像素的重绘,但是坐标变换还是要做的)。

结语


3D空间中常用的剔除方式基本算是介绍完毕了,上面那些策略的理论基础都还是比较容易的,然而实际处理中就要考虑诸多情况了。由于个人经历以及能力所限,本人不能对这些策略结合实际应用做进一步的介绍,不能不说是一个遗憾(只是实现过松散OctTree)。这篇文章先将这些方法简要整理一下,以后有机会的话会对其中某一个话题进行进一步的探讨。

写在前面


最近打算把之前写过的一个简易的仿DX固定管线的Software Renderer整理整理发到Git上面去。由于之前写这个东西的时候更偏重于理解学习,因此代码烂的一塌糊涂,实在是不能直接就丢上去——况且原来的版本只包含了能够正确显示3D物体的必要步骤(坐标变换、透视纹理以及1/Z深度缓冲,外加一个背面消隐),除此之外并没有更多的features,甚至没有裁剪(只是保证相机位置使得投影可以正常工作)。

这次重新整理编码,打算给这个玩意儿再加上一些features,比如光照,Alpha混合之类。如果情况允许的话会想尝试把DX中可编程管线部分使用脚本来模拟一下,不过这个不属于核心部分,可能不会去做,不过至少作为一个想法,这里先记一下吧。

这次重新编码的过程会被记录在一个文件中,这个文件也一并上传到了Git,通过这个记录可以大致知晓编码过程中自己是怎么思考的,另外在Blog上也会写一些东西来记录一些重点作为笔记,尽量写得易懂一些,以便帮助理解以及日后复习。

就“把复杂的东西讲得容易接受”这个事情吐个槽吧:有时候好像不加一些数学公式跟证明就显得很Low,显得没下功夫一样……难道必须加一些晦涩的东西才能显出有深度?

正文


这次先来看图形学中关于矩阵的问题。这里我们不讨论为何图形学要采用矩阵来作为主要的计算方式(这一点再全文讨论完毕之后我们自然就有答案了),只关心那些个矩阵是如何得来的。

矩阵乘法

无论是向量乘以矩阵还是矩阵乘以向量,我们都可以认为是矩阵乘以矩阵——向量本身可以视为是一个行维度或者列维度为1的特殊矩阵。

假设我们有两个矩阵 Mat01[M][N]与Mat02[L][K],注意高亮部分的NL,只有当N == L的时候这两个矩阵才可以相乘。另外,矩阵乘法是不支持交换律的。

现在时Mat01[M][N]与Mat02[N][K]相乘,那么我们将得到一个M*K的矩阵Mat03[M][K],后面来我们使用小写字母来代表矩阵元素的下标,也就是Mat[m][k]代表矩阵Mat的第m行,第k列位置的元素。有下面的规则:

Result[m][k] = MatLeft[m][0] * MatRight[0][k] + MatLeft[m][1] * MatRight[1][k] + ... + MatLeft[m][n-1] * MatRight[n-1][k]

上面的规则也说明了矩阵乘法为何没有交换律这个问题了。

在有了矩阵乘法这个概念之后,我们就可以进一步进行后续的说明了。

平移矩阵

假定我们想要移动空间中的某一个点,假如我们想要把位于[x, y, z]的点平移[l,m,n]到达[x', y' z'],也就是[x', y', z'] = [x, y, z] + [l, m, n] => {x' = x + l, y' = y + m, z' = z + n}。我们可以通过一个四维向量和一个4x4的矩阵来进行这一计算:

1
2
3
4
5
//trans
| 1 0 0 0 |
| x, y, z, 1 | * | 0 1 0 0 | = | x + l, y + m, z + n, 1 |
| 0 0 1 0 |
| l m n 1 |

上面计算过程中的矩阵就是我们需要的平移矩阵了。

旋转矩阵

旋转矩阵分为两类:固定围绕X/Y/Z轴旋转的矩阵和围绕任意一个向量旋转的矩阵。

围绕X/Y/Z轴旋转的矩阵

三者其是是相通的,这里只给围绕X轴的矩阵的推出过程,剩下的两个可以自行类比推导。

我们先来看以三角函数的方式该如何解决一个点[x, y, z]围绕X轴旋转的问题:

  • 围绕X轴进行旋转,那么也就是说,这个点的X坐标是永远不会发生改变的。
  • 在已知X坐标不会受到影响的情况下,我们可以将3D坐标投影到YOZ平面上对该问题进行求解。

(X,Y,Z)ProjOn YOZ

注意Y轴Z轴的位置关系,这里是从X轴负方向看去的。

上图中的点A: (y,z)围绕X轴转过了θ角度得到新的点A’: (y’,z’)。注意到旋转过程中点到X轴的距离不变,记其长度为L,再记OA’Y轴的夹角为α,那么OAY轴的夹角为(θ+α):那么可以有如下关系:

1
2
3
4
5
L * sinα = y’ ----------------->  ①
L * cosα = z’ -----------------> ②

L * sin(θ+α) = y --------------> ③
L * cos(θ+α) = z --------------> ④

其中③和④经由三角形和差倍角公式可以展开如下:

1
2
y = L * sinθ * cosα + L * cosθ * sinα ---------> ⑤
z = L * cosθ * cosα - L * sinθ * sinα ---------> ⑥

将①和②代入⑤和⑥:

1
2
y = y’ * cosθ + z’ * sinθ
z = z’ * cosθ - y’ * sinθ

解上面的方程,可得A’: (y’,z’)的结果(左手系):

1
2
y’ = y * cosθ - z * sinθ
z’ = y * sinθ + z * cosθ

将过程用矩阵表示出来:

1
2
3
4
5
6
7
//Rotate on X with θ
| 1 0 0 0 |
| 0 cosθ sinθ 0 |
| 0 -sinθ cosθ 0 |
| x, y, z, 1 | * | 0 0 0 1 | =

| x + l, y * cosθ - z * sinθ, y * sinθ + z * cosθ , 1 |

上面的矩阵就是左手坐标系下绕X轴旋转的矩阵了。另外两个矩阵这里不再推导,直接给出结果:

1
2
3
4
5
6
7
8
9
10
11
//Rotate on Y with θ
| cosθ 0 -sinθ 0 |
| 0 1 0 0 |
| sinθ 0 cosθ 0 |
| 0 0 0 1 |

//Rotate on Z with θ
| cosθ sinθ 0 0 |
| -sinθ cosθ 0 0 |
| 0 0 1 0 |
| 0 0 0 1 |

绕任意轴旋转的矩阵

计算围绕任意轴旋转的矩阵是一个比较繁琐的过程,然而通过上面的讨论我们知道投影是解决3D空间问题的一个有力武器,于是自然应当联想到能否通过将向量投影到与旋转轴垂直的平面上来解决问题。

来看看我们现有的已经条件:一条任意的旋转轴N: [Nx, Ny, Nz],一条待旋转的向量V: [x, y, z],需要到旋转θ角度后得到的向量V’: [x’ y’ z’]。我们将基于以下过程对问题进行求解:

  • 首先确向量N自身的坐标系:通过corss(N,V),我们可以的得到与N垂直的一个向量L,然后再经过cross(L,N)可以得到另外一个与L和N都互相垂直的向量M,于是我们可以得到一个小的坐标系统L,M,N

  • VL-M平面上投影,可以得到向量P;将V在向量N上投影,可以得到向量Q = N * Dot(V,N)。注意这两个向量是互相垂直的。于是我们只需要计算出P旋转后的结果P’,就可以由P’ + Q得到结果了。

  • 此时对于P’的求值就类似于前面的投影计算了。现在,问题集中在如何计算向量P上:由于V分为两个正交的投影向量PQ,那么V = P + Q,从而我们得到了P = V - Q。由此便可以计算出最终的V’了。

上面的过程看似需要计算(L,M,N),而实际上,可以直接通过cross(N,V)来得到这个小坐标系——V的平面投影PQcross(N,V)就是这个坐标系的三个正交轴。其中Q和我们的旋转轴N是相互重合的。

从而我们可以得到’P’ = P cosθ + cross(N,V) sinθ’,则最终结果为:
P = (V - N * Dot(V,N)) * cosθ + cross(N,V) * sinθ + N * Dot(V,N)

我们使V = [1, 0, 0],N = [Nx, Ny, Nz]来计算最终结果的X系数,这里不再写计算过程了:

1
[Nx*Nx(1-cosθ)+cosθ, Nx*Ny(1-cosθ)+Nz*sinθ, Nx*Nz(1-cosθ)-Ny*sinθ]

同理可以得到Y系数与Z系数:

1
2
3
4
// Y系数
[Nx*Ny(1-cosθ)-Nz*sinθ, Ny*Ny(1-cosθ)+cosθ, Ny*Nz(1-cosθ)+Nx*sinθ]
// Z系数
[Nx*Nz(1-cosθ)+Ny*sinθ, Ny*Nz(1-cosθ)-Nx*sinθ, Nz*Nz(1-cosθ)+cosθ]

合起来就是:

1
2
3
4
| Nx*Nx(1-cosθ)+cosθ     Nx*Ny(1-cosθ)+Nz*sinθ    Nx*Nz(1-cosθ)-Ny*sinθ    0 |
| Nx*Ny(1-cosθ)-Nz*sinθ Ny*Ny(1-cosθ)+cosθ Ny*Nz(1-cosθ)+Nx*sinθ 0 |
| Nx*Nz(1-cosθ)+Ny*sinθ Ny*Nz(1-cosθ)-Nx*sinθ Nz*Nz(1-cosθ)+cosθ 0 |
| 0 0 0 1 |

缩放矩阵

相较旋转矩阵而言,缩放矩阵就容易多了,一目了然:

1
2
3
4
5
// Scale
| Sx 0 0 0 |
|x, y, z, 1| * | 0 Sy 0 0 | = | x * Sx, y * Sy, z * Sz, 0 |
| 0 0 Sz 0 |
| 0 0 0 1 |

最基本的三种变换说完了,最后需要说明的是:使用矩阵进行变换之前必须先想好变换的顺序:先平移后旋转和先旋转后平移完全不一样。

相机矩阵

前面三个矩阵都可以用来表示物体自身的坐标系与时间坐标系之间的关系,在图形管线中World Transform之后就是View Transform,作用是将世界空间坐标系中的坐标变换到相机坐标系中。

与所有的坐标系一样,相机坐标系也有三个正交的轴,分别记为Up, Look, Right(这里略去基于欧拉角的相机系统,直接说明UNV相机系统)。相机本身又处于世界坐标系下,因此还需要一个世界中的位置。

那么如何将世界坐标系中的坐标转换到相机坐标系中呢?

我们需要把世界坐标系的远点[0, 0, 0]相机的的位置[Px, Py, Pz],然后将世界坐标系的三个正交轴与相机的三个正交轴通过旋转重合起来,当然在这个过程中所有处于世界坐标系中的点都进行了这个变换。

平移的话很简单,我们前面已经讨论过了,关键在于旋转。

回忆这样一个事实:当我们说一个点N的坐标为[x, y, z]的时候,实际上是指向量ONX, Y, Z三个轴上的投影的长度为x, y, z,并且我们默认了这样的一个前提:点N与坐标轴X, Y, Z同属于一个坐标系,且ON的起点为原点。再来看当前的情况,定义相机坐标系的UP, Look, Right三个轴实际上存在于世界坐标中,那么我们所谓的相机坐标中的位置,实际上就是所要变换的点在这三个作为坐标轴的向量上的投影长度。这样一来我们就可以不去考虑令人头痛的旋转问题了:在坐标平移之后直接进行投影即可,而投影实际上是一个点乘,我们记需要变换的坐标为V: [x, y, z],相机的位置为P: [Px, Py, Pz],三个轴的坐标分别为Up: [Ux, Uy, Uz], Right: [Rx, Ry,Rz], Look: [Lx, Ly, Lz],于是我们有:

1
2
3
4
5
6
7
8
9
y’ = Dot([x - Px, y - Py, z - Pz], [Ux, Uy, Uz])
= (x * Ux - Px * Ux) + (y * Uy - Py * Uy) + (z * Uz - Pz * Uz)
= Dot(V, Up) - Dot(P, Up)

x’ = Dot([x - Px, y - Py, z - Pz], [Rx, Ry, Rz])
= Dot(V,Right) - Dot(P, Right)

z’ = Dot([x - Px, y - Py, z - Pz], [Lx, Ly, Lz])
= Dot(V,Look) - Dot(P, Look)

上面之所以使用y->x->z的顺序是因为相机的Up实际上是要和左手坐标系下世界的Y轴重合的。
用矩阵表示就是:

1
2
3
4
|         Rx            Ux             Lx   0 |
| Ry Uy Ly 0 |
| Rz Uz Lz 0 |
| -dot(P,Up) -dot(P,Right) -dot(P,Look) 1 |

透视矩阵 ==> 裁剪矩阵

投影矩阵负责把相机坐标系中的点投规整到一个标准化方体之中,这个方体……呃,还是借用一张图吧:
PerspectiveTrans
注意到方体的“深度”值为0.0f ~ 1.0f,这与以后将会提到的深度缓冲有关。在采用右手坐标系的OpenGL中,这个区间是[-1.0f ~ 1.0f]。理论上方体大小可以任意指定,但实际中通常将方体限制在1,主要是为了方便计算。同时注意到投影中心点的位置以及坐标轴的方向,这两者在后面视口变换的过程中需要注意的点。

现在让我们开始讨论透视矩阵是如何得到并进行工作的:
PersVertex
图中是从Y轴负方向看去的XOZ平面。我们要将空间中的一点V投影到距离相机原点距离为d的投影平面上。

根据三角形相似关系,我们可以得到V在投影平面上的投影的X坐标:

1
d/Z0 = X’/X0 => X’ = (d*X0)/Z0

同理也可以得到V在投影平面上的投影的Y坐标:

1
d/Z0 = Y’/Y0 => Y’ = (d*Y0)/Z0

于是我们得到了2D平面上的一个点(X’, Y’)——在最终结果上,我们舍弃了原有的Z值,我们的矩阵应该满足:

1
| X0, Y0, Z0, 1 | * Matrix? = | (d*X0)/Z0, (d*Y0)/Z0, Z?, W? |

由于我们不关心最终结果的Z和W分量,因此可以暂时置为M = 1 和 N = 1,然后我们将X和Y分量的分母乘上去:

1
2
3
4
5
6
7
| X0, Y0, Z0, 1 | * Matrix? = | (d*X0), (d*Y0), 1*Z0, 1*Z0 |
推出矩阵:
```cpp
| d 0 0 0 | | 1 0 0 0 |
| 0 d 0 0 | ==> | 0 1 0 0 |
| 0 0 1 1 | | 0 0 1/d 1/d |
| 0 0 0 0 | | 0 0 0 0 |

貌似上面的矩阵是有问题的,得到的结构并非我们想要的结果,但实际上,4D向量需要进行一个齐次化的操作,也就是向量的所有分量都除以其W分量,就可以得到我们想要的X 和 Z了。
看起来我们只需要一个投影平面到XOY的距离d就可以了,那么我们回头看看OpenGL或者DirectX,发现创建一个透视矩阵需要的几个参数:视角Fov、投影屏幕宽高比例Aspect以及远近裁剪面FarZ与NearZ——似乎没有投影面这一说法,这又是怎么回事呢?

再次声明,我们不关心最终结果的Z和W分量,这使得我们对于能够影响到这两个分量的运算(矩阵的后面两列)可以进行适当的调整,以便于我们进行后续计算。OpenGLDirectX管线中最难理解的部分的就在这里了。

接下来引入视角。观察前面给出的草图,其中的两条斜的射线之间的夹角就是我们的视角Fov,为了分别,分别记为FovXFovY,前图中是以FovX示例。

引入视角的原因在于,进行投影的时候是不能将所有的点都进行投影的(理论上完全可以这样做),这样很有可能导致流水线被压垮,需要一定的手段进行筛选,包括远近裁剪平面也是如此,这形成了一个视锥体;再一个,也是最重要的,就是产生缩放效果的需要。这里的缩放指的是:当视角比较大的时候,我们看到的3D空间中的物体也就越多,由于屏幕的大小是固定的,因此视角越大,显示的物体就越小,从而能够显示更多的物体,视角越小则反之。我们以DirectX下的应用程序来示例一下,以帮助理解:下面是对于同样的一个显示一个柱体的程序,改变Fov的大小所带来的结果:

1
2
// 视角广,Fov = 90度
D3DXMatrixPerspectiveFovLH( &matProj, D3DX_PI/2, 1.0f, 1.0f, 100.0f );

GreatFov

1
2
// 视角窄小,Fov = 45度
D3DXMatrixPerspectiveFovLH( &matProj, D3DX_PI/4, 1.0f, 1.0f, 100.0f );

Fov

上面的对比可以清楚看到视角大小的选取对显示效果的影响。

在明白了为何需要视角之后,我们来看看视角Fov是怎样被放置到矩阵中的。很明显地,当Fov小时,投影变大,反之投影变小。由于我们唯一可以确定的就是投影平面的视距d,因此可以使用余切cot来反应这种反比例关系:cot(FOV/2),将其带入到矩阵中:

1
2
3
4
| cot(FovX/2)            0    0    0 |  
| 0 cot(FovY/2) 0 0 |
| 0 0 1/d 1/d |
| 0 0 0 0 |

然后,我们需要处理的只剩下远近裁剪平面FarZ和NearZ了。引入它们的目的前面提到一个裁剪,这里还有另外一个目的:将所有通过裁剪的点的坐标的Z值压缩到区间[0.0 ~ 1.0]中。很显然这个压缩需要除法,坐标的Z值将会除以(FarZ - NearZ)

1
2
3
4
| cot(FovX/2)            0                   0    0 |  
| 0 cot(FovY/2) 0 0 |
| 0 0 1/(d*(FarZ-NearZ)) 1/d |
| 0 0 0 0 |

等等,d的值呢?事实上d可以是任意值!回顾我们对这个矩阵的思考过程就可以发现这一事实,为了便于计算,我们令 d = 1

1
2
3
4
| cot(FovX/2)            0                   0    0 |  
| 0 cot(FovY/2) 0 0 |
| 0 0 1/(FarZ-NearZ) 1 |
| 0 0 0 0 |

但是这样得到的Z值将会是Z0/(FarZ-NearZ),很显然是不对的,应该为:(Z0-NearZ)/(FarZ-NearZ) = Z0/(FarZ-NearZ) - NearZ/(FarZ-NearZ)

1
2
3
4
| cot(FovX/2)            0                    0    0 |  
| 0 cot(FovY/2) 0 0 |
| 0 0 1/(FarZ-NearZ) 1 |
| 0 0 -NearZ/(FarZ-NearZ) 0 |

接下来又有一个但是了:经过上面的变化之后,Z’ = (Z0-NearZ)/(FarZ-NearZ),而W = Z0,其中Z’已经处于是[0.0f ~ 1.0f]之内了,齐次化之后不是更小么,显然不利于进行深度比较。于是再做一次下面的变化:

1
2
3
4
| cot(FovX/2)            0                           0    0 |  
| 0 cot(FovY/2) 0 0 |
| 0 0 FarZ/(FarZ-NearZ) 1 |
| 0 0 -(NearZ*FarZ)/(FarZ-NearZ) 0 |

现在来分别取 Z0 = NearZZ0 = FarZ来验证一下上面的矩阵:

1
2
3
4
5
6
7
8
9
// Z0 = NearZ:
Z’ = (FarZ * (NearZ-NearZ))/(FarZ-NearZ) = 0
W’ = NearZ
齐次化后 Z’ = 0

// Z0 = FarZ:
Z’ = (FarZ * (FarZ-NearZ))/(FarZ-NearZ) = FarZ
W’ = FarZ
齐次化后 Z’ = 1

至此,我们的裁剪矩阵已经完成了,还真是一段艰苦的旅程。至于其中的FovYFovX是可以相互替换的,方法就是通过d和Aspect

1
2
cot(FovX/2) = d/(Width/2) = d/((Height/2)*Aspect) = (d/(Height/2)) / Aspect = cot(FovY/2) / Aspect
cot(FovY/2) = d/(Height/2) = d/((Width/Aspect)/2) = d/(Width/2) * (Aspect) = cot(FovX/2) * Aspect

最终结果的形式取决于Fov代表水平视角还是垂直视角,这里采用了垂直视角:

1
2
3
4
5
// Persective Matrix
| cot(Fov/2)/Aspect 0 0 0 |
| 0 cot(Fov/2) 0 0 |
| 0 0 FarZ/(FarZ-NearZ) 1 |
| 0 0 -(NearZ*FarZ)/(FarZ-NearZ) 0 |

视口变换

我们终于来到了最后一个关卡了,前面的一系列变换确实有点复杂,写了几个小时脑袋都发昏了!相比之下,视口比变换就容易理解的多了,它的目的就是:得到正确的X与Y然后将其显示屏幕上。

经过前面的裁剪矩阵,我们的所处的坐标系就成了裁剪坐标系,也就是前面图中的那个方块,我们的X和Y都处于[-1.0f ~ 1.0f]之内。我们首先需要视口的宽Width与高Height,然后用它们各自的一半分别乘以X和Y,于是我们就得到了一个这样的投影结果:

PortNoOffset

对应的变换如下:

1
2
3
4
5
6
7
X’ = X*Width/2 + Width/2
Y’ = Y*Height/2 - Height/2

| width/2 0 0 0 |
| 0 height/2 0 0 |
| 0 0 1 0 |
| 0 0 0 1 |

然而考虑一下我们实际采用的窗口坐标系统,它的原点在客户区左上角的位置,并且Y轴d的方向是朝下的。这就要求我们再做另外一次也是最后一次的坐标转换:到屏幕坐标去。

首先将原来的坐标系原点移动到投影区域的左上角位置:

1
2
3
4
5
6
7
X’ = X*Width*0.5 + Width/2
Y’ = Y*Height*0.5 - Height/2

| width*0.5 0 0 0 |
| 0 height*0.5 0 0 |
| 0 0 1 0 |
| width*0.5 -height*0.5 0 1 |

然后将Y坐标翻转:

1
2
3
4
5
6
7
X’ = X*Width*0.5 + Width/2
Y’ = (-1)*Y*Height*0.5 - Height/2

| width*0.5 0 0 0 |
| 0 -height*0.5 0 0 |
| 0 0 1 0 |
| width*0.5 -height*0.5 0 1 |

上面的变换是假设Port的左上角为[0, 0]的情况下做出的变换,如果考虑视口左上角在篡改口客户区的位置,还需要加上对应的偏移量:

1
2
3
4
5
6
7
X’ = X*Width*0.5 + Width/2 + ViewPort.Left
Y’ = (-1)*Y*Height*0.5 - Height/2 + ViewPort.Top

| width*0.5 0 0 0 |
| 0 -height*0.5 0 0 |
| 0 0 1 0 |
| width*0.5+left -height*0.5+top 0 1 |

这便是最终的视口变换矩阵了。

结语

矩阵推导的部分终于结束了,重新编码的Software Renderer的数学运算部分也基本完成,这些应该是足以支撑整个程序了。上面的推导过程比较长,但这些东西是图形学的入门基础。虽然直接跟着DX或者OpenGL的教程走能写出来东西,但如果不彻底搞明白这些的话,但后面学习肯定会遇到瓶颈,以至于不得不回过头来再过一遍(我自己就属于这种情况……)。虽然这些东西在很多书上都会提及到,但是相当一部分也只是列出结果,对于思考过程语焉不详,即便有心也不一定知道如何下手,从而导致图形学上手困难。这里把我自己在推导这些矩阵的思考过程集中记录下来,以便在以后需要的时候能够帮助自己复习,同时也希望能够帮助一些刚开始学习图形学的人。

除了数学过程,Software Renderer推翻了原定使用GDI的打算,与《3D大师编程技巧》一样,使用DDraw来实现了,原因么,我只是自己想做做看而已。目前已搭好了基本的框架,可以在屏幕上点亮一些像素了:D

写在前面


在写算法的时候,我们应该会特别注意算法执行过程中每一步的“状态”如何,以及下一步“状态”将会发生什么样的改变——在DP、traceback的过程中我们尤其强调这一点。事实上,“状态”这个词可以推而广之,不仅局限在算法的范围内——一个被离散化的连续过程的每一个“切片”都有一个“状态”。这篇文章就从状态着手,来说明一些编程上的问题。

正文


状态驱动

Code Talks :)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
// 声明
struct Cond;
struct State;
struct Character;

// 条件
struct Cond
{
std::string m_condDesc;
bool operator<(const Cond& right) const
{
return strcmp(this->m_condDesc.c_str(), right.m_condDesc.c_str()) < 0;
}
bool operator==(const Cond& right) const
{
return strcmp(this->m_condDesc.c_str(), right.m_condDesc.c_str()) == 0;
}
void thisCond() { printf("%s!\n",m_condDesc.c_str()); }
Cond(const char* desc) { m_condDesc = desc; }
virtual ~Cond() {}
};
// 状态
struct State
{
std::map<Cond, State* > m_nextStates;
std::string m_stateDesc;
virtual void handleState(Character* c) = 0;
};
// 状态接受者
struct Character
{
std::string m_name;
State* m_preState;
State* m_curState;
void ReceiveCond(Cond* c)
{

auto it = m_curState->m_nextStates.find(* c);
auto end = m_curState->m_nextStates.end();
if (it != end)
{
m_preState = m_curState;
m_curState = it->second;
c->thisCond();
}
}
void HandleCurState()
{

m_curState->handleState(this);
}
void SetStart(State* s)
{

m_preState = NULL;
m_curState = s;
}
Character(const char* name) :m_name(name) {}
};
//Character的一些具体状态
struct fullOfPowerState : public State
{
void handleState(Character* c) { printf("Mr.%s is full of power!\n", c->m_name.c_str()); }
};

struct hungryState : public State
{
void handleState(Character* c) { printf("Mr.%s is hungry now, bread him.\n", c->m_name.c_str()); }
};

struct thirstyState : public State
{
void handleState(Character* c) { printf("Mr.%s is thirsty now, water him.\n", c->m_name.c_str()); }
};

struct tiredState : public State
{
void handleState(Character* c) { printf("Mr.%s is tired now, let him have a rest.\n", c->m_name.c_str()); }
};

int main(int argc, char const * argv[]) {
//定义一些状态
fullOfPowerState power;
hungryState hungry;
thirstyState thirsty;
tiredState tired;
//定义状态转移条件
Cond hungryCond("hungry");
Cond thirstyCond("thirsty");
Cond tiredCond("tired");
Cond breadCond("bread");
Cond waterCond("water");
Cond restCond("rest");
//构建状态转移图
power.m_nextStates.insert(std::pair<Cond,State* >(hungryCond,&hungry));
power.m_nextStates.insert(std::pair<Cond, State* >(thirstyCond, &thirsty));
power.m_nextStates.insert(std::pair<Cond, State* >(tiredCond, &tired));

hungry.m_nextStates.insert(std::pair<Cond, State* >(breadCond, &power));

thirsty.m_nextStates.insert(std::pair<Cond, State* >(waterCond, &power));

tired.m_nextStates.insert(std::pair<Cond, State* >(restCond, &power));
printf("\nNow Mr.K's turn\n");
//角色Mr.K
Character Mr_K("K");
Mr_K.SetStart(&power);
//Let Mr.K go
Mr_K.HandleCurState();
//饥饿
Mr_K.ReceiveCond(&hungryCond);
Mr_K.HandleCurState();
Mr_K.ReceiveCond(&breadCond);
Mr_K.HandleCurState();
//干渴
Mr_K.ReceiveCond(&thirstyCond);
Mr_K.HandleCurState();
Mr_K.ReceiveCond(&waterCond);
Mr_K.HandleCurState();
//疲累
Mr_K.ReceiveCond(&tiredCond);
Mr_K.HandleCurState();
Mr_K.ReceiveCond(&restCond);
Mr_K.HandleCurState();
printf("\nNow Mr.A's turn\n");
//角色Mr.A
Character Mr_A("A");
Mr_A.SetStart(&hungry);
//Let Mr.A go
Mr_A.HandleCurState();
//补充面包
Mr_A.ReceiveCond(&breadCond);
Mr_A.HandleCurState();

system("pause");
return 0;

从代码中可以看到,基于状态要求把“主体”的一部分行为提取出来放到与它相关的状态之中,根据当前“外界”传递给“主体”的条件来控制其对应的行为。
上面程序的状态转移图如下,还是比较简单的:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

rest cond +-----------+
+------------------------- | |
| +----------------------> | tired |
| | tired cond | |
| | +-----------+
| |
↓ |
+----------+ hungry cond +-----------+
| | ----------------> | |
|full power| | hungry |
| | <---------------- | |
+----------+ bread cond +-----------+
↑ |
| |
| | +-----------+
| | tirsty cond | |
| +-----------------------> | thirsty |
+-------------------------- | |
water cond +-----------+

接下来,我们来看一下另外一张图:
01
将上图中的if判断语句单独提取出来,放置到代码块的转移边上,再次体会一下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
             +-----------+
| i = 1 | <----+
+-----------+ |
| | (v > 20)
(u < v) | (u >= v) |
+---------+ | | |
| i = 2 | <--+ +-> +--------+
| u += 1 | ------------> | v -= 1 |
+---------+ +--------+
|
(v <= 20)
+-----------+ |
| j = 1 | <--+
+-----------+

对比之下我们可以看出即便是一个普通的过程式代码,一样可以被认为是由状态组成的。而对于函数调用,我们可以这样理解:将函数视作一个节点,而caller执行call和callee执行ret都是转移的一条边。这个模型还可以进一步套用到OS的工作过程中:cpu的每次trap都有一个条件(外部中断或者异常)来触发,进而使得操作系统在各个状态之间来回切换来完成工作。

回忆一下前面表驱动与面向对象编程一文中提到的“OO中的类以及其继承关系可以认为是一个转发表”,我们可以看出实际上“State”与“Policy”在本质上其实并无二致,实质上都是对不同过程的组织而已,只不过它们在设计方法上不同,以满足不同场景的需求而已。

其实谈到这里我们其实应当具有了一些函数式编程的基础:在函数式编程中,函数作为第一等实体,跟我们前面经常谈到的“过程”是一脉同出,最大的区别在于两点:

  • 类型、term(或者说symbol),两者在函数式编程中是分离的;而在C/C++ like的语言中,绝大部分时候两者之间是视为一体的。
  • C++中,任何时候的操作都是在一个值上面,这使得很多时候人们会将值与过程分开看待;而在函数式编程中函数本身也被视为是一个值(其实可以换一个角度去理解:对值的操作其实是对一段内存的操作;代码本身就需要被放置在代码段,自然可以被视为是一个值。当然由于代码段是只读的,那么我们可以对它的一个副本进新操作,操作这个副本在计算机看来跟操作一个值并没有什么区别)。

状态机


提到状态就不能不提状态机,我们最常听到的就是图形库以及词法/语法分析时候的状态机。实际上,前面的C++代码就是一个微型的状态机的例子。这里再举例简单的例子说明一个我个人入门时用的基于tire树的状态机,该结构可以用来解决多模板的字符串匹配问题(该问题可以在刘汝佳《算法竞赛入门经典——训练指南》中“实用数据结构——Aho-Corasick自动机”部分找到)。这里用它来说明词法分析的大致过程(嗯,我需要学习一下怎么制图了,用ASCII作画好低效劣质……)。

假设我们有6个关键字[if, ifel, int, for, foreach, float],那么我们的Token状态可以有[TK_IF, TK_IFEL, TK_FOR, TK_FOREACH, TK_ID, TK_INT, TK_FLOAT, TK_START, TK_END]。

对于6个关键字我们可以首先将其组织为一个tire树:

tire

然后把状态加到节点上,把tire中的字符作为转移条件放到边上:

tire_state

OK,这样一个简单的有限状态机就构建好了,注意对于不存在的字符边,状态被导入到TK_START以便状态机重新运作。一段这样的字符串:if float for int forint将会得到下面的Tokens:

1
2
3
4
5
lexeme: if, state: TK_IF
lexeme: float, state: TK_FLOAT
lexeme: for, state: TK_FOR
lexeme: int, state: TK_INT
lexeme: forint, state: TK_ID

当然这只是简单的示例,实际代码的情况要处理一些边界情况才能得到满意的结果。

后记


回头看了一下最近写的三篇,结果发现高中毕业不写作文之后表达能力严重退化,连个重点都没有了…还是先考虑考虑怎么把话说清楚吧。之前写过的一些代码也该整理整理发上去了——这些代码基本都将会被重写一遍,其过程将在后面的博文中体现出来。

写在前面


现在回顾起来,我对于CPP的偏爱应该是有那么点“斯德哥尔摩综合症”在里面的:我所在本科学校的编程语言教学是直接从CPP上手的,大概学校是考量到了CPP囊括了C的原因。于是我在并未体会到C的简洁优雅的情况下就被动去接触繁复浩杂的CPP(对,我也不是主动去接触CPP的)——在没有一点编程的认知下便开始了“所谓”的OO概念的灌输,结果可想而知:写CPP程序期间各种百思不解,特别是CPP里面资源管理的复杂性使我一败涂地,甚至用不了模板出手。我那可怜的在没有任何框架帮助下的本科小毕设大概也就那么3k行不到就使我有了那种“无力感”——我很清楚的记得最后程序跑起来时给我的那种无法言喻的 “WTF, It works?!” 的感觉。后面我又断断续续的写一些东西,然而无论怎么写,一旦没有了Framework的支撑,我对于代码掌控能力就降了不止一个档次。

本科期间也不是没有去接触Java和C#(那时候大数据没有火起来,Python还不够普及),然而我比较死硬,死心眼子要跟CPP刚正面,也就没转去别的语言。当时有搞C#的大神说:“C#真是门好语言,能让你确实理解OO,理解设计模式”——当然现在看来大神之所以那样说极有可能大神本身就理解了OO,C#只是更容易表达他的OO思想。至于我么,说出来让大家见笑,设计模式这个东西我当时磕了一年,终究还是没弄明白这么个东西是到底怎么提出来的——“what, why, how”,三点一个都不明白。

现在来写这篇东西(可能不止这一篇),也算是我最近一段时间的对编程的理解了吧,就从OO,设计模式以及C/CPP开始。

正文


多态

这是一个基础的不能在基础的概念了,之前去实习面试的时候不止一次被问到过,然而这又不是一个很容易就能答好的问题。我当时的在听到这个问题的时候思考了一会,可能我的肢体语言跟面部表情让面试官觉得我对多态的概念比较模糊,于是直接转去问语言层面是怎么表达了(比如该加哪些关键字之类的……)。

对于这个问题我个人的看法是:多态是程序动态特性的集中体现。正如其字面意思:接口的不同实现形式即为多态,接口不同的实现方式即为“采取的策略”。在C/CPP中的多态其实是有两种形式的:静态和动态。

静态多态

亦即是常说的编译期多态,最常见于CPP模板编程中。这个概念我第一次是见于《Modern Design C++》这本书,书内着重介绍了一种基于模板的策略编程,应用场景以及举例如下:

  1. 场景:
    假如有两个模块A和B,程序要求当A模块调用B模块中的接口时,需要先保存当前程序的状态,然后在调用返回的时候恢复为原来的状态。即使是模块A对模块B的同一个接口进行调用,也有可能根据不同的需要保存不同的状态。
  2. 场景代码
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    struct StatePolicyA
    {
    StatePolicyA()
    {
    // do save policy A
    printf("Save states with Policy A\n");
    }
    ~StatePolicyA()
    {
    // do restore
    printf("Restore states with Policy A\n");
    }
    };

    struct StatePolicyB
    {
    StatePolicyB()
    {
    // do save policy B
    printf("Save states with Policy B\n");
    }
    ~StatePolicyB()
    {
    // do restore
    printf("Restore states with Policy B\n");
    }
    };

    struct ModuleB
    {
    void interface1()
    {

    printf("Module B interface 1\n");
    }
    void interface2()
    {

    printf("Module B interface 2\n");
    }
    };

    struct ModuleA
    {
    // other functions
    ModuleB m_moduelB;

    template<typename StatePolicy>
    void callModuleB_interface01()
    {

    StatePolicy state;
    m_moduelB.interface1();
    }

    template<typename StatePolicy>
    void callModuleB_interface02()
    {

    StatePolicy state;
    m_moduelB.interface2();
    }
    };


    int main(int argc, char* argv[]) {

    ModuleA ma;
    printf("With Policy A-->\n");
    ma.callModuleB_interface01<StatePolicyA>();
    printf("\nWith Policy B-->\n");
    ma.callModuleB_interface01<StatePolicyB>();
    system("pause");
    return 0;
    }

如果有看过前文View CPP Template 的话,就会发现这个技巧其实就是利用模板向函数中传递了一个类型而,当模板被实例化的时候,类型将会被具现,从而完成任务。

我们知道CPP中引入模板的初衷就是为了提供一个比C中的宏更有力和更安全的“代码生成机制”,也就是说模板的一部分功能可以由C的宏来实现,这个技巧最有代表性的应该就是Linux中的“侵入式链表”,这里不对其做详细介绍,只举一个简单的例子来说明这种技巧,下面的例子中,如果将宏展开来,就可以看到生成了包含有int以及float两个不同类型的链表节点结构:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
#define List(type) struct List##_##type { type data; List##_##type * next;}

typedef List(int) ListInt;
typedef List(float) ListFloat;
typedef List(ListInt) ListListInt;

int main(int argc, char* argv[]) {

ListInt ihead;
ihead.data = 1;
ihead.next = NULL;
printf("Int Head: data = %i, next = %p\n",ihead.data, ihead.next);

ListFloat fhead;
fhead.data = 3.5f;
fhead.next = NULL;
printf("Float Head: data = %f, next = %p\n", fhead.data, fhead.next);

ListListInt llhead;
llhead.data = ihead;
llhead.next = NULL;
printf("ListListInt Head: data = %i, next = %p\n", llhead.data.data, llhead.next);

system("pause");
return 0;
}

事实上静态的多态比OO语言中动态的多态更容易看清楚问题的本质:实际上多态只是不同的“过程”而已,这个不同的过程被通过一个统一的接口进行描述,使得“外界”看起来一样而已。

在这里我无意探讨宏于模板之间的异同以及模板的机理,原因很简单,我对语言理论的理解尚浅,还未有够格。在这里只引用一小段话来说明C++模板的特殊之处:

1
2
3
C++ templates are essentially language-aware macros.
Each instance generates a different refinement of the same code.
—— from notes of 'Programming Languages and Translators', Fall 2012, Columbia University.

动态多态

动态多态就非常常见了,我们熟知的具有OO性质的语言都提供了这种特性,一般是通过函数覆盖来实现的,比如C++中提供了关键字virtual,而在Java中则是将所有的成员函数都视为是可以覆盖的虚函数。这里就不再赘述了,当然一同提供的还有“向上、向下”的类型转换机制。这里我们就不提供相关代码了,只大致说明一下单继承情况下的虚函数机制,多继承可以以此为根基进行拓展。

在单继承情况下,每个类都包含了一个指向虚函数表vtable的指针,这个表中记录了当前类中那些成员函数是虚函数(在C++中,如果一个根基类没有声明virtual关键字,那么是不会包含这个vtable指针的):

1
2
3
4
5
6
7
8
9
10
11
12
13
BaseClass:
+---------+
| vtable -|-------+
| | |
| ... | |
+---------+ |
|
虚函数表: |
+-----------+ <---+
| vb func 1 |
| vb func 2 |
| |
+-----------+

如果这个类有一个继承它的子类,并且它重写覆盖了基类的virtual func 2,那么它的情况将会是下面这样:

1
2
3
4
5
6
7
8
9
10
11
12
13
DeriverClass : BaseClass
+---------+
| vtable -|-------+
| | |
| ... | |
+---------+ |
|
虚函数表: |
+-----------+ <---+
| vb func 1 | 这里复制了基类的成员函数的地址
| vd func 2 | 这里更新了自己定义的成员函数的地址
| |
+-----------+

这样继承类的实例对象就能够正确找到自己新定义的函数了。当然类的析构函数与普通的函数是不同的,析构函数会经过一个链式结构从根基类的析构函数开始,按照继承链执行每一个析构函数,直到这个类自己的析构函数执行完毕结束。
当然C++标准并没有规定虚函数表的指针应该放在类的哪个位置,这要视编译器的实现而定。

表驱动

通过前面对多态的描述我们可以看到,OO的过程实际上就是对不同的“process”通过一个统一的interface进行调用的过程。process这个概念放到C里面就是C函数;放在C++里面可以是C风格的函数,也可以是用class封装而成的Functor,在C++11中更可以是一个lamdba;至于在Java和C#中就是class或者lambda,在函数式编程中就是一个Closure,或者一个Continuation。并且通过前面的论述,我们成功在动态多态中引入了表的话题。

我们先来说结论:OO语言中的类,其实是一个大的函数转发表。

这个结论其实很好理解:除却上文提到的虚成员函数之外,对于类的非虚成员函数,我们知道它是类相关的。而在C语言中,是没有类这跟概念的,这一点我们可以通过一个小小的中间层的帮助来体会一下它们之间的异同:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
void func0(void) { printf("c style func 0\n"); };
void func1(void) { printf("c style func 1\n"); };
void func2(void) { printf("c style func 2\n"); };
void func3(void) { printf("c style func 3\n"); };

typedef void(* vvfunc)(void);
vvfunc table[4] = {func0,func1,func2,func3};

void printCfunc(vvfunc* tb, int index)
{

if (index < 0 || index >= 4) {
printf("Exception : index out of range!\n");
return;
}
tb[index](); // call func
}

int main(int argc, char* argv[]) {
printCfunc(table,0);
printCfunc(table,1);
system("pause");
return 0;
}

上面是C代码,通过一个函数表,和一个索引来调用我们的函数。注意这里仅仅是为了简化说明使用了数组结构,事实上使用一个map,用函数名来索引将会更加贴近我们编写实际看到的情况(当然编译器很有可能性还是使用一个offset来索引)。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
struct Funcs
{
void func0(void) { printf("oo style func 0\n"); };
void func1(void) { printf("oo style func 1\n"); };
void func2(void) { printf("oo style func 2\n"); };
void func3(void) { printf("oo style func 3\n"); };
};

Funcs oo_tb;

int main(int argc, char* argv[]) {
oo_tb.func0();
oo_tb.func1();
system("pause");
return 0;
}

上面是C++代码,通过类直接调用了函数。对比一下前面C风格的代码,这里的Funcs oo_tb就相当于vvfunc table[4]main函数中直接对成员函数的调用就相当于C代码部分由printCfunc通过一个函数表指针(C++中就是那个众所周知的this了)和一个索引(在C++中是无法看到的,这一步由编译器帮忙做了)来调用需要的函数。

上面的例子还是比较简单的,如果加入函数重载的特性,那么就需要另外考虑到语义的问题了(函数类型),不适合用来做一个简单明了的说明。但无论加不加重载,最基本的道理是相同的,都是通过转发表来实现。如果想更加了解表驱动编程,推荐阅读《代码大全》相关章节。

结语


这边文章的中心意思基本讲述完毕了,最后让我们以两个设计模式的例子来加深一下对文章的理解:

策略模式

  • C++代码

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    struct Plan
    {
    virtual void proc(void) = 0;
    };

    struct Plan_A : public Plan
    {
    virtual void proc(void) { printf("Plan A is choosed\n"); }
    };

    struct Plan_B : public Plan
    {
    virtual void proc(void) { printf("Plan B is choosed\n"); }
    };

    struct TheMan
    {
    void choose(Plan& plan) { plan.proc(); }
    };

    int main(int argc, char * argv[]) {
    TheMan man;
    Plan_A a;
    Plan_B b;
    man.choose(a);
    man.choose(b);
    system("pause");
    return 0;
    }
  • 对应的C代码

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    void PlanA(void) { printf("Plan A is choosed\n"); }
    void PlanB(void) { printf("Plan B is choosed\n"); }

    void man_choose(void(* pfun)(void))
    {
    pfun();
    }

    int main(int argc, char * argv[]) {
    man_choose(PlanA);
    man_choose(PlanB);
    system("pause");
    return 0;
    }

这个还是比较好理解的,只用了一次转发。

访问者模式

  • C++代码

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    // 写完才发现今天是12.25...
    struct Accepter_A;
    struct Accepter_B;

    //////////////////////////////////////////////////////////////////////////
    // visitor
    struct Visitor
    {
    virtual void visit(Accepter_A& accepter) = 0;
    virtual void visit(Accepter_B& accepter) = 0;
    };

    struct Visitor_SelfIntro : public Visitor
    {
    void visit(Accepter_A& accepter) { printf("I am A\n"); }
    void visit(Accepter_B& accepter) { printf("I am B\n"); }
    };

    struct Visitor_Greetings : public Visitor
    {
    void visit(Accepter_A& accepter) { printf("Hey~\n"); }
    void visit(Accepter_B& accepter) { printf("Hi~\n"); }
    };

    //////////////////////////////////////////////////////////////////////////
    // accepter
    struct Accepter
    {
    virtual void accept(Visitor& v) = 0;
    };

    struct Accepter_A
    {
    virtual void accept(Visitor& v) { v.visit( * this); }
    };

    struct Accepter_B
    {
    virtual void accept(Visitor& v) { v.visit( * this); }
    };

    int main(int argc, char * argv[]) {
    Accepter_A a;
    Accepter_B b;
    Visitor_SelfIntro intro;
    Visitor_Greetings greeting;

    a.accept(greeting);
    a.accept(intro);

    b.accept(greeting);
    b.accept(intro);

    system("pause");
    return 0;
    }
  • visitor模式被认为是GoF书中较为难以理解的模式,主要在于它使用了双重转发。我们先来梳理它是怎么工作的,然后将其换位对应的C代码:

  1. 过程:
    每一个accpter接受一个visitor,并根据visitor的具体实例来做出对应的动作。
    每个visitor定义所有的accpeter的具体动作,实际工作的时候只调用accpter来接受不同的visitor就可以达到目的了。所谓的双重转发是指:第一次转发在接受visitor时的动态转发,第二次转发在visitor根据接受它的accpter的类型进行一次静态转发(函数重载)。
  2. 对应的C代码
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    typedef void(* vvfunc)(void);
    void visitor_selfintro_A() { printf ("I am A\n");}
    void visitor_selfintro_B() { printf ("I am B\n");}
    vvfunc visitor_selfintro[2] = {visitor_selfintro_A,visitor_selfintro_B};

    void visitor_greetings_A() { printf ("Hey~\n");}
    void visitor_greetings_B() { printf ("Hi~\n");}
    vvfunc visitor_greetings[2] = {visitor_greetings_A,visitor_greetings_B};

    void accepter_A(vvfunc * vtb) {vtb[0]();}//注意,这里的转发过程直接写进去了
    void accepter_B(vvfunc * vtb) {vtb[1]();}//注意,这里的转发过程直接写进去了

    int main(int argc, char * argv[]) {
    accepter_A(visitor_greetings);
    accepter_A(visitor_selfintro);
    accepter_B(visitor_greetings);
    accepter_B(visitor_selfintro);
    system("pause");
    return 0;
    }

尝试把OO设计模式用C表述一下,结果发现蛮好玩的,后面应该还会继续做,争取把常用的几个模式给过一遍。

收工,写篇文章真心花时间……

写在前面


写文一向是一件我不太拿手的事情,大多数时候单是起头就能让我大半天没法下笔:比如这次,大概十多分钟才写了这么两句。

决定写文的原因很简单:看到简历之后别人几乎都会问那么一句:“有没有上传到Github上去?”。经历了几次阵痛之后我终于觉得这应该是个比较重要的问题了:毕竟跟有Git加持的人相比,在话题主导方面我显然不占优势,况且我也不是个会聊天的人。于是乎,在注册了Git一年之后我终于下决心要在上面PO文章交代码了。

至于为什么用这样一篇分析报告这个作为第一篇博文,应该出于:选一个稍微高上一点的,自己多少能看懂的。另外,毕竟CPP用惯了,更亲近一些。

正文


本文的主线是一篇Matt教授的文章。Matt教授的这篇文章主要讲述了C++模板的图灵完备性(Turing-completeness),并以用模板构建了一个“高阶函数式编程语言”来证明之。之所以在Matt教授主页上面的C++文章中选择这样的一篇文章,一方面是因为我个人为这是搞懂C++模板,进而踏入C++泛型编程领域的关键所在;另一方面是也因为这篇文章涉及到的代码量只有200行左右,不至于使得读者过于专注代码细节而影响对整体思想的理解,并且这也是一篇代码分析报告所能承受的量。

当然本人水平有限,文中不免有错误之处,希望有心的读者来邮件(solaxu@outlook.com)批评指正。

模板特化

模板和模板特化赋予了C++模板图灵完备性,这使得C++模板能够像无类型的term-rewriting语言(例如ML、Lisp)一样工作。下面是一个模板和它特化的例子:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
// ListNode模板类
template <typename A>
struct ListNode
{
A data;
ListNode<A>* next;
};

// ListNode模板类的一个特殊化形式,注意该形式实现了对ListNode<uinsigned int>类型的屏蔽
template <>
struct ListNode<unsigned int>
{
int data;
ListNode<int>* next;
};

对于模板实例化这里就不再赘述,我们只看上面模板特化部分的代码:当我们的客户端程序在使用ListNode希望实现一个数据类型为uinsigned int类型的ListNode的时候,它将会被实际实例化为一个数据类型为int的ListNode,这样就实现了对ListNode类型的屏蔽。

一个求阶乘的例子

这个例子几乎出现在任意一本讲述C++模板编程的的书籍中:

1
2
3
4
5
6
7
8
9
10
11
template <int N>
struct Factorial
{
enum {value = N * Factorial<N-1>::value};
};

template <>
struct Factorial<0>
{
enum {value = 1};
};

假定我们需要求Factorial<3>,已知结果为123=6,那么上面代码的工作方式可以表示如下(编译期工作):

1
2
3
4
5
6
7
Factorial<3>
{
value = 3 * (Factorial<2>)
= 3 * ( 2 * (Factorial<1>) )
= 3 * ( 2 * ( 1 * (Factorial<0>) ) )
= 3 * ( 2 * ( 1 * (1) ) )
}

在上面的代码中,每次对 Factorial<N> 的递归都先将N“压入” Factorial<N-1>的 “工作环境”之中,这个过程并没有计算(如果对Closure比较熟悉的话这应该不是一个难以理解的点)。也就是说,Factorial<N>每次计算都生成了一个“闭包”,该闭包接收一个参数N-1,直到N=0时,遇到特化模板Factorial<0>结束。最终Factorial<3>的展开式为3 * 2 * 1 * 1,这就可由编译器的常数优化来进行直接求结果了。

数据结构与类型匹配

这个小节的标题原文是“Encoding and pattern-matching data structures”,但是译作“编码与模式匹配数据结构”这个怎么感觉都不对味,所以做主张调整了标题。本小节建议诸位读者大人打开原文,与本文两厢对比来看。

在C++元编程中,“类型”(types)被作为数据结构的代指,而“特化”则用来“在某种程度上破坏”这种代指(参见第一小节对unsigned int的屏蔽)。请看下面代码:

1
2
3
4
5
6
7
8
9
10
struct Zero
{
enum {value = 0};
};

template <typename N>
struct Succ
{
enum {value = N::value + 1};
}

有了第二小节的铺垫,上面代码不难理解。在上述代码中,使用类型Zero来代指自然数0,类型Succ<N>来表示自然数N后面的一个数。这里需要注意的一点是:类模板其实是一个“类型运算函数”,它接收一个或者多个类型(也就是模板参数),最后得到一个新的类型。于是自然数1可以表示为Succ<Zero>自然数2可以表示为Succ< Succ<Zero> >…依次类推就可以表示所有的自然数了。

1
2
3
4
5
6
7
8
9
10
11
template <typename N>
struct MatchOne
{
enum{value = 0};
}

template <>
struct MatchOne<Succ<Zero> >
{
enum{value = 1};
}

上面的代码可以这样工作:

1
2
3
cout << MatchOne<Zero>::value << endl; //prints 0
cout << MatchOne<Succ<Zero> >::value <<endl; // prints 1
cout << MatchOne<Succ<Succ<Zero> > >::value << endl; // prints 0

打印结果使用模板特化的知识很容易理解:只有特定类型才能产生特定的输出,这里是1;其他类型都产出0

C++模板的图灵完备性

前面铺垫了许久,终于到了重头戏。
众所周知,λ演算是图灵完备的,那么C++的模板只需要能够支持完整的λ演算,就可以证明它的图灵完备性了,该小节将会通过200行左右的代码来展示模板的这种特性,Matt教授的代码可以在这里下载到。

接下来的代码中Environment可以认为是一个“运行环境”,这里的“运行环境不同于 经常提到的“程序的运行环境”。这里特指“一小段代码”(例如一个Closure),可以回顾一下前面提到的Factorial的“运行环境”。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

struct EmptyEnv;

template <int Name, typename Value, typename Env>
struct Binding{};

// EnvLookup 的工作方式与前面的 Factorial类似
template <int Name, typename Env>
struct EnvLookup{};

template <int Name>
struct EnvLookup<Name,EmptyEnv>{};// Name Not found, 特化以终止递归

// 针对Binding类型特化
template <int Name, typename Value, typename Env>
struct EnvLookup <Name, Binding<Name,Value,Env> >
{
Value typedef result;
};

template <int Name, int Name2, typename Value2, typename Env>
struct EnvLookup <Name, Binding<Name2,Value2,Env> >
{
typename EnvLookup<Name,Env> :: result typedef result;
} ;

我们先来看看main函数里面对EnvLookup的使用:

1
2
3
4
5
6
7
8
// Testing [2 => 1, 3 => 0](3):
int v = EnvLookup<3, Binding<2,
Succ<Zero>,
Binding<3, Zero,EmptyEnv>
>
> :: result :: value ;

assert(v == 0) ;

首先从Binding入手。
Binding的作用是:在一个Environment中,把一个Value绑定到一个Name上。那么前面的代码首先把一个Zero绑定到名字为3的引用上,而这个引用位于环境EmptyEnv中,于是在进行了Binding操作之后,我们得到了一个新的Environment,假设我们对其命名为EmptyEnvWith3=0,接下来我们再把Succ<Zero>这个值再绑定到EmptyEnvWith3=0中的名字为2的引用上。这个过程可以表述如图:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
EmptyEnv:
+--------------------+
| ref 1 = Null |
| ref 2 = Null |
| ref 3 = Null |
| ... |
| ref 4 = 3 + 2 |
| ... |
+--------------------+
这里的EmptyEnv只是名字如此,并不意味着它就一定是一个“空”的,另外这些值为 Null 的引用代表着它们没有相对应的值,只“占位”而已。
而且这里的Environment虽然看起来似乎是一个结构体,但是实际上并不是,可以在某种程度上认为它是一段代码。
(可以查找一下Continuation的相关信息帮助理解,它与Closure是近亲)。

经过第一次 Binding 操作之后:
EmptyEnv:(在上面我们给他另外起名为 "EmptyEnvWith3=0" 实际上还是 EmptyEnv,临时变量,右值而已)
+--------------------+
| ref 1 = Null |
| ref 2 = Null |
| ref 3 = Zero |
| ... |
| ref 4 = 3 + 2 |
| ... |
+--------------------+

经过第二次 Binding 操作之后:
EmptyEnv:
+--------------------+
| ref 1 = Null |
| ref 2 = Succ<Zero> |
| ref 3 = Zero |
| ... |
| ref 4 = 3 + 2 |
| ... |
+--------------------+

于是最后轮到EnvLookup出手了:它使用与前面求Factorial一样的过程来求出两次Binding操作之后EmptyEnv对应名字的result,然后通过这个result来获取到具体的值。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
// Core syntax
template <int FormalName, typename Body>
struct Lambda{};

template <typename Fun, typename Arg>
struct App{};

template <int Name>
struct Ref{};

// Sugar:
template <typename Cond, typename Then, typename Else>
struct If{};

template <typename T>
struct Lit{};

// Values
template <typename Lam, typename Env>
struct Closure {} ;

struct True {} ;
struct False {} ;

// Eval procs
template <typename Exp, typename Env>
struct Eval {} ;

template <typename Proc, typename Value>
struct Apply {} ;

template <typename T, typename Env>
struct Eval <Lit<T>, Env>
{
T typedef result ;
} ;

template <int Name, typename Env>
struct Eval <Ref<Name>, Env>
{
typename EnvLookup<Name, Env> :: result typedef result ;
} ;

template <int Name, typename Body, typename Env>
struct Eval <Lambda<Name,Body>, Env>
{
Closure<Lambda<Name, Body>, Env> typedef result ;
} ;

template <typename Fun, typename Arg, typename Env>
struct Eval<App<Fun, Arg> , Env> {
typename Apply<typename Eval<Fun,Env> :: result ,
typename Eval<Arg,Env> :: result > :: result
typedef result ;
} ;

template <typename Then, typename Else, typename Env>
struct Eval<If<True,Then,Else>,Env> {
typename Eval<Then,Env> :: result typedef result ;
} ;

template <typename Cond, typename Then, typename Else, typename Env>
struct Eval<If<Cond,Then,Else>,Env> {
typename Eval<If<typename Eval<Cond,Env> :: result,
Then,
Else>,
Env> :: result
typedef result ;
} ;

template <int Name, typename Body, typename Env, typename Value>
struct Apply<Closure<Lambda<Name,Body>, Env>, Value> {
typename Eval<Body, Binding<Name,Value,Env> > :: result
typedef result ;
} ;

好长一段代码……
还是先看看main中是怎么使用的吧:

1
2
3
4
5
6
7
8
9
10
11
12
// Testing ((lambda (x) x) 2):
enum { X } ;
int x = Eval <
App <
Lambda<X,Ref<X> >,
Lit <
Succ < Succ<Zero> > // ===> 2
>
>,
EmptyEnv
> :: result :: value ;
assert(x == 2) ;

首先把代码中的Succ <Succ <Zero> >简化为2
然后是Lambda<Name,Body>,相当于一个这样的表达:

1
2
lambda x :
return x

于是这里的App<Lambda,Lit>实际上就是上面第49行代码中的App<Fun, Arg>,其作用是:将参数Arg应用到一个方法Fun中,最后得到依旧是一个方法(这一点需要特别注意,这个方法并没有直接得到结果,而是形成了一个闭包被返回了)。
那么Eval<App,EmptyEnv>就成了:将上面得到的闭包放置在EmptyEnv这个“环境”之中进行求值从而得到结果。至于为什么会得到2,请自行按照模板展开的方式一层一层剥离下去就可以了,注意模板特化带来的影响,这里不再赘述。

再来看main中的另外一段代码:

1
2
3
4
5
6
7
8
9
10
// Testing (if #f 0 1):
int y = Eval <
If <
Lit<False>,
Lit<Zero>,
Lit<Succ<Zero> >
>,
EmptyEnv
> :: result :: value ;
assert(y == 1) ;

通过前面的分析,我们可以很容易的就知道这段代码的目的了:把If放置在EmptyEnv中执行一次,结果是Succ<Zero>,也就是1

结语


Matt教授的文章跟代码分析完了,通过分析我们可以看到C++模板是图灵完备的,这使得C++模板元编程成为了可能。
在文章的结束部分,Matt教授给出了一个Exercise,这个Exercise需要我们自己在给出的框架上填写一段代码来实现自然数的加法,其实就是把两个代表不同数的Succ合并起来成为一个Succ,然后传递给Eval

1
2
3
4
5
6
7
8
9
10
// 模板特化Add,感觉有点不太对劲……虽然Pass了
template <>
struct Add
<
Succ<Succ<Zero> >,
Succ<Succ<Succ<Zero> > >
>
{
Succ<Succ<Succ<Succ<Succ<Zero> > > > > typedef result;
};