123005_基于MATLAB的图像恢复算法研究.doc

中北大学
课 程 设 计 说 明 书

学生姓名： 学 号：
学生姓名： 学 号：
学生姓名： 学 号：
学 院： 信息商务学院
专 业： 电子信息工程
题 目： 信息处理综合实践:
基于MATLAB的图像恢复算法研究

指导教师： 职称:

年 月 日

中北大学
课程设计任务书

13/14 学年第 一 学期

学 院： 信息商务学院
专 业： 电子信息工程
学 生 姓 名： 学 号：
学 生 姓 名： 学 号：
学 生 姓 名： 学 号:
课程设计题目： 信息处理综合实践:
于MATLAB的图像恢复算法研究
起 迄 日 期：
课程设计地点： 电子信息科学与技术专业实验室

指 导 教 师：
系 主 任：

下达任务书日期: 年月日
课 程 设 计 任 务 书
1．设计目的：
1、通过本课程设计的学习，学生将复习所学的专业知识，使课堂学习的理论知识应用于实践，通过本课程设计的实践使学生具有一定的实践操作能力；
2、掌握Matlab使用方法，能熟练运用该软件设计并完成相应的信息处理；
3、通过图像处理实践的课程设计，掌握设计图像处理软件系统的思维方法和基本开发过程。
2．设计内容和要求（包括原始数据、技术参数、条件、设计要求等）：
(1)创建一个仿真运动/均值模糊PSF来模糊一幅图像（图像自选）；
(2)针对退化设计出复原滤波器，对退化图像进行复原(复原的方法自定)；
(3)对退化图像进行复原，显示复原前后图像，对复原结果进行分析，并评价复原算法；
(4)要求每位学生进行查阅相关资料，并写出自己的报告。注意每个学生的报告要有所侧重，写出自己所做的内容。

3．设计工作任务及工作量的要求〔包括课程设计计算说明书(论文)、图纸、实物样品等〕：

每个同学独立完成自己的任务，每人写一份设计报告，在课程设计论文中写明自己设计的部分，给出设计结果。

课 程 设 计 任 务 书
4．主要参考文献：
1. 阮秋琦等.数字图像处理（第三版）.北京：电子工业出版社.2011
2. 冈萨雷斯等.数字图像处理（MATLAB版）.北京：电子工业出版社.2001
3. 杨杰等.数字图像处理及MATLAB实现：学习与实验指导.北京：电子工业出版社.2010
4. 刘卫国等.MATLAB程序设计与应用.北京：高等教育出版社.2006
5. 许国根等.模式识别与智能计算的MATLAB实现.北京：北京航空航天大学出版社.2012

5．设计成果形式及要求：
毕业设计说明书
仿真结果
6．工作计划及进度：
2013年
12月16日 ~ 12月19日：查资料；
12月19日 ~ 12月24日：在指导教师指导下设计方案；
12月25日 ~ 12月27日：撰写课程设计说明书；
12月27日：答辩
系主任审查意见：

签字：
年 月 日

目 录
摘　要： 5
1．图像复原的概念 6
1.1图像复原的定义 6
1.2 图象恢复与图象增强的异同 6
1.3 图象退化的原因 6
1.4 维纳滤波的研究历史 6
1.5图象退化举例 7
2．退化模型 7
2.1图象退化模型概述 7
2.2连续函数退化模型 8
2.3离散函数退化模型 8
3.图象复原技术 9
3.1无约束恢复 9
3.2 逆滤波 9
3.3 维纳（Wiener）滤波器基本原理 9
3.4维纳滤波复原法 11
3.5图像复原例图 12
4．图像复原的MATLAB实现实例 13
5.结束语 13
参考文献： 14
附录： 14
（1）.维纳滤波复原源代码： 14
（2）.规则化滤波复原程序源代码： 14
（3）.Lucy-Richardson复原滤波源代码： 15
（4）.盲目去卷积复原源代码： 15

摘　要：
图像复原是图象处理的一个重要课题。图像复原也称图象恢复，是图象处理中的一大类技术。它的主要目的是改善给定的图像质量。当给定了一幅退化了的或者受到噪声污染了的图像后，利用退化现象的某种先验知识来重建或恢复原有图像是复原处理的基本过程。可能的退化有光学系统中的衍射，传感器非线性畸变，光学系统的像差，摄影胶片的非线性，大气湍流的扰动效应，图像运动造成的模糊及几何畸变等等。噪声干扰可以由电子成像系统传感器、信号传输 过程或者胶片颗粒性造成。各种退化图像的复原都可归结为一种过程，具体地说就是把退化模型化，并且采用相反的过程进行处理，以便恢复出原图像。文章介绍了图象退化的原因，几种常用的图像滤波复原技术，以及用MATLAB实现图像复原的方法。

1．图像复原的概念
1.1图像复原的定义
图像复原也称图象恢复，是图象处理中的一大类技术。所谓图像复原，是指去除或减轻在获取数字图像过程中发生的图像质量下降（退化）这些退化包括由光学系统、运动等等造成图像的模糊，以及源自电路和光度学因素的噪声。图像复原的目标是对退化的图像进行处理，使它趋向于复原成没有退化的理想图像。成像过程的每一个环节（透镜，感光片，数字化等等）都会引起退化。在进行图像复原时，既可以用连续数学，也可以用离散数学进行处理。其次，处理既可在空间域，也可在频域进行。
1.2 图象恢复与图象增强的异同
相同点：改进输入图像的视觉质量 。
不同点：图象增强目的是取得较好的视觉结果(不考虑退化原因)； 图象恢复根据相应的退化模型和知识重建或恢复原始的图像(考虑退化原因)。
1.3 图象退化的原因
图象退化指由场景得到的图像没能完全地反映场景的真实内容，产生了
失真等问题。其原因是多方面的。如：
透镜象差/色差
聚焦不准（失焦，限制了图像锐度）
模糊（限制频谱宽度）
噪声（是一个统计过程）
抖动（机械、电子）
1.4 维纳滤波的研究历史
维纳是著名的数学家，后来被誉为信息理论家。维纳的著作不仅是一个很好的创见，而且具有结合工程的实际意义，是线性滤波理论研究的一个重要的开端.
在第二次世界大战中，由于雷达的发明以及防空炮火控制的任务，把大量有修养的数学家和物理学家都动员到信息科学这个研究领域中来了，这个时候人们活跃于这个领域，并有许多重大的科学创造。数学家维纳对于滤波理论的研究成果，就是这时候重大的科学创见之一。
通讯与控制中的滤波问题，指的是从获得的信号与干扰中尽可能地滤除干扰，分离出所期望的信号，或者说，是通过对一系列带有误差的实际测量数据的处理，得出所期望数据的估计值川。维纳的工作是从研究处在统计平衡的时间序列开始的，维纳证明:在一定条件下，处在统计平衡的时间序列的时间平均等于相平均。维纳正是基于这点提出了他著名的滤波和预测理论。滤波问题就是尽可能地恢复一个被噪声干扰了的信号的问题。实质上，就是预测一个被噪声干扰了的时间序列的问题，因此，滤波问题也可以视为一个预测问题。数学上讲，预测就是从一个时间序列的过去的数据估算整个序列的统计参数。
工程上的滤波问题也是理论上的一类统计估计问题，最佳线形滤波是最佳线性估计的方法之一，在最佳估计中最小均方误差估计是最有现实意义的。估计理论的课题是众多的，最小均方误差估计只是估计理论的一个小的分支。然而，它却是最重要又最富有实际意义的一个分支，对系统所加的线性条件起初是为了简化理论分析，非线性滤波问题是在理论处理上比线性滤波问题要困难和复杂的多，但是后来证明:在一定条件下，在最小均方误差准则下得到的最佳线性系统是所有系统中的最佳者。
近代滤波理论的发展对于信息科学的发展是有重大贡献的，它概括了通讯与控制中信息过滤的统计本质。这是由于滤波理论与通讯和控制中的许多课题有密切的联系，从而赋予了滤波理论以极大的生命力，滤波理论本来是一个小的研究领域，但是它联系着许多大的广泛的研究领域，因此它的价值己经超出了它起源时自身的价值，也就是它能够继续活跃地向前发展的保证。
几十年来滤波理论已经发展成了一个广阔的研究领域，可以有许多不同的方法来介绍它的内容，有的可以选择不同的重点。本文主要是关于维纳滤波的，介绍维纳滤波的基本概念以及讲其维纳滤波的应用。
从数学的观点来说滤波理论是统计学中的估计理论的一个重要分支，从工程的观点来看它又是系统工程研究的一个重要组成部分。
1.5图象退化举例
如图1所示是两个图象退化的例子。

图1 退化图像与原始图像

2．退化模型
2.1图象退化模型概述
图像复原处理的关键问题在于建立退化模型。在用数学方法描述图像时，它的最普遍的数学表达式为

这样一个表达式可以代表一幅活动的、彩色的立体图像。当研究的是静止的、单色的、平面的图像时，则其数学表达式就简化为

基于这样的数学表达式，可建立如图2所示的退化模型。由图2的模型可见，一幅纯净的图像是由于通过了一个系统H及加性噪声而使其退化为一幅图像的。

图2 图像退化模型

图像复原可以看成是一个估计过程。如果已经给出了退化图像并估计出系统参数H，从而可近似地恢复。这里，是一种统计性质的噪声信息。当然，为了对处理结果做出某种最佳的估计，一般应首先明确一个质量标准。根据图像的退化模型及复原的基本过程可见，复原处理的关键在于对系统H的基本了解。就一般而言，系统是某些元件或部件以某种方式构造而成的整体。退化模型可分为连续函数退化模型和离散函数退化模型。
2.2连续函数退化模型

假定系统H对坐标为(a,b)处的冲激函数d(x-a,y-b)的冲激响应为h(x,a,y,b)，则

此式说明，如果系统H对冲激函数的响应为已知，则对任意输入的响应可用上式求得，即，线性系统H完全可以由冲激响应来表征。图像中冲激响应也称为点扩散函数。
在有噪音的情况下：

2.3离散函数退化模型
对和进行均匀取样后，就可引伸出离散函数的退化模型。用一维的来说明。如果f (x)和h(x)周期分别A和B的序列，为避免卷积周期重叠需要对它们进行周期扩展为周期为M ≥ A + B – 1。

f(x) 0 ≤ x ≤ A-1 h(x) 0 ≤ x≤ B-1
fe(x)= he(x)=
0 A-1≤ x ≤ M-1 0 B-1< x ≤ M-1

那么它们的时域离散卷积可定义为下式：

显然，上式也是具有周期M的序列。

如果用矩阵来表示上述离散退化模型，可写成下式之形式:

退化过程为：

图像f(x,y)被线性操作h(x,y)所模糊，并叠加上噪声n(x,y)，构成了退化后的图像g(x,y)。退化后的图像与复原滤波器卷积得到复原的f(x,y)图像。

图3 图像的退化/复原过程模型

3.图象复原技术
3.1无约束恢复
由退化模型得：
最小均方误差准：
在最小二乘方意义上说，希望找到一个使下式的值最小:

3.2 逆滤波
设M = N ，则：

退化函数H (u, v)与F (u, v)相乘为退化过程，用H (u, v)去除G (u, v) 是复原过程，称其为逆滤波。可描述为：

记M (u, v)为复原转移函数，则其等于1 / H (u, v).

3.3 维纳（Wiener）滤波器基本原理
维纳(Wiener)滤波是用来解决从噪声中提取信号问题的一种过滤(或滤波)的方法。实际上这种线性滤波问题，可以看成是一种估计问题或一种线性估计问题。
一个线性系统，如果它的单位样本响应为，当输入一个随机信号，且

其中表示信号，表示噪声，则输出为

我们希望通过线性系统后得到的尽量接近于,因此称为的估计值，用表示，即

=+

图4.1 维纳滤波器的输入一输出关系

如图4-1所示。这个线性系统称为对于的一种估计器。
实际上，式(2.2)的卷积形式可以理解为从当前和过去的观察值, ,…,…来估计信号的当前值。因此，用进行过滤的问题可以看成是一个估计问题。由于我们现在涉及的信号是随机信号，所以这样一种过滤问题实际上是一种统计估计问题。
一般，从当前的和过去的观察值,,…估计当前的信号值称为过滤或滤波;从过去的观察值，估计当前的或将来的信号值 称为预测或外推;从过去的观察值，估计过去的信号值称为平滑或内插。因此维纳过滤又常常被称为最佳线性过滤与
预测或线性最优估计。这里所谓最佳与最优是以最小均方误差为准则的。这里只讨论过滤与预测问题。
如果我们以:与分别表示信号的真值与估计值，而用表示它们之间的误差，即

显然，可能是正的，也可能是负的，并且它是一个随机变量。因此，用它的均方值来表达误差是合理的，所谓均方误差最小即它的平方的统计平均值最小:

采用最小均方误差准则作为最佳过滤准则的原因还在于它的理论分析比较简单，不要求对概率的描述。并且在这种准则下导出的最佳线性系统对其它很广泛一类准则而言也是最佳的。
3.4维纳滤波复原法
维纳滤波法是由Wiener首先提出的，应用于一维信号处理，取得了很好的效果。之后，维纳滤波法被用于二维信号处理，也取得了不错的效果，尤其在图像复原领域，由于维纳滤波计算量小，复原效果好，从而得到了广泛的应用和发展。
维纳滤波器寻找一个使统计误差函数

最小的估计。是期望值操作符，是未退化的图像。该表达式在频域可表示为

其中，
表示退化函数

表示的复共轭
表示噪声的功率谱
表示未退化图像的功率谱
比率称为信噪功率比。在中维纳滤波使用函数deconvwnr来实现的。

它一种最小均方误差滤波器。

设 Rf 是 f 的相关矩阵：

Rf 的第 ij 元素是E{fi fj}，代表 f 的第 i 和第 j 元素的相关。

设 Rn是n 的相关矩阵：

根据两个象素间的相关只是它们相互距离而不是位置的函数的假设，可将Rf和Rn都用块循环矩阵表达，并借助矩阵W来对角化：

fe(x, y)的功率谱，记为Sf (u, v) ；ne(x, y)的功率谱，记为Sn(u, v)。D是1个对角矩阵，D(k, k) = l(k)，则有：

定义：

代入：

两边同乘以W –1，有：

最后整理得：

3.5图像复原例图

以下的几幅图是用MATLAB软件根据不同的复原方法进行的图像复原。根据图4例图可看出不同复原方法的区别。

原图 退化图像

全逆滤波 半径受限逆滤波 维纳滤波结果

图4复原例图

4．图像复原的MATLAB实现实例

维纳滤波复原 规则化滤波复原

Lucy-Richardson复原 盲目去卷积复原

图5 图像复原实例
5.结束语
本文简要介绍了图像退化的原因，图像退化的模型，图像复原的概念，几种常用的图像复原的方法，以及利用MATLAB实现图像复原的几个例子。简单的讲述了MATLAB在图像复原中的应用。
参考文献：
[1] 阮秋琦编著.—2版。北京：电子工业出版社，2007.2.
[2] (美)卡斯尔曼(castleman，k．R)著；朱志刚等译.数字图像处理.北京：电子工业出版社，2002.2.
[3] 孙家广等主编.计算机图形学.第3版.北京：清华大学出版社，1998.9.
[4] 罗军辉等主编.MATLAB7.0在图像处理中的应用.第1版.北京：机械工业出版社，2007.7.
附录：
（1）.维纳滤波复原源代码：
I=checkerboard(8); noise=0.1*randn(size(I));
PSF=fspecial('motion',21,11);
Blurred=imfilter(I,PSF,'circular');
BlurredNoisy=im2uint8(Blurred+noise);
NP=abs(fftn(noise)).^2;
NPOW=sum(NP(:)/numel(noise));
NCORR=fftshift(real(ifftn(NP)));
IP=abs(fftn(I)).^2;
IPOW=sum(IP(:)/numel(noise));
ICORR=fftshift(real(ifftn(IP)));
ICORR1=ICORR(:,ceil(size(I,1)/2));
NSR=NPOW/IPOW;
subplot(221);imshow(BlurredNoisy,[]);
title('模糊和噪声图像');
subplot(222);imshow(deconvwnr(BlurredNoisy,PSF,NSR),[]);
title('deconbwnr(A,PSF,NSR)');
subplot(223);imshow(deconvwnr(BlurredNoisy,PSF,NCORR,ICORR),[]);
title('deconbwnr(A,PSF,NCORR,ICORR)');
subplot(224);imshow(deconvwnr(BlurredNoisy,PSF,NPOW,ICORR1),[]);
title('deconbwnr(A,PSF,NPOW,ICORR_1_D)');
（2）.规则化滤波复原程序源代码：
I=checkerboard(8);
PSF=fspecial('gaussian',7,10);
V=.01;
BlurredNoisy=imnoise(imfilter(I,PSF),'gaussian',0,V);
NOISEPOWER=V*numel(I);
[J LAGRA]=deconvreg(BlurredNoisy,PSF,NOISEPOWER);
subplot(221);imshow(BlurredNoisy);
title('A=Blurred and Noisy');
subplot(222);imshow(J);
title('[J LAGRA]=deconvreg(A,PSF,NP)');
subplot(223);imshow(deconvreg(BlurredNoisy,PSF,[],LAGRA/10));
title('deconvreg(A,PSF,[],0.1*LAGRA)');
subplot(224);imshow(deconvreg(BlurredNoisy,PSF,[],LAGRA/10))
title('deconvreg(A,PSF,[],10*LAGRA');
（3）.Lucy-Richardson复原滤波源代码：
I=checkerboard(8);
PSF=fspecial('gaussian',7,10);
V=.0001;
BlurredNoisy=imnoise(imfilter(I,PSF),'gaussian',0,V);
WT=zeros(size(I));
WT(5:end-4,5:end-4)=1;
J1=deconvlucy(BlurredNoisy,PSF);
J2=deconvlucy(BlurredNoisy,PSF,20,sqrt(V));
J3=deconvlucy(BlurredNoisy,PSF,20,sqrt(V),[],WT);
subplot(221);imshow(BlurredNoisy);
title('A=Blurred and Noisy');
subplot(222);imshow(J1);
title('deconvlucy(A,PSF)');
subplot(223);imshow(J2);
title('deconvlucy(A,PSF,NI,DP)');
subplot(224);imshow(J3);
title('deconvlucy(A,PSF,NI,DP,[],WT)');
（4）.盲目去卷积复原源代码：
I=checkerboard(8);
PSF=fspecial('gaussian',7,10);
V=.0001;
BlurredNoisy=imnoise(imfilter(I,PSF),'gaussian',0,V);
WT=zeros(size(I));WT(5:end-4,5:end-4)=1;
INITPSF=ones(size(PSF));FUN=inline('PSF+P1','PSF','P1');
[J P]=deconvblind(BlurredNoisy,INITPSF,20,10*sqrt(V),WT,FUN,0);
subplot(221);imshow(BlurredNoisy);title('A=Blurred and Noisy');
subplot(222);imshow(PSF,[]);title('True PSF');
subplot(223);imshow(J);title('Deblured Image');
subplot(224);imshow(P,[]);
16