融合雙邊濾波算子的多相水平集圖像分割算法

史 娜， 孔慧華， 秦 鵬

(中北大學理學院數學系， 太原 030051)

超聲檢查是一種無輻射的安全檢查，對人體沒有創傷、沒有損害，現已被廣泛用于醫學臨床診斷。乳腺超聲圖像的準確快速分割對于婦科疾病的臨床診斷具有重要意義，有助于醫生提高診斷的正確率。乳腺超聲圖像具有低對比度的特點，其斑點噪聲較多，腫瘤區域的邊緣模糊且形狀不規則。對于惡性乳腺腫瘤而言，待分割的目標邊緣形狀復雜、內部灰度較不均勻、目標與背景的對比度較低，因此，準確得到乳腺腫瘤病灶區域的邊緣輪廓，是一項有難度且富有挑戰的工作。

目前，針對醫學超聲圖像分割問題的活動輪廓模型主要有以下3種類型:參數活動輪廓模型、幾何活動輪廓模型、無邊緣活動輪廓模型等。

根據曲線演化的不同，1988年，Kass等[1]提出了演化曲線參數化表示的活動輪廓模型，也稱為蛇(snake)模型，該模型利用初始動態閉合曲線在感興趣區域特征的約束下進行演化，當能量泛函達到極小值時，目標區域的邊緣輪廓就是參數化表示的曲線位置。雖然蛇模型可以提供光滑閉合的邊緣曲線，但無法自動處理運動界面的拓撲結構變化和局部形變。

1988年，Osher等[2]提出經典的幾何活動輪廓模型，該算法將N維曲面或曲線的演化問題，轉化為N+1維水平集函數的隱函數來求解，避免曲線的重新參數化表示，使曲線的演化更加靈活，具有自適應演化曲線拓撲結構變化的優勢。2001年，Chan等[3]以水平集方法為基礎，利用簡化的Mumford-Shah泛函，建立分段常值的水平集分割模型(CV模型)；1997年，Caselles等[4]提出測地線活動輪廓模型；2002年，Vese等[5]在考慮全局區域信息的CV模型基礎上，提出分片光滑和分段常值的多相水平集分割模型。

由于CV模型不能很好地解決灰度非均勻圖像的分割問題，2007年，Li等[6]提出一種基于圖像局部灰度信息擬合的活動輪廓(local binary fitting, LBF)模型，該模型以高斯函數為核函數對圖像進行平滑處理；2009年，Wang等[7]結合全局能量泛函，將LBF模型改進為全局和局部能量信息相結合的圖像分割模型。

近期，中外學者在超聲圖像分割領域，還提出基于多相水平集的活動輪廓模型[8-13]，王潔等[14]提出基于Vese-Chan模型的連續最大流方法，劉磊等[15]提出一種基于加權指數平均比率算子改進的CV模型，蘭紅等[16]利用改進的CV模型對磁共振圖像進行分割。

針對乳腺超聲圖像中腫瘤區域的結構特征，提出一種融合雙邊濾波算子的多相水平集圖像分割模型。改變以往將高斯函數作為核函數，引入雙邊濾波算子來進行邊緣平滑，作為非線性濾波算子，它能有效保護邊緣信息。結合超聲圖像的全局和局部灰度能量信息，以及多相水平集函數，使分割算法更適用于確定乳腺腫瘤區域的邊界。

1 傳統的無邊緣活動輪廓模型

1.1 基于全局信息的Chan-Vase模型

傳統的無邊緣活動輪廓CV模型[3]，其算法以簡化Mumford-Shah泛函為基礎，利用初始演化曲線內外部圖像灰度的分段常值，去逼近目標和背景區域。將曲線長度懲罰項和區域面積的懲罰項加入能量泛函，并用φ表示水平集函數，C表示演化曲線，對于待分割圖像Ω(x，y)，所定義的能量泛函為

(1)

式(1)中:μ、ν、λ1、λ2分別為取值為正的參數，演化曲線C內部區域為Ω1，其平均灰度為c1；外部區域為Ω2，其平均灰度為c2，表達式為

(2)

(3)

式中：H(φ)通常為Heaviside函數；ε為很小的正數，即

(4)

δ(φ)=H′(φ)

(5)

通過變分法最小化能量泛函，得到水平集演化的偏微分方程，采用有限差分法進行數值近似計算，梯度下降流為

(6)

式(6)中：面積參數ν≥0；尺度參數μ>0。

1.2 基于局部信息的LBF模型

2007年，針對CV模型無法正確分割灰度非均勻圖像的問題， Li等[6]提出一種局部二值能量擬合的活動輪廓(local binary fitting model，LBF)模型，該模型提供了有效的解決方案。

對于定義域為Ω的灰度圖像I(x，y)，給定點x，y∈Ω，待分割圖像的局部灰度信息的能量泛函為

|I(y)-f1(x)|2H[φ(y)]dy+λ2×

{1-H[φ(y)]}dy

(7)

式(7)中:H(·)為Heaviside函數;λ1和λ2分別為取值為正的參數;f1和f2這兩個平滑函數分別為輪廓曲線C內部和外部的圖像局部結構信息，而Gσ為標準差為σ高斯核函數，表達式為

(8)

f2(x)]dx+νL(C)

(9)

根據變分水平集方法，求解能量泛函式(9)的極小值，可得f1和f2的表達式為

(10)

最終，LBF模型的曲線演化方程為

(11)

式(11)中:κ為演化曲線的曲率;e1和e2滿足:

(12)

2 融合雙邊濾波算子的多相水平集分割模型

2.1 雙邊濾波算子

雙邊濾波算子(bilateral filtering)是一種依賴像素點的空間鄰近相似度和灰度差異度的非線性濾波器[17]。雙邊濾波在高斯濾波的基礎上，引入像素的灰度值進行局部加權平均，針對圖像的斑點噪聲進行平滑時，能夠較好地保持邊緣結構特征。

雙邊濾波算子的濾波核權重ωi在空間近鄰函數Gs基礎上，增加灰度相似函數Gr，表達式為

(13)

(14)

(15)

(16)

(17)

式中:i為像素索引；j為窗口ωi覆蓋區域的像素索引；x為像素空間位置坐標；I為像素灰度值；參數σs為空間近鄰因子；參數σr為灰度相似因子；σs越大平滑效果越好，σr越小保邊效果越好。

2.2 多相水平集方法

2002年，針對分段常值和分片光滑圖像，Chan等[5]采用N個水平集函數劃分2N個目標區域，建立了多相水平集的圖像分割模型。

首先，引入兩個水平集函數φ1和φ2，待分割區域被劃分為四個互不相交、互不重疊的部分(圖1):

圖1 2個水平集函數劃分4個圖像分割區域

Heaviside函數和Dirac函數的表達式如下:

針對2個水平集函數劃分4個圖像區域的特征函數分別為M1、M2、M3、M4，表達式如下:

2.3 基于全局或局部的活動輪廓模型

針對乳腺超聲圖像邊界模糊、感興趣區域對比度較低、斑點噪聲較大的情況，在CV模型和局部二值擬合模型的基礎上，綜合考慮乳腺超聲圖像的全局和局部區域的灰度能量信息，建立全局和局部二值擬合模型(local and global binary fitting model， LGBF)，其能量泛函的水平集函數表達式為

ELGBF(φ，f1，f2，c1，c2)=(1-α)EL(φ，f1，f2)+

αEG(φ，c1，c2)+ER

(18)

H[φ(y)]}dy]dx

(19)

H[φ(x)]}dx

(20)

式中:輸入圖像為I:Ω→R；點x，y∈Ω；C為閉合演化曲線；λ1、λ2≥0；圖像區域Ω被曲線C分為目標(C內部)和背景(C外部)兩個區域;c1為曲線C內部平均灰度值;c2為曲線C外部的平均灰度值;φ為水平集函數;H為Heaviside函數。局部能量泛函為EL(φ，f1，f2);全局能量泛函為EG(φ，c1，c2);α為圖像的全局信息參數;GBF為雙邊濾波算子。

為保持乳腺超聲圖像演化曲線的光滑性，同時，避免重新初始化水平集函數，在總的能量泛函ELGBF中引入懲罰項ER，包括:長度懲罰項L(φ)和能量懲罰項P(φ)。當總能量泛函ELGBF達到最小值時，長度約束項L(φ)使得演化曲線C保持盡可能短且光滑，能量約束項P(φ)使得水平集函數φ在曲線演化過程中保持為符號距離函數。懲罰項ER的表達式為

(21)

式(21)中:長度懲罰項參數ν>0;能量懲罰項參數μ>0。

為更準確地獲取乳腺腫瘤病灶區域的細節信息，可將ELGBF模型推廣到多相水平集情況，滿足乳腺腫瘤超聲圖像的多目標分割需求。具體的能量泛函定義為

ELGBF(φ1，φ2，f1，f2，f3，f4，c1，c2，c3，c4)=

ν[Lε(φ1)+Lε(φ2)]+μ[P(φ1)+P(φ2)]

(22)

(23)

式(23)中：f1、f2、f3和f4分別為演化曲線C1和C2所分割的四個圖像區域的局部灰度二值擬合能量函數;c1、c2、c3和c4分別為演化曲線C1和C2所分割的4個圖像區域中的全局平均灰度值。

由變分法，固定水平集函數φ1和φ2，最小化總的能量泛函ELGBF，得到滿足Euler-Lagrange方程的局部灰度二值擬合的能量函數f1、f2、f3和f4，以及全局灰度平均值c1、c2、c3和c4，具體表達式為

(24)

(25)

對應最小化φ1的梯度下降流方程為

(26)

式(26)中:div(·)為散度算子；δε(·)為正則化的Dirac函數，擬合能量為F13=Hε(φ2)(e3-e1);F24=[1-Hε(φ2)](e4-e2)。

同理，可得φ2的梯度下降流方程為

(27)

在乳腺腫瘤超聲圖像的分割過程中，演化曲線在全局和局部區域能量泛函的相互作用下，內外力相互補充，當演化曲線接近病灶區域的邊界時，反映乳腺超聲圖像局部結構信息的能量函數占主導地位；當演化曲線遠離病灶區域時，表征乳腺超聲圖像全局結構信息的能量函數起主要作用。

3 實驗結果與分析

3.1 實驗數據集

使用的真實病例的乳腺超聲圖像進行實驗，包括67幅惡性乳腺腫瘤的超聲圖像，圖像的長寬均在300～400像素。該數據集中的所有乳腺腫瘤超聲圖像，均配有醫生手工標記的病灶區域。在實際圖像分割實驗時，對每幅圖像邊緣的機器信息和圖像參數信息進行剪裁。

3.2 實驗結果的評價函數

乳腺腫瘤超聲圖像具有以下結構特征:①病灶區域的灰度值比周圍正常組織小很多，并且灰度值相對均勻，但是病灶邊緣不規則且拓撲結構復雜；②由于超聲成像機制所限，可檢測到的乳腺腫瘤區域的直徑通常大于5 mm，若小于該閾值，可視為背景區域；③在對病人進行超聲檢查時，可將病灶區域限制在超聲圖像中心區域的窗口內。

為了定量化評價各模型的分割效果，采用的評價指標分別是:Dice相似系數、Jaccard相似度、假陽性率RFP及假陰性率RFN。具體計算公式為

(28)

式(28)中:N(·)為區域中像素的個數；Sg為醫生手工標注的腫瘤區域；Sm為分割算法所獲得的分割區域；O為Sg和Sm的公共區域。Dice相似系數和Jaccard相似度都能夠表示算法分割結果Sm與醫生手動標記Sg結果的重疊率，Dice相似系數和Jaccard相似度的值越接近于1，分割越準確。假陽性率RFP的值表示醫生手工標記的病灶區域去除公共區域的占比，而假陰性率RFN的值表示分割算法所得分割區域去除公共區域的占比，RFP和RFN的值越接近于0，分割效果越好。

3.3 結果與分析

為驗證本文算法的可行性，首先對仿真圖像進行分割實驗。將本文中提出的LGBF模型與CV模型和LBF模型的分割結果進行對比。然后，對乳腺超聲圖像數據集進行分割實驗，利用評價指標對各算法的分割效果做出定量評價。最后，將融合雙邊濾波算子的LGBF模型推廣到多相水平集情況，能夠得到乳腺腫瘤內部的病灶細節信息。

實驗環境為Intel Core i7-9700CPU@ 3.00 GHz，內存為8 GB，64 bit操作系統，實驗軟件為MATLABR2014b。實驗中常用的參數設置為ν=0.001×2552，μ=1，λ1=1，λ2=1，ε=1，Δt=0.1。

在仿真圖像的分割實驗中，參數設置如下:全局結構信息的權重參數α=0.4，雙邊濾波算子的窗寬ω=4，空間近鄰因子σs=4，灰度相似因子σr=25。將傳統的CV模型、LBF模型與本文LGBF模型進行對比分割實驗[分辨率如圖2(a)左下方所示]。活動輪廓模型的初始演化曲線如圖2(b)所示。由圖2(c)可見，CV模型只考慮全局灰度信息，故無法正確分割灰度不均勻的圖像。針對第1、2幅圖像，LBF模型能夠有效分割灰度不均勻圖像，然而，在第3幅圖像的分割實驗中，由于只注重局部灰度信息，在成功分割染色體的同時，出現對染色體過度分割的現象(如圖2中第3行第4列所示)。

本研究提出的LGBF模型充分考慮了圖像的局部和全局灰度能量信息，采用雙邊濾波算子對圖像進行預處理，能夠有效分割出目標區域，分割效果良好，如圖2(e)所示。

圖2 灰度仿真圖像的分割結果

然后，針對乳腺惡性腫瘤的超聲圖像數據集進行分割實驗。在數據集中，病灶區域的邊緣輪廓可以分為兩類，一類是腫瘤區域的邊界模糊,如圖3中第1、2幅圖像所示;另一類是腫瘤區域的邊界較光滑，如圖3中第3、4幅圖像所示。

本文的融合雙邊濾波算子的LGBF模型，在進行乳腺超聲圖像進行分割實驗時，根據乳腺超聲圖像的結構特征，實驗參數設置如下:全局結構信息的權重參數α=0.3，雙邊濾波算子的窗寬ω=9，空間近鄰因子σs=4，灰度相似因子σr=15。

圖像分辨率和初始演化曲線如圖3(a)和圖3(b)所示。從圖3可以看出，CV模型在對乳腺病灶區域進行分割時，演化曲線受全局灰度能量的作用，將圖像灰度信息較一致的區域進行分割，導致病灶區域和正常組織區域混淆在一起，不能得到正確的病灶區域，如圖3(c)所示。而傳統的LBF模型過分注重超聲圖像的局部能量信息，在分割出病灶區域的同時，目標區域和背景區域出現了過度分割的現象，如圖3(d)所示。

本文融合雙邊濾波算子的LGBF模型，在曲線演化過程中，利用雙邊濾波算子對超聲圖像進行平滑處理后，演化曲線受乳腺腫瘤超聲圖像全局和局部灰度能量結構信息的綜合作用，能夠準確逼近醫生手工標記的病灶區域邊界,如圖3(e)所示。最后，將本文所提算法推廣到多相水平集情況，得到乳腺超聲圖像病灶區域的細節分割效果,如圖3(f)所示。利用評價函數，對不同算法的實驗結果進行定量分析，分別對比Jaccard相似度、Dice系數、假陽性率RFP以及假陰性率RFN，發現本文算法的指標分別為0.938 7、0.945 1、0.074 2、0.083 3，均優于傳統的CV模型和LBF模型，如表1所示。

表1 三種算法的分割結果對比

圖3 乳腺腫瘤超聲圖像的分割結果

4 結論

融合雙邊濾波算子對乳腺腫瘤超聲圖像進行平滑處理，同時，根據乳腺超聲圖像中腫瘤區域與其他人體組織結構的灰度信息有相似性特點，所建立的能量泛函要綜合考慮超聲圖像的全局和局部灰度能量信息。因此，本文多相LGBF模型優于傳統的CV模型和LBF模型，能夠準確得到病灶區域的邊緣輪廓，分割曲線的平滑性優于其他算法，分割結果與醫生手動標記的病灶區域相吻合。