基于栅格的全局光子图重建

蔡 鹏, 孔德慧, 尹宝才

（北京工业大学计算机学院多媒体与智能软件技术北京市重点实验室，北京 100124）

基于栅格的全局光子图重建

蔡 鹏, 孔德慧, 尹宝才

（北京工业大学计算机学院多媒体与智能软件技术北京市重点实验室，北京 100124）

基于光子图的光子映射算法能产生高质量的照片级图像。对于光照复杂的场景，光子图需要存储大量光子以提高生成图像的质量，这不仅占用大量的内存空间，而且光照估计的时间长。论文提出基于栅格的全局光子图重建的算法，即在光子包围盒被栅格化后，其非空栅格中一定比例的光子被用来重建新的光子图，并保证重建前后栅格内光子能量和守恒，这使得重建前后光子图的光照估计的效果相近。通过增加特定栅格中的重建光子数目，能有效减少由几何偏差引起的光照估计误差，增强直接聚焦（焦散）和间接聚焦光照的绘制效果；并使用简单方法检测生成图像中少量噪声，增加少量采样即可有效减少相应的噪声。全局光子图重建算法的计算成本低，并保持生成图像的视觉独立性。

光子映射；全局光子图重建；几何偏差；焦散

真实感图形学在虚拟现实、动画影视和系统仿真等领域中得到广泛应用。光线跟踪[1-2]能计算阴影、光的反射和折射以及其他全局光照效果。但对于辉映和聚焦（或通过镜面所间接看到的聚焦）光照效果的计算，光线跟踪是非常低效的。Jensen等[3-6]提出基于光子图的光子映射技术，对于光照复杂的场景，它需要存储大量光子在光子图中，以使光照估计的效果接近真实光照。这将占用大量的内存空间，并且它需要在所存储的光子中定位最近的若干光子，故其计算成本很高。

本文提出基于栅格的全局光子图重建方法，相对于使用相同光子数目的光子图所生成的图像，利用重建后的全局光子图进行光照估计所生成的图像，具有更少的噪声，并减少几何偏差引起的光照估计的误差，增强直接聚焦（焦散）和间接聚焦的绘制效果；通过简单的方法，检测出生成图像易产生噪声的区域，并增加少量采样即可减少相应的噪声。本文的全局光子图的重建时间少，并保持生成图像的视觉独立性。

1 以前的工作

Jensen提出的具有两个子阶段的标准的光子映射方法[6]。

第一个子阶段：利用光子跟踪建立光子图，即从场景的光源(L)出发跟踪大量的光子，构造高密度的焦散光子图和低密度的全局光子图。当光子击中漫反射表面(D)或略微光滑表面(G)时，存储光子信息到相应的光子图中，并决定光子是否被吸收或反射。焦散光子图存放的光子至少被镜面反射或折射(S)一次，即它经历的路径为L(S)+(D︱G)。全局光子图存储所有经历各种路径(L(D︱G|S)*(D︱G))的光子（+、*、︱分别代表≥1次、≥0次、或）。

光子图存储的光子i，具有能量或通量Φi、输入方向ωi和位置pi等信息。场景中光照估计位置x的反射方向ω的光照Lr，可以通过在光子图中定位距x最近的M个光子来估计：

其中，fr是双向反射分布函数(BRDF)，rM(x)是x与M个光子中距x最远的某个光子的距离。这对应最近邻居密度估计[7]，通常在球体范围内查找最近的若干光子。为了有效地定位最近的光子，光子图采用平衡二叉树作为光子搜索的数据结构。

第二个子阶段：使用光子图进行光照估计，即从图像每个像素出发的光线，在场景中经过多次或零次的镜面反射或折射后，击中漫反射表面时，基于式(1)计算交点位置的光照。算法以两种方式重建光照，第一种方式是先用局部光照公式计算直接光照，焦散光子图计算焦散光照，全局光子图计算间接光照，然后进行综合而得出实际的光照值；第二种方式用全局光子图计算场景中所有的全局光照，当光子数目越多，则光照计算的精度越高。本文只用全局光子图进行光照计算。

Peter和 Pietrek[8]提出基于重要性驱动的光子图构造算法，即用重要性图引导更多的光子到场景中视觉重要的区域。基于这种重要性驱动的光子跟踪，构造一个包含场景中全局光照信息的光子图。然而，在产生的光子图中，高能量的光子与低能量的光子相混合，造成光照估计中光照不一致现象增多。

Suykens和Willems[9]介绍了光子密度控制的方法，避免光子存储在过于密集的或不重要的区域；并提出基于视点重要性的方法，但全局光子图的视觉独立性被丢失，且计算成本高（包括上面Peter等的方法）。

Schregle[10]提出基于光子图的光子密度的光照估计的偏差补偿方法。他给出了最小带宽（即对应的光照估计的偏差最小）和最大带宽（即噪声最小），根据噪声的高斯分布的概率统计，在最小带宽和最大带宽之间，利用二分折半搜索得到最佳带宽（光照估计中的最近的光子数目），以使光照估计的偏差和噪声均最小。

Herzog等[11]在光线空间内，使用光线Splatting技术对具有漫反射和略微光滑表面的场景进行光子密度估计，以有效地计算间接光照。此方法能很好地减少由几何偏差引起的光照误差。

Spencer等[12]提出层次性的光子映射，他们通过光子的能量传递，将平衡二叉树的下层的光子能量逐层传递到上层的光子能量中，以提高光照估计的效率。他们用光线的足迹计算光照估计的面积，而不是用光子密度得到它；并根据光照估计的区域的特征，动态地调整光线的足迹，大大减少光照估计的噪声。

Hachisuka等[13]提出渐进式的光子映射，它是一个多道处理过程，即光线跟踪子阶段和许多光子跟踪子阶段。在有限内存的空间内，渐进的光子映射能使图像的像素光照值收敛到真实的光照值。其缺点是当切换视点时，需要重新进行时间很长的光子跟踪过程。

2 基于栅格的全局光子图重建

2.1 基本方案

在光子映射的光子跟踪阶段，场景中的光源发出大量的光子，其中与漫反射表面交互的光子被保存在全局光子图（以下简称为光子图）中。

本文先计算光子图中所有光子的包围盒，将其均匀分割成一定数目的栅格，并按光子产生的顺序将其存储到相应的栅格中。

如图1所示，在每个非空栅格中取出前面的若干光子（按光子的存放顺序）以用来重建光子图，并按重建前后栅格内光子能量和守恒的原则计算光子能量。设定的光子包围盒的栅格数目为Ngd，则计算光子包围盒的栅格结构如下：

其中，V为光子包围盒的体积，dx、dy、dz分别为包围盒的棱长，Nx、Ny、Nz分别为包围盒的x、y、z轴上的栅格数目。

图1 基于栅格的光子图重建的基本方案

基于栅格的全局光子图重建的基本方案是在重建前的每个非空栅格中，按光子的存放顺序以相同比例取出前面的若干光子进行光子图重建，即：

其中，Nold为重建前光子的数目，Nset为设定的重建光子的数目，Nnew为实际的重建光子的数目。Nold[i]和Nnew[i]分别为第i个非空栅格的重建前和重建后光子的数目，Ngrid为非空栅格的数目。若Nnew[i]为0，则Nnew[i]置为1。重建后的栅格中Nnew[i]个光子对应重建前相应栅格的Nold[i]个光子中的前Nnew[i]个（按光子的存放顺序）。通常，Nold远大于Nset和Nnew，Nnew接近Nset。

在第i个非空栅格中，记光子的红、绿、蓝3个通道的能量，在重建前分别为Φr0[j]、Φg0[j]、Φb0[j]，在重建后分别为Φr1[j]、Φg1[j]、Φb1[j]。对于(Nold[i]-Nnew[i])个非重建光子，其相应通道的能量平均分配到Nnew[i]个重建光子的相应通道的能量上，以保证栅格内光子能量和守恒，具体计算过程, 如式(5)所示。通过设置适当的栅格尺寸，大多数相邻的非空栅格的光子的能量变化不大。因此，用其中很少的光子来重建光子图，其光照估计的效果与基于重建前光子图的光照估计的效果相近。

2.2 优化方案

论文在基于栅格的全局光子图重建的基本方案基础上，为达到减少几何偏差引起的光照误差和增强焦散光照的绘制效果的目的，提出相应的优化方案。

对于曲率大的漫反射表面的光照估计，在球体范围内查找最近的若干光子，容易从不属于光照估计范围的表面引入光子，如图2(a)和图2(b)的左部。通常的做法是在光照估计的表面法线方向上将光照估计范围的球体压成椭球体或圆盘，如图2(a)和图2(b)的右上部，但存在其投影面积不能正确表示光子所覆盖面积的问题，通过凸包算法计算面积或利用滤波减小相应的光照误差[6]。

图2 在场景中曲率大的漫反射表面某位置光照估计

上述做法在某种程度上减小光照估计中由几何偏差引起的光照误差，但同时大大增加光照估计的计算成本。与之不同，论文增加对应场景中漫反射表面的曲率大的栅格中的重建光子，如图2(a)和图2 (b)的右下部所示，则曲率大的位置的光子密度变大，在光照估计的球体范围内包含同样光子数目的情况下，则球体范围变小，即从不属于光照估计范围的表面引入光子的数目变小，从而达到减小此类误差的目的，其优点是在光照估计中不增加额外的计算成本。特别是，对于光子泄露[11]，传统的椭球体估计范围包含垂直的不透明挡板另一侧的光子，增加了光照估计的误差，如图 2(b)上部，文献[11]也不能解决此问题；本文的方法能有效解决光子泄露的情况，如图2(b)下部，原因如前所述。

在场景中，聚焦光照比其它光照对视觉影响大。通常将光子图分成高密度的焦散光子图和低密度的全局光子图的做法[6]，只能增强直接聚焦(焦散)的光照效果，而不能增强间接聚焦的光照效果（即被漫反射表面反射的光子经镜面反射或折射后而产生的聚焦光照）。论文通过增加含聚焦光子（所有经镜面反射或折射后击中漫反射表面的光子）较多的栅格的重建光子数目，能增强直接聚焦和间接聚焦的光照估计的效果。

如上所述，将式(3)修改如下：

其中，Kj为栅格j中漫反射表面曲率的度量值，Sj为栅格j中聚焦光照增强的权值。如式(6)所示，若相应栅格的Kj和Sj值越大，则重建后相应栅格中的光子越多。

图3 场景中三角网格被分成小三角网格，并且小三角网格通过其中心被定位在相应的栅格中

如图3所示，以栅格边长的若干分之一为单位长度，将场景中对应的漫反射表面的每个三角网格拆分成一些相等的小三角网格。若某个小三角网格的中心在某个栅格内，即此栅格包含了该小三角网格。令njk为栅格j和与之邻近栅格的对应空间中的第k个小三角网格，m为小三角网格个数，由式(7)计算漫反射表面的曲率的度量值Kj。设ε为大于零的参数值（通常取ε= 0.03），它能调节漫反射表面的曲率度量值大的栅格(Kj＞ε)和曲率度量值为零的栅格（Kj=ε，即对应栅格空间中漫反射表面接近平面）的Kj值的相对大小。

由式(8)计算聚焦光照增强的权值Sj，其中Nold[j]和Nspe[j]分别为第j个非空栅格中重建前的光子和聚焦光子的数目，满足0≤Nspe[j]≤Nold[j]。Nold和Nspe分别为重建前的所有非空栅格中光子和聚焦光子的数目和，满足0≤Nspe≤Nold。在通常的场景中，Nspe[j]和Nold[j]均远小于Nspe和Nold。第j个非空栅格的聚焦光子数目Nspe[j]的系数比非聚焦光子数目(Nold[j]-Nspe[j])的系数大，即对Sj影响大。因此，光子数目Nold[j]越大，特别是聚焦光子数目Nspe[j]越大，则此栅格的聚焦光照增强的权值Sj越大。

2.3 图像的滤波处理

对于漫反射表面的曲率度量值大或聚焦光子比例大的栅格，基于栅格的全局光子图重建的的优化算法增加了这些栅格中重建光子数目。通过式(6)和式(5)计算，对应这些区域的光子能量与周围区域的光子能量相差比较大，在光照估计时容易产生少量的噪声。

3 实验研究与结果分析

本文给出基于栅格的全局光子图重建的实验研究与结果分析。以下所有实验是在 3GB内存、2.66GHz Intel (R)CPU的硬件平台上实现，生成图像的像素均为400×400。

采取分别用基于栅格的光子图重建的基本方案和优化方案做实验，3个场景分别为：场景1（含玻璃球、镜面球）、场景2（含两个玻璃正方体）、场景3（含具有漫反射表面的波动几何体和挡板）。3个场景的重建前光子图中均存储107个光子，重建后光子图的光子数目均稍小于5×105个。

本文设置场景的栅格数目Ngd为107，得到相应的光子包围盒的结构（X、Y和Z轴上的栅格数目）：216×216×216。3个场景的非空栅格数目Ngrid分别为：231540、231528和329572。

每个场景的光子图的重建前光子数目均为107个，依光子产生的顺序存储到相应场景的栅格中。如表1所示，列出3个场景的两种方案的相关参数（重建前的光子数目Nold、设定的重建光子数目Nset）与实际的重建光子数目Nnew和光子图重建时间（包含建立平衡二叉树的时间）。

表1 两种重建方案的参数和结果（时间：毫秒）

对于重建前具有107个光子的光子图，设定光照估计的光子包围球最多含 320个最近的光子。对于在光子数目（5×105个）大致相同的条件下进行对比的光子图和本文的两种方案的重建光子图，设定光子包围球最多含 80个最近的光子。这4种光子图的光子搜索的数据结构均为平衡二叉树。

如表2所示，列出每个场景的相应光子图的光照估计时间（即图像的生成时间）。因为，重建前光子图的光子数目为重建后的 20倍，且光子包围球的最大光子数目为重建后的4倍，所以表2中重建前光子图的光照估计时间远大于重建光子图的相应时间。

表2 场景的4种光子图的光照估计时间（时间：毫秒）

以上4种光子图的光照估计的效果，如图4～图6所示，对于这3个场景，重建前的光子数目（107个光子）比两种方案的重建后的光子数目（约为5×105个光子）的20倍稍多，但光照估计的效果差别很小。相对于本文的基本方案，其优化方案更接近重建前光子图的光照估计的效果，这是因为由几何偏差引起的光照估计误差被减小，且聚焦光照的绘制效果得到增强。图6中的场景3的左局部表明优化方案有效地减少了光子泄露（由特殊几何偏差引起的光照误差，具体解释见本文2.2小节）。在图4～图6中，基于空间分布的光子图重建的两种方案生成的光子图远比作为对比的光子图（具有大致相同的光子数目，即5×105个）的光照估计的效果好。

图4 场景1对应的4种光子图的光照估计的效果

图5 场景2对应的4种光子图的光照估计的效果

图6 场景3对应的4种光子图的光照估计的效果

图7 利用优化算法生成的光子图进行光照估计时，检测出图像中易产生噪声的区域(白色)。

在 2.3小节，提到优化方案生成的光子图，进行光照估计而生成的图像有少量的噪声，如图7所示，显示了用2.3小节的简单方法检测出的这些易产生噪声的区域。图4～图6是由优化方案生成的光子图进行光照估计而产生的图像，并在易产生噪声图像区域增加了少量的采样点，所以能有效地减少这些少量的噪声，仅比基本方案增加了少量的光照估计时间，如表2所示。

4 结 论

本文提出基于栅格的全局光子图重建方法，相对于相同光子数目的光子图的生成图像，利用重建后的光子图进行光照估计生成的图像，具有更少的噪声，并且减少由几何偏差引起的光照估计的误差，增强直接聚焦（焦散）和间接聚焦的绘制效果；以及通过简单的方法，检测生成图像中容易产生噪声的区域，并增加少量采样以有效减少这些噪声。

[1]Lafortune E P, Willems Y D. A theoretical framework for physically based rendering [J]. Computer Graphics Forum, 1994, 13(2)：97-107.

[2]Veach E, Guibas L J. Optimally combining sampling techniques for monte carlo rendering [J]. Computer Graphics Proceedings, Annual Conference Series,1995：419-428.

[3]Jensen H W. Global illumination using photon maps [J].Eurographics Rendering Workshop. New York City：1996：21-30.

[4]Jensen H W. Rendering caustics on non-lambertian surfaces [J]. Computer Graphics Forum, 1997, 16(1)：57-64.

[5]Jensen H W, Christensen P H. Efficient simulation of light transport in scenes with participating media using photon maps [J]. Computer Graphics, 1998：311-320.

[6]Jensen H W. Realistic image synthesis using photon mapping [M]. Peters A K, 2001.

[7]Silverman B W. Density Estimation for statistics and data analysis [M]. New York：Chapmann and Hall,1986.

[8]Peter I, Pietrek G. Importance driven construction of photon maps [J]. Proceedings of Eurographics Rendering Workshop. New York：1998：269-280.

[9]Suykens F, Willems Y D. Density control for photon maps [J]. Proc. Eurographics Workshop Rendering(Rendering Techniques ’00), 2000：23-34.

[10]Schregle R. Bias compensation for photon maps [J].Computer Graphics Forum, 2003, 22 (4)：729-742.

[11]Herzog R, Havran V, Kinuwaki S, et al. Global illumination using photon ray splatting [J].Eurographics, 2007, 26(3)：503-513.

[12]Spencer B, Jones M W. Hierarchical photon mapping [J].IEEE Transactions on visualization and computer graphics, 2009, 15 (1)：49-61.

[13]Hachisuka T, Ogaki S, Jensen H W. Progressive photon mapping [J]. ACM Trans. Graph, 2008, 27(5)：1-8.

A Grid-based Reconstruction of Global Photon Maps

Cai Peng, Kong Dehui, Yin Baocai
( Beijing Key Laboratory of Multimedia and Intelligent Software Technology, College of Computer Science and Technology,Beijing University of Technology, Beijing 100124, China )

Photon mapping based on photon maps can efficiently render high-quality photorealistic images. It stores many photons in photon map to improve rendering qualities of complex illumination, leading to high memory overhead and long reconstruction time. This paper presents a grid-based reconstruction of global photon maps. We divide bounding box of photons into a number of grids. Based on the number of photons, we take a certain rate of photons in grids to reconstruct global photon map, and ensure conservation of total energy of photons in each grid before and after reconstruction, which make estimated illumination effects of photon maps before and after reconstruction almost same. By adding photons in special grids, we effectively reduce illumination error due to all kinds of geometry biases and enhance rendering effects of direct and indirect caustic. This paper checks a few noises in rendering pictures by simple method, and adds a few samples to effectively reduce these noises. Our reconstruction method of photon maps takes fewer time and ensures viewpoint independence of rendering pictures.

photon mapping; reconstruction of global photon map; geometry bias; caustic

TP 391

A

2095-302X (2013)04-0010-06

2012-09-02；定稿日期：2013-03-06

国家自然科学基金资助项目（U0935004,61227004）

蔡 鹏（1977-），男，江苏句容人，博士研究生，主要研究方向为计算机图形学、真实感图形。E-mail：caipeng568@163.com

孔德慧（1968-），女，内蒙古包头人，教授，博士，博士生导师，主要研究领域为虚拟现实与图形学。E-mail：kdh@bjut.edu.cn