朱家明，許 諾，張倚天，包 清，杭 俊

(揚(yáng)州大學(xué) 信息工程學(xué)院，江蘇 揚(yáng)州 225127)

0 引言

隨著科技的發(fā)展，很多醫(yī)生把醫(yī)學(xué)影像作為疾病診斷的重要依據(jù)之一。由于核磁共振成像(MRI)圖像具有對(duì)含水分的軟組織較為敏感的特性，經(jīng)常用于對(duì)腦部組織進(jìn)行解剖結(jié)構(gòu)的分析與檢測(cè)[1]。但是MRI圖像的成像也有缺點(diǎn),例如由于受到設(shè)備、環(huán)境和不均勻的射頻信號(hào)影響，圖像大多存在噪聲和灰度不均的問(wèn)題[2]。

在醫(yī)學(xué)圖像處理領(lǐng)域，圖像分割是一個(gè)重要的技術(shù)[3]，其任務(wù)是從圖像中提取有診斷價(jià)值的目標(biāo)，是用于醫(yī)生臨床診斷的計(jì)算機(jī)輔助系統(tǒng)的重要組成部分。對(duì)于圖像分割，各國(guó)學(xué)者提出了多種方法，例如：閾值法、模式識(shí)別法、區(qū)域生長(zhǎng)法、神經(jīng)網(wǎng)絡(luò)法、小波變換法、模糊分割法、水平集方法[4]以及遺傳算法等。其中水平集方法是曲線演化問(wèn)題的一種有效的求解方法[5-7]，其最初由Osher和Sethian[8]提出，用于分析火苗的外形變化過(guò)程。在1989年，Mumford和Shah[9]提出了M-S模型。Chan[10]于2001年提出了單水平集CV模型，Vese[11]于2002年提出了多水平集CV模型，Li和Xu[12]于2010年提出了距離正則化水平集模型。為了解決MRI圖像的分割問(wèn)題，本文提出了一種結(jié)合閾值法預(yù)分割的水平集分割方法?？紤]到MRI圖像具有的高斯噪聲，先對(duì)其作低頻濾波處理。理論分析和實(shí)驗(yàn)結(jié)果都驗(yàn)證了該方法的有效性。

1 頻域?yàn)V波去噪

1.1 高斯噪聲與頻域?yàn)V波

作為圖像分割的前提條件，要求對(duì)噪聲先作有效抑制。高斯噪聲是醫(yī)學(xué)圖像中廣泛存在的一種噪聲，通常是由于電子元件的敏感性以及各種外界的電磁干擾造成的。消除高斯噪聲常用的方法有空間濾波與頻域?yàn)V波。因?yàn)樵肼曅盘?hào)一般都是高頻信號(hào)，當(dāng)采用低通濾波時(shí)，選擇合適的截止頻率就能夠?yàn)V除大部分噪聲，同時(shí)保留圖像的邊緣和細(xì)節(jié)。而采用低通濾波的前提就是先對(duì)圖像進(jìn)行傅里葉變換，使之轉(zhuǎn)換成頻域模型；然后針對(duì)頻域模型設(shè)計(jì)低通濾波器進(jìn)行濾波；最后對(duì)濾波后的頻域模型進(jìn)行傅里葉反變換，得到去噪的圖像。

計(jì)算機(jī)處理的都是離散化的系統(tǒng)，一幅灰度圖像可以表示為二維的離散函數(shù)，離散函數(shù)的二維傅里葉變換與傅里葉反變換公式分別為：

(1)

(2)

1.2 低通濾波

低通濾波方法有巴特沃斯濾波、切比雪夫?yàn)V波以及理想低通濾波等。其中理想低通濾波器屬于假想的低通濾波器，對(duì)于高頻(即高于截止頻率)信號(hào)完全截止，而對(duì)于低頻信號(hào)完全無(wú)失真?zhèn)鬏?。由于理想低通濾波去噪的效果要優(yōu)于其他低通濾波，因此本文采用理想低通濾波去除噪聲，如圖1所示。

(a) FFT頻譜

(b) FFT頻譜

(c) ButterworthLPF頻譜

(d) 理想LPF頻譜

(e) 原始圖像

(f) 噪聲圖像

(g) ButterworthLPF圖像

(h) 理想LPF圖像

2 改進(jìn)的水平集分割方法

2.1 閾值法的預(yù)分割

用閾值法先預(yù)分割，即設(shè)定水平集函數(shù)初值。水平集分割法步驟：先定一個(gè)初值(初始輪廓)，再按照演化方程進(jìn)行演化,即不斷更新水平集函數(shù)取值。相應(yīng)改變了零水平集的輪廓曲線，完成了圖像分割。當(dāng)達(dá)到設(shè)定標(biāo)準(zhǔn)，就停止演化并輸出結(jié)果。因?yàn)樗郊某踔禃?huì)對(duì)分割結(jié)果有一定的影響，所以選取合適的初值是有必要的。本文采用閾值法預(yù)分割(定水平集初值)，用整幅圖像灰度區(qū)間的黃金分割數(shù)作為閾值，即取權(quán)值α=0.382。閾值的計(jì)算公式如下：

Ith=Ilow+0.382(Ihigh-Ilow)，

(3)

式中，Ith為閾值，Ilow為整幅圖像的灰度值的下限，Ihigh為灰度值的上限。其中默認(rèn)目標(biāo)灰度值比背景灰度值小，反之如果目標(biāo)灰度值比背景灰度值大，則取α=0.618，計(jì)算公式如下：

Ith=Ilow+0.618(Ihigh-Ilow)。

(4)

閾值法確定的初始輪廓圖如圖2所示。

(a) α=0.382

(b) α=0.618

(c) α=0.85

2.2 水平集模型描述

水平集方法是依據(jù)水平集函數(shù)φ的取值將圖像劃分為目標(biāo)和背景兩相。水平集函數(shù)定義為符號(hào)距離函數(shù)。即圖上一點(diǎn)(x,y)與目標(biāo)輪廓曲線C的距離。表達(dá)式如下：

(5)

由式(5)可見(jiàn)，水平集函數(shù)取值為零的點(diǎn)的集合(零水平集)即為目標(biāo)的邊緣輪廓。水平集函數(shù)原理如圖3所示，從圖中可見(jiàn)水平集函數(shù)φ將整個(gè)圖像分割成互不重合的兩個(gè)部分。

圖3 水平集區(qū)域劃分圖Fig.3 Level set area division diagram

水平集模型主要包括兩大類模型：一類是基于邊緣信息的模型，如Snake模型；另一類是基于區(qū)域信息的模型，如M-S模型、C-V模型和變分水平集等。本文考慮結(jié)合邊緣信息和區(qū)域信息的模型。

2.3 邊界擴(kuò)展函數(shù)

本文設(shè)計(jì)一種新型的邊界擴(kuò)展函數(shù)，并將它結(jié)合到傳統(tǒng)的能量函數(shù)中，用于解決目標(biāo)邊界劃分細(xì)節(jié)不準(zhǔn)確的問(wèn)題?；驹O(shè)想如下：在目標(biāo)輪廓線附近，當(dāng)目標(biāo)外部灰度接近目標(biāo)內(nèi)部的灰度時(shí)，就使得目標(biāo)輪廓線擴(kuò)張，把外面這部分包含到目標(biāo)內(nèi)部來(lái)。該邊界擴(kuò)展函數(shù)設(shè)計(jì)如下：

(1-H(φ))dxdy，

(6)

式中，u(x,y)表示目標(biāo)外部某一點(diǎn)的灰度，cδ表示該點(diǎn)附近的目標(biāo)局部的平均灰度。

。

(7)

水平集函數(shù)φ的演化方程為：

(8)

2.4 能量函數(shù)與演化方程

主動(dòng)輪廓模型的能量泛函定義為：

E(C)=∮Ωg(|?I(C(s))|)ds。

(9)

對(duì)式(9)引入Hesviside函數(shù)，結(jié)合格林公式，可改寫成：

Es(φ)=?Ωg(x,y)‖?H(φ)‖dxdy，

(10)

式中，φ表示水平集函數(shù)，g(x,y)為反映邊緣特征的函數(shù)，表達(dá)式如下：

p≥1。

(11)

曲線演化的目標(biāo)是使能量泛函取得極小值。根據(jù)E-L方程和梯度下降法，解得演化方程為：

(12)

C-V模型的能量泛函定義為：

E(c)=μ·length(C)+

λ1?Ω1(u(x,y)-c1)2dxdy+

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

(13)

式中，μ、λ1、λ2為正的常數(shù)，u(x,y)為圖像灰度函數(shù)。引入Hesviside函數(shù)后能量泛函為：

Ecv(φ)=μ?Ω‖?H(φ)‖dxdy+

λ1?ΩH(φ)(u(x,y)-c1)2dxdy+

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

(14)

對(duì)應(yīng)的演化方程為：

(15)

將泛函(14)對(duì)c1求偏導(dǎo)，令它等于0可求得c1。同理可求得c2,公式如下：

(16)

(17)

為了避免在演化過(guò)程中出現(xiàn)重新初始化的問(wèn)題[15]，在能量泛函(14)中增加懲罰項(xiàng)Ep(φ):

(18)

式中，ν為正的常數(shù)。綜上所述，能量函數(shù)定義為：

E(φ)=Ecv(φ)+Ep(φ)+Es(φ)+Eex(φ)=

μ?Ω‖?H(φ)‖dxdy+

λ1?ΩH(φ)(u(x,y)-c1)2dxdy+

λ2?Ω(1-H(φ))(u(x,y)-c2)2dxdy+

α?Ωg(x,y)‖?H(φ)‖dxdy+

。(19)

對(duì)應(yīng)的演化方程為：

λ1(u-c1)2+λ2(u-c2)2+

(20)

式中，Δ表示拉普拉斯算子。

2.5 算法流程

算法流程為：① 加載原始圖像；② 用理想低頻濾波消除噪聲；③ 用閾值法進(jìn)行預(yù)分割；④ 用水平集方法進(jìn)行圖像分割；⑤ 輸出分割結(jié)果圖。

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

實(shí)驗(yàn)環(huán)境：硬件Huawei MateBook X Pro,RAM 8.0 G，CPU CORE i7-8550U 1.80 GHz; 軟件 Windows10中文版,MATLAB 7.0。選取一幅腦部MRI圖像進(jìn)行去噪預(yù)處理，迭代次數(shù)為200次。參數(shù)設(shè)置為：α=3,μ=5,v=1.5,λ1=λ1=10-6,采樣時(shí)間間隔Δt=0.02。分別采用C-V算法和本文算法來(lái)分割去噪后的MRI圖像，分割結(jié)果如圖4所示。由圖可見(jiàn)，本文的算法分割的效果優(yōu)于C-V算法。

(a) 原始圖像

(b) C-V算法

(c) LSE算法

采用Jaccard Similarity(JS)指標(biāo)來(lái)對(duì)分割結(jié)果進(jìn)行評(píng)價(jià)[13]：

(21)

指標(biāo)值越大代表分割效果越好，表1為C-V算法與本文算法的JS指標(biāo)，數(shù)據(jù)顯示本文算法的分割效果優(yōu)于C-V算法。

表1 兩種算法的JS指標(biāo)

4 結(jié)束語(yǔ)

本文提出一種結(jié)合新型邊界擴(kuò)展函數(shù)的水平集圖像分割方法。首先用理想低頻濾波對(duì)圖像進(jìn)行去噪處理，同時(shí)兼顧去噪效果和保留細(xì)節(jié)兩個(gè)目標(biāo)，然后用閾值法確定目標(biāo)的初始輪廓。提出新型邊界擴(kuò)展函數(shù)，并整合到能量函數(shù)中。在其中增加懲罰項(xiàng)，避免了在演化過(guò)程中的重新初始化。實(shí)驗(yàn)結(jié)果表明，該方法能夠準(zhǔn)確地分割含有高斯噪聲的灰度不均的MRI圖像。