基于MATLAB的图像恢复算法研究

1、PAGE PAGE 28中北大学课 程 设 计 说 明 书学生姓名： 学 号： 学生姓名： 学 号： 学生姓名： 学 号： 学 院： 信息商务学院 专 业： 电子信息工程 题 目： 信息处理综合实践: 基于MATLAB的图像恢复算法研究 指导教师： 职称: 年 月 日中北大学课程设计任务书 13/14 学年第 一 学期学 院： 信息商务学院 专 业： 电子信息工程 学 生 姓 名： 学 号： 学 生 姓 名： 学 号： 学 生 姓 名： 学 号: 课程设计题目： 信息处理综合实践: 于MATLAB的图像恢复算法研究 起 迄 日 期： 课程设计地点： 电子信息科学与技术专业实验室 指 导 教 师

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

3、前后图像，对复原结果进行分析，并评价复原算法；(4)要求每位学生进行查阅相关资料，并写出自己的报告。注意每个学生的报告要有所侧重，写出自己所做的内容。3设计工作任务及工作量的要求包括课程设计计算说明书(论文)、图纸、实物样品等：每个同学独立完成自己的任务，每人写一份设计报告，在课程设计论文中写明自己设计的部分，给出设计结果。课 程 设 计 任 务 书4主要参考文献：阮秋琦等.数字图像处理（第三版）.北京：电子工业出版社.20112. 冈萨雷斯等.数字图像处理（MATLAB版）.北京：电子工业出版社.20013. 杨杰等.数字图像处理及MATLAB实现：学习与实验指导.北京：电子工业出版社.20

4、104. 刘卫国等.MATLAB程序设计与应用.北京：高等教育出版社.20065. 许国根等.模式识别与智能计算的MATLAB实现.北京：北京航空航天大学出版社.20125设计成果形式及要求：毕业设计说明书仿真结果6工作计划及进度：2013年12月16日 12月19日：查资料；12月19日 12月24日：在指导教师指导下设计方案；12月25日 12月27日：撰写课程设计说明书； 12月27日：答辩系主任审查意见： 签字： 年 月 日目 录 TOC o 1-3 h z u HYPERLINK l _Toc375851342 摘要： PAGEREF _Toc375851342 h 5 HYPERL

5、INK l _Toc375851343 1图像复原的概念 PAGEREF _Toc375851343 h 6 HYPERLINK l _Toc375851344 1.1图像复原的定义 PAGEREF _Toc375851344 h 6 HYPERLINK l _Toc375851345 1.2 图象恢复与图象增强的异同 PAGEREF _Toc375851345 h 6 HYPERLINK l _Toc375851346 1.3 图象退化的原因 PAGEREF _Toc375851346 h 6 HYPERLINK l _Toc375851347 1.4 维纳滤波的研究历史 PAGEREF _

6、Toc375851347 h 6 HYPERLINK l _Toc375851348 1.5图象退化举例 PAGEREF _Toc375851348 h 7 HYPERLINK l _Toc375851349 2退化模型 PAGEREF _Toc375851349 h 7 HYPERLINK l _Toc375851350 2.1图象退化模型概述 PAGEREF _Toc375851350 h 7 HYPERLINK l _Toc375851351 2.2连续函数退化模型 PAGEREF _Toc375851351 h 8 HYPERLINK l _Toc375851352 2.3离散函数退化

7、模型 PAGEREF _Toc375851352 h 8 HYPERLINK l _Toc375851353 3.图象复原技术 PAGEREF _Toc375851353 h 9 HYPERLINK l _Toc375851354 3.1无约束恢复 PAGEREF _Toc375851354 h 9 HYPERLINK l _Toc375851355 3.2 逆滤波 PAGEREF _Toc375851355 h 9 HYPERLINK l _Toc375851356 3.3 维纳（Wiener）滤波器基本原理 PAGEREF _Toc375851356 h 9 HYPERLINK l _To

8、c375851357 3.4维纳滤波复原法 PAGEREF _Toc375851357 h 11 HYPERLINK l _Toc375851358 3.5图像复原例图 PAGEREF _Toc375851358 h 12 HYPERLINK l _Toc375851359 4图像复原的MATLAB实现实例 PAGEREF _Toc375851359 h 13 HYPERLINK l _Toc375851360 5.结束语 PAGEREF _Toc375851360 h 13 HYPERLINK l _Toc375851361 参考文献： PAGEREF _Toc375851361 h 14

9、HYPERLINK l _Toc375851362 附录： PAGEREF _Toc375851362 h 14 HYPERLINK l _Toc375851363 （1）.维纳滤波复原源代码： PAGEREF _Toc375851363 h 14 HYPERLINK l _Toc375851364 （2）.规则化滤波复原程序源代码： PAGEREF _Toc375851364 h 14 HYPERLINK l _Toc375851365 （3）.Lucy-Richardson复原滤波源代码： PAGEREF _Toc375851365 h 15 HYPERLINK l _Toc3758513

10、66 （4）.盲目去卷积复原源代码： PAGEREF _Toc375851366 h 15摘要：图像复原是图象处理的一个重要课题。图像复原也称图象恢复，是图象处理中的一大类技术。它的主要目的是改善给定的图像质量。当给定了一幅退化了的或者受到噪声污染了的图像后，利用退化现象的某种先验知识来重建或恢复原有图像是复原处理的基本过程。可能的退化有光学系统中的衍射，传感器非线性畸变，光学系统的像差，摄影胶片的非线性，大气湍流的扰动效应，图像运动造成的模糊及几何畸变等等。噪声干扰可以由电子成像系统传感器、信号传输 过程或者胶片颗粒性造成。各种退化图像的复原都可归结为一种过程，具体地说就是把退化模型化，并且

11、采用相反的过程进行处理，以便恢复出原图像。文章介绍了图象退化的原因，几种常用的图像滤波复原技术，以及用MATLAB实现图像复原的方法。1图像复原的概念1.1图像复原的定义图像复原也称图象恢复，是图象处理中的一大类技术。所谓图像复原，是指去除或减轻在获取数字图像过程中发生的图像质量下降（退化）这些退化包括由光学系统、运动等等造成图像的模糊，以及源自电路和光度学因素的噪声。图像复原的目标是对退化的图像进行处理，使它趋向于复原成没有退化的理想图像。成像过程的每一个环节（透镜，感光片，数字化等等）都会引起退化。在进行图像复原时，既可以用连续数学，也可以用离散数学进行处理。其次，处理既可在空间域，也可在

12、频域进行。1.2 图象恢复与图象增强的异同相同点：改进输入图像的视觉质量 。不同点：图象增强目的是取得较好的视觉结果(不考虑退化原因)； 图象恢复根据相应的退化模型和知识重建或恢复原始的图像(考虑退化原因)。1.3 图象退化的原因 图象退化指由场景得到的图像没能完全地反映场景的真实内容，产生了失真等问题。其原因是多方面的。如： 透镜象差/色差 聚焦不准（失焦，限制了图像锐度） 模糊（限制频谱宽度） 噪声（是一个统计过程）抖动（机械、电子）1.4 维纳滤波的研究历史维纳是著名的数学家，后来被誉为信息理论家。维纳的著作不仅是一个很好的创见，而且具有结合工程的实际意义，是线性滤波理论研究的一个重要的

13、开端. 在第二次世界大战中，由于雷达的发明以及防空炮火控制的任务，把大量有修养的数学家和物理学家都动员到信息科学这个研究领域中来了，这个时候人们活跃于这个领域，并有许多重大的科学创造。数学家维纳对于滤波理论的研究成果，就是这时候重大的科学创见之一。 通讯与控制中的滤波问题，指的是从获得的信号与干扰中尽可能地滤除干扰，分离出所期望的信号，或者说，是通过对一系列带有误差的实际测量数据的处理，得出所期望数据的估计值川。维纳的工作是从研究处在统计平衡的时间序列开始的，维纳证明:在一定条件下，处在统计平衡的时间序列的时间平均等于相平均。维纳正是基于这点提出了他著名的滤波和预测理论。滤波问题就是尽可能地恢

14、复一个被噪声干扰了的信号的问题。实质上，就是预测一个被噪声干扰了的时间序列的问题，因此，滤波问题也可以视为一个预测问题。数学上讲，预测就是从一个时间序列的过去的数据估算整个序列的统计参数。 工程上的滤波问题也是理论上的一类统计估计问题，最佳线形滤波是最佳线性估计的方法之一，在最佳估计中最小均方误差估计是最有现实意义的。估计理论的课题是众多的，最小均方误差估计只是估计理论的一个小的分支。然而，它却是最重要又最富有实际意义的一个分支，对系统所加的线性条件起初是为了简化理论分析，非线性滤波问题是在理论处理上比线性滤波问题要困难和复杂的多，但是后来证明:在一定条件下，在最小均方误差准则下得到的最佳线性

15、系统是所有系统中的最佳者。 近代滤波理论的发展对于信息科学的发展是有重大贡献的，它概括了通讯与控制中信息过滤的统计本质。这是由于滤波理论与通讯和控制中的许多课题有密切的联系，从而赋予了滤波理论以极大的生命力，滤波理论本来是一个小的研究领域，但是它联系着许多大的广泛的研究领域，因此它的价值己经超出了它起源时自身的价值，也就是它能够继续活跃地向前发展的保证。 几十年来滤波理论已经发展成了一个广阔的研究领域，可以有许多不同的方法来介绍它的内容，有的可以选择不同的重点。本文主要是关于维纳滤波的，介绍维纳滤波的基本概念以及讲其维纳滤波的应用。 从数学的观点来说滤波理论是统计学中的估计理论的一个重要分支，

16、从工程的观点来看它又是系统工程研究的一个重要组成部分。1.5图象退化举例如图1所示是两个图象退化的例子。图1 退化图像与原始图像2退化模型2.1图象退化模型概述图像复原处理的关键问题在于建立退化模型。在用数学方法描述图像时，它的最普遍的数学表达式为这样一个表达式可以代表一幅活动的、彩色的立体图像。当研究的是静止的、单色的、平面的图像时，则其数学表达式就简化为基于这样的数学表达式，可建立如图2所示的退化模型。由图2的模型可见，一幅纯净的图像是由于通过了一个系统H及加性噪声而使其退化为一幅图像的。图2 图像退化模型图像复原可以看成是一个估计过程。如果已经给出了退化图像并估计出系统参数H，从而可近似

17、地恢复。这里，是一种统计性质的噪声信息。当然，为了对处理结果做出某种最佳的估计，一般应首先明确一个质量标准。根据图像的退化模型及复原的基本过程可见，复原处理的关键在于对系统H的基本了解。就一般而言，系统是某些元件或部件以某种方式构造而成的整体。退化模型可分为连续函数退化模型和离散函数退化模型。2.2连续函数退化模型假定系统H对坐标为(,)处的冲激函数(x-,y-)的冲激响应为h(x,y,)，则此式说明，如果系统H对冲激函数的响应为已知，则对任意输入的响应可用上式求得，即，线性系统H完全可以由冲激响应来表征。图像中冲激响应也称为点扩散函数。在有噪音的情况下：2.3离散函数退化模型对和进行均匀取样

18、后，就可引伸出离散函数的退化模型。用一维的来说明。如果f (x)和h(x)周期分别A和B的序列，为避免卷积周期重叠需要对它们进行周期扩展为周期为M A + B 1。 f(x) 0 x A-1 h(x) 0 x B-1fe(x)= he(x)=0 A-1 x M-1 0 B-10.2)=0;Y=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,

19、Facecolor,interp,Edgecolor,none,Facelighting,phong); 二、理想高通滤波器IA=imread(lena.bmp);f1,f2=freqspace(size(IA),meshgrid);Hd=ones(size(IA);r=sqrt(f1.2+f2.2);Hd(r0.2)=0;Y=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(u

20、int8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); Butterworth低通滤波器IA=imread(lena.bmp);f1,f2=freqspace(size(IA),meshgrid);D=0.3;r=f1.2+f2.2;n=4;for i=1:size(IA,1) for j=1:size(IA,2) t=r(i,j)/(D*D); Hd(i,j)=1/(tn+1); endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift

21、(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); Butterworth高通滤波器IA=imread(lena.bmp);f1,f2=freqspace(size(IA),meshgrid);D=0.3;r=f1.2+f2.2;n=4;for i=1:size(IA,1) for j=1:size(IA,2) t=(D*D)/r

22、(i,j); Hd(i,j)=1/(tn+1); endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 高斯低通滤波器IA=imread(lena.bmp);IB=imread(babarra.bmp);f1,f2=

23、freqspace(size(IA),meshgrid);D=100/size(IA,1);r=f1.2+f2.2;Hd=ones(size(IA);for i=1:size(IA,1) for j=1:size(IA,2) t=r(i,j)/(D*D); Hd(i,j)=exp(-t); endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);fi

24、guresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 高斯高通滤波器IA=imread(lena.bmp);IB=imread(babarra.bmp);f1,f2=freqspace(size(IA),meshgrid);%D=100/size(IA,1);D=0.3;r=f1.2+f2.2;for i=1:size(IA,1) for j=1:size(IA,2) t=r(i,j)/(D*D); Hd(i,j)=1-exp(-t); endendY=fft2(double(IA);Y=fftshift(Y);Ya=

25、Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 梯形低通滤波器IA=imread(lena.bmp);IB=imread(babarra.bmp);f1,f2=freqspace(size(IA),meshgrid);%D=100/size(IA,1);D0=0.1;D1=0.4;r=sqr

26、t(f1.2+f2.2);Hd=zeros(size(IA);Hd(r=D0 & r(i,j)=D1 Hd(i,j)=(D1-r(i,j)/(D1-D0); end endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong)

27、; 梯形高通滤波器IA=imread(lena.bmp);IB=imread(babarra.bmp);f1,f2=freqspace(size(IA),meshgrid);%D=100/size(IA,1);D0=0.1;D1=0.4;r=sqrt(f1.2+f2.2);Hd=ones(size(IA);Hd(r=D0 & r(i,j)=D1 Hd(i,j)=(D0-r(i,j)/(D0-D1); end endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplo

28、t(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 用其他方法编写的理想低通、理想高通、Butterworth低通、同态滤波程序理想低通i1=imread(lena.bmp);i2=imnoise(i1,salt & pepper,0.1);f=double(i2);k=fft2(f);g=fftshift(k);N1,N2=size(g);d0=50;u0=floor(N1/2)+1;v0=flo

29、or(N2/2)+1;for i=1:N1 for j=1:N2 d=sqrt(i-u0)2+(j-v0)2); if d=d0 h=1; else h=0; end y(i,j)=g(i,j)*h; endendy=ifftshift(y);E1=ifft2(y);E2=real(E1);figuresubplot(2,2,1),imshow(uint8(i1);subplot(2,2,2),imshow(uint8(i2);subplot(2,2,3),imshow(uint8(E2); 理想高通i1=imread(lena.bmp);i2=imnoise(i1,salt & pepper,0.1);f=double(i2);k=fft2(f);g=fftshift(k);N1,N2=size(g)