数字图像处理上机作业五12页word资料.docx
《数字图像处理上机作业五12页word资料.docx》由会员分享，可在线阅读，更多相关《数字图像处理上机作业五12页word资料.docx（10页珍藏版）》请在冰豆网上搜索。
数字图像处理上机作业五12页word资料
数字图像第四讲作业
1.设计一个程序对受到高斯白噪声及椒盐噪声干扰的图像进行3x3,5x5邻域的平均平滑以及中值滤波.（添加噪声参看imnoise函数,空域卷积可用imfilter2函数实现）。
分析：
1.邻域平均平滑可以采用imfilter函数，选择正确的卷积核就可以进行相应的邻域平均平滑操作了。
3x3的卷积核为：
H1=1/8*[111
101
111];
5x5的卷积核为：
H2=1/24*[11111
11111
11011
11111
11111];
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）;
fori=2:
m+1
forj=2:
n+1
I0（i,j）=I（i-1,j-1）;
end
end
同理，5*5的中值平滑也可以同样操作，只不过是在外围补上两圈零而已。
代码及注释如下：
主函数：
clear
I=imread（'Lenna.bmp'）;
J=imnoise（I,'gaussian'）;
K=imnoise（I,'salt&pepper'）;
%H1为3*3邻域平滑的卷积核，H2为5*5邻域平滑的卷积核
H1=1/8*[111
101
111];
H2=1/24*[11111
11111
11011
11111
11111];
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邻域平滑'）;
figure
subplot（131）;imshow（J）;title（'高斯白噪声'）;
subplot（132）;imshow（J3）;title（'高斯白噪声的3*3中值平滑'）;
subplot（133）;imshow（J4）;title（'高斯白噪声的5*5中值平滑'）;
figure
subplot（131）;imshow（K）;title（'椒盐噪声'）;
subplot（132）;imshow（K1）;title（'椒盐噪声的3*3邻域平滑'）;
subplot（133）;imshow（K2）;title（'椒盐噪声的5*5邻域平滑'）;
figure
subplot（131）;imshow（K）;title（'椒盐噪声'）;
subplot（132）;imshow（K3）;title（'椒盐噪声的3*3中值平滑'）;
subplot（133）;imshow（K4）;title（'椒盐噪声的5*5中值平滑'）;
中值平滑子函数zhongzhi（）如下：
functionJ=zhongzhi（I,k）
[m,n]=size（I）;
ifk==3%3*3的中值平滑
I0=zeros（m+2,n+2）;
fori=2:
m+1
forj=2:
n+1
I0（i,j）=I（i-1,j-1）;%将到操作的图像矩阵I外围不上0
end
end
fori=2:
m+1
forj=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列
end
end
elsek==5%5*5的中值平滑
I0=zeros（m+4,n+4）;
fori=3:
m+2
forj=3:
n+2
I0（i,j）=I（i-2,j-2）;
end
end
fori=3:
m+2
forj=3:
n+2
a=I0（i-2:
i+2,j-2:
j+2）;
b=a（1:
25）;
J（i-2,j-2）=median（b）;
end
end
end
J=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.设计一个程序对受到高斯白噪声及椒盐噪声干扰的图像在频域内分别采用理想低通和2阶butterworth滤波器进行平滑处理.
分析：
1）.在频域进行低通滤波，可以先通过fft2得到图形的频响，然后通过fftshift将零频点移到中心位置，，再将频率响应的相应点利用数组乘法乘以传递函数即可。
最后再由频域转换回空域即可。
2）.理想低通的传递函数为
在程序中我用I0来表示，I0的生成过程如下：
fori=1:
m
forj=1:
n
if（i-i0）^2+（j-j0）^2<=60^2
I0（i,j）=1;
else
I0（i,j）=0;
end
end
end
3）巴特沃斯的传递函数为：
在程序中我用H来表示，H的生成过程如下：
fori=1:
m
forj=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右下方有一点，表示数组乘法而不是矩阵乘法。
最后由频域还原回空域，就得滤波后图像。
程序及注释如下：
clear
I=imread（'Lenna.bmp'）;
J=imnoise（I,'gaussian'）;%加高斯白噪声
K=imnoise（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）;
%****************理想低通******************
fori=1:
m
forj=1:
n
if（i-i0）^2+（j-j0）^2<=60^2
I0（i,j）=1;
else
I0（i,j）=0;
end
end
end
J3=ifft2（fftshift（J2.*I0））;
K3=ifft2（fftshift（K2.*I0））;
subplot（121）;imshow（J）;title（'高斯白噪声'）
subplot（122）;imshow（uint8（J3））;title（'高斯白噪声的理想低通滤波'）
figure
subplot（121）;imshow（K）;title（'椒盐噪声'）
subplot（122）;imshow（uint8（K3））;title（'椒盐噪声的理想低通滤波'）
%*************二阶巴特沃斯滤波****************
d0=125;%截止频率为125Hz
k=2;%二阶
fori=1:
m
forj=1:
n
d=sqrt（（i-i0）^2+（j-j0）^2）;
H（i,j）=1/（1+0.414*（d/d0））^（2*k）;
end
end
J4=ifft2（fftshift（J2.*H））;
K4=ifft2（fftshift（K2.*H））;
figure
subplot（121）;imshow（J）;title（'高斯白噪声'）
subplot（122）;imshow（uint8（J4））;title（'高斯白噪声的二阶巴特沃斯低通滤波'）
figure
subplot（121）;imshow（K）;title（'椒盐噪声'）
subplot（122）;imshow（uint8（K4））;title（'椒盐噪声的二阶巴特沃斯低通滤波'）
运行结果如下：
结论：
理想低通滤波后图像会有振铃现象，而二阶巴特沃斯滤波后，基本上没有振铃现象。
3.用egde函数提取一幅图像的边缘（sobel算子，canny算子，prewitt算子，LOG算子）
分析：
直接调用edge即可。
程序及注释入下：
clear
I=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'）
运行结果如下：
4根据提供的数据实现CT图像的重建.data的列向量是0~180度的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的矩阵，每一行的值相同;
F=F.*H;%|R|·F（R,θ）
B=abs（ifft（F））;
subplot（1,2,2）;imagesc（B）;
f=zeros（a,a）;
D=zeros（a,a）;
forn=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;
end
D=D/180;
figure;imagesc（D）;
运行结果如下：
重建图：
希望以上资料对你有所帮助，附励志名3条：
1、积金遗于子孙，子孙未必能守；积书于子孙，子孙未必能读。
不如积阴德于冥冥之中，此乃万世传家之宝训也。
2、积德为产业，强胜于美宅良田。
3、能付出爱心就是福，能消除烦恼就是慧。
升级会员 