基于小波域稀疏最优的图像修复方法

王　斌,胡辽林,曹京京,薛瑞洋,刘光飞

(西安理工大学机械与精密仪器工程学院,陕西西安 710048)



基于小波域稀疏最优的图像修复方法

王斌,胡辽林,曹京京,薛瑞洋,刘光飞

(西安理工大学机械与精密仪器工程学院,陕西西安 710048)

摘要：由模糊和噪声引起的图像退化属于非线性病态逆问题,修复比较困难.由于小波的稀疏表示能力较强,为提高修复质量,提出利用正交小波作为稀疏基,以小波系数的稀疏性为先验构造凸函数,最小化后得到修复图像;并提出将优化问题转化为逼近算子形式,利用不动点理论求解;证明了只需对构造出来的迭代形式的解析解反复迭代就可以得到最优解.对方法的构造过程、收敛性和复杂度进行了细致的分析,给出了迭代解,并结合加速方法提高了算法速度.仿真表明,本文方法具有较强的修复能力,收敛速度较快,能够有效去除模糊和噪声,保留图像的边缘和细节信息.

关键词:模糊;小波;稀疏性;凸函数;加速

1引言

图像在拍摄、传输和储存中会降低质量,图像修复很重要.图像降质的主要原因有模糊、噪声、像素丢失、有损压缩等.造成模糊的主要原因包括相机振动(运动模糊)、对焦不准(高斯模糊)、目标高速运动(运动模糊)、大气湍流(湍流模糊)、人为因素(均值模糊).目前,去模糊方法多种多样,大致可以分为四大类:空域方法、概率类方法、变换域方法、空域与变换域相结合方法.

空域方法中常用的是三维块匹配(BM3D,Block Matching 3D),首先把图像分块,利用块之间的相似性,把结构相似的块组合起来形成三维数组,再用联合滤波的方法处理三维数组,最后通过逆变换,把结果返回到原图像中,得到修复图像.该方法峰值信噪比较高,视觉效果也很好,但是,用三维模型处理二维问题,且中间含有图像分块、相似性检测、联合滤波等高复杂度运算,因此总体复杂度很高,计算较慢.

概率类方法以贝叶斯推断和最大后验概率(MAP,Maximum A Posterior)为基础,通过确立退化图像和原始图像的概率分布来恢复图像,对退化过程的知识要求少,鲁棒性强.但是由于MAP需要较强的先验概率信息作为推断基础,因此其适用范围有限;贝叶斯推断要假设图像服从高斯分布,往往导致修复结果过于平滑,无法保持边缘信息.此类方法中较有代表性的是Richardson-Lucy(RL)方法、最小均方差法(MMSE,Minimum Mean Square Error)、变分贝叶斯法,前两种方法本质上是等价的.

常用的变换域方法包括:小波与超小波域、训练字典域等.小波与超小波域方法以变换域系数的分布特点为先验,对系数进行操作[1,2],此类方法具有强大的多分辨特点,但系数操作很复杂且存在“振铃”现象,需要引入消除“振铃”的辅助措施.字典学习是预先学习一个具有强稀疏表示能力的字典,寻找待修复图像在字典域的最稀疏表示,从而去除模糊[3],其优点在于学习得到的字典能够克服小波的一些痼疾,缺点是准备工作比较难,且目标函数通常为非凸函数,只有局部解,没有全局解,对参数和初始条件的依赖较强;此外,由于此类方法利用图像分块的思想确立像素点与其临近点的信息关联,当图像块重新拼接成完整图像时存在信息不连贯、重叠等问题.

空域与变换域相结合方法的基本思路是寻找图像在空域或变换域的先验信息,利用Tikhonov正则化方法建立能量泛函模型,然后利用优化方法求解模型,得到修复后的图像.常用的空域先验为全变差域的稀疏性,常用的变换域先验为图像在变换域的稀疏性.空域先验常用的模型包括TV-l2模型[4],TV-l0模型[5],TV-mix norm模型[6],此类方法最大的优点是可以有效保留边缘细节信息.但是,由于模糊核与全变差范数很难分离,需引入特殊的分裂方法,如分裂-Bregman方法[7,8]、对偶范数分裂法[4]及Fermat’s rule分裂法[9],而且此类方法还会导致过度光滑.变换域先验常用的模型包括小波稀疏模型[10～12],Contourlet模型,Curvelet模型.文献[12]提出分裂增广拉格朗日收缩法求解小波稀疏模型,可以推广到全变差模型.

针对各种方法存在的问题,本文提出了利用正交小波作为稀疏基,以小波系数的稀疏性作为先验约束项,噪声图像和当前值的差的l2范数作为保真项,构造凸函数,最小化此凸函数就得到修复后的图像;并将凸函数优化问题转化为逼近算子形式,利用不动点理论来求解;对方法的构造过程、收敛性和复杂度进行了细致的分析,给出了迭代解,并结合加速方法提高了算法速度[23].相比于其他方法,本文方法的突出优点是:迭代流程清晰简单,对参数的依赖性小,容错性强,收敛速度快.仿真表明,本文方法具有较强的修复能力,收敛速度较快,能够保留图像的边缘和细节信息.此外,还比较了小波基和离散余弦变换(DCT)基的优劣,并比较了其他方法.小波稀疏模型发源于稀疏表示,本文是对小波稀疏模型的深化和发展.由于小波对二维图像的稀疏表示能力有限,修复质量不是理论上的最优,下一步将结合更强的稀疏表示工具,构造效果更好的方法.

2图像修复模型及求解

2.1修复模型

图像降质模型可表述为

d*x+w=b

(1)

x为n×n的无噪声图像,d为模糊核,w和b为噪声和退化图像,*表示卷积.每种模糊核都对应一个Toeplitz矩阵,式(1)可表示为

Dx+w=b

(2)

D为d对应的Toeplitz矩阵.但Toeplitz矩阵太大,难以直接处理,故选择文献[24]中的方法构造与图像大小相同的模糊矩阵D.

正常图像经小波变换后,非零系数较少,而受到模糊和噪声干扰的图像则不是这样.通过约束小波域的稀疏性,可得到清晰图像,即求解l0范数最小化问题,但无法直接求解.文献[13]中,Donoho证明了在字典冗余度较低的情况下,l0与l1范数最小化问题是等价的.本文采用正交小波作为字典,符合以上条件.所以图像修复的数学表达式为

(3)

W()表示正交小波变换,WT( )表示反变换,F(x)为目标函数.上式的解x*为

(4)

设xw是x对应的小波系数,则问题(3)归结为

(5)

易知式(3)和(5)两种形式是等价的,第一项都是保真项,第二项为约束项,l1范数用于约束小波域的稀疏度.λ为正则化常数,表示‖W(x)‖1或‖xw‖1所占权重,多采用人工调试的办法选取λ,第一项都为可微凸函数,第二项都为不可微凸函数,无法直接求解.下面只针对式(5)求解.

2.2问题转化

(6)

(7)

利用函数微分的性质,可推导出凸函数显然还满足

f(y)+〈▽f(y),x-y〉≤f(x)≤f(y)+〈▽f(y),

(8)

(9)

(10)

+‖xw‖1

(11)

定义逼近算子为

(12)

令P(xw)=xw-2(Lf)-1W(DT(DWT(xw)-b)),显然式(11)可写为

proxλ(Lf)-1‖xw‖1(P(xw))

(13)

将式(13)写成迭代形式即

(14)

其中xw,k表示第xw经第k次迭代得到的值.

2.3问题求解

正常的下半连续凸函数的次微分为对任意x,y∈domf都满足

f(x)-f(y)-〈p,x-y〉≥0

(15)

0∈∂(f(y)+g(y))∈▽f(y)+∂g(y)

(16)

引理1稳定非扩张映射proxλ(Lf)-1‖xw‖1(P(xw))的不动点也是令F(xw)取最小值的x值.

(17)

两边同乘ω变为

(18)

移项可得

(19)

由于proxωg(·)=(I+ω∂g)-1(·)[14],则式(19)等价于

(20)

剩下的问题就是如何求解xw,k=proxλ(Lf)-1‖xw‖1(P(xw,k-1)).事实上利用shrinkage-operator[15]可得

(21)

3算法流程及复杂度分析

3.1算法流程

标准逼近迭代算法流程如算法1所示.

3.2加速

算法1虽然可以解决图像修复问题,但是收敛速度较慢,处理较大尺寸的图像时比较困难.因此有必要采用加速策略,以提高收敛速度.本文主要应用文献[15]中的加速方法,如算法2所示.该方法拥有O(1/k2)的总体收敛速度,接近一阶算法的收敛速度极限,而且该方法鲁棒性强,迭代过程较简单.

3.3复杂度分析

算法1包括三大部分:初始化,迭代计算xw,m共k次,第k+1步做小波逆变换.对于式(21),根据文献[15],其时间复杂度一般认为是O(1/(k+1)),而本算法中的shrinkage-operator包含大量内部运算,主要是小波变换、逆变换各一次和两次矩阵相乘.设每次小波变换时间复杂度为O(n2),矩阵相乘的时间复杂度为O(n3),则本算法总时间复杂度为O(2n2(n+1)/(k+1)).算法2的收敛速度为O(1/(k+1)2).由于算法2的每次迭代也包含小波变换、逆变换各一次和两次矩阵相乘,故其总时间复杂度为O(2n2(n+1)/(k+1)2).

4仿真及分析

仿真采用256×256标准PGM格式的cameraman图像,稀疏基采用2维正交小波变换(Daubechie小波和Symlet小波),采用算法2.假设图像先被模糊干扰,后被噪声干扰.实际上,若图像先被噪声干扰,后被模糊干扰,也不会影响结论,因为大多数模糊核,都对图像有平滑作用.若图像先被噪声干扰,然后被模糊干扰,那么模糊会减轻噪声的干扰,所以处理前一种情况的难度更大一些.首先,添加方差为0.001的弱高斯噪声和掩膜为5×5、标准差Sigma=4的高斯模糊,如图1所示.小波基为Haar小波,分解层数为2层,经过100次迭代,峰值信噪比(PSNR)为32.47dB,从图1(d)可以看出,人脸、相机、手套等细节恢复得很好,基本没有模糊残留.若噪声方差增加到0.1,则修复图像的PSNR为29.97dB.

然后,添加方差为0.001的弱高斯噪声和掩膜为9×9、标准差Sigma=4的强高斯模糊,如图2所示.显然,随着模糊的增强,修复的质量有所下降,PSNR为27.35dB.但本文算法依然能够较好地修复图像,并较好地保持局部细节.若噪声方差增加到0.1,其他条件不变,则修复图像的PSNR为26.21dB.

接下来,验证离焦模糊.对原始图像添加方差为0.001的高斯噪声和掩膜为5×5、标准差Sigma=4的离焦模糊,如图3所示,修复后图像的PSNR为33.13dB.从图3的(c)和(d)可以看出,无论是整体轮廓还是人脸、相机、手套等细节都恢复得很好,修复质量较好.若噪声方差增加到0.1,其他条件不变,则修复图像的PSNR为30.14dB.再对原始图像添加方差为0.001的高斯噪声和掩膜为9×9、标准差Sigma=4的离焦模糊,如图4所示,修复后图像的PSNR为31.97dB,说明本文算法在去除离焦模糊时具有较强的稳定性.若噪声方差增加到0.1,其他条件不变,则修复图像的PSNR为29.88dB.由于图像尺寸为256×256,选择的掩膜尺寸都是奇数,这是因为高斯模糊和离焦模糊的点扩散函数是高斯型曲线,选择奇数尺寸比较合适.

图5为不同掩膜尺寸下、不同类型模糊修复后的PSNR.可以看出,当掩膜尺寸较小时,四种模糊的修复效果相差无几.随着掩膜尺寸的增大,四种模糊中运动模糊的PSNR下降最平缓,修复效果最好;另外三种相差不大,高斯模糊修复效果最差.

图6为不同稀疏基的修复效果比较,比较对象为小波和DCT.可以看出,无论是处理高斯模糊还是离焦模糊,小波都要优于DCT,这是因为小波的稀疏能力强于DCT.

本文提出的方法可以适用于高斯噪声和模糊同时存在的情况,即便噪声强度很大,也可以有效地修复图像.之所以选取较小的噪声方差,主要是因为我们更侧重于去模糊,因为去模糊的难度更大,只要把模糊核去掉,本文方法就变成了一个去噪方法.

表1为使用不同小波基的PSNR与结构相似度(SSIM)比较,其中SSIM值越大表示效果越好,其最大值为1.统一添加5×5的高斯模糊.表1中Daub4和Sym4表示4层分解.一般来说,消失矩越高的小波效果越好,但事实并非如此.因为小波对数据的稀疏能力,不仅与消失矩有关,也与小波基的平滑性及待处理数据与小波基的相似性紧密相关.例如Haar小波(即1阶Daub小波),虽然平滑性欠佳,但是由结构上的相似性,它对指纹图像进行稀疏表示时就强于其它小波.而且,消失矩越高,支撑长度越大,会影响稀疏表示能力,所以表1中数据的起伏是正常的.

表1　不同小波基的PSNR(dB)/SSIM比较

由于图像修复方法众多,理论、速度、预设条件相差很多,选择几种新方法做比较,比较对象分别为GPRS-wavelet[16]、l0-abs[17]、PLMMSE[18]、ADMDU-DEB[3]、BM3D五种算法.其中,GPRS-wavelet和l0-abs属于空域与变换域相结合的算法,是该方向比较有代表性的算法;PLMMSE属于贝叶斯方法;ADMDU-DEB属于字典学习类方法;BM3D属于空域方法.这五个对比对象属于较有代表性的算法,具有一定的说服力.表2中纵向为不同图像,分别为标准格式的cameraman、pepper、barbara和lena,横向为不同算法,统一添加9×9的高斯模糊后得到的PSNR/SSIM.可以看出,本文算法与ADMDU-DEB和BM3D较为接近,明显强于其它三种算法.在具体指标上,本文算法与ADMDU-DEB和BM3D互有高低,但差距很小,说明修复能力相当.但是,ADMDU-DEB和BM3D复杂度很高,ADMDU-DEB需要预先采集大量样本图片训练字典,训练过程较复杂,参数选取也需要相当的技巧;而BM3D以三维模型来处理二维问题,算法复杂度很高.本文算法对比ADMDU-DEB的优势在于简易性,因为本文算法只需要小波基即可.在计算方面,推导过程看似繁琐,但迭代过程很简单,只包含基本运算如矩阵相乘和相加,且迭代时没有图像分块(ADMDU-DEB)、相似性检测(BM3D)等高复杂度运算,因此在简单性上优于ADMDU-DEB和BM3D.

表2　不同算法的PSNR(dB)/SSIM比较

5结论

为解决图像修复问题,提出了利用优化后的正交小波作为稀疏基,以小波变换系数的稀疏性作为先验信息,再选择当前值与原始值之差的l2范数作为保真项,从而构造出一个凸函数,最小化此凸函数可以得到修复后的图像.此凸函数具有proximal 算子的形式,可以利用不动点理论求解.对方法的构造过程、收敛性和算法复杂度进行了细致分析.由于算法1收敛速度不够快,为了获得较高的收敛速度,利用加速方法对算法进行加速.仿真表明,此方法具有较强的修复能力,广泛适用于不同的模糊模型,且收敛速度较快,能较好地保留高频细节和轮廓,峰值信噪比(PSNR)较高.今后的工作是更深入的研究不同稀疏基的稀疏性能,包括超小波基、冗余字典基、学习得到的冗余字典基,并将其引入本文方法中.

参考文献

[1]RAM I,COHEN I,ELAD M.Patch-ordering-based wavelet frame its use in inverse problem[J].IEEE Transactions on Image Processing,2014,23(7):2779-2792.

[2]ZHANG Yan,HIRAKAWA K.Blur processing using double discrete wavelet transform[A].Proceedings of the 2013 IEEE Conference on Computer Vision and Pattern Recognition[C].Portland,Oregon,USA:IEEE,2013.1091-1098.

[3]LIU Qie-gen,LIANG Dong,SONG Ying,LUO Jian-hua,ZHU Yue-min,LI Wen-shu.Augmented lagrangian-based sparse representation method with dictionary updating for image deblurring[J].Siam Journal on Imaging Sciences,2013,6(3):1689-1718.

[4]BECK A,TEBOULLE M.Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems[J].IEEE Transactions on Image Processing,2009,18(11):2419-2434.

[5]XU Li,ZHENG Shi-cheng,JIA Jia-ya.Unnaturall0sparse representation for natural image deblurring[A].Proceedings of the 2013 IEEE Conference on Computer Vision and Pattern Recognition[C].Portland,Oregon,USA:IEEE,2013.1107-1114.

[6]LEFKIMMIATIS S,BOURQUARD A,UNSER M.Hessian-based norm regularization for image restoration with biomedical applications[J].IEEE Transactions on Image Processing,2012,21(3):983-995.

[7]王静,吕科,何宁,王茜.基于分裂Bregman方法的全变差图像去模糊[J].电子学报,2012,40(8):1503-1508.

WANG Jing,LV Ke,HE Ning,WANG Qian.Total variation image deblurring based on split Bregman method[J].Acta Electronica Sinica,2012,40(8):1503-1508.(in Chinese)

[8]徐焕宇,孙权森,李大禹,宣丽.基于投影的稀疏表示与非局部正则化图像复原方法[J].电子学报,2013,2014,42(7):1299-1304.

XU Huan-yu,SUN Quan-sen,LI Da-yu,XUAN Li.Projection-based image restoration via sparse representation and non-local regularization[J].Acta Electronica Sinica,2014,42(7):1299-1304.(in Chinese)

[9]LI Qia,MICCHELLI C A,SHEN Li-xin,XU Yue-sheng.A proximity algorithm accelerated by Gauss-Seidel iterations for L1/TV de-noising models[J].Inverse Problem,2012,28(9):095003.

[10]FIGUEIREDO M A,NOWAK R,WRIGHT S J.Gradient projection for sparse reconstruction:Application to compressed sensing and other inverse problems[J].IEEE Journal on Selected Topics Signal Process,2007(1):586-597.

[11]BIOUCAS-DIAS J M,FIGUEIREDO M A.A new twist:two-step iterative shrinkage/thresholding algorithms for image restoration[J].IEEE Transactions on Image Processing,2007,16(12):2992-3004.

[12]AFONSO M V,BIOUCAS-DIAS J M,FIGUEIREDO M A.Fast image recovery using variable splitting and constrained optimization[J].IEEE Transactions on Image Processing,2010,19(9):2345-2356.

[13]DONOHO D.For most large underdetermined systems of linear equations the minimall1-norm solution is also the sparsest solution[J].Communications on Pure and Applied Mathematics,2006,59(6):797-829.

[14]MOREAU J J.Proximité et dualité dans un espace hilbertien[J].Bulletin de la Société Mathématique de France,1965,93:273-299.

[15]BECK A,TEBOULLE M.A fast iterative shrinkage-thresholding algorithm for linear inverse problems [J].SIAM Journal on Imaging Sciences,2009,2 (1):183-202.

[16]FIGUEIREDO M A,NOWAK R,WRIGHT S J.Gradient projection for sparse reconstruction:Application to compressed sensing and other inverse problems[J].IEEE Journal on Selected Topics Signal Process,2007,1(4):586-597.

[17]PORTILLA J.Image restoration through l0analysis-based sparse optimization in tight frames[A].Proceedings of the 16th IEEE International Conference on Image Processing[C].Cairo,Egypt:IEEE,2009.3909-3912.

[18]MICHAELI T,SIGALOV,ELDAR Y C.Partially linear estimation with application to sparse signal recovery from measurement pairs[J].IEEE Transactions on Signal Processing,2012,60(5):2125-2137.

王斌男,1989年10月生于陕西铜川市.现为西安理工大学研究生.研究方向为稀疏表示与压缩感知.

E-mail:1208030272@stu.xaut.edu.cn

胡辽林男,1968年5月生于四川岳池县.现为西安理工大学副教授、博士.主要从事信号处理方面的研究.

E-mail:huliaolin@163.com

Image Restoration Based on Sparse-Optimal Strategy in Wavelet Domain

WANG Bin,HU Liao-lin,CAO Jing-jing,XUE Rui-yang,LIU Guang-fei

(FacultyofMechanicalandPrecisionInstrumentEngineering,Xi’anUniversityofTechnology,Xi’an,Shaanxi710048,China)

Abstract:Due to the nonlinear ill-posing characteristic of image degradation caused by blur and noise,image restoration is generally difficult.To improve image quality,orthogonal wavelet was taken as sparse basis and the sparsity of wavelet coefficients as prior.By constructing a convex function and minimizing this function,we obtained the recovered image.The minimization problem was transferred into a proximal operator which is solved by the fixed point theory.We proved that the optimal solution can be obtained through repeatedly iterating an analytical formula.The construction process,convergence rate and complexity of the method were discussed,and an accelerated algorithm was presented.Simulation results indicate that our method can remove blur and noise,and keep detail information meanwhile.

Key words:blur;wavelet;sparsity;convex function;acceleration

作者简介

DOI:电子学报URL:http://www.ejournal.org.cn10.3969/j.issn.0372-2112.2016.03.016

中图分类号:TP391.4;TP751.1

文献标识码：A

文章编号：0372-2112 (2016)03-0600-07

基金项目:陕西省自然科学基金(No.2014JM7273)

收稿日期:2014-09-03;修回日期:2015-01-04;责任编辑:覃怀银