本节数学符号较多，高能预警

引言

上一节我们完成了全局光照路径追踪算法的实现，并制作了三种最基础的材质实例。这一节我们主要实现抽象光源以及对光源的直接采样。

Part0. bug修正

在上一节的末尾，我对若干场景进行了测试，整体效果较好，但是有一些明显的bug需要修正。首先是在Raygen shader中，若光线弹射次数超过了预设值，应该直接break循环，而不需要给radiance置0。如果给radiance置0，那么路径中自发光物体的光照贡献也会被一并清除，自发光就失去意义了。

for (int bounces = 0; ; ++bounces)
{
     if (bounces >= optixLaunchParams.maxBounce) {
     //radiance = 0.0f;
     break;
}

另一个问题是导入模型的问题，需要对knownvertices调整一下位置，否则在多材质的场景中可能会出现异常的三角面：

for (int materialID : materialIDs) {
     if (materialID == -1)
          continue;
     std::map<tinyobj::index_t, int> knownVertices;
...

Part1. 直接光采样理论

之前的场景基本上都是开放或半开放的场景。现在我们考虑一下全封闭的场景，全封闭意味着外面的环境贴图无法提供任何光照，那么此时我们就只能以封闭场景中物体的自发光作为光源。然而当我们创建了一个自发光强度极强的球体放在封闭场景中央，整个场景依旧漆黑无比：

究其原因是，我们场景中大部分物体是漫反射，其散射方向完全随机，这就导致K次采样中只有寥寥无几的几次（甚至0次）击中到自发光上。换言之，正常随机采样很有可能每次都求交不到有自发光的物体上，而且光源越小出现这种悲剧的概率越大，这样最终场景就会非常昏暗。

这里有一个解决思路：将直接光（光源）和间接光（其他物体弹射的光源）分开采样。直接光采直接光的，非直接光采非直接光的，最后将两者加在一起就成了最终的结果。这样就可以保证每次采样不会错过直接光的贡献。

至于如何分开采样，我们依旧回到渲染方程：

$$L(x→out)=Emission + \int_\omega^\omega f_r(x,in→out)L(x←in)cos(N,in)d\omega$$

根据蒙特卡洛积分原理，可以将积分写成采样/pdf的形式：

$$L(x→out)=Emission + \frac{1}{K}\sum \frac{f_r(x,in→out)L(x←in)cos(N,in)}{pdf}$$

这里的L(x←in)指的是从外界发射进来的所有光线的采样，那么我们就可以将其拆解为直接光和间接光部分：

$$L(x→out)=E + \frac{1}{K}\sum (\frac{f_r(x,in→out)L(x←光源)cos(N,in)}{pdf} + \frac{f_r(x,in→out)L(x←表面)cos(N,in)}{pdf})$$

由于路径追踪算法中，当光线击中某个表面后，不需要反复采样取平均（否则会出现光线数量指数级爆炸），因此可以将求和取平均的符号去掉，那么此时计算公式就变成了：

$$L(x→out)=E + (\frac{f_r(x,in→out)L(x←直接光)cos(N,in)}{pdf} + \frac{f_r(x,in→out)L(x←间接光)cos(N,in)}{pdf})$$

这里埋藏了一个坑，那就是：两项的pdf是否是同一个pdf？按照咱们的推导来看，pdf确实是相同的，但是我们实际采样的时候，光源会有一套采样方案（即根据光源外表形状决定的采样方案），而非直接光的采样又会用另外一套方案，因此此时两者的pdf就是不同的了。

现在看这个公式有点太臃肿了，我们令：

$$F_{C_1} = \frac{f_C * cos(N,直接光)}{pdf}，F_{C_2} = \frac{f_C * cos(N,间接光)}{pdf}$$

则有：

$$L_C=E_C + F_{C_1} * L_{直接光} + F_{C_2}*L_{间接光}$$

所以我们的思路就变成了：当求交到一个物体表面后，先去场景中找直接光源求交，获得直接光照后*bsdf*cosine/pdf_光源；再去场景中找物体表面求交，获得间接光照后*bsdf*cosine/pdf_bsdf；将两者加在一起，就是本次采样得到的总能量了。（相当于每次求交后要采样两次，一次直接光，一次间接光）

对于递归型path-tracing程序而言，算法其实很明确，这在games101的ray-tracing第四讲有展示过：

但是对于我们的迭代型path-tracing就又不一样了。迭代程序的算法该是什么样的呢？我们依然思考一个例子，光照D发出光线打中C，C反弹到B，B反弹到A，A反弹进相机。如果我们考虑直接光照采样，那么上述事件就变成了：光照D发出光线打中C，C反弹到B的同时D也发出光线打中B，B反弹到A的同时D也发出光线打中A，A反弹进相机。

好，我们从D到A依次分析出射光线的能量。首先，D射向C的光线就是E_D，C射向B的光线就是：

$$L_C出射 = E_C + F_{C_1}*E_{(C←D)}$$

然后我们考虑了直接光采样，因此，B所接收的光线不止来源于C，还来源于D，因此：

$$L_B出射 = E_B + F_{B_1} * E_{(B←D)} + F_{B_2} * (E_C + F_{C_1}*E_{(C←D)})$$

$$L_A出射 = E_A + F_{A_1} * E_{(A←D)} + F_{A_2} * (E_B + F_{B_1} * E_{(B←D)} + F_{B_2} * (E_C + F_{C_1}*E_{(C←D)}))$$

$$L_A出射 = E_A + F_{A_1} * E_{(A←D)} + F_{A_2} * E_B + F_{A_2} * F_{B_1} * E_{(B←D)} + F_{A_2} * F_{B_2} * E_C + F_{A_2} * F_{B_2} * F_{C_1} * E_{(C←D)}$$

如果还是看不出规律，我就再做一下高亮处理：

$$L_A出射 = 1 * E_A + 1 * F_{A_1} * E_{(A←D)} + F_{A_2} * E_B + F_{A_2} * F_{B_1} * E_{(B←D)} + F_{A_2} * F_{B_2} * E_C + F_{A_2} * F_{B_2} * F_{C_1} * E_{(C←D)}$$

我们从迭代的角度来分析一下这个过程。初始化radiance依旧是0，accum依旧是1（accum就是红色部分）；第一次我们求交到了A，那么首先radiance+=E_A * accum；然后我们去采样直接光，得到直接光后radiance+=F_A1 * 直接光的亮度贡献 * accum；然后我们采样下一次弹射方向，根据方向我们能计算出F_A2的值，给accum*=F_A2；

接下来求交到了B点，首先B点可能有自发光，radiance+=E_B * accum；然后采样直接光，得到直接光后radiance+=F_B1*直接光的亮度贡献*accum；采样下一次弹射方向，根据方向计算出F_B2的值，给accum*=F_B2，依此类推……

所以具有直接光采样的光路迭代伪代码就是：

radiance = 0
accum = 1
wo = the first direction of the camera ray
o = the camera position
for i in range(0,depth,1)
    intersect from a ray whose direction is 'wo' and origin is 'o'
    get intersection, the intersection point call 'p'
    radiance += E_p * accum
    sample the direct light, we assum the vector from p to light is wi1
    radiance += (bsdf(wo,wi1) * cosine / pdf_light) * L_direct * accum
    sample a indirect direction, we assum this direction is wi2
    accum *= (bsdf(wo,wi2) * cosine / pdf_pq)
    wo = wi2
    o = p

同时注意一点：直接光采样的时候，需要判断是否有物体遮挡在直接光和当前交点之间，因此一般会制作一个阴影射线来检查可见性。这个是后话了。

如此，我们便完成了直接光的采样，这样在path-tracing的过程中，每一次弹射都不会错过光源的影响，那么我们就解决了“随机采样大概率碰不到光源导致场景偏暗”的问题了。（这部分介绍的内容相对比较晦涩，需要对上面的公式推导有一定掌握，然后才能对算法的过程进行理解。）

当然细心的同学也会发现一个问题：我们本质上是将完全随机的采样空间劈成了两部分，一部分采样空间完全指向直接光，而剩下的采样空间自然是完全不能指向直接光的。然而我们在采样“剩下那部分空间”的时候依然用的是完全随机弹射，这就意味着它仍然可能弹射到光源上，同时说明我们劈成的两部分截然不同的采样空间发生了重叠，这必然是不能被允许的事情。

那么怎么办呢？这时候就会提出一个解决办法：将直接光源抽象化。

什么是抽象化？就是尽管它有体积，但是它无法被正常求交得到，它只是存在于场景某处的一个没有外形的发光体。如此设计，我们依然可以知道直接光的位置、光强等信息从而完成直接光采样，同时又能避免非直接光采样时不小心采到光源。当然，这就意味着我们要定义一种新的结构体来描述光源。

Part.2 矩形光源定义

光源这部分将是一个庞大的工程，因为我们考虑的事情非常多：场景如果有多个光源如何采样？不同形状类型的光源如何采样？对应的pdf如何计算？如何检测光源的可见性？所以这一部分要花不少时间来推导。不过作为工程嘛，都是从简入繁，这里我们先提出一个最最简单的光源：矩形光源。

有人会说其实无限小点光源才是最简单的光源，然而在光追算法中无限小点光源是个很不友好的选择，因为无限小意味着对其采样的pdf无法计算，那么就无法正确拟合蒙特卡洛积分的值，所以选取矩形光源作为最基础的光源。

首先我把原来的material_def改名成了message_def，用来定义材质、光源等各种基础信息。在该文件中我们做如下补充：

enum light_kind
{
    REC_L, POINT_L
};

struct light_mes {
    light_kind lgt_kind;
    vec3f position = 0.0f;
    vec3f u = vec3f(1, 0, 0);
    vec3f v = vec3f(0, 1, 0);
    float radiance = 0.0f;
    float emitter = 200.0f;
};

REC_L就是我们要定义的矩形光源，后面的POINT_L就是有体积的点光源，第一步我们只创建矩形光源；

同时我们在launchParams.h中定义light_mes light，用来将光源传入gpu空间；

然后我们将SampleRenderer构造函数的传入参数种添加一个light_mes，将场景的光源信息绑定进来。在main函数中，我们定义一个矩形光源并输入数据，将light传入：

light_mes light = { REC_L,vec3f(0.0f),vec3f(100,0,0),vec3f(0,100,0),0.0f,200.0f };
const float worldScale = length(model->bounds.span());

SampleWindow* window = new SampleWindow("PL Tracer", model, camera, light, worldScale);

上面这部分内容就是注册光照用的，不多赘述。此时，我们在shader中已经可以成功获取到场景中的光源信息了。

在光路迭代的过程中，我们首先要对直接光进行采样，那么第一件事就是：如何采样？这里我们必须给出一个采样方案（即给定一个pdf分布），如此我们才能根据pdf发射射线，并在最终结果中除以pdf。为了简化问题，我们限定场景中只有一个矩形光源，那么我们就要考虑如何对矩形光源采样。

其实采样方案非常简单，那就是：对矩形面均匀采样。因为矩形光源上处处光强都相等，没有哪个区域和其他区域有显著不同，因此我们直接均匀对对矩形上的每个区域进行采样。

假设矩形面积是S，所以我们的pdf就是1/S了吗？这里是一个很大的陷阱，我们再来看一下渲染方程的直接光部分：

$$L_{直接光}=\int_\omega^\omega f_r(x,in→out)L(x←in)cos(N,in)d\omega$$

注意：最后一项是dw，也就是我们的渲染方程是对立体角积分，而不是对光源的面积积分！这意味着我们采样也必须是对dw采样而非对光源上的dS采样（换句话说，你对dS均匀采样，转化成对dw的采样肯定是不均匀的，所以对dw采样的pdf没法算了）。既然如此，我们不妨将dw转化成dS的形式，这样采样对象也就自然转化到光源上的dS了。

那么如何找dw和dS间的关系呢？首先我们先明晰立体角的定义：“以观测点为球心，构造一个单位球面；任意物体投影到该单位球面上的投影面积，即为该物体相对于该观测点的立体角。”换句话说就是，我们在观测点创建一个半径为1的球体，dw就是光源上的dS在球体上的投影面积。

既然是计算投影，首先我们要乘以一个cos(L,N光源)计算出矩形光源面 正对着单位球体的面积部分；然后我们可以根据相似三角形的性质，得到：

$$\frac{投影面积}{dS}=\frac{1^2}{distance^2}$$

根据立体角定义，投影面积就是dw；同时dS是向光线方向的投影面积，因此：

$$\frac{d\omega}{dS*cos\theta}=\frac{1}{distance^2}$$

$$d\omega=\frac{cos\theta}{distance^2}*dS$$

假设出射点是O，在光源上的采样点是P，那么公式就变成了：

$$d\omega=\frac{cos(N_{光源},PO)}{(length(PO))^2}*dS$$

那我们代入到渲染方程的直接光部分，就是：

$$L_{直接光}=\int_\omega^\omega f_r(x,in→out)L(x←in)cos(N,in)\frac{cos(N_{光源},PO)}{(length(PO))^2}*dS$$

这就是大名鼎鼎的渲染方程-区域光公式。这里我们的积分对象变成了dS，此时我们就可以对光源上的面进行均匀采样了。这意味着当我们采样直接光的时候，在矩形光源上均匀采样，得到的光强先要乘一个dot(N光源,PO)，然后再除以距离平方，得到的值才是真正采样到的光强。现在回到之前那个抽象的采样式子：

$$L_C=E_C + F_{C_1} * L_{直接光} + F_{C_2}*L_{间接光}$$

我们只关注中间直接光采样的部分:

$$F_{C_1} * L_{矩形光} = \frac{f_C * cos(N_{平面},OP)}{pdf} * \frac{cos(N_{光源},PO)}{(length(PO))^2}*L_{矩形光强}$$

这个就和闫老师的课程上贴的那句代码吻合了：

这个式子就指引了我们，在shader中要如何做直接光采样。首先我们发射一根阴影射线，查找交点距离是不是比光源远，如果比光源远说明光源可见；在可见的前提下，在矩形光源上均匀采样，得到一个点P，然后计算该点的光强贡献：光源光强*cos(N光源，PO)/(|PO|^2)，这里得到的是光强贡献，接下来要乘对应的bsdf值、乘cos(N平面，OP)，最后除以pdf也就是除以(1/S)。