吳建偉 ，王子牛

(1. 貴州大學(xué)計(jì)算機(jī)科學(xué)與信息學(xué)院，貴州貴陽(yáng)550025;2. 貴州大學(xué)網(wǎng)絡(luò)與信息化管理中心中心，貴州貴陽(yáng)550003)

醫(yī)學(xué)圖像三維重建就是利用二維醫(yī)學(xué)圖像序列(如CT 切片，即計(jì)算機(jī)斷層掃描切片)重建出三維圖像模型，圖像三維重建的常用算法主要包括體繪制和面繪制。體繪制能表現(xiàn)更多結(jié)構(gòu)上的細(xì)節(jié)，但要達(dá)到滿意的精度則需要很大的運(yùn)算量，運(yùn)行和存儲(chǔ)性能難以得到滿足。于是，面繪制成為醫(yī)學(xué)圖像三維重建的主流算法。MC(Marching Cubes)算法是面繪制中的經(jīng)典算法，由W.E.Lorenson 和H.E.Cline 在1987 年提出［1］。MC 算法較好的解決了在標(biāo)量數(shù)據(jù)場(chǎng)中抽取等值面的問(wèn)題，適合在不規(guī)則的醫(yī)學(xué)體數(shù)據(jù)場(chǎng)中進(jìn)行等值面重建，并以大量的三角面片逼近真實(shí)的醫(yī)學(xué)三維圖像。但是，MC 算法也有一些缺點(diǎn)，主要在于抽取等值面時(shí)存在拓?fù)浣Y(jié)構(gòu)二義性［2-5］，處理的數(shù)據(jù)量大時(shí)對(duì)處理速度和存儲(chǔ)空間要求很高。隨著計(jì)算機(jī)硬件和并行計(jì)算的發(fā)展，后一問(wèn)題得到了一定程度的緩和，而本文主要關(guān)注前一問(wèn)題。作為MC 算法的發(fā)展，MT(Marching Tetrahedrons)算法［6］沒有等值面抽取的二義性問(wèn)題，但計(jì)算量是MC 算法的好幾倍。為了在不明顯增加計(jì)算量的前提下，較好的解決MC 算法的二義性問(wèn)題，本文結(jié)合兩種算法的優(yōu)點(diǎn)，提出一種改進(jìn)算法，在CT 圖像三維重建的實(shí)驗(yàn)中進(jìn)行了驗(yàn)證，并得到了較好的結(jié)果。

1 MC 算法

1.1 MC 算法的基本原理

假設(shè)有一個(gè)CT 切片數(shù)據(jù)集，可以看作是一個(gè)三維標(biāo)量數(shù)據(jù)場(chǎng)，其中每一個(gè)切片其實(shí)是一個(gè)二維數(shù)組，數(shù)組中的某一位置的值即該標(biāo)量數(shù)據(jù)場(chǎng)中該位置的標(biāo)量值。在相鄰兩層切片中分別各取4 個(gè)數(shù)據(jù)點(diǎn)，構(gòu)成一個(gè)小立方體，稱為立方體素，而整個(gè)標(biāo)量數(shù)據(jù)場(chǎng)可由這些立方體素的集合表示。MC算法就是對(duì)這些立方體素一個(gè)一個(gè)的進(jìn)行處理，抽取等值面。首先設(shè)定一個(gè)等值面值，對(duì)每一個(gè)立方體素的8 個(gè)頂點(diǎn)，若其頂點(diǎn)的標(biāo)量值:

(1)小于等值面值，則該頂點(diǎn)在等值面內(nèi)部;

(2)等于等值面值，則該頂點(diǎn)在等值面上;

(3)大于等值面值，則該頂點(diǎn)在等值面外部。

可將在等值面內(nèi)部的頂點(diǎn)標(biāo)記為1，其余頂點(diǎn)標(biāo)記為0，然后將8 個(gè)標(biāo)記值按一定順序排列，可由一個(gè)字節(jié)表示。由于每個(gè)頂點(diǎn)的標(biāo)記值有兩種可能，故每個(gè)立方體素的頂點(diǎn)狀態(tài)共有28=256 種可能。而每一個(gè)頂點(diǎn)狀態(tài)確定了立方體素各邊是否與等值面相交，從而確定立方體素三角面片的拓?fù)錉顟B(tài)。通過(guò)旋轉(zhuǎn)對(duì)稱性和互補(bǔ)對(duì)稱性可將一個(gè)立方體素的拓?fù)錉顟B(tài)歸納為以下15 種情形，如圖1 所示:

圖1 立方體素15 種基本拓?fù)錉顟B(tài)

1.2 立方體素的邊與等值面的交點(diǎn)

在確定立方體素的拓?fù)浣Y(jié)構(gòu)后，應(yīng)該計(jì)算出立方體素與等值面的交點(diǎn)?？梢约僭O(shè)數(shù)據(jù)場(chǎng)中各點(diǎn)的值沿立方體素邊界呈線性變化，這樣就可以通過(guò)線性插值的方法求出立方體素的邊與等值面的交點(diǎn)，即三角面片的頂點(diǎn)。設(shè)立方體素某邊與等值面相交于一點(diǎn)P，該邊的兩個(gè)端點(diǎn)坐標(biāo)為(x1，y1，z1)，(x2，y2，z2)，標(biāo)量值為V1，V2，等值面值為V，那么P點(diǎn)的坐標(biāo)(x，y，z)為:

1.3 三角面片的頂點(diǎn)法向量

在OpenGL 中繪制三角面片需要提供三角面片頂點(diǎn)的法向量，從而可以使用光照效果。由于數(shù)據(jù)場(chǎng)中各數(shù)據(jù)點(diǎn)的標(biāo)量值的梯度垂直于等值面，所以數(shù)據(jù)場(chǎng)中某點(diǎn)P 的法向量可用該點(diǎn)處的梯度代替。設(shè)P 點(diǎn)處的法向量為n(x，y，z)，其計(jì)算公式為:

其中，V(i，j，k)代表點(diǎn)(i，j，k)處的標(biāo)量值，Δx，Δy，Δz 為三個(gè)方向上的立方體素的邊長(zhǎng)。

三角面片的頂點(diǎn)是通過(guò)對(duì)邊上兩端點(diǎn)坐標(biāo)進(jìn)行插補(bǔ)得到的，如式(1)。同樣，在已經(jīng)求得各數(shù)據(jù)點(diǎn)的法向量后，可以通過(guò)對(duì)邊上兩端點(diǎn)的法向量進(jìn)行插補(bǔ)來(lái)得到三角片頂點(diǎn)的法向量。頂點(diǎn)法向量的計(jì)算公式仍為式(1)，但(x1，y1，z1)，(x2，y2，z2)分別表示兩端點(diǎn)的法向量，(x，y，z)為要求的頂點(diǎn)法向量。

1.4 MC 算法的二義性

在圖1 中，有些頂點(diǎn)狀態(tài)下，可以有兩種不同的交點(diǎn)連接方式，也就會(huì)得到不同的三角面片，即存在二義性問(wèn)題。MC 算法的二義性問(wèn)題包括立方體素面二義性和體二義性。

(1)立方體素面二義性

在立方體某一面的兩個(gè)對(duì)角頂點(diǎn)的標(biāo)量值小于等值面值，而另外兩個(gè)對(duì)角頂點(diǎn)的標(biāo)量值大于等值面值時(shí)，就會(huì)有兩種交點(diǎn)連接方式，如圖2 所示:

圖2 立方體素表面的不同交點(diǎn)連接方式

(2)立方體素體二義性

正如立方體素表面存在不同的交點(diǎn)連接方式，立方體素內(nèi)部也有類似情況，比如當(dāng)立方體素的兩個(gè)體對(duì)角頂點(diǎn)在等值面內(nèi)部，而其余頂點(diǎn)均在等值面外部時(shí)，有如圖3 所示兩種連接方式:

圖3 立方體素內(nèi)部的不同交點(diǎn)連接方式

2 MT 算法

MT 算法與MC 算法的不同在于它把一個(gè)立方體再細(xì)分為若干個(gè)四面體，常見的劃分方法將一個(gè)立方體劃分為5 個(gè)，6 個(gè)或24 個(gè)四面體。四面體是最簡(jiǎn)單的多面體，其等值面拓?fù)浣Y(jié)構(gòu)可歸結(jié)為3種，如圖4 所示:

圖4 四面體素的3 種拓?fù)錉顟B(tài)

類似于MC 算法，只要確定了四面體的4 個(gè)頂點(diǎn)的標(biāo)記狀態(tài)，就可以確定等值面與四面體邊的交點(diǎn)的連接方式，也就得到了用來(lái)逼近等值面的三角面片。重要的是，對(duì)于一個(gè)四面體來(lái)說(shuō)，其邊上的交點(diǎn)的連接方式是唯一的，即對(duì)于單個(gè)四面體來(lái)說(shuō)其交點(diǎn)的連接不存在二義性。但是若相鄰兩個(gè)立方體的剖分為四面體的方式不一致，則會(huì)產(chǎn)生剖分二義性。如圖5 所示，這兩個(gè)立方體剖分后，在其相接的面上兩條剖分線并不重合，這可能使得兩個(gè)立方體中生成的三角面片不能以一致的方式進(jìn)行連接，導(dǎo)致空洞的產(chǎn)生。

圖5 相鄰兩立方體剖分不一致

3 基于MC 和MT 算法的改進(jìn)算法

為了避免MC 算法的面和體二義性，需先判斷立方體是否具有面和體二義性，這可以用一個(gè)二義性表來(lái)枚舉所有256 種情況。若立方體不具有面和體二義性，可以用MC 算法進(jìn)行處理;否則，對(duì)立方體執(zhí)行MT 算法，把立方體剖分為多個(gè)四面體進(jìn)行處理［7］。然而，對(duì)立方體進(jìn)行剖分會(huì)產(chǎn)生剖分二義性問(wèn)題:一方面，相鄰的兩個(gè)剖分立方體之間可能拓?fù)洳灰恢拢缜懊娴膱D5 所示;另一方面，剖分立方體也可能與相鄰的未剖分的立方體拓?fù)洳灰恢?。如圖6 所示，兩個(gè)相鄰的立方體中的三角面片沒有嚴(yán)密的連接起來(lái)，從而形成了中間白色區(qū)域的空洞。

圖6 剖分立方體與相鄰的未剖分立方體拓?fù)洳灰恢?/p>

下面將分別討論這兩種情形。

3.1 相鄰兩剖分立方體之間拓?fù)洳灰恢?/h3>

圖7 立方體6 -剖分的一種方式

MT 算法常用的立方體剖分方式有5 -剖分，6-剖分和24 -剖分。為了避免相鄰兩個(gè)剖分立方體之間的拓?fù)洳灰恢滦裕梢允褂靡环N6 -剖分方式，如圖7 所示:在這種剖分方式下，立方體的左右兩面的剖分線方向相同，因此在橫向方向上的兩個(gè)相鄰剖分立方體之間是拓?fù)湟恢碌?。同理，在另外兩個(gè)方向上，兩個(gè)相鄰剖分立方體之間也是拓?fù)湟恢碌?。所以，只要相鄰的剖分立方體的剖分都使用如圖7 所示的方式，就可以避免它們之間的剖分二義性問(wèn)題。

3.2 剖分立方體與相鄰的未剖分立方體拓?fù)洳灰恢?/h3>

對(duì)于相鄰兩個(gè)立方體中一個(gè)是有二義性需要剖分成四面體，另一個(gè)是沒有二義性不需要剖分成四面體的情形?？梢詫]有二義性的立方體在每個(gè)面上作一條對(duì)角線，該對(duì)角線的方向需與圖7 中的一致。這樣，剖分立方體和未剖分立方體在其相接的表面上的拓?fù)涫峭耆恢碌?，也就不存在剖分二義性了，如圖8 所示:

圖8 兩個(gè)不同立方體的拓?fù)湟恢滦?/p>

需要注意的是，當(dāng)一個(gè)立方體在每個(gè)面上各增加一條對(duì)角線之后，立方體的邊由12 條增加到18條。因此，MC 算法原有的相交邊索引表不再適用，需要重新生成。同樣，原有的三角面片索引表也需要根據(jù)新的立方體結(jié)構(gòu)重新生成。

4 實(shí)驗(yàn)結(jié)果及分析

本文的實(shí)驗(yàn)環(huán)境為:CPU 為AMD Athlon II X4 640 Processor，800 MHz × 4;內(nèi) 存2G;顯 卡 為NVIDIA GeForce GTS 450;操作系統(tǒng)為32 位Ubuntu14.04;編譯器為gcc 4.6。

使用的CT 數(shù)據(jù)為:大腦，200 ×160 ×160，8位;腹部及骨盆，512 ×512 ×174，8 位。傳統(tǒng)的MC算法，MT 算法與改進(jìn)后的算法的三維重建效果如圖9，圖10 所示。

實(shí)驗(yàn)結(jié)果數(shù)據(jù)的比較見表1。

圖9 大腦表面三維重建圖像

圖10 腹部及骨盆部位三維重建圖像

表1 傳統(tǒng)算法與改進(jìn)算法的實(shí)驗(yàn)結(jié)果數(shù)據(jù)

對(duì)比圖9 右邊的局部放大圖可以看出，MC 算法生成了不正確的三角面片，而MT 算法和本文的改進(jìn)算法均生成了正確的三角面片。另外，MT 算法和本文改進(jìn)算法的重建圖像顯示了更多的細(xì)節(jié)，在重建精度上比MC 算法高。從表1 可知，本文改進(jìn)算法的執(zhí)行時(shí)間介于MC 算法和MT 算法之間，并且更接近于MC 算法。

5 結(jié)語(yǔ)

針對(duì)MC 算法二義性和MT 算法計(jì)算量大的問(wèn)題，提出了一種結(jié)合兩種算法優(yōu)點(diǎn)的改進(jìn)算法。該算法解決了MC 算法的二義性問(wèn)題，三維重建效果更接近于MT 算法，而計(jì)算量更接近于MC 算法，從整體上優(yōu)于這兩種算法。本文實(shí)現(xiàn)了傳統(tǒng)的MC 算法，MT 算法和本文提出的改進(jìn)算法，實(shí)驗(yàn)結(jié)果與預(yù)想的效果一致。

［1］Lerensen W E，Cline H E. Marching cubes:A high resolution 3D surface construction algorithm［J］. Computer Graphics，1987，21(3):163 -169.

［2］T S Newman，H Li.A survey of the marching cubes algorithm［J］.Computers Graphics，2006，30(5):854 -879.

［3］王曉，何建英.考慮MC 算法二義性后的三角化方法［J］.湖北工業(yè)大學(xué)學(xué)報(bào)，2007，22(6):65 -68.

［4］楊吉宏，薛凌燕，張民，等.基于Double Marching Cubes 的表面重建算法［J］.計(jì)算機(jī)工程與設(shè)計(jì)，2010，31(4):795 -797.

［5］A Gueziec，R Hummel.Exploiting Triangulated Surface Extraction using Tetrahedral Decomposition［J］.IEEE Transactions on Visualisation and Computer Graphics，1995(4):328 -342.

［6］Nielson G M，Hamann B. The asymptotic decider:Resolving the ambiguity in marching cubes［C］//Proceedings of Visualization ’91，San Diego:IEEE，1991:83 -91.

［7］Huiyan Jiang，Ye Zhang，Yudong Zhao.A medical image 3D reconstruction method based on improved MC and MT algorithms［J］.Industrial Electronics and Applications，2009(25 -27):3765 -3768.