随着计算机技术的兴起与不断发展，计算机的性能在其中得到不断提升，利用计算机来渲染模拟真实生活中的场景受到越来越多人的关注。近年来计算机图形学有诸多的研究成果，这些理论成果及其实际应用使得计算机通过图像表达现实生活影像的能力越来越强，图像达到了很高的拟真度。在计算机图形学中，高品质渲染作为近些年的关注焦点之一，其主要是以基于光线追踪思想衍生出的各种算法予以实现的。而基于光线追踪思想的路径追踪算法则是本文研究的重点。

计算机图形学中关于真实图形渲染的实现手法基本上是以通过离线或实时计算机渲染呈现现实生活场景来实现的，呈现品质的好坏直接影响着人眼观察时的感受，因此在未来的高拟真虚拟现实的应用领域有很重要的地位。而图像品质的提升伴随着实时交互性能和静帧渲染速度的下降，使得所有图形学算法必须在上述二者之间找到平衡。同时，由于近十年来计算机硬件性能的飞速发展，尤其是图形处理器(GPU)浮点运算性能与显存带宽的指数型增长趋势，使得基于GPU的并行运算的思想越来越受广大研究人员的欢迎，许多基于串行算法的应用程序在更改为基于GPU型的并行算法后其效率有几十倍的提升，因此在如此背景下，使得在GPU平台上实现路径追踪等大规模与大计算量的算法成为可能。多年前国外已有基于NVIDIA GPU的路径追踪渲染器成品面世^[1]，其渲染效率之高，甚至在某些低多边形数量的场景中能够达到基础级别的交互帧率，从而作为替代传统光栅图形渲染管线的另外一种渲染机制并在游戏中得以应用，且近些年不断发展，成为成熟的游戏实时图形渲染引擎之一。

近几十年来，有不少研究者已经对有关全局光照的诸多算法做出研究，并得到了不错的效果。在1986年，Kajiya在论文“The rendering equation”^[2]中实际上推导得出了全局光照方程应具有普遍形式

$ \begin{array}{l} {L_o}(x, {\mathit{\boldsymbol{\omega }}_o}) = {L_e}(x, {\mathit{\boldsymbol{\omega }}_o}) + \\ {\smallint _{\mathit{\Omega}} }{f_r}(x, {\mathit{\boldsymbol{\omega }}_i}, {\mathit{\boldsymbol{\omega }}_o}){L_i}(x, {\mathit{\boldsymbol{\omega }}_i})({\mathit{\boldsymbol{\omega }}_i}\cdot\;n){\rm{d}}{\mathit{\boldsymbol{\omega }}_i} \end{array} $

即

$ {L_o} = {L_e} + {L_r} $

式中，${L_e}$为自身光照强度贡献值，${L_r}$为环境光照强度贡献值。

上式说明了光子在漫布场景且其能量分布趋于平衡时所满足的关系。对于全局光照，理论上需要进行半球面积分来得到该点环境光强的贡献。考虑到积分在实际计算时的困难，近几十年出现了许多算法来近似求解${L_r}$项，求解${L_r}$项的不同方法实质上造就了不同的全局光照算法，路径追踪算法即是其中之一。

路径追踪算法是由Kajiya在1986年发表全局光照渲染方程之后连带提出的渲染概念。其主要思想为，从用户定义的虚拟摄影机作为采样光线的出射点，以像素为单位对屏幕空间内的所有像素进行采样，每条经摄影机“发出”的采样光线经过以下几个过程：

1) 光线与场景物体求交；

2) 获取交点处的物体纹理与材质信息；

3) 根据物体表面材质信息确定出射光方向；

4) 具有新方向的采样光线递归与场景求交，并重复步骤2)—4)，直到满足最大递归深度为止；

5) 收集采样光线在遍历场景过程中的信息，得到该像素在一次采样后的颜色；

6) 对同一个像素进行重复采样，得到的大量样本取平均值，得到逐渐收敛的像素颜色。整帧收敛的像素组合成一个完整的渲染结果图。

在上述过程中，步骤1)占据了路径追踪算法的绝大部分时间，由于采样光线的方向随机，初始点随机(分布在场景物体表面)，场景几何信息复杂，若要做到快速的路径追踪则一般需要借助空间剖分算法。相关研究人员对此做出了不少的研究，并得出了不少有效的剖分方式。如利用最小包围盒来以物体为单位划分空间^[3]；又如二叉空间剖分(BSP)树^[4]、基于层次包围盒(BVH)^[5]的系列算法与多级均匀网格^[6]技术等，极大地提高了三维网格的遍历效率。本程序中使用的是简单且快速的基于表面积启发式剖分(SAH)的空间剖分算法^[7]。

在步骤2)3)中，经典全局光照方程^[2]中对其描述的是积分形式，具体是通过双向反射分布函数(BRDF)来完成对新的出射光线方向的确定。利用BRDF进行对物体材质表面的渲染时均能够体现极佳的拟真度。但传统BRDF仍是以积分形式体现，且在进行采样时单一的材质表面需要大量数据集来拟合，其典型的数据集合体现为枚举参数变化而成的大型数据表格。参数的细致化会让最终的渲染图像有极佳的拟真度，但这些因素使得传统的基于积分与数据采样方式得到的渲染图像的生成效率较低，且当所模拟的材质表面参数发生更改时，需要重新测量目标表面BRDF参数，且测量过程复杂，测量种类有限，此类情况在一些要求实时更改材质表面参数并进行快速预览的条件下就显得较为局限。因此在实际应用中，相关领域研究人员基于一系列测得的实验数据的条件下，分别对BRDF采样方法进行函数拟合近似，得到了不错的效果；文献[8]研究了BRDF的重要性与应用范围，重点考虑了其在基于GPU的实时渲染软件中的应用，配以金属车漆进行实例分析，将BRDF的拟合方式划归为两类，一类是经验拟合型，一种是实验拟合型，总结了上述两种拟合的优缺点，最终得出较为理想的BRDF所应具有的形式特点：真实、简单与低内存开销；文献[9]通过极少的真实各向同性BRDF采样值，重建出BRDF在全半球面的分布情况，使得采样速度获得巨大提升，并与参考图像对比后相差无几；文献[10]通过对常见材质进行分类，将150种分类后的材质视为模板，对每种材质模板，生成分布特征不同的采样点，利用极少的采样值最大限度还原原本材质的BRDF特性，使得最终渲染出的物体表面具有较高拟真度的同时，也能够通过不同模板的参数合成来近似模拟未知表面的BRDF特性，具有较高的灵活性；文献[11]根据两组与物体成角(一组以较小入射角拍摄，一组以较大入射角拍摄)的近距拍摄图，获得各向同性的平面材质样本的BRDF分布特性；综合上述文献，发现现今对BRDF采样的模拟方式仍存在复杂或较高内存开销的缺点，进一步地，出于对通用性与实用性的考虑，概括出经函数拟合的BRDF应具有以下几个特征：

1) 参数少且易于理解和控制，参数具有一定的物理意义；

2) 能够近似模拟生活中的大部分材质表面；

3) 采样过程能够足够快，使得该采样算法不会成为全部路径追踪算法中的性能瓶颈；

4) 极低硬件存储开销(主要为内存存储开销)。

基于上述要求，国外诸多研究人员分别通过不同的拟合方式提出了可用的实际算法或概率分布模型。例如Cook-Torrance微表面算法模型^[12]，其假定光线射入粗糙物体时均发生理想反射，但粗糙物体表面的法线分布方向不均，因而造成非理想反射的效果，可较好地模拟粗糙微表面的光学特性；又如Schlick模型^[13]，用于模拟高光反射时的菲涅尔光学特性；以及He-Torrance-Sillion-Greenberg模型^[14]等。但该模型由于表达式过于复杂而使得实现难度较大，未被多数研究者采用。在国内，杨健等人^[15]提出了一种基于光线跟踪的表面微结构渲染方法，将传统的基于图像的表面细微结构渲染转为基于高度图的方法，使得原始的上百MB的数据量压缩为5 MB以内，但仍未给出较为简化的采样形式；李敏等人^[16-17]分别基于光线追踪技术与Blinn-Phong算法，利用多项式公式拟合的方法大幅改善了红外图像中对高光部位的仿真效果与渲染速度，但未给出可见光波段BRDF的简化采样方法。

上述诸多国外模型均从不同角度对物体表面的特性进行建模，得出了较为理想的渲染结果；且国内已有诸多研究者对BRDF进行了研究，但均未得到简化的全局光照方程的形式以及简易的路径追踪算法。如果能够通过经验公式直接生成出射光线的方向，则可直接利用相应的路径追踪算法执行物体的渲染过程而不必关心计算密集度和复杂度较高的采样环节，这对于整个路径追踪算法效率的提升是十分有益的。且如果经验公式中的参数少且易于控制，通用性强，就可以在仅改变极少参数的情况下快速预览物体表面特性的变化。考虑到理想漫反射表面反射环境光时所具有的特性^[18]与镜面反射时的特点，本文提出一种算法用于解决传统的基于BRDF采样的渲染方式中高硬件存储开销与速度较慢的缺点，亦或作为上述国外诸多已有研究成果的替代方案。具体思路为通过函数值的分布特性生成具有某种分布特性的权重值，并将两个参考向量的方向根据该权重值进行线性插值并规范化后，得到最终具有某种分布特性的出射光向量。当函数选取合适时，所模拟出的出射光的分布能够很好地模拟该表面实际对光的散射情况(图 1)。实际结果表明，当采用余弦与指数函数作为权重生成函数时，该算法能够真实模拟生活中的磨砂金属材质表面，且该算法能够通过控制参数，分别得到高光反射，理想漫反射与镜面反射的分布特性。

由于现在的基于真实度的图形渲染基本都需要求解Kajiya提出的全局光照方程，而传统的全局光照方程又引入了BRDF的概念，且由于传统的BRDF采样方法具有高内存开销与灵活性较差等缺点，因此用原始方法实现Kajiya所提出的全局光照算法在实际中会遇到极大的困难。

本文目的在于简化全局光照方程中有关BRDF采样的方法，提出一种新的采样方式来解决传统BRDF采样时的高硬件存储开销与形式复杂的缺点。为此，需要先对全局光照方程的形式进行简化，并根据简化后的方程形式设计新采样方式，随后给出与之配套的一整套路径追踪算法和实现的关键伪代码。最终利用MERL BRDF数据库(https://www.merl.com/brdf/)所测得的实际数据进行原始BRDF算法的渲染，再与该算法的渲染结果进行对比，得出该算法的优势与不足。

受Kajiya全局光照方程^[2]的启发，为简化光照方程对其做简单分析，首先观察光照方程形式为

$ \begin{array}{l} {L_o}(x, {\mathit{\boldsymbol{\omega }}_o}) = {L_e}(x, {\mathit{\boldsymbol{\omega }}_o}) + \\ {\int_{\mathit{\Omega}} }{f_r}(x, {\mathit{\boldsymbol{\omega }}_i}, {\mathit{\boldsymbol{\omega }}_o}){L_i}(x, {\mathit{\boldsymbol{\omega }}_i})({\mathit{\boldsymbol{\omega }}_i}\cdot\mathit{\boldsymbol{n}}){\rm{d}}{\mathit{\boldsymbol{\omega }}_i} \end{array} $

(1)

式中，${\mathit{\boldsymbol{\omega }}_i}$为入射光方向，${\mathit{\boldsymbol{\omega }}_o}$为出射光方向，$\mathit{\boldsymbol{n}}$为表面法向，$x$为入射点，其示意图如图 2所示。

令第2项积分式中

$ \begin{array}{l} \;\;\;\;\;\;\;\;{L_r}(x, {\mathit{\boldsymbol{\omega }}_o}) = \\ {\smallint _{\mathit{\Omega}} }{f_r}(x, {\mathit{\boldsymbol{\omega }}_i}, {\mathit{\boldsymbol{\omega }}_o}){L_i}(x, {\mathit{\boldsymbol{\omega }}_i})({\mathit{\boldsymbol{\omega }}_i}\cdot\;\mathit{\boldsymbol{n}}){\rm{d}}{\mathit{\boldsymbol{\omega }}_i} \end{array} $

(2)

则式(1)可写为

$ {L_o}(x, {\mathit{\boldsymbol{\omega }}_o}) = {L_e}(x, {\mathit{\boldsymbol{\omega }}_o}) + {L_r}(x, {\mathit{\boldsymbol{\omega }}_o}) $

(3)

全局光照方程表达的是光能在全场景中传播流动时，在物体表面的每个点$x$处所满足的能量关系，在式(3)中，${L_o}(x, {\mathit{\boldsymbol{\omega }}_o})$表示在点$x$处沿出射方向${\mathit{\boldsymbol{\omega }}_o}$出射的光能总量，也是该方程所要求解的最终量，${L_o}(x, {\mathit{\boldsymbol{\omega }}_o})$是在点$x$处沿出射方向${\mathit{\boldsymbol{\omega }}_o}$自身发射的光能总量，${L_r}(x, {\mathit{\boldsymbol{\omega }}_o})$表示在点$x$处沿着方向${\mathit{\boldsymbol{\omega }}_o}$该物体所反射的环境光照强度。

式(1)反映出的思想即为：物体上任一点沿指定出射方向射出的光子能量，不仅取决于该物体本身在该点处沿该方向射出的能量，也取决于环境光对该物体的照明贡献。

针对${L_r}$项，由式(2)的组成特点可分为以下几个关键参数：

1) 反射系数项${f_r}$；

2) 沿任意方向入射的光照强度${L_i}$；

3) 漫反射余弦项$({\mathit{\boldsymbol{\omega }}_i}\cdot\mathit{\boldsymbol{n}})$。

${f_r}$即上述所提及的双向反射分布函数，它的含义是从入射方向${\mathit{\boldsymbol{\omega }}_i}$射入物体表面的光子能量与从出射方向${\mathit{\boldsymbol{\omega }}_o}$射出物体表面的光子能量的比例关系。

式(1)的含义为：一个具有一定能量的光束沿着${\mathit{\boldsymbol{\omega }}_i}$射入物体表面的某点$x$，由于物体表面吸收了部分能量，导致能量有部分衰减，这部分衰减后的能量与该物体自身发出的能量(例如该物体为光源，则会自身发出能量，若为普通物体则这部分忽略)的总和会根据该物体表面的材质特点重新分布，并按照一定的规律沿着不同方向出射，其中要求解的就是沿着出射方向${\mathit{\boldsymbol{\omega }}_o}$射出的那部分能量的大小，此即式(1)等号左侧的${L_o}(x, {\mathit{\boldsymbol{\omega }}_o})$。而对于能量在射入时的衰减，则利用${f_r}(x, {\mathit{\boldsymbol{\omega }}_i}, {\mathit{\boldsymbol{\omega }}_o})$来描述。${f_r}(x, {\mathit{\boldsymbol{\omega }}_i}, {\mathit{\boldsymbol{\omega }}_o})$是有关入射点$x$，入射方向${\mathit{\boldsymbol{\omega }}_i}$与出射方向${\mathit{\boldsymbol{\omega }}_o}$的函数。

经典全局光照方程在数学上清晰阐述了光能在整个空间中达到动态平衡时在场景任一点所满足的能量关系，但在实际算法应用中却仍需要进一步细化，以便得到可用的全局光照算法。以下是几个细化上述思路的步骤：

1) 式(1)为了得到统一的方程形式，将光源与普通物体一并而论，实际中应将二者的概念区分开，将光源仅视为光能流动的起始处，而将一般物体视为受光物，其仅能接收来自光源和经由其他物体反射到该物体的能量；

2) 在步骤1)的假设前提下，实质上就依据物体表面是否为发光物进行分类。对于非发光物的表面，其受光情况又可分为两类：一种是受光源直接照射，即从光源发出的光不与场景其他物体交互，直接射入该物体表面的光；一种是受到环境光照射，指的是直接光经过多次场景物体的反射作用后，最终射向该物体表面的光，如图 3所示。

在上述细化步骤完成后，再进行下述简化：

1) 光能在传输时自身携带能量，由于复合色的光含有多个波长的光的能量叠加，为了简化，这里将光的能量视为红绿蓝三原色三分量的叠加，每个分量的能量用浮点数近似表示，则一束光的能量可由一个三分量的向量唯一表示，即$\mathit{\boldsymbol{C}}(R, G, B)$。

2) 一束具有能量${\mathit{\boldsymbol{C}}_1}({R_1}, {G_1}, {B_1})$的光在与物体进行交互时，由于每个物体表面具有漫反射颜色${\mathit{\boldsymbol{C}}_2}({R_2}, {G_2}, {B_2})$，且${R_2}、{G_2}、{B_2} < 1$。该束光在经过表面反射后剩余的能量为$C = {C_1} + {C_2} = ({R_1} + {R_2}, {G_1} + {G_2}, {B_1} + {B_2})$。

3) 一束具有能量${\mathit{\boldsymbol{C}}_1}({R_1}, {G_1}, {B_1})$与另一束具有能量${\mathit{\boldsymbol{C}}_2}({R_2}, {G_2}, {B_2})$的光同时射入一个物体的同一个位置时，该位置接收到的总光强为

$ \mathit{\boldsymbol{C}} = {\mathit{\boldsymbol{C}}_1} + {\mathit{\boldsymbol{C}}_2} = ({R_1} + {R_2}, {G_1} + {G_2}, {B_1} + {B_2}) $

4) 基于基本的光路可逆原理，若一束光从光源发出并经过多次场景物体表面的反射作用最终射入摄影机内从而被感知，则该光的通路也可以反向的从摄影机“发出”，经多次场景物体的反射后射向光源，此即“采样光线”的基本思想。

当细化与简化步骤完成后，再次分析此时的全局光照方程所具有的形式。

如图 4所示，假设在非发光物体表面上，其所受的直接光简写为$\mathit{\boldsymbol{D}}$，环境光简写为$\mathit{\boldsymbol{I}}$，物体表面的漫反射颜色简写为$\mathit{\boldsymbol{R}}$，该表面朝某方向出射并直接射入摄影机的光能为$\mathit{\boldsymbol{C}}$，在排除该物体表面为光源的条件下，有

$ \mathit{\boldsymbol{C}} = (\mathit{\boldsymbol{D}} + \mathit{\boldsymbol{I}})\mathit{\boldsymbol{R}} $

(4)

在实际的路径追踪算法中，由概述部分的内容，应从摄影机处发射“采样光线”，并将该光线与场景求交，在求得与场景的第一个交点后，根据表面材质参数确定出射光的方向，再以该交点处为新的采样光线的原点，以新的方向再次与场景求交，重复此过程直到满足程序规定的最大深度为止，收集采样过程中的信息，得到一个像素在一次迭代过程后的颜色。如图 5所示，摄影机首先发射一条深度为1的采样光线，与场景物体求交后，由交点位置与出射方向确定深度为2的采样光线，然后再次与场景求交。图 5画出了最大采样深度为3的示例。

而对于式(4)中环境光照强度$I$的计算，则与式(4)完全类似，有

$ {I_n} = ({D_{n + 1}} + {I_{n + 1}}){R_{n + 1}} $

(5)

式中，角标$n$的含义是当前采样光线所在的采样深度。可见式(5)是一个递推式。合并式(4)(5)并将式(4)添加角标$n$，可得到方程组

$ \left\{ \begin{array}{l} C = ({D_1} + {I_1}){R_1}\\ {I_n} = ({D_{n + 1}} + {I_{n + 1}}){R_{n + 1}} \end{array} \right. $

(6)

将式(6)展开为

$ \begin{array}{l} \;\;\;\;\;\;C = ({D_1} + {I_1}){R_1} = \\ {\rm{ }}\;\;\;({D_1} + ({D_2} + {I_2}){R_2}){R_1} = \\ \;\;({D_1} + ({D_2} + ({D_3} + {I_3}){R_3}){R_2}){R_1} = \\ \;({D_1} + ({D_2} + ({D_3} + ({D_4} + \\ ( \ldots ({D_n} + {I_n}){R_n} \ldots )){R_4}){R_3}){R_2}){R_1} \end{array} $

可得

$ C = \sum\limits_{i = 1}^n {\left( {\prod\limits_{j = 1}^i {{R_j}} } \right)} {D_i} + \left( {\prod\limits_{i = 1}^n {{R_i}} } \right){I_n} $

(7)

此即全局光照方程的简化形式，式中$n$为最大采样深度，$C, I, R, D$均为前述的三分量颜色向量，其加与乘的计算方式已经在简化步骤的2)与3)中说明。当$n$较大时(一般取812为合适)，忽略第2项值较小的$\left({\prod\limits_{i = 1}^n {{R_i}} } \right){I_n}$，因此得到简化后的$C$的最终表达式为

$ C = \sum\limits_{i = 1}^n {\left( {\prod\limits_{j = 1}^i {{R_j}} } \right)} {D_i} $

(8)

$ C = \mathop {\lim }\limits_{n \to \infty } \frac{1}{n}\sum\limits_{i = 1}^n {{C_n}} $

(9)

若采样点与光源之间存在障碍，则可按照图 6所示方法，将该深度的直接光照强度${D_n}$置为0。

根据式(8)就可以得到单个像素在一次迭代下的采样值。因路径追踪算法的特点之一即是路径的不确定性(图 7)，摄影机向同一位置出射的两条光线在先后采样的过程中出现了截然不同的路径，但由于物体表面参数的原因，使得离开表面的光线方向不定，但却有一定的分布特征，虽然每次采样结束后同一个像素的颜色都不同，但大量采样后，对同一个像素的大量采样值取平均，其平均颜色就会收敛到一个理想值附近，该像素的颜色也就最终得以确定，即式(9)。具体实现算法可参见后面的部分。

在1.1节中已经得到了简化的全局光照方程具有的形式，但并未给出采样光线在遇到物体表面时其下一个深度的采样光线的方向。考虑到一般粗糙物体表面的BRDF分布主要集中在以理想镜面反射方向处，而理想漫反射表面的BRDF则均匀分布在半球面内，根据上述分析，鉴于余弦函数与指数函数的分布特性，提出以下方法，从而可以快速计算出射光线的采样方向，其示意图如图 8所示。

1) 当具有方向$\mathit{\boldsymbol{d}}$的光线射到表面时，计算其镜面反射方向$\mathit{\boldsymbol{r}}$；

2) 设随机变量$σ$服从$\left({0, \frac{{\rm{ \mathsf{ π} }}}{2}} \right)$的均匀分布，$ε$为给定的一个正数，计算权重值为

$ w = {({\rm{cos}}\sigma )^{\frac{1}{{{\varepsilon ^2}}}}} $

3) 在半球面内取均匀分布的随机方向$\mathit{\boldsymbol{t}}$，将方向向量$\mathit{\boldsymbol{t}}$与$\mathit{\boldsymbol{r}}$做权重为$w$的线性插值得到出射方向

$ \mathit{\boldsymbol{o}}{\rm{'}} = w\mathit{\boldsymbol{r}} + \left( {1-w} \right)\mathit{\boldsymbol{t}} $

4) 将$\mathit{\boldsymbol{o}}{\rm{'}}$规范化得到单位向量$\mathit{\boldsymbol{o}}$，则最终得到的$\mathit{\boldsymbol{o}}$为出射光线的方向。

上述方法中的$ε$可称为“聚拢因子”，其大小决定了出射光线$\mathit{\boldsymbol{o}}$与理想镜面反射方向$\mathit{\boldsymbol{r}}$的偏离程度的大小。为避免除法耗时，可以预先计算$\varepsilon ' = \frac{1}{{{\varepsilon ^2}}}$。

在上述算法中，由于$\mathit{\boldsymbol{t}}$的方向是随机给定的，则据此生成的$\mathit{\boldsymbol{o}}$也有一定的随机性，从而模拟了光线在射入粗糙表面后所产生的散射现象。且$\mathit{\boldsymbol{o}}$的分布范围受权重值$w$的控制，从而达到了通过给定$w$来控制表面粗糙度，进而控制$\mathit{\boldsymbol{o}}$的分布情况的目的。

下面通过改变$ε$的值，按上述采样方法，令$\mathit{\boldsymbol{r}} = \mathit{\boldsymbol{n}}$即表面法向向量，随机产生1 000个方向并绘图，得到实际的分布特性如图 9所示。

观察以上$ε$的值与其产生的随机样本内的关系，可以发现当$ε→∞$时产生的随机样本靠近表面法线的概率很高，而当$ε→0$时，随机方向基本与物体表面法矢无关，在上半球面内呈随机分布。当$\mathit{\boldsymbol{r}}$指向上半球面内的任意方向时，该算法能够保证生成的所有样本均分布在上半球面内，从而能够通过控制$ε$的值来控制反射光线的分布特性，间接模拟光线在射向一般漫反射表面时产生的模糊现象, 如图 10所示。实际渲染结果表明，该算法能够较好地模拟现实生活中的大部分表面。

为了使得最终的渲染器能够较为逼真地渲染物体，必须引入分层材质的概念。其基本思路为：视物体表面是由两层组成，一层是表皮层，光子在其上直接进行反射；另一层是内层，该层的特点是较为粗糙，近似认为光子在其上进行的是理想漫反射。光在击中表面时，可根据概率决定其在击中表面时的行为(高光反射或漫反射)。这个模型的典型代表是车漆等材质(图 11)，其表面同时具有粗糙的颗粒感与高光质感。在模拟金属等材质时，也需要利用该种模型以用于更好地控制金属材质的质感(图 11、图 12)。

根据上面所述内容，并由文献[19-20]中所陈述的路径追踪渲染器的实现代码，可以整理归纳得到本文所涉及的路径追踪的主要算法。这里光线$\mathit{\boldsymbol{R}}$被定义为空间中具有发射原点$\mathit{\boldsymbol{o}}$和发射方向$\mathit{\boldsymbol{d}}$的射线；所有颜色均由具有R、G、B三分量的向量表示，且$R、G、B < 1$。并假设该算法在单次迭代后得到对应像素的颜色值为$\mathit{\boldsymbol{C}}$，该算法结束后实质上实现了前文所述的式(8)的计算。向量相乘运算法则定义为

$ \begin{array}{l} \mathit{\boldsymbol{A}}({x_1}, {x_2}, \ldots, {x_n}) \times \mathit{\boldsymbol{B}}({y_1}, {y_2}, \ldots, {y_n}) = \\ \;\;\;\;\;\;\;\;\;({x_1}{y_1}, {x_2}{y_2}, \ldots, {x_n}{y_n}) \end{array} $

主要算法具体如下：

1) 初始化变量。对于屏幕中的每一个像素，计算由摄影机中心向屏幕该像素位置发射的采样光线$\mathit{\boldsymbol{Ray}} = \left({\mathit{\boldsymbol{o}}, \mathit{\boldsymbol{d}}} \right)$；定义颜色积聚量$\mathit{\boldsymbol{A}}$=(0, 0, 0)；定义颜色遮罩值$\mathit{\boldsymbol{M}}$=(1, 1, 1)；定义初始递归深度$i=1$与最大递归深度${i_{\max }}$。

2) 对于当前递归深度$i$，进行以下操作：

(1) 将当前光线$\mathit{\boldsymbol{Ray}} = \left({\mathit{\boldsymbol{o}}, \mathit{\boldsymbol{d}}} \right)$与场景求交，若无交点，则根据当前采样光线的方向$\mathit{\boldsymbol{d}}$对HDR环境贴图进行采样，得到的颜色值记为$\mathit{\boldsymbol{H}}$，令$\mathit{\boldsymbol{C}} = \mathit{\boldsymbol{A}} + \mathit{\boldsymbol{H}}$，算法结束；若有交点，则将得到的交点位置记为$\mathit{\boldsymbol{o}}{\rm{'}}$，并继续进行下一小步的操作；

(2) 获取被射中物体的表面参数(聚拢因子$ε$)与被射中物体的表面漫反射颜色$\mathit{\boldsymbol{R}}$，并根据前文所述的新采样算法，得到反射光线的出射方向$\mathit{\boldsymbol{d}}$，从而确定下个递归深度的采样光线$\mathit{\boldsymbol{Ray}}{\rm{'}} = \left({\mathit{\boldsymbol{o}}{\rm{'}}, \mathit{\boldsymbol{d}}{\rm{'}}} \right)$；

(3) 计算得到该点的直接光照强度$\mathit{\boldsymbol{D}}$；

(4) 令$\mathit{\boldsymbol{M}} = \mathit{\boldsymbol{M}} \times \mathit{\boldsymbol{R}}$，从而更新$\mathit{\boldsymbol{M}}$的值；并令$\mathit{\boldsymbol{A}} = \mathit{\boldsymbol{A}} + \mathit{\boldsymbol{M}} \times \mathit{\boldsymbol{D}}$，更新$\mathit{\boldsymbol{A}}$的值；

(5) 若$i = = {i_{{\rm{max}}}}$，则令$\mathit{\boldsymbol{C}} = \mathit{\boldsymbol{A}}$，算法结束；否则令$i = i + 1$，以便进入下一个递归深度。

3) 重复步骤2)，当算法退出时，得到式(8)中所述$\mathit{\boldsymbol{C}}$值。

上述步骤实质上是式(8)的实现，若要得到该像素的最终颜色，则由式(9)所述，应多次重复步骤1)得到${\mathit{\boldsymbol{C}}_1}, {\mathit{\boldsymbol{C}}_2}, \ldots, {\mathit{\boldsymbol{C}}_n}$。并令$C = \mathop {{\rm{lim}}}\limits_{n \to \infty } \frac{1}{n}\sum\limits_{i = 1}^n {{C_n}} $，随着$n$的增加，最终的$\mathit{\boldsymbol{C}}$会不断趋于稳定。

在图 13的渲染结果中，每组的左图是基于MERL BRDF数据库渲染出的实际参考结果。中图是基于本论文中的方法所渲染出的近似结果，右图是两张图片的差值图像。为保证严谨性且能更好地对比渲染差别，这里统一了摄影机位置、摄影机参数与光源位置，将物体表面的纹理贴图信息去除，保留主要颜色，并调节适当参数(聚拢因子，表面漫反射颜色等)使模拟结果尽可能贴近理想图像。其中MERL数据库中每种材质所测得的BRDF数据量均为34 MB，而使用本文方法可以使其降为仅几个浮点数的存储量，其存储开销可忽略不计。实验定量数据见表 1，其中相对误差估计的值是对差值图像内所有像素颜色取统计平均后，与原图像素颜色均值的比值。从图 13的差值图像可以看出，该算法在模拟大部分材质表面时，具有良好的拟真度。但其精简后的材质参数在模拟现实材质时仍具有不理想的地方，具体体现即为以下几个方面。

1) 对于菲涅尔反射特性的模拟。实际的材质表面，在光子以大角度入射时会有高反射率的特点，而该算法渲染出的实际结果却不能良好的体现该特性，如图 14(a)中给出的细节对比图。

2) 对于材质细节的表现能力。如图 14(b)，在物体表面反射率产生突变时，模拟结果就产生了较大的误差，这是因为采样过程中参数过少，所能模拟的分布特性较为单一所造成的。

图 15与图 16是利用本文算法编制出的程序的实际渲染结果。图 15说明了当参数$ε$连续发生变化时，被渲染的物体表面可由理想漫反射到理想镜面反射之间光滑过渡，从而统一了漫反射、高光反射与镜面反射三者的采样形式。图 16模拟的材质为镀金材质，聚拢因子的值已在图中标注。可对比发现其光泽度不断变高，且具有良好的视觉体验。

在上述算法中，通过入射光的方向和入射点处的表面法向向量计算出镜面反射向量，再根据表面法向所确定的半球面中随机均匀生成一个漫反射方向向量，通过随机变量与函数结合，生成权重，根据权重将漫反射方向与镜面反射方向线性插值后规范化，得到最终光线的出射方向，根据该思路，可以得到基于该算法的几个扩展，同时也可作为对当前算法局限性的补偿方案。

1) 在生成符合某种特性的权重分布时，引入了描述物体表面光滑程度的“聚拢因子”$ε$，用来控制光线出射方向相对于理想镜面反射方向的散射程度。在上述算法中，$ε$一旦确定后即为常数，这在一定程度上限制了算法的灵活性。如果将$ε$视为入射角的函数即$\varepsilon = \varepsilon ({\theta _i})$，就可以模拟更加复杂的材质表面，$\varepsilon ({\theta _i})$的具体形式应在达到最佳视觉拟真效果的前提下视具体情况而定。

2) 由渲染结果中对该快速算法失真性的讨论，可以得出如下的解决方案：因新的采样方式中将原本几十兆的存储数据量缩减为几个浮点数来存储，其存储开销可以忽略不计，所以可采用多层的材质表面拟合出原来的真实表面。例如，若要模拟一个既具有明显高光且又在非高光区域有粗糙漫反射特性的材质表面，可设置两层参数，其中一层的聚拢因子值较高，用于模拟高光特性，另一层的聚拢因子值较低，用于模拟粗糙特性，再将这两层按照一定的权重比例进行混合，当光子射入该混合参数表面时，可按照之前设定的概率，随机进入多层材质中的一层，再按照本文方法进行快速采样，即模拟了混合参数表面的材质特性。其具体的参数设置等细节可作为以后的研究方向继续展开深入研究。

3) 本文方法可类比到具有透射特性的表面的渲染算法中。即计算理想折射光的方向，再在透射物体的内表面处随机生成方向，将折射光方向与随机方向进行线性插值，这样即可模拟类似毛玻璃的模糊透射效果。

由渲染结果中所示图片以及实验结果可以看出，本文所述的反射光线快速采样方法能够较为真实地模拟现实生活中的各向同性表面，且通过控制聚拢因子$ε$可以模拟具有不同粗糙度的材质表面，$ε$的值可使得表面反射特性从镜面反射到漫反射连续变化，在需要对场景进行不失真实度的快速可视化渲染时，可借助该方法来完成路径追踪中对出射光线方向的快速采样。本文方法在模拟金属表面、漫反射表面和镜面等材质时具有高拟真特性，且硬件存储开销极小。在模拟其他可能出现的材质表面时，需对出射光线采样算法稍微做出更改，但主要思路与文中的思路相同。日后的研究方向为：利用文中所述的快速采样方法，解决对复杂材质表面进行多层拟合时的参数设置问题；以及利用该算法对各向异性表面进行渲染时所需要做出的改进。