董伟胜 中国西安电子科技大学电子工程学院 wsdong@mail.xidian.edu.cn 张磊 香港理工大学计算机系 cslzhang@comp.polyu.edu.hk 石光明 中国西安电子科技大学电子工程学院 gmshi@xidian.edu.cn
摘 要
本文对于图像恢复任务提出了一种新的称为集中稀疏表示(CSR)的稀疏表示模型。为了重建高还原度的图像，通过给定的字典，退化图像的稀疏编码系数预计应该尽可能接近那些未知的原始图像。然而，由于可用的数据是原始图像的退化版本(如噪声、模糊和/或者低采样率)，正如许多现有的稀疏表示模型一样，如果只考虑局部的稀疏图像，稀疏编码系数往往不够准确。为了使稀疏编码更加准确，通过利用非局部图像统计，引入一个集中的稀疏性约束。为了优化，局部稀疏和非局部稀疏统一到一个变化的框架内。大量的图像恢复实验 验证了我们的CSR模型在以前最先进的方法之上取得了令人信服的改进。
1、介绍
图像恢复（IR）目的是为了从，比如说通过一个低端摄像头或者在有限条件下得到图像的图像退化版本（例如噪声、模糊和/或者低采样率）， 来恢复一副高质量的图像。对于观察的图像y，IR问题可以表示成：
y = Hx + v (1) 其中H是一个退化矩阵，x是原始图像的矢量，v是噪声矢量。由于IR的病态特性，尝试把观察模型和所需解决方案的先验知识合并到一个变分公式的正则化技术，已经被广泛地研究。对于正则方法，对自然图像适当的先验知识进行寻找和建模是最重要的关注点之一，因此学习自然图像先验知识的各种方法已经被提出来了【25,5,6,12】。
近年来，对于图像恢复基于建模的稀疏表示已经被证明是一种很有前途的模型【9,5,13,20,16,21,27,15,14】。在人类视觉系统【23,24】的研究中，已经发现细胞感受区域使用少量的从一个超完备的编码集中稀疏选出的结构化基元来编码自然图像。在数学上，一个x ∈ R的信号可以表示为一个字典中的几个原子的线性组合，例如，X ≈ ,用|0 最小化：
N
x = argmin ||||0, s.t.||x − Φ||2 < , (2)
α
其中||||0 对中的非零系数进行统计，是一个很小的平衡稀疏性和近似误差的常数。在实践中，为了高效的凸优化【13,9】，|0范数通常被|1 范数所替代。另一方面，对照分析设计字典（比如小波字典）、学习字典【1，21,19】，例如图像分块可以更好地适应信号x和描述图像结构特征。 在IR的方案中，我们有原始图像x的一些退化观察值y，即y = Hx + v。为了从y中重建x，首先y通过求解下列最小化问题从Φ进行稀疏编码：
y = argmin ||||1, s.t.||y − H Φ||2 < , α
(3)
然后重构的x，由xˆ表示，得到xˆ=Φy。显然，式（3）中的y预期上与x非常接近，ˆ可以与真实图像x非常接近。不幸的是，因为y受噪声干扰，是所以估算的图像x模糊和不完整的，式（3）中的编码矢量y可能与所需矢量x偏差很大，导致图像x的恢复不准确。换句话说，模型式（3）可以确保式y是稀疏的，但不能确保y与x无限接近。 在本文中，我们介绍稀疏噪音编码（SCN）的概念来方便问题讨论。y的SCN定义为： Va = x - y, (4) 我们可以看到在给定的字典Φ中，IR的结果取决于SCN Va的水平，因为图像重建误差Vx= xˆ-x≈Φy-Φx=ΦVa。SCN Va的定义还表明了一种提高IR质量的方法，即降低Va的水平。 常规的稀疏表示模型，比如在式（2）或式（3）中的，主要利用了自然图像（基于分块）的局部稀疏性。每一个图像块独立编码，不考虑其他块。尽管如此，稀疏编码系数不是随机分布的，因为图像局部分块是非局部相关的。非局部平均值（NLM）方法，其目的是利用图像非局部冗余，已经成功地运用在许多图像处理应用中，特别是在去噪方面【6】。这意味着可以利用图像的非局部相似性来减少SCN，因此，复原图像的质量可以得到改善。事实上，最近的一些研究工作，比如【10】和【20】，都是基于这样的考虑的。例如，在【20】中提出了一组能够同时编码相似的块的稀疏编码方案，并取得了很好的去噪效果。 在本文中，我们提出了一个能够有效减少SCN的集中系数表示（CSR）模型，从而提高了基于稀疏的IR性能。其基本思想是将图像的局部稀疏约束（即一个局部块可以被一些从一个字典稀疏选择的原子进行编码）和集中稀疏约束（即稀疏编码系数应该接近其平均值）整合到一个统一的变化的优化框架中。具体而言，除了要求每个局部块的编码系数是稀疏的之外，我们也通过利用非局部相似性引起的稀疏性来强制让稀疏编码系数具有小的SCN，这一点可以用|1 范数来表示其特征。通过进行大量的IR实验，实验结果表明，所推荐的CSR算法明显优于许多最顶尖的IR方法。 2、集中稀疏表示模型
2.1 图像恢复的稀疏编码噪音
按照【16】中所用的符号，我们用x∈R 表示原始图像，用xi=xRi表示在一个大小为n×n的图像块中位置为i，其中Ri在i位置上的从x提取的矩阵块xi，给定一个字典Φ∈Rn×M，n 其中x表示所有i的连接。上面的公式只是说明整体图像已经通过平均xi的每一个重建块完成了重建。 （a） (b) (c) (d) 图1. Lena图像（a）是有噪音和模糊时的SCN分布；（b）是低采样率时的。（c）和（d）分别表示在log域中的（a）和（b）的相同分布。 在IR的应用中，x是不能编码的，我们所有的只是观察的退化的y=Hx + v。x的稀疏编码是基于y通过最小化 y = argmin { || y - H Φ◦ ||2 || ||1 }, (6) 2 + α 然后图像就被重建成xˆ =Φ◦y。正如我们在式（4）中所定义和讨论的，系数y将偏ˆ的图像恢复质量。 离x，并且稀疏编码噪声（SCN）Va = y-x决定了x在此我们进行了一些实验来调查SCN Va的统计。我们使用Lena图像作为例子。原始图像x首先被模糊（用一个标准差为1.6的高斯模糊内核），再加上一个标准差为2的高斯白噪声来获得一个有噪声和模糊的图像y。然后我们分别通过最小化 Φ◦ ||2 + || ||1 }， (7) x = argmin { || x - α 和式（6）计算x和y。DCT字典在这个实验被采用。然后通过Va=y-x计算SCN。在图1（a）中，我们绘出对应字典中的第四个原子（其他原子的分布类似）对应的Va分布。在图1（b）中，当观察的数据y刚开始模糊（由一个标准差为1.6的高斯模糊内核）然后采样后绘出Va分布。我们可以看到SCN Va的实验分布在零处达到最高值，拉普拉斯算子函数可以很好地表示其特点，但高斯函数有很多更大的拟合误差。在图1（c）和（d）中我们log域显示这些分布以便更好地观察拟合的尾巴。这个观察促使我们使用之前的拉普拉斯算子来对SCN建模，这些将在第三节介绍。 2 2.2 集中稀疏表示 ˆ，然而困难在于系数矢量x是未知 很明显，抑制SCN Va有助于改善IR的输出xˆ 的，所以不能直接测量Va。但是，如果我们能进行一些对x的合理估计，用αx表示，然后就可以用αy − αˆ x表示x的一个近似值。直观的说，为了抑制x和提高y的精确度，有一个新的稀疏编码模型： 其中 是一个常数，lp （p可以是1或2）范数用来测量与αˆ x 间的距离。与式（6）相比，式（8）使αy 接近αˆ x（因此SCN Va能被抑制）从而使αy是稀疏的，因此致使稀疏编码比通过求解式（6）得来的更为可取。如果αˆ x=0并且p=1，式（8）的模型将简化为式（6）的常规模型。 现在问题转化为如何找到一个未知矢量x的合理估计。通过将x看做一个随机变x的一个很好的无偏估计自然是它的平均值；ˆ 矢量，也就是说，我们可以令αx=E[x]。事实上，我么可以通过假设SCN的Va几乎是零均值（请参阅2.1节中的经验观察值）ˆ 来由E[y]近似E[x]，并有αx=E[x]≈E[y]。然后式（8）可以转化为 我们把上面的模型称为集中稀疏表示（CSR），因为它使稀疏编码系数强制接近其分布中心（即平均值）。 对于每一个图像块i的稀疏编码i，如果我们有足够的i的样本，E[i]是可以近 似计算的。幸运的是，给定的块i在自然图像中通常有许多非局部相似块。然后E[i]可以通过与块i相关的非局部相似块（包括块i）的那些稀疏编码矢量的加权平均数来计算。为此，我们可以对每一个块i通过块匹配而形成一个簇，用Ci表示，然后计算每个簇里稀疏编码的平均值。用i,j表示找到的与块i相似的块j的稀疏编码。然后我们用所有i,j的加权平均值来近似于E[i]，即 i = wj∈Cii,jai,j （10） 其中wi,j 是权重。与非局部平均方法【6】中相似，wi,j可以被设置成与块i与块j间距离成反比： wi,j = exp( - || xˆi - xˆi,j ||2 2/h) / W， （11） 其中xˆi = Φαˆ i 和xˆi,j = Φαˆ i,j 是块i和块j的估计，W是一个归一化因数，h是一个预定的标量。 通过用 i 作为E[i]的估计，式（9）中的CSR模型可以写成： 从式（12）中我们可以更清楚地看到CSR模型结合局部稀疏性（即||||1）和非局部相似性诱导稀疏性（即||i- i||lp）到一个统一的变化公式。通过利用局部和非局部的冗余，更好的IR结果是可以预期的。 实际上式（12）表明了一个CSR模型的迭代最小化方法。我们初始化 i为0，即 i(-1) =0。然后从一些初始的稀疏编码结果，表示为y(0)，我们可以得到x的初始估计，表示为x(0)，通过x(0)= Φ◦y(0)。在x(0)的基础上，我们可以找到每一个局部块i的相似块，并因此每个块的非局部平均编码矢量，即i，就可以根据y(0) 通过式（10）和式（11）来更新。更新后的平均值由 i (0) 表示，将会在CSR模型的下一轮稀疏编码处理中用到。这个步骤被重复使用直到收敛为止。在第j次迭代中，稀疏编码由下式执行： 从上面的讨论中可以看出，稀疏编码和非局部聚类步骤在拟定的CSR方案中交替执行。在迭代过程中，稀疏编码y(j)的准确性是逐步提高的，这又反过来提高了非局部聚类的准确性，非局部聚类准确性的提高也进一步提高了稀疏编码的准确性。最后，当连接稀疏编码和非局部聚类过程达到一个局部最小值就得到了所需的稀疏编码y。请注意，式（12）的模型不是凸的，但当均值是 i固定的时则是凸的。也就是说，式（13）中的稀疏编码步骤当p≥1时是凸的。 3、CSR算法 3.1 和参数定义 在拟定的CSR模型中，有两个参数和，两者分别平衡局部冗余引起的稀疏性和非局部冗余引起的稀疏性。这两个参数可以使用训练集通过经验来设置。然而，一个更为合理和适合的对两者的设置，不仅能提高收敛速度，还能在很大程度上提高IR质量【7】。在本章节中，我们提供了一个CSR模型稀疏编码步骤的贝叶斯解释，这也给我们一种明确的方式来确 定参数和。事实上，在基于小波的图像去噪【26】中，小波表示法和贝叶斯框架之间的联系已经很好地建立起来了。这样的连接有助于调节确定方法和概率方法之间的差异。 从2.2节的讨论中我们知道CSR模型有两个输出，稀疏编码y和均值i，然而我们只对前者感兴趣，因为图像是从y中恢复来的。换句话说， i或者E[y]是用于优化过程中的隐藏变量。为了方便下面的发展，我们让=-E[]。然后给定观察值y，和的最大后验概率（MAP）估计可以按下式计算： 按照贝叶斯公式，我们可以得到： 这两个式子分别对应于可能性和前一个式子。根据式（1）的观察模型，可能性式子被表征为： 其中n 是添加的高斯噪声标准差。 我们由经验发现和几乎是不相关的。例如，我们在去模糊实验所用的九个图像（请 参考4.1节），和间的相关系数范围是从0.039到0.153。因此，本文中我们假设和是相互独立的并且是独立同分布随机变矢量。如图1所示，拉普拉斯分布可以很好的表征SCN。因此，我们可以用独立同分布拉普拉斯算子分布给SCN信号建模。同时，稀疏稀疏在文献上用独立同分布拉普拉斯算子分布来表征也是可以接收的。因此，式（15）中的后验概率可以表示为： 其中i和i分别是和的在第i次的值，i和i分别是i和i的标准差。 把式（16）和（17）代入式（15），我们可以得到： 通过比较式（18）和式（9），很明显应该选择l1 范数（即p=1）来表征集中稀疏术语，使得的贝叶斯最优估计得以实现。这仅仅是因为分布可以很好地由拉普拉斯算子建模。因此，式（9）中的CSR模型可以指定为 最终通过比较式（18）和（19），我们得到 在实现过程中，i和i是从收集到的非局部相似块的i和i集合中估计出来的。这种估计比那些只使用局部块的更为强大。然后，在每次迭代中（或在多次迭代中）i和i随着和的更新而更新。 3.2 字典Φ的选择 在前面的章节中，我们提出CSR模型和相关算法是基于假设字典Φ已经给定。提议的CSR模型是一个总模型并且Φ的选择是不同的。例如，可以使用小波字典，或者一个通过使用像奇异值分解KSVD【1】算法的示例图像获得的学习字典。然而，分析设计的小波字典和学习型KSVD字典是通用字典；也就是说，它们可以用来表示任何的图像块，但是它们可能没有足够的灵活性来稀疏表示一个给定的局部块。在提出的CSR算法中，许多与给定块类似的块被收集起来。这促使我们对每一个块使用一个与之相适应的字典。我们了解的是每一个块的或每一个相似块簇的字典，而不是一个通用字典。具体地说，我们对每一个簇应用PCA来获得一个PCA学习字典，并用这个字典来编码这个簇中的块。这样一种PCA策略已经在【28,8,15,14】中的IR中应用。 3.3 算法总结 从前面章节的分析我们可以看到，式（9）和式（12）中的CSR模型可以迭代求解。从一些初始化工作中可以看出，一旦解决了y的稀疏编码，就可以计算出非局部平均值i，然后可以通过式（19）来更新y，等等。式（19）有两个l1 范数限制，并且可以通过增广拉格朗日乘数（ALM）方法大体求解。在本文中，我们扩展了一个在【13】中的从一个l1 约束到两个l1 约束的迭代收缩（IS）算法。尽管IS算法需要多次迭代才收缩，但是每一次迭代都非常简单。更重要的是，式（20）中的参数i和i可以直接被转换到IS算法中的阈值。由于篇幅限制，我们忽略这里的细节，有兴趣的读者可以参考【13】。 本文提出的基于CSR的图像恢复算法的主要步骤在算法1中进行总结。 算法1：CSR图像恢复 初始化：利用标准稀疏模型【13】计算初始估计值xˆ； 外循环：在l = 1,2,3，…….，L上迭代 用PCA对每一个相似块的簇更新字典； 用式（20）更新正则化参数（和）； 从稀疏编码y(l-1) 计算非局部平均值i(l-1)； 由扩展迭代收缩算法【13】通过求解式（19）计算y(l) 。 4、实验结果 我们进行了大量的IR实验来证明所提出的CSR模型的性能。在图像去噪方面，CSR可以取得与BM3D【10】和分组稀疏编码方法【20】非常相似的结果。由于篇幅限制，在本节中我们只展示图像去模糊和超分辨率的结果。在我们的实验中块的大小为6×6,。CSR方法的源代码和更多的实验结果可以在http://www.comp.polyu.edu.hk/∼cslzhang/CSR.htm中找到。 4.1 图像去模糊 CSR模型的去模糊性能在模拟的模糊图像和真实的模糊图像上进行了验证。为了模拟一个模糊图像，原始图像由一个模糊内核进行模糊，然后加入标准差n =2和n = 2的加性高斯噪声。用于模拟的两个模糊内核：一个9×9的均衡内核和一个标准差为1.6的高斯模糊内核。对于真正的运动模糊图像，我们借用了【17】中的内核估计方法来估计模糊内核。 我们比较了CSR模糊方法和几个最近发展的去模糊方法，例如受约束的电视去模糊方法（表示为FISTA）【2】、SA-DCT去模糊方法【18】以及BM3D去模糊方法【11】。需要注意的是FISTA是最近提出的基于电视的去模糊方法，能很好地重建分段光滑区域。SA-DCT和BM3D是两个著名的图像恢复方法，往往会有最顶尖的图像去模糊效果。 表1是一组9幅摄影图像的PSNR（峰值信噪比）结果报告。从表1中，我们可以得到 这样的结论：对于均匀模糊和高斯模糊，提出的CSR去模糊方法显著优于其他竞争方法。图2是去模糊方法的视觉比较，从中我们可以看到，与其他方法相比，CSR模型使图像边缘和纹理更加整洁与清晰。在补充材料中有更多的例子。 表1 去模糊图像的PSNR(dB)结果 Images FISTA [2] SA-DCT [18] BM3D [11] CSR Butterfly 28.37 27.50 27.21 29.75 Boats 29.04 29.25 29.93 31.10 C. Man 26.82 27.02 27.30 28.55 √9×9 均匀模糊, σn = 2 Parrot 29.11 30.07 30.50 32.09 Lena 28.33 28.43 28.97 29.95 Starfish 27.75 28.11 28.61 30.30 Barbara 25.75 26.58 27.99 27.93 Peppers 28.43 28.36 28.60 29.64 Leaves 26.49 27.04 27.45 29.97 Average 27.79 28.04 28.51 29.95 FISTA [2] SA-DCT [18] BM3D [11] CSR 27.73 26.46 26.56 28.66 27.93 28.14 29.21 30.48 26.13 26.11 26.61 27.68 9×9均匀模糊,σn = 2 27.50 27.14 27.97 29.00 28.88 29.10 29.75 30.57 27.40 27.58 28.34 29.23 25.24 25.75 27.26 27.42 27.57 28.02 29.05 26.03 25.86 26.60 28.64 27.14 27.08 27.81 28.94 FISTA [2] SA-DCT [18] BM3D [11] CSR 30.36 29.85 29.01 30.75 29.36 30.28 30.63 31.40 27.20 √标准差为1.6的高斯模糊，σn = 2 29.65 30.84 30.71 32.31 31.23 32.46 32.22 33.44 29.47 30.43 30.69 31.23 25.03 27.00 28.19 27.81 26.80 27.44 27.46 28.24 29.42 29.22 29.03 30.17 29.33 29.70 29.67 31.44 28.97 29.69 29.74 30.70 FISTA [2] SA-DCT [18] BM3D [11] CSR 29.67 29.42 28.56 30.14 27.89 29.88 30.21 31.19 标准差为1.6的高斯模糊，, σn = 2 25.94 26.99 27.08 27.81 29.18 30.04 30.23 31.47 30.74 31.79 31.72 32.60 28.00 29.96 30.28 30.94 24.54 26.08 27.02 26.53 27.94 28.90 28.73 30.03 28.62 29.16 29.10 30.56 28.07 29.13 29.21 30.09 图2. 海星图像（9×9的均匀模糊，n =2）的图像去模糊算法性能比较。（a）噪声与模糊；（b）FISTA【2】（PSNR=27.75dB）；（c）BD3M【11】（PSNR=28.61dB）；（d）CSR（PSNR=30.30dB）。 我们也将提出的CSR模型应用到一些真实运动模糊图像，但真正的模糊内核是不知道的。由于模糊内核估计是一项重要的任务并且超出了本文的讨论范围，所以我们借用了【17】中的模糊内核估计方法来估计模糊内核。然后将估计的模糊内核送入提出的CSR方法。在图3中，通过【17】中的盲目去模糊方法和提出的CSR，我们给出了一个真实模糊图像的去模糊后的图像。我们可以看出，通过我们的方法恢复出的图像更加清晰并且有更多的细节。 图3. 用【17】中的内核估计方法估计出的模糊内核对一幅真实模糊图像去模糊性能的比较。（a）输入的模糊图像；（b）【17】方法的去模糊图像；（c）CSR方法的去模糊图像；（d）特写镜头。 4.2 超分辨率单一图像 在超分辨率单一图像中，观测到的低分辨率（LR）图像是由第一个模糊内核模糊和一个比例因子采样得到的。因此，从单一的LR图像恢复出高分辨率（HR）图像比图像去模糊更加不确定。在本章节中，我们用提出的CSR方法和其他与之相竞争的方法来进行超分辨率单一图像实验。观察到的LR图像是由一个模糊内核模糊一个HR图像，比如说一个标准差为1.6的7×7的高斯滤波器，然后在水平和垂直方向对模糊图像进行比例因子为3的下采样。LR图像也添加了标准差为5的加性高斯噪声，使得图像恢复问题更加有挑战性。由于人的视觉系统对亮度变化更为敏感，所以我们只对亮度成分进行图像恢复方法并且对色差成分使用简单的双三次内插器。 我们对提出的基于CSR的IR方法和一些最近开发的图像超分辨率方法进行比较，例如【12】中的softcut方法、【22】中的基于TV的方法和【27】中的基于稀疏表示的方法。由于【27】中的基于稀疏表示的方法不能同时提高分辨率和去模糊，所以正如【27】中所建议的，我们用迭代反投影来对方法【27】的输出进行去模糊。 表2是几种方法对一组9幅自然图像的PSNR结果报告。从表2中我们可以得出结论：本文提出的CSR方法明显优于其他方法。这显示出CSR模型对于解决图像逆问题的优越性。 图4显示的是CSR和其他方法的主观比较。我们可以看出【22】中的基于TV的方法容易产生分段不变结构；【12】中的softcut方法产生过于平滑的图像局部结构；【27】中的基于稀疏表示方法重建的图像边缘包含一些可见的伪造。显然，通过CSR的图像重建提供了最佳的视觉质量。重建的边缘相比其他三种方法都更加清晰，并且恢复出更多的图像精细结构。 表2 重建HR图像PSNR(dB)结果（亮度成分） Images TV [22] Softcut [12] Sparsity [27] Proposed CSR Butterfly 26.56 24.74 24.70 28.19 flower 27.51 27.31 27.87 29.54 Girl 31.24 31.82 32.87 33.68 Pathenon 26.00 25.95 26.27 27.23 Noiseless Parrot 27.85 27.99 28.70 30.68 Noisy TV [22] Softcut [12] Sparsity [27] Proposed CSR 25.49 24.53 23.61 26.84 26.57 26.98 26.60 28.07 29.86 31.30 30.71 32.03 25.35 25.72 25.40 26.40 27.01 27.69 27.15 29.47 26.74 27.48 27.22 28.03 23.11 22.91 22.45 23.78 28.13 29.13 28.31 29.96 29.70 30.57 29.57 31.73 26.88 27.37 26.78 28.48 Raccoon 27.54 27.82 28.51 29.29 Bike 23.66 23.15 23.23 24.72 Hat 29.20 29.50 29.63 31.33 Plants 31.34 31.19 31.55 34.00 Average 27.88 27.72 28.15 29.85 图4. 植物图片（比例因子为3，n=0）的超分辨率性能比较。（a）LR图像；（b）基于TV的方法【22】（PSNR=31.34dB）；（c）基于稀疏性方法【27】（PSNR=31.55dB）；（d）CSR方法（PSNR=34.00dB）。 5、总结 图像恢复（IR）在图像处理和计算机视觉应用中是一个基本的主题，并且已经被广泛研究。在本文中，我们研究了基于稀疏编码技术的IR。为了更好地了解IR稀疏编码的有效性，我们引入了稀疏编码噪声（SCN）的概念，并且根据经验发现SCN服从拉普拉斯算子分布。为了抑制SCN从而改善IR质量，利用图像非局部相似性提出了集中稀疏表示（CSR）模型。 除了局部稀疏，我们也执行稀疏系数来获得小的SCN，即接近它们的分布中心。局部和非 局部冗余引起的稀疏，都由l1 范数表示，统一成一个变分公式。然后一个CSR模型的贝叶斯解释可以精确地确定正则化参数。IR实验结果表明，CSR图像恢复方法显著优于其他IR方法。 6、致谢 部分支持这项研究工作的有国家自然科学基金（No.61033004，60736043，61070138,61072104），中国中央高校基础研究基金（K50510020003），香港研究资助局研究基金（香港理工大学5375/09E）。 7、参考文献 [1] M. Aharon, M. Elad, and A. Bruckstein. K-svd: an al- gorithm for designing overcomplete dictionaries for sparse representation. IEEE Trans. Signal Process., 54(11):4311– 4322, Nov. 2006. 1, 4 [2] A. Beck and M. Teboulle. Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems. IEEE Trans. On Image Process., 18(11):2419– 2434, Nov. 2009. 5, 6, 7 [3] D. Bertsekas, A. Nedic, and A. Ozdaglar. Convex Analysis and Optimization. Athena Scientific, 2003. 5 [4] J. Biemond, R. Lagendijk, and R. Mersereau. Iterative methods for image deblurring. Proceedings of the IEEE, 78(5):856–883, 1990. 1 [5] A. M. Bruckstein, D. L. Donoho, and M. Elad. From sparse solutions of systems of equations to sparse modeling of sig- nals and images. SIAM Review, 51(1):34–81, Feb. 2009. 1 [6] A. Buades, B. Coll, and J. M. Morel. A non-local algorithm for image denoising. In IEEE Int. Conf. on CVPR, volume 2, pages 60–65, 2005. 1, 2, 3 [7] E. Candes, M. B. Wakin, and S. P. Boyd. Enhancing sparsity by reweighted l1 minimization. Journal of Fourier Analysis and Applications, 14:877–905, Dec. 2008. 4 [8] P. Chatterjee and P. Milanfar. Clustering-based denoising with locally learned dictionaries. IEEE Trans. Image Pro- cessing, 18(7):1438–1451, July 2009. 5 [9] S. Chen, D. Donoho, and M. Saunders. Atomic decomposi- tions by basis pursuit. SIAM Review, 43:129–159, 2001. 1,2 [10] K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian. Image denoising by sparse 3-d transform-domain collaborative fil- tering. IEEE Trans. Image Process., 16(8):2080–2095, Aug. 2007. 2, 5 [11] K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian. Image restoration by sparse 3d transform-domain collaborative fil- tering. In SPIE Conference Series, 2008. 5, 6, 7 [12] S. Dai, M. Han, W. Xu, Y. Wu, Y. Gong, and A. K. Katsagge- los. Softcuts: a soft edge smoothness prior for color image super-resolution. IEEE Trans. Image Process., 18(5):969– 981, May 2009. 1, 6, 8 [13] I. Daubechies, M. Defriese, and C. DeMol. An iterative thresholding algorithm for linear inverse problems with a sparsity constraint. Commun. Pure Appl. Math., 57:1413– 1457, 2004. 1, 2, 5 [14] W. Dong, G. Shi, L. Zhang, and X. Wu. Super-resolution with nonlocal regularized sparse representation. In Proc. of SPIE Vis. Comm. and Image Process. (VCIP), volume 7744, pages 77440H–1–77440H–10, 2010. 1, 5 [15] W. Dong, L. Zhang, G. Shi, and X. Wu. Image deblur- ring and super-resolution by adaptive sparse domain selec- tion and adaptive regularization. IEEE Trans. Image Pro- cessing, 20(7):1838–1857, July 2011. 1, 5 [16] M. Elad and M. Aharon. Image denoising via sparse and redundant representations over learned dictionaries. IEEE Trans. Image Process., 15(12):3736–3745, Dec. 2006. 1, 2 [17] R. Fergus, B. Singh, A. Hertzmann, S. T. Roweis, and W. T. Freeman. Removing camerashake from a single image. ACM Trans. Graph. (SIGGRAPH), pages 787–794, 2006. 5, 6 [18] A. Foi, V. Katkovnik, and K. Egiazarian. Pointwise shape- adaptive dct for high-quality denoising and deblocking of grayscale and color images. IEEE Trans. Image Process., 16(5):1395–1411, May 2007. 5, 7 [19] J. Mairal, F. Bach, J. Ponce, and G. Sapiro. Online dictionary learning for sparse coding. In Proc. of the 26th Annual Int.Conf. on Machine Learning, pages 689–696, 2009. 1 [20] J. Mairal, F. Bach, J. Ponce, G. Sapiro, and A. Zisserman. Non-local sparse models for image restoration. In Proc. of the IEEE ICCV, Tokyo, Japan, 2009. 1, 2, 5 [21] J. Mairal, M. Elad, and G. Sapiro. Sparse representation for color image restoration. IEEE Trans. on Image Processing, 17(1):53–69, Jan. 2008. 1 [22] A. Marquina and S. J. Osher. Image super-resolution by tv-regularization and bregman iteration. J. Sci. Comput., 37:367–382, 2008. 6, 7, 8 [23] B. Olshausen and D. Field. Emergence of simple-cell re- ceptive field properties by learning a sparse code for natural images. Nature, 381:607–609, 1996. 1 [24] B. Olshausen and D. Field. Sparse coding with an overcom- plete basis set: a strategy employed by v1? Vision Research, 37:3311–3325, 1997. 1 [25] L. Rudin, S. Osher, and E. Fatemi. Nonlinear total variation based noise removal algorithms. Phys. D, 60:259–268, 1992. 1 [26] L. Sendur and I. W. Selesnick. Bivariate shrinkage func- tions for wavelet-based denoising exploiting interscale de- pendency. IEEE Trans. Signal Process., 50(11):2744–2756,Nov. 2002. 4 [27] J. Yang, J. Wright, Y. Ma, and T. Huang. Image super- resolution as sparse representation of raw image patches. In Proc. of the IEEE CVPR, Jun. 2008. 1, 6, 7, 8 [28] L. Zhang, W. Dong, D. Zhang, and G. Shi. Two-stage image denoising by principal component analysis with local pixel grouping. Pattern Recognition, 43:1531–1549, Apr. 2010. 5 代码及注释 //添加噪声 function im_o = Add_noise( im, v ) seed = 0; randn( 'state', seed ); noise = randn( size(im) ); im_o = double(im) + v*noise; //模糊函数 function z = Blur (mode,x,psf) ws = size(psf); t = (ws-1)/2; if mode == 'fwd', //正变换模式 s = x; //输入图像 se = [s(:,t:-1:1,:), s, s(:,end:-1:end-t+1,:)]; se = [se(t:-1:1,:,:); se; se(end:-1:end-t+1,:,:)]; if size(x,3)==3 //输出z = Hx z(:,:,1) = conv2(se(:,:,1),psf,'valid'); z(:,:,2) = conv2(se(:,:,2),psf,'valid'); z(:,:,3) = conv2(se(:,:,3),psf,'valid'); else z = conv2(se,psf,'valid'); end elseif mode == 'trn', //伴随变换模式 y = x; ye = [y(:,t:-1:1,:), y, y(:,end:-1:end-t+1,:)]; ye = [ye(t:-1:1,:,:); ye; ye(end:-1:end-t+1,:,:)]; if size(ye,3)==3 //输出z = H’x z(:,:,1) = conv2(ye(:,:,1),psf,'valid'); z(:,:,2) = conv2(ye(:,:,2),psf,'valid'); z(:,:,3) = conv2(ye(:,:,3),psf,'valid'); else z = conv2(ye,psf,'valid'); end end //图像去模糊 clc; clear; addpath('Utilities'); blur_type = 1; //统一高斯模糊内核 if blur_type == 1 //blur_type = 1, blur_par表示内核大小 blur_type = 2, blur_par表示高斯内核标准差 blur_par = 9; //统一模糊内核大小默认为9 else blur_par = 1.6; //高斯模糊内核的标准偏差默认为3 end nSig = sqrt(2); //加性高斯噪声的标准差 out_dir = 'Results\\Deblurring_results\\'; im_dir = 'Data\\Deblurring_test_images\\'; im_name = 'Butterfly.tif'; if blur_type==1 if nSig<2 par.t1 = 0.06*nSig^2; par.c1 = 0.07*nSig^2; par.t2 = 0.19*nSig^2; par.c2 = 0.56*nSig^2; par.eps2 = 0.31; else par.t1 = 0.04*nSig^2; par.c1 = 0.06*nSig^2; par.t2 = 0.17*nSig^2; par.c2 = 0.45*nSig^2; par.eps2 = 0.32; end else par.t1 = 0.03*nSig^2; par.c1 = 0.04*nSig^2; par.t2 = 0.12*nSig^2; par.c2 = 0.36*nSig^2; par.eps2 = 0.31; end par.nSig = nSig; //高斯噪声的标准差 par.iters = 160; par.nblk = 12; par.sigma = 1.4; par.eps = 0.3; par.I = double( imread(fullfile(im_dir, im_name)) ); [par.bim par.fft_h] = Generate_blur_image(par.I, blur_type, blur_par, par.nSig); [im PSNR SSIM] = CSR_Deblurring( par ); fname = strcat('CSR_', im_name); imwrite(im./255, fullfile(out_dir, fname)); disp( sprintf('%s: PSNR = %3.2f SSIM = %f\\n\\n', im_name, PSNR, SSIM) ); //CSR方法图像去模糊 function [im_out PSNR SSIM] = CSR_Deblurring( par ) time0 = clock; bim = par.bim; [h w ch] = size(bim); par.step = 2; par.win = 6; par.h = h; par.w = w; par.cls_num = 64; par.hp = 80; par.s1 = 22; dim = uint8(zeros(h, w, ch)); ori_im = zeros(h,w); if ch == 3 b_im = rgb2ycbcr( uint8(bim) ); dim(:,:,2) = b_im(:,:,2); dim(:,:,3) = b_im(:,:,3); b_im = double( b_im(:,:,1)); if isfield(par, 'I') ori_im = rgb2ycbcr( uint8(par.I) ); ori_im = double( ori_im(:,:,1)); end else b_im = bim; if isfield(par, 'I') ori_im = par.I; end end disp(sprintf('The PSNR of the blurred image = %f \\n', csnr(b_im(1:h,1:w), ori_im, 0, 0) )); d_im = Deblurring(b_im, par, ori_im, b_im, 0, 4); d_im = Deblurring(b_im, par, ori_im, d_im, 1, 2); if isfield(par,'I') [h w ch] = size(par.I); PSNR = csnr( d_im(1:h,1:w), ori_im, 0, 0 ); SSIM = cal_ssim( d_im(1:h,1:w), ori_im, 0, 0 ); end if ch==3 dim(:,:,1) = uint8(d_im); im_out = double(ycbcr2rgb( dim )); else im_out = d_im; end disp(sprintf('Total elapsed time = %f min\\n', (etime(clock,time0)/60) )); return; //图像去模糊函数 function d_im = Deblurring(b_im, par, ori_im, d_im0, flag, K) d_im = b_im; [h1 w1] = size(ori_im); b2 = par.win*par.win; lam = zeros(0); gam = zeros(0); fft_h = par.fft_h; Y_f = fft2(b_im); A_f = conj(fft_h).*Y_f; H2_f = abs(fft_h).^2; cnt = 0; if flag==1 d_im = d_im0; end for k = 1:K Dict = KMeans_PCA( d_im, par, par.cls_num ); Dict.D0 = dctmtx(b2); [blk_arr wei_arr] = Block_matching( d_im, par); if flag==1 [lam, gam] = Sparsity_estimation( d_im, par, Dict, blk_arr, wei_arr ); end Reg = @(x)CSR_Regularization(x, par, Dict, blk_arr, wei_arr, lam, gam, flag ); f = d_im; for iter = 1 : par.iters cnt = cnt + 1; f_pre = f; if (mod(cnt, 40) == 0) if isfield(par,'I') //去模糊图像的PSNR(dB)值 PSNR = csnr( f(1:h1,1:w1), ori_im, 0, 0 ); fprintf( 'CSR deblurring: iter. %d : PSNR = %f\\n', cnt, PSNR ); end end for i = 1 : 3 im_f = fft2((f_pre)); Z_f = im_f + (A_f - H2_f.*im_f)./(H2_f + par.eps2); z = real(ifft2((Z_f))); f_pre = max(min(z,255),0); end f = Reg( f_pre ); end d_im = f; end //稀疏估计函数 function [lam, gam] = Sparsity_estimation( im, par, Dict, blk_arr, wei_arr ) b = par.win; s = par.step; b2 = b*b; [h w] = size(im); PCA_idx = Dict.cls_idx; PCA_D = Dict.PCA_D; s_idx = Dict.s_idx; seg = Dict.seg; A = Dict.D0; N = h-b+1; M = w-b+1; r = [1:s:N]; r = [r r(end)+1:N]; c = [1:s:M]; c = [c c(end)+1:M]; X0 = zeros(b*b, N*M); X_m = zeros(b*b,length(r)*length(c),'single'); N = length(r); M = length(c); L = N*M; k = 0; for i = 1:b for j = 1:b k = k+1; blk = im(i:end-b+i,j:end-b+j); X0(k,:) = blk(:)'; end end idx = s_idx(seg(1)+1:seg(2)); set = 1:size(X_m,2); set(idx) = []; for k = 1 : L i = set(k); cls = PCA_idx(i); P = reshape(PCA_D(:, cls), b2, b2); coe = P*X0(:, blk_arr(i,:)); v0(:,i) = mean(coe.^2, 2); coe = P*(X0(:, blk_arr(i,:)) - repmat(X_m(:, i), 1, par.nblk) ); vu0(:,i) = mean(coe.^2, 2); end v0 = max(0, v0-par.nSig^2); vu0 = max(0, vu0-par.nSig^2); lam = (par.c1*sqrt(2)*par.nSig^2)./(sqrt(v0) + par.eps); gam = (par.c2*sqrt(2)*par.nSig^2)./(sqrt(vu0) + par.eps); return; 运行结果 去模糊前 去模糊后 普通去模糊方法 去模糊前 去模糊后 CSR方法 因篇幅问题不能全部显示，请点此查看更多更全内容