基于光流场动态参与介质绘制算法

王元龙

摘要：为了解决动态参与介质场景连续帧的实时绘制问题，提出了一种基于光流场动态参与介质场景的绘制算法。首先，采用区域匹配的方法计算关键帧之间的光流场；然后，通过插值的方法计算中间帧之间的光流场，采用帧间光流连贯性函数表示帧与帧之间光流一致性的度量，保证帧与帧之间介质运动不会发生突变；最后，按照所产生的光流场，绘制连续动态变化的参与介质场景。在连续5帧动态参与介质场景绘制中，所提算法比基于径向基函数模型（RBF）参与介质的光子映射算法效率提高了近3倍，能够达到连续帧的实时绘制，且绘制质量比较高。

关键词：多散射；参与介质；光子映射；并行计算；光流场

中图分类号：TP391 文献标志码：A

Abstract： In order to achieve the realtime rendering of continuous frame for participating media scene， a rendering algorithm based on optical flow was proposed. First， the regional matching method was used to calculate the field of optical flow between key frames. Then the field of optical flow between intermediate frames was calculated by the interpolation method， and the optical coherence function between frames was used to denote the consistency of optical flow to guarantee that the media motion wont be suddenly changed. Finally， the dynamic scene of continuous frames was rendered according to the field of optical flow. In the participating media scene rendering for 5 continuous frames， the efficiency of the proposed algorithm increase nearly five times than that of based on Radial Basis Function （RBF） model， realtime rendering of consecutive frames is implemented and rendering quality is relatively high.

Key words：multiple scattering； participating media； photon mapping； parallel computing； optical flow field

0 引言

绘制介质参与的场景时，通常需要考虑光照在物体表面或介质中的散射现象。因此，在真实感图像渲染领域，散射效果的模拟（像烟雾、云和大部分非金属表面等等的绘制）成为越来越被研究人员关注的课题。基于光线跟踪的蒙特卡罗方法能够渲染高质量的介质参与场景，模拟非均匀介质的不同属性[1-3]。但是，由于在介质中散射到处发生，需要跟踪多次散射路径直到光能被吸收或路径离开参与介质，因此该方法渲染效率低，渲染一幅简单场景图像需要几分钟。

研究者提出了光子映射的方法[4-7]来提高介质参与场景的绘制效率。光子映射方法最初被应用于全局光照下没有参与介质的场景渲染中[4]。Jensen等[5]把全局光照的光子映射方法扩展到有参与介质的场景中，与文献[4]中的方法一样，用一个球形邻域来搜索影响该点的光子。该方法实现了各种全局光照效果，但要达到好的渲染结果需要十几万的光子，不但计算量比较大，而且容易出现搜索球叠加或存在搜索不到区域的情况。Jarosz等[6]对文献[5]中的辐射度估计作了改进，提出一种用光子束（Photon Beam， PB）来估计内散射辐射度的方法，他们根据光子辐射度估计方法重新推出了在介质中的光线传输等式，整体估计每一条从视点发出光线的辐射度，用可变核（Variable Kernel， VK）的方法来搜索最近的光子和用层次包围盒（Bounding Box Hierarchy， BBH）构建光子贴图。Jarosz等[7]把PB的概念不仅应用到估计辐射度的查询操作上，还应用到光子的数据表示中，提高了渲染图像的质量。但是，PB方法需要考虑光束中的所有光子，所以计算量仍然很大，且对光束宽度采用确定的值，造成了大的偏差，虽然可以通过减小宽度来减小偏差，但增加了噪声。最近，Zhang 等[8]提出一种预处理的辐射度传输方法来实现介质参与场景全局光照效果的实时渲染。该方法能够快速地估计视点光线的辐射度，但该方法需要的预计算时间较长（73～395s）。文献[9]也提出了基于径向基函数（Radial Basis Function， RBF）模型参与介质的光子映射算法。该算法用存储RBF模型来代替存储光子从而减少存储空间，实现了实时渲染。但在整个光子映射过程中需要的预计算时间较长（0.81～5.61s）。目前，光子映射算法存在预计算时间较长的问题，难以实现动态介质参与场景连续帧的实时绘制。

本文在文献[9]的基础上，将光流场概念扩展并应用到动态介质参与场景连续帧的绘制中，提出一种基于光流场动态介质参与场景的绘制算法。首先，采用区域匹配的方法计算关键帧之间的光流场，通过最小化残差找到帧与帧之间特征对的光流矢量；然后，通过插值的方法来计算中间帧之间的光流场，采用帧间光流连贯性函数来表示帧与帧之间的光流一致性的度量，保证帧与帧之间介质运动不会发生突变；最后，按照所产生的光流场，绘制连续动态变化的介质参与场景。另外，在CPU+GPU的框架下实现了本文算法，进一步加速动态参与介质的绘制。

1 相关工作

文献[9]中提出了基于RBF模型参与介质的光子映射算法。在构建RBF模型之前，首先需要计算基函数的中心，由于聚类个数不确定，因此采用自组织映射无监督聚类方法[10]来计算基函数的中心。在辐射度估计阶段，使用构建好的RBF模型来计算穿过参与介质的每一条视点光线辐射度。用式（1）计算视点光线的内散射辐射度：

2 基于光流场的动态介质参与场景绘制算法

大部分光流算法的研究是针对动态场景，从图像序列中得到运动目标的运动矢量[11-12]。光流场是将三维动态场景投影到二维图像平面上，利用图像序列的灰度变化、空域变化和相关性等信息计算得到的运动矢量场。本文把光流场的概念应用到绘制介质参与场景的光子映射算法中，来模拟介质的动态变化效果。首先介绍算法流程，然后分别阐述算法中采用的关键技术。在介绍算法之前，定义两个关于帧的定义。

关键帧 在动态显示场景时，场景中对象变化较大的帧。

中间帧 在关键帧之间，对象变化比较平滑的帧。

2.1 算法流程

本文算法经过并行设计，使用CPU+GPU的框架完成，充分利用了GPU并行计算能力，其主要步骤如下：

步骤1 分别计算关键帧场景的RBF模型。本文采用文献[17]的方法计算参与介质的RBF模型。

步骤2 计算关键帧之间光流场。将关键帧的RBF模型作为基本元素，通过区域匹配的方法计算关键帧之间光流场。本文把RBF的中心和大小作为特征进行匹配，找到帧与帧之间特征对，根据特征对计算得到光流场。

步骤3 采用插值的方法计算中间帧之间的光流场。本文采用帧间光流连贯性函数来表示帧与帧之间的光流一致性度量。

步骤4 渲染最终的动态介质参与场景连续帧。

注 本文在GPU显存中存储第一帧场景的RBF模型及所有帧间光流场，第一帧以外其他帧场景的RBF模型通过帧间光流场计算得到。

2.2 关键帧之间光流场的计算

按照数学方法划分，计算光流场的方法分为4种：微分法、区域匹配法、基于能量的方法和基于相位的方法，由于区域匹配法计算简单、实时性强等特点，得到了较为广泛的应用[13]。本文采用区域匹配法计算关键帧间的光流场。区域匹配法认为上一帧中某个区域，在下一帧中会形成对应的模式，光流场定义为两帧中相应区域块之间产生最佳拟合的位移矢量，实质上是在前帧与后帧图像之间实施的一种对应任务，这类方法是将光流定义为不同时刻图像区域的位移量d=（dx， dy），使得不同时刻的图像区域匹配最佳。匹配图像区域常常等价为最小化残差的过程，给定两帧图像I和J，需要找到位移量d= （dx， dy）使残差最小：

如果用聚类算法将像素点匹配转化为聚类特征点会大大减少计算量[14]。本文在光子映射阶段采用自组织映射无监督聚类方法把光子的分布用RBF来模拟（详见文献[9]），因此将RBF看作区域块来搜索相应的位移矢量场。在连续关键帧I和J匹配过程中，第I、J帧分别用RBF的半径作为匹配值，记为RI（i， j，k）、RJ（i， j，k），其中（i， j，k）为RBF的中心位置。以此中心建立一个搜索窗口W（wi，wj，wk），搜索窗口不能太小也不能太大，基于统计的实验结果，本文采用大小为6×6×6的搜索窗口。与计算图像光流场不同（光流场通常表示物体点的三维运动在二维图像上的表现，而本文计算三维空间中介质运动的光流场），本文在J帧中（i， j，k）位置附近搜索与I帧（i， j，k）对应的RBF，三维空间中计算光流场相应的残差定义为ε（d）=（RI（i， j，k）-RJ（i+di， j+dj，k+dk））2，通过最小化残差的方法计算得到相应的光流矢量，d（di，dj，dk）就是需要计算得到的光流矢量。计算光流场如算法1所示。

算法1 从关键帧动态介质中进行光流场的计算。

2.4 GPU并行实现

本文算法使用计算统一设备架构（Compute Unified Device Architecture， CUDA）来完成，算法流程中的步骤都经过了并行设计，可以并行地在GPU上运行。设计CUDA程序首先需要尽可能地发掘GPU体系结构的特点：减少分支操作、减少读写操作、充分利用线程资源等。

在计算关键帧之间光流场阶段，本文先将表示光子分布的RBF数据绑定到GPU上的纹理存储器。因为GPU上的纹理硬件有一个很有用的特性，它提供了一个纹理缓存，因此纹理存储器的平均访问速度比全局存储器更快。存储RBF采用层次包围盒（BBH）结构，加速光线与RBF求交的过程（详见文献[9]）。对于光流场计算，式（4）中要求在搜索窗口W（wi，wj，wk）范围查找匹配的RBF，本文设计了在BBH结构中搜索对应RBF的方法，如图3所示为在BBH中搜索具体过程的示例。首先在后一帧中找到前一帧RBF对应的位置，该位置可能存在RBF，也可能不存在，因此在图3中用虚线框表示，如果存在RBF，用式（4）计算前一帧与后一帧中该RBF的残差；然后在搜索窗口 W（wi，wj，wk）范围中搜索使残差最小的RBF，图3中w-1、w+1等表示向前和向后搜索，其中“-”表示向前搜索，“+”表示向后搜索。搜索过程为：1）如果位置（i， j，k）为其父节点的右（左）节点，那么向前（后）搜索其左（右）兄弟节点；2）如果没有达到搜索范围继续搜索其父节点的兄弟节点，直到达到搜索范围停止搜索。

在辐射度估计阶段（即绘制阶段），估计每条视点光线的辐射度，需要求得与该视点光线相交的所有RBF，这个过程称为光线遍历。由于每条视点光线是彼此独立的，因此遍历方法实现的一个重要特点是高度并行性（在GPU上并行的跟踪每一条光线，设置的视平面大小为800×600，而光线的设置为从视点分别到视平面的每个像素的方向，所以跟踪的光线也为800×600条，从而充分利用线程资源）。每一帧的辐射度估计，先从视点并行地向成像平面上每个像素发射视点光线，然后用当前帧的RBF模型来估计辐射度，在最后显示时，把辐射度值转化为颜色值，再将其绑定到纹理对象上，通过OpenGL缓存技术将其最终的图像显示出来。但每帧得到RBF模型的方式不同，在绘制第一帧的过程中直接使用光子跟踪阶段得到的RBF模型来进行辐射度估计，而在后续帧的绘制过程中需要先通过光流场计算得到当前帧的RBF模型。

3 实验结果及讨论

本文实验使用：PIV 2.83GHz 四核CPU，4GB主存，NVIDIA GeForce GTX470显卡的PC机，编译环境为VS2008和CUDA2.3，所有实验结果图像的分辨率均为800×600。介质数据与文献[9]相同，介质中每个位置上的密度值随介质位置的不同而不同，能很好地模拟不均匀介质的散射。光流场在GPU纹理存储器中采用一个四维数据类型float4表示（在GPU中， float4的读取操作比4个float分开读取获得更高的效率），其中前三维用来表示光流场的矢量的方向，第四维用来表示RBF半径的变化率（即矢量的大小），计算光流场的平均时间大约为1s。

为了验证搜索窗口大小对光流场的影响，本文计算不同大小搜索窗口得到光流场，以及通过光流场绘制了后一帧介质的散射效果，如图4所示。图4（a）、4（b）分别为文献[9]中基于RBF光子映射算法绘制的前后两帧介质参与场景，图4（c）、4（d）、4（e）分别是不同大小搜索窗口计算得到的光流场以及采用该光流场绘制后一帧的介质散射结果。从图4中可以看出，对于基于光流场算法绘制后续帧的介质散射效果算法，采用大小为6×6×6的搜索窗口能够达到与采用预处理算法[9]效果类似的绘制结果，但效率提高了大约20倍。表1为基于RBF 光子映射算法和基于光流场算法绘制后一帧介质参与场景的运行统计数据， 其中，基于RBF 光子映射算法需要在每一帧进行预处理操作，而基于光流场算法在绘制过程中使用光流场计算后一帧相于对前一帧的RBF， 因此可以提高连续帧的绘制效率。正如图4和表1所示，搜索窗口越大绘制效果越好，同时绘制时间越长。采用大小为6×6×6的搜索窗口就能够达到与文献[9]类似的较好效果，且能够实现连续动画的效果，因此本文采用大小为6×6×6的搜索窗口实现连续帧的实时绘制。

图5为在关键帧之间插值生成中间帧光流场的示意图，其中，第1帧和第5帧视为关键帧，第2、3、4帧为插值计算得到的中间帧，其中xy表示从第x帧到第y帧之间的光流场。图6所示为通过图5光流场产生的连续5帧动态介质参与场景的绘制结果，第1帧采用文献[9]中基于RBF光子映射算法绘制，后面的连续帧采用基于光流场算法绘制。从图中可以看出，基于光流场算法绘制的连续多帧介质散射效果质量比较高，而且可以实现帧间实时的绘制。其中基于BRF光子映射算法[9]，在每一帧都需要进行预处理操作，每一帧平均时间为5.31s，5帧共26.55s。而基于光流场算法仅仅需要在第1帧进行预处理操作（5.31s），后续4帧通过光流场计算各帧RBF模型（0.26×4=1.04s），总计6.35s。因此在连续5帧动态介质场景绘制中比基于RBF光子映射算法效率提高了近3倍。

4 结语

本文针对动态介质场景连续帧的绘制，提出一种基于光流场动态参与介质绘制算法，其借用了运动分析中光流场的方法进行帧与帧之间相关性的分析，并且采用CPU+GUP的并行框架，实现了连续帧的实时绘制。从实验结果可以看出，用本文算法模拟的效果与利用基于RBF光子映射算法绘制的效果相当，但本文算法提高了绘制连续帧的效率。

本文在插值计算中间帧光流场的过程中，需要尝试不同的迭代值找出最佳的迭代次数，不能应用于所有介质场景的自动绘制。将来要重点解决的问题是考虑设置自动判别迭代次数的机制。

参考文献：

[1]KAJIYA J， HERZEN B. Ray tracing volume densities[C]// Proceedings of the 11th Annual Conference on Computer Graphics and Interactive Techniques. New York： ACM， 1984： 165-174.

[2]MAX N. Efficient light propagation for multiple anisotropic volume scattering[C]// Proceedings of Eurographics Workshop on Rendering. Heidelberg： Springer， 1994： 87-104.

[3]JENSEN H， MARSCHNER S， LEVOY M， et al. A practical model for subsurface light transport[C]// Proceedings of the 28th Annual Conference on Computer Graphics and Interactive Techniques. New York： ACM， 2001： 511-518.

[4]JENSEN H W. Global illumination using photon maps[C]// Proceedings of the 7th Eurographics Workshop on Rendering Techniques. Berlin： Springer， 1996：21-30.

[5]JENSEN H W， CHRISTENSEN P H. Efficient simulation of light transport in scenes with participating media using photon maps[C]// Proceedings of the 25th Annual Conference on Computer Graphics and Interactive Techniques. New York： ACM， 1998：311-320.

[6]JAROSZ W， ZWICKER M， JENSEN H W. The beam radiance estimate for volumetric photon mapping[J]. Computer Graphics Forum， 2008， 27（2）：557-566.

[7]JAROSZ W， NOWROUZEZAHRAI D， SADEGHI I， et al. A comprehensive theory of volumetric radiance estimation using photon points and beams[J]. ACM Transactions on Graphics， 2011， 30（1）：1-19.

[8]ZHANG Y B， ZHAO D， MA LK. Realtime volume rendering in dynamic lighting environments using precomputed photon mapping[J]. IEEE Transactions on Visualization and Computer Graphics， 2013， 19（8）：1317-1330.

[9]WANG Y L， GUO P. Rendering scenes with participating media based on RBFs for photon mapping using graphics hardware[C]// Proceedings of the 12th ACM SIGGRAPH International Conference on VirtualReality Continuum and Its Applications in Industry. New York： ACM， 2013： 167-176.

[10]VESANTO J， ALHONIEMI E. Clustering of the selforganizing map[J]. IEEE Transactions on Neural Networks， 2000， 11（3）： 586-600.

[11]IKETANI A， NAGAI A， KUNO Y， et al. Realtime surveillance system detecting persons in complex scenes[J]. RealTime Imaging， 2001， 7（5）： 433-446.

[12]潘金山，苏志勋，王伟.运动细节估计的光流场方法[J].计算机辅助设计与图形学学报，2011， 23（8）：1433-1441.（PAN J S， SU Z X， WANG W. Optical flow method for motion details estimation[J]. Journal of ComputerAided Design & Computer Graphics，2011，23（8）：1433-1441.）

[13]艾海舟，苏延超. 图像处理、分析与机器视觉[M]. 北京：清华大学出版社， 2003： 542-548.（AI H Z， SU Y C. Image Processing Analysis and Machine Vision[M]. Beijing： Tsinghua University Press， 2003： 542-548.）

[14]林强，田桂兰. 聚类匹配法光流检测中共同运动匹配策略的研究[J].电子学报，1997， 25（1）：58-61.（LIN Q， TIAN G L. An investigation on the matching tactics of common motion in detecting flow by cluster algorithm[J].Acta Electronica Sinica，1997， 25（1）：58-61.）