CG 变换与软光栅化器核心算法的工程实现

世界空间变换链#

\underbrace{\text{模型坐标}}_{\text{Local}} \xrightarrow{\mathbf{M}} \underbrace{\text{世界坐标}}_{\text{World}} \xrightarrow{\mathbf{V}} \underbrace{\text{观察坐标}}_{\text{View}} \xrightarrow{\mathbf{P}} \underbrace{\text{裁剪坐标}}_{\text{Clip}} \xrightarrow{\div w} \underbrace{\text{NDC}}_{[-1,1]^3} \xrightarrow{\text{视口}} \underbrace{\text{屏幕坐标}}_{\text{Screen}}

复合矩阵的顺序

\vec{v}_{clip} = \mathbf{P} \cdot \mathbf{V} \cdot \mathbf{M} \cdot \vec{v}_{local}

Vertex Shader 的 out 变量到 Fragment Shader 已经是透视校正插值过的结果。

齐次坐标系统#

齐次坐标的分类：

点的表示： $(x, y, z, 1)$ — 第四个分量为 1
向量的表示： $(x, y, z, 0)$ — 第四个分量为 0
无穷远点： $(x, y, z, 0)$ 且 $(x, y, z) \neq (0, 0, 0)$

Eigen 和 GLM 里某些写法的细微区别

向量表示 Eigen：Vector3f result_point = transformed_point.head<3>() / transformed_point.w();
GL：vec3 result_point = transformed_point.xyz / transformed_point.w;

这里是在将裁剪空间（齐次坐标）转换成归一化设备坐标（NDC）向量的 w 分量为0，不需要进行透视除法

矩阵表示 Eigen：

1
Eigen::Matrix4f create_scaling_matrix(float sx, float sy, float sz)
2
{
3
    Eigen::Matrix4f scaling = Eigen::Matrix4f::Identity();
4
    scaling(0, 0) = sx;
5
    scaling(1, 1) = sy;
6
    scaling(2, 2) = sz;
7
    return scaling;
8
}

GLM：

1
glm::mat4 create_scaling_matrix(float sx, float sy, float sz)
2
{
3
glm::mat4 scaling = glm::mat4(1.0f);
4
scaling[0][0] = sx;  // GLM 用 [列][行] 索引
5
scaling[1][1] = sy;
6
scaling[2][2] = sz;
7
return scaling;
8
}

MVP 变换链的整体实现

MVP 整体实现

1
Matrix4f mvp = projection * view * model;
2

3
for (auto& vertex : vertices) {
4
    Vector4f clip_coord = mvp * Vector4f(vertex.x, vertex.y, vertex.z, 1.0f);
5

6

7
    // 其实在 OpenGL 中，后面的固定管线阶段都由 OpenGL 自动完成
8
    // 透视除法（齐次坐标到 NDC）
9
    Vector3f ndc_coord = clip_coord.head<3>() / clip_coord.w();
10

11
    // 视口变换
12
    Vector3f screen_coord = viewport_transform * ndc_coord;
13
}

几何变换数学原理#

绕任意轴旋转—Rodrigues 旋转公式

R(\vec{n}, \theta) = I\cos\theta + [\vec{n}]_{\times}\sin\theta + \vec{n}\vec{n}^T(1-\cos\theta)

仿射变换 平移不是线性变换。线性变换必须满足 $T(\vec{0}) = \vec{0}$

几个重要矩阵的数学推导*#

法向量变换矩阵#

法向量变换时不能用模型矩阵，这会使法向量扭曲，产生产错误光照。应当使用模型矩阵的逆转置

逆转置矩阵的完整推导

设切向量 $\vec{t}$ 在切平面内，法向量 $\vec{n}$ 与切平面垂直：

\vec{n} \cdot \vec{t} = 0

切向量沿着曲面方向、和顶点差向量一样按 $\mathbf{M}$ 变换：

\vec{t}' = \mathbf{M} \vec{t}

设变换后的法向量为 $\vec{n}' = \mathbf{N} \vec{n}$ ，其中 $\mathbf{N}$ 是我们要求的法向量变换矩阵。变换后仍需保持垂直：

\vec{n}' \cdot \vec{t}' = 0

写成矩阵形式（点积写作行向量乘列向量）：

(\mathbf{N} \vec{n})^T (\mathbf{M} \vec{t}) = 0 \implies \vec{n}^T \mathbf{N}^T \mathbf{M} \vec{t} = 0

原始条件 $\vec{n} \cdot \vec{t} = 0$ 写成矩阵形式为 $\vec{n}^T \vec{t} = 0$ 。两式对任意切向量 $\vec{t}$ 都要成立，意味着：

\mathbf{N}^T \mathbf{M} = \mathbf{I}

解得：

\boxed{\mathbf{N} = (\mathbf{M}^{-1})^T = (\mathbf{M}^T)^{-1}}

观察矩阵 View Transformation#

视图矩阵的作用是将世界空间变为观察空间，更加详细的推导参阅 Learn OpenGL CN

简单推导 [点击展开]

相机的三个属性：位置向量，视线方向，向上向量，向右向量。构建局部坐标系：

1
Vector3f camera_direction   = (eye - target).normalized();   // camera视线方向的反方向（+z 轴，即指向我们自己）
2
Vector3f      = up.cross(camera_direction).normalized(); // 右方向（+x 轴），up 叉乘 foward
3
Vector3f camera_up = camera_direction.cross(camera_right);           // 上方向（y 轴）

几何上：把向量变换到新坐标系里，等价于投影。
因为摄像机的三个轴 R,U,D 是互相垂直的单位向量（正交矩阵），所以转置等于逆。
矩阵乘法顺序是从右往左，所以 LookAt 矩阵先进行的平移，再旋转。

LookAt = \begin{bmatrix} R_x & R_y & R_z & 0 \\ U_x & U_y & U_z & 0 \\ D_x & D_y & D_z & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \cdot \begin{bmatrix} 1 & 0 & 0 & -P_x \\ 0 & 1 & 0 & -P_y \\ 0 & 0 & 1 & -P_z \\ 0 & 0 & 0 & 1 \end{bmatrix}

我们可以将 LookAt 矩阵当做 view 矩阵。幸运的是，GLM已经提供了这些支持。glm::LookAt ，便可以构建一个 view 矩阵

标准透视投影矩阵#

投影矩阵的作用是将观察空间转变为裁剪空间

完整推导点击展开

视锥体参数：

$\text{fov}$ ：垂直视场角（Field of View）
$\text{aspect}$ ：宽高比 $= \text{width} / \text{height}$
$n$ ：近裁剪面距离
$f$ ：远裁剪面距离

在近平面 $z = -n$ 处，视锥体尺寸为：

\text{top} = n \cdot \tan\left(\frac{\text{fov}}{2}\right), \quad \text{bottom} = -\text{top}

\text{right} = \text{top} \cdot \text{aspect}, \quad \text{left} = -\text{right}

坐标下标说明：
- e（eye）：相机坐标系坐标（视图坐标）
- c（clip）：裁剪坐标，投影矩阵输出，齐次形式 $(x_c, y_c, z_c, w_c)$
- n（ndc）：裁剪坐标除以 $w_c$ 的结果，即 $(x_n, y_n, z_n)$
已知视锥体参数（对称视锥，近平面 $z=-n$ ，远平面 $z=-f$ ）： $\text{top} = n \cdot \tan(\frac{\text{fov}}{2}),\quad \text{bottom} = -\text{top}$ $\text{right} = \text{top} \cdot \text{aspect},\quad \text{left} = -\text{right}$
1. 相似三角形投影到近平面
  视锥体内点 $(x_e, y_e, z_e)$ 在 $z=-n$ 上的投影： $x_p = \frac{n \cdot x_e}{-z_e}, \qquad y_p = \frac{n \cdot y_e}{-z_e}$
1. 投影坐标映射到 NDC（ $[-1,1]$ ）
  对称性使线性映射简化： $x_n = \frac{x_p}{\text{right}} = \frac{x_p}{\text{top} \cdot \text{aspect}}$ $y_n = \frac{y_p}{\text{top}}$
1. 构造裁剪坐标（引入 $w_c = -z_e$ ， $w_c$ 记录了点距离相机的深度，因为 OpenGL 相机朝向 -z 轴，所以 $z_e < 0$ ，即 $w_c > 0$ ）
  将 $x_p$ 代入，并设 $w_c = -z_e$ （目的：让 GPU 后续自动除以 $-z_e$ ）： $x_c = x_n \cdot w_c = \frac{n}{\text{top} \cdot \text{aspect}} \cdot x_e$ $y_c = y_n \cdot w_c = \frac{n}{\text{top}} \cdot y_e$ 得到矩阵第 1、2 行系数： $m_{00} = \frac{n}{\text{top} \cdot \text{aspect}} = \frac{1}{\text{aspect} \cdot \tan(\text{fov}/2)}$ $m_{11} = \frac{n}{\text{top}} = \frac{1}{\tan(\text{fov}/2)}$
1. Z 映射到 NDC
  $z_e \in [-n, -f]$ 映射到 $z_n \in [-1, 1]$ 。设 $z_c = A z_e + B$ ， $w_c = -z_e$ ，则 $z_n = z_c / w_c = -A - B / z_e$ 。
  代入边界条件：
  $\begin{cases} z_e = -n \to z_n = -1: & -A + \frac{B}{n} = -1 \\ z_e = -f \to z_n = 1: & -A + \frac{B}{f} = 1 \end{cases}$ 解得： $B = -\frac{2nf}{f-n},\qquad A = -\frac{f+n}{f-n}$ 因此： $z_c = -\frac{f+n}{f-n} z_e - \frac{2fn}{f-n}$ $w_c = -z_e$ 意味着矩阵第 4 行为 $(0,0,-1,0)$ 。
1. 完整投影矩阵（列主序）
  $\begin{pmatrix} \frac{1}{\text{aspect} \cdot \tan(\frac{\text{fov}}{2})} & 0 & 0 & 0 \\ 0 & \frac{1}{\tan(\frac{\text{fov}}{2})} & 0 & 0 \\ 0 & 0 & -\frac{f+n}{f-n} & -\frac{2fn}{f-n} \\ 0 & 0 & -1 & 0 \end{pmatrix}$
注：深度范围若需 $[0,1]$ （Vulkan/DX），第三行系数需调整；本推导基于 OpenGL 默认 $[-1,1]$ 。

代码实现 Eigen

1
Eigen::Matrix4f get_projection_matrix(float eye_fov, float aspect_ratio,
2
                                      float zNear, float zFar)
3
{
4
    // Students will implement this function
5

6
    Eigen::Matrix4f projection = Eigen::Matrix4f::Identity();
7

8
    // TODO: Implement this function
9
    // Create the projection matrix for the given parameters.
10
    // Then return it.
11

12
    float rad = eye_fov * MY_PI / 180.0f;
13

14
    float tanHalffov = std::tan(rad / 2);
15

16
    projection(0, 0) = 1.0f / (aspect_ratio * tanHalffov);
17
    projection(1 ,1) = 1.0f / tanHalffov;
18
    projection(2, 2) = -((zFar + zNear)/(zFar - zNear));
19
    projection(2, 3) = -((2*zFar*zNear) / (zFar - zNear));
20
    projection(3, 2) = -1.0f;
21
    projection(3, 3) = 0.0f;
22

23
    return projection;
24
}

TBN 矩阵#

详细数学推导过程

软光栅化器的实现*#

齐次空间裁剪#

齐次空间裁剪发生在 projection 矩阵之后，透视除法之前。（裁剪空间内）
为什么要裁剪？ 如果一个三角形是跨视锥边界的，那么这个三角形必须被切开，于是要找三角形边与近平面的交点来组成新的三角形。而相机后的点最大特征便是 w < 0。
当然上面的这一切都可以交给 GPU 自己完成，无需自己编码。

光栅化核心#

光栅化在屏幕空间进行，下面是光栅化器的实现方式

根据三角形顶点计算最小包围盒
在包围盒中遍历（扫描转换）：找出哪些像素落在这个三角形内部。
计算重心坐标：算出三角形内部每一个像素点对应的重心坐标 (α,β,γ)
属性插值：根据这些坐标，把顶点着色器输出的颜色、法线、纹理坐标等属性，透视校正插值给覆盖到的每一个像素。
深度测试，更新深度缓冲区，并设置颜色。

光栅化遵循 Top-Left 规则： 只有左边（从下到上的边）和顶边（水平、从左到右的边）上的像素算作”内部”，其他边上的像素算作”外部”

Bresenham 画线算法#

推导过程

设线段起点、终点分别为：

(x_0,y_0),\ (x_1,y_1)

定义：

dx=x_1-x_0,\qquad dy=y_1-y_0

以下推导假设斜率满足：

0\le k=\frac{dy}{dx}\le1

因此每一步固定令 $x=x+1$ ，只需决定 $y$ 是否加 $1$ ，即在 E（下）和 NE（上）两个候选像素之间进行选择。

采用直线的一般形式：

F(x,y)=dy(x-x_0)-dx(y-y_0)

其中：

$F(x,y)=0$ ：点在线上；
$F(x,y)>0$ ：点在线的上侧；
$F(x,y)<0$ ：点在线的下侧。

为了判断应选择哪个像素，取 E 与 NE 的中点：

M=\left(x+1,\ y+\frac12\right)

将其代入直线方程，得到决策变量：

d=F(M)

由于包含 $\frac12$ ，为了避免浮点运算，将其整体乘以 $2$ ：

D=2d= 2dy(x + 1 - x_0) - 2dx(y - y_0) - dx

于是初始决策变量为：

D_0=2dy-dx

之后无需重新计算，只需递推更新：

若选择 E：（x++, y 不变，将新的 x’ = x + 1 代入 D 的表达式，然后对比新旧两式便会发现新式在旧式基础上只多了一个 2dy。同理选 NE 的情况）

D \leftarrow D+2dy

若选择 NE：（x++, y++）

D \leftarrow D+2(dy-dx)

因此整个算法只包含整数的加减运算，无需乘法、除法和浮点计算，这也是 Bresenham 算法高效的原因。

叉积同向法判断点是否在三角形内部#

对于三角形的三条边，分别计算：

\overrightarrow{AB}\times\overrightarrow{AP},\quad \overrightarrow{BC}\times\overrightarrow{BP},\quad \overrightarrow{CA}\times\overrightarrow{CP}

其中：

$P$ 为待判断的点；
二维叉积定义为：

\mathbf{a}\times\mathbf{b}=a_xb_y-a_yb_x

若三个叉积结果同号（全为正或全为负），说明点 $P$ 始终位于三条边的同一侧，因此 $P$ 在三角形内部（或边界上）；若结果有正有负，则点位于三角形外部。

本质：利用叉积判断点是否始终位于三角形三条边的同一侧。

代码实现：
先明确像素坐标的概念，像素坐标一般指这个像素左下角的顶点坐标。但我们判断像素是否在三角形中时，传入的应该是采样点坐标，即（x + 0.5, y + 0.5）

1
static float cross(const Vector3f& a, const Vector3f& b)
2
{
3
    return a.x() * b.y() - a.y() * b.x();
4
}
5

6
static bool insideTriangle(int x, int y, const Vector3f* _v)
7
{
8
    Vector3f p(x, y, 1);
9

10
    Vector3f AB = _v[1] - _v[0];
11
    Vector3f BC = _v[2] - _v[1];
12
    Vector3f CA = _v[0] - _v[2];
13

14
    Vector3f AP = p - v[0];
15
    Vector3f BP = p - v[1];
16
    Vector3f CP = p - v[2];
17

18
    float c1 = cross(AB, AP);
19
    float c2 = cross(BC, BP);
20
    float c3 = cross(CA, CP);
21

22
    return (c1 >= 0 && c2 >= 0 && c3 >= 0) ||
23
           (c1 <= 0 && c2 <= 0 && c3 <= 0);
24
}

重心坐标的计算#

在 GPU 里用 z= 0 的叉积计算面积：

叉积运算：

(x_1, y_1) \times (x_2, y_2) = x_1y_2 - y_1x_2

于是我们设总面积的两倍为 D：

D = (B - C) \times (A - C) = (x_1 - x_2)(y_0 - y_2) - (y_1 - y_2)(x_0 - x_2)

进而：

\alpha = \frac{(x_B - x_P)(y_C - y_P) - (x_C - x_P)(y_B - y_P)}{D}

\beta = \frac{(x_C - x_P)(y_A - y_P) - (x_A - x_P)(y_C - y_P)}{D}

\gamma = 1 - \alpha - \beta

软光栅化器的代码实现

1
std::tuple<float, float, float> barycentric_2d(float x, float y, const Eigen::Vector3f v[3])
2
{
3
    float x0 = v[0].x(), y0 = v[0].y();
4
    float x1 = v[1].x(), y1 = v[1].y();
5
    float x2 = v[2].x(), y2 = v[2].y();
6

7
    float denom = (y1 - y2) * (x0 - x2) + (x2 - x1) * (y0 - y2);
8
    float a = ((y1 - y2) * (x  - x2) + (x2 - x1) * (y  - y2)) / denom;
9
    float b = ((y2 - y0) * (x  - x2) + (x0 - x2) * (y  - y2)) / denom;
10
    float c = 1.0f - a - b;
11
    return {a, b, c};
12
}

透视校正插值#

（透视校正插值便是光栅化最关键的一步），由于透视除法的影响，插值计算不能只是单纯的线性插值，而是需要用到透视校正插值。例如对于属性 uv（或深度）：

GPU 先将每个顶点的 uv 处理成 uv_A/w_A，uv_B/w_B，uv_c/w_c；同时单独处理 w：1/w_A, 1/w_B, 1/w_C
然后根据像素 P 的重心坐标 $(\alpha, \beta, \gamma)$ （P 在三角形里的比例关系），对 uv/w 和 1/w 进行线性插值。（之所以每个属性都要除以 w，是因为化解透视投影的影响，透视投影会对远处的坐标产生压缩，进而导致不能直接进行线性插值。而最后 uv/w 还要通过除以 1/w 来还原真实的 uv，所以我们还需要对 1/w 进行插值计算）

I = \alpha \frac{w v_A}{w_A} + \beta \frac{w v_B}{w_B} + \gamma \frac{w v_C}{w_C}

W = \alpha \frac{1}{w_A} + \beta \frac{1}{w_B} + \gamma \frac{1}{w_C}

最后恢复:

uv_P = \frac{I}{W}

这里的 $uv_p$ 便是 GPU 传给 Fragment Shader 的 uv 值，不过以上的过程在 OpenGL 里都是 GPU 里自动完成的。

软光栅化器的代码实现：

1
// alpha, beta, gamma: 屏幕空间重心坐标
2
// z_a/b/c: 三个顶点在相机空间的深度（正值）
3
// f_a/b/c: 三个顶点上待插值的属性
4
template <typename T>
5
T perspective_correct(float alpha, float beta, float gamma,
6
                      float z_a, float z_b, float z_c,
7
                      const T& f_a, const T& f_b, const T& f_c)
8
{
9
    float inv_z = alpha / z_a + beta / z_b + gamma / z_c;
10
    T     num   = alpha * f_a / z_a + beta * f_b / z_b + gamma * f_c / z_c;
11
    return num / inv_z;
12
}

软光栅化器的核心代码实现#

展开代码

1
// 判断点是否在三角形内部
2
static bool insideTriangle(int x, int y, const Vector3f* _v)
3
{
4
    Vector3f p(x, y, 1);
5

6
    Vector3f AB = _v[1] - _v[0];
7
    Vector3f BC = _v[2] - _v[1];
8
    Vector3f CA = _v[0] - _v[2];
9

10
    Vector3f AP = p - v[0];
11
    Vector3f BP = p - v[1];
12
    Vector3f CP = p - v[2];
13

14
    float c1 = cross(AB, AP);
15
    float c2 = cross(BC, BP);
16
    float c3 = cross(CA, CP);
17

18
    return (c1 >= 0 && c2 >= 0 && c3 >= 0) ||
19
           (c1 <= 0 && c2 <= 0 && c3 <= 0);
20
}
21

22
// 计算重心坐标，返回三个权重
23
std::tuple<float, float, float> barycentric_2d(float x, float y, const Vector3f v[3])
24
{
25
    float x0 = v[0].x(), y0 = v[0].y();
26
    float x1 = v[1].x(), y1 = v[1].y();
27
    float x2 = v[2].x(), y2 = v[2].y();
28

29
    float denom = (y1 - y2) * (x0 - x2) + (x2 - x1) * (y0 - y2);
30
    float a = ((y1 - y2) * (x  - x2) + (x2 - x1) * (y  - y2)) / denom;
31
    float b = ((y2 - y0) * (x  - x2) + (x0 - x2) * (y  - y2)) / denom;
32
    float c = 1.0f - a - b;
33
    return {a, b, c};
34
}
35

36
// 透视校正插值
37
template <typename T> T perspective_correct(float alpha, float beta, float gamma,
38
                      float z_a, float z_b, float z_c,
39
                      const T& f_a, const T& f_b, const T& f_c)
40
{
41
    float inv_z = alpha / z_a + beta / z_b + gamma / z_c;
42
    T     num   = alpha * f_a / z_a + beta * f_b / z_b + gamma * f_c / z_c;
43
    return num / inv_z;
44
}
45

46
// 根据像素的二维坐标得到深度缓冲区对应的数据下标（一维数组）
47
int rst::rasterizer::get_index(int x, int y)
48
{
49
    return (height-1-y)*width + x;
50
}
51

52
// 将一个三角形对象的顶点变成一个长度为3的数据，每个数组元素是4维向量
53
std::array<Vector4f, 3> Triangle::toVector4() const
54
{
55
    std::array<Eigen::Vector4f, 3> res;
56
    std::transform(std::begin(v), std::end(v), res.begin(), [](auto& vec) { return Eigen::Vector4f(vec.x(), vec.y(), vec.z(), 1.f); });
57
    return res;
58
}
59

60
// 屏幕空间光栅化
61
void rst::rasterizer::rasterize_triangle(const Triangle& t)
62
{
63
    auto v = t.toVector4();
64

65
    // bounding box
66
    float minX = std::min({t.v[0].x(), t.v[1].x(), t.v[2].x()});
67
    float maxX = std::max({t.v[0].x(), t.v[1].x(), t.v[2].x()});
68
    float minY = std::min({t.v[0].y(), t.v[1].y(), t.v[2].y()});
69
    float maxY = std::max({t.v[0].y(), t.v[1].y(), t.v[2].y()});
70

71
    // 拓宽上界和下界，防止因浮点数引发边缘漏点
72
    int x_min = floor(minX);
73
    int x_max = ceil(maxX);
74
    int y_min = floor(minY);
75
    int y_max = ceil(maxY);
76

77
    // 遍历每一个像素
78
    for (int i = minX; i <= maxX; i++)
79
    {
80
        for (int j = minY; j <=maxY; j++)
81
        {
82
            // 判断是否在三角形中
83
            if (insideTriangle(i + 0.5f, j + 0.5f, t.v))
84
            {
85
                // 计算三顶点权重
86
                auto [alpha, beta, gamma] = barycentric_2d(i + 0.5f, j + 0.5f, t.v);
87
                // 深度值插值
88
                float z = perspective_correct(alpha, beta, gamma, v[0].z(), v[1].z(), v[2].z(), v[0].w(), v[1].w(), v[2].w());
89
                // 深度测试
90
                int index = get_index(i, j);
91
                if (z < depth_buf[index])
92
                {
93
                    // 更新深度值
94
                    depth_buf[index] = z;
95
                    // 插值其余顶点属性
96
                    auto interpolated_color = perspective_correct(alpha, beta, gamma, t.color[0], t.color[1], t.color[2], v[0].w(), v[1].w(), v[2].w());
97
                    auto interpolated_normal = perspective_correct(alpha,beta,gamma,t.normal[0],t.normal[1],t.normal[2],v[0].w(),v[1].w(),v[2].w());
98
                    auto interpolated_texcoords = perspective_correct(alpha, beta, gamma, t.tex_coords[0], t.tex_coords[1], t.tex_coords[2], v[0].w(), v[1].w(), v[2].w());
99
                    auto interpolated_shadingcoords = perspective_correct(alpha,beta,gamma,view_pos[0],view_pos[1],view_pos[2],v[0].w(),v[1].w(),v[2].w());
100

101

102
                    set_pixel(Eigen::Vector3f(i, j, 1), t.getColor());
103
                }
104
            }
105
        }
106
    }
107
}

萘Naie_Char