徐思敏,金正猛,閔莉花,王 皓,郭小亞

(南京郵電大學 理學院,江蘇 南京 210023)

0 引 言

圖像分割是根據(jù)某種均勻性或一致性原則把圖像分成具有特定性質(zhì)區(qū)域的過程,是圖像處理和計算機視覺領(lǐng)域的基本任務(wù)之一,在機器視覺、醫(yī)學成像、自動駕駛、對象檢測和交通控制系統(tǒng)等領(lǐng)域有廣泛應(yīng)用。在自然或醫(yī)學圖像成像過程中,由于光照或設(shè)備偏移場的影響,圖像會出現(xiàn)灰度異質(zhì)問題,從而嚴重影響圖像分割的精度。該文主要關(guān)注針對灰度異質(zhì)圖像的分割方法。

近年來,活動輪廓圖像分割模型備受關(guān)注?，F(xiàn)有的活動輪廓模型可分為：基于邊緣的分割模型和基于區(qū)域的分割模型。基于邊緣的分割模型主要依賴圖像的邊緣信息捕捉目標物體的邊界,如：著名的Snake模型[1]和測地活動輪廓模型(Geodesic Active Contour,GAC)[2]等。基于區(qū)域的分割方法,主要利用圖像中物體呈分片區(qū)域的特點來引導初始輪廓的演化,如：經(jīng)典的Mumford-Shan(MS)模型[3]。由于MS模型的求解過于復(fù)雜,Chan等人[4]將不同區(qū)域的灰度值近似為分片常數(shù),并用水平集函數(shù)[5]來表示輪廓線,提出經(jīng)典的Chan-Vese(CV)模型。近年來,為了提高邊緣處的分割精度,Bresson等人[6]將GAC模型中的測地輪廓長度項引入CV模型中,提出基于測地輪廓的CV分割模型(簡稱g-CV模型),該模型能更好地捕捉邊緣信息,提高對目標物體的分割精度。由于CV和g-CV模型都是在圖像灰度均勻的前提下建立的,對分割灰度異質(zhì)圖像,其分割精度難以得到保證。

為了有效分割灰度異質(zhì)圖像,Zosso等人[7]根據(jù)Retinex理論[8-9]將圖像分解為光照部分和反射部分之和,提出CVB模型用于圖像分割和偏移場矯正。該模型能較好分割灰度異質(zhì)圖像,但仍然存在以下問題：一方面,CVB模型沒有整合目標物體邊緣信息,容易出現(xiàn)過度分割或欠分割問題。另一方面,文獻[7]使用結(jié)合交替極小化[10]和閾值動力學[11-12]的算法求解CVB模型,計算過程中需要求解復(fù)雜輪廓演化的偏微分方程,求解速度并不可觀。

為了提高CVB模型的分割精度,該文在CVB模型的基礎(chǔ)上引入測地輪廓長度項,提出新的變分分割模型(簡稱g-CVB模型)。此外,受輪廓長度近似公式[13]的啟發(fā),該文提出用特征函數(shù)表示測地輪廓長度項,并結(jié)合交替極小化算法和迭代卷積閾值法設(shè)計新模型的求解算法。該算法的優(yōu)勢體現(xiàn)在：(1)可以通過簡單的迭代卷積閾值法進行交替求解,計算效率高。(2)在交替求解過程中,每個子問題的求解都是穩(wěn)定的。最后,實驗結(jié)果表明：該算法不僅能有效分割灰度異質(zhì)圖像,其收斂速度還有明顯提升。

1 相關(guān)工作

文中,Ω?R2表示具有Lipschitz邊界的有界圖像域,I：x∈Ω為輸入圖像,Γ為封閉邊界曲線,*表示卷積運算符,?表示梯度算子。

1.1 基于特征函數(shù)的GAC模型

傳統(tǒng)GAC模型使用水平集表示能量泛函,在求解時存在計算效率低和數(shù)值不穩(wěn)定的問題。為了解決上述問題,Ma等人[13]使用特征函數(shù)表示能量泛函,結(jié)合迭代卷積閾值法(Iterative Convolution Threshold Method,ICTM)求解模型,其方法概括如下。定義特征函數(shù)u(x)：

(1)

其中,Γ表示待分割對象的邊界,ΩΓ表示Γ內(nèi)部的區(qū)域。根據(jù)文獻[14],目標的測地線長度近似表示為：

(2)

其中：τ>0,邊緣檢測函數(shù)為：

(3)

其中,γ為大于0的參數(shù),Gτ和Gσ均為高斯核函數(shù)。

由測地線長度項結(jié)合面積項得到基于特征函數(shù)的GAC模型的能量泛函近似為：

文獻[13]中Ma等人基于上式使用迭代卷積閾值法進行求解,相較于傳統(tǒng)水平集方法,該算法求解速度更快,收斂更穩(wěn)定。

1.2 CVB模型

2017年,Zosso等人[7]將Retinex理論應(yīng)用到CV圖像分割模型中,提出CVB模型。該模型不僅能有效分割灰度異質(zhì)圖像,且能較好地對輸入圖像進行偏置矯正。

設(shè)i(x,y)是灰度異質(zhì)圖像,Retinex理論將其分解為光照部分b(x,y)和反射部分s(x,y)的乘積：i(x,y)=b(x,y)·s(x,y),兩邊同時進行對數(shù)變換得：

logi(x,y)=logb(x,y)+logs(x,y)

(4)

令I(lǐng)=log(i),B=log(b),S=log(s),則式(4)可以簡化為：

I=B+S

(5)

其中,B被認為是光滑函數(shù),表示圖像中光照偏移場部分,S被認為是分片常數(shù)函數(shù),表示反射部分。假設(shè)輸入函數(shù)I滿足式(5),得到如下分割灰度異質(zhì)圖像的CVB模型：

其中,λ1,λ2,α,β>0為參數(shù)。使用水平集函數(shù)φ的零水平來表示Γ,引入Heaviside函數(shù)[5]H(φ),令u：=H(φ)且凸松弛為u：Ω→[0,1],得到下式：

s.t.I=B+S

模型第一項是偏移場B的平滑項,第二、三項是數(shù)據(jù)保真項,最后一項是長度正則項。

在文獻[7]中,Zosso等人結(jié)合交替極小化和閾值動力學算法來求解CVB模型,并結(jié)合相場方法[11]和MBO方案[12]設(shè)計基于閾值動力學的算法來求解關(guān)于u的子問題,但因為該算法需要求解復(fù)雜的輪廓演化的偏微分方程,求解速度較慢。此外,上述CVB模型中的長度正則項是以分割結(jié)果的長度項最小為目標,在分割尖銳邊界時容易出現(xiàn)欠分割問題。

2 文中模型及算法

2.1 文中模型與算法

為了充分利用圖像邊緣信息,該文利用GAC模型中的測地長度項取代CVB模型中原有的長度項,提出如下的g-CVB模型：

s.t.I=B+S

(6)

用式(1)中的特征函數(shù)u(x)隱式表示Γ,可將g-CVB模型轉(zhuǎn)化為：

s.t.I=B+S

(7)

其中,τ、λ1、λ2、α為正參數(shù),M：={u∈BV(Ω,R)|u={0,1} },BV(Ω,R)為有界變差函數(shù)空間[15]。

結(jié)合交替極小化和迭代卷積閾值法,該文設(shè)計一種高效的數(shù)值求解算法。將式(7)轉(zhuǎn)換為無約束極值問題,其對應(yīng)的增廣拉格朗日函數(shù)為：

(8)

其中,ρ是罰參數(shù),θ是拉格朗日乘子。首先固定c1,c2：

接下來,分別對各子問題進行求解：

(1)求解關(guān)于u的子問題。

這里記：

(10)

注意到能量泛函最小化問題式(10)的可行集M是非凸的,直接求解式(10)是困難的,因此將M松弛到它的凸包K：= {u∈BV(Ω,R)|u∈[0,1] }上,得到問題式(10)松弛后的極小化問題：

(11)

該文在2.2節(jié)中給出求解式(10)與式(11)的等價性證明。

接下來,通過迭代卷積閾值法來求解關(guān)于u的問題式(11)。很容易證明,Eτ(u)是一個凹泛函,而凹泛函的圖像總是低于它的線性逼近,故可將求Eτ(u)的最小值問題近似等價于求Eτ(u)的線性逼近的最小值問題[17]。具體來說,計算Eτ(u)在第k次迭代uk處的一階泰勒展開式：

其中：

繼而,通過求解下列線性化問題得到k+1次迭代uk+1：

對于?x∈Ω,可解出：

由于凸集上的線性泛函的最小值必在邊界處達到,因此：

(2)求解關(guān)于B的子問題。

〈θk,I-B-Sk〉Ω

其對應(yīng)的歐拉-拉格朗日方程為：

ρBk+1-2αΔβk+1=θk+ρ(I-Sk)

通過快速傅里葉變換,求解得：

(3)求解關(guān)于S的子問題。

|I-Bk+1-S|2dx+〈θk,I-Bk+1-S〉Ω

解得：

其中：

Q1=2[λ1uk+1c1+λ2(1-uk+1)c2]+

ρ(I-Bk+1)+θk

Q2=2[λ1uk+1+λ2(1-uk+1)]+ρ

(4)更新θ。

θk+1=θk+ρ(I-Bk+1-Sk+1)

依次迭代上述求解過程直到滿足收斂條件。具體流程如算法1所示。

2.2 等價性證明

下面的引理1證明了求解式(10)與式(11)的等價性。

引理1：求解原問題式(10)與求解其松弛后的問題式(11)是等價的,即若u*是式(10)的解,則它同樣是式(11)的解,反之亦然。

對于|t|

二階導為：

算法1：文中算法的求解流程。

1.初始賦值：S0=I,B0=(0,0),θ0=0,u0取決于初始輪廓.

2.迭代：更新uk+1,Bk+1,Sk+1如下：

2.3 穩(wěn)定性分析

在本節(jié)中,定理1證明了對于任意的τ>0,Eτ(u)在迭代過程中逐漸減小。即算法1是無條件穩(wěn)定的。

定理1 令uk(k=0,1,2,…)為算法1中u的第k次迭代,對任意τ>0,有Eτ(uk+1)≤Eτ(uk)。

證明：對式(9)做線性化處理可得：

將uk帶入Eτ(u)得到：

進一步,

因為uk+1是序列線性規(guī)劃的一個解,故：

Lτ(uk+1,uk)≤Lτ(uk,uk)

進而：

Eτ(uk+1)-Eτ(uk)=Lτ(uk+1,uk)-Lτ(uk,uk)+ζ其中：

可得：

因此Eτ(uk+1)-Eτ(uk)≤0。證畢。

3 實驗結(jié)果及分析

本節(jié)中,對多幅不同類型的灰度異質(zhì)的圖像進行分割,并同CV模型[4]、Cai模型[18]、CVB模型[7]、ICTM-GAC算法[13]以及ICTM-CV算法[19]的分割結(jié)果進行對比,以檢驗文中方法的有效性。為了定量地評價不同方法的分割性能,這里采用Dice相似系數(shù)(Dice Similarity Coefficient,DSC)[20]和Hausdorff距離(Hausdorff Distance,HD)[21]作為評價指標。它們的定義分別是：

(12)

HD(A,B)：=

(13)

其中,S1代表分割后目標物體的區(qū)域,S2代表對應(yīng)的真實區(qū)域,A為分割后的二值圖像,B是相應(yīng)的真實二值圖像,n表示分割目標邊界集的總數(shù)。DSC值越高,HD值越低,表示分割結(jié)果越精確。在提出的模型(7)中,設(shè)置參數(shù)τ=0.4,λ1=λ2,平滑參數(shù)α=9 000,不同圖像對應(yīng)的λ1、λ2的值列于表1中。在算法1中固定懲罰參數(shù)ρ=92,設(shè)置停止準則tol=10-5。

表1 參數(shù)λ1、λ2的取值

3.1 合成圖像與MR圖像的數(shù)值實驗結(jié)果

首先,圖1展示了不同模型與算法對三幅合成圖像的分割結(jié)果。其中,第一列為帶有初始輪廓(白色線)的合成圖像,后面依次是CV模型、Cai模型、CVB模型、ICTM-GAC算法、ICTM-CV算法和文中算法的分割結(jié)果。從圖中可以看出,輸入圖像中的偏移場嚴重影響了CV模型和Cai模型的分割性能,這兩種模型都未能對目標進行整體分割。CVB模型在分割灰度異質(zhì)圖像時,能對輸入圖像進行偏移場矯正,但不能分割出一些尖角邊緣。ICTM-GAC算法分割效果不佳,ICTM-CV算法受偏移場影響分割出錯誤邊緣,然而文中算法能精確地分割出尖角邊緣,得到最優(yōu)的分割結(jié)果。進一步,觀察表2中各個模型分割結(jié)果的DSC、HD值,發(fā)現(xiàn)文中算法有較高的DSC值與較低的HD值,分割結(jié)果較好。

圖1 不同算法對合成圖像的分割結(jié)果

表2 圖1中分割結(jié)果的DSC、HD值、迭代收斂步數(shù)及收斂時間

其次,對加有不同偏移場的灰度異質(zhì)MR圖像中的腦白質(zhì)進行分割,實驗結(jié)果如圖2所示?？梢园l(fā)現(xiàn),文中算法能準確識別細微的腦灰質(zhì)得到更精確的分割結(jié)果。進一步,通過表3的數(shù)據(jù)進行定量分析,文中算法取得的DSC和HD值要優(yōu)于其他算法,且其收斂速度更快。

圖2 不同算法對MR腦圖像白質(zhì)的分割結(jié)果

表3 圖2中分割結(jié)果的DSC、HD值、迭代收斂步數(shù)及收斂時間

3.2 原始圖像數(shù)值實驗結(jié)果

為驗證文中方法對原始圖像分割的有效性,圖3給出了文中方法對來自魏茨曼分割數(shù)據(jù)集[22]上的8幅灰度異質(zhì)原始自然圖像的分割結(jié)果,并與其他方法作對比。文中方法能精準分割出目標物體邊緣且在一些細節(jié)邊緣表現(xiàn)較好。表4給出了各算法對6幅原始圖像分割結(jié)果的DSC值,文中算法的DSC均值最大,且方差較小,穩(wěn)定性較強?？梢哉f明,文中算法對原始圖像分割精度較高。

圖3 不同算法對原始圖像的分割結(jié)果

表4 圖3中分割結(jié)果的DSC值 %

4 結(jié)束語

在CVB模型的基礎(chǔ)上,提出基于測地輪廓長度和特征函數(shù)的灰度異質(zhì)圖像分割方法。通過對灰度異質(zhì)的合成圖像、MR圖像和原始圖像進行分割,結(jié)果表明該方法分割精度較高,且收斂速度明顯提升。