圖像去噪LOT模型的分裂Bregman方法＊

龐志峰,楊余飛?,林 玲

(1.湖南大學(xué)數(shù)學(xué)與計(jì)量經(jīng)濟(jì)學(xué)院,湖南長(zhǎng)沙 410082;2.廣東外語(yǔ)外貿(mào)大學(xué)信息科學(xué)技術(shù)學(xué)院,廣東廣州 510420)

由于圖像在生成、存儲(chǔ)或者傳遞過(guò)程中經(jīng)常被噪聲污染,因此圖像去噪是圖像處理中的一個(gè)最基本的問(wèn)題.隨著計(jì)算機(jī)技術(shù)的不斷進(jìn)步,最近20年來(lái),基于偏微分方程(PDE)的圖像去噪方法得到了很快的發(fā)展[1-2].傳統(tǒng)的PDE去噪方法主要是濾除圖像的高頻成分,然而由于圖像的細(xì)節(jié)也分布在高頻區(qū)域,所以總是在濾除噪聲的同時(shí)模糊了圖像的邊緣.為了能在濾除噪聲的同時(shí)也能保持圖像的邊緣,Rudin等[3]提出了如下的全變分(Total Variation)去噪模型(ROF模型):

式中:Ω∈R2,u為待修復(fù)的原始圖像,f為噪聲圖像,為正則化因子.然而,全變分正則化模型經(jīng)常在圖像恢復(fù)過(guò)程中出現(xiàn)階梯現(xiàn)象.為了克服階梯現(xiàn)象,文[4-6]考慮下面的高階PDE去噪模型(LLT模型):

式中:為正則化因子.但是由于高階PDE演化圖像的邊緣比全變分模型快,所以經(jīng)常在圖像的邊緣引起模糊現(xiàn)象.

為了在去除噪聲的同時(shí)不但保持圖像的邊緣,而且避免階梯現(xiàn)象,Lysaker等[7]建議下面的兩步算法(LOT模型):

第1步:令n=(n1,n2)擬合噪聲圖像的法向量,即:

第2步:令n0=(n0x,n0y)為式(3)的解,求解下面的優(yōu)化問(wèn)題:

事實(shí)上,式(4)和式(1)以及式(2)含有的L1項(xiàng)使得數(shù)值求解并不容易,從而尋求有效的數(shù)值解法解此類問(wèn)題是當(dāng)前的一個(gè)重要研究?jī)?nèi)容.最近,基于Bregman迭代算法,Goldstein和Osher[8]建議用分裂Bregman方法解帶有L1項(xiàng)的圖像恢復(fù)問(wèn)題.由于此算法具有編程簡(jiǎn)單、數(shù)值求解過(guò)程比較穩(wěn)定、在計(jì)算過(guò)程中保持正則化參數(shù)為一個(gè)常數(shù)、占有內(nèi)存小且具有較快的計(jì)算速度和收斂速度等優(yōu)勢(shì),廣泛應(yīng)用于圖像恢復(fù)問(wèn)題[8-11].因此,基于分裂Bregman方法能有效地處理含有L1項(xiàng)的優(yōu)化問(wèn)題,本文提出用該方法來(lái)解LOT去噪模型的第2步,也就是優(yōu)化問(wèn)題(4).

1 圖像去噪LOT模型的分裂Bregman算法

分裂Bregman方法由于能有效處理含有L1項(xiàng)的優(yōu)化問(wèn)題并且在數(shù)值求解過(guò)程中比較穩(wěn)定,因此最近被廣泛應(yīng)用到圖像復(fù)原、圖像修補(bǔ)、圖像壓縮等領(lǐng)域.本節(jié)首先回顧文獻(xiàn)[8]中的分裂Bregman方法,然后提出用該方法來(lái)解LOT去噪模型的第2步式(4),并給出具體的算法.

1.1 分裂Bregman算法

式中:ω為正則化參數(shù).很明顯,問(wèn)題(5)在一定條件下等價(jià)于下面的約束優(yōu)化問(wèn)題:

考慮下面的優(yōu)化問(wèn)題

對(duì)于式(5),由于含有兩個(gè)變量u,d,因此利用文[9]中的思想,引入一個(gè)輔助變量c可以通過(guò)下面的交替算法求解:

所以,我們有下面的分裂Bregman算法解優(yōu)化問(wèn)題(7).

算法1 分裂Bregman算法

步0 取初值u0=f,d0=0,c0=0,令k:=0;

步1 計(jì)算式(7)求得(uk+1,dk+1,ck+1);

步2 若終止條件滿足,則停止;否則,令k:=k+1,轉(zhuǎn)步1.

注 在算法1中,一般情況下H(u)為連續(xù)可微的,因此我們可以直接求出式(7)的第1個(gè)子問(wèn)題的歐拉方程,然后利用一般的數(shù)值方法求解.然而雖然式(7)的第2個(gè)子問(wèn)題不可微,但是該問(wèn)題是嚴(yán)格凸問(wèn)題并且含有一個(gè)不可微項(xiàng)|d|L1(Ω),因此在求歐拉方程時(shí),我們需要用到廣義導(dǎo)數(shù).經(jīng)過(guò)簡(jiǎn)單的計(jì)算并引入閾值算子有:

其中:

1.2 LOT模型的分裂Bregman方法

由于LOT模型的第1步(3)中含有1個(gè)等式約束|n|=1,使得求解并不容易,因此文[7]中引入一個(gè)輔助變量θ,使n=(nx,ny)=(cos θ,sin θ).另一方面,由于|n|=|θ|,因此若令n:=(,)=則問(wèn)題(3)可轉(zhuǎn)化為下面的無(wú)約束問(wèn)題:

進(jìn)一步,式(8)對(duì)應(yīng)的歐拉方程為:

利用下面的事實(shí):

文[7]中建議解下面的穩(wěn)態(tài)方程

現(xiàn)在考慮LOT模型的第2步,若令n0=(nx,ny)且設(shè)H(u)=‖u-f‖(Ω),并引入輔助變量d=u,類似于1.1節(jié)的分裂Bregman迭代算法,優(yōu)化問(wèn)題(4)可寫(xiě)成如下類似式(7)的形式:

事實(shí)上,式(12)和式(7)的第2個(gè)優(yōu)化問(wèn)題有所不同,很明顯,在式(12)中,(n0,d)L2(Ω)+‖u-dk-ck‖關(guān)于d是可微的.因此結(jié)合式(10)和式(11)以及算法1,對(duì)圖像去噪LOT模型,我們提出下面的分裂Bregman算法:

算法2 圖像去噪LOT模型的分裂Bregman算法.

第1步:利用半隱式時(shí)間演化法式(10)和式(11).

步0 取初值n0=(,),令k:=0;

其中Δt為時(shí)間步長(zhǎng)

第2步:利用分裂Bregman算法(12).

為獲得T-Map的3維空間域，依次將規(guī)范重心坐標(biāo)λF表達(dá)式中限定自由度方向的和設(shè)定為“0”，實(shí)現(xiàn)T-Map維數(shù)由4維降為3維。

步1 取初值u0=f,d0=0和c0=0,令k:=0;

步2 通過(guò)解式(15)～式(19),計(jì)算(uk+1,dk+1,ck+1):

其中:

步3 若停止準(zhǔn)則滿足,則算法停止,輸出復(fù)原圖像u;否則,令k:=k+1,轉(zhuǎn)步2.

注:事實(shí)上,算法2中的式(18)對(duì)應(yīng)式(15)的第1個(gè)式子且對(duì)應(yīng)的歐拉方程為:

其中:I和Δ分別為單位算子和Laplace算子.因此在離散情況下,我們可以用Gauss-Seidel迭代:

具體算法流程圖如圖1所示.

2 實(shí)驗(yàn)結(jié)果及其評(píng)價(jià)

本節(jié)進(jìn)行相關(guān)的數(shù)值試驗(yàn)來(lái)說(shuō)明所提出算法的有效性.實(shí)驗(yàn)平臺(tái)為:Intel core(TM)2 DUO CPU,2.20 GHz,2G內(nèi)存,Windows XP操作系統(tǒng),Matlab7.5軟件.

圖1 算法2流程圖Fig.1 The flow diagam of algorithm 2

首先我們比較LOT模型與經(jīng)典的ROF模型以及LLT模型的去噪效果.考慮添加標(biāo)準(zhǔn)差為10的高斯白噪聲的Lena圖像面部(圖2(a)和圖2(b)).

圖2 Lena圖像面部Fig.2 The face of lema image

由于Lena圖像的面部含有平滑區(qū)域,因此在使用ROF模型去噪時(shí)經(jīng)常出現(xiàn)階梯現(xiàn)象.另一方面,由于圖像的邊緣區(qū)域出現(xiàn)圖像像素值的跳躍,因此在使用LLT模型去噪時(shí)這些區(qū)域經(jīng)常出現(xiàn)模糊.為了更好地比較這3種模型,對(duì)于ROF模型(1)以及LLT模型(2),我們建議用文[3,6]的半隱式梯度下降法,其中相關(guān)的正則化參數(shù)和時(shí)間步長(zhǎng)t依次為:λ=0.08,tROF=0.02和α=0.275,tLLT=0.075.對(duì)于LOT模型的第1步,我們?nèi)≌齽t化參數(shù)β=1.125和時(shí)間步長(zhǎng)t=0.08;在第2步中,我們建議用分裂Bregman算法,其中γ=0.115,μ=0.002.若設(shè)定迭代次數(shù)為200次,則從圖2中可以看出LOT模型復(fù)原的圖像(e)在一定程度上克服了ROF模型復(fù)原圖像(c)中的階梯現(xiàn)象,又能避免LLT模型復(fù)原圖像(d)中的邊緣模糊現(xiàn)象,尤其是在圖像中帽檐區(qū)域.為了更好地理解這3種模型的復(fù)原效果,我們繪出了復(fù)原圖像與噪聲圖像的差.很明顯,如圖2所示(e1)同時(shí)具有(d1)和(c1)的特征.

下面我們比較文[7]中的算法和本文提出的分裂Bregman算法,考慮添加方差為12的白色高斯噪聲的Cameraman圖像(圖3).Cameraman圖像中不但含有像素跳躍區(qū)域(如:相機(jī)支架),而且也含有圖像漸變區(qū)域(如:天空).對(duì)于LOT模型的第一步:法向擬合步,我們?nèi)匀挥梦腫7]中的半隱式梯度下降法,其中相關(guān)的參數(shù)依次設(shè)定為β=1.125和時(shí)間步長(zhǎng)t=0.1,迭代次數(shù)為300次.現(xiàn)在考慮復(fù)原圖像步,即:LOT模型第2步,其中正則化參數(shù)γ=0.135.設(shè)定文[8]中的用梯度下降法(GD)解式(4)的時(shí)間步長(zhǎng)t=0.1,用分裂Bregman方法解(4)的正則化參數(shù)μ=0.004,若相鄰兩次迭代結(jié)果滿足8.25×10-4,則第2步迭代終止.

圖3 Cameraman圖像Fig.3 Cameraman image

如圖3所示,兩種算法幾乎有著同樣的復(fù)原結(jié)果.另外,從表1中可以看出這兩種算法在滿足終止標(biāo)準(zhǔn)時(shí),文[7]中的原始梯度下降法需要57次迭代,大約耗時(shí)2.890 6 s,而分裂Bregman迭代算法僅僅需要27次迭代,耗時(shí)大約1.828 1 s.另外,從收斂曲線圖中(圖4)可以看出分裂Bregman方法有著較快的收斂速度.同時(shí)還可以看出此時(shí)分裂Bregman方法有著較好的修復(fù)效果.

表1 實(shí)驗(yàn)Cameraman的相關(guān)結(jié)果Tab.1 The pelated results of cameraman experiments

圖4 收斂曲線Fig.4 Convergence carves

3 結(jié) 語(yǔ)

由于分裂Bregman方法能有效地處理含有L1項(xiàng)的優(yōu)化問(wèn)題,因此本文利用分裂Bregman方法來(lái)處理LOT模型的第2步.實(shí)驗(yàn)結(jié)果表明,與傳統(tǒng)的梯度下降法相比較,該方法不但有著較快的收斂速度,而且也能保持圖像的平滑區(qū)域和邊界區(qū)域的信息.另外,由于LOT模型的第1步含有一個(gè)等式約束,因此值得進(jìn)一步研究有效地處理第1步的數(shù)值方法.

[1] CHAN T,SHEN J.Image processing and analysis:Variational,PDE,wavelet and stochastic methods[M].SIAM,2005:145-198.

[2] A UBERT G,KORNPROBST P.Mathematical problems in image processing:partial differential equations and the calculus of variations[M].New York:Springer Verlag,2002:67-136.

[3] RUDIN L,OSHER S,FAT EMI E.Nolinear total variation based noise removal algorithms[J].Physica D,1992,60:259-268.

[4] STEIDI G.A note on the dual treatment of higher order regularization functionals[J].Computing,2006,76:135-148.

[5] SETZER S,STEIDI G.Variational methods with higher-order derivatives in image processing[C]//Brentwood Approximation,Nashboro Press,2008,12:360-386.

[6] LYSAKER M,LUNDERVOLD A,LUNDERVOLD,et al.Noise removal using fourth-order partial differential equation with applications to medical magnetic resonance images in space and time[J].IEEE T rans Image Process,2003,12:1579-1590.

[7] LYSAKER M,OSHER S,TAI X.Noise removal using smoothed normals and surface fitting[J].IEEE T rans Image Proce,2004,13(10):1345-1357.

[8] GOLDSTEIN T,OSHER S.The split Bregman for L1 regularized problems[J].SIAM Jornal on Imaging Sciences,2009,2(2):323-343.

[9] CAI Ji,OSHER S,SHEN Z.Split Bregman methods and frame based image restoration[R].UCLA CAM:Report,2009:9-28.

[10]TAI X,WU C.Augmented lagrangian method,dual methods and split Bregman iteration fo r ROF model[R].UCLA CAM:Report,2009.

[11]SETZER S.Split bregman algorithm,Douglas-Rachford splitting and frame shrinkage[C]//Scale Space and Variational M ethods in Computer Vision,Berlin:Springer,2009,5567:464-476.