总想写一下体渲染技术的历史根源，虽然现在各种渲染方程层出不穷，但追根溯源还是最开始的光传输理论等方法。做了这么久的体渲染了，也没有好好把最初的算法进行一下回顾，所以这里我要开始写一篇很长很长的博客，把体渲染最古老的论文都研究研究，写一个尽量通俗易懂的包含介绍性和知识性的文章。
Table of Contents
怎么来模拟烟雾？
Blinn 在1982年图形学还刚起步不久的时候，发表了一篇论文，在论文中阐述了自己的模型，为了计算光在烟雾中的反射。
大家一看就能知道是最简单的反射问题。
光与烟雾作用的模拟主要分为两部分：
- 设计一个函数，来模拟不同位置烟雾的浓度（烟雾粒子的密度）
- 根据这个函数来模拟光的交互过程
但是因为这个模型并不是体渲染当前的主要思想，所以我暂且不再进行表述，如果以后有需要我再来补充这一节。
渲染方程
Kajiya发布了一篇标题只有三个单词的文章：渲染方程。该文章不但综合了以前的多种模型，还发布了一种比较新的技术。不过其实这篇文章主要还是讲述的光在不同表面如何进行散射的问题。注意这种散射与粒子中的体渲染散射还是不一样的，但是基本思路差不太多。
我们先看一下这个渲染方程：
如果之前看过体渲染方程，对这个肯定非常熟悉。但是这个方程并不是体渲染的方程，而是比较简单的表面渲染。几何项可以理解为一个缩放关系，表示从x'到x点的图中被其他微粒表面遮挡多少。整个方程可以理解为：
从x'到x的光 = 从x'自发出的光 + 全部从四面八方射到x'然后再散射到x方向的光
从x''射到x'上再反射到x方向的光有多少，这个比例就是
。
渲染方程深入
我们定义能量：
表示的是每个时间微元，从 x 微面到 x' 微面传递的能量。
只看 x' 点发射的能量：
这里的1/r^2意思是能量随着距离 r 变长，呈r^2倍衰减。
我们只看从 x'' 处发射到 x' 然后散射到 x点的能量：
辐射度强度定义为单位时间单位立体角中投射到单位面积上的能量：
立体角可以这么进行计算： 表示的是光投影到的面的面积。
然后对前面的公式进行替换，最后可以得到：
即 
。
还有自发射能量 ： 
渲染方程在多个表面的散射
就是光从 x'' 这里发出来，经过表面 x' 反射到 x 上的能量。这里面主要是由 
来控制传输到 x 方向上的能量比例的。这个函数叫做BRDF双向反射分布函数（在论文中称作BRF，双向反射函数，估计是那个年代BRDF理论还不是很完备）
优化（简化）方案
我们把一开始的方程进行变形：
变形为：
M是一个线性操作符，对应着积分。
然后再进行整理：
无穷级数收敛的一个条件是算子M的谱半径小于1。对该展开式的物理解释很有吸引力。给出了 x 点与 x' 点间辐射传输的最终密度，即直接项（
）、一次散射项（
）、二次散射项（
）等的总和。
我们叫 Utah 近似（来自Utah 大学）：
注意这个近似中，M只作用域后面的
。
还有一种近似方案，是基于辐射度的近似——近似BRDF积分，但是因为也与后面的其他研究意义不是很大，所以就不再介绍了。
虽然近似的方法有很多，包括使用后面的MC蒙特卡洛来计算积分，以及分层MC计算积分，但是万变不离其宗，都是研究的如何近似积分的，因为与渲染的关系并不是很大，这里就不再详细展开。
结论
这篇文章介绍的技术，主要是光从一个表面散射到另外一个表面的过程的方程及其计算方案（简化和优化计算），其实，所有的光照模型无非都是这种思想：
- 光跟物体接触以后是穿透还是反射（散射）还是两者都有？
- 光散射的能量是怎么分布的
理解了这个目标，就能很清楚的知道每一个论文具体研究的是什么内容了。
3D标量场可视化
我们迎来了重头戏——大名鼎鼎的 density emitted model（密度发射模型）。
Sabella在图形学顶级期刊SIGGRAPH上发表的关于体渲染的论文，提出了光线追踪渲染3D标量场的技术，该算法计算沿采样光线的场的特性，如衰减强度（attenuated intensity）、峰值密度（the peak density）和重心（center of gravity）。
渲染体数据的奥义是能够可视化里面任何的数据区，比如比较高的标量值，比如“热点区域”，还比如场梯度变化（即类似标量场的数据形成了怎样的表面）。这个论文里把体空间（标量场）设为一个可变的密度发射对象。
该模型的三个特点：
- 颜色映射用于可视化各个属性。
- 使用相位函数来可视化场梯度。
- 重新设计了Kajiya和Von Herzen光线跟踪算法，以提高计算效率，消除阴影，同时保持遮挡效果。
两种已经存在的渲染模型（老方案）
这两种模型都是表示的光交互模型。
Cross-section rendering
简单描述一下这个方法：这个方法就是“基于纹理映射的体渲染算法”的最初版本，把三维图像比作一堆二维图像的叠加，然后控制只显示哪些部分，来得到最终的图像。
图片来自于《Real Time Volume Graphics》
例如，MIP最大密度投影，选择每个成像采样光线上最大的一个值作为最终结果。
Threshold rendering
主要用于医学图像显示方面，只对里面一定阈值内的数据感兴趣，比如人体CT图像中的骨组织，可以使用该方法。有三种阈值渲染的方法：
- 从后到前遍历
- 光线追踪
- 表面重建（建立三维面模型）
阈值渲染模型会丢弃三维图像中保存的大部分数据。除了像医学成像这样的应用程序（其中关注的是少量的密度范围），它就像横截面渲染一样，因为它需要许多视野，这些视图必须在我们心里进行组合，以便用户完全感知该标量场。在像地球物理成像这样的应用中，三维地震图像并不像医学图像那样描述连贯性。因此，阈值渲染并不令人满意，因为它会生成碎片而不是粘合曲面。
散射模型
i） 当光线被靠近观察者的云层散射时，模型的一部分会被遮挡。
ii）根据光源的位置创建阴影。
iii）颜色变化是由于不同波长的散射量不同而产生的分离
密度定义
这里的密度定义是无维度关系的，而且因为是随时变化的，所以要用微分形式：
是每个粒子的体积，V是包含粒子的当前微元的体积。
表示的是当前空间微元的粒子数。
亮度方程
我们令密度场表示为
在采样光线上的单个微元中的粒子发光强度是：
（单个粒子微元发光的亮度）
（采样光线沿途中所有的粒子数）
表示的是到人的眼睛的光强相当于粒子发出的光，同时还在概率P(0;V)下衰减。我们假设概率符合泊松分布：
所以单个微元发光亮度到最终图像上为：
总亮度就是：
论文里的这个公式看起来不是很丑（主要是论文太老，翻新以后显示的其他公式实在是太丑了），所以就直接截图了。
我们用常量来代替
，然后令
来标准化整个方程（毕竟我们是在计算机程序中应用，所以不需要完全按照物理的辐射度来优化，只需要找到一个对应关系就好了），得到：
看着是不是超级像体渲染的方程？要注意一点，这里的衰减虽然是按照的概率分布计算的，但是因为概率分布与粒子浓度有关，所以最后得到的方程也是与“发射吸收方程”一样。
因为 的值在 0-1之间，所以我们采用转换：
我们总体介绍一下这个方程：
高的 表示衰减更快，高的
会增强密集部分相对于较扩散区域的外观，低的
表示更弥散：
离散化的方法就不再介绍，与渲染的研究没有什么关系。
颜色映射
现在要注意到，我们以前解决的问题只是明暗度，但是如果我们想显示彩色图呢？（例如下图，左图是明暗度图，右图是上色了的图）
但是因为这种颜色映射方式比较古老，与现在的“传输函数”方面的研究没有太大关系，所以我们就暂且不提了。
结论
与后面的渲染方法大同小异，考虑到了粒子的衰减和粒子的发光，虽然本文的研究中没有阴影效果（没有体现出散射），但是仅靠发光和衰减渲染出的地震波图像效果还是可以的。
光传输理论在3D标量场可视化中的应用
一个里程碑式的文章，诞生于1990年，由Krüger发表于IEEE Vis会议。
以前很多研究主要是将3D标量场里面的所有体素都当成粒子，然后研究光传输是如何与这些粒子交互的。这篇论文乍看一眼其实和最早的研究很相近，但是实际上它们之间还是有区别的，它是假设有一些虚拟的粒子，这些粒子穿过3D标量场，与场内发生作用，最后投到成像面上的过程。
渲染方程
s·D可以理解为 求微分。 表示的是吸收，
表示的是总的散射，
表示衰减系数。
即可以理解为：（注意减号表示的是衰减的正值，相当于“加上负值”）
光强度在x点的变化率 = 光在x点的衰减值（包括被吸收和被散射） + x点自己发出的光 + 抑制能量变化过大的项（Sin乘以能量变化率，例如如果能量上升，就减去一些，上升越大减的值就越大） + 从四面八方散射到x点上的光的总和
做一下积分，就能够得到：
我们可以把这个公式理解为：
在x点（已经传输了距离R）的总能量 = 原始能量Is到x点以后的衰减 + 沿途中每个点（连续点）的能量之和（包含被吸收和被周围散射过来的光增强以及自发射能量这三部分）
数据特征映射
源项
这里的模型叫做“源——衰减 模型”，其实也就是“吸收发射模型”。
我们定义3D标量场为 F(x) ：
但是这种方程适用性不好，不能显示深度信息和场覆盖作用（遮挡），所以使用下面的方程：
s表示采样光线（积累）的方向， 表示的是当前位置的梯度（法向量）。
表示的是场梯度的绝对值。
这种方程很适用于医学影像体渲染，因为可以很好的区分边界。
衰减项
还是分了两种情况，根据场标量值决定的衰减值以及根据场梯度决定的衰减值。
散射项
是散射系数：
以及散射相位函数：
注意该相位函数描述了从前面看当前粒子的透明度 以及 从后面看当前粒子的透明度
。
颜色映射
通过简单映射的方法：
可以实现的效果是：
加入与物理有关的“停止能量”项：
可以得到：
增强重要区域
我们更换新的数据场：
在原始数据场F1上加一个偏离项。偏离场 
可以用方差，自相关来进行描述：
括号里面代表空间平均值，Cr表示自相关长度，与空间网格大小有关。粒子密度 I 也需要矫正：
H 和 N取决于映射。密度方差可以表示为：
增强显示出一些点：
但是因为描述原因，并不是很清楚里面的更深层次的设计细节，所以这里也只是简单提一下。
结论
给个结论吧。上个世纪的老论文了，很多描述都不清不楚的，但是总体思想其实就是把成像当成粒子穿过体空间，然后粒子与体空间的粒子进行交互，然后最终结果作为成像结果的过程。
体渲染数学方法和模型
Hege等人在1993年发布的文章中具体阐述了体渲染的数学方法，凡是做过体渲染的人不可能不知道这套体系。其实说实话，从内容来看，接下来的这篇文章更像是前面的技术的一个总结。不过这篇文章确实阐述的是比较完备的体渲染算法，它还是基于线性传输理论。
理论介绍
光具有波动性，进而有多种状态，比如干涉衍射等，还有不同的波长，这里我们不考虑，而是通过应用几何光学来研究渲染方程（即线性传输理论）。在几何光学中，使用频率而不是波长，因为当折射率改变时，频率会保持不变。而且有个小细节：光传播需要时间，模拟光的弛豫是可以的，但是没有必要。我们假设一开始就是光的稳定态。
虽然前面已经介绍了这个方程，但是为了全面这里再粘贴一下：
如果用单位体积单位时间dv在x点上的光子密度
表示，光子的速度是 c，一个光子的能量为 （普朗克常量，物理学知识）：
辐射能：
辐射能密度（点光源）：
（单位立体角的辐射能）。
射出的能量：
射入的能量：
吸收发射和散射
吸收项：
表示在当前微元吸收的能量，
表示的是衰减（吸收）系数。da表示面积微元。
注意到吸收项与角 没有什么关系。
被称作频率为 的光子在材料中的“平均自由路径（mean free path）”，它定义了特征长度规模（characteristic length scale）。
发射项：
区分真实的或热的吸收和发射过程以及散射过程是很重要的。在前一种情况下，从光束中移除的能量转化为材料热能，能量分别以材料能量为代价发射到光束中。相反，在散射过程中，光子与散射中心相互作用，并从沿不同方向移动的事件中出现，通常也具有不同的频率。如果频率不变，我们就说弹性散射，否则就说非弹性散射。
散射项：
用 
表示吸收系数，用 
表示散射系数，这样就能得到总的吸收（衰减）系数：
表示为 albedo，当albedo为1的时候表示，射过来的光会被全部反射掉，
同理，发射系数也应该包含源项发射 和散射到表面的光 ：
我们定义散射外出的能量方程：
该方程表示从 n 方向到 n' 方向，以及从频率 到频率
' 的能量。
其中 p 表示相位函数：（这里的4π别的文章中说过，并没有什么含义，只是为了抵消掉球面积分因子，让最后积分结果为1）
以及散射到x点的散射能量方程：
弹性散射中，
。注意其实大部分应用里都是弹性散射，所以不需要考虑频率的变化。同时，很多相位函数只关心 cosθ ，而不是n和n'，这叫做各向同性材料。最简单的各向同性材料是：
表示的是光向各个方向散射量都是一样的。
- Rayleigh散射
满足球面积分为1，该函数对于前向散射和后向散射是没有什么区分的（仅仅依赖于 θ）。
- Henvey - Greenstein相位函数
改变参数k提供了前向散射（k>0）、各向同性散射（k=0）和后向散射（k<0）之间的连续介质。
传输方程的微分形式
x+dx的能量与x处的能量之间的差值 = 发射出的能量与吸收的能量之间的差值，所以：
同时，
这就是《Real-Time Volume Graphics》第一章给出的方程。
同时注意散射包含一个积分项：
因为一般是弹性散射，我们不考虑频率的变化：
边界条件
所谓边界条件，就是微积分方程中的一些限定，比如遇到未定义区间应该怎么处理。
对于不透明表面，边界条件很容易指定。我们假设表面法线N总是指向辐射场所在的体。
explicit boundary conditions：
辐射到体中的强度值由表面发射函数E给出。
表示指向体的所有方向的集合。
implicit or reflecting boundary conditions：
surface scattering kernel ：k ： 
explicit 以及 reflective boundary conditions 也可以结合起来。
考虑到表面的时候，边界条件就需要区分正表面
以及负表面
：
在这个定义中，我们选择了所有输出的辐射。类似地，我们也可以选择入射到表面元素上的所有辐射：
我们把这些都结合起来：
即 x点处的能量 = x点发出的能量 + 其他方向射到x方向的能量
在没有参与介质的情况下，该方程简化为众所周知的渲染方程，该方程是所有表面绘制的基础。在这里，在更一般的情况下，它只是作为传递方程的边界条件。
渲染方程的积分形式
表示 沿着 
上的方向积分。
在
和 
之间的光学深度
回忆一下，
表示的是一个频率为的光子的平均自由路径。
应用到方程中：
光学深度可以分解为：
所以得到体渲染方程：
该方程就是大家特别熟悉的体渲染（光学）方程。
渲染方程
我们再进行一些简化，就可以得到著名的渲染方程：
以及如下定义：
即：
同样，后面的解决思路不再描述，因为解决方案无非就是离散化之类的，我们会在最后一篇survey里面详细介绍一下。
DVR光学模型
终于到封顶之作了。虽然其实在Hege的模型已经把很多方面都阐述的很清楚了，但是我们仍然需要这最后一篇论文来继续强化体渲染的模型。
1995年Max在图形学顶会之一的TVCG上发表了一篇文章，详细的介绍了体渲染的算法细节的方方面面。这是一篇survey，这里面的模型包括：
- 只有体素吸收光的体渲染方程
- 只有体素发射光的体渲染方程
- 体素既能发射光又能吸收光的体渲染方程
- 外部照明单散射没有阴影
- 外部照明单散射有阴影
- 多重散射
这篇survey的讲述布局也是按照的这个顺序。我们从头到尾把这个survey讲述下。
只有吸收的体渲染方程
最简单的参与介质有完全黑色的冷粒子，它们吸收所有被拦截的光，而不散射或发射任何光。
设粒子是个球体，半径为 r ，投射面积 ，令
作为单位体的粒子数（密度），设有一个圆柱形面板：
厚度为 ，面积为 E ，所以体积就是
，总共包含
个粒子。
如果圆柱横截面积足够小，则粒子投到这里的覆盖率就会很低，所以 ，光能通过这个圆柱体的分数表示为：
趋近于0的时候，重叠率也会趋近于0，得到衰减方程：
I(s)表示在s点的光密度。
叫做衰减系数。
其中 表示的是在s=0位置的光强度。
表示的是从0到s位置的穿透率（透明度）。
在体渲染中， 是衰减系数，定义从 0 到 l 位置的不透明度（opacity）：
后面介绍了一坨关于传输函数的问题，很简单，这里懒得再讲了。
只有发射的体渲染方程
除了介质会吸收光，介质也会发射光或者接受外部照明，这里我们只看介质的发光：
当然真实粒子既吸收也发射光，但在粒子大小或数密度接近零的极限下，当发射以补偿的方式趋于无穷大时，我们可以忽略吸收。这是非常热的稀薄气体的情况，它发光但几乎是透明的。我们将通过假设上面讨论的小球形粒子是透明的来模拟这种情况，然后在下一节中我们将包括它们的吸收。
如果粒子是透明的，密度为C（每个单位目标面积）， 的发光流为
，每单位面积的发光流为
包含了发射和吸收的体渲染方程
实际云中的粒子阻挡入射光，并自己也发光。因此，实际的微分方程应包括源项和衰减项：
（后面我们要介绍一个特别的
，但是这里我们假设这个源项是自己通过传输函数设置的）
我们通过一系列的解方程手段得到下面的方程：
然后设：
表示从位置 s 处一直到D处的衰减：
这个公式应该这么理解：最后得到的光强度 = 初始光强从初始位置一直衰减到D的强度 + 每个位置的光源发出的光强从当前位置一直衰减到D的强度的累加。
计算方法
对于某些传递函数和插值方法，可以解析地计算上式中的积分。但是，对于更一般的情况，需要数值积分。积分的最简单数值逼近是黎曼和：
变为：
后面的论文因为都是关于黎曼和的，就不再解释了，直接贴图：
粒子模型
（相当于在颜色C上预乘了一个透明度）
C设为常数：
粒子模型对应于发光粒子的物理情况，但有时定义g(s)与(s) 无关是很方便的。例如，可以直接用标量场 f 给出 g，甚至可以用与确定
无关的标量场表示g。这比粒子模型的灵活性稍大，因为它允许g即使在
为零时也不为零，从而允许完全透明的发光气体，不需要一个无限明亮的C(s)。
有散射无阴影的模型
还记得之前说过的Utah吗，就是指外部光源没有阻挡的打到要渲染的物体上（所以没有阴影）。
是从 射到X上的光源，r表示BRDF。
表示粒子的反照率，
代表散射而不是吸收，p是相位函数。
对于球形粒子或者随机反射粒子来说，只关心 与
之间的夹角。对于漫反射散射来说，比如粒子比波长要大很多，可以使用：
根据表面法向的渲染：通过求梯度得到表面法向
，
注意我们如何应用散射项。方法是在源项中包含两个部分，自发射项和散射项：
有散射有阴影的模型
方向的光源射到X点上的强度。实际中，积分并不是到无穷，到体的边界就好了。
令
，重写
得到：
用散射和遮挡光来替换，得到：
表示阴影光线，与递归光线跟踪中使用的“阴影触角”相对应，但阴影触角沿着主光线发送到每个点X-tw处的光源并返回分数透明度。当 为常数且存在不透明多边形对象时，请考虑多个光源。
多重散射
上面的方法只有在反照率或密度较低时才有效，因此不太可能发生多次散射。这在大气云中通常不是这样，所以在图6中，远离光源的云端看起来不自然地黑暗。为了纠正这一点并解释多重散射，人们可以应用最初在热辐射热传输领域中开发的“光能传递”方法。多次散射计算对于高反照率参与介质的真实感渲染非常重要，但是在计算机时间上很昂贵，并且对于大多数科学可视化应用来说，这是一种过度消耗。
多重散射涉及到方向相关的光流，因此有必要求出 I(X，w)，即每个光流方向上每个点X的强度。与光流相反，沿X的观察光线距离s处的点是 。将 4π 单位球体上所有可能入射方向
的光在
处的散射积分，得到源项：
计算方法不再赘述，就是如何近似这个积分方程的方法罢了。
参考文献
J. Blinn. Light Reflection Functions for Simulation of Clouds and Dusty Surfaces. In Proc. of ACM SIGGRAPH, pages21–29, 1982a.
J. Kajiya. The Rendering Equation. In Proc. of ACM SIGGRAPH, pages 143–150, 1986.
P. Sabella. A Rendering Algorithm for Visualizing 3D Scalar Fields. In Proc. of ACM SIGGRAPH, pages 51–58, 1988.
W. Krüger. The Applicaton of Transport Theory to Visualization of 3-D Scalar Data Fields. In Proc. of IEEE Visualization,pages 273–280, 1990a.
H. Hege, T. Höllerer, and D. Stalling. Volume Rendering: Mathematical Models and Algorithmic Aspects. TechnicalReport ZIB-93-07, ZIB, 1993.
N. Max. Optical Models for Direct Volume Rendering. IEEETransactions onVisualization and Computer Graphics, 1(2):99–108,1995.
NAKAO,Megumi. Real-time Volume Graphics[J]. Medical Imaging Technology, 2007, 25.