MATLAB数字图像处理实验.doc

MATLAB数字图像处理实验

———————————————————————————————— 作者：
———————————————————————————————— 日期：

13

个人收集整理 勿做商业用途
数字图像处理实验
图 像 处 理 实 验(一）直 方 图
灰度变换是图像增强的一种重要手段，使图像对比度扩展，图像更加清晰，特征更加明显。
灰度级的直方图给出了一幅图像概貌的描述，通过修改灰度直方图来得到图像增强。
1、 灰度直方图
（1） 计算出一幅灰度图像的直方图
clear
close all
I=imread('004.bmp’）；
imhist（I)
title（’实验一(1) 直方图');
(2) 对灰度图像进行简单的灰度线形变换,
figure
subplot（2,2，1）
imshow(I）;
title（’试验2—灰度线性变换’);
subplot(2，2，2）
histeq（I）；
(3） 看其直方图的对应变化和图像对比度的变化。
原图像 f(m,n) 的灰度范围 ［a，b] 线形变换为图像 g（m，n)，灰度范围[a',b’］
公式： g（m，n）=a'+（b’—a’)＊ f(m,n) /（b—a）
figure
subplot（2，2,1)
imshow(I)
J=imadjust（I,［0。3，0。7］,[0,1]，1）；
title（’ 实验一（3）用 g(m，n）=a’+(b’-a’)＊ f（m,n) /（b—a)进行变换 ');
subplot（2,2,2）
imshow(J）
subplot(2，2,3)
imshow（I)
J=imadjust(I，[0.5 0.8],[0,1］，1）；
subplot（2,2，4）
imshow（J）
（4) 图像二值化 (选取一个域值，(5） 将图像变为黑白图像）
figure
subplot(2，2，1)
imshow(I）
J=find（I<150);
I(J)=0;
J=find（I>=150);
I（J)=255；
title（' 实验一（4）图像二值化 （ 域值为150 ) ’）;
subplot(2，2，2）
imshow（I)
clc;
I=imread（’14499.jpg');
bw=im2bw（I，0。5)；%选取阈值为0.5
figure；
imshow(bw） ％显示二值图象
图象处理变换（二）
1． 傅立叶变换
熟悉其概念和原理,实现对一幅灰度图像的快速傅立叶变换，并求其变换后的系数分布。
2． 离散余弦变换
熟悉其概念和原理,实现对一幅灰度和彩色图像作的离散余弦变换,选择适当的DCT系数阈值对其进行DCT反变换。
％ 图象的FFT变换
clc；
I=imread（’005.bmp’）；
subplot（1,2，1)
imshow（I）;
title（’原图’）；
subplot(1,2,2）
imhist(I)；
title('直方图')；
colorbar
J=fft2(I);
figure；
subplot(1,2，1)
imshow（J）；
title(’FFT变换结果')；
subplot(1,2，2)
K=fftshift（J);
imshow(K）；
title(’零点平移’）;
figure；
imshow(log(abs（K))，［]),colormap(jet(64）),colorbar;
title('系数分布图'）;
% 图象的DCT变换
RGB=imread（’005。bmp’)；
figure;
subplot(1,2,1）
imshow(RGB）;
title(’彩色原图')；
a=rgb2gray(RGB);
subplot（1,2，2)
imshow（a)；
title（’灰度图');
figure;
b=dct2(a);
imshow(log(abs(b）），［])，colormap(jet(64）)，colorbar;
title(’DCT变换结果'）；
figure；
b(abs（b）〈10）=0；
% idct
c=idct2（b)/255;
imshow（c);
title（'IDCT变换结果’）；
图象处理变换(三）小波变换
实验内容： 熟悉小波变换的概念和原理，熟悉matlab小波工具箱主要函数的使用。利用二维小波分析对一幅图象作2层小波分解，并在此基础上提取各层的低频信息实现图像的压缩。
程序如下：
clc
close all
clear
a=imread（’005。bmp’）;
subplot（1，2，1）；
imshow(a);
title（’原始图象')；
I=rgb2gray（a);
subplot(1,2,2);
imshow(I）；
title('原始图象的灰度图’）;
％ 进行二维小波变换
[a，b] = wavedec2(I, 2, ’bior3.7')；
% 提取各层低频信息
figure;
c = appcoef2( a， b， ’bior3。7', 1 )；
subplot（1，2,1)；
imshow（c， ［])；
title(’一层小波变换结果'）;
d = appcoef2（ a, b， ’bior3.7’, 2 );
subplot(1，2,2）；
imshow（d, ［]）;
title(’二层小波变换结果’)；
图象处理实验（四） 模板运算
一、实验内容：
（1）平滑：平滑的目的是模糊和消除噪声。平滑是用低通滤波器来完成，在空域中全是正值。
（2)锐化:锐化的目的是增强被模糊的细节。锐化是用高通滤波器来完成，在空域中，接近原点处为正,在远离原点处为负。
利用模板进行图象增强就是进行模板卷积。
1、 利用二个低通邻域平均模板(3×3和9×9）对一幅图象进行平滑，验证模板尺寸对图象的模糊效果的影响。
2、 利用一个低通模板对一幅有噪图象（GAUSS白噪声）进行滤波，检验两种滤波模板（分别使用一个5×5的线性邻域平均模板和一个非线性模板:3×5中值滤波器)对噪声的滤波效果。
3、 选择一个经过低通滤波器滤波的模糊图象，利用sobel和prewitt水平边缘增强高通滤波器（模板)对其进行高通滤波图象边缘增强，验证模板的滤波效果。
4、 选择一幅灰度图象分别利用 一阶Sobel算子和二阶Laplacian算子对其进行边缘检测，验证检测效果。
二、实验步骤：
1、利用低通邻域平均模板进行平滑
I=imread('girl.bmp’）；
subplot（1，3，1)；
imshow(I);
title('原图’)；
J=fspecial('average'）；
J1=filter2(J,I)/255;
subplot（1，3，2）;
imshow（J1)；
title（’3＊3滤波')；
K=fspecial（'average’，9);
K1=filter2（K,I）/255；
subplot(1，3,3）；
imshow(K1）；
title（'9＊9滤波')；
2、中值滤波和平均滤波
I=imread(’girl。bmp')；
J=imnoise(I,'gaussian’,0,0.01);
subplot（2,2,1)；
imshow（I);
title('原图');
subplot(2,2,2);
imshow（J）；
title(’noise')；
K=fspecial（’average'，5);
K1=filter2（K,J）/255;
subplot(2，2，3);
imshow(K1）;
title（'average’）；
L=medfilt2（J,[3 5］);
subplot(2,2,4);
imshow（L）；
title('medium’）；
3、高通滤波边缘增强
I=imread（’girl。bmp'）；
subplot（2，2，1);
imshow(I）;
title(’original pic')；
J=fspecial（'average’，3）;
J1=conv2（I，J）/255；
％J1=filter2（J,I)/255；
subplot（2，2，2）;
imshow(J1）;
title（’3*3lowpass’）;
K=fspecial（’prewitt’）;
K1=filter2（K，J1)*5;
subplot（2,2，3)；
imshow（K1);
title(’prewitt');
L=fspecial(’sobel'）;
L1=filter2（L,J1）*5;
subplot（2，2，4）；
imshow(L1）；
title('sibel'）；
4、边缘检测
分别用sobel和laplacian算子来进行，程序如下：
I=imread(’girl。bmp'）;
subplot(1,3,1）；
imshow（I);
title（’original pic’)；
K=fspecial(’laplacian',0.7）；
K1=filter2（K，I）/100；
subplot(1，3,2）；
imshow(K1);
title('laplacian');
L=fspecial（’sobel')；
L1=filter2（L,I)/200；
subplot（1,3，3)；
imshow(L1）；
title（'sibel'）；
图像处理实验（五)图像分割
实验目的：1、 学习边缘检测
2、 学习灰度阀值分割
实验内容:
1、分别用sobel、Laplacian—Gaussian方法对一幅灰度图像进行边缘提取，2、给出对比结果
i=imread(’eight.tif’)；
figure;
subplot(2，2,1)；
imshow(i）;
title(’原始图像'）；
subplot（2，2,3)；
imshow(i)；
title（'原始图像’);
i1=edge(i,'sobel')；
subplot(2,2，2);
imshow(i1)；
title（'sober方法提取的边缘’)；
i2=edge（i,’log')；
subplot(2，2,4）;
imshow(i2）；
title（’Laplacian—Gaussian方法提取的边缘'）；
比较提取边缘的效果可以看出，sober算子是一种微分算子，对边缘的定位较精确，但是会漏去一些边缘细节。而Laplacian—Gaussian算子是一种二阶边缘检测方法,它通过寻找图象灰度值中二阶过零点来检测边缘并将边缘提取出来，边缘的细节比较丰富。通过比较可以看出Laplacian-Gaussian算子比sober算子边缘更完整，效果更好。
3、利用双峰法对一幅灰度图像进行灰度分割处理
i=imread(’eight。tif')；
subplot(1,2,1);
imhist（i）；
title（'原始图像直方图’）;
thread=130/255；
subplot（1,2,2）；
i3=im2bw(i，thread）;
imshow（i3）;
title('分割结果’)；
根据原图像的直方图，发现背景和目标的分割值大约在130左右，并将灰度图像转为二值图像,分割效果比较理想.
图像处理实验（六）图像压缩与编码
实验目的： 学习JPEG压缩编码
实验内容：
一．实现基本JPEG的压缩和编码分三个步骤:
1． 首先通过DCT变换去除数据冗余；
2． 使用量化表对DCT系数进行量化；
3． 对量化后的系数进行Huffman编码.
具体源程序由主程序及两个子程序（DCT量化、Huffman编码）组成：
1．主程序
I=imread('autumn.tif'）;
yiq=rgb2ntsc（I)；
my=[16 11 10 16 24 40 51 61；12 12 14 19 26 58 60 55；
14 13 16 24 40 57 69 56；14 17 22 29 51 87 80 62；
18 22 37 56 68 109 103 77；24 35 55 64 81 104 113 92；
49 64 78 87 103 121 120 101;72 92 95 98 112 100 103 99］；
miq=［17 18 24 47 99 99 99 99；18 21 26 66 99 99 99 99;
24 26 56 99 99 99 99 99；47 66 99 99 99 99 99 99；
99 99 99 99 99 99 99 99；99 99 99 99 99 99 99 99;
99 99 99 99 99 99 99 99;99 99 99 99 99 99 99 99]；
I1=yiq（：,:,1);I2=(：，：,2)；
［m n］=size（I1）；
t1=8;ti1=1；
while(t1
t1=t1+8;ti1=ti1+1;
end
t2=8；ti2=1；
while(t2
t2=t2+8;ti2=ti2+1;
end
times=0;
for k=0:ti1—2
for j=0:ti2-2
dct8x8（I1（k*8+1:k*8+8,j＊8+1:j＊8+8），my,times＊64+1）;
dct8x8(I2(k*8+1：k＊8+8,j*8+1:j*8+8),miq，times*64+1)；
times=times+1;
end
block(I2（k*8+1:k*8+8,j*8+1：t2），[8 8], ’dctmtx(8）'）；
end
for j=0:ti2—2
dct8x8（I1（k*8+1:t1，j＊8+1：j＊8+8)，times＊64+1）；
times=times+1;
end
dct8x8(I1(k＊8+1：t1,j*8+1：t2)，times*64+1）;
2．function dct8x8（I，m，s) %定义DCT量化子程序
T=inline(’dctmtx（8）')；
y=blkproc(I,［8 8],T)；
y=round(y。/m)；
p=1；te=1；
while(p<=64)
for q=1：te
y1（s+p)=y（te-q+1,q);p=p+1；
end
for q=te:-1：1
y1（s+p）=y（te—q+1,q);p=p+1；
end
end
f=haffman（y1）;
c（s：s+64，1）=f（:，1）;c（s：s+64,2）=f（:,2)；c（s:s+64，3）=f(：,3）
3．function c=haffman(I） %定义Huffman编码子程序
[m,n]=size(I）;
p1=1;s=m＊n;
for k=1:m
for h=1：n

f=0;
for b=1:p1-1
if（c(b,1）==I(k，h)) f=1;break；end
end
if（f==0) c（p1,1）=I(k，h）；p1=p1+1;end
end
end
for g=1:p1-1
p(g）=0；c（g,2）=0；
for k=1:m
for h=1:n
if(c（g，1）==I（k,h)) p（g）=p（g）+1；end
end
end
p（g）=p（g）/s;
end
pn=0；po=1；
while(1）
if(pn〉=1.0） break;
else
[pm p2］=min(p(1:p1-1））；p(p2)=1.1;
［pm2，p3]=min(p(1:p1-1))；p(p3)=1.1；
pn=pm+pm2；p（p1）=pn;
tree(po,1）=p2;tree(po，2）=p3；
po=po+1;p1=p1+1；
end
end
for k=1：po—1
tt=k；m1=1；
if（or（tree（k,1）<9，tree（k,2)<9))
if(tree(k,1)〈9)
c（tree(k，1),2）=c（tree（k,1）,2）+m1;
m2=1;
while(tt
m1=m1*2；
for h=tt:po—1
if(tree（h，1)==tt+g）
c(tree（k，1）,2)=c(tree(k,1)，2)+m1;
m2=m2+1;tt=h;break；
elseif（tree（h,2)==tt+g)
m2=m2+1；tt=h;break；
end
end
end
c（tree(k，1）,3）=m2;
end
tt=k;m1=1;
if（tree(k，2）〈9）
m2=1；
while(tt
m1=m1＊2；
for h=tt:po—1
if（tree（h,1）==tt+g）
c（tree（k，2)，2)=c（tree（k,2)，2)+m1；
m2=m2+1；tt=h；break;
elseif(tree(l，2）==tt+g)
m2=m2+1,tt=h；break；
end
end
end
c（tree（k,2),3）=m2;
end
end
end
二．JPEG2000采用小波变换编码，
小波变换压缩编码实现程序为
load wbarb；
subplot（2,2,1),image(X)；colormap（map)
title(’原始图象')；
［c,s］=wavedec2（X,2, ’bior3.7'）；
thr=20；
ca1=appcoed2(c，s， 'bior3.7',1)；
ch1=detcoef2（’h’，c，s，1);
cv1=detcoef2('v’,c，s，1）;
cd1=detcoef2(’d’，c,s,1）；
a1=wrcoef2（'a’,c,s， 'bior3。7’,1）；
h1=wrcoef2(’h'，c，s， ’bior3。7'，1)；
v1=wrcoef2（'v',c,s， ’bior3。7'，1);
d1=wrcoef2('d’,c，s, ’bior3.7’，1);
c1=[a1,h1,v1，d1］;
ca1=appcoed2（c，s， 'bior3。7',1）；
ca1=wcodemat(ca1,440, ’mat',0）;
ca1=0.5＊ca1
subplot(2,2,2),image(ca1）
title（’压缩图象一')
ca2=appcoed2（c，s, 'bior3。7’，2）；
ca2=wcodemat(ca2，440, 'mat’，0)；
ca2=0.5*ca2;
subplot(2,2,3),image(ca2）
title('压缩图象二’）
图象处理实验(七）应用KL变换进行图象的特征提取
一、实验要求：应用KL变换进行图象的特征提取。熟悉MATLAB的相关命令。
二、实验目的：掌握如何应用KL变换进行图象的特征提取.
三、实验内容：选择一幅稍大的灰度图象（最好用纹理图象)，按下面步骤进行实验：
（1）应用9×9的窗口对上述图象进行随机抽样，共抽样200块子图象;
（2）将所有子图象按列相接变成一个81维的行向量；
（3）对所有200个行向量进行KL变换，求出其对应的协方差矩阵的特征向量和特征值，按降序排列特征值以及所对应的特征向量；
（4)选择前40个最大特征值所对应的特征向量作为主元，将原图象块向这40个特征向量上投影，所获得的投影系数就是这个子块的特征向量。
（5)求出所有子块的特征向量。
四、实验结果：
源程序如下：
clear
close all
clc
M=rand（200,200）；
I=imread（’a。bmp'）;
［mx，my]=size（I）;
imshow(I）；
title（'原始图象’）;
for i=1:200
for j=1:199
if(ceil（M（i，j）*mx)
x(i)=ceil(M（i,j)*mx）;
y(i)=ceil（M（i，j+1)＊my);
end
end
end
for i=1：200
I1(：,：,i）=imcrop(I,[x（i）,y(i），8,8］）;
I2（:,i）=reshape（I1（:,:,i),1，81);
end
I2=double(I2）；
C=I2'＊I2/200;
[a,s]=eig(C）；
A=a(161：200，：）;
U=A*I2'；
figure;
imshow（U);
title('特征向量')；