光子映射学习#

主要按[Realistic Image Synthesis Using Photon Mapping]这个书来学习光子映射
[A Practical Guide to Global Illumination using Photon Maps]大致是其精简版,也包括一些不同信息可以参考。

持续更新


传统光线追踪至少有两个问题。
- diffuse inter-reflection
没理解他说的是颜色不准还是指noise问题
  • caustics问题

这两个问题都是间接光问题。光子映射可以利用预计算的方法来解决。
是光线追踪的一个补充,能更有效地解决问题。
光子映射是一个two-pass方法。目的是计算diffuse面的间接光。
第一遍从光源发出光子,打到非镜面时存下信息。所有光子都在diffuse面上。
第二遍用统计方法,根据第一遍的光子映射信息算出radiance。

一个重点是光子映射与场景的geo信息是解耦的。所以可以用于非常复杂的场景。

1. Photon emission#

每个光子带相等的flux(辐射功率)。

点光源,光子的发射是均匀分布。
方相光源,从场景外均匀发射。
面光源?向半球发射,cos分布,normal方向概率大。

一个光源的总功率是\(P_{light}\),那么发射\(n_e\)个光子,每个光子的功率

\({\Large P_{photon} = \frac{P_{light}}{n_e}}\)

多光源#

每个光源都要发射光子。

并不是多光源就需要更多光子,维持一个总量即可,各个光源按能量分配光子数。
还可使用importance sampling,对某些光源侧重。

2. Photon tracing#

photon tracing与ray tracing类似。
区别:光子追踪是散播flux,光线追踪是收集radiance
因此与材质的交互也有不同。

一个例子是对于radiance,折射由折射率决定。光子不是这样算。 光子打到物体上,可能反射,传播,或吸收。概率由material的反射率决定。

这里和ray tracing有巨大的区别。
我一开始没有仔细看这一节,而是直接用ray tracing的方法计算光子经过brdf处理后的能量,这样不断往下传播。也能出似是而非的图,但是有一系列问题。
论文实际需要用Russian roulette来决定对光子的操作。
一来可以有效的去除不重要的光子,减少计算量。
二来可以使每个光子的能量差不多。后续计算能量时效果更好。brdf处理后能量可能差距巨大。

反射还是吸收?#

设r为我们取的反射概率

p = randint() # [0-1]的随机数
if p < r:
    以原始能量反射
else:
    吸收
为什么这样操作是对的。比如材质r=0.7,有100个光子。
0.7的两种解读: 1. 反射100个光子,每个的能量为0.7。 2. 反射70个光子,每个的能量不变。

两种都是对的。我们就可以使用第二种解读。保持正确性的前提下,可以有效地减少计算。

反射,传播,还是吸收?#

以单色为例,设d为diffuse反射概率,s为specular反射概率。 可以按三个概率区间决定光子如何传播。

[0, d] diffuse反射
[d, s + d] specular反射
[s + d, 1] 吸收

对于多色。

Photon storing#

只有光子碰到diffuse面时才储存。储存在photon map里。
每个光子在路径中可存多次。
存的数据为[位置,能量,入射角度]。

扩展到participating media#

使用光子映射能更容易地处理participating media。比如烟雾。

三个photon map#

为了效率,把光子分到三个map里。具体原因后续会讲。

  • Caustic photon map
    在碰到diffuse面之前至少有过一次镜面反射的光子。\({\Large LS^+D}\)
  • Global photon map
    所有diffuse面的GI。\({\Large L\{S|D|V\}^*D}\)
  • Volume photon map
    participating media的间接光。\({\Large L\{S|D|V\}^+V}\)

上面用了[Heckbert90]的描述光子路径的语法。

  • L 从光源发射

  • S 镜面反射或传输

  • D 非镜面反射或传输

  • V volume散射

  • {x|y|z} xyz中的一个

  • \(x^+\) 1个或以上x

  • \(x^*\) 任意个x

kd-tree#

需要频繁地查找一个光子附近的光子。需要适应非常不均匀分布的物体。
这里最常用的数据结构是平衡kd-tree。需要自己实现。

3. 计算面上的辐射率#

photon map可以看成输入功率的表示,计算radiance就是要对其积分。

\(x\)上的reflected radiance:

\({\Large L_r(x,\omega) =\int_{\Omega_x}f_r(x,\omega', \omega)L_i(x,\omega')|cos\theta_i|d\omega_i'}\)

\(\omega\)为出射方向,\(\omega'\)为入射方向,\(L_i\)为入射radiance。
半球上所有方向(\(\omega'\))的入射radiance,经过brdf和cos处理后积起来。

需要把radiance和flux的关系展开。

\({\Large E = \frac{d\Phi}{dA}}\)

\({\Large L=\frac{dE}{cos\theta dω} = \frac{d^2\Phi}{cos\theta dA d\omega}}\)

\({\Large L_i(x,\omega')=\frac{d^2\Phi_i(x, \omega')}{cos\theta_i d\omega_i'dA_i}}\)

代入\(L_r\)

\({\Large L_r(x,\omega) =\int_{\Omega_x}f_r(x,\omega', \omega)\frac{d^2\Phi_i(x, \omega')}{d\omega_i'dA_i}d\omega_i'}\)

\({\Large L_r(x,\omega) =\int_{\Omega_x}f_r(x,\omega', \omega)\frac{d^2\Phi_i(x, \omega')}{dA_i}}\)

也可以写成

\({\Large L_r(x,\omega) =\int_{\Omega_x}f_r(x,\omega', \omega)dE_i(x, \omega')}\)

下面做近似。
输入功率\(\Phi_i\)用距离点\(x\)最近的n个光子来计算。每个光子的功率是\({\Large \Delta \Phi_p(\omega_p)}\)
得到

\({\Large L_r(x,\omega) \approx \sum_{p=1}^{n}f_r(x,\omega_p', \omega)\frac{\Delta \Phi_p(x, \omega_p)}{\Delta A}}\)

从点x生成一个球,不断变大直到包含n个光子。用此时的半径算A。如图8所示。

\({\Large L_r(x,\omega) \approx \frac{1}{\pi r^2} \sum_{p=1}^{n}f_r(x,\omega_p', \omega) \Delta \Phi_p(x, \omega_p)}\)

这个近似包含很多假设,准确度依赖于photon map和这个公式中采用的光子数量。
数量越多越准确。
球体在物体边缘或墙角可能包含错误的光子,造成结果错误。
这些错误区域也是跟光子数量相关。
如果n趋近于无穷大,那么近似可以写成等号。

可以用其他形状代替球体,各有好坏。 球体比较好算。

从我测试来看,例如发射20000个光子,生成100000个左右的diffuse光子,4000左右的caustic光子。
4000其实挺少的。计算能量的时候某个交点周围会取n个caustic光子,这个n如何取值?
如果定的比较大,那么会到较大范围去找齐光子,那么r可能会很大,就会有严重问题。
跑老远找光子,贡献又不大,但会把r拖大,造成把整体能量变的很小。
如果定的很小,那么图像就类似青一块紫一块。
还有边界问题。比如在一个突出的墙角,可能会穿过墙壁找光子。

3.1 Filtering#

如果photon map中的光子数量太少,边缘的radiance就会模糊。
这种artifact有时有益,有时有害。

要减少这种模糊,估算的radiance需要做过滤。给越靠近点x的光子增加越多权重。

收集光子用的是球体,但过滤并不是3d的。收集的光子都看作是在同一微小表面。

cone filter#

\({\Large \omega_{pc}=1-\frac{d_p}{kr}}\)

加上过滤和normal

\({\LARGE L_r(x,\omega) \approx \frac{\sum_{p=1}^{n}f_r(x,\omega_p', \omega) \Delta \Phi_p(x, \omega_p)\omega_{pc}}{(1-2/3k)\pi r^2}}\)

\(k\ge 1\)

Gaussian filter#

\({\Large \omega_{pg}=\alpha(1-\frac{1-e^{-\beta \frac{d_p^2}{2r^2}}}{1-e^{-\beta}})}\)

\(\alpha = 0.918\)

\(\beta = 1.953\)

\({\Large L_r(x,\omega) \approx \frac{1}{\pi r^2} \sum_{p=1}^{n}f_r(x,\omega_p', \omega) \Delta \Phi_p(x, \omega_p)\omega_{pg}}\)

3.2 rendering#

photon map可以独立使用。有几种方法进行可视化。
暂不考media。

总的输出radiance为自发的radiance加反射出的radiance。

\({\Large L_o(x,\omega) =L_e(x,\omega) + L_r(x,\omega)}\)

反射radiance为半球上所有方向的入射radiance,经过cos项和brdf处理,积分。

\({\Large L_r(x,\omega) =\int_{\Omega_x}f_r(x,\omega', \omega)L_i(x,\omega')|cos\theta_i|d\omega_i'}\)

\(L_r\)可用蒙特卡洛光线追踪算出。即我们之前实现的路径追踪。
这些方法比较耗时。而光子映射更高效。

brdf分成specular和diffuse两部分

\({\Large f_r(x,\omega', \omega) = f_{r,s}(x,\omega', \omega) + f_{r,d}(x,\omega', \omega)}\)

\(L_i\)分成三部分

\({\Large L_{i,l}(x,\omega')}\) 光源的直接光

\({\Large L_{i,c}(x,\omega')}\) caustics。间接光。至少经过一次specular反射或折射,再打到diffuse面。

\({\Large L_{i,d}(x,\omega')}\) 间接光。至少一次diffuse反射。

\({\Large L_{i}(x,\omega') = L_{i,l}(x,\omega') + L_{i,c}(x,\omega') + L_{i,d}(x,\omega')}\)

最终的\(L_r\)由4项组成

\({\Large L_r(x,\omega) =\int_{\Omega_x}f_r(x,\omega', \omega)L_i(x,\omega')|cos\theta_i|d\omega_i'=}\)

\({\Large \int_{\Omega_x}f_r(x,\omega', \omega)L_{i,l}(x,\omega')|cos\theta_i|d\omega_i' + }\)

\({\Large \int_{\Omega_x}f_{r,s}(x,\omega', \omega)(L_{i,c}(x,\omega')+L_{i,d}(x,\omega'))|cos\theta_i|d\omega_i' + }\)

\({\Large \int_{\Omega_x}f_{r,d}(x,\omega', \omega)L_{i,c}(x,\omega')|cos\theta_i|d\omega_i' + }\)

\({\Large \int_{\Omega_x}f_{r,d}(x,\omega', \omega)L_{i,d}(x,\omega')|cos\theta_i|d\omega_i' }\)

其中

直接光#

\({\Large \int_{\Omega_x}f_r(x,\omega', \omega)L_{i,l}(x,\omega')|cos\theta_i|d\omega_i'}\)

最重要的部分。

  • 准确计算
    用蒙特卡洛光线追踪。如果交点和光源之间有其他物体,那么直接剔除。
    另有一种shadow photons方法。
  • 近似计算
    用global photon map计算。

Specular和glossy反射#

\({\Large \int_{\Omega_x}f_{r,s}(x,\omega', \omega)(L_{i,c}(x,\omega')+L_{i,d}(x,\omega'))|cos\theta_i|d\omega_i' }\)

这里不必使用任何photon map。因为这个值由\({\large f_{r,s}}\)主导,只在镜面反射方向的狭小范围内起作用。

也可以用pm优化,但需要大量光子。不合算。

Caustics#

\({\Large \int_{\Omega_x}f_{r,d}(x,\omega', \omega)L_{i,c}(x,\omega')|cos\theta_i|d\omega_i'}\)

  • 准确计算
    需要从caustics photon map计算功率。光子数量越大结果越准。
  • 近似计算
    从global photon map计算。

Multiple Diffuse Reflections#

\({\Large \int_{\Omega_x}f_{r,d}(x,\omega', \omega)L_{i,d}(x,\omega')|cos\theta_i|d\omega_i' }\)

计算至少有一次diffuse反射的光路。因为有diffuse,所以这个值比较“柔和”。

  • 准确计算
    要用蒙特卡洛光线追踪。
    通常蒙特卡洛计算diffuse间接光比较低效,但是这个积分有一些特殊性质,最终效果还行。
    间接光很柔和,因为我们已经把caustics分出去了。caustics通常是蒙特卡洛光追中noise的主因。
  • 近似计算
    用global photon map。

可以看到所有近似计算都可以用global photon map。global photon map本质上包含了所有光照!


4. Participating Media#

之前的讨论中,都是考虑光子与mesh表面碰撞,严格来说这只发生在真空中。
即使是很纯净的空气也会改变传播方向。
在更复杂或要求更高的场景中,就不能忽略这个效果,比如烟雾、云等等。
这类介质叫做participating media。
另一类participating media是有些透明度的材质,比如大理石、皮肤、植物等。
涉及到subsurface。要解决这些问题,得解决基本的participating media。

光子映射很适合解决这类问题。是第一个完全模拟出subsurface的算法。

4.1 Light Scattering in Participating Media#

当一个光子进入一种介质,可以不受影响继续往前走,或者在某个位置与介质交互。

交互的结果:要么被吸收,要么被发散出去(scattered)。
发散因子\({\Large \sigma_s}\)(scattering coefficient)
吸收因子\({\Large \sigma_a}\)(absorption coefficient)
这些因子的单位是百分比每单位距离。比如选定0.2/米,可以理解为光子每移动1米,有20%概率被吸收。
或者说100个光子移动1米,有20个会被吸收。
或者说光子在这个介质中,被吸收前,平均能移动1/0.2=5米。

当光穿过介质,可以看作是radiance的连续改变。

在某一点x处,\(\omega\)方向
对于发散(out-scattering),改变为

\({\Large (\omega\cdot\nabla)L(x, \omega)=-\sigma_s(x) L(x, \omega)}\)

对于吸收,改变为

\({\Large (\omega\cdot\nabla)L(x, \omega)=-\sigma_a(x) L(x, \omega)}\)

合起来,对于\({\Large \sigma_t = \sigma_s + \sigma_a}\)

\({\Large (\omega\cdot\nabla)L(x, \omega)=-\sigma_t(x) L(x, \omega)}\)

穿过介质的过程中,存在光的in-scattering,所以会有个gain。

\({\Large (\omega\cdot\nabla)L(x, \omega)=\sigma_s(x) \int_{\Omega_4\pi} p(x, \omega', \omega) Li(x, \omega')d\omega'}\)

\(\omega'\)是入射方向。跟渲染方程类似,积分是求所有方向的入射光的贡献。
其中p是phase function,描述光的分布。

另外如果介质本身会发光,又有一个gain。

\({\Large (\omega\cdot\nabla)L(x, \omega)=\sigma_a(x) L_e(x, \omega)}\)

一共4个改变,两个加两个减。合起来有

\({\Large (\omega\cdot\nabla)L(x, \omega) =}\)

\({\Large dL(x, \omega) = \sigma_a(x) L_e(x, \omega) +\sigma_s(x) \int_{\Omega_4\pi} p(x, \omega', \omega) Li(x, \omega')d\omega' -\sigma_t(x) L(x, \omega)}\)

4.2 体渲染方程#

要解出\({\Large L(x, \omega)}\)有一堆数学上的推导,对于我们数学不熟的人比较跳跃。
论文里直接就给答案,没有推导,只说了一句,对两边同时积分一段距离s。
首先对于 \({\Large dL(x, \omega)=-\sigma_a(x) L(x, \omega)}\)
L在某个x处的改变量,等于该处的L值乘一个系数。
本质是一个微分方程。我们先要学一下基本的微分方程,一般微积分书里都有。
比如托马斯微积分第16章。
另我的文档一阶微分方程基础学习。其中也包含了一些微积分的重要但容易忽略的知识。
L的这个问题本质和微积分教材里的人口问题、阻力问题是一样的。
都是变化率与实时值相关,然后包含初始值问题。

对于人口问题的微分方程

\({\Large \frac{dP}{dt}= kP,P(0) = P_0}\)

\({\Large P = P_0e^{kt}}\)

变量是时长t。但我们的L变量是任意点x,而且形式不一样,能否转成一个标准的微分方程?

对于人口来说是t从0开始,经过一段时间之后P的情况。
类似可以用距离s的思想来重写,表示成光子从0开始 移动一段距离s之后L的情况。
那么

\({\Large dL(x, \omega)=-\sigma_a(x) L(x, \omega)}\)

转变成微分方程

\({\Large \frac{dL(s)}{ds}=-\sigma_a(s) L(s)}\)

这样就舒服了,可以用这个思想重写整个式子。

\({\Large \frac{dL(s)}{ds} = \sigma_a(s) L_e(s) +\sigma_s(s) \int_{\Omega_4\pi} p(s, \omega') Li(s, \omega')d\omega' -\sigma_t(s) L(s)}\)

写成微分方程标准形式

\({\Large L'(s) + \sigma_t(s) L(s) = \sigma_a(s) L_e(s) +\sigma_s(s) \int_{\Omega_4\pi} p(s, \omega') Li(s, \omega')d\omega'}\)

\({\Large L'(s) + p(s) L(s) = q(s)}\)

此时p和q可以说都是确定的,只需解微分方程求L。 可参考上述的一阶微分方程基础学习


只包含吸收的情况#

如果简化一下,只考虑吸收,那么q为0。有

\({\Large L'(s) + p(s) L(s) = 0}\)

\({\Large L(s) = Ce^{-P(s)}}\)

给定初始值L(0),可以算出C,进而算出L(s)。

\({\Large L(0) = Ce^{-P(0)}}\)

\({\Large C = L(0) e^{P(0)}}\)

\({\Large L(s) = L(0) e^{P(0)}e^{-P(s)} = L(0) e^{P(0)-P(s)}}\)

也可以写成

\({\Large L(s) = L(0) e^{\int_{s}^{0}p(t)dt}}\) (t只是一个变量符号,可以随意换,写成t避免冲突)

甚至可以不写t

\({\Large L(s) = L(0) e^{\int_{s}^{0}p}}\)

写成从0积到s更直观,加个负号。

\({\LARGE L(s) = L(0) e^{-\int_{0}^{s}p_s(t)dt}}\)

这里的s是经过的距离,即从初始位置积到某个距离s。
如果取L(0)为1,L(s)就是积到s时,对于初始值的倍数!

如果是均匀的介质,Ps为常数k,有

\({\LARGE L(s) = e^{-ks}}\)

也可以按具体的点考虑,那么就是从点x1积到点x2。

\({\LARGE L(s) = L(x_1) e^{-\int_{x_1}^{x_2}p_x(t)dt}}\)

也就是上一节说的两种思想,一个是针对距离,一个是针对某个坐标。
那么p相应地应该有两个版本\({\Large p_s}\)\({\Large p_x}\)。最后的结果肯定一样。

完整的情况#

\({\Large L'(s) + p(s) L(s) = q(s)}\)

\({\LARGE L(s) = \int_{a}^{s} e^{-\int_{t}^{s}p(t')dt'} q(t)dt + Ce^{-\int_{a}^{s}p(t')dt'} }\)

按[一阶微分方程学习]最后的讨论。选定a=0,结合初始值L(s=0)确定C

\({\Large L(0) = C }\)

代回原方程得到最终的体渲染方程

\({\LARGE L(s) = \int_{0}^{s} e^{-\int_{t}^{s}p(t')dt'} q(t)dt + L(0)e^{-\int_{0}^{s}p(t')dt'} }\)

意思就是一个初始能量为L(0)的光子,在介质中移动s距离后的最终能量。


4.3 Phase Function#

phase function描述散射的光如何在介质中分布。

p定义再某点上某个方向w对某个方向w’的贡献值。

p必须满足 给定一个方向w,所有方向w’射过来的能量对w方向的贡献,积起来值为1。

\({\Large \int_{\Omega_4\pi}^{}p(x, \omega', \omega) d\omega' = 1}\)

p和brdf相似,但是有两个大不同: p没有单位,且是normalized。

可以用入射光和出射光的夹角\(\theta\)简写

各向同性#

是个常数,各个方向均匀分布

\({\Large p(\theta) = \frac{1}{4\pi}}\)

Henyey-Greenstein#

最常用的phase function。最初被用来解释星尘的散射,后来被用来描述海洋、云、皮肤、石头等等。

\({\LARGE p(\theta) = \frac{1-g^2}{4\pi(1+g^2-2gcos\theta)^{1.5}}}\)

正的g为前向散射,负的g为后向散射,0为各向同性。
\(\theta=0\)为前向,\(\theta=\pi\)为反向。
可以画出不同g值对应的图形。可以看到g=0时,是一个圆形,任何角度都是\(1/4\pi\)
\(g^2\)越趋近于1,圆形越扁,越分布在正向和反向附近。
g值正地越多越往正向集中,负地越多越往负向集中。

4.4 Ray Marching#

体渲染方程,除了一些简单场景之外,只能用数值积分来计算。
L是关于距离s的函数,很自然想到模拟光子不断前进,就是所谓的Ray Marching
把前进路线分成很多小步,每步算一下体渲染方程。必要时做一些简化。

那么就是一个近似的解,细分越多越精确。

啊照这样说,那是不是根本就不用推导那一大片微分方程?每一步直接算最初的dL和更新L不就行了?


4.4.1 自适应ray matching#

上面的ray marching时平均步长。如果是非均匀的材质,或者有本地波动(比如阴影),最好用自适应步长。 例如当材质特性发生较大变化时,可采用较小的步长,以更好地反映出细节。

方法1 可以以递归的形式来操作,从两个端点开始,从中间破开,左右递归下去。
直到距离达到设定的限制,或者两个端点的差别小到一定程度(不用继续破开)。

方法2 根据sigma_t来调整步长。这个对于未知的材质可能很有效。

取均匀分布的随机数

貌似涉及均匀采样


4.5 Photon Tracing#

4.5.1 Photon Scattering#

两种情况 光子与介质交互,或发散或被吸收。 光子与某表面交互,反射或传输或吸收。

光子进入介质时,不看作与表面交互,而是直接往里面走,进行与介质的交互。

之前已提到sigma的意义,光子进行下一个交互前平均移动的距离。

光子如果发生scattering,对phase函数进行采样,得到新的方向。
推导见(从0开始学习PBRT)。

4.5.2 photon storing#

每当光子与介质进行交互,就要存光子,无论是吸收还是散射。

对于介质,我们起一个独立的volume photon map。 因为介质交互中的radiance的计算方式与表面交互不一样。

一个优化是,只存至少散射过一次的光子。 这样省去了直接光的部分。直接光可以用传统ray tracing高效算出,不需要费劲用光子来算。

4.5.3 photon emission#

光子也可以从介质中发出。


4.6 radiance的计算#


4.7 渲染#


4.8 次表面散射#

次表面散射存在于所有非金属材质。

一般会简化掉,做成光子碰到表面后立即从交点往外反射离开,而不是往里面散射。 这样简化对于皮肤、玉石、牛奶等材质来说,效果就很差。

要准确地处理translucent材质,就需要遵循事实:光从某点进入,会经过次表面散射,从其他点出去。

photon mapping是第一种完全实现这个原理的算法。

photon tracing#

同上述的photon tracing

渲染#



结果#

还有较多的细节没仔细看。

目前只做了点光源。

光子能量不知道怎么取值。大概加了个scale,比如5000,能得到还算正常的图。

difffuse光子数量通常是caustic数量的5倍。

光子数量问题,较难调整。数量少的话随机性太大。多的话很耗时。
等换了新电脑再好好试试。
knn中k的取值,是取周围5个光子,10个还是100个500个?
没有明确的说法。一般是说越大越准。但有一系列问题。
比如墙角,如果光子密度小,又要强行取一批光子,那么可能越过墙角,明显是错的。
也就是光子取的太多,范围可能太大,中间可能有遮挡。
取的少么,又可能不够准。
最好就是发射巨量的光子,使得每个交点周围的很小范围内都能找到需要的光子。

它这个渲染本质上是看光子的密度。那么限定k和限定r其实是一回事?


资料#

[Realistic Image Synthesis Using Photon Mapping]

https://www.scratchapixel.com/lessons/3d-basic-rendering/volume-rendering-for-developers/volume-rendering-summary-equations

这个资料总的很有帮助,有较多推导。但是有几个符号和字母小错,而且有的推导貌似不严密。

论文学习#

`Progressive Photon Mapping <http://graphics.ucsd.edu/~henrik/papers/progressive_photon_mapping/progressive_photon_mapping.pdf>`__

基本光子映射的一个问题是需要大量的光子来达到较高的精确度,相应地需要大量的memory。

ppm是一个多pass算法,可以不断迭代以增加准确度,同时不必维护每次迭代的photon map。从而不必无限制地增加memory用量。

基本光子映射的\(L_r\)

\({\Large L_r(x,\omega) \approx \frac{1}{\pi r^2} \sum_{p=1}^{n}f_r(x,\omega_p', \omega) \Delta \Phi_p(x, \omega_p)}\)

这只是一个近似。只有当n趋近于无穷大时,才会收敛到一个正确的无偏结果。

ppm算法第一遍是ray tracing,后续每个pass都是photon tracing。
每次photo tracing都增加了gi的准确度。所以称为progressive。

ray tracing pass#

和蒙特卡洛光追类似,从眼睛向每个pixel发光,直到碰到第一个非specular面。

对于每个hitpoint,存下列信息。

  • 交点位置

  • normal

  • 入射光方向

  • brdf

  • pixel坐标

  • pixel权重

  • 当前统计半径

  • 光子累积数量

  • 累积能量

photon tracing pass#

任务是对ray tracing的结果数据中的能量进行累积。

每一遍可以追踪一定数量的光子,建立photon map。

每一遍结束时,遍历所有交点,找到交点周围的光子。
用这些光子更新能量。
这样形成一个不断完善gi的过程。直到图像达到满意的质量。
上一遍用过的光子以后就不需要了,可以扔掉。
每一遍过后都可以更新图像。

能量的估算#

基本光子映射的\(L_r\)

\({\Large L_r(x,\omega) \approx \frac{1}{\pi r^2} \sum_{p=1}^{n}f_r(x,\omega_p', \omega) \Delta \Phi_p(x, \omega_p)}\)

依赖于光子的密度

\({\Large d(x)=\frac{n}{\pi r^2}}\)

其中假定n个光子分布在交点所处的微小disk面上,半径为r。

如果再生成一个新的map,又可以算一个新的密度

\({\Large d_{new}(x)=\frac{n_{new}}{\pi r^2}}\)

r不变。两个密度平均一下能得到更准确的密度。
要注意这个方法并不是consistent的,且不会准确收敛在x上。
它只是在r上算平均值。没法体现r之内的细节。且准确度依赖于每一遍的map。
按这个方法处理一系列photon map就会收敛到正确结果。
其中的一个关键点是,减小估算半径,同时累加光子
半径消减#

给x定义一个半径\({\large R(x)}\)。目标是增加其中的光子数\({\large N(x)}\)同时减小半径。

假设已有\({\large N(x)}\)个光子。此时经过一个pass又找到\({\large M(x)}\)个光子。可计算新密度

\({\Large d'(x)=\frac{N(x)+M(x)}{\pi R(x)^2}}\)

接着要缩小半径

\({\Large R'(x) = R(x) - R_{cut}}\)

假设\({\large R(x)}\)中的光子密度处处相等。那么新的半径中的光子数量为

\({\Large N'(x) =\pi R'(x)^2d'(x) = \pi (R(x) - R_{cut})^2d'(x) }\)

另一方面,为了结果能收敛,每一遍我们必须增加光子。
引入\(\alpha\)来控制实际的加入量。可以随便自定这个量,并不是定死M
M是计算密度改变这个事实,而\(\alpha\)是用于估算。

\({\Large N'(x) = N(x) + \alpha M(x)}\)

几个方程联立可解出\(R_{cut}\)\(R'(x)\)

\({\Large \pi (R(x) - R_{cut})^2 \frac{N(x)+M(x)}{\pi R(x)^2} = N(x) + \alpha M(x)}\)

\({\Large R(x) - R_{cut} = R(x) \sqrt{\frac{N(x)+\alpha M(x)}{N(x)+ M(x)}} = R'(x)}\)

这样就可以用\(\alpha\)来控制每次的光子添加量和半径的缩减量。
直接指定半径缩减量,反推\(\alpha\)也是可以的。用\(\alpha\)更方便。

注意新的半径\(R'(x)\)对于每个交点都要算。

功率的矫正#
设某交点现有光子的功率为\(F_N\),一遍过后新增功率\(F_M\)
做半径缩减后,新的功率简单按面积的比例计算即可

\({\Large F'_N = Flux_{new} = (F_N+F_M)\frac{\pi R'(x)^2}{\pi R(x)^2}}\)

\({\Large = (F_N+F_M)\frac{ (R(x) \sqrt{\frac{N(x)+\alpha M(x)}{N(x)+ M(x)}})^2}{ R(x)^2}}\)

\({\Large = (F_N+F_M) \frac{N(x)+\alpha M(x)}{N(x)+ M(x)}}\)

辐射率估算#

基本光子映射方程

\({\Large L_r(x,\omega) \approx \frac{1}{\pi r^2} \sum_{p=1}^{n}f_r(x,\omega_p', \omega) \Delta \Phi_p(x, \omega_p)}\)

其中求和即与\(F\)等价。

\({\Large L_r(x,\omega) \approx \frac{1}{\pi R(x)^2} F_N(x)}\)

可以看出本质就是把执行一次的大量光子计算(pm),转化为了一遍遍增加光子,同时增强画质(ppm)。

最后要normalize

\({\Large L_r(x,\omega) \approx \frac{F_N(x)}{\pi R(x)^2 \times 总光子数} }\)

总结#

与pm有不少区别。

每个hitpoint设定一个初始r值。
每次添加光子,比如2000个。计算每个hitpoint在r范围内新增的光子数。
如果这样那么knn就是不必要的?因为必须算出r范围内所有新增光子,而不是自定义取最近的k个光子。

论文学习#

`Stochastic Progressive Photon Mapping <http://graphics.ucsd.edu/~henrik/papers/sppm/stochastic_progressive_photon_mapping.pdf>`__

ppm#

ppm是一个多pass算法,用统计的方法累积光子,来解决渲染方程。

radiance可以写成

\({\Large L_r(x,\omega) \approx \frac{F_i(x)}{N_e(i) \pi R_i(x)^2} }\)

其中的i为当前pass数,每个pass中F,N,R都得到更新,结果变得更好。
当i趋近于无穷大时L达到准确值(收敛)。

\({\Large L(x,\omega) = \lim_{i \to \infty} \frac{F_i(x)}{N_e(i) \pi R_i(x)^2} }\)

sppm#

sppm的动机是,计算一块区域内的准确平均radiance。
比如motion blur需要算平均值。depth-of-field也需要算平均值。