基于非局部總變差的圖像分割活動輪廓模型

摘 要:在一般活動輪廓模型的連續全局極小化方法基礎上，利用四種非局部總變差，給出了一種具有連續全局極小解的非局部活動輪廓模型。由于該模型的非局部特性，在分割過程中能有效地去除圖像中的噪聲，同時保留那些重復的精細結構。數值實驗證明，該模型能將圖像中的主體結構和精細結構很好地分割出來，而標準活動輪廓模型的分割結果中則丟掉了許多小的精細結構。

關鍵詞:圖像分割; 活動輪廓模型; 水平集方法; 連續全局極小化方法; 非局部總變差

中圖分類號:TP 391文獻標志碼:A

文章編號:1001-3695(2010)06-2373-04

doi:10.3969/j.issn.10013695.2010.06.108

Image segmentation active contour model based on nonlocal total variation

ZHANG Wenjuan1，2， FENG Xiangchu2

(1.Dept. of Mathematics Physics， Xi’an Technological University， Xi’an 710032， China; 2.School of Science， Xidian University， Xi’an 710071， China)

Abstract:Based on the continuous global minimizing method for the general active contour model， this paper proposed a nonlocal active contour model having continuous global minimizing solutions by using four nonlocal total variation. Due to the nonlocal property of the model proposed here， noise could be effectively removed in the course of segmentation， and the repeated fine structures could be preserved simultaneously. Numerical experiment shows this model can segment the main structures and the fine structures very well， however many small fine structures are lost in the results of the standard active contour model.

Key words:image segmentation; active contour model; level set method; continuous global minimizing method; nonlocal total variation

0 引言

近年來空域非局部信息在圖像去噪中得到了廣泛應用，非局部濾波及非局部變分模型比以往很多經典的變分模型和PDE(偏微分方程)方法產生了更好的去噪效果。在非局部方法中，用圖表示圖像，圖像的像素點定義為圖的節點，對具有相似圖像特征(強度、顏色、紋理等)的像素點，即使它們的空間距離很大，也認為是相近的，在去噪過程中，這些像素點會被保留下來。由于以往經典的變分模型和PDE方法是建立在圖像局部信息基礎上的，對分片光滑圖像的去噪問題來說，這已經足夠了，而當圖像中含有類似于噪聲的紋理或其他精細結構時，利用這些方法會導致這些小結構丟失或部分丟失。基于圖像的非局部信息的去噪方法中，影響較大的是Buades等人[1]提出的非局部平均算法，它是一種圖像濾波模型，其濾波器是基于圖像各像素點之間的鄰域塊距離計算出來的，Gilboa等人[2]給出了非局部平均算法的變分形式，以及相應的非局部偏微分方程。雖然關于非局部去噪模型方面的研究還有很多，但基于非局部信息的圖像分割模型還很少。

圖像分割是針對圖像的不同特性把一幅給定圖像分割為幾個子區域的過程，這些特性包括圖像的強度、顏色以及紋理等。由于分割所依據的圖像特性的不同，形成了不同的圖像分割方法，這里所研究的內容屬于變分圖像分割模型中最普遍的活動輪廓模型。變分模型的主要思想是求與原始圖像相關的某個特定的目標泛函的極小解，其中活動輪廓模型是在算法過程中使得圖像的邊界曲線活動，達到穩定態后將圖像分為幾個子區域，從而達到分割圖像的目的，活動輪廓模型中較普遍的是ChanVese活動輪廓模型[3]、測地線活動輪廓模型[4]和MumfordShah模型[5]。求解活動輪廓模型最常用的是由Osher等人[6]提出的水平集方法，它用一個高維函數的零水平集來表示圖像的邊界曲線，將邊界曲線的活動轉換為水平集函數的變化過程。利用水平集函數不僅可以用其零水平集表示邊界曲線，而且可以表示邊界曲線的內部和外部，此外，還可以處理圖像區域的拓撲變化，是近年來圖像分割及其他領域一種非常重要的方法。但是，水平集方法極小化過程很慢，且需要不斷重新定義水平集函數，而且利用水平集方法求解活動輪廓模型所得為局部極小解，即要取得好的圖像分割效果，初始輪廓的位置必須靠近精確解，盡管有很多算法可以提高活動輪廓的速度，但遺憾的是這些算法得到的均為局部極小解。文獻[7，8]針對前面所述的三種活動輪廓模型提出了一種全局連續極小化方法，利用此方法求得的目標函數的極小解不依賴于初始輪廓的位置，事實上，這種全局極小化方法可以推廣到任意活動輪廓模型[9]。本文目的是在文獻[7，8]的活動輪廓模型的全局連續極小化方法基礎上，將模型中的總變差正則項推廣到非局部情形，給出了四種非局部總變差的形式[10～12]。實驗證明，利用非局部總變差比標準總變差能更好地分割精細結構，這對圖像分割(尤其是醫學圖像分割)研究是很有意義的。

1 非局部總變差

首先給出一些基本的非局部算子的定義，這里用圖表示圖像，灰度圖像u(x):ΩR2→R的每個像素點為圖的節點，u(x)在點x∈Ω的非局部梯度向量NLu(x)的定義如下:

(NLu)(x，y)=(u(y)-u(x))ω(x，y):Ω×Ω→R(1)

其中:0<ω(x，y)<∞為點x和y之間的權系數，計算權系數的方法有很多種，參見文獻[2]，這里的向量概念不同于以往，是關于x和y的函數。兩個非局部向量p1，p2:Ω×Ω→R在點x∈Ω的點積(p2#8226;p2)(x)定義為

(p1#8226;p2)(x)=∫Ω p1(x，y)p2(x，y)dy:Ω→R(2)

由此，非局部向量p:Ω×Ω→R在點x∈Ω的模的大小為

|p|(x)=∫Ω p(x，y)2dy:Ω→R(3)

非局部散度定義為非局部梯度的共軛，因此對向量p:Ω×Ω→R，可推出在點x∈Ω有

(divNLp)(x)=∫Ω(p(x，y)-p(y，x))ω(x，y)dy:Ω→R(4)

從而，非局部拉普拉斯為

ΔNLu(x)=12divNL(NLu(x))=∫Ω(u(y)-u(x))(5)

ω(x，y)dy:Ω→R

下面給出四種非局部正則化泛函，首先，Gilboa等人[10]給出了兩種正則化非局部泛函:第一種是基于上面的非局部梯度;第二種是基于差分，設為一正函數，滿足(0)=0，基于非局部梯度的泛函為

J1(u)=∫Ω(|NLu|2)dx=∫Ω(|ωu|2)dx=∫Ω∫Ω(u(y)-u(x))2ω(x，y)dydx(6)

當取(s)=s時得到第一種非局部總變差，且此時泛函為凸:

NLTV1(u)=∫Ω∫Ω(u(y)-u(x))2ω(x，y)dy dx(7)

泛函式(6)對u的導數(EulerLagrange)為

uJ1(u)=-2∫Ω(u(y)-u(x))ω(x，y)(′(|NLu|2(x))+′(|NLu|2(y)))dy(8)

基于差分的非局部泛函為

J2(u)=∫Ω×Ω(u(y)-u(x))2ω(x，y)dydx(9)

同樣，取(s)=s時得到另一種非局部總變差:

NLTV2(u)=∫Ω×Ω|u(x)-u(y)|ω(x，y)dydx(10)

泛函式(9)對u的導數為

uJ2(u)=-4∫Ω(u(y)-u(x))ω(x，y)′((u(y)-u(x))2ω(x，y))dy(11)

注意到，當(s)=s時，泛函式(6)和(9)相等，非局部總變差式(7)和(10)分別對應于局部二維情形下的各向同性總變差和各向異性總變差，可以通過常用的最速下降法對泛函式(6)(9)極小化。對泛函式(9)，還可以利用基于圖—切割技術的快速算法，對于以式(7)為正則項的非局部ROF、非局部TVL1以及非局部TVG等模型還可以通過投影算法求解。

Gilboa等人[11]給出了第三種正則化非局部泛函:

J3(u)=12∫Ω×Ω(|u(x)-u(y)|)ω(x，y)dxdy(12)

其中:為凸的正函數，同樣要滿足(0)=0。對L1型的泛函，進一步假設lims→∞ (s)/s=1，當(s)=s時，得到第三種非局部總變差:

NLTV3(u)=12∫Ω×Ω|u(x)-u(y)|ω(x，y)dxdy(13)

泛函式(12)對u的EulerLagrange為

uJ3(u)=∫Ω′(|u(x)-u(y)|)u(x)-u(y)|u(x)-u(y)|ω(x，y)dxdy(14)

此外，Kindermann等人[12]給出了第四種非局部正則化泛函:

J4(u)=∫Ω×Ω(|u(x)-u(y)|2)ω(x，y)dxdy(15)

其中:為非負可微函數，則第四種非局部總變差為

NLTV4(u)=∫Ω×Ω|u(x)-u(y)|ω(x，y)dxdy(16)

泛函式(15)對u的EulerLagrange為

uJ4(u)=2(∫Ω′(|u(x)-u(y)|2)(u(x)-u(y))ω(x，y)dy)(17)

可以利用最速下降法和圖—切割技術求泛函式(12)和(15)的極小解。上述四種非局部正則化泛函之間的聯系以及在何種情況下用哪種泛函效果更好，目前尚在研究當中，就本文給出的圖像分割實驗來說，利用上面四種泛函的分割效果相差無幾。

2 基于非局部總變差的活動輪廓模型

圖像分割中，一種成功的變分模型是由Kass等人最先提出的活動輪廓模型，該模型可以將圖像的各種信息，如邊界、區域以及形狀等結合起來，二相位活動輪廓模型的一般形式為

EAC(C)=∫Cgbds+λ∫Cinginrdx+λ∫CoutgoutrdxC(s):[0，1]→R2(18)

其中:C為參變量s∈[0，1]的閉曲線(輪廓);gb:Ω→R+為邊緣函數;Cin，outΩ分別表示曲線C的內部和外部;gin，outr為曲線C的內部和外部區域檢測函數;λ為一正常數。取ginr=goutr=0，邊緣檢測函數gb=gb(|u0(C(s))|)，且gb(t):R+→R+單調下降，則模型式(18)變為測地線活動輪廓模型EGAC(C)=∫Cgb(|u0(C(s))|)ds，這個模型是良態的且可以通過水平集方法求解，所以得到了廣泛應用;若取gb=1，ginr=(μin-I)2 ，goutr=(μout-I)2，其中，I為要分割的圖像，則得到另一種廣泛應用的活動輪廓模型，即ChanVese活動輪廓模型:

ECV(C)=∫Cds+λ∫Cin(μin-I)2dx+λ∫Cout(μout-I)2dx

它是MumfordShah模型的二相位分片常數逼近。MumfordShah模型的二相位分片光滑逼近為

EVC(C)=∫Cds+λ∫Cin|sin|2+(sin-I)2dx+λ∫Cout|sout|2+(sout-I)2dx

事實上，就是在泛函式(18)中取gb=1，ginr=|sin|2+(sin-I)2，goutr=|sout|2+(sout-I)2。另外，許多其他的活動輪廓模型也都可以表示為式(18)的形式。

求解活動輪廓模型常用的是水平集方法，式(18)的水平集形式為

ELSM()=∫Ωgb|H()|dx+λ∫ΩginrH()dx+λ∫Ωgoutr(1-H())dx，:Ω→R(19)

其中:即為水平集函數。輪廓曲線C可以用的零水平集表示，即C={x|x∈Ω，(x)=0}，利用水平集函數還可以表示C的內部和外部:Cin={x|x∈Ω，(x)>0}，Cout={x|x∈Ω，(x)<0}，H為Heaviside函數，當≥0時，H()=1，<0時，H()=0。可以看出，式(19)與(18)是等價的，式(19)中的第一項∫Ωgb|H()|dx=∫Ωgb||δ()dx=∫Cgbds，第二項等價于λ∫Cinginrdx，同理，第三項等價于λ∫Coutgoutrdx。對能量泛函式(19)極小化可以通過標準的EulerLagrange方程求解，水平集方法是求解活動輪廓模型式(18)的一種非常重要的方法，它將圖像邊界曲線的活動轉換為水平集函數的迭代過程，可以處理拓撲的自然變化，而且利用水平集方法可以在分割過程中很方便地引入圖像的各種信息，如邊界、區域和形狀等;另外，水平集方法的數值實現是基于雙曲守恒定律及迎風格式，可以保證算法的穩定性和準確率。但是，這里要假設函數H、δ的正則性，且必須不斷重新定義水平集函數為帶符號的距離函數(這意味著解Eikonal方程||=1)以保證表示輪廓的零水平集準確而光滑地活動，而且利用水平集方法所求得的為局部極小解，也就是說，解的好壞依賴于初始輪廓的位置，初始輪廓要盡量靠近精確解，才能對圖像準確地實施分割。

Chan等人[7]給出了ChanVese模型全局連續極小解的存在性定理以及求解方法;文獻[8]中Bresson等人將ChanVese模型的全局連續極小化方法推廣到測地線活動輪廓模型和二相位分片光滑MumfordShah模型，也給出了計算這兩種模型全局連續極小解的數值方法。利用文獻[7，8]中的方法不僅使得模型的解不依賴于初始輪廓的位置，在運算速度上也比水平集方法有了很大提高，而且還避免了水平集方法需要重新定義水平集函數和函數H、δ的正則化等諸多問題。上面所述的全局連續極小化方法可以應用于式(18)所示的一般的活動輪廓模型，其變分形式如下

EGMAC(u)=∫Ωgb|u|dx+λ∫Ωginrudx+λ∫Ωgoutr(1-u)dxu:Ω→[0，1](20)

下面的定理說明了式(18)定義的變分問題全局連續極小解的存在性及求解方法。

定理1 設gb:Ω→R+，對任意給定的ginr，goutr:Ω→R及λ∈R+，若u*為EGMAC的任一極小解，則對任意的v∈R，集合ΩC(v)={x∈Ω:u(x)>v}的特征函數(這里C為集合ΩC的邊界)為式(20)和(18)定義的EGMAC和EAC的全局極小解。

本文將式(20)中的總變差正則項推廣到非局部情形，因為若利用標準總變差作為正則項，則意味著水平集的曲率越大，正則化程度就越高。當對圖像(尤其是醫學圖像)中有用的精細結構進行分割時，由于這些結構的曲率較大，故會被光滑掉或部分光滑掉，而若利用非局部總變差作為正則項，其非局部特性會使得分割過程中去除圖像噪聲的同時很好地保留圖像中重復的小結構。這里，要求如下能量泛函的極小解:

ENLGMAC(u)=gbNLTV(u)+λ∫Ωginrudx+λ∫Ωgoutr(1-u)dxu:Ω→[0，1](21)

其中:NLTV(u)可以取前面給出的四種形式。下面求解該變分問題，首先引入一個新的函數v:Ω→R，使得

ENLGMAC(u，v)=gbNLTV(u)+λ∫Ω(ginr-goutr)vdx+12θ∫Ω(u-v)2dx(22)

其中:λ，θ>0是兩個常數。顯然，θ要足夠小，因為泛函(22)關于u，v均為凸，所以式(22)分別關于u，v求極小，然后迭代直到收斂，即可求得式(21)的極小解:

a)固定u，泛函式(22)關于v求極小，即求

infv∈[0，1]λ∫Ω(ginr-goutr)vdx+12θ∫Ω(u-v)2dx(23)

容易求得其解為v=min(max(u-λθ(ginr-goutr)，0)，1)。

b) 固定v，泛函式(22)關于u求極小，即求

infugbNLTV(u)+12θ∫Ω(u-v)2dx(24)

此為非局部ROF，利用最速下降法求解，令p(u)=uNLTV(u)，式(8)(11)(14)(17)分別給出了四種非局部總變差相應的p(u)的計算方法，只要取合適的函數即可。式(24)相應的EulerLagrange方程為

gbp(u)+1θ(u-v)=0(25)

取初始值u|t=0=v，迭代ut=-gbp(u)+1θ(v-u)即可求得非局部ROF式(24)的極小解。

可以證明，在圖像塊流形上非局部總變差收斂于標準的總變差，所以只要在流形上利用Coarea公式，就可以把定理1推廣到由強度塊流形定義的圖像表示上來[9]。如果u*為ENLGMAC的任意極小解，在流形上對u*閾值意味著在Ω上對u*閾值， 則對任意的v∈[0，1]，集合ΩC(v)={x∈Ω:u*(x)>v}的特征函數為ENLGMAC的一個全局極小解，同時也是活動輪廓模型式(18)的全局極小解。

3 數值實驗

下面以ChanVese模型為例檢驗本文方法的圖像分割效果，在式(21)中取gb=1，ginr=(μin-I)2，goutr=(μout-I)2。其中:I為要分割的圖像，泛函(21)分別關于μin、μout求極小可得μin=∫ΩI#8226;u(x)dx∫Ωu(x)dx，μout=∫ΩI#8226;(1-u(x))dx∫Ω(1-u(x))dx，即為原始圖像I在Cin，Cout內的強度均值。圖1給出了要分割的人腦醫學圖像，圖像大小為256×256，圖中的矩形框為初始輪廓的位置，同時也是圖2(d)和圖3(d)所對應的部分。圖2為利用標準總變差為正則項的具有全局連續極小解的ChanVese模型的圖像分割結果。圖3為以非局部總變差為正則項的具有全局連續極小解的ChanVese模型的分割結果。圖2(a)及圖3(a)分別為式(20)和(21)的全局連續極小解，利用第2章中所述的方法求得式(20)和(21)的任一極小解u*，任取v∈[0，1]，則u(x)=1 u*(x)>v0 u*

4 結束語

本文利用非局部凸泛函在圖像處理(尤其是圖像去噪)中所起的突出作用，將非局部總變差作為具有全局連續極小解的圖像分割活動輪廓模型的正則項，給出了一種非局部圖像分割方法，利用這種方法可以根據實際情況任意選取初始輪廓的位置，均可以得到好的分割效果，而且克服了水平集方法關于Heaviside函數H和δ函數的正則性問題和需要不斷重新定義水平集函數的問題。此外，由非局部總變差的特性及數值試驗證明，這種非局部圖像分割模型可以將圖像中小的精細結構很好地分割出來，這是以往標準的活動輪廓模型實現不了的。關于本文給出的四種非局部總變差之間的聯系和區別，以及本文方法理論方面的研究還有待在以后的工作中繼續進行。

參考文獻:

[1]BUADES A， COLL B， MOREL J M. A review of image denoising algorithms， with a new one[J]. SIAM Multiscale Modeling and Simulation， 2005，4(2):490-530.

[2]GILBOA G， OSHER S. Nonlocal linear image regularization and supervised segmentation[J]. SIAM Multiscale Modeling and Simulation，2007，6(2):595-630.

[3]CHAN T F， VESE L A. Active contours without edges[J]. IEEE Trans on Image Processing，2001，10(2):266-277.

[4]CASELLES V， KIMMEL R， SAPIRO G. Geodesic active contours[J]. International Journal of Computer Vision，1997，22(1):61-79.

[5]MUMFORF D， SHAH J. Optimal approximations of piecewise smooth functions and associated variation problems[J]. Communications on Pure and Applied Mathematics，1989，42:577-685.

[6]OSHER S， SETHIAN J A. Fronts propagating with curvaturedependent speed: algorithms based on HamiltonJacob formulations[J]. Journal of Computational Physics，1988，79(1):12-49.

[7]CHAN T F， ESEDOGLU S， NIKOLOVA M. Algorithms for finding global minimizers of image segmentation and denoising models[J]. SIAM Journal on Applied Mathematics，2006，66(5):1632-1648.

[8]BRESSON X， ESEDOGLU S， VANDERGHEYNST P， et al. Fast global minimization of the active contour/snake models[J]. Journal of Mathematical Imaging and Vision，2007，28(2):151-167.

[9]BRESSON X， CHAN T F. Nonlocal unsupervised variational image segmentation models[R]. [S.l.]:UCLA， 2008.

[10]GILBOA G， OSHER S. Nonlocal operators with applications to image processing[J]. SIAM Multiscale Modeling and Simulation，2008，7(3):1005-1028.

[11]GILBOA G， DARBON J， OSHER S， et al. Nonlocal convex functionals for image regularization[R]. [S.l.]:UCLA， 2006.

[12]KINDERMANN S， OSHER S， JONES P W. Deblurring and denoising of images by nonlocal functionals[J]. SIAM Multiscale Modeling and Simulation，2005，4(4):1091-1115.