郭立文,廖永忠
1(陜西國防工業(yè)職業(yè)技術(shù)學院 計算機與軟件學院,西安 710300) 2(湖南第一師范學院 信息科學與工程學院,長沙 410205)
成像過程中的相機晃動是圖像模糊的一個重要成因,去除這種運動模糊是圖像去模糊研究的一個重要研究方向,圖像去模糊過程本質(zhì)是對一類病態(tài)的反卷積方程的求解,如何高效、準確的求解此類方程是圖像去模糊研究的難點問題.
圖像的盲去模糊是指在圖像的模糊核函數(shù)未知情況下的圖像去模糊方法,盲去模糊算法是當前圖像恢復研究十分活躍的一個領域,一些學者在近幾年取得很多高水平的研究成果,相關的期刊和會議上有很多此類研究的文章[1-4].MIT的Levin提出一種基于邊緣概率密度函數(shù)的模糊核估計算法[5],用高斯混合分布模型表示自然圖像的先驗分布,然后采用期望最大化(EM)算法交替迭代實現(xiàn)模糊圖像的恢復.Taeg Sang Cho在文獻[6]中利用圖像直線邊緣在模糊過程中發(fā)生變化,其Radon變換與模糊核函數(shù)之間有確定的對應關系,從而確定模糊核函數(shù)的方法,然后用最大后驗概率(MAP)完成模糊圖像的恢復.Bae H提出一種快速的模糊圖像恢復算法[7],其認為由于自然圖像梯度的稀疏分布特性,圖像模糊信息都包含在圖像的邊界中,對圖像進行分塊,然后提取一些包含較多信息邊界的塊,再把它們進行拼接,利用這個模糊圖像塊完成模糊核的估計,該算法能降低計算量.Lin Zhong 提出強噪聲下大尺寸模糊核的圖像去模糊圖像算法,先利用方向濾波去噪,然后Radon變換重構(gòu)模糊核函數(shù)[8].文獻[9]以多組代價函數(shù)的優(yōu)化為目標,實現(xiàn)圖像去模糊效果且極大的提高了算法的實時性.與此同時基于傳統(tǒng)的參數(shù)估計的圖像盲去模糊算法也被提出[10,11].
這些算法的計算量十分巨大,且對計算機內(nèi)存要求很高.本文提出了一種新的交替快速迭代算法,能極大的提高圖像去模糊的速度,使得圖像盲去模糊算法的實時實現(xiàn)成為可能.
一般來說,運動模糊圖像的數(shù)學模型為:
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)一般假定為零均值的高斯噪聲.
對一個二維圖像f(x,y),其Radon變換可以看成圖像f(x,y)沿與某一直線的一維投影,或者說是圖像f(x,y)定義在與某一直線的線積分.所以f(x,y)的Radon變換數(shù)學表達式為[12]:
(2)
式中的直線l的方程為:b-xcosθ-ysinθ=0,利用沖激函數(shù)δ(·)的特性,Radon變換的數(shù)學表達式可以寫成如下表達式:
(3)
式中b表示直線l到坐標原點的距離.
經(jīng)驗告訴我們,如果一幅自然圖像中有直線輪廓,模糊后的結(jié)果是直線輪廓邊緣會發(fā)生明顯的變化,產(chǎn)生毛邊,即邊緣的模糊,而這些邊緣的模糊與模糊核函數(shù)是有直接關系,這里通過利用圖像邊緣的模糊方程來推導其與模糊核函數(shù)之間的關系.
前面已經(jīng)定義了一個二維圖像f(x,y)的Radon變換,如式(3)所示,根據(jù)二維函數(shù)的卷積定義,不考慮噪聲的影響,模糊圖像的卷積方程寫成積分方程的形式為:
(4)
假設二維圖像f(x,y)是一個理想的直線,其與x軸夾角為θ′,且θ′=θ+π/2,直線方程可以寫成如下表達式:
f(x,y)=δ(xsinθ′-ycosθ′)
(5)
又因為f(x1-x,y1-y)是f(x,y)的線性變換,所以也是直線,用沖激函數(shù)表達此直線為:
(6)
對于一個理想直線圖像,模糊后所得到的模糊圖像假設為gL(x,y),其表達式可以重新寫為:
(7)
代入θ′=θ+π/2,有:
(8)
根據(jù)Radon變換定義,有:
(9)
gL(x1,y1)=Rh(ρxy,θ)
(10)
由此可見,當f(x,y)是一個理想直線輪廓時,模糊后的直線法向截面輪廓等于模糊核函數(shù)的Radon變換,這給我們提供了一個思路,如果找到一幅圖像中的模糊直線邊緣,就可以通過求解反Radon變換來確定模糊核函數(shù).
這里選用了文獻[5]的一個模糊核函數(shù),如圖1所示,其與一個只包含一條理想直線輪廓(不同角度)的黑白圖像互相卷積來進行實驗,對比結(jié)果如圖2所示.
圖1 測試用模糊核函數(shù)Fig.1 Kernel of testing
圖2 測試結(jié)果Fig.2 Result of testing
為簡化起見,這里分別用f,g,h表示f(x,y),g(x,y),h(x,y),根據(jù)圖像的貝葉斯恢復理論有:p(f,h|g)∝p(g|f,h)p(f)p(h),對每一個像素來說,模糊圖像的梯度值和清晰圖像的梯度值是相近的,這里假定其差值服從均值為零,方差為σ2的高斯分布[13],故有:
(11)
考慮已知f,h后,一方面,模糊圖像g的分布由噪聲n分布規(guī)律來決定,一般假定噪聲為零均值的高斯分布,用pg表示,另一方面,為引入正則項,限制解的范圍,引入一個新的條件,前面已經(jīng)推導模糊圖像的邊界與模糊核函數(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變換,采用了共軛梯度法進行優(yōu)化求解:
(20)
式中粗體分別表示相應的算子矩陣,可得方程的解為:
(21)
圖像濾波是為了去噪,由于本算法利用直線邊緣來估計模糊核,即模糊核函數(shù)只與圖像直線邊緣附近像素的變化有關,因此希望直線邊緣附近圖像在去噪同時又能保持圖像邊緣.一些簡單的去噪方法會導致圖像直線邊緣發(fā)生變異,變得更加模糊,破壞了直線邊緣所攜帶的原始信息,因此,需要選擇合適的濾波算子,在實現(xiàn)圖像去噪的同時,進一步平滑光滑區(qū)域部分的圖像,保持直線邊緣附近的圖像信息,這里采用引導濾波來去噪.
引導濾波是近些年提出的圖像濾波方法[14],其核心思想是假設該輸出圖像與引導圖像(或是輸入圖像)在一個二維圖像窗口wk內(nèi)滿足線性關系.
qi=akIi+bk,?i∈wk
(22)
式中,qi是輸出圖像,Ii是輸入圖像,i和k是索引,在窗口wk內(nèi),一般假定ak和bk是常數(shù).對上式兩邊取梯度有:qi=akIi,這表明其梯度特性沒有發(fā)生變化,邊緣細節(jié)保持不變.設輸入圖像pi,通過讓輸入圖像和輸出濾波圖像的差值最小來確定兩個參數(shù)akbk的值,如下式所示:
(23)
(24)
直接采用輸入濾波圖像作為引導圖像的時候,即Ii=pi,如果=0,則是式(6)的解,輸出等于輸入,引導濾波器無效.如果>0,在圖像平滑的區(qū)域,ak的值較小或等于0,而bk的值近似于等于此時引導濾波器為加權(quán)均值濾波器;而在圖像邊緣特征比較豐富的圖像區(qū)域,ak近似于1,bk近似于0,引導濾波器也基本無效,從而保持原來圖像邊緣的基本特性.故引導濾波器的濾波效果在窗口尺寸固定的情況下,基本是由參數(shù)確定.引導濾波另外一個更重要的特性是其計算時間復雜度,引導濾波的時間復雜度為O(No)(No為待處理圖像的像素數(shù)目).圖3為對LENA圖像加噪后,分別用雙邊濾波和引導濾波去噪后效果圖像.
圖3 圖像去噪Fig.3 Image denoising
實驗中發(fā)現(xiàn),由于不同圖像的尺寸大小的差異,以及不同圖像所包含的直線輪廓數(shù)量也各不相同,故而算法的計算時間和內(nèi)存消耗相差懸殊,特別是在圖像尺寸較大的情況下,計算機內(nèi)存不足和計算時間過長的問題特別突出.我們已經(jīng)知道圖像的直線輪廓像素數(shù)目在整幅圖像中是一個較小的值,這啟發(fā)我們可以通過提取圖像主要輪廓邊界,降低圖像尺寸來完成模糊核的估計,實現(xiàn)模糊圖像的恢復.
圖像的模糊核函數(shù)依賴于圖像中直線輪廓法向截面輪廓的反Radon變換,而反Radon變換的知識告訴我們,獲得在不同角度下圖像Radon變換的向量越多,提供的信息量就越大,使用最大后驗概率的獲得精確解的概率就越大.
因此,本文首先對圖像進行分塊,然后找到包含不同角度最多直線邊界的幾個塊進行拼合(本文選用4個),拼合時候要考慮模糊核尺寸的大小,利用掩膜消除塊邊界的影響,然后利用拼合的圖像完成模糊核函數(shù)的重構(gòu),從而降低計算量和減少內(nèi)存的消耗.圖4為對一幅圖像采用本方法先分塊,然后選擇合適的圖像塊,最后進行拼合后的結(jié)果.
圖4 圖像塊選擇Fig.4 Image patch choosing
這里采用了文獻[6]方法檢測邊界,因為直線輪廓的法向截面輪廓與模糊核函數(shù)的Radon變換相等,需要對提取的直線輪廓進行篩選,選擇在直線邊界附近沒有其他輪廓干擾的直線邊界作為候選邊界.
這里選用三幅自然圖像作為測試圖像,如圖5所示,從文獻[5]選擇三個模糊核函數(shù)作測試,如圖6所示,其相互卷積加2%的高斯白噪聲,得到9幅含噪的模糊圖像,分別采用本文的算法與文獻[6]的最大后驗概率法(MAP Radon),進行對比試驗(計算機配置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
分別選用文獻[6]提供的測試圖像和文獻[15]提供的一副測試圖像進行實驗,實驗的結(jié)果如圖7所示,計算時間如表2所示,從表中可見,本文算法的時間是遠遠小于MAPRadon算法,且在運行過程中,不會出現(xiàn)該算法最常見的內(nèi)存不夠的現(xiàn)象,在圖像的恢復質(zhì)量上,本文算法恢復的圖像,在視覺上的恢復效果絲毫不遜色于對比算法所獲得的恢復效果,在一些細節(jié)上,本文算法更具有優(yōu)勢.
表2 計算時間對比Table 2 Comparison of the computation time(s)
本文提出一種基于圖像信息塊的單幅模糊圖像盲恢復算法,通過檢測圖像的直線邊界,在圖像中找出包含最多信息邊界的圖像塊進行拼接,以此拼接圖像來估計模糊核函數(shù),完成模糊圖像的盲恢復.實驗結(jié)果證明,相對于MAPRadon算法,本文算法在不損失恢復圖像質(zhì)量的前提下,明顯提高了計算速度并降低了內(nèi)存的占用量,降低了計算的時間和對硬件的要求.