CN104361624A - 一种计算机三维模型中全局光照的渲染方法 - Google Patents

Abstract

本发明公开了一种计算机三维模型中全局光照的渲染方法，包括了以下步骤：1)读入模型文件，并将其点和三角面信息存储于kd-tree划分的数据结构中。2)通过采样生成若干条从光源出发的光源路径。3)将最终渲染的图像平面分批采样，每次采样一部分点，在这些像素点处生成从视点出发的视点路径，并从2)中选择一条光源路径，采用双向路径跟踪的路径连接方法与之做连接计算得到当前像素的颜色值，重复该步骤直到整个图像上的像素点都被采样到。4)重复之前的2)和3)直到整个渲染结果在可接受的范围内。

Description

一种计算机三维模型中全局光照的渲染方法

技术领域

本发明涉及到一种计算机三维模型中全局光照的渲染方法，属于计算机图形学等领域。

背景技术

近年来，计算机技术飞速发展，而应用于现代娱乐生活的游戏，电影中的图形学也是高速发展。在电影中，离线渲染技术中的全局光照渲染已经可以制造出以假乱真的效果。全局光照计算是计算机图形学中一个核心的组成部分，它的任务是使预先组织好的三维场景生成完全真实感的图像。为了达到这一目的，我们必须模拟光线在场景中进行传播时所发生的各种物理现象，如相互的反射，焦散，阴影等。这不但需要我们对整个场景中的各个物体的材质以及几何特性等有一个正确的描述，还需要我们对由于多次反射和折射所带来的一个无穷维的积分进行求解。由于全局光照的数学模型过于复杂，在巨大的计算量面前，算法的时间成本很难降下来。为了提高渲染算法效率，降低时间成本，很多全局光照算法被提出来。

在实际应用中，当前传统的方法是使用基于点的全局光照算法或者路径跟踪一类的算法，基于点的全局光照能够很快速地处理漫反射场景，但是对于镜面，透明面的场景没有很好的处理方法；另一方面，路径跟踪(包括双向路径跟踪)一系列的方法能很好地处理镜面，透明面等各种材质，但其算法效率比较低。

发明内容

发明目的：本发明所要解决的技术问题是针对现有方法效率的不足，提供一种更为高效的计算机三维模型中全局光照的渲染方法。

技术方案：本发明公开了一种计算机三维模型中全局光照的渲染方法，该方法的特征在于能够在较短的时间内对一个给定的三维场景绘制出具有真实感的图像，并且能够在最终结果出来之前多次看到其计算的中间结果，以便更好地调节渲染参数。具体包括以下步骤：

1.加载模型：读入三维模型场景的点坐标，三角形序列，点的法向量，点的纹理；用kd-tree(k-dimensional树的简称，是一种分割k维数据空间的数据结构)对三维模型场景空间进行划分，并将三角面存放于kd-tree中，同时从与该三维模型场景对应的材质文件中读入材质属性；

步骤1-1，加载模型文件：读入三维模型场景文件，获取点坐标v，点法向量vn，点纹理vt，三角形f，以及材质对应项(描述对应于材质文件中的材质名称)；其中点坐标v和点法向量vn是一个三维向量，代表点在空间的位置，点纹理vt是一个二维向量，代表该点的纹理在纹理贴图中的坐标，三角形f是一个三维向量，包含了一个三角面中三个点的信息；

步骤1-2，加载材质文件：加载与步骤1-1中的三维模型场景对应的材质文件，读取每个面所用材质属性，包括漫反射系数，镜面反射系数，透明度以及纹理路径；

步骤1-3，计算包围盒：首先计算出整个模型空间的最大范围，用一个轴对齐的长方体将整个模型包围起来。这个长方体的计算公式如下：

X_max＝max{v0_x,v1_x,…vn_x},X_min＝min{v0_x,v1_x,…vn_x},

Y_max＝max{v0_y,v1_y,…vn_y},Y_min＝min{v0_y,v1_y,…vn_y},

Z_max＝max{v0_z,v1_z,…vn_z},Z_min＝min{v0_z,v1_z,…vn_z}.

其中长方体在x，y，z方向上的坐标为X,Y,Z，而X_max和X_min分别代表长方体在x方向上的最大值和最小值，Y_max和Y_min分别代表长方体在y轴方向上的最大值和最小值，Z_max和Z_min分别代表长方体在z轴方向上的最大值和最小值。vi_x，vi_y，vi_z分别表示点集v中第i的点的x，y，z坐标，其中i＝0,1,……,n，(n是点集v所包含的点的数量)。

步骤1-4，建立kd-tree并存储数据：将整个长方体空间利用kd-tree划分，并将之前读取三维模型场景文件得到的三角形数据保存于kd-tree中。kd-tree在k维空间中用轴对齐的子平面一次次地划分整个空间，最后得到一个二叉树的数据结构。Kd-tree的划分基于SAH策略，SAH策略全称Surface Area Heuristic策略，是基于表面积启发的策略，其核心思想是在划分子空间的时候尽量以最小的子空间表面积包围更多的三角形。假设对于空间S选取剖分面P，把该空间划分为左右两个子空间S_L和S_R，并且同时把场景中的三角形面片划分成两个集合N_L和N_R，上述的SAH优化函数公式为：

${\mathrm{SAH}}_{\cos t}\left(\mathrm{node}\right)=C_{\mathrm{walk}}+\frac{\mathrm{SA}\left(S_L\right)}{\mathrm{SA}\left(S\right)}{N}_L{C}_{\mathrm{hit}}+\frac{\mathrm{SA}\left(S_R\right)}{\mathrm{SA}\left(S\right)}{N}_R{C}_{\mathrm{hit}}$

其中node是指当前空间节点，SA(S)表示空间S的表面积，C_hit表示光线与每个三角形的求交代价。C_walk表示光线遍历当前node的代价，在划分每个节点的时候采用使得SAH_cost(node)最小的划分方式。且在满足以下任一条件的时候终止划分，包括：

1)当前节点的深度超出一个预先设定好的阈值Deep(1≤Deep<∞,且为整数)；

2)当前节点中的三角形个数少于一个预先设定的阈值TN(0≤TN<20)；最终得到一个包含了所有场景三角形信息的kd-tree数据结构，这个kd-tree只有叶子节点有三角形数据，非叶子节点的三角形数据为空。

2.采样光源：使用与双向路径跟踪中采样光源一样的方法(详情请见论文<ROBUSTMONTE CARLO METHODSFOR LIGHT TRANSPORT SIMULATION>Eric Veach,chapter 10)进行采样得到若干条光源路径。在传统的双向路径跟踪算法中，在每次计算一个像素时所产生的光源路径都是互相独立的，这样使得基于CPU的双向路径跟踪计算速度变得很慢，而本发明的方法的不同之处在于不同像素之间重复利用光源路径，以减少光源路径的生成的数量，达到提升效率的目的。虽然这种相关性的采样方式会使得整个算法的方差增大，但是仍能够保持其无偏性，而节省下来的时间完全能够弥补方差增大所带来的噪声的影响，使得算法的总体效率提升。由于产生的光源路径数量较少，所以光源路径的分布情况对噪声的影响会更明显，本发明使用了分层采样的方法使得这些光源路径在路径空间(所有采样路径的集合)中分布更加均匀，这样也可以使最终的求解收敛速度更快。具体包括以下步骤：

步骤2-1，从光源上采样路径的起始点：从光源上随机采样一个点P₀作为光源路径的起始点(因为一条光源路径可以看成由若干条直线相连而成，在采样这样一条光源路径的时候，我们相当于是在采样若干个点，然后将这些点按顺序相连，得到光源路径)，如果是点光源，那么路径起始点就是光源位置，如果是面光源或者是光源体，那么路径的起始点是这个面光源或者光源体上的随机采样点；

步骤2-2，从路径起点P₀随机产生路径：依次采样路径的下一个点，直到光线射出整个三维模型场景，或者被三维模型场景表面吸收，或者是路径长度到达某个阈值d。产生路径的方式类同路径跟踪算法，即光线在场景空间传播，遇镜面发生镜面反射，遇漫反射面随机在表面法线方向上半球面取下一反射方向，遇透明表面发生折射(详细算法过程参见<The rendering equation>ACM SIGGRAPH Computer Graphics Volume 20Issue 4,Aug 1986,Pages 143-150)。这时将产生出一条采样得到的光源路径，计算这条光源路径的贡献度。其计算公式为：

${\mathrm{\&alpha;}}_0^L=1,$

${\mathrm{\&alpha;}}_1^L=\frac{L_e^{\left(0\right)}\left(y_0\right)}{P_A\left(y_0\right)},$

${\mathrm{\&alpha;}}_r^L=\frac{f_s(y_{r-3}\&RightArrow;y_{r-2}\&RightArrow;y_{r-1})}{P_{{\mathrm{\&sigma;}}^{\&perp;}}(y_{r-2}\&RightArrow;y_{r-1})}{\mathrm{\&alpha;}}_{r-1}^L,(r\&GreaterEqual;2)$

其中表示光源路径上光源位置的颜色贡献度，表示光源路径上第1个点位置的颜色贡献度，表示光源路径上第r个点位置的颜色贡献度，P_A(y₀)表示光源上采样到y₀点的概率，表示光源在y₀的亮度值，f_s(y_r-3→y_r-2→y_r-1)则表示在已得到点y_r-2的情况下，采样到点y_r-1的概率。其中而P_σ(y_r-2→y_r-1)表示采样y_r-2→y_r-1方向的概率，θ表示采样方向与y_r-2点的表面法向量的夹角。根据分层采样的原理(其原理请参见<Advanced Global Illumination>SecondEdition,Page 71)，产生出若干条这样的光源路径。

3.采样视点并计算像素颜色值：对最终渲染结果的图像平面进行划分，每次根据该划分选出一组像素。对平面划分的意义在于，可以使每次选择到的像素能够均匀地分布到图像平面的各个位置，而如果每一组选到的像素点过于集中的话，就会导致最终结果的噪声增大。像素的选择准则是：1)每一批选择的像素点应当尽量地均匀分布在整个像素平面；2)每一批选择的像素点应当尽量少地存在像素相邻的情况。选定好一批像素点后，利用透视原理，以视点到这些像素的方向作为初始方向采样路径，得到一批视点路径，从步骤2中选择一条光源路径，与这些视点路径相连接，得到完整的从光源到视点的路径。计算整个路径对像素的颜色贡献值。重复这个过程直到对图像平面的采样点覆盖整个图像平面。最终会得到一次迭代后的渲染结果图像。具体包括以下步骤：

步骤3-1，划分图像平面，假设图像平面一共有N个像素，我们每次不重复地取出M个像素点，M<N，用以产生视点路径，总共分轮取完，不同轮次取出的像素点也不能有重复的。具体的像素取法如下：

3-1-1.以A*B的block(一个block是指最终渲染结果图像平面上的一个长方形区域)划分整个像平面，假设图像平面长L，高为H，那么它将被分成份，其中表示X向上取整。在选择A与B的值的时候应当尽量保证L能整除A，H能够整除B，这样可以使算法效率更高；

3-1-2.每一批像素的选取是从每个block中分别随机选出一个像素，但是在之前轮已经选择过的像素不能再被选取，一共M个像素。为了不用每一次都进行这个随机，我们可以预先生成一组随机序列，即假设一个block中有个像素，那么我们就生成一个1-C一共C个数的随机排列，每一轮选取的时候按照这个排列来选择。

步骤3-2，从视点向选中的像素点发射光线，作为初始路径，依次产生出视点路径，路径的产生方式与步骤2-2中产生光源路径的方式类似，产生出的视点路径的贡献度计算公式为：

${\mathrm{\&alpha;}}_0^E=1,$

${\mathrm{\&alpha;}}_1^E=\frac{W_e^{\left(0\right)}\left(z_0\right)}{P_A\left(z_0\right)},$

${\mathrm{\&alpha;}}_{r^{\&prime;}}^E=\frac{f_s(z_{r^{\&prime;}-3}\&RightArrow;z_{r^{\&prime;}-2}\&RightArrow;z_{r^{\&prime;}-1})}{P_{{\mathrm{\&sigma;}}^{\&perp;}}(z_{r^{\&prime;}-2}\&RightArrow;z_{r^{\&prime;}-1})}{\mathrm{\&alpha;}}_{r^{\&prime;}-1}^E,(r^{\&prime;}\&GreaterEqual;2)$

其中表示视点路径上视点位置的颜色贡献度，表示视点路径上第1个点位置的颜色贡献度，表示视点路径上第r′个点位置的颜色贡献度，PA(z0)表示在视点棱镜表面上采样到z0点的概率(如果视点是一个理想的点，那么该值为1)，表示z0点在视点棱镜表面的权值；

步骤3-3，路径连接：对每一个选中的像素点都产生一条视点路径后，从步骤2-2中选择一条光源路径与它们连接，即连接视点路径和光源路径上的最后一个点，使之成为一条完整的从光源到视点的路径。依照如下公式求得该完整路径的贡献度：其中i是指光源路径的长度，j是指视点路径的长度，与可分别由步骤2-2以及步骤3-2中的公式得到，c_i,j的计算公式如下：

c_0，j＝L_e(z_j-1→z_j-2)，c_i，0＝W_e(y_i-2→y_i-1)，

$c_{i,j}=f_i(y_{i-2}\&RightArrow;y_{i-1}\&RightArrow;z_{j-1})G(y_{i-1}\&LeftRightArrow;z_{j-1})f_i(z_{j-2}\&RightArrow;z_{j-1}\&RightArrow;y_{i-1}),i,j>0$

其中i是指光源路径的长度，j是指视点路径的长度，G是形状因子函数，定义为：其中cosθ_o表示视点路径和光源路径上最后一个点相连接得到的直线L与视点路径上最后一个点所在的三角形表面法向量的余弦值，cosθ′表示该直线L与光源路径上最后一个点所在的三角形表面法向量的余弦值；

V则是可见性函数，如果点x和x′在场景中可见，那么反之 $V(x\&LeftRightArrow;x^{\&prime;})=0.$

步骤3-4，计算加权贡献度：加权贡献度的计算公式为：

其中路径权值而比值p_a/p_i可以利用以下公式求得：

$\frac{p_1}{p_0}=\frac{P_A\left(x_0\right)}{P_{{\mathrm{\&sigma;}}^{\&perp;}}(x_1\&RightArrow;x_0)G(x_1\&LeftRightArrow;x_0)},$

$\frac{p_{a+1}}{p_a}=\frac{P_{{\mathrm{\&sigma;}}^{\&perp;}}(x_{a-1}\&RightArrow;x_a)G(x_{a-1}\&LeftRightArrow;x_a)}{P_{{\mathrm{\&sigma;}}^{\&perp;}}(x_{a+1}\&RightArrow;x_a)G(x_{a+1}\&LeftRightArrow;x_a)},(0<a<i)$

$\frac{p_{i+1}}{p_i}=\frac{P_{{\mathrm{\&sigma;}}^{\&perp;}}(x_{i-1}\&RightArrow;x_i)G(x_{i-1}\&LeftRightArrow;x_i)}{P_A\left(x_i\right)},$

步骤3-5，计算整个图像平面的像素颜色值：将步骤3-4中求出的路径加权贡献度作为当前像素颜色值存放起来，然后重复进行步骤3-1到3-4，直到整个图像平面上的像素的颜色值都被计算过一次为止，注意步骤3-3中选取光源路径的时候，之前计算中已经选择过的光源路径不能被重复选取。正如步骤3-2-2中所说，我们可以事先对每一个block生成这样一组随机的序列，这样可以避免每一轮选择像素的时候都要随机一次。此步骤最终会生成一张中间渲染图像。

4.产生最终结果图像：使用渐进式的渲染方式，重复步骤2和步骤3直到渲染结果的误差小于一个预先设定的值为止。这一步中涉及到误差判断，但是在实际应用中，一般是使用渐进式渲染的轮数，或者总时间来作为标准，在本方法中，使用了渲染轮数作为结束条件，每一轮就相当于进行一次步骤2到步骤3的操作，当总的渲染轮数达到某一阈值的时候，结束渲染得到最终结果图像，具体包括以下步骤：

步骤4-1，重复步骤2和步骤3，每次计算一整张新的平面颜色值，并与之前的结果相融合。假设当前是第T轮计算，而前面T-1轮中得到的某像素点Pi的颜色值为f_T-1(Pi)＝a，第T轮得到的加权路径贡献度为Co，那么第T轮结束后，该点的颜色值被更新为：当T＝1的时候，f_T(Pi)＝Co。

步骤4-2，评价渲染结果：用渲染轮数，或者时间来进行评价，当渲染总轮数达到某个阈值N(根据场景的复杂度来定，通常设定范围为:渲染轮数20-5000,渲染时间10分钟-10小时)的时候，我们认为渲染结果达标。结束渲染，得到最终结果图像。

在本发明中的全局光照算法属于路径跟踪系列的方法，但是由于本方法重用了光路，并且采用渐进式的渲染框架，所以比起传统的方法效率更高，以及能在运行过程中看到粗略的效果图，可以更快地看到渲染的大体效果。

附图说明

图1为本发明方法的基本流程图。

图2为整个算法计算单个像素颜色值的时候路径连接示意图。

图3为一个典型的3维kd-tree空间结构图。

图4为对平面的三种采样方式。

图5为cornell box场景计算一个循环后的结果。

图6为cornell box场景计算八个循环后的结果。

图7为cornell box场景计算三十二个循环后的结果。

图8为本发明方法与标准双向路径跟踪方法在sponza场景下的结果对比。

图9为本发明方法渲染的sponza场景与reference的差异放大5倍后的结果。

图10为标准bpt方法渲染的sponza场景与reference的差异放大5倍后的结果。

具体实施方式

下面结合附图和具体实施方式对本发明做更进一步的具体说明。

本方法的流程图如图1所示，是一个迭代式的过程：首先是场景模型文件导入，并采用kd-tree对整个场景进行划分，合理的划分可以大大减少计算的时间；紧接着是从光源采样产生光源路径，首先需要从光源上随机采样一个点，把这个点作为路径初始点，然后通过后续一系列的采样产生若干条子路径，最后得到一整条光源路径，重复这个过程直到得到若干条光源路径；然后需要从视点采样路径，但是由于上一步的光源路径是要被重用的，这样会由于采样的相关性导致噪声的增大，所以需要对最终渲染结果的图像平面进行划分，每次选择不同的像素采样路径并与光源路径进行连接计算，以加快该方法的收敛速度。从视点生成路径的时候采用的办法与从光源产生光源路径相同。每产生一批视点路径，就不重复地从光源路径集合中选取一条光源路径与它们连接成为一条完整的从光源到视点的路径，作为该像素点最终颜色的一个估计值，在多次估计后，得到的结果就能更加接近该像素的真实颜色值。当整个图像平面都完成一次颜色值的估计后，一次迭代过程结束。最后一直重复这个计算的过程，经过多次计算估值后，可以让计算得到的像素颜色值足够接近真实值，这时我们认为渲染结果达标，输出计算结果图像。

具体地说，如图1所示，本发明公开了一种计算机三维模型中全局光照的渲染方法，主要包括以下几个步骤：

步骤1，对于一个场景渲染的第一步是需要得到这个场景的三维信息，即场景文件。本发明用的是通用的.obj模型(OBJ格式文件是Alias|Wavefront公司为它的一套基于工作站的3D建模和动画软件"Advanced Visualizer"开发的一种标准3D模型文件格式)，需要读入模型的点坐标，三角形序列，点的法向量，点的纹理。并且用kd-tree进行空间划分，将三角面存放于kd-tree中，同时从与模型文件对应的.mtl文件(Material LibraryFile的缩写，描述的是物体的材质信息)中读入材质属性。

步骤2，采样光源：使用与双向路径跟踪方法中采样光源一样的方法(详情请见论文<ROBUST MONTE CARLO METHODSFOR LIGHT TRANSPORTSIMULATION>Eric Veach,chapter 10)进行采样得到若干条光源路径。在传统的双向路径跟踪算法中，在每次计算一个像素时所产生的光源路径都是互相独立的，这样使得基于CPU的双向路径跟踪算法的计算速度变得很慢，而本发明的方法的不同之处在于不同像素之间重复利用光源路径，以减少光源路径生成的数量，达到提升效率的目的。虽然这种相关性的采样方式会使得整个算法的方差增大，但是仍能够保持其无偏性，而节省下来的时间完全能够弥补方差增大所带来的噪声的影响，使得算法的总体效率提升。由于产生的光源路径数量较少，所以光源路径的分布情况对噪声的影响会更明显，本发明使用了分层采样的方法使得这些光源路径在路径空间中分布更加均匀，这样也可以使最终的求解收敛速度更快。

步骤3，采样视点并计算像素颜色值：对最终渲染结果的图像平面进行划分，每次选出一组像素。对平面划分的意义在于，避免每次选择到的像素能够均匀地分布到图像平面的各个位置，而如果每一组选到的像素点过于集中的话，就会导致最终结果的噪声增大。像素的选择准则是：1)每一组选择的像素点应当尽量地均匀分布在整个像素平面；2)每一组选择的像素点应当尽量少地存在像素相邻的情况。选定好一组像素点后，利用透视原理，以视点到这些像素的方向作为初始方向采样路径，得到一组视点路径，从步骤2中选择一条光源路径，与这些视点路径相连接，得到完整的从光源到视点的路径。计算整个路径对像素的颜色贡献值。步骤2与3中产生光源路径与视点路径并相连的这个过程如图2所示。重复这个过程直到对图像平面的采样点覆盖整个图像平面。最终会得到一次迭代后的渲染结果图像。

步骤4，产生最终结果图像：使用渐进式的渲染方式，重复步骤2和步骤3直到渲染结果的误差小于一个预先设定的值为止。这一步中涉及到误差判断，但是在实际应用中，一般是使用渐进式渲染的轮数，或者总时间来作为标准，在本方法中，使用了渲染轮数作为结束条件，每一轮就相当于进行一次步骤2到步骤3的操作，当总的渲染轮数达到某一阈值(阈值根据场景的复杂度不同而不同，范围在)的时候，结束渲染得到最终结果图像。

步骤1，加载模型的详细步骤如下：

步骤1-1，加载模型文件：读入.obj格式文件，获取点坐标信息v，点法向量信息vn，点纹理数据vt，三角形信息f，以及材质对应项(描述对应于材质文件中的材质名称)。其中点坐标v和点法向量vn是一个三维向量，代表点在空间的位置，点的纹理数据vt是一个二维向量，代表该点的纹理在纹理贴图中的坐标，f是一个三维向量，描述了一个三角面中三个点在点坐标集v中的索引值。

步骤1-2，加载材质文件：读入与步骤1-1中的模型场景对应的.mtl文件，读取每个面所用材质属性，包括漫发射系数kd，镜面反射系数ks，透明度kt，以及表面的纹理贴图文件路径。

步骤1-3，计算包围盒：首先计算出整个模型空间的最大范围，用一个轴对齐的长方体将整个模型包围起来。这个长方体的计算公式如下：

X_max＝max{v0_x,v1_x,…vn_x},X_min＝min{v0_x,v1_x,…vn_x},

Y_max＝max{v0_y,v1_y,…vn_y},Y_min＝min{v0_y,v1_y,…vn_y},

Z_max＝max{v0_z,v1_z,…vn_z},Z_min＝min{v0_z,v1_z,…vn_z}.

其中长方体在x，y，z方向上的坐标为X,Y,Z，而X_max和X_min分别代表长方体在x方向上的最大值和最小值，Y_max和Y_min分别代表长方体在y轴方向上的最大值和最小值，Z_max和Z_min分别代表长方体在z轴方向上的最大值和最小值。vi_x，vi_y，vi_z分别点集v中第i的点的x，y，z坐标，其中i＝0,1,……,n，(n是点集v包含的点的数量)。

步骤1-4，建立kd-tree并存储数据：将整个长方体空间利用kd-tree划分，并将之前读取.obj文件得到的三角形数据保存于kd-tree中。Kd-tree是指在k维空间中用轴对齐的子平面一次次地划分整个空间，最后得到一个二叉树的数据结构。Kd-tree的划分基于SAH策略，SAH(Surface Area Heuristic)策略的基本思想是竟可能把更多的三角面包含到更小的子空间中，这样可以使得ray在遍历整个kd-tree的时候有更高的效率。假设对于空间S选取剖分面P，把该空间划分为左右两个子空间S_L和S_R，并且同时把场景中的所有三角形面片划分成两个集合N_L和N_R，上述的SAH优化函数便可以写成公式：

${\mathrm{SAH}}_{\cos t}\left(\mathrm{node}\right)=C_{\mathrm{walk}}+\frac{\mathrm{SA}\left(S_L\right)}{\mathrm{SA}\left(S\right)}{N}_L{C}_{\mathrm{hit}}+\frac{\mathrm{SA}\left(S_R\right)}{\mathrm{SA}\left(S\right)}{N}_R{C}_{\mathrm{hit}}$

其中node是指当前空间节点，SA(S)表示空间V的表面积，C_hit表示光线与每个三角形的求交代价。C_walk表示光线遍历当前node的代价，在划分每个节点的时候采用使得SAH_cost(node)最小的划分方式。且在满足以下任一条件的时候终止划分，包括：

1)当前节点的深度超出一个预先设定的阈值Deep(1≤Deep<∞,且为整数)；

2)当前节点中的三角形个数少于一个预先设定的阈值TN(0≤TN<20)；

最终得到一个包含了所有场景三角形信息的kd-tree数据结构，这个kd-tree只有叶子节点有三角形数据，非叶子节点的三角形数据为空。一个典型的3维kd-tree如图3所示。

步骤2，采样光源的详细步骤如下：

步骤2-1，从光源上采样路径的起始点：从光源上随机采样一个点P₀作为光源路径起始的点，如果光源是点光源，那么路径起始点就是光源位置，如果是面光源或者是光源体，那么路径的起始点是这个面光源或光源体上的随机采样点。

步骤2-2，从路径起点P₀随机产生路径：依次采样路径的下一个点，直到光线射出整个场景，或者被场景表面吸收，或者是路径长度到达某个阈值值d。产生路径的方式同路径跟踪算法，即光线在空间传播，遇镜面发生镜面反射，遇漫反射面随机在表面法线方向上半球面取下一反射方向，遇透明表面发生折射(详细算法过程参见<Therendering equation>ACM SIGGRAPH Computer Graphics Volume 20Issue 4,Aug 1986,Pages 143-150)。这时将产生出一条光源路径，计算这条光源路径的贡献度。其计算公式为：

${\mathrm{\&alpha;}}_0^L=1,$

${\mathrm{\&alpha;}}_1^L=\frac{L_e^{\left(0\right)}\left(y_0\right)}{P_A\left(y_0\right)},$

${\mathrm{\&alpha;}}_r^L=\frac{f_s(y_{r-3}\&RightArrow;y_{r-2}\&RightArrow;y_{r-1})}{P_{{\mathrm{\&sigma;}}^{\&perp;}}(y_{r-2}\&RightArrow;y_{r-1})}{\mathrm{\&alpha;}}_{r-1}^L,(r\&GreaterEqual;2)$

其中表示光源路径上光源位置的颜色贡献度，表示光源路径上第1个点位置的颜色贡献度，表示光源路径上第r个点位置的颜色贡献度，P_A(y₀)表示光源上采样到y₀点的概率，表示光源在y₀的亮度值，f_s(y_r-3→y_r-2→y_r-1)则表示在已得到点y_r-2的情况下，采样到点y_r-1的概率。其中而P_σ(y_r-2→y_r-1)表示采样y_r-2→y_r-1方向的概率，θ表示采样方向与y_r-2点的表面法向量的夹角。最后根据分层采样的思想，产生出若干条这样的光源路径。分层采样的基本思想是将采样空间均分成N份，然后再在这均分的每一份中随机采样一次，一共得到N个样本。关于分层采样在图形学中的应用，更详细的内容请参见<AdvancedGlobal Illumination>Second Edition,Page 71。

步骤3，采样视点并计算像素颜色值详细步骤如下：

步骤3-1，划分图像平面，假设图像平面一共有N个像素，我们每次不重复地取出M个像素点，用以产生视点路径，总共分轮取完，不同轮次取出的像素点也不能有重复的，如图4的上图。而同一轮所采样到的像素中也要尽可能地不包含有相邻的像素，因为同一轮中过多的像素挤到一起会对最后结果造成不好的影响(如果不划分图像平面，就可能出现像素点分布不均匀，以及堆积的情况，如图4的中图，同时也要保证随机性，不能规则地在每个块中按顺序取点，如图4的下图，不然最后的结果会有带样式的噪声出现(pattern noise)影响视觉感受，具体的像素取法如下：

3-1-1.以A*B的block(一个block是指最终渲染结果图像平面上的一个长方形区域)划分整个像平面，假设图像平面长L，高为H，那么它将被分成份，其中表示X向上取整。尽量保证L能整除A，H能够整除B，这样可以使算法效率更高；

3-1-2.每一批像素的选取是从每个block中分别随机选出一个像素，但是在之前轮已经选择过的像素不能再被选取，一共M个像素。为了不用每一次都进行这个随机，我们可以预先生成一组随机序列，即假设一个block中有个像素，那么我们就生成一个1-C一共C个数的随机排列，每一轮选取的时候按照这个排列来选择。随机序列的生成基本思想是数组乱序，即对一个有序的数组采用一定的操作产生一个乱序数组，但要保证每一种排列得到的概率都一样，本文采用的方法基于《计算机程序设计艺术》上提供的数组乱序的算法，其C语言实现代码如下：

该方法能够保证所产生的随机序列是在整个数组排列空间中均匀分布的，其证明过程具体请参考《计算机程序设计艺术》第35章。该算法的时间复杂度为O(n)，本文改进了该方法，在固定划分好图像平面后，只生成一次随机序列，并保存下来，每一轮都采用这个随机序列，但是每一轮的开始位置不一样。

步骤3-2，从视点向选中的像素点发射光线，作为初始路径，依次产生出视点路径，路径的产生方式与步骤2-2中产生光源路径的方式类似，产生出的视点路径的贡献度计算公式为：

${\mathrm{\&alpha;}}_0^E=1,$

${\mathrm{\&alpha;}}_1^E=\frac{W_e^{\left(0\right)}\left(z_0\right)}{P_A\left(z_0\right)},$

${\mathrm{\&alpha;}}_{r^{\&prime;}}^E=\frac{f_s(z_{r^{\&prime;}-3}\&RightArrow;z_{r^{\&prime;}-2}\&RightArrow;z_{r^{\&prime;}-1})}{P_{{\mathrm{\&sigma;}}^{\&perp;}}(z_{r^{\&prime;}-2}\&RightArrow;z_{r^{\&prime;}-1})}{\mathrm{\&alpha;}}_{r^{\&prime;}-1}^E,(r^{\&prime;}\&GreaterEqual;2)$

其中表示视点路径上视点位置的颜色贡献度，表示视点路径上第1个点位置的颜色贡献度，表示视点路径上第r′个点位置的颜色贡献度，PA(z0)表示在视点棱镜表面上采样到z0点的概率(如果视点是一个理想的点，那么该值为1)，表示z0点在视点棱镜表面的权值；

步骤3-3，路径连接：对每一个选中的像素点都产生一条视点路径后，从步骤2-2中选择一条光源路径与它们相连接，即连接视点路径和光源路径上的最后一个点，使之成为一条完整的从光源到视点的路径。依照如下公式求得该完整路径的贡献度： $C_{i,j}^{*}={\mathrm{\&alpha;}}_i^L{c}_{i,j}{\mathrm{\&alpha;}}_j^E.$ 其中c_i,j的定义如下：

c_0,j＝L_e(z_j-1→z_j-2),c_i,0＝W_e(y_i-2→y_i-1),

$c_{i,j}=f_i(y_{i-2}\&RightArrow;y_{i-1}\&RightArrow;z_{j-1})G(y_{i-1}\&LeftRightArrow;z_{j-1})f_i(z_{j-2}\&RightArrow;z_{j-1}\&RightArrow;y_{i-1}),i,j>0$

其中i是指光源路径的长度，j是指视点路径的长度，G是形状因子函数，定义为：其中cosθ_o表示视点路径和光源路径上最后一个点相连接得到的直线L与视点路径上最后一个点所在的三角形表面法向量的余弦值，cosθ′表示该直线L与光源路径上最后一个点所在的三角形表面法向量的余弦值。

V是可见性函数，如果点x和x′在场景中可见，那么反之 $V(x\&LeftRightArrow;x^{\&prime;})=0.$

步骤3-4，计算加权贡献度：在之前计算得到的路径贡献度如果直接作为最终的颜色贡献度的话会导致噪声比较大，用多重重要性采样的方法可以极大地缓解这个问题，具体来说就是同一条路径的采样方式有很多种，但是通过不同的采样方法得到的几率是不一样的，在具体实现中需要加大一些采样概率低的路径的权重。权贡献度的计算公式为：

其中i是指光源路径的长度，j是指视点路径的长度，路径权值而比值p_a/p_i可以利用以下公式求得：

$\frac{p_1}{p_0}=\frac{P_A\left(x_0\right)}{P_{{\mathrm{\&sigma;}}^{\&perp;}}(x_1\&RightArrow;x_0)G(x_1\&LeftRightArrow;x_0)},$

$\frac{p_{a+1}}{p_a}=\frac{P_{{\mathrm{\&sigma;}}^{\&perp;}}(x_{a-1}\&RightArrow;x_a)G(x_{a-1}\&LeftRightArrow;x_a)}{P_{{\mathrm{\&sigma;}}^{\&perp;}}(x_{a+1}\&RightArrow;x_a)G(x_{a+1}\&LeftRightArrow;x_a)},(0<a<i)$

$\frac{p_{i+1}}{p_i}=\frac{P_{{\mathrm{\&sigma;}}^{\&perp;}}(x_{i-1}\&RightArrow;x_i)G(x_{i-1}\&LeftRightArrow;x_i)}{P_A\left(x_i\right)},$

步骤3-5，计算整个图像平面的像素颜色值：将步骤3-4中求出的路径加权贡献度作为当前像素颜色值存放起来，然后重复进行步骤3-1到3-4，直到整个图像平面上的像素的颜色值都被计算过一次为止，注意步骤3-3中选取光源路径的时候，之前计算中已经选择过的光源路径不能被重复选取。正如步骤3-2-2中所说，我们可以事先对每一个block生成这样一组随机的序列，这样可以避免每一轮选择像素的时候都要随机一次。最终会生成一张中间渲染图像。

步骤4，产生最终结果图像详细步骤如下：

步骤4-1，重复步骤2和步骤3，每次计算一整张新的平面颜色值，并与之前的结果相融合。假设当前是第T轮计算，而前面T-1轮中得到的某像素点Pi的颜色值为f_T-1(Pi)＝b，第T轮得到的加权路径贡献度为Co，那么第T轮结束后，该点的颜色值被更新为：当T＝1的时候，f_T(Pi)＝Co；

步骤4-2，评价渲染结果：用渲染轮数，或者时间来进行评价，当渲染总轮数达到某个阈值N(根据场景的复杂度来定，通常设定范围为:渲染轮数20-5000,渲染时间10分钟-10小时)的时候，我们认为渲染结果达标。结束渲染，得到最终结果图像。图5-图7展示了同一个场景在渲染轮数不同的情况下的结果，其中图5是渲染1轮后的结果，图6是渲染8轮后的结果，图7是渲染了32轮的结果。可以看出渲染轮数越多，得到的结果就越好。

实施例

本实施例的实验硬件环境是：Intel DualE2200 2.2GHz，2G内存，MicrosoftWindows7操作系统旗舰版，编程环境是visual studio 2012,directx9.0，测试用的场景模型文件来源于网上免费模型库，主要来自http://graphics.cs.williams.edu/data/meshes.xml。

实验场景包括sponza，cornell box等场景以及其变形。每一轮的采样率设置为2，光源路径与视点路径的最大长度均设置为6。划分图像平面时候的block设置为4*4的大小，这样保证了相当高的光源路径复用率，实际上，我们专门测试过block设置得更小的情况，如果将block设置到1*1的话，也就是说每一轮只用产生一条light path，用于整张图像中像素颜色值的计算，这样会产生一些难以消除的错误。

与标准bpt(双向路径跟踪)的比较显示，本发明的方法在相同时间内能得到更好的结果，图8显示的是两种方法在sponza场景上比较直观的比较(渲染时间均为1小时)。图9显示的是标准bpt(双向路径跟踪)方法在sponza场景上的与reference(参考结果)的差异，图10显示的是本发明的方法在sponza场景上的与reference(参考结果)的差异，可见在相同时间内本发明方法得到的结果略好。在这个场景下，本发明方法的MSE(平均平方误差)值为0.366021，标准bpt(双向路径跟踪)方法的MSE值为0.474184，可见本发明方法在计算收敛速度上要优于标准的bpt(双向路径跟踪)方法。

本发明在标准bpt(双向路径跟踪)方法的基础上使用了重用光源路径的思想提升渲染效率，并且通过分批采样图像空间产生视点路径的方法来缓解由于相关性采样所带来的噪声。总之本发明具有效率高，算法无偏，适用面广的特点。

本发明提供了一种计算机三维模型中全局光照的渲染方法，具体实现该技术方案的方法和途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。本实施例中未明确的各组成部分均可用现有技术加以实现。

Claims (5)

1.一种计算机三维模型中全局光照的渲染方法，其特征在于，包括以下步骤：

步骤1，加载模型：读入三维模型场景的点坐标，三角形序列，点的法向量，点的纹理；用kd-tree对三维模型场景空间进行划分，并将三角面存放于kd-tree中，同时从与该三维模型场景对应的材质文件中读入材质属性；

步骤2，采样光源：进行采样，得到若干条光源采样路径，使用分层采样的方法使得这些路径在路径空间中分布均匀；

步骤3，采样视点并计算像素颜色值：对最终渲染结果的图像平面进行划分，每次根据该划分选出一批像素，以视点到这些像素的方向作为初始方向采样路径，得到一组视点路径；从步骤2的结果中选择一条光源路径，与所述视点路径相连接，得到完整的从光源到视点的路径，计算整个路径对像素的颜色贡献值；重复直到对图像平面的采样点覆盖整个图像平面，得到一次迭代后的渲染结果；

步骤4，使用渐进式的渲染方式，重复步骤2和步骤3直到渲染结果的误差小于一个预先设定的值为止，生成最终的渲染结果。

2.如权利要求1所述的一种计算机三维模型中全局光照的渲染方法，其特征在于，步骤1包括以下步骤：

步骤1-1，加载模型文件：读入三维模型场景文件，获取点坐标v，点法向量vn，点纹理vt，三角形f，以及材质对应项；其中点坐标v和点法向量vn是一个三维向量，代表点在空间的位置，点纹理vt是一个二维向量，代表该点的纹理在纹理贴图中的坐标，三角形f是一个三维向量，包含了一个三角面中三个点的信息；

步骤1-2，加载与步骤1-1中的三维模型场景对应的材质文件，读取每个面所用材质属性；

步骤1-3，计算包围盒：首先计算出整个三维模型空间的最大范围，用一个轴对齐的长方体将整个模型包围起来，这个长方体的计算公式如下：

X_max＝max{v0_x,v1_x,…vn_x},X_min＝min{v0_x,v1_x,…vn_x},

Y_max＝max{v0_y,v1_y,…vn_y},Y_min＝min{v0_y,v1_y,…vn_y},

Z_max＝max{v0_z,v1_z,…vn_z},Z_min＝min{v0_z,v1_z,…vn_z}.

其中长方体在x，y，z轴方向上的坐标为X,Y,Z，而X_max和X_min分别代表长方体在x轴方向上的最大值和最小值，Y_max和Y_min分别代表长方体在y轴方向上的最大值和最小值，Z_max和Z_min分别代表长方体在z轴方向上的最大值和最小值；vi_x，vi_y，vi_z分别表示点集v中第i的点的x，y，z轴坐标，其中i＝0,1,……,n，n是点集v包含的点的数量；

步骤1-4，建立kd-tree并存储数据：将整个长方体空间利用kd-tree划分，并将之前读取三维模型场景文件得到的三角形保存于kd-tree中，kd-tree在k维空间中用轴对齐的子平面一次次地划分整个空间，最后得到一个二叉树的数据结构；

Kd-tree的划分基于SAH策略，假设对于空间S选取剖分面P，把该空间划分为左右两个子空间S_L和S_R，并且同时把三维模型场景中的三角形面片划分成两个集合N_L和N_R，上述的SAH优化函数公式为：

${\mathrm{SAH}}_{\cos t}\left(\mathrm{node}\right)=C_{\mathrm{walk}}+\frac{\mathrm{SA}\left(S_L\right)}{\mathrm{SA}\left(S\right)}{N}_L{C}_{\mathrm{hit}}+\frac{\mathrm{SA}\left(S_R\right)}{\mathrm{SA}\left(S\right)}{N}_R{C}_{\mathrm{hit}},$

其中node是指当前空间节点，SA(S)表示空间S的表面积，C_hit表示光线与每个三角形的求交代价，C_walk表示光线遍历当前node的代价，在划分每个节点的时候采用使得SAH_cost(node)最小的划分方式，且在满足以下任一条件的时候，终止划分，包括：

1)当前节点的深度超出一个预先设定的阈值Deep，1≤Deep<∞,且为整数；

2)当前节点中的三角形个数少于一个预先设定的阈值TN，0≤TN<20；

最终得到一个包含了所有三维模型三角形的kd-tree数据结构，所述kd-tree只有叶子节点有三角形数据，非叶子节点的三角形数据为空。

3.如权利要求2所述的一种计算机三维模型中全局光照的渲染方法，其特征在于，步骤2包含下列步骤：

步骤2-1，从光源上采样路径的起始点：从光源上随机采样一个点P₀作为光源路径的起始点，如果是点光源，那么路径起始点就是光源位置，如果是面光源或者是光源体，那么路径的起始点是这个面光源或者光源体上的随机采样点；

步骤2-2，从路径起点P₀随机产生路径：依次采样路径的下一个点，直到光线射出整个三维模型场景，或者被三维模型场景表面吸收，或者是路径长度到达阈值d；产生路径的方式为：光线在场景空间传播，遇镜面发生镜面反射，遇漫反射面随机在表面法线方向上半球面取下一反射方向，遇透明表面发生折射，这时将产生出一条采样得到的光源路径，计算这条光源路径的贡献度，计算公式为：

${\mathrm{\&alpha;}}_0^L=1,$

${\mathrm{\&alpha;}}_1^L=\frac{L_e^{\left(0\right)}\left(y_0\right)}{P_A\left(y_0\right)},$

${\mathrm{\&alpha;}}_r^L=\frac{f_s(y_{r-3}\&RightArrow;y_{r-2}\&RightArrow;y_{r-1})}{P_{{\mathrm{\&sigma;}}^{\&perp;}}(y_{r-2}\&RightArrow;y_{r-1})}{\mathrm{\&alpha;}}_{r-1}^L,r\&GreaterEqual;2,,$

其中表示光源路径上光源位置的颜色贡献度，表示光源路径上第1个点位置的颜色贡献度，表示光源路径上第r个点位置的颜色贡献度，P_A(y₀)表示光源上采样到y₀点的概率，表示光源在y₀的亮度值，f_s(y_r-3→y_r-2→y_r-1)则表示在已得到点y_i-2的情况下，采样到点y_r-1的概率；

$P_{{\mathrm{\&sigma;}}^{\&perp;}}(y_{r-2}\&RightArrow;y_{r-1})=\frac{P_{\mathrm{\&sigma;}}(y_{r-2}\&RightArrow;y_{r-1})}{\cos \mathrm{\&theta;}},$

其中P_σ(y_r-2→y_r-1)表示采样y_r-2→y_r-1方向的概率，而θ表示采样方向与y_r-2点的表面法向量的夹角。

4.如权利要求3所述的一种计算机三维模型中全局光照的渲染方法，其特征在于，步骤3包含下列步骤：

步骤3-1，划分图像平面，假设图像平面一共有N个像素，每次不重复地取出M个像素点，M<N，用以产生视点路径，总共分轮取完，每次取出的M个像素点和之前已经取出的不重复；步骤3-2，从视点向选中的像素点发射光线，作为初始路径，依次产生出视点路径，视点路径的贡献度计算公式为：

${\mathrm{\&alpha;}}_0^E=1,$

${\mathrm{\&alpha;}}_1^E=\frac{W_e^{\left(0\right)}\left(z_0\right)}{P_A\left(z_0\right)},$

${\mathrm{\&alpha;}}_{r^{\&prime;}}^E=\frac{f_s(z_{r^{\&prime;}-3}\&RightArrow;z_{r^{\&prime;}-2}\&RightArrow;z_{r^{\&prime;}-1})}{P_{{\mathrm{\&sigma;}}^{\&perp;}}(z_{r^{\&prime;}-2}\&RightArrow;z_{r^{\&prime;}-1})}{\mathrm{\&alpha;}}_{r^{\&prime;}-1}^E,r^{\&prime;}\&GreaterEqual;2,$

其中表示视点路径上视点位置的颜色贡献度，表示视点路径上第1个点位置的颜色贡献度，表示视点路径上第r′个点位置的颜色贡献度，P_A(z₀)表示在视点棱镜表面上采样到z₀点的概率，表示z₀点在视点棱镜表面的权值；

步骤3-3，路径连接：对每一个选中的像素点都产生一条视点路径后，从步骤2-2中选择一条光源路径与它们连接，即连接视点路径和光源路径上的最后一个点，从而形成一条完整的从光源到视点的路径；

依照如下公式求得该完整路径的贡献度：

$C_{i,j}^{*}={\mathrm{\&alpha;}}_i^L{c}_{i,j}{\mathrm{\&alpha;}}_j^E,$

其中i是指光源路径的长度，j是指视点路径的长度，与分别由步骤2-2以及步骤3-2得到，c_i,j的计算公式如下：

c_0,j＝L_e(z_j-1→z_j-2),c_i,0＝W_e(y_i-2→y_i-1),

$c_{i,j}=f_i(y_{i-2}\&RightArrow;y_{i-1}\&RightArrow;z_{j-1})G(y_{i-1}\&LeftRightArrow;z_{j-1})f_i(z_{j-2}\&RightArrow;z_{j-1}\&RightArrow;y_{i-1}),i,j>0,$

其中i是指光源路径的长度，j是指视点路径的长度，G是形状因子函数，定义为：其中cosθ_o表示视点路径和光源路径上最后一个点相连接得到的直线L与视点路径上最后一个点所在的三角形表面法向量的余弦值，cosθ′表示该直线L与光源路径上最后一个点所在的三角形表面法向量的余弦值；

V是可见性函数，如果点x和x′在场景中可见，那么反之 $V(x\&LeftRightArrow;x^{\&prime;})=0;$

步骤3-4，计算加权贡献度：加权贡献度的计算公式为：

$C_{i,j}=w_{i,j}{C}_{i,j}^{*},$

其中i是指光源路径的长度，j是指视点路径的长度，路径权值比值p_a/p_i根据以下公式求得：

$\frac{p_1}{p_0}=\frac{P_A\left(x_0\right)}{P_{{\mathrm{\&sigma;}}^{\&perp;}}(x_1\&RightArrow;x_0)G(x_1\&LeftRightArrow;x_0)},$

$\frac{p_{a+1}}{p_a}=\frac{P_{{\mathrm{\&sigma;}}^{\&perp;}}(x_{a-1}\&RightArrow;x_a)G(x_{a-1}\&LeftRightArrow;x_a)}{P_{{\mathrm{\&sigma;}}^{\&perp;}}(x_{a+1}\&RightArrow;x_a)G(x_{a+1}\&LeftRightArrow;x_a)},(0<a<i)$

$\frac{p_{i+1}}{p_i}=\frac{P_{{\mathrm{\&sigma;}}^{\&perp;}}(x_{i-1}\&RightArrow;x_i)G(x_{i-1}\&LeftRightArrow;x_i)}{P_A\left(x_i\right)},$

步骤3-5，计算整个图像平面的像素颜色值：将步骤3-4中求出的路径加权贡献度作为当前像素颜色值存放起来，然后重复进行步骤3-1到3-4，直到整个图像平面上的像素的颜色值都被计算过一次为止，步骤3-3中选取光源路径的时候，之前计算中已经选择过的光源路径不重复选取，最终生成了一张中间渲染图像。

5.如权利要求3所述的一种计算机三维模型中全局光照的渲染方法，其特征在于，步骤4包含下列步骤：

步骤4-1，重复步骤2和步骤3，每次计算一整张新的平面颜色值，并与之前的结果相融合；假设当前是第T轮计算，而前面T-1轮中得到的某像素点Pi的颜色值为f_T-1(Pi)＝b，第T轮得到的加权路径贡献度为Co，那么第T轮结束后，该点的颜色值被更新为：当T＝1的时候，f_T(Pi)＝Co；

步骤4-2，评价渲染结果：当渲染总轮数或时间达到阈值的时候，渲染结束。