实验名称：图像退化与复原实验目的1. 了解光电图像的退化原因；2. 掌握和理解基本的噪声模型，并能对图像进行加噪处理；3. 了解点扩展函数（PSF）与光学传递函数（OTF）的关系，熟悉几种经典的退化模型的模拟试验和OTF估计方法；4. 熟悉和掌握几种经典的图像复原方法及其基本原理；5. 能熟练利用MATLAB或C/C++工具进行图像的各种退化处理，并能编程实现退化图像的复原。
实验原理光电成像系统出现图像退化的过程是复杂多变的，为了研究的需要，通常情况下都把退化简化为化为一个线性移不变过程，见下图1所示。
f（工,v）I I 厂、冒（工J）了=0退化函数H = --------------- U + t——复原滤波器I------- 曝声V I ------------------- II "（3）I II 」I I! 降质过程I 屯原!图i光电图像退化与复原原理图因此，在空域中退化过程可以表示如下:g (x,y) = f (x,y) * h(x,y) + h(x,y)只有加性噪声不存在情况下，退化过程可以模型化如下表达式:g(x,y) = f (x,y) + h(x,y)其频域表达式为:G(u,v) =F(u,v)+N(u, v)针对这种退化图像的复原，除了周期噪声以外，通常都可以采用空间域滤波的方法进行图像复原，此时图像复原与图像增强几乎是没有区别的。
常见的空间 域滤波方法有均值滤波器和统计排序滤波器。
当退化图像存在线性移不变退化时， 图像的复原不能采用简单空间域滤波器 来实现，要实现线性移不变退化图像的复原， 必须知道退化系统的退化函数，即 点扩展函数h(x,y)。
在点扩展函数已知的情况下，常见图像复原方法有逆滤波 和维纳滤波两种。
在考虑噪声的情况下，逆滤波的原理可以表示如下: F (")二^^= F(u,v)+^) H(u,v) ' / H(u,v) 通常情况下，N(u,v)是未知的，因此即使知道退化模型也不能复原图像 此外，当H (u,v )的任何元素为零或者值很小时，N (u,v )/H (u,v )的比值决定 着复原的结果，从而导致图像复原结果出现畸变。
对丁这种情况， 通常采用限制 滤波频率使其难以接近原点值，从而减少遇到零值的可能性。
维纳滤波则克服了逆滤波的缺点，其数学模型表示如下： F(u ,v)顼|H *(u,v)「 2(u,v)H (u,v) + S h(u,v)/S f(u,v)G(u,v)然而，为退化图像的功率谱很少是已知的，因此常常用下面表达式近似:F(u,v) = & 1 |H *(u,v)|2:扩声 |H (u,v)|2+ kG(u,v)因此，本实验的内容就是利用上述经典图像复原的原理，对降质退化图像进 行复原。
四.实验步骤本次实验主要包括光电图像的退化模型和复原方法实现两大部分内容。
(一)图像的退化图像 1、大气湍流的建模1) 湍流引起图像退化的光学传递函数(OTF生成。
已知湍流退化模型的OTF表达式如下：H (u,v) = exp[-k(u2 + v2)5/6 ] (7)其中，k为一个常数，反映了大气湍流的严重程度。
(u,v)分别代表了(x,y )方向的频率坐标。
为了生成中心化的OTF可以考虑将式(7)改写为：, 2 2\5/6H (u,v)= exp[-k((u- M /2) + (v- N/2)2) ] (8)其中，M , N为图像的长和宽。
2) 读入一幅灰度图像，设定式(8)中k= 0.0025 ,进行退化试验。
分别显示原始图像、退化模型和退化图像。
3) 设定k= 0.001、0.00025重复上一步的试验。
读入原始图像对原图形进行灰度处理显示原图像将上述结果进行fft处理得到FP设计退化湍流模型为H显示传递函数由FP与H进行相关处理,得到退化图像结束图2大气湍流的退化过程2、运动模糊的图像退化试验1) 匀速直线运动引起图像退化的光学传递函数( OTF)生成。
已知相机匀速直线运动的OTF表达式如下:H (u,v) = ------- T------ sin[p(ua+ vb)]e- jp(ua+vb)p(ua + vb)其中，T为相机曝光时间，a,b分别表示（x,y）方向的速度；（u,v）分别对应（x, y）方向的频率坐标。
2）读入一幅灰度图像，设定式（9）中T = 1.0, a=b=0.1,编写MATLAB代码进行模糊退化试验。
要求分别显示原始图像、退化模型和退化图像。
3）设定不同的值，a, b值，重复上一步的试验。
4）利用数字显微镜或其它图像采集设备，移动物体过程中，采集图像。
读入原始图像对原图形进行灰度处理显示原图像（二）图像复原试验1、逆滤波将上述结果进行fft处理得到FP设计运动模糊模型为H由FP与H进行相关处理,得到退化图像结束图3运动模糊的图像退化显示传递函数显示退化图像1）根据试验（一）设计一幅退化图像（包括噪声污染+模糊退化两部分），其中模糊退化可选高斯模糊、大气湍流模糊或运动模糊（方向可任意指定，如10度、20度、45度等），噪声模型可自行设定。
2) 利用MATLAB 编程实现利用全逆滤波方法对退化图像的复原。
要求在同一个窗口下显示原始退化图像、复原结果及复原结果与理想图像的差值图共 3个图，并对复原 结果进行必要的分析。
逆滤波复原公式如下：其中，G(u,v)为退化图像的傅立叶变换，H(u,v)为退化系统的光学传递函数(OTF 。
3) 伪逆滤波：为了防止逆滤波中H(u,v)过小，使得复原后的图像数据过大和放大噪声，可采用频谱半径(阈值)限制下的逆滤波方法，即其中，R 为中心化频谱 H(u, v)中某点到原点(零频)的距离或半径。
另一种替代方法是直其中，(T 为一个阈值，用丁限制频谱的幅度值。
这种方法被称为伪逆滤波。
实验要求利用式(11)方式的伪逆滤波重复实验步骤内容 2)所涉及的图像。
2、Wiener 滤波1)针对以上逆滤波设计的退化图，编程实现利用Wiener 滤波对其进行复原。
滤波原理如下：F (u,v)=|5虫冬 G(u,v)U (u ,v)H(u,v) + kF(u,v) 二G(u,v)H(u,v)(10)R ?2V+2 接限制H (u,v)的值，即P(u,v) =0,H (u, v) £ s H (u, v) > s(12)(13)其中，G(u,v)为退化图像的傅立叶变换，H(u,v)为退化系统的光学传递函数(OTF, k 为一个与信噪比有关的调节因子。
要求在同一个窗口下显示理想图像(退化前)、退化图像、复原结果等共3个图，并对复原结果进行必要的分析。
2)改变k 值，重复试验内容1 )。
以上应根据原理自行编写代码，不允许直接调用MATLAfi 带的deconvwnr()函数分别显示三种 方式下的复原结果F1结束图4全逆，伪逆，wiener 滤波复原过程五.实验结果及分析1、大气湍流的建模显示原始 退化图像将复原结果F1与理想图像做减法，得到差值F2分别显示三种 方式下的差值F2分析：由上述结果可知，大气湍流会使图像变得模糊，而k值越大，其模糊效果越明显。
2、运动模糊的图像退化试验分析：由上述结果可知，随着a, b的值变大，图像模糊变得明显，人眼看起来好像是由丁运动速度过快造成的模糊。
3、图像复原试验分析：图10,图11,图12分别为全逆，伪逆，wiener对运动模糊（a=b=0.1）滤波的结果，从中可以看出wiener是三者中对运动模糊复原效果最好的滤波方式，且wiener中k值越小复原效果越好。
六.实验心得体会和建议心得体会：通过这次实验使我了解了图像退化的原因，以及相关的退化模型，并学会以matlab为平■台利用退化模型对图像进行退化处理以及退化图像的复原处理。
建议：可以要求利用C或C+辿行图像的退化与复原处理。
七.程序源代码% title : atmosphere% explain :本程序利用大气湍流模型对理想图像进行退化f=imread('3.jpg');figure(1)subplot(131),imshow(f),title('f=rgb2gray(f);Fp=fft2(f);原始图彳')[m,n]=size(f);% 绘制网格点[v,u]=meshgrid(1:n,1:m);u=u-floor(m/2);v=v-floor(n/2);k=0.00025;Duv=u.A2+v.A2;H=exp(-k.*Duv.A(5/6));G=H.*fftshift(Fp);f1=abs(ifft2(G));nchar = num2str(k);ltext = strcat('k=', nchar);% 标题注释subplot(132),imshow(H),title([' subplot(133),imshow(f1,[]),title(' 传递函数',ltext]);退化图像')；% title : move% explain :本程序利用运动模糊模型对理想图像进行退化f=imread('3.jpg');figure(1)subplot(131),imshow(f),title('f=rgb2gray(f);原始图彳')[m,n]=size(f);[v,u]=meshgrid(1:n,1:m);u=u-floor(m/2);v=v-floor(n/2);T=1.0;a=0.3;b=0.3;% a=0.1,b=0.1;% a=0.01,b=0.01;z=pi*(u*a+v*b)+eps;H=T./z.*sin (z) .*exp(-1j*z);Fp=fft2(f);G=H.*fftshift(Fp);f1=abs(ifft2(G));nchar = num2str(a);ltext = strcat('a=b=', nchar);subplot(132),imshow(H),title([' 传递函数',ltext]);subplot(133),imshow(f1,[]),title(' 退化图像');% title : recovery% explain :本程序利用运动模糊对加入椒盐噪声的理想图像进行退化，之后分别用全逆滤波%伪逆滤波，wiener滤波对设计的退化图进行处理，观察三种滤波的复原效果。
clc,clear all,close all;f=imread('4.jpg');f=rgb2gray(f);figure,imshow(f),title(' 原始图像')f1=imnoise(f,'salt & pepper',0.02);% 加入椒盐噪声后的图像[m,n]=size(f);[v,u]=meshgrid(1:n,1:m);% 画网格点u=u-floor(m/2);v=v-floor(n/2);%==================================================================== %运动模糊%==================================================================== % T=1.0,a=0.1,b=0.1;% T=1.0,a=0.01,b=0.01;T=1.0; a=0.1 ; b=0.1;z=pi*(u*a+v*b)+eps;H=T./z.*sin (z) .*exp(-1j*z);G=H.*fftshift(fft2(f1)); % 退化图像的频域 f2=real(ifft2(ifftshift(G)));% 加了椒盐噪声的退化图像 G1=fftshift(fft2(f2));% 求原始退化图像的频域表示%====================================================================以逆滤波%==================================================================== F=G1./H;F1=real(ifft2(ifftshift(F)));以反中心化，反傅立叶变换取实部得到复原结果图(包含噪声) F2=F1-double(f);%复原结果与理想图像的差值figure,subplot(131),imshow(f2,[]),title(' 原始退化图像'); subplot(132),imshow(F1,[]);title(' 复原结果'); subplot(133),imshow(F2,[]);title(' 差值图');%=======================%»逆滤波%======================= Duv=u.A 2+v.A 2; if Duv.A(1/2)<=1000 P=1./H; else P=0; endF=G1.*P;F1=real(ifft2(ifftshift(F))); F2=F1-double(f); figure,subplot(131),imshow(f2,[]),title(' subplot(132),imshow(F1,[]);title('subplot(133),imshow(F2,[]);title(' %==================================================================== %wiener 滤波%==================================================================== k1=0.000025;F=((1./H).*(abs(H)A2)./((abs(H)A2)+k1)).*G1; F1=real(ifft2(ifftshift(F)));% 反中心化，反傅立叶变换取实部得到复原结果图(包含噪声)F2=F1-double(f); figure,原始退化图像'); 复原结果'); 差值图')；subplot(131),imshow(f2,[]),title(' subplot(132),imshow(F1,[]);title(' subplot(133),imshow(F2,[]);title('八.思考题1. 简要叙述图像退化的原因。
