自適應尺度參數的非均勻圖像快速水平集分割

龔進慧 張貴倉

(西北師范大學數學與統計學院 甘肅 蘭州 730070)

0 引 言

圖像分割在計算機視覺和圖像處理領域仍然占有重要的地位。圖像分割面臨的主要挑戰是如何分割具有噪聲、強度不均勻和紋理復雜的圖像。最近幾十年來，學者們提出了許多圖像分割方法。主動輪廓模型( Active Contour Model，ACM)[1]自提出以來一直受到人們的廣泛關注。ACM算法能夠提供平滑、封閉的輪廓線，以亞像素的精度覆蓋目標邊界。水平集方法是由Osher等[2]提出的活動輪廓模型族的成員之一。水平集方法用較高維函數的零水平集來表示輪廓，通常稱為水平集函數，并且將輪廓的運動轉換為水平集函數的演化。根據能量泛函中嵌入的圖像性質，現有的水平集方法可分為兩類：基于邊緣的模型[3-5]和基于區域的模型[6-10]。

基于邊緣的模型通常使用局部圖像的梯度信息來迫使活動輪廓朝目標的期望邊界移動，這對于分割具有清晰邊界的圖像特別有效。然而，這些模型對噪聲和初始輪廓都很敏感，在梯度值比較小的弱邊界上很容易發生邊緣泄漏。基于區域的模型利用圖像統計信息來推動活動輪廓覆蓋目標邊界。它們對噪聲和輪廓初始化不太敏感，因為它們利用區域圖像信息來驅動水平集演化。此外，它們還可以分割具有弱邊界甚至無邊界的圖像。最流行的基于區域的模型之一是Chan-Vese(CV)模型[6]，這是眾所周知的Mumford-Shah(MS)分割模型[11]的分段常量情況。文獻[9，12-13]提出了一種基于區域的壓力函數，并引入了一種快速的水平集演化公式，以使模型更快地收斂。隨后，將局部灰度信息廣泛應用于主動輪廓模型中，提高了噪聲和強度不均勻性對圖像分割的精度[14-18]。文獻[14-15]將局部二元擬合(LBF)能量引入基于區域的活動輪廓模型的能量函數，以獲得具有強度不均勻性的圖像的滿意分割結果。然而，基于局部區域的模型對初始輪廓也很敏感，容易陷入局部極小值。

強度不均勻性是現實世界圖像中的常見現象，尤其是在磁共振(MR)圖像中，其通常由成像裝置的缺陷或照明變化引起。具有強度不均勻性的圖像通常被描述為一個真實的分段常數圖像乘以一個空間變化的光滑場[18]。這種空間變化的光滑場稱為偏置場。偏置校正方法已被廣泛研究，以處理具有強度不均勻的圖像[19]。許多非均勻圖像分割方法[20-23]都是基于偏置校正方法的。Li等[20]提出了一種新穎的水平集模型(LIC)，該模型考慮局部圖像強度信息以演化水平集函數，并利用局部加權K均值聚類方法來評估偏置場。Zhang等[22]提出了一種局部統計活動輪廓模型(LSACM)，它引入了聚類方差信息。然而，這些模型對初始輪廓的位置也很敏感。此外，Li等[20]指出核函數的尺度參數應根據強度不均勻性的程度適當選擇。然而，在LIC模型中卻沒有提到如何估計強度不均勻性的程度。Min等[24]提出了一種新的雙極小化模型(DM)，該模型采用多層結構自適應地確定每個像素的尺度參數。

本文提出一種新的水平集非均勻圖像分割方法。從非均勻圖像模型出發，推導出最優分割平面，為每個像素提供最優分割。在平面上，提出一種新的基于區域的壓力函數，并利用該函數在整個圖像區域的水平集公式中定義了一個能量泛函。該方法通過對能量泛函進行最小化，在對非均勻圖像進行分割的同時，對偏置場進行估計。此外，還引入了一個新的偏置場初始化，以提高初始輪廓的魯棒性，這可以確保分割過程首先聚焦在有很大變化的區域，然后延伸到整個圖像區域，而不管初始輪廓的位置如何。為了準確估計偏置場，針對核函數設計了一種新的自適應尺度參數。此外，在數值實現中采用了一種簡單快速的水平集進化公式來提高效率。

1 背景模型

1.1 CV模型

Chan等[6]提出的CV模型假定物體和背景的強度屬于兩個可分離的區域。該模型的主要概念是找出圖像的具體劃分，將圖像劃分為兩個聚類。對于域Ω中給定的圖像I，通過最小化以下能量泛函來建立CV模型：

E=γ·length(C)+v·Area(inside(C))+

(1)

式中：γ≥0、v≥0、λ1>0、λ2>2是固定參數;length(C)是閉合輪廓線C的長度；Area是C的內部區域面積；c1和c2是輪廓C內部和外部的平均強度。在水平集方法中，等高線C、inside(C)和outside(C)的不相交區域表示為：

通過將式(1)相對于水平集函數φ最小化，相應的梯度流方程為：

v-λ1(I-c1)2+λ2(I-c2)2]

(2)

式中：δ是狄利克雷函數；▽是梯度算子，div()是散度算子。

同樣，c1和c2是通過以下方式獲得的：

(3)

(4)

式中：H(φ(x))是Heaviside函數。

CV模型對灰度均勻圖像分割效果良好。此外，它對初始輪廓的位置不敏感。但是，當輪廓C內外的圖像強度不均勻時，常數c1和c2可能不能準確地描述圖像。因此，CV模型可能無法對灰度不均勻的圖像進行分割。

1.2 張氏模型

張氏模型[9]將基于區域的模型與基于梯度的圖像分割模型相結合，在該模型中引入基于區域的壓力函數S。基于區域的模型的能量泛函被定義為：

(5)

式中：α是固定參數；μ是對象的灰度級下限。

下限μ可以作為圖像的最優分割平面。因此，關于水平集函數φ最小化式(5)，并促使輪廓附著到灰度級大于μ的區域。相應的梯度流動方程為：

(6)

在CV模型中，當忽略前兩項并且λ1=λ2=1時，能量泛函的梯度流動方程簡化為：

(7)

1.3 LIC模型

為了正確分割強度不均勻的圖像，Li等[20]提出了LIC模型。該方法基于一個常用的圖像模型，其定義為：

I=bJ+n

(8)

式中：I是所測量的圖像強度；b是偏置場；J是真實圖像；n是加性噪聲。在文獻[20]中，關于真實圖像J和偏置場b的兩個假設如下：

(1) 偏置場b在圖像域Ω中變化緩慢。

(2) 真實圖像J可以近似地以在每個不相交區域中的常數來表示，即J(x)≈ci,x∈Ωi，其中{Ωi}i=1，2，…，N是圖像域Ω的分區。

能量方程被定義為：

(9)

式中：K(·)是核函數。

當考慮LIC模型中的兩相情況并引入水平集函數φ時，能量函數可以寫成：

(10)

式中：M1(φ(x))=H(φ(x))，M2(φ(x))=1-H(φ(x))是每個集群的成員函數。

最佳的ci和b可通過以下方法獲得：

(11)

(12)

式中：*代表卷積運算。

對于固定的c和b，通過求解以下梯度流方程，可以得到E相對于φ的極小值。

(13)

LIC模型采用加權K-均值聚類方法對偏置場進行評價。該方法通過偏置校正，可以分割出一些強度不均勻的圖像。然而，LIC模型對初始輪廓位置的魯棒性不強，核函數的尺度參數是經驗性選擇的，這是LIC模型的固有缺陷。另外，在LIC模型中采用有限差分格式進行數值實現。有限差分格式的時間步長T應滿足CFL條件。為了保持數值穩定性，CFL條件將時間步長Δt限制為很小的長度，從而導致收斂緩慢。

2 本文方法

2.1 基于區域的壓力函數設計

在LIC模型中，從式(13)導出以下方程：

(14)

卷積運算可以看作是一個獲得加權平均值的過程。在假設(1)的基礎上，可以得到：

|K*b2-(K*b)2|<ξ

式中：ξ是一個小的正常數。因此，可以假設K*b2≈(K*b)2，上面的方程可以表示為：

(15)

式中：α=2(K*b)(c1-c2)。

(16)

式中：K是具有標準差σb的高斯核函數。

在考慮N=2并引入水平集函數φ的情況下，可以將能量泛函寫成：

(17)

類似于LIC模型，可以得到最優ci和b：

(18)

(19)

Heaviside函數H(φ)被正則化如下：

(20)

相應的梯度流動方程是：

(21)

2.2 新的偏置場初始化

當初始化偏置場b(x)=1,x∈Ω時，LIC模型在第一次迭代中退化為CV模型。在第一次迭代之后，零水平集可能遠離目標的實際邊界，這可能導致錯誤的偏置場估計和分割結果。因此，LIC模型對初始輪廓的位置敏感。為了克服這個缺點，引入了一個新的偏置場，它可以確保分割過程首先聚焦在物體真實邊界附近的區域，然后擴散到整個圖像區域。

由于整體圖像強度變化的趨勢與偏置場的趨勢相似，因此通過保持這一趨勢可以更好地估計偏置場。根據圖像模型I=bJ+n和假設a，即在圖像域Ω中，偏置場b變化緩慢，新的偏置場初始化定義為：

b0=K*(I/N0)

(22)

式中：N0是圖像強度的平均值。

通過偏差校正，校正后的圖像為：

(23)

S=I-K*(K*I)

(24)

由于目標邊界位于圖像強度變化較大的區域，新的偏置場初始化可以保證分割過程集中在真實邊界附近的區域。不考慮初始輪廓的位置，因此可以先分割目標實際邊界附近的圖像。隨著水平集的演化，遠離目標真實邊界的區域可以被分割。新的偏置場初始化可以提高初始輪廓定位的魯棒性。

2.3 自適應尺度參數

核函數的尺度參數在偏置場估計中起著重要作用。小尺度參數可能使算法陷入局部極小值，而較大的參數可能無法正確估計偏置場。為了克服這一缺點，設計了一種新的自適應尺度參數，該參數依賴于局部圖像的不均勻程度。根據圖像模型I=bJ+n和假設a和b，可以看出，圖像強度在圖像邊界附近的區域中變化很大，而在其他區域中變化緩慢。圖像強度的變化是由圖像在遠離真實物體邊界的平滑區域內的不均勻性所引起的。因此，可以使用局部區域中的方差來測量圖像不均勻程度。可以根據平滑區域中局部圖像不均勻性的程度來獲得比例參數。

在橋臺搭板進行設置期間，需要充分考慮到其坡度值，保證搭板長度的合理性。不僅要確保車輛行駛的安全性，更是要求橋體荷載力達到預期設計標準，能夠承受更多的車輛負重，這樣一來才可以將沉降控制在合理范圍之內。雖然該項施工技術的應用有著較強的科學性，但還是不能夠有效地解決橋頭跳車問題，并且還可能導致路基壓實度無法滿足相關要求，嚴重的情況下還會導致搭板出現斷裂的問題，引發嚴重的交通安全事故。所以，施工單位在施工之前需要充分的考慮沉降段的實際情況，確保搭板寬度和橋面寬度相適應，搭板邊緣在緣石邊緣的內部不能大于50cm，同時還要確保搭板與配筋設計方案的合理性，從而有效提高工程施工質量，給人們出行安全提供保障。

首先，采用中值濾波和雙邊濾波對圖像I進行預處理，以減小噪聲的影響，得到I1。然后，需要得到遠離目標邊界的光滑區域Ωh。利用拉普拉斯算子檢測目標邊界。由于大部分像素位于光滑區域，因此使用區域Ωh0粗略地表示遠離目標邊界的區域。

Ωh0=Ω(x)F(x)

(25)

式中：F=K1*I1；T=τ·min(max(|F|),T0)是閾值；T0是恒定的閾值；K1是八領域的拉普拉斯算子；τ是時間步長。

局部區域中的圖像強度變化越大，|F|的值越大。由于圖像強度在遠離物體真實邊界的區域內變化緩慢，并且大部分像素位于光滑區域，所以F(x)在區域Ωh0附近，服從高斯分布。為了去除邊界和噪聲區域，在Ωh1中更精確地表示平滑區域。

Ωh1=Ω(x) max(cF-1.5σF,-T)<
F(x)

(26)

式中：cF、σF是F(x)在Ωh0區域的均值和標準差。為了使平滑區域更加精確，應該刪除區域Ωh1的外部邊緣，并獲得平滑區域Ωh。

在獲得平滑區域Ωh之后，得到以x為中心的(2ω+1)×(2ω+1)的正方形領域Nω(x)中的圖像I1的局部方差D(x)。對于均勻圖像，平滑區域的局部方差為零，尺度參數為無窮大。尺度參數與局部圖像不均勻程度成反比。核函數的尺度參數可以設計為：

(27)

式中：β是一個常數；cD是圖像I1在區域Ωh中D(x)的平均值。從式(27)中可以清楚地看到，圖像不均勻程度越大，尺度參數越小。可以根據圖像自適應地選擇適當的比例參數。對于σb>(length(I)+width(I))/4，(length(I)表示圖像I的長度，width(I)表示圖像I的寬度),圖像強度是均勻的，可以用張氏模型對圖像進行分割。

2.4 算法實現

2.4.1算法流程圖

本文算法流程如圖1所示。

圖1 本文算法流程

2.4.2算法過程描述

該算法的主要過程概括如下：

步驟1分別用式(22)和式(27)初始化偏置場b0和σb，并將水平集函數φ0初始化為以下二元函數：

式中:Ω0是圖像閾Ω中的一個子集;?Ω0是Ω0的分界線。

步驟2K=1∶m。

(1) 利用式(18)計算ci。

(2) 根據式(30)演化水平集函數φ。

(3) 用高斯濾波器對水平集函數進行正則化。

(4) 如果φ>0，則φ=1；否則，φ=-1。

步驟3利用式(19)計算b。

步驟4檢查水平集函數是否滿足停止條件，如果不，返回步驟2。

事實上，δ(φ)可以被|▽φ|取代[13]，那么式(21)的穩態等效于：

(28)

基于觀察，水平集函數在演化過程中必須保持為符號距離函數，為了防止它在界面附近太陡或者太平[26]，水平集函數必須滿足條件|▽φ|=1。因此，式(28)可以簡化為：

(29)

為了保持數值實現的穩定性，在式(29)的迭代過程中，需要對水平集函數進行正則化。利用高斯濾波過程對水平集函數進行正則化[12]。使用加法運算符分裂(AOS)方案。水平集演化過程可以分為兩部分，第一部分是求解式(29)，另一部分是對水平集函數進行正則化。采用帶標準偏差σR的高斯濾波器對水平集函數進行正則化，并采用顯式有限差分格式對式(29)進行離散，如下：

(30)

式中：Δt表示時間步長；n表示迭代次數。

由于式(30)右邊對于φ是獨立的，可以選擇方程的相對大的時間步長。因此，所提出的模型比LIC模型更快。實際上，每次迭代時對偏置場的估計是耗時的主要因素。當偏置場固定時，僅在一次迭代后，水平集不能達到穩定狀態。因此，利用n次迭代來演化所提出的用于固定偏置場的方法中的水平集函數。

在實驗中，當水平集函數的變化保持小于給定的閾值η時，迭代將自動停止，即，如果式(31)被滿足，迭代將停止。

(31)

2.5 部分實驗編程

//初始化水平集函數

initialLSF = c0*ones(size(Img));

initialLSF(30:90,50:90) = -c0;

u=initialLSF;

hold on;

contour(u,[0 0],′r′);

title(′Initial contour′);

//初始輪廓

imagesc(Img,[0, 255]); colormap(gray); axis off; axis

equal

epsilon=1;

b=ones(size(Img));

//初始化偏置場

K=fspecial(′gaussian′,round(2*sigma)*2+1,sigma);

//高斯核

KI=conv2(Img,K,′same′);

//圖像卷積

KONE=conv2(ones(size(Img)),K,′same′);

//高斯核與常值1函數的卷積，在計算LBF時會用到

[row,col]=size(Img);

N=row*col;

for n=1:iter_outer

[u, b, C]= lse_bfe(u,Img, b, K,KONE, nu,timestep,mu,epsilon, iter_inner);

//水平集迭代

if mod(n,2)==0

pause(0.001);

imagesc(Img,[0, 255]); colormap(gray); axis off; axis equal;

hold on;

contour(u,[0 0],′r′);。

iterNum=[num2str(n), ′iterations′];

title(iterNum);

hold off;

end

end

3 實驗結果與分析

本文方法的性能在合成圖像和真實圖像上得到了廣泛的評價。所有實驗都是在MATLAB R2016a編程環境下進行的，使用的計算機型號為Lenovo G40，其配置為CPU 2.10 GHz、內存8 GB和Windows 10操作系統。實驗中使用的圖像像素各不相同，圖2中第一幅像素為84×84,第二幅為98×100，第三幅為98×98，第四幅是MATLAB自帶的一幅圖像。圖3中圖像來自蘭州第一人民醫院眼科影像資源庫，其像素為490×480。圖4是在加入噪聲的情況下的分割結果。圖5(d)是合成圖像，像素為103×93。經過多次實驗，最終選擇參數為：β=5.2，ω=1，T0=5，σR=0.9，Δt=1.0，η=1，m=5。其中時間步長可以取得比傳統水平集方法大一些，本文實驗過程中就分別取到了從0.1到100的寬幅范圍。使用大的時間步長雖然可以加快迭代速度，但是選得太大可能會在邊緣的地方引起錯誤。即大的時間步長和精確的邊緣定位是矛盾的，通常情況下選擇時間步長都小于10。

(a) 初始輪廓的原始圖像 (b) 最終分割結果圖2 本文方法在四幅強度不均勻圖像上的實驗結果

(a) 原始圖像(b) 最終分割結果(c) 估計偏置場 (d) 偏差校正圖像初始輪廓 圖3 本文方法在兩幅強度不均勻的真實圖像中的應用

(a) 泊松噪聲(b) 0.01的椒鹽噪聲(c) 0.03的椒鹽噪聲圖4 加噪聲圖像分割結果

圖5 三幅真實圖像和一幅強度不均勻的合成圖像的分割

3.1 實驗結果

在這一部分中，將本文模型與CV模型、LBF模型、LIC模型在強度不均勻的合成圖像和真實圖像上進行了比較。

圖2為本文方法對四幅強度不均勻圖像的分割。在這四幅圖像中，強度不均勻性明顯。結果表明，本文方法能夠對灰度不均勻的圖像進行分割，并取得了滿意的分割效果。

與LIC模型相似，本文模型也能夠分割灰度不均勻的圖像，并同時估計偏置場。該方法應用于兩幅強度不均勻的真實醫學圖像，如圖3所示。可以看出，該模型得到了準確的分割結果，并且在偏置校正圖像中，輪廓內部和外部區域的圖像強度變得相當均勻。實驗結果表明，該模型在圖像分割和偏置場估計方面具有良好的性能。

圖4是在加入噪聲的情況下本文算法的圖像分割結果。第1列的圖像中加入了泊松噪聲，第2列的圖像中加入了密度為0.01的椒鹽噪聲，第3列的圖像中加入了密度為0.03的椒鹽噪聲，第2行分別是各種噪聲下的圖像分割結果。實驗結果表明，在加入上述各種噪聲的情況下，本文算法仍可以較好地分割目標圖像。

3.2 實驗分析

為了進一步驗證本文算法的圖像分割性能，并與其他算法進行比較，本文采用Jaccard Similarity(JS)等圖像分割性能指標[25-28]，利用下面兩個圖像分割性能評價指標進行評價：

1) 面積重疊誤差。

(32)

式中：A和M分別表示水平集算法分割區域以及手動分割區域的像素數量。

當面積重疊誤差指標為零時，表示最精確分割，為100%時表示算法分割的結果與真實結果完全不重合，即錯誤分割。

2) 邊界平均距離。設C是算法分割出的目標邊界，對于邊界C上的每一點Pi(i=1,2,…,N)計算Pi與真實目標邊界S的距離dist(Pi,S)，然后求平均值。

(33)

如圖5和表1分別給出了本文算法與其他經典算法的實驗對比結果。在圖5中給出了利用CV模型、LBF模型、LIC模型，以及本文方法對三幅真實圖像和一幅強度不均勻合成圖像的分割結果, 對于最后一幅圖像，σR設置為0.5。CV模型采用全局強度信息的均值，如圖5(b)所示，無法分割出強度不均勻的圖像。LBF模型利用局部強度的均值來擬合被測圖像，從而可以對灰度不均勻的圖像進行分割。然而，該模型在能量泛函中引入了許多局部極小值，在某些情況下導致輪廓卡在不正確的位置上，如圖5(c)所示。圖5(d)和圖5(e)分別是LIC模型與本文算法分割結果，因為它們在水平集演化過程中分割圖像并同時估計了偏置場，因此比其他兩種模型對于灰度不均勻圖像的分割效果更好。但LIC模型對于灰度嚴重不均勻的合成圖像的分割效果不是很理想，而本文模型不管是對真實圖像還是對灰度嚴重不均勻的合成圖像都取得了很好的分割效果。

表1 各算法的迭代次數、運算時間、面積重疊誤差、邊界平均距離對比

可以看出，本文算法在面積重疊誤差、邊界平均距離、迭代次數、運算時間上均優于其他幾種算法。

為了進一步分析新偏置場初始化的影響，給出LIC模型與本文模型對于不同的初始輪廓分割結果的比較，如圖6所示。結果表明，當初始輪廓位于物體內部時，LIC模型得到了正確的分割結果。然而，當初始輪廓同時包含目標和背景時，LIC模型無法對目標進行分割。而本文模型不管輪廓在什么位置，都取得了很好的分割結果，因此本文模型對初始輪廓具有較強的魯棒性。

(a) 原始圖像及初始輪廓 (b) LIC模型 分割結果 (c) 本文模型 分割結果 圖6 LIC模型與本文模型對于初始輪廓位置的比較

4 結 語

本文提出一種新的水平集分割方法。在非均勻圖像模型的基礎上，推導出圖像域的最優分割平面。在平面上，提出一種新的基于區域的壓力函數，并在整個圖像區域的水平集公式中定義了一個能量泛函。通過對能量泛函進行最小化來實現圖像分割和偏置場估計。此外，引入了新的偏置場初始化以提高對初始輪廓的魯棒性，并且為核函數設計了新的自適應尺度參數來準確地估計偏置場。在人工合成圖像與真實圖像上的仿真實驗結果表明，該方法能夠迅速、準確地對非均勻圖像進行分割，并且對初始輪廓具有較強的魯棒性。該方法分割結果能夠輔助人工判讀，用于輔助臨床實習教學，對診斷提供一定參考信息。近幾年許多學者基于水平集的圖像分割方法做了大量的工作，該方法已經成功用于醫學圖像中病灶區域的分割，比如對于心臟的分割，大腦灰白質的分割，未來它還將會被應用于軍事、航天、人臉識別等更多領域。