出投笔记 -
图像复原技术及其MATLAB实现
摘 要
图像复原的目的是从观测到的退化图像重建原始图像，它是图像处理、模式识别、机器视觉等的基础，在天文学、遥感成像、医疗图像等领域获得了重要应用。运动模糊图像的复原是图像复原的重要组成部分。
由运动模糊图像复原出原图像关键问题是获取点扩展函数,模糊方向和长度的鉴别至关重要。本文通过对运动模糊图像的频域幅度图的黑带条纹(即图像零点个数)分析,计算出运动模糊PSF的参数。
获得PSF的参数后，本文主要采用了逆滤波法、维纳滤波法、最小线性二乘法、Richardson-Lucy算法对模糊图像进行复原，并对各种复原方法的结果进行了分析与对比。
关键词：图像复原；运动模糊；模糊方向；模糊长度
引 言
图像复原是图像处理领域一个具有现实意义的课题。运动模糊图像的研究越来越受到关注，这种模糊是成像过程中普遍存在的问题，其复原在许多领域都有广泛的应用。
实际上，图像复原设计三个方面的内容：退化图像的成像模型，图像复原算法和复原图像的评价标准。不同的成像模型、问题空间、优化规则和方法都会导致不同的图像复原算法，适用于不同的应用领域。现有的复原方法概括为以下几个类型：去卷积复原算法、线性代数复原、图像盲反卷积算法等，其他复原方法多是这三类的衍生和改进。其中，去卷积方法包括维纳去卷积、功率谱平衡与几何平均值滤波等，这些方法都是非常经典的图像复原方法。但是需要有关于原始图像、降质算子较多的先验信息和噪声平衡性的假设只适合于不变系统及噪声于信号不相关的情形，特别是降质算子病态的情况下，图像复原结果还不太理想。
由于图像复原技术在图像处理中占有重要的地位，已经形成了一些经典的常用图像复原算法，如无约束最小二乘法、有约束最小二乘方法、逆滤波、维纳、最大熵复原等，至今还被广泛使用。但这些复原算法都是假设系统的点扩散函数PSF(即系统对图像中点的脉冲响应，是导致图像退化的原因)为已知，实际情况是系统的点扩散函数由于大气扰动、光学系统的相差、相机和对象之间的相对运动等多种因素的影响，往往是未知的。这就需要人们用某种先验知识在系统的点扩散函数未知时进行估计，然而这种先验知识并不容易取得也不够精确，这就需要我们在对己模糊图像分析和处理的基础之上估计最逼近的PSF。
在运动模糊方向的鉴别方面，由于匀速直线运动的点扩散函数是矩形函数，其模糊图像对应的频域上有周期性的零值条纹，运动方向与零值条纹方向相垂直，本文就是借用此法获取模糊图像的PSF参数。本文主要针对运动模糊图像的复原进行研究，讨论分析了匀速直线运动模糊的退化模型，研究了运动方向和模糊尺度的估计，介绍了常用的几种图像复原方法。对模糊图像用几种复原算法分别进行了复原，根据复原结果，讨论分析了各算法的优缺点及适用的恢复环境。
第1章 绪论
1.1 研究背景
图像复原是数字图像处理中的一个重要课题。它的主要目的是改善给定的图像质量并尽可能恢复原图像。图像在形成、传输和记录过程中，受多种因素的影响，图像的质量都会有不同程度的下降，典型的表现有图像模糊、失真、有噪声等，这一质量下降的过程称为图像的退化。图像复原的目的就是尽可能恢复被退化图像的本来面目。
在成像系统中，引起图像退化的原因很多。例如，成像系统的散焦，成像设备与物体的相对运动，成像器材的固有缺陷以及外部干扰等。成像目标物体的运动，在摄像后所形成的运动模糊。当人们拍摄照片时，由于手持照相机的抖动，结果像片上的景物是一个模糊的图像。由于成像系统的光散射而导致图像的模糊。又如传感器特性的非线性，光学系统的像差，以致在成像后与原来景物发生了不一致的现象，称为畸变。再加上多种环境因素，在成像后造成噪声干扰。人类的视觉系统对于噪声的敏感程度要高于听觉系统，在声音传播中的噪声虽然降低了质量，但时常是感觉不到的。但景物图像的噪声即使很小都很容易被敏锐的视觉系统所感知。图像复原的过程就是为了还原图像的本来面目，即由退化了的图像恢复到能够真实反映景物的图像。
在交通系统、刑事取证中图像的关键信息至关重要，但是在交通、公安、银行、医学、工业监视、军事侦察和日常生活中常常由于摄像设备的光学系统的失真、调焦不准或相对运动等造成图像的模糊，使得信息的提取变得困难。但是相对于散焦模糊，运动模糊图像的复原在日常生活中更为普遍，比如高速运动的违规车辆的车牌辨识，快速运动的人群中识别出嫌疑人、公安刑事影像资料中提取证明或进行技术鉴定等等，这些日常生活中的重要应用都需要通过运动模糊图像复原技术来尽可能地去除失真，恢复图像的原来面目。因此对于运动模糊图像的复原技术研究更具有重要的现实意义。
1.2 相关领域的研究现状及存在的问题
图像恢复是数字图像处理中的一个重要分支，它研究的是如何从所得的退化图像中以最大的保真度复原出真实图像。成像系统的缺陷，传播媒介中的杂质，以及图像记录装置与目标之间的相对运动等因素，都不可避免地造成了图像的某些失真和不同程度的降质。然而在众多的应用领域中，又需要清晰的、高质量的图像，因此，图像恢复问题具有重要的意义。
与图像增强相似，图像复原的目的也是改善图像的质量。图像复原可以看作图像退化的逆过程，是将图像退化的过程加以估计，建立退化的数学模型后，补偿退化过程造成的失真，以便获得未经干扰退化的原始图像或图像的最优估计值，从而改善图像质量。图像复原是建立在退化的数学模型基础上的，且图像复原是寻求在一定优化准则下的原始图像的最优估计，因此，不同的优化准则会获得不同的图像复原，图像复原结果的好坏通常是按照一个规定的客观准则来评价的。
运动模糊图像的恢复是图像恢复中的重要课题之一，随着科学技术的不断发展，它在各个领域中的应用越来越多，要求也越来越高，可广泛应用于天文、军事、道路交通、医学图像、工业控制及侦破领域，具有重要的现实意义。
图像复原作为图像处理的一个重要领域，对于该问题国内外展开了诸多关键技术的研究。实际上，图像复原涉及三个方面的内容：退化图像的成像模型，图像复原算法和复原图像的评价标准。不同的成像模型、问题空间、优化规则和方法都会导致不同的图像复原算法。适用于不同的应用领域。现有的复原方法概括为以下几个类型：去卷积复原算法、线性代数复原、图像盲反卷积算法等，其他复原方法多是这三类的衍生和改进。其中，去卷积方法包括维纳去卷积、功率谱平衡与几何平均值滤波等，这些方法都是非常经典的图像复原方法。但是需要有关于原始图像、降质算子较多的先验信息和噪声平衡性的假设只适合于不变系统及噪声于信号不相关的情形，特别是降质算子病态的情况下，图像复原结果还不太理想。
线性代数复原技术是基于已知降质算子和噪声的统计特征，从而利用线性代数原理的复原技术，它为复原滤波器的数值提供了一个统一的设计思路和较透彻的解释。但是当降质函数有接近零的特征值时，复原的结果对噪声特别敏感，且该方法是把整幅图像一并处理，计算量大，同时也没有考虑纹理、边界等高频信号与噪声的区别，这将使纹理、边界等重要特征在图像复原过程中被破坏。针对这些问题，国外主要在改进算法的效率上做了许多工作，如全局最小二乘法、约束总体最小二乘法和正则化约束总体最小二乘法。
图像盲反卷积是图像复原的另一个重要的手段，它针对没有或少有关于降质函数和真实信号灯先验知识的复原问题，直接根据退化图像来估计降质函数和真实信号。目前有以下几种算法：零叶面分离法、预先确定降质函数法、三次相关法、迭代盲反卷积法等。这些算法在先验信息不足的情况下对降质图像进行复原，由于原始图像以及点扩展函的先验知识只是部分已知的，造成图像复原的解往往不唯一，而且解的好坏与初始条件的选择以及附加的图像假设等直接有关。同时，由于加性噪声的影响使得图像的盲目复原成病态。即若对点扩展函数直接求逆进行复原，通常会带来高频噪声放大的问题导致算法性能的恶化，所以当图像的信噪比水平较低时获得的结果往往不太理想。
正则化方法作为一种解决病态反问题的常用方法，通常用图像的平滑性作为约束条件，但是这种正则化策略通常导致复原图像的边缘模糊。为了克服边缘退化问题，最近几年，不少学者对各种“边缘保持”的正则化方法进行了比较深入的研究，提出了一些减少边缘退化的正则化策略，这些策略通常需要引入非二次正则化泛函，从而使问题的求解成为一个非线性问题。沿着这一思路，Geman和Yang提出了“半二次正则化”的概念来解决这种策略中出现的非线性优化问题。其后，Charbonni等人在此基础上研究了一种新的半二次正则化方法。从而可以利用确定性算法来得到问题的最优解。另一个较新的发展使Vogel等人提出的基于全变差模型的图像复原算法。尽管这些算法都在一定意义上提高了复原图像的质量，但边缘模糊的问题并未得到理想的解决。
另外，近年来小波的理论得到迅速发展，并光法应用于图像复原中。基于小波变换的迭代正则化图像复原算法，兼顾抑制噪声的增长和保留图像的重要边界。具有噪声估计能力的图像恢复正则化方法。Belge等人以广义高斯模型作为小波系数的先验分布，提出了一种小波域边缘保持正则化的方法。同时给出了小波域图像复原的一般框架，但其复原方法相对于传统复原方法提高的并不显著，赵书斌等人以混合高斯模型逼近小波系数的分布，并引入小波域隐马尔可夫模型作为自然图像的先验概率模型对图像超分辨率复原问题进行正则化，复原效果不错，但该方法还是不能避免计算量过大的缺点。从图像复原的Bayesian框架出发，小波域局部高斯模型的线性图像复原方法，该方法较好的再现了图像的各种边缘信息，取得不错的复原效果。
2 逆滤波复原
图5-1逆滤波复原过程图，图(a)是选取的原始图像，图(b)是利用MATLAB对原始图像进行运动模糊和加噪声仿真而生成的仿真图像，模糊长度为10个像素。经过逆滤波复原图像为图（c）。
（a）原图像 （b）模糊加噪图像
（c）复原图像
图5-1 逆滤波复原过程
从恢复出来的图5-1（c）可以看出复原效果不理想，出现较大的振铃现象。从理论分析上看是由于一般情况下，逆滤波复原不能正确估计退化函数
的零点，因此必须采用一个折中的方法进行解决。实际上，逆滤波不用
，而是采用另外一个关于
的函数
。函数
称为转移函数。改进方法为
(5-1)
式中K和d均为小于1的常数。
采用逆滤波对运动模糊加噪图像进行复原。在噪声相同情况下，参数K分别选取0.0001、0.001、0.005、0.01、0.1和0.9。图中（a）-（f）为对应参数下的复原图像。通过转移函数替代原退化模型得到图5-2逆滤波减小振铃现象复原图像。
（a）k=0.0001 (b) k=0.001
（c）k=0.005 (d) k=0.01
(e) k=0.1 (f) k=0.9
图5-2 逆滤波减小振铃现象复原图像图
从复原结果图5-2可以看出随着K值的增大复原效果逐渐变差，K=0.001到0.01之间复原效果较好。
从理论上分析，逆滤波方法不能正确估计退化函数
的零点。采用一个折衷的方法加以解决。一般情况，可以将图像的退化过程视为一个具有一定带宽的带通滤波器，随着频率的升高，该滤波器的带通特性很快下降，即
的幅度随着
平面原点的距离的增加而迅速下降，而噪声项
的幅度变化是比较平缓的。在远离平面的
的值就会变得很大，而对于大多数图像来说，
却变小，在这种情况下，噪声反而占优势，自然无法满意的恢复出原图像。这一规律说明，应用逆滤波时仅在原点领域采用
方能有效。
5.3 有约束最小二乘方复原的实现
通过MATLAB仿真来实现有约束的最小二乘方复原，图5-3是有约束的最小二乘方复原图。分别取参数为0、1、0.001、0.01、0.05、0.1对应图5-3里面的(a)-(f)。
图5-3 有约束最小二乘方在不同参数下的恢复情况图
5.4 维纳滤波复原的实现
图(a)是选取的原始图像，图(b)是利用MATLAB对原始图像进行运动模糊和加噪声仿真而生成的仿真图像，模糊长度为5个像素。
采用维纳滤波恢复算法对模糊图像进行恢复，在加噪声的情况下，参数k分别选取0．0001、0．001、0．005、0．01、0.1和1。各种图中(c)-(h)为对应参数下的恢复图像。图5-4有噪声下维纳滤波在不同参数下的恢复情况。
（a）原始图像 （b）含噪声运动模糊图像d=5，v=0.001
(c) K=0.0001 (d) K=0.001
(e) K=0.005 (f) K=0.01
(g) K=0.1 (h) K=1
图5-4 有噪声下维纳滤波在不同参数下的恢复情况
可以看出，恢复图像还是都有一定的振铃现象。K=0.0001时，图像振铃效应比较小，但其噪声很大。k=0.00l相对前一幅恢复图像振铃效应明显一点，但噪声有所减少。k=0.005和k=0.01的恢复效果也是看起来区别不明显，虽然它们的噪声都减少了，但图像整体都相对前面有明显的模糊，且振铃效应明显。K=0.1的模糊程度比较大，而k=1时，图像最模糊，且亮度很暗。
总的看来，主观评价认为k=0.001时的恢复效果最清楚，恢复质量最好。K=0.0001时次之，k=1时的效果最差。
利用公式（5-1）与（5-2）计算出平均平方误差，如表5-1所示：
表5- 1平均平方误差客观评价方法得分
参数k 0.0001 0.001 0.005 0.01 0.1 1
平均平方误差（M） 2861 6463 6937 6576 2861 2513
从表5-1可以看出，采用平均平方误差准则时K=0.1的平均平方误差和K=0.0001一样，但是其对应的图像很模糊。对于非迭代方法的维纳滤波恢复法，k值的选取对图像恢复质量有很重要的影响。从上面实验结果可以看出，虽然对于每幅特定图像的评价得分不尽相同，但基本上当k值在[0.001，0.1]的范围取值时，恢复图像质量最好。
从理论上分析，维纳滤波复原的图像，在图像的频率特征和附加噪声已知的情况下，采用维纳滤波去卷积比较有效。维纳滤波复原法不存在极点，即当
很小或变为零时，分母至少为K，而且
的零点也转换成了维纳滤波器的零点，抑制了噪声，所以它在一定程度上克服了逆滤波复原方法的缺点。
5.5 Richardson-Lucy复原的实现
图(a)是原始图像，图(b)是对原图进行运动模糊仿真而生成的仿真图像，模糊长度为10个像素，模糊方向为水平方向。
采用Richardson-Lucy恢复算法对模糊图像进行恢复，迭代次数参数分别选取20次、50次、100次、150次、200次和300次。所有图的(c-h)为对应迭代次数下的复原图像。
(a) 原图像 (b) 水平运动10像素加噪声图像
(c) 迭代20次 (d) 迭代50次
(e) 迭代100次 (f) 迭代150次
(g) 迭代200次 (h) 迭代300次
图5-5 R-L算法在不同参数下的复原图像
从图5-5可以看出，恢复的图像整体差别不大。从图像人物、背景等分辨，质量随着迭代次数增大而提高。迭代100次以后恢复效果区别不大，仔细辨认，迭代200次和300次更好一些．所以主观评价认为200次或300次时复原质量最好。
本文通过MATLAB编程，利用公式（5-1）、（5-2）计算出图5-5：（c）-(h)各恢复图像的平均平方误差。通过计算平均平方误差的倒数（M）来做客观分析。客观分析如表5-2：
表5- 2复原图像与原图像的平均平方误差
迭代次数 20 50 100 150 200 300
平均平方误差（M） 4348 4762 5000 5263 5263 5263
从表5-2可以看出，采用均方误差准则评价时，平均平方误差差别不大。总体上随着迭代次数增加分数增大，迭代次数为150次后平均平方误差一样。
上述分析表明，在R-L恢复算法下，对于这幅图像，传播波方程客观评价准则和平均平方误差准则的评价结果基本一致，并和主观评价结果吻合。
从理论上分析，Richardson-Lucy算法能够按照泊松噪声统计标准求出与给定PSF卷积后，最有可能成为输入模糊图像的图像。PSF已知但是图像噪声信息未知时，也可以使用这个函数进行有效的工作。随着复原迭代的次数增加，可以提高复原图像的似然性，最终将会收敛在具有最大似然性处。
结论与展望
图像复原是图像处理领域一个具有现实意义的课题。运动模糊图像的研究越来越受到关注，这种模糊是成像过程中普遍存在的问题，其复原在许多领域都有广泛的应用。
图像复原需要根据相应的退化模型知识重建或恢复原始的图像。也就是说，图像复原技术就是要将图像退化的过程模型化，并由此采取相反的过程以得到原始图像。运动模糊是由于在拍摄过程中相机与景物之间相对运动而产生，因此对于匀速直线运动造成的运动模糊图像来说，图像退化模型的两个重要参数相对运动的方向和运动模糊尺度的估计就成了图像复原的关键问题。本文以匀速直线运动造成的模糊图像为基础，研究退化函数的参数估计方法，所做的工作及创新之处总结如下：
论文的工作总结
（1）论文研究了模糊图像尤其是水平方向运动模糊图像的退化模型，任意方向的匀速直线运动模糊图像只需要通过坐标旋转至水平方向，其图像特征的描述可由水平匀速直线运动模糊图像类推得出。
（2）论文研究了运动方向和模糊尺度的估计，通过对运动模糊图像的频域幅度图的黑带条纹(即图像零点个数)分析来估算出运动模糊PSF的参数运动模糊方向和运动模糊长度)的，同时通过查阅文献获得另一种对模糊尺度的估算即对模糊图像进行一阶微分，然后进行自相关运算，可得到一条鉴别曲线，曲线上会出现对称的相关峰，峰值为负，两相关峰之间的距离等于运动模糊长度。
（3）对于运动模糊图像的恢复，介绍分析了逆滤波、有约束的最小二乘方、维纳滤波和Richardson-Lucy四种常用的恢复方法。并且采用Richardson-Lucy迭代算法和维纳滤波方法在选取不同参数的情况下对运动模糊图像进行了恢复。利用逆滤波方法进行恢复时，复原图像的效果整体不是很好，存在着较明显的振铃效应，加噪情况下复原图像的噪声也比较严重。本文通过理论分析及仿真，探索出减小振铃现象的一些方法，但还不够完善，今后还需要继续深入研究如何改进算法、减少振铃效应和噪声，以提高复原图像的质量。针对有约束的最小二乘方、维纳滤波复原方法，本文主要通过参数变化来控制复原效果，最终选出最优准则。Richardson-Lucy迭代算法从理论上看是迭代次数越大，复原效果越好。考虑到程序的有效性，本文采用了最大为迭代300次。从主客观评价对复原图像的评价来看迭代次数超过150次以后效果就基本一样。
展望
由于本人的能力有限，对图像复原技术的研究还不够系统、不够深入，无论在理论上，还是在工程应用中，还需要做大量深入、细致的研究工作。因此在这方面的研究还只是个开始，很多地方都需要改进与提高，例如：
（1）运动模糊图像的复原大多是对整幅图像进行全局的复原，然而在实际应用中并非完全如此。例如，由于物体运动而产生的相对运动，其运动模糊只出现在物体运动的轨迹上，而背景是清晰的。在这种情况下就不能对全局进行处理，应首先分割出运动模糊区域，然后再进行参数估计，图像复原。如何分割运动模糊区域，分割的依据如何等将成为以后研究工作的一部分。
（2）本文研究的运动模糊图像参数估计算法仅限于匀速直线运动造成的模糊，而缺乏对非匀速的、轨迹为曲线的运动模糊研究，且得到的参数还具有一定的误差。
参考文献
[1] 杨帆，等.数字图像处理与分析（第二版）[M].北京：航空航天大学出版社，2007.
[2] 黄爱民，等.数字图像处理分析基础[M].北京：中国水利水电出版社，2005.
[3] 孙兆林.MATLAB 6.x图像处理[M].北京:科学出版社，2003.
[4] 贾永红.计算机图像处理与分析[M].武汉：武汉大学出版社，2001.
[5] 姚敏.数字图像处理[M].北京：机械工业出版社，2006.
[6] 孟永定,马佳.基于MATLAB实现数字图像恢复[J].电脑学习,2007,1(1),30-32.
[7] 刘红岩,徐志鹏.基于MATLAB的数字图像恢复[J].科技信息(学术研究).2008,3(12),23-26.
[8] 孟昕,张燕平.运动模糊图像恢复的算法研究与分析[J].计算机技术与发展,2007,17(8):74-76.
[9] 孟昕,周琛琛,郝志廷.运动模糊图像恢复算法相关研究发展概述[J].安徽电子信息职业技术学院学报,2008,7(6),38-41.
[10] 曾志高,谭骏珊.匀速直线运动模糊图像的恢复技术研究[J].陕西理工学院学报(自然科学版),2006,22(2),36-38.
[11] 李云浩,王建设.匀速直线运动模糊图像的退化数学模型试验研究[J].江西理工大学学报,2006,27(4),28-30.
[12] 谢伟,秦前清.基于倒频谱的运动模糊图像PSF参数估计[J].武汉大学学报(信息科学版).2008,5(02),30-32.
[13] Banham Mark M,Katsaggelos A K.Digital image restoration. IEEE Signal Processing Magazine . 1997.
[14] Mccallum B C.Blind deconvolution by simulated annealing. Optics Communitions . 1990.
[15] Hardie C and Boncelet C G.LUM filters: A class rank order based filters for smoothing and sharpening. IEEE Transactions on Signal Processing . 1993.
附录C 主要源程序
(1)模糊图像的傅里叶频谱获取
j=imread("车牌.jpg");
figure(1);
imshow(j);
title("原图像");
len=20;
theta=30;
psf=fspecial("motion",len,theta);
j1=imfilter(j,psf,"circular","conv");
figure,imshow(j1);
title("PSF模糊图像");
J=rgb2gray(j);
K=fft2(J); %傅里叶变换
M=fftshift(K); %直流分量移到频谱中心
N=abs(M); %计算频谱幅值
P=(N-min(min(N)))/(max(max(N))-min(min(N)))*225;%归一化
Figure；
imshow(P);
title("原图像的傅里叶变换频谱");
J1=rgb2gray(j1);
K1=fft2(J1); %傅里叶变换
M1=fftshift(K1); %直流分量移到频谱中心
N1=abs(M1);%计算频谱幅值
P1=(N1-min(min(N1)))/(max(max(N1))-min(min(N1)))*225;%归一化
Figure；
imshow(P1);
title("模糊图像的傅里叶变换频谱");
(2)模糊长度获取程序
f1=rgb2gray(j1);
f1=im2double(f1);
h = fspecial("Sobel"); %Sobel算子
J = conv2(f1,h,"same"); %Sobel算子微分
IP=abs(fft2(J)); %图像能量谱密度
S=fftshift(real(ifft2(IP)));
figure,plot(S);
title