一種針對MLT的采樣分布改進方法*

賀懷清， 陳 帥， 劉浩翰， 計 瑜(中國民航大學 計算機科學與技術學院， 天津 300300)

信息科學與工程

一種針對MLT的采樣分布改進方法*

賀懷清， 陳 帥， 劉浩翰， 計 瑜
(中國民航大學 計算機科學與技術學院， 天津 300300)

為了解決Metropolis光線追蹤算法隨機采樣分布的稀疏性導致的生成圖像隨機噪聲問題，提出了一種基于采樣分布的突變策略.改進方法將采樣過程分為兩個階段對Metropolis光線追蹤算法進行改進，通過一階段采樣分布矩陣對二階段突變策略進行影響，在當前采樣點的8鄰域像素位置進行小范圍采樣.在滿足細節平衡條件的同時，設置多個采樣點的標量貢獻函數平均值比值為接收概率.結果表明，在相同渲染時間下改進算法生成的圖像較原方法噪聲更少，與目前較好的改進算法效果相當.通過改進突變策略的方式解決了間接光照場景中噪聲較大的問題.

Metropolis光線追蹤； 采樣分布； 全局光照； 馬爾科夫鏈蒙特卡羅； 真實感渲染； 突變策略； 接收概率； 路徑選擇

照片級真實感渲染的快速發展為電影和游戲等應用帶來了與真實世界相同的效果[1]，同時其在民航飛行模擬的視景仿真中也有重要應用.離線渲染技術以全局光照為基礎，不僅考慮光源直接照射物體的直接光照，也考慮場景中其他物體相互作用、反射和透射等間接光照，在渲染時間增加的同時也極大地提升了場景渲染的真實性.使用Monte Carlo光線追蹤對光傳遞渲染方程的解空間進行離散采樣，可以計算場景中真實光照的近似結果，從而生成具有真實感光照的場景生成圖像.基于Monte Carlo光線追蹤思想的算法主要有以下幾種：路徑追蹤(path tracing，PT)、雙向路徑追蹤(bidirectional path tracing，BDPT)、光子映射(photon mapping)、Metropolis光線追蹤(Metropolis light transport，MLT)[2].相比之下，MLT算法在具有復雜光照場景的全局光照渲染中具有高效性和魯棒性.MLT應用MCMC(Markov chain Monte Carlo)采樣方法在局部范圍內探索空間中的光線傳輸路徑，其特點為對于提議路徑只由當前路徑決定而與歷史狀態無關.因此，在找到一條有效路徑后在其局部繼續探索可能的有效路徑就可以在復雜光照場景中提高采樣效率，快速找到并生成有效路徑.后續研究者們研究MLT算法的改進，主要工作集中在突變策略[3]、重要性采樣[4]和啟動偏差[5]等方面，由此得到兩種典型的改進算法：1)PSSMLT(primary sample space Metropolis light transport)通過改變原始突變策略，優點在于使用空間變換的方法使提議樣本的生成轉化為產生一組隨機數向量[6]，但重要性采樣的選取較為單一；2)MMLT(multiplexed Metropolis light transport)結合多重重要性采樣(multiple importance sampling，MIS)，其優點在于可以通過MIS自動地自適應構建路徑[7]，并限制了子路徑的連接且算法易于實現，但其對多次反射和焦散等復雜路徑的采樣效率不高[8].

分析MLT算法發現，通過不斷地接收和拒絕提議樣本對場景進行渲染，可在生成過程中在亮度值較高的區域拒絕多次而使收斂速度降低，表現為不斷拒絕提議樣本并連續使用當前樣本.而在提議樣本的接收和拒絕過程中，高亮度區域多次重復采樣使生成的圖像噪聲較大.為此，本文改進MLT采樣分布為兩個階段來生成圖像.由于PSSMLT算法的突變策略與本文的改進突變策略無相關性，一階段由PSSMLT算法采樣并生成采樣矩陣，通過采樣矩陣影響二階段圖像生成.同時，改變小范圍擾動作為對當前采樣點8鄰域像素位置的比例采樣，并將區域的標量貢獻函數平均值比值作為接收概率.實驗證明，在具有間接光照場景中可以有效降低生成圖像噪聲，使得相同渲染時間下生成的圖像質量更好.

1 Metropolis光線追蹤算法

1.1 Monte Carlo積分

Monte Carlo積分可以通過N個樣本向量yi∈Ω(i∈N)來估計積分的值，其樣本的生成符合概率分布函數(probability density function，PDF)p(yi).重要性采樣技術使p(yi)在一定程度上成比例于g(yi)，可以降低積分估計的方差，其表達式為

(1)

(2)

1.2 路徑積分

Kajiya首次提出了渲染方程，奠定了基于物理的真實感渲染的基礎.渲染方程描述了光線在場景中反射的數學模型，通過物理方程描述場景中光線與幾何表面的相互作用.在對場景的渲染過程中，其本質是對路徑積分的求解，并由路徑積分計算最終生成圖像的輻射亮度.對于生成圖像中的每個像素j，其輻射亮度Ij可以表示為

(3)

式中：Ω′(M)為路徑空間中所有路徑的集合；z∈Ω′(M)為從光源到傳感器生成路徑的頂點向量，z=(z0，z1，z2，…，zq)，q為路徑長度，q≥1；hj(z)為每個像素j的過濾函數，其值非零；f(z)為路徑貢獻函數，由路徑向量z中各個頂點的反射算子乘積構成；dz為生成路徑的面密度.

1.3 光線追蹤算法

Veach and Guibas將MCMC采樣方法應用于光線追蹤算法，通過生成一系列采樣點xi對場景進行采樣，而新生成的提議樣本xi對當前樣本xi-1的突變生成只與當前樣本相關，與采樣過程中的歷史狀態無關.將式(2)應用于路徑積分可得

(4)

(5)

式中，Q(xi-1→y)為轉移函數，描述樣本xi-1突變為樣本y的概率密度函數.為使采樣過程快速收斂于穩態分布，要在容許的情況下盡可能增大接收概率，因此可以將接收概率較大者的值設為1.由此，提議路徑被接收的概率為a，若被接收則xi=y，否則xi=xi-1，這也被稱作Metropolis-Hastings更新規則.同時，為了保證采樣過程可以收斂，還需要滿足細節平衡條件，即

f*(xi-1)Q(xi-1→y)a(xi-1→y)=

f*(y)Q(y→xi-1)a(y→xi-1)

(6)

2 基于采樣分布的改進

2.1 改進之處

MLT算法通過設置每個像素平均采樣數量與生成圖像像素的乘積作為總采樣數量.在采樣過程中，需要在已有路徑小范圍擾動生成提議樣本，尋找路徑貢獻值更大的路徑樣本，由此，MLT算法可以通過復用已有路徑探索更為有效的路徑樣本.其中，接收概率在MLT算法中主要決定提議樣本的接收與拒絕，這使得當前樣本的選取只與當前樣本相關，而與采樣歷史狀態無關.但是，對于單一的鏡頭采樣只衡量了單一像素點的信息，沒有考慮到像素周圍的相似像素點信息，這樣也可能會導致亮度值較大的像素點被多次采樣，占用過多的采樣數量，從而減少場景中其他部分的采樣數量形成噪聲.本文算法將渲染過程分為兩個階段，第一階段為總采樣數量的3/4，第二階段為剩余采樣數量.一階段采樣采用PSSMLT算法中的突變策略，即大范圍突變采用隨機值，而小范圍突變使用指數分布函數對樣本進行小范圍擾動，指數分布函數的表達式為

yi=xi-1+te-log(t/s)ξ

(7)

式中：ξ為標量隨機變量；提議樣本yi是通過對當前樣本xi-1進行指數分布擾動生成的，其取值范圍為[s，t].同時，在一階段采樣過程中，通過生成并記錄一階段鏡頭采樣的采樣矩陣來記錄一階段的采樣分布.采樣矩陣的行列分別為生成圖像的水平像素數和垂直像素數.在二階段采樣中，大范圍突變仍然采用隨機值.在一次大范圍突變后，將小范圍突變中鏡頭突變樣本y0設置為對當前像素周圍8鄰域像素位置，其余樣本(y1，y2，…，yk)仍采用式(7)中的指數分布函數進行擾動.若當前像素的坐標為l0(m，n)，則其8鄰域像素位置可以表示為

Ω= {l0(m-1，n-1)，l0(m-1，n)，

l0(m-1，n+1)，l0(m，n-1)，

l0(m，n+1)，l0(m+1，n-1)，

l0(m+1，n)，l0(m+1，n+1)}

(8)

對于周圍8鄰域像素位置的采樣次數，則由一階段形成的采樣矩陣來確定，其值為鄰域像素位置對應矩陣值整除當前像素對應矩陣值，如圖1a所示.若當前像素位置對應矩陣值為0，為了避免除數為0，設置其值為1.在隨后的采樣過程中，每次確定一個提議樣本后，則按照二階段突變策略對提議樣本周圍進行采樣，在達到采樣數量后由接收概率決定新的提議樣本是否接收或拒絕，從而進行下一次的提議樣本及其周圍樣本的采樣過程，如圖1b所示.

圖1 本文算法采樣圖示Fig.1 Sampling diagram of present algorithm

(9)

(10)

因此，本文算法的接收概率將受當前樣本像素位置及周圍像素位置共同影響.

2.2 算法描述

結合采樣分布對Metropolis光線追蹤算法進行改進[9]，計算新的提議樣本的前提是已經解析完場景文件的參數信息并進入渲染流程，算法描述如下：

1) 初始化N個樣本，并計算其輻射亮度和標量因子b，用于估計整個圖像亮度和參與路徑積分計算；

2) 一階段令初始樣本x0為xi，并以樣本中的交點為頂點形成有效路徑；

3) 利用式(4)計算相應f(x)函數值來估計當前樣本xi的輻射亮度，并計算當前樣本對圖像的貢獻Icurrent；

4) 對當前樣本進行突變，生成新的提議樣本y，并按照步驟3)中的方式計算提議樣本y所形成路徑的輻射亮度，并計算提議樣本對圖像的貢獻Iproposed；

5) 利用式(5)計算接收概率a來決定是否接受提議樣本，但每個樣本被拒絕的次數最多不能超過M，否則接受該提議樣本；

6) 利用式(4)和標量因子b計算最終當前樣本輻射亮度并通過期望方式進行累加；

7) 如果樣本被接受，則將提議樣本作為當前樣本，如果被拒絕，保持當前樣本不變，重復步驟4)～6)，直至達到由場景中設置的一階段最大采樣數量；

8) 二階段對當前樣本進行突變，生成新的提議樣本y，并根據一階段采樣矩陣生成8鄰域像素位置樣本，按照步驟3)中的方式計算提議樣本y及8鄰域像素位置所形成路徑的輻射亮度，并利用式(9)計算提議樣本對圖像的共同貢獻SumIproposed；

9) 利用式(10)計算接收概率a來決定是否接受提議樣本，但每個樣本被拒絕的次數最多不能超過M，否則接受該提議樣本；

10) 利用式(4)和標量因子b計算二階段當前樣本輻射亮度并通過期望方式對二階段輻射亮度值進行累加；

11) 如果樣本被接受，則將提議樣本作為當前樣本，將SumIproposed的值賦給SumIcurrent，如果被拒絕，保持當前樣本不變，重復步驟8)～10)，直至達到由場景中設置的最大采樣數量；

12) 將得到的輻射亮度值記錄在虛擬膠片中，并對輻射亮度值進行值映射，記錄在exr格式的圖像中.

2.3 時空復雜度分析

設生成圖像的大小為u×v，每像素平均采樣數量為w，則總采樣數量為u×v×w.MLT算法進行隨機采樣的過程實際為對總采樣數量進行突變、連接和值映射的過程，MLT算法的時間復雜度為O(uvw).其中，對于MLT算法需要的空間為圖像像素大小，令每像素值映射存儲空間為c，其空間復雜度可以表示為S(uvc).本文算法兩個階段采樣數量相比MLT算法沒有改變，總采樣數量為u×v×w，則本文算法的時間復雜度為O(uvw).對于空間復雜度，本文算法在MLT算法基礎上增加了采樣矩陣的存儲，其存儲類型為整形，大小為生成圖像像素個數，令每像素存儲空間為r，其空間復雜度可以表示為S(uvc+uvr).

3 實驗結果與分析

3.1 實驗環境及度量標準

為了驗證改進效果，將本文算法應用于PSSMLT算法中，并與PSMLT(path space Metropolis light transport)、PSSMLT和MMLT進行比較，其中，PSMLT為原始MLT方法.由于本文算法的改進側重于提議樣本的突變策略，所以在重要性函數選取較為單一的PSSMLT算法上進行實驗改進.由于MMLT算法使用多重重要性采樣對樣本進行采樣，是目前較為優秀的基于重要性采樣的改進算法，因此，將這種算法也作為對比算法之一.使用3個經典場景進行渲染實驗，分別為具有直接光照的Cornell場景、具有間接光照的Room場景和共同具有直接光照和間接光照的Veach場景.實驗環境為：Intel Core I7-6700 CPU@3.4 GHz，內存4 GB，Win7專業版，編程環境為Visual Studio 2010.

對生成的圖像進行質量評價，采用兩種全參考圖像評價方法[10-12].第一種為峰值信噪比(peak signal to noise ratio，PSNR)，生成的結果中PSNR的值越大，則生成圖像質量越好.第二種為結構相似度(structural similarity，SSIM)，SSIM的范圍為[0，1]，且滿足最大值唯一性，當參考圖像和生成圖像完全相同時，SSIM=1.

3.2 采樣分布變化

由于本文算法對馬爾科夫鏈所形成采樣分布的影響，本文算法會使圖像的采樣分布發生變化.圖2為采樣密度比較.圖2中使用Veach場景，PSSMLT和本文算法突變概率為0.01，平均每像素采樣個數為128，二階段起始位置設置為96倍像素個數.參考圖像使用PSSMLT，每像素平均采樣個數為32 768，突變概率為0.01.

圖2 采樣密度比較Fig.2 Comparison in sampling density

通過采樣密度比較可以發現，相較于PSSMLT中高亮區域的高密度采樣，本文算法將采樣更多地分布于場景中，減少高亮度區域重復采樣的同時，將部分采樣分布于場景中的陰影部分，從而減少了亮度較低區域的噪聲.

3.3 實驗結果對比

場景渲染過程中的參數需要進行設置，同時也對一些其他參數設置了默認值，即為防止連續拒絕提議樣本，提議樣本最大拒絕次數設為512.在光線追蹤中，光線路徑的最大長度設為7.其中，基于PSSMLT的基礎場景，突變概率設為0.25；基于MMLT的場景，突變概率設為0.1.由于PSSMLT和MMLT兩種算法的執行效率不同，在渲染過程中對不同場景分別使用統一的時間對場景進行渲染，但相同時間中PSSMLT及其改進算法執行較慢，為了保證采樣的遍歷性，對不同場景中的每像素平均采樣個數分別進行設置，如表1所示.

表1 場景平均每像素樣本值Tab.1 The scene averaged samples per pixel

三種場景的參考圖像均使用每像素平均采樣個數32 768進行了10 h以上的渲染，從而保證了生成圖像的收斂性和遍歷性.為保證實驗準確性，不同場景的多種算法的渲染時長均設置為相同值，三種場景的渲染時間如表2中第二列所示.對于本文算法所增加的空間主要集中在采樣矩陣占用的空間中，其數據類型為整形，則三種場景所增加的空間大小如表2中第三列所示.

表2 場景渲染時間Tab.2 Rendering time of scene

實驗一中，對具有直接光照的Cornell場景進行渲染，如圖3所示，其整幅圖像的PSNR和SSIM數據對比如表3所示，其中分別使用PSMLT、PSSMLT、MMLT和本文算法進行比較.由實驗可知，在直接光照場景中，PSMLT算法的PSNR值優于本文算法，這是由于其突變策略可以較好地覆蓋場景中的路徑.而相比于PSMLT、PSSMLT和MMLT，本文算法在陰影部分有明顯的噪聲消除效果，且在亮度值較大的光源部分沒有受到采樣點減少的影響.但是在直接光照場景中，場景中的有效光線路徑分布較為平均，這也使得場景采樣過多導致PSNR值較低，而SSIM值4種方法則基本相同.

圖3 Cornell場景生成圖像對比Fig.3 Comparison in image of Cornell scene

表3 Cornell場景評價指標對比Tab.3 Comparison in evaluationindexes of Cornell scene

實驗二中，采用具有間接光照的Room場景進行渲染，如圖4所示，其整幅圖像的PSNR和SSIM數據對比如表4所示.

由實驗可知，在間接光照場景中，本文算法分別較PSMLT、PSSMLT和MMLT算法減少了壺身和陰影部分的噪點.相比于直接光照場景，間接光照場景中更難以找到場景中的有效路徑，體現在陰影和遮擋部分需要經過光線多次彈射而形成路徑的部分，PSMLT和PSSMLT在光線較少的陰影部分的提議樣本會被多次拒絕，從而導致大量噪聲.本文算法與MMLT的生成圖像質量相當，而相對于PSMLT和PSSMLT算法有明顯提升.

圖4 Room場景生成圖像對比Fig.4 Comparison in image of Room scene

表4 Room場景評價指標對比Tab.4 Comparison in evaluationindexes of Room scene

實驗三中，采用具有直接光照和間接光照的Veach場景進行渲染，如圖5所示，其整幅圖像的PSNR和SSIM數據對比如表5所示.

圖5 Veach場景生成圖像對比Fig.5 Comparison in image of Veach scene表5 Veach場景評價指標對比Tab.5 Comparison in evaluationindexes of Veach scene

方法PSNR/dBSSIMPSMLT27 5854500 638878PSSMLT28 1288700 587478MMLT30 7687800 789231本文算法29 4300910 653658

由實驗可知，在混合光照場景中，本文算法在PSMLT、PSSMLT算法基礎上減少了場景中的部分噪點，尤其表現在陰影部分和金屬材質部分.本文算法的PSNR值和SSIM值略低于MMLT算法，對于桌面反射投影至墻面的光線有部分噪聲，但與MMLT效果基本相同.

4 結 論

本文對MLT算法的渲染過程和突變策略進行了改進.本文算法一階段以PSSMLT作為基礎，通過一階段形成的采樣矩陣影響二階段的突變策略，并通過新的接收概率決定提議樣本的接收和拒絕，在圖像生成過程中接受更好的提議路徑，使采樣分布在場景中更平均，從而減小生成圖像噪聲.由上述實驗可知，在三個經典場景中，相比于PSMLT和PSSMLT，本文算法質量較好，相比于MMLT，本文算法對于間接光照場景有部分提升，有效降低圖像生成過程中產生的隨機噪聲，提高了生成圖像的質量.對于兩個階段起始位置選取判斷的標準，也是判斷馬爾科夫鏈收斂的標準，從而使兩個階段的渲染具有自適應性.將本文的突變策略與具有多重重要性采樣的MMLT算法相結合將會成為后續研究的重點.

[1] 牛連強，張丹，陶峰.直線的光柵轉換算法與快速反走樣繪制技術 [J].沈陽工業大學學報，2012，34(1)：73-78.

(NIU Lian-qiang，ZHANG Dan，TAO Feng.Raster-conversion algorithm and fast anti-aliased drawing technique for line [J].Journal of Shenyang University of Technology，2012，34(1)：73-78.)

[2] Veach E，Guibas L J.Metropolis light transport [C]//Proceedings of the 24th Annual Conference on Computer Graphics and Interactive Techniques.New York，USA，1997：65-76.

[3] Jakob W，Marschner S.Manifold exploration：a Markov chain Monte Carlo technique for rendering scenes with difficult specular transport [J].ACM Transactions on Graphics，2012，31(4)：1-13.

[4] Hoberock J，Hart J C.Arbitrary importance functions for Metropolis light transport [J].Computer Graphics Forum，2010，29(6)：1993-2003.

[5] Szirmay K L，Dornbach P，Purgathofer W.On the start-up bias problem of metropolis sampling [C]//Proceedings of Winter School of Computer Graphics’ 99.Plzen，eskrepublika，1999：273-280.

[6] Kelemen C，Szirmay K L，Antal G，et al.A simple and robust mutation strategy for the metropolis light transport algorithm [J].Computer Graphics Forum，2002，21(3)：531-540.

[7] Hachisuka T，Kaplanyan A S，Dachsbacher C.Multiplexed metropolis light transport [J].ACM Transactions on Graphics，2014，33(4)：1-10.

[8] Li T M，Lehtinen J，Ramamoorthi R，et al.Anisotropic Gaussian mutations for Metropolis light transport through Hessian-Hamiltonian dynamics [J].ACM Transactions on Graphics，2015，34(6)：1-13.

[9] Pharr M，Humphreys G.Physically based rendering：from theory to implementation [M].San Francisco：Morgan Kaufmann，2004：26-27.

[10]Chu J，Chen Q，Yang X.Review on full reference image quality assessment algorithms [J].Application Research of Computers，2014，31(1)：13-22.

[11]Lu L，Liu H，Huang X Y.Deeply research on diversity evaluation mode of “fundamentals of programming” [J].Springer London，2012，154：1595-1601.

[12]Keller A.Quasi-Monte Carlo image synthesis in a nutshell [C]//Proceedings in Mathematics & Statistics.Berlin，Germany，2013：213-249.

AnimprovedalgorithmforsamplingdistributionbasedonMLT

HE Huai-qing, CHEN Shuai, LIU Hao-han, JI Yu
(College of Computer Science and Technology, Civil Aviation University of China, Tianjin 300300, China)

In order to solve the random noise problem in the image generation process caused by the sparsity of random sampling distribution in Metropolis light tracking (MLT) algorithm, a mutation strategy based on sampling distribution was proposed. The sampling process was divided into two stages to improve the MLT algorithm, and the mutation strategy in the second stage was affected by the sampling distribution matrix generated from the first stage. In addition, a small scale of sampling was performed at the positions of eight-neighbor pixels of current sampling point. While satisfying the detailed balance condition, the average ratio of scalar contribution function of multiple sample points was set as the acceptance probability. The results show that the improved algorithm can generate smaller noise in the image compared with the original method at the same time, and the effect is equivalent to the better improved algorithms. The problem of big noise in the indirect illumination scene can be solved through improving the mutation strategy.

Metropolis light tracking; sampling distribution; global illumination; Markov chain Monte Carlo; photorealistic rendering; mutation strategy; acceptance probability; path selection

2017-02-23.

國家自然科學基金民航聯合研究基金資助項目(U1333110).

賀懷清(1969-)，女，吉林白山人，教授，博士，主要從事圖形圖像與可視化分析等方面的研究.

* 本文已于2017-06-21 21∶21在中國知網優先數字出版. 網絡出版地址： http：∥www.cnki.net/kcms/detail/21.1189.T.20170621.2121.022.html

10.7688/j.issn.1000-1646.2017.06.09

TP 391.41

A

1000-1646(2017)06-0646-08

(責任編輯：鐘 媛 英文審校：尹淑英)