孫海鵬， 余偉巍， 席 平

?

基于Level Set的交互式快速分割算法

孫海鵬， 余偉巍， 席 平

（北京航空航天大學(xué)機(jī)械工程及自動化學(xué)院，北京100191）

三維醫(yī)學(xué)圖像數(shù)據(jù)量大，并且受噪聲、邊界模糊等原因的影響，致使三維分割過程消耗時間較長，容易產(chǎn)生欠分割或過度分割。針對以上問題，提出一種基于Level Set的三維快速分割算法，采用Fast Marching獲取二維分割區(qū)域，優(yōu)化輪廓邊界，利用直線數(shù)值微分算法（Digital Differential Analyzer，DDA）提取輪廓像素；進(jìn)一步引入掃描線種子填充思想，實現(xiàn)醫(yī)學(xué)圖像的三維快速分割。實驗結(jié)果表明，上述算法能夠快速準(zhǔn)確地分割出感興趣區(qū)域。

計算機(jī)應(yīng)用；醫(yī)學(xué)圖像三維分割；Level Set算法；數(shù)值微分算法；掃描線種子填充

醫(yī)學(xué)圖像分割技術(shù)，是圖像分割領(lǐng)域中一個重要分支，自上世紀(jì)90年代起一直受到學(xué)術(shù)界的廣泛重視，是虛擬手術(shù)系統(tǒng)中不可缺少的一個重要組成部分，為結(jié)構(gòu)分析、運動分析和三維可視化等提供基礎(chǔ)數(shù)據(jù)來源，分割結(jié)果的有效性直接影響最終虛擬手術(shù)結(jié)果的有效性。

目前主流醫(yī)學(xué)圖像分割算法大多針對二維醫(yī)學(xué)圖像，如：M Kass提出的Snakes算法，它是用能量最小化作為框架，通過定義內(nèi)能和外能來模擬物理上的力學(xué)原理最終實現(xiàn)醫(yī)學(xué)圖像的分割，但snake算法受能量函數(shù)的初始狀態(tài)影響較大；Osher和Sethian提出的Level Set算法，它將曲線演化問題轉(zhuǎn)化為偏微分方程數(shù)值求解問題，具有很強(qiáng)的處理拓?fù)浣Y(jié)構(gòu)變化的能力，但Level Set的計算復(fù)雜度較高。針對Level Set算法計算復(fù)雜度較高的問題，提出了快速水平集算法(Fast Marching Level Set Method)；此外，還有一些基于結(jié)合特定理論的分割方法，如：基于數(shù)學(xué)形態(tài)的圖像分割，基于模糊技術(shù)的圖像分割，基于神經(jīng)網(wǎng)絡(luò)的圖像分割，基于知識的圖像分割等。但為了實現(xiàn)構(gòu)造三維醫(yī)學(xué)圖像需要多次進(jìn)行二維醫(yī)學(xué)圖像分割，比較耗時，并且每次二維醫(yī)學(xué)圖像分割結(jié)果都極大地影響到三維醫(yī)學(xué)圖像的精確性。

然而，現(xiàn)有的直接針對三維體數(shù)據(jù)分割，例如基于6聯(lián)通區(qū)域增長法，算法耗時長，且容易出現(xiàn)欠分割或過度分割等情況。

針對以上問題，本文在分析和融合多種分割算法的基礎(chǔ)上，面向三維醫(yī)學(xué)圖像分割，提出一種基于Level Set的三維快速分割算法，采用Fast Marching 快速步進(jìn)獲取二維分割區(qū)域，優(yōu)化輪廓邊界，利用直線掃描數(shù)值微分算法（DDA）方式提取輪廓像素，引入掃描線種子填充思想，實現(xiàn)醫(yī)學(xué)圖像的三維快速分割，有效解決了三維醫(yī)學(xué)圖像分割過程中所耗時間長、易產(chǎn)生欠分割和過度分割的問題。

1 Level Set算法

Level Set算法的基本思想是將平面閉合曲線隱含地表達(dá)為二維曲面函數(shù)的水平集，即具有相同函數(shù)值的點集，通過水平集函數(shù)曲面的進(jìn)化隱含地求解曲線的運動。盡管這種轉(zhuǎn)化使得問題在形式上變得復(fù)雜，但在問題求解上帶來很多優(yōu)點，其最大的優(yōu)點在于曲線的拓?fù)渥兓軌虻玫胶茏匀坏靥幚?，而且可以獲得唯一的滿足熵條件的弱解。

水平集函數(shù)的演化滿足如下基本方程

Φ(x, t)為水平集函數(shù)，Φ(x, t=0)為初始設(shè)置的演化曲線，其零水平集表示目標(biāo)輪廓線，即

（2）

Level Set算法中最典型的一種算法是FastMarching算法。Fast Marching算法只考慮一種界面運動的特殊情況，即界面的運動速度>0。假定是界面經(jīng)過一個指定點(x, y)的時間，這樣，就滿足如下的方程

式(3)簡單地說明了到達(dá)時間的梯度和界面的運動速度成反比。廣義上說，有兩種方法可以用來求解近似運動界面隨時間變化的位置：一種是通過迭代和數(shù)字近似來獲取式(1)中的微商；另一種是構(gòu)建式(3)中到達(dá)時間T的解決方案。而快速水平集算法依賴于后一種方法。

要得到式(3)中的到達(dá)時間，等價于求解下面的二次方程

（4）

由于式(4)差分法的結(jié)構(gòu)決定演化曲線傳遞的信息是單向的，也就是時間由小到大的傳遞過程。因此，求解式(4)采取從最小的時間向外求解過程。

Fast Marching算法的過程分為初始化和循環(huán)。

初始化：

（1）活動點 是指所有網(wǎng)格點中時間固定的點。在此算法中即用戶指定的種子點P，表示時間T(x, y)=0，如圖2黑色方格所示；

（2）窄帶點 所有在窄帶中的點叫做窄帶點。在此算法中就是所有種子點的4-鄰接的點，時間T(x, y)=1/F(x, y)，如圖2灰色方格所示；

（3）遠(yuǎn)離點 除了活動點和窄帶點以外，所有其他的網(wǎng)格點都是遠(yuǎn)離點，T(x, y)=TIME_ MAX，如圖2白色方格所示。

(a)?????(b) ?????(c)

循環(huán)：

（1）開始循環(huán)，假定點P是窄帶中具有最小時間T的點；

（2）標(biāo)定點P為活動點，并從窄帶中刪除；

（3）標(biāo)定點P的4-聯(lián)通點：如果點P的鄰接點為活動點，則不改變時間；如果其鄰接點為窄帶點，則按照公式(4)更新鄰接點的時間；如果其鄰接點為遠(yuǎn)離點，則標(biāo)定該鄰接點為窄帶點，同時按照公式(4)更新該鄰接點的時間；

（4）如果某一點的到達(dá)時間超過指定的限值，則循環(huán)結(jié)束，否則跳到(1)。

依據(jù)Fast Marching算法對關(guān)鍵層圖像進(jìn)行分割，分割結(jié)果如圖3所示。

(a) 所選關(guān)鍵層醫(yī)學(xué)圖像?????(b) 關(guān)鍵層分割結(jié)果

2 輪廓追蹤與輪廓精簡

為了得到圖像的外部特征，關(guān)鍵層圖像在Level Set算法分割過程中進(jìn)行了二值化處理，得到目標(biāo)輪廓的單值區(qū)域，通過對該單值區(qū)域的輪廓追蹤，得到目標(biāo)輪廓的外部輪廓點。

2.1 輪廓追蹤

本文提出的輪廓追蹤方法，其基本原理如下：

首先將輪廓上最左下方的點作為輪廓搜索的起始點。其次按照下述的“探測準(zhǔn)則”來尋找目標(biāo)輪廓上的其它像素。

探測準(zhǔn)則：從第1個邊界點開始，定義初始的搜索方向為左上方；如果左上方的點是黑點，則為這個邊界點的新的邊界點，否則，搜索方向順時針旋轉(zhuǎn)45°。直到找到第1個黑點為止。之后把這個黑點作為新的邊界點，在當(dāng)前方向上逆時針旋轉(zhuǎn)90°，繼續(xù)同樣的方法搜索下一個黑點,直到返回最初的邊界點。如圖4所示。

追蹤算法對圖3(b)進(jìn)行追蹤結(jié)果如圖5所示。

圖4 輪廓追蹤法示意圖 (箭頭代表搜索方向)

圖5 輪廓追蹤結(jié)果

2.2 輪廓精簡

在輪廓追蹤過程中，目標(biāo)輪廓的每一個像素點都被存儲，因此造成較大的數(shù)據(jù)冗余。為方便以后輪廓曲線的調(diào)整，必須對追蹤結(jié)果進(jìn)行精簡。目前提出的輪廓精簡算法主要有等距采樣法和曲率采樣法兩種，但是前者會導(dǎo)致大量特征點的丟失，而后者因為曲率計算公式會涉及2 級導(dǎo)數(shù)的計算，比較復(fù)雜。因此，本文采用“弦差法”進(jìn)行輪廓精簡，“弦差法”的原理如圖6所示。

圖6 “弦差法”原理

首先定義距離閾值ΔD描述輪廓精簡的精度，然后取a作為輪廓精簡的起點，a為終點，計算a到aa的距離D，如果D<ΔD，說明點a在控制弦差范圍之內(nèi)，終點加1，變?yōu)閍，計算點a, a到弦aa的距離，比較D，D與距離閾值ΔD之間的大小，依次循環(huán)，直到起點a和終點a之間的點到弦aa的距離D>ΔD為止，此時說明a和a之間的點已經(jīng)超出了精度控制范圍，因此可將a作為輪廓精簡點予以保留，而a和a之間的點可以舍棄，把a(bǔ)作為起始點，重復(fù)以上的步驟，尋找下一精簡點。

依次循環(huán)，直到遍歷所有輪廓點，精簡結(jié)束。圖7是運用上述方法取距離閾值ΔD= 0.6時對圖5進(jìn)行輪廓精簡后的輸出圖像，精簡前輪廓點數(shù)為429，精簡后輪廓點數(shù)為35，從精簡結(jié)果可以看出，精簡后圖像與精簡前圖像形狀極其一致，而點的數(shù)據(jù)量大大減少，因此精簡效果十分明顯。

圖7 輪廓精簡結(jié)果

3 構(gòu)造輪廓曲面

3.1 生成Bézier輪廓曲線和投影

3.1.1 Bézier曲線

Bézier曲線是1962年貝齊爾提出的一種參數(shù)多項式類型的構(gòu)造曲線的方法。

由于處理作為頂點的絕對矢量比Bézier多邊形邊的相對矢量方便；故本文使用英國的弗里斯特改寫B(tài)ézier曲線方程如下

（6）

3.1.2 構(gòu)造Béizer輪廓曲線及投影

如直接利用在2.2節(jié)中的輪廓精簡方法得到的一系列的輪廓點構(gòu)造Bézier曲線會造成曲線次數(shù)較高，調(diào)整輪廓形狀不方便等不良效果。因此，將精簡后的兩兩相鄰的輪廓點作為一段Bézier曲線的起始控制頂點和終端控制頂點，分段構(gòu)造輪廓Bézier曲線。

兩個控制頂點僅僅可以構(gòu)造1次Bézier曲線，為了得到更好的光順性，需人為地在始末控制頂點間添加兩個控制頂點，滿足構(gòu)造3次Bézier曲線的要求。添加控制頂點及構(gòu)造Bézier輪廓曲線過程如圖8所示。

圖8 添加Bézier曲線控制頂點示意圖

其中、、、、、為精簡后的控制頂點，取與的連線，在連線上取點，使點滿足以下條件

取與的連線，在連線上取點，使點滿足以下條件

（8）

再取與的連線，在連線上取點，使點滿足以下條件

取與的連線，在連線上取點，使點滿足以下條件

（10）

以此類推可以人為地在精簡后的兩個相鄰控制頂點間添加兩個新控制頂點，以如(,,,)，(,,,)…(,,,)等為控制頂點分段構(gòu)造3次Bézier曲線。

如圖9所示為構(gòu)造的3次Bézier樣條曲線示意圖

圖9 構(gòu)造3次Bézier樣條曲線示意圖

圖10為運用上述構(gòu)造輪廓曲線方法對圖7精簡輪廓后的控制頂點構(gòu)造Bézier輪廓曲線的結(jié)果圖。

圖10 構(gòu)造Bézier輪廓曲線

由于待分割目標(biāo)輪廓在某截面方向上輪廓大體相似，故將關(guān)鍵層的輪廓曲線沿與關(guān)鍵層面相垂直的軸向投影便可得到待分割目標(biāo)的大致體輪廓。本文以橫切面和方向為例，以CT切片組的方向的切片間距為增量，將之前得到的Bézier輪廓曲線沿方向在整個CT切片組空間中進(jìn)行投影，得到一組Bézier輪廓曲線族。如圖11所示。

圖11 Bézier輪廓曲線投影

3.2 DDA輪廓插值及構(gòu)造輪廓曲面

雖然各截面上目標(biāo)輪廓大體相似，但每層截面上目標(biāo)輪廓存在位置偏差，需調(diào)整各輪廓曲線的位置，完全包括待分割目標(biāo)。移動Bézier曲線的控制頂點調(diào)整此前生成的Bézier輪廓曲線族到待分割區(qū)域，如圖12所示。

圖12 Bézier曲線族調(diào)整圖

為了滿足體分割中空間封閉的條件，需要得到各輪廓曲線上所有的點在空間中的位置坐標(biāo)以此來構(gòu)造空間封閉輪廓曲面。但由于在反求Bézier曲線上點坐標(biāo)存在著一定的困難和計算量較大等因素，現(xiàn)簡化為通過對每條Bézier輪廓曲線上各控制頂點進(jìn)行DDA插值的方法來得到近似的輪廓曲線上各點在空間中的位置坐標(biāo)。

將得到的每條輪廓線上的點按照下述方法構(gòu)造成空間封閉輪廓曲面。

（1）設(shè)定垂直于關(guān)鍵層面的軸向的正負(fù)方向，對輪廓曲線族中各輪廓曲線進(jìn)行編號，使各輪廓曲線編號沿軸向正向依次增加；對每條輪廓曲線上的點沿著逆時針方向進(jìn)行編號。

（2）從第一層輪廓曲線開始，沿點序號遞增方向，取當(dāng)前輪廓曲線中相鄰的兩個輪廓點、以及上層輪廓曲線中與點的編號相同的輪廓點。將、、三點連成封閉三角形。取與點同層且點序號增加1的點，將、、三點連成封閉三角形。

（3）依次類推，直到循環(huán)到最上面一層輪廓曲線終止。

根據(jù)封閉的輪廓曲面對待分割目標(biāo)數(shù)據(jù)進(jìn)行輪廓曲面內(nèi)外標(biāo)定，本文將輪廓曲面外的體數(shù)據(jù)點的灰度值置為-3000。

如圖13所示為根據(jù)調(diào)整后的輪廓曲線族構(gòu)造的輪廓曲面以及區(qū)域內(nèi)外判斷。

圖13 輪廓曲面構(gòu)造圖

4 掃描線種子填充思想

體素是基本圖元，是組成體數(shù)據(jù)的基本單位。體素又分為表面體素和內(nèi)部體素。如果定義一個至少與其他體素6鄰接或位于體數(shù)據(jù)邊緣的體素為其表面體素，那么其表面體素可以組成一個18連通的封閉表面，同時它也可以表示為內(nèi)部體素6連通的物體。

掃描線種子填充思想正是根據(jù)內(nèi)部體素的6連通屬性，從未被填充鄰居區(qū)域中，選擇新的種子，利用遞歸函數(shù)來實現(xiàn)。由于體數(shù)據(jù)的數(shù)據(jù)量大，而且遞歸函數(shù)在實現(xiàn)時要消耗大量的內(nèi)存空間和花費大量時間。1985年Fishkin和Basky 通過分配給較短的像素跨度較高的優(yōu)先級的方法提高效率。1995年Albert等為了減少內(nèi)存的消耗和改善效率，提出了沿著某個方向連續(xù)填充體素的方法來代替僅僅填充相鄰的一個體素的方法。1998年Feng和Soon從4鄰域跨度并沿著一定的方向和堆疊一個或者多個新種子來填充一個跨度的體素。

但以上方法依然有兩個缺點。首先，即使只有一個種子需要填充，但是不止一個的同樣跨度的種子會被壓入堆棧。其次，即使填充的區(qū)域沒有包含新的種子點，但是尋找種子過程依然會在填充和未被填充的區(qū)域中執(zhí)行的。

為了避免以上兩個問題的發(fā)生，本文使用改進(jìn)后的三維種子填充算法來實現(xiàn)分割。

對填充的過程進(jìn)行記錄就可以阻止多余種子的填充。即在尋找新種子點之前，需要比較相鄰行和當(dāng)前填充跨度的范圍，之后根據(jù)比較的結(jié)果，判斷是否進(jìn)行新種子尋找。

相鄰行與當(dāng)前填充跨度的范圍的比較可以被分作下述5種情況，如圖14所示：

第1種情況 沒有相鄰的跨度與當(dāng)前正在填充的跨度重疊。那么在相鄰行的范圍內(nèi)尋找新的種子。

第2種情況 相鄰的跨度包含了當(dāng)前正在填充的跨度。不需要尋找新的種子點。

第3種情況 相鄰的跨度的左半部份與當(dāng)前正在填充的跨度相重疊。那么僅僅需要在相鄰的跨度的右半部分尋找新的種子點。

第4種情況 相鄰的跨度的右半部分與當(dāng)前正在填充的跨度相重疊。那么僅僅需要在相鄰的跨度的左半部份尋找新的種子點。

第5種情況 相鄰的跨度完全包含在當(dāng)前正在填充的跨度中。那么僅僅需要在沒有重疊的相鄰的跨度內(nèi)尋找新的種子點。

圖14 相鄰行與當(dāng)前填充跨度范圍比較分類圖

5 實驗結(jié)果分析

本文在INTEL的酷睿2雙核1.66G，2G內(nèi)存筆記本上進(jìn)行三維分割算法實驗，實驗數(shù)據(jù)為Ⅰ:512×512×24和Ⅱ:512×512×57的兩組CT切片圖像。圖15為Ⅰ、Ⅱ組圖像的未分割組織圖，圖16為對兩組數(shù)據(jù)使用Level Set算法三維分割結(jié)果，圖17為本文算法三維分割結(jié)果。

(Ⅰ) ?????(Ⅱ)

(Ⅰ) ?????(Ⅱ)

(Ⅰ) ?????(Ⅱ)

表1為本文算法和Level Set算法分割的比較結(jié)果?？杀葏?shù)為分割時間。

表1 本文算法與Level Set算法三維分割結(jié)果比較

由于三維種子填充算法存在回溯性，所以僅使用三維種子填充算法不能進(jìn)行相連通部分的分割，否則會出現(xiàn)欠分割的情況。本文算法對（Ⅰ）組數(shù)據(jù)的三維分割時間是利用Level Set算法三維分割的時間13.4%，對（Ⅱ）組數(shù)據(jù)的三維分割時間是利用Level Set算法三維分割的時間15.8%，大大提高了分割效率。

6 總 結(jié)

本文提出一種新型的醫(yī)學(xué)圖像分割算法。該算法第一，利用Level Set算法對關(guān)鍵層圖像（該層圖像具有一組醫(yī)學(xué)圖像中與待分割組織的輪廓大體相似的特性）進(jìn)行分割；第二，通過輪廓追蹤得到分割后的關(guān)鍵層的外輪廓點，但由于在輪廓追蹤過程中，目標(biāo)輪廓的每一個像素點都被存儲，造成數(shù)據(jù)冗余，不利于進(jìn)行輪廓曲線投影后的輪廓調(diào)節(jié)。所以對輪廓點使用“弦差法”進(jìn)行輪廓精簡；第三，將輪廓點作為控制頂點生成Bézier曲線作為新的輪廓線；第四，將Bézier輪廓曲線進(jìn)行投影，得到一族輪廓曲線；第五，調(diào)整每條輪廓曲線上的控制頂點使其達(dá)到恰當(dāng)位置，利用DDA對調(diào)整到位的控制頂點進(jìn)行插值，生成新的封閉的輪廓邊界，利用輪廓曲線族構(gòu)造輪廓曲面；第六，在封閉的輪廓邊界內(nèi)運用掃面線種子填充思想進(jìn)行三維分割；最后完成分割任務(wù)。

[1] KASS M, WITKIN A, TERZOPOULOS D. Snakes: active contour models [J]. International Journal of Computer Vision, 1987, 1(4): 321-331.

[2] Malladi R, Sethian J A, Vemuri B. Shape modeling with front propagation: a level set approach [J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 1995, 17(2): 158-174.

[3] Sethian J A. A fast marching level set method for monotonically advancing fronts [C]//Proceedings of the National Academy of Sciences, 1996: 9389-9392.

[4] 張海波. 醫(yī)學(xué)CT圖像的三維分割與可視化研究[D]. 濟(jì)南: 山東師范大學(xué), 2005.

[5] Feng L, Soon S H. An effective 3D seed fill algorithm [J]. Comput Graph, 1998, 22(5): 641-644.

[6] 郭開波, 周 鋼. 一種基于斷層測量圖片的實體輪廓提取算法[J]. 計算機(jī)輔助工程, 2001, (4): 50-54.

[7] 施法中. 計算機(jī)輔助幾何設(shè)計與非均勻有理B樣條[M]. 北京: 高等教育出版社, 2001. 115-165.

The Three-Dimensional Fast Segmentation Algorithm Based on Level Set Method

SUN Hai-peng, YU Wei-wei, XI Ping

( School of Mechanical Engineering and Automation, Beijing University of Aeronautics and Astronautics, Beijing 100191, China )

Because of the large volume of medical image data, the impact of noise, blurred boundaries and other reasons, the three-dimensional segmentation process is time-consuming, and easily produces less or over segmentation. To solve the above problems, this paper proposes a three-dimensional fast segmentation algorithm based on Level Set, using Level Set Fast Marching Method to obtain two-dimensional segmental region, optimizing the boundary contour, using the Digital Differential Analyzer method to extract contour pixels, finally introducing the idea of the Scan Line Seed-filling to achieve the three-dimensional fast segmentation. The actual clinical CT images of vertebral segmentation experiment result shows that this method can quickly and accurately separate out the interested area.

computer application;three-dimensional medical image segmentation; Level Set method; DDA; scan line seed-filling

TP 391.41

A

1003-0158(2011)03-0045-07

2009-11-02

孫海鵬（1985-），男，內(nèi)蒙古多倫人，碩士研究生，主要研究方向為計算機(jī)可視化、醫(yī)學(xué)圖像處理。