賀江琳，劉世偉，王星月，岳晴，王遠(yuǎn)軍

上海理工大學(xué)醫(yī)學(xué)影像工程研究所，上海200093

前言

腦膠質(zhì)瘤通常呈彌漫性浸潤(rùn)生長(zhǎng)，范圍廣泛，邊界不清。受累區(qū)域的腦組織腫脹，溝變淺或消失，腦室變小。病變?cè)缙冢X膠質(zhì)瘤占位效應(yīng)常不明顯，中線(xiàn)結(jié)構(gòu)常沒(méi)有移位；病變中、晚期可表現(xiàn)出占位效應(yīng)，若病變偏一側(cè)，占位效應(yīng)征象可較早出現(xiàn)。腫瘤細(xì)胞多侵犯大腦半球兩個(gè)或兩個(gè)以上的部位，皮層及皮層下白質(zhì)均可受累［1-2］。

由于腦膠質(zhì)瘤的邊界模糊，其磁共振（Magnetic Resonance,MR）圖像分割有很大的難度。因此，為提高腫瘤分割的準(zhǔn)確性，研究人員提出上百種圖像分割的方法，文獻(xiàn)［3］詳細(xì)討論了學(xué)者們提出的經(jīng)典分割方法。近年來(lái)，國(guó)內(nèi)外研究人員關(guān)于腦腫瘤分割也提出不少新穎的方法，包括K均值聚類(lèi)和模糊C均值聚類(lèi)的分割方法［4］、基于梯度幅度的分水嶺分割法［5］、形態(tài)學(xué)閾值分割法［6-8］、基于Hough 變換定位與遺傳算法［9］、基于灰度分布的分割方法［10］、基于學(xué)習(xí)的方法［11-13］等。實(shí)踐證明，這些方法對(duì)于腦腫瘤的分割有明顯效果，在臨床上發(fā)揮了一定的作用，但是現(xiàn)在已有的任何一種分割算法在分割圖像時(shí)，都難以取得比較滿(mǎn)意的成果，結(jié)合多種算法的分割方法是目前研究的主要方向。1988年Kass［14］基于偏微分的方法提出的活動(dòng)輪廓模型對(duì)圖像分割領(lǐng)域產(chǎn)生了深遠(yuǎn)的影響，使得近年來(lái)，基于偏微分方程的圖像分割方法成為研究的熱點(diǎn)。后來(lái)隨著水平集方法的提出，幾何活動(dòng)輪廓模型與水平集方法相結(jié)合的曲線(xiàn)演化方法成為圖像分割領(lǐng)域的焦點(diǎn)，其中分割的快速、精確、自適應(yīng)和魯棒性成為研究注重的重點(diǎn)。

本研究主要對(duì)基于活動(dòng)輪廓模型的分割算法展開(kāi)研究與討論，并提出一種基于邊緣檢測(cè)的活動(dòng)輪廓（Geodesic Active Contours,GAC）模型與局部圖像擬合（Local Image Fitting,LIF）模型相結(jié)合的混合水平集分割方法。

1 基于水平集方法的分割模型

1.1 水平集函數(shù)

水平集的方法最早是由0sher 和Sethian［15］提出的，其基本思想是將平面閉合曲線(xiàn)隱含地表達(dá)成三維連續(xù)函數(shù)曲面的一個(gè)具有相同函數(shù)值的同值曲線(xiàn)，用連續(xù)演化曲線(xiàn)來(lái)表達(dá)目標(biāo)輪廓，并定義一個(gè)能量泛函使得其自變量包括曲線(xiàn)，從而實(shí)現(xiàn)將圖像分割過(guò)程轉(zhuǎn)換為求解能量泛函的最小值的過(guò)程。

給定一個(gè)常數(shù)c，定義水平集函數(shù)為?(x,y,t )，c，c稱(chēng)為水平集函數(shù)值，當(dāng)c=0 時(shí)，稱(chēng){(x,y)|?(x,y,t)}=0為零水平集。水平集函數(shù)常初始化為符號(hào)距離函數(shù)，即?(x,y)=±d(x,y)，其中，d(x,y)表示點(diǎn)(x,y)到輪廓線(xiàn)的距離，并且符號(hào)距離函數(shù)滿(mǎn)足 |??|=1。任意給定平面一條封閉曲線(xiàn)，相應(yīng)的符號(hào)距離函數(shù)如圖1所示。

圖1 符號(hào)距離函數(shù)Fig.1 Symbol distance function

1.2 C-V（Chan-Vese）模型

Mumford 和Shah［16］于1989年提出M-S 模型，該模型充分考慮了圖像數(shù)據(jù)、初始估計(jì)、目標(biāo)輪廓以及基于先驗(yàn)知識(shí)的約束等因素，在適當(dāng)初始化后可以自動(dòng)收斂，且對(duì)初始輪廓位置不敏感。但由于函數(shù)表達(dá)式本身的限制，該模型求取極小值比較困難，不能很好應(yīng)用于實(shí)際。據(jù)此，Chan 和Vese［17］針對(duì)M-S模型的缺點(diǎn)，進(jìn)行改進(jìn)和簡(jiǎn)化，提出基于區(qū)域信息的活動(dòng)輪廓模型：C-V模型。

C-V 模型假設(shè)待分割圖像I(x,y) 可以被閉合曲線(xiàn)C分為兩種同質(zhì)區(qū)域，即目標(biāo)區(qū)域以及背景區(qū)域。設(shè)這兩個(gè)區(qū)域的平均灰度值分別為c1和c2，該模型對(duì)應(yīng)的擬合能量函數(shù)如下：

其中，H(?) 為Heaviside 函數(shù)，Ω 為待分割灰度圖像I(x,y) 的定義域。當(dāng)閉合輪廓線(xiàn)處于兩個(gè)同質(zhì)區(qū)域的邊界時(shí)，ECV達(dá)到最小值。根據(jù)變分法和Euler-Lagrange方程，對(duì)式（1）求解最小值，可以得到水平集函數(shù)的演化方程為：

其中，c1和c2按如下方式進(jìn)行更新：

C-V模型充分利用圖像的全局灰度信息，具有全局特征，可以自動(dòng)檢測(cè)圖像的內(nèi)外邊界，并且有很好的抗噪能力。但在現(xiàn)實(shí)中，大量真實(shí)圖像的灰度分布是不均勻的，因此，傳統(tǒng)的C-V模型無(wú)法得到滿(mǎn)意的結(jié)果。

1.3 LIF模型

為了解決灰度不均勻圖像的分割問(wèn)題，Zhang等［18］提出LIF模型，該模型是一種嵌入圖像局部信息的區(qū)域活動(dòng)輪廓模型，利用分片光滑函數(shù)來(lái)近似待分割圖像I(x,y)。對(duì)于圖像區(qū)域Ω 中的一點(diǎn)(x,y)，引入水平集函數(shù)?(x,y)，LIF的能量泛函形式如下：

其中，ILFI(x,y) 為局部擬合項(xiàng)（Local Fitted Image,LFI），定義如下：

其中，m1、m2是圖像在點(diǎn)(x,y)處的兩個(gè)局部平均灰度值，即：

其中，Wk(x,y)為矩形窗口函數(shù)，通常采用截?cái)喔咚勾翱赪σ(x,y)，其大小為(4k+1)×(4k+1)，k為小于σ（標(biāo)準(zhǔn)差）的最大整數(shù)。通過(guò)變分法與最速下降法最小化式（4）可得相應(yīng)的水平集函數(shù)演化方程：

為了考慮水平集函數(shù)的正則性，獲得更好的分割效果，需要考慮重新初始化的問(wèn)題。Li 等［19］提出改進(jìn)的變分水平集方法，可以避免重新初始化，也可以通過(guò)高斯濾波正則水平集的方法，對(duì)水平集函數(shù)施以高斯濾波處理，能夠光滑水平集函數(shù)中的正則項(xiàng)，從而允許靈活地初始化水平集，實(shí)現(xiàn)快速分割［18,20］。因此，在LIF模型中，通常采用高斯濾波的方法來(lái)實(shí)現(xiàn)對(duì)水平集符號(hào)距離函數(shù)的正則化過(guò)程以及保證水平集演化的穩(wěn)定性。

當(dāng)水平集函數(shù)固定時(shí)，m1和m2的取值是一種與σ大小有關(guān)的局部量。這種局部特征使得該模型在分割灰度不均勻的圖像時(shí)可以取得較好的效果，其分割結(jié)果優(yōu)于C-V 模型。此外，相比于LBF 模型，LIF 模型的計(jì)算量大大縮減，提高計(jì)算效率，不足的是，LIF 模型對(duì)初始輪廓線(xiàn)的位置非常敏感，其分割結(jié)果依賴(lài)于初始輪廓線(xiàn)的位置、大小以及形狀。

1.4 GAC模型

GAC模型是在Snake模型的基礎(chǔ)上提出來(lái)的，該模型是一種典型的基于邊緣檢測(cè)的能量模型，其基本思想是利用黎曼空間中的測(cè)地線(xiàn)概念，把尋找圖像中的邊界線(xiàn)問(wèn)題轉(zhuǎn)化為尋找一條加權(quán)弧長(zhǎng)的最小值［21］。由于水平集的引入，GAC模型消除了Snake模型對(duì)自由參數(shù)的依賴(lài)，在曲線(xiàn)演化過(guò)程中能夠自動(dòng)處理拓?fù)浣Y(jié)構(gòu)的變化，并檢測(cè)到有明顯輪廓的物體邊界。

對(duì)于給定的待分割圖像I(x,y) ，GAC 模型的能量泛函定義為活動(dòng)輪廓線(xiàn)C的加權(quán)長(zhǎng)度，具體形式如下：

其中，L(C)表示閉合曲線(xiàn)的弧長(zhǎng)，s表示閉合曲線(xiàn)的弧長(zhǎng)參數(shù)，g表示邊緣檢測(cè)函數(shù)。引入水平集函數(shù)?(x,y,t)并代替曲線(xiàn)C，最小化式（8）可得相應(yīng)的梯度

下降流：

其中，k為曲率，c 是一個(gè)可選常數(shù)。c 的取值決定著式（9）中恒定速度項(xiàng)g(??) c 的運(yùn)動(dòng)方向：其方向總是恒定指向曲線(xiàn)的內(nèi)部（c ＞0）或外部（c ＜0），GAC 模型中添加此項(xiàng)有利于曲線(xiàn)在演化時(shí)，能夠在一定程度上緩解曲率為負(fù)的情況下對(duì)曲線(xiàn)運(yùn)動(dòng)帶來(lái)的阻礙作用，從而有利于對(duì)凹陷邊界的正確檢測(cè)。

相比于Snake模型，GAC模型不再依賴(lài)于曲線(xiàn)參數(shù)，并且可以自動(dòng)適應(yīng)曲線(xiàn)的拓?fù)浣Y(jié)構(gòu)變化，能夠有效分割邊界清晰的物體。但是GAC模型在分割弱邊緣或離散邊緣的物體時(shí)容易產(chǎn)生邊界泄漏的問(wèn)題，并且該方法對(duì)初始輪廓線(xiàn)的位置要求較高，一般選擇在待檢測(cè)物體的附近，否則分割結(jié)果將不太理想。此外，該模型中恒定速度項(xiàng)c 的取值必須為先驗(yàn)參數(shù)，c 值過(guò)大將導(dǎo)致曲線(xiàn)在圖像較為平滑區(qū)域產(chǎn)生錯(cuò)誤的分割結(jié)果，因此參數(shù)c 的取值需要慎重對(duì)待。

2 本文模型

由于圖像的邊界和區(qū)域的重要性，在分割過(guò)程中既要保證圖像區(qū)域信息的完整，也要保證圖像邊界信息不丟失，所以結(jié)合LIF模型與GAC模型，本文提出一種新的混合水平集模型。

本文的水平集能量模型是LIF模型與GAC模型的一個(gè)線(xiàn)性組合，具體形式如下：

其水平集演化偏微分方程為：

式中，α的取值決定著模型在分割時(shí)，LIF 模型與GAC模型各自發(fā)揮作用的大小。當(dāng)α=0 時(shí)，該模型就成為GAC 模型；當(dāng)α=1 時(shí)，該模型則變成LIF 模型。合理選擇α值的大小，本文模型可以改進(jìn)GAC模型中對(duì)弱邊緣以及離散邊緣產(chǎn)生的邊界泄漏的問(wèn)題；并且綜合運(yùn)用圖像的局部區(qū)域信息，發(fā)揮LIF 模型的優(yōu)點(diǎn)，減少模型的計(jì)算量，從而提高曲線(xiàn)演化的效率。此外，常數(shù)速度項(xiàng)cg( ?I)δε(?)里c 的取值一般為先驗(yàn)參數(shù)，但由于c 值的大小不易具體確定，若c 值過(guò)大將導(dǎo)致曲線(xiàn)在圖像較為平滑區(qū)域產(chǎn)生錯(cuò)誤的分割結(jié)果。因此，本文采用變系數(shù)V(I)來(lái)代替常數(shù)c［22］。V(I)的定義如下：

其中，β1與β2為各項(xiàng)的權(quán)重系數(shù)且β1=ω?β2，H為Heaviside函數(shù)，的定義如下：

采用變系數(shù)的方法可以提高模型檢測(cè)圖像凹陷邊界的能力，更好地保護(hù)弱邊界，抑制分割曲線(xiàn)滲透弱邊界的問(wèn)題，并且提高曲線(xiàn)的演化速度。變系數(shù)V(I)只改變?cè)?shù)c 的大小，不改變其正負(fù)。

另外，在傳統(tǒng)水平集方法中，為了保持符號(hào)距離函數(shù) ||??=1 的優(yōu)良性質(zhì)，我們有必要使水平集函數(shù)在演化的過(guò)程中始終近似于符號(hào)距離函數(shù)。早期往往采用周期性的重新初始化方法，此種方法極大地提高了模型的計(jì)算復(fù)雜度，而且重新初始化的操作時(shí)間和執(zhí)行頻率難以具體確定。因此，本文采用Zhang 等［18］提出的一種高斯濾波正則水平集函數(shù)法來(lái)實(shí)現(xiàn)水平集函數(shù)的正則化，即通過(guò)對(duì)水平集函數(shù)前一次迭代的結(jié)果進(jìn)行高斯濾波來(lái)實(shí)現(xiàn)正則化過(guò)程，并且為了簡(jiǎn)化符號(hào)距離函數(shù)，水平集函數(shù)初始化為如下形式：

其中，d 為大于零的常數(shù)，且d＞2ε是定義的正則的Dirac函數(shù)δε的寬度。

本文模型的算法可歸納如下：（1）采用C-V 模型預(yù)處理腦膠質(zhì)瘤MR圖像，提取腦組織；（2）根據(jù)經(jīng)驗(yàn)設(shè)置迭代次數(shù)n，參數(shù)的具體取值見(jiàn)本文第3 節(jié)；（3）初始化水平集函數(shù)?(x,y,t) ，具體形式如式（15）；（4）通過(guò)式（12）來(lái)演化水平集函數(shù)?，每次演化結(jié)束后對(duì)水平集函數(shù)?進(jìn)行高斯濾波；（5）若程序的迭代次數(shù)等于n，算法則停止，否則重復(fù)（4）。

實(shí)驗(yàn)過(guò)程中，由于MR圖像中非腦組織區(qū)域的部分灰度信息會(huì)對(duì)腦腫瘤分割產(chǎn)生干擾，因此本文在預(yù)處理時(shí)采用C-V模型與數(shù)學(xué)形態(tài)學(xué)相結(jié)合的方法提取腦組織，為后續(xù)的腦腫瘤分割實(shí)驗(yàn)消除一定的影響因素。此外，迭代次數(shù)n的取值十分重要，本文通過(guò)多次實(shí)驗(yàn)驗(yàn)證后，從中選擇最優(yōu)的數(shù)值作為標(biāo)準(zhǔn)的迭代次數(shù)。為避免初始輪廓線(xiàn)位置對(duì)分割實(shí)驗(yàn)結(jié)果帶來(lái)的影響，本文在實(shí)驗(yàn)中采用手動(dòng)設(shè)定初始輪廓線(xiàn)圓心位置的方式。

3 實(shí)驗(yàn)結(jié)果與分析

本文實(shí)驗(yàn)數(shù)據(jù)為一組格式為DICOM文件的MR圖像，其數(shù)據(jù)矩陣大小為512×512，數(shù)據(jù)類(lèi)型為16 位無(wú)符號(hào)整型，該序列共有248 幅MR 圖像，本文選取其中可以看到腦腫瘤的30幅MR圖像進(jìn)行實(shí)驗(yàn)。實(shí)驗(yàn)環(huán)境的操作系統(tǒng)為64 位的Windows 10，處理器為Inter(R) Core(TM) i3-2350M CPU(2.30 GHz)，內(nèi)存為8 GB；采用的軟件為MATLAB R2014a，本文腦腫瘤分割算法的所有代碼均在此平臺(tái)上運(yùn)行。

為定量評(píng)價(jià)分割結(jié)果，本文使用有經(jīng)驗(yàn)的臨床醫(yī)生手動(dòng)分割結(jié)果作為評(píng)價(jià)分割準(zhǔn)確性的標(biāo)準(zhǔn)。本文使用的定量評(píng)價(jià)指數(shù)分別為Jc（Jaccard Index,Jc）指數(shù)與Dice相似性系數(shù)（Dice Similarity Index,DSI）。

Jc 指數(shù)的定義為：

其中，Rdoc為臨床醫(yī)生手動(dòng)分割的結(jié)果（評(píng)價(jià)標(biāo)準(zhǔn)），Rpro為算法分割的結(jié)果（評(píng)價(jià)對(duì)象）。算法分割結(jié)果越準(zhǔn)確，則Jc 值越接近1。

DSI定義為：

式中，Rdoc與Rpro表示的意義與式（16）一樣，算法分割結(jié)果與評(píng)價(jià)標(biāo)準(zhǔn)重合度越高，則DSI越接近1，表示分割結(jié)果越精確。

在實(shí)驗(yàn)過(guò)程中，本文依據(jù)式（12）和式（13）演化水平集函數(shù)，相應(yīng)參數(shù)的取值為：α=0.8，ε=1.5，β2=0.28，ω=0.1 ，迭代的步長(zhǎng)Δt=0.1，迭代次數(shù)n=200。此外，用于計(jì)算LIF 模型中局部均值的截?cái)喔咚勾翱赪(x,y)σ的σ取值為3，k為1，即采用窗口大小為5×5的高斯窗口。除參數(shù)α、迭代次數(shù)n與迭代步長(zhǎng)Δt，其余參數(shù)取值均保持不變的情況下，分別采用本文的模型（α=0.8，n=200）、LIF 模型（α=1，n=600）與GAC 模型（α=0，n=100，Δt=3.0，式（12）中常數(shù)c=0.2）對(duì)預(yù)處理后的MR 圖像進(jìn)行腦腫瘤分割。本文主要展示實(shí)驗(yàn)中4 幅腦腫瘤MR 圖像的分割結(jié)果，圖像相關(guān)實(shí)驗(yàn)結(jié)果如圖2所示。

由以上實(shí)驗(yàn)結(jié)果可知，在分割同一幅圖像時(shí)，LIF模型需要的迭代次數(shù)大于本文模型，且分割結(jié)果稍差于本文模型，在圖像弱邊緣處無(wú)法得到精確的分割結(jié)果；而GAC 模型分割圖像時(shí)，速度較慢，需要增加迭代的步長(zhǎng)，從而減少迭代次數(shù)，并且在分割凹陷邊界時(shí)，僅采用GAC 模型往往不能得到正確的分割結(jié)果。本文模型有效地結(jié)合GAC模型與LIF模型的優(yōu)點(diǎn)，在曲線(xiàn)演化過(guò)程中，可以減少算法的計(jì)算量，改善LIF 模型與GAC 模型在弱邊緣處產(chǎn)生的邊界泄漏問(wèn)題，提高算法的分割精度。因此，在相同迭代次數(shù)的情況下，本文模型的分割結(jié)果優(yōu)于GAC 模型與LIF模型。

圖2 3種模型分割結(jié)果圖Fig.2 Segmentation results of 3 models

本文采用醫(yī)生的手動(dòng)分割結(jié)果作為評(píng)價(jià)標(biāo)準(zhǔn)。對(duì)本文的實(shí)驗(yàn)分割結(jié)果做定量驗(yàn)證分析，相關(guān)數(shù)據(jù)如表1所示。從表1可以看出，在兩個(gè)定量評(píng)價(jià)指標(biāo)中，本文模型的Jc 的平均值為91.00%，DSI的平均指標(biāo)為95.27%。兩者的數(shù)據(jù)皆接近于1，高于GAC 模型與LIF模型的定量評(píng)價(jià)指標(biāo)，說(shuō)明本文的分割結(jié)果較為精確，即本文模型在分割圖像時(shí)具有一定的有效性。

4 結(jié)論

本文基于混合水平集算法分割腦腫瘤展開(kāi)研究與探討。首先，對(duì)傳統(tǒng)的腦腫瘤分割方法進(jìn)行簡(jiǎn)單的介紹。然后，以水平集算法為核心，分別介紹C-V模型、GAC模型與LIF模型，并將GAC模型與LIF模型的各能量項(xiàng)按一定加權(quán)比例結(jié)合起來(lái)，實(shí)現(xiàn)混合水平集算法。通過(guò)實(shí)驗(yàn)，我們得到了較好的分割結(jié)果，實(shí)驗(yàn)中的不足之處在于混合水平集模型對(duì)初始輪廓線(xiàn)的位置較敏感，需要手動(dòng)確定其圓心位置，在以后的研究中，可以對(duì)此做進(jìn)一步的改進(jìn)。

表1 本文算法分割精度的定量評(píng)價(jià)（%）Tab.1 Quantitative evaluation of the segmentation accuracy(%)