基于快速k近鄰的光子映射算法研究

王海波

摘要：光輻射強度估算是光子映射算法一個關(guān)鍵技術(shù)，傳統(tǒng)使用簡單、有效的k近鄰（kNN）算法，但kNN具有計算復雜度高，內(nèi)存需求量的缺點，新算法針對kNN的缺點，改進kNN搜索光子的方式，先將空間分割為多個固定長度的立方體，每個立方體體包含一定數(shù)量的光子數(shù)，通過測試各個立方體與光線接觸點之間的位置搜索接觸點周圍的k個最近鄰光子，進而估算光輻射強度，實驗表明新算法搜索光子的速度更快，而且圖形清晰度更高。

關(guān)鍵詞：光子映射;光輻射強度估算;kNN;空間網(wǎng)格

中圖分類號：TP3? ?文獻標識碼：A

文章編號：1009-3044（2019）28-0286-03

光子映射是較好的真實感圖形渲染方式之一，與光線跟蹤方法比較，光子映射能較好地渲染輝映、焦散效果[1]。光子映射分為兩個階段第一個階段發(fā)射光子，跟蹤光子，建立光子圖;第二階段利用光子圖估算光照，從圖形像素的角度發(fā)出光線，如果遇到反射或折射后，記錄接觸點，搜索接觸點周圍的光子，對接觸點進行光輻射強度估算，渲染圖形。因此光輻射強度估算是光子映射算法的第二階段的一個關(guān)鍵技術(shù)。

為了準確估算光輻射強度，需要尋找到對接觸點光輻射強度有影響的光子，傳統(tǒng)的光子映射算法采用K-近鄰算法對光輻射強度進行估算，k-近鄰算法具有簡單、有效的優(yōu)點，但發(fā)射光子量較大時k-近鄰算法的計算時間比較久。

為了改進計算精度，提高計算效率Garcia采用常核函數(shù)方法，并且放棄第K個光子輻射通量的計算，以避免偏差[2][3];PerChristensen主張對第K個光子計算1/2的光輻射通量[4];Reinhard Klein主張根據(jù)第k-1與第k個光子的平均距離來計算第k個光子對接觸點的影響[5];對于核函數(shù)可選擇圓錐核函數(shù)、高斯核函數(shù)、Epanechnikov核函數(shù)、Silverman核函數(shù)[6][7][8][9]。為實現(xiàn)光輻射強度估算無無偏差性，Hachisuka等提出漸進式光子映射算法，可從理論上對光輻射強度的計算達到無偏差計算，但需要循環(huán)發(fā)射大量光子[10][11]。

kNN算法具有簡單、有效的特點，因此如何利用kNN算法的優(yōu)點，提高光輻射強度估算也是新算法研究的方向之一。本算法先將空間分割為多個固定長度的立方體，每個立方體包含一定數(shù)量的光子數(shù)，通過測試各個球體與接觸點之間的位置速搜索測試點周圍的k個最近鄰光子，進而估算光輻射強度。

1 光子映射

光子映射是一種全局光照，分兩階段。第一階段從場景的光源發(fā)出大量的光子，跟蹤光子建立光子圖，建立密度低的全局光子圖及密度高的焦散光子圖。焦散光子圖要求光子至少被鏡面反射或折射一次，全局光子圖存儲各種路徑的光子，當光子擊中漫反射表面或略微光滑表面時，光子信息都會被存儲到相應(yīng)的光子圖中。光子圖存儲每個光子，并記錄光子能量或通量、輸入方向和位置等信息，光子圖采用平衡二叉樹作為光子搜索的數(shù)據(jù)結(jié)構(gòu)。

第二階段是利用光子圖計算光輻射強度，即從圖形的每個像素發(fā)出光線，光線在場景中經(jīng)過若干次鏡面反射、折射后、漫反射后，記錄接觸點x既估算點，然后利用全局光子圖及焦散光子圖，搜索離估算點最近的若干光子如圖1所示，并計算估算點的光輻射強度如式（1）所示

其中，x是估算點，fr是雙向反射分布函數(shù)（BRDF），[Φi]是第i個光子的光通量，[rk]是x與 k個光子中距x最遠的某個光子的距離。當光子數(shù)目越多，則光輻射強度的計算精度越高。因此光子的發(fā)射量影響到光輻射強度的計算，如何快速查找光子成為人們研究的方向。

2 快速K近鄰模型

K近鄰（KNN）算法是非常有效的基于距離的算法，被廣泛應(yīng)用于回歸模型，主要思想是：待測樣本的特征屬性等于最近鄰的k個樣本特征屬性的平均值。假設(shè)訓練集合A包含s個樣本，每個樣本又含有t個屬性既：{Ai=（Ai1，Ai2，…，Ait），i=1，2，…，s };待測樣本為B，屬性為（B1，B2，…，Bt）;則求得樣本B的特征屬的步驟如下：

3 快速kNN

3.1 構(gòu)建新算法

新算法使用空間網(wǎng)格的方法先將空間劃分為若干立方體，將所有光子置入到這些立方體中。立方體網(wǎng)格單元太大太小，都會對整個查找過程產(chǎn)生不良影響，若立方體網(wǎng)格單元太小，會增加存儲量，降低效率;若立方體網(wǎng)格單元太大，則每個立方體單元會包含過多面片，對求交造成困難，因此新算法劃分立方體的個數(shù)為p=n/t其中n是光子總數(shù)，t是未知數(shù)，一般取值為20，立方體的邊長為L，滿足[L3=n/t]。如果有多個不包含光子的立方體連續(xù)在一起，則合并為一個立方體，保證生成的空間單元數(shù)不超過O（n），如算法1所示。接著搜索估算點周圍的k個立方體，并從k個立方體中搜索最近的k個光子，完成對估算點光輻射強度的估算如算法2，圖2所示。

3.2 算法分析

設(shè)發(fā)射光子數(shù)n，k為kNN算法參數(shù)，m為估算點的個數(shù)。傳統(tǒng)KNN的時間復雜度為O（mnk），表示每個估算點要計算同n個光子的距離，同時為求出k個最近鄰的光子，內(nèi)存中要維持一個k長度的插入排序表，在排序表中，每插入一個新值，對表的最大操作次數(shù)為k。新算法中設(shè)立方體的總數(shù)為p，時間復雜性包含求k個最小立方體及其所包含樣本中的k個最近鄰樣本。每個立方體平均包含的光子數(shù)為n/p，因此時間復雜性為O（mpk+mnk2/p）=O（mk（p+nk/p）），故當p+nk/pp2/（p-k）時，新算法的時間復雜性低于傳統(tǒng)kNN算法，由于k取值同立方體數(shù)相比要小得多，故當n>p>k時，實際取值p<=n/20，因此n>p>k條件是成立的，新算法的運行效率優(yōu)于傳統(tǒng)的kNN算法。

4 算法實現(xiàn)

采用vs2013和OpenGL的編程環(huán)境，在一臺配置為Intel（R） Core（TM） i5，8GB內(nèi)存，NVIDIA GeForce 610M顯卡，win7下進行實驗。

實驗所用的測試場景如圖1，圖2，圖3 所示，表1給出了各圖的參數(shù)及渲染時間，其中圖1的發(fā)射的光子數(shù)為10M，圖2發(fā)射的光子數(shù)為15M，圖3發(fā)射的光子數(shù)為12M。由表中可以看出：新算法與傳統(tǒng)kNN的算法發(fā)射的光子數(shù)及需搜索的光子數(shù)是一樣的，新算法的每個立方體包含的光子數(shù)為20，從表中可以看出新算法的渲染時間比傳統(tǒng)kNN算法要快，圖1快24.7%，圖2快33.0%，圖3快23.1%，圖2增快的幅度最大是因為空間分割最緊密。從圖1、圖2、圖3的圖形渲染質(zhì)量方面來看新算法也比傳統(tǒng)的kNN算法更清晰。

5 總結(jié)

新算法針對kNN搜索光子的計算量大的缺點，改進了kNN搜索光子的方式，先將空間分割為多個固定長度的立方體，每個立方體體包含一定數(shù)量的光子數(shù)，通過測試各個球體與接觸點之間的位置速搜索測試點周圍的k個最近鄰光子，進而估算光輻射強度，實驗表明新算法搜索光子的速度更快，而且圖形清晰度更高。

參考文獻：

[1] Per H Christensen， Henrik W Jensen， Toshi Kato， Frank Suykens. A Practical Guide to Global Illumination Using Photon Mapping [R].USA： Siggraph， 2002

[2]Garcia， R.， Urena， C.， Sbert， M.. Description and solution of an unreported intrinsic bias in photon mapping density estimation with constant kernel[J]. Comput. Graph. Forum， 2012，31（1）：33-41.

[3]Garcia， R.， Urena， C.， Poch， J.， et al.， 2014. Overestimation and underestimation biases in photon mapping with non-constant kernels[J]. IEEE Trans. Visual. Comput. Graph.， 20（10）：1441-1450.

[4]Jensen H W， Christensen P. High quality rendering using ray tracing and photon mapping[C]// Siggraph 07： Acm Siggraph Courses， Acm. 2007.

[5] R. Garc??a， C. Ure ? na， and M. Sbert， “Description and solution of an unreported intrinsic bias in photon mapping density estimation

with constant kernel.” Presented at Eurographics 24th Symposium on Rendering. Invited paper from Computer GraphicsForum， 2013.

[6] Roland， S. Bias compensation for photon maps Comput[J]. Graph. Forum， 2003，22（4）：729-742.

[7]Myszkowski， K.， 1997. Lighting reconstruction using fast and adaptive density estimation techniques. Proc. Eurographics Workshop on Rendering Techniques， p.251-262.

[8] Schjoth， L.， 2009. Anisotropic Density Estimation in Global Illumination. PhD Thesis， University of Copenhagen， Denmark

[9]Schjoth， L.， Sporring， J.， Olsen， O.F.， 2008. Diffusion based photon mapping.

[10] T. Hachisuka， S. Ogaki， and H. W. Jensen， “Progressivephoton mapping，” ACM Trans. Graph.， vol. 27， pp. 130：1–130：8， December 2008.

[11] C. Knaus and M. Zwicker， “Progressive photon mapping： A probabilistic approach，” ACM Trans. Graph.， vol. 30，no. 3， pp. 25：1–25：13， May 2011.

【通聯(lián)編輯：梁書】