光子映射学习#
[Realistic Image Synthesis Using Photon Mapping]这个书来学习光子映射[A Practical Guide to Global Illumination using Photon Maps]大致是其精简版，也包括一些不同信息可以参考。持续更新
caustics问题
一个重点是光子映射与场景的geo信息是解耦的。所以可以用于非常复杂的场景。
1. Photon emission#
每个光子带相等的flux(辐射功率)。
一个光源的总功率是\(P_{light}\)，那么发射\(n_e\)个光子，每个光子的功率
\({\Large P_{photon} = \frac{P_{light}}{n_e}}\)
多光源#
每个光源都要发射光子。
2. Photon tracing#
散播flux，光线追踪是收集radiance。一个例子是对于radiance，折射由折射率决定。光子不是这样算。 光子打到物体上，可能反射，传播，或吸收。概率由material的反射率决定。
反射还是吸收？#
设r为我们取的反射概率
p = randint() # [0-1]的随机数
if p < r:
以原始能量反射
else:
吸收
两种都是对的。我们就可以使用第二种解读。保持正确性的前提下，可以有效地减少计算。
反射，传播，还是吸收？#
以单色为例，设d为diffuse反射概率，s为specular反射概率。 可以按三个概率区间决定光子如何传播。
对于多色。
Photon storing#
扩展到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 mapparticipating media的间接光。\({\Large L\{S|D|V\}^+V}\)
上面用了[Heckbert90]的描述光子路径的语法。
L 从光源发射
S 镜面反射或传输
D 非镜面反射或传输
V volume散射
{x|y|z} xyz中的一个
\(x^+\) 1个或以上x
\(x^*\) 任意个x
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'}\)
需要把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')}\)
\({\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)}\)
可以用其他形状代替球体，各有好坏。 球体比较好算。
3.1 Filtering#
要减少这种模糊，估算的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#
总的输出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'}\)
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#
光子映射很适合解决这类问题。是第一个完全模拟出subsurface的算法。
4.1 Light Scattering in Participating Media#
当一个光子进入一种介质，可以不受影响继续往前走，或者在某个位置与介质交互。
当光穿过介质，可以看作是radiance的连续改变。
\({\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'}\)
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 \frac{dP}{dt}= kP，P(0) = P_0}\)
\({\Large P = P_0e^{kt}}\)
变量是时长t。但我们的L变量是任意点x，而且形式不一样，能否转成一个标准的微分方程？
\({\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}}\)
如果是均匀的介质，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}}\)
完整的情况#
\({\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}}}\)
4.4 Ray Marching#
Ray Marching。那么就是一个近似的解，细分越多越精确。
啊照这样说，那是不是根本就不用推导那一大片微分方程？每一步直接算最初的dL和更新L不就行了？
4.4.1 自适应ray matching#
上面的ray marching时平均步长。如果是非均匀的材质，或者有本地波动(比如阴影)，最好用自适应步长。 例如当材质特性发生较大变化时，可采用较小的步长，以更好地反映出细节。
方法2 根据sigma_t来调整步长。这个对于未知的材质可能很有效。
取均匀分布的随机数
貌似涉及均匀采样
4.5 Photon Tracing#
4.5.1 Photon Scattering#
两种情况 光子与介质交互，或发散或被吸收。 光子与某表面交互，反射或传输或吸收。
光子进入介质时，不看作与表面交互，而是直接往里面走，进行与介质的交互。
之前已提到sigma的意义，光子进行下一个交互前平均移动的距离。
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倍。
它这个渲染本质上是看光子的密度。那么限定k和限定r其实是一回事？
资料#
[Realistic Image Synthesis Using Photon Mapping]
论文学习#
`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趋近于无穷大时，才会收敛到一个正确的无偏结果。
ray tracing pass#
和蒙特卡洛光追类似，从眼睛向每个pixel发光，直到碰到第一个非specular面。
对于每个hitpoint，存下列信息。
交点位置
normal
入射光方向
brdf
pixel坐标
pixel权重
当前统计半径
光子累积数量
累积能量
photon tracing pass#
任务是对ray tracing的结果数据中的能量进行累积。
每一遍可以追踪一定数量的光子，建立photon map。
能量的估算#
基本光子映射的\(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}}\)
减小估算半径，同时累加光子。半径消减#
给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) }\)
可以随便自定这个量，并不是定死M。\({\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)}\)
注意新的半径\(R'(x)\)对于每个交点都要算。
功率的矫正#
\({\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有不少区别。
论文学习#
`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} }\)
\({\Large L(x,\omega) = \lim_{i \to \infty} \frac{F_i(x)}{N_e(i) \pi R_i(x)^2} }\)