体渲染(Volume Rendering)公式推导
前言
由于这是个人对图形学领域的知识积累，可能会直接保留参考文章中深刻见解的语句，并结合自我的理解提出一些个人看法，因此本文不具有版权保护，仅作分享交流。
体渲染 (Volume Rendering)
体渲染 (Volume Rendering) 是图形学中最重要的技术之一，大名鼎鼎的 NeRF 的核心就是基于此利用神经网络优化各个参数实现的。
所谓计算机图形学，就是让计算机模拟出一个真实的世界。而渲染的作用，是将这个虚拟的世界投影成图像，正如自然界中的各种光线经过碰撞后，投影到我们视网膜上的样子。
体渲染属于整个渲染技术的分支，它的目的主要是为了解决云、烟、果冻这类非刚性物体的渲染建模，可以简单理解为是为了处理密度较小的非固体的渲染。当然进一步推广到固体的渲染也说的通，而 NeRF 也是这样做的。为了建模这种非刚性物体的渲染，体渲染把气体等物质抽象成一团飘忽不定的粒子群，这里可以从物理的角度理解，就是从分子的角度看待气体，将气体物质建模为一系列的气体分子。光线在穿过这类物体时，其实就是光子在跟粒子发生碰撞的过程。
如上图所示，光沿直线方向穿过一堆粒子 (粉色部分)，如果能计算出每根光线从最开始发射，到最终打到成像平面上的辐射强度，我们就可以渲染出投影图像（最终获取的还是经过粒子群后的最终辐射强度），而体渲染的任务就是对这个过程进行建模。为了简化计算，这里假设光子只跟它附近的粒子发生作用，这个范围就是图中圆柱体大小的区间。
关于辐射、辐射强度、渲染
- 辐射是物理学中的一个术语，指的是能量以波或是次原子粒子移动的型态，在真空或介质中传送，最熟悉的就是电磁波向四周辐射出磁场或电场，而可见光就是一种电磁波，并且具有波粒二象性，一束光也即是一系列光子向着某个方向不断前进的过程。
- 辐射强度在近代物理中指点辐射源在某方向上单位立体角内传送的辐射通量，通俗地理解可以认为是单位时间内垂直于光线传播方向的某个体积元内正在传送的光子数量，表示光子传送的分布密度。当某点传送的光子数量越多，表示其辐射强度越高。
- 渲染，计算最终辐射强度可以理解为计算到达屏幕上各点的光子密度，只不过在渲染出图像这个过程中，每个光子都包含有其他的属性，例如颜色、密度等。渲染将会在初步建立的图像上面添加位图纹理或者程序纹理、照明、凸凹纹理映射以及相对于其它物体的位置。
作用类型
体渲染把光子与粒子发生作用的过程，进一步细化为四种类型：
- 吸收 (absorption)：光子被粒子吸收，导致入射光通过该点后的辐射强度减弱；
- 放射 (emission)：粒子本身可能发光，比如气体分子加热到一定程度，会诱导原子跃迁释放光子，这进一步增大辐射强度；
- 外散射 (out-scattering)：光子在撞击到粒子后会发生碰撞，可能导致方向发生偏移，这会减弱当前方向的入射光强度；
- 内散射 (in-scattering)：其他方向的光子在撞到粒子后，可能和当前方向上的光子重合，从而增强当前光路上的辐射强度。
如下图所示的四种情况，$L_i$ 表示入射光，$L_o$ 表示出射光。
在只有这四种情况的作用时，计算出射光与入射光之间的变化量，可以认为是这四种情况的叠加，即
$$
L_o - L_i = d L(x, w) = emission + inscattering - outscattering - absorption
$$
其中 $x$ 表示光线上某个位置点， $w$ 表示光线发射方向。
吸收 (absorption)
微观分析
How to model the absorption?
如上图所示吸收的过程，我们首先要明确一点：吸收的过程是指粒子吸收了光子，使得光子不会在转移、传送，这个过程其实可以看作是光子被粒子遮挡，一旦遮挡就永远不会达到（理想遮挡），那么只需要计算遮挡比例，然后将其作为对入射光的衰减比例即可。每一点的位置对应的遮挡比例不同该怎么办？积分来解决。在表示每一点的遮挡比例后，可以将其进行积分已得到最终经过吸收后的光照强度。
假设粒子均为半径 $r$ 的球体，每个粒子的投影面积（对光线的遮挡面积）为 $A = \pi r^2$，圆柱体中粒子的体密度是 $\rho$，圆柱体底面积为 $E$。
当圆柱体足够薄 (薄到跟粒子一样厚) 的时候，可以认为粒子之间不会互相重叠 (也就是粒子都平铺在圆柱体一个横截面上)。
假定这个厚度是 $\Delta s$，那么在这个厚度内，圆柱体体积为 $E\Delta s$，包含的例子总数为 $\rho E\Delta s$。这些粒子遮挡的面积为 $\rho E\Delta s A$，占整个底面积的比例为 $\frac{\rho E\Delta s A}{E} = \rho \Delta s A$，即，当一束光通过这个圆柱体的时候，有 $\rho \Delta s A$ 的概率会被遮挡。
无限延展
当在圆柱体一段均匀发射无数光线（假设都朝向相同的方向），在另一端接收，被遮挡（吸收）的光线无法通过。根据建模原理，可以将遮挡比例看作是光线的衰减系数，即接收到的光的强度的均值是入射光的 $\rho \Delta s A$ 倍。（均值是指接收到的光线还要考虑接受面积，更确切地说是单位面积上的光强，这样才能与入射光进行比较）这可以形式化地表示为，
$$
I_o(s) - I_i(s) = \Delta I =-\rho(s) A \Delta(s)I_i(s) \tag{1} \label{eq1}
$$
这里粒子密度 $\rho, I_o(s), I_i(s)$ 都是关于位置 $s$ 的函数，因为每个区域的密度、光照强度都不同，而且微分量表示的是光线经过该位置介质后被遮挡的那部分。
等式 $\eqref{eq1}$ 可以进一步转换为常微分方程：
$$
\frac{dI}{ds} = \frac{\Delta I}{\Delta s} = -\rho(s) A I(s)=-\tau_a(s)I(s) \tag{2} \label{eq2}
$$
上式中，$I_i(s)$ 简写为了 $I(s)$ 也就是说最后解微分方程得到的是入射光的强度，并且将 $\rho(s)A$ 合并为 $\tau_a(s)$，根据定义 $\rho(s)$ 表示某个位置 $s$ 的体密度，而 $A$ 表示粒子的投影面积，因此 $\tau_a(s)$ 表示的对应于粒子投影面积的面密度，衡量了粒子对光线的吸收能力。
解方程 $\eqref{eq2}$，得到
$$
I(s) = I_0 \exp(-\int_{0}^{s} \tau_a(\sigma) d\sigma) \tag{3} \label{eq3}
$$
其中，$I_0$ 表示起始光线的强度，是常微分方程中的常数项，或者说计算的起点，$I_0 = I_i \exp(-\int_{i}^{0} \tau_a(\sigma) d\sigma)$.
物理意义
对公式 $\eqref{eq3}$ 物理意义的探讨，如果介质（粒子群）是均匀的，即 $\tau_a(t)$ 处处相等，那么解得：
$$
I(s) = I_0 \exp(-\tau_a s) \tag{4} \label{eq4}
$$
这表示入射光经过介质（粒子群）后，辐射强度会呈指数衰减，称为比尔-朗伯定律（Beer-Lambert Law）。
根据此可以定义透射比（transmittance）：出射光与入射光光强之比
$$
T(s) = \frac{I_o}{I_i} = \exp(-\int_i^o\tau_a(t)dt) \tag{5} \label{eq5}
$$
在物理意义上，它表示粒子群某一点的透明度，数值越大，说明粒子群越透明，光线衰减的幅度就越小。而透明度本身是关于 $\tau_a(t)$ 的方程，$\tau_a(t)$ 越大，$T(s)$ 就越小。由公式 $\eqref{eq2}$ 可以知道，，它是由粒子密度和投影面积决定的。这在直觉上也很好理解，如果粒子密度大，粒子本身也比较大，那么遮住光线的概率也会相应提升，自然透明度也就下降了。
$\tau_a(t)$ 也被称为光学厚度（optical depth），可以从这个图理解面密度的表示含义：
放射 (emission)
How to model the emission?
理解下面公式的关键在于对粒子跃迁释放光子的过程进行建模，在足够小的位置面上，粒子所释放的光子具有相同的辐射强度，粒子的密度以及粒子的大小越大，对应粒子释放的光子数也就越多。所建模的过程就是对某个位置面前后光线辐射强度的增量变化进行量化，这样就可以通过积分计算出任意位置面的放射强度。
对于某个位置面 $s$ 的粒子，对应的当个光子的放射强度表示 $I_e(s)$（同一位置面释放的光子辐射强度相同，不同位置面释放的光子辐射强度不同），计算位置面 $s$ 对应粒子的密度以及粒子大小（因为假设粒子在该平面上没有重叠，所以可以将衡量粒子密度以及粒子大小的程度当作粒子的覆盖比例）。类似与遮挡比例，可以计算出对应于某个位置面 $s$ 的覆盖比例，即 $\frac{\rho A E \Delta s}{E} = \rho A\Delta s$，因此，在位置面 $s$ 能够接收到粒子放射的平均辐射强度为 $I_e(s)\rho A\Delta s$，形式化地表示为：
$$
\frac{dI}{ds} = \frac{\Delta I}{\Delta s} = \rho(s) A I_e(s) = \tau_a(s)I_e(s) \tag{6} \label{eq6}
$$
其中， $I_e(s)$ 是一个关于 $s$ 的函数，表示圆柱体不同位置的粒子，对应所释放光子的辐射强度不同，而且微分量表示的是该位置面前后因粒子放射产生的平均辐射强度。 对于 $\tau_a(s)$ 而言，主要的考虑是粒子放射的光强与吸收类似，均是表示粒子的面密度（光学厚度）。
外散射 (out-scattering)
粒子除了吸收光子，也可能会弹射光子，这个过程称为外散射，即光子被弹射出原本的光路，导致光线强度减弱。同吸收一样，外散射对光线的「削弱」程度，也跟光学厚度相关，不过过程相对吸收来说又复杂一些，因此我们用 $\tau_s$ 来表示外散射对光线的削弱比例，以区别于 $\tau_a$，$I(s)$ 表示入射光的强度。同样地，这一过程可以表示为：
$$
\frac{dI}{ds} = \frac{\Delta I}{\Delta s} = -\tau_s(s)I(s) \tag{7} \label{eq7}
$$
内散射 (in-scattering)
光子可以被弹射走，自然就有其他光路的光子被弹射到当前光路，这一过程就是内散射。内散射的过程比外散射又更加复杂，因为弹射到当前光路的光子可能来自多条不同的光路，因此需要综合考虑其他光路的辐射强度以及各种弹射角度。
为了简化计算，我们认为其他光路的辐射强度为 $I_s$，而弹射到当前光路的能量损失比为 $\tau_s$ (注意这里和外散射使用的是相同的系数，既然都是散射，那有些性质上自然是共通的)。内散射的过程可以表示为：
$$
\frac{dI}{ds} = \frac{\Delta I}{\Delta s} = \tau_s(s)I_s \tag{8} \label{eq8}
$$
体渲染方程
经过以上四个过程的结合，我们可以得到对光线经过介质传播后到屏幕上建立的模型，形式化表述为：
$$
\begin{align}
\frac{dI}{ds} &= -\tau_a(s)I(s) -\tau_s(s)I(s) + \tau_a(s)I_e(s) + \tau_s(s)I_s(s) \newline
&= -\tau_t(s)I(s) + \tau_a(s)I_e(s) + \tau_s(s)I_s(s)
\end{align}
\tag{9} \label{eq9}
$$
其中，为了简化计算，定义 $\tau_t(s) = \tau_a(s) + \tau_s(s)$。吸收和外散射都会削弱光线的辐射强度，并且由于它们都和入射光有关，因此它们共同构成了体渲染中的衰减项 (attenuation item)，而放射和内散射都来自独立的光源，会增强当前光路上的辐射强度，因此被称为源项 (source item)。
根据以下一阶线性非齐次微分方程求解回顾，求解微分方程 $\eqref{eq9}$：
- 首先将其变换为如下形式，其中 $I’ = \frac{dI}{ds}$，$P(s) = \tau_t(s)$，$Q(s) = \tau_a(s)I_e(s) + \tau_s(s)I_s(s)$：
$$
\begin{align}
\frac{dI}{ds} + \tau_t(s)I(s) &= \tau_a(s)I_e(s) + \tau_s(s)I_s(s) \newline
I’ + P(s)I &= Q(s)
\end{align}
\tag{10} \label{eq10}
$$ - 根据一阶线性非齐次微分方程的通解公式 $\eqref{eq22}$，得到如下的渲染方程：
$$
\begin{align}
I &= e^{-\int P(u)du}(\int Q(t)e^{\int P(u)du}dt +C)\newline
& = e^{-\int_0^s P(u)du}(\int_0^s Q(t)e^{\int_0^t P(u)du}dt +C) \newline
& = e^{-\int_0^s \tau_t(u)du}(\int_0^s [\tau_a(t)I_e(t) + \tau_s(t)I_s(t)]e^{\int_0^t \tau_t(u)du}dt +C)\newline
& = \int_0^s [\tau_a(t)I_e(t) + \tau_s(t)I_s(t)]e^{-\int_t^s \tau_t(u)du}dt +Ce^{-\int_0^s \tau_t(u)du}
\end{align}
\tag{11} \label{eq11}
$$ - 当 $t = 0$ 时，记此时的光强为 $I_0$ 求解出系数 $C = I_0$，故得到最终的体渲染方程:
$$
I(s) = \int_0^s [\tau_a(t)I_e(t) + \tau_s(t)I_s(t)]e^{-\int_t^s \tau_t(u)du}dt +I_0e^{-\int_0^s \tau_t(u)du}
\tag{12} \label{eq12}
$$
其中，$I_0$ 表示入射光强度，在 NeRF 中被当作是背景光。公式 $\eqref{eq12}$ 表明，出射光的强度主要由入射光 (背景光) 和源项 (放射以及内散射) 等构成的，而这两者在传递过程中都伴随着相同的衰减系数 ($\tau_t$)。
回顾一阶线性非齐次微分方程求解
给定一阶线性非齐次微分方程：$y’ +P(x)y = Q(x)$，求 $y(x)$ 的通解以及 $x=s$ 点的解。
- 本质原理推导（最广泛的求解思路——积分因子）
- 先求解不定积分的通解，引入积分因子 $\mu(x) = e^{\int P(x)dx}$，并将原方程两边同时乘以积分因子 $\mu(x)$，得到：
$$
\mu(x)y’ + \mu(x)P(x)y = \mu(x)Q(x)\tag{13} \label{eq13}
$$
由于 $\mu(x)’ = \mu(x)P(x)$，因此上式可以进一步变形为：
$$
\frac{d}{dx}(\mu(x)y) = \mu(x)Q(x)\tag{14} \label{eq14}
$$
对上式两边关于 $x$ 积分得到：
$$
\mu(x)y = \int \mu(x)Q(x)dx + C\tag{15} \label{eq15}
$$
最后两边同时除以 $\mu(x)$ 即可得到关于 $y$ 的通解：
$$
\begin{align}
y(x) &= \frac{1}{\mu(x)}\int \mu(x)Q(x)dx + \frac{C}{\mu(x)}\newline
&= e^{\int -P(x)dx}\int Q(x)e^{\int P(x)dx}dx + Ce^{\int -P(x)dx}
\end{align}\tag{16} \label{eq16}
$$
上式即是我们最熟悉的一阶线性非齐次微分方程的通解公式，在高等数学课程中始终牢记的模板公式。但是如果我们想进一步求 $y(s)$ 公式中的积分上下限又该如何处理呢？
下面先根据高等数学中常见的做法，给一个直观的求解，直接将积分下限换成起点 $t$，积分上限 $s$：
$$
\begin{align}
y(s) &= e^{\int_t^s -P(x)dx}\int_t^s Q(x)e^{\int_t^x P(u)du}dx + Ce^{\int_t^s -P(x)dx}
\end{align}\tag{17} \label{eq17}
$$
其中，由于 $\int_t^s Q(x)e^{\int_t^x P(u)du}dx$ 的 $e^{\int_t^x P(u)du}$ 是关于 $x$ 的积分上限函数，因此其上限一定是 $x$。整理上式得到：
$$
\begin{align}
y(s) &= \int_t^s Q(x)e^{-\int_x^s P(u)du}dx + Ce^{\int_t^s -P(x)dx}
\end{align}\tag{18} \label{eq18}
$$2. 深入理解怎样确定积分上下限
重新审视积分因子，积分因子 $\mu(x)=e^{\int P(x)dx}$ 需要满足哪些条件才能作为辅助函数，使得其可以结合 $y$ 和 $y’$。1.首先一定需要以指数函数构造；2.导数中的另一部分一定是 $P(x)$。
那么，在添加积分上下限后，$\mu(x) = e^{\int_a^x P(u)du}$，即是我们学过的积分上限函数。
回到公式 $\eqref{eq14}$，我们需要对公式两边同时从 $t$ 积分到 $s$，即得到：
$$
\begin{align}
\int_t^s\frac{d}{dx}(\mu(x)y) &= \int_t^s \mu(x)Q(x)dx + C\newline
[\mu(x)y]_t^s &= \int_t^s \mu(x)Q(x)dx + C\newline
\mu(s)y(s) - \mu(t)y(t) &= \int_t^s \mu(x)Q(x)dx + C
\end{align}
\tag{19} \label{eq19}
$$
$x=t$ 表示计算的起点，一般都当作常数处理，因此可以进一步化简为：
$$
\begin{align}
\mu(s)y(s) &= \int_t^s \mu(x)Q(x)dx + C + \mu(t)y(t) \newline
y(s) &= \frac{1}{\mu(s)}\int_t^s \mu(x)Q(x)dx + \frac{C}{\mu(s)} + \frac{\mu(t)y(t)}{\mu(s)}
\end{align}
\tag{20} \label{eq20}
$$
对于上式，非常有意思，我们实际上带入积分上下限想计算 $x=s$ 点的解，但是 $x=s$ 点可以表示任意的点，也即当作变量处理，这样就可以得到始于 $x=t$ 点的任一点的解。重新整理一下，得到：
$$
\begin{align}
y(s) &= e^{\int_a^s-P(u)du}\int_t^s Q(x)e^{\int_a^xP(u)du}dx + e^{\int_a^s-P(u)du}(C+ \mu(t)y(t))\newline
y(s) &= e^{\int_a^s-P(u)du}\int_t^s Q(x)e^{\int_a^xP(u)du}dx + e^{\int_a^s-P(u)du}C’\newline
y(s) &= \int_t^s Q(x)e^{-\int_x^sP(u)du}dx + e^{\int_t^s-P(u)du}e^{\int_a^t-P(u)du}C’\newline
y(s) &= \int_t^s Q(x)e^{-\int_x^sP(u)du}dx + e^{\int_t^s-P(u)du}C’’
\end{align}
\tag{21} \label{eq21}
$$
上式中，对于 $C’’ = e^{\int_a^t-P(u)du}C’$ 的处理，因为我们并不关心 $\int_a^t$ 到底是怎样计算的，而是只需要关注于我们的计算区域 $\int_t^s$ 中的值计算，其他参数都可以通过初值来求得。
通过上式就可以计算出关于任一点 $s$ 的解，初始点 $x = t$，如果 $t=0$ 进一步可以表示为：
$$
\begin{align}
y(s) &= \int_0^s Q(x)e^{-\int_x^sP(u)du}dx + e^{\int_0^s-P(u)du}C’’
\end{align}
\tag{22} \label{eq22}
$$
- 高等数学中的常数变异法
齐次方程的通解为：$y = Ce^{-\int P(x)dx}$，利用常数变易法对齐次方程的解求导：
$$
y’ = C’e^{-\int P(x)dx} - CP(x)e^{-\int P(x)dx}
$$
这里 $C’$ 是关于 $x$ 的函数，因此保留导数表示，将其带入原方程可得
$$
\begin{aligned}
& C’e^{-\int P(x)dx} - CP(x)e^{-\int P(x)dx} + P(x) Ce^{-\int P(x)dx} = Q(x) \newline
& \Longrightarrow C’e^{-\int P(x)dx} = Q(x)
\end{aligned}
$$
求解出系数 $C$ 得到：
$$
C = \int Q(x)e^{\int P(x)dx}dx + C_1
$$
故此，非齐次微分方程的通解为：
$$
\begin{align}
y &= e^{-\int P(x)dx}(\int Q(x)e^{\int P(x)dx}dx + C_1)\newline
&= e^{-\int P(x)dx}\int Q(x)e^{\int P(x)dx}dx + e^{-\int P(x)dx}C_1
\end{align}
$$
其中，$e^{-\int P(x)dx}C_1$ 表示齐次方程的通解，$e^{-\int P(x)dx}\int Q(x)e^{\int P(x)dx}dx$ 表示非齐次方程的特解，两者相加即为非齐次微分方程的通解。
上式是高等数学中采用的常数变异法的求解过程，关于如何将其引入积分上下限可以参见公式 $\eqref{eq21}$。
参考链接
Scratchapixel. From the Radiative Transfer Equation to the Volume Rendering Equation.
Jermmy, AI小男孩. NeRF入门之体渲染 (Volume Rendering).
体渲染(Volume Rendering)公式推导