A Recap on Volumetric Path Tracing

半年前退坑的我大概也想不到终于有一天我会遗忘这些基础技艺并且会重返$\text{Rendering}$

$\text{0. 闲话：重拾Volumetric Rendering的一些碎碎念}$

这学期选了图形学课，自然而然也要准备大作业。我的想法自然是结合曾经学过的那些东西，为自己的$\text{Rendering}$之路画上一个句号。这其中就包括那个我没有爬上去的门槛——双向路径追踪，半年前彻底地读了$\text{Veach}$写的理论与实现，但是由于各种原因，最后自己并没有实现。

考虑到双向路径追踪其实做的是对路径空间的操作，在这个理论框架下体积渲染和表面渲染是可以统一的，所以这个学期我最终考虑实现双向的体积路径追踪，一方面填上$\text{BDPT}$这个坑，另一方面也把学过的体积路径追踪集成进来，不过目前来看，这个目标真的遥遥无期。。。

到现在，如果在中文互联网上搜索“体渲染”这个关键词，自然会搜出一大堆$\text{NeRF}$相关的讲解，然而$\text{NeRF}$所谓的“体渲染”堪称粗糙，在这个语境下作出了大量简化并且散射也直接被忽略了。。。这样自然也就略去了参与介质的很多有趣的光学性质，多少还是很遗憾。这种体渲染能叫体渲染吗（暴论）

$\text{Update: } 失败啦，没做双向$

$\text{1. 光传输的算子表述}$

在初读$\text{A Null-Scattering Path Integral Formulation of Light Transport}$时，这篇论文里的形式化表述非常让我恼火写的什么玩意儿。

现在再来看的时候发现这篇神棍文默认你读过$\text{Veach}$的博士论文，因为太过有名这里就不说叫啥名字了相信能点开这篇文章的人不会一无所知，并且采用了那里面的形式化表述。所以首先得回去阅读这篇开山之作对此有一个完全的理解。

一个小插曲是之前读的$\text{BDPT}$那篇文章其实是$\text{Veach}$的论文的其中一章，不得不佩服上古时代的肝帝，那些年头的人真的是在拿应用数学对待$\text{Rendering}$。

所以为了方便后面的学习，在这里好好学习一下$\text{Veach}$中的算子表述。

$\text{1.1. 光线空间}$

在基于辐射度学的渲染框架中，我们认为$\text{Radiance}$沿着光线传播，最终形成稳定分布的光场，而实际上$\text{Radiance}$（以下用中文统称为辐射度）这个物理量可以被看成是光路的函数。实际上光路还可以承载很多的物理量以及数学量，为了方便表述，我们首先形式化地定义光路的空间。

定义光线空间$\mathscr{R=M}\times \mathscr{S^2}$，其中$\mathscr{M}$为整个场景的表面（这篇文章中尚未考虑到参与介质），$\mathscr{S}$则是球面，这还是非常好理解的，于是一条光线就可以定义为$\bold{r}(\bold{x}, \omega)$，一个表面的点和一个方向。

除此以外，也可以有其他的表示，例如我们也可以定义这样的路径空间:

$$\mathscr{R=M}^2$$

这样就是用光线的起点和终点来表示。

$\text{1.2. 通量测度}$

接下来我们定义空间$\mathscr{R}$上的测度$\mu$,这被称为通量测度,这个自然就非常直观，就是面积的测度乘上立体角的测度虽然背后的数学严谨性我不清楚，但是谁爱关心谁关心去。

$$d\mu(\bold{r}) = dA(\bold{x})d\sigma_\bold{x}^{\perp}(\omega)$$

于是我们就可以在这个框架下定义辐射度，当然这是辐射度的的定义，但是我们有了一个紧凑的表示

$$L(\bold{r}) = \frac{d\Phi(\bold{r})}{d\mu(\bold{r})}$$

回想上面说的，光线空间也可以用$\mathscr{M}^2$进行表示，在这种情况下，通量测度就定义为

$$d\mu(\bold{r}(\bold{x}, \bold{y})) = V(\bold{x}, \bold{y})\frac{\cos(\theta_\bold{x})\cos(\theta_\bold{y})}{||\bold{x}-\bold{y}||^2}dA(\bold{x})dA(\bold{y})$$

这也就是我们经常在路径追踪代码中看到的测度转换的细节。

$\text{1.3. 光线空间的函数}$

这个没啥好说的，我们可以在空间$\mathscr{R}$上定义函数，例如辐射度$L(\bold{r})$，后面还会引入重要度$W(\bold{r})$。

这些自然还可以定义其范数，我们不会像数学系的那些同学们一样追求各种稀奇古怪的反例，我们只考虑对于某个给定的$p$范数有界的那些函数构成的函数空间$L_p(\mathscr{R})$。可以证明我不会，这个空间是完备的（也就是说所有的柯西列都收敛），因此是一个巴拿赫空间。这个函数空间上自然还可以定义内积等等概念。

好吧，我承认上面这些扯淡都是为了强行拔高本文的姿势水平

$\text{1.4. 散射与传播算子}$

函数空间的存在允许我们对于光线定义物理量了，接下来要操作这些物理量自然就要定义函数空间上的算子。

散射算子。定义如下算子$\bold{K}$为散射算子：

$$(\bold{K}h)(\bold{x}, \omega_o) = \int_{\mathscr{S}^2}f_s(\bold{x},\omega_i\rightarrow \omega_o)h(\bold{r})d\sigma_\bold{x}^\perp(\omega_i)$$

当然，$\text{Veach}$的这篇文章里只考虑了表面渲染，因此只需要在同一个点对于不同的立体角积分。

于是这样我们就可以用$L_o = \bold{K}L_i$这样的式子十分紧凑地表达光的散射。

传播算子。$\text{Veach}$的论文里用了一段来形式化：

首先定义$\bold{x}_\mathscr{M}(\bold{r}(\bold{x},\omega))$为沿着$\bold{x}$点$\omega$方向出发第一次遇到的表面$\mathscr{M}$方向上的点，当然如果遇不到那就是没定义。原文把这个定义也形式化了但是我觉得真没那个必要。

然后可以定义传播算子$\bold{G}$如下：

$$(\bold{G}h)(\bold{x},\omega)=h(\bold{x}_\mathscr{M}(\bold{x},\omega), -\omega)$$

当然，这个算子也很好理解，就是描述了物理量从一个点传播到另一个点，在体渲染中这个算子还可以把衰减项也考虑进来。

这些定义并不能算是$\text{Veach}$的原创，$\text{Arvo}$在$1995$年就提出过类似的记号，和$\text{Veach}$的略有不同。$\text{Veach}$说这种定义的好处在于，当$f_s$是对称函数的时候，两个算子都是自伴算子，而在$\text{Arvo}$的定义中没有这种好性质，因此他不得不引入另一个算子$\bold{H}$,考虑$\bold{H}$以后就是等价的的定义了。

算子自伴意味着如果我们把光路逆过来也会得到相同的物理量。

在上述两个算子的基础上就可以定义传输算子

$$\bold{T}=\bold{K}\bold{G}$$

如果把场景中的所有发光体发出的辐射度定义为函数$L_e$，我们熟悉的表面渲染的渲染方程就可以在上述记号下简洁地写为算子方程

$$L = L_e+\bold{T}L$$

更多分析学的讨论这里就不继续了，这个方程并不总是有解的，我们只需要考虑有解的情形就行——实际上，对于物理世界里正常的场景，都是有解的，这个时候我们可以定义其解算子$\bold{S}$为

$$\bold{S}=(\bold{I}-\bold{T})^{-1}$$

于是就有$L=\bold{S}L_e$，$L$被称为平衡辐射度函数。

$\text{1.5. 传感器和测量方程}$

我们输出图片的方式自然是像素，但是拿像素来形式化显然不够本质，毕竟我们可以用不同的形式来表述一张图片。如果使用光线追踪/光栅化等等方法，那么我们自然可以把图像表示为一个颜色序列表示每个像素的颜色。而如果我们使用有限元，那么输出图片就是一个系数组合表示每个基函数对应的系数。当然，我们渲染出来的东西甚至也可以不是颜色。

更一般的，对于线性的传感器，我们可以定义“响应函数”如下

$$W_e(\bold{x}, \omega) = \frac{dS(\bold{x}, \omega)}{d\Phi(\bold{x}, \omega)}$$

$S$是传感器观测到的物理量，在图形学中，$W_e$被称为重要度函数，$W_e$大多数时候都是$0$，因为我们考虑的是正常的相机，大部分时候只关注一个平面上收到的辐射通量。这里的$W_e$考虑了你可能在相机上玩的所有花活儿——透镜系统，滤波器等等都在这个框架之下。

如何描述一张图片呢？这就是测量方程：

$$I(\bold{x},\omega)=\int_\mathscr{R}W_e(\bold{r})L_i(\bold{r})d\mu(\bold{r})$$

也可以通过点积紧凑地写为

$$I=<W_e, L_i>$$

这是一般性的方程，对整个光线空间积分，实际上因为很多情况下我们只关注胶片上的图片长什么样子，只需要算出胶片上的$I$即可。

$L_i$和$L$的关系可以通过传输算子联系起来（让辐射度传输到相机上），于是我们可以把测量方程写为

$$I=<W_e, L_i>=<W_e, \bold{G}L>=<W_e, \bold{GS}L_e>$$

$\text{1.6. 通过伴随算子的重要度传输}$

回忆一下，对于一个算子$\bold{H}$，我们可以定义其伴随算子$\bold{H}^*$为

$$<\bold{H}^*f, g>=<f, \bold{H}g>$$

使用伴随算子的记号，又可以把测量方程写为

$$I=<(\bold{GS})^*W_e, L_e>$$

当然，实际使用我们还是得把算子展开，这也就需要把$\bold{G}$和$\bold{K}$的伴随算子写出来，$\bold{G}$在任何场景下都是自伴的，但是$\bold{K}$仅在场景中$\text{BSDF}$均对称时才是自伴的，我们首先定义伴随的$\text{BSDF}$。

$$f_s^*(\bold{x}, \omega_o\rightarrow \omega_i)=f_s(\bold{x}, \omega_i\rightarrow \omega_o)$$

然后我们可以定义$\bold{K}$的伴随算子

$$(\bold{K}^*h)(\bold{x}, \omega_o)=\int_{\mathscr{S}^2}f_s^*(\bold{x},\omega_i\rightarrow \omega_o)h(\bold{x}, \omega_i)d\sigma_\bold{x}^\perp(\omega_i)$$

然后论文里很不负责任地假设$\text{BSDF}$都是对称的。。。于是$\bold{K}$就是自伴的了。

在这种情况下可以证明$\bold{GS}$也是自伴算子，我们可以把测量方程里的伴随去掉，于是就有

$$I=<(\bold{GS})^*W_e, L_e>=<\bold{GS}W_e, L_e>$$

所以前面的讨论约等于放了一发空炮，但是这也是一个非常重要的提醒：一定要小心场景中的不对称性，对于不对称的$\text{BSDF}$要使用其伴随。这种情况并不罕见，实际上如果稍微写过一点儿表面材料就知道透明绝缘介质的折射会导致这种情况。

由于$\bold{GS}$是自伴的（在我们美好的假设里），这个时候优雅的对称性就允许我们跟光传输一样一样地写出重要度传输，也就是。

$$W=W_e+\bold{T}W$$

$$W=\bold{S}W_e$$

这也就是为什么很多时候我们从相机出发计算路径追踪也是对的。

不过$\bold{K}$不自伴也不是什么难事，将$\bold{K}$替换为$\bold{K}^*$定义$\bold{T}_W$然后把上面的$\bold{T}$替换为$\bold{T}_W$即可。

最后的最后，我们重新回顾四个重要的物理量，入射辐射度$L_i$，出射辐射度$L_o$，入射重要度$W_i$，出射重要度$W_e$，这四个量之间的关系可以用下面的式子表示：

$$L_i=\bold{G}L_o$$

$$W_i=\bold{G}W_o$$

$$L_o=\bold{K}L_i$$

$$W_o=\bold{K^*}W_i$$

对于整个场景，大多数时候我们最初始能拿到的条件就是$W_e$和$L_e$。少数时候用入射的$W_e$和$L_e$考虑则更为自然，论文里说这通常是在使用有限元法解渲染方程的时候，不过考虑到本文的目的，这里就不展开了。

扯句闲话就是，这一套理论和物理学与工程学中的许多问题也是相关的，很多问题和光传输也有着类似的性质，这一套框架部分也是从那些领域中来的，甚至可以追溯至上世纪七十年代的中子传输问题。某种意义上可以说图形学的老祖宗还真就是计算物理

$\text{2. 体积渲染的算子表述}$

$\text{2.1. Recap: 体积辐射传输}$

接下来时间快进到$2019$年$\text{null-scattering}$的论文。体积渲染相比于表面渲染的本质区别在于光线空间可以取遍整个场景，不再局限于表面。

体积辐射传输方程在此就不再叙述了，这些属于基础中的基础，可以参考$\text{CSE-272}$作业笔记。

考虑所谓的“空散射”之后，体积辐射传输方程可以写为

$$(\omega\cdot\nabla)L(\bold{x},\omega)=-\bar{\mu}_tL(\bold{x},\omega)+\bar{\mu}_tL^{\bar{m}}(\bold{x},\omega)$$

其中

$$L^{\bar{m}}(\bold{x},\omega)=\frac{\mu_a(\bold{x})}{\bar{\mu}_t(\bold{x})}L_e^m(\bold{x},\omega)+\frac{\mu_s(\bold{x})}{\bar{\mu}_t(\bold{x})}\int_{\mathscr{S}^2}\rho_m(\bold{x},\omega'\rightarrow\omega)L(\bold{x},\omega')d\omega'+\frac{\mu_n(\bold{x})}{\bar{\mu}_t(\bold{x})}L(\bold{x},\omega)$$

上式中第一项为发光项，第二项为实散射项，第三项为空散射项。

第一个方程还可以进一步$\text{ODE}$得到

$$L(\bold{x},\omega)=\int_0^z\bar{T}(\bold{x},\bold{y})\bar{\mu}_t(\bold{y})L(\bold{y}, \omega)dy+\bar{T}(\bold{x},\bold{z})L^s(\bold{z},\omega)$$

$L^s$顾名思义就是散射的辐射度。在本文中，所有由于考虑空散射而造成的“假想”的物理量都有一个上划线，可以据此区分。$\bar{T}$就是考虑空散射后的传输项。

$$T(\bold{x},\bold{y})=\exp(-\bar{\tau}(\bold{x},\bold{y}))=\exp(-\int_0^y\bar{\mu}_t(\bold{x}-s\omega)ds)$$

$\text{2.2. 路径空间}$

你当然需要注意路径空间和光线空间的区别。此外，在本文中，为了方便表达$\delta$分布，还引入了$\text{Dirac}$测度$\delta_\bold{x}$。

定义路径空间为$$\mathscr{P}=\bigcup\limits_{k=0}^\infty (\mathscr{A}\cup\mathscr{V}\cup\mathscr{V}_\delta)^k$$

这个非常直观，一共有三种路径节点：在表面上，在参与介质中且改变了传播方向，在参与介质中并且没有改变传播方向。

$$d\bar{\bold{x}}=\prod\limits_i d\bold{x}_i$$

这里$d\bar{\bold{x}}$表示路径的测度，$\bold{x}_i$表示路径上的第$i$个节点，$d\bold{x}_i$则需要根据节点类型进行不同处理。

这里唯一较为特殊的地方大概是$\mathscr{V}_\delta$中的节点，其测度定义为

$$d\bold{x}_i=d\delta_{\bold{x}^-_{i}\leftrightarrow\bold{x}^+_i}(\bold{x}_i)$$

这很好理解，$\text{Dirac}$测度意味着我们把积分限制在了这一个方向上。

接下来，为了统一介质与表面，我们引入一般化的函数，这些函数根据点的不同类型取不同的值，有了一致的函数表达以后就可以使用算子表述一统天下了。

定义：

$$G(\bold{x},\bold{y})=\frac{D(\bold{x},\omega_{\bold{xy}})V(\bold{x},\bold{y})D(\bold{y},\omega_{\bold{yx}})}{||\bold{x}-\bold{y}||^2}$$

其中对于表面上的点$D(\bold{x}, \omega)=|n(\bold{x})\cdot \omega|$而对于体积上的点$D(\bold{x}, \omega)=1$。

这个$G$函数说白了就是用来将立体角测度转化为节点测度的。

对于表面上的点，$$L_e(\bold{x},\omega)$$就是发光的辐射度，而对于参与介质上的点，这个值还要乘上一个吸收系数$\mu_a(\bold{x})$，不过我好奇介质自己发光自己也会吸收么。。。

然后是一般性的散射函数（这个名字我自己起的）$\rho$，对于表面上的点，就是$\text{BSDF}$，对于参与介质上的点，就是$\text{Phase Function}$乘上散射系数（乘散射系数这真的是个写$\text{VPT}$特别容易翻车的点），而对于$\mathscr{V}_\delta$上的点，$\rho$的值是$\mu_n(\bold{x})H(\omega\cdot\omega')$，$H$是阶跃函数。

使用上面的表达就可以改写体积渲染方程：

$$L(\bold{x},\omega)=\int_0^z\bar{T}(\bold{x},\bold{y})L_o(\bold{y}, \omega)dy+\bar{T}(\bold{x},\bold{z})L_o(\bold{z},\omega)$$

$$L_o(\bold{x},\omega)=L_e(\bold{x},\omega)+\int_{\mathscr{S}^2}\rho(\bold{x},\omega'\rightarrow\omega)L(\bold{x},\omega')D(\bold{x},\omega')d\omega'+\mu_n(\bold{x})L(\bold{x},\omega)$$

$\text{2.3. 算子表述}$

这就是和$1.$类似的手法了，定义如下四个算子

$$(\mathbf{R}_\text{m}h)(\mathbf{x},\omega)=\int_{\mathcal{S}^2}\int_0^z\rho(\omega,\mathbf{x},\omega^{\prime})\overline{T}(\mathbf{x},\mathbf{y})D(\mathbf{x},\omega^{\prime})h(\mathbf{y},\omega^{\prime})\mathrm{d}y\mathrm{d}\omega^{\prime}$$

$$(\mathbf{R_s}h)(\mathbf{x},\omega)=\int_{\mathcal{S}^2}\rho(\omega,\mathbf{x},\omega^{\prime})\overline{T}(\mathbf{x},\mathbf{z})D(\mathbf{x},\omega^{\prime})h(\mathbf{z},\omega^{\prime})\mathrm{d}\omega^{\prime}$$

$$(\mathbf{N_m}h)(\mathbf{x},\omega)=\mu_\mathbf{n}(\mathbf{x})\int_0^z\overline{T}(\mathbf{x},\mathbf{y})h(\mathbf{y},\omega)\mathrm{d}y$$

$$(\mathbf{N_s}h)(\mathbf{x},\omega)=\mu_\mathbf{n}(\mathbf{x})\bar{T}(\bold{x},\bold{z})h(\mathbf{z},\omega)$$

于是可以定义$\mathbf{R}=\mathbf{R_m}+\mathbf{R_s},\mathbf{N}=\mathbf{N_m}+\mathbf{N_s} $

$\mathbf{R}$表示“实”，而$\mathbf{N}$表示“空”。

这样一来，把前面两个渲染方程中的$L$干掉，留下$L_o$，其积分方程就写成如下紧凑的算子形式

$$L_o=L_e+(\mathbf{R}+\mathbf{N})L_o$$

同样定义测量函数，在本文中，就默认使用像素表达了。