基于改進的雙水平集的MRI圖像快速分割方法

朱家明，李祥健，徐婷宜

(揚州大學 信息工程學院，江蘇 揚州 225127)

0 引言

隨著科技的不斷進步，醫學影像成為醫生后期診斷的重要依據。常用于醫學領域的影像學成像技術有：計算機斷層掃描成像技術(CT)、核磁共振成像(MRI)、單光子輻射斷層成像(SPECT)、正電子放射斷層成像技術(PET)、光學斷層纖維成像、核磁共振顯微成像以及正電子斷層顯微成像等。其中MRI圖像由于對含水分的軟組織比較敏感，可以用于對腦部等組織進行解剖結構分析檢測的場合[1]。然而MRI圖像的成像質量也存在不足：受到設備、環境以及射頻信號不均勻的影響，圖像經常存在噪聲和灰度不均的現象[2]。

醫學圖像分割是醫學圖像處理的一個關鍵技術[3]，是服務于臨床醫生的計算機輔助診斷的重要組成部分。其任務是從醫學圖像中提取感興趣的目標，即解剖組織。國內外學者提出了很多圖像分割方法，如閾值法、區域生長法；結合特定理論工具的方法，如模式識別法、人工神經網絡法等；基于模糊分割的方法，如小波變換法；基于遺傳算法的方法以及水平集方法等[4]。其中水平集方法是解決曲線演化問題的一種經典方法[5-7]。水平集方法最初由Osher和Sethian[8]提出，用于解決遵循熱力學方程的火苗外形變化過程。Mumford和Shah[9]于1989年提出了M-S模型。在此基礎上，Chan和Vese[10-11]分別于2001年和2002年提出了單水平集CV模型和多水平集CV模型。Li和Xu[12]于2010年提出了距離正則化水平集模型。

1 結合小波變換與中值濾波去噪

1.1 小波變換去噪

有效抑制噪聲是進行圖像分割的前提條件。小波變換[13]以其較好的時頻特性在圖像去噪領域受到了廣泛的關注。小波變換方法有2個關鍵步驟，分別是閾值函數的選取和閾值的選取。

閾值函數直接影響到信號去噪的效果。常用的閾值函數包括：硬閾值法、軟閾值法及模平方法等[14]。

設ω為小波系數，wλ為加入閾值后小波系數的大小，λ為閾值。

硬閾值函數如下所示：

(1)

軟閾值函數如下所示：

(2)

模平方函數如下所示：

(3)

硬閾值法去除高斯噪聲效果很好，但是圖像灰度變化較大的鄰域出現了振蕩，即振鈴現象。軟閾值法去除高斯噪聲效果較好，但是圖像邊緣出現了模糊。模平方法效果約等于硬閾值法與軟閾值法的折中。

針對上述閾值函數的不足，本文提出S型分段閾值函數，其特點是保持函數值的連續與光滑。在保持較好的去噪性能前提下，既避免了硬閾值函數的振鈴現象，又避免了軟閾值函數的邊緣模糊現象。

S型分段閾值函數如下：

(4)

函數圖像如圖1所示。

圖1 S型閾值函數圖像Fig.1 S-type threshold function

閾值的選取對于去噪效果同樣很關鍵。本文采用Donoho和Johnstorn[15]提出的統一閾值：

(5)

式中，δ為噪聲標準方差，N為信號的尺寸或長度。

1.2 中值濾波去噪

中值濾波去噪法是一種對圖像局部進行平均平滑的技術[16]，屬于一種非線性濾波。中值濾波的基本原理是將鄰域中的所有像素值按灰度級進行排列，取隊列的中間值為輸出像素值。它處理圖像的效果主要取決于2個要素：鄰域的空間范圍和模板所覆蓋的像素數。時間序列x(t)的中值濾波可以表示如下：

y(t)=MED{x(t-k),…,x(t),…,x(t+k)}，

式中，中值濾波器的長度為2k+1。改變濾波器的長度或者模板可改變濾波器的性能。

2 改進的雙水平集分割方法

2.1 單水平集模型描述

單水平集是依據水平集函數φ的取值將圖像劃分為目標和背景兩相。其適用于分割單個目標和灰度均勻目標的任務。圖2為單水平集模型原理圖。可見水平集函數φ將圖像分割成2個互不重合的部分。

圖2 單水平集區域劃分圖Fig.2 Single level set area division diagram

曲線擬合能量泛函定義為：

E(c1,c2,c)=μ·Length(c)+λ1?Ω1(u(x,y)-c1)2dxdy+

λ2?Ω2(u(x,y)-c2)2dxdy，

(6)

式中，μ,λ1,λ2為正常數，表示各能量項的權重系數，采用水平集φ后能量泛函為：

E(c1,c2,φ)=μ·Length(φ)+λ1?Ω1H(φ)(u(x,y)-c1)2dxdy+λ2?Ω2(1-H(φ))(u(x,y)-c2)2dxdy，

(7)

式中，H(*)為Hesviside函數。

曲線演化的目標是使得上述能量泛函取得極小值。根據多元函數求極值的理論，求泛函對c1的一階偏導數，并令它為0，可得c1，同理可得c2。公式如下：

(8)

根據求函數極值的梯度下降法，令水平集函數φ的導數等于能量泛函的負梯度得到水平集函數的演化方程如下：

(9)

2.2 雙水平集模型描述

為了將水平集方法推廣到分割多目標與灰度不均目標的任務，Chan和Vese提出了多相CV模型，即可通過n條水平集函數將圖像分割成2n個互不重合的部分。

圖3為四相CV模型原理圖。水平集函數φ1和φ2將圖像分割成Ω11，Ω12，Ω21，Ω22四個區域，它們之間互不重合。

圖3 四相CV區域劃分圖Fig.3 Division diagram of four phase CV region

定義能量函數如下：

(10)

式中，c11，c12，c21，c22分別對應Ω11，Ω12，Ω21，Ω22各區域的灰度均值，它們的計算方法如下：

(11)

考慮到演化過程中出現的重新初始化[15]的問題，在原來的能量函數中增加懲罰項p(φ)，表達式如下：

(12)

式中，ν為能量懲罰項的系數。

新的能量函數更新如下：

Ep=E(c,φ)+p(φ)。

(13)

針對傳統的水平集函數演化方程在演化后期容易出現演化速度放緩停滯等問題，提出了改進的演化方程。增加一組加速因子K1,K2作為乘法因子，自適應地調整加快演化速度。

(14)

圖4 加速因子Ki的函數圖像Fig.4 Function image of acceleration factor Ki

其中，

(15)

根據Euler-Lagrange方程和梯度下降法，最終的水平集函數φ1,φ2的演化方程為：

(16)

(17)

2.3 算法流程

① 首先加載原始圖像；

② 其次用小波變換去除高斯噪聲；

③ 然后用中值濾波去除椒鹽噪聲；

④ 最后用多水平集快速分割算法對圖像進行分割；

⑤ 輸出分割效果圖。

3 實驗結果與分析

本文實驗環境：軟件 Windows10, MATLAB 7.0;硬件Huawei MateBook X Pro, CPU CORE i7-8550U 1.80 GHz，RAM 8.0 G。

為了比較中值濾波和本文算法對高斯噪聲圖像的去噪效果，在腦部MRI圖像中加入高斯噪聲。分別用中值濾波和本文算法進行去噪，由圖5(b)和圖5(c)可以看出，中值濾波只能去除部分高斯噪聲，而本文算法能夠完全去除噪聲。

圖5 中值濾波與本文算法去噪效果對比圖Fig. 5 Comparison of denoising effect between median filter and the algorithm in this paper

為了比較小波變換和本文算法對椒鹽噪聲的去噪效果，分別用兩種算法對含有椒鹽噪聲的圖6(a)進行去噪，通過對比分割效果圖可以看出本文算法可以完全去除圖像中的椒鹽噪聲。需要說明：如果椒鹽噪聲密度很高,那么一次中值濾波可能無法消除全部噪聲。此時應該進行兩次或多次中值濾波。

圖6 小波變換與本文算法去噪效果對比圖Fig.6 Comparison of denoising effect between wavelet transform and the algorithm in this paper

圖7 傳統水平集分割效果圖Fig.7 Traditional level set segmentation renderings

圖8 本文模型分割效果圖Fig.8 Model segmentation renderings in this paper

通過對比圖7和圖8可知，對于傳統的雙水平集算法，當迭代到30次時，分割的效果只顯示出部分輪廓；當迭代到100次時，分割的效果基本理想，顯示出比較準確的輪廓。而對于快速水平集算法，當迭代到第30次時，分割的效果已經達到理想效果。繼續迭代到100次時，分割的結果也沒有明顯變化。實驗結果驗證了快速算法的有效性。

對于分割結果采用Jaccard Similarity(JS)指標進行評價[17]：

(18)

JS值越大表明表示分割效果越好，表1為本文算法與傳統水平集算法的JS指標，可看出本文算法的分割效果更好。

表1 2種算法的JS指標
Tab.1 JS index of two algorithms

圖像組織分割算法傳統水平集本文算法腦部MRI圖像1白質0.780.86灰質0.730.85

4 結束語

本文結合小波變換與中值濾波2種去噪算法的特點，對MRI圖像進行去噪預處理，兼顧了較好的去噪效果和保留圖像細節的目標。在傳統水平集方法中增加了自適應的加速因子以加快圖像分割速度，并且通過在能量函數中增加懲罰項解決曲線在演化過程中重新初始化的問題。實驗表明，本文的方法可以快速有效地分割含有高噪聲、灰度不均的MRI圖像。