一種自適應(yīng)的全變分圖像去噪算法

一種自適應(yīng)的全變分圖像去噪算法

1非線性圖像去噪方法圖像噪聲去除是數(shù)據(jù)處理的重要部分，在實(shí)際的圖像采集、傳輸和存儲過程中不可避免地會引入各種噪聲。噪聲有很多類型和原因。在大多數(shù)情況下，需要處理圖像以消除噪聲，以便將處理后的圖像適合分析和提取信息。有很多方法可以去除噪聲。在頻域處理中，有以高濾波為代表的線性濾波和以中值濾波為代表的非線性濾波。由于其簡單、復(fù)雜、低噪聲計(jì)算，高濾波是一種廣泛應(yīng)用的噪聲去除方法，但由于噪聲和信號的處理，它可以去除噪聲，削弱圖像的亮度，并在圖像邊緣移動(dòng)。因此，在這個(gè)問題上，提出了許多非線性噪聲去除方法。其中，基于微分法的圖像消噪是一種具有代表性的方法。由于這一方法因?yàn)槠x差分方程是數(shù)學(xué)的一個(gè)重要分支，因此形成了一個(gè)完整的理論體系和數(shù)值方法。同時(shí)，該方法從新的角度（如能量擴(kuò)散、橫向發(fā)展等）來解釋圖像噪聲的去除過程。其中，這種方法的代表性是pm方程和全變分算法（圖像遺傳，模型）。本文主要討論的是TV算法,該方法自1992年由Rudin、Osher和Fatime首次提出后,因其去噪效果明顯而引起廣泛關(guān)注,但這種方法也有自身的缺點(diǎn),如階梯效應(yīng)、需要已知圖像的噪聲方差等.針對這些問題,本文提出一種自適應(yīng)全變分圖像去噪算法,該算法主要是修改傳統(tǒng)TV算法中的逼近項(xiàng),用受噪聲影響較小的模糊圖像代替原始的含噪聲圖像,使得在整個(gè)圖像處理的計(jì)算過程中不需要已知圖像的噪聲方差,便可較好地對圖像進(jìn)行去噪處理,達(dá)到顯著的去噪效果,同時(shí)依據(jù)邊緣信息對模型中拉格朗日因子進(jìn)行加權(quán)處理,使得處理后不僅不會模糊圖像的邊緣,而且較好地抑制傳統(tǒng)TV算法所產(chǎn)生的階梯效應(yīng).2算法的設(shè)計(jì)和實(shí)現(xiàn)2.1tv模型中的抽象函數(shù)傳統(tǒng)TV算法就是最小化能量泛函:E=∫Ω(|?u|+12λ(u-u0)2)dxdy,(1)E=∫Ω(|?u|+12λ(u?u0)2)dxdy,(1)其中,∫Ω|?u|dxdy∫Ω|?u|dxdy和∫Ω12λ(u-u0)2dxdy∫Ω12λ(u?u0)2dxdy分別為TV模型的正則項(xiàng)(RegularizationTerm)和逼近項(xiàng)(FidelityTerm),Ω為圖像區(qū)域,λ為拉格朗日因子,u0為初始含有噪聲的圖像,u為經(jīng)處理后的圖像,式(1)的歐拉-拉格朗日(Euler-Lagrange,E-L)方程如下:div(?u|?u|)+λ(u0-u)=0.(?u|?u|)+λ(u0?u)=0.由最速下降法(SteepestDescentMethod)可解得圖像以時(shí)間為演化參數(shù)的演化公式如下:{ut=div(?u|?u|)+λ(u0-u)u|t=0=u0λ=1σ2|Ω|∫Ωdiv(?u|?u|)(u-u0)dxdy?????????????ut=div(?u|?u|)+λ(u0?u)u|t=0=u0λ=1σ2|Ω|∫Ωdiv(?u|?u|)(u?u0)dxdy其中,|Ω||Ω|為圖像面積,λ也被稱為是全局的尺度因子.從上式可看出,TV模型中關(guān)于逼近項(xiàng)的演化公式存在著2個(gè)問題:1)λ的求解需要已知圖像噪聲方差σ2,并且λ是一個(gè)全局變量,但圖像中不同部分(如邊緣、背景等)的信噪比是不同的,用相同的λ作為逼近項(xiàng)解前的系數(shù)顯然不是一個(gè)理想的選擇;2)當(dāng)原圖像噪聲較大(σ2較大)時(shí),u0及其相關(guān)的逼近項(xiàng)就會存在較大誤差,同時(shí)σ2與λ成反比,λ值相對較小,使得演化公式主要由正則項(xiàng)的解div(?u|?u|)(?u|?u|)決定,所以在處理邊緣時(shí)很容易出現(xiàn)階梯效應(yīng).實(shí)際情況也是當(dāng)圖像噪聲越大時(shí),階梯效應(yīng)越明顯.2.2拉格朗日因子求解本文解決上述問題的出發(fā)點(diǎn)是逼近項(xiàng)中的u0,因?yàn)橛?jì)算λ需要的噪聲方差是由u0決定,當(dāng)噪聲較大時(shí),u0還會引起較大的誤差,所以本文對u0進(jìn)行一次高斯濾波得到一個(gè)較模糊的圖像ug=u0*G(G為高斯核),用ug代替u0,這樣逼近項(xiàng)就變成∫Ω12λ(u-ug)2dxdy∫Ω12λ(u?ug)2dxdy,新的能量泛函如下:Enew=∫Ω(|?u|+12λ(u-ug)2)dxdyEnew=∫Ω(|?u|+12λ(u?ug)2)dxdy,其歐拉-拉格朗日方程如下:div(?u|?u|)+λ(ug-u)=0,這樣得到新的演化公式:ut=div(?u|?u|)+λ(ug-u).(2)如2.1節(jié)小結(jié)所述,將λ作為一個(gè)全局變量不是一個(gè)理想的選擇,λ應(yīng)該根據(jù)圖像局部信息的不同而不同,所以將上式兩邊同時(shí)乘上(ug-u)并在局部圖像區(qū)域Ω0(Ω0是以當(dāng)前像素點(diǎn)為中心點(diǎn)大小為(2m+1)×(2m+1)的窗口)進(jìn)行積分,當(dāng)達(dá)到穩(wěn)定解時(shí),ut將趨于0,所以得0=∫Ω0div(?u|?u|)(ug-u)dxdy+λ∫Ω0(ug-u)2dxdy.(3)下面的重點(diǎn)是如何求解拉格朗日因子λ,而圖像局部信息主要是邊緣信息,所以接下來本文分邊緣區(qū)域和非邊緣區(qū)域進(jìn)行討論.1)邊緣區(qū)域.不失一般性,假設(shè)u0組成形式如下:u0=u*+n,其中,u*為無噪聲有效信號,n為均值為0的噪聲信號,由于高斯濾波的去噪作用,且在邊緣區(qū)域,圖像的有效信號功率通常要遠(yuǎn)大于噪聲信號功率,即P(u*)?P(n),所以可得ug=u0*G≈u**G.(4)同時(shí)對圖像信號的高斯濾波可用能量擴(kuò)散的角度來解釋,其過程同樣可用以時(shí)間為參數(shù)的演化公式來表示:u*t=u*xx+u*yy,(5)其中,u*t表示u*對時(shí)間的導(dǎo)數(shù),u*xx和u*yy表示u*對圖像坐標(biāo)方向的二階導(dǎo)數(shù),由式(4)和式(5)可得到一個(gè)ug的近似表達(dá)式:ug≈u**G=u*+u*t=u*+u*xx+u*yy,將上式代入式(3)等號右邊的第二項(xiàng),則第二項(xiàng)變?yōu)棣恕姚?(u*+u*xx+u*yy-u)2dxdy,將上式分解為注意到u*為無噪聲有效信號,顯然u演化的目標(biāo)便是不斷地趨近于u*,所以當(dāng)達(dá)到穩(wěn)定解的時(shí)候,新的模型中存在一個(gè)約束條件:∫Ω0(u-u*)2dxdy=0,上式的條件表示當(dāng)ut將趨于0時(shí),u將趨于u*,等價(jià)對于任一常數(shù)C,∫Ω0C(u-u*)dxdy=0,此時(shí),式(6)的第一項(xiàng)和第三項(xiàng)被消除,同時(shí)當(dāng)達(dá)到穩(wěn)定解的時(shí)候,u可近似代替u*,所以式(6)可近似等價(jià)為λ∫Ω0(uxx+uyy)2dxdy,將上式替換式(3)等號右邊的第二項(xiàng),可得到求解λ表達(dá)式為λ=∫Ω0div(?u|?u|)(u-ug)dxdy∫Ω0(uxx+uyy)2dxdy+ε?(7)其中,ε是一個(gè)大于0的常數(shù),引入該常數(shù)是為了避免出現(xiàn)∫Ω0(uxx+uyy)2dxdy取值接近0的情況.2)非邊緣區(qū)域.非邊緣區(qū)域是圖像的背景和平坦區(qū)域,其信號主要是由噪聲組成,而且不需要考慮邊緣區(qū)域存在的階梯效應(yīng).在非邊緣區(qū)域,并不能按式(7)計(jì)算拉格朗日因子λ,因?yàn)槭?7)的成立依賴于式(4)的近似u0*G≈u**G,而這一近似表達(dá)成立的前提條件是在邊緣區(qū)域圖像的有效信號功率通常要遠(yuǎn)大于噪聲信號功率,即P(u*)?P(n),而在非邊緣區(qū)域這一條件很可能是不成立的,尤其是在噪聲較大的時(shí)候,所以在非邊緣區(qū)域圖像的演化公式應(yīng)該有所不同.本文從演化公式中正則項(xiàng)與逼近項(xiàng)的作用來考慮非邊緣區(qū)域的圖像演化.圖像噪聲的去除是通過正則項(xiàng)來實(shí)現(xiàn),它現(xiàn)實(shí)對噪聲信號的平滑.雖然在平滑過程中正則項(xiàng)容許邊緣的不連續(xù),但逼近項(xiàng)的作用則是要求經(jīng)過演化后的圖像近似地趨于指定圖像ug(針對演化式(2)而言),從某種程度來說,正則項(xiàng)與逼近項(xiàng)作用是相反的,兩項(xiàng)經(jīng)過演化后以拉格朗日因子λ為比例權(quán)重到達(dá)平衡.在非邊緣區(qū)域中逼近項(xiàng)的作用實(shí)際上是多余的,因?yàn)樵谶@些區(qū)域中并不存在圖像的邊緣信息和紋理信息,所以使這些區(qū)域中的像素點(diǎn)與其周圍像素點(diǎn)共同趨于某一相同的灰度值是合理的,完全沒有必要趨于指定圖像ug對應(yīng)的灰度值.因此在圖像的非邊緣區(qū)域,演化公式可完全由正則項(xiàng)部分(div(?u|?u|))決定,這樣不僅可以避免拉格朗日因子λ不能由式(7)計(jì)算的問題,而且消除逼近項(xiàng)在演化公式中的作用,從而加速圖像的演化,有利于噪聲的去除,同時(shí)也不必考慮邊緣區(qū)域才會出現(xiàn)的階梯效應(yīng).2.3保證邊緣信號響應(yīng)的區(qū)域因?yàn)榍懊媸欠诌吘墔^(qū)域與非邊緣區(qū)域討論演化公式,所以如何區(qū)分邊緣區(qū)域與非邊緣區(qū)域就成了一個(gè)需要解決的問題.但是否需要準(zhǔn)確區(qū)分兩個(gè)區(qū)域,答案是否定的.首先前面所討論的演化公式并不是基于準(zhǔn)確區(qū)分兩個(gè)區(qū)域的,所以準(zhǔn)確區(qū)分沒有必要;其次,如果準(zhǔn)確區(qū)分兩個(gè)區(qū)域會直接導(dǎo)致兩個(gè)問題:兩個(gè)區(qū)域的演化公式不能統(tǒng)一;兩個(gè)區(qū)域的相鄰區(qū)域很可能存在相鄰且灰度值相近的像素點(diǎn)卻按照差異很大的演化公式進(jìn)行演化,從而使演化的結(jié)果差異很大,顯然這也是不合理的.所以本文需要一種度量來表征圖像中某一像素點(diǎn)屬于邊緣區(qū)域或非邊緣區(qū)域的程度,在無噪聲情況下,這種度量必須滿足2個(gè)性質(zhì):1)單調(diào)性,它需要保證其對邊緣信號的響應(yīng)在邊緣區(qū)域大于(或小于)在非邊緣區(qū)域;2)連續(xù)性,即當(dāng)從邊緣區(qū)域過渡到非邊緣區(qū)域(或相反方向),該度量在數(shù)值是連續(xù)變化的.在有噪聲的情況下,這種度量即使不能嚴(yán)格滿足這兩個(gè)性質(zhì),也需要近似滿足這兩個(gè)性質(zhì).從本文后面的結(jié)果可看出,采用滿足這兩個(gè)性質(zhì)的度量對演化公式進(jìn)行修正,可最好地統(tǒng)一邊緣區(qū)域與非邊緣區(qū)域的演化公式,同時(shí)也可避免相鄰且灰度值相近的像素點(diǎn)因演化公式的巨大差異而產(chǎn)生差異很大的演化結(jié)果.本文用局部方差作為上述的度量,顯然局部方差在邊緣區(qū)域?qū)吘壭盘柕捻憫?yīng)要大于非邊緣區(qū)域,且其也滿足上述所說的連續(xù)性,當(dāng)然為了最大程度避免噪聲的影響,同時(shí)基于我們不需要準(zhǔn)確區(qū)分邊緣區(qū)域與非邊緣區(qū)域這一前提,本文并不對u0計(jì)算局部方差,而對模糊圖像ug計(jì)算局部方差,計(jì)算結(jié)果為V,計(jì)算局部方差的區(qū)域是以當(dāng)前像素點(diǎn)為中心點(diǎn)大小為(2l+1)×(2l+1)的窗口,如果當(dāng)前像素點(diǎn)的圖像坐標(biāo)為(x,y),則V(x,y)計(jì)算表達(dá)式如下:V(x,y)=1(2l+1)2i=l∑i=-lj=l∑j=-l[ug(x+i,y+j)-1(2l+1)2i=l∑i=-lj=l∑j=-lug(x+i,y+j)]2.(8)但不同圖像的局部方差在數(shù)值上存在明顯的差異,所以本文對上述的局部方差V作了歸一化處理,使V線性映射到[01]的區(qū)間內(nèi),其最終結(jié)果為—V.顯然當(dāng)歸一化后的局部方差越趨于0,說明當(dāng)前位置的像素點(diǎn)越可能屬于非邊緣區(qū)域,反之當(dāng)該值越趨于1,說明當(dāng)前位置的像素點(diǎn)越可能屬于邊緣區(qū)域.所以用—V乘以演化公式中的逼近項(xiàng)部分,便得到了統(tǒng)一的演化公式:{ut=div(?u|?u|)+λ—V(ug-u)u|t=0=u0λ=∫Ω0div(?u|?u|)(u-ug)dxdy∫Ω0(uxx+uyy)2dxdy+ε(9)從上式可看出,在邊緣區(qū)域時(shí),—V→1,上式由式(2)和式(7)決定,符合邊緣區(qū)域的演化過程,在非邊緣區(qū)域時(shí),—V→0,上式由演化公式的正則項(xiàng)部分決定,這也同樣符合非邊緣區(qū)域的演化過程.3邊緣區(qū)域犯罪性檢驗(yàn)為驗(yàn)證本文算法,本文用兩幅8位的灰度圖像(圖1(b))進(jìn)行了實(shí)驗(yàn)((b)分別是在(a)上加入均值為0,方差為50的高斯噪聲).并將本文算法與傳統(tǒng)TV算法在整幅圖像與相同位置的局部放大圖像的結(jié)果進(jìn)行對比,如圖1所示.需要說明的是,在本文的實(shí)驗(yàn)中式(9)中局部圖像區(qū)域Ω0與式(8)計(jì)算局部方差的區(qū)域均選擇相同大小為5×5的窗口,式(9)中常量ε取值為1,并且實(shí)驗(yàn)的計(jì)算過程中本文算法與傳統(tǒng)TV算法均選取相同的迭代規(guī)則和終止條件,在同一個(gè)λ下迭代10次,迭代步長為0.2,當(dāng)圖像平均每一個(gè)像素的灰度值改變量小于0.08時(shí)終止迭代,各結(jié)果圖λ的迭代次數(shù)如下所示.圖1(c)左圖λ的迭代次數(shù)為76,右圖為74;圖1(d)左圖λ的迭代次數(shù)為77,右圖為70.假設(shè)圖像的大小為M×N,則傳統(tǒng)TV算法與本文算法每迭代一次λ所對應(yīng)的計(jì)算復(fù)雜度為O(MN)和O((2m+1)2MN).從實(shí)驗(yàn)結(jié)果圖可很顯然地看出本文算法在不需要已知圖像噪聲方差的情況,與傳統(tǒng)TV算法(如果不知道噪聲方差,傳統(tǒng)TV算法是無法達(dá)到如圖1(c)所示的去噪效果)一樣可較好去除圖像的噪聲,不會像線性濾波那樣模糊圖像邊緣,較好地保持邊緣的對比度,同時(shí)還可看出傳統(tǒng)TV算法結(jié)果的階梯效應(yīng)要比本文算法明顯的多,這一點(diǎn)從圖1中兩者的局部放大圖(e)和(f)可很明顯的看出.所以本文算法從已知條件約束和最終的結(jié)果圖來看要明顯地優(yōu)于傳統(tǒng)TV算法.接下來本文對圖像處理的結(jié)果進(jìn)行評價(jià),顯然信噪比(SignaltoNoiseRatio,SNR)是評價(jià)圖像去噪效果的一個(gè)重要指標(biāo),但僅用這一個(gè)指標(biāo)是不夠的,因?yàn)镾NR并不能反應(yīng)出圖像邊緣區(qū)域的圖像質(zhì)量,如果去噪算法產(chǎn)生較明顯的階梯效應(yīng),就會產(chǎn)生不規(guī)則的人造邊緣,這樣圖像質(zhì)量明顯是不好的,但由于階梯效應(yīng)的產(chǎn)生,使得這些邊緣處會產(chǎn)生更大的灰度值變化,因此它的信噪比也會更大.所以本文需要一個(gè)輔助指標(biāo),在圖像噪聲得到明顯的改善,SNR有了大幅度提升,圖像邊緣沒有被模糊的前提下,來評價(jià)圖像邊緣的平滑情況,因?yàn)殡A梯效應(yīng)越明顯,邊緣平滑性越差.以∫Ω|?u|2dxdy作為一幅圖像平滑性的度量是一種合理的選擇,因?yàn)樗怀湓S邊緣有跳變,實(shí)際如果以∫Ω|?u|2dxdy作為圖像的能量泛函并使其最小化,那么圖像的演化公式就是線性擴(kuò)散方程,如式(5),即高斯濾波,所以∫Ω|?u|2dxdy的值越小說明圖像越平滑.但本文需要考慮的是邊緣區(qū)域的平滑情況,所以需要對∫Ω|?u|2dxdy進(jìn)行加權(quán)平均處理,邊緣區(qū)域權(quán)重高,非邊緣區(qū)域權(quán)重低,2.3節(jié)所提到的歸一化局部方差ˉV滿足這權(quán)重關(guān)系.結(jié)合上面的敘述,本文提出邊緣區(qū)域的平滑性(SmoothnessofEd