基于邊力函數能量泛函的圖像分割算法

唐銀敏 李自納 胡 靜

1(鄭州大學信息工程學院 河南 鄭州 450000) 2(南陽農業職業學院信息工程學院 河南 南陽 473000)

0 引 言

圖像分割在醫學、定位、計算機視覺等領域發揮了極其重要的作用[1]。然而,由于低分辨率、高噪聲和模糊的邊界,它具有極大的挑戰性。隨著網絡時代的不斷發展,實時性以及準確性要求使得實現快速高精度的圖像分割成為研究發展的熱點[2]。

活動輪廓模型(Active Contour Model,ACM)的基本思想是利用輪廓曲線提取目標,輪廓曲線包括兩類:顯式ACM和隱式ACM[3-4]。顯式ACM在求解過程中使用參數方程顯式表示演化曲線[5]。隱式采用有符號距離函數(Signed Distance Functions,LSF)代替參數方程來表示演化曲線[6]。實際上,噪聲或弱邊界往往存在于真實圖像中。由于噪聲的影響和灰度值的不連續性,標準的基于邊緣的最小二乘法無法檢測出圖像的邊界。因此,Li等[7]提出了一種距離正則化水平集演化方法,以保持LSF的正則性,并將等高線曲線推向所需的位置。基于區域的最小二乘法利用圖像的區域統計信息來定位目標邊緣。比較基于邊緣的最小二乘法,這些模型對噪聲不敏感,例如區域可縮放擬合模型[8]、局部圖像擬合模型[9]和魯棒噪聲區域ACM[10]。特別的,Chan等[11]提出了一種基于區域的Chan-Vese模型,它能成功地分割出亮度均勻的圖像。局部圖像擬合模型是另一種基于區域的LSM模型,它以截斷高斯函數作為局部區域約束來分割圖像,可以分割出強度不均勻的圖像,但對強噪聲圖像敏感,RSF模型也能分割出非均勻圖像,但當圖像是嚴重噪聲或非均勻時,分割效果較差,而且該模型會使集函數過于平坦,不利于水平集函數的正則性。

為解決上述方法中存在的問題,提出一種基于有序統計濾波能量驅動的魯棒主動輪廓模型。通過實驗結果分析可知提出的方法能夠更加快速、準確地分割圖像,并且對初始輪廓和參數具有較強的魯棒性。

1 基礎理論

1.1 幾何主動輪廓模型

幾何主動輪廓(Geodesic Active Contour,GAC)模型是一種基于邊緣指示函數的主動輪廓模型[12]。該模型提出了如下能量函數:

(1)

式中:Eint代表內部能量項,即長度項;Eext代表外部能量項,即面積項。

Eint(C)=α|C′(q)|2

(2)

Eext(C)=λg|▽I(C′(q))|2

(3)

式中:α和λ是常數;▽表示梯度算子;g是邊緣指示函數。g定義如下:

(4)

式中:Gσ是具有標準差σ的高斯濾波函數;*表示卷積。可以得出梯度下降流方程:

(5)

1.2 距離正則化水平集演化模型

距離正則化水平集演化模型(Distance Regularized Level Set Evolution,DRLSE)的能量函數如下:

(6)

式中:μ、λ和α是常數,分別是距離正則項、長度項和面積項的系數。式(6)中的g表示邊緣指示函數:

(7)

式中:Gσ是具有標準差σ的高斯濾波函數。式(6)中的p表示雙阱勢函數,可以寫成:

(8)

Hε和δε函數表達式為:

(9)

(10)

可以通過變化的微積分獲得以下梯度下降流方程:

αgδε(φ)

(11)

式中:dp(x)=p′(x)/x表示演化速度函數;p′(x)表示雙阱勢函數的導數函數。

1.3 CV模型

設I:Ω→R為輸入圖像,C為閉合曲線,提出了如下能量函數:

(12)

式中:λ1、λ2和v為常數;outside(C)代表輪廓C以外的區域;inside(C)代表輪廓C內的區域;c1和c2是兩個常數,分別近似于區域outside(C)和inside(C)中的強度平均值。式(12)中的前兩項是數據項,它們用于驅動曲線停止在目標邊界處。第三項是長度項,用于平滑和縮短曲線。從式(12)可知,只有當輪廓C位于目標邊界時,能量ECV才能達到最小值。

在水平集方法中,輪廓C?Ω由Lipschitz函數φ:Ω?R的零水平集表示,稱為水平集函數,水平集函數如式(13)所示。當x點分別位于曲線C的外部和內部時,則水平集函數φ分別為正和負;否則,水平集函數φ等于零。

(13)

能量函數ECV重寫如下:

(14)

式中:Hε(x)和δε(x)分別是近似Heaviside函數和近似Dirac函數[13]。

(15)

(16)

使用標準的梯度下降方法來最小化式(14)中的能量函數,可以得到以下梯度下降流方程:

(17)

其中c1和c2是:

(18)

最后,用迭代法φk+1=φk+Δt×(?φ/?t)求出水平集函數φ。c1和c2涉及圖像強度的全局屬性,它們不包含任何局部信息,因此,在對強度不均勻的圖像進行分割時,分割結果將是錯誤的。

2 順序統計濾波模型

2.1 二階微分特性

GAC模型和DRLSE模型通常使用面積項αgδε(φ)來加速零水平集的演化速度。但由面積項驅動的演化過程是單向的,這意味著水平集函數不能自適應地選擇演化方向,只能根據設定方向向內收縮或向外擴張,而不能根據圖像信息調整演化方向。因此,本文的目的之一是尋找一種既能準確定位邊界位置又能指導曲線演化方向的邊緣信息。

以CV模型為例,CV模式下的數據擬合項為-[λ1|I(y)-c1(x)|2-λ2|I(y)-c2(x)|2],c1和c2分別代表輪廓C外部和內部區域中強度的平均值。并且設置λ1=λ2=1,則CV模型中的數據擬合項可以重寫為FCV(φ,c1,c2)=-(λ1|I-c1|2-λ2|I-c2|2)=(I-c2)2-(I-c1)2。I(x,y)表示圖像中的一個點,因此(I-c1)2表示c1和點I(x,y)之間的差的平方,而(I-c2)2表示c2和點I(x,y)之間的差的平方。(I-c1)2和(I-c2)2都代表一次微分計算。之后,(I-c2)2-(I-c1)2表示二次微分計算。因此,數據擬合項FCV(φ,c1,c2)也可以稱為二階微分數據項。二階微分的零點是灰度值劇烈變化的地方,即目標邊界的位置。

2.2 邊力函數

在本文中,利用順序統計濾波來構造具有二階微分特性的邊力函數(Edge Force Function,EFF)。令Ω→Rd為圖像域,當d=1時,它代表灰度圖像,而d=3時,它代表彩色圖像。首先,定義一個以x為中心、寬度為ω的正方形鄰域Ωx。

其次,通過使用MATLAB中的二維順序統計濾波器,將局部窗口Ωx中的像素值從低到高排序,然后獲取圖像強度的最大值和圖像強度的最小值分別為f1和f2。

f1(x)=max[I(y)|y∈Ωx]f2(x)=min[I(y)|y∈Ωx]

(19)

式中:I(y)表示局部窗口Ωx中的圖像強度值。對于寬度為ω的某點x,可以直接根據式(19)計算f1(x)和f2(x)的值。

然后,介紹了以下能量函數:

(20)

式中:C是閉合輪廓;inside(C)和outside(C)分別表示輪廓C內部和外部的區域;λ1和λ2是常數。當式(20)中的能量函數被最小化時,輪廓C將位于邊界上。

曲線C可以由Lipschitz函數φ的零水平集代替,因此,式(20)中的能量函數可以重新定義如下:

(21)

δε(φ)在式(22)中是Hε(φ)的導數。Hε(x)和δε(x)表示如下:

(22)

通過用最速下降法最小化式(21)中的能量函數,可以得到梯度下降流方程[14]:

(23)

式中:α是可以控制OSF模型分割速度的正常數;參數s表示為控制OSF模型分割范圍而增加的圖像的標準差。ei(x)表示如下:

(24)

然后,將EFF定義如下:

fEFF(x)=α(λ1e1(x)-λ2e2(x))

(25)

演化水平集函數:

φn+1=φn+Δt·(?φ/?t)

(26)

式中:Δt是時間步長。(Ai+1-Ai)/Ai<10-5是曲線演化的停止標準,其中Ai和Ai+1分別表示迭代前后輪廓所包圍的區域。

最后,使用優化的正則化函數φL和優化的長度函數φC分別對水平集函數進行正則化并平滑和縮短曲線。

(27)

φC=Gk*φL

(28)

式中:φL是優化的正則化函數,其目的是改善過零區域的斜率,并限制兩端的斜率,以確保水平集函數具有很大的正則性;參數β用于控制過零區域的斜率,以確保分割效率;Gk是具有標準偏差k且大小為round(2×k)×2+1的高斯核。這兩個函數的加入不僅減少了計算量,而且消除了長度和懲罰項的系數。

該模型的優勢體現在以下幾個方面:首先,由于f1(x)和f2(x)與大小受Ωx控制的局部區域中的圖像強度有關,因此,OSF模型可以分割強度不均勻的圖像。其次,由于EFF的二階微分特性,曲線可以在演化過程中自適應地選擇演化方向。第三,EFF是在迭代之前計算的,因此,OSF模型僅需要在迭代期間更新δε(φ),即可大大提高分割速度。第四,OSF模型大大提高了初始化和參數的魯棒性。最后,由于添加了優化的正則化函數φL和優化的長度函數φC,極大地降低了計算復雜度,并且消除了長度和懲罰項的系數。

3 實驗與結果分析

3.1 實驗設置

OSF模型的步驟如下:

步驟1設置各種參數并初始化水平集函數。

步驟2根據式(19)計算f1(x)和f2(x)。

步驟3根據式(24)和式(25)計算EFF。

步驟4根據式(26)更新水平集函數。

步驟5分別根據式(27)和式(28)計算φL和φC。

步驟6如果達到迭代次數,則停止,并輸出分割結果,否則返回步驟4。

在OSF模型中使用以下參數:Δt=0.1,ε=1,k=2.5,β=10,α=1.5,ω=7,λ1=λ2=1,C0=1。初始水平集函數φ0初始化為二進制步進函數,在零水平集外部和內部分別取1和-1。OSF模型是在2.6 GHz Inter Core i5個人計算機上的MATLAB R2017b中實現的。

3.2 初始輪廓的魯棒性

在本節中,將選擇五幅強度不均勻或邊界較弱的不同圖像,以證明OSF模型對初始輪廓具有魯棒性。對于每幅圖像,用三種不同的形狀設置初始輪廓,分別是矩形、圓形和三角形。每個形狀的初始輪廓在大小和位置上都不同。

初始輪廓和分割結果如圖1所示。因為EFF是在迭代之前計算的,它與輪廓無關。因此,可以選擇不同的初始輪廓,在不同的位置嵌入、包含或遠離目標。

圖1 不同形狀、大小和位置的初始輪廓的分割結果

在圖1中,白虛線代表初始輪廓,可以看到初始輪廓具有不同的位置、形狀和大小。灰線代表最終的分割結果,可以看到,在初始輪廓具有不同形狀、大小和位置的情況下,OSF模型可以精確地分割每幅圖像。因此,可以得出一個結論,即OSF模型對初始輪廓具有魯棒性,能清楚顯示出不同形狀、大小位置的輪廓,分割結果明確,OSF模型的魯棒性效果明顯。

3.3 分割實驗

在本節中,將利用十幅強度不均勻、邊界弱、背景不均勻或對比度低的典型圖像來進行OSF模型的分割實驗。初始輪廓、演化過程和最終分割結果如圖2所示。在圖像d和j中,設置β=9,在圖像g中,設置α=1.0只是為了獲得更好的分割效果。

圖2 圖像的分割結果

在圖2中,規則的線代表初始輪廓,不規則的線代表演化曲線。在每個組中,第一列代表原始圖像和初始輪廓。第二和第三列表示演化過程,最后一列表示準確且令人滿意的最終分割結果。因此,可以得出結論,OSF模型可以正確地分割強度不均勻或邊界較弱的圖像。

3.4 對比實驗

在本節中,將比較OSF模型的分割結果與一些經典的基于區域模型的分割結果,包括RSF模型、LIF模型、RSF&LoG模型和LPF模型。選擇圖3中的圖像k-圖像r來證明OSF模型的良好分割速度和分割精度。圖3顯示了具有相同初始輪廓的OSF模型、RSF模型、LIF模型、RSF&LoG模型和LPF模型之間的比較。根據花費的時間和每個模型的分割精度,將評估每個模型的性能。并根據戴斯相似系數(Dice Similarity Coefficient,DSC)[15]估算結果:

圖3 對比實驗結果

(29)

式中:ST表示通過手動獲得的真實區域;SE表示通過實驗獲得的目標區域。DSC的值越接近1,則分割精度越高。分析結果列于表1。

表1 每個圖像的迭代次數、花費的時間和DSC

在圖3中,第一列代表原始圖像和初始輪廓。第二到最后一列分別代表RSF模型、LIF模型、RSF&LoG模型、LPF模型和OSF模型的分割結果。表1顯示了每幅圖像的迭代次數、花費的時間和DSC。

在表1中,可以注意到OSF模型花費的時間比其他模型少得多,因為EFF是在曲線演化以及優化的正則項和長度項之前計算的。RSF模型和RSF&LoG模型花費的時間相對較長,因為它們每次迭代都具有四個卷積。LPF模型的分割時間也很長,這是因為要計算預擬合函數fs(x)和fl(x)。由于每次迭代的計算成本低,LIF模型的分割時間相對較短。此外,所提出模型的DSC在圖像k、l、m、o、p、q、r中排名第一,在圖像n中排名第二。因此,可以得出一個結論:與其他四個基于區域的經典模型相比,OSF模型可以更準確、更快地分割圖像。

4 討 論

4.1 邊力函數

圖4(a)顯示了圖像的EFF的計算結果,可以看出圖中有明顯的黑白環,黑環和白環相交的位置是目標邊界所在的位置;圖4(b)顯示了橫截面上具有二階微分特性的邊緣力能量EEFF的值。可以注意到,邊緣力能量的值在目標邊界外部和附近的區域為正,而邊緣力能量的值在目標邊界內部和附近的區域為負。下行中所示的二階微分的零點是點A、B、C、D,它們表示目標邊界所在的位置。因此,通過以上分析可以知道,可以通過計算EFF的二階微分零點來檢測目標邊界。

圖4 EFF的計算結果和邊緣力能量的值

4.2 自適應演化曲線

EFF在本文中的另一個作用是吸引演化曲線以自適應地演化到目標邊界。圖5顯示了一次迭代后的合成圖像的分割結果。規則的線表示初始輪廓,不規則的線表示一次迭代后的輪廓。如圖5所示,點a在目標內部,位于初始輪廓上。因此,根據圖5可以知道點a的邊緣力能量值為負,并且根據梯度下降流方程,點a應該移動到能量較高的區域。并且,由于能量值在初始輪廓之外的區域中為正,因此,點a向外擴展到點a1,移動速度為|?φ/?t|Δt,其中Δt是時間步長。

圖5 曲線的演化方向

點b也在初始輪廓上,但是它位于目標邊界之外。因此,根據圖5,可以知道點b的邊緣力能量值為正,因此根據梯度下降流方程,點b應該移動到能量較低的區域。并且,因為能量值在初始輪廓內的區域中為負,所以點b向內收縮到點b1,移動速度為|(?φ/?t)|Δt。

點c位于目標邊界上,因此根據圖5,可以知道點c的邊緣力能量的值等于零,即|(?φ/?t)|=0,則運動停止。

4.3 參數β、k、α和ω的影響

4.3.1 關于參數β

式(27)對本文中的水平集函數進行了正則化。利用式(27)中的參數β控制過零區域的斜率,以保證水平集函數的正則性。首先選擇圖6中的圖像s、t、u、v、w、x來顯示參數β的有效性。圖6中的第一行代表原始圖像和初始輪廓;第二行表示沒有參數β的分割結果;第三行表示有參數β的最終分割結果。可以明顯地注意到,沒有參數β的最終分割結果是錯誤的,而具有參數β的分割結果是正確的。

圖6 在沒有參數β和有參數β的情況下獲得的分割結果

在分割實驗中,參數β在圖像d和j中進行了微調,以實現更好的分割。因此,在本節中,將選擇圖6中的圖像v來討論參數β的魯棒性。圖7(a)顯示了原始圖像和初始輪廓。圖7(b)-圖7(h)表示當參數β等于7到13時獲得的最終分割結果。當參數β等于8到12時,最終的分割結果令人滿意。因此,可知OSF模型對參數β具有較強魯棒性。

圖7 具有不同參數β值的分割結果

4.3.2 關于參數k

式(28)用于平滑和縮短演化過程中的曲線。但是,高斯濾波會模糊圖像的邊界,并會導致丟失少量信息。因此,式(28)中的參數k應根據圖像信息設置。在本節中,將設置k=2.0、k=2.5、k=3.0,以證明在OSF模型中選擇的k=2.5是最合適的。在圖8中,選擇圖6中的圖像u、v、w、x和圖像y進行比較實驗。第一行表示原始圖像和初始輪廓。第二行到最后一行分別代表k=2.0、k=2.5、k=3.0的最終的分割結果。表2記錄了每幅圖像的迭代次數和花費的時間。

圖8 使用參數k的不同值進行分割的結果

表2 不同參數k下圖像u到y的迭代次數和花費的時間

從圖8和表2中可以發現,當k=2.0時,圖像v和y均未成功分割。同時,當k=3.0時,圖像v和w分割失敗。可以使用參數k的不同值成功分割圖像u和x。從表2可以看出,當k=2.5時,迭代次數和花費的時間最少。當k=2.0和k=2.5時,可以成功分割圖像w。從表2可以看出,當k=2.5時,迭代次數和花費的時間都更少。類似的,當k=2.5和k=3.0時,可以將圖像y分割,當k=2.5時,迭代次數和花費的時間會更少。因此,可以得出結論,在OSF模型中選擇k=2.5是最有利的。

4.3.3 關于參數α

式(23)中使用的參數α可以控制OSF模型的分割速度。在分割實驗中,為了更好地分割,在圖像g中對參數α進行了微調。因此,選擇圖8中的圖像x來證明OSF模型對于參數α是魯棒的。

如圖9所示,第一列表示原始圖像和初始輪廓,而第二列到最后一列顯示最終的分割結果,其中α=0.5、α=1.5、α=2.5、α=3.5、α=4.5和α=5.5。可以看到,只有當α=0.5時,OSF模型才不能正確分割圖像。相反,當α=1.5、α=2.5、α=3.5、α=4.5和α=5.5時,最終的分割結果都是正確的。因此,可以得出一個結論,即OSF模型對參數α具有魯棒性。

圖9 不同參數α值的分割結果

4.3.4 關于參數ω

局部窗口Ωx中的參數ω決定了所提出的模型一次要處理的像素數。因此,參數ω可以在一定程度上控制所提出模型的分割速度。以圖像y為例,通過更改參數ω的值并保持其他參數固定來說明參數ω的有效性。如圖10所示,當ω=4時,最終的分割結果沒有捕獲到完整的目標。相反,當ω=8時,最終分割結果包括非目標部分。但是,當ω=5、ω=6和ω=7時,最終結果都是正確的,這表明所提出的模型對參數ω具有魯棒性。此外,當ω等于5、6和7時,表3中記錄了所提出模型所需的迭代次數和分割時間。從表3中可以發現ω=5、ω=6和ω=7都可以得到理想的結果,但是ω的值越大,所提出的模型所需的迭代次數和分割時間就越少,這可以證明所提方法中使用的參數ω可以控制所提模型的分割速度。因此,通過以上分析,可以得出關于參數ω的以下結論:(1) 所提出的模型在一定范圍內對參數ω具有魯棒性。(2) 所提出的方法中使用的參數ω會影響所提出模型的分割速度,這是因為參數ω可以一次控制所提出的模型要處理的像素數量,因此ω的值越大,則一次處理更多像素。

圖10 不同參數ω值的分割結果

表3 不同參數ω下圖像y的迭代次數和花費的時間

5 結 語

針對傳統圖像分割算法中存在的圖像分割速度慢以及初始輪廓和參數魯棒性差等問題,提出一種基于邊力函數能量泛函的圖像分割算法。通過實驗結果分析可得如下結論:

(1) 提出的OSF模型能夠快速、準確地獲取目標邊界,實現精確的圖像分割。

(2) 由于能量函數中的f1和f2僅與圖像強度有關,所以OSF模型對初始輪廓具有較強的魯棒性。

(3) 正則項的引入大大地減少計算量,并消除了長度系數和懲罰項系數,使得算法對參數選擇也具有一定的魯棒性。