采样与误差

Participating media的采样

根据我们在上一篇软光栅下体积云理论与实现(二)得到的散射光路，为了得到着色点的能量，我们把任务主要分为了下面几个步骤：

1. 判断场景中交点
2. 计算交点的光照能量情况
3. 计算下一条光路的方向
4. 重复上述流程若干次后停止

在上文的伪代码中，需要首先判断当前交点是否在participating media中，并且要知道与media的交点情况。
为了解决上面这些问题，需要引入对于media的采样方法。

最显而易见的思路就是，每隔一段距离，做一次采样来判断是否处在media的环境中，类似于每个箭头处，做一次采样，只要采样够多，肯定可以获得这个信息。

这种方法就像做光栅化一样，为了判断交点，可以用一个“AABB”的包围盒来加速。
虽然在我们的例子中是一个cube box ，实际项目中会用噪声进行腐蚀，但是为了加速完全可以拟定一个粗略的包围盒作为物体的geometry。
这样就不需要时刻都进行采样，而是进入了物体的geometry之后，再进行采样，类似于在外部采样数为0，进入后才开始采样。

进入geometry 之后再采样

这样就解决了第一个问题：光线有没有进入media。

第二个问题是：计算交点的光照能量情况。
进入了participating media之后，需要知道光线与media的交点情况，因为即使进入了media也并非一定会对光线产生散射。
例如有一片很薄的云，薄到比散射的平均自由程都要薄，或者说光线进入了物体的斜角处，都会有还没来得及产生散射效果，就要从media中传播出去的情况。
反之如果media对于光线产生的散射效果，为了计算光照能量，也需要知道光线传播过程中media的beam transmittance，这里的采样方式可以用伪代码表示为：

media采样（ray） {
    采样的距离 = 最小自由程 * （0,1]的随机数
    采样传播的时间 = 采样的距离 / ray的模
    
    if (采样传播的时间 < ray交点的传播时间) {
        存在media的交点
        获得对应点TR
    }
}

这里对于采样的距离，使用了最小自由程 $$1 / \sigma_t $$ 的若干倍，是因为在前期的假设中，已经把participating media划分为了全部均匀的voxel，在每个voxel的长度为最小自由程，那么在其中采样自然是(0,1]的范围中。

接下来判断了ray传播的交点与voxel中交点的时间，如果ray的下一个交点（从participating media出去的点）都比voxel要小，那么可以假设这段media并没有造成scattering effect。否则按照上节的TR计算方式，对于当前点计算TR。

采样理论与优化

在participating media的采样中，涉及到了对于voxel的采样，伪代码中绝大部分都直观易懂。
但是对于采样的随机数这里，我想再多聊两句。

在渲染中，采样采的好可以事半功倍，采样偏差的大能让算力白忙活。采样就像武侠中的内功一样，也许成像算法一模一样，在不同的人手中，也有天差地别的表现。
在surface scattering的渲染中，就在一些情况下会使用对于光源直接采样来代替对于BRDF的采样，这样会让采样效率更高。
对于采样效率，有时均匀的采样并不一定是最好的，例如对于类似sin函数的均匀采样，很有可能重建的结果就是一条直线。

也很容易想到优化方法，在值比较高的地方可以多采样一些（这些地方对于整体贡献比较大），对于值比较低的地方，可以少采样一些。

但是该如何知道哪些地方值高呢？
要是表达式知道的情况下还是选择，有的时候正是因为表达式没有解析解，才需要采样，这样的话岂不是进入了是鸡生蛋，还是蛋生鸡的问题。

例如涉及brdf的积分表达式，断是不可能知道哪里的值高的

这种情况下，经常会采用montecarlo的方法，根据某一种PDF对于表达式进行采样，获得最后的平均，不一定是简单的算术平均，本质是采样的结果 / PDF的算术平均 。

假设有一个积分 $$\int^{1}_{0} 3x^4$$其实一眼能看出来他的结果是1，如果用上面所说的monte carlo的方法来解，首先可以假设某种采样方式，可以获得采样值$$x$$ ,$$x$$ 存在一种概率密度关系$$f(x) = 1/N$$,而这里的pdf积分也是1。

即对于每个采样样本是均匀分布的，是一个常数。从而得到$$f(x)/p(x) = f(x)$$,那么对于montecarlo的结果是可以获得其积分的值的。只要样本足够大，肯定最后都会收敛到期望值。

这样我们现在得到了如何求某个积分的方法

1. 对于积分估计一个合适的pdf
2. 根据pdf的分布进行montecarlo采样

但是问题又来了：

1. 什么叫合适的pdf呢？
2. 应该如何制造出符合pdf分布的采样？

Importance sampling

合适的PDF就是收敛速度快的PDF,假设有两个PDF:PDF1\PDF2,对于同一个sample function进行采样，以收敛到同一个标准差为基准，PDF1需要几千次sample，PDF2需要几万次sample，那么肯定是PDF1更好。

那么怎么样能让PDF收敛的更快呢，importance sampling的方法给出了一种答案：PDF约接近sample function，收敛速度就越快。相对而言，当然可以均匀的采样，不过就要接受效果不好、收敛慢的弊端。

当然重要性采样也并不是唯一的方法，回到之前问题，“当不知道sample function分布的情况下，如何选择采样方式”，还有一种方法会先使用数量较小的样本进行采样，确定大致的图形情况，再在重要的区域安排较多的采样。 这种方法就是adaptive sampling。

这样就获得了什么是合适的pdf的定义，接下来问题就是应该如何制造出符合pdf分布的采样？

Inversion method

在工程实践中，最好获得的随机分布就是0~1的均匀分布了。
但是根据上面的要求，对于不同的pdf，我们需要某种方法来生产符合对应pdf的分布。
那么如果能找到某种转化关系，让工程中最好实现的[0,1]均分分布，能转化成PDF对应的分布，就最好了，类似于存在这种函数：

int seed = 均匀的0-1随机数
pdf(seed) {
// do something

return 符合pdf规范的样本
}

这种函数的构建，一般会采用inversion method 方法来处理：

1. 假设我们有一组随机变量x = [1,2,3,4]; 对应的pd函数为 $$p(x)$$，其密度分布如下：

2. PDF的积分为1。也就是说 $$\int_1^4 p(x) dx = 1$$即，全部概率之和为1.
3. 对于一个选定的随机量$$\xi$$某段pdf被选中的概率，与pdf自己的值（bar的高度）相同。

4. 有了y的值（$$\xi$$）现在可以反求x的值，也就是$$X = F^{-1}(\xi)$$ ，即样本为p3,密度分布pdf（p3）

这样就完成了，根据任意0-1的均分分布，获得符合pdf分布的样本的方法，对于积分的情况，可以把上图的bar无限细分就可以得到积分的情况。

当然inversion method并不是唯一的方法还有“rejection method”，由于在工程中使用的是inversion method，这里就不再多说。

Multiple Importance Sampling

结合importance sampling 的思路和inversion method的方法，终于可以对于sample function获得比较好的采样结果。

但是在实际工程中，经常会针对具体情况，结合使用，例如BRDF模型中，一般由两种方法用于估计反射光线的方向：1. 直接对于BRDF采样 2.直接对于光源采样。

这两种方法都有各自适合的场景：

1. Sampling BRDF

在光源依次减小的情况下，下面的材质依次增加粗糙，可以看到在大光源、光滑材质的情况下lambert BRDF表现不错，但是小光源粗糙材质效果就很差了。
从逻辑上来解释，作为均分分布的BRDF，在大光源的时候，哪怕是均匀采样，很容易就能找到光路，所以效果不错。

2. Sampling lighting
   

同理，对于采样光源的情况，小光源表现很好，对于小光源，如果是均匀采样，则会很难构成光路，但是如果直接对于光源采样，光线利用率就会好很多。

人都是贪心的，能不能结合两者的优点，互相弥补缺点？

当然，这里可以使用 multiple importance sampling 的方法来把两种 importance sampling下的采样方法结合起来，来完成variance redution 的目的。

对于采样的结果，假如颜色比较黑、接近于没有能量的地方，就减弱对于整体的影响；反之如果能量比较多，则说明对于整体光照结果的表现影响更多，则可以增加一些权重，通过这样动态的调节，就可以避免特定的某一种采样方法造成的误差。

在具体实现上，可以通过某种系数来控制每个采样结果对于整体的影响：

例如sampling light 的PDF系数为 $$f$$,则对于光照的权重为：$$\frac{f}{f+g}$$

则sampling BRDF 的PDF系数为 $$g$$,光照的权重则为：$$ 1 - \frac{f}{f+g}$$

在实际工程中，也会通过增加一个指数系数来让效果更加明显。例如：$$\frac{ f^2} {f^2 + g^2}$$

这样就可以让BRDF在无论光源、材质粗糙程度的情况下，都可以获得比较好的效果。

浮点数误差

无论是采样的精度，还是建模的精度，老天爷都有无穷的精度可以使用，而优势也不止于此，计算机天然还有一种精度的劣势：浮点数精度。

在实际的工程实践中，会存在一个无法回避的问题：
由于我们的path tracing存在着诸多的相交判断，无论是光线与材质的，还是光线在voxel传播中，都会有多次的相交判断，并且在participating media中存在对于最小自由程的采样，涉及到许多精细的浮点数判断，这些都埋藏着隐藏的暗雷。

无尽的边界

在participating media的传播中，对于传入、传出media的边界时，一定会涉及到光线在交点的再一次传输。

在误差下，对于交点的判断可能会导致光线无法传出meida的边界

在理想的情况下，光线与media的边界存在交点p,以p为起点沿原本的方向会继续传输，构建一条光路（如上图黑色射线），但是如果是存在误差情况下，起点p可能距离边界还有一点距离，这个时候构建新的射线再去判断交点，又会判断为还在media的范围中，造成无限的loop；

显而易见的方法就是我们对于交点增加一个极小的offset，就像处理zfighting一样，这样就不会有这个问题了。

解决了，但…没有完全解决。
offset的问题是，到底要给多少呢？对于类似垂直的情况可以通过给极小的offset解决，但是如果是接近平行的情况

这种情况下，就需要非常大的offset才能解决，但是对于较大的offset，对于画面的质量势必会有影响。

所以我们需要另外一种方法，更加合理的计算offset，可以通过结合media 边界的信息，给出一个更加合理的offset方向：

1. 获得一个交点p,同时可以得到对应geometry的normal
2. 不直接验证射线方向，而是根据normal的方向，设计出x\y\z三个微小的分量
3. 选择某一个分量作为offset，例如选择一个距离p点最近的offset
4. 得到新的光线追踪起始点，保证和p点在同一侧，继续下一轮追踪

经过这样的调整，由于浮点数精度带来的无限循环问题就会得到解决。

辐射方程的应用

聊完了采样理论，解决了一些精度上的问题，终究还是要落到实际应用中，归根结底我们所需要解决的问题实际就是一个渲染方程：

我们要解决一个问题：看到的这个点，到底是什么颜色。
接下来，就要利用所有前文提到的内容，来完成这道渲染方程
由两个部分组成:

1. 物体的自发光
   对于participating media而言，也有对应的自发光情况。在“散射效果对能量的影响”，中提到了关于自发光的部分：$$L(p,\omega) = \sigma_a§ L_e(p,\omega )$$

2. 物体的散射光
   从四面八方各个方向散射的光，并在观察方向上的光，由于散射的影响，又三个部分组成：负责减益相关的： aborption、out-scattering 与 增益相关的 In-scattering：$$L(p,\omega) = -(\sigma_aL +\sigma_sL) + \sigma_sL_s$$

对于in-scattering term 与表面反射非常相似，也是由类似BRDF的phase function 和 incident radiance组成：$$L_s = \int_{s^2} f(x,\omega, \omega’)Ld\omega’$$

最终,将上面所有的部分组合起来就可以得到 participating media的微分渲染方程，即radiative transfer equation：

这就解决了特定点的着色问题，而针对具体的in-scattering的散射情况，即$$L_s = \int_{s^2} f(x,\omega, \omega’)Ld\omega’$$

方程的求解，其中的incident radiance又是由其他部分的participating media所得到的，将其展开递归进行求解，就可以得到对应的volme rendering equation：

结合radiative transfer equation和volme rendering equation 我们还设计了对应的光路：

1. 第一次bounce，打到box cloud所在的“物体表面”p0，直接透射进box
2. 对于participating medium做最小自由程的采样，看能不能采样到medium，得到与medium中粒子的交点p1
3. 采样到了之后计算p0~p1的transmittance，同时计算光源采样\bsdf采样下的radiance,在box以内的 朝向光源的光路计算transmittance
4. 根据phase function算下一次bounce的方向，重复2-4步，得到p2\p3
5. 当第二步没采样到medium之后，从box出来按照表面散射处理，获得p4点的randiance
6. 累加各个点的radiance 结束

在具体的实现中，我们使用regular tracing的方法来采样homogeneous medium，这种方法类似 ray matching，不过同时也有和ray matching一样的缺点：在统计上是有偏的。

不过同时，对于上面的问题和heterogeneous medium 也可以使用 delta tracing的方法来继续优化。或者直接采用BSSRDF的方法直接对于效果进行拟合。

而且在实现的过程中，工程并不像数学公式那样一帆风顺。
采样理论部分我们解决了采样优化上的问题，通过多重重要性采样解决了不同情况下的效果问题。
其次还有计算机本身引入的浮点数问题，通过结合交点周边情况的信息，增加对应的offset以解决了计算机本身带来的浮点数问题。
最终终于到了最后的效果，这是楼主自己渲染的不带任何后处理的结果，和大家分享一下：

如果有文章中不清楚的地方，也欢迎直接看：对应的代码

希望这三篇文章(体积云一、体积云二)能让大家感受到乐趣与收获。

这…就结束了吗？

明白一种算法是容易的，但是难得是在复杂系统中进行集成，落地到实际的工程实践中。

cs2中的体积云可以和子弹交互，影响了传统的拉烟策略

能否像cs2一样把可交互的体积云带入到piccolo中，后面继续加油带给大家，希望能让大家喜欢。

如果能点个免费的赞支持下楼主就更好了~ 谢谢大家。