《数字图像处理上机作业五.doc》由会员分享，可在线阅读，更多相关《数字图像处理上机作业五.doc（12页珍藏版）》请在三一文库上搜索。
1、数字图像第四讲作业 数字图像第四讲作业 1. 设计一个程序对受到高斯白噪声及椒盐噪声干扰的图像进行3x3,5x5邻域的平均平滑以及中值滤波. (添加噪声参看imnoise函数, 空域卷积可用imfilter2函数实现)。 分析：1.邻域平均平滑可以采用imfilter函数，选择正确的卷积核就可以进行相应的邻域平均平滑操作了。3x3的卷积核为： H1=1/8*1 1 1 1 0 1 1 1 1;5x5的卷积核为：H2=1/24*1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 ; 2.中值平滑可以先编写中值平滑子函数zhongzhi()，然后在主
2、函数中调用即可。以3*3中值平滑为例来分析其操作过程，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+1 for j=2:
3、n+1 I0(i,j)=I(i-1,j-1); end end 同理，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 1 1 0 1 1 1 1;H2=1/24*1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 ;J1=imfilter(J,H1); %高斯白噪声的3*3邻域平
4、滑J2=imfilter(J,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邻
5、域平滑);subplot(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邻域平滑);subplot(133);imshow(K2);title(椒盐噪
6、声的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,k)m,n=size(I);if k=3 %3*3的中值平滑 I0=zeros(m+2,n+2);for i=2:m+1 for j=2:n+1 I0(i,j)=I(i-1,j-1); %将到操作的图像矩阵I外围不上0 endendfo
7、r i=2:m+1 for j=2:n+1 a=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+2 for j=3:n+2 I0(i,j)=I(i-2,j-2); end end for i=3:m+2 for j=3:n+2 a=I0(i-2:i+2,j-2:j+2); b=a(1:25); J(i-2,j-2)=media
8、n(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).理想低通的传递函数为 在程序中我用I0来表示，I0的生成过程如下： for i=1:m for j=1:n if (i-i0)2+(j-j0)2=602 I0(i,j)=1; else I0(i,j)=0; end endend3)巴特沃斯的传递函数为：在程序中我用H来表示，H的生成过程如下：for i=1:m fo
10、r j=1:n d=sqrt(i-i0)2+(j-j0)2); H(i,j)=1/(1+0.414*(d/d0)(2*k); end end 4). 将频率响应的相应点利用数组乘法乘以传递函数就可对频响进行低通滤波处理，如J2.*H，注意J2右下方有一点，表示数组乘法而不是矩阵乘法。最后由频域还原回空域，就得滤波后图像。程序及注释如下：clearI = imread(Lenna.bmp);J=imnoise(I,gaussian); %加高斯白噪声K=imnoise(I,salt & pepper);%加椒盐噪声J1=fft2(single(J);J2=fftshift(J1);K1=fft2
11、(single(K);K2=fftshift(K1);m,n=size(I);i0=round(m+1)/2);j0=round(n+1)/2);%*理想低通*for i=1:m for j=1:n if (i-i0)2+(j-j0)2=602 I0(i,j)=1; else I0(i,j)=0; end endendJ3=ifft2(fftshift(J2.*I0);K3=ifft2(fftshift(K2.*I0);subplot(121);imshow(J);title(高斯白噪声)subplot(122);imshow(uint8(J3);title(高斯白噪声的理想低通滤波)figu
12、resubplot(121);imshow(K);title(椒盐噪声)subplot(122);imshow(uint8(K3);title(椒盐噪声的理想低通滤波)%*二阶巴特沃斯滤波*d0=125; %截止频率为125Hzk=2; %二阶for i=1:m for j=1:n d=sqrt(i-i0)2+(j-j0)2); H(i,j)=1/(1+0.414*(d/d0)(2*k); end endJ4=ifft2(fftshift(J2.*H);K4=ifft2(fftshift(K2.*H);figuresubplot(121);imshow(J);title(高斯白噪声)subpl
13、ot(122);imshow(uint8(J4);title(高斯白噪声的二阶巴特沃斯低通滤波)figuresubplot(121);imshow(K);title(椒盐噪声)subplot(122);imshow(uint8(K4);title(椒盐噪声的二阶巴特沃斯低通滤波)运行结果如下：结论：理想低通滤波后图像会有振铃现象，而二阶巴特沃斯滤波后，基本上没有振铃现象。3. 用egde函数提取一幅图像的边缘（sobel算子，canny算子，prewitt算子， LOG算子）分析：直接调用edge即可。程序及注释入下：clearI = imread(Lenna.bmp);BW1 = edge(
14、I,sobel);BW2 = edge(I,canny);BW3 = edge(I,prewitt);BW4 = edge(I,log);figure;imshow(BW1);title(sobel)figure; imshow(BW2);title(canny)figure; imshow(BW3);title(prewitt)figure; imshow(BW4);title(log)运行结果如下：4 根据提供的数据实现CT图像的重建. data的列向量是0180度的ct扫描数据(投影数据)。见附件Data.mat分析：可以先对data各列求fft，求出N个方向上投影集合的傅立叶变换即不同
15、的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的矩阵，每一行的值相同;F=F.*H; %|R|F(R,)B=abs(ifft(F);subplot(1,2,2);imagesc(B);f=zeros(a,a);D=zeros(a,a);for n=1:b f=B(:,n)*ones(1,a); i=imrotate(f,n); %旋转n度 p,q=size(i); f=i(floor(p/2-a/2+1):floor(p/2+a/2),floor(q/2-a/2+1):floor(q/2+a/2); %在旋转后的图形中截取和D同样大小的一部分图象 D=D+f;endD=D/180;figure;imagesc(D);运行结果如下：重建图： Page 12 of 12