数字图像处理上机作业五

数字图像第四讲作业 Page 1 of 12数字图像第四讲作业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 1 1 1 1 1 11 1 0 1 11 1 1 1 1 1 1 1 1 1 ;2.中值平滑可以先编写中值平滑子函数 zhongzhi()，然后在主函数中调用即可。以 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+1for j=2:n+1I0(i,j)=I(i-1,j-1);endend同理，5*5 的中值平滑也可以同样操作，只不过是在外围补上两圈零而已。数字图像第四讲作业 Page 2 of 12代码及注释如下：主函数：clearI = imread(Lenna.bmp);J=imnoise(I,gaussian);K=imnoise(I,salt %H1 为 3*3 邻域平滑的卷积核，H2 为 5*5 邻域平滑的卷积核H1=1/8*1 1 11 0 11 1 1;H2=1/24*1 1 1 1 1 1 1 1 1 11 1 0 1 11 1 1 1 1 1 1 1 1 1 ;J1=imfilter(J,H1); %高斯白噪声的 3*3 邻域平滑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 邻域平滑);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(椒盐噪声的 5*5 邻域平滑 );figuresubplot(131);imshow(K);title(椒盐噪声);subplot(132);imshow(K3);title(椒盐噪声的 3*3 中值平滑 );subplot(133);imshow(K4);title(椒盐噪声的 5*5 中值平滑 );中值平滑子函数 zhongzhi()如下：数字图像第四讲作业 Page 3 of 12function J=zhongzhi(I,k)m,n=size(I);if k=3 %3*3 的中值平滑I0=zeros(m+2,n+2);for i=2:m+1for j=2:n+1I0(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);J(i-2,j-2)=median(b);endendendJ=uint8(J);运行结果如下：1）加高斯白噪声后图像，及 3*3、5*5 邻域平滑数字图像第四讲作业 Page 4 of 122）加高斯白噪声后图像，及 3*3、5*5 中值滤波3）加椒盐噪声后图像，及 3*3、5*5 邻域平滑4）加椒盐噪声后图像，及 3*3、5*5 中值平滑数字图像第四讲作业 Page 5 of 12结论：平滑滤波和中值滤波对噪声都有一定的抑制作用，且阶数越高滤波效果越好，中值滤波对椒盐噪声的抑制效果特别明显，中值滤波效果比平滑滤波好一些，轮廓比较清晰。另外中值平滑在 matlab 中有现成的函数medfilt2，经过与 medfilt2 函数对比，已验证 zhongzhi()效果与medfilt2 完全相同。2. 设计一个程序对受到高斯白噪声及椒盐噪声干扰的图像在频域内分别采用理想低通和 2 阶 butterworth 滤波器进行平滑处理 .分析：1). 在频域进行低通滤波，可以先通过 fft2 得到图形的频响，然后通过fftshift 将零频点移到中心位置， ，再将频率响应的相应点利用数组乘法乘以传递函数即可。最后再由频域转换回空域即可。2).理想低通的传递函数为 02,01),(DvuvuH在程序中我用 I0 来表示，I0 的生成过程如下：for i=1:mfor j=1:nif (i-i0)2+(j-j0)2=602I0(i,j)=1;elseI0(i,j)=0;endend数字图像第四讲作业 Page 6 of 12end3)巴特沃斯的传递函数为： nDvuvuH20),(1),(在程序中我用 H 来表示，H 的生成过程如下：for i=1:mfor j=1:nd=sqrt(i-i0)2+(j-j0)2);H(i,j)=1/(1+0.414*(d/d0)(2*k);end end4). 将频率响应的相应点利用数组乘法乘以传递函数就可对频响进行低通滤波处理，如 J2.*H，注意 J2 右下方有一点，表示数组乘法而不是矩阵乘法。最后由频域还原回空域，就得滤波后图像。程序及注释如下：clearI = imread(Lenna.bmp);J=imnoise(I,gaussian); %加高斯白噪声K=imnoise(I,salt %加椒盐噪声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-i0)2+(j-j0)2=602I0(i,j)=1;elseI0(i,j)=0;endendendJ3=ifft2(fftshift(J2.*I0);K3=ifft2(fftshift(K2.*I0);subplot(121);imshow(J);title(高斯白噪声)subplot(122);imshow(uint8(J3);title(高斯白噪声的理想低通滤波 )数字图像第四讲作业 Page 7 of 12figuresubplot(121);imshow(K);title(椒盐噪声)subplot(122);imshow(uint8(K3);title(椒盐噪声的理想低通滤波)%*二阶巴特沃斯滤波*d0=125; %截止频率为 125Hzk=2; %二阶for i=1:mfor j=1:nd=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(高斯白噪声)subplot(122);imshow(uint8(J4);title(高斯白噪声的二阶巴特沃斯低通滤波 )figuresubplot(121);imshow(K);title(椒盐噪声)subplot(122);imshow(uint8(K4);title(椒盐噪声的二阶巴特沃斯低通滤波)运行结果如下：数字图像第四讲作业 Page 8 of 12数字图像第四讲作业 Page 9 of 12结论：理想低通滤波后图像会有振铃现象，而二阶巴特沃斯滤波后，基本上没有振铃现象。3. 用 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)figure; imshow(BW3);title(prewitt)figure; imshow(BW4);title(log)数字图像第四讲作业 Page 10 of 12运行结果如下：4 根据提供的数据实现 CT 图像的重建. data 的列向量是 0180度的 ct 扫描数据( 投影数据 )。见附件 Data.mat分析：可以先对data各列求fft，求出N 个方向上投影集合的傅立叶变换即不同的F。求出 |R|F(R, )乘积的ifft变换B。B 每一列的值扩展成为一张图，即数字图像第四讲作业 Page 11 of 12各行值相同各列为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(