郭立文，廖永忠

1(陜西國防工業(yè)職業(yè)技術(shù)學(xué)院 計算機(jī)與軟件學(xué)院，西安 710300) 2(湖南第一師范學(xué)院 信息科學(xué)與工程學(xué)院，長沙 410205)

1 引 言

成像過程中的相機(jī)晃動是圖像模糊的一個重要成因，去除這種運動模糊是圖像去模糊研究的一個重要研究方向，圖像去模糊過程本質(zhì)是對一類病態(tài)的反卷積方程的求解，如何高效、準(zhǔn)確的求解此類方程是圖像去模糊研究的難點問題.

圖像的盲去模糊是指在圖像的模糊核函數(shù)未知情況下的圖像去模糊方法，盲去模糊算法是當(dāng)前圖像恢復(fù)研究十分活躍的一個領(lǐng)域，一些學(xué)者在近幾年取得很多高水平的研究成果，相關(guān)的期刊和會議上有很多此類研究的文章[1-4].MIT的Levin提出一種基于邊緣概率密度函數(shù)的模糊核估計算法[5]，用高斯混合分布模型表示自然圖像的先驗分布，然后采用期望最大化(EM)算法交替迭代實現(xiàn)模糊圖像的恢復(fù).Taeg Sang Cho在文獻(xiàn)[6]中利用圖像直線邊緣在模糊過程中發(fā)生變化，其Radon變換與模糊核函數(shù)之間有確定的對應(yīng)關(guān)系，從而確定模糊核函數(shù)的方法，然后用最大后驗概率(MAP)完成模糊圖像的恢復(fù).Bae H提出一種快速的模糊圖像恢復(fù)算法[7]，其認(rèn)為由于自然圖像梯度的稀疏分布特性，圖像模糊信息都包含在圖像的邊界中，對圖像進(jìn)行分塊，然后提取一些包含較多信息邊界的塊，再把它們進(jìn)行拼接，利用這個模糊圖像塊完成模糊核的估計，該算法能降低計算量.Lin Zhong 提出強(qiáng)噪聲下大尺寸模糊核的圖像去模糊圖像算法，先利用方向濾波去噪，然后Radon變換重構(gòu)模糊核函數(shù)[8].文獻(xiàn)[9]以多組代價函數(shù)的優(yōu)化為目標(biāo)，實現(xiàn)圖像去模糊效果且極大的提高了算法的實時性.與此同時基于傳統(tǒng)的參數(shù)估計的圖像盲去模糊算法也被提出[10，11].

這些算法的計算量十分巨大，且對計算機(jī)內(nèi)存要求很高.本文提出了一種新的交替快速迭代算法，能極大的提高圖像去模糊的速度，使得圖像盲去模糊算法的實時實現(xiàn)成為可能.

2 原理與模型

2.1 圖像模糊方程

一般來說，運動模糊圖像的數(shù)學(xué)模型為：

g(x，y)=f(x，y)?h(x，y)+n(x，y)

(1)

式中g(shù)(x，y)表示降晰后的模糊圖像，f(x，y)表示降晰前的清晰圖像，h(x，y)為模糊核函數(shù)(PSF)，n(x，y)一般假定為零均值的高斯噪聲.

2.2 Radon變換

對一個二維圖像f(x，y)，其Radon變換可以看成圖像f(x，y)沿與某一直線的一維投影，或者說是圖像f(x，y)定義在與某一直線的線積分.所以f(x，y)的Radon變換數(shù)學(xué)表達(dá)式為[12]：

(2)

式中的直線l的方程為：b-xcosθ-ysinθ=0，利用沖激函數(shù)δ(·)的特性，Radon變換的數(shù)學(xué)表達(dá)式可以寫成如下表達(dá)式：

(3)

式中b表示直線l到坐標(biāo)原點的距離.

3 模糊核的估計

經(jīng)驗告訴我們，如果一幅自然圖像中有直線輪廓，模糊后的結(jié)果是直線輪廓邊緣會發(fā)生明顯的變化，產(chǎn)生毛邊，即邊緣的模糊，而這些邊緣的模糊與模糊核函數(shù)是有直接關(guān)系，這里通過利用圖像邊緣的模糊方程來推導(dǎo)其與模糊核函數(shù)之間的關(guān)系.

3.1 模糊核與模糊直線邊緣關(guān)系

前面已經(jīng)定義了一個二維圖像f(x，y)的Radon變換，如式(3)所示，根據(jù)二維函數(shù)的卷積定義，不考慮噪聲的影響，模糊圖像的卷積方程寫成積分方程的形式為：

(4)

假設(shè)二維圖像f(x，y)是一個理想的直線，其與x軸夾角為θ′，且θ′=θ+π/2，直線方程可以寫成如下表達(dá)式：

f(x，y)=δ(xsinθ′-ycosθ′)

(5)

又因為f(x1-x，y1-y)是f(x，y)的線性變換，所以也是直線，用沖激函數(shù)表達(dá)此直線為：

(6)

對于一個理想直線圖像，模糊后所得到的模糊圖像假設(shè)為gL(x，y)，其表達(dá)式可以重新寫為：

(7)

代入θ′=θ+π/2，有：

(8)

根據(jù)Radon變換定義，有：

(9)

gL(x1，y1)=Rh(ρxy，θ)

(10)

由此可見，當(dāng)f(x，y)是一個理想直線輪廓時，模糊后的直線法向截面輪廓等于模糊核函數(shù)的Radon變換，這給我們提供了一個思路，如果找到一幅圖像中的模糊直線邊緣，就可以通過求解反Radon變換來確定模糊核函數(shù).

3.2 模糊核實驗測試

這里選用了文獻(xiàn)[5]的一個模糊核函數(shù)，如圖1所示，其與一個只包含一條理想直線輪廓(不同角度)的黑白圖像互相卷積來進(jìn)行實驗，對比結(jié)果如圖2所示.

圖1 測試用模糊核函數(shù)Fig.1 Kernel of testing

圖2 測試結(jié)果Fig.2 Result of testing

3.3 模糊核的估計與圖像的恢復(fù)

為簡化起見，這里分別用f，g，h表示f(x，y)，g(x，y)，h(x，y)，根據(jù)圖像的貝葉斯恢復(fù)理論有：p(f，h|g)∝p(g|f，h)p(f)p(h)，對每一個像素來說，模糊圖像的梯度值和清晰圖像的梯度值是相近的，這里假定其差值服從均值為零，方差為σ2的高斯分布[13]，故有：

(11)

考慮已知f，h后，一方面，模糊圖像g的分布由噪聲n分布規(guī)律來決定，一般假定噪聲為零均值的高斯分布，用pg表示，另一方面，為引入正則項，限制解的范圍，引入一個新的條件，前面已經(jīng)推導(dǎo)模糊圖像的邊界與模糊核函數(shù)的Radon變換相等，這里定義為：

p(g|f，h)∝pg(g|f，h)*pe(g|f，h)

(12)

假定圖像噪聲為零均值的高斯噪聲，有：

(13)

(14)

定義代價函數(shù)：

E(f，h)=-log(p(f，h|g))

(15)

為簡化計算，加快計算速度，取p(h)為高斯分布，有：

E(f，h)=λ1‖f?h-g‖2+λ2‖f-g‖2+
λ3‖h‖2+λ4∑q‖gLθq-Rθqh‖2

(16)

圖像的去模糊過程轉(zhuǎn)變?yōu)檫@個方程的優(yōu)化求解，交替優(yōu)化這個代價函數(shù)就能得到模糊核函數(shù)和清晰圖像.

優(yōu)化這類代價函數(shù)，采用交替迭代的算法，分為兩步，第一步，固定h，求解f：

(17)

式(17)包含兩個二次項，利用傅里葉變換快速求解：

(18)

第二步，固定f，求解h：

(19)

方程的三項都是二次項，但是第三項包含有Radon變換，采用了共軛梯度法進(jìn)行優(yōu)化求解：

(20)

式中粗體分別表示相應(yīng)的算子矩陣，可得方程的解為：

(21)

4 圖像塊的處理與圖像塊的選擇

4.1 引導(dǎo)濾波

圖像濾波是為了去噪，由于本算法利用直線邊緣來估計模糊核，即模糊核函數(shù)只與圖像直線邊緣附近像素的變化有關(guān)，因此希望直線邊緣附近圖像在去噪同時又能保持圖像邊緣.一些簡單的去噪方法會導(dǎo)致圖像直線邊緣發(fā)生變異，變得更加模糊，破壞了直線邊緣所攜帶的原始信息，因此，需要選擇合適的濾波算子，在實現(xiàn)圖像去噪的同時，進(jìn)一步平滑光滑區(qū)域部分的圖像，保持直線邊緣附近的圖像信息，這里采用引導(dǎo)濾波來去噪.

引導(dǎo)濾波是近些年提出的圖像濾波方法[14]，其核心思想是假設(shè)該輸出圖像與引導(dǎo)圖像(或是輸入圖像)在一個二維圖像窗口wk內(nèi)滿足線性關(guān)系.

qi=akIi+bk，?i∈wk

(22)

式中，qi是輸出圖像，Ii是輸入圖像，i和k是索引，在窗口wk內(nèi)，一般假定ak和bk是常數(shù).對上式兩邊取梯度有：qi=akIi，這表明其梯度特性沒有發(fā)生變化，邊緣細(xì)節(jié)保持不變.設(shè)輸入圖像pi，通過讓輸入圖像和輸出濾波圖像的差值最小來確定兩個參數(shù)akbk的值，如下式所示：

(23)

(24)

直接采用輸入濾波圖像作為引導(dǎo)圖像的時候，即Ii=pi，如果=0，則是式(6)的解，輸出等于輸入，引導(dǎo)濾波器無效.如果>0，在圖像平滑的區(qū)域，ak的值較小或等于0，而bk的值近似于等于此時引導(dǎo)濾波器為加權(quán)均值濾波器；而在圖像邊緣特征比較豐富的圖像區(qū)域，ak近似于1，bk近似于0，引導(dǎo)濾波器也基本無效，從而保持原來圖像邊緣的基本特性.故引導(dǎo)濾波器的濾波效果在窗口尺寸固定的情況下，基本是由參數(shù)確定.引導(dǎo)濾波另外一個更重要的特性是其計算時間復(fù)雜度，引導(dǎo)濾波的時間復(fù)雜度為O(No)(No為待處理圖像的像素數(shù)目).圖3為對LENA圖像加噪后，分別用雙邊濾波和引導(dǎo)濾波去噪后效果圖像.

圖3 圖像去噪Fig.3 Image denoising

4.2 圖像塊的選擇與提取

實驗中發(fā)現(xiàn)，由于不同圖像的尺寸大小的差異，以及不同圖像所包含的直線輪廓數(shù)量也各不相同，故而算法的計算時間和內(nèi)存消耗相差懸殊，特別是在圖像尺寸較大的情況下，計算機(jī)內(nèi)存不足和計算時間過長的問題特別突出.我們已經(jīng)知道圖像的直線輪廓像素數(shù)目在整幅圖像中是一個較小的值，這啟發(fā)我們可以通過提取圖像主要輪廓邊界，降低圖像尺寸來完成模糊核的估計，實現(xiàn)模糊圖像的恢復(fù).

圖像的模糊核函數(shù)依賴于圖像中直線輪廓法向截面輪廓的反Radon變換，而反Radon變換的知識告訴我們，獲得在不同角度下圖像Radon變換的向量越多，提供的信息量就越大，使用最大后驗概率的獲得精確解的概率就越大.

因此，本文首先對圖像進(jìn)行分塊，然后找到包含不同角度最多直線邊界的幾個塊進(jìn)行拼合(本文選用4個)，拼合時候要考慮模糊核尺寸的大小，利用掩膜消除塊邊界的影響，然后利用拼合的圖像完成模糊核函數(shù)的重構(gòu)，從而降低計算量和減少內(nèi)存的消耗.圖4為對一幅圖像采用本方法先分塊，然后選擇合適的圖像塊，最后進(jìn)行拼合后的結(jié)果.

圖4 圖像塊選擇Fig.4 Image patch choosing

這里采用了文獻(xiàn)[6]方法檢測邊界，因為直線輪廓的法向截面輪廓與模糊核函數(shù)的Radon變換相等，需要對提取的直線輪廓進(jìn)行篩選，選擇在直線邊界附近沒有其他輪廓干擾的直線邊界作為候選邊界.

5 圖像恢復(fù)

5.1 合成模糊圖像的實驗結(jié)果

這里選用三幅自然圖像作為測試圖像，如圖5所示，從文獻(xiàn)[5]選擇三個模糊核函數(shù)作測試，如圖6所示，其相互卷積加2%的高斯白噪聲，得到9幅含噪的模糊圖像，分別采用本文的算法與文獻(xiàn)[6]的最大后驗概率法(MAP Radon)，進(jìn)行對比試驗(計算機(jī)配置CPU為3.1GHz，內(nèi)存為10GB，軟件MATLAB2012a)，實驗的結(jié)果如表1所示.從表可知，本文算法的計算時間和PSNR值均優(yōu)于MAP Radon算法，同時，本文算法克服了MAP Radon算法占用內(nèi)存大的缺點.

圖5 測試用圖像Fig.5 Image of testing

圖6 測試用模糊核函數(shù)Fig.6 Kernel of testing

表1 與算法比較PSNR值(db)和計算時間(s)Table 1 Comparison of PSNR values and the computation times

圖7 去模糊結(jié)果Fig.7 Image deblurring results

5.2 實際自然模糊圖像恢復(fù)

分別選用文獻(xiàn)[6]提供的測試圖像和文獻(xiàn)[15]提供的一副測試圖像進(jìn)行實驗，實驗的結(jié)果如圖7所示，計算時間如表2所示，從表中可見，本文算法的時間是遠(yuǎn)遠(yuǎn)小于MAPRadon算法，且在運行過程中，不會出現(xiàn)該算法最常見的內(nèi)存不夠的現(xiàn)象，在圖像的恢復(fù)質(zhì)量上，本文算法恢復(fù)的圖像，在視覺上的恢復(fù)效果絲毫不遜色于對比算法所獲得的恢復(fù)效果，在一些細(xì)節(jié)上，本文算法更具有優(yōu)勢.

表2 計算時間對比Table 2 Comparison of the computation time(s)

6 結(jié) 論

本文提出一種基于圖像信息塊的單幅模糊圖像盲恢復(fù)算法，通過檢測圖像的直線邊界，在圖像中找出包含最多信息邊界的圖像塊進(jìn)行拼接，以此拼接圖像來估計模糊核函數(shù)，完成模糊圖像的盲恢復(fù).實驗結(jié)果證明，相對于MAPRadon算法，本文算法在不損失恢復(fù)圖像質(zhì)量的前提下，明顯提高了計算速度并降低了內(nèi)存的占用量，降低了計算的時間和對硬件的要求.