基于最近離散點的光線跟蹤

蔡 鵬， 尹寶才， 孔德慧

（北京工業(yè)大學(xué)計算機(jī)學(xué)院多媒體與智能軟件技術(shù)北京市重點實驗室，北京 100124）

基于點云模型的光線跟蹤算法得到很多的關(guān)注，因為它在全局光照方面提供了高質(zhì)量的繪制效果。目前，點云模型的光線跟蹤算法的點元，除了包含三維坐標(biāo)，還包含法向量或半徑等空間信息[1-3]。但通過激光掃描等方法獲得的原始數(shù)據(jù)，即點云數(shù)據(jù)，只含有三維坐標(biāo)的空間信息，點云模型的離散點的法向量和半徑等屬性是預(yù)先計算出來的，這需要很多時間，特別就大規(guī)模的點云模型而言；而且對于曲率變化大或支離破碎的模型，計算的離散點的法向量極有可能不一致，即造成光照計算的錯誤（不一致）。本文提出基于最近離散點的光線跟蹤算法，直接對原始的點云模型進(jìn)行光線跟蹤，主要貢獻(xiàn)有以下3項：

1）對原始的點云模型進(jìn)行光線跟蹤，不需計算離散點的法向量和半徑等信息，因此，減少了內(nèi)存開銷。對于無法加載的大規(guī)模的點云模型（含法向量）或網(wǎng)格模型到內(nèi)存進(jìn)行繪制的情況，本文算法在內(nèi)存的容量限制下僅加載上述模型的坐標(biāo)信息，并在較少的時間內(nèi)，完成光線跟蹤，同時保證光線與離散點的局部表面交點的法向量的一致性（即使模型的曲率變化大或者是支離破碎的模型），從而產(chǎn)生光照一致的高質(zhì)量的圖像。

2）通過改變光線跟蹤的參數(shù)（最近離散點的數(shù)目），即可達(dá)到漸進(jìn)地多分辨率顯示原始的點云模型的目的。對于噪聲多的原始的點云模型，設(shè)置較大的最近離散點的數(shù)目，以有效地減少其繪制的噪聲；對于噪聲少的原始的點云模型，設(shè)置較小的最近離散點的數(shù)目，以更多地顯示其局部幾何特征。

3）通過平衡二叉樹在設(shè)定范圍內(nèi)搜索離光線迭代點最近的N個離散點，并用柵格的加速結(jié)構(gòu)避免不必要的迭代搜索計算，提高了光線跟蹤的時間效率。

1 以前的工作

目前，許多學(xué)者提出有效的方法來計算光線與點云表面的交點[4-8]，我們列出一些主要的點云模型的光線跟蹤的方法如下：

Schaufler和Jenson[4]將光線想象成具有一定半徑的的圓柱，通過計算圓柱前面的一些點的位置和法向量的平均值得到交點，此方法需要半徑更大的圓柱來避免空洞。這將需要更多的交點計算時間，并且降低了繪制質(zhì)量。

Adamson和Alexa[5]提出基于光線上的迭代點到多項式擬合的逼近表面投影的高度差的求交算法。這個方法能獲得非常高的繪制質(zhì)量，但它的計算成本很高。

Wand和Strasser[6]將光線考慮成各向異性的圓錐，通過光線的圓錐體與點云的多分辨率層次中預(yù)濾波的采樣點相交來計算交點，算法提供了沒有走樣的圖像，但計算時間很長。

Wald和Seidel[7]提出基于點云的隱式表面的光線跟蹤的框架。他們使用在光線迭代點周圍的Splats（具有圓心、半徑和法向量屬性的圓平面）的圓心和法向量的平均來定義局部平面逼近函數(shù)，計算光線上等間隔的點的函數(shù)值。通過在不同符號的函數(shù)值的兩個迭代點插值計算得到交點。

Linsen等人[8]在Splat（定義同上）產(chǎn)生時，計算Splat的法向場。當(dāng)光線與表面的求交時，通過Splat的法向場插值得到交點的法向量。此方法能得到光滑的顯示效果。

以上點云模型的光線跟蹤算法都增加了某些點元的屬性，如法向量、半徑、包圍球內(nèi)點云的擬合表面的多項式系數(shù)、點云的層次信息和Splat的法向場等，這將需要更多的內(nèi)存空間和預(yù)處理時間。本文直接對原始的點云模型進(jìn)行光線跟蹤，因此，減少了內(nèi)存開銷，并且計算時間較少。

Hart[9]提出光線迭代點到隱式表面的距離作為高度差和步長的基于球體范圍的光線跟蹤算法。Adamson[5]等提出基于點云模型的多項式逼近表面的高度差的求交算法。本文方法也提出高度差和步長的概念，與之不同的是，當(dāng)光線迭代點在離它實際最近的N′個離散點的包圍球內(nèi)時，本文定義迭代點到這些離散點的局部平面的垂直距離為高度差，并取其一半作為步長，而且具體迭代求交的過程和局部平面的計算也不同。

2 基于最近離散點的光線跟蹤

2.1 基于最近離散點的光線跟蹤算法的框架

本文首先給出基于最近離散點的光線跟蹤算法的基本思想如下：

1）計算原始的點云模型的平衡二叉樹的存儲結(jié)構(gòu)和柵格的空間信息。

2）光線迭代點逐步向點云模型的局部表面逼近，并使光線直接穿過沒有離散點的空柵格。

3）用平衡二叉樹在設(shè)定范圍內(nèi)搜索離光線迭代點最近的N個離散點，并計算實際離迭代點最近的N′個離散點的局部平面。若光線與局部平面滿足相交條件，得到交點；否則，計算下一個迭代點。

2.2 點云模型的局部平面的定義

當(dāng)光線上的迭代點Ai逼近點云局部表面時，以Ai為圓心、R′為半徑的包圍球包含N′個離散點，即點云模型中離Ai最近的N′個離散點。計算N′個離散點的包圍球的球心為離散點k的坐標(biāo)；并計算o到N′個離散點中最遠(yuǎn)點的距離作為包圍球的半徑r，如圖1所示。

圖1 離迭代點Ai最近的N′個離散點的局部表面(a)和對應(yīng)的局部平面(b)

設(shè)離迭代點Ai最近的N′(N′≥3)個離散點對應(yīng)的局部平面過包圍球的球心o，其法向量為n。⊿oPiPj是以球心o和N′個離散點中的兩個點Pi和Pj為頂點的三角形，其中j= (i+1) modN′，即⊿oPiPj的頂點Pi和Pj的順序是按N′個離散點在最大堆（搜索最近離散點時用的堆結(jié)構(gòu)）的結(jié)點順序。

其中，SΔoPiPj和nΔoPiPj分別為⊿oPiPj的面積和法向量。令noAi為o到迭代點Ai的法向量，即，使nΔoPijP的方向與noiA一致，即：

＜,＞為向量內(nèi)積運算，以下均同。SΔoPiPj由海倫公式求得，nΔoPiPj由⊿oPiPj兩條邊的單位向量的叉積求得的向量經(jīng)單位化得到。迭代點Ai沿光線方向逐步向點云模型的局部表面逼近（詳見下節(jié)），確保Ai總是在點云表面的一側(cè)，即法向量noAi的方向相對于點云表面是一致性的。SΔoPiPj越大，則nΔoPiPj對局部平面的法向量n的計算貢獻(xiàn)越大，因此，取SΔoPiPj為nΔoPiPj的權(quán)值。通過式（1）和式（2），計算的局部平面的法向量n相對于點云表面是一致性的，并且是連續(xù)的，從而保證了高質(zhì)量的繪制效果。

2.3 基于最近離散點的光線跟蹤算法

本文采用平衡二叉樹搜索最近的N個離散點（光子映射中多用平衡二叉樹搜索空間中某點的最近若干光子[10]）。平衡二叉樹中的每個結(jié)點對應(yīng)點云模型的一個離散點，按結(jié)點的垂直于某坐標(biāo)軸的分裂面的值（即離散點的對應(yīng)坐標(biāo)的值），相應(yīng)空間的離散點被分成左、右兩個子樹，即左子樹中離散點的對應(yīng)坐標(biāo)的值小于該分裂面的值，而右子樹中離散點則相反。這樣從根結(jié)點開始，搜索離迭代點Ai最近的N個離散點時，相對于分裂面，在每個結(jié)點可以排除一半不必要的搜索的點。因此，平衡二叉樹的搜索效率很高，同時，初始的搜索半徑R越小，則搜索效率越高[10]。

搜索算法使用最大堆存放離迭代點Ai實際最近的N′個離散點。最大堆采用平衡二叉樹的結(jié)構(gòu)，其中每個結(jié)點的關(guān)鍵值（Ai與相應(yīng)離散點的距離）比左、右兩子結(jié)點大。搜索算法如下：

//d初始為最大搜索半徑R，堆的結(jié)點個數(shù)為N；

// 調(diào)用LocatePoints(1)，即從根結(jié)點開始搜索離迭// 代點Ai最近的N個離散點；

// 返回離迭代點Ai實際最近的N′個離散點在堆中。

LocatePoints(p)

{

if (2p+1 ＜ 點云的離散點的總數(shù)目)

{

令δ為迭代點Ai到結(jié)點(離散點)p的分裂面的距離;

if (δ＜ 0) {

LocatePoints(2p); // 迭代點Ai在左子樹,

// 先遍歷左子樹;

if (δ2＜d2) // 檢查右子樹;

LocatePoints(2p+ 1); }

else {

LocatePoints(2p+ 1); // 迭代點Ai在右子樹,

// 先遍歷右子樹;

if (δ2＜d2) // 檢查左子樹;

LocatePoints(2p); }

}

δ2為迭代點Ai到結(jié)點(離散點)p的平方距離;

if (δ2＜d2) {將結(jié)點(離散點)p插入最大堆中;d2= 迭代點Ai到最大堆根結(jié)點的平方距離;}

}

對點云模型的包圍盒進(jìn)行柵格化，每個柵格對應(yīng)一個布爾變量，即當(dāng)柵格中有離散點時設(shè)為真，否則設(shè)為假。這樣可以減少算法在點云空間的迭代搜索的計算量，因為光線上的迭代點Ai直接跳過空柵格（沒有離散點的柵格），不用在空的柵格中搜索最近的離散點，故能提高效率。

圖2 光線上的迭代點Ai逐步逼近點云局部表面(a)，Ai不在(b)和在(c)最近N′個離散點的包圍球內(nèi)的情況

當(dāng)光線上的迭代點Ai在某一非空柵格（有離散點的柵格）中，對于點云模型的所有離散點，通過平衡二叉樹中在最大半徑R范圍內(nèi)搜索離Ai最近的N個離散點，得到實際最近的N′個離散點。若N′為零，以最大搜索半徑R為步長沿光線方向前進(jìn)，直至逼近點云局部表面，如圖2(a)所示。當(dāng)N′大于零時，即光線上的迭代點Ai接近點云局部表面，分為兩種情況。

第1種情況，如圖2(b)所示，迭代點Ai不在N′個離散點的包圍球（球心半徑r為o到N′個離散點中最遠(yuǎn)點的距離，Pk為離散點k的坐標(biāo)）內(nèi)，即光線不與此包圍球相交。因此，光線不與包圍球內(nèi)的點云局部表面相交，計算下一個迭代點Ai+1。

第2種情況，如圖2(c)所示，迭代點Ai在N′個離散點的包圍球內(nèi)，即光線與此包圍球相交，判斷是否滿足光線與點云局部表面相交的條件。若滿足相交的條件，得到相應(yīng)的交點；否則，計算下一個迭代點Ai+1。

光線與點云局部平面相交的情況如圖3所示，設(shè)光線與過球心o、法向量為n（n的計算見上節(jié)）的局部平面相交于B。若Ai與B的距離和半徑r的比值小于閾值α?xí)r，或Ai到局部平面的垂直距離與r的比值小于閾值β時，則光線與點云局部表面相交，即相交的判斷條件為：

圖3 光線與點云局部平面（過點o、法向量為n）相交于B點，h為Ai到局部平面的垂直高度

若光線與原始的點云模型相交，布爾變量intersect被設(shè)置為TRUE；否則，它被設(shè)置為FALSE。具體求交算法如下：

輸入：原始的點云模型、平衡二叉樹、柵格、最大的搜索半徑R、最近的離散點個數(shù)N、閾值α和β

輸出：布爾變量intersect，交點的坐標(biāo)和法向量

步驟1計算光線與點云包圍盒的交點，若有交點，即光線迭代點A0，轉(zhuǎn)到步驟2；否則，返回intersect= FALSE。

步驟2計算光線迭代點Ai所在的柵格，若柵格在點云包圍盒外，返回intersect= FALSE；否則，判斷當(dāng)前柵格是否有離散點，若無，計算光線與鄰近柵格的交點作為下一個迭代點Ai+1，轉(zhuǎn)到步驟2，若有，轉(zhuǎn)到步驟3。

步驟3對于點云模型的所有離散點，利用平衡二叉樹在最大搜索半徑R的范圍內(nèi)搜索離迭代點Ai最近的N個離散點，得到實際的N′個最近的離散點（見上面的LocatePoints(p)算法）。

步驟4若N′= 0，則Ai以步長R沿光線方向計算下一個迭代點Ai+1，轉(zhuǎn)到步驟2；否則，轉(zhuǎn)到步驟5。

步驟5若N′＜ 3，則計算Ai到N′個離散點中最近點的距離，取其一半作為步長，沿光線方向計算下一個迭代點Ai+1，轉(zhuǎn)到步驟2；否則，轉(zhuǎn)到步驟6。

步驟6若迭代點Ai不在N′個離散點的包圍球內(nèi)，則計算Ai到其中最近點的距離，取其一半作為步長，沿光線方向計算下一個迭代點Ai+1，轉(zhuǎn)到步驟2，如圖2(b)所示；否則，轉(zhuǎn)到步驟7，如圖2(c)所示。

步驟7由式（1）和式（2）計算過包圍球圓心的局部平面的法向量n，并計算光線與局部平面的交點B（見圖3）。若滿足條件式（3），則得到相應(yīng)的交點即為B，其法向量即為n，返回intersect= TRUE；否則，以迭代點Ai到局部平面的垂直距離h的一半作為步長，沿光線方向計算下一個迭代點Ai+1，轉(zhuǎn)到步驟2。

上述算法通過改變光線跟蹤的參數(shù)（最近離散點的數(shù)目），即可漸進(jìn)地多分辨率顯示原始的點云模型。對于噪聲多的原始的點模型，設(shè)置較大的最近離散點的數(shù)目，以有效地減少其繪制的噪聲；對于噪聲少的原始的點模型，設(shè)置較小的最近離散點的數(shù)目，以更多地顯示其局部幾何特征。

3 實驗設(shè)計和結(jié)果

本文所有實驗都在單個的Intel (R) 2.66GHz CPU平臺上執(zhí)行的。所有的點云模型都是從www.dirdim.com網(wǎng)站下載的免費的點云數(shù)據(jù)文件，即只包含三維坐標(biāo)信息。點云模型的離散點的數(shù)目如表1所示。

表1 4個點云模型的離散點的數(shù)目

實驗的設(shè)置參數(shù)：柵格數(shù)目NGd、離散點的最大搜索半徑R、最近的離散點個數(shù)N、光線迭代點與局部平面的光線交點的距離和包圍球半徑比值的閾值α、迭代點到局部平面的垂直距離與包圍球半徑比值的閾值β（詳見2.3節(jié)）。設(shè)L為點云模型的包圍盒的X、Y和Z軸長度和。

實驗的結(jié)果數(shù)據(jù)：實際最近離散點的平均非零個數(shù)N′、平衡二叉樹的建立時間TSetTree、柵格的建立時間TSetGd、點云模型的光線跟蹤時間TRay、計算時間的總和TTotal(即TTotal=TSetTree+TSetGd+TRay)。

本文對點云模型Demo Head和Facility Scanning的實驗參數(shù)N分別設(shè)置不同的值，得到相應(yīng)的實驗結(jié)果。所有點云模型的實驗參數(shù)的設(shè)置如表2所示，其光線跟蹤實驗的結(jié)果統(tǒng)計如表3所示。

表2 點云模型的光線跟蹤實驗的設(shè)置參數(shù)

在圖4-6中，原始的點云模型的光線跟蹤生成圖像的分辨率都為400×500像素。通過設(shè)置迭代點的最近離散點的不同的個數(shù)N，即當(dāng)N逐漸變大（即4、8、12、16），點云模型Demo Head的繪制效果變得更光滑，同時，點云表面的一些細(xì)節(jié)部分變得模糊，相應(yīng)的光線跟蹤時間也增多（見表3）。如圖5所示，當(dāng)N逐漸變大（即8、16、24、32），點云模型Facility Scanning的繪制的噪聲逐漸減少。由圖4和圖5可見，對于噪聲少的模型（如Demo Head），可以設(shè)置較小的最近離散點個數(shù)，以更多地保留其局部幾何特征；對于噪聲多的模型（如 Facility Scanning），設(shè)置較大的最近離散點個數(shù)，以有效地減少其噪聲。

圖4 原始的點云模型的光線跟蹤的繪制效果（400×500像素）

圖6 原始的點云模型的光線跟蹤的繪制效果（400×500像素）

表3 點云模型的光線跟蹤實驗的結(jié)果統(tǒng)計（時間單位：毫秒）

光線迭代點的實際最近離散點的個數(shù)N′略小于設(shè)置的值N，這是由于迭代點逐漸逼近點云表面并且最大搜索半徑R的范圍的限制（詳見2.3節(jié)）。點云模型Facility Scanning的離散點的數(shù)目遠(yuǎn)大于其它的點云模型，因此，相應(yīng)的計算時間也更多（見表3）。由圖4-6可見，對于原始的點云模型，在較少的時間內(nèi)（見表3），本文的算法能產(chǎn)生高質(zhì)量的繪制效果。

4 結(jié) 論

本文提出了基于最近離散點的光線跟蹤算法，即在較少的時間內(nèi)，直接對僅包含三維坐標(biāo)的點云模型進(jìn)行光線跟蹤，不需計算離散點的法向量和半徑等信息。算法保證了光線與局部平面的交點的法向量計算的一致性，從而產(chǎn)生十分滿意的繪制效果。

本文算法通過改變光線跟蹤的參數(shù)（最近離散點的數(shù)目），即可達(dá)到漸進(jìn)地多分辨率顯示原始的點云模型的目的。對于噪聲多的原始的點模型，設(shè)置較大的最近離散點的數(shù)目，以有效地減少其繪制的噪聲；對于噪聲少的原始的點模型，設(shè)置較小的最近離散點的數(shù)目，以更多地顯示其局部幾何特征。

今后，計劃將本文的算法擴(kuò)展到GPU上，達(dá)到實時交互的目的。

[1]Levoy M, Whitted T. The use of points as display primitives [J]. Technical Report 85-022, Computer Science Department, University of North Carolina at Chapel Hill, 1985.

[2]Kobbelt L, Botsch M. A survey of point-based techniques in computer graphics [J]. Computers &Graphics, 2004, 28(6): 801-814.

[3]Sainz M, Pajarola R. Point-based rendering techniques [J].Computers & Graphics, 2004, 28: 869-879.

[4]Schaufler G, Jensen H W. Ray tracing point sampled geometry [C]//Proceedings of the 11th Eurographics Workshop on Rendering, 2000: 319-328.

[5]Adamson A, Alexa M. Ray tracing point set surfaces [C]//Proceedings of the Shape Modeling International,2003: 272-279.

[6]Wand M, Strasser W. Multi-resolution point-sample raytracing [C]//Graphics Interface, 2003: 139-148.

[7]Wald I, Seidel H P. Interactive ray tracing of point-based models [C]//Eurographics Symposium on Point-Based Graphics, 2005: 1-8.

[8]Linsen L, M¨uller K, Rosenthal P. Splat-based ray tracing of point clouds [J]. Journal of WSCG, 2007,15(2): 51-58.

[9]Hart J C. Sphere tracing: a geometric method for the antialiased ray tracing of implicit surfaces [J]. The Visual Computer, 1996, 12(9): 527-545.

[10]Jensen H W. Realistic image synthesis using photon mapping [M]. A.K.Peters; 2001.