写在前面
最近打算把之前写过的一个简易的仿DX固定管线的Software Renderer整理整理发到Git上面去。由于之前写这个东西的时候更偏重于理解学习，因此代码烂的一塌糊涂，实在是不能直接就丢上去——况且原来的版本只包含了能够正确显示3D物体的必要步骤（坐标变换、透视纹理以及1/Z深度缓冲，外加一个背面消隐），除此之外并没有更多的features，甚至没有裁剪（只是保证相机位置使得投影可以正常工作）。
这次重新整理编码，打算给这个玩意儿再加上一些features，比如光照，Alpha混合之类。如果情况允许的话会想尝试把DX中可编程管线部分使用脚本来模拟一下，不过这个不属于核心部分，可能不会去做，不过至少作为一个想法，这里先记一下吧。
这次重新编码的过程会被记录在一个文件中，这个文件也一并上传到了Git，通过这个记录可以大致知晓编码过程中自己是怎么思考的，另外在Blog上也会写一些东西来记录一些重点作为笔记，尽量写得易懂一些，以便帮助理解以及日后复习。
就“把复杂的东西讲得容易接受”这个事情吐个槽吧：有时候好像不加一些数学公式跟证明就显得很Low，显得没下功夫一样……难道必须加一些晦涩的东西才能显出有深度？
正文
这次先来看图形学中关于矩阵的问题。这里我们不讨论为何图形学要采用矩阵来作为主要的计算方式（这一点再全文讨论完毕之后我们自然就有答案了），只关心那些个矩阵是如何得来的。
矩阵乘法
无论是向量乘以矩阵还是矩阵乘以向量，我们都可以认为是矩阵乘以矩阵——向量本身可以视为是一个行维度或者列维度为1的特殊矩阵。
假设我们有两个矩阵 Mat01[M][N]与Mat02[L][K]，注意高亮部分的N和L，只有当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平面上对该问题进行求解。
注意Y轴和Z轴的位置关系，这里是从X轴负方向看去的。
上图中的点A: (y,z)围绕X轴转过了θ角度得到新的点A’: (y’,z’)。注意到旋转过程中点到X轴的距离不变，记其长度为L，再记OA’与Y轴的夹角为α，那么OA与Y轴的夹角为(θ+α)：那么可以有如下关系：1
2
3
4
5L * sinα = y’ -----------------> ①
L * cosα = z’ -----------------> ②
L * sin(θ+α) = y --------------> ③
L * cos(θ+α) = z --------------> ④
其中③和④经由三角形和差倍角公式可以展开如下：1
2y = L * sinθ * cosα + L * cosθ * sinα ---------> ⑤
z = L * cosθ * cosα - L * sinθ * sinα ---------> ⑥
将①和②代入⑤和⑥：1
2y = y’ * cosθ + z’ * sinθ
z = z’ * cosθ - y’ * sinθ
解上面的方程，可得A’: (y’,z’)的结果（左手系）：1
2y’ = 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 | //Rotate on Y with θ |
绕任意轴旋转的矩阵
计算围绕任意轴旋转的矩阵是一个比较繁琐的过程，然而通过上面的讨论我们知道投影是解决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。将
V在L-M平面上投影，可以得到向量P；将V在向量N上投影，可以得到向量Q = N * Dot(V,N)。注意这两个向量是互相垂直的。于是我们只需要计算出P旋转后的结果P’，就可以由P’ + Q得到结果了。此时对于
P’的求值就类似于前面的投影计算了。现在，问题集中在如何计算向量P上：由于V分为两个正交的投影向量P和Q，那么V = P + Q，从而我们得到了P = V - Q。由此便可以计算出最终的V’了。
上面的过程看似需要计算(L,M,N)，而实际上，可以直接通过cross(N,V)来得到这个小坐标系——V的平面投影P，Q，cross(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]的时候，实际上是指向量ON在X, 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
9y’ = 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 |
透视矩阵 ==> 裁剪矩阵
投影矩阵负责把相机坐标系中的点投规整到一个标准化方体之中，这个方体……呃，还是借用一张图吧：
注意到方体的“深度”值为0.0f ~ 1.0f，这与以后将会提到的深度缓冲有关。在采用右手坐标系的OpenGL中，这个区间是[-1.0f ~ 1.0f]。理论上方体大小可以任意指定，但实际中通常将方体限制在1，主要是为了方便计算。同时注意到投影中心点的位置以及坐标轴的方向，这两者在后面视口变换的过程中需要注意的点。
现在让我们开始讨论透视矩阵是如何得到并进行工作的：
图中是从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分量，这使得我们对于能够影响到这两个分量的运算（矩阵的后面两列）可以进行适当的调整，以便于我们进行后续计算。OpenGL与DirectX管线中最难理解的部分的就在这里了。
接下来引入视角。观察前面给出的草图，其中的两条斜的射线之间的夹角就是我们的视角Fov，为了分别，分别记为FovX，FovY，前图中是以FovX示例。
引入视角的原因在于，进行投影的时候是不能将所有的点都进行投影的（理论上完全可以这样做），这样很有可能导致流水线被压垮，需要一定的手段进行筛选，包括远近裁剪平面也是如此，这形成了一个视锥体；再一个，也是最重要的，就是产生缩放效果的需要。这里的缩放指的是：当视角比较大的时候，我们看到的3D空间中的物体也就越多，由于屏幕的大小是固定的，因此视角越大，显示的物体就越小，从而能够显示更多的物体，视角越小则反之。我们以DirectX下的应用程序来示例一下，以帮助理解：下面是对于同样的一个显示一个柱体的程序，改变Fov的大小所带来的结果：1
2// 视角广，Fov = 90度
D3DXMatrixPerspectiveFovLH( &matProj, D3DX_PI/2, 1.0f, 1.0f, 100.0f );
1
2// 视角窄小，Fov = 45度
D3DXMatrixPerspectiveFovLH( &matProj, D3DX_PI/4, 1.0f, 1.0f, 100.0f );
上面的对比可以清楚看到视角大小的选取对显示效果的影响。
在明白了为何需要视角之后，我们来看看视角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 = NearZ 和 Z0 = 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
至此，我们的裁剪矩阵已经完成了，还真是一段艰苦的旅程。至于其中的FovY 与 FovX是可以相互替换的，方法就是通过d和Aspect：1
2cot(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,于是我们就得到了一个这样的投影结果：
对应的变换如下：1
2
3
4
5
6
7X’ = 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
7X’ = 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
7X’ = 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 | X’ = X*Width*0.5 + Width/2 + ViewPort.Left |
这便是最终的视口变换矩阵了。
结语
矩阵推导的部分终于结束了，重新编码的Software Renderer的数学运算部分也基本完成，这些应该是足以支撑整个程序了。上面的推导过程比较长，但这些东西是图形学的入门基础。虽然直接跟着DX或者OpenGL的教程走能写出来东西，但如果不彻底搞明白这些的话，但后面学习肯定会遇到瓶颈，以至于不得不回过头来再过一遍（我自己就属于这种情况……）。虽然这些东西在很多书上都会提及到，但是相当一部分也只是列出结果，对于思考过程语焉不详，即便有心也不一定知道如何下手，从而导致图形学上手困难。这里把我自己在推导这些矩阵的思考过程集中记录下来，以便在以后需要的时候能够帮助自己复习，同时也希望能够帮助一些刚开始学习图形学的人。
除了数学过程，Software Renderer推翻了原定使用GDI的打算，与《3D大师编程技巧》一样，使用DDraw来实现了，原因么，我只是自己想做做看而已。目前已搭好了基本的框架，可以在屏幕上点亮一些像素了:D