UCSD CSE272-Advanced Image Synthesis-Homework

Posted on Edited on
A record of doing the lab for CSE-272.

关于为什么要在这里开个新坑。。。因为我发现这门课作业并不是很多。UCSD的这门神课冷门得令人吃惊，但是它涵盖的内容也极其之广泛，可以说涵盖了离线渲染的相当大一部分内容（虽然有的可能讲得不深），值得一学。

另一方面，我自己的渲染器最近很长时间debug无果，估计是遇到瓶颈了，来学习学习Tzumao Li大佬的代码，把这门课为数不多的作业做了，实现一些新的算法然后继续完善自己的渲染器应该还是OK的。

update:\text{update:}实际上我发现作业量也不算小。。。虽然只有两次编程作业但每次任务量都很大。不过指南写得超级详细贴心所以做起来也不慢。

Hw-0: Lajolla Renderer\text{Hw-0: Lajolla Renderer}

起步作业是安装此课程指定的渲染器Lajolla\text{Lajolla}。照着任务书上的命令走就行，不过可能需要安装一下HDRView\text{HDRView}来查看渲染结束得到的exr\text{exr}文件，这些都是相当简单的事情。HDRView\text{HDRView}一次make可能会失败，多试几次也就成功了。

这次任务书还简要回顾了一些基于物理渲染的基础入门知识：渲染方程，重要性采样，多重重要性采样和路径空间积分等，PBR书上都有，就算没写过也不难理解要是这都理解不了就别学了。

简要说一下对于Lajolla\text{Lajolla}代码的初印象：不同于GAMES101\text{GAMES101}最后的光追框架和pbrt\text{pbrt}的框架，Lajolla\text{Lajolla}很多地方是通过std::variant\text{std::variant}来实现多态的，然后相应的运算封装成一个算子结构体（其实不是数学上的那个算子啦），重载括号运算写成仿函数来完成。换做是我如果用纯粹OOP\text{OOP}的想法（也是我自己渲染器的写法）可能会写一个虚基类然后搞一堆派生类出来实现一堆接口。。。不过实话说我也不明白这两种写法有何利弊。但是初来乍到搞不清楚std::variant\text{std::variant}的话还是容易不明白Lajolla\text{Lajolla}在干什么。

Hw-1: Disney Principled BSDF\text{Hw-1: Disney Principled BSDF}

实话说，这更像是一些hack\text{hack}的感觉而不是真正的基于物理的有坚实科学依据和数学依据的东西。。。但是这些东西写了确实能丰富渲染器的材质，而且也是艺术家很常用的，所以写一写也无妨。作业里的思考题也可以想一想。但是这一部分作业不会全做，主要做一些我认为能学到新东西的部分以及思考题。

1. Diffuse\text{1. Diffuse}

1.1. Coding\text{1.1. Coding}

按照任务书上给出的公式开始写代码。实现两种BSDF\text{BSDF}分别描述普通的漫反射和次表面散射然后线性插值混合。这里就要注意Lajolla\text{Lajolla}代码框架的特点了。

首先实现Disney principled diffuse\text{Disney principled diffuse}材质。照着公式写倒是没啥难的，但是注意传进仿函数eval_op\text{eval\_op}的bsdf\text{bsdf}——可能上来就会把这玩意儿的base_color\text{base\_color}成员和M_1_PI\text{M\_1\_PI}乘在一起然后喜提报错。这里稍微解释一下，也提醒自己注意。

bsdf\text{bsdf}是DisneyDiffuse\text{DisneyDiffuse}类的实例，这个类的定义如下：

注意这里的成员都是Texture<Spectrum>\text{Texture<Spectrum>}类的——不是Spectrum\text{Spectrum}类的。稍微翻一翻代码就能发现：Texture\text{Texture}是一个std::variant\text{std::variant}的别名，后者中包括了三种纹理，从名字上看分别是纯色纹理，图片纹理和方格纹理。这三种纹理类也都是提供了一致的接口。

然后真正在访问这些纹理的信息的时候，需要使用material.h\text{material.h}中提供的函数eval\text{eval}来获取。后者需要传入一堆参数，只需要将bsdf\text{bsdf}中的成员和一些eval_op\text{eval\_op}类本身的vertex\text{vertex}成员的一些成员传进去就可以。

对于重要性采样，我们学习抄袭lambertian.inl\text{lambertian.inl}的余弦重要性采样就差不多，文档里也写了你可以随便复制粘贴，在sample_bsdf_op\text{sample\_bsdf\_op}里最后返回函数传参里的roughness\text{roughness}传入bsdf\text{bsdf}算出来的roughness\text{roughness}就行。另外一个算pdf\text{pdf}的甚至改都不用改。

上代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
#include "../vector.h"
Spectrum eval_op::operator()(const DisneyDiffuse &bsdf) const {
    if (dot(vertex.geometry_normal, dir_in) < 0 ||
            dot(vertex.geometry_normal, dir_out) < 0) {
        // No light below the surface
        return make_zero_spectrum();
    }
    // Flip the shading frame if it is inconsistent with the geometry normal
    Frame frame = vertex.shading_frame;
    if (dot(frame.n, dir_in) < 0) {
        frame = -frame;
    }

    // Homework 1: implement this!
    const TVector3 H = normalize(dir_in + dir_out);
    Real tmp = dot(H, dir_out);
    Real Roughness = eval(bsdf.roughness, vertex.uv, vertex.uv_screen_size, texture_pool);
    Real NdotWin = abs(dot(vertex.geometry_normal, dir_in));
    Real NdotWout = abs(dot(vertex.geometry_normal, dir_out));
    Real FD90 = 0.5 + Real(2) * Roughness * tmp * tmp;
    Real FD_wi = 1.0 + (FD90 - 1.0) * pow(1.0 - NdotWin, 5);
    Real FD_wo = 1.0 + (FD90 - 1.0) * pow(1.0 - NdotWout, 5);
    Spectrum BaseColor = eval(bsdf.base_color, vertex.uv, vertex.uv_screen_size, texture_pool);
    Spectrum f_baseDiffuse = BaseColor * M_1_PI * FD_wi * FD_wo * NdotWout;
    Real FSS90 = Roughness * tmp * tmp;
    Real FSS_wi = 1.0 + (FSS90 - 1.0) * pow(1 - NdotWin, 5);
    Real FSS_wo = 1.0 + (FSS90 - 1.0) * pow(1 - NdotWout, 5);
    Spectrum f_subsurface = 1.25 * BaseColor * M_1_PI * (FSS_wi * FSS_wo * (1.0 / (NdotWin + NdotWout) - 0.5) + 0.5) * NdotWout;
    Real Subsurface = eval(bsdf.subsurface, vertex.uv, vertex.uv_screen_size, texture_pool);
    return (1.0 - Subsurface) * f_baseDiffuse + Subsurface * f_subsurface;
}

Real pdf_sample_bsdf_op::operator()(const DisneyDiffuse &bsdf) const {
    if (dot(vertex.geometry_normal, dir_in) < 0 ||
            dot(vertex.geometry_normal, dir_out) < 0) {
        // No light below the surface
        return 0;
    }
    // Flip the shading frame if it is inconsistent with the geometry normal
    Frame frame = vertex.shading_frame;
    if (dot(frame.n, dir_in) < 0) {
        frame = -frame;
    }
    
    // Homework 1: implement this!
    return fmax(dot(frame.n, dir_out), Real(0)) / c_PI;
}

std::optional<BSDFSampleRecord> sample_bsdf_op::operator()(const DisneyDiffuse &bsdf) const {
    if (dot(vertex.geometry_normal, dir_in) < 0) {
        // No light below the surface
        return {};
    }
    // Flip the shading frame if it is inconsistent with the geometry normal
    Frame frame = vertex.shading_frame;
    if (dot(frame.n, dir_in) < 0) {
        frame = -frame;
    }
    
    // Homework 1: implement this!
    return BSDFSampleRecord{
            to_world(frame, sample_cos_hemisphere(rnd_param_uv)),
            Real(0) /* eta */,  eval(bsdf.roughness, vertex.uv, vertex.uv_screen_size, texture_pool)/* roughness */};
}

TextureSpectrum get_texture_op::operator()(const DisneyDiffuse &bsdf) const {
    return bsdf.base_color;
}

1.2. 思考题\text{1.2. 思考题}

1. 将这两种BSDF与Lambertian漫反射对比, 有什么区别？为什么？\text{1. 将这两种BSDF与Lambertian漫反射对比, 有什么区别？为什么？}

2. 修改插值的比例并比较区别。为什么会出现这样的区别？在什么光照下，两种BSDF的差别最为显著？\text{2. 修改插值的比例并比较区别。为什么会出现这样的区别？在什么光照下，两种BSDF的差别最为显著？}

2. Metal\text{2. Metal}

1. Coding\text{1. Coding}

这个是老熟人了，经典的Torrance-Sparrow\text{Torrance-Sparrow}模型配合经典的Trowbridge-Reitz\text{Trowbridge-Reitz}微表面法线分布。不过菲涅尔项并没有使用基于物理的公式，而是使用了Schlick\text{Schlick}近似来计算。

总结一下就是一些耳熟能详的东西：Torrance-Sparrow\text{Torrance-Sparrow}模型，Schlick\text{Schlick}近似，Trowbridge-Reitz\text{Trowbridge-Reitz}微表面法线分布，Smith\text{Smith}几何遮蔽。

但是这里要求我们使用一种更好的重要性采样方法。

理解起来有很多细节，但是如果纯粹是为了完成作业的话可以直接抄slide\text{slide}（逃）。

这里有个重要的问题需要澄清。就是关于指南里的BRDF\text{BRDF}和概率密度问题。

1.1. 不一样的BRDF\text{1.1. 不一样的BRDF}

指南里的BRDF\text{BRDF}是这么写的：

f=FDG4∣n⋅ωi∣f=\frac{FDG}{4|n\cdot \omega_i|}

其他地方绝大多数的BRDF\text{BRDF}相比与这里的都会在分母上多一个cos\text{cos}项。这个其实好理解，只是写法问题，因为这个分母上的cos\text{cos}在渲染方程里会和积分项中用来表示投影的cos\text{cos}自然而然地抵消掉。这么考虑的话进来，指南里的这个写法也没啥毛病（注释里也解释了这一点）。

1.2. pdf也不一样\text{1.2. pdf也不一样}

这个问题就大了，我在学习这种可视法线重要性采样的时候参考了44个地方，概率密度出现了33种写法，看起来挺吓人的，实际上都是一样的。

以下统一用ωi\omega_i表示入射光的方向。

  1. PBR\text{PBR}和Heitz\text{Heitz}的slide\text{slide}如是说：按照可视法线的分布函数采样的概率密度为

Dωi(ωh)=D(ωh)G1(ωi)max⁡(0,∣ωi⋅ωh∣)cos⁡θiD_{\omega_i}(\omega_h)=\frac{D(\omega_h)G_1(\omega_i)\max(0, |\omega_i\cdot \omega_h|)}{\cos \theta_i}

  1. 这篇文章写的是：采样的概率密度是

p(ωi,ωo)=D(ωh)G1(ωi)max⁡(0,∣ωi⋅ωh∣)4∣n⋅ωi∣∣ωh⋅ωo∣p(\omega_i, \omega_o)=\frac{D(\omega_h)G_1(\omega_i)\max(0, |\omega_i\cdot \omega_h|)}{4|n\cdot \omega_i||\omega_h\cdot \omega_o|}

  1. Lab-1\text{Lab-1}的指南上这么写：用这种方法重要性采样

D(ωh)G1(ωi)4∣n⋅ωi∣\frac{D(\omega_h)G_1(\omega_i)}{4|n\cdot \omega_i|}

写成这样其实也很容易看出来了，这几种形式的关联是这样的。

PBR\text{PBR}给出的是半程向量的概率密度，为了转换为出射光的概率密度，乘上Jacobian\text{Jacobian}也就是14∣ωh⋅ωo∣\frac{1}{4|\omega_h\cdot \omega_o|}就是2中写的形式。

而对于2中的形式，由于是半程向量，分子上的∣ωh⋅ωi∣|\omega_h\cdot\omega_i|和分子上的∣ωh⋅ωo∣|\omega_h\cdot\omega_o|实际上相等可以消掉，这就是指南上的形式。

我们实际写代码的时候都以指南为准就不会错了。实在不知道怎么写其实可以看看Lajolla\text{Lajolla}本身的接口和函数实现，比如可视法线采样其实在Lajolla\text{Lajolla}里实现了一个各向同性的版本当然我是写完了改bug的时候才发现。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
#include "../microfacet.h"

Spectrum eval_op::operator()(const DisneyMetal &bsdf) const {
    if (dot(vertex.geometry_normal, dir_in) < 0 ||
            dot(vertex.geometry_normal, dir_out) < 0) {
        // No light below the surface
        return make_zero_spectrum();
    }
    // Flip the shading frame if it is inconsistent with the geometry normal
    Frame frame = vertex.shading_frame;
    if (dot(frame.n, dir_in) < 0) {
        frame = -frame;
    }
    // Homework 1: implement this!
    Spectrum baseColor = eval(bsdf.base_color, vertex.uv, vertex.uv_screen_size, texture_pool);
    Real Roughness = eval(bsdf.roughness, vertex.uv, vertex.uv_screen_size, texture_pool);
    Real Anisotropic = eval(bsdf.anisotropic, vertex.uv, vertex.uv_screen_size, texture_pool);
    Real aspect = sqrt(1.0 - 0.9 * Anisotropic);
    Real alphax = max(0.0001, Roughness * Roughness / aspect);
    Real alphay = max(0.0001, Roughness * Roughness * aspect);
    Vector3 H = normalize(dir_in + dir_out);
    Vector3 Hl = to_local(frame, H);
    Spectrum F = baseColor + (Spectrum(1.0, 1.0, 1.0) - baseColor) * pow(1 - abs(dot(H, dir_out)), 5);
    Real tmp = Hl.x * Hl.x / (alphax * alphax) + Hl.y * Hl.y / (alphay * alphay) + Hl.z * Hl.z;
    Real D = M_1_PI / (alphax * alphay * tmp * tmp);
    Vector3 wi_l = to_local(frame, dir_in);
    Vector3 wo_l = to_local(frame, dir_out);
    Real G_in = 2.0 / (1.0 + sqrt(1.0 + (pow(wi_l.x * alphax, 2) + pow(wi_l.y * alphay, 2)) / (wi_l.z * wi_l.z)));
    Real G_out = 2.0 / (1.0 + sqrt(1.0 + (pow(wo_l.x * alphax, 2) + pow(wo_l.y * alphay, 2)) / (wo_l.z * wo_l.z)));
    Real G = G_in * G_out;
    Spectrum ret = 0.25 * F * D * G / abs(wi_l.z);
    return ret;
}

Real pdf_sample_bsdf_op::operator()(const DisneyMetal &bsdf) const {
    if (dot(vertex.geometry_normal, dir_in) < 0 ||
            dot(vertex.geometry_normal, dir_out) < 0) {
        // No light below the surface
        return 0;
    }
    // Flip the shading frame if it is inconsistent with the geometry normal
    Frame frame = vertex.shading_frame;
    if (dot(frame.n, dir_in) < 0) {
        frame = -frame;
    }
    // Homework 1: implement this!
    Real Roughness = eval(bsdf.roughness, vertex.uv, vertex.uv_screen_size, texture_pool);
    Real Anisotropic = eval(bsdf.anisotropic, vertex.uv, vertex.uv_screen_size, texture_pool);
    Real aspect = sqrt(1.0 - 0.9 * Anisotropic);
    Real alphax = max(0.0001, Roughness * Roughness / aspect);
    Real alphay = max(0.0001, Roughness * Roughness * aspect);
    Vector3 Hl = to_local(frame, normalize(dir_in + dir_out));
    Real tmp = Hl.x * Hl.x / (alphax * alphax) + Hl.y * Hl.y / (alphay * alphay) + Hl.z * Hl.z;
    Real D = M_1_PI / (tmp * tmp * alphax * alphay);
    Vector3 wi_l = to_local(frame, dir_in);
    Real G_in = 2.0 / (1.0 + sqrt(1.0 + (pow(wi_l.x * alphax, 2) + pow(wi_l.y * alphay, 2)) / (wi_l.z * wi_l.z)));
    Real ret = 0.25 * D * G_in / max(0.0001, abs(wi_l.z));
    return ret;
}

std::optional<BSDFSampleRecord>
        sample_bsdf_op::operator()(const DisneyMetal &bsdf) const {
    if (dot(vertex.geometry_normal, dir_in) < 0) {
        // No light below the surface
        return {};
    }
    // Flip the shading frame if it is inconsistent with the geometry normal
    Frame frame = vertex.shading_frame;
    if (dot(frame.n, dir_in) < 0) {
        frame = -frame;
    }
    // Homework 1: implement this!
    Real Roughness = eval(bsdf.roughness, vertex.uv, vertex.uv_screen_size, texture_pool);
    Real Anisotropic = eval(bsdf.anisotropic, vertex.uv, vertex.uv_screen_size, texture_pool);
    Real aspect = sqrt(1.0 - 0.9 * Anisotropic);
    Real alphax = max(0.0001, Roughness * Roughness / aspect);
    Real alphay = max(0.0001, Roughness * Roughness * aspect);
    Vector3 wi_l = to_local(frame, dir_in); // map the incoming light to local frame
    Vector3 wi = normalize(Vector3(wi_l.x * alphax, wi_l.y * alphay, wi_l.z));
    Vector3 T1 = normalize(abs(wi.z) > 0.9999 ?
                        Vector3(1, 0, 0) : cross(Vector3(0, 0, 1), wi));
    Vector3 T2 = cross(wi, T1);
    Real R = sqrt(rnd_param_uv[0]);
    Real phi = c_TWOPI * rnd_param_uv[1];
    Real t1 = R * cos(phi);
    Real t2 = R * sin(phi);
    t2 = 0.5 * (1.0 - wi.z) * sqrt(1.0 - t1 * t1) + 0.5 * (wi.z + 1.0) * t2;
    Vector3 Hl = t1 * T1 + t2 * T2 + sqrt(max(0.0, 1.0 - t1 * t1 - t2 * t2)) * wi;// in local frame but scaled
    Hl = normalize(Vector3(Hl.x * alphax, Hl.y * alphay, max(0.0, Hl.z))); // trans to local frame
    Vector3 wo_l = 2.0 * dot(Hl, wi_l) * Hl - wi_l;
    return BSDFSampleRecord{to_world(frame, wo_l), Real(0), Roughness};
}

TextureSpectrum get_texture_op::operator()(const DisneyMetal &bsdf) const {
    return bsdf.base_color;
}

1.3. 基于物理的公式实际上是不基于物理的？\text{1.3. 基于物理的公式实际上是不基于物理的？}

这里指南中提到了一个问题，算Fresnel\text{Fresnel}项我们用的是Schlick\text{Schlick}近似式而非准确的公式。这么做一部分原因是因为算得更快，另一个原因则是尽管后者似乎会比前者更准确，但是对于用RGB\text{RGB}颜色表示的光（也就是Graphics\text{Graphics}中用的最多的）来说实际上是不准确的。背后的原因其实主要跟颜色背后的物理原理以及色度学有些关系。

指南里给了一份参考，等我看完了就更新。

2. 思考题\text{2. 思考题}

1. 修改roughness参数，除了粗糙度不一样，你还观察到了什么，为什么？\text{1. 修改roughness参数，除了粗糙度不一样，你还观察到了什么，为什么？}

不用实验也猜得到，粗糙度越来越大，表面会越来越暗，闫令琪在GAMES202\text{GAMES202}里提过，这是因为我们没有考虑微表面间多次反射因此这个模型其实是有能量损失的，越粗糙这个损失就越大，表面显得也就越暗，GAMES202\text{GAMES202}讲了一种需要离线预计算的近似算法(当然，说预计算是因为这个预计算就能用在实时渲染里)——Kulla-Conty\text{Kulla-Conty}近似，当然我没做，此处略去。

update on 2022.11.17:\text{update on 2022.11.17:}

之前重新回顾了一下GAMES202\text{GAMES202}相关内容并且看了一下本课程的Microflake Model\text{Microflake Model}部分。Kulla-Conty\text{Kulla-Conty}近似只是常用于工业界实时渲染的近似解决方案。实际上真正解决这个问题的是Heitz\text{Heitz}等人于20162016年提出的Microflake Model\text{Microflake Model}（这相当于是体积渲染中的Microfacet Model\text{Microfacet Model}）。他们证明了Microfacet Model\text{Microfacet Model}的反射行为均可由使用Microflake Model\text{Microflake Model}的光滑平面来模拟，借用体积渲染的方法解决了这个问题。这又是个新的坑了。。。

作业剩下的部分很多就是套公式，没有什么新东西，而且大部分也是hack\text{hack}，就不写了。直接进入重点——体积路径追踪。

Hw-2: Volumetric Path Tracing\text{Hw-2: Volumetric Path Tracing}

体积辐射度传输方程\text{体积辐射度传输方程}

首先看体积渲染的辐射度传输方程：

ddtL(p⃗(t),ω)=−(σa(p⃗(t))+σs(p⃗(t)))L(p⃗(t),ω)+Le(p⃗(t),ω)+σs(p⃗(t))∫S2ρ(p⃗(t),ω,ω′)L(p⃗(t),ω′)dω′\frac{\text{d}}{\text{d}t}L(\vec p(t),\omega) = -(\sigma_a(\vec p(t))+\sigma_s(\vec p(t)))L(\vec p(t), \omega) + L_e(\vec p(t),\omega) + \sigma_s(\vec p(t))\int_{\mathcal{S}^2}\rho(\vec p(t),\omega,\omega')L(\vec p(t), \omega')\text{d}\omega'

这么长，能看出来个锤子啊

其实公式含义非常明确，σa\sigma_a和σs\sigma_s分别是吸收系数和散射系数。LeL_e是（体积的）出射的Radiance\text{Radiance}，ρ\rho是相函数，其作用类似于平面反射中的BRDF\text{BRDF}。

这么一说，其实大体上讲相比于普通的渲染方程也就多了开头那一项。

按照每一项分析一下这个式子在说什么。

考虑从p⃗(t)\vec p(t)点到p⃗(t+dt)\vec p(t+\text{d}t)点，沿ω\omega方向的Radiance\text{Radiance}会有哪些变化。

  1. Radiance\text{Radiance}可能会按照一定比例被吸收：这就是σa(p⃗(t))L(p⃗(t))\sigma_a(\vec p(t))L(\vec p(t))，因为是吸收了，所以符号为负。

  2. Radiance\text{Radiance}可能会按照一定比例被散射到别的方向去：这就是σs(p⃗(t))L(p⃗(t))\sigma_s(\vec p(t))L(\vec p(t))，既然散射到别的方向去了，那沿着ω\omega方向的自然就减少了，所以符号也是负号。

  3. 介质自己会发光：这就是Le(p⃗(t),ω)L_e(\vec p(t),\omega)

  4. 从其他方向的Radiance\text{Radiance}散射到ω\omega方向：我们对整个球面积分，统计所有方向贡献到ω\omega方向的Radiance\text{Radiance}，乘上系数就是最后一项。

加起来就是整个方程了。考虑到前面有一项同时计算了σa+σs\sigma_a+\sigma_s对Radiance\text{Radiance}的削弱，我们可以记σt=σa+σs\sigma_t=\sigma_a+\sigma_s，后面会用。

接下来，我们会由易到难，分66步完成整个体积路径追踪。

2.1. Lajolla中参与介质的数据结构与接口\text{2.1. Lajolla中参与介质的数据结构与接口}

在开始之前先要熟悉Lajolla\text{Lajolla}里怎么表示介质。

其实没啥好说的。。。就是这里有个新的神奇C++\text{C++}语法糖叫std::visit\text{std::visit}

2.2. 仅吸收Radiance的单色均匀介质\text{2.2. 仅吸收Radiance的单色均匀介质}

2.2.1. Coding\text{2.2.1. Coding}

从最简单的开始，这里首先作出44个假设:

  1. 只有一坨介质，并且是均匀的。

  2. 这坨介质不散射光，即σs=0\sigma_s=0。

  3. 场景中的表面只发光，不反射也不折射。

4.这坨介质是单色的， σa\sigma_a的RGB\text{RGB}三通道都是相同的。

于是辐射度传输方程得到了简化：

ddtL(p⃗(t),ω)=−σaL(p⃗(t),ω)\frac{d}{dt}L(\vec p(t), \omega)=- \sigma_aL(\vec p(t), \omega)

这个方程好解:

L(p⃗(t),ω)=L0exp⁡(−σat)L(\vec p(t), \omega)=L_0\exp(-\sigma_a t)

连积分都没有，甚至不需要蒙特卡洛。而且任务书里把伪代码都给你写好了，照着写就行。

事实上也推荐照着伪代码写，因为这里给出的伪代码是可以用于第一个测试场景的，这个场景非常简单，不需要想很多的复杂情况，仅用伪代码的流程就可以完成。

其实翻翻代码看看什么地方用什么接口就行。不过注意如下几个坑：

  1. 生成光线可以用sample_primary\text{sample\_primary}函数，但是注意这个函数传入参数中的向量都是在[0,1][0,1]之间，也就是说对于像素坐标得除以相机的长宽再传。

  2. 光线打到的是灯，外面才是介质，所以应该用exterior_medium\text{exterior\_medium}。

  3. 调用emission\text{emission}函数来获取发光的Radiance\text{Radiance}，但是注意传入参数中光线方向要取反。

代码不难

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
Spectrum vol_path_tracing_1(const Scene &scene,
                            int x, int y, /* pixel coordinates */
                            pcg32_state &rng)
{
    // Homework 2: implememt this!
    Real X = (next_pcg32_real<Real>(rng) + x) / scene.camera.width;
    Real Y = (next_pcg32_real<Real>(rng) + y) / scene.camera.height;
    const Ray ray = sample_primary(scene.camera, Vector2(X, Y));
    std::optional<PathVertex> opt_path_vertex = intersect(scene, ray);
    Spectrum L(0.0, 0.0, 0.0);
    if(opt_path_vertex)
    {
        PathVertex path_vertex = opt_path_vertex.value();
        Real t = distance(ray.org, path_vertex.position);
        Spectrum Le = emission(path_vertex, -ray.dir, scene);
        if(path_vertex.exterior_medium_id >= 0)
        {
          const Medium in_medium = scene.media[path_vertex.exterior_medium_id];
          L += Le * exp(-get_sigma_a(in_medium, path_vertex.position) * t);
        }
    }
    return L;
}

2.2.2. 思考题\text{2.2.2. 思考题}

1. 修改absorption参数，结果有什么不同？为什么？\text{1. 修改absorption参数，结果有什么不同？为什么？}

这简直是送，不用实验也知道absorption\text{absorption}越大图片会越来越暗。。。

2.3. 吸收和散射Radiance的单色均匀介质\text{2.3. 吸收和散射Radiance的单色均匀介质}

2.3.1. Coding\text{2.3.1. Coding}

还是先给出这个简化模型的假设。

  1. 只有一坨介质，并且是均匀的。

  2. 光在这坨介质中只散射一次。

  3. 场景中的表面只发光，不反射也不折射。

  4. 这坨介质是单色的， σa\sigma_a和σs\sigma_s的RGB\text{RGB}三通道都是相同的。

此时辐射度传输方程就变成原始的形式了。但是不慌，积分项还是不难算的。

Lsc1(p⃗(t),ω)=∫Mρ(p⃗(t),ω,ω′)Le(p⃗′,ω′)exp⁡(−σt∥p⃗(t)−p⃗′∥)∣ω′⋅n⃗p⃗′∣∥p⃗(t)−p⃗′∥2visible(p⃗(t),p⃗′)dp⃗′L_{sc_1}(\vec p(t), \omega)=\int_\mathcal M \rho(\vec p(t), \omega, \omega')L_e(\vec p', \omega')\exp(-\sigma_t\parallel \vec p(t)-\vec p' \parallel)\frac{|\omega'\cdot \vec n_{\vec p'}|}{\parallel \vec p(t)-\vec p' \parallel^2}\text{visible}(\vec p(t),\vec p')\text{d}\vec p'

这TM甚至比原来还要难看仔细一看这个积分其实就非常好。

因为我们做出了光线只会散射一次的假设，所以对所有入射方向辐射度的积分可以直接回溯其来源。它们只会散射这一次，也就是说在此之前都是沿直线传播，这个直线传播的过程就回到了我们在2.22.2中解决的模型。

因此，忽略积分项中的相函数和辐射度和指数，最后剩下的部分其实是对发光表面重要性采样的Jacobian\text{Jacobian}，或者更通俗的说，是我们把积分从对球面变成对所有发光表面时的转换系数。

这么一说，这个积分应该怎么计算就很显而易见了，跟路径追踪差不多。

辐射度传输方程于是变为如下：

ddtL(p⃗(t),ω)=−σtL(p⃗(t),ω)+σsLsc1(p⃗(t),ω)\frac{\text{d}}{\text{d}t}L(\vec p(t), \omega)=-\sigma_t L(\vec p(t), \omega) + \sigma_s L_{sc_1}(\vec p(t),\omega)

这东西是个ODE\text{ODE}，但是没有封闭形式了。不过不慌，作为一个ODE\text{ODE}来讲要解这个东西还是不算难。两边乘以exp⁡(tσt)\exp(t\sigma_t)。

exp⁡(tσt)[ddtL(p⃗(t),ω)+σtL(p⃗(t),ω)]=exp⁡(tσt)σsLsc1(p⃗(t),ω)\exp(t\sigma_t)[\frac{\text{d}}{\text{d}t}L(\vec p(t), \omega) +\sigma_t L(\vec p(t), \omega)]= \exp(t\sigma_t)\sigma_s L_{sc_1}(\vec p(t),\omega)

两边积分

exp⁡(tσt)L(p⃗(t),ω)=σs∫0texp⁡(t′σt)Lsc1(p⃗(t′),ω)dt′+L(p⃗(0),ω)\exp(t\sigma_t)L(\vec p(t), \omega)=\sigma_s\int_0^t\exp(t'\sigma_t) L_{sc_1}(\vec p(t'),\omega)\text{d}t'+L(\vec p(0), \omega)

L(p⃗(t),ω)=σs∫0texp⁡((t′−t)σt)Lsc1(p⃗(t′),ω)dt′+exp⁡(−σtt)L(p⃗(0),ω)L(\vec p(t), \omega)=\sigma_s\int_0^t\exp((t'-t)\sigma_t) L_{sc_1}(\vec p(t'),\omega)\text{d}t'+\exp(-\sigma_t t)L(\vec p(0), \omega)

我们之前一直没有说的事情是光的传播方向：虽然单向路径追踪是从相机发光去找光源，实际上光的传输应该是从光源到相机。也就是说，这个方程里t=0t=0表示的是光源。

考虑单向路径追踪里我们是逆向去计算的，为了更加清楚直观，不妨反过来，设相机是t=0t=0点，而光源对应的点为p⃗(thit)\vec p(t_{hit})。这也很清楚，相当于我们从相机发出光线并且打到光源。

因此，我们将上面的方程反过来写

L(p⃗(t),ω)=σs∫tthitexp⁡(−t′σt)Lsc1(p⃗(t′),ω)dt′+exp⁡(−σt(thit−t))L(p⃗(thit),ω)L(\vec p(t), \omega)=\sigma_s\int_{t}^{t_{hit}}\exp(-t'\sigma_t) L_{sc_1}(\vec p(t'),\omega)\text{d}t'+\exp(-\sigma_t (t_{hit}-t))L(\vec p(t_{hit}), \omega)

可能不容易看出来干了什么，其实我们是做了一个简单的换元。虽然数学过程比较隐晦，但是物理含义还是相当直观。

我们要算的其实就是t=0t=0的时候（摄像机收到的Radiance\text{Radiance}）

因此在上面令t=0t=0

L(p⃗(0),ω)=σs∫0thitexp⁡(−t′σt)Lsc1(p⃗(t′),ω)dt′+exp⁡(−σtthit)L(p⃗(thit),ω)L(\vec p(0), \omega)=\sigma_s\int_{0}^{t_{hit}}\exp(-t'\sigma_t) L_{sc_1}(\vec p(t'),\omega)\text{d}t'+\exp(-\sigma_t t_{hit})L(\vec p(t_{hit}), \omega)

这就得到了任务书里的形式（其实直接在原始式子里代入thitt_{hit}然后反过来是一样的），不过那里在积分项里写的是Lsc1(p⃗,ω)L_{sc_1}(\vec p,\omega)而没有写tt应该取什么。

这里我们就又有机会用我们最喜欢的Monte-Carlo\text{Monte-Carlo}积分了。我们来重要性采样exp⁡(−tσt)\exp(-t\sigma_t)。即我们需要采样一个分布p(t)p(t)使得p(t)∝exp⁡(−tσt)p(t) \propto \exp(-t\sigma_t)，这里t>0t>0。

运用inversion method\text{inversion method}这并不难做。设p(t)=cexp⁡(−tσt)p(t)=c\exp(-t\sigma_t)，则有

c∫0+∞exp⁡(−tσt)dt=1c\int_0^{+\infty}\exp(-t\sigma_t)dt=1

轻松解出c=σtc=\sigma_t

故p(t)=σtexp⁡(−tσt)p(t)=\sigma_t \exp(-t\sigma_t)，积分得其累积分布函数：

P(t)=∫0+∞σtexp⁡(−tσt)dt=1−exp⁡(−tσt)P(t)=\int_0^{+\infty}\sigma_t \exp(-t\sigma_t)dt=1-\exp(-t\sigma_t)

对其取逆，得P−1(u)=−ln⁡(1−u)σtP^{-1}(u) = -\frac{\ln(1-u)}{\sigma_t}。

因此我们均匀采样一个分布uu，然后求t=−ln⁡(1−u)σtt=-\frac{\ln(1-u)}{\sigma_t}就可以按照我们希望的指数分布采样tt了。

但是这里需要注意一点。实际上只有在[0,thit][0,t_{hit}]上采样的tt才有意义。对于那些t≥thitt \ge t_{hit}，就直接取光源上的的对应点并且算进发光项。发生t≥thitt\ge t_{hit}的概率就是exp⁡(−thitσt)\exp(-t_{hit}\sigma_t)，容易发现通过这种方式计算出的发光项的期望和我们上面方程里的exp⁡(−σtthit)L(p⃗(thit),ω)\exp(-\sigma_t t_{hit})L(\vec p(t_{hit}), \omega)这一项完全相同。于是我们在上面得到的方程里的发光项刚好就不用算了，直接用这采样方式带来的added bonus\text{added bonus}来替代就行。

代码并不难写，看代码使用已经写好的接口就可以轻松完成，这个过程不复杂，退一万步还有伪代码可以看。

不过伪代码里有一个计算散射Radiance\text{Radiance}的函数需要自己实现，也很简单，使用一些已有的接口对场景中的光源上的点进行采样并计算其pdf\text{pdf}，然后把蒙特卡洛积分式搬进去就行，注意相函数的方向。

不过还是有两个坑要注意：

  1. 在计算进入照相机的Radiance\text{Radiance}时，如果光线与场景没有交点，不能想当然就返回00。因为我们还需要考虑散射至相机的Radiance\text{Radiance}。这时就相当于thitt_{hit}变为正无穷大，一样地进行采样和计算。介质可以通过scene.camera.medium_id\text{scene.camera.medium\_id}获取。

  2. 注意光线的方向！无论是在相函数里还是在计算散射光里这个都很重要。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
Spectrum L_scatter(const Scene& scene, const Vector3& pos, const Vector3& dir_out, pcg32_state& rng)
{
    Real u = next_pcg32_real<Real>(rng);
    int light_id = sample_light(scene, u);
    Vector2 uv(next_pcg32_real<Real>(rng), next_pcg32_real<Real>(rng));
    u = next_pcg32_real<Real>(rng);
    PointAndNormal light_src = sample_point_on_light(scene.lights[light_id], pos, uv, u, scene);
    Ray ray;
    ray.org = pos;
    ray.dir = normalize(light_src.position - pos);
    ray.tnear = Real(0);
    ray.tfar = infinity<Real>();
    std::optional<PathVertex> opt_path_vertex = intersect(scene, ray);
    Spectrum L(0.0, 0.0, 0.0);
    if(opt_path_vertex)
    {
        PathVertex path_vertex = opt_path_vertex.value();
        Real t = distance(pos, path_vertex.position);
        Spectrum Le = emission(path_vertex, pos - light_src.position, scene);//emission() has already considered visibility
        if(path_vertex.exterior_medium_id >= 0)
        {
            Medium medium = scene.media[path_vertex.exterior_medium_id];
            PhaseFunction phase = get_phase_function(medium);
            Real sigma_t = get_sigma_a(medium, path_vertex.position).x
                          + get_sigma_s(medium, path_vertex.position).x;
            Spectrum rho = eval(phase, ray.dir, dir_out);
            Spectrum transmit = rho * Le * exp(-sigma_t * t)
                              * abs(dot(ray.dir, light_src.normal)) / (t * t);
            Real pdf = pdf_point_on_light(scene.lights[light_id], light_src, pos, scene)
                            * light_pmf(scene, light_id);
            L += transmit / pdf;
        }
    }
    return L;
}
Spectrum vol_path_tracing_2(const Scene &scene,
                            int x, int y, /* pixel coordinates */
                            pcg32_state &rng) {
    // Homework 2: implememt this!
    Real X = (next_pcg32_real<Real>(rng) + x) / scene.camera.width;
    Real Y = (next_pcg32_real<Real>(rng) + y) / scene.camera.height;
    const Ray ray = sample_primary(scene.camera, Vector2(X, Y));
    std::optional<PathVertex> opt_path_vertex = intersect(scene, ray);
    Spectrum L(0.0, 0.0, 0.0);
    if(opt_path_vertex)
    {
        PathVertex path_vertex = opt_path_vertex.value();
        Real t_hit = distance(ray.org, path_vertex.position);
        Real u = next_pcg32_real<Real>(rng);
        if(path_vertex.exterior_medium_id >= 0)
        {
            const Medium ext_medium = scene.media[path_vertex.exterior_medium_id];
            Real sigma_s = get_sigma_s(ext_medium, path_vertex.position).x;
            Real sigma_t = get_sigma_a(ext_medium, path_vertex.position).x + sigma_s;
            Real t = -log(1.0 - u) / sigma_t;
            if(t < t_hit)
            {
                Real pdf = sigma_t * exp(-sigma_t * t);
                Spectrum L_sc = L_scatter(scene, ray.org + t * ray.dir, -ray.dir, rng);
                L += sigma_s * exp(-sigma_t * t) * L_sc / pdf;
            }
            else
            {
                Spectrum Le = emission(path_vertex, -ray.dir, scene);
                L += Le;
            }
        }
    }
    else
    {
        Real u = next_pcg32_real<Real>(rng);
        if(scene.camera.medium_id >= 0)
        {
            const Medium& medium = scene.media[scene.camera.medium_id];
            Real sigma_s = get_sigma_s(medium, ray.org).x;
            Real sigma_t = get_sigma_a(medium, ray.org).x + sigma_s;
            Real t = -log(1.0 - u) / sigma_t;
            Real pdf = sigma_t * exp(-sigma_t * t);
            Spectrum L_sc = L_scatter(scene, ray.org + t * ray.dir, -ray.dir, rng);
            L += sigma_s * exp(-sigma_t * t) * L_sc / pdf;
        }
    }
    return L;
}

2.3.2. 思考题\text{2.3.2. 思考题}

1. 修改absorption和scattering参数，结果有什么不同？为什么？\text{1. 修改absorption和scattering参数，结果有什么不同？为什么？}

2.3.3. 加分项\text{2.3.3. 加分项}

指南里说，自己渲染得到的图片可能会比较noisy\text{noisy}，这很正常，有可能是因为没有实施等角采样。

但是为感觉我的图也不noisy\text{noisy}啊，而且我也不知道啥是等角采样啊。。。

好在指南里说推荐我们全都搞完了再来折腾这个加分项，那就先鸽着吧。

2.4. 吸收且多次散射的多重单色均匀介质，仅用相函数采样\text{2.4. 吸收且多次散射的多重单色均匀介质，仅用相函数采样}

这个就有点复杂咯。

2.4.1. Coding\text{2.4.1. Coding}

惯例，先给出这个模型的假设。

  1. 有多坨介质，并且是均匀的。

  2. 场景中的表面只发光，不反射也不折射。

  3. 这些介质是单色的， σa\sigma_a和σs\sigma_s的RGB\text{RGB}三通道都是相同的。

  4. 光可以散射多次，但是我们只根据相函数采样。

先把上面解完ODE\text{ODE}得到的积分式抄下来。

L(p⃗(0),ω)=σs∫0thitexp⁡(−t′σt)Lsc(p⃗(t′),ω)dt′+exp⁡(−σtthit)L(p⃗(thit),ω)L(\vec p(0), \omega)=\sigma_s\int_{0}^{t_{hit}}\exp(-t'\sigma_t) L_{sc}(\vec p(t'),\omega)\text{d}t'+\exp(-\sigma_t t_{hit})L(\vec p(t_{hit}), \omega)

但是现在 LscL_{sc}的形式就几乎退化成原始的形式了，说“几乎”是因为我们的几个系数还是常数，不随空间变化。

Lsc(p⃗(t),ω)=σs∫S2ρ(p⃗(t),ω,ω′)L(p⃗(t),ω′)dω′L_{sc}(\vec p(t),\omega)=\sigma_s\int_{\mathcal{S}^2}\rho(\vec p(t),\omega,\omega')L(\vec p(t), \omega')\text{d}\omega'

我们的采样策略如下：

首先采样一个距离tt，如果t≥thitt\ge t_{hit}那么我们采样光源并且计算光源的贡献，这个很好做。

否则我们计算散射光的贡献，首先根据相函数采样一个方向ω′\omega'，为了计算Radiance\text{Radiance}我们还需要采样一个距离tt再计算Radiance\text{Radiance}，为了计算这个Radiance\text{Radiance}我们再重复这个流程。。。直到我们达到了一个发光表面的时候就终止并且返回。

这个做法是OK的，原因还是我们说过的那个added bonus\text{added bonus}，因为采样距离超过thitt_{hit}的概率刚好是LL前面的系数。我们其实这样采样不仅是在采样距离，也是在以一定概率来采样是否进行相函数采样。

发现了吗？这个过程其实和普通路径追踪的递归计算特别像。

然后就没什么理论上的东西了。指南里还给了几个代码细节上的说明，都说的很明白。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
  Spectrum vol_path_tracing_3(const Scene &scene,
                            int x, int y, /* pixel coordinates */
                            pcg32_state &rng)
{
    // Homework 2: implement this!
    Spectrum cur_path_throughput(1.0, 1.0, 1.0);
    Real X = (next_pcg32_real<Real>(rng) + x) / scene.camera.width;
    Real Y = (next_pcg32_real<Real>(rng) + y) / scene.camera.height;
    Ray ray = sample_primary(scene.camera, Vector2(X, Y));
    int cur_medium_id = scene.camera.medium_id, depth = 0;
    Spectrum L(0.0, 0.0, 0.0);
    while(true)
    {
        if(scene.options.max_depth >= 0 && depth >= scene.options.max_depth)break;
        if(scene.options.rr_depth >= 0 && depth >= scene.options.rr_depth)//first do a Russian Roulette
        {
            Real P_RR = min(cur_path_throughput.x, 0.95);
            Real rnd_rr = next_pcg32_real<Real>(rng);
            if(P_RR < rnd_rr) return L / P_RR;
        }
        std::optional<PathVertex> opt_path_vertex = intersect(scene, ray);
        if(cur_medium_id >= 0)// if we are now in a medium
        {
            Real sigma_s = get_sigma_s(scene.media[cur_medium_id], ray.org).x;
            Real sigma_t = get_sigma_a(scene.media[cur_medium_id], ray.org).x + sigma_s;
            Real u = next_pcg32_real<Real>(rng);
            Real t = -log(1.0 - u) / sigma_t;
            Real pdf = sigma_t * exp(-sigma_t * t);
            if (opt_path_vertex) //if intersect
            {
                PathVertex path_vertex = opt_path_vertex.value();
                Real t_hit = distance(ray.org, path_vertex.position);
                if (t >= t_hit) // we calc the radiance emitted and continue;
                {
                    if(path_vertex.material_id == -1)//simply go through the interface
                    {
                        depth++;
                        ray.org += t * ray.dir;
                        //cur_path_throughput *= exp(-sigma_t * t_hit) / pdf;
                        update_medium(path_vertex, ray, cur_medium_id);
                        //if(cur_medium_id >= 0)
                        //{
                        //    sigma_t = get_majorant(scene.media[cur_medium_id], ray).x;
                        //    cur_path_throughput *= exp(-sigma_t * (t - t_hit));
                        //}
                        continue;
                    }
                    else
                    {
                        L += cur_path_throughput * emission(path_vertex, -ray.dir, scene);
                        break;
                    }
                }
            }//else we sample a direction and calc the scattered radiance
            cur_path_throughput *= exp(-sigma_t * t) / pdf;
            ray.org += t * ray.dir;//now we go to this new place to calc the radiance scattered there
            PhaseFunction phaseFunction = get_phase_function(scene.media[cur_medium_id]);
            Vector2 rnd_uv(next_pcg32_real<Real>(rng), next_pcg32_real<Real>(rng));
            std::optional<Vector3> opt_dir = sample_phase_function(phaseFunction, ray.dir, rnd_uv);
            if (opt_dir)
            {
                Vector3 dir = opt_dir.value();
                Real phase_pdf = pdf_sample_phase(phaseFunction, dir, -ray.dir);
                cur_path_throughput *= sigma_s * eval(phaseFunction, dir, -ray.dir) / phase_pdf;
                ray.dir = dir;
            }
        }
        else
        {
            if(!opt_path_vertex) break;
            else
            {
                PathVertex path_vertex = opt_path_vertex.value();
                L += cur_path_throughput * emission(path_vertex, -ray.dir, scene);
            }
        }
        depth++;
    }
    return L;
}

2.4.2. 思考题\text{2.4.2. 思考题}

1. 修改absorption，scatter和max_depth参数，结果有什么不同？为什么？\text{1. 修改absorption，scatter和max\_depth参数，结果有什么不同？为什么？}

2.5. 用多重重要性采样结合相函数采样和next event estimation\text{2.5. 用多重重要性采样结合相函数采样和next event estimation}

Lab-2\text{Lab-2}进度条过半警告。。。

2.5.1. 多重重要性采样的理解\text{2.5.1. 多重重要性采样的理解}

解释一下，所谓next event estimation\text{next event estimation}，在普通路径追踪里就是我们直接采样一次光源作为直接光。在这里也差不多。

这里相比于2.4\text{2.4}其实只多了两件事：一个是进行next event estimation\text{next event estimation}，另一个是使用多重重要性采样(MIS)\text{(MIS)}将其与相函数采样结合。

要进行MIS\text{MIS}，首先要理解MIS\text{MIS}是在干什么。实际上，重要性采样可以认为是根据一种分布，通过换元的方法给出方差更小的无偏的统计量。而多重重要性采样则是根据多种不同的分布，通过换元并且加权给出方差更小的无偏的统计量。

可以看这里的高赞解答

不管是体积路径追踪还是普通路径追踪，我们其实一般都是分开计算光源的直接光和来自其他表面的间接光。也就是

总辐射度=直接辐射度(来自光源)+间接辐射度\text{总辐射度=直接辐射度(来自光源)+间接辐射度}

首先看普通路径追踪的情况。对于直接光，我们有两种策略来采样，一种是根据BRDF对立体角重要性采样，另一种是对光源表面重要性采样。这是因为直接辐射度在光源以外的表面都是没有的。

而对于间接光采样，我们只有一种重要性采样方式可以使用，那就是按照BRDF对立体角重要性采样。因为对表面进行重要性采样需要我们对于场景中所有表面的辐射度都有所了解，但是我们实际上要算的就是所有表面的辐射度，这是未知的，所以不能重要性采样。

所以，间接光照就好算，而直接光照有两种采样策略。积分项里有辐射度和BRDF的乘积。一个经典的场景是光源很小而BRDF很接近镜面，这种情况下，仅使用BRDF采样和仅使用光源表面采样的统计量都会导致较大的方差（一个直观的理解是，假如只使用BRDF重要性采样，那么大部分情况下很难采样到很小的光源，统计量的值为0。与之相对应的，如果我们很幸运地采样到光源，统计量f(x)p(x)\frac{f(x)}{p(x)}的数量级则会和光源的辐射度相当，统计量的值不“稳定”。仅使用光源重要性采样同理）。

因此，在普通路径追踪里，我们做的事情是

总辐射度=按BRDF采样的直接辐射度(加权)+按光源表面采样的直接辐射度(加权)+间接辐射度\text{总辐射度}=\text{按BRDF采样的直接辐射度(加权)}+\text{按光源表面采样的直接辐射度(加权)}+\text{间接辐射度}

加权是必须加权的，不然同样的直接辐射度我们用两种采样方式估计了两次，结果翻倍了，不是无偏的。

使用NEE\text{NEE}的一种理解是我们先把直接光考虑进来，用BRDF采样出新的方向后就只计算间接光不算直接光了，这种情况下，如果使用BRDF采样新的方向又打到了光源，我们会直接返回00，这相当于给BRDF采样的直接光赋权重为00，而给光源表面采样的直接光赋权重为11 。这样的估计是无偏的，但是可能会导致较大的方差。而在多重重要性采样里，这个权重不再是00和11，而是设计好的相加为11的两个权重，使得估计是无偏的，而且方差能够控制得较小。

那么在体积路径追踪里也是这个道理：

总辐射度=直接辐射度+散射的间接辐射度\text{总辐射度}=\text{直接辐射度}+\text{散射的间接辐射度}

总辐射度=按相函数采样的直接辐射度(加权)+按光源表面采样的直接辐射度(加权)+间接辐射度\text{总辐射度}=\text{按相函数采样的直接辐射度(加权)}+\text{按光源表面采样的直接辐射度(加权)}+\text{间接辐射度}

然后有一些地方的直接光不需要进行MIS\text{MIS}（同样也是不需要采样的地方），一种是第一次从摄像头直接打到光源，另一种是BRDF\text{BRDF}为理想镜面反射或者理想折射的时候，或者光源为点光源的时候。岔开话题说一嘴，这里其实表明δ\delta函数在场景里面是很麻烦的，总是要特殊处理。。。所以lajolla\text{lajolla}中并没有理想的镜面或者点光源这种东西。

2.5.2. 伪代码的理解\text{2.5.2. 伪代码的理解}

有了上面的想法，再来看体积光追里应该如何运用MIS\text{MIS}就更容易理解了。我们知道，MIS\text{MIS}算权重时，用的是使用不同采样方式采样同一样本的pdf\text{pdf}。因此一个pdf\text{pdf}是对光源表面采样的pdf\text{pdf}，另一个pdf\text{pdf}是使用相函数采样这个样本的pdf\text{pdf}。

这里就要注意了，因为中间有介质，我们如果要采到这个点，还需要在介质中按指数分布采样来光线步进。也就是说，如果我们通过相函数采样采这个点，我们需要如下两步：

  1. 首先采样一个方向，刚好对着这个采样点。

  2. 然后沿着这个方向进行指数分布采样，要保证每一步采样都能“跨过”介质的界面。因为按照我们之前的采样方法，假如不能保证“跨过”界面，那么就会采样散射光，这就不是我们要采样的直接光了。

所以我们在光线步进的同时还需要把第二步的概率也累乘起来，这才能在做了NEE\text{NEE}的时候正确加权。

对于NEE\text{NEE}，我们使用的是对光源表面采样来估计直接光，MIS\text{MIS}告诉我们如果使用power heuristics\text{power heuristics}其权重是pnee2pphase2+pnee2\frac{p^2_{nee}}{p^2_{phase}+p^2_{nee}}。

重点分析一下指导中对于NEE\text{NEE}部分的MIS\text{MIS}的伪代码

1
2
3
4
5
6
7
8
9
contrib = T_light * G * f * L / pdf_nee
# Multiple importance sampling: it’s also possible
# that a phase function sampling + multiple exponential sampling
# will reach the light source.
# We also need to multiply with G to convert phase function PDF to area measure.
pdf_phase = pdf_sample_phase(phase_function, dir_view, dir_light) * G * p_trans_dir
# power heuristics
w = (pdf_nee * pdf_nee) / (pdf_nee * pdf_nee + pdf_phase * pdf_phase)
return w * contrib

根据代码前文，p_trans_dir累乘了每一步采样穿过介质界面的概率，因此，我们先计算出使用相函数采样到该方向的概率，乘上p_trans_dir，在已知pdf_nee以及Jacobian\text{Jacobian}项G的情况下就可以算出来采样的权重。

使用另一种采样方式（也就是相函数采样）采样到直接光时应当赋予的权重是pphase2pphase2+pnee2\frac{p^2_{phase}}{p^2_{phase}+p^2_{nee}}。注意这和前面NEE\text{NEE}采样的pp不是一个东西，因为我们做NEE\text{NEE}采的样本和我们用相函数采样的样本不是同一个样本。

那么对于相函数采样，我们也需要计算赋予的权重，这就需要知道pphasep_{phase}和pneep_{nee}，前者是在采样后或者采样时就可以算出来的，但是后者不然。

当我们使用相函数采样打到一个光源的时候，我们又需要计算直接光的光照，这时只可能是相函数采样采到了此光源。我们要计算使用NEE\text{NEE}采样采到此点的概率密度。

与普通路径追踪不同的是，在体积路径追踪中，我们不能直接用上一次的采样点来算NEE\text{NEE}的概率密度，因为我们进行距离采样而没有进行相函数采样（也就是穿过了介质的界面）的时候是没有进行NEE\text{NEE}的，这个时候光线的传播方向没有改变，仅仅只是辐射度衰减了。因此，如果我们要考虑当前采到的光源相对于哪个点是直接光照，那自然是上一次发生相函数采样的点，这个点可能存在于几次采样之前，因此我们需要将其记录下来，才能在采样到光源时准确计算使用NEE\text{NEE}采样直接光的概率密度。

2.5.3. Coding\text{2.5.3. Coding}

一个是对光源进行采样的next event estimation\text{next event estimation}的概率密度，这个其实就是所有光源可见面积的倒数，可以用接口直接算。

相函数采样的概率密度也有接口可以直接算。至于穿过介质的概率，这个在前面已经推导过了，是个简单的指数函数。

所以说我们的思路是非常清晰的。。。但是写好以后debug这个程序花了我半辈子的功夫，前前后后检查path tracing\text{path tracing}部分和next event estimation\text{next event estimation}的部分耗时不下几十个小时，甚至还动用了github上面已有的正确代码来对着查错。。。

最后的最后，发现这两个部分都没错，错的是更新介质的函数update_medium搞错了，又一次忘记了光线的方向要反过来。。。感觉自己这种大聪明可以换个地球生活了下次再犯这种错误就别说自己是搞Graphics\text{Graphics}的了，Graphics\text{Graphics}没有这样的制杖（

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
void update_medium(const PathVertex& pathVertex, const Ray& ray, int& cur_medium_id)
{
    if(pathVertex.exterior_medium_id != pathVertex.interior_medium_id)
    {
        if (dot(ray.dir, pathVertex.geometry_normal) >= 0)
            cur_medium_id = pathVertex.exterior_medium_id;
        else cur_medium_id = pathVertex.interior_medium_id;
    }
}
Real Geometry(const PointAndNormal& pointAndNormal, const Vector3& pos)
{
    Vector3 dif = pos - pointAndNormal.position;
    Real dist = length(dif);
    dif /= dist;
    return max(0.0, dot(dif, pointAndNormal.normal)) / (dist * dist);
}
Spectrum next_event_estimation(const Ray& ray, int depth, int cur_medium_id, const Scene& scene, pcg32_state& rng)
{
    int light_id = sample_light(scene, next_pcg32_real<Real>(rng));
    const PointAndNormal& light_PN = sample_point_on_light(scene.lights[light_id], ray.org,
                                                    Vector2(next_pcg32_real<Real>(rng), next_pcg32_real<Real>(rng)), next_pcg32_real<Real>(rng), scene);
    Ray nee_ray;
    nee_ray.org = ray.org;
    nee_ray.dir = normalize(light_PN.position - ray.org);
    nee_ray.tnear = get_intersection_epsilon(scene);
    nee_ray.tfar = infinity<Real>();
    const PhaseFunction& phaseFunction = get_phase_function(scene.media[cur_medium_id]);
    Real nee_throughput = eval(phaseFunction, -ray.dir, nee_ray.dir).x;
    int nee_dep = 0;
    Real pdf_nee = pdf_point_on_light(scene.lights[light_id], light_PN, ray.org, scene) * light_pmf(scene, light_id);
    Real p_march = 1.0;//reduce as we pass through different media
    Real G = Geometry(light_PN, ray.org);
    while(true)//we do a next event estimation
    {
        if(scene.options.max_depth >= 0 && depth + nee_dep >= scene.options.max_depth)
            return make_zero_spectrum();
        std::optional<PathVertex> opt_pathVertex = intersect(scene, nee_ray);
        if(opt_pathVertex)
        {
            PathVertex pathVertex = opt_pathVertex.value();
            Real t = distance(pathVertex.position, nee_ray.org);
            if(cur_medium_id >= 0)
            {
                Real sigma_a = get_sigma_a(scene.media[cur_medium_id], nee_ray.org).x;
                Real sigma_s = get_sigma_s(scene.media[cur_medium_id], nee_ray.org).x;
                Real sigma_t = sigma_a + sigma_s;
                nee_throughput *= exp(-sigma_t * t);
                p_march *= exp(-sigma_t * t);
            }
            nee_ray.org = pathVertex.position;
            if(pathVertex.material_id >= 0 && get_area_light_id(scene.shapes[pathVertex.shape_id]) != light_id)
                return make_zero_spectrum();
            if(distance(nee_ray.org, light_PN.position) < get_intersection_epsilon(scene)) // we reach the light source, and we calc the contrib
            {
                Spectrum L_nee = nee_throughput * emission(pathVertex, -nee_ray.dir, scene) * G / pdf_nee;
                Real pdf_phase = pdf_sample_phase(phaseFunction, -nee_ray.dir, ray.dir)
                                    * G * p_march;
                Real w = pdf_nee * pdf_nee / (pdf_nee * pdf_nee + pdf_phase * pdf_phase);
                return L_nee * w;
            }
            update_medium(pathVertex, nee_ray, cur_medium_id);
            nee_dep++;
        }
        else return make_zero_spectrum();
    }
}
Spectrum vol_path_tracing_4(const Scene &scene,
                            int x, int y, /* pixel coordinates */
                            pcg32_state &rng) {
    // Homework 2: implement this!
    Spectrum cur_path_throughput(1.0, 1.0, 1.0);
    Real X = (static_cast<Real>(x) + next_pcg32_real<Real>(rng)) / scene.camera.width;
    Real Y = (static_cast<Real>(y) + next_pcg32_real<Real>(rng))/ scene.camera.height;
    Ray ray = sample_primary(scene.camera, Vector2(X, Y));
    int cur_medium_id = scene.camera.medium_id, depth = 0;
    Spectrum L(0.0, 0.0, 0.0);
    Vector3 pos_nee_cache;
    Real p_march = 1.0;
    Vector3 dir_in_cache;
    int medium_id_cache = -1;
    bool never_scatter = true;
    Real p_path = 1.0;
    while(true)
    {
        if(scene.options.max_depth >= 0 && depth >= scene.options.max_depth)break;
        if(scene.options.rr_depth >= 0 && depth >= scene.options.rr_depth)//first do a Russian Roulette
        {
            Real rnd_rr = next_pcg32_real<Real>(rng);
            Real P_RR = std::min(0.95, cur_path_throughput.x / p_path);
            if(rnd_rr > P_RR) break;
            else cur_path_throughput /= P_RR;
        }
        std::optional<PathVertex> opt_path_vertex = intersect(scene, ray);
        if(cur_medium_id >= 0)// if we are now in a medium
        {
            Real t_hit;
            PathVertex path_vertex;
            if(!opt_path_vertex) t_hit = infinity<Real>();
            else
            {
                path_vertex = opt_path_vertex.value();
                t_hit = distance(ray.org, path_vertex.position);
            }
            Real sigma_s = get_sigma_s(scene.media[cur_medium_id], ray.org).x;
            Real sigma_t = get_sigma_a(scene.media[cur_medium_id], ray.org).x + sigma_s;
            Real u = next_pcg32_real<Real>(rng);
            Real t = -log(1.0 - u) / sigma_t;
            if (t >= t_hit)
            {
                p_march *= exp(-t_hit * sigma_t);
                p_path *= exp(-t_hit * sigma_t);
                if(path_vertex.material_id == -1)//simply go through the interface
                {
                    depth++;
                    ray.org = path_vertex.position;
                    update_medium(path_vertex, ray, cur_medium_id);
                    continue;
                }
                else // perform MIS
                {
                    PointAndNormal light_PN = {path_vertex.position, path_vertex.geometry_normal};
                    int light_id = get_area_light_id(scene.shapes[path_vertex.shape_id]);
                    if(never_scatter && light_id != -1)
                        L += cur_path_throughput * emission(path_vertex, -ray.dir, scene);
                    else if(medium_id_cache != -1 && light_id != -1)
                    {
                        Real G = Geometry(light_PN, pos_nee_cache);
                        PhaseFunction phaseFunction = get_phase_function(scene.media[medium_id_cache]);
                        Real pdf_dir = pdf_sample_phase(phaseFunction, dir_in_cache, ray.dir)
                                       * p_march * G;
                        Real pdf_nee = pdf_point_on_light(scene.lights[light_id], light_PN,
                                                          pos_nee_cache, scene) * light_pmf(scene, light_id);
                        Real w = pdf_dir * pdf_dir / (pdf_dir * pdf_dir + pdf_nee * pdf_nee);
                        L += cur_path_throughput * emission(path_vertex, -ray.dir, scene) * w;
                    }
                    else break;
                }
            }
            else//we perform nee and phase function sampling
            {
                p_path *= sigma_t;
                cur_path_throughput *= sigma_s / sigma_t;
                ray.org += t * ray.dir;//now we go to this new place to calc the radiance scattered there
                Spectrum L_nee = next_event_estimation(ray, depth, cur_medium_id, scene, rng);
                L += cur_path_throughput * L_nee;
                medium_id_cache = cur_medium_id;
                dir_in_cache = -ray.dir;
                pos_nee_cache = ray.org;
                never_scatter = false;
                p_march = 1.0;
                PhaseFunction phaseFunction = get_phase_function(scene.media[cur_medium_id]);
                Vector2 rnd_uv(next_pcg32_real<Real>(rng), next_pcg32_real<Real>(rng));
                std::optional<Vector3> opt_dir = sample_phase_function(phaseFunction, ray.dir, rnd_uv);
                if (opt_dir)
                {
                    Vector3 dir = opt_dir.value();
                    Real phase_pdf = pdf_sample_phase(phaseFunction, dir, -ray.dir);
                    p_path *= phase_pdf;
                    cur_path_throughput *= eval(phaseFunction, dir, -ray.dir) / phase_pdf;
                    ray.dir = dir;
                }
                else break;
            }
        }
        else
        {
            if(!opt_path_vertex) break;
            PathVertex path_vertex = opt_path_vertex.value();
            PointAndNormal light_PN = {path_vertex.position, path_vertex.geometry_normal};
            Real t_hit = distance(ray.org, path_vertex.position);
            int light_id = get_area_light_id(scene.shapes[path_vertex.shape_id]);
            if(light_id >= 0)
            {
                if(never_scatter)
                    L += cur_path_throughput * emission(path_vertex, -ray.dir, scene);
                else if(medium_id_cache != -1)
                {
                    Real G = Geometry(light_PN, pos_nee_cache);
                    PhaseFunction phaseFunction = get_phase_function(scene.media[medium_id_cache]);
                    Real pdf_dir = pdf_sample_phase(phaseFunction, dir_in_cache, ray.dir)
                                   * p_march * G;
                    Real pdf_nee = pdf_point_on_light(scene.lights[light_id], light_PN,
                                                      pos_nee_cache, scene) * light_pmf(scene, light_id);
                    Real w = pdf_dir * pdf_dir / (pdf_dir * pdf_dir + pdf_nee * pdf_nee);
                    L += cur_path_throughput * emission(path_vertex, -ray.dir, scene) * w;
                }
                else break;
            }
            ray.org += t_hit * ray.dir;
            update_medium(path_vertex, ray, cur_medium_id);
        }
        depth++;
    }
    return L;
}

2.5.4. 思考题\text{2.5.4. 思考题}

1. 什么时候使用next event estimation会比相函数采样更为有效？\text{1. 什么时候使用next event estimation会比相函数采样更为有效？}

 在我们的测试用例里，哪一个更加高效？为什么？\text{ 在我们的测试用例里，哪一个更加高效？为什么？}

个人猜测，当光源面积较小，或者整个场景的σt\sigma_t都比较小，使用next event estimation\text{next event estimation}更有效，前者是因为光源小的时候，使用相函数采样很难命中光源，后者是因为散射系数很大的时候会频繁发生散射。

在我们的测试用例里面，第一个是小光源但是吸收系数比较正常，全空间都有介质，第二个是小光源并且吸收系数大得离谱，大部分空间是真空。自己测试了一下，控制两个场景的spp\text{spp}相同，然后。。。目测一下渲染结果的噪音大不大。目测感觉第二个更高效，如果要解释原因的话。。。可能是因为大部分空间都是真空，小部分空间是σt\sigma_t极大的球，这种情况下，相函数采样的贡献非常之小（因为pneep_{nee}几乎为00），MIS\text{MIS}差不多退化为完全的next event estimation\text{next event estimation}，单独的这一种采样方式会比结合两种采样方式更为有效…吧。

2. 我们渲染了一个非常致密的介质球，它和直接采用Lambertian模型渲染的球有何区别？\text{2. 我们渲染了一个非常致密的介质球，它和直接采用Lambertian模型渲染的球有何区别？}

 他们之间为什么会有相似或者不同?\text{ 他们之间为什么会有相似或者不同?}

这点很好理解，介质非常致密(σt\sigma_t非常大)的时候，辐射度进入介质内部很快就衰减到非常小，因此，光在介质内部的行进几乎可以忽略不计，主要的辐射度都来自于在介质分界面处一个薄层的散射，这就相当于是表面反射一样，相函数就相当于BSDF。不过，在原理上讲，体积光照还是比Lambertian\text{Lambertian}模型要复杂很多，可以认为Lambertian\text{Lambertian}模型是这种体积渲染的极限情况。

3. Jim Kajiya曾于1991年预言，10年以内所有的渲染都会是体积渲染。你认为他的理由是什么？\text{3. Jim Kajiya曾于1991年预言，10年以内所有的渲染都会是体积渲染。你认为他的理由是什么？}

 为什么他的预言没有实现？\text{ 为什么他的预言没有实现？}

第二个问题其实部分解释了这个事情，表面光照其实可以看成体积渲染简化到极致的情形，体积渲染算法会比表面光照算法复杂很多（无论是实时还是离线），需要的数据量和运算量也大很多。不过，如果硬件和软件条件允许，采用体积渲染的图像肯定更加准确（表面光照只考虑了表面的情况，没有考虑物体内部的光传输情况，采用体积渲染，这个问题就直接解决了）。不过现在的硬件和软件都还没有达到这个能够比较普遍地支撑体积渲染的地步，所以目前还是表面光照居多。

至于为什么他的预言尚未实现。。。可能是现在的硬件和软件都不大行？我也不知道我TM要是知道我还会在这里学这个课啊，你TM最好先把我拉进实验室再来问我怎么理解

不过某种意义上也实现了，Jim Kajiya当年肯定没想到Vision人能搞出来一个叫Nerf的东西

2.6. 加入表面光照\text{2.6. 加入表面光照}

有了前面的铺垫，再来考虑这个就非常好想了，我们不过就是把前面的辐射传输方程中涉及到光源辐射度的部分都换成场景表面的辐射度，而场景表面的辐射度可以用渲染方程来算。指南里没有给出具体细节，不过这个其实也并不难想就是了。

2.6.1. Coding\text{2.6.1. Coding}

如果我说的是“不难想”而不是“不难”的话，这意味着很难写

当然啊，实际操作起来细节极多，要把几个测试场景都跑出来非常不容易。。。特别是测试场景5−25-2以及茶壶场景都有透明绝缘界面包裹着介质的情况，在计算表面光照（无论是做NEE\text{NEE}还是采样）之前都记得要更新当前介质。

而且写这份代码的时候还遇到了很大的数值精度问题。。。用于测试的Cornell Box\text{Cornell Box}天花板不知道为什么会有两个矩形，光线相交就出现了误差，删掉一个矩形天花板才亮起来。

下面的代码写得非常草率，实话说很多代码都肉眼可见地可以精简，但是本着能跑就不要动其实是懒的精神，就原模原样地放上来咯。

以后要是有机会重写的话再写一份更好的吧（

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
Spectrum next_event_estimation_bsdf(const Ray& ray, int depth, const PathVertex& path_vertex, int cur_medium_id,
                                    const Scene& scene, pcg32_state& rng)
{
    int light_id = sample_light(scene, next_pcg32_real<Real>(rng));
    PointAndNormal light_PN = sample_point_on_light(scene.lights[light_id], ray.org,
                                                    Vector2(next_pcg32_real<Real>(rng), next_pcg32_real<Real>(rng)),
                                                    next_pcg32_real<Real>(rng), scene);
    Ray nee_ray;
    nee_ray.org = ray.org;

    nee_ray.dir = normalize(light_PN.position - ray.org);
    nee_ray.tnear = get_intersection_epsilon(scene);
    nee_ray.tfar = infinity<Real>();
    const Material& mat = scene.materials[path_vertex.material_id];
    Spectrum nee_throughput = eval(mat, -ray.dir, nee_ray.dir, path_vertex, scene.texture_pool);
    int nee_dep = 0;
    Real pdf_nee = pdf_point_on_light(scene.lights[light_id], light_PN, ray.org, scene) * light_pmf(scene, light_id);
    Real p_march = 1.0;//reduce as we pass through different media
    Real G = Geometry(light_PN, ray.org);
    update_medium(path_vertex, nee_ray, cur_medium_id);
    while(true)//we do a next event estimation
    {
        if(scene.options.max_depth >= 0 && depth + nee_dep >= scene.options.max_depth)
            return make_zero_spectrum();
        std::optional<PathVertex> opt_pathVertex = intersect(scene, nee_ray);
        if(opt_pathVertex)
        {
            PathVertex pathVertex = opt_pathVertex.value();
            Real t = distance(pathVertex.position, nee_ray.org);
            if(cur_medium_id >= 0)
            {
                Spectrum sigma_a = get_sigma_a(scene.media[cur_medium_id], nee_ray.org);
                Spectrum sigma_s = get_sigma_s(scene.media[cur_medium_id], nee_ray.org);
                Spectrum sigma_t = sigma_a + sigma_s;
                nee_throughput *= exp(-sigma_t.x * t);
                p_march *= exp(-sigma_t.x * t);
            }
            if(pathVertex.material_id >= 0 && get_area_light_id(scene.shapes[pathVertex.shape_id]) != light_id)
                return make_zero_spectrum();
            nee_ray.org = pathVertex.position;
            if(distance(nee_ray.org, light_PN.position) < get_intersection_epsilon(scene)) // we reach the light source, and we calc the contrib
            {
                Spectrum L_nee = nee_throughput * emission(pathVertex, -nee_ray.dir, scene) * G / pdf_nee;
                Real pdf_bsdf = pdf_sample_bsdf(mat, -ray.dir, nee_ray.dir, path_vertex, scene.texture_pool)
                                    * G * p_march;
                Real w = pdf_nee * pdf_nee / (pdf_nee * pdf_nee + pdf_bsdf * pdf_bsdf);
                return L_nee * w;
            }
            update_medium(pathVertex, nee_ray, cur_medium_id);
            nee_dep++;
        }
        else return make_zero_spectrum();
    }
}
// The fifth volumetric renderer: 
// multiple monochromatic homogeneous volumes with multiple scattering
// with MIS between next event estimation and phase function sampling
// with surface lighting
Spectrum vol_path_tracing_5(const Scene &scene,
                            int x, int y, /* pixel coordinates */
                            pcg32_state &rng) {
    // Homework 2: implement this!
    Spectrum cur_path_throughput(1.0, 1.0, 1.0);
    Real X = (static_cast<Real>(x) + next_pcg32_real<Real>(rng)) / scene.camera.width;
    Real Y = (static_cast<Real>(y) + next_pcg32_real<Real>(rng))/ scene.camera.height;
    Ray ray = sample_primary(scene.camera, Vector2(X, Y));
    ray.tnear = get_intersection_epsilon(scene);
    int cur_medium_id = scene.camera.medium_id, depth = 0;
    Spectrum L(0.0, 0.0, 0.0);
    Vector3 pos_nee_cache;
    Real p_march = 1.0;
    Vector3 dir_in_cache;
    int medium_id_cache = -1;
    PathVertex pv_nee_cache;
    pv_nee_cache.material_id = -1;
    bool never_scatter_or_reflect = true;
    Real p_path = 1.0;
    while(true)
    {
        if(scene.options.max_depth >= 0 && depth >= scene.options.max_depth)break;
        if(scene.options.rr_depth >= 0 && depth >= scene.options.rr_depth)//first do a Russian Roulette
        {
            Real rnd_rr = next_pcg32_real<Real>(rng);
            Real P_RR;
            if(0.95 * p_path <= cur_path_throughput.x) P_RR = 0.95;
            else P_RR = cur_path_throughput.x / p_path;
            if(rnd_rr > P_RR) break;
            else cur_path_throughput /= P_RR;
        }
        std::optional<PathVertex> opt_path_vertex = intersect(scene, ray);
        if(cur_medium_id >= 0)// if we are now in a medium
        {
            Real t_hit;
            PathVertex path_vertex;
            if(!opt_path_vertex) t_hit = infinity<Real>();
            else
            {
                path_vertex = opt_path_vertex.value();
                t_hit = distance(ray.org, path_vertex.position);
            }
            Spectrum sigma_s = get_sigma_s(scene.media[cur_medium_id], ray.org);
            Spectrum sigma_t = get_sigma_a(scene.media[cur_medium_id], ray.org) + sigma_s;
            Real u = next_pcg32_real<Real>(rng);
            Real t = -log(1.0 - u) / sigma_t.x;
            if (t >= t_hit)
            {
                p_march *= exp(-t_hit * sigma_t.x);
                p_path *= exp(-t_hit * sigma_t.x);
                if(path_vertex.material_id == -1)//it is a index-matching surface, we go through the interface
                {
                    depth++;
                    ray.org = path_vertex.position;
                    update_medium(path_vertex, ray, cur_medium_id);
                    continue;
                }
                else // if it is a light then perform MIS, else perform surface lighting
                {
                    PointAndNormal surface_PN = {path_vertex.position, path_vertex.geometry_normal};
                    int light_id = get_area_light_id(scene.shapes[path_vertex.shape_id]);
                    if(light_id != -1)
                    {
                        if(never_scatter_or_reflect) L += cur_path_throughput * emission(path_vertex, -ray.dir, scene);
                        else
                        {
                            Real pdf_dir, pdf_nee;
                            if(medium_id_cache != -1)// perform MIS for phase function sampling in scattering
                            {
                                Real G = Geometry(surface_PN, pos_nee_cache);
                                PhaseFunction phaseFunction = get_phase_function(scene.media[medium_id_cache]);
                                pdf_dir = pdf_sample_phase(phaseFunction, dir_in_cache, ray.dir)
                                               * p_march * G;
                                pdf_nee = pdf_point_on_light(scene.lights[light_id], surface_PN,
                                                                  pos_nee_cache, scene) * light_pmf(scene, light_id);
                            }
                            else if(pv_nee_cache.material_id != -1)// perform MIS for bsdf sampling in reflecting
                            {
                                Real G = Geometry(surface_PN, pv_nee_cache.position);
                                const Material& mat = scene.materials[pv_nee_cache.material_id];
                                pdf_dir = pdf_sample_bsdf(mat, dir_in_cache, ray.dir, pv_nee_cache, scene.texture_pool)
                                          * p_march * G;
                                pdf_nee = pdf_point_on_light(scene.lights[light_id], surface_PN, pv_nee_cache.position, scene)
                                          * light_pmf(scene, light_id);
                            }
                            else break;
                            Real w = pdf_dir * pdf_dir / (pdf_dir * pdf_dir + pdf_nee * pdf_nee);
                            L += cur_path_throughput * emission(path_vertex, -ray.dir, scene) * w;
                            break;
                        }
                    }
                    else //not light or index-macthing surface, so it must be a solid surface, perform nee
                    {
                        never_scatter_or_reflect = false;
                        ray.org = path_vertex.position;
                        Spectrum L_nee = next_event_estimation_bsdf(ray, depth, path_vertex, cur_medium_id, scene, rng);
                        L += cur_path_throughput * L_nee;
                        medium_id_cache = -1;
                        pv_nee_cache = path_vertex;
                        dir_in_cache = -ray.dir;
                        pos_nee_cache = ray.org;
                        p_march = 1.0;
                        const Material& mat = scene.materials[path_vertex.material_id];
                        Vector2 rnd_uv(next_pcg32_real<Real>(rng), next_pcg32_real<Real>(rng));
                        std::optional<BSDFSampleRecord> opt_bRec = sample_bsdf(mat, -ray.dir, path_vertex, scene.texture_pool,
                                                                               rnd_uv, next_pcg32_real<Real>(rng));
                        if (opt_bRec)
                        {
                            BSDFSampleRecord bRec = opt_bRec.value();
                            Real bsdf_pdf = pdf_sample_bsdf(mat, -ray.dir, bRec.dir_out, path_vertex, scene.texture_pool);
                            if(bsdf_pdf == 0.0) break;
                            cur_path_throughput *= eval(mat, -ray.dir, bRec.dir_out, path_vertex, scene.texture_pool) / bsdf_pdf;
                            p_path *= bsdf_pdf;
                            ray.dir = bRec.dir_out;
                            update_medium(path_vertex, ray, cur_medium_id);
                        }
                        else break;
                    }
                }
            }
            else// we perform nee and phase function sampling
            {
                p_path *= sigma_t.x;
                cur_path_throughput *= sigma_s.x / sigma_t.x;
                ray.org += t * ray.dir;//now we go to this new place to calc the radiance scattered there
                Spectrum L_nee = next_event_estimation(ray, depth, cur_medium_id, scene, rng);
                L += cur_path_throughput * L_nee;
                medium_id_cache = cur_medium_id;
                pv_nee_cache.material_id = -1;
                dir_in_cache = -ray.dir;
                pos_nee_cache = ray.org;
                never_scatter_or_reflect = false;
                p_march = 1.0;
                PhaseFunction phaseFunction = get_phase_function(scene.media[cur_medium_id]);
                Vector2 rnd_uv(next_pcg32_real<Real>(rng), next_pcg32_real<Real>(rng));
                std::optional<Vector3> opt_dir = sample_phase_function(phaseFunction, -ray.dir, rnd_uv);
                if (opt_dir)
                {
                    Vector3 dir = opt_dir.value();
                    Real phase_pdf = pdf_sample_phase(phaseFunction, -ray.dir, dir);
                    p_path *= phase_pdf;
                    cur_path_throughput *= eval(phaseFunction, -ray.dir, dir) / phase_pdf;
                    ray.dir = dir;
                }
                else break;
            }
        }
        else
        {
            if(!opt_path_vertex) break;
            PathVertex path_vertex = opt_path_vertex.value();
            PointAndNormal surface_PN = {path_vertex.position, path_vertex.geometry_normal};
            int light_id = get_area_light_id(scene.shapes[path_vertex.shape_id]);
            if(light_id >= 0)
            {
                if(never_scatter_or_reflect) L += cur_path_throughput * emission(path_vertex, -ray.dir, scene);
                else
                {
                    Real pdf_dir, pdf_nee;
                    if(medium_id_cache != -1)// perform MIS for phase function sampling in scattering
                    {
                        Real G = Geometry(surface_PN, pos_nee_cache);
                        PhaseFunction phaseFunction = get_phase_function(scene.media[medium_id_cache]);
                        pdf_dir = pdf_sample_phase(phaseFunction, dir_in_cache, ray.dir)
                                    * p_march * G;
                        pdf_nee = pdf_point_on_light(scene.lights[light_id], surface_PN, pos_nee_cache, scene)
                                                        * light_pmf(scene, light_id);
                    }
                    else if(pv_nee_cache.material_id != -1)// perform MIS for bsdf sampling in reflecting
                    {
                        Real G = Geometry(surface_PN, pv_nee_cache.position);
                        const Material& mat = scene.materials[pv_nee_cache.material_id];
                        pdf_dir = pdf_sample_bsdf(mat, dir_in_cache, ray.dir, pv_nee_cache, scene.texture_pool)
                                    * p_march * G;
                        pdf_nee = pdf_point_on_light(scene.lights[light_id], surface_PN, pv_nee_cache.position, scene)
                                    * light_pmf(scene, light_id);
                    }
                    else break;
                    Real w = pdf_dir * pdf_dir / (pdf_dir * pdf_dir + pdf_nee * pdf_nee);
                    L += cur_path_throughput * emission(path_vertex, -ray.dir, scene) * w;
                }
                break;
            }
            else if(path_vertex.material_id >= 0)
            {
                never_scatter_or_reflect = false;
                ray.org = path_vertex.position;
                Spectrum L_nee = next_event_estimation_bsdf(ray, depth, path_vertex, cur_medium_id, scene, rng);
                L += cur_path_throughput * L_nee;
                medium_id_cache = -1;
                pv_nee_cache = path_vertex;
                dir_in_cache = -ray.dir;
                pos_nee_cache = ray.org;
                p_march = 1.0;
                const Material& mat = scene.materials[path_vertex.material_id];
                Vector2 rnd_uv(next_pcg32_real<Real>(rng), next_pcg32_real<Real>(rng));
                std::optional<BSDFSampleRecord> opt_bRec = sample_bsdf(mat, -ray.dir, path_vertex, scene.texture_pool,
                                                                       rnd_uv, next_pcg32_real<Real>(rng));
                if (opt_bRec)
                {
                    BSDFSampleRecord bRec = opt_bRec.value();
                    Real bsdf_pdf = pdf_sample_bsdf(mat, -ray.dir, bRec.dir_out, path_vertex, scene.texture_pool);
                    if(bsdf_pdf == 0.0) break;
                    cur_path_throughput *= eval(mat, -ray.dir, bRec.dir_out, path_vertex, scene.texture_pool) / bsdf_pdf;
                    p_path *= bsdf_pdf;
                    ray.dir = bRec.dir_out;
                }
                else break;
            }
            else ray.org = path_vertex.position;
            update_medium(path_vertex, ray, cur_medium_id);
        }
        depth++;
    }
    return L;
}

2.6.2. 思考题\text{2.6.2. 思考题}

1. 改变绝缘界面材料的折射率，这如何影响渲染结果？为什么？\text{1. 改变绝缘界面材料的折射率，这如何影响渲染结果？为什么？}

2. 在茶壶场景中，茶壶是用玻璃包裹介质渲染的，这和用纯玻璃茶壶渲染的区别是什么？\text{2. 在茶壶场景中，茶壶是用玻璃包裹介质渲染的，这和用纯玻璃茶壶渲染的区别是什么？}

都不知道，以后再说

2.7. 非均匀多色介质\text{2.7. 非均匀多色介质}

终于来到最后一阶段了。。。也是最复杂的。

原本我们一直认为σa,σs,σt\sigma_a,\sigma_s,\sigma_t都是不随空间变化的标量。而现在，一方面我们不能认为σt\sigma_t不随空间变化（这样我们不能再把传输项里的积分写成封闭形式了），另一方面由于σt\sigma_t对于各色光分量不同，我们也不能将多色光视为同一束光来考虑了。有了前面的所有铺垫，下面我们就来解决这两个问题。

2.7.1. Null Scattering\text{2.7.1. Null Scattering}

考虑传输项的定义为

τ(p⃗(0),p⃗(t))=exp⁡(∫0tσt(p⃗(t′))dt′)\tau(\vec{p}(0), \vec{p}(t))=\exp(\int_0^t\sigma_t(\vec{p}(t'))dt')

之前，由于σt\sigma_t是常数，括号内部的积分是有封闭形式的，有了封闭形式我们才有后面的重要性采样啊什么的。但是现在括号内部积分写不出封闭形式了，而且由于E[exp⁡(X)]≠exp⁡[E(X)]E[\exp(X)] \neq \exp[E(X)]，我们无法再用我们喜欢的蒙特卡洛积分来估计这个传输项。需要另外想办法。

如果介质在空间中的分布不均匀，我们可以认为这是介质的粒子在空间中不同位置的密度不同导致的。我们把这个新问题转化成我们解决过的问题，可以给介质里填充“空粒子”，给不均匀的介质里“掺水”得到均匀介质。

如果形式化地表达这个过程，我们首先把之前的微分方程写下来。

ddtL(p⃗(t),ω)=−σt(p⃗(t))L(p⃗(t),ω)+σs(p⃗(t))∫S2ρ(ω,ω′)L(p⃗(t),ω′)dω′\frac{\text{d}}{\text{d}t}L(\vec{p}(t),\omega)=-\sigma_t(\vec{p}(t))L(\vec{p}(t), \omega) + \sigma_s(\vec{p}(t))\int_{S^2}\rho(\omega,\omega')L(\vec{p}(t), \omega')\text{d}\omega'

边界条件就是

L(p⃗(thit),ω)=Lsurface(p⃗(thit))L(\vec{p}(t_{hit}),\omega)=L_{surface}(\vec{p}(t_{hit}))

向介质中“掺水”也就是在这里减一项加一项σnp⃗(t)\sigma_n{\vec{p}(t)}

ddtL(p⃗(t),ω)=−(σt(p⃗(t))+σn(p⃗(t)))L(p⃗(t),ω)+σn(p⃗(t))L(p⃗(t),ω)+σs(p⃗(t))∫S2ρ(ω,ω′)L(p⃗(t),ω′)dω′\frac{\text{d}}{\text{d}t}L(\vec{p}(t),\omega)=-(\sigma_t(\vec{p}(t))+\sigma_n(\vec{p}(t)))L(\vec{p}(t), \omega) +\sigma_n(\vec{p}(t))L(\vec{p}(t), \omega) + \sigma_s(\vec{p}(t))\int_{S^2}\rho(\omega,\omega')L(\vec{p}(t), \omega')\text{d}\omega'

σn\sigma_n就代表了填充的空粒子的密度，或者说空粒子“贡献”的系数，我们只需要让系数σt+σn\sigma_t+\sigma_n在当前行进方向上是常数就行。

我们如果要把介质“掺水”到一定密度，那么“掺水”到什么密度比较好呢？一个常用的选择是当前光线行进方向上系数的上界σm\sigma_m=sup⁡t∈[0,thit){σt(p⃗(t))}\sup\limits_{t\in[0,t_{hit})}\{\sigma_t(\vec{p}(t))\}。

这样我们积分一下就能算出来

L(p⃗(0),ω)=∫0thitτm(p⃗(0),p⃗(t′))(σn(p⃗(t′))L(p⃗(t′),ω)+σs(p⃗(t′))∫S2ρ(ω,ω′)L(p⃗(t),ω′)dω′)L(\vec{p}(0), \omega)=\int_0^{t_{hit}}\tau_m(\vec{p}(0),\vec{p}(t'))(\sigma_n(\vec{p}(t'))L(\vec{p}(t'),\omega)+\sigma_s(\vec{p}(t'))\int_{S^2}\rho(\omega,\omega')L(\vec{p}(t), \omega')\text{d}\omega')

+τm(p⃗(0),p⃗(thit))Lsurface(p⃗(thit))+\tau_m(\vec{p}(0),\vec{p}(t_{hit}))L_{surface}(\vec{p}(t_{hit}))

而τm(p⃗(0),p⃗(t))=exp⁡(−σmt)\tau_m(\vec{p}(0),\vec{p}(t))=\exp(-\sigma_m t)，这就是我们熟悉的东西。

这个式子很长，但是可以看出来没有积分在指数上了，传输项是我们熟悉的样子，就可以放心地设计蒙特卡洛积分算法。

  1. 首先按照指数分布采样一个tt

  2. 如果t>thitt>t_{hit}，那么我们自然而然地进行表面光照的计算。

  3. 如果t<thitt<t_{hit}，那么有两种可能的情况会发生。

第一种情况是我们估计散射的辐射度，这就对应于计算σs(p⃗(t′))∫S2ρ(ω,ω′)L(p⃗(t),ω′)dω′)\sigma_s(\vec{p}(t'))\int_{S^2}\rho(\omega,\omega')L(\vec{p}(t), \omega')\text{d}\omega')。

另一种情况是我们继续采样并估计L(p⃗(t′),ω)L(\vec{p}(t'),\omega)，这对应于计算σn(p⃗(t′))L(p⃗(t′),ω)\sigma_n(\vec{p}(t'))L(\vec{p}(t'),\omega)。

实际上，我们也可以这样来理解。两种情况就相当于光线行进到此处，要么撞到了介质粒子，发生了散射，要么撞到了空粒子，继续前进。我们需要分配这两个概率以使得方差较小（因为理论上来讲只要采样方法和样本概率能保持一致，蒙特卡洛积分就是无偏的）。一个比较好的方案是σtσm\frac{\sigma_t}{\sigma_m}和σnσm\frac{\sigma_n}{\sigma_m}。这很好理解，因为σt\sigma_t是介质粒子密度，σn\sigma_n是空粒子密度。相当于在所有粒子中有σnσm\frac{\sigma_n}{\sigma_m}的粒子是空粒子。

指导书这里提到，实际上更好的方案是以σsσs+σn\frac{\sigma_s}{\sigma_s+\sigma_n}的概率进行散射，但是考虑到σs+σn\sigma_s + \sigma_n可能会为00，为了简单我们不考虑这种情况。

由此，我们就得到了经典的delta tracking\text{delta tracking}或者Woodcock tracking\text{Woodcock tracking}。

2.7.2. delta tracking for next event estimation\text{2.7.2. delta tracking for next event estimation}

前面是重要性采样的方案。但是我们在做next event estimation\text{next event estimation}的时候，我们不需要重要性采样，我们需要的是配合我们前面的null scattering\text{null scattering}计算采样概率以及计算传输项。

τ(p⃗(0),p⃗(t))=exp⁡(−∫0tσt(p⃗(t′))dt′)\tau(\vec{p}(0),\vec{p}(t))=\exp(-\int_0^t\sigma_t(\vec{p}(t'))\text{d}t')

我们还是不喜欢这个东西，因此把它写成微分方程，并且尝试把σt\sigma_t换掉。

ddtτ=(−σm+σn(p⃗(t)))τ\frac{\text{d}}{\text{d}t}\tau=(-\sigma_m+\sigma_n(\vec{p}(t)))\tau

边界条件为τ(0)=1\tau(0)=1

很容易想到两边同时乘上exp⁡(−σmt)\exp(-\sigma_mt)，写出(实际上也就是Duhamel\text{Duhamel}原理)

τ(p⃗(t))=exp⁡(−σmt)+∫0texp⁡(σm(t′−t))σn(p⃗(t′))τ(p⃗(t′))dt′\tau(\vec{p}(t))=\exp(-\sigma_m t)+\int_0^t\exp(\sigma_m(t'-t))\sigma_n(\vec{p}(t'))\tau(\vec{p}(t'))\text{d}t'

注意左右两边都有待求的τ\tau，我们不能直接写出τ\tau的表达式。这个式子看起来很复杂，但它是能够使用蒙特卡洛积分处理的形式，之前的表达式因为积分式在指数上，不能准确估计，这里就没有这个问题了。

指导书里提到了一嘴，事实上仅仅使用含有−σt-\sigma_t的表达式进行上述过程，也是可以得到一个能用蒙特卡洛积分处理的形式，但是结果表明这种方法非常非常慢。

用来求解上述积分式的蒙特卡洛方法就是经典的ratio tracking\text{ratio tracking}方法。注意到，这个东西和渲染方程一样，都是Fredholm\text{Fredholm}积分方程，那么我们可以用类似于渲染方程的路径追踪的解法来解它。

具体而言，这个积分方程中的“核”是exp⁡(σm(t′−t))σn(p⃗(t′))\exp(\sigma_m(t'-t))\sigma_n(\vec{p}(t'))，我们要对其进行重要性采样，很自然的，我们就应该采样exp⁡(−σmt′)\exp(-\sigma_m t')，然后把积分式带入并估计采样的新点的τ\tau值，最后除以概率密度作为估计。不过我们还是有exp⁡(−σmt)\exp(-\sigma_m t)的可能性采到比tt更大的数，这个概率刚好就是前面的常数项，因此当发生这种情况时，贡献为11。这样我们就用蒙特卡洛方法解决了这个积分方程。

2.7.3. 多色介质处理和MIS\text{2.7.3. 多色介质处理和MIS}

我们考虑颜色的不同通道（大多数情况下为RGB\text{RGB}三通道）。很有可能对应三个通道的σ\sigma不一样，这样我们在做距离采样的时候就不知道用哪个。

一个解决方法是，首先等概率随机采样一个通道，然后使用该通道的σ\sigma进行距离的采样。这样做的话，概率密度就是每个通道下采样该距离的平均值。

和之前的采样方法相比，这里有一个很重要的微妙的区别。由于现在是用蒙特卡洛方法估计传输项，我们的样本不再是一个点，而是一系列的点。同时为了能在MIS\text{MIS}中准确计算权重，我们同时需要ratio tracking\text{ratio tracking}和delta tracking\text{delta tracking}采样的概率密度，这就要同时计算两种概率密度。

这里解释一下指导书中涉及这一部分的伪代码。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
if t < dt: # haven’t reached the surface
# sample from real/fake particle events
    real_prob = sigma_t / majorant
    if next(rng) < real_prob[channel]:
# hit a "real" particle
        scatter = true
        transmittance *= exp(-majorant * t) / max(majorant)
        trans_dir_pdf *= exp(-majorant * t) * majorant * real_prob / max(majorant)
# don’t need to account for trans_nee_pdf since we scatter
    break
    else:
# hit a "fake" particle
        transmittance *= exp(-majorant * t) * sigma_n / max(majorant)
        trans_dir_pdf *= exp(-majorant * t) * majorant * (1 - real_prob) / max(majorant)
        trans_nee_pdf *= exp(-majorant * t) * majorant / max(majorant)
else: # reach the surface
    transmittance *= exp(-majorant * dt)
    trans_dir_pdf *= exp(-majorant * dt)
    trans_nee_pdf *= exp(-majorant * dt)
    break

需要解释的一点是，这里所有除以max(majorant)的操作都是为了数值鲁棒性，因为我们算概率密度或者概率的时候很可能要乘上σ\sigma，而σm\sigma_m搞不好会比较大，不归一化的话几轮迭代下来就喜提数值爆炸了。对于概率和概率密度做这个事情当然会改变数值————不过不要紧，因为我们并不实际使用其数值，概率和概率密度的数值往往是作为齐次分式的分母出现的，我们只要把分子（例如transimittance）也除以max(majorant)就可以抵消了。

首先看外层的else，这个很容易理解，我们采样到了表面，那么无论是中间的传输项（用τm\tau_m计算的“掺水”介质的传输项，有显式解），还是使用delta tracking\text{delta tracking}方法或者ratio\text{ratio}采样的概率（注意不是概率密度），都是e−σmte^{-\sigma_mt}。

再看内层的if，这表示的是我们撞到介质粒子进行散射的情况。在delta tracking\text{delta tracking}的积分式中，积分项是exp⁡(σmt)\exp(\sigma_mt)，这就是transmittance。

我们使用delta tracking\text{delta tracking}采样到这个散射事件的概率密度则是exp⁡(σmt)σm×σtσm=exp⁡(σmt)σt\exp(\sigma_mt)\sigma_m\times \frac{\sigma_t}{\sigma_m}=\exp(\sigma_mt)\sigma_t。这很好理解，第一项表示采样当前距离的概率密度，第二项则表示采样到散射事件的概率，这自然就是trans_dir_pdf。NEE\text{NEE}使用的ratio tracking\text{ratio tracking}方法的采样概率不用计算，因为我们马上要进行散射，再算也没意义。

最后看内层的else。这表示我们采样到了Null Scattering\text{Null Scattering}事件。根据前面的式子，积分项即transmittance就是exp⁡(σmt)σn(p⃗(t),t)\exp(\sigma_mt)\sigma_n(\vec{p}(t),t)，此时我们可以继续前进估计直接辐射度。在该点采样到Null Scattering\text{Null Scattering}事件的概率是exp⁡(σmt)σm×σnσm=exp⁡(σmn)σn\exp(\sigma_mt)\sigma_m\times \frac{\sigma_n}{\sigma_m}=\exp(\sigma_mn)\sigma_n，而根据前文我们提到的ratio tracking\text{ratio tracking}使用的采样方法，ratio tracking\text{ratio tracking}采样到该点的概率是exp⁡(σmt)σm\exp(\sigma_mt)\sigma_m，因为那里不需要考虑Null Scattering\text{Null Scattering}所以这就是trans_nee_pdf的值。

最后需要注意的是，我们这里所有的概率和概率密度都是三个通道的向量。但是任何时候更新throughput\text{throughput}用的都是三个通道的算术平均值，这是因为按照我们对通道的采样方式（等概率采样），采样到任意通道完成上述过程的概率就是pR+pG+pB3\frac{p_R+p_G+p_B}{3}。当然，多重重要性采样的权重计算也是如此。

对于NEE\text{NEE}使用的delta tracking\text{delta tracking}，计算其实是差不多的。只不过这里如果纯粹看指导书上的内容可能会和代码不太能完全对应，稍微解释一下。

如果看我们前面在ratio tracking\text{ratio tracking}里给出的积分方程，会发现它有exp⁡(σm(t′−t))\exp(\sigma_m(t'-t))项，但是给出的伪代码里不是这样的。

这其实涉及到一个定义不同的问题，可以按照我们在2.3.1\text{2.3.1}里做的那样，“反过来”，得到的就是伪代码里的东西。

这一块儿仅用口舌和仅用代码都很难说清楚，更加推荐读一读指导书的第一篇参考论文以及Tzumao Li本人在课上提到的A Null-Scattering Path Integral Formulation of Light Transport。后一篇论文里扯了很多很形式化的东西。。。Tzumao Li说不看这些推导你也能学会，谁TM关心路径空间的测度啊

2.7.4. Coding\text{2.7.4. Coding}

作为最后的绿宝石水花任务，这次的代码一如既往地细节多。个人觉得主要是这方面的实现公说公有理婆说婆有理，五花八门，不好说什么样的写法才是最标准的。。。比如指导书的讲解并没有特别说明index-matching surface\text{index-matching surface}应该如何处理。。。按照自己的理解写出来了一份（至少看起来是）无偏的代码，至于方差如何，实在很难衡量了。希望自己有机会遇到更强的大佬并且彻底把这些问题研究清楚吧。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
Spectrum nee(const Ray& ray, int depth, int cur_medium_id, const Scene& scene, pcg32_state& rng)
{
    int light_id = sample_light(scene, next_pcg32_real<Real>(rng));
    const PointAndNormal& light_PN = sample_point_on_light(scene.lights[light_id], ray.org,
                                                           Vector2(next_pcg32_real<Real>(rng), next_pcg32_real<Real>(rng)), next_pcg32_real<Real>(rng), scene);
    Ray nee_ray;
    nee_ray.org = ray.org;
    nee_ray.dir = normalize(light_PN.position - ray.org);
    nee_ray.tnear = get_intersection_epsilon(scene);
    nee_ray.tfar = infinity<Real>();
    const PhaseFunction& phaseFunction = get_phase_function(scene.media[cur_medium_id]);
    Spectrum nee_throughput = eval(phaseFunction, -ray.dir, nee_ray.dir);
    int nee_dep = 0;
    Real pdf_light = pdf_point_on_light(scene.lights[light_id], light_PN, ray.org, scene) * light_pmf(scene, light_id);
    Real G = Geometry(light_PN, ray.org);
    Real pdf_delta = 1, pdf_ratio = 1.0;
    while(true)//we do a next event estimation
    {
        if(scene.options.max_depth >= 0 && depth + nee_dep >= scene.options.max_depth)
            return make_zero_spectrum();
        std::optional<PathVertex> opt_pathVertex = intersect(scene, nee_ray);
        if(opt_pathVertex)
        {
            PathVertex pathVertex = opt_pathVertex.value();
            if(cur_medium_id >= 0)
            {
                Real t_hit = distance(pathVertex.position, nee_ray.org);
                int channel = floor(next_pcg32_real<Real>(rng) * 3);
                Spectrum sigma_m = get_majorant(scene.media[cur_medium_id], nee_ray);
                int null_collision_cnt = 0;
                Spectrum T = make_const_spectrum(1), trans_pdf_delta = make_const_spectrum(1),
                            trans_pdf_ratio = make_const_spectrum(1);
                while(true)
                {
                    Real u = next_pcg32_real<Real>(rng);
                    if(sigma_m[channel] <= 0.0) break;
                    Real t = -log(1.0 - u) / sigma_m[channel];
                    if (t >= t_hit)
                    {
                        if(pathVertex.material_id == -1)//it is an index-matching surface, we go through the interface
                        {
                            nee_dep++;
                            T *= exp(-sigma_m * t_hit);
                            trans_pdf_delta *= exp(-sigma_m * t_hit);
                            trans_pdf_ratio *= exp(-sigma_m * t_hit);
                        }
                        break;
                    }
                    else// in ratio tracking, we just perform null scattering event
                    {
                        Spectrum sigma_s = get_sigma_s(scene.media[cur_medium_id], nee_ray.org);
                        Spectrum sigma_t = get_sigma_a(scene.media[cur_medium_id], nee_ray.org) + sigma_s;
                        nee_ray.org += t * nee_ray.dir;
                        Spectrum sigma_n = sigma_m - sigma_t;
                        T *= exp(-sigma_m * t) * sigma_n / Max(sigma_m);
                        trans_pdf_delta *= exp(-sigma_m * t) * sigma_n / Max(sigma_m);
                        trans_pdf_ratio *= exp(-sigma_m * t) * sigma_m / Max(sigma_m);
                        null_collision_cnt++;
                        if(null_collision_cnt >= scene.options.max_null_collisions) break;
                        if(Max(T) <= 0.0) return make_zero_spectrum();
                    }
                    t_hit -= t;
                }
                nee_throughput *= (T / avg(trans_pdf_ratio));
                pdf_ratio *= avg(trans_pdf_ratio);
                pdf_delta *= avg(trans_pdf_delta);
            }
            nee_ray.org = pathVertex.position;
            if(pathVertex.material_id >= 0 && get_area_light_id(scene.shapes[pathVertex.shape_id]) != light_id)
                return make_zero_spectrum();
            if(distance(nee_ray.org, light_PN.position) < get_intersection_epsilon(scene)) // we reach the light source, and we calc the contrib
            {
                Real pdf_phase = pdf_sample_phase(phaseFunction, -nee_ray.dir, ray.dir)
                                 * G * pdf_delta;
                Real pdf_nee = pdf_light * pdf_ratio;
                Spectrum L_nee = nee_throughput * emission(pathVertex, -nee_ray.dir, scene) * G / pdf_light;
                Real w = pdf_nee * pdf_nee / (pdf_nee * pdf_nee + pdf_phase * pdf_phase);
                return L_nee * w;
            }
            update_medium(pathVertex, nee_ray, cur_medium_id);
            nee_dep++;
        }
        else return make_zero_spectrum();
    }
}
Spectrum nee_bsdf(const Ray& ray, int depth, const PathVertex& path_vertex, int cur_medium_id,
                                    const Scene& scene, pcg32_state& rng)
{
    int light_id = sample_light(scene, next_pcg32_real<Real>(rng));
    const PointAndNormal& light_PN = sample_point_on_light(scene.lights[light_id], ray.org,
                                                           Vector2(next_pcg32_real<Real>(rng), next_pcg32_real<Real>(rng)), next_pcg32_real<Real>(rng), scene);
    Ray nee_ray;
    nee_ray.org = ray.org;
    nee_ray.dir = normalize(light_PN.position - ray.org);
    nee_ray.tnear = get_intersection_epsilon(scene);
    nee_ray.tfar = infinity<Real>();
    const Material& mat = scene.materials[path_vertex.material_id];
    Spectrum nee_throughput = eval(mat, -ray.dir, nee_ray.dir, path_vertex, scene.texture_pool);
    int nee_dep = 0;
    Real pdf_light = pdf_point_on_light(scene.lights[light_id], light_PN, ray.org, scene)
                        * light_pmf(scene, light_id);
    Real G = Geometry(light_PN, ray.org);
    Real pdf_delta = 1, pdf_ratio = 1;
    update_medium(path_vertex, nee_ray, cur_medium_id);
    while(true)//we do a next event estimation
    {
        if(scene.options.max_depth >= 0 && depth + nee_dep >= scene.options.max_depth)
            return make_zero_spectrum();
        std::optional<PathVertex> opt_pathVertex = intersect(scene, nee_ray);
        if(opt_pathVertex)
        {
            PathVertex pathVertex = opt_pathVertex.value();
            if(cur_medium_id >= 0)
            {
                Real t_hit = distance(pathVertex.position, nee_ray.org);
                int channel = floor(next_pcg32_real<Real>(rng) * 3);
                Spectrum sigma_m = get_majorant(scene.media[cur_medium_id], nee_ray);
                int null_collision_cnt = 0;
                Spectrum T = make_const_spectrum(1), trans_pdf_delta = make_const_spectrum(1),
                                trans_pdf_ratio = make_const_spectrum(1);
                while(true)
                {
                    Real u = next_pcg32_real<Real>(rng);
                    if(sigma_m[channel] <= 0.0) break;
                    Real t = -log(1.0 - u) / sigma_m[channel];
                    if (t >= t_hit)
                    {
                        if(pathVertex.material_id == -1)//it is a index-matching surface, we go through the interface
                        {
                            depth++;
                            T *= exp(-sigma_m * t_hit);
                            trans_pdf_delta *= exp(-sigma_m * t_hit);
                            trans_pdf_ratio *= exp(-sigma_m * t_hit);
                            nee_ray.org = pathVertex.position;
                        }
                        break;
                    }
                    else// in ratio tracking, we just perform null scattering event
                    {
                        Spectrum sigma_s = get_sigma_s(scene.media[cur_medium_id], nee_ray.org);
                        Spectrum sigma_t = get_sigma_a(scene.media[cur_medium_id], nee_ray.org) + sigma_s;
                        nee_ray.org += t * nee_ray.dir;
                        Spectrum sigma_n = sigma_m - sigma_t;
                        T *= exp(-sigma_m * t) * sigma_n / Max(sigma_m);
                        trans_pdf_delta *= exp(-sigma_m * t) * sigma_n / Max(sigma_m);
                        trans_pdf_ratio *= exp(-sigma_m * t) * sigma_m / Max(sigma_m);
                        null_collision_cnt++;
                        if(null_collision_cnt >= scene.options.max_null_collisions) break;
                        if(Max(T) <= 0.0) return make_zero_spectrum();
                    }
                    t_hit -= t;
                }
                nee_throughput *= (T / avg(trans_pdf_ratio));
                pdf_ratio *= avg(trans_pdf_ratio);
                pdf_delta *= avg(trans_pdf_delta);
            }
            nee_ray.org = pathVertex.position;
            if(pathVertex.material_id >= 0 && get_area_light_id(scene.shapes[pathVertex.shape_id]) != light_id)
                return make_zero_spectrum();
            if(distance(nee_ray.org, light_PN.position) < get_intersection_epsilon(scene)) // we reach the light source, and we calc the contrib
            {
                Real pdf_bsdf = pdf_sample_bsdf(mat, -ray.dir, nee_ray.dir, path_vertex,
                                                     scene.texture_pool) * G * pdf_delta;
                Real pdf_nee = pdf_light * pdf_ratio;
                Spectrum L_nee = nee_throughput * emission(pathVertex, -nee_ray.dir, scene) * G / pdf_light;
                Real w = pdf_nee * pdf_nee / (pdf_nee * pdf_nee + pdf_bsdf * pdf_bsdf);
                return L_nee * w;
            }
            update_medium(pathVertex, nee_ray, cur_medium_id);
            nee_dep++;
        }
        else return make_zero_spectrum();
    }
}

// The final volumetric renderer: 
// multiple chromatic heterogeneous volumes with multiple scattering
// with MIS between next event estimation and phase function sampling
// with surface lighting
Spectrum vol_path_tracing(const Scene &scene,
                          int x, int y, /* pixel coordinates */
                          pcg32_state &rng) {
    // Homework 2: implement this!
    Spectrum cur_path_throughput(1.0, 1.0, 1.0);
    Real X = (static_cast<Real>(x) + next_pcg32_real<Real>(rng)) / scene.camera.width;
    Real Y = (static_cast<Real>(y) + next_pcg32_real<Real>(rng))/ scene.camera.height;
    Ray ray = sample_primary(scene.camera, Vector2(X, Y));
    ray.tnear = get_intersection_epsilon(scene);
    int cur_medium_id = scene.camera.medium_id, depth = 0;
    Spectrum L(0.0, 0.0, 0.0);
    Real pdf_ratio = 1, pdf_delta = 1;// p_march for nee
    PathVertex pv_nee_cache;
    Real pdf_dir_cache = 0;
    bool never_scatter_or_reflect = true;
    while(true)
    {
        if(scene.options.max_depth >= 0 && depth >= scene.options.max_depth)break;
        if(scene.options.rr_depth >= 0 && depth >= scene.options.rr_depth)//first do a Russian Roulette
        {
            Real rnd_rr = next_pcg32_real<Real>(rng);
            Real P_RR = 0.95;
            if(rnd_rr > P_RR) break;
            else cur_path_throughput /= P_RR;
        }
        std::optional<PathVertex> opt_path_vertex = intersect(scene, ray);
        bool reach_light = false;
        bool scatter = false;
        bool surface_lighting = false;
        Spectrum sigma_s, sigma_t;
        Real t_hit;
        PathVertex path_vertex;
        PointAndNormal surface_PN;
        int light_id, channel;
        if(!opt_path_vertex) t_hit = infinity<Real>();
        else
        {
            path_vertex = opt_path_vertex.value();
            surface_PN = {path_vertex.position,  path_vertex.geometry_normal};
            t_hit = distance(ray.org, path_vertex.position);
        }
        if(cur_medium_id >= 0)// if we are now in a medium
        {
            channel = floor(next_pcg32_real<Real>(rng) * 3);
            Spectrum sigma_m = get_majorant(scene.media[cur_medium_id], ray);
            int null_collision_cnt = 0;
            Spectrum T = make_const_spectrum(1), trans_pdf_delta = make_const_spectrum(1),
                        trans_pdf_ratio = make_const_spectrum(1);
            while(true)// perform free-flight sampling
            {
                Real u = next_pcg32_real<Real>(rng);
                if(sigma_m[channel] <= 0.0) break;
                Real t = -log(1.0 - u) / sigma_m[channel];
                if (t >= t_hit)
                {
                    trans_pdf_delta *= exp(-t_hit * sigma_m);
                    trans_pdf_ratio *= exp(-t_hit * sigma_m);
                    T *= exp(-t_hit * sigma_m);
                    ray.org = path_vertex.position;
                    if(path_vertex.material_id == -1)//it is a index-matching surface, we go through the interface
                    {
                        depth++;
                        update_medium(path_vertex, ray, cur_medium_id);
                    }
                    else
                    {
                        light_id = get_area_light_id(scene.shapes[path_vertex.shape_id]);
                        if(light_id != -1) reach_light = true;
                        else surface_lighting = true;
                    }
                    break;
                }
                else// we sample to decide whether it is a null collision or scattering event
                {
                    ray.org += t * ray.dir;
                    sigma_s = get_sigma_s(scene.media[cur_medium_id], ray.org);
                    sigma_t = get_sigma_a(scene.media[cur_medium_id], ray.org) + sigma_s;
                    if(next_pcg32_real<Real>(rng) < sigma_t[channel] / sigma_m[channel])
                    {
                        scatter = true;
                        T *= exp(-sigma_m * t) / Max(sigma_m);
                        trans_pdf_delta *= exp(-sigma_m * t) * sigma_t / Max(sigma_m);
                        break;
                    }
                    else
                    {
                        Spectrum sigma_n = sigma_m - sigma_t;
                        T *= exp(-sigma_m * t) * sigma_n / Max(sigma_m);
                        trans_pdf_delta *= exp(-sigma_m * t) * sigma_n / Max(sigma_m);
                        trans_pdf_ratio *= exp(-sigma_m * t) * sigma_m / Max(sigma_m);
                        if(Max(T) <= 0.0) break;
                        null_collision_cnt++;
                        if(null_collision_cnt >= scene.options.max_null_collisions) break;
                    }
                }
                t_hit -= t;
            }
            cur_path_throughput *= (T / avg(trans_pdf_delta));
            pdf_delta *= avg(trans_pdf_delta);
            pdf_ratio *= avg(trans_pdf_ratio);
        }
        else
        {
            if(!opt_path_vertex) break;
            light_id = get_area_light_id(scene.shapes[path_vertex.shape_id]);
            ray.org = path_vertex.position;
            if(light_id >= 0) reach_light = true;
            else if(path_vertex.material_id >= 0) surface_lighting = true;
            else update_medium(path_vertex, ray, cur_medium_id);
        }
        if(scatter)
        {
            cur_path_throughput *= sigma_s;
            Spectrum L_nee = nee(ray, depth, cur_medium_id, scene, rng);
            L += cur_path_throughput * L_nee;
            pv_nee_cache = path_vertex;
            never_scatter_or_reflect = false;
            PhaseFunction phaseFunction = get_phase_function(scene.media[cur_medium_id]);
            Vector2 rnd_uv(next_pcg32_real<Real>(rng), next_pcg32_real<Real>(rng));
            std::optional<Vector3> opt_dir = sample_phase_function(phaseFunction, -ray.dir, rnd_uv);
            if (opt_dir)
            {
                Vector3 dir = opt_dir.value();
                Real phase_pdf = pdf_sample_phase(phaseFunction, -ray.dir, dir);
                pdf_dir_cache = phase_pdf;
                pdf_ratio = pdf_delta  = 1.0;
                cur_path_throughput *= eval(phaseFunction, -ray.dir, dir) / phase_pdf;
                ray.dir = dir;
            }
            else break;
        }
        if(reach_light)
        {
            if(never_scatter_or_reflect) L += cur_path_throughput * emission(path_vertex, -ray.dir, scene);
            else
            {
                Real G = Geometry(surface_PN, pv_nee_cache.position);
                Real pdf_dir = pdf_delta * G * pdf_dir_cache;
                Real pdf_nee = pdf_point_on_light(scene.lights[light_id], surface_PN, pv_nee_cache.position,
                                           scene) * light_pmf(scene, light_id) * pdf_ratio;
                Real w = pdf_dir * pdf_dir / (pdf_dir * pdf_dir + pdf_nee * pdf_nee);
                L += cur_path_throughput * emission(path_vertex, -ray.dir, scene) * w;
            }
            break;
        }
        if(surface_lighting)
        {
            never_scatter_or_reflect = false;
            Spectrum L_nee = nee_bsdf(ray, depth, path_vertex, cur_medium_id, scene, rng);
            L += cur_path_throughput * L_nee;
            pv_nee_cache = path_vertex;
            const Material& mat = scene.materials[path_vertex.material_id];
            Vector2 rnd_uv(next_pcg32_real<Real>(rng), next_pcg32_real<Real>(rng));
            std::optional<BSDFSampleRecord> opt_bRec = sample_bsdf(mat, -ray.dir, path_vertex, scene.texture_pool,
                                                                   rnd_uv, next_pcg32_real<Real>(rng));
            if (opt_bRec)
            {
                BSDFSampleRecord bRec = opt_bRec.value();
                Real bsdf_pdf = pdf_sample_bsdf(mat, -ray.dir, bRec.dir_out, path_vertex, scene.texture_pool);
                if(bsdf_pdf == 0.0) break;
                cur_path_throughput *= eval(mat, -ray.dir, bRec.dir_out, path_vertex, scene.texture_pool) / bsdf_pdf;
                pdf_dir_cache = bsdf_pdf;
                pdf_delta = pdf_ratio = 1.0;
                ray.dir = bRec.dir_out;
                update_medium(path_vertex, ray, cur_medium_id);
            }
            else break;
        }
        depth++;
    }
    return L;
}

这，就是，本章的最后的代码了！

2.7.5. 加分项\text{2.7.5. 加分项}

研究一下表达介质的文件格式，试着生成属于自己的介质并且渲染它！

对我而言这就是下一个阶段的事情了，本次作业使用的介质文件是Wenzel Jakob\text{Wenzel Jakob}怎么哪哪儿都有这个佬啊写的烟雾模拟生成的这个佬怎么什么都会啊。下一步学习物理模拟，来生成这些东西吧。

接下来还有些不怎么容易的思考题。。。。这些问题值得长期思考，经过这么长时间体积渲染的折磨我也想不出个啥了能回答的回答一下，不能回答的一并坑在这里吧（

UPD: \text{UPD: }最后啊，一年之后，写了个GPU\text{GPU}加速的烟雾模拟，效果比Wenzel Jakob\text{Wenzel Jakob}的难看多了，但是再烂也是自己写的。

2.7.6. 思考题\text{2.7.6. 思考题}

怎样的密度分布使得上面的渲染算法效率最低/高？\text{怎样的密度分布使得上面的渲染算法效率最低/高？}

个人觉得当密度分布的变化频率大并且波动很大的时候会导致效率很低。。。但是只是感觉，说也不好说啊。。。

你能想到针对恶劣环境的改进办法吗？\text{你能想到针对恶劣环境的改进办法吗？}

想不到，谢谢。

如何让Null Scattering方法适用于发光介质？\text{如何让Null Scattering方法适用于发光介质？}

猜一下，加入发光粒子的选项？

为什么无偏的渲染介质的算法如此重要？\text{为什么无偏的渲染介质的算法如此重要？}

啊这。。。搞学术的不该追求无偏吗。。。

用Ravi Ramamoorthi\text{Ravi Ramamoorthi}的话来说就是

It is easy to render a beautiful image, but it is hard to render a correct one.

设计有偏但是更快的算法是否明智？\text{设计有偏但是更快的算法是否明智？}

个人猜测影视工业界应该会比较喜欢稍微有偏但是更快的算法吧，毕竟电影看着让人舒服就行，如果为了追求一个无偏多花了一大堆电费和时间最后观众也看不出来就很吃亏嘛（

3. 总结\text{3. 总结}

终于结束了，这个过程中学到了很多，渲染算法的巨量细节耗费的debug\text{debug}时间极其长，还遗留了太多问题要解决。当然，也有更多好玩有趣的话题值得探索。

2023\text{2023}年春天的CSE 272\text{CSE 272}课程正在进行中。今年的Final Project\text{Final Project}是从给定的若干项目中（例如BDPT,BSSRDF,BCSDF for Hair...\text{BDPT,BSSRDF,BCSDF for Hair...}）选择一个并且复现。甚至也有一些研究项目可供选择。当然精力不够了这些就慢慢学吧，先不急着实现了。只能说实在羡慕国外能有如此充实的渲染课程和如此高质量的作业虽然大部分时间几乎都是在和体积渲染作斗争。当然也希望小破校的渲染科研水平以后越来越好，能够开出这样优秀的课程吧以小破校目前选修课的深度我看是没可能了。

GAMES303\text{GAMES303}之前，Rendering\text{Rendering}就到此为止吧。下一个就是你了，承太郎simulation\text{simulation}！