蒙特卡洛路径追踪

  • 摘要
  • 1 蒙特卡洛积分(Monte Carlo Integration)
  • 2 蒙特卡洛路径追踪(Monte Carlo Path Tracing)
  • Reference

(本篇文章同步发表于知乎专栏:https://zhuanlan.zhihu.com/p/146714484 欢迎三连关注)

摘要

在上一篇文章中,我们通过对辐射度量学当中一系列概念的定义,引入了渲染方程,一个正确的光线传播模型,但并没有去涉及如何解出该渲染方程,或者说如何通过该渲染方程计算出屏幕坐标上每一个坐标的像素值。在本篇文章中会利用蒙特卡洛路径追踪来完成这个目标。

1 蒙特卡洛积分(Monte Carlo Integration)

首先让我们先搞懂蒙特卡洛路径追踪的这个“蒙特卡洛”的前缀到底指什么。

蒙特卡洛积分的目的: 当一个积分很难通过解析的方式得到答案的时候可以通过蒙特卡洛的方式近似得到积分结果,如下图所示:

显然对于这样一个函数,很难去用一个数学式子去表示,因此无法用一般解析的方法直接求得积分值,而这时候就可以采用蒙特卡洛的思想了。

蒙特卡洛积分的原理及做法: 对函数值进行多次采样求均值作为积分值的近似

该做法十分容易理解,想象一下如果对上图这个函数值进行均匀采样的话,其实就相当于将整个积分面积切成了许许多多个长方形,然后将这些小长方形的面积全部加起来。没错,该做法其实就与黎曼积分的想法几乎一致。但蒙特卡洛积分更加的general,因为它可以指定一个分布来对被积分的值进行采样,定义如下:

如图所示,我们希望求出一个函数f(x)f(x)f(x)在积分域[a,b][a,b][a,b]上的积分值,选定一个采样的分布p(x)p(x)p(x),通过对该分布来进行多次的函数值采样,最后估计的值如图中最下方式子所示。

这里对该式子进行一个简单的推导。相信大家都知道,求均值的做法其实也是对期望的逼近,因此:
1N∑i=1Nf(Xi)p(Xi)≈Ex∼p(x)(f(x)p(x))\frac{1}{N} \sum_{i=1}^{N} \frac{f\left(X_{i}\right)}{p\left(X_{i}\right)} \approx \mathbb{E}_{x\sim p(x)}(\frac{f\left(x\right)}{p\left(x\right)})N1​i=1∑N​p(Xi​)f(Xi​)​≈Ex∼p(x)​(p(x)f(x)​)
那么对于这样一个服从某一分布的期望的计算套公式直接计算得:
EXi∼p(x)(f(Xi)p(Xi))=∫abf(x)p(x)p(x)dx=∫abf(x)dx\begin{aligned} \mathbb{E}_{X_i\sim p(x)}(\frac{f\left(X_{i}\right)}{p\left(X_{i}\right)}) &=\int_{a}^{b} \frac{f(x)}{p(x)}p(x)d x\\ &=\int_{a}^{b} f(x)d x\end{aligned}EXi​∼p(x)​(p(Xi​)f(Xi​)​)​=∫ab​p(x)f(x)​p(x)dx=∫ab​f(x)dx​
通过以上推导即可明白蒙特卡洛的近似正是对积分值的一个无偏估计。

但在本文中为了方便,所有的采样都使用均匀采样,因此很容易推出:

因此,蒙特卡洛在此来说就是一个帮助求得困难积分值的方法

2 蒙特卡洛路径追踪(Monte Carlo Path Tracing)

回顾一下上篇文章中所得到的渲染方程:
Lo(p,ωo)=Le(p,ωo)+∫Ω+Li(p,ωi)fr(p,ωi,ωo)(n⋅ωi)dωiL_{o}\left(p, \omega_{o}\right)=L_{e}\left(p, \omega_{o}\right)+\int_{\Omega^{+}} L_{i}\left(p, \omega_{i}\right) f_{r}\left(p, \omega_{i}, \omega_{o}\right)\left(n \cdot \omega_{i}\right) \mathrm{d} \omega_{i}Lo​(p,ωo​)=Le​(p,ωo​)+∫Ω+​Li​(p,ωi​)fr​(p,ωi​,ωo​)(n⋅ωi​)dωi​

要想解出以上方程的解主要有两个难点:

  1. 积分的计算
  2. 递归形式

而解决这些难点自然就要利用上节中所提到的蒙特卡洛积分方法了。

在进入具体计算之前,对渲染方程做出一点小修改,即舍弃一下自发光项(因为除了光源其他物体不会发光), 以方便进行计算推导:
Lo(p,ωo)=∫Ω+Li(p,ωi)fr(p,ωi,ωo)(n⋅ωi)dωiL_{o}\left(p, \omega_{o}\right)=\int_{\Omega^{+}} L_{i}\left(p, \omega_{i}\right) f_{r}\left(p, \omega_{i}, \omega_{o}\right)\left(n \cdot \omega_{i}\right) \mathrm{d} \omega_{i}Lo​(p,ωo​)=∫Ω+​Li​(p,ωi​)fr​(p,ωi​,ωo​)(n⋅ωi​)dωi​

从具体例子出发,首先仅仅考虑直接光照:

Lo(p,ωo)=∫Ω+Li(p,ωi)fr(p,ωi,ωo)(n⋅ωi)dωiL_{o}\left(p, \omega_{o}\right)=\int_{\Omega^{+}} L_{i}\left(p, \omega_{i}\right) f_{r}\left(p, \omega_{i}, \omega_{o}\right)\left(n \cdot \omega_{i}\right) \mathrm{d} \omega_{i}Lo​(p,ωo​)=∫Ω+​Li​(p,ωi​)fr​(p,ωi​,ωo​)(n⋅ωi​)dωi​

观察该修改过之后的方程其实就只是一个单纯的积分计算了,其物理含义为着色点p到摄像机或人眼的Radiance值。

回想第一章所提的,对于一个困难积分只要选定一个被积分变量的采样分布即可通过蒙特卡洛的方法得到积分结果的近似值,而此时的被积分值为ωi\omega_iωi​,选定ωi∼p(ωi)\omega_i \sim p(\omega_i)ωi​∼p(ωi​),不难得出积分近似结果如下:
Lo(p,ωo)≈1N∑i=1NLi(p,ωi)fr(p,ωi,ωo)(n⋅ωi)p(ωi)L_{o}\left(p, \omega_{o}\right) \approx \frac{1}{N} \sum_{i=1}^{N} \frac{L_{i}\left(p, \omega_{i}\right) f_{r}\left(p, \omega_{i}, \omega_{o}\right)\left(n \cdot \omega_{i}\right)}{p\left(\omega_{i}\right)}Lo​(p,ωo​)≈N1​i=1∑N​p(ωi​)Li​(p,ωi​)fr​(p,ωi​,ωo​)(n⋅ωi​)​

正如一开始所说,先单独考虑直接光照,因此只有当采样的方向ωi\omega_iωi​击中光源的时候,光源才会对该着色点有贡献,计算伪代码如下:

显而易见的,单独仅仅考虑直接光照自然是不够的,还需要间接光照,即当采样的ωi\omega_iωi​方向碰撞到了别的物体,如下图所示:

此时采样的光线碰撞到了另一个物体的Q点,那么该条路径对着色点P的贡献是多少呢?自然是在点Q的直接光照再乘上反射到该方向上的百分比了!显然这是一个类似光线追踪的递归过程,不同在于该方法通过对光线方向的采样从而找出一条条可行的路径,这也正是为什么叫路径追踪的原因,伪代码如下:

至此,我们成功通过蒙特卡洛的方式解出了渲染方程的积分值,也通过考虑直接光照与间接光照解决了递归的问题。但该方法至此有一个非常致命的问题:

我们通过每次对光线方向的采样从而解出方程,假设每次采样100条,那么从人眼出发的第一次采样就是100条,在进行第二次反射之后就是10000条,依次类推,反射越多次光线数量便会爆炸增长,计算量会无法负担,那么如何才能使得光线数量不爆炸增长呢?唯有每次只采样一个方向!N=1

每次如果只采样一个方向那么所带来的问题也是显而易见的,积分计算的结果会非常的noisy,虽然蒙特卡洛积分是无偏估计,但样本越少显然偏差越大。但该问题很好解决,如果每次只去寻找一条路径结果不好,那么重复多次寻找到多条路径,将多条路径的结果求得平均即可!如下图所示:

改良之后的Path Tracing伪代码如下:

通过对经过像素的光线重复采样,每次在反射的时候只按分布随机选取一个方向,解决了只对经过像素的光线采样一次,而对反射光线按分布采样多次所导致的光线爆炸问题。

那么现在所有的问题都解决了吗?还没有!因为shade函数的递归没有出口,永远不会停下。
但这里并不没有采用类似光线追踪当中设定反射深度显示的给出递归出口的方法,而是非常精妙的采用了俄罗斯轮盘赌(Russian Roulette)

给你一把左轮,两发子弹,你不知道哪一发会真正的射出子弹,因此拿这把左轮射自己,你有4/6的概率活下来,这就是俄罗斯轮盘赌的概念。

将其应用在路径追踪当中,首先设定一个概率PPP, 有P的概率光线会继续递归并设置返回值为Lo/PL_o /PLo​/P,有1−P1-P1−P的概率光线停止递归,并返回0。这样巧妙的设定之下光线一定会在某次反射之后停止递归,并且计算的结果依然是无偏的,因为Radiance的期望不变,证明如下:
E=P⋆(Lo/P)+(1−P)⋆0=L0E=P^{\star}(L o / P)+(1-P)^{\star} 0=L_{0}E=P⋆(Lo/P)+(1−P)⋆0=L0​

shade函数的伪代码变更如下,使得可以停止递归了:

至此,我们的路径追踪算法已经完成大半,只差最后一个小问题!现在的路径追踪效率非常的低下,如图所示:

在每次计算直接光照的时候,通过均匀采样任选一个方向,但很少会的光线可以hit光源,尤其当光源较小的时候,这种现象越明显,大量采样的光线都被浪费了。

因此在计算直接光照的时候改进为**直接对光源进行采样!**这样所有采样的光线都一定会击中光源(如果中间没有别的物体),没有光线再会被浪费了。假设光源的面积为A,那么对光源进行采样的 pdf=1/Apdf = 1/Apdf=1/A (因为∫pdf⁡dA=1\int \operatorname{pdf} d A=1∫pdfdA=1),但原始的渲染方程:
Lo(p,ωo)=∫Ω+Li(p,ωi)fr(p,ωi,ωo)(n⋅ωi)dωiL_{o}\left(p, \omega_{o}\right)=\int_{\Omega^{+}} L_{i}\left(p, \omega_{i}\right) f_{r}\left(p, \omega_{i}, \omega_{o}\right)\left(n \cdot \omega_{i}\right) \mathrm{d} \omega_{i}Lo​(p,ωo​)=∫Ω+​Li​(p,ωi​)fr​(p,ωi​,ωo​)(n⋅ωi​)dωi​
很明显是对光线方向ωi\omega_iωi​进行积分的,如果想要对光源进行采样的并依然使用蒙题卡洛的方法,那么一定要将其修改为对光源面积 dA的积分,换言之就是需要找到dA与dωi\omega_iωi​的关系即可。如下图所示:

关系式中的cosθ′cos\theta^{\prime}cosθ′是为了计算出光源上微分面积元正对半球的面积,之后再按照立体角的定义dω=dAr2\mathrm{d} \omega=\frac{\mathrm{d} A}{r^{2}}dω=r2dA​,除以着色点x与光源采样点x’距离的平方即可。于是根据图中二者的关系可将渲染方程改写如下:
Lo(x,ωo)=∫Ω+Li(x,ωi)fr(x,ωi,ωo)cos⁡θdωi=∫ALi(x,ωi)fr(x,ωi,ωo)cos⁡θcos⁡θ′∥x′−x∥2dA\begin{aligned} L_{o}\left(x, \omega_{o}\right) &=\int_{\Omega^{+}} L_{i}\left(x, \omega_{i}\right) f_{r}\left(x, \omega_{i}, \omega_{o}\right) \cos \theta \mathrm{d} \omega_{i} \\ &=\int_{A} L_{i}\left(x, \omega_{i}\right) f_{r}\left(x, \omega_{i}, \omega_{o}\right) \frac{\cos \theta \cos \theta^{\prime}}{\left\|x^{\prime}-x\right\|^{2}} \mathrm{d} A \end{aligned}Lo​(x,ωo​)​=∫Ω+​Li​(x,ωi​)fr​(x,ωi​,ωo​)cosθdωi​=∫A​Li​(x,ωi​)fr​(x,ωi​,ωo​)∥x′−x∥2cosθcosθ′​dA​
这样便成功从ωi\omega_iωi​积分转到了对光源面积A的积分,就可以利用蒙特卡洛的方法对光源进行采样从而计算直接光照的积分值了,对于间接光照,依然采用先前的方法进行光线方向的均匀采样。最终伪代码如下,分直接光照和间接光照两部分计算:

tips:计算直接光照的时候还需要判断光源与着色点之间是否有物体遮挡,该做法也很简单,只需从着色点x向光源采样点x’发出一条检测光线判断是否与光源之外的物体相交即可,如图所示:

最后以一张闫老师课程的path tracing作业的截图作为结束!

如果本文对你有帮助的话求点赞求收藏求一个大大的关注

计算机图形学十五:基于物理的渲染(蒙特卡洛路径追踪)相关推荐

  1. [图形学] 基于物理的渲染(PBR)

    reference :<real-time rendering 4> BRDF 概述 基于物理的渲染即计算沿着视线进入相机的光照辐射.如果不考虑吸收或散射的介质,则进入相机的辐射等于离开相 ...

  2. 【基于物理的渲染(PBR)白皮书】(五)几何函数相关总结

            本文由@浅墨_毛星云 出品,首发于知乎专栏,转载请注明出处           文章链接: https://zhuanlan.zhihu.com/p/81708753 在基于物理的渲染 ...

  3. 图形学基础 | 基于物理的渲染理论(PBR)

    转载自: https://learnopengl.com/PBR/Theory Learn OpenGL PBR Theory PBR 基于物理的渲染(Physically Based Renderi ...

  4. (二十)unity shader之——————基于物理的渲染技术(PBS):下篇(PBS技术拓展:全局光照、伽马校正、HDR)

    前面两篇文章我们介绍了PBS实现的数学和理论基础,和standard shader的原理和实现,还有一些其他的渲染相关的unity技术.其中有些概念和技术没有讲的很详细,现在对这些重要的概念进行更深入 ...

  5. (十九)unity shader之——————基于物理的渲染技术(PBS):中篇(Unity 5中的Standard Shader的实现和使用)

    一.unity 5中的standard shader 在unity5中新创建一个模型或是新创建一个材质时,默认使用的着色器都是一个名为standard 的着色器.这个standard shader使用 ...

  6. PBR:基于物理的渲染(Physically Based Rendering)+理论相关

    一: 关于能量守恒 出射光线的能量永远不能超过入射光线的能量(发光面除外).如图示我们可以看到,随着粗糙度的上升镜面反射区域的会增加,但是镜面反射的亮度却会下降.如果不管反射轮廓的大小而让每个像素的镜 ...

  7. 【基于物理的渲染(PBR)白皮书】(二) PBR核心理论与渲染光学原理总结

            本文由@浅墨_毛星云 出品,首发于知乎专栏,转载请注明出处           文章链接: https://zhuanlan.zhihu.com/p/56967462 这是[基于物理的 ...

  8. 3D纹理映射(噪点,基于物理的渲染)

    3D纹理映射(噪点,基于物理的渲染) 1.简介 •2D纹理颜色/扩散映射粘贴粘贴图片(纹理)到对象的表面上 •在3D对象的表面上粘贴2D纹理有多种方法 •3D纹理映射将对象放置在3D模式中 •3D模式 ...

  9. 基于物理的渲染(PBR)白皮书 | 迪士尼原则的BRDF与BSDF相关总结

    基于物理的渲染(Physically Based Rendering , PBR)技术,自迪士尼在SIGGRAPH 2012上提出了著名的"迪士尼原则的BRDF(Disney Princip ...

最新文章

  1. spring加载properties文件
  2. 厉害了!Antiilatency推出移动位置追踪器!
  3. Eclipse 配置 maven 的两个 settings 文件
  4. javascript写的关于静态页面获取URL传递参数的函数[原创]
  5. vue1升级到vue2的问题
  6. Python键鼠操作自动化库PyAutoGUI简介
  7. 面板大小调整_3天学会premiere完全自学教程-了解时间轴面板
  8. 【LeetCode】剑指 Offer 68 - I. 二叉搜索树的最近公共祖先
  9. 基于JAVA+Servlet+JSP+MYSQL的网上订餐管理系统
  10. 移动端上下拖动调整顺序效果_HTML5 移动端的上下左右滑动问题
  11. PgSQL · 特性分析 · 金融级同步多副本分级配置方法
  12. 计算机英语(王艺)译文(unit6-unit12)
  13. 解释一下什么是vue实例
  14. 如何给网页设置自己的字体
  15. STC8H开发(二): 在Linux VSCode中配置和使用FwLib_STC8封装库(图文详解)
  16. 可变剪接分析流程(rMATS)
  17. 缓存(Cookie,SessionStorage,localStorage)详解
  18. 我的中国“芯”——资深后端工程师成长分享——“胡”说IC工程师完美进阶
  19. 积分墙广告的七个真相(触控软文)
  20. Excel的常用快捷键

热门文章

  1. 云计算的前世今生(上)
  2. RAD Studio 10.3.x RIO 常规快捷键操作
  3. 计算机本科应届生薪资大多是多少?外行人18k垫底25k人均水平
  4. 2014年4月7日 再次深入激辩余额宝
  5. 一些著名的软件都用什么语言编写?C/C++ yyds
  6. 免费开源智慧农业物联网云平台 V3.0.1.2含源码
  7. 《Google软件测试之道》—第2章2.4节与工具开发工程师Ted Mao的访谈
  8. 通信协议原理及应用CAN基础知识
  9. QNX BSP包目录结构
  10. java基础知识总结大全(8万多字)