基于循環(huán)和對稱邊界的圖像反卷積快速算法

徐夢娜，蔡 信，金文清，韓玉兵

（南京理工大學(xué) 電子工程與光電技術(shù)學(xué)院，江蘇 南京210094）

0 引 言

在現(xiàn)實生活中有很多因素會導(dǎo)致圖像降質(zhì)，如果把圖像降質(zhì)過程用一個卷積來描述，那么圖像復(fù)原就可以用反卷積來實現(xiàn)。但這通常是個反問題，反問題所共有的一個問題是病態(tài)性，需要正則化來修改最終的解。常用的正則化參數(shù)選取方法主要有經(jīng)驗選取方法、廣義交叉驗證方法GCV （generalize cross validation）、L－Curve方法，但這些方法通常需要圖像先驗知識或是計算量很大。近年來一種動態(tài)自適應(yīng)的參數(shù)選取方法應(yīng)用越來越廣泛，Eldar將這種方法推廣到圖像反卷積［1］中，R.Giryes在此基礎(chǔ)上又結(jié)合了小波的概念［2］。但在這些反卷積算法中多需要寫出具體的表達(dá)式，而大部分的算法是非線性的，無法寫出具體的表達(dá)式。我們的方法通過給出迭代表達(dá)式而不是具體算式規(guī)避了這些缺點，并且由于模糊矩陣的高度結(jié)構(gòu)化，通過邊界條件的選擇可以大大縮小時間效益［3－5］。

1 卷積快速算法

一般地，圖像的降質(zhì)模型可表示為

其中，F(xiàn) 為原圖像，G 為觀測到的降質(zhì)圖像，H 為降質(zhì)模型的確定性部分，一般假設(shè)為一個線性算子，代表圖像獲取過程中的各種扭曲、模糊及降采樣等，ξ 為加性噪聲部分。反卷積就是從觀測圖像G 中盡可能復(fù)原出F，這顯然是一個病態(tài)反問題，本文采用著名的total variation （TV）正則化以及Split Bregman方法［6］來求解，則有迭代算法

對于式 （2），采用文獻(xiàn) ［7］中的TV 去噪方法來求解，正則化參數(shù)λ的選取采用文獻(xiàn) ［8］中方法。TV 模型中，最重要的結(jié)構(gòu)，模糊HF 及其轉(zhuǎn)置操作HTG，在適當(dāng)?shù)倪吔鐥l件下可以被看成離散卷積，如何使這個結(jié)構(gòu)達(dá)到最好的效益就是我們所關(guān)心的［9］。一種有效的降低計算成本的方法就是設(shè)計快速算法來實現(xiàn)卷積操作，本文將主要討論周期邊界和對稱邊界條件下的兩種計算的快速算法。

1.1 基于循環(huán)邊界的DFT快速算法

循環(huán)邊界條件是指超出圖像邊界的像素數(shù)據(jù)是未超出部分的完整復(fù)制，基于這種邊界條件，模糊矩陣是一個塊循環(huán)矩陣，這些矩陣能被離散傅里葉矩陣對角化，因此通過傅里葉變換就很容易找到其逆矩陣［10］。為簡單起見，先從一維信號開始，假設(shè)一維原信號f’＝ （…，f－m＋1，…，f0，f1，…，fn，fn＋1，…，fn＋m，…）T，模糊函 數(shù)為h＝ （…，0，0，h－m，h－m＋1，…，h0，…，hm－1，hm，0，0，…）T，則卷積可表示為

由式 （3）可知，g 不僅僅由 （f1，…，fn）T決定，還由 （f－m＋1，…，f0）T和 （fn＋1，…，fn＋m）T決定。在加上邊界條件之前，我們先將式 （3）重新寫為［11］

其中

當(dāng)設(shè)定循環(huán)邊界條件，在式 （3）對所有j有fj＝fnj，則式 （4）可寫為

其中， （0｜Tl）和 （Tr｜0）是對Tl和Tr分別加上 （nm）列0后得到的n×n Toeplitz矩陣。加上循環(huán)邊界條件后最大的優(yōu)勢就是H 是一個循環(huán)矩陣，可以被離散傅里葉矩陣對角化，從而大大地減低運(yùn)算問題。令

由矩陣乘法可得Hwk＝λkwk，k＝0，1，…，m－1，再由矩陣?yán)碚摽芍?，λk是循環(huán)矩陣H 的特征根，而wk是該矩陣H 對應(yīng)于此特征根的特征向量，將m 個wk排成矩陣W＝ ［w0w1w2… wm－1］，就可以將H 表達(dá)成特征分解的形式H＝WΛW－1，式中Λ＝diag ｛λ0，λ1，…，λm－1｝為對角矩陣。利用以上的式子，容易得出

上述結(jié)論可以推廣到二維函數(shù)的分塊循環(huán)矩陣中［10，12］，與一維函數(shù)的循環(huán)矩陣對角化思路一樣，可使二維函數(shù)的分塊循環(huán)矩陣對角化。設(shè)H 是一個m1m2×m1m2塊循環(huán)矩陣，先將卷積函數(shù)h補(bǔ)零使尺寸擴(kuò)展到m1×m2，然后作離散傅里葉變換，記為B（μ，ν），易驗證H 有特征分解式H＝WΛW－1。其中W 是H 的特征向量排成的矩陣，是非奇異的，Λ 是對角矩陣，它是通過把B（μ，ν）的m1×m2個元素從第一行起，逐個順序填入m1m2×m1m2對角陣構(gòu)成。H 的特征矩陣W 和W－1之前的關(guān)系可以驗證得到

由此可采用離散傅里葉變換來處理，不需要實際地存放和處理大尺度的矩陣H。

根據(jù)以上理論，我們給出一個簡單的算例

采用MATLAB 內(nèi)置函 數(shù)imfilter（X，h，’circular’）以及DFT 算法得到卷積X＊h的結(jié)果如式 （6）所示

采用編寫的Matlab程序ImfilterTranspose（）函數(shù)以及DFT 算法得到HTX 的結(jié)果如式 （7）所示

由式 （6）、式 （7）可知，由卷積方法得到的結(jié)果與DFT 方法得到的結(jié)果一致，證明基于循環(huán)邊界的卷積可以由DFT 代替。

1.2 基于對稱邊界的DCT快速算法

對稱邊界條件是指超出圖像邊界的像素數(shù)據(jù)是未超出部分的鏡像對稱映射，基于這種邊界條件，模糊矩陣不再是一個塊循環(huán)矩陣，而是塊Toeplitz加Hankel矩陣 （block Toeplitz－plus－Hankel matrix），盡 管 這 種 矩 陣 的 結(jié) 構(gòu) 更 復(fù)雜，但也能被離散余弦矩陣對角化［4，11］，因此通過離散余弦變換也很容易找到其逆矩陣。

為簡單起見，也從一維反卷積開始，當(dāng)設(shè)定對稱邊界條件，在式 （3）中有

基于以上理論，則式 （4）可寫為

其中，J 是n×n逆矩陣。

首先定義C 為n×n離散余弦變換矩陣

其中，δi，j為克羅內(nèi)克函數(shù)，C 是正交的，則有CTC＝CCT＝I，同時，對于任何n維向量v，Cv和CTv 都能通過DCT 的實部操作來得到。定義Γ＝ ｛CTΛC｜Λ 為n×n對稱陣｝是空間包含所有能被C 對角化的矩陣，令Q＝CTΛC∈Γ，對CQ＝ΛC 兩邊同乘e1＝ （1，0，…，0）T，可得到Q 的特征值

由式 （9），Q 的特征值能夠通過對Q 的第一列做一次DCT 得到，事實上，在Γ 里的所有矩陣都是由其第一列決定的。再定義向量v＝ （v0，v1，…，vn－1）T的轉(zhuǎn)移為σ（v）≡ （v1，v2，…，vn－1）T，T（v）是 以v 為 第 一 列 的n×n對稱Toeplitz矩陣，H（x，y）是以x 為第一列向量，以y為最后一列向量的n×n Hankel矩陣［11］。下面再給出幾個定理：

引理1［13］設(shè)Γ 是所有能被離散余弦變換C 對角化的矩陣集，則

從引理1可以看出能被C 對角化的矩陣是一些特殊的Toeplitz加Hankel矩陣。

定理1［11］模糊函數(shù)h 是對稱的，則對所有i有hi＝h－i，則式 （8）中H 能被重新寫為

其中，u＝ （h0，h1，…，0，…，0）T。

由定理1就能得到式 （8）中的解f 由f＝CΛ－1CTg 得到，其中Λ 是由H 特征值組成的對角矩陣，可由一次DCT變換得到。因此f 能由3次DCT 變換得到。

一維的結(jié)果可以很容易地拓展到二維，在對稱邊界條件下，H 是塊Toeplitz和Hankel矩陣。

定理2［11］假如模糊函數(shù)hi，j是對稱的，則對所有i，j有hi，j＝hi，－j＝h－i，j＝h－i，－j，則H 能被二維離散余弦變換矩陣CC 對角化，其中為張量積。

令P 為置換 矩 陣 滿 足 ［PTΛP］i，j；k，l＝ ［Λ］k，l；i，j，1≤i，j≤n，1≤k，l≤n，其中Λ 中第 （k，l）個塊中的第 （i，j）個數(shù)被第 （i，j）個塊中的第 （k，l）個數(shù)置換，則有PT（C I）P＝（IC）及PTΛP ＝，是由塊矩陣構(gòu)成的對角陣，每個都和式 （8）中的H 有相同的形式，則對所有j有C＝，是一個對角陣，則有

根據(jù)以上理論，基于MATLAB7.9.0.529版本，我們也給出一個簡單的算例：

采用MATLAB內(nèi)置函數(shù)imfilter（X，h，'symmetric'）以及DCT 算法得到卷積X＊h的結(jié)果如式 （11）所示

采用編寫的Matlab程序ImfilterTranspose（）函數(shù)以及DCT 算法計算HTX 得到的結(jié)果如式 （12）所示

由式 （11）、式 （12）可知，卷積方法得到的結(jié)果與DCT 方法得到的結(jié)果一致，證明基于對稱邊界的卷積可以由DCT 代替。

2 實驗研究

將以上算法應(yīng)用在TV 正則化的圖像反卷積中，模糊算子采用標(biāo)準(zhǔn)差0.8的5×5高斯模糊核，噪聲方差為0.05（灰度值歸一化為 ［0，1］），對圖像進(jìn)行反卷積重建，實驗平臺基于Matlab7.9.0.525 版本，CPU 速度為2GHZ，內(nèi)存為1GHZ 的32 位操作系統(tǒng)。選用大小為256×256 的Cameraman及Peppers圖像，給定循環(huán)邊界條件，在迭代算法中應(yīng)用DFT 快速算法，結(jié)果分別如圖1及圖2所示。選用大小為256×256的Boats及Elaine圖像，給定對稱邊界的條件，在迭代算法中應(yīng)用DCT 快速算法，結(jié)果分別如圖3及圖4所示。

由圖1－圖4可以看出，反卷積都取得了比較好的效果，較大程度地恢復(fù)出了原圖像。在此基礎(chǔ)上，我們再來比較這兩種快速算法的時間效益以及圖像的信噪比提高情況，信噪比計算公式采用下式

結(jié)果如表1，表2所示，其中計算時間是對應(yīng)一個正則化參數(shù)的完整的反卷積過程。

表1 卷積方法與DFT 方法對比

表2 卷積方法與DCT 方法對比

由表1，表2可以看出，不論是基于循環(huán)邊界條件還是對稱邊界條件，圖像輸出信噪比得到明顯提高，基于這兩種邊界的快速算法在計算中也節(jié)省了大量的時間。

3 結(jié)束語

在圖像反卷積算法中，采用正則化的方法來尋找解，但由于模糊矩陣十分龐大，在實際計算時不易顯式得到和存儲，尤其是在計算矩陣轉(zhuǎn)置的時候，很難實現(xiàn)，但是在多數(shù)迭代算法中又需要用到。為了減少運(yùn)算時間，根據(jù)圖像卷積時特殊邊界的特性，對基于循環(huán)邊界條件和對稱邊界條件的圖像，在需要計算模糊矩陣及其轉(zhuǎn)置矩陣的時候，可以分別用DFT 及DCT 快速算法來實現(xiàn)，既能節(jié)約計算時間，又能準(zhǔn)確地保證圖像復(fù)原質(zhì)量。實驗結(jié)果已經(jīng)表明該算法適用不同的圖像并都能取得良好的效果。

［1］Eldar Y.Generalized SURE for exponential families：Applications to regularization ［J］.IEEE Trans Signal Process，2009，57 （2）：471－481.

［2］Giryes R，Elad M，Eldar Y C.The projected GSURE for automatic parameter tuning in iterative shrinkage methods［J］.Appl Comput Harmon Anal，2010，3 （30）：407－422.

［3］Shi Y，Chang Q.Acceleration methods for image restoration problem with different boundary conditions ［J］.Appl Numer Math，2008，58 （5）：606－614.

［4］Lv X G，Song Y Z，Wang S X.Image restoration with a highorder total variation minimization method ［J］.Applied Mathematical Modelling，2013，37 （16－17）：8210－8224.

［5］Ji H，Wang K.Robust image deblurring with an inaccurate blur kernel［J］.IEEE Trans Image Process，2012，21 （4）：1624－1634.

［6］Goldstein T，Osher S.The split bregman method for L1－regularized problems［J］.SLAM Journal on Imaging Science，2009，2 （2）：323－343.

［7］Getreuer P.Rudin－osher＿fatemi total variation denoising using split bregman ［OL／J］.Image Processing on Line，http：／／dx.doi.org／10.5201／ipol.2012.g－tvd，2012.

［8］Ramini S，Blu T，Unser M.Monte－carlo sure：A black－box optimization of regularization parameters for general denoising algorithms［J］.IEEE Trans Image Process，2008，17 （9）：1540－1554.

［9］Wang Yulun，Yang Junfeng，Yin Wotao，et al.A new alternating minimization algorithm for total variation image reconstruction ［J］.SIAM J Imag Sci，2008，1 （3）：248－272.

［10］Gonzalez R，Woods R.Digital image processing ［M］.New York：Addison Wesley，1992.

［11］Michael K，Raymond H，Wun Cheung.A fast algorithms for deblurring models with neumann boundary conditions［J］.SIAM J Sci ComPut，1999，21 （3）：851－866.

［12］MIAO Qing.The research and application of regularization in image restoration ［D］.Changsha：National Defense Science and Technology University，2005：17－19 （in Chinese）.［苗晴.圖像復(fù)原中正則化方法的研究及應(yīng)用 ［D］.長沙：國防科學(xué)技術(shù)大學(xué)，2005：17－19.］

［13］Chan R，Chan T，Wong C.Cosine transform based preconditioners for total variation deblurring ［J］.IEEE Trans Image Processing，1999，8 （10）：1472－1478.

［14］Getreuer P.Total variation deconvolution using split bregman［J／OL］.Image Processing on Line，http：／／dx.doi.org／10.5201／ipol.2012.g－tvd，2012.