基于測地輪廓和特征函數的灰度異質圖像分割

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

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

0 引 言

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

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

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

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

1 相關工作

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

1.1 基于特征函數的GAC模型

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

(1)

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

(2)

其中：τ>0,邊緣檢測函數為：

(3)

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

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

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

1.2 CVB模型

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

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

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

(4)

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

I=B+S

(5)

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

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

s.t.I=B+S

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

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

2 文中模型及算法

2.1 文中模型與算法

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

s.t.I=B+S

(6)

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

s.t.I=B+S

(7)

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

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

(8)

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

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

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

這里記：

(10)

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

(11)

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

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

其中：

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

對于?x∈Ω,可解出：

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

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

〈θk,I-B-Sk〉Ω

其對應的歐拉-拉格朗日方程為：

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

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

(3)求解關于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 穩定性分析

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

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

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

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

進一步,

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

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 實驗結果及分析

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

(12)

HD(A,B)：=

(13)

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

表1 參數λ1、λ2的取值

3.1 合成圖像與MR圖像的數值實驗結果

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

圖1 不同算法對合成圖像的分割結果

表2 圖1中分割結果的DSC、HD值、迭代收斂步數及收斂時間

其次,對加有不同偏移場的灰度異質MR圖像中的腦白質進行分割,實驗結果如圖2所示。可以發現,文中算法能準確識別細微的腦灰質得到更精確的分割結果。進一步,通過表3的數據進行定量分析,文中算法取得的DSC和HD值要優于其他算法,且其收斂速度更快。

圖2 不同算法對MR腦圖像白質的分割結果

表3 圖2中分割結果的DSC、HD值、迭代收斂步數及收斂時間

3.2 原始圖像數值實驗結果

為驗證文中方法對原始圖像分割的有效性,圖3給出了文中方法對來自魏茨曼分割數據集[22]上的8幅灰度異質原始自然圖像的分割結果,并與其他方法作對比。文中方法能精準分割出目標物體邊緣且在一些細節邊緣表現較好。表4給出了各算法對6幅原始圖像分割結果的DSC值,文中算法的DSC均值最大,且方差較小,穩定性較強。可以說明,文中算法對原始圖像分割精度較高。

圖3 不同算法對原始圖像的分割結果

表4 圖3中分割結果的DSC值 %

4 結束語

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