張璐鵬　賈世宇

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

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

中圖分類號： TP391.41文獻標識碼： A

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

1相關(guān)工作

1.1距離場

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

A.Fuhrmann等人[5]提出通過分析三角網(wǎng)格法線方向產(chǎn)生的棱柱來快速計算局部距離場的方法；Rong G等人[6]提出在常數(shù)時間內(nèi)計算輸入網(wǎng)格近似距離場的jump flooding算法；S.P.Mauch[7]根據(jù)掃描轉(zhuǎn)換解決靜態(tài)HamiltonJacobi方程，進而計算精確距離場的CSC算法等。上述算法全部建立在封閉流形幾何的基礎(chǔ)上，無法計算非流形幾何體的距離場。由于沒有明確定義劃分非流形幾何體的內(nèi)部、外部和邊界，所以在計算距離場時，通常通過從三角網(wǎng)格多個方向進行投影，生成二進制體數(shù)據(jù)[8]，如果三角網(wǎng)格包含孔洞或自交，將會體現(xiàn)在掃描轉(zhuǎn)換后的結(jié)果中。雖然可以通過對每個體素的掃描轉(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)學者們提出了大量的體素化算法，J.Huang等人[10]提出的保證拓撲一致性的體素化方法；Fang S等人[11]利用硬件加速方法實現(xiàn)體素化算法，但這些算法通常需要大量內(nèi)存，并且輸出一個非縮放的三維體素網(wǎng)格；M.Schwarz等人[12]提出一種基于八叉樹的稀疏固體體素化算法，雖然內(nèi)存要求大大降低，但是該方法仍受限于可用內(nèi)存。

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

2八叉樹體模型生成算法

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

2.1計算距離場

距離場是一個標量場，表示場內(nèi)每個點到給定對象上最近點的距離，該距離有符號，表示該點位于模型的內(nèi)部或外部。定義為

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

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

2.1.1提取無符號距離場

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

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

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

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

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

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

2.1.2計算符號距離場

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

引理1如果x為偏移曲面Eσ外的一點，則符號距離場dS定義為

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

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

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

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

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

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

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

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

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

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

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

2.3體素化

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

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

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

3實驗結(jié)果與分析

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

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

4結(jié)束語

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

參考文獻：

[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]. 圖學學報， 2005， 26（4）： 17.