基于Metropolis光線跟蹤的組合濾波器

吳熙　徐慶　卜紅娟　王征

摘要：

蒙特卡羅方法是計(jì)算全局光照的基礎(chǔ)，目前已經(jīng)有很多基于蒙特卡羅的全局光照算法，但大多數(shù)算法在渲染時(shí)間上都有一定局限性。在蒙特卡羅方法基礎(chǔ)上，結(jié)合Metropolis光線跟蹤算法和組合濾波器，提出一種新的全局光照算法。該算法分為兩個(gè)部分，首先使用多組不同尺度的濾波器對(duì)圖像進(jìn)行處理，然后將多組濾波器處理后的結(jié)果組合成最終的結(jié)果。該算法使用相對(duì)均方根誤差作為選擇濾波尺度的依據(jù)，在采樣和重建過(guò)程中自適應(yīng)地為每個(gè)像素選擇合適的濾波器，以最大化降低誤差，得到更好的重建結(jié)果。實(shí)驗(yàn)結(jié)果表明，該算法相對(duì)于傳統(tǒng)Metropolis算法在效率和圖像質(zhì)量上都有較大提高。

關(guān)鍵詞：

Metropolis算法；蒙特卡羅方法；光線跟蹤；組合濾波器；全局光照算法

中圖分類號(hào)：

TP391.41

文獻(xiàn)標(biāo)志碼：A

Abstract：

The Monte Carlo method is the basis of calculating global illumination. Many Monte Carlobased global illumination algorithms have been proposed. However， most of them have some limitations in terms of rendering time. Based on the Monte Carlo method， a new global illumination algorithm was proposed， combining the Metropolis ray tracing algorithm with an integrated filter. The algorithm is composed of two parts. In the first part， multiple sets of filters with different scales were used to smooth the image； in the second part， filtered images were combined into the final result. Relative Mean Squared Error （RMSE） was used as a basis for the selection of filtering scale， and an appropriate filter was adaptively selected for each pixel during the process of sampling and reconstruction， aiming to reduce the errors to a minimum degree and gain better reconstruction results. Experimental results show that the proposed method outperforms many traditional Metropolis algorithms in terms of both efficiency and image quality.

英文關(guān)鍵詞Key words：

Metropolis algorithm； Monte Carlo method； ray tracing； integrated filter； global illumination algorithm

0引言

目前Monte Claro渲染熱門(mén)研究方向大致分為兩個(gè)。第一個(gè)方向是通過(guò)尋找更加有效的光徑來(lái)提高圖像合成質(zhì)量。Kelemen等[1]提出一種新的突變策略來(lái)提高M(jìn)etropolis效率，將路徑突變用一個(gè)空間里面的隨機(jī)數(shù)來(lái)代替，能很好地提高M(jìn)etropolis算法效率。Lehtinen等[2]提出一種基于圖像梯度的Metropolis算法，直接計(jì)算圖像梯度，通過(guò)求解Poisson方程直接從圖像梯度重建圖像。Hachisuka等[3]通過(guò)將多重重要性采樣和Metropolis算法相結(jié)合，提出一種多路復(fù)用Metropolis算法，用視點(diǎn)路徑的長(zhǎng)短來(lái)計(jì)算多重重要性的權(quán)值，增加路徑利用率，提高合成圖像效率。

第二個(gè)方向是對(duì)圖像進(jìn)行自適應(yīng)渲染濾波。Rousselle等[4]提出一種自適應(yīng)的NonLocal Means算法，使用雙緩存的方法，在合成圖像過(guò)程中生成兩個(gè)采樣數(shù)目相同的緩存，并分別計(jì)算濾波系數(shù)，交叉處理，去除噪聲對(duì)真實(shí)像素的影響。Kaplanyan等[5]提出一種局部自適應(yīng)的Photon Mapping算法，該算法通過(guò)平衡渲染過(guò)程中的噪聲和偏差來(lái)使渲染結(jié)果錯(cuò)誤率最小。Delbracio等[6]提出一種基于光線顏色直方圖來(lái)促進(jìn)圖像渲染的方法，該方法通過(guò)統(tǒng)計(jì)每個(gè)像素采樣點(diǎn)的顏色分布來(lái)判斷像素之間是否共享光線。Liu等[7]提出通過(guò)梯度來(lái)對(duì)圖像進(jìn)行聚類分析，然后通過(guò)多項(xiàng)式擬合圖像的方式進(jìn)行自適應(yīng)渲染。Zwicker等[8]對(duì)自適應(yīng)渲染方法作了一個(gè)詳盡的概述。

為進(jìn)一步提高蒙特卡羅渲染效率，本文將這兩個(gè)方向相結(jié)合，提出一種新的蒙特卡羅渲染方法。該算法相比傳統(tǒng)渲染算法，不是通過(guò)暴力增加采樣點(diǎn)來(lái)提高圖像質(zhì)量，而是通過(guò)尋找更有效光徑和更適宜的濾波尺度來(lái)提高渲染質(zhì)量。

1Metropolis原理

假設(shè)定義一個(gè)狀態(tài)空間Ω和一個(gè)非負(fù)函數(shù)f：Ω→R+，Metropolis算法通過(guò)f產(chǎn)生一系列采樣點(diǎn)xi。當(dāng)選定起始點(diǎn)x0，下一個(gè)采樣點(diǎn)x1將通過(guò)x0突變產(chǎn)生，之后每一個(gè)采樣點(diǎn)xi都是由它的前一個(gè)采樣點(diǎn)xi-1通過(guò)突變產(chǎn)生。這種每個(gè)采樣點(diǎn)（也稱為狀態(tài)）xi只與它的前一個(gè)采樣點(diǎn)xi-1有關(guān)的隨機(jī)過(guò)程叫作馬爾可夫鏈。每次突變都可能被接收或者是被拒絕，但無(wú)

論從哪一個(gè)采樣點(diǎn)x0開(kāi)始，最后得到的采樣點(diǎn)xi都是服從和函數(shù)f成比例的分布，這個(gè)分布就叫作平穩(wěn)分布。

為得到服從平穩(wěn)分布的采樣點(diǎn)，需要使用特定方式來(lái)接受或拒絕突變。Metropolis算法采用細(xì)節(jié)平衡的方法來(lái)保證馬爾可夫過(guò)程能夠趨向于平穩(wěn)分布。

1.1細(xì)節(jié)平衡

假設(shè)當(dāng)前時(shí)刻狀態(tài)為xi，通過(guò)式（1）可以得到狀態(tài)xi+1。在生成狀態(tài)xi+1時(shí)，首先根據(jù)預(yù)先制定的突變策略得到臨時(shí)狀態(tài)x′。這里突變函數(shù)定義為M（|），該函數(shù)表示從狀態(tài)轉(zhuǎn)移到狀態(tài)的概率，函數(shù)M（|）被稱為臨時(shí)轉(zhuǎn)移函數(shù)。在狀態(tài)轉(zhuǎn)移過(guò)程中，每次突變得到的臨時(shí)狀態(tài)x′都可能被接收或拒絕，將取決于接收概率r（|）。生成新?tīng)顟B(tài)表達(dá)式：

xi+1=x′，概率為r（x′|xi）

xi，概率為1-r（x′|xi） （1）

如何定義接收概率r（|）是馬爾可夫鏈快速達(dá)到平穩(wěn)分布的關(guān)鍵。假設(shè)已經(jīng)達(dá)到平穩(wěn)分布，那么必須定義轉(zhuǎn)移概率T，使已經(jīng)達(dá)到的平穩(wěn)分布保持?？紤]狀態(tài)空間中兩個(gè)狀態(tài)和，從宏觀上來(lái)看，它們達(dá)到平穩(wěn)分布。然而在微觀上，它們一直進(jìn)行著相互轉(zhuǎn)移，狀態(tài)之間的相互轉(zhuǎn)移一直沒(méi)有停止，它們之間處于一個(gè)動(dòng)態(tài)平衡，即從狀態(tài)到狀態(tài)的采樣點(diǎn)和從狀態(tài)到狀態(tài)的采樣點(diǎn)一定相等，這就稱之為細(xì)節(jié)平衡[9]，即：

f（）M（|）r（|）=f（）M（|）r（|） （2）

可以證明當(dāng)滿足這個(gè)條件時(shí)，平穩(wěn)分布可以得到保持。通過(guò)式（2）可以得出計(jì)算接收概率r（|）的方法，由于f（）已知，而M（|）由任意選取的突變策略決定，所以可以計(jì)算出r（|）/r（|）。為盡快達(dá)到平穩(wěn)分布，需要一個(gè)較大的接收概率，因此接收概率r（|）定義為：

r（|）=min1，f（）M（|）f（）M（|）（3）

這樣，就能夠使得r（|）和r（|）中較大的被置為1，也就表示兩個(gè)狀態(tài)之間的轉(zhuǎn)移，一定會(huì)有一個(gè)方向被接收，而另一個(gè)方向的轉(zhuǎn)移有一定概率會(huì)被接收。這樣設(shè)定接收概率，既可以保證細(xì)節(jié)平衡也可以保證馬爾可夫過(guò)程以較快的速度達(dá)到平穩(wěn)分布。

1.2Metropolis渲染

在圖像渲染過(guò)程中，渲染方程[10-11]描述能量從光源到達(dá)人眼的過(guò)程，對(duì)它求解是圖像合成的關(guān)鍵問(wèn)題。在物體表面上某點(diǎn)x沿方向α出射的輻射亮度的定義如式（4）所示：

L（x→α）=Le（x→α）+

∫Ωxfr（x， β α）L（x ← β）cos β dωβ（4）

這里Ωx是點(diǎn)x所有半球方向的集合，Le（x→α）是點(diǎn)x自身所發(fā)出的輻射亮度， fr是雙向散射分布函數(shù)（Bidirectional Scattering Distribution Function， BSDF），L（x←β）是點(diǎn)x從β方向得到的入射輻射度。

式（4）是渲染方程最基本形式，由它可得到多種變形，比如路徑積分形式[12]等，它們最終目的都是方便求取出射輻射亮度。Metropolis光線傳輸[9，11]將式（4）改寫(xiě)成蒙特卡羅積分形式：

Ij≈CN∑ihj（xi）L*（xi）L（xi）（5）

其中：Ij是像素j處的輻射亮度，C是使用蒙特卡羅方法求得L*的積分值，N是通過(guò)馬爾可夫鏈產(chǎn)生的隨機(jī)路徑個(gè)數(shù)，hj是濾波權(quán)值系數(shù)，L*是圖像的貢獻(xiàn)函數(shù)，馬爾可夫鏈生成的路徑就是和L*成比例的隨機(jī)數(shù)。

2組合濾波器

本文提出的方法中，組合濾波器指在渲染過(guò)程中，使用多個(gè)尺度不同但類型相同的濾波器對(duì)圖像進(jìn)行處理，處理后會(huì)得到多組圖像濾波結(jié)果，然后根據(jù)圖像之間的特征來(lái)選取不同尺度的濾波器進(jìn)行組合得到最后結(jié)果。因?yàn)樵诿商乜_渲染中，圖像不同區(qū)域噪聲水平不同，例如場(chǎng)景邊界、陰影等比較復(fù)雜的區(qū)域渲染時(shí)會(huì)有較多噪聲，而一些平滑區(qū)域則噪聲較少。所以不同區(qū)域需要不同尺度的濾波器以取得最好效果。使用組合濾波器可以在噪聲水平不同的區(qū)域使用不同尺度濾波器。在這個(gè)方法中，有三個(gè)關(guān)鍵問(wèn)題需要考慮：首先是濾波器種類的選取，其次為濾波器之間尺度差異的選取，最后是組合結(jié)果時(shí)依據(jù)的標(biāo)準(zhǔn)。本文算法基本流程如圖1所示。

2.1組合濾波

本節(jié)將詳細(xì)介紹組合濾波器原理。濾波器的選擇上，本文選擇較為簡(jiǎn)單的高斯濾波器，因?yàn)楦咚篂V波器能夠快速地對(duì)圖像進(jìn)行濾波，在渲染過(guò)程中能夠滿足快速合成圖像的需求。本文提出的算法選取m組高斯濾波器，定義它們尺度為ai（i=1，2，…，m），它們滿足a1

組合濾波器關(guān)鍵點(diǎn)是如何選取組合圖像濾波結(jié)果的依據(jù)。圖像渲染過(guò)程中，渲染后處理在很多情況下都選取方差作為濾波器的輸入信息，Rousselle等[13]使用均方根誤差（Mean Squared Error， MSE）作為選擇濾波器的依據(jù)，Li等[14]選擇MSE的一個(gè)無(wú)偏估計(jì)量（Steins Unbiased Risk Estimator， SURE）作為選擇雙邊濾波器的依據(jù)。本文提出的方法中使用相對(duì)均方根誤差（Relative Mean Squared Error， RMSE）作為選擇濾波器的依據(jù)。

在渲染過(guò)程中，因?yàn)闊o(wú)法直接知道渲染的真實(shí)結(jié)果，也就無(wú)法直接計(jì)算MSE，所以本文使用RMSE選擇濾波器。在渲染過(guò)程中，可以估算偏差和方差，并以此估算RMSE。在計(jì)算RMSE時(shí)有三個(gè)關(guān)鍵點(diǎn)需要注意。首先，濾波尺度從小到大對(duì)應(yīng)著圖像的結(jié)果是越來(lái)越平滑；其次，對(duì)于大多數(shù)像素，隨著濾波尺度從小到大變化其偏差呈單調(diào)遞增趨勢(shì)，而方差則是呈單調(diào)遞減趨勢(shì)；最后，偏差和方差之間需要一個(gè)權(quán)衡系數(shù)使二者達(dá)到一個(gè)相對(duì)平衡。RMSE計(jì)算公式如下：

RMSE（ai→ai+1）=MSE[ai+1]-MSE[ai]=

ΔVar+ΔBias=ξ*（Var[ai+1]-Var[ai]）+Bias[ai+1]2-Bias[ai]2 （6）

其中：RMSE表示相對(duì)均方根誤差，ai是濾波器的尺度，MSE指在各自尺度下的均方根誤差。本文所定義的RMSE是在不知道圖像真實(shí)值得情況下估算的結(jié)果。為了能夠消去真實(shí)值，將RMSE定義為兩個(gè)不同尺度的MSE之差。RMSE定義分為方差和偏差兩部分，本文提出參數(shù)ξ來(lái)均衡方差和偏差之間的差距。

2.2參數(shù)計(jì)算

在未知圖像真實(shí)值得情況下估算圖像偏差是本文算法很重要的模塊。本文用ai和ai+1分別表示小尺度和大尺度的濾波器，pi和pi+1分別表示使用小尺度和大尺度濾波器濾波后的像素值。文獻(xiàn)[15]表明濾波尺度和偏差之間的關(guān)系：

偏差可以通過(guò)式（10）求得，那么求取RMSE就只需解出ΔVar部分即可。對(duì)于各個(gè)尺度下的方差求取，本文參考Rousselle等[13]中的做法。本文和Rousselle方法的不同之處在于計(jì)算RMSE所需的ΔVar需要在兩個(gè)相鄰尺度下的方差之差上乘上權(quán)衡系數(shù)ξ，系數(shù)ξ的選取直接影響濾波效果。因?yàn)镽MSE是由偏差和方差兩部分組成，并且當(dāng)濾波尺度由小變大時(shí)，由于方差和偏差變化的單調(diào)性，通常情況下會(huì)有ΔVar<0，ΔBias>0，所以RMSE的值就會(huì)在0附近變化。因?yàn)榉讲詈推钔?jì)算出來(lái)后不在同一個(gè)數(shù)量級(jí)，所以需要一個(gè)系數(shù)ξ來(lái)平衡它們。通過(guò)大量實(shí)驗(yàn)得出，方差和偏差之間相差的數(shù)量級(jí)約為10000，因此ξ取值為10000。為了使該算法適應(yīng)不同的場(chǎng)景，在不同的場(chǎng)景中都能得到最佳的效果，針對(duì)不同的場(chǎng)景還需要對(duì)ξ進(jìn)行細(xì)微調(diào)整。為實(shí)現(xiàn)算法方便引入一個(gè)新的系數(shù)n用于微調(diào)系數(shù)ξ，在具體實(shí)現(xiàn)時(shí)，直接將兩個(gè)系數(shù)相乘即可。通過(guò)實(shí)驗(yàn)分析，n的取值在1～10的大多數(shù)場(chǎng)景都能取得較好的渲染效果。

本文通過(guò)RMSE來(lái)選擇所需要的濾波器尺度。ΔVar和ΔBias分別表示同一像素較大尺度和較小尺度之間的方差和偏差變化。通常在尺度較小時(shí)方差變化會(huì)大于偏差變化，所以RMSE為負(fù)數(shù)。本文選擇的濾波器是使RMSE絕對(duì)值最小時(shí)的濾波器，這樣能夠很好地平衡方差和偏差的變化得到最優(yōu)的結(jié)果。在實(shí)際操作中，隨著濾波器的尺度由小變大，當(dāng)RMSE>0時(shí)，就選擇較小尺度的濾波器，否則就繼續(xù)遍歷。

3實(shí)驗(yàn)結(jié)果與分析

本文算法基于PBRTV2系統(tǒng)[16-17]實(shí)現(xiàn)，所有結(jié)果都是在Intel Core i52400 CPU@ 3.10GHz，內(nèi)存4GB上實(shí)現(xiàn)。本算法易實(shí)現(xiàn)，本文選用8組高斯濾波器進(jìn)行組合。為選擇合適的系數(shù)ξ，進(jìn)行了大量實(shí)驗(yàn)。由實(shí)驗(yàn)得出，大多數(shù)圖像首先需要一個(gè)改變方差數(shù)量級(jí)的系數(shù)10000，才能夠細(xì)微調(diào)整方差和偏差之間的平衡系數(shù)n，經(jīng)過(guò)實(shí)驗(yàn)這個(gè)系數(shù)可以選擇1～10。

圖2是平衡系數(shù)n為5時(shí)得到的Toasters場(chǎng)景每個(gè)像素分布16個(gè)采樣點(diǎn)的尺度選擇圖像。圖中由黑到白分別對(duì)應(yīng)著濾波尺度1～8。由圖中可以看出，在三角面中心區(qū)域，圖像比較平滑，該區(qū)域圖像的像素值較為相似。所以選擇的濾波尺度往往比較大，這樣可以使圖像更加平滑，去除渲染所帶來(lái)的噪聲；而在三角面邊緣區(qū)域圖像像素值變化較大，所以選擇的濾波尺度往往比較小，這樣有利于保持圖像的特征，使圖像不至于變得過(guò)于模糊。

圖3中給出cornellmlt和Toasters兩個(gè)場(chǎng)景的渲染結(jié)果。

圖3（a）是cornellmlt場(chǎng)景渲染結(jié)果，圖3（b）是toasters場(chǎng)景渲染結(jié)果。Lehtinen等[2]和Hachisuka等[3]所提出的對(duì)Metropolis方法的改進(jìn)側(cè)重點(diǎn)和本文有所不同，所以僅選擇將本文提出的算法和經(jīng)典Metropolis方法[1]以及只用單個(gè)的高斯濾波進(jìn)行對(duì)比。cornellmlt場(chǎng)景渲染生成的圖像大小是250×250，平衡系數(shù)n=4，每個(gè)像素分布128個(gè)采樣點(diǎn)，首先計(jì)算場(chǎng)景中的直接光照，采用雙向路徑追蹤渲染圖像。合成圖像大概需要30s，合成速度很快，在可以接受范圍內(nèi)。Toasters場(chǎng)景中，圖片的原始大小是400×400，每個(gè)像素分布16個(gè)采樣點(diǎn)，合成圖像所用時(shí)間大約需要10s。圖3中，左邊是原始Metropolis所合成的圖像，中間是采用單個(gè)高斯濾波器所合成的圖像，右邊是本文算法合成圖像結(jié)果?？偟膩?lái)說(shuō)，本文算法能夠根據(jù)圖像特征來(lái)自適應(yīng)選擇濾波器尺度，提高渲染效率。

從圖3的結(jié)果可以直觀看出，本文提出的組合濾波器算法渲染結(jié)果優(yōu)于經(jīng)典的Metropolis算法和單個(gè)的高斯濾波器算法。為更加進(jìn)一步證明本文的方法確實(shí)有效，又對(duì)PBRT官方網(wǎng)站[17]上的多個(gè)經(jīng)典場(chǎng)景進(jìn)行了實(shí)驗(yàn)。以每個(gè)像素2048個(gè)采樣點(diǎn)所生成的圖像作為參考圖像，計(jì)算圖像的峰值信噪比（Peak SignaltoNoise Ratio， PSNR）和MSE，各個(gè)場(chǎng)景的渲染結(jié)果計(jì)算的結(jié)果如表1所示。通過(guò)表1中MSE和PSNR的比較，可以看出本文的組合濾波器相比原始的Metropolis算法和單一的高斯濾波效果有明顯的提高。

4結(jié)語(yǔ)

本文基于蒙特卡羅方法提出一種新的渲染方法，該算法將組合濾波器和Metropolis算法相結(jié)合，大大提高了蒙特卡羅渲染效率，在較短時(shí)間內(nèi)便可以得到高質(zhì)量的真實(shí)感圖像。后續(xù)研究中，我們的主要工作是利用更少的采樣點(diǎn)，尋找更有效路徑來(lái)得到高質(zhì)量的合成圖像；同時(shí)，將進(jìn)一步改進(jìn)本算法，尋找自適應(yīng)的方法來(lái)設(shè)置該算法中需要的平衡系數(shù)，并且嘗試將本算法應(yīng)用到動(dòng)畫(huà)的合成過(guò)程中，以提高其渲染的效率。

參考文獻(xiàn)：

[1]

KELEMEN C， SZIRMAYKALOS 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.

[2]

LEHTINEN J， KARRAS T， LAINE S， et al. Gradientdomain metropolis light transport [J]. ACM Transactions on Graphics， 2013， 32（4）： Article No. 95.

[3]

HACHISUKA T， KAPLANYAN A S， DACHSBACHER C. Multiplexed metropolis light transport [J]. ACM Transactions on Graphics， 2014， 33（4）： Article No. 100.

[4]

ROUSSELLE F， KNAUS C， ZWICHER M. Adaptive rendering with nonlocal means filtering [J]. ACM Transactions on Graphics， 2012， 31（6）： Article No. 195.

[5]

KAPLANYAN A S， DACHSBACHER C. Adaptive progressive photon mapping [J]. ACM Transactions on Graphics， 2013， 32（2）： Article No. 16.

[6]

DELBRACIO M， MUS P， BUADES A， et al. Boosting Monte Carlo rendering by ray histogram fusion [J]. ACM Transactions on Graphics， 2014， 33（1）： Article No. 8.

[7]

LIU X D， ZHENG C W. Adaptive cluster rendering via regression analysis [J]. Visual Computer， 2015， 31（1）： 105-114.

[8]

ZWICKER M， JAROSZ W， LEHTINEN J， et al. Recent advances in adaptive sampling and reconstruction for Monte Carlo rendering [J]. Computer Graphics Forum， 2015， 34（2）： 667-681.

[9]

鮑世強(qiáng).基于雙向路徑追蹤的 Metropolis 全局光照算法的研究[D].天津：天津大學(xué)，2006. （BAO S Q. The study of bidirectional path tracing based metropolis global illumination method [D]. Tianjin： Tianjin University， 2006.）

[10]

KAIJYA J T. The rendering equation [J]. ACM SIGGRAPH Computer Graphics， 1986， 20（4）： 143-150.

[11]

VEACH E， GUIBAS L J. Metropolis light transport [C] //Proceedings of the 24th Annual Conference on Computer Graphics and Interactive Techniques. New York： ACM Press/AddisonWesley Publishing Company， 1997： 65-76.

[12]

PHARR M， HUMPHREYS G. Physically Based Rendering： From Theory to Implementation [M]. San Francisco， CA： Morgan Kaufmann， 2004： 652-660.

[13]

ROUSSELLE F， KNAUS C， ZWICKER M. Adaptive sampling and reconstruction using greedy error minimization [J]. ACM Transactions on Graphics， 2011， 30（6）： Article No. 159.

[14]

LI T M， WU Y T， CHUANG Y Y. SUREbased optimization for adaptive sampling and reconstruction [J]. ACM Transactions on Graphics， 2012， 31（6）： Article No. 194.

[15]

SILVERMAN B W. Density Estimation for Statistics and Data Analysis [M]. Boca Raton， FL： CRC press， 1986：1-22.

[16]

Implementation of pbrtv2 [S/OL]. [20151205]. https：//github.com/mmp/pbrtv2/wiki.

[17]

Scenes and data for rendering [S/OL]. [20151205]. http：//www.pbrt.org.