数字图像处理上机作业五

1、数字图像第四讲作业1.设计一个程序对受到高斯白噪声及椒盐噪声干扰的图像进行3x3,5x5邻域的平均平滑以及中值滤波. (添加噪声参看imnoise 函数, 空域卷积可用imfilter2函数实现。分析:1.邻域平均平滑可以采用imfilter函数,选择正确的卷积核就可以进行相应的邻域平均平滑操作了。3x3的卷积核为:H1=1/8*1 1 11 0 11 1 1;5x5的卷积核为:H2=1/24*1 1 1 1 11 1 1 1 11 1 0 1 11 1 1 1 11 1 1 1 1 ;2.中值平滑可以先编写中值平滑子函数zhongzhi(,然后在主函数中调用即可。以3*3中值平滑为例来分析其

2、操作过程,3*3中值平滑就是将以各项素为中心的9个像素值的中间值作为平滑后的新的像素值赋给该像素。因此可以通过I(i-1:i+1,j-1:j+1得到对应于I(i,j点的九个像素值,然后在由median函数可求出这九个值的中值,赋给新矩阵的(i,j点即可。注意I(i-1:i+1,j-1:j+1操作可能会有i-1=0,j-1=0或i+1、j+1大于矩阵最大行列数的情况,从而出现错误。在这里我的处理是在I矩阵的外围补上一圈0,即出现上述情况时像素值以0来代替。具体代码为:I0=zeros(m+2,n+2;for i=2:m+1for j=2:n+1I0(i,j=I(i-1,j-1;endend同理,

3、5*5的中值平滑也可以同样操作,只不过是在外围补上两圈零而已。代码及注释如下:主函数:clearI = imread('Lenna.bmp'J=imnoise(I,'gaussian'K=imnoise(I,'salt & pepper'%H1为3*3邻域平滑的卷积核,H2为5*5邻域平滑的卷积核H1=1/8*1 1 11 0 11 1 1;H2=1/24*1 1 1 1 11 1 1 1 11 1 0 1 11 1 1 1 11 1 1 1 1 ;J1=imfilter(J,H1; %高斯白噪声的3*3邻域平滑J2=imfilter(J

4、,H2; %高斯白噪声的5*5邻域平滑K1=imfilter(K,H1; %椒盐噪声的3*3邻域平滑K2=imfilter(K,H2; %椒盐噪声的5*5邻域平滑J3=zhongzhi(J,3; %高斯白噪声的3*3中值平滑J4=zhongzhi(J,5; %高斯白噪声的5*5中值平滑K3=zhongzhi(K,3; %椒盐噪声的3*3中值平滑K4=zhongzhi(K,5; %椒盐噪声的5*5中值平滑subplot(131;imshow(J;title('高斯白噪声'subplot(132;imshow(J1;title('高斯白噪声的3*3邻域平滑'subp

5、lot(133;imshow(J2;title('高斯白噪声的5*5邻域平滑'figuresubplot(131;imshow(J;title('高斯白噪声'subplot(132;imshow(J3;title('高斯白噪声的3*3中值平滑'subplot(133;imshow(J4;title('高斯白噪声的5*5中值平滑'figuresubplot(131;imshow(K;title('椒盐噪声'subplot(132;imshow(K1;title('椒盐噪声的3*3邻域平滑'subplo

6、t(133;imshow(K2;title('椒盐噪声的5*5邻域平滑'figuresubplot(131;imshow(K;title('椒盐噪声'subplot(132;imshow(K3;title('椒盐噪声的3*3中值平滑'subplot(133;imshow(K4;title('椒盐噪声的5*5中值平滑'中值平滑子函数zhongzhi(如下:function J=zhongzhi(I,km,n=size(I;if k=3 %3*3的中值平滑I0=zeros(m+2,n+2;for i=2:m+1for j=2:n+1I

7、0(i,j=I(i-1,j-1; %将到操作的图像矩阵I外围不上0endendfor i=2:m+1for j=2:n+1a=I0(i-1:i+1,j-1:j+1;b=a(1:9; %将3*3的矩阵化成1*9的矩阵,便于median操作J(i-1,j-1=median(b; %取中值,保存为平滑后矩阵J的i-1行、j-1列endendelse k=5 %5*5的中值平滑I0=zeros(m+4,n+4;for i=3:m+2for j=3:n+2I0(i,j=I(i-2,j-2;endendfor i=3:m+2for j=3:n+2a=I0(i-2:i+2,j-2:j+2;b=a(1:25;

8、J(i-2,j-2=median(b;endendendJ=uint8(J;运行结果如下:1加高斯白噪声后图像,及3*3、5*5邻域平滑 2加高斯白噪声后图像,及3*3、5*5中值滤波 3加椒盐噪声后图像,及3*3、5*5邻域平滑 4加椒盐噪声后图像,及3*3、5*5中值平滑 结论:平滑滤波和中值滤波对噪声都有一定的抑制作用,且阶数越高滤波效果越好,中值滤波对椒盐噪声的抑制效果特别明显,中值滤波效果比平滑滤波好一些,轮廓比较清晰。另外中值平滑在matlab 中有现成的函数medfilt2,经过与medfilt2函数对比,已验证zhongzhi(效果与medfilt2完全相同。2. 设计一个程序

9、对受到高斯白噪声及椒盐噪声干扰的图像在频域内分别采用理想低通和2阶butterworth 滤波器进行平滑处理.分析:1.在频域进行低通滤波,可以先通过fft2得到图形的频响,然后通过fftshift将零频点移到中心位置,再将频率响应的相应点利用数组乘法乘以传递函数即可。最后再由频域转换回空域即可。2.理想低通的传递函数为>+=022022,0,1,(D v u D v u v u H在程序中我用I0来表示,I0的生成过程如下:for i=1:mfor j=1:nif (i-i02+(j-j02<=602I0(i,j=1;elseI0(i,j=0;endendend3巴特沃斯的传递函

10、数为:n D v u D v u H 20,(12(11,(-+=在程序中我用H 来表示,H 的生成过程如下:for i=1:mfor j=1:nd=sqrt(i-i02+(j-j02;H(i,j=1/(1+0.414*(d/d0(2*k;endend4. 将频率响应的相应点利用数组乘法乘以传递函数就可对频响进行低通滤波处理,如J2.*H ,注意J2右下方有一点,表示数组乘法而不是矩阵乘法。最后由频域还原回空域,就得滤波后图像。程序及注释如下:clearI = imread('Lenna.bmp'J=imnoise(I,'gaussian' %加高斯白噪声K=i

11、mnoise(I,'salt & pepper'%加椒盐噪声J1=fft2(single(J;J2=fftshift(J1;K1=fft2(single(K;K2=fftshift(K1;m,n=size(I;i0=round(m+1/2;j0=round(n+1/2;%*理想低通*for i=1:mfor j=1:nif (i-i02+(j-j02<=602I0(i,j=1;elseI0(i,j=0;endendendJ3=ifft2(fftshift(J2.*I0;K3=ifft2(fftshift(K2.*I0;subplot(121;imshow(J;ti

12、tle('高斯白噪声'subplot(122;imshow(uint8(J3;title('高斯白噪声的理想低通滤波'figuresubplot(121;imshow(K;title('椒盐噪声'subplot(122;imshow(uint8(K3;title('椒盐噪声的理想低通滤波'%*二阶巴特沃斯滤波*d0=125; %截止频率为125Hzk=2; %二阶for i=1:mfor j=1:nd=sqrt(i-i02+(j-j02;H(i,j=1/(1+0.414*(d/d0(2*k;endendJ4=ifft2(fftsh

13、ift(J2.*H;K4=ifft2(fftshift(K2.*H;figuresubplot(121;imshow(J;title('高斯白噪声'subplot(122;imshow(uint8(J4;title('高斯白噪声的二阶巴特沃斯低通滤波' figuresubplot(121;imshow(K;title('椒盐噪声'subplot(122;imshow(uint8(K4;title('椒盐噪声的二阶巴特沃斯低通滤波' 运行结果如下: 结论:理想低通滤波后图像会有振铃现象,而二阶巴特沃斯滤波后,基本上没有振铃现象。3.

14、用egde函数提取一幅图像的边缘(sobel算子,canny算子,prewitt算子,LOG算子分析:直接调用edge即可。程序及注释入下:clearI = imread('Lenna.bmp'BW1 = edge(I,'sobel'BW2 = edge(I,'canny'BW3 = edge(I,'prewitt'BW4 = edge(I,'log'figure;imshow(BW1;title('sobel'figure; imshow(BW2;title('canny'figu

15、re; imshow(BW3;title('prewitt'figure; imshow(BW4;title('log'运行结果如下: 4根据提供的数据实现CT图像的重建. data的列向量是0180度的ct扫描数据(投影数据。见附件Data.mat分析:可以先对data各列求fft,求出N个方向上投影集合的傅立叶变换即不同数字图像第四讲作业 的F。求出|R|·F(R,乘积的ifft变换B。B每一列的值扩展成为一张图，即各 行值相同各列为B，程序表示为f=B(:,n*ones(1,a，然后对各图进行相应的旋转 变换i=imrotate(f,n。最后将各图叠加就完成了重建。 代码及注释如下： subplot(1,2,1;imagesc(data; F=fft(data; %对 data 每一列做 fft 变换 a,b=size(F; H=zeros(a,b; H=abs(-(a-1/2:(a-1/2'*ones(1,b; %H 为 a*b 的矩阵，每一行的值相同; %|R|·F(R, F=F.*H; B=abs(ifft(F; subplot(1,2,2;imagesc(B; f=zeros(a,a; D=zeros