【笔记】Progressive Photon Mapping渐进式光子映射

渐进式光子映射Progressive Photon Mapping

感觉有点难找到合适的教程。

直接看pbrt有点头大，于是选择了从论文开始。这里的内容主要是论文翻译

PPM相对于其他无偏的离线渲染方法（主要指PT，BPT，MLT）和一些有偏的方法（PM），解决的是SDS的问题（SDS中的caustic现象）。

光子映射Photon mapping

是一个2pass 的方法

  • 1st pass
    • photon tracing
  • 2nd pass
    • rendering using photon map

对于一个photon map，任何一个表面位置x的exitant radiance估计为
$$
L(x,\overrightarrow w)\approx\sum_{p=1}^n\frac{f_r(x,\overrightarrow w,\overrightarrow w_p)\phi_p(x_p,\overrightarrow w_p)}{\pi r^2}
$$
$n$ :邻近的光子数量，用来估计incoming radiance

$\phi_p$ :第p个光子的flux（power）

这种估计假设了局部的光子代表了x接收到的radiance,并且x周围的表面是平坦（flat）的。

这也是Photon mapping方法偏差的来源。photon tracing过程本身是无偏的，但是photon分布的结果在radiance estimate的过程中进行平均（blurred）。

光子密度（photon density）增加，radiance的估计也会收敛到正确的结果，这说明光子映射是一致的(consistent)。

为了保证收敛到正确结果，需要在photon map中储存无限的光子，并且radius半径需要收敛到0。我们可以通过一种操作满足这些要求：

在photon map 中使用$N$个光子，但是只有$N^\beta$ ($\beta\in [0,1]$)个光子用来进行radiance 估计，当$N$趋近无穷时，$N$和$N^\beta$ 都趋近无穷，但是$N^\beta$是N的高阶无穷小，保证了r收敛到0，以下会把它称为辐射度估计方程。
$$
L(x,\overrightarrow w)\approx \lim {N\to\infin}\sum{p=1}^{N^\beta}\frac{f_r(x,\overrightarrow w,\overrightarrow w_p)\phi_p(x_p,\overrightarrow w_p)}{\pi r^2}
$$
在标准的光子映射中，这个结果是理论正确的，但是光子储存在内存中，这使得无法获得精确的结果。

渐进式光子映射将解决这个问题。

渐进式光子映射Progressive Photon Mapping

刚才又去浏览了一下PBRT，我知道我为什么看起来难受了，Literate Programming是什么鬼啊……文学编程……从零开始直接跳去看PM的内容着实有点难受了。还是接着从论文开始吧。

PPM是多pass的算法，先是ray tracing，然后子序列的pass用来做photon tracing，每一次的photon tracing pass 都会提升全局光照结果的准确性。

Ray Tracing Pass

使用标准的ray tracing通过图像中的每个像素找到场景中所有可见表面，每一条ray path都包含了所有specular的反弹，直到遇到第一个non-specular的表面。

场景中specular表面比较多时，也可以用俄罗斯轮盘赌（Russian Roulette）来停止。而如果击中的表面BRDF有non-specular的部分，对于每个ray path我们都储存路径上所有的hit points（对于这句话我的理解是，如果不是完全镜面的表面，就会有能量的吸收，有一部分光子停留在这里）

1
2
3
4
5
6
7
8
9
10
11
struct  hitpoint {
    position x;//击中位置
    normal   n;//x所在的法线
    vector   w;//入射光线方向
    integer BRDF;//BRDF的index
    float    x,y;//像素位置
    color    wgt;//像素权重
    float    R;//当前的光子查找半径
    integer  N;//累计光子数
    color    t;//累计反射的flux能量
}

Photon Tracing Pass

这个步骤是用来累计光子能量的。可以分成很多的pass去做，每个pass追踪一系列光子，每个 photon tracing pass结束后，就去查看所有hit points，找到半径区域的光子。使用新加入的光子来修正光照计算。光子的贡献一经记录，就可以把光子丢掉了，然后再去处理下一个photon tracing pass。直到累积了足够数量的光子。

我们甚至可以在每个PTP后渲染一遍场景，累积的光子越多，场景的质量就越高。

Progressive Radiance Estimate渐进的辐射度估计

传统的PM算法估计着色点的局部光子密度
$$
d(x) = \frac{n}{\pi r^2}
$$
这个估计的假设是周围是平面，在一个半径为r的圆盘上去估计。如果我们在新的光子图上，要产生新的光子，在同样的圆盘上去估计密度，那就是
$$
d’(x)=\frac{n’}{\pi r^2}
$$
把$d(x)$ 和$d’(x)$ 进行平均，我们可以获得半径 r上更准确的估计。这种方法可以获得更加平滑的radiance估计，但是最终结果由于平均计算会失去很多细节。并且平均过程会破坏一致性，使得无法收敛到正确结果。

渐进式的辐射度估计结合多个光子图，能够收敛到正确结果，也能解决细节问题。关键方法是在每个hit point的辐射度估计中，随着光子数量的累计减少半径。这有效地保证了光子密度在极限估计趋于无穷。

接下来会描述光子密度是如何渐进式增长的。我们在ray tracing pass生成的每个hit points上都会计算辐射度估计。初始化时，x对应的半径R(x)会设置一个非零值，比如说对应像素的footprint（虎书中对像素footprint的解释是，屏幕空间像素映射到纹理空间的形状，这里可以看作世界空间）。也可以第一次photon tracing pass后，通过使用光子图来估计半径。

Radius Reduction半径缩减

每个hit point都有一个半径R(x)，我们的目标是，半径内的光子数量累计增加时，减少半径。

image-20220515211027302

hit point x的密度d(x)，使用上面的公式就可以算，假设已经做了一些photon tracing了，在x处累计了N(x)的光子，如果这个时候在做一个photon tracing pass，并且在R(x)范围内又找到M(x)个光子，我们可以把新的M(x)个光子加上去
$$
\hat d(x)=\frac{N(x)+M(x)}{\pi R(x)^2}
$$
下一步就是用dR(x)来减少半径R(x)，如果我们假设半径R(x)内的光子密度是常数，我们可以算出新的圆盘半径$\hat R(x)=R(x)-dR(x)$ 内光子总数
$$
\hat N(x) = \pi\hat R(x)^2\hat d(x) = \pi(R(x)-dR(x))^2\hat d(x)
$$
为了满足辐射度估计方程中的条件 ，每一次迭代的光子总数都需要有增加的($\hat N(x)>N(x)$，为了简便，使用了一个系数$\alpha \in [0,1]$ 来控制光子的比例
$$
\hat N(x) = N(x) + \alpha M(x)
$$
也就是说，我们每次迭代可以把$\alpha M(x)$ 个新的光子加上去，可以计算出对应需要减少的半径$dR(x)$
$$
\pi(R(x)-dR(x))^2\hat d(x) = \hat N(x)
\\Leftrightarrow \pi(R(x)-dR(x))^2\frac{N(x)+M(x)}{\pi R(x)^2} = N(x) + \alpha M(x)
\\Leftrightarrow dR(x) = R(x) - R(x)\sqrt{\frac{N(x) + \alpha M(x)}{N(x)+M(x)}}
$$
所以更新的$\hat R(x)$
$$
\hat R(x)=R(x)-dR(x)=R(x)\sqrt{\frac{N(x) + \alpha M(x)}{N(x)+M(x)}}
$$
注意，在每个hit point，这个公式都是独立计算的。

Flux Correction能量修正

当hit point接收到新的M(x)个光子，我们还需要加上这些光子所携带的能量。还需要把前面计算的半径缩减考虑在内。每个hit point储存接收到的BRDF预乘后未归一化的能量。把它叫做$\tau(x,\vec w)$ ,对于N(x)个光子
$$
\tau_N(x,\vec w) = \sum_{p=1}^{N(x)}f_r(x,\vec w,\vec w_p)\phi_p’(x_p,\vec w_p)
$$
$\vec w$是hit point的入射光线的方向，$\vec w_p$是入射光子的方向，$\phi_p’(x_p,\vec w_p)$ 是光子p携带的未归一化的能量。注意！这个阶段的能量在标准光子映射中，是没有被发出光子的数量除掉的。

同样的，新的M(x)个光子提供的能量
$$
\tau_M(x,\vec w) = \sum_{p=1}^{M(x)}f_r(x,\vec w,\vec w_p)\phi_p’(x_p,\vec w_p)
$$
如果半径是常数，那我们可以干嘛，直接把这两个能量加起来了。但是半径减少了，我们还要考虑到已经变成在半径外面的那些光子。

一种方法是，维护一个圆盘内所有光子的列表，半径衰减后，不在圆盘内的，就把它们移出去。 但是这个方法不实用，因为光子列表消耗太多内存了。因此，我们假设圆盘内的光照和光子密度是常数，会有以下的结果
$$
\tau_{\hat N}(x,\vec w) = (\tau_N(x,\vec w)+\tau_M(x,\vec w))\frac{\pi \hat R(x)^2}{\pi R(x)^2}
\=\tau_{N+M}(x,\vec w)\frac{\pi(R(x)\sqrt{\frac{N(x)+\alpha M(x)}{N(x)+M(x)}})^2}{\pi R(x)^2}
\=\tau_{N+M}(x,\vec w)\frac{N(x)+\alpha M(x)}{N(x)+M(x)}
$$
$\tau_{\hat N}$ 就是半径缩减后$\hat N$ 个光子相关的缩减后的能量。最开始的时候，假设的是半径内的光子密度和光照是常数，这可能不正确，但是随着半径越来越小，这个结果会变得越来越正确，除了恰好位于照明不连续处的点。但这不构成问题，因为不连续的光照是未定义的，击中点恰好在不连续的位置的概率是0。

Radiance Evaluation辐射度估计

每一次光子追踪后，我们都可以估计击中点的辐射度。回调储存的数据，包括当前半径、当前的乘以BRDF后的截断能量。估计的辐射度要乘以对应像素的权重，再加到对应像素上。

为了估计辐射度，我们还需要知道发射出的光子总数$N_{emitted}$ 用来归一化$\tau (x,\vec w)$

辐射度估计如下
$$
L(x,\vec w) = \int_{2\pi}f_r(x,\vec w,\vec w’)L(x,\vec w’)(\vec n\cdot\vec w’)dw’
\\approx\frac{1}{\Delta A}\sum_{p=1}^nf_r(x,\vec w,\vec w’)\Delta\phi_p(x_p,\vec w_p)
\=\frac{1}{\pi R(x)^2}\frac{\tau(x,\vec w)}{N_{emitted}}
$$
和正常的光子映射相似，这个公式没有限制为Lambertian材质，因为我们要把能量预先乘上BRDF再储存为$\tau(x,\vec w)$ 。

如果R（x）定义的圆盘位于未照明区域内，则半径R（x）不会减小（因为M（x）= 0 ）。虽然这种情况看起来破坏了一致性，它仍然会收敛到正确的结果$L(x,\vec w) = 0$ ，因为随着$N_{emitted}\to \infin, \tau(x,\vec w)也不会增加，L(x,\vec w)\to0$ ，

总之论文根据实验数据给出N(x)和R(x)是正确收敛的，半径衰减至0，光子数量增长至无穷，辐射度也会正确收敛。渐进式的辐射度估计保证了每次迭代中每个击中点光子密度的增加，和辐射度估计方程是一致的。

再往后就是论文对不同方法的效果的比较了，差不多到这里就结束了。