西交大数字图像处理第六次作业.docx

1、西交大数字图像处理第六次作业 数字图像处理第六次作业 姓名： 班级： 学号： 提交日期：2015年4月27日1、在测试图像上产生高斯噪声lena图-需能指定均值和方差；并用滤波器（自选）恢复图像；（1）问题分析：图像退化模型：图像退化过程被建模为一个退化函数和一个加性噪声项，对一幅输入图像f(x,y)进行处理，产生一副退化后的图像g(x,y)。给定g(x,y)和关于退化函数H的一些知识以及关于加性噪声项(x,y)的一些知识后，图像复原的目的就是获得原始图像的一个估计。如果H是一个线性的、位置不变的过程，那么空间域中的退化图像可由下式给出： g(x,y)=h(x,y)f(x,y)+(x,y)其中

2、，h(x,y)是退化函数的空间表示；符号”表示空间卷积.等价的频率域表示： G(u,v)=H(u,v)F(u,v)+N(u,v) 高斯噪声： 所谓高斯噪声是指它的概率密度函数服从高斯分布（即正态分布）的一类噪声。一个高斯随机变量z的PDF可表示为： P（z）= 其中z代表灰度，u是z的均值，是z的标准差。高斯噪声的灰度值多集中在均值附近。 算术均值滤波器：令Sxy表示中心在点(x,y)处，大小为mn的矩形子图像窗口的一组坐标。算术均值滤波器在Sxy定义的区域中计算被污染的图像g(x,y)的平均值。在点(x,y)处复原图像的值： 这个操作可以使用大小为mn的一个空间滤波器来实现，其所有系数均为其

3、值的1/mn。均值滤波器平滑一幅图像中的局部变化，虽然模糊了结果，但降低了噪声。中值滤波器： 统计排序滤波器是空间域滤波器，空间滤波器的响应是基于由该滤波器包围的图像区域中的像素值的排序。中值滤波器使用一个像素邻域中灰度级的中值来替代该像素值，即： 编程思路： 首先利用imnoise函数给图像添加高斯噪声，之后分别利用算术平均值滤波器和中值滤波器进行滤波并观察效果。（2）MATLAB函数： g=imnoise(f,type,parameters) 函数功能：使用函数imnoise来用噪声污染一幅图像、 调用格式： J = imnoise(I,type) J = imnoise(I,type,p

4、arameters) 参数Type对应的噪声类型如下： gaussian高斯白噪声 localvar0均值白噪声 poisson泊松噪声 salt & pepper盐椒噪声 speckle乘性噪声（3）处理结果：1）添加高斯噪声 lena.bmp原始图像 lena添加高斯噪声（u=0，s2=0.01） lena添加高斯噪声（u=0，s2=0.05） lena添加高斯噪声（u=0，s2=0.1） lena添加高斯噪声（u=0.1，s2=0.01） lena添加高斯噪声（u=0.5，s2=0.01） 2）图像恢复（选取被均值为0，方差为0.01的高斯噪声污染的图像为例）利用算术均值滤波器恢复图像(

5、5x5模板) lena添加高斯噪声（u=0，s2=0.01） 算术均值滤波的结果 利用中值滤波器恢复图像 (5x5模板) lena添加高斯噪声（u=0，s2=0.01） 中值滤波的结果 （4）结果分析及总结：首先通过imnoise函数分别产生了被不同均值和方差的高斯噪声污染的图像。当高斯噪声均值不变为0时，随着方差增加，图像噪声越严重；当高斯噪声方差不变时，均值会影响到整个图像的灰度值，使整个图像变亮。与理论上均值和方差对图像的影响一致。分别使用算术均值滤波器和中值滤波器对加噪图像进行恢复。两种方法在一定程度上都可以降低噪声。算术均值滤波器降低噪声的同时也模糊了图像。2、推导维纳滤波器并实现下

6、边要求；(a) 实现模糊滤波器如方程Eq. (5.6-11).(b) 模糊lena图像：45度方向，T=1；(c) 再模糊的lena图像中增加高斯噪声，均值= 0 ，方差=10 pixels 以产生模糊图像；(d)分别利用方程 Eq. (5.8-6)和(5.9-4)，恢复图像。（1）问题分析：1）维纳滤波器的推导： 图像的退化模型为： （1） 其中，s(n1,n2)为原始图像，b(n1,n2)为退化函数，w(n1,n2)为噪声函数，x(n1,n2)为退化的图像。并假设s与w不相关，w为0均值的平稳随机过程。图像的复原模型为： （2） 其中，为恢复的图像，为恢复滤波器。 误差度量为： （3） 基

7、于正交性原理，若要求误差最小，则必有下式成立： （4） 将（3）式带入（4）式有： （5） 即 （6） 换元得： （7） 等式两端同时取傅里叶变换得： （8） 即 （9） 公式（8）中（10）公式（10）两端同时取傅里叶变换得： （11） 公式（8）中（12） 公式（12）两端同时取傅里叶变换： （13） 将（11）式和（13）式带入（8）式得 （14） 将符号化成与书中一致的表示 （15） 故表达式由下式给出 （16） 2）约束最小二乘方滤波 对于约束最小二乘方滤波，期望是找一个最小准则函数C，定义如下： （17） 其约束为 （18） 其中，是欧几里得向量范数，是未退化图像的估计。 这个最佳

8、问题在频率域中的解决由下面的表达式给出： （19）其中，是一个参数，必须对它进行调整以满足式（18）的条件，P(u,v)是函数 （20） 的傅里叶变换。（2）MATLAB函数：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用于指定在滤波过程中是使用“相关”还是

9、“卷积”。boundary_options用于处理边界充零问题，边界的大小由滤波器的大小确定。具体参数选项见下表：选项描述filtering_mode corr通过使用相关来完成，该值为默认。conv通过使用卷积来完成。boundary_options X输入图像的边界通过用值X（无引号）来填充扩展，其默认值为0。replicate图像大小通过复制外边界的值来扩展。symmetric图像大小通过镜像反射其边界来扩展。circular图像大小通过将图像看成是一个二维周期函数的一个周期来扩展size_optionsfull输出图像的大小与被扩展图像的大小相同same输出图像的大小与输入图像的大小相

10、同。这可通过将滤波掩模的中心点的偏移限制到原图像中包含的点来实现，该值为默认值。2）fspecial功能：fspecial函数用于建立预定义的滤波算子。用法：h = fspecial(type)h = fspecial(type，para)其中type指定算子的类型，para指定相应的参数；选项描述typeaverageaveraging filter为均值滤波，参数为hsize代表模板尺寸，默认值为3 3diskcircular averaging filter为圆形区域均值滤波，参数为radius代表区域半径，默认值为5.gaussianGaussian lowpass filter为高斯

11、低通滤波，有两个参数，hsize表示模板尺寸，默认值为3 3，sigma为滤波器的标准值，单位为像素，默认值为0.5.laplacianlaplacian filter为拉普拉斯算子，参数alpha用于控制算子形状，取值范围为0，1，默认值为0.2。logLaplacian of Gaussian filter为拉普拉斯高斯算子，有两个参数，hsize表示模板尺寸，默认值为3 3，sigma为滤波器的标准差，单位为像素，默认值为0.5.motionmotion filter为运动模糊算子，有两个参数，表示摄像物体逆时针方向以theta角度运动了len个像素，len的默认值为9，theta的默认

12、值为0；prewittPrewit thorizontal edge-emphasizing filter用于边缘增强，大小为3 3，无参数sobelSobel horizontal edge-emphasizing filter用于边缘提取，无参数unsharpunsharp contrast enhancement filter为对比度增强滤波器。参数alpha用于控制滤波器的形状，范围为0，1，默认值为0.2.（3）处理结果：（a）实现模糊滤波器如方程Eq. (5.6-11). 模糊滤波器的频域表达式为： 故实现该滤波器，只需先将输入图像进行傅里叶变换并移至图像中心，之后将图像的傅里叶变

13、换和模糊滤波器的傅里叶变换进行阵列相乘，将得到的结果经过傅里叶反变换返回到空间域即可实现该滤波器。具体程序见附件mohu_filter.m（b）模糊lena图像：45度方向，T=1；（a=0.1，b=0.1；T=1） lena.bmp原始图像 lena运动模糊的结果（a=0.1，b=0.1；T=1） lena.bmp原始图像 lena运动模糊的结果（调用MATLAB中的函数） （c）再模糊的lena图像中增加高斯噪声，均值= 0 ，方差=10 pixels 以产生模糊图像；lena运动模糊的结果（a=0.1，b=0.1；T=1） 添加高斯噪声的结果（均值为0，方差为0.01） lena运动模糊

14、的结果（MATLAB版） 添加高斯噪声的结果（MATLAB版） （d）分别利用方程 Eq. (5.8-6)和(5.9-4)，恢复图像。维纳滤波的结果： lena运动模糊+高斯噪声.bmp 维纳滤波结果（K=0.06） 维纳滤波结果（K=0.01） 维纳滤波结果（K=0.5） lena运动模糊+高斯噪声（MATLAB版） 维纳滤波结果（MATLAB版） 约束最小二乘方滤波的结果： lena运动模糊+高斯噪声（MATLAB版） 约束最小二乘滤波结果（MATLAB版） （3）结果分析及总结：首先分别通过自己编写的模糊函数和MATLAB中提供的imfilter和fspecial函数配合使用对图像len

15、a进行了模糊滤波。发现套用书上的公式图像是斜向下45度运动模糊，而MATLAB中函数模糊的结果是斜向上45度运动模糊，不过这不是重点可以接受，模糊的基本效果还是一致的。之后调用imnoise函数对两幅图像加入高斯噪声，得到第二问的结果。之后分别使用自己编写的函数和MATLAB中提供的deconvwnr函数进行维纳滤波。调用MATLAB中函数滤波后的图像得到了一定的改善，运动模糊的影响基本被消除，但噪声的影响仍然较大，导致图像质量下降；对于自己编写的维纳滤波函数，难点在于寻找令信噪比最大的K值，报告中显示了部分K值对应的滤波结果，其中K=0.06，为信噪比最大时的滤波结果，从结果看，视觉上的效果

16、并不是很理想，要想达到更好的效果可能需要寻找更加合适的K值。最后采用MATLAAB提供的deconvreg函数进行约束最小二乘方滤波。从滤波后的结果看，约束最小二乘方滤波得到了比维纳滤波更好的结果，尤其是对噪声的滤除。由于实在是没看懂书上的公式应该怎么实现，所以只好调用函数了。尤其实P(u,v)的维数和其他几个矩阵的维数不同，怎么遍历啊？附录【参考文献】1 冈萨雷斯.数字图像处理（第三版）北京：电子工业出版社，20112 周品.MATLAB数字图像处理 北京：清华大学出版社，20123 杨杰.数字图像处理及MATLAB实现 北京：电子工业出版社，2010【源代码】add_gaussian_no

17、ise.m(题目1添加高斯噪声)I=imread(lena.bmp);figure(1);imshow(I);title(源图像lena.bmp);imwrite(I,lena原始图像.bmp);I2=imnoise(I,gaussian,0.5,0.01);figure(2);imshow(I2);title(加入gaussian噪声后的lena.bmp);imwrite(I2,lena加入gaussian噪声后(u=0.5，s2=0.01).bmp);suanshujunzhi_filter.m(题目1算术均值滤波恢复图像)I=imread(lena_gaussian_noise.bmp)

18、;figure(1);imshow(I);title(lena加入gaussian噪声后(u=0，s2=0.01).bmp);imwrite(I,lena加入gaussian噪声后(u=0，s2=0.01).bmp);n=5; h=1/n2.*ones(n,n);I2=conv2(I,h,same); I2=uint8(I2);figure(2);imshow(I2); title(算术均值滤波的结果（5x5）);imwrite(I2,算术均值滤波的结果（5x5）.bmp);median_filter.m(题目1中值滤波恢复图像)I=imread(lena_gaussian_noise.bmp

19、);figure(1);imshow(I);title(lena加入gaussian噪声后(u=0，s2=0.01).bmp);imwrite(I,lena加入gaussian噪声后(u=0，s2=0.01).bmp);figure(2);n=5; a=ones(n,n);p=size(I); x1=double(I);x2=x1;for i=1:p(1)-n+1 for j=1:p(2)-n+1 c=x1(i:i+(n-1),j:j+(n-1); e=c(1,:); for u=2:n e=e,c(u,:); end mm=median(e); x2(i+(n-1)/2,j+(n-1)/2)

20、=mm; endendI2=uint8(x2);imshow(I2); title(中值滤波的结果（5x5）);imwrite(I2,中值滤波的结果（5x5）.bmp);mohu_filter.m(题目2运动模糊滤波器)I=imread(lena.bmp);figure(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:M for v=1:N H(u,v)=(T/(pi*(u*a+v

21、*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);figure(2);imshow(g);title(运动模糊化lena.bmp);imwrite(g,lena运动模糊的结果.bmp); Wiener_filter.m(题目2维纳滤波器自编版)I=imread(lena运动模糊+高斯噪声.bmp);figure(1);imshow(I);title(lena运动模糊+

22、高斯噪声);imwrite(I,lena运动模糊+高斯噪声.bmp);g=double(I); G=fft2(g); G=fftshift(G); M,N=size(G);a=0.1;b=0.1;T=1;K=0.0259;for u=1:M for v=1:N H(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/H(u,v)*(abs(H(u,v)2/(abs(H(u,v)2+K)*G(u,v); endendF=ifftshift(F);f=ifft2(F);f=256.*f./max(m

23、ax(f);f=uint8(real(f);figure(2);imshow(f);title(维纳滤波的将结果);imwrite(f,维纳滤波的结果(K=0.0259).bmp);wiener_filter_k.m(题目2维纳滤波器寻找信噪比最大的K值)I=imread(lena运动模糊+高斯噪声.bmp);figure(1);imshow(I);title(lena运动模糊+高斯噪声);imwrite(I,lena运动模糊+高斯噪声.bmp);g=double(I); G=fft2(g); G=fftshift(G); M,N=size(G);a=0.1;b=0.1;T=1;i=1;for

24、mat longfor k=0.01:0.01:0.11 for u=1:M for v=1:N H(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/H(u,v)*(abs(H(u,v)2/(abs(H(u,v)2+k)*G(u,v); end end F=ifftshift(F); f=ifft2(F); f=256.*f./max(max(f); f=uint8(real(f); figure; imshow(f); title(维纳滤波的结果); e=f-uint8(g); SNR(

25、i)=sum(sum(g.2)/sum(sum(e.2); i=i+1;endidx=find(max(SNR)wiener_filter_matlab.m(题目2维纳滤波器MATLAB版)I=imread(lena.bmp);H=fspecial(motion,50,45);I1=imfilter(I,H,circular,conv);figure(1);imshow(I1);title(运动模糊后的lena.bmp(角度为45度);imwrite(I1,lena运动模糊（调用matlab中的函数）.bmp);I2=imnoise(I1,gaussian,0,0.01);figure(2)i

26、mshow(I2);title(加噪并模糊的lena.bmp);imwrite(I2,lena运动模糊+高斯噪声（调用matlab中的函数0.bmp);figure(3);noise=imnoise(zeros(size(I),gaussian,0,0.01);NSR=sum(noise(:).2)/sum(im2double(I(:).2);I3=deconvwnr(I2,H,NSR);imshow(I3);title(维纳滤波的结果);imwrite(I3,lena维纳滤波的结果（调用MATLAB中的函数）.bmp);yueshuzuixiaercheng_filter_matlab.m（

27、题目2约束最小二乘滤波（MATLAB版）I=imread(lena.bmp);h=fspecial(motion,50,45);I1=imfilter(I,h,circular,conv);I2=imnoise(I1,gaussian,0,0.01);figure(1);imshow(I2);title(lena运动模糊+高斯噪声);imwrite(I2,lena运动模糊+高斯噪声（MATLAB版）.bmp);V=0.0001;NoisePower=V*prod(size(I);g,LAGRA=deconvreg(I1,h,NoisePower);figure(2);imshow(g);title(约束最小二乘滤波的结果（MATLAB版）);imwrite(g,约束最小二乘滤波的结果（MATLAB版）.bmp);