体传输

Date: 2018-11-06 12:41
Categories: 渲染

草稿，未完…
【x*】代表相对重要可以优先看的参考文献。

基础

直观认识

吸收：黑色的烟
散射（外散射，内散射）：云
发射：爆炸产生的红色发光烟雾

理论

对于一个微分距离：
$$\frac{dL(x,\omega )}{dz}=-{\mu }_{t}L(x,\omega )+{\mu }_a(x)L_e(x,\omega )+{\mu }_s(x)L_s(x,\omega )$$

其中${\mu }_t(x)={\mu }_a(x)+{\mu }_s(x)$，前半截是lossed，后半截是gains.

注意：只有两个系数，吸收和散射系数，吸收分为自吸收($-{\mu }_{a}L$)和发射吸收(${\mu }_a{L}_e$)，散射分为自外散射($-{\mu }_{s}L$)和外来内散射(外来光线通过phase函数:${\mu }_s{L}_s$)。

所以要确定一个介质只需要确定三个参数：吸收系数，散射系数，phase函数（包括g）.

把这个式子换成参数式得到一个一阶线性常微分方程，简单写为：

$$L^{\prime }(s)+{\mu }_t(s)L(s)=Z(s)$$

从0积分到z(所求的x点)有：

$$\frac{dL(s)}{ds}=-{\mu }_t(s)L(s)$$
两端对s积分从0到z(定积分完成后就不再有s)：
$${\int }_0^z\frac{dL(s)}{L(s)}={\int }_0^z-{\mu }_t(s)ds$$
$$\to L(z)=C{e}^{-{\int }_0^z{\mu }_t(s)ds}$$
变C为$u(z)$有（）：
$$L(z)=u(z)e^{-{\int }_0^z{\mu }_t(s)ds}$$
两端对z求导：
$$\frac{dL(z)}{dz}=\frac{du(z)}{dz}{e}^{-{\int }_0^z{\mu }_t(s)ds}-u(z){\mu }_t(z)e^{-{\int }_0^z{\mu }_t(s)ds}$$
注意到：
$$u(z){\mu }_t(z)e^{-{\int }_0^z{\mu }_t(s)ds}={\mu }_t(z)L(z)$$
全部带回原方程并改变一下变量名字（方便之后定积分）得到：
$$\frac{du(s)}{ds}=Z(s)e^{{\int }_0^s{\mu }_t(r)dr}$$
移项后对两端做0到z的定积分有：
$${\int }_0^{z}du(s)={\int }_0^{z}Z(s)e^{{\int }_0^s{\mu }_t(r)dr}ds$$
$$u(z)-u(0)={\int }_0^{z}Z(s)e^{{\int }_0^s{\mu }_t(r)dr}ds$$
带回得到：
$$L(z)=e^{-{\int }_0^z{\mu }_t(r)dr}[u(0)+{\int }_0^{z}Z(s)e^{{\int }_0^s{\mu }_t(r)dr}ds]$$
通过()计算可知：$L(0)=u(0)$，于是：
得到最终解：
$$L(z)=e^{-{\int }_0^z{\mu }_t(r)dr}[L(0)+{\int }_0^{z}Z(s)e^{{\int }_0^s{\mu }_t(r)dr}ds]$$

解微分方程的时候一般有个目标的上下限，如果确定使用定积分解方程，所有积分步骤都要使用定积分（这一步会更替自变量），关键问题在于：原本微分式子，做定积分后原本的自变量失效（因为定积分是常数！），现在的自变量替换为积分上或下限。而如果使用不定积分来解决问题则一切照旧。

上面的积分形式称为体渲染方程(VRE)其中$e^{-{\int }_0^z{\mu }_t(r)dr}L(0)$称为背景辐射(background radiance)。
直观地讲：不仅线上每个点都要辐射，而且每个点都要对应地衰减。

蒙特卡洛体传输

三阶矩代表对称度，0代表分布对称，正数代表左倾斜，负数代表右倾斜
四阶矩越大越拖尾(-<+)

Diffusion Approximation

Delta-Eddington Approximation【1】结合一个delta函数和一个two-term近似解决了Eddington近似在高度不对称情况下的poor accuracy. 同时这个近似取消了HG的分母。

DA的推导【2】【3】为了简化使用Delta-Eddington Approximation带入辐射度微分方程的phase函数，把delta部分积分出来与展开的衰减系数合并得到一个与原来形式近似的辐射度方程。令：

$${\mu }_s^{\prime }={\mu }_s(1-f)$$
$${\mu }_t^{\prime }={\mu }_a+{\mu }_s^{\prime }$$

整理一下得到（注意这里提前假设source term是0）：

$$\frac{dL(r,\omega )}{dz}=-{\mu }_t^{\prime }L(r,\omega )+\frac{{\mu }_s^{\prime }}{4\pi }{\int }_{4\pi }L(r,{\omega }^{\prime })(1+3g^{\prime }\omega \cdot {\omega }^{\prime })d{\omega }^{\prime }$$

其中${\mu }_s^{\prime }$和${\mu }_t^{\prime }={\mu }_a+{\mu }_s^{\prime }$分别叫做reduced scattering coefficient和reduced extinction coefficient 。 不散射的部分被$1-f$reduced了, and this light is scattered out of the collimated beam and into the diffuse portion of the radiance。Similarity Theory指出，对于一个改变的参数集，可以得到RTE近似一样的解（ By the principle of similarity, we use these reduced coefficients and assume isotropic scattering without significantly altering the solution【6】）。这个可以解释为：
Physically, we are lengthening or shortening the distance traveled by scattered light by 1− f, which after many scattering events is similar to the effect of an anisotropic phase function.
参考【7】。

Similarity Theory : 研究条件相似与解相似关系的理论，相似三角形是一个几何应用。RTE上的应用可以参考【4】【5】【6*】The authors derived a set of relations between two sets of scattering parameters so that the resulting RTEs have identical solution radiance fields when their directional frequencies are bounded. A highly simplified order-1 form of this theory is used extensively in diffusion methods and has been applied for accelerating Monte Carlo simulation of light transport.

delta-Eddington函数第一项代表前方的strong散射，第二项则代表更加漫反射type的散射。f代表前向峰值程度，g’则代表漫反射部分的对称性。
radiance分为collimated和diffuse两部分，前者代表入射beam方向还未散射的radiance，后者代表其他方向散射到beam方向的漫反射radiance：

$$L=L_{coll}+L_d$$

单独考察$L_{coll}$在Reduced Equation中的衰减，但是注意直接使用系数${\mu }_t$衰减会比Reduced情形下衰减的更多，因为Reduced Equation情形下的衰减有漫反射补偿（delta-Eddington函数的delta部分的加入补偿了原来的衰减），所以应该使用${\mu }_t^{\prime }$来进行衰减，于是有：
$$\frac{d{L}_{coll}}{ds}=-{\mu }_t^{\prime }{L}_{coll}$$

再注意微分方程的初始条件为：$L_{coll}$介质边界来的一道beam，应为菲涅尔穿越介质边界时会被减弱，于是：
$$L_{coll}=T\pi F_0(r)\delta (\omega \cdot z-{\omega }_0\cdot z)$$

于是解得：
$$L_{coll}=T\pi F_0(r)\delta (\omega \cdot z-{\omega }_0\cdot z)e^{-{\mu }_t^{\prime }z/{\mu }_0}$$

这些方程全部会带到辐射度方程化简移项得到：
$$\frac{d{L}_d}{ds}=-{\mu }_t^{\prime }{L}_d+\frac{{\mu }_s^{\prime }}{4\pi }{\int }_{4\pi }{L}_d({\omega }^{\prime })[1+3g^{\prime }(\omega \cdot {\omega }^{\prime })]d{\omega }^{\prime }+$$
$$\frac{{\mu }_s^{\prime }}{4\pi }T\pi F_0(r)e^{-{\mu }_t^{\prime }z/{\mu }_0}[1+3g^{\prime }(\omega \cdot z)({\omega }_0\cdot z)]$$

如果假设phase函数对称并使用算子来表示方向导数则上述公式进一步简化为：
$$(\omega \cdot ▽)L_d(r,\omega )=-{\sigma }_t^{\prime }{L}_d(r,\omega )+\frac{{\sigma }_s^{\prime }}{4\pi }{\int }_{4\pi }{L}_d(r,{\omega }^{\prime })d{\omega }^{\prime }+\frac{{\sigma }_s^{\prime }}{4\pi }T\Phi (r)e^{{\sigma }_t^{\prime }z/{\mu }_0}$$

为了使得解这个方程更加容易，下面需要对$L_d$做近似展开，此处可以看做是球谐展开【2 P37】（实际上此处也是以1,2阶矩来表示函数）：
$$L_d(r,\omega )\approx \frac{1}{4\pi }\varphi (r)+\frac{3}{4\pi }E(r)\cdot \omega$$

讲展开代入化简后的方程，并注意到：
$$\varphi (r)={\int }_{4\pi }{L}_d(r,\omega )d\omega$$
$$E(r)={\int }_{4\pi }{L}_d(r,\omega )\omega d\omega$$

注意上式被积函数不是点乘，所以结果是一个向量。
第一个式子叫做fluence，代表各向同性flux，可以考虑为光入射一个medium中小体积的平均量；第二个叫做vector irradiance或者vector flux，代表某方向的各向异性flux，代表某一点的净光子流（net photon current）.

讲得到的结果两端做全球积分，以及注意：

$${\int }_{4\pi }\omega \cdot Ad\omega =0$$
$${\int }_{4\pi }(\omega \cdot A)(\omega \cdot B)=\frac{4\pi }{3}(A\cdot B)$$

算子可以当成一个独立的向量进行计算

于是到的结果：
$$▽\cdot E(r)=-{\sigma }_a\varphi (r)+{\sigma }_s^{\prime }T\Phi (r)e^{{\sigma }_t^{\prime }z/{\mu }_0}$$

此表达式说明diffuse flux取决于scalar flux的吸收以及collimated光的内散射。

联系两个量的方程

如果我们把公式两边乘以一个$\omega$再做全球积分，可以得到：

$$▽\varphi (r)=-3{\sigma }_t^{\prime }E(r)$$

我们最终要得到表面的漫反射radiance，等价于找到表面法方向的radiant exitance（单位面积出射的辐射流），如果考虑到flux强度$E^{\prime }(r)=E(r)\cdot n$（乘以cos），这个量也可以考虑为射出材质的光子流，于是有：

$$E(r)=-D▽\varphi (r)$$

这意味着有净光子量流从高能量区域流向底能量区域。

Diffusion equation

如果方程两端再取散度得到：

$$▽\cdot E(r)=-\frac{1}{3{\sigma }_t^{\prime }}{▽}^2\varphi (r)$$

再讲两端直接积分的方程代入其中，消掉E®，于是可以得到classic diffusion equation：

$${▽}^2\varphi (r)-{\sigma }_{tr}^2(r)=-3{\sigma }_s^{\prime }{\sigma }_t^{\prime }T\Phi (r)e^{{\sigma }_t^{\prime }z/{\mu }_0}$$

其中${\sigma }_{tr}^2=3{\sigma }_t^{\prime }{\sigma }_a$，叫做effective transport coefficient .
注意，方程右边代表源函数，这里使用的是collimated源，如果改用其他源应该修改这个源函数。

求解Diffusion方程，讲结果代入联系方程，就能得到最终的radiant exitance .

求解的方法很多地方有写，在此不表，可以参考PBRT等 .

注意得到radiant exitance之后带回原来的方程得到$L_d$，注意可以直接使用$E\cdot n$来直接计算标量E，这和对$L_d$积分求标量E是一样的（因为展开的时候E就已经是L的一阶矩）。

Jensen的bssrdf两边（出射和入射）都上了菲涅尔，然后中间的$R_d$就不再考虑菲涅尔，而PBRT的bssrdf考虑到$L_d$只在体内部有效，出射计算radiant的时候考虑的菲涅尔，所以pbrt的bssrdf公式只有一个菲涅尔term (PBRT P927 15.30).

注意：$$n\cdot \int \omega \cdot L(\omega )d\omega =\int n\cdot \omega L(\omega )d\omega =\int cos\theta L(\omega )d\omega$$

Sampling light transport paths

Basic Methods

单项方法中，直接均匀采样光源面积，效率较低，下列技术可以uniformly 采样立体角：
Stratified sampling of spherical triangles（siggraph）
An area-preserving parametrization for spherical rectangles（EGSR 2013）
Solid angle sampling of disk and cylinder lights（EGSR 2016）
Area-preserving parameterizations for spherical ellipses（EGSR 2017）
或者采样polygon明度：
Stratified sampling of 2-manifolds（sig course 01）

In fact, for point lights the variance is unbounded as individual contribu-
tions tend to infinity when the last free vertex（也就是倒数第二个顶点） in the medium approaches (by random chance) the point light—the geometry term along the shadow ray tends to infinity in such cases （点光源那个点和倒数第二个点距离过近）.（对于一般的光源也有这个问题）
一种避免这个情况的方法就是，先采样灯光上那个点，然后根据这个点的情况来采样之前的几个点，有两种方法：
Equiangular sampling： Importance sampling techniques for path tracing in participating media（EGSR 2012）
简单来说，这种方法，现在光源上选择一个点，然后沿着$x_{k-2}$到$x_{k-1}$采样一$x_{k-1}$，采样正比于$1/d^2$.
[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-8csyQS8U-1605237030673)(./vol_1.png)]
PDF和CDF在此不表，可以查阅原论文.

对于路径积分的形式来说，总的采样概率应该正比于path的contribution，但是实践中如果要从路径积分的pdf中来采样，就要对每个随机变量边缘化他的pdf，这是无法分析性计算出来的，也是不可行的。
所以通常的采样方法总是使用local采样（phase函数等），这并没有从joint分布上来控制采样，就可能会造成高方差，特别是对于光源在高度各向异性的media中。
For the special case of double scattering in media, Georgiev et al. [GKH ∗ 13] have found a way to sample from a prescribed joint distribution proportional to the product of geometry and phase function terms along a sequence of four path vertices—three segments.
*Joint importance sampling of low-order volumetric scattering（siggraph Asia 2013）

在有色media中波长dependency往往也有很大影响。比如： when the extinction varies with wavelength, there exist a few interesting variants of importance sampling.
The straightforward way to handle such situations is to just trace a single wavelength (or one channel of a color space in RGB rendering). Doing so, however, introduces chromatic noise and is wasteful because there is usually a large amount of correlation of the path contribution between different wavelengths. We could potentially save cost by calculating more wavelength contributions for every generated path.
Hero wavelength sampling：
Hero wavelength spectral sampling（EGSR 2014）
An improved technique for full spectral rendering（WSCG 2009）

【6.5.2. Spectral tracking with null-collision methods】略

Advanced methods

Bidirectional path tracing starts random walks to form transport paths both at the sensor and the light sources. In the presence of participating media, this proceeds just the same as in the vacuum case, only that propagation distances need to be sampled to simulate interactions with the medium and transmittance must be evaluated along deterministic connections.

其中transmittance解算是对称的，它与方向无关，但是pdf却不是：the normalization depends on the extinction coefficient at the target point µ t (t).

Reference

【1】The Delta-Eddington Approximation for Radiative Flux Transfer
【2】Towards Realistic Image Synthesis of Scattering Materials
【3】A radiative transfer framework for rendering materials with anisotropic structure
【4】Similarity relations for
anisotropic scattering in monte carlo simulations of deeply penetrating neutral
particles
【5】Similarity relations for the interaction parameters in radiation transport
【6】High-Order Similarity Relations in Radiative Transfer
【7】LIGHT TRANSPORT IN TISSUE