在此回顾记录2019年12月实习期间编写Soft Renderer时的思考。
主要模块
以下为我实现Soft Renderer时划分的主要模块
- 坐标变换MVP + 视口变换
- Bresenham划线算法
- Scanline-Zbuffer多边形填充
- Phong光照模型 + Phong着色
- 纹理映射（透视校正插值）
管线流程
模型读取
顶点的数据结构设定如下，其中pos，normal，texcoord从.obj文件中读取：
class Vertex |
三角面的数据结构如下：
class Triangle |
坐标变换MVP + 视口变换
本项目中采用的是Direct3D标准，即所有向量为行向量，乘法顺序为从左向右。
Model Matrix
模型读取完成后，首先对模型本身的平移(Translate)、旋转(Rotate)、缩放(Scale)三种变换进行处理。
TRS中，必须保证T要在最后相乘。
$$
Model Matrix = Scale \times Rotation \times Translate
$$
void Transform::UpdateModelMatrix() |
Scale
Scale矩阵的设置为：
$$
\left[
\begin{matrix}
scaleX & 0 & 0 & 0 \\
0 & scaleY & 0 & 0 \\
0 & 0 & scaleZ & 0 \\
0 & 0 & 0 & 1
\end{matrix}
\right]
$$
void Transform::SetScale(float x, float y, float z) |
Rotation
欧拉旋转的顺序参照Unity标准，即ZXY。
$$
Rotation = Rz \times Rx \times Ry
$$
其中绕不同轴的Rotation矩阵的设置为：
$$
Rx =
\left[
\begin{matrix}
1 & 0 & 0 & 0\\
0 & \cos x & \sin x & 0\\
0 & -\sin x & \cos x & 0 \\
0 & 0 & 0 & 1
\end{matrix}
\right]
Ry =
\left[
\begin{matrix}
\cos y & 0 & -\sin y & 0 \\
0 & 1 & 0 & 0 \\
\sin y & 0 & \cos y & 0 \\
0 & 0 & 0 & 1
\end{matrix}
\right]
Rz =
\left[
\begin{matrix}
\cos z & \sin z & 0 & 0 \\
-\sin z & \cos z & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & 0 & 1
\end{matrix}
\right]
$$
void Transform::SetRotation(float x, float y, float z) |
其中D2R函数含义为Degree to Radian
const float PI = 3.14159265358f; |
Rotation矩阵必然正交，原理为：
$$
\left[
\begin{matrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & 0 & 1
\end{matrix}
\right]
\text {（原坐标轴）}
\times R =
R\text {（新坐标轴）}
$$
经旋转后的新坐标轴必然保持正交。
Translate
Translate矩阵的设置为：
$$
\left[
\begin{matrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 1 & 0 \\
TranslateX & TranslateY & TranslateZ & 1
\end{matrix}
\right]
$$
void Transform::SetTranslate(float x, float y, float z) |
View Matrix
接下来，需要将顶点由World Space转入View Space，思路为：
$$
任一点 \xrightarrow{Follow} 相机 \xrightarrow{Translate} 原点 \xrightarrow{Rotate} 将uvw旋转至与xyz重合
$$
相机的数据结构为：
class Camera { |
求出View Space的坐标向量uvw：
Vector3 forward = (camera.lookAt - camera.pos).Normalize(); |
$$
u = right \\v = up \\w = forward
$$
接下来的目标是将View Space的坐标轴uvw乘以View Matrix变换至与World Space的坐标轴xyz重合，即：
$$
u \rightarrow x \\v \rightarrow y \\w \rightarrow z
$$
而由上述思路可知：
$$
Matrix_{view} = M_T \times M_R
$$
首先易知：
$$
M_T =\left[\begin{matrix}1 & 0 & 0 & 0 \\0 & 1 & 0 & 0 \\0 & 0 & 1 & 0 \\-Camera_{posX} & -Camera_{posY} & -Camera_{posZ} & 1\end{matrix}\right]
$$
乘以$M_T$后，此时两坐标系原点已重合，然后考虑 $M_R$，首先考虑其逆过程$xyz \rightarrow uvw $：
（uvw为方向向量，不会因上一步的平移而改变）
$$
[1, 0, 0] \times R = [u_x, u_y, u_z] \\ [0, 1, 0] \times R = [v_x, v_y, v_z] \\ [0, 0, 1] \times R = [w_x, w_y, w_z]
$$
则对于$uvw \rightarrow xyz $
$$
M_R = R^{-1}
$$
由旋转矩阵必正交可得：
$$
R^{-1} = R^T
$$
则
$$
M_R = R^T =\left[\begin{matrix}u_x & v_x & w_x & 0 \\u_y & v_y & w_y & 0 \\u_z & v_z & w_z & 0 \\0 & 0 & 0 & 1\end{matrix}\right]
$$
最终代码为：
Matrix44 t = Matrix44::Identity(); |
Projection Matrix
接下来进行透视变换，透视变换的过程为：
$$
视锥体 \xrightarrow{透视矩阵} Clip Space \xrightarrow{透视除法} NDC
$$
其中，Direct3D的NDC(Normalized Device Coordinate)范围为 $[-1, 1] [-1, 1] [0, 1]$
定义屏幕宽高比为aspect，Y方向的视场角为fovy，相机到近平面与远平面的距离分别为$ N $和$ F $
现在，我们试图将视锥体内的顶点线性映射至NDC空间，可得：
$$
x’ = x \frac {1}{z\tan \frac{fovx}{2}} \\
y’ = y \frac {1}{z\tan \frac{fovy}{2}} \\
z’ = \frac{z - N}{F - N}
$$
想要将以上形式写成矩阵乘法形式非常困难，主要原因是$x’$和$y’$之中多出一个$\frac1z$，而且$z’$中的$z$无法提出为因式，于是这里我们需要使用一个十分重要的trick：
由于z值的最终用途是比较顶点深度，并不需要精确映射，仅需满足在$(N, 0)$点到$(F, 1)$点之间保持单调即可。同时，我们先从$x’$和$y’$中提出公因子$\frac1z$，等到矩阵乘法之后再做处理，于是我们可以先将投影矩阵写为：
$$
\left[
\begin{matrix}
\cot(fovx/2) & 0 & 0 & 0 \\
0 & \cot(fovy/2) & 0 & 0 \\
0 & 0 & a & 1 \\
0 & 0 & b & 0
\end{matrix}
\right]
$$
其中$a$和$b$为待求未知数，而其中的1巧妙地将变换前的$z$保存在变换后的$w’$中，此时：
$$
x’ = x \cot (fovx/2) \\
y’ = y \cot (fovy/2) \\
z’ = az + b \\
w’ = z
$$
此时进行透视除法(Perspective provider)，即除以$w’$，得到如下结果：
$$
x’ = \frac 1z x \cot (fovx/2) \\
y’ = \frac 1z y \cot (fovy/2) \\
z’ = \frac {az + b}{z} \\
w’ = 1
$$
其中$x’$和$y’$完全符合要求，$z’ = \frac {az + b}{z}$ 为反比例函数，在$[N, F]$区间内必然单调，符合要求。
将$(N, 0)$和$(F, 1)$点代入该式，得：
$$
\begin{cases}
\frac {aN + b}{N} = 0 \\
\frac {aF + b}{F} = 1
\end{cases}
$$
解得：
$$
\begin{cases}
a = \frac {F}{F - N} \\
b = \frac {-NF}{F - N}
\end{cases}
$$
设$Q = \frac {F}{F - N}$，则最终：
$$
\left[
\begin{matrix}
\cot(fovx/2) & 0 & 0 & 0 \\
0 & \cot(fovy/2) & 0 & 0 \\
0 & 0 & Q & 1 \\
0 & 0 & -NQ & 0
\end{matrix}
\right]
$$
void Transform::SetprojectionMatrix(float fovy, float w, float h, float zn, float zf) |
注意，这一步中最终深度信息z并非线性映射，其函数二阶导小于0，即z经映射后更接近于远平面。
Viewport Matrix
最后一步，视口变换本质上是将顶点从NDC空间映射至视口空间，即为：
$$
NDC \rightarrow Viewport \\
x: [-1, 1] \rightarrow [x_0, x_0 + width] \\
y: [-1, 1] \rightarrow [y_0, y_0 + height] \\
z: [0, 1] \rightarrow [z_{min}, z_{max}]
$$
首先进行区间缩放，缩放矩阵如下，其中Viewport坐标系Y轴竖直向下，需要与之前取反：
$$
\left[
\begin{matrix}
\frac{width}2 & 0 & 0 & 0 \\
0 & -\frac{height}2 & 0 & 0 \\
0 & 0 & z_{max} - z_{min} & 0 \\
0 & 0 & 0 & 1
\end{matrix}
\right]
$$
缩放后的区间为：
$$
x: [-1, 1] \rightarrow [-\frac{width}2, \frac{width}2] \\
y: [-1, 1] \rightarrow [-\frac{height}2, \frac{height}2] \\
z: [0, 1] \rightarrow [0, z_{max} - z_{min}]
$$
之后我们需要平移矩阵与Viewport区间进行对齐：
$$
\left[
\begin{matrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 1 & 0 \\
x_0 + \frac{width}2 & y_0 + \frac{height}2 & z_{min} & 1
\end{matrix}
\right]
$$
两矩阵相乘则得到Viewport Matrix如下：
$$
\left[
\begin{matrix}
\frac{width}2 & 0 & 0 & 0 \\
0 & -\frac{height}2 & 0 & 0 \\
0 & 0 & z_{max} - z_{min} & 0 \\
x_0 + \frac{width}2 & y_0 + \frac{height}2 & z_{min} & 1
\end{matrix}
\right]
$$
void Transform::SetViewportMatrix(float x, float y, float w, float h, float zmin = 0, float zmax = 1) |
Bresenham划线算法
该算法主要用于模型的线框显示，主要思路是通过增量绘制而减少乘除，同时通过乘2的方式避免浮点数出现，最终提高绘制效率。
具体算法细节暂略。
Scanline-Zbuffer多边形填充
该阶段主要进行三角面内部颜色的填充。
首先根据y值对顶点进行排序为v0, v1, v2。
然后分别计算三条线段上x线性插值系数随y增加而累加的步长
float step01 = (v1.pos.y - v0.pos.y) == 0 ? 0 : 1 / (v1.pos.y - v0.pos.y); |
第一段，从v0出发，对线段v0v1和v0v2之间区域进行水平扫描线填充：
int scanlineY = (int)v0.pos.y; |
第二段，从v1处继续，对线段v1v2和v0v2之间区域进行水平扫描线填充：
for (; scanlineY <= v2.pos.y; scanlineY++) |
其中填充像素时，需要将待绘制像素的z值与zBuffer中z值比较：
void RenderBuffer::FillPixel(int x, int y, float z, COLORREF c) |
Phong光照模型 + Phong着色
下面采用Phong光照模型为物体添加光照，以正确显示其在Directional Light下的颜色。
Phong着色与Gourand着色对应，指的是在片元着色器阶段而不是在顶点着色器阶段计算光照。
Directional Light的数据结构为：
class DirectionalLight |
Phong光照模型中
$$
Phong = Ambient + Diffuse + Specular
$$
其中与Directional Light无关的环境光Ambient非常简单：
$$
Ambient = AmbientColor \times AmbientStrength
$$
漫反射Diffuse和镜面反射Specular则需要分别计算其因子，然后与光照颜色和片元颜色相乘：
$$
Diffuse = LightColor \times Fragment Color \times DiffuseFactor \\
Specular = LightColor \times Fragment Color \times SpecularFactor
$$
漫反射因子公式为：
$$
DiffuseFactor = L \cdot N
$$
其中，入射光方向$L$为-DirectionalLight.direction，法向量方向$N$的计算稍后讲解。
镜面反射因子公式为，其中$shininess$越大，高光部分越密集：
$$
SpecularFactor = (R \cdot V)^{shininess}
$$
其中，画图计算易得反射光方向R与视线方向V：
$$
R = 2N(L \cdot N) - L \\
V = CameraPos - FragmentWorldPos
$$
以上代码如下，注意颜色相乘时需要使用哈达玛积(Hadamard Product)使对应分量相乘：
Vector3 DirectionalLight::Apply(Vertex v, Camera camera) |
于是，光照计算中最重要的难点就是法向量变换，由于我们是在WorldSpace下进行光照计算，所以需要通过一个NormalMatrix将模型的法向量信息从LocalSpace变换到WorldSpace。
首先考虑切向量$T$的变换，取三角面两顶点$V1$，$V2$，则有：
$$
T = V1 - V2
$$
记ModelMatrix为$M$，则变换后的切向量：
$$
T’ = V1M - V2M = (V1 - V2)M = TM
$$
即切向量的变换矩阵与顶点完全相同。
而法向量与切向量必然始终垂直，记NormalMatrix为$M_N$则有：
$$
N \cdot T = 0 \\
N’ \cdot T’ = NM_N(TM)^T = NM_NM^TT^T = 0
$$
可得：
$$
M_NM^T = I
$$
即：
$$
M_N = (M^T)^{-1}
$$
由之前的顶点变换矩阵公式，且缩放矩阵$S^T = S$，旋转矩阵正交，法向量平移无效，可知：
$$
M = SRT \\
M_N = (T^TR^TS^T)^{-1} = S^{-1}R
$$
代码如下：
void Transform::UpdateNormalMatrix() |
Matrix44 Matrix44::InverseDiagonal() |
纹理映射（透视校正插值）
之前在模型读取时，我们读取了三角形顶点所对应的纹理坐标texcoord，然而在求三角形内部纹理坐标时，如果直接在屏幕坐标系内对三角形内部的点进行重心插值（或双线性插值，结果相同），则会得到不正确的结果，这主要由之前的透视变换造成，需要进行透视校正插值(Perspective-Correct Interpolation)。
画出一条线段$A(x_0, z_0)B(x_1, z_1)$透视投影在距离相机$d$的近平面上线段$a(u_0, d)b(u_1, d)$的图，其中待插值点$C(x, z)$投影在点$c(u, d)$，易得：
$$
\frac {x_0}{z_0} = \frac {u_0}{d} \\
\frac {x_1}{z_1} = \frac {u_1}{d} \\
\frac {x}{z} = \frac {u}{d}
$$
设线段$ab$上线性插值系数为$s$，线段$AB$上线性插值系数为$t$，则有：
$$
u = u_0 + s(u_1 - u_0) \\
x = x_0 + t(x_1 - x_0) \\
z = z_0 + t(z_1 - z_0)
$$
我们已知系数$s$，需要求正确的插值系数$t$，解得：
$$
t = \frac{sz_0}{sz_0 + (1 - s)z_1}
$$
进一步得：
$$
\frac 1z = (1 - s) \frac 1{z_0} + s \frac 1{z_1}
$$
若$AB$两点分别有属性$I_0I_1$，则插值后得到$C$点属性为：
$$
\frac Iz = (1 - s) \frac {I_0}{z_0} + s \frac {I_1}{z_1}
$$
接下来问题的关键是，如何正确获取顶点在ViewSpace中的$z$值，此处不能直接使用透视除法后的顶点$z$值，因为其仅代表非线性映射后的深度值（我在这个问题上掉坑并找了三天BUG）。前文透视投影部分中提到，变换前的$z$保存在变换后的$w’$中，所以在进行透视除法之前，一定要注意保存一份RHW(Reciprocal of Homogeneous W, $\frac 1{w’}$)值备用。
此时，纹理坐标即为一个需要透视校正插值的属性$I$，于是代码如下，此处$\frac 1z$即等于RHW：
float z = 1 / LerpFloat(v0.rhw, v1.rhw, k); |
全文完，由纸质实习笔记整理而成。