基于光流場動態參與介質繪制算法

王元龍

摘要：為了解決動態參與介質場景連續幀的實時繪制問題，提出了一種基于光流場動態參與介質場景的繪制算法。首先，采用區域匹配的方法計算關鍵幀之間的光流場；然后，通過插值的方法計算中間幀之間的光流場，采用幀間光流連貫性函數表示幀與幀之間光流一致性的度量，保證幀與幀之間介質運動不會發生突變；最后，按照所產生的光流場，繪制連續動態變化的參與介質場景。在連續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.）