基于MATLAB复原离焦模糊图像
摘 要 图像在获取、传输和存储过程中会受到如模糊、失真、噪声等原因的影响，这些原因会使图像的质量下降。因此，我们需要采取一定的方法尽可能地减少或消除图像质量的下降，恢复图像的本来面目，这称为图像复原。通过阅读图像复原技术相关资料，本文主要探讨了维纳(Wiener)过滤器，限制最小二乘滤波算法，Lucy-Richardson使用相关工具箱函数的算法和盲解卷积算法deconvwnr函数、deconvreg函数、deconvlucy函数、deconvblind模拟函数。此外，本文模拟了上述算法，并分析了四种算法的实验结果。
关键词 图像恢复；维纳滤波恢复；限制最小二乘滤波恢复；Lucy-Richardson恢复；盲解卷积恢复
前 言
在现实生活中，人们应该接触到许多图像和图像。在景物成像过程中，可能会出现模糊、失真或混合噪声，最终导致图像质量下降。这种现象被称为图像退化。因此，我们可以采取一些技术手段，尽量减少甚至消除图像质量的下降，恢复图像的原始面貌，即在预定义的意义上改善给定的图像，即图像恢复。虽然图像增强和图像恢复之间存在重叠部分，但前者主要是主观处理，而图像恢复大多是客观处理。恢复试图通过使用退化现象的先验知识重建或恢复退化图像。因此，恢复技术倾向于模拟退化，并使用相反的处理来恢复原始图像，即考虑使用模糊函数来消除图像的模糊。造成图像模糊的原因有很多，如运动、高斯噪声、斑点噪声、胡椒盐噪声等。
本文主要研究离焦模糊图像的恢复。离焦模糊图像是指拍摄过程中场景和相机的相对运动引起的离焦 ，或成像区域不同深度的对象引起的不同程度的离焦 ，以及自动调焦系统混淆引起的照片离焦。因此，本文的研究应用MATLAB模型退化现象，使用几种常用的滤波方法MATLAB模拟实现为人们在不同的应用和图像数据条件下选择不同的恢复算法提供了一定的依据。
1 图像退化/恢复处理的模型
由于退化是线性的，空间不变的退化函数可以模型化为卷积。同样，恢复有时也称为反卷积。
2 噪声特性
图像噪声是图像在摄入或传输过程中受到的随机信号干扰，是图像中阻碍人们接受信息的各种因素。大多数情况下，图像噪声被视为一个多维随机过程，因此描述噪声的方法可以借用随机过程的描述，即概率分布函数和概率密度分布函数。在许多情况下，这些函数难以测量和描述，甚至无法获得，因此常用的统计特征，如平均值、方差和相关函数。
3噪声的分类
1、按干扰源分类
根据其干扰源，图像噪声可分为外部噪声和内部噪声。
外部噪声，即指系统外部干扰以电磁波或经电源串进系统内部而引起的噪声。如电气设备，天体放电现象等引起的噪声。
内部噪声：一般又可分为以下四种：
- 由光和电的基本性质所引起的噪声。如电流的产生是由电子或空穴粒子的集合，定向运动所形成。因这些粒子运动的随机性而形成的散粒噪声；导体中自由电子的无规则热运动所形成的热噪声；根据光的粒子性，图像是由光量子所传输，而光量子密度随时间和空间变化所形成的光量子噪声等。
- 电器的机械运动产生的噪声。如各种接头因抖动引起电流变化所产生的噪声；磁头、磁带等抖动或一起的抖动等。
- 器材材料本身引起的噪声。如正片和负片的表面颗粒性和磁带磁盘表面缺陷所产生的噪声。随着材料科学的发展，这些噪声有望不断减少，但在目前来讲，还是不可避免的。
- 系统内部设备电路引起的噪声。如电源引入的交流噪声；偏转系统和位移电路引起的噪声。
2.根据噪声和信号之间的关系，
乘性声音与图像信号有关，通常随着图像信号的变化而变化，如飞点扫描图像中的声音、电视扫描光栅、胶片颗粒等。由于载送每个象素信息载体的变化而产生的噪声由信息本身调节。在某些情况下，如果信号变化很小，噪音就不大。为了便于分析和处理，乘性噪声通常被视为加性噪声，并且总是假设信号和噪声是相互统计和独立的。
3. 按概率密度函数分为
这更重要，主要是因为引入数学模型，这有助于使用数学手段去除噪音。
(1)白噪音（White Noise）：常量功率谱。白噪声的一个特例是高斯噪声（Gaussian Noise）。在空间域和频域中，由于高斯噪声在数学上的易处理性，这种噪声（也称为正态噪声）模型在实践中经常使用。事实上，这种易处理性非常方便，因此高斯模型通常适用于街道条件。其方形图曲线服从一维高斯分布：
(2)椒盐噪声（Pepper Noise）：胡椒盐噪声是由图像传感器、传输通道、解码处理等引起的黑白亮点噪声，通常由图像切割引起。胡椒盐噪声是指两种噪声，一种是盐噪声（salt noise），另一种是胡椒噪音（pepper noise）。盐=白色，椒=黑色。前者为高灰度噪声，后者为低灰度噪声。一般来说，两种噪声同时出现，图像上呈现黑白杂点。该噪声在图像中更为明显，对图像分割、边缘检测、特征提取等后续处理具有严重的破坏性。
(3)冲击噪声（Impulsive Noise）：指图像被个别噪声像素破坏，这些噪声像素的亮度明显不同于其领域。冲击噪声突然，常由外部因素引起；其噪声范围可能相当大，不能通过提高信噪比来避免，是传输中的主要错误。
(4)量化噪声（Quatization Noise）：是指不同量化级别的噪声。例如，当图像的亮度水平降低一半时，就会出现伪轮廓。
4直接逆滤波
恢复退化图像最简单的方法是构成以下形式的估计：
然后使用傅里叶逆变器来获得相应的图像估计。这种方法变成了逆滤波器。我们可以在前面讨论的模型中表示估计
这个简单的表达告诉我们，如果我们及时准确地知道它，我们就无法恢复它，因为噪声重量是一个随机函数，他的傅里叶变换是未知的。此外，在实践中，有许多零的情况也是一个问题。即使项目可以忽略不计，它也将控制零值的恢复估计。
使用逆滤波器时，典型的方法是形成比例，然后将频率范围限制在接近原点的频率上，以获得逆滤波器。概念是，零不太可能出现在接近原点的地方，因为变换值通常是该区域的最高值。
5维纳滤波
维纳滤波（N.Wiener1942年提出的第一种方法)是最早、最著名的线性图像恢复方法。维纳滤波器寻找一个统计误差函数
我们感兴趣的两个量平均噪声功率和平均图像功率，分别定义为
在这种情况下，即使真实的比率未知，交互式地变化常量并观察复原的结果的实验就变成了一件简单的事。当然，假设函数为常量是一种粗糙的的近似。在前述滤波器方程中，用一个常量数组来代替就产生所谓的参数维纳滤波器。
在IPT中，维纳滤波是使用函数deconvwnr来实现的，函数deconvwnr有三种可能的语法形式。在所有的这些形式中，g代表退化图像，fr是复原图像。第一种语法形式
fr=deconvwnr（g，PSF）
假设信噪比为零。从而，维纳滤波器的这种形式就是前文中的逆滤波器。语法
fr=deconvwnr（g，PSF，NSPR）
假设噪信功率比已知，或是个常量或是个数组；函数接受其中的任何一个。这是用于实现参数维纳滤波器的语法，在这种情况下，NSPR可以是一个交互的标量输入。最后，语法
fr=deconvwnr（g，PSF，NACORR,FAVORR）
假设噪声和未退化图像的自相关函数NACORR和FAVORR是已知的。注意deconvwnr的这种形式使用和的自相关来代替这些函数的功率谱。从相关理论我们可知
其中，“”表示相关操作，表示傅立叶变换。这个表达式说明了对于deconvwnr的使用，通过计算功率谱的傅里叶变换，可以得到自相关函数。噪声的自相关有类似的注释。
若复原图像呈现出由算法中使用的离散傅立叶变换所引入的振铃，则它往往会提醒我们在调用
J=edgetaper(I,PSF)
该函数利用点扩散函数PSF模糊了输入图像I的边缘。这样，输出图像J就是图像I和I的模糊版本的加权和。这个由PSF的自相关函数所决定的加权数组在它的中心区域取J等于I，而在接近边缘的地带等于I的模糊版本。
使用deconvwnr函数复原模糊噪声图像：
f=imread('cameraman.jpg');
LEN=30;
THETA=45;
PSF=fspecial('motion',LEN,THETA);
MF=imfilter(f,PSF,'circular','conv');
wnr=deconvwnr(MF,PSF);
subplot(2,2,1);imshow(f);title('原始图像');
subplot(2,2,2);imshow(MF);title('模糊后的图像');
subplot(2,2,3);imshow(wnr);title('恢复后的图像');
经过仿真，如图所示，我们得到的结果虽然仍有一些噪声存在，但是已经和原图很接近了。因为原图像和噪声函数都是已知的，所以可以正确地估算参量。在实践中，当这些量之一（或更多）未知时，挑战便是在试验中智能地选择所用的函数，知道获得可接受的结果为止。
6 约束的最小二乘方滤波
约束最小二乘滤波恢复函数deconvreg的调用格式：deconvreg(I，PSF，NP，LRANGE，REGOP)，其中，I表示输入像，PSF表示点扩散函数，NP、LRANGE（ 输入） 和 REGO是可选参数，分别表示图像的噪声强度、拉氏算子的搜索范围和约束算子，同时，该函数也可以在指定的范围内搜索最优的拉氏算子。下面通过程序事例说明约束最小二乘滤波恢复图像：
I=imread ('cameraman.jpg');
PSF=fspecial('gaussian',10,4);
Blurred=imfilter(I,PSF,'conv');
V=.03;
BN=imnoise(Blurred,'gaussian',0,V);
NP=V*prod(size(I));
[reg LAGRA]=deconvreg(BN,PSF,NP);
Edged=edgetaper(BN,PSF);
reg2=deconvreg(Edged,PSF,NP/1.2);
reg3=deconvreg(Edged,PSF,[],LAGRA);
figure
subplot(2,3,1);imshow (I);title('原始图像');
subplot(2,3,2);imshow (BN);title('加入高斯噪声的图像');
subplot(2,3,3);imshow (reg);title('恢复后的图像');
subplot(2,3,4);imshow(reg2);title('振铃抑制图像');
subplot(2,3,5);imshow(reg3);title('拉格朗日算子恢复图像');
图6约束最小二乘恢复对比图
利用振铃抑制恢复图像是几种中恢复效果最好的，其他几种方法也可以恢复但是比较模糊，效果不是很明显。
7 使用Lucy-Richardson算法的迭代非线性复原
Lucy-Richardson(LR)算法假设图像服从Poission 分布，采用最大似然法进行估计，是一种基于贝叶斯分析的迭代算法。其最优估计以最大似然准则作为标准，即要使概率密度函数
大噪声的缺陷。因此，处理噪声项是 LR 算法应用于低信噪比图像复原的关键。
在IPT中，L-R算法是由名为deconvlucy的函数完成的，该函数的语法为
fr=deconvlucy(g,PSF,NUMIT,DAMPAR,WEIGHT)
其中，fr代表复原的图像，g代表退化的图像，PSF是点扩散函数，NUMIT为迭代次数（默认为10次），DAMPAR和WEIGHT定义如下。
DAMPAR是一个标量，它指定了结果图像与原图像g之间的偏离阈值。当像素偏离原值的范围在DAMPAR之内时，就不用再迭代。这既抑制了这些像素上的噪声，又保存了必要的图像细节。默认值为0（无衰减）。
WEIGHT是一个与g同样大小的数组，它为每一个像素分配一个权重来反映其重量。当用一个指定的PSF来模拟模糊时，WEIGHT
可以从计算像素中剔除那些来自图像边界的像素点，因此，PSF造成的模糊是不同的。若PSF的大小为，则在WEIGHT中用到的零边界的宽度是ceil（n/2）。默认值是同输入图像g同等大小的一个单位数组。若复原图像呈现出由算法中所用的离散傅里叶变换所引入的振铃，则在调用函数deconvlucy之前，要利用函数edgetaper。
下面通过程序事例说明 Lucy-Richardson算法恢复图像：
I=imread('lajiao.jpg');
PSF=fspecial('gaussian',5,5) ;
Blurred=imfilter(I,PSF,'symmetric','conv');
V=.003;
BN=imnoise(Blurred,'gaussian',0,V);
luc=deconvlucy(BN,PSF,5);
figure
subplot(2,2,1);imshow(I);title('原始图像');
subplot(2,2,2);imshow (Blurred);title('模糊后的图像');
subplot(2,2,3);imshow (BN);title('加噪后的图像');
subplot(2,2,4);imshow (luc);title('恢复后的图像');
图7.1 Lucy-Richardson恢复对比图
从图中我们可以看到，经多次迭代，尤其是在低信噪比情况下，重建图像可能会出现一些斑点，这些斑点并不代表图像的真实结构，是输出图像过于逼近噪声所产生的结果。
8 盲去卷积
通常图像恢复方法均在成像系统的点扩展函数PSF已知下进行, 实际上它通常是未知的. 在 PSF未知的情况下, 盲去卷积是实现图像恢复的有效方法。因此，把那些不以PSF知识为基础的图像复原方法统称为盲去卷积算法。
在过去的20年里，一种盲去卷积的方法已经受到了人们的极大重视，它是以最大似然估计(MLE)为基础，即一种用被随机噪声所干扰的量进行估计的最优化策略。简要的说，关于MLE方法的一种解释就是将图像数据看成随机量，它们与另外一族可能的随机量之间有着某种似然性。似然函数用、和来加以表达，然后，问题就变成了寻求最大似然函数。在盲去卷积中，最优化问题规定的约束条件并假定收敛时通过迭代来求解，得到的最大和就是还原的图像和PSF。
工具箱通过函数deconvblind来执行盲去卷积，它有如下语法：
[f,PSFe]=deconvblind(g,INITPSF)
其中，g代表退化函数，INITPSF是点扩散函数的出事估计。PSFe是这个函数最终计算到的估计值，fr是利用估计的PSF复原的图像。 用来去的复原图像的算法是L-R迭代复原算法。PSF估计受其初始推测尺寸的巨大影响，而很少受其值的影响。
若复原图像呈现出由算法中使用的离散傅里叶变换所引入的振铃，则我们在调用函数deconvblind值钱，通常要使用函数edgetaper。下面我们使用函数deconvblind估计PSF：
I=imread('lajiao'.jpg');
PSF=fspecial('motion',10,30);
Blurred=imfilter(I,PSF,'circ','conv') ;
INITPSF=ones(size(PSF));
[J P]=deconvblind (Blurred,INITPSF,20);
figure
subplot(2,2,1);imshow (I);title('原始图像');
subplot(2,2,2);imshow (Blurred);title('模糊后的图像')
subplot(2,2,3);imshow (J);title('初步恢复后的图像');
该算法优点是,同时恢复了图像和点扩张函数,在对失真情况毫无先验知识的情况下, 仍能实现对模糊图像的恢复操作。利用 MATLAB实现的图像恢复, 并对恢复图像的失真情况做了改善。在进行图像恢复时,重建 PSF,对图像进行重建, 得到恢复的图像。
总 结
本文介绍了图像退化的原因并且简要介绍了当前主流的图像复原方法，并通过对各种复原方法的仿真，了解了各种方法的优劣性，为我们在实际生活提供依据。但是无论是哪一种方法都有所局限性，我们应该努力致力于研究新型的优秀的图像复原方法，来获得更好的图像复原效果。同时，我们知道算法利用的信息越多信息的准确性越高，则复原图像的质量就越高。而且采用 MATLAB实现图像恢复，通过几条简单的MATLAB命令就可完成一大串高级计算机语言才能完成的任务,简捷明快。大多数图像处理模型是可以通过使用MATLAB的基本函数通过编程实现的。
参考文献
[1]张燕妮.离焦模糊复原方法研究[D]．大连：大连理工大学，2005
[2]徐飞.施晓红. MATLAB应用图像处理[M]．西安电子科技大学出版社，2002.
[3]路敬袜，王冬梅. 基于MATLAB的图像复原算法的研究[D]．硕士学位论文.
[4]阮秋琦等译. 数字图像处理（MATLAB版）[M]．北京：电子工业出版社，2005.
[5]张德丰,张葡青.维纳滤波图像恢复的理论分析与实现[J].中山大学学 报,2006,45(6)44-47.
[6]14曹茂永，孙农亮，郁道银．离焦模糊图像清晰度评价函数的研究[J]．仪器仪表学报，2007，
22(3)：259-268．
[7]章毓晋．图像处理和分析[M]．北京：清华大学出版社，1999．
[8]苏开亮.运动模糊图像的恢复及恢复质量评价[D]．西安：西安电子科技大学出版社，2010，21-25．
[9]龚声蓉，刘纯平，王强. 数字图像处理与分析[M]，清华大学出版社.2006.
[10]Edmund Y L. Digital restoration of defocus image in the wavelet domain[J].Applied Optics.2002,141(126):4806-4811.
[11]Tnoa Zznia.Blur identification analysis in blind image deconvolution using Markov random fields[J]
.Pattern Recognition and Image Analysis,2001,11(4):699-700.
[12]Williamk,Faramarzd.Fast computational techniques for Pseudoinverse and Wiener image restoration [J].IEEE Trans on computer,1997,26(6)571-580.