作业7-光线追踪
1 总览
在之前的练习中，我们实现了 Whitted-Style Ray Tracing 算法，并且用 BVH 等加速结构对于求交过程进行了加速。在本次实验中，我们将在上一次实验的基础上实现完整的 Path Tracing 算法。至此，我们已经来到了光线追踪版块的最后一节内容。
2 调通框架
2.1 框架的变化
相比上一次实验，本次实验对框架的修改较大，主要在以下几方面：
- 修改了 main.cpp，以适应本次实验的测试模型 CornellBox
- 修改了 Render，以适应 CornellBox 并且支持 Path Tracing 需要的同一 Pixel 多次 Sample
- 修改了 Object，Sphere，Triangle，TriangleMesh，BVH，添加了 area 属性与Sample 方法，以实现对光源按面积采样，并在 Scene 中添加了采样光源的接口 sampleLight
- 修改了 Material 并在其中实现了 sample, eval, pdf 三个方法用于 Path Tracing 变量的辅助计算
2.2 需要迁移的内容
需要从上一次编程练习中直接拷贝以下函数到对应位置：
Triangle::getIntersectionin Triangle.hpp: 将你的光线-三角形相交函数粘贴到此处，请直接将上次实验中实现的内容粘贴在此。IntersectP(const Ray& ray, const Vector3f& invDir,const std::array<int, 3>& dirIsNeg)in the Bounds3.hpp: 这个函数的作用是判断包围盒 BoundingBox 与光线是否相交，请直接将上次实验中实现的内容粘贴在此处，并且注意检查 t_enter = t_exit 的时候的判断是否正确。getIntersection(BVHBuildNode* node, const Ray ray)in BVH.cpp: BVH查找过程，请直接将上次实验中实现的内容粘贴在此处.
2.3 编译运行
基础代码只依赖于 CMake，下载基础代码后，执行下列命令，就可以编译这个项目：
cmake -B build
cmake --build build
之后可以通过
./build/Raytracing
来执行程序。请务必确保程序可以正确编译之后，再进入下一节内容
3 开始实现
3.1 代码框架
在本次实验中只需要修改下面这一个函数:
// Implementation of Path Tracing
Vector3f Scene::castRay(const Ray &ray, int depth) const
castRay(const Ray ray, int depth)in Scene.cpp: 在其中实现 Path Tracing 算法
3.2 可能用到的函数
可能用到的函数有：
3.2.1 intersect()
Intersection Scene::intersect(const Ray &ray) const
intersect(const Ray ray)in Scene.cpp: 求一条光线Ray与场景的交点，这个函数实际上就是调用Scene类的成员变量bvh的Intersect()方法
Ray结构体包含起点、方向等成员变量，有()和<<两个符号重载，()用于把t代入参数方程进行计算，<<用于打印结构体变量的信息。
struct Ray{
//Destination = origin + t*direction
Vector3f origin;
Vector3f direction, direction_inv;
double t;//transportation time,
double t_min, t_max;
Vector3f operator()(double t) const;
friend std::ostream &operator<<(std::ostream& os, const Ray& r)
};
返回的类型是Intersection结构体，这个结构体的定义如下
struct Intersection
{
Intersection(){
happened=false;
coords=Vector3f();
normal=Vector3f();
distance= std::numeric_limits<double>::max();
obj =nullptr;
m=nullptr;
}
bool happened;
Vector3f coords;
Vector3f tcoords;
Vector3f normal;
Vector3f emit;
double distance;
Object* obj;
Material* m;
};
| 成员变量 | 变量类型 | 含义 |
|---|---|---|
bool happened |
布尔值 | 标记光线是否与物体发生相交（true = 相交，false = 未相交） |
Vector3f coords |
三维向量 | 相交点的空间坐标（x,y,z） |
Vector3f tcoords |
三维向量 | 相交点的纹理坐标（Texture Coords），用于材质贴图映射 |
Vector3f normal |
三维向量 | 相交点处物体表面的法向量（用于计算光照、反射 / 折射） |
Vector3f emit |
三维向量 | 相交点处物体的自发光颜色（如光源的发光颜色，RGB） |
double distance |
双精度浮点数 | 从光线起点到相交点的距离（用于筛选最近的相交物体） |
Object* obj |
指针 | 指向发生相交的那个物体（如球体、平面等） |
Material* m |
指针 | 指向相交物体的材质（如漫反射、镜面反射材质） |
| #### 3.2.2 sampleLight() | ||
|
||
sampleLight(Intersection &pos, float &pdf) in Scene.cpp: 在场景的所有光源上按面积 uniform 地 sample 一个点，并计算该 sample 的概率密度。 |
||
函数通过参数pos返回采样点，通过参数pdf返回对应的采样概率密度。 |
概率密度函数与概率密度
概率密度函数 (Probability Density Function, PDF) 是描述连续型随机变量在某个取值点附近 “可能性” 相对大小的函数。
设 \(X\) 是连续型随机变量，若存在非负可积函数 \(f(x)\)，使得对任意实数 \(a \le b\)，都有：\(P(a \le X \le b) = \int_{a}^{b} f(x)\,dx\)则称 \(f(x)\) 为 \(X\) 的概率密度函数。
\(f(x_0)\)被称为随机变量\(X\)在某一点\(x_0\)处的概率密度。可以通俗理解为
\(\mathrm{d}P\)表示随机变量\(X\)落在\(x\)附近\(\mathrm{d}x\)范围内的概率
因此说到概率密度，就要明确随机变量是什么。需要注意的是这个随机变量可以是多维的。
void Scene::sampleLight(Intersection &pos, float &pdf) const
{
float emit_area_sum = 0;
for (uint32_t k = 0; k < objects.size(); ++k) {
if (objects[k]->hasEmit()){
emit_area_sum += objects[k]->getArea();
}
}
float p = get_random_float() * emit_area_sum;
emit_area_sum = 0;
for (uint32_t k = 0; k < objects.size(); ++k) {
if (objects[k]->hasEmit()){
emit_area_sum += objects[k]->getArea();
if (p <= emit_area_sum){
objects[k]->Sample(pos, pdf);
break;
}
}
}
}
- 遍历场景中所有物体，筛选出带自发光属性的物体（
hasEmit()为 true），累加它们的面积得到emit_area_sum。 get_random_float()生成[0,1)的随机浮点数，乘以总发光面积，得到一个落在[0, 总发光面积)范围内的随机值p，用于后续 “命中” 具体的发光物体。- 再次遍历发光物体并累加面积：一旦累加面积大于p，说明随机命中了该发光物体，调用其自身的
Sample方法：在该物体表面随机采样一个点，将采样点的相交信息写入pos，并计算该采样的pdf，随后跳出循环。
以Sphere::Sample()为例
void Sample(Intersection &pos, float &pdf){
float theta = 2.0 * M_PI * get_random_float(), phi = M_PI * get_random_float();
Vector3f dir(std::cos(phi), std::sin(phi)*std::cos(theta), std::sin(phi)*std::sin(theta));
pos.coords = center + radius * dir;
pos.normal = dir;
pos.emit = m->getEmission();
pdf = 1.0f / area;
}
一方面，Sample()会在光源表面随机采样一个点，把这个点的信息通过
pos返回，另一方面Sample()会计算光源表面积的倒数即1/area，把结果通过pdf返回。
3.2.3 sample()
Vector3f Material::sample(const Vector3f &wi, const Vector3f &N)
sample(const Vector3f wi, const Vector3f N) in Material.hpp: 按照该材质的性质，给定入射方向与法向量，用某种分布采样一个出射方向Vector3f Material::sample(const Vector3f &wi, const Vector3f &N){
switch(m_type){
case DIFFUSE:
{
// uniform sample on the hemisphere
float x_1 = get_random_float(), x_2 = get_random_float();
float z = std::fabs(1.0f - 2.0f * x_1);
float r = std::sqrt(1.0f - z * z), phi = 2 * M_PI * x_2;
Vector3f localRay(r*std::cos(phi), r*std::sin(phi), z);
return toWorld(localRay, N);
break;
}
}
}
对于DIFFUSE材质，需要在半球面上按面积均匀分布地随机采样一个方向。上述代码先在 \([0,1]\) 范围内随机取样一个z值，在 \([0,2\pi]\) 随机采样一个角度 phi，由z和phi就能在半球面上确定一个方向。
这种方式得到的方向在球面上是按面积平均的,粗略证明如下：
(注意：这里为方便取\(\theta\)为仰角，而不是天顶角)
首先有公式
所以
由此说明均匀分布的 \(z\) 等价于均匀分布的 \(S\)
3.2.4 pdf()
float Material::pdf(const Vector3f &wi, const Vector3f &wo, const Vector3f &N)
pdf(const Vector3f wi, const Vector3f wo, const Vector3f N) in Material.hpp: 给定一对入射、出射方向与法向量，计算 sample 方法得到该出射方向的概率密度float Material::pdf(const Vector3f &wi, const Vector3f &wo, const Vector3f &N){
switch(m_type){
case DIFFUSE:
{
// uniform sample probability 1 / (2 * PI)
if (dotProduct(wo, N) > 0.0f)
return 0.5f / M_PI;
else
return 0.0f;
break;
}
}
}
由于DIFFUSE材质向由法线确定的半球面反射是均匀分布的，因此概率密度为定值。在前面我们已经提到“说到概率密度就要明确随机变量是什么”。在这里随机变量是方立体角。立体角的定义是面积除以半径平方。因此半球面的立体角为
半球上的概率密度为
另一种稍加繁琐的推导方法
以法线\(\boldsymbol{N}\)为极轴建立球坐标系，任意方向与 N 的夹角为极角\(\theta\)（上半球中\(\theta \in [0, \frac{\pi}{2}]\)），绕 N 的旋转角为方位角\(\phi\)（\(\phi \in [0, 2\pi]\)），则立体角微元为：
球面坐标系
球坐标系 - 维基百科，自由的百科全书
在计算机图形学、光学的球面坐标系中，有两个角度参数，分别是方位角\(\phi\)和天顶角\(\theta\)。
需要注意的是，物理学中常定义方位角\(\phi\)和天顶角\(\theta\)，而数学中常定义方位角\(\theta\)和天顶角\(\phi\)，为了将其看作平面极坐标系的拓展。
上半球内的 PDF 是常数，设该常数为\(C\)（即\(p(\boldsymbol{w_i}) = C\)），
其中
解得
3.2.5 eval()
Vector3f Material::eval(const Vector3f &wi, const Vector3f &wo, const Vector3f &N)
eval(const Vector3f wi, const Vector3f wo, const Vector3f N) in Material.hpp: 给定一对入射、出射方向与法向量，计算这种情况下的 f_r 值
eval是英文单词evaluate的缩写，音标 /ɪˈvæljueɪt/，本义是评估、求值、计算，这是编程中极其通用的缩写形式，而在这段 PBR 渲染的代码里，eval有贴合光照计算的专属含义。
Vector3f Material::eval(const Vector3f &wi, const Vector3f &wo, const Vector3f &N){
switch(m_type){
case DIFFUSE:
{
// calculate the contribution of diffuse model
float cosalpha = dotProduct(N, wo);
if (cosalpha > 0.0f) {
Vector3f diffuse = Kd / M_PI;
return diffuse;
}
else
return Vector3f(0.0f);
break;
}
}
}
代码中使用了“朗伯 漫反射模型”
朗伯漫反射（Lambertian Diffuse）是最基础、最贴合现实中哑光材质的光散射模型，由光学中的 朗伯余弦定律 推导而来，专门描述理想漫反射表面的光照表现 —— 这类表面会把入射的光能量均匀地向半球空间（表面前方的所有方向）散射，没有任何方向的偏向性。
\(K_d\) 是材质的漫反射系数，Vector3f类型（RGB），取值范围是[0,1]，比如红色哑光墙的\(K_d=(1,0,0)\)，白色纸的\(K_d=(1,1,1)\)，黑色布的\(K_d=(0,0,0)\)，直接决定材质的基础颜色
为什么BRDF的漫反射项要除以π？ - 知乎
下面是我参考这篇文章，对于为什么要除以\(\pi\)的理解。
一般的渲染方程
\(L_e\) 部分为自发光，积分部分为反射光；\(\omega_i\)为指向光来向的向量；\(p\) 可以认为是反射点表面的性质；
反照率(albedo)
反照率（albedo）是行星物理学中用来表示天体反射本领的物理量，定义为物体的辐射度（radiosity）与辐照度（irradiance）之比。射是出，照是入，出射量除以入射量，得到无量纲量。反照率的公式定义如下
其中\(\sum L_o\)表示出射量(o表示out)，\(\sum L_i\)表示入射量(i表示in)
绝对黑体（black body）的反照率是0。煤炭呈黑色，反照率 接近0，因为它吸收了投射到其表面上的几乎所有可见光。镜面将可见光几乎全部反射出去，其反照率接近1。
albedo 翻译成反照率，与 reflectance（反射率）是有区别的。反射率用来表示某一种波长的反射能量与入射能量之比；而反照率用来表示全波段的反射能量与入射能量之比。BRDF 的 R 是 reflectance，方程仅关注一种波长。
由于不同波长的光能用 RGB 三原色表示，所以 albedo 也用0到1区间的 vec3 向量表示。当只存在漫反射，不存在镜面反射的情况下，albedo 才等于 diffuse，即等于漫反射系数 \(K_d\)
在不考虑自发光，仅有漫反射的情况下，\(f_r(p,\omega_i,\omega_o)\) 为常数，记为
所以有
证明\(\int _{\Omega^{+}} \cos\theta \ d\omega=\pi\)
首先有
所以
$$
\begin{align}
\int_{\Omega^+}\cos \theta \mathrm{d}\omega &= \int_{\Omega^+}\cos \theta\sin \theta\, \mathrm{d}\theta\, \mathrm{d}\phi \
&= \int_{0}^{2\pi}\mathrm{d}\phi \int_{0}^{\frac{\pi}{2}}\sin\theta\cos\theta\,\mathrm{d}\theta \
&= 2\pi \int_{0}^{\frac{\pi}{2}}\sin\theta\,\mathrm{d}(\sin\theta)\
&= 2\pi\cdot \left. \frac{\sin^2\theta}{2}\right|_{0}^{\frac{\pi}{2}}\
&= 2\pi\cdot \frac{1}{2} =\pi
\end{align}
$$
证毕
因为
所以
注意
\(L_i(p, \omega_i)\) 应该是关于 \(\theta\) 和 \(\phi\) 的函数吧，这里把它提出积分其实不太正确，但是原文章就是这样算的，我会在最后进行分析
其中
关于\(\sum L_i\)和\(\sum L_o\)的疑问
有人可能会问：
1. 我们之前在计算面元接收到来自\(\omega_i\) 方向的radiance时，需要乘以\(\cos\theta\)。这里的\(\sum L_i\)表示的不就是面元接收到的光吗，为什么这里不需要乘以\(\cos\theta\)？
- 还有就是\(\sum L_o\)表示的是面元向\(\omega_o\)方向发出的光(radiance)，根据光路可逆，入射和出射应该是等价的呀，为什么这里也不用乘以\(\cos\theta\)？
首先回答第一个问题。 \(\sum L_i\)表示的不是面元接收到的光，而是周围环境照射给面元的光。
因为我们把中间的微元看作面元，满足朗伯余弦定理，因此二者并不相等，从“环境照射的光”到“面元接受的光”需要乘以一个余弦值。如果中间的微元是当作微小球体的粒子，二者就能相等了。
再来回答第二个问题。 虽然光路是可逆的，但是radiance这个概念的定义其实是单向的。
Radiance 是对光线传播中的度量，是每单位立体角单位面积上的功率；
\(\mathrm{d}A\cos\theta\) 为单位面积在传播方向上的有效面积
在光的传播中我们可以通俗地分为“发射方”和“接收方”。这里的单位立体角和单位面积都指的都是接收方的立体角和面积，与发射方无关。因此我称radiance这个概念是单向的。这就导致计算接收光的时候要加 \(\cos\theta\)，而计算出射光的时候不加 \(\cos\theta\)。
因为反照率(albebo)公式
在只存在漫反射，不存在镜面反射的情况下，albedo 等于 diffuse，即等于漫反射系数\(K_d\) ,所以
由此，我们计算得到了朗伯漫反射的 \(f_r\)
现在再来回顾前文提到的问题。首先
由于 \(L_i(p, \omega_i)\) 不能提出来，等式左右乘以 \(\mathrm{d}\omega_o\) 并关于 \(\omega_o\) 积分，之后做除法，得
要让等号成立，我们需要让
这个公式并不恒成立，比如 \(L_i(p, \omega_i)=\cos\theta\) 时，则积分变为：
对 \(\phi\) 积分：\(\int_{0}^{2\pi}d\phi = 2\pi\)
对 \(\theta\) 积分（同之前计算）：\(\int_{0}^{\frac{\pi}{2}} \cos\theta(1-2\cos\theta)\sin\theta d\theta = -\frac{1}{6}\)
最终结果：\(2\pi \times (-\frac{1}{6}) = -\frac{\pi}{3} \neq 0\)，直接证明公式不恒成立。
3.3 可能用到的变量
可能用到的变量有：
RussianRoulette in Scene.cpp: 也即P_RR, 表示Russian Roulette 的概率
3.4 实现RayTracing
课程中介绍的 Path Tracing 伪代码如下 (为了与之前框架保持一致，wo 定义与课程介绍相反):
shade(p, wo)
Uniformly sample the light at xx (pdf_light=1/A)
shoot a ray from p to x
if the ray is not blocked in the middle
L_dir = L_i * f_r * cos_theta * cos_theta_x / |x-p|^2
L_indir = 0.0
Test Russian Rouletee with probability P_RR
Uniformly sample the hemisphere toward wi (pdf_hemi=1/2pi)
Trace a ray r(p, wi)
If ray r hit a non-emitting object at q
L_indir = shade(q, wi) * f_r * cos_theta / pdf_hemi / P_RR
Return L_dir + L_indir
按照本次实验给出的框架，我们进一步可以将伪代码改写为：
shade(p, wo)
sampleLight(inter, pdf_light)
Get x, ws, NN, emit form inter
Shoot a ray from p to x
If the ray is not blocked in the middle
L_dir = emit * eval(wo, ws, N) * dot(ws, N) * dot(ws, NN) / |x-p|^2 / pdf_light
L_indir = 0.0
Test Russian Roulette with probability RussianRoulette
wi = sample(wo, N)
Trace a ray r(p, wi)
If ray r hit a non-emitting object at q
L_indir = shade(q, wi) * eval(wo, wi, N) * dot(wi, N) / pdf(wo, wi, N) / RussianRoulette
Return L_dir + L_indir
为什么这个伪代码是可行的（指光源和非光源分别采样）
4 遇到的问题
GAMES101/Assignment7/code/Bounds3.hpp at master · Foggy-whale/GAMES101 · GitHub
4.1 无法正确处理相交
我实现了castRay()之后,运行程序发现结果总是完全黑色。最开始我以为是castRay()对于光源相关的计算和判断出了问题。
但是经过调试发现，intersect()函数计算得到的相交物体hasEmit()始终返回的是false，于是我直接调试查看由摄像机发出的光线能否射到光源上，最终发现程序没有检测到摄像机发出的光线射到光源，这与实际不符。
为了更方便进行调试，我修改castRay(),让它生成深度图。
Vector3f Scene::castRay(const Ray& ray, int depth ) const {
// 场景深度范围配置（根据你的场景实际尺寸调整）
const float NEAR_DEPTH = 0.0f; // 最近深度（射线起点）
const float FAR_DEPTH = 1500.0f; // 最远深度（超出此范围视为最远）
const float EPSILON = 1e-4f; // 射线偏移量，解决自相交
const float GAMMA = 2.2f; // 伽马校正系数（固定2.2适配人眼）
// 1. 检测射线与场景的交点
Intersection inter = this->intersect(ray);
// 2. 未命中/超出最远深度 → 返回纯黑（最大深度值0.0）
if (!inter.happened) {
return Vector3f(1.0f,0.0f,0.0f);
}
// 3. 计算原始深度：射线起点到交点的欧氏距离
float raw_depth = (inter.coords - ray.origin).norm();
// 4. 步骤1：钳制深度到场景有效范围（避免异常值）
float clamped_depth = std::clamp(raw_depth, NEAR_DEPTH, FAR_DEPTH);
// 5. 步骤2：核心非线性归一化（近=1、远=0，符合常规深度图）
float normalized_depth = clamped_depth / FAR_DEPTH;
// 6. 步骤3：伽马校正（拉伸近处细节，优化视觉效果）
// 公式：1 - (归一化深度)^(1/γ) → 近处细节差异更明显
float gamma_corrected_depth = 1.0f - std::pow(normalized_depth, 1.0f / GAMMA);
// 7. 鲁棒性处理：自相交（深度过近）→ 强制纯白（最近深度值1.0）
if (raw_depth < EPSILON) {
gamma_corrected_depth = 1.0f;
}
// 8. 返回RGB三通道相同值（直接输出为深度图）
return Vector3f(gamma_corrected_depth);
}
深度图使用伽马矫正来拉伸近处细节，优化视觉效果。
先深度钳制，再非线性归一化，最后伽马矫正。无法射中物体时则直接着色为红色。
从结果可以发现很多面消失了，比如正方体的上表面，floor和right。
最开始我以为可能是面朝向的问题，我把模型导入到blender里面，打开视图的面朝向选项，此时面的正面会被着色为蓝色，背面会被着色为红色。由下图中blender显示的结果可以发现，面朝向都是正确的。因此只可能是相交判断出了问题。
将我的代码和网上的实现进行比较，最终发现我Bounds3.hpp中的函数Bounds3::IntersectP没有正确处理边界情况。这个函数负责判断光线是否与包围盒相交，我原本写的判断条件为
return t_enter < t_exit && t_exit >= 0;
这忽视了包围盒为一个平面的情况，这种情况下
t_enter==t_exit，我们也应该认为光线和包围盒相交了。于是，我把这行代码修改为
return t_enter <= t_exit && t_exit >= 0;
这样，代码的逻辑才是正确的。最终能够得到正确的深度图。
5 我的实现
核心的两个公式如下
L_dir = emit * m->eval(wo, ws, N) * dotProduct(ws, N) * dotProduct(-ws, NN)
/ (x_p_distance * x_p_distance) / pdf_light;
L_indir = castRay(r, depth + 1) * m->eval(wo, wi, N) * dotProduct(wi, N)
/ m->pdf(wo, wi, N) / RussianRoulette;
- wo表示从着色点到观察者的矢量，即出射光
- ws表示从着色点到光源上的采样点的矢量
- N表示着手色处的法线矢量
- NN表示光源上采样点的法线矢量
- wi表示射向着色点的矢量，即入射光
需要特别注意第一个公式中dotProduct(-ws, NN)ws的符号
最开始渲染的结果下图所示
之后添加一个if判断，处理在depth=0直接打到光源的情况
else if (depth==0 && inter.obj->hasEmit() == true) {
return inter.m->getEmission();
}
这样就能正确渲染光源
完整代码如下：
const int MAX_RECURSION_DEPTH = 5; // 最大递归深度，避免栈溢出
// Implementation of Path Tracing;
Vector3f Scene::castRay(const Ray& ray, int depth) const {
// TO DO Implement Path Tracing Algorithm here
if (depth > MAX_RECURSION_DEPTH) {
return Vector3f(0.0f);
}
// get the intersection of the ray to scene
Intersection inter = this->intersect(ray);
if (inter.happened == false) {
return Vector3f(0.0);
}
if (depth>0 && inter.obj->hasEmit() == true) {
return Vector3f(0.0);
} else if (depth==0 && inter.obj->hasEmit() == true) {
return inter.m->getEmission();
}
Vector3f p = inter.coords;
Vector3f wo = -ray.direction;
Vector3f N = inter.normal;
Material *m = inter.m;
Intersection light_sample;
float pdf_light;
sampleLight(light_sample, pdf_light);
// Get x, ws, NN, emit from light_sample
Vector3f x = light_sample.coords;
Vector3f ws = (x - p).normalized();
const float EPSILON = 1e-4f;
Ray inter2light_ray = Ray(p + ws * EPSILON, ws);
Vector3f NN = light_sample.normal;
Vector3f emit = light_sample.emit;
Intersection light_inter;
light_inter = this->intersect(inter2light_ray);
Vector3f L_dir = Vector3f(0.0);
const double x_p_distance = (x - p).norm();
// if(light_inter.happened && light_inter.obj->hasEmit()) {
// printf("light_inter.happened==true\n");
// }
if (std::fabs(light_inter.distance - x_p_distance) < 0.001) {
L_dir = emit * m->eval(wo, ws, N) * dotProduct(ws, N) * dotProduct(-ws, NN)
/ (x_p_distance * x_p_distance) / pdf_light;
}
Vector3f L_indir = Vector3f(0.0);
if (get_random_float() < RussianRoulette) {
Vector3f wi = m->sample(wo, N);
Ray r = Ray(p, wi);
L_indir = castRay(r, depth + 1) * m->eval(wo, wi, N) * dotProduct(wi, N)
/ m->pdf(wo, wi, N) / RussianRoulette;
}
return L_dir + L_indir;
}
6 多线程加速
6.1 openmp环境配置
我使用的是MSYS2 UCRT64工具链。MSYS2提供了多种工具链比如Clang32，Clang64，MinGW64等，我使用的是UCRT64工具链环境。
首先更换MSYS2 UCRT64的pacman下载源为国内源，参考下面的链接
msys2更换国内源（多个文件（不是3个文件的版本！））msys2国内源-CSDN博客
打开MSYS2软件内的\etc\pacman.d\ 的 mirrorlist.ucrt64，其他环境的下载源也推荐更换，对应 mirrorlist.mingw64 等文件
在 mirrorlist.ucrt64 的开头添加
## Primary
Server = https://mirrors.tuna.tsinghua.edu.cn/msys2/mingw/ucrt64/
Server = http://mirrors.ustc.edu.cn/msys2/mingw/ucrt64/
Server = https://mirrors.ustc.edu.cn/msys2/mingw/ucrt64/
Server = http://mirrors.aliyun.com/msys2/mingw/ucrt64/
Server = https://mirrors.aliyun.com/msys2/mingw/ucrt64/
打开MSYS2 UCRT64终端，执行
pacman -Syu
来更新UCRT64的环境
之后，在UCRT64终端中参考下面的链接继续配置
pkg-helper/VSCode环境搭建基于C++的OpenMP+MPI的并行程序的Windows运行环境.md at main · LckOrLck/pkg-helper · GitHub
通过以下命令安装 OpenMP 支持：
pacman -S mingw-w64-ucrt-x86_64-llvm-openmp
通过以下命令安装 MS MPI（虽然在本实验中不需要用到MPI，但还是建议安装一下，可以锻炼自己配环境的能力）：
pacman -S mingw-w64-ucrt-x86_64-msmpi
将MSYS2的安装路径下的ucrt64文件夹中的bin、lib文件夹所在路径添加到本机环境变量中。
下载并安装 MS MPI（Windows 环境）
打开 MS MPI 下载页面, 下载并安装以下两个文件：
msmpisetup.exe：主安装程序，用于安装 MPI 运行时环境。msmpisdk.msi：MPI 开发工具包，包含头文件、库文件和示例程序。（可选）
安装好Microsoft MPI之后要在cmd中运行以配置环境变量
set msmpi
这是尝试在cmd中执行mpiexec指令，如果终端输出mpiexec的help信息（解释如何运行mpiexec），说明正常。否则如果终端无法识别mpiexec指令，需要把MS MPI的Bin路径添加到系统环境变量的Path中。如果MS MPI使用的是默认安装路径，则添加下述路径即可
C:\Program Files\Microsoft MPI\Bin
之后，应该终端应该就能识别mpiexec指令了。这部分可以参考 Windows 配置 MPI 环境 - 知乎 以及文章的评论区。
使用 文章 中给的测试程序hello.cpp
// 这个程序使用 MPI 来获取进程的编号，并使用 OpenMP 来打印每个线程的信息。
#include <iostream>
#include <omp.h>
#include <mpi.h>
int main(int argc, char* argv[]) {
// 初始化 MPI
MPI_Init(&argc, &argv);
int world_size, world_rank;
MPI_Comm_size(MPI_COMM_WORLD, &world_size);
MPI_Comm_rank(MPI_COMM_WORLD, &world_rank);
omp_set_num_threads(4);
// 使用 OpenMP 打印每个进程的信息
#pragma omp parallel
{
int thread_id = omp_get_thread_num();
int num_threads = omp_get_num_threads();
std::cout << "Hello from rank " << world_rank
<< " of " << world_size
<< " with thread " << thread_id
<< " of " << num_threads << " threads" << std::endl;
}
// 结束 MPI
MPI_Finalize();
return 0;
}
注意编译选项的设置
{
"tasks": [
{
"type": "cppbuild",
"label": "C/C++: g++.exe 生成活动文件",
"command": "C:/msys64/ucrt64/bin/g++.exe",
"args": [
"-fdiagnostics-color=always",
"-g",
"${file}",
"-fopenmp",
"-lmsmpi",
"-o",
"${fileDirname}\\${fileBasenameNoExtension}.exe",
""
],
"options": {
"cwd": "C:/msys64/ucrt64/bin"
},
"problemMatcher": [
"$gcc"
],
"group": {
"kind": "build",
"isDefault": true
},
"detail": "调试器生成的任务。"
}
],
"version": "2.0.0"
}
在powershell终端执行
mpiexec -n 4 hello.exe
在命令行设置openmp线程数
set OMP_NUM_THREADS=3
在编译选项中同时加入-openmp和-lomp会导致运行程序时得不到预期结果，会多次输出
Hello from rank 0 of 1 with thread 0 of 1 threads
因此只能加上
-openmp
6.2 修改源码并编译
在 CMakeLists.txt 中 add_exectable 前面添加下述语句
# Enable OpenMP if available
find_package(OpenMP)
if(OpenMP_CXX_FOUND)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
endif()
之后修改函数
Render() in Renderer.cpp，void Renderer::Render(const Scene& scene)
{
std::vector<Vector3f> framebuffer(scene.width * scene.height);
float scale = tan(deg2rad(scene.fov * 0.5));
float imageAspectRatio = scene.width / (float)scene.height;
Vector3f eye_pos(278, 273, -800);
// change the spp value to change sample ammount
int spp = 16;
std::cout << "SPP: " << spp << "\n";
for (uint32_t j = 0; j < scene.height; ++j) {
// parallelize inner loop over columns; keep rows sequential to update progress safely
#pragma omp parallel for schedule(dynamic)
for (uint32_t i = 0; i < scene.width; ++i) {
// generate primary ray direction
float x = (2 * (i + 0.5f) / (float)scene.width - 1) * imageAspectRatio * scale;
float y = (1 - 2 * (j + 0.5f) / (float)scene.height) * scale;
Vector3f dir = normalize(Vector3f(-x, y, 1));
int idx = j * scene.width + i;
Vector3f pixel_color(0.0f);
for (int k = 0; k < spp; k++){
pixel_color += scene.castRay(Ray(eye_pos, dir), 0) / (float)spp;
}
framebuffer[idx] = pixel_color;
}
UpdateProgress(j / (float)scene.height);
}
UpdateProgress(1.f);
// save framebuffer to file
FILE* fp = fopen("binary.ppm", "wb");
(void)fprintf(fp, "P6\n%d %d\n255\n", scene.width, scene.height);
for (auto i = 0; i < scene.height * scene.width; ++i) {
static unsigned char color[3];
color[0] = (unsigned char)(255 * std::pow(clamp(0, 1, framebuffer[i].x), 0.6f));
color[1] = (unsigned char)(255 * std::pow(clamp(0, 1, framebuffer[i].y), 0.6f));
color[2] = (unsigned char)(255 * std::pow(clamp(0, 1, framebuffer[i].z), 0.6f));
fwrite(color, 1, 3, fp);
}
fclose(fp);
}
主要是在变量i的for循环前面添加预编译指令
#pragma omp parallel for schedule(dynamic)
告诉编译器处理这个for循环时使用动态分配的多线程
由于openmp是安装在UCRT64环境中的，如果使用windows安装的cmake编译项目会导致找不到openmp，因此我在UCRT64环境中安装cmake。由于UCRT64会自动配置openmp的环境变量，因此这样在UCRT64终端运行cmake时就能正确找到openmp。
打开MSYS2 UCRT64终端，运行下述指令
# 安装 UCRT64 版本的 CMake
pacman -S --needed mingw-w64-ucrt-x86_64-cmake
该命令会自动将 CMake 安装到
/ucrt64/bin/ 目录（这个目录默认在 UCRT64 终端的 PATH 中）。
参考链接 Visual Studio Code 中的终端配置文件_Vscode中文网 新建一个UCRT64的终端配置
点击配置终端设置，在配置中随便找一个点击“在settings.json设置”进入settings.json，
在"terminal.integrated.profiles.windows"中添加"MSYS2_UCRT64"(自定义的名称)
{
"terminal.integrated.profiles.windows": {
"MSYS2_UCRT64": {
"path": "C:\\msys64\\msys2_shell.cmd",
"args": [
"-defterm",
"-here",
"-no-start",
"-ucrt64"
]
}
}
}
之后在vscode中新建一个UCRT64的终端，在里面执行cmake相关指令，就能正确编译。
6.3 加速效果
在联想拯救者y7000p 2024上运行
处理器 Intel(R) Core(TM) i7-14650HX，2200 Mhz，16 个内核，24 个逻辑处理器
已安装的物理内存(RAM) 16.0 GB
默认参数下，不使用多线程需要二十多分钟，使用多线程后cpu占用到达90%以上，运行时间缩短到2分钟以内。
评论
如果你已登录 GitHub，就可以直接在这里评论。