MatLab代码大全.doc

第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)
shading interp; %平滑三维曲面上的小格
axis off; %关闭坐标系
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),[-1 5],'notruesize'); %显示图像d的傅里叶变换谱如3.5(b)所示
例3.3 MATLAB图像及其傅里叶变换谱
MATLAB程序：
figure(1);
load imdemos saturn2; %装入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.2 Radon变换的MATLAB实现
例3.8真彩图像的Radon变换
MATLAB程序如下：
RGB=imread('image2.jpg'); %装入真彩图像
GRAY=rgb2gray(RGB); %将真彩图像转换为灰度图像
figure(2);
imshow(GRAY); %显示灰度图像如图3.16(a)
[R,xp] = radon(GRAY,[0 45]); %计算变换角度为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);
%进行线性灰度变换
for i=1:M
for j=1:N
if I(i,j)<=30
I(i,j)=I(i,j);
elseif I(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);
for i=1:M
for j=1:N
if I(i,j)<=50
I(i,j)=40;
elseif I(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);
for k=1:8
J=zeros(M,N);
for i=1:M
for j=1:N
temp=I(i,j);
s1=0;s2=0;
range=[k:-1:1];
for d=range
s1=2^(8-d)+s1;s2=2^(8-d+1);
if temp>=s1 & temp<s2;
J(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);
for i=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([0 260 0 5000]);

4.2.3 图像间运算

例：用图像平均减少随机噪声
I=imread('tire.tif');
[M,N]=size(I);
II1=zeros(M,N);
for i=1:16
II(:,:,i)=imnoise(I,'gaussian',0,0.01);
II1=II1+double(II(:,:,i));
if or(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,[5 5]); %进行5×5模板中值滤波
k3=medfilt2(J,[7 7]); %进行7×7模板中值滤波
k4=medfilt2(J,[9 9]); %进行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;
for i=1:M
for j=1:N
d=sqrt((i-n1)^2+(j-n2)^2);
if d<=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)
for i=1:M
for j=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;
for i=1:M %进行理想高通滤波和理想高通加强滤波
for j=1:N
d=sqrt((i-n1)^2+(j-n2)^2);
if d>=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;
for i=1:M %进行巴特沃斯高通滤波和巴特沃斯高通加强滤波
for j=1:N
d=sqrt((i-n1)^2+(j-n2)^2);
if d==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);
for i=1:M
for j=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); %建立复原点扩散函数
[J P]= 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
0 0 1];
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,tform,'FillValue',0.5); % 修改双线性插值变换的背景色为灰色
figure(4); imshow(g3);

例5.9 利用“连接点”实施图像配准复原
MATLAB程序如下：
f=imread('text.jpg'); % 读入256´256原始图，如图5.17(a)
figure(1); imshow(f);
g=imread('textg.jpg'); % 读入几何失真图，，如图5.17(b)
figure(2); imshow(g);
%利用cpselect(g, f)函数交互选择“连接点”
base_points =[ 256.4000 256.1273; 1.5818 256.4182;
256.4182 1.0000; 200.5636 203.4727;
147.3273 183.9818; 96.4182 145.0000;
44.6364 35.0364; 157.5091 30.3818];
input_points =[ 280.0455 304.6182; 1.3455 255.2545;
255.2545 1.0000; 205.8545 225.8909;
145.7455 196.5273; 90.4727 146.7818;
38.6545 32.4364; 148.5091 31.0545];
tform = cp2tform(input_points, base_points, 'projective');
gp = imtransform(g, tform, 'XData', [1 256], 'YData', [1 256]);
figure(3); imshow(gp); % 显示复原图像，如图5.17(c)

第6章 彩色图像处理

例6.1 考虑生成一幅128*128的RGB图像，该图像左上角为红色，右上角为蓝色，左下角为绿色，右下角为黑色。其Matlab程序如下：
iR = zeros(128, 128); %生成一个128*128的零矩阵，作为R分量
iR(1:64, 1:64) = 1; %将左上角的64*64设置成1
iG = zeros(128, 128); %生成一个128*128的零矩阵，作为G分量
iG(65:128, 1:64) = 1; %将右下角的64*64设置成1
iB = zeros(128, 128); %生成一个128*128的零矩阵，作为B分量
iB(1:64, 65:128) = 1; %将右上角的64*64设置成1
I = cat(3, iR, iG, iB); %使用cat函数将三个分量组合
imshow(I) %显示生成的RGB图像，如图6.5所示

例6.2 将一幅RGB图像转换到CMY空间：
rgb_I = imread('peppers.bmp'); %载入一幅彩色图像
cmy_I = imcomplement(rgb_I); %函数imcomplement转换到CMY空间
imshow(I); %显示原图，如图6.7（a）所示
figure, imshow(I2); %显示转换后图，如图6.7（b）所示

例6.3 将一幅三原色图像从RGB空间转换到HSI空间，其结果见图6.9。
rgb = imread('三原色.bmp'); %载入一幅图像
imshow(rgb); %见图6.9（a）
rgb = im2double(rgb); %将图像转换成double类型
r = rgb(:, :, 1); %提取图像的r分量
g = rgb(:, :, 2); %提取图像的g分量
b = rgb(:, :, 3); %提取图像的b分量

I = (r + g + b)/3; %计算I分量
%计算S分量
tmp1 = min(min(r, g), b);
tmp2 = r + g + b;
tmp2(tmp2 == 0) = eps;
S = 1 - 3.*tmp1./tmp2;
%计算H分量
tmp1 = 0.5*((r - g) + (r - b));
tmp2 = sqrt((r - g).^2 + (r - b).*(g-b));
theta = acos(tmp1./(tmp2 + eps));
H = theta;
H(b > g) = 2*pi - H(b > g);
H = H/(2*pi);
H(S == 0) = 0;

hsi = cat(3, H, S, I);
figure, imshow(H); %见图6.9（b）
figure, imshow(S); %见图6.9（c）
figure, imshow(I); %见图6.9（d）

例6.4
rgb = imread(‘yellowRose.jpg’); %载入一幅图像
imshow(rgb); %显示，见图6.10（a）
R = rgb(:, :, 1); %提取图像的R、G、B分量
G = rgb(:, :, 2);
B = rgb(:, :, 3);
figure, imshow(R); %分别显示图像的R、G、B分量。见图6.10(b)。
figure, imshow(G); %见图6.10(c)
figure, imshow(B); %见图6.10(d)
m = fspecial(‘average’); %生成一个空间均值滤波器
R_filtered = imfilter(R, m); %分别对图像的R、G、B分量进行滤波。
G_filtered = imfilter(G, m);
B_filtered = imfilter(B, m);
rgb_filtered = cat(3, R_filtered, G_filtered, B_filtered);
figure, imshow(rgb_filtered); %见图6.10（e）

例6.5 彩色图像锐化。
I = imread(‘yellowRose.jpg’); %见图6.10（a）
imshow(I);
lapMatrix = [1 1 1; 1 –8 1; 1 1 1]; % Laplacian模板
I_tmp = imfilter(fb, lapMatrix, ‘replicate’); %滤波
I_sharped = imsubtract(I, I_