数字图像处理第六次作业摘要本次报告主要记录第六次作业中的各项任务完成情况。
本次作业以Matlab 2013为平台，结合matlab 函数编程实现对所给图像文件的相关处理：1.在测试图像上产生高斯噪声lena 图-需能指定均值和方差；并用滤波器（自选）恢复图像；2.实现下边要求；(a) 实现模糊滤波器如方程Eq. (5.6-11) ；(b) 模糊lena 图像：45度方向，T=1；(c) 在模糊的lena 图像中增加高斯噪声，均值= 0 ，方差=10 pixels 以产生退化图像；(d)分别利用方程 Eq. (5.8-6)和(5.9-4)，恢复图像。
以上任务完成后均得到了预期的结果。
1.在测试图像上产生高斯噪声lena 图-需能指定均值和方差；并用滤波器（自选）恢复图像； （1）实验原理与方法图像复原处理是建立在图像退化的数学模型基础上的，这个退化数学模型能够反映图像退化的原因。
图像的退化过程可以理解为施加于原图像上的运算和噪声两者联合作用的结果，图像退化模型如图1所示，可以表示为：g +(x,y )H [f (x,y )]n(x,y )f (x,y )*h(x,y )n(x,y )==+图1 图像退化模型高斯噪声是指它的概率密度函数服从高斯分布（即正态分布）的一类噪声。
一个高斯随机变量z 的PDF 可表示为：222(z u )P(z )]-=σ 其中z 代表灰度，u 是z 的均值，是z 的标准差。
高斯噪声的灰度值多集中在均值附近。
本文采用5×5模板的中值滤波器和高斯滤波器（σ=1.5） 作业四中已经介绍过，中值滤波器是使用一个像素邻域中灰度级的中值来替代该像素值，即，xy(s,t )S ˆf (x,y ){g(s,t )}median ∈=。
高斯滤波是一种根据高斯函数的形状来选择模板权值的线性平滑滤波方法，具体操作是：用一个模板（或称卷积）扫描图像中的每一个像素，用模板确定的邻域内像素的加权平均灰度值去替代模板中心像素点的值。
利用matlab 中imnoise 函数加入高斯噪声： g=imnoise(f,type,parameters) 调用格式：g = imnoise(I,type)g = imnoise(I,type,parameters) 参数Type 对应的噪声类型如下： 'gaussian'高斯白噪声 'localvar'0均值白噪声'poisson'泊松噪声'salt & pepper'盐椒噪声'speckle'乘性噪声滤波程序同作业四（2）处理结果原始图像lena.bmp加入gaussian噪声后的lena.bmp(u=0.5，s2=0.01)原始图像lena.bmp加入gaussian噪声后的lena.bmp(u=0，s2=0.01)原始图像lena.bmp加入gaussian噪声后的lena.bmp(u=0.5，s2=0.1)（3）结果分析通过imnoise 函数产生了被均值和方差可选的高斯噪声污染的图像。
当高斯噪声均值不变为0时，随着方差增加，图像噪声越严重；当高斯噪声方差不变时，均值会影响到整个图像的灰度值，使整个图像变亮。
与理论上均值和方差对图像的影响一致。
分别使用高斯滤波器和中值滤波器对加噪图像进行恢复。
两种方法在一定程度上都可以降低噪声。
高斯滤波器降低噪声的同时保存的图像细节更丰富，亮度比原噪声图像和中值滤波后图像暗更接近原始图像，中值滤波后图像亮度基本与原噪声图像相同。
2．实现下边要求：(a) 实现模糊滤波器如方程Eq. (5.6-11). (b) 模糊lena 图像：45度方向，T=1；(c) 在模糊的lena 图像中增加高斯噪声，均值= 0 ，方差=10 pixels 以产生退化图像；(d)分别利用方程 Eq. (5.8-6)和(5.9-4)，恢复图像；j (ua vb )TH(u,v )sin[(ua vb )]e (ua vb )-π+=π+π+（5.6-11）利用上式模糊lena 图像：45度方向，T=1，即使a=b=0.1，T=1。
对原始图像的图像矩阵做傅里叶变换并移至图像中心得到频域矩阵F ，使H 与F 相乘后反傅里叶变换到空域得到变换后图像。
维纳滤波综合了退化函数和噪声统计特性两个方面进行复原处理，其目标是寻找一个滤波器，使得复原后图像ˆf (x,y )与原始图像f (x,y )的均方误差最小：原始图像lena.bmplena 加入gaussian 噪声后(u=0.5，s 2=0.01).bmp中值滤波（5x5）高斯滤波5x5{}2ˆE f (x,y )f (x,y )min⎡⎤-=⎣⎦因此维纳滤波器又称为最小均方误差滤波器，在频率中用下式表达：221f H(u,v )ˆF(u,v )[]G(u,v )H(u,v )H(u,v )S (u,v )/S (u,v )η=+其中G(u,v)是退化图像的傅里叶变换，H(u,v)是退化函数。
2S (u,v )N(u,v )η= 为噪声功率谱，2f S (u,v )F(u,v )= 为未退化图像的功率谱。
式5.8-6为221H(u,v )ˆF(u,v )[]G(u,v )H(u,v )H(u,v )K=+属于维纳滤波 式5.9-4为22H *(u,v )ˆF(u,v )[]G(u,v )H(u,v )P(u,v )=+γ其中，γ是一个参数，必须对它进行调整以满足22ˆg Hf-=η的条件，属于约束最小二层方滤波利用以上两式恢复图像的流程与对模糊原始图像的流程相似：对原始图像的图像矩阵做傅里叶变换并移至图像中心得到频域矩阵F ，通过H 得到ˆF，使ˆF 与F 相乘后反傅里叶变换到空域得到变换后图像。
在实现 5.9-4时借助matlab 工具包以得到更好的效果。
用fspecial 和imnoise 函数得到45度方向，T=1的模糊lena 图像，并在此图像上实现维纳滤波和约束最小二乘方滤波。
1）imfilter功能：对任意类型数组或多维图像进行滤波。
用法：B = imfilter(A,H)B = imfilter(A,H,option1,option2,...)或写做g = imfilter(f, w, filtering_mode, boundary_options, size_options)其中，f 为输入图像，w 为滤波掩模，g 为滤波后图像。
filtering_mode 用于指定在滤波过程中是使用“相关”还是“卷积”。
boundary_options 用于处理边界充零问题，边界的大小由滤波器的大小确定。
2）fspecial功能：fspecial 函数用于建立预定义的滤波算子。
用法：h = fspecial(type) h = fspecial(type ，para)其中type 指定算子的类型，para 指定相应的参数。
（2）处理结果lena.bmp 原始图像运动模糊化lena.bmp运动模糊化lena.bmp模糊lena.bmp 加入高斯噪声(u=0，s 2=0.01)lena 运动模糊+高斯噪声维纳滤波的结果(K=0.06)(3)结果分析1.按照书上公式编写的模糊函数图像是斜向下45度运动模糊，matlab 函数是斜向上45度运动模糊的，公式的程序得到图像棱角比较分明边界比较明显。
2.使用自己编写的函数进行维纳滤波，难点在于寻找令信噪比最大的K 值，报告中显示了K=0.06时的滤波结果，从结果看，视觉上的效果并不是很理想，噪声依然很大，要想达到更好的效果可能需要寻找更加合适的K 值或者直接使用matlab 的deconvreg 函数实现。
3.最后采用MATLAAB 提供的deconvreg 函数进行约束最小二乘方滤波。
从滤波后的结果看，约束最小二乘方滤波得到了比维纳滤波更好的结果，噪声基本消除，图像变得模糊但是平滑。
附录： 参考文献：[1] Rafael C. Gonzalez., et al. 数字图像处理（第三版）, 电子工业出版社, 2011. [2] 周品. MATLAB 数字图像处理北京, 清华大学出版社, 2012源代码：1.img1.m[产生高斯噪声并用高斯滤波器和中值滤波器滤波]I=imread('lena.bmp'); figure(1); subplot(1,2,1)imshow(I);lena 运动模糊+高斯噪声约束最小二乘滤波的结果title('原始图像lena.bmp');imwrite(I,'原始图像lena.bmp');I2=imnoise(I,'gaussian',0.5,0.01);subplot(1,2,2)imshow(I2);title('加入gaussian噪声后的lena.bmp(u=0.5，s^2=0.01)');imwrite(I2,'加入gaussian噪声后的lena.bmp(u=0.5，s^2=0.01).bmp');figure(2);subplot(2,2,1)imshow(I);title('原始图像lena.bmp');subplot(2,2,2)imshow(I2);title('lena加入gaussian噪声后的（u=0.5，s^2=0.01).bmp');n=5;a=ones(n,n);p=size(I2);x1=double(I2);x2=x1;for i=1:p(1)-n+1for j=1:p(2)-n+1c=x1(i:i+(n-1),j:j+(n-1));e=c(1,:);for u=2:ne=[e,c(u,:)];endmm=median(e);x2(i+(n-1)/2,j+(n-1)/2)=mm;endendI3=uint8(x2);subplot(2,2,3)imshow(I3);title('中值滤波（5x5）');imwrite(I3,'中值滤波（5x5）.bmp');[I2,map]=imread('加入gaussian噪声后的lena.bmp(u=0.5，s^2=0.01).bmp'); k=1.5;Img=double(I2);n=5;n1=floor((n+1)/2);for i=1:nfor j=1:nb(i,j) =exp(-((i-n1)^2+(j-n1)^2)/(4*k))/(4*pi*k);endendImg1=conv2(Img,b,'same');d=uint8(Img1);subplot(2,2,4);imshow(d,map);title('高斯滤波5x5')2.img2.m[运动模糊+高斯噪声、维纳滤波、约束最小二乘滤波]I=imread('lena.bmp');figure(1);subplot(1,2,1)imshow(I);title('lena.bmp原始图像');imwrite(I,'lena原始图像.bmp');f=double(I);F=fft2(f);F=fftshift(F);[M,N]=size(F);a=0.1;b=0.1;T=1;for u=1:Mfor v=1:NH(u,v)=(T/(pi*(u*a+v*b)))*sin(pi*(u*a+v*b))*exp(-sqrt(-1)*pi*(u*a+v*b)); G(u,v)=H(u,v)*F(u,v);endendG=ifftshift(G);g=ifft2(G);g=256.*g./max(max(g));g=uint8(real(g));subplot(1,2,2);imshow(g);title('运动模糊化lena.bmp');imwrite(g,'lena运动模糊的结果.bmp');figure(2)subplot(1,2,1);imshow(g);title('运动模糊化lena.bmp');imwrite(g,' lena运动模糊的结果.bmp');I2=imnoise(g,'gaussian',0,0.01);subplot(1,2,2)imshow(I2);title('模糊lena.bmp加入高斯噪声(u=0，s^2=0.01)');imwrite(I2,'模糊lena.bmp加入高斯噪声(u=0，s^2=0.01).bmp');figure(3)I=imread('lena.bmp');h=fspecial('motion',50,45);I1=imfilter(I,h,'circular','conv');I2=imnoise(I1,'gaussian',0,0.01);subplot(1,2,1)imshow(I2);title('lena运动模糊+高斯噪声');imwrite(I2,'lena运动模糊+高斯噪声.bmp');g1=double(I2);G1=fft2(g1);G1=fftshift(G1);[M,N]=size(G1);a=0.1;b=0.1;T=1;K=0.06;for u=1:Mfor v=1:NH1(u,v)=(T/(pi*(u*a+v*b)))*sin(pi*(u*a+v*b))*exp(-sqrt(-1)*pi*(u*a+v*b)); F(u,v)=1/H1(u,v)*(abs(H1(u,v)))^2/((abs(H1(u,v)))^2+K)*G1(u,v);endendF=ifftshift(F);f=ifft2(F);f=256.*f./max(max(f));f=uint8(real(f));subplot(1,2,2)imshow(f);title('维纳滤波的结果(K=0.06)');imwrite(f,'维纳滤波的结果(K=0.06).bmp');I=imread('lena.bmp');h=fspecial('motion',50,45);I1=imfilter(I,h,'circular','conv');I2=imnoise(I1,'gaussian',0,0.01);figure(4);subplot(1,2,1)imshow(I2);title('lena运动模糊+高斯噪声');imwrite(I2,' lena运动模糊+高斯噪声.bmp');V=0.0001;NoisePower=V*prod(size(I));[g,LAGRA]=deconvreg(I1,h,NoisePower);subplot(1,2,2)imshow(g);title('约束最小二乘滤波的结果');imwrite(g,'约束最小二乘滤波的结果.bmp');（注：可编辑下载，若有不当之处，请指正，谢谢!）。