基于混合水平集的腦膠質瘤分割方法

賀江琳，劉世偉，王星月，岳晴，王遠軍

上海理工大學醫學影像工程研究所，上海200093

前言

腦膠質瘤通常呈彌漫性浸潤生長，范圍廣泛，邊界不清。受累區域的腦組織腫脹，溝變淺或消失，腦室變小。病變早期，腦膠質瘤占位效應常不明顯，中線結構常沒有移位；病變中、晚期可表現出占位效應，若病變偏一側，占位效應征象可較早出現。腫瘤細胞多侵犯大腦半球兩個或兩個以上的部位，皮層及皮層下白質均可受累［1-2］。

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

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

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

1.1 水平集函數

水平集的方法最早是由0sher 和Sethian［15］提出的，其基本思想是將平面閉合曲線隱含地表達成三維連續函數曲面的一個具有相同函數值的同值曲線，用連續演化曲線來表達目標輪廓，并定義一個能量泛函使得其自變量包括曲線，從而實現將圖像分割過程轉換為求解能量泛函的最小值的過程。

給定一個常數c，定義水平集函數為?(x,y,t )，c，c稱為水平集函數值，當c=0 時，稱{(x,y)|?(x,y,t)}=0為零水平集。水平集函數常初始化為符號距離函數，即?(x,y)=±d(x,y)，其中，d(x,y)表示點(x,y)到輪廓線的距離，并且符號距離函數滿足 |??|=1。任意給定平面一條封閉曲線，相應的符號距離函數如圖1所示。

圖1 符號距離函數Fig.1 Symbol distance function

1.2 C-V（Chan-Vese）模型

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

C-V 模型假設待分割圖像I(x,y) 可以被閉合曲線C分為兩種同質區域，即目標區域以及背景區域。設這兩個區域的平均灰度值分別為c1和c2，該模型對應的擬合能量函數如下：

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

其中，c1和c2按如下方式進行更新：

C-V模型充分利用圖像的全局灰度信息，具有全局特征，可以自動檢測圖像的內外邊界，并且有很好的抗噪能力。但在現實中，大量真實圖像的灰度分布是不均勻的，因此，傳統的C-V模型無法得到滿意的結果。

1.3 LIF模型

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

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

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

其中，Wk(x,y)為矩形窗口函數，通常采用截斷高斯窗口Wσ(x,y)，其大小為(4k+1)×(4k+1)，k為小于σ（標準差）的最大整數。通過變分法與最速下降法最小化式（4）可得相應的水平集函數演化方程：

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

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

1.4 GAC模型

GAC模型是在Snake模型的基礎上提出來的，該模型是一種典型的基于邊緣檢測的能量模型，其基本思想是利用黎曼空間中的測地線概念，把尋找圖像中的邊界線問題轉化為尋找一條加權弧長的最小值［21］。由于水平集的引入，GAC模型消除了Snake模型對自由參數的依賴，在曲線演化過程中能夠自動處理拓撲結構的變化，并檢測到有明顯輪廓的物體邊界。

對于給定的待分割圖像I(x,y) ，GAC 模型的能量泛函定義為活動輪廓線C的加權長度，具體形式如下：

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

下降流：

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

相比于Snake模型，GAC模型不再依賴于曲線參數，并且可以自動適應曲線的拓撲結構變化，能夠有效分割邊界清晰的物體。但是GAC模型在分割弱邊緣或離散邊緣的物體時容易產生邊界泄漏的問題，并且該方法對初始輪廓線的位置要求較高，一般選擇在待檢測物體的附近，否則分割結果將不太理想。此外，該模型中恒定速度項c 的取值必須為先驗參數，c 值過大將導致曲線在圖像較為平滑區域產生錯誤的分割結果，因此參數c 的取值需要慎重對待。

2 本文模型

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

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

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

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

其中，β1與β2為各項的權重系數且β1=ω?β2，H為Heaviside函數，的定義如下：

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

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

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

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

實驗過程中，由于MR圖像中非腦組織區域的部分灰度信息會對腦腫瘤分割產生干擾，因此本文在預處理時采用C-V模型與數學形態學相結合的方法提取腦組織，為后續的腦腫瘤分割實驗消除一定的影響因素。此外，迭代次數n的取值十分重要，本文通過多次實驗驗證后，從中選擇最優的數值作為標準的迭代次數。為避免初始輪廓線位置對分割實驗結果帶來的影響，本文在實驗中采用手動設定初始輪廓線圓心位置的方式。

3 實驗結果與分析

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

為定量評價分割結果，本文使用有經驗的臨床醫生手動分割結果作為評價分割準確性的標準。本文使用的定量評價指數分別為Jc（Jaccard Index,Jc）指數與Dice相似性系數（Dice Similarity Index,DSI）。

Jc 指數的定義為：

其中，Rdoc為臨床醫生手動分割的結果（評價標準），Rpro為算法分割的結果（評價對象）。算法分割結果越準確，則Jc 值越接近1。

DSI定義為：

式中，Rdoc與Rpro表示的意義與式（16）一樣，算法分割結果與評價標準重合度越高，則DSI越接近1，表示分割結果越精確。

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

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

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

本文采用醫生的手動分割結果作為評價標準。對本文的實驗分割結果做定量驗證分析，相關數據如表1所示。從表1可以看出，在兩個定量評價指標中，本文模型的Jc 的平均值為91.00%，DSI的平均指標為95.27%。兩者的數據皆接近于1，高于GAC 模型與LIF模型的定量評價指標，說明本文的分割結果較為精確，即本文模型在分割圖像時具有一定的有效性。

4 結論

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

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