数字图像处理上机实验报告(基于Matlab)

1、华北水利水电大学数字图像处理上机实验报告学院：数学与信息科学专业：信息与计算科学姓名：孙志攀学号：201216511指导老师：黄春艳日期：2015年5月4日第一次实验1. 图像的运算实验目的：图像的算术运算在图像处理中有着广泛的应用，它除了可以实现自身所需的算术操作，还能为许多复杂的图像处理提供准备。例如，图像减法就可以用来检测同一场景或物体生成的两幅或多幅图像的误差。我们可以使用MATLAB基本算术符(、等)来执行图像的算术操作，但是在此之前必须将图像转换为适合进行基本操作的双精度类型。图像处理工具箱包含了一个能实现所有数值数据的算术操作的函数集合。列举如下：函数名功能描述Imabsdiif

2、两幅图像的绝对差值Imadd两个图像的加法imcomplment补足一幅图像imdivide两个图像的除法imlincomb计算两幅图像的线形组合lmmultiply两个图像的乘法Imsubtract两个图像的减法使用图像工具箱中的图像代数运算函数无需再进行数据类型间的转换，这些函数能够接受uint8和uintl6数据，并返回相同格式的图像结果。图像的算术运算函数使用以下截取规则使运算结果符合数据范围的要求：超出数据范围的整型数据将被截取为数据范围的极值，分数结果将被四舍五入。例如，如果数据类型是uint8，那么大于255的结果(包括无穷大inf)将被设置为255。无论进行哪一种代数运算都要保

3、证两幅输入图像的大小相等，且类型相同。实验内容：图像的加法MATLAB中调用imadd函数实现图像相加，格式如下：Z=imadd(X,Y)；其中Z=XY。图像的减法图像减法也称为差分方法，MATLAB中调用imsubtract函数实现图像相减。调用格式如下：Z=imsubtract(X,Y)；其中Z=XY。图像的乘法两幅图像进行乘法运算可以实现掩模操作，即屏蔽掉图像的某些部分。一幅图像乘以一个常数通常被称为缩放。如果使用的缩放因数大于1，那么将增强图像的亮度，如果因数小于1则会使图像变暗。MATLAB中调用immultiply函数实现两幅图像相乘。调用格式如下：Zimmultiply(X，Y)

4、；其中Z=XY。图像的除法除法运算可用于校正成像设备的非线性影响。MATLAB中调用imdivide函数进行两幅图像相除。调用格式如下：Z=imdivide(X,Y)，其中Z=XY。图像四则运算因对uint8、uintl6数据，每步运算都要进行数据截取，将会减少输出图像的信息量。图像四则运算较好的办法是使用函数imlincomb。该函数按双精度执行所有代数运算操作，仅对最后的输出结果进行截取，该函数的调用格式如下：Z=imlincomb(A,X,B,Y,C)，其中，Z=AXBYCZ=imlincomb(A,X,C)，其中：ZAXCZ=imlincomb(A,X,B,Y)，其中：Z=AXBY1.

5、 图像的加法基于Matlab的代码如下：clear; I=imread(C:UsersAdministratorDesktop大海.jpg)J=imread(C:UsersAdministratorDesktop恋人.png)K=imadd(I,J);subplot(2,2,1);imshow(I);subplot(2,2,2);imshow(J);subplot(2,2,3);imshow(K);图片素材：运行结果窗口截图：2. 图像的减法基于Matlab的代码如下：sea=imread(C:UsersAdministratorDesktop恋人.png)background=imopen(

6、sea,strel(disk,20);sea2=imsubtract(sea,background);subplot(1,2,1),imshow(sea);subplot(1,2,2),imshow(sea2); 素材图片：运行结果窗口截图：3. 图像的乘法基于Matlab的代码如下：I=imread(C:UsersAdministratorDesktop恋人.png)J=immultiply(I,1.2);K=immultiply(I,0.6);subplot(1,3,1),imshow(I);subplot(1,3,2),imshow(J);subplot(1,3,3),imshow(K)

7、;素材图片：运行结果窗口截图：4. 图像的除法基于Matlab的代码如下：Love=imread(C:UsersAdministratorDesktop恋人.png)I=double(Love);J=0.73*I+100;Ip=imdivide(I,J);subplot(2,2,1),imshow(uint8(I);subplot(2,2,3),imshow(uint8(J);subplot(2,2,4),imshow(uint8(Ip),);运行结果窗口截图：5. 图像的逻辑运算（与、或、非、异或等运算）基于Matlab的代码如下：clear allA=zeros(128);A(40:70,

8、60:100)=1; %生成128*128的零矩阵A，将宽40-70，高60-100的值设为1B=zeros(128);B(50:80,40:70)=1; %生成128*128的零矩阵B，将宽50-80，高40-70的值设为1 C=A(:,:,1); %取单色% D=B(:,:,1);figure(Name,图像逻辑运算) res1=C&D; %C&D%res2=C|D; %C|D%res3=xor(C,D); %xor C D% res3_2=xor(D,C) %xor D C% res4=not(C)¬(D); %C% res4_2=not(C)|not(D); %D% subplo

9、t(2,3,1)subplot(2,3,1)imshow(res1);title(C&D); % &运算subplot(2,3,2)imshow(res2);title(C|D); % |运算subplot(2,3,3)imshow(res3);title(xor C D); % or运算subplot(2,3,4)imshow(res3_2);title(xor D C); subplot(2,3,5)imshow(res4);title(C&D); % 运算subplot(2,3,6)imshow(res4_2);title(C|D);运行结果窗口截图：2. 图像的滤波（空间、频域）实验目

10、的 掌握matlab程序设计语言，掌握matlab基本数据类型、核心函数及辅助函数的使用。掌握理想和线性平滑滤波器的设计方法。实验内容 利用理想和线性平滑滤波器实现图像的滤波；三、实验原理二维理想低通滤波器的传递函数为：D0是指定非负数值，D（u，v）是（u，v）点距频率中心的距离。如果要研究的图像尺寸为M X N，则它的变换也有相同的尺寸。在半径为D0的圆内，所有频率无衰减地通过滤波器，而在此半径之外的所有频率完全被衰减掉。高斯高通滤波器传递函数为：D（u，v）是距傅立叶变换中心原点的距离。D0是截止频率。高斯低通滤波器的傅立叶变换也是高斯的。二维理想高通滤波器的传递函数为：D0是从频率矩形

11、中点测得的截止频率长度，它将以D0为半径的圆周内的所有频率置零，而毫不衰减地通过圆周外的任何频率。但其物理上是不可实现的。高斯高通滤波器传递函数为：高通滤波器能够用高斯型低通滤波器的差构成。这些不同的滤波器有更多的参数，因此能够对滤波器的形状进行更多的控制。四实验设备和仪器1计算机 2. matlab开发平台空间域滤波（举例：线性平滑滤波器）基于Matlab的代码如下：I=imread(C:UsersAdministratorDesktop撑着油纸伞的姑娘.jpg);subplot(231)imshow(I)title(原始图像)I=rgb2gray(I);I1=imnoise(I,salt

12、& pepper,0.02);subplot(232)imshow(I1)title( 添加椒盐噪声的图像)k1=filter2(fspecial(average,3),I1)/255; %进行3*3模板平滑滤波k2=filter2(fspecial(average,5),I1)/255; %进行5*5模板平滑滤波k3=filter2(fspecial(average,7),I1)/255; %进行7*7模板平滑滤波k4=filter2(fspecial(average,9),I1)/255; %进行9*9模板平滑滤波subplot(233),imshow(k1);title(3*3 模板平滑

13、滤波);subplot(234),imshow(k2);title(5*5 模板平滑滤波);subplot(235),imshow(k3);title(7*7 模板平滑滤波);subplot(236),imshow(k4);title(9*9 模板平滑滤波);运行结果窗口截图：频域滤波（举例：理想低通滤波器）基于Matlab的代码如下：RGB = imread(C:UsersAdministratorDesktoplove.png);I0 = rgb2gray(RGB); subplot(2,3,1),imshow(I0);title(原图);I1 = imnoise(I0,gaussian)

14、; %对原图像加噪声subplot(2,3,2),imshow(I1);title(加入噪声后) %将灰度图像的二维不连续Fourier 变换的零频率成分移到频谱的中心s=fftshift(fft2(I1);subplot(2,3,3),imshow(log(1+abs(s),);title(fftshift); M,N=size(s); %分别返回s的行数到M中，列数到N中n1=floor(M/2); %对M/2进行取整n2=floor(N/2); %对N/2进行取整%ILPF滤波（程序中以d0=15为例）d0=50; %初始化d0for i=1:Mfor j=1:N d=sqrt(i-n1

15、)2+(j-n2)2); %点（i,j）到傅立叶变换中心的距离 if d=d0 %点（i,j）在通带内的情况 h(i,j)=1; %通带变换函数 else %点（i,j）在阻带内的情况 h(i,j)=0; %阻带变换函数ends(i,j)=h(i,j)*s(i,j); %ILPF滤波后的频域表示endends=ifftshift(s); %对s进行反FFT移动%对s进行二维反离散的Fourier变换后，取复数的实部转化为无符号8位整数s=uint8(real(ifft2(s); subplot(2,3,4),imshow(h);title(传递函数); %显示GHPF滤波器的传递函数subpl

16、ot(2,3,5),imshow(s); title(ILPF滤波（d0=50); %显示ILPF滤波后的图像运行结果窗口截图：第二次实验1.图像复原算法图像在形成、传输和记录过程中，由于受到多种原因的影响，图像的质量就会有所下降，典型的表现为图像模糊、失真、有噪声等，这一过程称为图像的退化。图像复原是试图利用退化过程的先验知识使已退化的图像恢复本来面目，即根据退化的原因，分析引起退化的环境因素，建立相应的数学模型，并沿着使图像降质的逆过程恢复图像。目的在于消除或减轻在图像获取以及传输的过程中造成的图像品质下降，恢复图像的本来面目。因此，复原技术就是把退化模型化，并采用相反的过程进行处理，以便

17、尽可能复原被退化图像的本来面目。广义上讲，图像复原是一个求逆问题，逆问题经常存在非唯一解，甚至无解。要想恢复全真的景物图像比较困难。为了得到逆问题的有用解，图像复原本身往往需要一个质量标准，即衡量接近全真景物图像的程度，或者说，对图像的估计是否达到最佳的程度。需要有先验知识以及对解的附加约束条件。典型的图像复原是根据图像退化的先验知识建立一个退化模型，以此模型为基础，采用各种逆退化处理的方法进行恢复，使图像质量得到改善。图像退化的数学模型一般来说，图像的生成可以简单地被描述为如下数学模型：f(x,y)是成像景物，H是综合退化因子，g(x,y)是退化图像。图像f(x,y)可以表示为：用卷积符号*

18、表示为：因此还有：式中，是像素点的特性函数，为冲击响应。假定成像系统是线性移不变系统：退化模型如图所示+不考虑加性噪声：考虑加性噪声：卷积等同于频域内乘积：逆滤波复原逆滤波复原法也叫做反向滤波法，其主要过程是首先将要处理的数字图像从空间域转换到傅里叶频域中，进行反向滤波后再由频率域转回到空间域，从而得到复原的图像信号5。1在不考虑噪声的情况下：上式两边进行傅里叶变换得则原始图像= 然后进行傅里叶逆变换，就可以得到原始图像。由此可看出，如果已知退化图像的傅里叶变换和“滤波”传递函数，则可以求得原始图像的傅里叶变换，经反傅里叶变换就可以求得原始图像f(x,y)，这就是逆滤波法的基本原理。但在实际中

19、用逆滤波法存在病态的情况：当H(u,v)=0时，或非常小的数值点上，F(u,v)将变成无穷大或非常大的数。2在有噪声的情况下：逆滤波原理可以写成：G(u,v)=F(u,v)H(u,v)+N(u,v)写成逆滤波的方式：F(u,v)=F(u,v)+ 但实际用逆滤波存在病态的情况：噪声存在，当H(u,v)很小或为零时，则噪声被放大。这意味着退化图像中小噪声的干扰在H(u,v)较小时，会对逆滤波恢复的图像产生很大的影响，有可能使恢复的图像和f(x,y)相差很大，甚至面目全非。实验证明，当退化图像的噪声较小，即轻度降质时，采用逆滤波复原的方法可以获得较好的结果。通常，在离频率平面原点较远的地方数值较小或

20、为零，因此图像复原在原点周围的有限区域内进行，即将退化图像的傅里叶频谱限制在没出零点而且数值又不是太小的有限范围内。维纳滤波复原逆滤波比较简单，但没有清楚地说明如何处理噪声，而维纳滤波综合了退化函数和噪声统计特性两个方面进行复原处理。维纳滤波是维纳在1949年提出的，并应用于一维平稳时间序列，获得了满意的结果。这是最早也是最著名的线性滤波技术。 采用维纳滤波是假设图像信号可以近似看成平稳随机过程的前提下，按照使f(x,y)和f(x,y)之间的均方误差达到最小的准则函数来实现图像复原的，即式中，E（）代表求期望值。因此维纳滤波又称为最小均方误差滤波器。维纳滤波需要假定下述条件成立：1、 系统为线

21、性空间移不变系统。2、 退化图像、原始图像、噪声都是均匀随机场，噪声的均值为零，且与图像不相关。维纳滤波的复原滤波函数，即滤波器的传递函数为： 没有噪声时，维纳滤波退化为逆滤波。有噪声时，维纳滤波利用信噪功率比对恢复过程进行修正，在信噪功率比很小的区域内，P(u,v)的值也很小，这使恢复图像较小地依赖于退化图像。在H(u,v)很小或等于零时，P(u,v)的分母不为零，维纳滤波没有病态问题7。在实际系统中，维纳滤波经常用下式近似：，K为特殊常数。 的维纳滤波要求未退化图像和噪声的功率必须是已知的。虽然用近似的方法能得到好的结果，但功率谱比常数K的估计一般没有合适的解。基于Matlab的代码如下：

22、clear; I=imread(C:UsersAdministratorDesktop明晓妹妹原始图.jpg)imshow(I); I=rgb2gray(I); %将原图像转化为黑白图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

23、(I0),gaussian,0,0.001) subplot(2,2,2);imshow(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

24、)+K); endendF2=H1.*F0; I3=ifft2(fftshift(F2); subplot(2,2,4);imshow(uint8(I3);窗口截图如下：经过仿真，如上图所示，可以看出逆滤波复原与维纳滤波复原的区别和联系。维纳滤波后虽然仍有一些噪声存在，但已经和原图很接近了。因为原图像和噪声函数都是已知的，可以正确的估算参量。使用Lucy-Richardson算法的迭代非线性复原L-R算法是一种迭代非线性复原算法，它是从最大似然公式印出来的，图像用泊松分布加以模型化的。当下面这个迭代收敛时模型的最大似然函数就可以得到一个令人满意的方程：*代表卷积，代表未退化图像的估计，g和h和

25、以前定义一样。这个算法的本质是显而易见的。它的非线性本质是在方程右边用来除产生的4。在IPT中，L-R算法是由名为deconvlucy的函数完成的，此函数的语法为fr= deconvlucy(g,PSF,NUMIT,DAMPAR,WEIGHT)其中，fr代表复原的图像，g代表退化的图像，PSF是点扩散函数，NUMIT为迭代次数（默认为10次），DAMPAR是一个标量，它指定了结果图像与原图像g之间的偏离阈值。WEIGHT是一个与g同样大小的数组，它为每一个像素分配一个权重来反映其重量。基于Matlab的代码如下：clear; I=imread(C:UsersAdministratorDeskt

26、op孙志攀帅比原始图.jpg)PSF=fspecial(gaussian,5,5) ;Blurred=imfilter(I,PSF,symmetric,conv); V=.003; BN=imnoise(Blurred,gaussian,0,V); luc=deconvlucy(BN,PSF,5); figuresubplot(2,2,1);imshow(I);title(原始图像);subplot(2,2,2);imshow (Blurred);title(模糊后的图像);subplot(2,2,3);imshow (BN);title(加噪后的图像);subplot(2,2,4);imsh

27、ow (luc);title(恢复后的图像);窗口截图如下：2.读入一幅彩色图像，显示各层图像。（方法一）基于Matlab的代码如下：I=imread(C:Documents and SettingsAdministrator桌面1.jpg);%如果是win7系统则读取路径为I=imread(C:UsersAdministratorDesktop1.jpg);IR=I;IG=I;IB=I;IR(:,:,2)=0;IR(:,:,3)=0;IG(:,:,1)=0;IG(:,:,3)=0;IB(:,:,1)=0;IB(:,:,2)=0;figure;subplot(2,2,1);imshow(I,)

28、;title(I);subplot(2,2,2);imshow(IR,);title(IR);subplot(2,2,3);imshow(IG,);title(IG);subplot(2,2,4);imshow(IB,);title(IB);（方法二）基于Matlab的代码如下：I=imread(C:Documents and SettingsAdministrator桌面1.jpg);R = I(:,:,1);G = I(:,:,2);B = I(:,:,3);figure;subplot(2,2,1);imshow(I,);title(I);subplot(2,2,2);imshow(R,

29、);title(R);subplot(2,2,3);imshow(G,);title(G);subplot(2,2,4);imshow(B,);title(B);3.图像的变换（小波或FT变换）clear allclose allA=imread(C:Documents and SettingsAdministrator桌面1.jpg);%读入并且显示出一个图像文件subplot(1,2,1);imshow(A);title(原始的图像);%显示原始图像作为对照if length(size(A)=3A=rgb2gray(A);endsubplot(1,2,2);imshow(A);title(灰度图像);%对灰度图像进行傅里叶变换并输出频谱A2=fft2(A);A2=fftshift(A2);%将图像进行二维傅里叶变换figure,imshow(log(abs(A2)+1),0,12);%显示傅里叶变换后的图像title(傅里叶变