CN105118085A - Kd-树与Voronoi图混合的辐射亮度计算方法 - Google Patents

Abstract

本发明涉及一种Kd-树与Voronoi图混合的辐射亮度计算方法，其步骤包括：1)构造场景拓扑结构及其Kd-树层次组织；2)以多边形面片为单位存储生成的光子图；3)遍历三维场景表面的每个多边形，收集属于每个多边形的光子并组织成二维Kd-树结构；4)对每个多边形表面构造Voronoi图划分，并计算Voronoi单元的辐射通量密度；5)计算从视点发射的射线与场景的相交点x；6)根据x所属的多边形面片，利用所保存的二维Kd-树结构收集临近的光子；7)根据Voronoi图划分找到收集到的光子所覆盖的精确面积，并计算x点的辐射亮度。本发明能够利用局部表面光子分布的几何特征进行更加精确的辐射亮度内核估计。

Description

Kd-树与Voronoi图混合的辐射亮度计算方法

技术领域

本发明属于计算机图形学技术领域，具体涉及一种Kd-树(Kd-tree)与Voronoi图混合的精确辐射亮度计算方法。

背景技术

具有高度真实感的绘制技术是计算机图形学中最震撼人心的研究领域之一，而其中精确细致的全局光照明模拟尤其重要。光子映射算法(PhotonMapping，简称PM)是一种基于光子图实现的全局光照明算法(JensenH.W.:Globalilluminationusingphotonmaps.InRenderingTechniques96.Springer,1996,pp.21–30)，由于光子映射相比光线追踪能够更好的模拟焦散现象和SDS，即经过镜面反射(折射)-漫反射-镜面反射(折射)所产生的现象，而且光线追踪算法在采样较低的条件下，使用蒙特卡洛采样精确计算初始的反射时会生成噪声，因此光子映射算法可以得到更加逼真自然的图像。

光子映射算法中光子图的生成也需要借助光线追踪的方法，光线追踪算法对镜面反射光线和透射光线的追踪要优于光子映射算法，在实际场景中往往选取二者的混合算法解决绘制问题，即使用光线追踪方法计算光源引起的直接光照，而使用光子图计算环境引起的间接漫反射光照。

光子映射算法中核心的辐射亮度估计的精度常常严重依赖于半径非常小的估算内核以及在估算内核上的光子分布。对于几何特征复杂的绘制场景，光子映射不能保证在绘制中可以有效地控制偏差和噪声错误的生成。光子映射算法是一个有偏算法，偏差产生的原因在于辐射亮度估算本质上是计算x点处光子密度的过程，由于x点是抽象点不存在体积和大小，只能近似的使用x点附近区域内的平均光子密度来代替x点的光子密度，这样做必然会引入不精确的偏差。间接光照计算在光子数目趋近于无穷大是可以趋近于正确的结果，但是实际上，由于光子映射的整个绘制方法中发射的光子数目不可能达到无穷大(否则计算机系统无法进行处理)，分布在辐射亮度估计区域内的光子通常是有限甚至是稀疏的，这样就造成了估算内核的半径在有限光子的条件下不能趋近于0，直观上，与相交点距离“较远”的光子也会用于计算相交点的辐射亮度，这一部分辐射量是不精确的。此外，传统的光子映射算法的辐射亮度计算存在两个假设条件：1.估算内核圆盘所在的局部物体表面是平坦的几何表面，即在相交点x附近的局部表面可以完整地包含估算内核圆盘；2.所有与相交点x的距离小于r的光子都会对x点产生有效的辐射量贡献(ToshiyaT.,OgakiS.,JensenH.W.:Progressivephotonmapping.ACMTrans.Graph.27,5(Dec.2008),130:1–130:8)。同时，在光子映射方法中，光子图往往存储在三维空间Kd-tree(Kd-树)层次组织结构中离散光子，已经完全丢弃了待绘制三维场景的几何信息，而在进行光子收集阶段由于对点x周围的空间几何信息完全未知，因此估算计算内核对整个圆盘范围内均有光子分布的假设在实际中也不能满足，即辐射亮度估算方程的分子和分母都是不精确的，这样导致的错误辐射亮度估算结果为偏差错误。偏差问题所导致的场景绘制结果的表现为模糊，即丢失了局部细节。

发明内容

针对上述问题，本发明提出一种Kd-tree与Voronoi图混合的精确辐射亮度计算方法(RadianceEstimate)，能够利用局部表面光子分布的几何特征进行更加精确的辐射亮度(Radiance)内核估计。

本发明采用的技术方案如下：

一种Kd-tree与Voronoi图混合的精确辐射亮度计算方法，其步骤包括：

1)针对待绘制的三维场景，构造场景拓扑结构及其Kd-tree层次组织；

2)光源发射光子并生成光子图，以多边形面片为单位存储该光子图，包括存储光子的基本属性以及存储光子所停留的多边形面片的唯一标识，同时每个多边形都记录落在其中的光子的唯一标识；

3)遍历三维场景表面的每个多边形，针对每个多边形，收集属于该多边形的光子并组织成二维Kd-tree结构；

4)对每个多边形表面构造Voronoi图划分，并计算划分出的Voronoi单元的辐射通量密度；

5)计算从视点发射的射线与场景的相交点x；

6)根据相交点x所属的多边形面片，利用该多边形面片所保存的二维Kd-tree结构收集临近的光子；

7)对收集到的光子，根据Voronoi图划分找到这些光子所覆盖的精确的面积，并利用辐射亮度计算公式计算出x点的辐射亮度。

进一步地，根据步骤7)计算出的辐射亮度合成出最终的真实感图像。

进一步地，步骤1)构造场景拓扑结构的方法是：首先为每个场景中的节点建立邻域网格结构，通过遍历一次场景中的网格结构，得到从网格的节点到邻域网格集合的对应关系表；之后为每个网格建立1-邻域网格集合，使用的方法是针对一个网格，遍历它的所有节点，选取每个节点邻域网格集合的并集作为该网格的1-邻域网格集合；待所有网格的1-邻域网格集合创建完成之后，递归地创建网格的2-邻域网格集合，依次类推。

进一步地，步骤4)构造Voronoi剖分时，首先利用Voronoi图计算每个采样光子所在Voronoi单元格的面积大小，然后利用快速邻近查找到以相交点为圆心以内核半径为半径的圆盘内的所有光子，并获取其有效的辐射度以及Voronoi单元格的面积。

进一步地，步骤6)使用栈结构收集临近的光子。

进一步地，步骤7)采用如下公式计算辐射亮度：

$L_s(x,\stackrel{\→}{\mathrm{\ω}})\≈\frac{\underset{p=1}{\overset{n}{\Σ}}{f}_r(x,\stackrel{\→}{\mathrm{\ω}},{\stackrel{\→}{\mathrm{\ω}}}_p){\mathrm{\Δ\φ}}_p(x,{\stackrel{\→}{\mathrm{\ω}}}_p)}{\underset{p=1}{\overset{n}{\Σ}}{\mathrm{\ΔA}}_p},$

其中，x是视点光线与物体表面相交点；n是相交点x附近用于辐射亮度估计的光子数目；是出射方向；是光子入射方向；是双向反射分布函数；是光子p所携带的有效辐射通量；△A_p是包含光子p的Voronoi多边形的面积。

本发明提出了一种新的具备几何感知的精确辐射亮度计算方法，能够利用局部表面光子分布的几何特征进行更加精确的辐射亮度估算；为了实现更精确的内核估算，本发明提出每个光子占据的局部表面的Voronoi图剖分单元的辐射度概念，利用该辐射度进行精确的辐射亮度估计计算。与标准的光子映射算法不同，本发明能够绘制中消除由于非精确估算内核导致的可见偏差。

附图说明

图1是本发明的辐射亮度计算方法的步骤流程图。

图2是全局光子与焦散光子示意图。

图3是Kd-tree建树过程示意图。

图4是光子及其占据的Voronoi剖分单元的示意图。

图5是相交点x的估算内核范围内的Voronoi多边形示意图。

图6是传统光子映射方法和本发明方法的绘制效果对比图。

图7是本发明绘制的场景的真实感效果图。

具体实施方式

为使本发明的上述目的、特征和优点能够更加明显易懂，下面通过具体实施例和附图，对本发明做进一步说明。

本发明提出了一种新的具备几何感知的精确辐射亮度计算方法，采用基于Kd-tree与Voronoi图混合结构的算法，能够利用局部表面光子分布的几何特征进行更加精确的辐射亮度内核估计。本发明的整个新的光子映射方法(辐射亮度计算方法)的流程如图1所示：

1)针对待绘制的三维场景，构造场景拓扑结构及其Kd-tree层次组织；

2)光子跟踪过程开始：光源发射光子生成光子图，光子图存储以多边形面片为单位，即每个光子一旦停止跟踪停留在某个场景的表面时，除了存储其光子的基本属性外，还存储该光子所停留(即所属)多边形面片的唯一标识(ID)，同时每个多边形都记录落在其中的光子的唯一标识(ID)；

3)遍历三维场景表面的每个多边形，针对每个多边形，收集属于该多边形的光子并组织成二维Kd-tree结构；

4)对每个多边形表面构造Voronoi图划分，并计算划分出的Voronoi单元的辐射通量密度；

5)光线跟踪过程开始：计算从视点发射的射线与场景的相交点x；

6)根据相交点x所属的多边形面片，利用该多边形面片所保存的二维Kd-tree结构收集临近的光子；

7)对收集到的光子，根据Voronoi图划分找到这些光子所覆盖的精确的面积，并利用辐射亮度计算公式计算出x点的辐射亮度；

8)合成出最终的真实感图像。

上述过程中，构造场景拓扑结构和Kd-tree层次组织(步骤1)→光子跟踪→光线跟踪→合成最终真实感图像(步骤8)，这四步基本步骤与现有的光子映射方法流程一致，不具备创新性，但是光子跟踪过程中的2、3、4三个步骤的方法与以往的光子跟踪方法不同，光线跟踪过程中的6、7步骤与以往的光线跟踪方法以及辐射亮度计算方法不同，具有创新性。

下面具体描述本发明方法的各个步骤。表1说明了本发明所使用的符号的含义。

表1.主要符号对照表

对于绘制算法，辐射亮度是一类非常重要的计算，涉及到多种辐射量的概念，以下给出定义：

辐射能Q是一种电磁波能量，与光照相关的辐射能量来源于光源。

辐射通量φ(又被称为辐射量)定义为辐射能Q每秒通过物体表面的数量，即φ＝dQ/dt。

辐射度E(即物体表面的辐射通量密度)定义为单位表面积内的辐射通量，即E＝dφ/dA。

辐射亮度L(光源的辐射亮度又被称为辐射强度)定义为每单位立体角、每单位面积上的辐射通量，即L＝dE/dt。

由于材质具有线性光学特性，入射辐射量正比于出射辐射量，定义两者之间的关系为双向反射分布函数(BidirectionalReflectanceDistributionFunction，简称BRDF)，BRDF描述了光线与物体表面之间的反射机制，其形式为f_r(x,ω_i,ω_o)。其中x点为入射光线与反射光线的交点，ω_o为出射方向，ω_i为入射方向。经典的Phong光照模型和Lambert光照模型都是对BRDF的一种近似表达。双向透射率分布函数(BidirectionalTransmittanceDistributionFunction，简称BTDF)将BRDF扩展到透射表面的光照计算，BRDF和BTDF统称为双向散射分布函数(BidirectionalScatteringDistributionFunction，简称BSDF)，为了表达方便，本发明中使用的辐射亮度计算公式仅考虑BRDF的场景，BTDF场景可以类比得出相应的辐射亮度计算公式。

步骤一：

1-1.模型场景组织

场景使用自定义的树型结构组织。场景是由多个模型组成的，一个根节点保存各个模型的节点信息。本发明可以处理各种网格表示的模型，包括任意多边形网格、三角形网格、四边形网格等，为便于描述，后续统一用三角形网格或者三角形面网格作为本发明具体实现方法进行论述。在每个模型的节点数据中，指定的模型顶点数据被保存在连续的顶点缓冲区中，另一部分三角面数据被放入与网格相关的缓冲区中，利用以上数据在程序进行绘制之前的预处理工作中将每个模型中的数据组织为空间加速结构(如Kd-tree结构)，该Kd-tree结构为了加速计算也可以采用在GPU上建立的方式。本发明按照层次场景结构来组织三角面数据，将这些场景数据打包发给光线跟踪和光子跟踪模块，与光线追踪相关的过程如产生光线、相交测试、相交处理等分别封装为GPU端内核函数，使得单个光线的追踪过程成为上述多个内核函数的组合序列。CPU端需要分别为光线追踪过程和光子追踪过程中产生的多种光线指定其光线生成、相交测试以及相交处理等多个内核函数入口，使得每种光线在被发射执行前其将要进行的追踪过程都是明确且有效的。

1-2.构建场景拓扑结构

本方法一个重要预处理工作就是建立场景网格的拓扑结构，它与原始模型的网格输入不同，我们需要使用邻接表的数据结构实现快速查找。首先我们的算法为每个场景中的节点建立邻域三角形结构，仅需要遍历一次场景中的网格结构，便可得到从网格的节点到邻域网格集合的对应关系表。之后我们为每个网格建立1-邻域网格集合，使用的方法是针对一个网格，遍历它的所有节点，选取每个节点邻域网格集合的并集作为该网格的1-邻域网格集合。待所有网格的1-邻域网格集合创建完成之后，可以递归的创建网格的2-邻域网格集合，依次类推，尽管这部分计算时间可能随着邻域增加而呈现指数型增长，但是作为预处理工作的时间是可以接受的。

1-3.三角面场景优化

针对三角面网格组成的场景输入，假设场景物体为较为光滑的表面，即不存在高频噪声，当相交点在高频噪声上时，其所在的切平面具有极大的随机性。实际场景中高频噪声对于间接光照不会产生显著的影响，对于包含高频噪声的场景输入，我们的算法采用表面滤波的方法除去噪声部分，即假设场景输入为不含高频噪声的表面光滑三角面网格。在非平面的物体表面完成精确估算内核的计算需要将非平面表面投影到相交点所在的切平面。当相交点x位于三角面内部时，x点的法相使用顶点法相插值得到。此时，x所在三角面投影到切平面的面积计算可以使用三角面所表示的曲面投影得到，这种方法在已知物体表面为均匀曲率的条件下可以得到快速准确的估算面积。

步骤二：

2-1.光子图概念

将环境引起的间接漫反射光照分为两种，并使用两个光子图分别计算这两部分间接光照计算中。从光源发射的光子，经过多次镜面反射或者规则透射最终落在一个具备漫反射材质的物体表面，这部分光子被存储在焦散光子图中；同样从光源发射的光子，但是经过至少一次漫反射而最终落在一个具备漫反射材质的物体表面，这部分光子被存储在全局光子图中。全局光子图和焦散光子图分别用于计算环境引起的漫反射中包含漫反射路径和不包含漫反射路径的光照分量。将两种光照来源分开的原因在于，焦散光子图用于在漫反射或者具有轻微光泽表面生成焦散，全局光子图用于生成柔和的环境间接光照。焦散现象是高密度的光子聚集形成的，算法从光源向镜面反射或者规则透射表面发射光子建立焦散光子图，而全局光子图作为场景中均匀辐射光照的粗略近似，算法从光源向场景中所有的物体发出光子建立全局光子图。由于两种光子图的特性不同，可以将两种光子图建立为不同的加速结构，两种光子图的辐射亮度计算分开进行可以提高算法的精确度和整体计算性能。

图2中展示了全局光子和焦散光子的区别。图中光子从光源L出发经过透明材质的圆球1最后停留在具有漫反射属性的地面2上形成焦散光子图；光子从光源出发经过墙面3或者地面2最后停留在具有漫反射属性的地面或者墙面上形成全局光子图。

2-2.建立光子图

本方法中的光子追踪过程主要完成生成光子图工作。光子追踪过程的光线从光源出发，按照光源半球空间采样，每一个光线对应于半球空间上一个立体角。光子被设定为该光线与场景相交的具备漫反射属性的表面，一旦光线与表面相交，那么相交点信息将被记录下来称为光子，如果该物体表面同时具备反射和透射材质，则使用俄罗斯轮盘赌算法计算并确定究竟在该表面发生反射、投射还是吸收，从而会产生新的反射或者投射光线用于下一个光子追踪过程。

2-3.俄罗斯轮盘赌

俄罗斯轮盘赌是借助赌博游戏实现重要性采样的方法，在光子映射算法中应用俄罗斯轮盘赌的场景主要是在创建光子图阶段。由于理想的镜面反射和漫反射物体现实中是不存在的，实际的物体表面往往是多种材质的复合体。基于俄罗斯轮盘赌的光子追踪中实现类似如下过程，针对一个较为复杂的场景，通常假设物体同时具有透明反射和漫反射材质，两者的重要性比例分别为β_t和β_d，同时还有一定的概率光线会被物体吸收β_a，β_t+β_d+β_a＝1。产生介于0-1之间的随机数λ₁，首先检测λ₁<β_a是否成立，如果成立，表明光子的能量被完全吸收，追踪结束。接着检测λ₁<β_a+β_d是否成立，如果成立，表明光子落在具备漫反射属性的表面，当前光子的追踪过程结束，此时的光子将被记录在光子图中留作后须辐射亮度计算使用。如果检测不成立，那么表明光线入射到具有透明材质的表面，此时，依照入射场景将使用Fresnel双向反射中的一种反射与投射分布曲线，仍然需要产生领域一个介于0-1之间的随机数λ₂，如果λ₂<k_r则该光子是反射，否则该光子折射，经过上述过程，能完全确定当前光子的传输形态。

所有被记录下来的光子最终就形成了光子图，光子图存储以多边形面片为单位，即每个光子一旦停止跟踪停留在某个场景的漫反射表面时，除了存储其光子的基本属性外，还存储该光子所停留(即所属)多边形面片的唯一标识(ID)，同时每个多边形都记录落在其中的光子的唯一标识(ID)。

步骤三：光子遍历构建Kd-tree

遍历三维场景中的每个三角形。对每个三角形，我们均采用二维Kd-tree组织落在该三角形表面的光子，然后在其上实现方便的二维KNN(k-NearestNeighbor)搜索，由于从三维退化到二维相比较原始光子映射算法的三维Kd-tree，KNN搜索的效率得到了进一步的提高，一旦给出相交点x的位置，就可以在x所在表面网格的二维光子Kd-tree中快速的定位到邻近空间内与其相距r以内的光子，使用这些光子进行辐射亮度估计，同时利用到光子所携带的Voronoi多边形面积计算精确估算内核。

Kd-tree(k-dimensional树的简称，或者Kd树，或者k-d树，或者kd-tree)是一种分割k维数据空间的数据结构，是目前加速空间遍历过程的有效手段，Kd-tree沿着坐标轴方向递归的将场景中所有采样光子划分为子空间，每一个子空间存储包含在其内的采样光子。Kd-tree算法分为建树过程和搜索树过程。

构建Kd-tree空间结构需要分层对场景进行划分，具体的划分方法是选择与坐标轴平面平行的分割面将当前光子空间分为左右两部分，将空间中的光子对象对应插入到左右两子空间中，两个子节点的子空间没有交集。从整个场景空间开始，用上述方法递归的划分子空间，直到当前空间中光子对象数目少于阈值为止。Kd-tree空间分割面的选取有多种方法，重心法，SAH评估法等(ZhouK,HouQ,WangR,etal.Real-timeKD-treeconstructionongraphicshardware.ACMTransactionsonGraphics(TOG),2008,27(5):126.)。本实施例使用简单的中位数方法选取采样光子作为分割点。如图3以2D情况下的Kd-tree为例，其中(b)图展示了(a)图中场景的Kd-tree建树结果，树的内部节点既包含采样光子，同时指示空间中的分割平面，树的叶节点中仅包含采样光子。

步骤四：

4-1.Voronoi图划分

在后续的每个光子追踪过程中需要遍历KNN邻近搜索到的光子，并获得这些邻近光子的准确的Voronoi多边形面积。因此，需要对这些光子所占据的三角形表面进行Voronoi剖分。我们首先利用Voronoi图计算每个采样光子所在Voronoi单元格的面积大小，然后利用快速邻近查找到以相交点为圆心以内核半径为半径的圆盘内的所有光子，并获取其有效的辐射度以及Voronoi单元格的面积。

Voronoi图是计算几何中重要的几何结构，Voronoi图的求解点集或其他几何对象与距离相关的问题时起着重要作用。Voronoi图又叫泰森多边形或Dirichlet图，它是由一组由连接两邻点直线的垂直平分线组成的连续多边形组成。Voronoi图的基本概念如下：

设p₁，p₂是平面上的两点，L是的垂直平分线，L将平面分为L_L和L_R两部分，使用H(p₁,p₂)表示半平面L_L，H(p₂,p₁)表示半平面L_R。那么给定平面上n个点的集合S＝{p₁,p₂,...,p_n}。定义V(p_i)＝∩_i≠j(p_i,p_j)，由定义可知，V(p_i)是一个不多于n-1条边的多边形区域，称之为Voronoi多边形。对于S中的每个点都可以做一个Voronoi多边形，这样n个Voronoi多边形所组成的图即为Voronoi图，即为Vor(S)(周培德：计算几何—算法设计与分析(第4版).清华大学出版社.2011年)。

按照Voronoi图的定义，对于位于一个Voronoi多边形内的x，S集合中距离x最近的点一定是位于该Voronoi多边形内的那个点。实质上Voronoi图是对二维平面的一个完整划分，将平面的每部分区域划分到每一个Voronoi多边形内，即划分为集合S的每一个点区域。可以得出，Voronoi图是一种包含集合S点分布以及表面区域信息的重要结构，在几何中通常会使用Voronoi图来表达和分析离散点的几何分布特征。

4-2.近似构造Voronoi图

为网格中的光子集合构造Voronoi图，最经典的方法是使用增量构造方法构造Voronoi图，但是其时间复杂度达到了Ο(n²)，其中n为输入的网格上光子数目，这一时间复杂度相对空间加速结构如Kd-tree的代价过高，另一种较快的构建算法是分支法，分支法通过与Kd-tree类似的二分策略构造Voronoi图可以是构造时间复杂度缩减到Ο(nlog(n))，但是这一时间复杂度仍然严重依赖于网格中光子的数目n，我们的算法在实现时考虑使用跳跃式泛洪算法(RongG,TanTS.JumpfloodinginGPUwithapplicationstoVoronoidiagramanddistancetransform.Proceedingsofthe2006symposiumonInteractive3Dgraphicsandgames.ACM,2006:109-116.)，泛洪算法是一类基于计算机的近似算法，与严格得出Voronoi图的拓扑结构不同的是，泛洪算法尝试在构造纹理缓冲区模拟网格表面的面积，使用像素的概念统计计算每个Voronoi多边形的面积大小。

使用泛洪方法求解Voronoi图，该算法的核心思想是以可变步长l从原始的光子位置出发不断的给周围“黑色”区域染色，在第一轮迭代中步长取为l＝e/2，之后每个迭代中将步长减半，直到最后将正方形表面所有区域“染色”，即得到正方形区域对于输入8个采样光子的Voronoi划分，统计标记不同颜色的点的数目，可以近似得到每个Voronoi多边形的面积大小。

关于洪泛算法的时间复杂度的额外说明，上述染色过程实际上包含两部分内容，即1.计算采样点与要染色像素点的距离；2.比较该距离与已保存最小距离，如果产生更小距离则更新颜色，否则维持像素点原来的颜色。从上述泛洪算法求解步骤可以很容易看出，采样光子中的颜色可以在log(e²)轮的洪泛算法中到达正方形纹理缓冲区域所有的像素格点。在每一轮的“染色”过程中，需要“染色”的像素数目以4为等比不断增长，在倒数第二轮时缓冲区域所有的像素格点均至少被“染色”一次。而且最后一轮进行l＝1的“染色”时最多会进行4e²个像素的染色过程。由于每个像素的染色过程可以近似为Ο(1)的计算复杂度，所以，洪泛算法的计算复杂度为Ο(e²)。

使用上述泛洪算法求解Voronoi图可以实现与采样光子数目无关的时间复杂度快速计算出Voronoi多边形的面积大小，Voronoi多边形的面积精确度取决于使用纹理缓冲区的分辨率大小。在实验中最大选取缓冲区大小为128*128，完全可以满足所有网格对面积计算精确度的需求。该算法的另一个优势在于可以方便的在GPU上并行实现，使得程序的性能得到进一步的提升，在保证精确度的基础上同时使得算法整体的计算时间可以得到有效控制。

如图4所示，图中黑点为光子，每个光子所占据的Voronoi剖分单元用不同颜色(图中显示为不同灰度)的区域，其中Voronoi多边形不同的颜色来表征不同的光子辐射量密度。

步骤五：

光线追踪过程主要完成生成相交点x和计算直接光照两个工作。光线追踪过程的光线从视点出发，按照屏幕空间采样，每一个光线对应于屏幕上一个像素。一个光线的有效相交点x被设定为该光线与场景相交的第一个具备漫反射属性的表面，一旦光线与表面相交，那么相交点信息将被记录下来用于之后收集光子，同时该物体表面的材质信息也会被记录以计算相交点的BRDF，如果该物体表面同时具备反射和投射材质，那么就会产生新的反射或者投射光线用于下一个光线追踪过程。

与记录相交点x不同的是计算直接光照过程在光线与场景中的每一次相交时都会被执行，直接光照的计算通过向光源发射采样阴影光线的方法完成，阴影光线如果在相交点与光源的路径上不存在其它物体那么采样光源信息成功，否则视为物体遮挡，采样成功的阴影光线越多，则光源对相交点产生的直接光照便越强。这里发射的阴影光线数目越多，最终得到的阴影和直接光子便更加柔和，我们的算法中通常为每个相交点发射3×3＝9条阴影光线，即可获得比较柔和自然的直接光照。

光线追踪算法中遇到最多的采样技术是针对半球空间和平面的采样。针对物体表面的面采样以正方形为例，最简单的采样为随机采样，方法非常简单，假设生成n个位于(0，0)和(1.0，1.0)正方形内的采样点，可以随机的产生2n个0-1之间的随机数，每两个配对称为一个采样点，但是这样的方法分布过于随机，相比较一维随机采样，直接将随机采样扩展到二维得到的采样分布并不均匀。一种简单的采样方式为抖动采样。抖动采样与随机采样不同的是加入了分层机制，首先第一层将原始正方形等分为n个Voronoi多边形，接着在这n个Voronoi多边形内分别利用随机采样选取一个采样点，这样得到抖动采样结果相比较随机采样拥有更好的采样点分布。上述介绍的采样方法主要用于光线跟踪算法中对于面光源的模拟，从相交点向面光源发射多条阴影检测光线，那么阴影检测光线的朝向如何选取既可以利用抖动采样在面光源上采样多个点之后确定多条检测光线，这样做会大大削弱由于采样光线不足导致的阴影锯齿效果。

针对半球空间的采样技术通常选取二维随机参数的半球采样方法，假设相交点为x，其所在平面的法向为在相交点x的法向空间建立三维坐标系其中输入的二维随机参数为(sample.x，sample.y)，其中sample.x和sample.y都位于0-1之间。选取sample.x作为单位方向向量在平面圆盘上的转角弧度值，而sample.y作为单位方向投影在平面圆盘上的半径大小，可以得到半球空间的采样方向为dir，好处在于可以获得在在半球空间内均匀的采样方向，即当方向的贡献越大时，相应的在和上可选的采样就越少，符合半球的空间特征。上述方法主要用于入射光与具备理想漫反射相交时产生的随机二次光线，具体公式如下：

phi＝sample.x·2π

$\begin{array}{c}dir=\sqrt{sample.y}\&CenterDot;\cos \left(phi\right)\&CenterDot;\stackrel{\&RightArrow;}{U}\\ +\sqrt{sample.x}\&CenterDot;\cos \left(phi\right)\&CenterDot;\stackrel{\&RightArrow;}{V}\\ +\max (0,1-sample.y)\&CenterDot;\stackrel{\&RightArrow;}{W}\end{array}$

步骤六：

利用构建的Kd-tree空间结构实现邻域查找需要使用到栈结构，起始栈结构为空，搜索过程将从相交点x所在的多边形中的包含的所有光子作为根节点开始检测，如果当前检测节点包含光子p_c，那么将根据该检测光子p_c所在子树划分平面l，判断相交点x落在了划分平面l的左边还是右边：假设相交点x落在了划分平面的左边，这时计算相交点x与划分平面l的距离d，比较d与内核半径r的大小，如果d<r说明相交点x非常接近划分平面l，那么划分平面的左右两边都需要继续递归搜索下去，即将检测光子p_c的两个子节点全部压栈；反之说明相交点x远离划分平面l，这时候仅需递归搜索检测光子p_c的左子节点。当相交点x落在了划分平面的右边时也可以利用同样的策略进行递归搜索。如果当前检测光子p_c不包含子节点，那么就从栈中弹出首个节点进行递归，在检测到每个落在内核中光子时将其保存下来，整个递归过程持续直到栈中不存在节点，待搜索过程结束便可搜集了所有位于估算内核中的光子。如果估算内核半径超过了单个多边形的边界，则从估算内核半径所包含的临近多边形所对应的kd-树组织的光子集合按照上述递归搜索规则进行查找，直到找到所有估算半径之内满足条件的光子的集合。

平衡Kd-tree是紧致而高效的，能够确保使用Kd-tree在n个光子中定位m个光子的时间复杂度为Ο(mlog(n))。由于这m个光子处在临近的空间内，实际中寻找m个光子的时间会进一步减少，从而保证光子映射算法的计算效率。如果发射的光子数目趋近于无穷大，辐射亮度估算中用于计算的辐射亮度估算面积就无限趋近于0，此时的辐射亮度计算会趋近于正确的结果。

步骤七：基于几何的精确辐射亮度估计

7-1.辐射度计算

首先，寻找光子映射中携带有局部几何信息的物理量。尽管光子密度本身携带有局部几何信息，但是计算光子密度本身是不精确的过程。虽然在光子密度层面不能保证精确，但是在比光子密度更“精细”的程度上考虑几何信息是可以保证精确的，局部表面上的辐射度E就是这样的物理量，辐射度的定义决定了它受到辐射通量这一表征光子的量以及局部面积这一表征相交点几何特征的量共同作用。

本发明重新定义辐射度E，它强烈的依赖于相交点附近的光子分布，同时能够近似携带相交点所在平面的拓扑与边缘信息。由Voronoi多边形的概念，在光子追踪过程中，我们为分布在同一网格上的光子建立Voronoi图的几何结构。按照Voronoi图的定义，每个Voronoi多边形中的光子是距离Voronoi多边形内点的最近的光子，其抽象意义为，该光子控制了Voronoi多边形范围的物体表面。可以近似地将相交点x的辐射密度量表现为相交点所处Voronoi多边形内的辐射度，定义Voronoi单元内辐射度E_(i,p)的概念如下：

$E_{(i,p)}=\frac{{\mathrm{\&Delta;\&phi;}}_p(x,{\stackrel{\&RightArrow;}{\mathrm{\&omega;}}}_p)}{{\mathrm{\&Delta;A}}_p}$

这里，是光子p所携带的有效辐射通量，是光子p的入射方向，而△A_p是包含光子p的Voronoi多边形的面积。基于Voronoi的划分策略，辐射通量E_(i,p)是局部几何特征的间接表达。而且，E_(i,p)的用处不仅局限于获得精确的估算内核面积。

如图5所示，图中相交点x的估算内核范围内的Voronoi多边形被设置为浅灰色，精确估算内核也应该选择浅灰色Voronoi多边形的面积之和。

本发明方法是将光子图中离散的辐射能量转化为连续均匀的物体表面的辐射亮度分布。辐射度E_(i,p)是一种离散的能量分布表达，如果为相交点x引入更多的辐射通量，可以得到更加柔和的辐射亮度结果，如图5所示。图中处在相交点x的估算内核范围内的Voronoi多边形被设置为浅灰色，那么相交点x的辐射亮度应该选择浅灰色区域所附加光子的辐射量来进行辐射亮度估计，并且相交点x的精确估算内核也应该选择浅灰色Voronoi多边形的面积之和，从图中可明显看出精确估算内核的面积是要小于完整圆盘的面积πr²。

连续均匀的辐射度可以由如下公式计算：

$\tilde{E}=\frac{{\mathrm{\&Sigma;}}_{p=1}^n{\mathrm{\&Delta;\&phi;}}_p(x,{\stackrel{\&RightArrow;}{\mathrm{\&omega;}}}_p)}{{\mathrm{\&Sigma;}}_{p=1}^n{\mathrm{\&Delta;A}}_p}$

公式中n个Voronoi多边形的面积和为该估算面积比内核估算面积πr²更加精确。如果发射足够多的光子，可以在相对较小的局部区域中收集到大量的辐射通量。因此，该方法可以获得均匀柔和的辐射亮度估计结果。

7-2.辐射亮度计算

引入了对分布在物体表面光子的Voronoi图划分之后，我们将Voronoi多边形的信息附加在光子上，一个光子不仅包含了所携带的辐射通量，光子所在的空间位置信息和入射方向等PM算法中包含的数据，同时也会附加其所在物体表面Voronoi多边形的信息，对于相交点x附近使用半径r收集到的光子集合，可以方便地取出其对应的Voronoi多边形面积，这样光子的辐射量与光子所在精确估算内核面积之间具备一对一的关系，精确面积计算可以在光子收集的过程中以非常低的时间消耗完成。在精确估算面积的基础上，本发明所提出的辐射亮度估计公式：

$L_s(x,\stackrel{\&RightArrow;}{\mathrm{\&omega;}})\&ap;\frac{\underset{p=1}{\overset{n}{\&Sigma;}}{f}_r(x,\stackrel{\&RightArrow;}{\mathrm{\&omega;}},{\stackrel{\&RightArrow;}{\mathrm{\&omega;}}}_p){\mathrm{\&Delta;\&phi;}}_p(x,{\stackrel{\&RightArrow;}{\mathrm{\&omega;}}}_p)}{\underset{p=1}{\overset{n}{\&Sigma;}}{\mathrm{\&Delta;A}}_p}$

其中，n是相交点x附近用于辐射亮度估计的光子数目，是出射方向；是光子入射方向；是双向反射分布函数。与传统光子映射方法相比，Voronoi多边形面积之和比内核面积πr²更加精确。给出了基于局部表面几何特征的辐射亮度估计解决方案，可以很好地适应具有复杂局部特征的表面的光照明计算。

步骤八：

经过上述的计算光线跟踪相交点的辐射亮度之后，采用把亮度转换为相应的色彩值的计算，就可以生成最终的真实感图像。

本发明采用了多个国际常用的光子映射绘制算法的测试场景进行实验，其中既包含了以漫反射为主的场景，也包含了以焦散为主的场景，以及多种形式混合的场景。无论何种场景，实验结果均表明，本发明算法的辐射亮度计算都可以得到稳定的数值并合成出真实感的绘制结果。图6为经典CornelBox绘制效果对比，其中(a)图为传统光子映射方法，(b)图为本发明方法的结果。图7为本发明绘制的场景的真实感效果图。

以上实施例仅用以说明本发明的技术方案而非对其进行限制，本领域的普通技术人员可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明的精神和范围，本发明的保护范围应以权利要求书所述为准。

Claims (6)

1.一种Kd-树与Voronoi图混合的精确辐射亮度计算方法，其步骤包括：

1)针对待绘制的三维场景，构造场景拓扑结构及其Kd-树层次组织；

2)光源发射光子并生成光子图，以多边形面片为单位存储该光子图，包括存储光子的基本属性以及存储光子所停留的多边形面片的唯一标识，同时每个多边形都记录落在其中的光子的唯一标识；

3)遍历三维场景表面的每个多边形，针对每个多边形，收集属于该多边形的光子并组织成二维Kd-树结构；

4)对每个多边形表面构造Voronoi图划分，并计算划分出的Voronoi单元的辐射通量密度；

5)计算从视点发射的射线与场景的相交点x；

6)根据相交点x所属的多边形面片，利用该多边形面片所保存的二维Kd-树结构收集临近的光子；

7)对收集到的光子，根据Voronoi图划分找到这些光子所覆盖的精确的面积，并利用辐射亮度计算公式计算出x点的辐射亮度。

2.如权利要求1所述的方法，其特征在于：根据步骤7)计算出的辐射亮度合成出最终的真实感图像。

3.如权利要求1所述的方法，其特征在于：步骤1)构造场景拓扑结构的方法是：首先为每个场景中的节点建立邻域网格结构，通过遍历一次场景中的网格结构，得到从网格的节点到邻域网格集合的对应关系表；之后为每个网格建立1-邻域网格集合，使用的方法是针对一个网格，遍历它的所有节点，选取每个节点邻域网格集合的并集作为该网格的1-邻域网格集合；待所有网格的1-邻域网格集合创建完成之后，递归地创建网格的2-邻域网格集合，依次类推。

4.如权利要求1所述的方法，其特征在于：步骤4)构造Voronoi剖分时，首先利用Voronoi图计算每个采样光子所在Voronoi单元格的面积大小，然后利用快速邻近查找到以相交点为圆心以内核半径为半径的圆盘内的所有光子，并获取其有效的辐射度以及Voronoi单元格的面积。

5.如权利要求1所述的方法，其特征在于：步骤6)使用栈结构收集临近的光子。

6.如权利要求1所述的方法，其特征在于：步骤7)采用如下公式计算辐射亮度：

$L_s(x,\stackrel{\&RightArrow;}{\mathrm{\&omega;}})\&ap;\frac{\underset{p=1}{\overset{n}{\&Sigma;}}{f}_r(x,\stackrel{\&RightArrow;}{\mathrm{\&omega;}},{\stackrel{\&RightArrow;}{\mathrm{\&omega;}}}_P){\mathrm{\&Delta;\&phi;}}_p(x,{\stackrel{\&RightArrow;}{\mathrm{\&omega;}}}_p)}{\underset{p=1}{\overset{n}{\&Sigma;}}{\mathrm{\&Delta;A}}_p},$

其中，x是视点光线与物体表面相交点；n是相交点x附近用于辐射亮度估计的光子数目；是出射方向；是光子p的入射方向；是双向反射分布函数；是光子p所携带的有效辐射通量；△A_p是包含光子p的Voronoi多边形的面积。