剖面數據場中的三維表面建模

摘 要：在三維表面建模技術中，Marching Cubes算法是應用最為廣泛的方法之一。該算法簡單高效，但與此同時，研究人員也發現它存在一些不足。在構造等值面時，Marching Cubes算法要把所有體素全部檢測一遍，即使有些體素沒有和等值面相交，這影響了算法效率；此外在這個過程中，Marching Cubes算法還會忽略掉一些本來在等值面上的點，降低了表面重建的精度。針對這些問題，本文對算法進行了改進。在構造等值面時，不檢測空的體素以提高算法的速度，并且把一些被忽略的等值點添加進來以提高算法的精度。

關鍵詞：三維表面建模；Marching Cubes算法；算法效率；逼近精度

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

Abstract：Marching Cubes algorithm is one of the most widely used methods in 3D surface modeling technologies.Although the algorithm is simple and efficient，researchers have found some shortcomings in it.When constructing the isosurface，the Marching Cubes algorithm detects all voxels，even if some voxels do not intersect with the isosurface，which affects the efficiency of the algorithm.In addition，the Marching Cubes algorithm ignores some points that are originally on the isosurface in this process，which reduces the accuracy of surface reconstruction.In order to solve these problems，this paper improves the algorithm.When constructing isosurface，it does not detect empty voxels to increase the efficiency of the algorithm，and adds some neglected equivalents to improve the accuracy of algorithm.

Keywords：3D surface modeling；Marching Cubes algorithm；algorithmic efficiency；approximation accuracy

1 引言（Introduction）

隨著科學計算可視化和計算機硬件技術水平的發展，三維表面建模技術越來越多地應用于生活、工程和科研領域[1]。在三維表面建模技術中，等值面抽取方法的重建效果比較好[2]，而Marching Cubes是其中比較經典的算法[3]。在醫學領域，醫生需要在核磁共振圖像（MRI）的幫助下，掌握患者身體內部的病變情況，對病人進行治療[4，5]；在地質領域，地質專家需要通過地質勘測圖像了解地層結構，指導礦藏開采或者預測地質災害[6，7]。這些圖像都可以由Marching Cubes算法重建獲得。然而Marching Cubes算法也存在一些不足之處，并且不斷有學者和研究人員進行研究和改進[8，9]。

本文以Marching Cubes算法為基礎，對算法進行了一定的改進，提高了算法的效率和逼近精度。利用三維礦體剖面輪廓線數據，對三維礦體表面進行重建，還原了礦體原貌。在通過實驗圖像驗證了算法可行性之后，對算法進行了時間復雜度分析，并且總結了算法的優缺點。

2 基本原理（The basic principle）

本文的基礎數據是礦山剖面輪廓線數據，在進行三維表面建模之前，先要對此數據進行處理以生成體數據。生成體數據之后，用改進的Marching Cubes算法構造等值面，完成建模工作。主要步驟有兩個：（1）構造體數據，這個步驟中排除了對空體素的檢測；（2）生成等值面，這個過程中增加了等值點。

2.1 構造體數據

體數據是體素級表面建模的基礎。本算法將每一層剖面都定義為一個Nx×Ny的二維網格，所有剖面垂直疊加就形成一個規模為Nx×Ny×Nz的三維空間區域，算法過程就在這個區域中進行。

在一系列二維平面上，定義場函數為f（x，y），對于平面上的某個坐標為（x，y）的網格點P的場函數：（1）如果該點在剖面輪廓線內部，其值為正；（2）如果該點位于剖面輪廓線的外部，其值為負；（3）如果該點落在剖面輪廓線上，其值為0。為了使得重建結果盡量平滑而且逼近物體真實表面，算法要選擇合適的場函數，否則會出現明顯的突變，造成繪制結果失真。本文采取歐氏距離函數作為場函數：

dis（x，y）表示在這個二維網格中，該網格點到該剖面上的輪廓線上點的最近距離。

在同樣大小的三維空間區域中，體數據的密度會影響Marching Cubes算法重建結果的逼近精度[10]，因此在構造體數據的過程中，合理劃分網格單元的大小對算法結果有直接的影響。在對算法處理時間和存儲空間允許的情況下，盡量將網格劃分的精細一些，可以得到更為精細的重建結果。

在Marching Cubes算法中，等值面與體素棱邊的交點用線性插值求取。通常在計算場函數值時，要對所有體素端點都進行計算，但是實際上與等值面相交的體素比起沒有與等值面相交的體素，數量要少很多[11]。本文是基于剖面數據進行表面重建，因此在剖面上考慮此問題，對原算法進行了改進。假設剖面是一個Nx×Ny的二維網格，在某一個剖面上，與等值線相交的網格數量占整個網格的比例很小。

如圖1所示，在一個剖面上，與等值線相交的網格數量所占比例很小，如果計算所有網格頂點的場函數值，無效計算很多，將會降低算法的效率。為了改善這一狀況，提高算法的效率，需要盡量減少無效的計算。為此要設置一個標識符，一個表示各頂點狀態值的數組，以及一個存放場函數值的數組，首先判斷頂點與等值面的關系，若在等值線內（包括在等值線上），則其狀態值為1，反之狀態值為0。當確定了一個網格四個頂點的狀態值后，如果狀態值全為0或者全為1，則說明該網格與等值面無交點，因此不需要計算其場函數值。步驟如下：

Step1：輸入剖面數據；

Step2：將剖面數據網格化；

Step3：計算每個網格四個頂點的狀態值，通過頂點狀態值確定是否與等值面相交，計算與等值面相交的網格的四個頂點的場函數值；

Step4：逐個處理每個剖面，構造出體數據。

2.2 生成等值面

基于體素級的三維表面建模，即基于等值面生成的表面建模，需要從大量的體數據中把近似表示物體表面的等值面提取出來。在Marching Cubes算法中，在一個體素中等值面的提取過程如下：首先計算體素每個頂點的狀態值，確定該體素是否是和等值面相交的體素；確定了體素和等值面的關系之后，對于和等值面相交的體素，計算其八個頂點的場函數值，然后根據已經給定的等值面閾值求取交點（即等值點），求出體素內所有的交點后，把這些交點按照一定的方式連接起來，再進行三角剖分，就生成了該體素內的等值面。Marching Cubes算法的前提是數據場沿著體素棱邊呈線性變化，求取等值點時根據其所在棱邊兩個端點的坐標進行插值即可求得。

體素的每個頂點都有兩種狀態，要么在等值面內，要么在等值面外，而8個頂點，則共有28種狀態，也就是說根據頂點狀態的不同，可以把體素類型分成256種，記錄256種狀態比較復雜，而且這些狀態中很多是可以相互轉換的，根據旋轉對稱性和互補對稱性進行簡化，只需要記錄15種等值面連接構型[12]。在求出體素與等值面的交點后，再根據這15種等值面構型連接三角形組成最終的等值多邊形。

將數據場內所有和等值面相交的體素處理完后，這些體素內部的等值面組合起來就構成了所要重現物體的表面網格模型，而可視化最終要將抽象數據轉變為容易理解的圖像信號，因此為了更好地逼真物體真實表面，需要對網格模型進行光照處理。處理方法為在體素內構造等值面的同時，也要計算每個剖分三角形三個頂點的法向量，然后求取平均值作為該三角面片的法向量，利用此法向量就可以計算該面片上的光照強度，最終得到逼近真實物體的繪制結果[13]。

Marching Cubes算法在提取等值面時，通過插值計算出體素棱邊與等值面的交點，然后再通過一定的方式把插值點連接起來形成多邊形近似表示物體表面。在實際情況中，等值面和體素表面的交線是雙曲線，Marching Cubes算法以連接插值點的方法，用直線段近似表示曲線段生成物體表面。在這個過程中，原算法并不考慮是否有原輪廓線上的點落入體素之中，因此插值后形成的等值輪廓線和原輪廓線之間存在一定的誤差，而原輪廓線上的點是物體表面上的點，因此使用等值輪廓線提取出的等值面與物體原表面也存在一定的誤差。當體數據場密度較大的時候，利用直線段來近似表示曲線段是可行的，因為當線段長度非常小時，直線段和曲線段的逼近程度會很高，小到一定程度的時候肉眼是無法分辨的。但是如果體數據場密度較小，插值求取交點后連接形成的直線段長度較大，此時其和實際的曲線段交線相比較，逼近程度就會差很多，這個逼近程度會隨著體數據場密度的減小而降低。

如圖3所示，網格中的閉合實曲線是原輪廓線，虛線是經過插值計算后連接插值點形成的等值輪廓線，實心黑點是原輪廓線上的數據點。從圖中可以看到原輪廓線上的點不全在等值線上，這樣就造成連接插值點形成的等值線和原輪廓線吻合程度較低。

為了降低逼近誤差，本文在構造等值面時，除了計算體素棱邊和等值面的交點之外，還考慮了原輪廓線上的點，對于有原輪廓線上的點落入的體素，將插值點和原輪廓線上的點以一定的方式連接起來，然后再通過三角剖分生成等值面。通過這一方法，可以使得剖面數據中原輪廓線上的點包含在等值線上，而且在連接的時候加入原輪廓線上的點，那么生成的等值線段就不再是是線段，而是由兩條直線段組成的折線段。如圖4為一個有原輪廓線點落入的網格中等值線的連接方式，曲線段P1P2是等值面和該體素表面的實際交線，直線段P1P2是傳統Marching Cubes算法經過插值計算出交點坐標后連接形成的等值線段，Q是落入該表面上的原輪廓線上的點，連接的時候如果將點Q考慮在內，那么最終以折線段P1QP2作為等值線，很明顯可以看到折線段P1QP2和點O圍城的多邊形面積更接近曲線段P1P2和點O圍城的扇形面積。

連接插值點形成等值線時，如果加入原輪廓線上的點，可能會改變原來的連接方式。對于體素的某個面來說，如果原輪廓線上的點落在該面的某條邊上，那么計算插值點的時候，插值點剛好就是這個原輪廓線上的點，此時即使將該原輪廓線上的點考慮在內，也不會改變該面上插值點的連接方式。但是如果原輪廓線上的點完全落入該面內，那么就會改變插值點的連接方式。如圖5所示，圖a和圖b分別是調整前后的等值點連接方式，調整后左邊和右邊兩個網格內等值點的連接方式都發生了改變，中間的網格由于沒有原輪廓線上的點落入，因此連接方式沒有改變。

根據算法前提假設，場函數沿著體素棱邊呈線性變化，因此體素的一個表面四條邊與等值面的交點數量有0、2和4這三種情況。當體素某個表面四條邊與等值面的交點數量為0的時候，必定沒有原輪廓線上的點落入；當交點數量為2的時候，如果有原輪廓線上的點落入，那么改變連接方式比較簡單，只需要依次連接三個點即可。如果有四個交點，該面是二義性面，在連接的同時也要考慮解決面二義性問題，不同的連接方式可能會生成不同的等值面構型，甚至可能會出現拓撲錯誤的情況[14]，比較經典的解決方法有四面體剖分法[15]和雙曲線漸近線法[16]，不同的方法各有特點，本文采取了一種基于插值點連線交點的方法解決面二義性問題[17]。

體素級表面建模在提取等值面的過程中，是以直線段近似代替曲線，而此處將原輪廓線上的點考慮進來，實際上是以兩段折線段近似代替曲線，逼近程度要比之前高一些。在圖6（a）中可以看到，有兩個交點的情況有兩種，在每種情況中，考慮原輪廓線上的點之后，連接方式是唯一的，將原輪廓線上的點分別和兩個交點連接起來即可。但是在二義性面上，因為內部等值直線段有兩條，因此可能的連接方式有兩種，如圖6（b）所示，兩種不同的連接方式會導致體素內的三角面片拓撲關系不同。在這種情況下，結合輪廓線和交點的關系，計算該輪廓線上點到兩個直線段的距離，取距離較近的那個直線段為需要調整的直線段，因為距離近的直線段其上面點的場函數值更接近等值面閾值，這樣更符合實際情況。

在調整了體素表面上的等值點連接方式之后，體素內部三角剖分形成等值面的方式也隨著會發生改變，不能再完全按照之前介紹的15種基本構型來構造等值面了。但是如果全部重新調整等值面構型會增加較大的計算量，而調整插值點連接方式的面上，調整方式都是相對簡單的，因此可以在原來15種構型的基礎上進行改進。

如圖7所示，以Marching Cubes算法中的一種等值面構型為例。圖（a）是原來的等值面構型，該體素中一個頂點在等值面內（外），其余七個頂點都在等值面外（內）。按照經典Marching Cubes算法原理，在該體素中有三條棱邊和等值面相交，插值計算出交點坐標后，連接這三個交點構成一個三角形面片，最終以該三角形面片作為該體素內的逼近等值面。但是當加入原輪廓線上的點之后，上表面的連接方式發生了改變，此時就不能再按照圖（a）的方式構造等值面了。此時要在原來構型的基礎上，如圖7（b）所示，點P1、P2和P3分別是插值計算出的交點，Q是原輪廓線上的點，做如下調整：首先確定該體素中插值點的連接方式，也就是說要按照傳統Marching Cubes算法的方法確定等值面構型，然后通過點P1和P2找到原構型中和這兩個點處于同一個三角形面片中的另外一個頂點P3，之后分別連接QP1、QP2和QP3，最后去掉P1和P2之間的連線，這樣就完成了調整。由圖可知加入原輪廓線上的點Q后，通過調整等值面構型，使用兩個三角形面片

△QP1P3和△QP2P3代替了原來的三角面片△P1P2P3，這樣就使得重建后的逼近精度有了一定的提高。

在原Marching Cubes算法15種基本構型中，有些體素構型中可能包含多個三角形等值面片，0號構型中沒有三角形面片，1號構型中有1個三角形面片，2、3、4、8號構型中都有2個三角形面片，5、6、7號構型中包含三個三角形面片，9—14號構型中分別包含四個三角形面片。當體素中三角形面片數量大于等于2的時候，就可能出現有公共邊的三角面片。但是通過觀察這些構型，可以發現在所有存在公共邊面片的體素中，公共邊只在體素內部，而在體素表面上不存在包含于兩個或兩個以上三角形面片中的等值直線段。本文是基于剖面輪廓線數據進行三維表面重建，在構造體數據的時候，要將相鄰兩層剖面進行體素化，因此剖面上的原輪廓線數據只分布在體素的上下表面上，由于體素表面上的等值直線段不存在被幾個三角形面片共用的情況，因此在連接插值點和原輪廓線上的點時，不需要考慮圖7（b）中P1P2被幾個三角形面片共用的情況。

經過上面的處理，在構造等值面的時候既考慮了體素棱邊與等值面的交點，也考慮了原輪廓線上點的分布，因此在一個剖面上，最終形成的等值輪廓線與原輪廓線之間的誤差相比之前有了一定的減小，一定程度上提高了算法的逼近精度，圖8是對圖3進行處理前后的對比，可以明顯看到二者的區別。

綜合上述步驟，得出構造等值面的步驟如下：

Step1：讀入一個體素數據；

Step2：根據體素八個頂點的場函數值插值計算棱邊上的交點；

Step3：根據頂點狀態索引值查找對應的邊狀態索引值，再根據邊狀態索引值查找該體素內的等值面構型，確定等值面連接方式；

Step4：確定該體素表面是否包含原輪廓線上的點，如果包含原輪廓線上的點，則執行Step5，否則執行Step8；

Step5：查找包含原輪廓線點Q的面，計算該面上插值求出的等值點個數n，如果n=4，則執行Step6，如果n=2則執行Step7；

Step6：按照Step3中已經確定的連接方式，計算點Q到該面上兩條等值線段的距離，取距離較近的等值線段進行調整；

Step7：查找等值線段P1P2所在三角形面片的另外一個頂點P3，分別連接QP1、QP2和QP3，并且去掉等值線段P1P2；

Step8：算法結束。

2.3 基于剖面數據進行三維表面重建的過程

根據之前對算法的介紹，基于剖面數據的改進MC表面建模流程如下：

Step1：讀入剖面數據；

Step2：將相鄰兩層剖面數據網格化，得到體素；

Step3：選中體素，求取體素頂點狀態值，判斷體素是否與等值面相交；

Step4：對于和等值面相交的體素，求取每個頂點的場函數值；

Step5：檢查體素是否還有二義性面，對二義性面進行處理；

Step6：根據頂點狀態值找到相應的邊狀態值，根據邊狀態值查找對應的等值面連接方式，對于有原輪廓線上的點落入的體素，對其中的等值面連接方式進行調整；

Step7：計算剖分三角形的法向量，計算光照，繪制該三角形對應的等值面片；

Step8：選擇新的體素轉入Step3繼續處理，直到所有體素都處理完畢；

Step9：結束。

3 實驗結果（Experimental result）

本文在對剖面數據進行體素化構造體數據的過程中，將剖面網格劃分為規模大小為Nx×Ny的網格，假設每個剖面上平均有m條礦體輪廓線剖面，每條輪廓線上有n個點。下面對算法效率和結果進行分析說明。

（1）算法的時間復雜度

算法主要時間用于構造體數據和生成等值面兩部分，這兩個部分中有大量的計算和比較。體數據的構造由判斷體素頂點與等值面關系、計算頂點場函數值兩部分組成。生成等值面時需要查詢邊狀態表和等值面結構表。

在構造體數據時，計算量主要集中網格頂點狀態值和場函數值的計算，計算頂點狀態值的次數為Nx×Ny，所以時間復雜度為O（NxNy）；計算場函數的次數為4×m×n，由于m<

O（NxNy）+O（N）=O（NxNy）

生成等值面的過程中，需要在每個體素內提取等值面，計算次數Nx×Ny；而查詢邊狀態表和等值面結構表的時間復雜度都是是線性的O（N），因此二者相加得到生成等值面等值面的總體時間復雜度：

O（NxNy）+O（N）+O（N）=O（NxNy）

算法的總體時間復雜度由構造體數據和生成等值面兩部分的時間復雜度相加所得，近似為O（NxNy）。

（2）算法實驗效果

本文以礦體剖面輪廓線數據為實驗數據，通過實驗對算法進行了仿真和驗證，編程環境為VS2010，結果如下。

使用本文算法對礦體剖面數據進行處理，正確重建了礦體的三維表面模型，消除了二義性帶來的孔洞現象，效果如圖9所示為兩種礦體剖面數據重建后的線框模型和表面模型。

在重建過程中，隨著網格規模的增加，圖形的逼近程度也會增加，但是由此產生的計算量增大，因此所耗時間也會增加。圖10是對同一種礦體的剖面數據采取不同規模的網格重建后的結果。

如表1所示，同樣的剖面數據，采用不同的網格規模進行重建，當網格規模為20×20時所需時間為0.823s，當網格規模為40×40時所需時間為1.452s，因此要得到效果較好的重建圖形并且盡量節省時間，需要合理選用網格規模。

4 結論（Conclusion）

本文在Marching Cubes算法的基礎上，改進了原算法，不再對空體素頂點場函數值進行計算，提升了算法效率，并且通過增加交點的方式，提高了算法的逼近精度，此外也存在一些不足，歸納如下：

減少了無效的場函數值計算，傳統Marching Cubes算法在尋找等值點時要計算所有的體素頂點，包括很多本來和等值面不相交的體素也進行了計算，本算法通過設置標志位只對和等值面相交的體素的頂點計算了場函數值，減少了無效計算，提升了算法效率。

在構造等值面的過程中，傳統Marching Cubes算法只考慮插值計算求出的交點，將這些交點按照體素對應的等值面構型連接起來表示等值面。為了提高重建結果的逼近精度，本文在構造等值面時，將原輪廓線上的點和插值點一起連接起來，使得原輪廓線上的點包含于生成的等值面之中，并且由于使用折線段代替直線段近似表示等值線，使得算法重建的逼近精度有了一定的提高。

當數據場數據密度特別小時，沒有解決逼近精度低的問題，傳統Marching Cubes算法在較高密度的數據場中可以得到很好的表面繪制結果，但是隨著數據場數據密度的降低，算法重建結果的逼近程度也隨著降低，本文在數據場數據密度特別小時沒有做相應的處理，將作為后續問題繼續研究。

參考文獻（References）

[1] Minho Chang，et al.Interactive marching cubes algorithm for intraoral scanners[J].The International Journal of Advanced Manufacturing Technology，2017，89（5-8）：2053-2062.

[2] 黃偉.三維表面建模方法的研究與實現[D].長沙：中南大學，2010.

[3] William E.Lorensen，Harvey E.Cline.Marching cubes：A high resolution 3D surface construction algorithm[J].ACM SIGGRAPH Computer Graphics，1987，21（4）：163-169.

[4] 于洋.肺部CT血管分割及三維重建[D].哈爾濱工業大學，2016.

[5] 姚翠萍，周湘連，王晶，等.納米金標記細胞的熒光壽命成像及其三維重建[J]. 西安交通大學學報，2016，50（04）：153-158.

[6] 劉致寧，宋承云，李志勇，等.基于多地震屬性融合分割的地質異常體三維模型構建[J].Applied Geophysics，2016（3）：519-528.

[7] 鄒艷紅，何建春.移動立方體算法的地質體三維空間形態模擬[J].測繪學報，2012，41（06）：910-917.

[8] Seungki Kim.A Novel Interpolation Scheme for Dual Marching Cubes on Octree Volume Fraction Data[J].Computers & Graphics，2017，66（8）：169-178.

[9] 畢碩本，陸源，曾曉文，等.Marching Cubes改進算法及其氣象三維模擬[J]系統仿真學報，2017，29（07）：1405-1410；1418.

[10] 孫偉，張彩明，楊興強.Marching Cubes算法研究現狀[J].計算機輔助設計與圖形學學報，2007，19（7）：947-952.

[11] 錢峰，馬秀麗，楊勝齊，等.移動立方體算法的研究和改進[J].計算機工程與應用，2010，46（34）：177-180.

[12] Adriano Lopes and Ken Brodlie.Improving the Robustness and Accuracy of the Marching Cubes Algorithm for Isosurface[J].IEEE Transactions on Visualization and Computer Graphics，2003，9（1）：l6-29.

[13] 顧耀林，呂理偉.移動立方體算法中的三角剖分[J].計算機工程與設計，2006，27（1）：120-123.

[14] 梁秀霞，張彩明.拓撲結構正確的三線性插值曲面的三角片逼近[J].計算機研究與發展，2006，43（3）：528-535.

[15] Hummel R. Exploiting Triangulated Surface Extraction Using Tetrahedral Decomposition[M].IEEE Educational Activities Department，1995，1（4）：328-342.

[16] Nielson G，Hamann B.The asymptotic decider：resolving the ambiguity in marching cubes[C].Proceedings of Visualization' 91，Los Alamitos CA，1991：83-91；413.

[17] 王錚，李瑞明.移動立方體算法面二義性問題研究[J].軟件工程，2017，20（09）：6-8；5.

作者簡介：

王 錚（1985-），男，碩士，助教.研究領域：計算機圖形學.