2016, Vol. 45
[PDF全文]
2. 深圳信息职业技术学院软件学院 广东 深圳 518172;
3. 南昌工程学院信息工程学院 南昌 330099
2. School of Software, Shenzhen Institute of Information Technology Shenzhen Guangdong 518172;
3. Department of Information Engineering, Nanchang Institute of Technology Nanchang 330099
图像的质量在获取和传输过程中不可避免地会退化。如何有效地恢复退化图像，以保证后续图像处理应用是天文、医学成像、遥感、视频通信等众多领域的重要问题，其目标就是从退化的图像中尽可能准确地恢复出原始图像。
图像恢复是图像处理中一类典型的不适定逆问题。图像恢复结果的好坏关键取决于对图像退化过程的先验知识掌握的精确程度。正则化是求解不适定问题的一类重要方法，其基本思路是通过对问题解施加正则性约束来限定解空间，使问题具有良好的特性，从而改善问题性质，获得唯一的稳定解。图像恢复中常用的正则化项包括Tikonov正则化[1]、全变差正则化[2]、小波正则化[3]等。近年来，稀疏正则化理论在数学和信号处理领域受到高度关注，并广泛应用于图像去噪[4]、图像超分辨率重建[5]、压缩感知[6, 7]等领域。稀疏正则化假定图像在某种变换空间中具有稀疏分布特性，并利用这一特性对图像不适定问题进行正则化处理。然而，稀疏正则化仅考虑图像的局域相似结构，易造成图像纹理和结构信息的损失。图像的非局部自相似性模型将传统的图像局部空间邻域扩展到更广义的几何意义下的结构邻域，在更大的区域考虑图像的相似性，更好地利用了图像之间的相似性，在纹理和结构保持方面都有很大的改善，在图像处理的各个领域都到了广泛应用[8, 9, 10]。
本文提出一种非局域自相似约束的Shearlet稀疏正则化图像恢复模型。该模型同时考虑图像的局域和非局域自相似结构，以更好地恢复图像中的结构和纹理等特征；与此同时，引入变量分裂增广拉格朗日法，提出了模型的数值求解方法。通过图像去模糊和图像修复实验，结果表明新的模型能够更有效地恢复降质图像中特征信息。
1 Shearlet变换Shearlet变换[11]由合成小波理论衍生而来，继承了曲线波和轮廓波的优点，是一种新型多尺度几何分析工具。它通过对基本函数缩放、剪切和平移等仿射变换，生成具有不同特征的剪切波函数。具有合成膨胀的仿射函数如下：
| $ \begin{matrix} {{A}_{\mathbf{AB}}}(\mathbf{\psi })= \\ \left\{ {{\mathbf{\psi }}_{j,l,k}}(x)={{\left| \det \mathbf{A} \right|}^{j/2}}\psi ({{B}^{l}}{{A}^{j}}x-k)\ :\ j,l\in Z,k\in {{Z}^{2}} \right\} \\ \end{matrix} $ | (1) |
式中，$\psi \in {L^2}({R^2})$；${\bf{A}}$和${\bf{B}}$是$2 \times 2$可逆矩阵, $\left| {\det {\kern 1pt} {\bf{B}}} \right|$ $ = 1 $。矩阵 ${A^j}$ 与尺度变换相关联，${B^l}$与保持面积不变的几何变换(旋转和剪切)相关联。
当${A_a} = \left[ {\begin{array}{*{20}{c}} a&0\\ 0&{\sqrt a } \end{array}} \right]$，${B_s} = \left[ {\begin{array}{*{20}{c}} 1&s\\ 0&1 \end{array}} \right]$时，产生的函数集为：
| $\left\{ {{\psi _{ast}}(x) = {a^{ - {\textstyle{3 \over 4}}}}\psi {\rm{(}}A_a^{ - 1}B_s^{ - 1}x - t),a \in {R^ + },s \in R,t \in {R^2}} \right\}$ | (2) |
称为Shearlet函数。由此，对任意给定的函数$f \in {L^2}({R^2})$，其Shearlet变换定义为：
| ${\rm{S}}{{\rm{H}}_\psi }f(a,s,t) = \left\langle {f,{\kern 1pt} {\psi _{ast}}} \right\rangle {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} ,{\rm{ }}a \in {R^ + },{\rm{ }}s \in R,{\rm{ }}t \in {R^2}$ | (3) |
假定图像的退化模型为：
| ${\bf{y}} = {\bf{Hx}} + {\bf{n}}$ | (4) |
式中，x为原始图像； H为降质函数； y为观测图像；n为噪声。正则化约束的图像恢复问题等同于求解以下能量泛函的变分问题：
| $ {\bf{\hat x}} = \mathop {\arg \min }\limits_{\bf{x}} \left\{ {\frac{1}{2}\left\| {{\bf{y}} - {\bf{Hx}}} \right\|_2^2 + \lambda \varphi ({\bf{x}})} \right\} $ | (5) |
式中， $ \left\| {y - Hx} \right\|_2^2 $为保真项；$\varphi (x)$为正则化项；λ为
正则化参数。常用的正则化项包括：Tikonov[1]、全变差[2, 12]、稀疏[3, 4, 5, 13]等。
在图像恢复模型中，正则化项通常用于描述图像先验。全变差[2, 12]和图像稀疏模型[3, 4, 5, 13]是当前最常用的正则化方法。全变差正则化假定图像分片光滑的前提下最小化图像的全变差。全变差正则化仅考虑像素层面的光滑与连续，没有考虑图像的结构和相似性等，易造成图像的过度平滑，只适用于含少量纹理和结构的图像。图像稀疏模型表明自然图像在合适的基函数下总存在稀疏的表示，并且非零大系数能够描述图像内在的结构与本质特征，为图像恢复提供了新的思路。然而，由于图像作为复杂的二维信号，通常包含平滑、边缘和纹理等不同的形态结构，单一的基函数很难同时对图像的多种结构进行最稀疏的表示[12, 13]。相比于单一的基函数如小波等，图像块稀疏模型[4]假定每一图像块可由给定的字典进行稀疏表示，在图像恢复领域表示出更好的优势。虽然基于图像块稀疏模型能较好地刻画图像的局部结构，但它仅考虑各像素点及其周围的局部结构信息，并且对每个块独立地进行约束，而没有考虑到不同图像块之间的相似性。另外，图像块稀疏模型算法复杂、运算时间较长，并不适用于图像恢复等图像预处理操作。
本文提出一种非局域自相似约束的Shearlet稀疏正则化图像恢复模型。一方面，模型引入Shearlet稀疏表示，它能更稀疏地表示图像的几何结构，同时由于其简单和严谨的数学框架，具有快速离散算法，可以解决单一基函数非稀疏和图像块稀疏算法复杂问题。另一方面，模型引入非局域自相似来刻画图像中不同位置的像素点或图像块之间的相似性，更好地保持了图像的整体结构特征。
2.2 新的模型Shearlet变换具有简单的数学结构，通过对图像频率空间的逐层剖分，很好地捕获图像的方向几何结构，并对其进行稀疏表示。这些优良特性使得Shearlet更加适用于图像的恢复。自然图像中存在大量重复相似结构信息，这种相似结构信息不仅包括平滑区域，而且在纹理区域和边缘部分，同样存在大量结构相似性。它们不仅像素值相同，图像块结构也基本相似。图像的这种非局域相似特征，将有利于提高图像恢复的性能。鉴于Shearlet稀疏表示和非局部自相似的图像处理中的优势，结合稀疏正则化和和非局部自相似，本文提出了非局域自相似约束的Shearlet稀疏正则化图像恢复模型。新的模型描述如下：
| $ {\bf{\hat x}} = \mathop {\arg \min }\limits_{\bf{x}} \left\{ {\frac{1}{2}\left\| {{\bf{y}} - {\bf{Hx}}} \right\|_2^2 + \lambda {{\left\| {\left\langle {{\bf{x}},{\psi _{{\rm{ast}}}}} \right\rangle } \right\|}_1} + \gamma \left\| {{\bf{x}} - {\bf{wx}}} \right\|_2^2} \right\} $ | (6) |
式中， $ {\left\| {\left\langle {{\bf{x}},{\psi _{{\rm{ast}}}}} \right\rangle } \right\|_1} $为Shearlet稀疏正则化项； $ \left\| {x - wx} \right\|_2^2 $为非局域自相似正则化项；γ为正则化参数；w为
非局域自相似权重矩阵。为简化，将式(6)写为：
| $ {\bf{\hat x}} = \mathop {\arg \min }\limits_{\bf{x}} \left\{ {\frac{1}{2}\left\| {{\bf{y}} - {\bf{Hx}}} \right\|_2^2 + \lambda {{\left\| {S{\bf{x}}} \right\|}_1} + \gamma \left\| {{\bf{x}} - {\bf{wx}}} \right\|_2^2} \right\} $ | (7) |
式中， $ S $表示Shearlet变换矩阵。
图像的非局域自相似性[8, 9, 10]描述为处于图像中不同位置的像素点或图像块往往表现出很强的相关性，如自然图像中的长边缘、纹理图像以及具有周期性结构的图像等，都表现出远距离的相关性。基于此，非局域相关性的强弱通常采用非局域自相似权重矩阵${\bf{w}}$度量，其定义为：
| $ {\bf{w}}(i,j) = \left\{ {\begin{array}{*{20}{l}} {{w_{ij}}}&{{\rm{1}} \le {\rm{j}} \le {\rm{J}}}\\ 0&{{\rm{其他 }}} \end{array}} \right. $ | (8) |
式中， $ {w_{i,j}} $表示图像块${{\bf{x}}_i}$与其相对应的第 j个图像块${\bf{x}}_i^j$的非局部自相似权重值，${w_{i,j}}$值越大，说明图像块之间越相似。其定义为：
| $ {w_{ij}} = \frac{1}{{{c_i}}}\exp ( - ||{{\bf{x}}_i} - {\bf{x}}_i^j||_2^2/{h^2}) $ | (9) |
式中， $h$表示自相似控制因子；${c_i}$为归一化因子。可知，某一图像块与当前图像块越相似，则分配的权重越大。
2.3 模型求解式(6)和式(7)是组合约束最小化能量泛函问题，对其进行数值求解是实现图像恢复的关键。增广拉格朗日算法框架是计算优化领域的典型方法。文献[12]在变量分裂的基础上提出了用增广拉格朗日函数优化方法，并用其来处理稀疏正则化问题。基于变量分裂增广拉格朗日法具有参数选取鲁棒、迭代收敛快等特性。因此，本文采用基于变量分裂增广拉格朗日法，对该新的图像恢复变分模型进行数值求解。
引入中间变量${\bf{u}}$，则无约束优化问题式(7)可变换为以下约束优化问题：
| $ \begin{array}{c} {\bf{\hat x}} = \mathop {\arg \min }\limits_{\bf{x}} \left\{ {\frac{1}{2}\left\| {{\bf{y}} - {\bf{Hu}}} \right\|_2^2 + \lambda {{\left\| {{\bf{Sx}}} \right\|}_1} + \gamma \left\| {{\bf{u}} - {\bf{wu}}} \right\|_2^2} \right\}\\ {\rm{s}}{\rm{.t}}{\rm{.}}\quad {\bf{u}} = {\bf{x}} \end{array} $ | (10) |
对应的增广拉格朗日函数为：
| $ \begin{array}{c} L({\bf{u}},{\bf{x}},{\bf{d}}) = \frac{1}{2}\left\| {{\bf{y}} - {\bf{Hu}}} \right\|_2^2 + \lambda {\left\| {{\bf{Sx}}} \right\|_1} + \\ \gamma \left\| {{\bf{u}} - {\bf{wu}}} \right\|_2^2 + \frac{\mu }{2}\left\| {{\bf{x}} - {\bf{u}} - {\bf{d}}} \right\|_2^2 \end{array} $ | (11) |
式中，$\mu > 0$为正则化参数； $ d/\mu $表示增广拉格朗日乘子。增广拉格朗日法通过寻找$L\left( {{\bf{u}},{\bf{x}},{\bf{d}}} \right)$的鞍点求
解式(10)。因此，约束问题式(10)的解可以通过迭代求解以下两式得到：
| $ ({{\bf{u}}_{k + 1}},{{\bf{x}}_{k + 1}}) = \mathop {\arg \min }\limits_{{\bf{u}},{\bf{x}}} L({\bf{u}},{\bf{x}},{\bf{d}}) $ | (12) |
| $ {{\bf{d}}_{k + 1}} = {{\bf{d}}_k} - ({{\bf{x}}_{k + 1}} - {{\bf{u}}_{k + 1}}) $ | (13) |
采用交替方向法，最优化问题式(12)可转换成两个次优化问题。固定变量 x ，式(12)等价为：
| $ \begin{array}{l} {{\bf{u}}_{k + 1}} = \mathop {\arg \min }\limits_{\bf{x}} \frac{1}{2}\left\| {{\bf{y}} - {\bf{Hu}}} \right\|_2^2 + \\ \gamma \left\| {{\bf{u}} - {\bf{wu}}} \right\|_2^2 + \frac{\mu }{2}\left\| {{\bf{x}} - {\bf{u}} - {\bf{d}}} \right\|_2^2 \end{array} $ | (14) |
显然，式(14)是二次优化问题，令导数等于零，其解可以直接求出，即：
| $ \begin{array}{l} {{\bf{u}}_{k + 1}} = [{{\bf{H}}^{\rm{T}}}{\bf{y}} - \mu ({{\bf{x}}_k} - {{\bf{d}}_k})]/\\ [{{\bf{H}}^{\rm{T}}}{\bf{H}} - 2\gamma {({\bf{I}} - {\bf{w}})^{^{\rm{T}}}}({\bf{I}} - {\bf{w}}) - \mu ] \end{array} $ | (15) |
同理，固定变量${\bf{u}}$，式(12)等价为：
| $ {{\bf{x}}_{k + 1}} = \mathop {\arg \min }\limits_{\bf{u}} \lambda {\left\| {{\bf{Sx}}} \right\|_1} + \frac{\mu }{2}\left\| {{\bf{x}} - {\bf{u}} - {\bf{d}}} \right\|_2^2 $ | (16) |
式(16)是典型的${\ell _1}$稀疏正则化优化问题，其解可用软阈值法得到，即：
| $ {{\mathbf{x}}_{k+1}}={{\mathbf{S}}^{-1}}(\text{soft(}{{\mathbf{S}}^{\text{T}}}{{\mathbf{u}}_{k+1}}-{{\mathbf{S}}^{\text{T}}}{{\mathbf{d}}_{k}},\,{\lambda }/{\mu }\;\text{)}) $ | (17) |
给定待恢复图像${\bf{y}}$，非局域自相似约束的Shearlet稀疏正则化图像恢复算法具体步骤描述如下：
1) 初始化：迭代次数$k = 0$，正则化参数$\gamma $和$\mu $、拉格朗日乘子${{\bf{d}}_0}$、中间变量 $ {{\bf{u}}_0} $、复图像 $ {{\bf{x}}_0} $及最大迭代次数$K$；
2) 模型迭代求解：分别由式(15)、式(17)和式(13)，计算得到 $ {{\bf{u}}_{k + 1}} $、 $ {{\bf{x}}_{k + 1}} $和 $ {{\bf{d}}_{k + 1}} $；
3) 更新迭代次数： $k = k + {\rm{1}}$ ；
4) 判定终止条件：如果$k \le K$，转到步骤2)迭代；否则，结束迭代，输出恢复图像${{\bf{x}}_{k + 1}}$。
3 实验结果与分析为了验证算法的有效性，使用标准的512×512像素的灰度图像Lena、Barbara、Boats和Bridge进行去模糊和修复实验。实验中各参数设定为$K = 100$，$\gamma = \mu = 1$。恢复效果采用峰值信噪比(PSNR)和结构自相似指标(SSIM)[14]衡量。
3.1 去模糊实验去模糊实验中，降质函数${\bf{H}}$选定方差为1的7×7像素高斯函数，并加入高斯白噪声，使降质图像的模糊信噪比为40 dB。表 1给出了各算法去模糊后图像的PSNR和SSIM。其中TwIST为两步迭代收缩算法[13]，C-SALSA为约束的分裂增广拉格朗日算法[12]。由表 1可以知，本文算法去模糊后的图像具有更高的PSNR和SSIM数值，有更好的去模糊性能。
| 表1 各算法去模糊后，恢复图像PSNR和SSIM值比较 |
 |
表 1 各算法去模糊后，恢复图像PSNR和SSIM值比较
图 1给出了3种算法对降质的Lena图像去模糊后恢复图像的局部细节。从图 1可以看出，TwIST和C-SALSA算法虽然恢复降质图像的整体效果，但边缘轮廓等细节仍然较模糊，并且出现了振铃效应。本文算法同时考虑稀疏表示特性和图像自身非局域自相似特性，能够更好地恢复边缘轮廓特征，并且具有更好的整体视觉效果。
 |
| 图1 3种算法对Lena图像复原后的局部视觉质量比较 |
图像修复实验中，降质图像分别通过随机丢失40%和60%的像素得到。表 2给出随机丢失40%像素后，各算法修复图像的PSNR和SSIM。由表 2可知，本文算法修复图像的PSNR和SSIM明显高于TwIST和C-SALSA算法。
| 表2 随机丢失40%像素后，各算法恢复图像PSNR和SSIM值比较 |
|
图 2给出了Barbara图像随机丢失60%像素后，3种算法修复的图像。由图 2b和图 2c可以看出，TwIST算法仅能恢复部分丢失的数据，而C-SALSA算法能较为完整地恢复丢失数据，但在边缘区域出现了严重的模糊。本文算法能够很好地恢复边缘等结构特征，获得更好的视觉效果。
 |
| 图2 丢失60%像素后，3种算法恢复图像的局部视觉质量比较 |
为了便于进行算法复杂性分析，假定TwIST和C-SALSA均采用Shearlet变换。Shearlet变换的实现包含快速傅里叶变换和拉普拉斯金字塔算法两部分，其复杂度为 $ O(N\log N) $，其中$N$为图像的像素数。3种算法均为迭代算法，非局域正则子计算可采用快速傅里叶的方式，其时间复杂度为$O\left( N \right)$，稀疏正则子计算复杂度仅由Shearlet变换的算法复杂度决定。由于TwIST采用双迭迭代收缩，需要进行一次Shearlet变换和两次收缩处理，每次迭代其算法复杂度为$O(N\log N)$。C-SALSA每次迭代进行一次Shearlet变换和收缩处理，每次迭代其算法复杂度为$O(N\log N)$。而本文算法每次迭代除进行一次Shearlet变换和收缩处理外，还进行一次非局域正则化处理，每次迭代其算法复杂度为$O(N\log N) + $$O(N)$。
 |
| 图3 丢失60%和40%像素后，3种算法恢复图像的MSE 随时间的变化曲线 |
为了验证算法的收敛速度，图 3给出了像素丢失率为40%和60%时Barbara恢复图像的均方误差(MSE)随运行时间的变化曲线图。从图 3中可以看出，随运行时间的增加，MSE不断下降。本文算法的MSE随运行时间下降最快，且比其他算法的MSE小，表明该算法收敛速度更快、恢复效果好。
4 结 束 语图像恢复是一类不适定逆问题。正则化是求解该类问题的热点方法之一。本文引入非局部自相似和Shearlet稀疏性正则化，提出一种混合正则化图像恢复变分模型。Shearlet稀疏性表征图像几何结构，非局部自相似表征图像的全局自相似。图像去模糊和图像修复实验结果表明，本文算法能够较完整地恢复图像整体，同时较好地恢复了图像的边缘等结构特征，具有更好的视觉效果，并获得了更高的PSNR和SSIM值。
| [1] | KLANN E, MAASS P, RAMLAU R. Tikhonov regularization with wavelet shrinkage for linear inverse problems[J]. Journal for Inverse and Ill-Posed Problems, 2006, 14(6): 583-609. |
| [2] | LI Wei-hong, LI Quan-li, GONG Wei-guo, et al. Total variation blind deconvolution employing split bregman iteration[J]. Journal of Visual Communication and Image Representation, 2012, 23(3): 409-417. |
| [3] | DAUBECHIES I, DEFRISE M, de MOL C. An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[J]. Communications on Pure and Applied Mathematics, 2004, 57(11): 1413-1457. |
| [4] | ELAD M, FIGUEIREDO M A T, MA Y. On the role of sparse and redundant representation in image processing[J]. Proceedings of the IEEE-Special Issue on Applications of Sparse Representation and Compressive Sensing, 2010, 98(6): 972-982. |
| [5] | 邓承志, 田伟, 汪胜前, 等. 近似稀疏正则化的红外图像 超分辨率重建[J]. 光学精密工程, 2014, 22(6): 1648-1654. DENG Cheng-zhi, TIAN Wei, WANG Sheng-qian, et al. Super-resolution reconstruction of approximate sparsity regularized infrared images[J]. Opt Precision Eng, 2014, 22(6): 1648-1654. |
| [6] | DONOHO D L. Compressed sensing[J]. IEEE Transactions on Information Theory, 2006, 52(4): 1289-1306. |
| [7] | 吴迪, 王奎民, 赵玉新, 等. 分段正则化正交匹配追踪算 法[J]. 光学精密工程, 2014, 22(5): 1395-1402. WU Di, WANG Kui-min, ZHAO Yu-xin, et al. Stagewise regularized orthogonal matching pursuit algorithm[J]. Opt Precision Eng, 2014, 22(5): 1395-1402. |
| [8] | 李家德, 张叶, 贾平. 采用非局部均值的超分辨率重构[J]. 光学精密工程, 2013, 21(6): 1576-1585. LI Jia-de, ZHANG Ye, JIA Ping. Super-resolution reconstruction using nonlocal means[J]. Opt Precision Eng, 2013, 21(6): 1576-1585. |
| [9] | WANG J, JEON G, JEONG J. Deinterlacing algorithm with an advanced non-local means filter[J]. Optical Engineering, 2012, 51(4): 47-49. |
| [10] | ZHANG K B, GAO X B, TAO D C, et al. Single image super-resolution with non-local means and steering kernel regression[J]. IEEE Trans Image Process, 2012, 21(11): 4544-4556. |
| [11] | EASLEY G R, LABATE D, LIM W. Sparse directional image representations using the discrete shearlet transforms[J]. Applied and Computational Harmonic Analysis, 2008, 25(1): 25-46. |
| [12] | AFONSO M, BIOUCAS-DIAS J, FIGURIREDO M. An augmented lagrangian approach to the constrained optimization formulation of imaging inverse problems[J]. IEEE Transactions on Image Processing, 2011, 20(3): 681-695. |
| [13] | BIOUCAS-DIAS J M, FIGUEIREDO M A T. A new TwIST: two-step iterative shrinkage/thresholding algorithm for image restoration[J]. IEEE Transactions on Image Processing, 2007, 16(12): 2992-3004. |
| [14] | BRUNET D, VRSCAY E R, WANG Z. On the mathematical properties of the structural similarity index[J]. IEEE Transactions on Image Processing, 2012, 21(4): 1488-1499. |