徐 娟,潘振寬,魏偉波,王加忠

(青島大學(xué) 計(jì)算機(jī)科學(xué)技術(shù)學(xué)院,山東 青島 266071)

0 概述

多相圖像分割是圖像處理、圖像分析、計(jì)算機(jī)視覺(jué)等領(lǐng)域中的研究熱點(diǎn)之一,而三維多相圖像分割則是其中的重點(diǎn)和難點(diǎn)。三維多相圖像分割在模式識(shí)別、醫(yī)學(xué)診斷、地球物理勘探、三維圖像重建等方面具有重要的應(yīng)用價(jià)值[1-2]。

在圖像分割領(lǐng)域中,人們提出許多關(guān)于分割問(wèn)題的模型。CHAN和VESE將分段常值Mumford-Shah模型[3-4]簡(jiǎn)化,結(jié)合水平集方法[5-6]建立了基于區(qū)域的分段常值兩相圖像分割模型,即Chan-Vese模型[7-8]。隨后,VESE和CHAN又提出了用n個(gè)水平集函數(shù)劃分2n個(gè)區(qū)域的區(qū)域競(jìng)爭(zhēng)策略,進(jìn)一步將該模型擴(kuò)展為分段常值和分段光滑的多相圖像分割模型,即多相圖像分割Chan-Vese模型[9-10]。隨著水平集個(gè)數(shù)增加,目標(biāo)區(qū)域的表達(dá)及其對(duì)應(yīng)的能量泛函求解方程越來(lái)越復(fù)雜,為此,文獻(xiàn)[11]根據(jù)劃分區(qū)域的編號(hào)與其二進(jìn)制表達(dá)的關(guān)系建立每個(gè)區(qū)域的特征函數(shù)的通用表達(dá)式。文獻(xiàn)[12]采用文獻(xiàn)[13]提出的區(qū)域劃分方案,用n個(gè)連續(xù)水平集函數(shù)的零水平集建立了n相圖像分割的變分水平集方法。文獻(xiàn)[14]用n個(gè)離散常值的水平集函數(shù),結(jié)合拉格朗日多項(xiàng)式差值,建立了分段常值多相圖像分割的變分模型。文獻(xiàn)[15]提出一種新的區(qū)域劃分方法,該方法采用幾何的方式引入n個(gè)水平集來(lái)劃分n+1個(gè)區(qū)域。但上述模型都存在以下不足:均需對(duì)多個(gè)函數(shù)求極值,計(jì)算效率不高;隨著水平集個(gè)數(shù)增加,目標(biāo)區(qū)域的表達(dá)及對(duì)應(yīng)的能量泛函求解方程越來(lái)越復(fù)雜。此外,多相水平集模型難以確定用于多相分割的水平集函數(shù)的個(gè)數(shù)。針對(duì)上述問(wèn)題,文獻(xiàn)[16]基于Chan-Vese多相圖像分割模型的區(qū)域表達(dá)策略,提出了用一個(gè)連續(xù)變化的水平集函數(shù)建立表達(dá)n相變分形式的分層水平集Chung-Vese模型,采用最少的水平集函數(shù)來(lái)表達(dá)多個(gè)相,無(wú)需增加約束條件,避免了重疊和漏分現(xiàn)象。但該模型只提供了一個(gè)區(qū)域分割方案,沒(méi)有給出通用區(qū)域分割的公式表達(dá),當(dāng)相數(shù)較多時(shí)將導(dǎo)致能量泛函和相關(guān)水平集函數(shù)演化方程表達(dá)復(fù)雜,為此,文獻(xiàn)[17]對(duì)Chung-Vese模型提出區(qū)域標(biāo)記的統(tǒng)一化表達(dá)。多層水平集Chung-Vese變分模型在二維多相圖像分割中得到充分的研究和應(yīng)用[18],但是如何將其應(yīng)用到三維多相圖像分割中,還未見相關(guān)研究報(bào)導(dǎo)。

文獻(xiàn)[13]提出的變分水平集方法將變分法和水平集方法相結(jié)合,由于其具有自適應(yīng)復(fù)雜拓?fù)浣Y(jié)構(gòu)變化、二維和三維圖像分割模型表達(dá)一致、數(shù)值計(jì)算方法穩(wěn)定以及具有集成多模型信息能力等特點(diǎn),已被廣泛用于圖像分割的研究領(lǐng)域中[19]。但變分水平集方法存在計(jì)算效率低等問(wèn)題,為此,文獻(xiàn)[20]將分裂算法與Bregman迭代[21]相結(jié)合,通過(guò)引入輔助變量和迭代乘子,提出了快速Split-Bregman方法。該方法具有較高的計(jì)算效率和迭代精度。

文獻(xiàn)[13]設(shè)計(jì)了一個(gè)函數(shù)在多層水平集標(biāo)記的模型實(shí)現(xiàn)對(duì)圖像不同區(qū)域的分割,但只驗(yàn)證了該模型對(duì)二維圖像的分割效果,并未指出其對(duì)三維圖像是否具有同樣的通用性和實(shí)用性。由于一個(gè)Lipschitz連續(xù)多層水平集函數(shù)能有效地對(duì)兩相和多相三維圖像分割,本文將文獻(xiàn)[13]的二維圖像分割模型向三維圖像進(jìn)行擴(kuò)展,設(shè)計(jì)滿足圖像不同區(qū)域互不重疊、劃分方案對(duì)稱、符合唯一劃分條件的特征函數(shù),在變分水平集理論框架下基于分段常值假設(shè)建立統(tǒng)一的三維圖像的兩相和多相變分模型,并采用Split-Bregman投影方法[22]對(duì)模型求解,即在Split-Bregman方法基礎(chǔ)上引入投影變量,將約束轉(zhuǎn)化為解析投影公式進(jìn)行求解。

1 多層水平集函數(shù)的變分方法

多相圖像分割是按照一定的區(qū)域劃分方案將圖像分割成多個(gè)相的過(guò)程,一個(gè)有效的多相圖像分割模型需要解決的問(wèn)題包括能量泛函的表達(dá)、多相區(qū)域的劃分策略以及演化方程的數(shù)值計(jì)算等,其中,區(qū)域劃分方案的選擇和區(qū)域特征函數(shù)的設(shè)計(jì)是多相圖像分割的核心。為防止區(qū)域之間的重疊和漏分,多相圖像分割必須解決好區(qū)域之間的競(jìng)爭(zhēng)策略。下文將介紹本文區(qū)域劃分方案對(duì)應(yīng)的區(qū)域特征函數(shù)表達(dá)和多層水平集函數(shù)的變分方法。

1.1 變分水平集方法

水平集方法[7]是一種有效的演化曲線/曲面隱式表示,即將分割邊界/表面隱式表示為在歐拉框架中演化的更高維函數(shù)的水平集,該高維函數(shù)即水平集函數(shù)。水平集方法自動(dòng)處理拓?fù)渥兓?為分割準(zhǔn)則的設(shè)計(jì)提供了很大的靈活性,以適應(yīng)對(duì)圖像及其結(jié)構(gòu)的各種假設(shè)[23],包括不同的外觀模型[23-25]和形狀先驗(yàn)[26]。水平集方法對(duì)二維和三維甚至更高維圖像分割具有相同的表達(dá)形式。因此,在三維圖像分割研究中,二維圖像分割的能量模型以及區(qū)域劃分方法依然適用。

變分水平集方法是文獻(xiàn)[13]針對(duì)閉合曲線C演化的能量泛函E(C)最小化問(wèn)題提出的一種新的水平集方法,通過(guò)引入嵌入函數(shù)φ(x)和Heaviside函數(shù),將E(C)改造成E(φ(x)),再利用變分法求得關(guān)于φ(x)的偏微分方程[27],其中要求保證水平集函數(shù)在演化過(guò)程中始終保持符號(hào)距離函數(shù)的特征|φ|=1。C={(x)|φ(x)=c0},C是滿足函數(shù)φ(x)等于常值c0的點(diǎn)集,稱為函數(shù)φ(x)的一個(gè)水平集。函數(shù)φ(x)為曲線C的嵌入函數(shù),稱為水平集函數(shù)。當(dāng)常數(shù)c0=0時(shí),C={(x)|φ(x)=0}被稱為水平集函數(shù)φ(x)的零水平集,如圖1所示。本文基于多層水平集函數(shù)的n層水平集在圖像中的曲面演化對(duì)圖像進(jìn)行分割,需要通過(guò)變分水平集方法解決曲面演化的能量泛函的極值問(wèn)題。

圖1 水平集函數(shù)和零水平集對(duì)應(yīng)的輪廓線

1.2 多層水平集方法

多層水平集模型[14]在Chan-Vese模型基礎(chǔ)上引入外延生長(zhǎng)島動(dòng)力學(xué),利用水平集函數(shù)的多層等高線的曲線演化來(lái)完成圖像分割,得到改進(jìn)的水平集圖像分割算法模型,以處理多相圖像分割問(wèn)題。該模型也被稱為Chung-Vese模型,比多相水平集模型簡(jiǎn)單,運(yùn)算速度快。多層水平集函數(shù)的2層嵌套水平集示意圖如圖2所示。

圖2 多層水平集函數(shù)的2層嵌套水平集示意圖

1.3 多層水平集模型的區(qū)域劃分方案

Ω1:χ1(φ)=H(φ-l0)(1-H(φ-l1)),l0≤φ(x)≤l1

Ω2:χ2(φ)=H(φ-l1)(1-H(φ-l2)),l1<φ(x)≤l2

?

Ωn:χn(φ)=H(φ-ln-1)(1-H(φ-ln)),ln-1<φ(x)≤ln

(1)

由φ(x)∈[l0,ln](0=l0

x∈Ωi:li-1<φ(x)≤li

φ-l0>0,φ-li≤0

φ-l1>0,φ-li+1≤0

? ?

φ-li-1>0,φ-ln≤0

(2)

當(dāng)χi(x)=1時(shí),有:

χ1(x)=0,χ2(x)=0,…,χi-1(x)=0

χi+1(x)=0,χi+2(x)=0,…,χn(x)=0

(3)

即:

(4)

2 多層水平集的Split-Bregman投影方法

設(shè)f為三維圖像的圖像強(qiáng)度,多相圖像分割變分水平集模型通過(guò)引入1個(gè)包含n個(gè)不同水平集(0=l0

(5)

滿足約束條件:

(6)

以保證水平集函數(shù)在演化過(guò)程中始終保持符號(hào)距離函數(shù)的特征。ui=(u1,u2,…,un)為f在區(qū)域Ω中的分段常值,其估計(jì)式為:

(7)

本文采用規(guī)整化的Heaviside函數(shù)和Dirac函數(shù)[28]實(shí)現(xiàn)區(qū)域劃分的函數(shù)表達(dá),即當(dāng)ε→0時(shí),Hε(φ)→H(φ)(ε趨于0時(shí),Hε(φ)近似Heaviside函數(shù)H(φ)),其表達(dá)式為:

(8)

(9)

其中,ε為數(shù)值較小的正數(shù),為書寫簡(jiǎn)潔,本文仍使用原符號(hào)表達(dá)。所以,邊緣項(xiàng)可以等價(jià)表示為:

(10)

因此,多層水平集圖像分割模型的能量泛函可以等價(jià)表示為:

(11)

本文采用高效、穩(wěn)定的Split-Bregman投影方法對(duì)能量泛函求解極值。引入輔助變量w和Bregman迭代參數(shù)b,將式(11)轉(zhuǎn)變?yōu)槿缦伦兎帜Ｐ偷袷?

(12)

s.t.|w|=1

(13)

其中,w0=0,b0=0,θ(θ>0)為懲罰參數(shù)。采用交替優(yōu)化方法分別得到關(guān)于φ的歐拉拉格朗日方程及關(guān)于w的廣義軟閾值公式如式(14)和式(15)所示。

(14)

(15)

其中:

(16)

(17)

為滿足約束φ(x)∈[l0,ln],對(duì)函數(shù)φ(x)添加一個(gè)約束項(xiàng)為:

φk+1=min(max(l0,φk+1),ln)

(18)

(19)

更新b得:

bk+1=bk+φk+1-wk+1

(20)

本文算法步驟如下:

步驟1初始化φ0為水平集函數(shù),w0=b0=0,k=0。

步驟2估計(jì)ui,i=1,2,…,n。

步驟3依次求解式(14)、式(15)、式(19)得到φk+1,wk+1。

步驟4更新bk+1=bk+φk+1-wk+1。

3 實(shí)驗(yàn)與結(jié)果分析

本文實(shí)驗(yàn)PC機(jī)的配置為:操作系統(tǒng)Win7 x64,處理器Intel(R)Core(TM)i5-4590 CPU @ 3.30 GHz,RAM 4.00 GB,編程運(yùn)行環(huán)境MATLAB R2018b。三維圖像共有3個(gè)維度,常見的三維圖像有真實(shí)人體CT掃描圖像和人工合成的三維圖像。本文實(shí)驗(yàn)參數(shù)設(shè)置為空間步長(zhǎng)h=1,時(shí)間步長(zhǎng)Δt=0.1。

3.1 三維圖像2相分割與重建

在圖像分割模型中,2相圖像分割是多相圖像分割的一種特殊情況。本實(shí)驗(yàn)取真實(shí)人體牙齒周邊部位CT掃描圖像序列,圖像規(guī)格為150×128×105。在實(shí)驗(yàn)中,將原本彩色圖像預(yù)處理為灰度圖像,將原本的灰色背景處理為白色以實(shí)現(xiàn)2相圖像分割,其中,牙齒和下頜骨整體作為分割目標(biāo)。預(yù)處理后的牙齒圖像序列的3個(gè)斷層圖像如圖3所示。該實(shí)驗(yàn)采用一個(gè)水平集函數(shù)的2個(gè)水平集l1=1和l2=3,相關(guān)參數(shù)ε=8,α=1,γ=1,θ=1.5,分割過(guò)程如圖4所示。

圖3 2相牙齒圖像序列的斷層圖像

圖4 2相圖像序列分割過(guò)程

3.2 三維圖像3相分割與重建

本實(shí)驗(yàn)使用規(guī)格為158×128×99的真實(shí)人體下頜部位CT掃描圖像序列和128×128×105的牙齒部位CT掃描圖像序列。該實(shí)驗(yàn)采用一個(gè)水平集函數(shù)的3個(gè)水平集l1=1、l2=3和l3=7,相關(guān)參數(shù)ε=8,α=1,γ=0.09×2552,θ=0.3。下頜部位和牙齒部位的圖像序列中3個(gè)斷層圖像分別如圖5和圖6所示,其分割過(guò)程如圖7和圖8所示。

圖5 下頜圖像序列的斷層圖像

圖6 牙齒圖像序列的斷層圖像

圖7 3相下頜圖像序列分割過(guò)程

圖8 3相牙齒圖像序列分割過(guò)程

3.3 三維圖像4相分割與重建

本實(shí)驗(yàn)使用人工合成的三維立體幾何圖像序列,圖像規(guī)格為l1=1,l2=3,l3=7。該實(shí)驗(yàn)采用一個(gè)水平集函數(shù)的4個(gè)水平集l1=1,l2=3,l3=7和l4=15,相關(guān)參數(shù)ε=8,α=1,γ=0.01×2552,θ=0.08。幾何立體圖像序列中3個(gè)斷層圖像如圖9所示,其分割過(guò)程如圖10所示。

圖9 立體幾何圖像序列的斷層圖像

圖10 4相圖像序列分割過(guò)程

3.4 噪聲對(duì)本文模型分割效果的影響

不同類型的圖像符合不同的概率分布,限于篇幅,本文僅給出Gauss噪聲對(duì)4相立體幾何圖像序列分割效果的影響。圖11為立體幾何圖像序列隨機(jī)加入噪聲值分別為0.0、0.3、0.5、0.7、0.9和1.0后相對(duì)應(yīng)的分割結(jié)果,從圖中可以看出,當(dāng)噪聲值小于1時(shí),噪聲對(duì)區(qū)域劃分影響不大,但是會(huì)影響圖像的分割效果,導(dǎo)致分割后圖像表面不光滑。

圖11 隨機(jī)加入不同噪聲值后對(duì)應(yīng)的分割結(jié)果

3.5 與多相圖像分割Chan-Vese模型效果的對(duì)比

多相圖像分割Chan-Vese模型(以下簡(jiǎn)稱多相Chan-Vese模型)劃分n個(gè)區(qū)域需要引入lbn個(gè)水平集函數(shù)。除本文模型以外,該模型是在表達(dá)多個(gè)相時(shí)使用的水平集函數(shù)個(gè)數(shù)最少的多相圖像分割模型,因此,本文將其選為對(duì)比模型。圖12～圖15依次為2相、3相(下頜)、3相(牙齒)及4相圖像序列分別采用本文模型和多相Chan-Vese模型最終三維分割效果的主觀對(duì)比。表1為2相、3相(下頜)、3相(牙齒)及4相圖像序列分別采用本文模型和多相Chan-Vese模型最終三維分割效果的客觀對(duì)比。

圖12 本文模型和多相Chan-Vese 模型的2相圖像分割效果對(duì)比

圖13 本文模型和多相Chan-Vese 模型的3相牙齒圖像分割效果對(duì)比

Fig.13 Effect comparison of three-phase tooth image segmentation between the proposed model and the multiphase Chan-Vese model

圖14 本文模型和多相Chan-Vese模型的3相下頜圖像分割效果對(duì)比

Fig.14 Effect comparison of three-phase mandibular image segmentation between the proposed model and the multiphase Chan-Vese model

圖15 本文模型和多相Chan-Vese模型的4相圖像分割效果對(duì)比

表1 本文模型與多相Chan-Vese模型的實(shí)驗(yàn)數(shù)據(jù)比較

從以上實(shí)驗(yàn)結(jié)果可以看出,本文模型可以有效地實(shí)現(xiàn)三維多相圖像分割,并且迭代步數(shù)較少,分割速度較快。

4 結(jié)束語(yǔ)

本文針對(duì)三維多相圖像分割提出一種新的變分水平集模型。利用基于一個(gè)連續(xù)水平集函數(shù)的多層水平集對(duì)圖像進(jìn)行區(qū)域劃分,采用Split-Bregman方法實(shí)現(xiàn)三維多相圖像分割。本文使用真實(shí)人體牙齒和下頜部位的CT掃描圖像序列進(jìn)行實(shí)驗(yàn),結(jié)果表明,該模型能有效地分割圖像。本文模型具有一定的實(shí)際應(yīng)用價(jià)值,例如醫(yī)生可以通過(guò)分割后重建的三維圖像看見病人牙齒的缺損情況。但本文實(shí)驗(yàn)未對(duì)包含噪聲的原始數(shù)據(jù)進(jìn)行預(yù)處理,導(dǎo)致分割得到的圖像不光滑,后續(xù)將對(duì)此改進(jìn),提高分割圖像的光滑度。