適用于柔性體切割仿真的八叉樹體模型生成算法

張璐鵬　賈世宇

摘要： 針對(duì)任何三角幾何模型的體素化需求，本文基于距離場(chǎng)和八叉樹結(jié)構(gòu)的特性，提出將任意三角幾何模型離散為體素模型的算法。提取任意三角幾何模型的無符號(hào)距離場(chǎng)的偏移表面，自動(dòng)刪除內(nèi)部表面，然后從偏移表面計(jì)算符號(hào)距離場(chǎng)，并調(diào)整閾值，得到與原模型接近的重構(gòu)表面，最后對(duì)重構(gòu)表面進(jìn)行體素化，得到適用于切割仿真的體素模型，為驗(yàn)證算法的有效性，本文對(duì)基于空間分割的體素化算法和本文提出的算法進(jìn)行對(duì)比測(cè)試。測(cè)試結(jié)果表明，該算法對(duì)于包含自相交和非連通的輸入三角幾何的體素化效果更好，更為貼近原模型，且時(shí)間和內(nèi)存成本較低，效果較好。該算法為仿真柔性體切割的體素模型提供一種有效方案。

關(guān)鍵詞： 距離場(chǎng)； 八叉樹； 非流形； 自相交； 三角幾何； 體素化

中圖分類號(hào)： TP391.41文獻(xiàn)標(biāo)識(shí)碼： A

近年來，隨著可編程圖形處理器和計(jì)算機(jī)圖形學(xué)的發(fā)展，醫(yī)學(xué)仿真在醫(yī)學(xué)領(lǐng)域的應(yīng)用越來越廣泛。醫(yī)學(xué)仿真中的手術(shù)仿真通過精確的人體組織器官模型，真實(shí)模擬組織器官在手術(shù)器械的外力交互作用下變形乃至被切割等操作過程，而為實(shí)現(xiàn)人體器官柔性體變形和切割的真實(shí)感，模型需要被體素化后才能夠?qū)崿F(xiàn)相應(yīng)的數(shù)學(xué)計(jì)算進(jìn)行仿真。目前，用于生成體素模型的方法有很多。Huang J等人[1]提出利用模型法線到表面函數(shù)作為距離標(biāo)準(zhǔn)，將模型離散為體素模型的方法；S.Oomes等人[2]采用尺度空間理論（Scalespace Theory）轉(zhuǎn)換模型為反鋸齒體素?cái)?shù)據(jù)；S.Laine[3]提出基于拓?fù)浣Y(jié)構(gòu)的體素化方法。以上研究適用于絕大部分幾何模型體素化的方法，雖然這些方法精度和速度較高，但體素化應(yīng)用對(duì)象只用于流形幾何，普適性較差。針對(duì)人體器官模型無定向、非流形或包含自相交的三角幾何網(wǎng)格的特點(diǎn)，本文提出基于距離場(chǎng)和表面重構(gòu)算法填補(bǔ)模型表面孔洞，界定模型在數(shù)學(xué)意義上的內(nèi)外之分，然后對(duì)模型進(jìn)行體素化處理，從而產(chǎn)生適用于切割仿真的體素模型。

1相關(guān)工作

1.1距離場(chǎng)

距離場(chǎng)的概念由Rosenfeld和Pfaltz在1966年與數(shù)字圖像處理相關(guān)的論文中提出[4]。距離場(chǎng)存在有無符號(hào)區(qū)分，符號(hào)指的是某個(gè)點(diǎn)位于模型內(nèi)部還是外部。只有當(dāng)幾何模型為密封的流形幾何體時(shí)，才能計(jì)算出其符號(hào)距離場(chǎng)，而當(dāng)幾何模型為非流形時(shí)，便無法明確符號(hào)的定義，也就無法從數(shù)學(xué)意義上界定模型的內(nèi)外。

A.Fuhrmann等人[5]提出通過分析三角網(wǎng)格法線方向產(chǎn)生的棱柱來快速計(jì)算局部距離場(chǎng)的方法；Rong G等人[6]提出在常數(shù)時(shí)間內(nèi)計(jì)算輸入網(wǎng)格近似距離場(chǎng)的jump flooding算法；S.P.Mauch[7]根據(jù)掃描轉(zhuǎn)換解決靜態(tài)HamiltonJacobi方程，進(jìn)而計(jì)算精確距離場(chǎng)的CSC算法等。上述算法全部建立在封閉流形幾何的基礎(chǔ)上，無法計(jì)算非流形幾何體的距離場(chǎng)。由于沒有明確定義劃分非流形幾何體的內(nèi)部、外部和邊界，所以在計(jì)算距離場(chǎng)時(shí)，通常通過從三角網(wǎng)格多個(gè)方向進(jìn)行投影，生成二進(jìn)制體數(shù)據(jù)[8]，如果三角網(wǎng)格包含孔洞或自交，將會(huì)體現(xiàn)在掃描轉(zhuǎn)換后的結(jié)果中。雖然可以通過對(duì)每個(gè)體素的掃描轉(zhuǎn)換，并使用奇偶規(guī)則判定重構(gòu)合理的體數(shù)據(jù)[9]，但仍然存在幾何體結(jié)構(gòu)化的異常值處理困難的問題。

1.2體素化

1987年，A.Kaufman[10]首次提出體素化概念，根據(jù)三維模型的網(wǎng)格數(shù)據(jù)，轉(zhuǎn)換為包含模型表面信息和內(nèi)部屬性的體數(shù)據(jù)集。近年來，相關(guān)學(xué)者們提出了大量的體素化算法，J.Huang等人[10]提出的保證拓?fù)湟恢滦缘捏w素化方法；Fang S等人[11]利用硬件加速方法實(shí)現(xiàn)體素化算法，但這些算法通常需要大量?jī)?nèi)存，并且輸出一個(gè)非縮放的三維體素網(wǎng)格；M.Schwarz等人[12]提出一種基于八叉樹的稀疏固體體素化算法，雖然內(nèi)存要求大大降低，但是該方法仍受限于可用內(nèi)存。

隨著可編程圖形處理器的快速發(fā)展，研究者們利用計(jì)算機(jī)圖形處理器（graphics processing unit，GPU）加速改進(jìn)體素化方法。C.Crassin等人[13]使用GPU光柵單元進(jìn)行稀疏體素化；M.Ptzold等人[14]提出一種非網(wǎng)格化（gridfree）、非核心（outofcore）的GPU體素化方法。這些方法雖然利用GPU并行特性進(jìn)行算法改進(jìn)，但是仍然存在算法判斷過程復(fù)雜，內(nèi)存消耗較大，對(duì)硬件要求較高的缺陷。

2八叉樹體模型生成算法

針對(duì)包含自相交、孔洞、無定向、不連通及重復(fù)的幾何體在內(nèi)的任何三角幾何模型的體素化需求，本文提出一種基于八叉樹的體模型生成算法。該方法分為三個(gè)階段：一是計(jì)算得到輸入三角網(wǎng)格的符號(hào)距離場(chǎng)，并重構(gòu)模型表面網(wǎng)格；二是利用保證拓?fù)浣Y(jié)構(gòu)的Marching Cubes算法修復(fù)重構(gòu)的表面網(wǎng)格，得到封閉的流形表面網(wǎng)格；三是利用八叉樹對(duì)模型空間進(jìn)行劃分和掃描，區(qū)分模型內(nèi)外部數(shù)據(jù)，進(jìn)而獲取體數(shù)據(jù)，實(shí)現(xiàn)模型體素化。

2.1計(jì)算距離場(chǎng)

距離場(chǎng)是一個(gè)標(biāo)量場(chǎng)，表示場(chǎng)內(nèi)每個(gè)點(diǎn)到給定對(duì)象上最近點(diǎn)的距離，該距離有符號(hào)，表示該點(diǎn)位于模型的內(nèi)部或外部。定義為

dist（x，Ω）=inf‖p-x‖， （x∈R3，p∈Ω）（1）

式中，dist（x，Ω）表示從點(diǎn)x到模型Ω的無符號(hào)距離場(chǎng)；p是屬于模型Ω表面網(wǎng)格上的點(diǎn)；inf‖p-x‖表示點(diǎn)p到x的距離下界；R3表示距離場(chǎng)內(nèi)的點(diǎn)的集合。符號(hào)距離場(chǎng)是在原有無符號(hào)距離場(chǎng)的基礎(chǔ)上，使用正負(fù)號(hào)表示點(diǎn)x位于三角幾何模型Ω的內(nèi)部或外部。通常情況下，判定曲面內(nèi)部點(diǎn)的符號(hào)為負(fù)，曲面外部點(diǎn)的符號(hào)為正，若用σ表示模型Ω的內(nèi)部距離場(chǎng)的點(diǎn)集，則符號(hào)距離場(chǎng)φx的定義為

2.1.1提取無符號(hào)距離場(chǎng)

計(jì)算無符號(hào)距離場(chǎng)前進(jìn)行空間劃分[16]，首先將輸入三角網(wǎng)格用平行于坐標(biāo)軸的包圍盒包含，隨后根據(jù)距離場(chǎng)偏移距離σ擴(kuò)展已有包圍盒，從而生成采樣空間。根據(jù)所定義的劃分分辨率h，將整個(gè)采樣空間均勻劃分，模型包圍盒用固定值均勻劃分。均勻空間劃分如圖1所示。

依次計(jì)算每個(gè)采樣點(diǎn)到包圍盒立方體的最小距離和最大距離，取最大距離中的最小值作為采樣點(diǎn)到三角網(wǎng)格的最小距離的參考值，隨后遍歷所有包圍盒立方體，比較參考值和采樣點(diǎn)到包圍盒立方體的最小距離。如果當(dāng)前立方體的最小距離大于參考值，則忽略當(dāng)前立方體；若小于參考值，則計(jì)算采樣點(diǎn)到當(dāng)前立方體中所有三角面的距離，并更新最小值，最終得到無符號(hào)距離場(chǎng)。以無符號(hào)距離場(chǎng)作為輸入，使用Marching Cubes算法產(chǎn)生一個(gè)流形的偏移曲面，偏移曲面定義為

Sσ={x∈R3|dU（x，Ω）=σ}（3）

式中，dU（x，Ω）表示從采樣點(diǎn)x到輸入幾何Ω的無符號(hào)距離場(chǎng)。生成的偏移曲面可能存在嵌套，其定義為

Sσ=（∏iEiσ）∏（∏iIiσ）（4）

式中，Eiσ代表外部組件；Iiσ表示內(nèi)部組件；∏表示不相交的集合。使用UnionFind算法檢測(cè)所有連通分支，去除所有內(nèi)部組件，得到外部偏移曲面。

2.1.2計(jì)算符號(hào)距離場(chǎng)

通過計(jì)算得到一個(gè)流形的偏移曲面Eσ和輸入幾何的無符號(hào)距離場(chǎng)，通過使用偽法線測(cè)試[15]計(jì)算偏移曲面Eσ的符號(hào)距離場(chǎng)。由于蠻力執(zhí)行偽法線測(cè)試運(yùn)行效率較低，所以本文利用文獻(xiàn)[16]證明所得的兩條引理。

引理1如果x為偏移曲面Eσ外的一點(diǎn)，則符號(hào)距離場(chǎng)dS定義為

dS（x，Eσ）=dU（x，Ω）-σ≥0（5）

引理2如果點(diǎn)x滿足dU（x，Ω）<σ，則點(diǎn)x一定位于偏移曲面Eσ的內(nèi)部。

引理1表明在計(jì)算偏移曲面Eσ的外部符號(hào)距離場(chǎng)時(shí)，再次使用之前所計(jì)算的符號(hào)距離場(chǎng)，極大程度地縮減計(jì)算量。但是對(duì)于內(nèi)部而言，當(dāng)點(diǎn)x在內(nèi)部無限接近于偏移曲面Eσ時(shí)，引理1的等式不成立，所以之前計(jì)算的無符號(hào)距離場(chǎng)不能應(yīng)用于內(nèi)部符號(hào)距離場(chǎng)的計(jì)算。針對(duì)這一問題，由引理2可知，網(wǎng)格上所有滿足dU（x，Ω）<σ的點(diǎn)均無需進(jìn)行偽法線測(cè)試。

首先對(duì)偏移曲面Eσ進(jìn)行均勻的空間分割及邊界檢測(cè) ，標(biāo)記所有和偏移曲面相交的立方體，然后再次計(jì)算偏移曲面的無符號(hào)距離場(chǎng)，并執(zhí)行S遍歷，S遍歷后可得到曲面網(wǎng)格點(diǎn)的符號(hào)，最后在計(jì)算出Eσ的符號(hào)距離場(chǎng)基礎(chǔ)上，將其偏移-σ產(chǎn)生相對(duì)于輸入三角網(wǎng)格的符號(hào)距離場(chǎng)。

2.2修復(fù)重構(gòu)表面

為實(shí)現(xiàn)任意幾何體的體素化，得到輸入三角網(wǎng)格的符號(hào)距離場(chǎng)后，執(zhí)行Marching Cubes算法可以得到輸入幾何體的重構(gòu)表面Fσ。當(dāng)輸入幾何趨于封閉流形的非流形幾何時(shí)，重構(gòu)表面Fσ是一個(gè)封閉的流形幾何，可以直接進(jìn)行體素化。但是，并非所有的非流形幾何體都趨近于流形，當(dāng)一個(gè)非流形幾何由許多不連通的三角網(wǎng)格組成時(shí)，初次重構(gòu)的表面效果不足以體素化。偏移曲面圖如圖2所示。

為解決上述問題，本文在已有重構(gòu)表面的基礎(chǔ)上進(jìn)行修復(fù)。已知重構(gòu)表面Fσ是以符號(hào)距離場(chǎng)作為等值面執(zhí)行Marching Cubes所得，當(dāng)輸入三角網(wǎng)格本身為不連通的情況時(shí)，由于計(jì)算所得為精確符號(hào)距離場(chǎng)，所以距離場(chǎng)也是不連通的，導(dǎo)致重構(gòu)表面效果較差。此時(shí)，精確的符號(hào)距離場(chǎng)定義為

dU（x，Ω）=σ， σ→0（6）

當(dāng)符號(hào)距離場(chǎng)能夠連通并覆蓋整個(gè)輸入幾何體時(shí)，Marching Cubes的重構(gòu)表面才能達(dá)到足以進(jìn)行體素化效果，所以通過調(diào)整偏移量，讓符號(hào)距離場(chǎng)覆蓋整個(gè)輸入幾何體，即

dU（x，Ω）=σ+δ（7）

式中，δ作為調(diào)整的偏移量由用戶輸入。調(diào)整距離場(chǎng)后，再次利用Marching Cubes生成封閉的流形重構(gòu)表面。

2.3體素化

現(xiàn)將封閉的流形重構(gòu)表面體素化。針對(duì)重構(gòu)表面再次構(gòu)成包圍盒，并在頂點(diǎn)（xmin，ymin，zmin）固定與歐式空間方向相同的坐標(biāo)系，此時(shí)包圍盒空間定義為體素空間V3，每個(gè)體素由一組唯一整數(shù)序列（i，j，k）確定，且體素空間內(nèi)也存在“6鄰接”、“18鄰接”和“26鄰接”三種拓?fù)潢P(guān)系[16]。體素空間確定后，基于八叉樹原理進(jìn)行空間劃分[17]，一般八叉樹劃分方法僅對(duì)灰節(jié)點(diǎn)進(jìn)行劃分，本文則對(duì)3種節(jié)點(diǎn)均進(jìn)行劃分操作，生成均勻的體素網(wǎng)格，且體素標(biāo)志位滿足如下關(guān)系

f（i，j，k）=-1， 體素位于模型邊界外0， 體素與模型邊界相交1， 體素位于模型邊界內(nèi)（8）

通過空間細(xì)分獲取模型的表面數(shù)據(jù)，隨后利用逐行掃描根據(jù)模型表面數(shù)據(jù)快速區(qū)分模型內(nèi)外部數(shù)據(jù)，從而獲得體數(shù)據(jù)集，實(shí)現(xiàn)模型體素化。

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

為驗(yàn)證算法的有效性，本文對(duì)基于空間分割的體素化算法和本文提出的算法進(jìn)行對(duì)比測(cè)試?；诳臻g分割的體素化算法實(shí)質(zhì)是基于MIN提出的使用空間分割快速計(jì)算距離場(chǎng)算法的基礎(chǔ)上重構(gòu)表面并體素化的方法，由于此方法對(duì)趨于流形幾何的非流形幾何體的體素化效果較好，因此用于作對(duì)比測(cè)試。斯坦福兔子模型、骷髏模型和棕櫚植物模型的無符號(hào)距離場(chǎng)和符號(hào)距離場(chǎng)的測(cè)試結(jié)果如表1所示。

所有符號(hào)距離場(chǎng)均使用八核計(jì)算，3種模型均包含自相交、孔洞和非連通特性的三角幾何。體素化效果對(duì)比如圖3所示。本算法合理解決了具有精細(xì)和超薄特性的非流形幾何體的符號(hào)距離場(chǎng)計(jì)算問題。圖3b和圖3c展示斯坦福兔子模型分別使用基于空間分割的體素化算法和本文算法的體素化效果圖，可見兩種算法在對(duì)趨近于封閉流形的幾何模型進(jìn)行體素化時(shí)，所產(chǎn)生的效果并無較大差別；圖3e～i展示骷髏模型和棕櫚樹模型的體素化效果對(duì)比，經(jīng)過表面修復(fù)后，本文算法對(duì)于包含自相交和非連通的輸入三角幾何的體素化效果更好，體素化后的模型更為貼近原模型，而基于空間分割的體素化算法的體素化模型則明顯存在失真問題，無法精確還原模型。

4結(jié)束語

基于距離場(chǎng)和八叉樹結(jié)構(gòu)的特性，本文提出將任意三角幾何模型離散為體素模型的算法，解決了非流形幾何難以體素化的問題，采用表面修復(fù)，改善了精確距離場(chǎng)對(duì)重構(gòu)“三角形湯”幾何偏移曲面不連通的缺陷。實(shí)驗(yàn)結(jié)果表明，對(duì)非流形幾何進(jìn)行體素化時(shí)，相對(duì)基于空間分割的體素化算法，本算法能夠?qū)崿F(xiàn)更好的效果，考慮到三角幾何結(jié)構(gòu)的不確定性問題，可以適用于任意三角幾何的體素化，對(duì)通用體素化算法的研究具有一定參考價(jià)值。本文采用Marching Cubes算法進(jìn)行表面重構(gòu)，雖然保證了拓?fù)浣Y(jié)構(gòu)的一致性，但仍存在等值面精度較低的問題，導(dǎo)致修復(fù)效果并不完美。下一步研究可致力于改善算法效率和提升算法精度，利用GPU加速運(yùn)算、改善Marching Cubes算法或者選用更為高效的修復(fù)算法，使體素化的效果更為完美。

參考文獻(xiàn)：

[1]Huang J， Yagel R， Kurzion Y. An Accurate Method To Voxelize Polygonal Meshes[C]∥IEEE Symposium on Volume Visualization. Chapel Hill， North Carolina， USA： IEEE， 1998： 119126.

[2]Oomes S， Snoeren P， Dijkstra T. 3D Shape Representation： Transforming Polygons into Voxels[J]. International Conference on Scalespace Theory in Computer Vision， 1997， 1252： 349352.

[3]Laine S. A Topological Approach to Voxelization[J]. Computer Graphics Forum， 2013， 32（4）： 7786.

[4]Rosenfeld A， Pfaltz J. Sequential Operations in Digital Picture Processing[J]. Journal of the ACM， 1966， 13（4）： 471494.

[5]Fuhrmann A， Sobotka G. Distance Fields for Rapid Collision Detection in Physically Based Modeling[C]∥Proceedings of Graphi Con 2003. Moscow， Russia： IEEE， 2003： 5865.

[6]Rong G， Tan T S. Variants of Jump Flooding Algorithm for Computing Discrete Voronoi Diagrams[C]∥4th International Symposium on. Voronoi Diagrams in Science and Engineering. Clamorgan， UK： IEEE， 2007： 176181.

[7]Mauch S P. Efficient Algorithms for Solving Static HamiltonJacobi Equations[M]. USA： Mauch Sean Patrick， 2003.

[8]Spillmann J， Wagner M， Teschner M. Robust Tetrahedral Meshing of Triangle Soups[C]∥Vision， Modeling， Visualization （VMV）. Oxford， UK： IEEE， 2006： 916.

[9]Kaufman A， Shimony E. 3D ScanConversion Algorithms for VoxelBased Graphics[C]∥Proceedings of the 1986 Workshop on Interactive 3D graphics. Chapel Hill， North Carolina， USA： IEEE， 1987： 4575.

[10]Huang J， Yagel R， Filippov V， et al. An Accurate Method for Voxelizing Polygon Meshes[C]∥IEEE Symposium on Volume Visualization. Chapel Hill， North Carolina， USA： IEEE， 1998： 119126.

[11]Fang S， Chen H. Hardware Accelerated Voxelization[J]. Computers & Graphics， 2000， 24（3）： 433442.

[12]Schwarz M， Seidel H P. Fast Parallel Surface and Solid Voxelization on GPUs[J]. Acm Transactions on Graphics， 2010， 29（6）： 110.

[13]Crassin C， Neyret F， Sainz M， et al. Interactive Indirect Illumination Using Voxel Cone Tracing[J]. Symposium on Interactive 3d Graphics & Game， 2011， 30（7）： 19211930.

[14]Ptzold M， Kolb A. GridFree OutofCore Voxelization to Sparse Voxel Octrees on GPU[C]∥Proceedings of the 7th Conference on HighPerformance Graphics. Dublin， Ireland： ACM， 2015： 95103.

[15]Nooruddin F S， Turk G. Simplification and Repair of Polygonal Models Using Volumetric Techniques[J]. IEEE Transactions on Visualization and Computer Graphics， 2003， 9（2）： 191205.

[16]Xu H， Barbi J. Signed Distance Fields for Polygon Soup Meshes[C]∥Proceedings of Graphics Interface 2014. Processing Society， Montreal， Q C， Canada： ACM. 2014： 3541.

[17]Min D， Jia S， Zhang X. Multithreaded Fast Distance Field Computation Using Spatial Partition[J]. Science & Technology Vision， 2015（8）： 153154.

[18]Brentzen J A， Aanaes H. Signed Distance Computation Using the Angle Weighted Pseudonormal[J]. IEEE Transactions on Visualization and Computer Graphics， 2005， 11（3）： 243253.

[19]Meagher D. Geometric Modeling Using Octree Encoding[J]. Computer Graphics and Image Processing， 1982， 19（2）： 129147.

[20]吳曉軍， 劉偉軍， 王天然. 基于八叉樹的三維網(wǎng)格模型體素化方法[J]. 圖學(xué)學(xué)報(bào)， 2005， 26（4）： 17.