河北工业大学数字图像处理之matlab实验报告.doc
实验一一、 实验要求1、熟悉MATLAB软件的开发环境、基本操作以及图像处理工具箱，为编写图像处理程序奠定基础。2、掌握二进制（黑白）图像、灰度图像、彩色图像读、写和显示的基本方法，分别选择以上几种类型的图像，观察其图像类型参数。3、说明以下程序的显示结果为何是一幅几乎全黑的图像。clear; close all;myi=zeros(20,20);myi(2:2:18,2:2:18)=1;myi=uint8(myi);figure, imshow(myi,'notruesize');编写程序将图一中orangutan_1.tif图片拉伸一倍后形成orangutan_2.tif所示的图片。4、应用MATLAB（或C）语言编写一幅灰度图像直方图统计程序，并选择一幅图像显示其直方图，将结果与MATLAB图像处理工具箱中提供的灰度直方图函数imhist的处理结果进行比较。5、利用以上编写的程序，估算图一所示图像iris.tif中的瞳孔半径（以像素为单位）。 orangutan_1.tif orangutan_2.tif iris.tif图一二、 实验内容1，通过安装matlab软件，熟悉MATLAB软件的开发环境、基本操作以及图像处理工具箱2,图像的读、写、显示和图像类型参数查看通过 I= imread(FILENAME),读取图像，注意FILENAME=图像文件所在的位置+文件名；1） 通过imwrite(a,FILENAME,fmt),可以把a图像写进FILENAME图像，即是替换掉FILENAME，但是文件图像格式不替换。2） 通过imshow(I),可以把I图像显示出来。3） 例如：在matlab软件中，运行如下命令： I = imread('E:orangutan_2.tif'); figure,imshow(I) I1 = imread('E:iris.tif'); figure,imshow(I1) imwrite(I,'E:iris.tif','tif') I2 = imread('E:iris.tif'); figure,imshow(I2) 原图： orangutan_2.tif iris.tif运行结果：可以看到orangutan_2.tif已成功写入 iris.tif中。4） 在matlab的可以查看各图像类型的参数先在workspace中导入一幅名为xx.jpg的真彩图，然后在命令窗口输入如下命令： I8= imread('E:iris.tif');I64=double (I8)/255;imshow(I64);将一幅名为iris.tif的灰度图像转换成二值图像，然后观察workspace中显示的二值图像、灰度图像、和真彩图像的参数如下图所示3,全黑图像的解释和拉伸图像程序1） 、解释 由程序myi=uint8(myi);可知此程序编写的图像在计算机中采用的是8位存储，即灰度值为0255，其中0为黑，255为白，中间值为由黑向白过度的值。而程序myi=zeros(20,20);构造了一个20行20列的0值矩阵，程序myi(2:2:18,2:2:18)=1;将1赋值给矩阵的2，4，6，8，10，12，14，16，18行和2，4，6，8，10，12，14，16，18列，即这些位置的值为1。因为0，1都在人眼识别的黑色范围内，固图像是全黑的，如果将myi(2:2:18,2:2:18)=1;改为myi(2:2:18,2:2:18)=255;就会发现图像是黑白相间的条格。2） 、纵向二倍拉伸图像代码T1 = maketform('affine',1 0 0; 0 2 0; 0 0 1);I = imread('E:orangutan_1.tif');I1 = imtransform(I,T1);imshow(I),figure,imshow(I1)4,灰度图像直方图统计程序1）编写的程序img=imread('E:iris.tif');imshow(img);M,N=size(img);img=double(img);h=zeros(256,256);for i=1:M, for j=1:N, f=img(i,j); h(f+1)=h(f+1)+1; endendfigure;bar(h);grid;2）直接调用MATLAB图像处理工具箱中提供的灰度直方图函数imhistI = imread('E:iris.tif');imhist(I)3） 显示分别如下图1和图2图1图25,估算图像iris.tif中的瞳孔半径（以像素为单位）方法：从图1中带网格线的图中估算像素具体方法：100200之间的单位个数为8，则瞳孔面积大概为80000像素，利用求半径公式S=3.14R2得到半径R=(S/3.14)=(80000/3.14)160像素实验二应用MATLAB工具箱自带的图像函数或自己编程进行图像处理，包括显示图像、格式转换、图像傅立叶变换、图像灰度调整、对比度增强、直方图均衡、邻域平均法、中值滤波法、图像代数运算、边缘增强、伪彩色增强、假彩色合成、二值化、边缘检测等；分析处理结果，掌握方法和原理，巩固所学内容。部分示例以及需要完成的具体实验题目：1）二维傅立叶变换示例：选择一幅图像，对其进行离散傅立叶变换并显示离散傅立叶谱。% 根据以下生成的图像矩阵数据进行图像的傅立叶变换并显示其频谱clear;close all;f=zeros(30,30);f(5:24,13:17)=1;figure,imshow(f,'notruesize');F=fft2(f);% 可视化二维傅立叶变换结果的常用方法是使用对数，以便更详细地观察在0附近区域的细节figure,imshow(log(abs(F),-1 5,'notruesize');colorbar;colormap(jet);pause% 以上可视化f的DFT振幅谱存在以下两个问题：% 一、傅立叶变换采样比较粗糙% 原因：快速傅立叶变换算法只能处理大小为2的幂次的矩阵，而本例中矩阵维数不是2的幂次% 解决办法：在计算DFT时可以通过对f进行零填充来获得较好的傅立叶变换采样% 二、零频率系数显示在图形的左上角而不是传统的中心位置% 原因：计算机系统在执行傅立叶变换算法时是以图像的左上角为坐标原点的% 解决办法：用函数fftshift通过将图像FF的四个象限进行交换，使零频率系数位于图形的中心FF=fft2(f,256,256); % 对f进行零填充，得到一个256X256矩阵，然后再计算离散傅立叶变换figure,imshow(log(abs(FF),-1 5,'notruesize');colorbar;colormap(jet);pauseFFF=fftshift(FF); % 使用函数fftshift对第二个问题进行修正解决figure,imshow(log(abs(FFF),-1 5,'notruesize');colorbar;colormap(jet);return代码：% 根据以下生成的图像矩阵数据进行图像的傅立叶变换并显示其频谱clear;close all;f=zeros(30,30);f(5:24,13:17)=1;figure,imshow(f,'notruesize');F=fft2(f);% 可视化二维傅立叶变换结果的常用方法是使用对数，以便更详细地观察在0附近区域的细节figure,imshow(log(abs(F),-1 5,'notruesize');colorbar;colormap(jet);pauseFF=fft2(f,256,256); % 对f进行零填充，得到一个256X256矩阵，然后再计算离散傅立叶变换figure,imshow(log(abs(FF),-1 5,'notruesize');colorbar;colormap(jet);pauseFFF=fftshift(FF); % 使用函数fftshift对第二个问题进行修正解决figure,imshow(log(abs(FFF),-1 5,'notruesize');colorbar;colormap(jet);return2）直方图均衡化：利用Matlab提供的可进行图像直方图修正的函数，自己选择几幅直方图不均匀的图像（如图二pout.tif），对这些图像进行直方图均衡处理，显示处理前后的图像以及它们的灰度直方图，体会直方图均衡化算法的特点。代码：I=imread('pout.tif');H=histeq(I);%对源图像的直方图进行均衡化subplot(2,2,1),imshow(I);%显示源图像title('源图像');subplot(2,2,2),imshow(H);%显示均衡化后的图像title('均衡化后的图像');subplot(2,2,3),imhist(I);%显示原始图像的灰度直方图title('源图像的直方图');subplot(2,2,4),imhist(H);%显示均衡化后图像的直方图title('均衡化后的图像直方图');3）灰度修正技术：按照教材64页上的公式（4.1.6） (4.1.6)以及以下公式（1），分别编程实现图像的分段线性灰度变换和非线性灰度变换。 (1)利用以上自己编写的代码对图三所示的实验图像（bottle.tif）进行灰度变换，使拍摄主体瓶子的花纹更加清晰地表现出来，改善视觉效果，显示增强前、后的图像以及它们的灰度直方图，并考虑利用Matlab提供的imadjust函数是否能完成同样的工作。代码：F=imread('bottle.tif');raw,col=size(F);%读取图像的行列数max=double(F(1,1)+1;%转换为double类型for i=1:raw%循环查找灰度的最大值 for j=1:col k=double(F(i,j)+1; if max<kmax=k;end endenddisp('灰度最大值：');max%将最大值输出disp('-首先进行线性变换-');fo=input('请输入源图像的灰度范围：');%输入格式：num1 num2a=fo(1);%将第一个数赋值给ab=fo(2);%将第二个数赋值给bif b>max%判断b是否超过max disp('输入源图像的范围不能大于最大值！');endre=input('请输入变换后的灰度范围：');c=re(1);d=re(2);if d>256%判断d是否超过256 disp('输入变换后的范围不能超过255！');endG=zeros(raw,col);%初始一个图像的行列数for i=1:raw%对源图像的每一个节点进行三段线性灰度变换 for j=1:col k=double(F(i,j)+1; if k<a G(i,j)=(c/a)*F(i,j); elseif k<b G(i,j)=(d-c)/(b-a)*(F(i,j)-a)+c; else G(i,j)=(255-d)/(max-b)*(F(i,j)-b)+d; end endendout=uint8(round(G-1);%将图像转换回uint8格式figure,imshow(F);%显示源图像title('源图像');figure,imshow(out);%显示线性变换后的图像title('线性变换后');figure,imhist(F);%源图像的直方图title('源图像的直方图');figure,imhist(out);%线性变换后的直方图title('线性变换后的直方图');disp('-下面进行非线性变换，其变换公式为g=c*fy-');y=input('请输入指数y：');c=input('请输入乘数c：');g=zeros(raw,col);%初始第二个图像的行列数for i=1:raw%对源图像的每一个节点进行非线性灰度变换 for j=1:colk=double(F(i,j)+1; g(i,j)=c*(ky);%根据公式求得各个节点的值 endendmax=g(1,1);for i=1:raw%循环查找g(i,j)的最大值 for j=1:col if max<g(i,j)max=g(i,j);end endendscale=max/256;for i=1:raw%将所求的g(i,j)等比例转换到1-256的像素值; for j=1:colg(i,j)=round(g(i,j)/scale); endendout2=uint8(round(g-1);figure,imshow(out2);%显示非线性变换后的图像title('非线性变换后');figure,imhist(out2);%非线性变换后的直方图title('非线性变换后的直方图');灰度最大值：max = 256-首先进行线性变换-请输入源图像的灰度范围：70 120请输入变换后的灰度范围：20 200Warning: Image is too big to fit on screen; displaying at 75% scale.> In truesize>Resize1 at 308 In truesize at 44 In imshow at 161Warning: Image is too big to fit on screen; displaying at 75% scale.> In truesize>Resize1 at 308 In truesize at 44 In imshow at 161-下面进行非线性变换，其变换公式为g=c*fy-请输入指数y：3请输入乘数c：3Warning: Image is too big to fit on screen; displaying at 75% scale.> In truesize>Resize1 at 308 In truesize at 44 In imshow at 161>>4）图像平滑（去噪）：（1）选择一幅图像（如图四lena512.jpg），叠加零均值高斯噪声，然后分别利用邻域平均法和中值滤波法对该图像进行滤波，显示滤波前后的图像，比较各滤波器的滤波效果。（2）选择一幅图像，叠加椒盐噪声，寻则合适的滤波器将噪声滤除。（3）选择一幅叠加了零均值高斯噪声的图像，设计一种处理方法，即能去噪声又能保持边缘清晰。（选作）（4）利用图像空间低通滤波法，对于给定的任意图像和低通滤波器（模板），编写一个较为通用的函数实现对输入图像的平滑处理。（选作）代码：（1）x=imread('lena512.jpg');subplot(2,2,1);imshow(x);%显示源图像title('源图像');y=imnoise(x,'gaussian',0);%添加0均值的高斯噪声subplot(2,2,2);imshow(y);%显示添加高斯噪声的图像title('添加高斯噪声');mo=fspecial('average');%产生一个3*3的方形平均滤波模板x1=imfilter(x,mo,'replicate');%对图像进行邻域平均法处理%'replicate'指卷积填充边缘时用复制边界的值来扩展subplot(2,2,3);imshow(x1);%显示邻域平均法处理后的图像title('邻域平均法处理');x2=medfilt2(y);%对高斯躁声中值滤波subplot(2,2,4);imshow(x2);%显示中值滤波后的图像title('对高斯噪声进行中值滤波');（2）M=imread('lena512.jpg');subplot(1,3,1); imshow(M); %显示原始图像 title('源图像');P=imnoise(M,'salt & pepper'); %加入椒盐躁声P=imnoise(M,'salt & pepper'); %两次叠加椒盐躁声subplot(1,3,2);imshow(P); %两次叠加椒盐躁声后显示图像 title('加入椒盐噪声');h=medfilt2(P); %对椒盐躁声中值滤波subplot(1,3,3);imshow(h);title('中值滤波后');5）边缘检测技术：（1）选择三种不同的边缘检测算子对一幅图像（如图五rice.png）进行边缘检测，显示检测结果。（2）对含噪图像进行边缘检测，比较各边缘检测算子对噪声的敏感性，并提出抗噪性能较好的边缘检测方法。（3）编制一个比较通用的边缘提取函数。通过输入不同的参数，能够实现Roberts算子、Prewitt算子、Sobel算子以及Laplacian算子边缘检测。（选作）代码：（1）I=imread('rice.png');imshow(I); %显示源图像title('源图像');A1=edge(I,'sobel'); %用SOBEL算子进行边缘检测figure,subplot(1,3,1);imshow(A1); %SOBEL算子图像显示title('SOBEL算子');A2=edge(I,'roberts'); %用Roberts算子进行边缘检测subplot(1,3,2);imshow(A2); %Roberts算子图像显示title('Roberts算子');A3=edge(I,'prewitt'); %用prewitt算子进行边缘检测subplot(1,3,3);imshow(A3); %prewitt算子图像显示title('prewitt算子'); （2）I=imread('rice.png');P=imnoise(I,'salt & pepper'); %加入椒盐噪声subplot(2,2,1);imshow(P); %椒盐噪声图像显示title('椒盐噪声图像');B1=edge(P,'sobel'); %用SOBEL算子进行边缘检测subplot(2,2,2);imshow(B1); %SOBEL算子图像显示title('SOBEL算子');B2=edge(P,'roberts'); %用Roberts算子进行边缘检测subplot(2,2,3);imshow(B2); %Roberts算子图像显示title('Roberts算子');B3=edge(P,'prewitt'); %用prewitt算子进行边缘检测subplot(2,2,4);imshow(B3); %prewitt算子图像显示title('prewitt算子');实验心得：其实刚开始我只是听老师课上听讲，但是终究无法深刻理解计算机图像具体的处理过程，通过这次上机实验，不会的上网百度，然后请教同学，才对图像处理有了更深的理解。
