基于GPU的MC算法

李丹 彭海欣

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

關(guān)鍵詞：MarchingCubes；醫(yī)學(xué)圖像；相鄰體元；GPU

中圖分類號：TP391.41 文獻(xiàn)標(biāo)識碼：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

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

1 MarchingCubes算法基本原理

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

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

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

1.2 提取等值面

分析立方體的256種狀態(tài)可以發(fā)現(xiàn)，立方體頂點(diǎn)狀態(tài)具有反對稱性和旋轉(zhuǎn)性，即旋轉(zhuǎn)不會(huì)改變等值面的拓?fù)浣Y(jié)構(gòu)，反對稱性就是所有頂點(diǎn)狀態(tài)相反不會(huì)改變等值面的連接方式。根據(jù)這種情況可以將256種情況化簡為15種情況，如圖1所示。

可以建立邊索引表，將每條邊進(jìn)行編號，邊索引中代表這256種狀態(tài)下每條邊的狀態(tài)，根據(jù)8個(gè)頂點(diǎn)的狀態(tài)可以得到一個(gè)頂點(diǎn)索引值，根據(jù)該索引值可以查看邊索引表，得到等值邊情況。同時(shí)，還可以建立三角形索引記錄256種狀態(tài)下每個(gè)三角面片的頂點(diǎn)情況（用頂點(diǎn)編號記錄頂點(diǎn)，每3個(gè)代表一個(gè)三角面片），用得到的邊索引值獲取三角形索引的三角面片情況。根據(jù)8個(gè)頂點(diǎn)的狀態(tài)獲取頂點(diǎn)索引值，根據(jù)索引值查找邊索引，即現(xiàn)在狀態(tài)下的12條邊情況，獲取相對應(yīng)邊索引值。再根據(jù)2個(gè)索引值查找三角索引，獲取該立方體的等值面情況。

1.3 計(jì)算等值點(diǎn)信息

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

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

2 基于GPU的改進(jìn)的MC算法

2.1 改進(jìn)的MC算法算法原理

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

若某一立方體的某一面存在著等值邊，又因?yàn)槟Ｐ褪窍噙B通的，則該面相鄰的立方體的同一面必定存在等值邊，除了該立方體正好是邊界的情況外，含有等值邊的立方體一定存在著等值面。因此，通過判斷立方體面上是否有等值邊，可判斷相鄰立方體中一定存在等值面的立方體。逐步遍歷這些非空立方體，就可避免對非空立方體的處理。又因?yàn)槟Ｐ褪沁B通的，所以基本不會(huì)遺漏主要的特征的三角網(wǎng)格。

在觀察立方體的狀態(tài)時(shí)可以發(fā)現(xiàn)，若一條邊的2個(gè)頂點(diǎn)狀態(tài)相反時(shí)必定存在一個(gè)等值點(diǎn)，當(dāng)立方體的一個(gè)面存在2個(gè)等值點(diǎn)時(shí)，必定存在等值邊。由此可以發(fā)現(xiàn)，當(dāng)某個(gè)面的4個(gè)頂點(diǎn)存在2個(gè)標(biāo)記、2個(gè)不標(biāo)記或者3個(gè)頂點(diǎn)狀態(tài)相同、另一頂點(diǎn)狀態(tài)相反時(shí)，該面必定存在等值邊，如圖2所示。

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

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

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

2.2 基于GPU的并行化處理

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

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

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

3 試驗(yàn)分析

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

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

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

4 結(jié)語

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

參考文獻(xiàn)：

[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的醫(yī)學(xué)圖像三維重建[J].計(jì)算機(jī)工程與應(yīng)用，2012（18）：154-158.

[3]高峰，付忠良.基于改進(jìn)移動(dòng)立方體的醫(yī)學(xué)圖像三維重建算法[J].計(jì)算機(jī)應(yīng)用，2013（S1）：201-203，213.

[4]]王旭初，王贊，牛彥敏，等.融合構(gòu)型查找表與鄰接查找子表的改進(jìn)MC方法[J].重慶大學(xué)學(xué)報(bào)，2012（12）：68-77，83.

[5]王溪，秦新強(qiáng)，黨發(fā)寧，等.醫(yī)學(xué)圖像三維重建的規(guī)則移動(dòng)立方體法[J].2009（4）：477-481.