基于GPU的MC算法

李丹 彭海欣

摘 要：MarchingCubes算法是醫學圖像三維重建獲取等值面常用的方法，該方法原理簡單，但在遍歷立方體的時候會遍歷許多空立方體，浪費時間。本文提出一種基于頂點狀態獲取相鄰非空體元的方法，根據頂點的狀態判斷某一面是否含有等值邊，含有等值邊的那面相鄰體元有很大可能含有等值邊，減少空立方體的遍歷。此外，本文還將算法在基于GPU的基礎上進行并行化處理。

關鍵詞：MarchingCubes；醫學圖像；相鄰體元；GPU

中圖分類號：TP391.41 文獻標識碼：A 文章編號：1003-5168（2017）03-0057-04

Abstract： MarchingCubes algorithm is a commonly used method for 3D reconstruction of medical images. The method is simple， but a lot of empty cubes will be traversed when traversing the cube， which wastes time. This paper presented a method based on vertex state to obtain adjacent non - empty body elements and determine whether a side contains the equivalent edge according to the state of the vertex to. It is very likely that the adjacent surface element with the same edge can contain the equivalent edge， which reduces the ergodic of the empty cube. In addition， the algorithm was parallelized based on GPU.

Keywords： MarchingCubes；medical image；adjacent voxel；GPU

近年來，三維重建技術迅速發展，其應用領域十分廣泛，其中在醫學方面有著重要的學術意義及實踐意義。移動立方體MarchingCubes算法[1]（簡稱MC算法）原理簡單，實現容易，逼近精度較高[2]。但是其仍存在不足之處，遍歷立方體時會存在空立方體，浪費大量精力，避免對空立方體的遍歷，可以提高三維重建速度。由高峰等[3]提出的選取種子體元后根據種子體元衍生出整個器官的等值面，王旭初[4]提出建立接查找子表約束體元搜索路徑并根據鄰接查找子表特點設計堆棧結構實現搜索算法。王溪[5]提出了表面判斷查找表尋找等值面，不需要判斷體數據中的每個立方體，但需耗費存儲空間存儲記錄256中每個狀態下每個三角形的等值邊與面的相對情況，以及計算每條邊的情況。本文在此基礎上直接通過頂點狀態實現判斷，而不需要建立查找表及計算等值邊。GPU從出現后發展飛速，并且被利用到并行處理中，而GPU強大的計算能力及高效處理效率，可以更快提高MC方法的處理速度。因此，基于GPU的MC方法可以實現更快的重建速度。

1 MarchingCubes算法基本原理

MC算法也就是移動立方體算法，就是利用一個個立方體來獲取等值面信息，而等值面則被化為三角面片，MC算法獲得的就是三角面片信息。立方體也稱體元，是由相鄰8個頂點構成，掃描每個體元的8個頂點信息，根據給定的閾值，若包含閾值范圍，則體元中必定存在等值面，等值面與體元相交的點成為等值點，相交的邊稱為等值邊。將等值點連接則生成所需的等值面，具體過程如下。

1.1 確定含有等值面的立方體

將三維離散數據場看作立體柵格系統，圖像信息讀入后，相鄰兩幅圖像構成一層，各自相對應相近的4個點構成一個立方體。由于數據場沿著立方體是線性變化的，所以，若一條邊上的2個頂點的灰度值分別大于、小于所規定的閾值，則該邊上有且僅有一個與等值面相交的點，即等值點。立方體8個頂點的灰度值分別與規定的閾值比較，以此來判斷等值面是否與該立方體相交。根據閾值可以將頂點標記為標記和不標記2種狀態：若頂點的閾值大于規定的閾值，則頂點在等值面外部，不標記該頂點；若頂點的閾值小于規定的閾值，則頂點在等值面內部，標記該頂點。每個頂點有2種狀態，則1個立方體的8個頂點有256種狀態。

1.2 提取等值面

分析立方體的256種狀態可以發現，立方體頂點狀態具有反對稱性和旋轉性，即旋轉不會改變等值面的拓撲結構，反對稱性就是所有頂點狀態相反不會改變等值面的連接方式。根據這種情況可以將256種情況化簡為15種情況，如圖1所示。

可以建立邊索引表，將每條邊進行編號，邊索引中代表這256種狀態下每條邊的狀態，根據8個頂點的狀態可以得到一個頂點索引值，根據該索引值可以查看邊索引表，得到等值邊情況。同時，還可以建立三角形索引記錄256種狀態下每個三角面片的頂點情況（用頂點編號記錄頂點，每3個代表一個三角面片），用得到的邊索引值獲取三角形索引的三角面片情況。根據8個頂點的狀態獲取頂點索引值，根據索引值查找邊索引，即現在狀態下的12條邊情況，獲取相對應邊索引值。再根據2個索引值查找三角索引，獲取該立方體的等值面情況。

1.3 計算等值點信息

獲取立方體等值面位置后，就可以計算等值點坐標及法向量等信息。等值點計算采取插值計算，常用的有線性插值法和中點選擇法。中點選擇法是取該等值邊的中點，若用p1、p2分別代表2個頂點的坐標，n1、n2代表2個頂點的法向量，則公式可表示為：

MC算法的算法步驟為：①逐步掃描兩層數據圖像，構建體元；②將體元的所有頂點的灰度值與規定的閾值做比較，建立頂點索引；③若體元含有等值面，利用灰度差分法計算立方體各頂點的梯度；④根據邊頂點索引查找在邊索引表，獲得存在等值點所在的邊；⑤根據該相交的邊的2個頂點，用線性插值法獲取等值點的坐標信息及法向量；⑥根據索引值查找三角面片索引表，確定該立方體內三角面片的頂點連接方式；⑦由三角面片構成等值面；⑧重復以上步驟直到所有的體元遍歷結束。

2 基于GPU的改進的MC算法

2.1 改進的MC算法算法原理

一些MC算法中為提高計算速度，通常采用中點法來計算等值點，其計算量比線性插值法少。此外，大多數是建立在減少空立方體的遍歷上。在獲取圖像信息時多數區域信息是無效的，因此會建立許多空立方體，若避免對這些立方體的遍歷，基本會減少1/2的時間。本文算法主要思想為確定一個種子體元，以此體元為為起點，遍歷相鄰體元。

若某一立方體的某一面存在著等值邊，又因為模型是相連通的，則該面相鄰的立方體的同一面必定存在等值邊，除了該立方體正好是邊界的情況外，含有等值邊的立方體一定存在著等值面。因此，通過判斷立方體面上是否有等值邊，可判斷相鄰立方體中一定存在等值面的立方體。逐步遍歷這些非空立方體，就可避免對非空立方體的處理。又因為模型是連通的，所以基本不會遺漏主要的特征的三角網格。

在觀察立方體的狀態時可以發現，若一條邊的2個頂點狀態相反時必定存在一個等值點，當立方體的一個面存在2個等值點時，必定存在等值邊。由此可以發現，當某個面的4個頂點存在2個標記、2個不標記或者3個頂點狀態相同、另一頂點狀態相反時，該面必定存在等值邊，如圖2所示。

若存在上述狀態即只要一個面的頂點狀態不全相同，就可判定該面含有等值邊，同時也就可以判定該面相鄰立方體為非空立方體，則可以遍歷相鄰該方向上的立方體。可以對立方體的頂點和邊進行編號，以及根據立方體6個面定義6個方向，如圖3所示。

若某一面存在等值邊，則該方向標記為1，否則為0，并將該方向的體元壓入棧中，而下一個遍歷的體元就從棧中提取。為了防止同一體元壓入棧中或者壓入處理過的體元，應對體元設置標記1代表不需要再壓入體元，否則為0，在每次壓入棧前都需查看該體元狀態標記。每個體元在處理過程中不僅需要獲取等值面信息，還需獲取相鄰未處理過的存在等值面的體元信息并將其壓入棧中，然后下一個體元不需要再去遍歷判斷是否是非空體元，直接從棧中讀取非空體元，在讀取體院同時刪除棧中該因素，同時該體元狀態更新為1。每遍歷一個體元，都可能壓入新的非空體元，需要不斷從棧中提取體元，直至棧中為空。

該算法的算法步驟為：①逐步掃面圖像數據構造體元；②從中間體元開始，按原MC算法找到第一個非空體元；③該體元獲取完等值面時，獲取6個面的方向標記索引；④若某個方向標記為1，并且該方向的體元的狀態0，則將該方向體元坐標壓入棧中；⑤下個處理體元從棧中讀取，讀取完后將該體元信息從棧中刪除，并將該體元的狀態更改為1；⑥以此讀取棧中體元，直到棧為空。

2.2 基于GPU的并行化處理

GPU（Graphics Processing Unit）是圖形處理器的簡稱，這個概念是由NVIDIA公司在發布GeForce256繪圖處理芯片時首先提出。GPU與CPU類似，為復雜的數學計算和幾何計算而制造，在浮點計算及并行計算等方面有著比CPU高幾十倍的性能。因此，對于GPU對于浮點計算的高性能，能對MC算法的插值計算等計算步驟進行更快的計算，并且利用GPU的海量線程進行并行計算，在基于GPU的MC方法上具有更高的效率。將原始的MC方法及本文中改進的MC方法在基于GPU的環境下進行并行。

2.2.1 原MC算法并行化。原始MC算法主要借用GPU的浮點計算能力，其算法流程如下：①并行的獲取醫學圖像數據；②根據獲取的醫學圖像構建體元；③根據頂點狀態判斷是否含有等值點，若沒有返回上一步，繼續遍歷下一個體元，若有等值面則進行下一步；④獲取該體元內的等值面信息；⑤GPU中計算，通過線性插值方法，計算出體元邊界與等值而的交點，然后利用中心差分法，求出體元各角點處的法向，再通過線性插值方法，求出三角形各頂點處的法向；⑥遍歷體元直至獲取完所有體元等值面；⑦利用三角面片信息繪制模型。

2.2.2 基于頂點狀態的MC算法并行化。改進的在讀入離散數據構造立方體后，首先找出第一個非空立方體后，分別對6個面進行是否有等值邊的狀態判斷，并進行標記，獲取完該體元的等值面后根據等值邊標記將標記為1的方向的相鄰體元壓入棧中。計算同樣是在GPU中進行的。其算法流程如下：①并行獲取兩層數據構建立方體；②并行的從立方體中找出第1個非空立方體，從中間開始查找；③對立方體按照正常MC方法進行并行獲取等值面信息；④對6個面的狀態進行判斷，獲取相鄰含等值面的體元信息并壓入棧中；⑤并行從棧中獲取待處理的體元，同時將該體元從棧中刪除以及將該體元狀態改為1；⑥多線程處理體元，重復從步驟③開始處理，在步驟④時，壓入棧前先判斷體元狀態是否為0，若為1，則不壓入棧中；⑦并行的重復處理棧中體元直至為空為止。

3 試驗分析

在visual studio 2013環境下，用OpenGL實現三維繪制，脊椎骨和頭顱閾值為60。由試驗結果可以看出，改進的MC算法與原MC算法相比，模型樣子基本不變，只是頭顱頸椎部分少了一部分邊緣瑣碎部分，這部分不影響模型本身，而脊椎圖像基本沒有什么變化（見圖4）。

由表1結果可以看出，該算法在原算法的基礎上節省了大約50%的時間，雖然可能忽略了一些三角面片，但原模型基本特征沒有發生變化。

并行化試驗結果見表2，并行化后的算法比原算法快了將近一倍的速度，很明顯比串行速度快；而改進后的MC算法時間并行化后比最初時間快了將近3倍，速度有了很大的提高。

4 結語

對于醫學圖像三維重建，本文在原來的MC算法基礎上設置立方體方向來描述相鄰立方體，用頂點狀態確定含有等值面的立方體方向，將其壓入棧中待處理，以此避免空立方體的遍歷。本文對其進行實驗研究，結果表明該算法比原算法提高了將近一半的速度。同時，還將算法進行并行化處理，可以發現并行化的繪制速度明顯提高。基于頂點狀態的MC算法是在圖像模型具有連通性的情況下，在一個模型中可能存在互相不連通的情況，需要建立多個種子節點。同時，算法需要三維數組記錄每個立方體的是否處理情況，由于三角面片數量多，所以三維數組占用空間內存較大。

參考文獻：

[1]Willian E Lorense，Harvey E Cline. Marching Cubes：A High Resolution 3d Surface Construction Algorithm[J].Computer Graphics，1987（4）：163-169

[2]王旭初，王贊.基于最近鄰Marching Cubes的醫學圖像三維重建[J].計算機工程與應用，2012（18）：154-158.

[3]高峰，付忠良.基于改進移動立方體的醫學圖像三維重建算法[J].計算機應用，2013（S1）：201-203，213.

[4]]王旭初，王贊，牛彥敏，等.融合構型查找表與鄰接查找子表的改進MC方法[J].重慶大學學報，2012（12）：68-77，83.

[5]王溪，秦新強，黨發寧，等.醫學圖像三維重建的規則移動立方體法[J].2009（4）：477-481.