MatLab代码大全.docx
《MatLab代码大全.docx》由会员分享，可在线阅读，更多相关《MatLab代码大全.docx（104页珍藏版）》请在冰豆网上搜索。
MatLab代码大全
第2章图像获取
2.3.2二维连续傅里叶变换
例2.2
figure
（1）;%建立图形窗口1
[u,v]=meshgrid（-1:
0.01:
1）;%生成二维频域网格
F1=abs（sinc（u.*pi））;
F2=abs（sinc（v.*pi））;
F=F1.*F2;%计算幅度频谱F=|F（u,v）|
surf（u,v,F）;%显示幅度频谱，如图2.3（b）
shadinginterp;%平滑三维曲面上的小格
axisoff;%关闭坐标系
figure
（2）;%建立图形窗口2
F1=histeq（F）;%扩展F的对比度以增强视觉效果
imshow（F1）;%用图像来显示幅度频谱，如图2.3（c）
第3章图像变换
3.4.4二维FFT的MATLAB实现
例3.2简单图像及其傅里叶变换
MATLAB程序：
%建立简单图像d并显示之
d=zeros（32,32）;%图像大小32⨯32
d（13:
20,13:
20）=1;%中心白色方块大小为8⨯8
figure
（1）;%建立图形窗口1
imshow（d,'notruesize'）；%显示图像d如图3.5（a）所示
%计算傅里叶变换并显示之
D=fft2（d）;%计算图像d的傅里叶变换，fft2（d）=fft（fft（d）.'）.'
figure
（2）;%建立图形窗口2
imshow（abs（D）,[-15],'notruesize'）;%显示图像d的傅里叶变换谱如3.5（b）所示
例3.3MATLAB图像及其傅里叶变换谱
MATLAB程序：
figure
（1）;
loadimdemossaturn2;%装入MATLAB图像saturn2
imshow（saturn2）;%显示图像saturn2如图3.6（a）所示
figure
（2）;
S=fftshift（fft2（saturn2））;%计算傅里叶变换并移位
imshow（log（abs（S））,[]）;%显示傅里叶变换谱如3.6（b）所示
例3.4真彩图像及其傅里叶变换谱
MATLAB程序：
figure
（1）;
A=imread（'image1.jpg'）;%装入真彩图像，见图1.1（b）
B=rgb2gray（A）;%将真彩图像转换为灰度图像
imshow（B）;%显示灰度图像如图3.7（a）所示
C=fftshift（fft2（B））;%计算傅里叶变换并移位
figure
（2）;
imshow（log（abs（C））,[]）;%显示傅里叶变换谱如3.7（b）所示
3.5.4离散余弦变换的MATLAB实现
例3.5计算并显示真彩图像余弦变换的MATLAB程序如下：
RGB=imread（'image2.jpg'）;%装入真彩图像
figure
（1）;
imshow（RGB）;%显示彩色图像
GRAY=rgb2gray（RGB）;%将真彩图像转换为灰度图像
figure
（2）;
imshow（GRAY）;%显示灰度图像如图3.10（a）所示
DCT=dct2（GRAY）;%进行余弦变换
figure（3）;
imshow（log（abs（DCT））,[]）;%显示余弦变换如图3.10（b）所示。
3.8.2Radon变换的MATLAB实现
例3.8真彩图像的Radon变换
MATLAB程序如下：
RGB=imread（'image2.jpg'）;%装入真彩图像
GRAY=rgb2gray（RGB）;%将真彩图像转换为灰度图像
figure
（2）;
imshow（GRAY）;%显示灰度图像如图3.16（a）
[R,xp]=radon（GRAY,[045]）;%计算变换角度为0°和45°的Radon变换
figure;plot（xp,R（:
1））;title（'R_{0^o}（x\prime）'）
%显示0°方向上的Radon变换如图3.16（b）
figure;plot（xp,R（:
2））;title（'R_{45^o}（x\prime）'）
%显示45°方向上的Radon变换如图3.16（c）
例3.9连续角度的Radon变换
对于一组连续角度的Radon变换通常用一幅图像来表示。
本例先建立一幅简单图像，然后令变换角度从0°以1°的增量变化到180°时的Radon变换情况。
其MATLAB程序如下：
I=zeros（100,100）;%建立简单图像如图3.17（a）
I（25:
75,25:
75）=1;
figure
（1）;imshow（I）;
theta=0:
180;%规定变换角度的范围
[R,xp]=radon（I,theta）;%计算Radon变换
figure
（2）;
imagesc（theta,xp,R）;%以图像方式显示变换结果R，
%其x轴和y轴分别为theta和xp
title（‘R_{\theta}（X\prime）’）;%显示图像标题
xlabel（‘\theta（degrees）’）;%显示x坐标“
”
ylabel（‘X\prime’）;%显示y坐标“
”
set（gca,’Xtick’,0:
20:
180）;%设置x坐标刻度
colormap（hot）;%设置调色板
colorbar；%显示当前图像的调色板
第4章图像增强
4.2.1直接灰度变换
Matlab程序实现图像求反：
I=imread（'cameraman.tif'）;
imshow（I）
I=double（I）
I=256-1-I
I=uint8（I）
figure
imshow（I）
例4.1
用Matlab程序实现线性灰度变换的图像增强：
%读入并显示原始图像
I=imread（'pout.tif'）;
imshow（I）;
I=double（I）;
[M,N]=size（I）;
%进行线性灰度变换
fori=1:
M
forj=1:
N
ifI（i,j）<=30
I（i,j）=I（i,j）;
elseifI（i,j）<=150
I（i,j）=（200-30）/（150-30）*（I（i,j）-30）+30;
else
I（i,j）=（255-200）/（255-150）*（I（i,j）-150）+200;
end
end
end
例4.2
I=imread（'lena.bmp'）;
figure;imshow（I）;
I=double（I）;
I2=41*log（1+I）;
I2=uint8（I2）;
figure;imshow（I2）;
例4.3
灰度切割变换的Matlab的程序如下：
I=imread（'007.bmp'）;
figure;imshow（I）;
I=double（I）
[M,N]=size（I）;
fori=1:
M
forj=1:
N
ifI（i,j）<=50
I（i,j）=40;
elseifI（i,j）<=180
I（i,j）=220;
else
I（i,j）=40;
end
end
end
I=uint8（I）;
figure;imshow（I）;
例4.4
具体Matlab程序如下：
I=imread（'lena.bmp'）;
imshow（I）;
I=double（I）;
[M,N]=size（I）;
fork=1:
8
J=zeros（M,N）;
fori=1:
M
forj=1:
N
temp=I（i,j）;
s1=0;s2=0;
range=[k:
-1:
1];
ford=range
s1=2^（8-d）+s1;s2=2^（8-d+1）;
iftemp>=s1&tempJ（i,j）=255;break;
end
end
end
end
J=uint8（J）;
figure;imshow（J）;
end
4.2.2直方图修正
例4.6直方图均衡化效果实例
用Matlab中的histeq函数实现直方图均衡化的程序如下：
I=imread（'circuit.tif'）;
figure
subplot（221）;imshow（I）;
subplot（222）;imhist（I）
I1=histeq（I）;
figure;
subplot（221）;imshow（I1）
subplot（222）;imhist（I1）
例4.8：
直方图规定效果实例
用matlab中的histeq函数实现直方图均衡化的程序如下：
I=imread（'circuit.tif'）;
[M,N]=size（I）;
fori=1:
8:
257
counts（i）=i;
end
Q=imread（'circuit.tif'）;
N=histeq（Q,counts）;
figure
subplot（221）;imshow（N）;
subplot（222）;imhist（N）;
axis（[026005000]）;
4.2.3图像间运算
例：
用图像平均减少随机噪声
I=imread（'tire.tif'）;
[M,N]=size（I）;
II1=zeros（M,N）;
fori=1:
16
II（:
:
i）=imnoise（I,'gaussian',0,0.01）;
II1=II1+double（II（:
:
i））;
ifor（or（i==1,i==4）,or（i==8,i==16））;
figure;
imshow（uint8（II1/i））;
end
end
4.3空域滤波增强
Matlab实现的邻域平均法抑制噪声的程序：
I=imread（'eight.tif'）;
J=imnoise（I,'salt&pepper',0.02）;
subplot（231）,imshow（I）;title（'原图像'）;
subplot（232）,imshow（J）;title（'添加椒盐噪声图像'）
k1=filter2（fspecial（'average',3）,J）;%进行3×3模板平滑滤波
k2=filter2（fspecial（'average',5）,J）;%进行5×5模板平滑滤波
k3=filter2（fspecial（'average',7）,J）;%进行7×7模板平滑滤波
k4=filter2（fspecial（'average',9）,J）;%进行9×9模板平滑滤波
subplot（233）,imshow（uint8（k1））;title（'3×3模板平滑滤波'）;
subplot（234）,imshow（uint8（k2））;title（'5×5模板平滑滤波'）;
subplot（235）,imshow（uint8（k3））;title（'7×7模板平滑滤波'）;
subplot（236）,imshow（uint8（k4））;title（'9×9模板平滑滤波'）
例4.10：
使用中值滤波降低图像噪声
I=imread（'eight.tif'）;
J=imnoise（I,'salt&pepper',0.02）;
subplot（231）,imshow（I）;title（'原图像'）;
subplot（232）,imshow（J）;title（'添加椒盐噪声图像'）
k1=medfilt2（J）;%进行3×3模板中值滤波
k2=medfilt2（J,[55]）;%进行5×5模板中值滤波
k3=medfilt2（J,[77]）;%进行7×7模板中值滤波
k4=medfilt2（J,[99]）;%进行9×9模板中值滤波
subplot（233）,imshow（k1）;title（'3×3模板中值滤波'）
subplot（234）,imshow（k2）;title（'5×5模板中值滤波'）
subplot（235）,imshow（k3）;title（'7×7模板中值滤波'）
subplot（236）,imshow（k4）;title（'9×9模板中值滤波'）
例4.11：
梯度锐化实例
I=imread（'cameraman.tif'）;
subplot（131）,imshow（I）
H=fspecial（'Sobel'）;
H=H';%Sobel垂直模板
TH=filter2（H,I）;
subplot（132）,imshow（TH,[]）;
H=H';%Sobel水平模板
TH=filter2（H,I）;
subplot（133）,imshow（TH,[]）
4.4图像频域增强
例4.12：
频域低通滤波所产生的模糊
%理想低通过滤波器所产生的模糊和振铃现象
J=imread（'lena.bmp'）;
subplot（331）;imshow（J）;
J=double（J）;
%采用傅里叶变换
f=fft2（J）;
%数据矩阵平衡
g=fftshift（f）;
subplot（332）;imshow（log（abs（g））,[]）,color（jet（64））;
[M,N]=size（f）;
n1=floor（M/2）;
n2=floor（N/2）;
%d0=5,15,45,65
d0=5;
fori=1:
M
forj=1:
N
d=sqrt（（i-n1）^2+（j-n2）^2）;
ifd<=d0
h=1;
else
h=0;
end
g（i,j）=h*g（i,j）;
end
end
g=ifftshift（g）;
g=uint8（real（ifft2（g）））;
subplot（333）;
imshow（g）;
例4.12：
用巴特沃斯低通滤波器去除图像中的盐椒噪声
%实现Butterworth低通过滤波器
I=imread（'saturn.tif'）;
J=imnoise（I,'salt&pepper',0.02）;%给原图像加入椒盐噪声，如图4.33（a）所示
subplot（121）;imshow（J）;
tilte（'含有盐椒噪声的图像'）
J=double（J）;
%采用傅里叶变换
f=fft2（J）;
%数据矩阵平衡
g=fftshift（f）
[M,N]=size（f）;
n=3;
d0=20
n1=floor（M/2）
n2=floor（N/2）
fori=1:
M
forj=1:
N
d=sqrt（（i-n1）^2+（j-n2）^2）
h=1/（1+（d/d0）^（2*n））;
g（i,j）=h*g（i,j）;
end
end
g=ifftshift（g）;
g=uint8（real（ifft2（g）））;
subplot（122）;
imshow（g）;%结果如图4.33（b）所示
例4.13：
频域高通滤波增强示例
J=imread（'lenabu.bmp'）;
imshow（uint8（J））;title（'模糊图像'）
J=double（J）;
f=fft2（J）;%采用傅里叶变换
g=fftshift（f）;%数据矩阵平衡
[M,N]=size（f）;
n1=floor（M/2）;
n2=floor（N/2）;
d0=20;
fori=1:
M%进行理想高通滤波和理想高通加强滤波
forj=1:
N
d=sqrt（（i-n1）^2+（j-n2）^2）;
ifd>=d0
h1=1;
h2=1+0.5;
else
h1=0;
h2=0.5;
end
g1（i,j）=h1*g（i,j）;
g2（i,j）=h2*g（i,j）;
end
end
g1=ifftshift（g1）;
g1=uint8（real（ifft2（g1）））;
subplot（221）;imshow（g1）;%显示理想高通滤波结果
title（'理想高通滤波结果'）;
g2=ifftshift（g2）;
g2=uint8（real（ifft2（g2）））;
subplot（222）;imshow（g2）;%显示理想高通加强滤波结果
title（'理想高通加强滤波结果'）;
n=2;
d0=20;
fori=1:
M%进行巴特沃斯高通滤波和巴特沃斯高通加强滤波
forj=1:
N
d=sqrt（（i-n1）^2+（j-n2）^2）;
ifd==0
h1=0;
h2=0.5;
else
h1=1/（1+（d0/d）^（2*n））;
h2=1/（1+（d0/d）^（2*n））+0.5;
end
gg1（i,j）=h1*g（i,j）;
gg2（i,j）=h2*g（i,j）;
end
end
gg1=ifftshift（gg1）;
gg1=uint8（real（ifft2（gg1）））;
subplot（223）;imshow（gg1）;%显示巴特沃斯高通滤波结果
title（'巴特沃斯高通滤波结果'）
gg2=ifftshift（gg2）;
gg2=uint8（real（ifft2（gg2）））;
subplot（224）;imshow（gg2）;%显示巴特沃斯高通加强滤波结果
title（'巴特沃斯高通加强滤波结果'）;
例4.14：
同态滤波的增强效果
J=imread（'eight.tif'）;%读入原图
subplot（121）;imshow（J）;
J=double（J）;
f=fft2（J）;%采用傅里叶变换
g=fftshift（f）;%数据矩阵平衡
[M,N]=size（f）;
d0=10;
rl=0.5;
rh=2
c=4；
n1=floor（M/2）;
n2=floor（N/2）;
fori=1:
M
forj=1:
N
d=sqrt（（i-n1）^2+（j-n2）^2）;
h=（rh-rl）*（1-exp（-c*（d.^2/d0.^2）））+rl;
g（i,j）=h*g（i,j）;
end
end
g=ifftshift（g）;
g=uint8（real（ifft2（g）））;
subplot（122）;imshow（g）;
第5章图像复原
例5.1
C=imread（'image3.jpg'）;%装入清晰图像
subplot（1,2,1）；%将图形窗口分成两个矩形平面
imshow（C）;%在第一个矩形平面中显示装入的图像
LEN=30;%设置运动位移为30个象素
THETA=45;%设置运动角度为45o
PSF=fspecial（'motion',LEN,THETA）;%建立二维仿真线性运动滤波器PSF
MF=imfilter（C,PSF,'circular','conv'）;%用PSF产生退化图像
subplot（1,2,2）,imshow（MFUZZY）;%在第二个矩形平面中显示模糊后的图像
imwrite（MF,'image3-MF.jpg'）;%将运动模糊后的图像保存起来备用
例5.2消除图5.4（b）的运动模糊，其MATLAB程序如下：
[MF,map]=imread（'image3-MF.jpg'）;%装入运动模糊图像
figure
（1）;
imshow（MF）;%显示模糊图像
LEN=30;
THETA=45;
INITPSF=fspecial（'motion',LEN,THETA）;%建立复原点扩散函数
[JP]=deconvblind（MF,INITPSF,30）;%去卷积
figure
（2）;imshow（J）;%显示结果图像如图5.6（a）
figure（3）;imshow（P,[],'notruesize'）;%显示复原点扩散函数如图5.6（b）
例5.3逆滤波与维纳滤波的比较
F=checkerboard（8）;%生成原始图像F
figure
（1）;
imshow（F,[]）;
PSF=fspecial（'motion',7,45）;%生成运动模糊图像MF
MF=imfilter（F,PSF,'circular'）;
noise=imnoise（zeros（size（F））,'gaussian',0,0.001）;%生成高斯噪声
MFN=MF+noise;%生成运动模糊+高斯噪声图像MFN
figure
（2）;
imshow（MFN,[]）;
NSR=sum（noise（:
）.^2）/sum（MFN（:
）.^2）;%计算噪信比
figure（3）;
imshow（deconvwnr（MFN,PSF）,[]）;%逆滤波复原
figure（4）;
imshow（deconvwnr（MFN,PSF,NSR）,[]）;%维纳滤波复原
例5.7顺序统计滤波器比较
f=imread（'image4g.jpg'）;
figure
（1）;imshow（f）;
title（'原始图像'）;
g=imnoise（f,'salt&pepper',0.2）;
figure
（2）;imshow（g）;
title（'椒盐噪声污染的图像'）;
g1=double（g）/255;
j1=medfilt2（g1,'symmetric'）;
figure（3）;imshow（j1）;
title（'中值滤波图像'）;
j2=ordfilt2（g1,median（1:
3*3）,ones（3,3）,'symmetric'）;
figure（4）;imshow（j2）;
title（'中点滤波图像'）;
j3=ordfilt2（g1,1,ones（3,3））;
figure（5）;imshow（j3）;
title（'最小值滤波图像'）;
j4=ordfilt2（g1,9,ones（3,3））;
figure（6）;imshow（j4）;
title（'最大值滤波图像'）;
例5.8简单图像的affine变换
MATLAB程序如下：
f=checkerboard（24）;%建立原始图像，如图5.16（a）
figure
（1）;imshow（f）;
s=0.7;
theta=pi/6;
T=[s*cos（theta）s*sin（theta）0%建立变换矩阵：
旋转与尺度变换
-s*sin（theta）s*cos（theta）0
001];
tform=maketform（'affine',T）;
g1=imtransform（f,tform,'nearest'）;%最近邻插值变换,如图5.16（b）
figure
（2）;imshow（g1）;
g2=imtransform（f,tform）;%双线性插值变换,如图5.16（c）
figure（3）;imshow（g2）;
g3=imtransform（f,tfor
升级会员 