数字图像处理实验指导书(2015.11.04)

1、数字图像处理 实验指导书（自编）实验1 空域增强-直接灰度映射实验2 空域增强平滑滤波器（邻域平均 vs 中值滤波）实验3 图像变换-傅里叶变换及离散余弦变换原理及应用实验4 频域增强理想低通滤器vs Butterworth低通滤波实验5 图像恢复逆滤波 vs 维纳滤波实验6 图像恢复-交互式恢复去除周期性噪声实验7 图像编码-算术编码方法实验8 图像编码-维游程编码方法实验1 空域增强-直接灰度映射授课教师肖红授课班级信息12 学时理论48学时+实验16学时实验环境Matlab7.1实验内容：1熟悉matlab工作界面2直接灰度映射的典型方法图像求反 对比度拉伸对数变换 灰度切分3.注意观察

2、几种方法对图像进行处理后的对比效果注：标签图像文件格式（Tagged Image File Format，简写为TIFF）是一种灵活的位图格式，主要用来存储包括照片和艺术图在内的图像。它最初由Aldus公司与微软公司一起为PostScript打印开发。TIFF与JPEG和PNG一起成为流行的高位彩色图像格式。1直接灰度映射的典型方法图像求反clear;close all;X1=imread('cameraman.tif');figure,imshow(X1);f0=0;g0=255;f1=255;g1=0figure,plot(f0,f1,g0,g1)axis tight,xl

3、abel('f'),ylabel('g')title('intensitytransformation')%绘制变换曲线r1=(g1-g0)/(f1-f0);%等于-1b1=g0-r1*f0;%等于255m,n=size(X1);X2=double(X1);for i=1:mfor j=1:nf=X2(i,j);g(i,j)=0;g(i,j)=r1*f+b1; %这其实就是用255去减f，255-fendendfigure, imshow(mat2gray(g)figure,imshow(imcomplement(X1),title('利

4、用imcomplement对图像求反的效果')（1）原始图像（2）映射曲线（3）求反后效果图1-1 图像求反示意图对比度拉伸clear;close all;X1=imread('cameraman.tif');figure,imshow(X1)f0=0;g0=0;f1=60;g1=50;f2=150;g2=240;f3=255;g3=255;figure,plot(f0,f1,f2,f3,g0,g1,g2,g3)axis tight,xlabel('f'),ylabel('g')title('intensitytransforma

5、tion')%绘制变换曲线r1=(g1-g0)/(f1-f0);%斜率b1=g0-r1*f0;%截距r2=(g2-g1)/(f2-f1);b2=g1-r2*f1;r3=(g3-g2)/(f3-f2);b3=g2-r3*f2;m,n=size(X1);X2=double(X1);for i=1:mfor j=1:nf=X2(i,j);g(i,j)=0;if (f>=f0)&(f<f1) g(i,j)=r1*f+b1;elseif(f>=f1)&(f<=f2) g(i,j)=r2*f+b2;elseif(f>=f2)&(f<=f3

6、) g(i,j)=r3*f+b3;endendendfigure,imshow(mat2gray(g)（1）原始图像（2）映射曲线（3）增强后效果图1-2 增强对比度示意图图1-3 imadjdemo调整对比度对数变换clearclose allN=100f=zeros(50,50);f(15:35,23:28)=1; %这样是一个长方形，很适合观察傅里叶变换的特征F=fft2(f,N,N);F2=fftshift(abs(F);%再以灰度图的形式显示一下figure;subplot(121);imshow(F2,);%直接显示的话，图像很黑，看不清楚subplot(122);imshow(0

7、.25*log10(1+F2),);%这条语句用于压缩傅里叶频谱的动态范围。 %因为高频区域衰减的很快，导致频谱中低亮度值的可视细节丢失。通过计算对数，频谱的动态范围会大大减小，便于观察和处理。 %注：此处求的是log10，即以10为底的(1+F2)的对数。图1-4 傅立叶频谱动态范围压缩灰度切分 clear;close all;X1=imread('cameraman.tif');figure,imshow(X1);f0=0;g0=30;f1=50;g1=30;f2=50;g2=200;f3=180;g3=200;f4=180;g4=30;f5=255;g5=30;figur

8、e,plot(f0,f1,f1,f2,f3,f4,f5,g0,g1,g1,g2,g3,g4,g5)axis(0 255 0 255),xlabel('f'),ylabel('g')title('intensitytransformation')%绘制变换曲线m,n=size(X1);X2=double(X1);for i=1:mfor j=1:nf=X2(i,j);if (f>=f2) & (f<=f3) g(i,j)=200;else g(i,j)=30;endendendfigure, imshow(mat2gray(g)

9、图1-5 灰度切分函数示意图实验2 空域增强平滑滤波器（邻域平均 vs 中值滤波）授课教师肖红授课班级信息12 学时理论48学时+实验16学时实验环境Matlab7.1实验内容：1熟悉matlab工作界面2理解空域增强技术的平滑滤波技术邻域平均 中值滤波器3.观察并对比两种平滑滤波器的滤波效果1邻域平均% 邻域平均的线性平滑滤波法实现降噪的源程序如下：clear all;I=imread('cameraman.tif');%读入预处理图像imshow(I)%显示预处理图像%向原始图像中添加噪声J= imnoise(I,'salt & pepper',0.

10、02);figure,imshow(J);title('salt & pepper noise');K1=filter2(fspecial('average',3),J);%此处如果不除以255，则下面imshow显示的时候应该用imshow(K1,)，也可以。%进行3*3均值滤波K2=filter2(fspecial('average',5),J);%进行5*5均值滤波K3=filter2(fspecial('average',7),J);%进行7*7均值滤波figure,imshow(K1,)figure,imshow

11、(K2,)figure,imshow(K3,) 图2-1邻域平均效果2中值滤波% 调用medfilt2(A,m,n)来实现中值滤波的源程序：clear all;close all;I=imread('cameraman.tif');%原图像读入subplot(1,2,1),imshow(I),title('原始图像');%原始图像显示J= imnoise(I,'salt & pepper',0.08);%向原始图像中添加噪声subplot(1,2,2),imshow(J);title('含噪声图像');K1=medfilt

12、2(J,3,3);%对噪声图像进行滤波处理figure;subplot(1,2,1),imshow(K1),title('3*3中值滤波结果');%使用3*3模板完成中值滤波K2=medfilt2(J,9,9);subplot(1,2,2),imshow(K2),title('9*9中值滤波结果');图2-2中值滤波效果3思考不使用MATLAB提供的函数filter2和medfilt2，请同学们试着自己编写邻域平均和中值滤波的程序？实验3 图像变换-傅里叶变换及离散余弦变换原理及应用授课教师肖红授课班级信息12学时理论48学时+实验16学时实验环境Matlab7

13、.1实验内容：1熟悉matlab工作界面2理解傅里叶变换及离散余弦变换傅里叶变换 离散余弦变换3.观察傅里叶变换频谱的特征4.理解离散余弦变换在图像压缩方面的应用1熟悉matlab工作界面2理解傅里叶变换及离散余弦变换关于傅里叶变换clearclose allN=100f=zeros(50,50);f(15:35,23:28)=1; %这样是一个长方形，很适合观察傅里叶变换的特征%f(15:35,15:35)=1;%这样就是一个正方形，则和教材45页的图3.2.1基本完全吻合，但其实正方形不利于观察傅里叶变换figure;subplot(221);imshow(f,'notruesiz

14、e'),title('原始图像')F=fft2(f,N,N);F2=fftshift(abs(F);%再以灰度图的形式显示一下subplot(222);imshow(F,),title('原始图像的傅里叶频谱')subplot(223);imshow(F2,),title('经过平移的傅里叶频谱')%直接显示的话，图像很黑，看不清楚subplot(224);imshow(0.25*log10(1+F2),),title('平移后经过动态范围压缩的傅里叶频谱');%这条语句用于压缩傅里叶频谱的动态范围。 %因为高频区域衰减的

15、很快，导致频谱中低亮度值的可视细节丢失。通过计算对数，频谱的动态范围会大大减小，便于观察和处理。 %注：此处求的是log10，即以10为底的(1+F2)的对数。 %反变换g=ifft2(F);figure;imshow(real(g(1:50,1:50),'notruesize'),title('反变换后得到的重构图像');%得到扩充之前的原始图像x=1:N;y=1:N;figure;mesh(x,y,F2(x,y);%三维网格曲面colormap(gray);colorbar关于离散余弦变换例1：clear;close all;RGB = imread(

16、9;autumn.tif');I = rgb2gray(RGB);J = dct2(I);imshow(log(abs(J),);,title('DCT结果') %对变换结果取绝对值，然后取自然对数J(abs(J) < 10) = 0;K = idct2(J);figure,imshow(I),title('原始图像')figure,imshow(K,0 255),title('重构图像')例2：I = imread('cameraman.tif');I = im2double(I);T = dctmtx(8);%产

17、生DCT变换矩阵B = blkproc(I,8 8,'P1*x*P2',T,T');%计算二维DCTmask = 1 1 1 1 0 0 0 0 1 1 1 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;%二值掩模，用来压缩DCT的系数B2 = blkproc(B,8 8,'P1.*x',mask);%只保留DCT的10个系数I2 = blkproc(B2,8 8,'P1*x*P2&

18、#39;,T',T);%逆DCT，用来重构图像imshow(I), figure, imshow(I2) 实验4 频域增强理想低通滤器 vs Butterworth低通滤波授课教师肖红授课班级信息12学时理论48学时+实验16学时实验环境Matlab7.1实验内容：1理解频域增强的原理2理解理想低通滤波和Butterworth低通滤波的原理 理想低通滤波 Butterworth低通滤波3.观察并对比两种低通滤波器的滤波效果1理想低通滤波% 5-3-1 理想低通滤波器实现例程如下：clear all;X=imread('D:肖红study matlabimagesourcetif

19、lena.tif');% I1 = grayslice(X,16);%由256灰度级量化为16个灰度级，目的是为了产生虚假轮廓I1=X;figure,subplot(121),imshow(I1,),title('原始图像');f=double(I1);g=fft2(f);%傅立叶变换g=fftshift(g);%转换数据矩阵N1,N2=size(g);n=1;d0=30; %可以通过修改d0来决定过滤掉多少高频，从而查看不同的理想低通的效果 n1=fix(N1/2);n2=fix(N2/2);for i=1:N1for j=1:N2d=sqrt(i-n1)2+(j-n

20、2)2);if d>d0 %距离频域中心距离超过d0的都被完全过滤，体现理想低通 h=0; else h=1; endresult(i,j)=h*g(i,j); %这里是相乘，不是卷积endendresult=ifftshift(result);X2=ifft2(result);% X3=uint8(real(X2);subplot(122),imshow(real(X2),),title('截断频率d0=30的理想低通滤波结果');%存在振铃效应2Butterworth低通滤波% 5-3-1 Buterworth低通滤波器实现例程如下：clear all;%先读入一幅图

21、像，然后将其由256灰度级量化为16个灰度级，或者更少的灰度级，造成虚假轮廓现象。%然后再利用低通滤波（理想低通滤波或Buterworth低通滤波），来消除虚假轮廓现象。查看其对比效果。X=imread('D:肖红study matlabimagesourcetiflena.tif'); I1 = grayslice(X,16);%由256灰度级量化为16个灰度级figure,subplot(121),imshow(I1,),title('灰度级为16，存在虚假轮廓效应的原始图像');f=double(I1);g=fft2(f);%傅立叶变换g=fftshift

22、(g);%转换数据矩阵N1,N2=size(g);n=1;d0=30;n1=fix(N1/2);n2=fix(N2/2);for i=1:N1for j=1:N2d=sqrt(i-n1)2+(j-n2)2);% h=1/(1+0.414*(d/d0)(2*n); h=1/(1+(d/d0)(2*n);%这里为了和教材上完全一致，所以没有用0.414%计算Butterworth低通转换函数result(i,j)=h*g(i,j);endendresult=ifftshift(result);X2=ifft2(result);% X3=uint8(real(X2);%这里不能再用uint8格式了，

23、那样的话等同于又要量化成256个灰度级了，重新出现了虚假轮廓效应。subplot(122),imshow(real(X2),0 15),title('阶为1，截断频率d0=30的butterworth低通滤波结果')%正确的显示滤波结果，这一步非常非常关键实验5 图像恢复逆滤波 vs 维纳滤波授课教师肖红授课班级信息12学时理论48学时+实验16学时实验环境Matlab7.1实验内容：1理解无约束恢复和有约束恢复的区别2理解逆滤波和维纳滤波恢复图像的原理 逆滤波 维纳滤波3.比较逆滤波恢复和维纳滤波恢复的效果1程序清单clear;I=imread(' E:1-肖红老婆盘

24、study matlabimagesourcetifrice.tif');figure;subplot(2,2,1);imshow(I);title('原图像');m,n=size(I);F=fftshift(fft2(I);k=0.0025;for u=1:m for v=1:n H(u,v)=exp(-k)*(u-m/2)2+(v-n/2)2)(5/6); endendG=F.*H;I0=real(ifft2(fftshift(G);I1=imnoise(uint8(I0),'gaussian',0,0.001);subplot(2,2,2);ims

25、how(uint8(I1);title('模糊退化且添加高斯噪声的图像');F0=fftshift(fft2(I1);F1=F0./H;I2=ifft2(fftshift(F1);subplot(2,2,3);imshow(uint8(I2);title('全逆滤波复原图');K=0.1; for u=1:m for v=1:n H(u,v)=exp(-k*(u-m/2)2+(v-n/2)2)(5/6); H0(u,v)=(abs(H(u,v)2; H1(u,v)=H0(u,v)/(H(u,v)*(H0(u,v)+K); endendF2=H1.*F0;I3=i

26、fft2(fftshift(F2);subplot(2,2,4);imshow(uint8(I3);title('维纳滤波复原图');2结果分析 如运行结果所示，其中退化函数中常数k=0.0025（剧烈湍流）时，全逆滤波的结果不可用了，而维纳滤波的结果比较接近源图像；而且经过模糊退化加高斯噪声后，维纳滤波能消除一些噪声影响，复原图的效果也比较好。改变常数k=0.00025(轻微湍流)时，全逆滤波结果有所改善。但总体来看维纳滤波的效果更好些。实验6 图像恢复-交互式恢复去除周期性噪声授课教师肖红授课班级信息12学时理论48学时+实验16学时实验环境Matlab7.1实验内容：1理

27、解图像恢复的原理2利用交互式恢复去除周期性噪声（实验中为正弦干扰模式）1理解图像恢复的原理图像在形成,传输和记录的过程中,由于受多种原因的影响,图像的质量会有下降,典型表现为图像模糊,失真,有噪声等.这一降质的过程称为图像的退化. 而图像复原试图利用退化现象的某种先验知识(即退化模型),把已经退化了的图像加以重建和复原.其目的就是尽可能地减少或去除在获取图像过程中发生的图像质量的下降(退化),恢复被退化图像的本来面目.2利用交互式恢复去除周期性噪声%Periodic noise removal close all;clear all;I=imread('cameraman.tif

28、9;);I=im2double(I);% 生成正弦噪声 %u=32;v=-32;d=0.3;u1=2*pi*u/256;v1=2*pi*v/256;for j=1:256 for k=1:256 A(j,k)=d*sin(u1*j+v1*k); endend% 叠加生成正弦噪声 %J=I+A;F1=fft2(I);% F1am=abs(F1);% 这条语句没用% F1ph=angle(F1); % 这条语句没用F11=fftshift(F1);F11am=abs(F11);F2=fft2(J);F2am=abs(F2);%当x为复数时，abs（x）返回的不是绝对值而是x的模，保存下来为了后边使

29、用，事实上显示叠加正弦噪声后图像的频谱是用不上这两条语句的。F2ph=angle(F2); %以弧度的形式返回相位角，保存下来为了后边使用F22=fftshift(F2);F22am=abs(F22);subplot(2,3,1),imshow(I);title('原图像');subplot(2,3,2),imshow(J);title('叠加正弦噪声后的图像');subplot(2,3,4),imshow(0.25*log10(1+F11am);subplot(2,3,5),imshow(0.25*log10(1+F22am);% 交互式复原图像 %F2am0

30、=F2am(1,1);F2am(1,1)=0;thresh=max(F2am(:);for j=1:1:256 for k=1:1:256 if F2am(j,k)>thresh/2 F2am(j,k)=0; end endendF2am(1,1)=F2am0;R=F2am.*exp(i*F2ph);%参照章毓晋教材44页公式3.2.4 F(u,v)=|F(u,v)|*expj*相位角R2=fftshift(R);R2am=abs(R2);subplot(2,3,3),imshow(real(ifft2(R);title('复原后的图像')subplot(2,3,6),i

31、mshow(0.25*log10(1+R2am);F3am=abs(fftshift(fft2(J);figure;imshow(0.25*log10(1+F3am);实验7 图像编码-算术编码方法授课教师肖红授课班级信息12学时理论48学时+实验16学时实验环境Matlab7.1实验内容：1理解算术编码的原理2编写程序代码实现算术编码1理解算术编码的原理（1）算术编码是一种从整个符号序列出发，采用递推形式连续编码的方法。（2）在算术编码中，源符号和码字之间的一一对应关系并不存在。一个算术码字要赋给整个信源符号序列，而码字本身确定0和1之间的一个实数区间。（3）例：要编码的是来自一个4-符号信

32、源a1,a2,a3,a4的由5个符号组成的符号序列：b1b2b3b4b5=a1a2a3a3a4。设已知各个信源符号的概率为：p(a1)=0.2,p(a2)=0.2,p(a3)=0.4,p(a4)=0.2。编完最后一个符号得到一个区间0.06752，0.0688，任何1个该区间内的实数，如0.068就可用来表示整个符号序列。function ssbm(str)% 算术编码过程实例ssbm.m%l=0;r=1;d=1;%初始间隔about=. '本实例说明：=>>' '字符串不能太长，程序不加判断，请注意溢出；' '本实例只限定少数字符串；' '实例只是说明一下算术编码过程。'disp(about);%程序限定字符为：a、b、c、dp=0.2 0.2 0.4 0.2;%字符的概率分布，sum(p)=1n=length(str);disp('a b c d ')disp(num2str(p)for i=1:n switch str(i)