周 徽
(中國電子科技集團(tuán)公司第二十八研究所,南京 210007)
數(shù)字圖像復(fù)原/恢復(fù)技術(shù)(以下簡稱復(fù)原技術(shù))是數(shù)字圖像處理的重要研究內(nèi)容,已在諸多領(lǐng)域取得了成功應(yīng)用[1]。
由于現(xiàn)有的圖像系統(tǒng)都存在某種程度的缺陷,呈現(xiàn)出各種退化因素,因而進(jìn)行圖像復(fù)原是必要的??傮w上退化因素可以歸結(jié)為2類。①空間退化,又稱為模糊。成像系統(tǒng)中某些元件失常、圖像傳感器與成像目標(biāo)之間的相對(duì)移動(dòng)、攝像機(jī)失焦、不良天氣影響等因素都會(huì)造成圖像模糊。②點(diǎn)退化,通常指噪聲。最常見的點(diǎn)退化來自電子元件或數(shù)字化噪聲。
圖像退化過程可被模型化為一個(gè)退化函數(shù)和一個(gè)加性噪聲項(xiàng)共同作用于原始圖像f(x,y),從而產(chǎn)生一幅退化的圖像g(x,y),如圖1所示。根據(jù)該模型,退化圖像的數(shù)學(xué)描述見式(1)。
其中:g(x,y)代表退化圖像;f(x,y)代表未退化圖像;n(x,y)代表噪聲項(xiàng);h(x,y)代表模糊算子,又稱點(diǎn)擴(kuò)展函數(shù)(point spread function,PSF);符號(hào)“*”代表卷積。為了便于描述,式(1)改寫為矩陣形式
其中:f、g和n為NM×1維向量,由二維圖像按行或者列堆積而成;H是由PSF生成的MN×MN塊循環(huán)矩陣。
圖1 基本圖像退化/復(fù)原模型
理論上,在給定f(x,y)以及退化函數(shù)H和噪聲n(x,y)的一些先驗(yàn)知識(shí)后,便可以獲得原始圖像的一個(gè)近似估計(jì)f^(x,y)。然而,式(2)具有病態(tài)特性,導(dǎo)致直接求解不容易實(shí)現(xiàn),只能先設(shè)法對(duì)式(2)施加一定的約束,將其病態(tài)轉(zhuǎn)為良態(tài)后再設(shè)法求其最優(yōu)估計(jì)解。這個(gè)估計(jì)應(yīng)能盡可能接近原始圖像。
針對(duì)該問題,人們提出過許多方法,如在已知g、H與n有關(guān)的先驗(yàn)知識(shí)的前提下,去求原圖像的逼近解。在這些方法當(dāng)中,Tikhonov正則化方法最具優(yōu)勢[2]。該方法是一個(gè)在L2空間的最小方差近似問題,能夠平衡精確性與平滑性[3],因而得到廣泛研究。它將式(2)變?yōu)閹Ъs束條件的優(yōu)化解,即
其中:ψdata(f,g)表示數(shù)據(jù)逼近項(xiàng);φreg(Qf)表示正則化項(xiàng);Q為高通濾波算子;d為正則化參數(shù),用以控制數(shù)據(jù)逼近項(xiàng)與正則化項(xiàng)的空間分布。理論上,在 L(α,f)最小化,即條件下的解就是原圖像f(x,y)的最優(yōu)逼近解。
在正則化方法應(yīng)用中確定正則化參數(shù)是非常關(guān)鍵的一步。人們先后提出過許多關(guān)于如何選擇正則化參數(shù)的方法。Galatsanos等[4]提出了基于均方誤差(MSE)評(píng)判標(biāo)準(zhǔn)的正則化參數(shù)選擇和噪聲變量估計(jì)方法。Nakano等[5]提出用權(quán)矩陣的方式優(yōu)化正則化效果,而局部正則化優(yōu)化權(quán)矩陣從迭代過程獲得。Reeves[6]則給出了在限制性條件下,用交叉有效性準(zhǔn)則判斷正則化最優(yōu)度的方法。Molina等[7]將等級(jí)貝葉斯方法應(yīng)用到圖像復(fù)原中,提出了一個(gè)3步方法。該方法在每個(gè)步驟中完成1個(gè)正則參數(shù)的確定,降低了對(duì)噪聲先驗(yàn)知識(shí)的依賴性。Leung[8]提出用L曲線作為選取正則化參數(shù)的依據(jù),通過實(shí)驗(yàn)驗(yàn)證了L曲線的曲率與正則化參數(shù)優(yōu)異性之間的相關(guān)性,其結(jié)論是:當(dāng)曲率最大化時(shí),取得的參數(shù)是最優(yōu)的。Demoment[9]綜述了正則化方法的詳細(xì)內(nèi)容。
上面提到的正則化參數(shù)選取方法都具有全局性限制,即參數(shù)在算法開始之前就已經(jīng)確定,導(dǎo)致過分依賴于先驗(yàn)信息,無法在算法迭代過程中自適應(yīng)修正。與上面的方法不同,Katsaggelos等[10]提出用正則化函數(shù)取代全局正則化參數(shù)的方法,使算法能夠自適應(yīng)地修改正則化參數(shù),大幅提高了算法的效率。本文的工作就是充分利用該思想,提出一種改進(jìn)的自適應(yīng)加權(quán)正則化迭代圖像復(fù)原方法,增強(qiáng)算法對(duì)強(qiáng)噪聲存在時(shí)的適應(yīng)性。
用含有自變量f的正則化函數(shù)α(f)代替全局正則化參數(shù)α,建立目標(biāo)平滑泛函式(4)。
目標(biāo)泛函建立之后,正則化參數(shù)的選擇至關(guān)重要。自適應(yīng)正則化函數(shù)α(f)的選取應(yīng)滿足3條性質(zhì)[10]:① α(f)為平滑泛函的函數(shù);② α(f)與目標(biāo)泛函之間需滿足一定的極值性質(zhì),即當(dāng)時(shí),α(f)→0;‖g - Hf‖時(shí),α(f)→∞;③ 在α(f)取值范圍內(nèi)不改變目標(biāo)泛函的凸函數(shù)性質(zhì)。
根據(jù)上面性質(zhì),假設(shè)存在正實(shí)數(shù)λ使得式(5)成立,則
正則化項(xiàng)用于懲罰解的粗糙性,然而過分的懲罰會(huì)造成解的過分平滑,抑制算法對(duì)高頻成分的恢復(fù),從而導(dǎo)致復(fù)原結(jié)果在圖像強(qiáng)度變化劇烈處產(chǎn)生振鈴效應(yīng),影響圖像視覺效果。因此,在復(fù)原濾波時(shí),必須將圖像局部細(xì)節(jié)特性考慮在內(nèi),采取加權(quán)形式增強(qiáng)算法的局部適應(yīng)能力,提高算法對(duì)振鈴的抑制能力,改善復(fù)原質(zhì)量[5]。
文獻(xiàn)[11]認(rèn)為,圖像中小的梯度值對(duì)應(yīng)噪聲,需對(duì)其進(jìn)行平滑,而大的梯度值則對(duì)應(yīng)圖像邊緣輪廓,應(yīng)給予保護(hù)。根據(jù)該理論,定義加權(quán)正則化項(xiàng)為
其中W是由wij構(gòu)成的矩陣,wij定義為
加權(quán)后的自適應(yīng)正則化函數(shù)如式(10)所示。
圖像復(fù)原過程等價(jià)于式(4)定義的平滑函數(shù)逼近最優(yōu)解的過程。由于目標(biāo)泛函的凸性,恰當(dāng)?shù)剡x取正則化參數(shù)能夠使式(4)存在唯一的全局最優(yōu)解。該最優(yōu)解在平滑函數(shù)梯度為0時(shí)取得,即
將式(11)對(duì)f求導(dǎo),得
因?yàn)棣?fk)是第k次迭代結(jié)束后計(jì)算出的參數(shù),用于第k+1次迭代,所以▽fα(f)=0,從而式(12)化簡為
建立逐步逼近的Van Cittert格式的迭代等式:
首先分析算法的收斂性,改寫式(14)為
對(duì)式(15)左右兩邊同時(shí)取范數(shù),并根據(jù)矩陣范數(shù)與向量范數(shù)的協(xié)調(diào)性,有
反復(fù)使用這個(gè)不等式得
下面分析算法的收斂條件。因?yàn)椤珹‖<1等價(jià)于 ρ(A)< 1,有
其中ρ(·)定義為“·”的譜半徑。根據(jù)三角不等式,上式等價(jià)于
通過規(guī)范化H、C可使ρ(HTH)和ρ(CTC)均不大于1,從而算法收斂的充分條件演變?yōu)棣?f)<1,即
算法流程描述:
步驟1 設(shè)置初始值f0=g,迭代計(jì)數(shù)值k=0,設(shè)定算法終止條件ε=10-7,開定時(shí)器。
步驟2 在第k步迭代中,自動(dòng)修正正則化參數(shù)αk(fk),并利用式(14)對(duì)圖像進(jìn)行復(fù)原濾波。更新fk+1。
步驟4 關(guān)定時(shí)器,計(jì)算復(fù)原評(píng)價(jià)指標(biāo)ISNR,顯示復(fù)原結(jié)果。
選擇2組實(shí)驗(yàn)驗(yàn)證本文提出算法的有效性。第1組是樣條圖像,用來驗(yàn)證算法對(duì)圖像邊緣恢復(fù)的效果;第2組是經(jīng)典的lena圖像,用來驗(yàn)證算法對(duì)圖像細(xì)節(jié)恢復(fù)的效果。
實(shí)驗(yàn)1 樣條圖像尺寸為158×300的灰度圖像,如圖2(a)所示。本實(shí)驗(yàn)在PC機(jī)(CPU為P42.4GHz,內(nèi)存512 Mbytes)上進(jìn)行,模糊算子H取為11×11的均勻模糊,隨機(jī)疊加信噪比為20 dB的高斯噪聲,C取2-D拉普拉斯算子。判據(jù)‖fk+1-fk‖/‖fk‖ <10-7作為算法終止迭代的條件,改進(jìn)信噪比(improvement signal-to-noise ratio,ISNR)作為算法性能的評(píng)價(jià)指標(biāo),ISNR定義為
圖2 樣條圖像復(fù)原效果
圖2(a)、(c)、(e)分別為原始樣條圖像、退化圖像和復(fù)原圖像。圖2(b)、(d)、(f)分別為穿越圖2(a)、(c)、(e)中部的水平掃描線,其橫坐標(biāo)為水平象素點(diǎn),縱坐標(biāo)為對(duì)應(yīng)象素的亮度值。
由圖2(e)可以看出,相對(duì)于圖2(c),算法恢復(fù)的視覺效果是比較理想的。圖2(f)給出的復(fù)原圖像中部的掃描線說明樣條交接的輪廓線基本得以復(fù)原,顯然較圖2(d)有明顯改善。實(shí)驗(yàn)獲得的具體參數(shù)見表1和圖4。
實(shí)驗(yàn)2 lena圖像尺寸為300×300的灰度圖像,如圖3所示。其他條件同實(shí)驗(yàn)1。圖3(a)、(b)分別為退化的lena圖像和復(fù)原后的lena圖像。圖3(c)、(d)分別為圖3(a)、(b)的局部放大圖像。實(shí)驗(yàn)獲得的具體參數(shù)見表1和圖4。
圖3 lena圖像恢復(fù)效果
表1 實(shí)驗(yàn)結(jié)果參數(shù)
圖4 收斂誤差曲線
圖4給出2組實(shí)驗(yàn)中每一迭代步驟對(duì)應(yīng)的誤差判據(jù)變化曲線,其中圖4(a)是樣條圖像,圖4(b)是lena圖像,橫坐標(biāo)對(duì)應(yīng)迭代次數(shù),縱坐標(biāo)對(duì)應(yīng)誤差判據(jù)值。由圖4可以看出,無論迭代初始條件如何,算法能夠以較少的迭代次數(shù)快速平穩(wěn)地趨向于最優(yōu)估計(jì)解。
本文提出一種改進(jìn)的自適應(yīng)加權(quán)正則化迭代圖像復(fù)原算法。該算法能夠自適應(yīng)地選擇并自動(dòng)修正正則化參數(shù)。算法對(duì)初始值不敏感,在每步迭代中,算法完成1次復(fù)原濾波,并更新正則化參數(shù),確保結(jié)果能夠快速趨向于最優(yōu)。討論了算法的收斂性與收斂條件,通過建立顯式迭代復(fù)原算法尋找全局最優(yōu)解。該算法具有占用內(nèi)存小、計(jì)算速度快等優(yōu)點(diǎn)。無需任何噪聲信息作為先驗(yàn)知識(shí),算法通過觀測圖像本身自適應(yīng)估計(jì)并調(diào)整正則化參數(shù)以適應(yīng)強(qiáng)噪聲存在下的復(fù)原操作。用2組實(shí)驗(yàn)驗(yàn)證了算法的有效性。結(jié)果顯示:該算法能提高算法對(duì)強(qiáng)噪聲存在的適應(yīng)性,抑制Gibbs振鈴波紋,具有較好的評(píng)價(jià)結(jié)果和較高的實(shí)時(shí)性。
[1]Banham M R,Katsaggelos A K.Digital image restoration[J].Signal Processing,1997,14(2):24-41.
[2]Bouhamid A,Jbilou K.Sylvester Tikhonov-regularization methods in image restoration[J].Journal of Computational and Applied Mathematics,2007,206(1):86-98.
[3]Leung C M,Lu W S.Optimal determination of regularization parameters and the stabilizing operator[J].Signal Processing,1995,201:403-406.
[4]Galatsanos N P,Katsaggelos A K.Methods for choosing the regularization parameter and estimating the noise variance in image restoration and their relation[J].IEEE Transactions on Image Processing,1992,1(3):322-336.
[5]Nakano K,Eguchi M,Toyota Y,et al.On regularization for image restoration problems from the viewpoint of a Bayesian information criterion [C]//IECON-93.Hawaiian,USA:[s.n.],1993:2257-2261.
[6]Reeves S J.Optimal regularized image restoration with constraints[C]//ICASSP-92.San Francisco,USA:[s.n.],1992:301-304.
[7]Molina R,Katsaggelos A K,Mateos J.Bayesian and regularization methods for hyperparameter estimation in image restoration.Image Processing[J].IEEE Transactions on Image Processing,1999,8(2):231-246.
[8]Leung C M,Lu W S.An L-curve approach to optimal determination of regularization parameter in image restoration.Electrical and Computer Engineering[J].Electrical and Computer Engineering,1993,2:1021-1024.
[9]Demoment G.Image reconstruction and restoration:overview of common estimation structures and problems[J].IEEE Transactions on Signal Processing,1989,37(12):2024-2036.
[10]Moon G K,Katsaggelos A K.General choice of the regularization functional in regularized image restoration[J].IEEE Transactions on Image Processing,1995,4(5):594-602.
[11]Charbonnier N P,Blanc F L.Deterministic edge-preserving regularization in computed imaging[J].IEEE Transactions on Image Processing,1997,6(2):298-311.