计算机图形学笔记(二)：光栅化渲染管线

计算机图形学笔记(二)：光栅化渲染管线#

这部分是整个图形学的核心，也是最实用的部分。做Assignment的时候深深感受到，理解渲染管线的每个步骤对写出正确的代码太重要了。从顶点数据到最终像素，每一步都有它存在的意义。

目录#

图形渲染管线概述 - 现代GPU管线架构
顶点处理与变换 - MVP变换链与顶点着色器
图元装配与裁剪 - 视锥体裁剪算法详解
光栅化算法详解 - 三角形光栅化与重心坐标
深度测试与隐藏面消除 - Z-Buffer算法与优化
光照模型与着色 - Phong模型到PBR的演进
纹理映射技术 - Mipmap、法线贴图与环境映射

图形渲染管线概述#

5.1 渲染管线的整体架构#

5.1.1 管线阶段划分#

应用阶段（Application Stage）：

场景管理：构建场景图
视锥体裁剪：剔除不可见物体
细节层次控制：根据距离选择模型精度
动画更新：骨骼动画、变形动画

几何阶段（Geometry Stage）：

顶点着色：变换顶点位置
投影：3D到2D的映射
裁剪：移除视锥体外的几何体
屏幕映射：NDC到屏幕坐标

光栅化阶段（Rasterization Stage）：

三角形设置：计算边方程
三角形遍历：确定覆盖的像素
像素着色：计算每个像素的颜色
合并：深度测试、混合等

5.1.2 坐标系变换流程的数学推导#

完整的变换链#

坐标空间序列： $\text{模型坐标} \xrightarrow{M} \text{世界坐标} \xrightarrow{V} \text{观察坐标} \xrightarrow{P} \text{裁剪坐标} \xrightarrow{\text{透视除法}} \text{NDC} \xrightarrow{S} \text{屏幕坐标}$

MVP变换的数学表示#

复合变换矩阵： $\mathbf{V}_{clip} = \mathbf{P} \cdot \mathbf{V} \cdot \mathbf{M} \cdot \mathbf{V}_{model}$

其中：

$\mathbf{M} \in \mathbb{R}^{4 \times 4}$ ：模型变换矩阵（Model Matrix）
$\mathbf{V} \in \mathbb{R}^{4 \times 4}$ ：视图变换矩阵（View Matrix）
$\mathbf{P} \in \mathbb{R}^{4 \times 4}$ ：投影变换矩阵（Projection Matrix）

透视除法（Perspective Division）#

齐次坐标到NDC的转换：

透视除法： $(x_{clip}, y_{clip}, z_{clip}, w_{clip}) \rightarrow (x_{clip}/w_{clip}, y_{clip}/w_{clip}, z_{clip}/w_{clip})$

$\mathbf{V}_{ndc} = (x_{ndc}, y_{ndc}, z_{ndc}) = (x_{clip}/w_{clip}, y_{clip}/w_{clip}, z_{clip}/w_{clip})$

NDC范围： $x_{ndc}, y_{ndc}, z_{ndc} \in [-1, 1]$

视口变换（Viewport Transform）#

从NDC到屏幕坐标：

视口变换矩阵 $\mathbf{S}$ 的结构：

第1行： $(\frac{w}{2}, 0, 0, \frac{w}{2})$
第2行： $(0, \frac{h}{2}, 0, \frac{h}{2})$
第3行： $(0, 0, \frac{d}{2}, \frac{d}{2})$
第4行： $(0, 0, 0, 1)$

$\mathbf{S} = \text{viewport transformation matrix}$

其中：

$w$ ：屏幕宽度（像素）
$h$ ：屏幕高度（像素）
$d$ ：深度缓冲区范围

最终屏幕坐标：

对于NDC坐标 $(x_{ndc}, y_{ndc}, z_{ndc}, 1)$ ，屏幕坐标为： $\mathbf{V}_{screen} = \left(\frac{w}{2}(x_{ndc} + 1), \frac{h}{2}(y_{ndc} + 1), \frac{d}{2}(z_{ndc} + 1), 1\right)$

变换的几何意义#

坐标系的右手/左手约定：

OpenGL：右手坐标系， $z_{ndc} \in [-1, 1]$
DirectX：左手坐标系， $z_{ndc} \in [0, 1]$

深度值的非线性分布：透视投影后，深度值在近平面附近密集，远平面附近稀疏： $z_{ndc} = \frac{-(f+n)z - 2fn}{z(f-n)}$

5.2 可编程着色器架构#

5.2.1 顶点着色器（Vertex Shader）#

主要功能：

顶点位置变换
法向量变换
纹理坐标传递
光照计算（Gouraud着色）

典型顶点着色器：

1
# version 330 core
2

3
layout (location = 0) in vec3 aPos;
4
layout (location = 1) in vec3 aNormal;
5
layout (location = 2) in vec2 aTexCoord;
6

7
uniform mat4 model;
8
uniform mat4 view;
9
uniform mat4 projection;
10

11
out vec3 FragPos;
12
out vec3 Normal;
13
out vec2 TexCoord;
14

15
void main() {
16
    FragPos = vec3(model * vec4(aPos, 1.0));
17
    Normal = mat3(transpose(inverse(model))) * aNormal;
18
    TexCoord = aTexCoord;
19

20
    gl_Position = projection * view * vec4(FragPos, 1.0);
21
}

5.2.2 片段着色器（Fragment Shader）#

主要功能：

像素颜色计算
纹理采样
光照计算（Phong着色）
特效处理

典型片段着色器：

1
# version 330 core
2

3
in vec3 FragPos;
4
in vec3 Normal;
5
in vec2 TexCoord;
6

7
uniform sampler2D texture_diffuse1;
8
uniform vec3 lightPos;
9
uniform vec3 viewPos;
10

11
out vec4 FragColor;
12

13
void main() {
14
    // 环境光
15
    float ambientStrength = 0.1;
16
    vec3 ambient = ambientStrength * vec3(1.0);
17

18
    // 漫反射
19
    vec3 norm = normalize(Normal);
20
    vec3 lightDir = normalize(lightPos - FragPos);
21
    float diff = max(dot(norm, lightDir), 0.0);
22
    vec3 diffuse = diff * vec3(1.0);
23

24
    // 镜面反射
25
    float specularStrength = 0.5;
26
    vec3 viewDir = normalize(viewPos - FragPos);
27
    vec3 reflectDir = reflect(-lightDir, norm);
28
    float spec = pow(max(dot(viewDir, reflectDir), 0.0), 32);
29
    vec3 specular = specularStrength * spec * vec3(1.0);
30

31
    vec3 result = (ambient + diffuse + specular) * texture(texture_diffuse1, TexCoord).rgb;
32
    FragColor = vec4(result, 1.0);
33
}

顶点处理与变换#

6.1 顶点属性与数据结构#

6.1.1 顶点数据组织#

基本顶点属性：

1
struct Vertex {
2
    Vector3f position;    // 位置
3
    Vector3f normal;      // 法向量
4
    Vector2f texCoord;    // 纹理坐标
5
    Vector3f tangent;     // 切向量
6
    Vector3f bitangent;   // 副切向量
7
    Vector4f color;       // 顶点颜色
8
};

顶点缓冲对象（VBO）：

1
// 创建并绑定VBO
2
unsigned int VBO;
3
glGenBuffers(1, &VBO);
4
glBindBuffer(GL_ARRAY_BUFFER, VBO);
5
glBufferData(GL_ARRAY_BUFFER, vertices.size() * sizeof(Vertex),
6
             vertices.data(), GL_STATIC_DRAW);

6.1.2 顶点数组对象（VAO）#

VAO的作用：

存储顶点属性配置
简化渲染调用
提高渲染效率

VAO配置示例：

1
unsigned int VAO;
2
glGenVertexArrays(1, &VAO);
3
glBindVertexArray(VAO);
4

5
// 位置属性
6
glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, sizeof(Vertex), (void*)0);
7
glEnableVertexAttribArray(0);
8

9
// 法向量属性
10
glVertexAttribPointer(1, 3, GL_FLOAT, GL_FALSE, sizeof(Vertex),
11
                      (void*)offsetof(Vertex, normal));
12
glEnableVertexAttribArray(1);
13

14
// 纹理坐标属性
15
glVertexAttribPointer(2, 2, GL_FLOAT, GL_FALSE, sizeof(Vertex),
16
                      (void*)offsetof(Vertex, texCoord));
17
glEnableVertexAttribArray(2);

6.2 变换矩阵的计算与优化#

6.2.1 法向量变换的数学推导#

问题的提出#

核心问题：法向量不能直接用模型变换矩阵进行变换，否则会破坏垂直关系。

数学推导过程#

平面的隐式表示：设平面方程为： $\vec{n} \cdot \vec{p} = d$

其中 $\vec{n}$ 是法向量， $\vec{p}$ 是平面上的点， $d$ 是常数。

变换后的约束条件：变换后的平面方程应为： $\vec{n}' \cdot \vec{p}' = d$

点的变换关系：已知点的变换为： $\vec{p}' = \mathbf{M}\vec{p}$

法向量变换的推导：为保持垂直关系，需要： $\vec{n}' \cdot (\mathbf{M}\vec{p}) = \vec{n} \cdot \vec{p}$

将其写成矩阵形式： $(\vec{n}')^T \mathbf{M}\vec{p} = \vec{n}^T \vec{p}$

由于这对所有 $\vec{p}$ 都成立，因此： $(\vec{n}')^T \mathbf{M} = \vec{n}^T$

转置得到： $\mathbf{M}^T \vec{n}' = \vec{n}$

解得： $\vec{n}' = (\mathbf{M}^T)^{-1}\vec{n} = (\mathbf{M}^{-1})^T\vec{n}$

法向量变换矩阵#

法向量变换矩阵： $\mathbf{N} = (\mathbf{M}^{-1})^T$

其中 $\mathbf{M}$ 是模型变换矩阵的左上角 $3 \times 3$ 子矩阵。

特殊情况：

正交变换：当 $\mathbf{M}$ 是正交矩阵时， $\mathbf{M}^{-1} = \mathbf{M}^T$ ，因此 $\mathbf{N} = \mathbf{M}$
均匀缩放：当 $\mathbf{M} = s\mathbf{I}$ 时， $\mathbf{N} = \frac{1}{s}\mathbf{I}$

工程实现#

1
// 计算法向量变换矩阵
2
Matrix3f compute_normal_matrix(const Matrix4f& model_matrix) {
3
    Matrix3f upper_left = model_matrix.block<3,3>(0,0);
4
    return upper_left.inverse().transpose();
5
}
6

7
// 应用法向量变换
8
Vector3f transform_normal(const Vector3f& normal, const Matrix4f& model_matrix) {
9
    Matrix3f normal_matrix = compute_normal_matrix(model_matrix);
10
    Vector3f transformed = normal_matrix * normal;
11
    return transformed.normalized();  // 重新归一化
12
}

6.2.2 变换矩阵的TRS分解#

TRS分解的数学理论#

分解定理：任何非奇异的仿射变换矩阵 $\mathbf{M}$ 都可以唯一分解为： $\mathbf{M} = \mathbf{T} \cdot \mathbf{R} \cdot \mathbf{S}$

其中：

$\mathbf{T}$ ：平移矩阵（Translation）
$\mathbf{R}$ ：旋转矩阵（Rotation）
$\mathbf{S}$ ：缩放矩阵（Scaling）

分解算法的数学推导#

矩阵结构分析：

变换矩阵的分块结构：

左上角 $3 \times 3$ 子矩阵 $\mathbf{A}$ ：包含旋转和缩放
右上角列向量 $\vec{t}$ ：平移向量
左下角行向量 $\vec{0}^T$ ：零向量
右下角标量：1

$\mathbf{M} = \text{block matrix structure}$

其中 $\mathbf{A} \in \mathbb{R}^{3 \times 3}$ 是线性变换部分， $\vec{t} \in \mathbb{R}^3$ 是平移向量。

第一步：提取平移 $\vec{t} = \mathbf{M}_{[0:3,3]}$

第二步：分解线性变换 对矩阵 $\mathbf{A}$ 进行极分解： $\mathbf{A} = \mathbf{R} \cdot \mathbf{S}$

缩放因子计算： $s_x = \|\mathbf{A}_{[:,0]}\|, \quad s_y = \|\mathbf{A}_{[:,1]}\|, \quad s_z = \|\mathbf{A}_{[:,2]}\|$

旋转矩阵提取：

旋转矩阵通过归一化列向量获得： $\mathbf{R} = [\frac{\mathbf{A}_{[:,0]}}{s_x}, \frac{\mathbf{A}_{[:,1]}}{s_y}, \frac{\mathbf{A}_{[:,2]}}{s_z}]$

其中每一列都是归一化后的旋转轴。

处理反射变换#

行列式检查：如果 $\det(\mathbf{A}) < 0$ ，说明包含反射变换： $s_z = -s_z$ $\mathbf{R}_{[:,2]} = -\mathbf{R}_{[:,2]}$

完整的分解实现#

1
struct Transform {
2
    Vector3f translation;
3
    Quaternionf rotation;
4
    Vector3f scale;
5
};
6

7
Transform decompose_matrix(const Matrix4f& matrix) {
8
    Transform result;
9

10
    // 1. 提取平移
11
    result.translation = matrix.block<3,1>(0,3);
12

13
    // 2. 提取线性变换部分
14
    Matrix3f A = matrix.block<3,3>(0,0);
15

16
    // 3. 计算缩放因子
17
    result.scale.x() = A.col(0).norm();
18
    result.scale.y() = A.col(1).norm();
19
    result.scale.z() = A.col(2).norm();
20

21
    // 4. 处理反射（负行列式）
22
    if (A.determinant() < 0) {
23
        result.scale.z() = -result.scale.z();
24
    }
25

26
    // 5. 提取旋转矩阵
27
    Matrix3f R;
28
    R.col(0) = A.col(0) / result.scale.x();
29
    R.col(1) = A.col(1) / result.scale.y();
30
    R.col(2) = A.col(2) / result.scale.z();
31

32
    // 6. 转换为四元数
33
    result.rotation = Quaternionf(R);
34

35
    return result;
36
}
37
    rotation_matrix.col(1) = upper_left.col(1) / scale.y();
38
    rotation_matrix.col(2) = upper_left.col(2) / scale.z();
39

40
    // 转换为四元数
41
    rotation = Quaternionf(rotation_matrix);
42
}

图元装配与裁剪#

7.1 图元装配过程#

7.1.1 图元类型#

基本图元：

点（Points）
线段（Lines）
三角形（Triangles）

扩展图元：

线条带（Line Strip）
三角形带（Triangle Strip）
三角形扇（Triangle Fan）

7.1.2 索引缓冲对象（EBO）#

作用：减少顶点数据冗余

示例：

1
// 顶点数据（4个顶点定义矩形）
2
float vertices[] = {
3
     0.5f,  0.5f, 0.0f,  // 右上
4
     0.5f, -0.5f, 0.0f,  // 右下
5
    -0.5f, -0.5f, 0.0f,  // 左下
6
    -0.5f,  0.5f, 0.0f   // 左上
7
};
8

9
// 索引数据（2个三角形）
10
unsigned int indices[] = {
11
    0, 1, 3,  // 第一个三角形
12
    1, 2, 3   // 第二个三角形
13
};
14

15
// 创建EBO
16
unsigned int EBO;
17
glGenBuffers(1, &EBO);
18
glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, EBO);
19
glBufferData(GL_ELEMENT_ARRAY_BUFFER, sizeof(indices), indices, GL_STATIC_DRAW);

7.2 裁剪算法#

7.2.1 视锥体裁剪的数学理论#

视锥体的数学定义#

标准化设备坐标（NDC）中的视锥体：在齐次裁剪坐标系中，视锥体由六个半空间定义：

六个裁剪平面：

左平面： $x \geq -w$
右平面： $x \leq w$
下平面： $y \geq -w$
上平面： $y \leq w$
近平面： $z \geq -w$
远平面： $z \leq w$

其中 $(x, y, z, w)$ 是齐次裁剪坐标。

点的裁剪测试#

数学条件：点 $\mathbf{P}_{clip} = (x, y, z, w)$ 在视锥体内当且仅当： $-w \leq x \leq w \quad \text{且} \quad -w \leq y \leq w \quad \text{且} \quad -w \leq z \leq w$

几何意义：这等价于透视除法后的NDC坐标满足： $-1 \leq \frac{x}{w} \leq 1, \quad -1 \leq \frac{y}{w} \leq 1, \quad -1 \leq \frac{z}{w} \leq 1$

工程实现#

1
bool is_inside_frustum(const Vector4f& clip_pos) {
2
    float w = clip_pos.w();
3

4
    // 检查w分量的有效性
5
    if (w <= 0) return false;  // 在相机后方
6

7
    return (clip_pos.x() >= -w && clip_pos.x() <= w &&
8
            clip_pos.y() >= -w && clip_pos.y() <= w &&
9
            clip_pos.z() >= -w && clip_pos.z() <= w);
10
}
11

12
// 计算点到裁剪平面的距离
13
float distance_to_plane(const Vector4f& point, int plane_index) {
14
    switch (plane_index) {
15
        case 0: return point.w() + point.x();  // 左平面
16
        case 1: return point.w() - point.x();  // 右平面
17
        case 2: return point.w() + point.y();  // 下平面
18
        case 3: return point.w() - point.y();  // 上平面
19
        case 4: return point.w() + point.z();  // 近平面
20
        case 5: return point.w() - point.z();  // 远平面
21
        default: return 0.0f;
22
    }
23
}

7.2.2 Sutherland-Hodgman裁剪算法#

算法的数学原理#

基本思想：对多边形逐个裁剪平面进行裁剪，每次裁剪产生一个新的多边形，直到所有裁剪平面都处理完毕。

数学模型：设多边形的顶点序列为 $\{V_0, V_1, \ldots, V_{n-1}\}$ ，裁剪平面方程为： $\pi: \mathbf{n} \cdot \mathbf{p} + d = 0$

其中 $\mathbf{n}$ 是平面法向量， $d$ 是距离参数。

点与平面的位置关系#

符号距离函数： $\text{dist}(\mathbf{p}, \pi) = \mathbf{n} \cdot \mathbf{p} + d$

位置判断：

$\text{dist}(\mathbf{p}, \pi) > 0$ ：点在平面正侧（内部）
$\text{dist}(\mathbf{p}, \pi) = 0$ ：点在平面上
$\text{dist}(\mathbf{p}, \pi) < 0$ ：点在平面负侧（外部）

线段与平面的交点计算#

参数方程：线段 $\overline{P_1P_2}$ 的参数方程为： $\mathbf{L}(t) = \mathbf{P_1} + t(\mathbf{P_2} - \mathbf{P_1}), \quad t \in [0, 1]$

交点参数求解：将参数方程代入平面方程： $\mathbf{n} \cdot [\mathbf{P_1} + t(\mathbf{P_2} - \mathbf{P_1})] + d = 0$

解得： $t = -\frac{\mathbf{n} \cdot \mathbf{P_1} + d}{\mathbf{n} \cdot (\mathbf{P_2} - \mathbf{P_1})} = -\frac{\text{dist}(\mathbf{P_1}, \pi)}{\mathbf{n} \cdot (\mathbf{P_2} - \mathbf{P_1})}$

交点坐标： $\mathbf{P}_{intersection} = \mathbf{P_1} + t(\mathbf{P_2} - \mathbf{P_1})$

算法实现#

1
class SutherlandHodgmanClipper {
2
private:
3
    struct Plane {
4
        Vector3f normal;
5
        float d;
6

7
        float distance(const Vector3f& point) const {
8
            return normal.dot(point) + d;
9
        }
10

11
        bool is_inside(const Vector3f& point) const {
12
            return distance(point) >= 0;
13
        }
14
    };
15

16
public:
17
    std::vector<Vector3f> clip_polygon(const std::vector<Vector3f>& polygon,
18
                                      const std::vector<Plane>& planes) {
19
        std::vector<Vector3f> input = polygon;
20

21
        for (const auto& plane : planes) {
22
            std::vector<Vector3f> output;
23

24
            if (input.empty()) break;
25

26
            Vector3f prev_vertex = input.back();
27
            bool prev_inside = plane.is_inside(prev_vertex);
28

29
            for (const auto& curr_vertex : input) {
30
                bool curr_inside = plane.is_inside(curr_vertex);
31

32
                if (curr_inside) {
33
                    if (!prev_inside) {
34
                        // 从外部进入：添加交点
35
                        Vector3f intersection = compute_intersection(prev_vertex, curr_vertex, plane);
36
                        output.push_back(intersection);
37
                    }
38
                    // 添加当前顶点
39
                    output.push_back(curr_vertex);
40
                } else if (prev_inside) {
41
                    // 从内部离开：只添加交点
42
                    Vector3f intersection = compute_intersection(prev_vertex, curr_vertex, plane);
43
                    output.push_back(intersection);
44
                }
45
                // 两点都在外部：不添加任何点
46

47
                prev_vertex = curr_vertex;
48
                prev_inside = curr_inside;
49
            }
50

51
            input = output;
52
        }
53

54
        return input;
55
    }
56

57
private:
58
    Vector3f compute_intersection(const Vector3f& p1, const Vector3f& p2, const Plane& plane) {
59
        Vector3f direction = p2 - p1;
60
        float denominator = plane.normal.dot(direction);
61

62
        if (std::abs(denominator) < 1e-6f) {
63
            return p1;  // 线段平行于平面
64
        }
65

66
        float t = -plane.distance(p1) / denominator;
67
        return p1 + t * direction;
68
    }
69
};

光栅化算法详解#

8.1 三角形光栅化基础#

8.1.1 扫描线算法的数学原理#

算法概述与数学基础#

基本思想：扫描线算法通过逐行扫描的方式填充三角形，是经典的多边形光栅化方法。

核心数学原理：

直线方程：利用直线的参数方程或隐式方程计算交点
区间填充：在每条扫描线上填充左右边界之间的像素
增量计算：利用相邻扫描线的相关性优化计算

直线方程与交点计算#

参数直线方程：连接点 $P_1(x_1, y_1)$ 和 $P_2(x_2, y_2)$ 的直线：

参数方程（ $t \in [0, 1]$ ）： $x(t) = x_1 + t(x_2 - x_1)$ $y(t) = y_1 + t(y_2 - y_1)$

扫描线交点计算：给定水平扫描线 $y = y_s$ ，求解参数 $t$ ： $y_s = y_1 + t(y_2 - y_1) \implies t = \frac{y_s - y_1}{y_2 - y_1}$

代入x方程得交点： $x_s = x_1 + \frac{y_s - y_1}{y_2 - y_1}(x_2 - x_1) = x_1 + \frac{(y_s - y_1)(x_2 - x_1)}{y_2 - y_1}$

增量形式优化：相邻扫描线的x坐标增量： $\Delta x = \frac{x_2 - x_1}{y_2 - y_1} = \frac{\Delta x}{\Delta y}$

因此： $x_{i+1} = x_i + \Delta x$

边方程方法（Half-Space Test）#

隐式直线方程：对于直线 $AB$ ，隐式方程为： $f_{AB}(x, y) = (y_A - y_B)x + (x_B - x_A)y + (x_A y_B - x_B y_A) = 0$

半空间测试：点 $P(x, y)$ 相对于有向直线 $AB$ 的位置：

$f_{AB}(P) > 0$ ：点在直线左侧
$f_{AB}(P) = 0$ ：点在直线上
$f_{AB}(P) < 0$ ：点在直线右侧

三角形内部判断：对于逆时针方向的三角形 $ABC$ ，点 $P$ 在三角形内当且仅当：

边函数定义：

$E_1(P) = (y_A - y_B)x_P + (x_B - x_A)y_P + (x_A y_B - x_B y_A) \geq 0$

$E_2(P) = (y_B - y_C)x_P + (x_C - x_B)y_P + (x_B y_C - x_C y_B) \geq 0$

$E_3(P) = (y_C - y_A)x_P + (x_A - x_C)y_P + (x_C y_A - x_A y_C) \geq 0$

优化的扫描线实现#

1
void optimized_scanline_triangle(const Triangle& t) {
2
    auto vertices = t.toVector4();
3

4
    // 按y坐标排序
5
    std::sort(vertices.begin(), vertices.end(),
6
              [](const Vector4f& a, const Vector4f& b) { return a.y() < b.y(); });
7

8
    Vector4f v0 = vertices[0], v1 = vertices[1], v2 = vertices[2];
9

10
    // 计算边的增量
11
    float dx01 = (v1.x() - v0.x()) / (v1.y() - v0.y());
12
    float dx02 = (v2.x() - v0.x()) / (v2.y() - v0.y());
13
    float dx12 = (v2.x() - v1.x()) / (v2.y() - v1.y());
14

15
    // 上半部分三角形 (v0 到 v1)
16
    float x_left = v0.x(), x_right = v0.x();
17
    for (int y = v0.y(); y <= v1.y(); y++) {
18
        fill_horizontal_line(x_left, x_right, y, t.getColor());
19
        x_left += dx01;
20
        x_right += dx02;
21
    }
22

23
    // 下半部分三角形 (v1 到 v2)
24
    x_left = v1.x();  // 重新设置左边界
25
    for (int y = v1.y(); y <= v2.y(); y++) {
26
        fill_horizontal_line(x_left, x_right, y, t.getColor());
27
        x_left += dx12;
28
        x_right += dx02;
29
    }
30
}

8.1.2 重心坐标系统（Barycentric Coordinates）#

重心坐标的数学基础#

定义：对于三角形 $\triangle ABC$ 和平面内任意点 $P$ ，重心坐标是一组权重 $(\alpha, \beta, \gamma)$ ，使得： $P = \alpha A + \beta B + \gamma C$ $\alpha + \beta + \gamma = 1$

几何解释：

$\alpha, \beta, \gamma$ 分别表示点 $P$ 相对于顶点 $A, B, C$ 的”重量”
当 $\alpha, \beta, \gamma \geq 0$ 时，点 $P$ 在三角形内部
重心坐标提供了三角形内任意点的唯一表示

重心坐标的计算方法#

面积比方法：重心坐标等于子三角形面积与总三角形面积的比值： $\alpha = \frac{S_{\triangle PBC}}{S_{\triangle ABC}}, \quad \beta = \frac{S_{\triangle APC}}{S_{\triangle ABC}}, \quad \gamma = \frac{S_{\triangle APB}}{S_{\triangle ABC}}$

向量叉积计算：利用叉积计算有向面积： $\alpha = \frac{(\vec{BP} \times \vec{BC}) \cdot \vec{n}}{(\vec{BA} \times \vec{BC}) \cdot \vec{n}}$ 其中 $\vec{n}$ 是三角形的法向量。

GAMES101 Assignment 2实现：

1
// 基于面积的重心坐标计算
2
static std::tuple<float, float, float> computeBarycentric2D(float x, float y, const Vector3f* v) {
3
    float c1 = (x*(v[1].y() - v[2].y()) + (v[2].x() - v[1].x())*y + v[1].x()*v[2].y() - v[2].x()*v[1].y()) /
4
               (v[0].x()*(v[1].y() - v[2].y()) + (v[2].x() - v[1].x())*v[0].y() + v[1].x()*v[2].y() - v[2].x()*v[1].y());
5
    float c2 = (x*(v[2].y() - v[0].y()) + (v[0].x() - v[2].x())*y + v[2].x()*v[0].y() - v[0].x()*v[2].y()) /
6
               (v[1].x()*(v[2].y() - v[0].y()) + (v[0].x() - v[2].x())*v[1].y() + v[2].x()*v[0].y() - v[0].x()*v[2].y());
7
    float c3 = (x*(v[0].y() - v[1].y()) + (v[1].x() - v[0].x())*y + v[0].x()*v[1].y() - v[1].x()*v[0].y()) /
8
               (v[2].x()*(v[0].y() - v[1].y()) + (v[1].x() - v[0].x())*v[2].y() + v[0].x()*v[1].y() - v[1].x()*v[0].y());
9
    return {c1, c2, c3};
10
}
11

12
// 优化版本：避免重复计算
13
Vector3f barycentric_coordinates(const Vector2f& A, const Vector2f& B,
14
                                const Vector2f& C, const Vector2f& P) {
15
    Vector2f v0 = C - A;
16
    Vector2f v1 = B - A;
17
    Vector2f v2 = P - A;
18

19
    float dot00 = v0.dot(v0);
20
    float dot01 = v0.dot(v1);
21
    float dot02 = v0.dot(v2);
22
    float dot11 = v1.dot(v1);
23
    float dot12 = v1.dot(v2);
24

25
    float inv_denom = 1.0f / (dot00 * dot11 - dot01 * dot01);
26
    float u = (dot11 * dot02 - dot01 * dot12) * inv_denom;
27
    float v = (dot00 * dot12 - dot01 * dot02) * inv_denom;
28

29
    return Vector3f(1.0f - u - v, v, u);  // $(\alpha, \beta, \gamma)$
30
}

重心坐标的几何意义#

面积解释：

$\alpha = \frac{\text{Area}(\triangle PBC)}{\text{Area}(\triangle ABC)}$ ：点P到边BC的”距离权重”
$\beta = \frac{\text{Area}(\triangle APC)}{\text{Area}(\triangle ABC)}$ ：点P到边AC的”距离权重”
$\gamma = \frac{\text{Area}(\triangle APB)}{\text{Area}(\triangle ABC)}$ ：点P到边AB的”距离权重”

边界情况：

$P = A$ ： $(\alpha, \beta, \gamma) = (1, 0, 0)$
$P = B$ ： $(\alpha, \beta, \gamma) = (0, 1, 0)$
$P = C$ ： $(\alpha, \beta, \gamma) = (0, 0, 1)$
$P$ 在边 $AB$ 上： $\gamma = 0$ ， $\alpha + \beta = 1$

重心坐标的应用#

1. 点在三角形内部判断：

1
bool is_inside_triangle(float alpha, float beta, float gamma) {
2
    return alpha >= 0 && beta >= 0 && gamma >= 0;
3
}

2. 属性插值：对于三角形顶点的任意属性（颜色、深度、纹理坐标等），可以用重心坐标进行插值：

1
// 深度插值
2
float interpolated_depth = alpha * depth_A + beta * depth_B + gamma * depth_C;
3

4
// 颜色插值
5
Vector3f interpolated_color = alpha * color_A + beta * color_B + gamma * color_C;
6

7
// 纹理坐标插值
8
Vector2f interpolated_uv = alpha * uv_A + beta * uv_B + gamma * uv_C;

8.2 属性插值#

8.2.1 透视校正插值#

问题：屏幕空间的线性插值在透视投影下不正确

解决方案：在齐次坐标空间进行插值

透视正确插值公式：

深度倒数插值： $\frac{1}{z} = \frac{\alpha}{z_1} + \frac{\beta}{z_2} + \frac{\gamma}{z_3}$

纹理坐标u分量插值： $\frac{u}{z} = \alpha\frac{u_1}{z_1} + \beta\frac{u_2}{z_2} + \gamma\frac{u_3}{z_3}$

纹理坐标v分量插值： $\frac{v}{z} = \alpha\frac{v_1}{z_1} + \beta\frac{v_2}{z_2} + \gamma\frac{v_3}{z_3}$

最终纹理坐标： $u = \frac{u/z}{1/z}, \quad v = \frac{v/z}{1/z}$

代码实现：

1
Vector2f perspective_correct_interpolation(const Vector3f& bary,
2
                                          const Vector2f& uv1, float z1,
3
                                          const Vector2f& uv2, float z2,
4
                                          const Vector2f& uv3, float z3) {
5
    float inv_z = bary.x() / z1 + bary.y() / z2 + bary.z() / z3;
6
    float u_over_z = bary.x() * uv1.x() / z1 +
7
                     bary.y() * uv2.x() / z2 +
8
                     bary.z() * uv3.x() / z3;
9
    float v_over_z = bary.x() * uv1.y() / z1 +
10
                     bary.y() * uv2.y() / z2 +
11
                     bary.z() * uv3.y() / z3;
12

13
    return Vector2f(u_over_z / inv_z, v_over_z / inv_z);
14
}

8.2.2 现代GPU并行光栅化算法#

Tile-Based光栅化#

基本思想：现代GPU将屏幕分割成小的tile（通常8×8或16×16像素），每个tile并行处理。

算法流程：

Tile分割：将屏幕分成规则的tile网格
三角形分配：确定每个三角形覆盖哪些tile
并行光栅化：每个tile独立进行光栅化
结果合并：将各tile的结果合并到最终图像

数学模型：对于tile $(i, j)$ ，其屏幕坐标范围为：

Tile坐标计算： $x_{min} = i \times \text{tile\_size}$ $x_{max} = (i + 1) \times \text{tile\_size} - 1$ $y_{min} = j \times \text{tile\_size}$ $y_{max} = (j + 1) \times \text{tile\_size} - 1$

三角形-Tile相交测试：使用分离轴定理（SAT）或包围盒测试快速判断三角形是否与tile相交。

GPU Warp/Wavefront并行模型#

SIMD执行模型： GPU以warp（NVIDIA）或wavefront（AMD）为单位执行，通常包含32个线程。

并行光栅化实现：

1
// GPU kernel伪代码
2
__global__ void rasterize_triangle_kernel(Triangle* triangles, int num_triangles,
3
                                         float* framebuffer, float* depth_buffer,
4
                                         int width, int height) {
5
    int pixel_x = blockIdx.x * blockDim.x + threadIdx.x;
6
    int pixel_y = blockIdx.y * blockDim.y + threadIdx.y;
7

8
    if (pixel_x >= width || pixel_y >= height) return;
9

10
    float min_depth = FLT_MAX;
11
    Vector3f final_color(0, 0, 0);
12

13
    // 遍历所有三角形
14
    for (int i = 0; i < num_triangles; i++) {
15
        Triangle& tri = triangles[i];
16

17
        // 计算重心坐标
18
        auto [alpha, beta, gamma] = computeBarycentric2D(pixel_x, pixel_y, tri.vertices);
19

20
        // 内部测试
21
        if (alpha >= 0 && beta >= 0 && gamma >= 0) {
22
            // 插值深度
23
            float depth = alpha * tri.vertices[0].z +
24
                         beta * tri.vertices[1].z +
25
                         gamma * tri.vertices[2].z;
26

27
            // 深度测试
28
            if (depth < min_depth) {
29
                min_depth = depth;
30

31
                // 插值颜色
32
                final_color = alpha * tri.colors[0] +
33
                             beta * tri.colors[1] +
34
                             gamma * tri.colors[2];
35
            }
36
        }
37
    }
38

39
    // 写入framebuffer
40
    int pixel_index = pixel_y * width + pixel_x;
41
    framebuffer[pixel_index * 3 + 0] = final_color.x;
42
    framebuffer[pixel_index * 3 + 1] = final_color.y;
43
    framebuffer[pixel_index * 3 + 2] = final_color.z;
44
    depth_buffer[pixel_index] = min_depth;
45
}

Early-Z优化技术#

Z-Prepass：在主渲染pass之前，先进行一次只写深度的pass：

1
// Z-Prepass阶段
2
void z_prepass(const std::vector<Triangle>& triangles) {
3
    for (const auto& tri : triangles) {
4
        rasterize_triangle_depth_only(tri);
5
    }
6
}
7

8
// 主渲染阶段
9
void main_render_pass(const std::vector<Triangle>& triangles) {
10
    for (const auto& tri : triangles) {
11
        rasterize_triangle_with_early_z(tri);
12
    }
13
}

Early-Z测试：在片段着色器执行前进行深度测试，可以大幅减少不必要的着色计算。

层次化Z-Buffer（Hi-Z）#

基本原理：构建深度缓冲区的mipmap层次结构，每个层级存储下一级的最小/最大深度值。

快速剔除：

1
bool hierarchical_z_test(int x, int y, float z, int level) {
2
    int tile_x = x >> level;
3
    int tile_y = y >> level;
4

5
    float max_z = hi_z_buffer[level][tile_y * (width >> level) + tile_x];
6

7
    return z <= max_z;  // 如果当前深度大于tile最大深度，可以剔除
8
}

8.2.3 Mipmap与纹理过滤#

Mipmap原理：预计算多级纹理，根据像素覆盖面积选择合适级别

级别计算：

1
float calculate_mipmap_level(const Vector2f& duv_dx, const Vector2f& duv_dy,
2
                            int texture_width, int texture_height) {
3
    float du_dx = duv_dx.x() * texture_width;
4
    float dv_dx = duv_dx.y() * texture_height;
5
    float du_dy = duv_dy.x() * texture_width;
6
    float dv_dy = duv_dy.y() * texture_height;
7

8
    float max_sqr = std::max(du_dx * du_dx + dv_dx * dv_dx,
9
                            du_dy * du_dy + dv_dy * dv_dy);
10

11
    return 0.5f * std::log2(max_sqr);
12
}

三线性过滤：

在两个相邻mipmap级别分别进行双线性过滤
在两个结果间进行线性插值

深度测试与隐藏面消除#

9.1 Z-Buffer算法#

9.1.1 算法原理#

基本思想：为每个像素维护一个深度值，只绘制最近的表面

算法步骤：

1
// 初始化
2
for each pixel (x, y):
3
    color_buffer[x][y] = background_color
4
    depth_buffer[x][y] = infinity
5

6
// 渲染
7
for each triangle:
8
    for each pixel (x, y) in triangle:
9
        z = interpolated_depth(x, y)
10
        if z < depth_buffer[x][y]:
11
            depth_buffer[x][y] = z
12
            color_buffer[x][y] = shaded_color(x, y)

9.1.2 深度值计算的数学分析#

线性深度与非线性深度#

线性深度（观察空间）：在观察空间中，深度值与距离成线性关系： $z_{linear} = \frac{z_{eye} - n}{f - n}$

其中 $z_{eye}$ 是观察空间中的深度， $n$ 和 $f$ 分别是近平面和远平面距离。

非线性深度（透视投影后）：经过透视投影和透视除法后，NDC空间中的深度为： $z_{ndc} = \frac{f + n}{f - n} + \frac{2fn}{(f - n) \cdot z_{eye}}$

深度精度分析#

精度分布函数：深度缓冲区精度定义为： $\text{precision}(z) = \frac{dz_{buffer}}{dz_{eye}} = \frac{2fn}{z_{eye}^2(f-n)}$

关键观察：

精度与 $z_{eye}^2$ 成反比
近平面附近精度最高
远平面附近精度急剧下降

Z-Fighting的数学原因#

浮点精度限制：对于24位深度缓冲区，最小可分辨深度差为： $\Delta z_{min} = \frac{1}{2^{24}} \approx 5.96 \times 10^{-8}$

临界距离计算：两个表面在距离 $z$ 处可分辨的最小间距为： $\Delta z_{world} = \frac{z^2(f-n)}{2fn} \cdot \Delta z_{min}$

深度精度优化策略#

1. 近远平面比值优化：最优比值关系： $\frac{f}{n} < 10^3$

2. 对数深度缓冲：使用对数分布改善精度： $z_{log} = \frac{\log(z_{eye}/n)}{\log(f/n)}$

3. 反向Z缓冲：将深度映射反转： $z_{reverse} = 1 - z_{ndc}$

利用浮点数在0附近精度更高的特性。

工程实现#

1
// 深度精度分析工具
2
class DepthPrecisionAnalyzer {
3
public:
4
    static float compute_precision(float z_eye, float near, float far) {
5
        return (2.0f * far * near) / (z_eye * z_eye * (far - near));
6
    }
7

8
    static float min_resolvable_distance(float z_eye, float near, float far, int depth_bits) {
9
        float depth_resolution = 1.0f / (1 << depth_bits);
10
        float precision = compute_precision(z_eye, near, far);
11
        return depth_resolution / precision;
12
    }
13

14
    static void analyze_depth_distribution(float near, float far, int samples = 100) {
15
        std::cout << "Depth Precision Analysis:\n";
16
        std::cout << "Near: " << near << ", Far: " << far << "\n";
17
        std::cout << "Ratio: " << far/near << "\n\n";
18

19
        for (int i = 0; i < samples; i++) {
20
            float t = static_cast<float>(i) / (samples - 1);
21
            float z_eye = near + t * (far - near);
22
            float precision = compute_precision(z_eye, near, far);
23
            float min_dist = min_resolvable_distance(z_eye, near, far, 24);
24

25
            std::cout << "z=" << z_eye << ", precision=" << precision
26
                      << ", min_dist=" << min_dist << "\n";
27
        }
28
    }
29
};
30

31
// 反向Z缓冲实现
32
Matrix4f create_reverse_z_projection(float fov, float aspect, float near, float far) {
33
    Matrix4f proj = Matrix4f::Zero();
34

35
    float tan_half_fov = std::tan(fov * 0.5f);
36

37
    proj(0, 0) = 1.0f / (aspect * tan_half_fov);
38
    proj(1, 1) = 1.0f / tan_half_fov;
39
    proj(2, 2) = near / (far - near);           // 反向映射
40
    proj(2, 3) = (far * near) / (far - near);
41
    proj(3, 2) = 1.0f;                         // 注意：这里是+1而不是-1
42

43
    return proj;
44
}

9.2 其他隐藏面消除算法#

9.2.1 画家算法#

原理：按深度从远到近绘制物体

优点：

简单易实现
支持透明度混合

缺点：

需要排序，复杂度高
无法处理循环遮挡

实现：

1
struct Triangle {
2
    Vector3f vertices[3];
3
    float avg_depth;
4
    Color color;
5
};
6

7
bool depth_compare(const Triangle& a, const Triangle& b) {
8
    return a.avg_depth > b.avg_depth;  // 远到近排序
9
}
10

11
void painters_algorithm(std::vector<Triangle>& triangles) {
12
    // 计算每个三角形的平均深度
13
    for (auto& tri : triangles) {
14
        tri.avg_depth = (tri.vertices[0].z() +
15
                        tri.vertices[1].z() +
16
                        tri.vertices[2].z()) / 3.0f;
17
    }
18

19
    // 排序并绘制
20
    std::sort(triangles.begin(), triangles.end(), depth_compare);
21

22
    for (const auto& tri : triangles) {
23
        rasterize_triangle(tri);
24
    }
25
}

9.2.2 BSP树算法#

原理：用二叉空间分割树预处理场景

构建过程：

选择分割平面
将多边形分为前、后两组
递归构建子树

遍历渲染：

1
void render_bsp_tree(BSPNode* node, const Vector3f& view_pos) {
2
    if (!node) return;
3

4
    float distance = node->plane.distance_to_point(view_pos);
5

6
    if (distance > 0) {
7
        // 观察者在平面前方
8
        render_bsp_tree(node->back, view_pos);   // 先画后面
9
        render_polygons(node->polygons);         // 再画平面上的多边形
10
        render_bsp_tree(node->front, view_pos);  // 最后画前面
11
    } else {
12
        // 观察者在平面后方
13
        render_bsp_tree(node->front, view_pos);
14
        render_polygons(node->polygons);
15
        render_bsp_tree(node->back, view_pos);
16
    }
17
}

光照模型与着色#

10.1 光照的物理基础理论#

10.1.1 辐射度量学的数学框架#

基本辐射量的定义#

辐射能量（Radiant Energy）： $Q \quad [\text{焦耳, J}]$ 表示电磁辐射携带的总能量。

辐射通量（Radiant Flux/Power）： $\Phi = \frac{dQ}{dt} \quad [\text{瓦特, W}]$ 表示单位时间内通过某个表面的辐射能量。

辐射强度（Radiant Intensity）： $I(\omega) = \frac{d\Phi}{d\omega} \quad [\text{W/sr}]$ 表示点光源在某个方向上单位立体角内的辐射通量。

辐照度（Irradiance）： $E(\mathbf{p}) = \frac{d\Phi}{dA} \quad [\text{W/m}^2]$ 表示单位面积接收到的辐射通量。

辐射度（Radiance）： $L(\mathbf{p}, \omega) = \frac{d^2\Phi}{d\omega \, dA \cos\theta} \quad [\text{W/(sr}\cdot\text{m}^2\text{)}]$

这是最重要的辐射量，其中 $\theta$ 是方向 $\omega$ 与表面法向量的夹角。

辐射度的重要性质#

1. 感知相关性：辐射度直接对应于人眼或相机传感器接收到的光强。

2. 传播不变性：在真空中，辐射度沿直线传播时保持不变： $L(\mathbf{p}_1, \omega_{12}) = L(\mathbf{p}_2, \omega_{21})$

3. 可加性：多个光源的辐射度可以直接相加： $L_{total} = \sum_i L_i$

10.1.2 双向反射分布函数（BRDF）理论#

BRDF的严格数学定义#

微分形式定义： $f_r(\mathbf{p}, \omega_i, \omega_r) = \frac{dL_r(\mathbf{p}, \omega_r)}{dE_i(\mathbf{p}, \omega_i)} = \frac{dL_r(\mathbf{p}, \omega_r)}{L_i(\mathbf{p}, \omega_i) \cos\theta_i \, d\omega_i}$

其中：

$L_r(\mathbf{p}, \omega_r)$ ：反射辐射度
$E_i(\mathbf{p}, \omega_i)$ ：入射辐照度
$L_i(\mathbf{p}, \omega_i)$ ：入射辐射度

BRDF的物理约束条件#

1. 非负性约束： $f_r(\mathbf{p}, \omega_i, \omega_r) \geq 0$

2. 能量守恒约束： $\int_{\Omega} f_r(\mathbf{p}, \omega_i, \omega_r) \cos\theta_r \, d\omega_r \leq 1$

这确保反射的能量不超过入射能量。

3. 亥姆霍兹互易性： $f_r(\mathbf{p}, \omega_i, \omega_r) = f_r(\mathbf{p}, \omega_r, \omega_i)$

这是基于光路可逆性的物理原理。

BRDF的几何解释#

立体角微分： $d\omega = \frac{dA \cos\theta}{r^2}$

其中 $dA$ 是投影面积， $r$ 是距离， $\theta$ 是与法向量的夹角。

半球积分： BRDF在上半球面上的积分给出了反射率： $\rho(\omega_i) = \int_{\Omega} f_r(\omega_i, \omega_r) \cos\theta_r \, d\omega_r$

10.2 经典光照模型的数学推导#

10.2.1 Lambert漫反射模型#

物理原理与假设#

理想漫反射表面的特征：

各向同性：反射在所有方向上均匀分布
完全漫反射：没有镜面反射分量
表面粗糙：微观结构导致光线随机散射

Lambert BRDF的数学推导#

基本假设：反射辐射度在所有方向上恒定 $L_r(\omega_r) = \text{常数}$

能量守恒约束： $\int_{\Omega} f_r \cos\theta_r \, d\omega_r = \rho_d$

其中 $\rho_d$ 是漫反射率（albedo）。

半球积分计算： $\int_{\Omega} \cos\theta_r \, d\omega_r = \int_0^{2\pi} \int_0^{\pi/2} \cos\theta_r \sin\theta_r \, d\theta_r \, d\phi = \pi$

Lambert BRDF推导：由于 $f_r$ 为常数，结合能量守恒： $f_r \cdot \pi = \rho_d$

因此： $f_r = \frac{\rho_d}{\pi}$

Lambert余弦定律#

反射辐射度计算： $L_r = \int_{\Omega} f_r L_i(\omega_i) \cos\theta_i \, d\omega_i$

对于平行光源： $L_r = f_r L_i \cos\theta_i = \frac{\rho_d}{\pi} L_i \cos\theta_i$

物理意义：

$\cos\theta_i$ 项体现了投影面积效应
入射角越大，有效照射面积越小

工程实现#

1
Vector3f lambert_diffuse(const Vector3f& light_dir, const Vector3f& normal,
2
                        const Vector3f& light_color, const Vector3f& albedo) {
3
    // 计算入射角余弦值
4
    float cos_theta = std::max(0.0f, normal.dot(light_dir));
5

6
    // Lambert BRDF: $\frac{albedo}{\pi}$
7
    // 最终颜色 = $BRDF \times 入射光 \times \cos(\theta) \times \pi$ (积分因子)
8
    // 简化为：$albedo \times light\_color \times \cos(\theta)$
9
    return albedo * light_color * cos_theta;
10
}
11

12
// 更严格的实现（包含$\pi$因子）
13
Vector3f lambert_brdf_strict(const Vector3f& light_dir, const Vector3f& normal,
14
                            const Vector3f& light_color, const Vector3f& albedo) {
15
    float cos_theta = std::max(0.0f, normal.dot(light_dir));
16

17
    // 严格的Lambert BRDF
18
    Vector3f brdf = albedo / M_PI;
19

20
    // 渲染方程的离散形式
21
    return brdf * light_color * cos_theta * M_PI;  // $\pi$来自立体角积分
22
}

10.2.2 Phong反射模型的数学推导#

Phong模型的三分量结构#

Phong模型将光照分解为三个独立的分量：

1. 环境光（Ambient）： $I_a = k_a \cdot I_{ambient}$ 模拟复杂的间接光照，是一个常数项。

2. 漫反射（Diffuse）： $I_d = k_d \cdot I_{light} \cdot (\mathbf{N} \cdot \mathbf{L})$ 使用Lambert余弦定律。

3. 镜面反射（Specular）： $I_s = k_s \cdot I_{light} \cdot (\mathbf{R} \cdot \mathbf{V})^n$ 模拟光滑表面的镜面高光。

镜面反射的数学推导#

反射向量的计算：根据反射定律，入射向量 $\mathbf{L}$ 关于法向量 $\mathbf{N}$ 的反射向量为： $\mathbf{R} = 2(\mathbf{N} \cdot \mathbf{L})\mathbf{N} - \mathbf{L}$

推导过程：设入射向量为 $\mathbf{L}$ ，将其分解为法向分量和切向分量： $\mathbf{L} = \mathbf{L}_{\parallel} + \mathbf{L}_{\perp}$

其中：

$\mathbf{L}_{\parallel} = (\mathbf{L} \cdot \mathbf{N})\mathbf{N}$ （法向分量）
$\mathbf{L}_{\perp} = \mathbf{L} - \mathbf{L}_{\parallel}$ （切向分量）

反射时，切向分量不变，法向分量反向： $\mathbf{R} = \mathbf{L}_{\perp} - \mathbf{L}_{\parallel} = \mathbf{L} - 2\mathbf{L}_{\parallel} = \mathbf{L} - 2(\mathbf{L} \cdot \mathbf{N})\mathbf{N}$

整理得： $\mathbf{R} = 2(\mathbf{N} \cdot \mathbf{L})\mathbf{N} - \mathbf{L}$

完整的Phong光照方程#

$I_{total} = k_a I_a + k_d I_l (\mathbf{N} \cdot \mathbf{L}) + k_s I_l (\mathbf{R} \cdot \mathbf{V})^n$

其中：

$k_a, k_d, k_s$ ：材质系数
$n$ ：光泽度指数（shininess）
$\mathbf{V}$ ：视线方向向量

工程实现#

1
Vector3f phong_lighting(const Vector3f& position, const Vector3f& normal,
2
                       const Vector3f& view_dir, const Vector3f& light_pos,
3
                       const Vector3f& light_color, const Material& material) {
4
    // 环境光分量
5
    Vector3f ambient = material.ambient * light_color;
6

7
    // 计算光线方向
8
    Vector3f light_dir = (light_pos - position).normalized();
9

10
    // 漫反射分量（Lambert）
11
    float diff = std::max(normal.dot(light_dir), 0.0f);
12
    Vector3f diffuse = material.diffuse * light_color * diff;
13

14
    // 计算反射向量
15
    Vector3f reflect_dir = 2.0f * normal.dot(light_dir) * normal - light_dir;
16
    reflect_dir.normalize();
17

18
    // 镜面反射分量
19
    float spec = std::pow(std::max(view_dir.dot(reflect_dir), 0.0f), material.shininess);
20
    Vector3f specular = material.specular * light_color * spec;
21

22
    return ambient + diffuse + specular;
23
}

10.2.3 Blinn-Phong模型的数学改进#

Blinn-Phong的理论基础#

改进动机： Phong模型在计算反射向量时存在效率和稳定性问题，Blinn-Phong通过引入半角向量来解决这些问题。

半角向量的数学定义： $\mathbf{H} = \frac{\mathbf{L} + \mathbf{V}}{\|\mathbf{L} + \mathbf{V}\|}$

其中 $\mathbf{L}$ 是光线方向， $\mathbf{V}$ 是视线方向。

几何意义与物理解释#

几何意义： $\mathbf{H}$ 是光线方向和视线方向的角平分线，表示能够产生镜面反射的理想微表面法向量。

物理解释：在微表面理论中，只有法向量与半角向量平行的微表面才会将光线从 $\mathbf{L}$ 方向反射到 $\mathbf{V}$ 方向。

Blinn-Phong镜面反射公式#

修改后的镜面反射项： $I_s = k_s I_l (\mathbf{N} \cdot \mathbf{H})^{n'}$

指数关系：为了获得与Phong模型相似的视觉效果，通常需要调整光泽度指数： $n' \approx 4n_{phong}$

数学优势分析#

1. 计算效率：

Phong： $\mathbf{R} = 2(\mathbf{N} \cdot \mathbf{L})\mathbf{N} - \mathbf{L}$ （6次乘法，3次减法）
Blinn-Phong： $\mathbf{H} = \frac{\mathbf{L} + \mathbf{V}}{\|\mathbf{L} + \mathbf{V}\|}$ （3次加法，1次归一化）

2. 数值稳定性：避免了反射向量在掠射角附近的数值不稳定问题。

3. 物理合理性：更好地符合微表面理论的物理基础。

代码实现：

1
Vector3f blinn_phong_specular(const Vector3f& light_dir, const Vector3f& view_dir,
2
                             const Vector3f& normal, const Vector3f& light_color,
3
                             const Vector3f& specular_color, float shininess) {
4
    Vector3f half_dir = (light_dir + view_dir).normalized();
5
    float spec = std::pow(std::max(normal.dot(half_dir), 0.0f), shininess);
6
    return specular_color * light_color * spec;
7
}

10.3 着色技术#

10.3.1 平面着色（Flat Shading）#

原理：每个三角形使用单一颜色

法向量计算：

1
Vector3f calculate_face_normal(const Vector3f& v1, const Vector3f& v2, const Vector3f& v3) {
2
    Vector3f edge1 = v2 - v1;
3
    Vector3f edge2 = v3 - v1;
4
    return edge1.cross(edge2).normalized();
5
}

特点：

计算简单
多面体外观明显
适合低多边形风格

10.3.2 Gouraud着色#

原理：在顶点计算光照，三角形内部插值

算法步骤：

计算顶点法向量（相邻面法向量平均）
在顶点进行光照计算
在三角形内部插值颜色

顶点法向量计算：

1
void calculate_vertex_normals(Mesh& mesh) {
2
    // 初始化顶点法向量为零
3
    for (auto& vertex : mesh.vertices) {
4
        vertex.normal = Vector3f::Zero();
5
    }
6

7
    // 累加相邻面的法向量
8
    for (const auto& face : mesh.faces) {
9
        Vector3f face_normal = calculate_face_normal(
10
            mesh.vertices[face.v1].position,
11
            mesh.vertices[face.v2].position,
12
            mesh.vertices[face.v3].position
13
        );
14

15
        mesh.vertices[face.v1].normal += face_normal;
16
        mesh.vertices[face.v2].normal += face_normal;
17
        mesh.vertices[face.v3].normal += face_normal;
18
    }
19

20
    // 归一化
21
    for (auto& vertex : mesh.vertices) {
22
        vertex.normal.normalize();
23
    }
24
}

10.3.3 Phong着色#

原理：插值法向量，在每个像素进行光照计算

算法步骤：

在顶点存储法向量
在三角形内部插值法向量
在每个像素进行光照计算

法向量插值：

1
Vector3f interpolate_normal(const Vector3f& bary,
2
                           const Vector3f& n1, const Vector3f& n2, const Vector3f& n3) {
3
    Vector3f interpolated = bary.x() * n1 + bary.y() * n2 + bary.z() * n3;
4
    return interpolated.normalized();
5
}

质量对比：

Flat < Gouraud < Phong（质量递增）
Flat > Gouraud > Phong（性能递减）

纹理映射技术#

11.1 纹理映射基础理论#

11.1.1 纹理坐标系统的数学基础#

UV坐标系的定义#

标准化纹理坐标：纹理坐标 $(u, v)$ 定义在单位正方形内： $(u, v) \in [0, 1] \times [0, 1]$

其中：

$u$ 轴：水平方向，对应纹理的宽度
$v$ 轴：垂直方向，对应纹理的高度

坐标变换#

纹理坐标到像素坐标的映射：

最近邻采样的像素坐标： $x_{pixel} = \lfloor u \cdot (W - 1) \rfloor$ $y_{pixel} = \lfloor v \cdot (H - 1) \rfloor$

其中 $W$ 和 $H$ 分别是纹理的宽度和高度（以像素为单位）。

连续坐标映射：对于需要插值的情况：

连续坐标计算： $x_{continuous} = u \cdot W - 0.5$ $y_{continuous} = v \cdot H - 0.5$

减去0.5是为了将采样点放在像素中心。

工程实现#

1
struct TextureCoordinate {
2
    float u, v;
3

4
    // 转换为像素坐标（整数）
5
    std::pair<int, int> to_pixel_coords(int width, int height) const {
6
        int x = static_cast<int>(u * (width - 1));
7
        int y = static_cast<int>(v * (height - 1));
8
        return {x, y};
9
    }
10

11
    // 转换为连续像素坐标（用于插值）
12
    std::pair<float, float> to_continuous_coords(int width, int height) const {
13
        float x = u * width - 0.5f;
14
        float y = v * height - 0.5f;
15
        return {x, y};
16
    }
17
};

11.1.2 纹理采样理论#

点采样（最近邻采样）#

数学定义：点采样选择距离采样点最近的纹理像素： $T_{nearest}(u, v) = T\left(\left\lfloor u \cdot W + 0.5 \right\rfloor, \left\lfloor v \cdot H + 0.5 \right\rfloor\right)$

其中 $T(i, j)$ 表示纹理在像素 $(i, j)$ 处的颜色值。

双线性过滤的数学推导#

问题设定：给定连续纹理坐标 $(u, v)$ ，计算对应的纹理值。

坐标变换：

连续坐标计算： $x = u \cdot W - 0.5$ $y = v \cdot H - 0.5$

四个邻近像素：

像素坐标定义： $(x_0, y_0) = (\lfloor x \rfloor, \lfloor y \rfloor)$ $(x_1, y_0) = (x_0 + 1, y_0)$ $(x_0, y_1) = (x_0, y_0 + 1)$ $(x_1, y_1) = (x_0 + 1, y_0 + 1)$

插值权重：

权重计算： $f_x = x - x_0$ $f_y = y - y_0$

双线性插值公式：

分步计算： $T_{bilinear}(u, v) = [1-f_y]((1-f_x)T_{00} + f_x T_{10}) + f_y[(1-f_x)T_{01} + f_x T_{11}]$

展开形式： $T_{bilinear}(u, v) = (1-f_x)(1-f_y)T_{00} + f_x(1-f_y)T_{10} + (1-f_x)f_y T_{01} + f_x f_y T_{11}$

其中 $T_{ij}$ 表示像素 $(x_i, y_j)$ 的颜色值。

几何解释：双线性插值等价于先在x方向进行两次线性插值，再在y方向进行一次线性插值：

分步插值过程：

x方向插值： $C_0 = (1-f_x)T_{00} + f_x T_{10}$ $C_1 = (1-f_x)T_{01} + f_x T_{11}$

y方向插值： $T_{bilinear} = (1-f_y)C_0 + f_y C_1$

工程实现#

1
class TextureSampler {
2
public:
3
    static Color sample_nearest(const Texture& texture, float u, float v) {
4
        int x = static_cast<int>(u * texture.width + 0.5f) % texture.width;
5
        int y = static_cast<int>(v * texture.height + 0.5f) % texture.height;
6
        return texture.get_pixel(x, y);
7
    }
8

9
    static Color sample_bilinear(const Texture& texture, float u, float v) {
10
        float x = u * texture.width - 0.5f;
11
        float y = v * texture.height - 0.5f;
12

13
        int x0 = static_cast<int>(std::floor(x));
14
        int y0 = static_cast<int>(std::floor(y));
15
        int x1 = x0 + 1;
16
        int y1 = y0 + 1;
17

18
        float fx = x - x0;
19
        float fy = y - y0;
20

21
        // 获取四个邻近像素（带边界处理）
22
        Color c00 = texture.get_pixel_wrap(x0, y0);
23
        Color c10 = texture.get_pixel_wrap(x1, y0);
24
        Color c01 = texture.get_pixel_wrap(x0, y1);
25
        Color c11 = texture.get_pixel_wrap(x1, y1);
26

27
        // 双线性插值
28
        Color c0 = lerp(c00, c10, fx);
29
        Color c1 = lerp(c01, c11, fx);
30
        return lerp(c0, c1, fy);
31
    }
32

33
private:
34
    static Color lerp(const Color& a, const Color& b, float t) {
35
        return a * (1.0f - t) + b * t;
36
    }
37
};

11.1.3 纹理寻址模式的数学定义#

重复模式（Repeat/Wrap）#

数学定义： $u_{repeat} = u - \lfloor u \rfloor = u \bmod 1$

性质：

周期性： $u_{repeat}(u + k) = u_{repeat}(u)$ 对任意整数 $k$
值域： $u_{repeat} \in [0, 1)$

镜像模式（Mirror/Reflect）#

数学定义：

镜像映射函数：

当 $t \leq 0.5$ 时： $u_{mirror} = 2t$

当 $t > 0.5$ 时： $u_{mirror} = 2(1-t)$

其中 $t = (u/2) \bmod 1$ 。

等价表示： $u_{mirror} = 1 - \left|2 \cdot \left(\frac{u}{2} - \left\lfloor \frac{u}{2} + 0.5 \right\rfloor\right)\right|$

边缘拉伸模式（Clamp to Edge）#

数学定义： $u_{clamp} = \max(0, \min(1, u))$

分段函数表示：

钳制映射函数：

当 $u < 0$ 时： $u_{clamp} = 0$

当 $0 \leq u \leq 1$ 时： $u_{clamp} = u$

当 $u > 1$ 时： $u_{clamp} = 1$

边界颜色模式（Border Color）#

数学定义：

边界纹理映射函数：

当 $(u, v) \in [0, 1]^2$ 时： $T_{border}(u, v) = T(u, v)$

其他情况： $T_{border}(u, v) = C_{border}$

其中 $C_{border}$ 是预定义的边界颜色。

工程实现#

1
enum class WrapMode {
2
    REPEAT,
3
    MIRROR,
4
    CLAMP,
5
    BORDER
6
};
7

8
class TextureAddressing {
9
public:
10
    static float apply_wrap_mode(float coord, WrapMode mode) {
11
        switch (mode) {
12
            case WrapMode::REPEAT:
13
                return coord - std::floor(coord);
14

15
            case WrapMode::MIRROR: {
16
                float t = (coord * 0.5f) - std::floor(coord * 0.5f);
17
                return (t <= 0.5f) ? (2.0f * t) : (2.0f * (1.0f - t));
18
            }
19

20
            case WrapMode::CLAMP:
21
                return std::max(0.0f, std::min(1.0f, coord));
22

23
            case WrapMode::BORDER:
24
                return coord;  // 边界检查在采样时处理
25

26
            default:
27
                return coord;
28
        }
29
    }
30

31
    static std::pair<float, float> apply_wrap_mode_2d(float u, float v,
32
                                                     WrapMode u_mode, WrapMode v_mode) {
33
        return {apply_wrap_mode(u, u_mode), apply_wrap_mode(v, v_mode)};
34
    }
35
};

11.2 高级纹理技术#

11.2.1 Mipmap技术详解#

Mipmap生成算法：

1
void generate_mipmaps(Texture& texture) {
2
    int width = texture.width;
3
    int height = texture.height;
4
    int level = 0;
5

6
    while (width > 1 || height > 1) {
7
        int new_width = std::max(1, width / 2);
8
        int new_height = std::max(1, height / 2);
9

10
        Texture mip_level(new_width, new_height);
11

12
        for (int y = 0; y < new_height; ++y) {
13
            for (int x = 0; x < new_width; ++x) {
14
                // 2x2像素区域平均
15
                Color sum(0, 0, 0, 0);
16
                int count = 0;
17

18
                for (int dy = 0; dy < 2; ++dy) {
19
                    for (int dx = 0; dx < 2; ++dx) {
20
                        int src_x = x * 2 + dx;
21
                        int src_y = y * 2 + dy;
22

23
                        if (src_x < width && src_y < height) {
24
                            sum += texture.get_pixel(src_x, src_y);
25
                            count++;
26
                        }
27
                    }
28
                }
29

30
                mip_level.set_pixel(x, y, sum / count);
31
            }
32
        }
33

34
        texture.mip_levels[++level] = mip_level;
35
        width = new_width;
36
        height = new_height;
37
    }
38
}

Mipmap级别计算：

1
float calculate_mipmap_level(const Vector2f& duv_dx, const Vector2f& duv_dy,
2
                            int texture_width, int texture_height) {
3
    float du_dx = duv_dx.x() * texture_width;
4
    float dv_dx = duv_dx.y() * texture_height;
5
    float du_dy = duv_dy.x() * texture_width;
6
    float dv_dy = duv_dy.y() * texture_height;
7

8
    float max_sqr = std::max(du_dx * du_dx + dv_dx * dv_dx,
9
                            du_dy * du_dy + dv_dy * dv_dy);
10

11
    return 0.5f * std::log2(max_sqr);
12
}

11.2.2 法线贴图（Normal Mapping）#

切线空间法向量：

1
Vector3f sample_normal_map(const Texture& normal_map, float u, float v) {
2
    Color normal_color = sample_bilinear(normal_map, u, v);
3

4
    // 从[0,1]映射到[-1,1]
5
    Vector3f normal;
6
    normal.x() = normal_color.r * 2.0f - 1.0f;
7
    normal.y() = normal_color.g * 2.0f - 1.0f;
8
    normal.z() = normal_color.b * 2.0f - 1.0f;
9

10
    return normal.normalized();
11
}

TBN矩阵构建：

1
Matrix3f build_tbn_matrix(const Vector3f& normal, const Vector3f& tangent) {
2
    Vector3f N = normal.normalized();
3
    Vector3f T = tangent.normalized();
4

5
    // Gram-Schmidt正交化
6
    T = (T - N.dot(T) * N).normalized();
7
    Vector3f B = N.cross(T);
8

9
    Matrix3f TBN;
10
    TBN.col(0) = T;  // Tangent
11
    TBN.col(1) = B;  // Bitangent
12
    TBN.col(2) = N;  // Normal
13

14
    return TBN;
15
}

切线空间到世界空间变换：

1
Vector3f transform_normal(const Vector3f& tangent_normal, const Matrix3f& TBN) {
2
    return (TBN * tangent_normal).normalized();
3
}

11.2.3 环境映射#

球面映射：

1
Vector2f sphere_mapping(const Vector3f& direction) {
2
    Vector3f d = direction.normalized();
3
    float u = 0.5f + std::atan2(d.z(), d.x()) / (2.0f * M_PI);
4
    float v = 0.5f - std::asin(d.y()) / M_PI;
5
    return Vector2f(u, v);
6
}

立方体映射：

1
struct CubemapSample {
2
    int face;
3
    Vector2f uv;
4
};
5

6
CubemapSample cube_mapping(const Vector3f& direction) {
7
    Vector3f abs_dir = direction.cwiseAbs();
8
    CubemapSample sample;
9

10
    if (abs_dir.x() >= abs_dir.y() && abs_dir.x() >= abs_dir.z()) {
11
        // X面
12
        sample.face = direction.x() > 0 ? 0 : 1;  // +X : -X
13
        sample.uv.x() = direction.x() > 0 ? -direction.z() : direction.z();
14
        sample.uv.y() = -direction.y();
15
        sample.uv /= abs_dir.x();
16
    } else if (abs_dir.y() >= abs_dir.z()) {
17
        // Y面
18
        sample.face = direction.y() > 0 ? 2 : 3;  // +Y : -Y
19
        sample.uv.x() = direction.x();
20
        sample.uv.y() = direction.y() > 0 ? direction.z() : -direction.z();
21
        sample.uv /= abs_dir.y();
22
    } else {
23
        // Z面
24
        sample.face = direction.z() > 0 ? 4 : 5;  // +Z : -Z
25
        sample.uv.x() = direction.z() > 0 ? direction.x() : -direction.x();
26
        sample.uv.y() = -direction.y();
27
        sample.uv /= abs_dir.z();
28
    }
29

30
    // 转换到[0,1]范围
31
    sample.uv = (sample.uv + Vector2f(1, 1)) * 0.5f;
32
    return sample;
33
}