数字图像处理
—及MATLAB实现
实验指导参考答案
目录
实验1 MATLAB软件入门及数字图像处理工具箱（2学时） .................................................. 2 实验3 MATLAB实现数字图像的变换技术（2学时） ............................................................ 20 实验4 MATLAB实现数字图像增强（2学时） ...................................................................... 24
实验5 MATLAB实现图像的复原（2学时） ................................................................... 36 实验6 MATLAB实现图像分割（2学时） ......................................................................... 39 实验7 MATLAB实现彩色图像处理（2学时） ................................................................. 55 实验8 MATLAB实现数字图像其它处理（2学时）.......................................................... 60
实验1 MATLAB软件入门及数字图像处理工具箱（2学时）
【实验内容】 一、二两项略
一、运行MATLAB，要求完成： （1)、了解MATLAB工作环境 （2)、运用MATLAB帮助系统，查询imread,imwrite,imshow这3个工具箱函数的功能，给出解释
（3)、在桌面上建立一个文件夹，以本人姓名拼音缩写为文件名，将该文件夹通过set path设置，放到MATLAB的搜索路径下。截图最后操作的界面，保存至该文件
夹下。
二、了解MATLAB的一些简单的工具箱函数使用，在MATLAB下完成一些图像处理操作。
1、在命令窗口下依次输入如下带下划线的英文命令，无需写注释，对比查看结果： （1）、读入并显示一幅图像
>> clear； close all； %首先清除 MATLAB 所有的工作平台变量，关闭已打开的图形窗口。 >> I=imread（′pout.tif′）；%然后使用图像读取函数 imread 来读取一幅图像。假设要读取图像 pout.tif，并将其存储在一个名为 I的数组中 >> imshow（I） %使用 imshow 命令来显示数组 I 观察显示结果
（2）、检查内存中的图像
>> Whos %使用 whos命令来查看图像数据 I是如何存储在内存中， MATLAB 做出的响应如下： Name Size Bytes Class
I291 ×240 69840 uint8 array
Grand totalis69840 elementsusing69840 bytes （3）、实现直方图均衡化
如pout.tif图像对比较低，为了观察图像当前状态下亮度分布情况，可以通过使用 imhist函数创建描述该图像灰度分布的直方图。 >> figure，imhist(I）；%首先使用 figure命令创建一个新的图像窗口，避免直方图覆盖图像数组 I的显示结果。 >> I2 =histeq（I）； figure，imshow（I2）；%可以通过调用 histeq函数将图像的灰度值扩展到整个灰度范围中，从而达到提高数组 I的对比度。此时修改过的图像数据保存在变量 I2中。然后，再通过调用imhist函数观察其拓展后的灰度值的分布情况。 （4）保存图像
将新调节后的图像 I2保存到磁盘中。假设希望将该图像保存为 PNG 格式图像文件，使用 imwrite函数并指定一个文件名，该文件的扩展名为.png。 >> imwrite（I2，′pout2.png′）； （5）检查新生成文件的内容
利用 imfinfo函数可以观察上述语句写了什么内容在磁盘上。值得注意的是：在 imfinfo函数语句行末尾不要加上分号，以保证 MATLAB 能够显示图像输出结果；另外，要保证此时的路径与调用 imwrite时的路径一致。 >> imfinfo（′pout2.png′）
2、对一幅灰度图像 rice.tif进行一些较为高级的操作，在文本编辑器下输入如下程序，保存成test.m文件，按要求做，调试成功，并观察显示结果： clear； close all；
I=imread（′rice.tif′）； imshow（I）；%此处加一个断点，观察结果 background =imopen（I，stre（′disk15））；%实现了背景亮度的估计 I2 =imsubtract（I，background）；%得到背景一致的图像 figure，imshow（I2）；%此处加一个断点，观察结果 I3 =imadjust（I2，stretchlim（I2）,［0,1］）；%修改后图很暗，调节对比度 figure，imshow（I3）；%此处加一个断点，观察结果 level=graythresh（I3）； bw =im2bw（I3，level）；%通过阈值法实现图像的二值化 figure，imshow（bw）；%此处加一个断点，观察结果 ［labeled，numobjects］=bwlabel（bw，4）； % labelcomponents %找出米粒的个数： numobjects=80 grain =imcrop（labeled）；%使用imcrop选择标记区域
RGB_label(grain,@spring,’c’,’shuffle’);%伪彩色显示标记矩阵 Figure,imshow(RGB_label) %此处加一个断点，观察结果 graindata=regionprops（labeled，′basic′） allgrains=［graindata.area］；
max（allgrains）%找最大的那粒米的尺寸
biggrain =find（allgrains==695）%返回最大的那粒米的标记号，%biggrain =68 mean（allgrains）%米粒的平均尺寸 hist（allgrains，20）；%绘制一个包含20柱的直方图来说明米粒大小的分布 三、编程综合练习
1、试对lena图像分别进行4和16倍减采样，查看其减采样效果，分析。
提示：二维图像处理，行每隔2，列每隔2，抽取，则实现4倍减采样；行每隔4，列每隔4，抽取，则实现16倍减采样 参考程序如下： clear
a=imread('lena.bmp'); b=rgb2gray(a); figure
subplot(1,3,3),imshow(a),title('原图像'); m=[2 4]; for mn=1:2
[width,height]=size(b);
quartimage=zeros(floor(width/(m(mn))),floor(height/(m(mn)))); k=1; n=1;
for i=1:(m(mn)):width
for j=1:(m(mn)):height
quartimage(k,n)=b(i,j); n=n+1; end k=k+1; n=1; end
subplot(1,3,mn)
imshow(uint8(quartimage)); mm=m(mn)*m(mn);
title([int2str(mm) ,'倍减采样图像']) end
图1.1 减采样显示结果
2、任选一副彩色图像，将其转换成256级灰度图像，64级灰度图像，16级灰度图像，2级灰度图像，显示并存储显示结果图，观察转换后图像的变换，并分析。（即灰度级的量化） clear
a=imread('lena.bmp'); b=rgb2gray(a); figure
subplot(2,3,1)
imshow(b),title('256级原图像') [width,hei]=size(b); img64=zeros(width,hei); img32=zeros(width,hei); img8=zeros(width,hei); img2=zeros(width,hei); for i=1:width for j=1:hei
img64(i,j)=floor(b(i,j)/4); end end
subplot(2,3,2)
imshow(uint8(img64),[0,63]);title('64级图像') for i=1:width for j=1:hei
img32(i,j)=floor(b(i,j)/16); end end
subplot(2,3,3)
imshow(uint8(img8),[0,7]);title('8级图像') for i=1:width for j=1:hei
img2(i,j)=floor(b(i,j)/128); end
end
subplot(2,3,4)
imshow(uint8(img2),[0,1]);title('2级图像')
【思考题】 1、数字图像在MATLAB中以数据矩阵方式表达，如何通过MATLAB语句截取图像的一部分矩形区域进行分析处理。 答：由于数字图像在MATLAB中以数据矩阵方式表达，因此可用截取矩阵的方式来实现对图像的截取，如原图像大小为A(256,256),则A(1:100,1:200)则实现了对图像截取（1-100行，1-200列）进行读取。 2、图像的数字化工作包括哪两个方面？
答：图像的数字化工作包括采样和量化两个方面，采样指的是空间域的离散化，量化指的是灰度(亮度)的离散化。
3、图像工程包含哪3个层次？你最感兴趣的是图像处理哪个方面的内容？ 答；图像工程包含由低到高图像处理，图像分析，图像理解3个层次。
实验2 数字图像的MATLAB基本操作（2学时）
【实验内容】
一、图像读取，显示，保存
1、利用imread( )函数从内存中读取一幅绘有花朵的tif图像，假设其名为flower.tif，存入一个数组I中；
利用whos 命令提取该读入图像flower.tif的基本信息； 利用imshow()函数来显示这幅图像，并添加颜色条；
利用imfinfo函数来获取图像文件的压缩，颜色等等其他的详细信息； 利用imwrite()函数来压缩这幅图象，将其保存为一幅压缩了像素的jpg文件,放在桌面上，并显示出来，设为flower.jpg；语法：imwrite(原图像，新图像，‘quality’,q), q取0-100。 程序段及显示结果如下：
（1）I=imread('flower.tif'); Imshow(I)
（2）whos
Name Size Bytes Class I 296x394x3 349872 uint8 （3）imfinfo('flower.tif') ans =
Filename: 'flower.tif'
FileModDate: '09-十月-2013 10:35:38' FileSize: 48384 Format: 'jpg' FormatVersion: '' Width: 394 Height: 296 BitDepth: 24
ColorType: 'truecolor' FormatSignature: '' NumberOfSamples: 3
CodingMethod: 'Huffman' CodingProcess: 'Sequential' Comment: {}
（4）imwrite(I,'flower.jpg','quality',10);
（5）imwrite(I,'flower.bmp','quality',10);
Attributes
2、用函数imread()读入图像forest.tif 和trees.jpg,使用imshow命令在同一图形窗中显示两者；再打开另一个图形窗，利用subimage函数实现在同一个图形窗口中显示该两幅图像，对比前后两者的显示效果，给出解释。 程序如下：
I1=imread(‘forest.tif’); I2=imread(‘trees.jpg’) figure,
subplot(2,1,1),imshow(I1); subplot(2,1,2),imshow(I2); figure,
subplot(2,1,1),subimage(I1); sublot(2,1,2),subimage(I2);
对比前后两者的显示效果，imshow显示的图片有问题，而subimage则显示正常。原因:subimage可以在一个界面上同时显示几幅不同调色板的图像 3、总结imshow命令的使用方式，并给出简单实例语句。（如何显示索引图像，RGB图像，二进制图像，非图像矩阵等） 总结如下： 1、imshow(I)%显示图像最基本命令,imshow(I,[low,high])%按照最大灰度范围显示图像
2、imshow(I,map)%按照调色板map,显示索引图像 3、imshow（I,N）%显示具有N个灰度级的灰度图像 4、imshow(BW)%显示二值图像BW 可参考教材(余成波)3.3节内容 二、图像类型转换
1、编写一个m文件实现如下功能：读入任意一幅RGB图像，选择合适的函数将其变换为灰度图像和二值图像，并在同一个窗口内分成三个子区域来分别显示RGB图像，灰度图像，二值图像。
程序及显示结果如下： RGB=imread('xianhua.jpg'); Y=rgb2gray(RGB); BW=im2bw(RGB,0.4);
subplot(1,3,1),imshow(RGB); subplot(1,3,2),imshow(Y); subplot(1,3,3),imshow(B);
三、图像运算 1、线性点运算
(1)、给出如下程序
rice=imread('rice.png'); I=double(rice); A=0.43; B=60; J=I*A+B;
rice2=uint8(J);
subplot(1,2,1),imshow(rice); subplot(1,2,2),imshow(rice2);
运行该程序，观察显示效果，修改以上线性点运算程序，改变A和B值，分析A的值大于1和小于1时，B的值正负时，对原图像的影响。 参考程序如下： 程序如下：
（1）rice=imread('rice.png'); （2）rice=imread('rice.png'); I=double(rice); J=60;
J=I*0.43+60; subplot(1,2,1),imshow(rice); rice2=uint8(J) ; subplot(1,2,2),imshow(rice2); subplot(1,2,1),imshow(rice);
subplot(1,2,2),imshow(rice2);
（3） J=I*0.22+60; （4）J=I*0.68+60;
（5）J=I*0.43+30; （6） J=I*0.43+100;
（2）、任选一副图像或下图进行点运算，实现图像变亮、变暗和负片效果，在同一个窗口内分成四个子窗口来分别显示，注上文字标题。
图2.10 例图 参考程序如下:
mei=imread('mei.jpg');
J=mei*0.43+100; K=mei*0.43+30; Z=255-mei;
subplot(1,4,1),imshow(mei); subplot(1,4,2),imshow(J); subplot(1,4,3),imshow(K); subplot(1,4,4),imshow(Z);
2、代数运算 （1）、选取两幅大小相同的灰度图像A,B，或将不同大小的图像截取成相同大小后，将两幅图像执行加法运算得到C。 （2）、选取以上相加后的图像C作为源图像，将该混合图像与背景图像A或B做减法,显示结果。 （3）、选取一副灰度图像A，设置一个掩模模板B，需要保留的区域，图像值为1，而在需要被抑制的区域，则设为0，完成A.*B,显示结果。 提示：模板构造时，大小与源图像一致，如：B=zeros（256,256）;B(40:200,40:200)=1,就构成一个模板。
根据以上结果，分析图像加法，减法，乘法运算的应用特点 (1)、
I=imread('A.jpg'); J=imread('B.jpg'); K=imadd(I,J);
subplot(1,3,1),imshow(I);title(‘原图1’); subplot(1,3,2),imshow(J); title(‘原图2’); subplot(1,3,3),imshow(K); title(‘相加’) 结果如下：
（2）、相减 %K为上题图
I=imread('A.jpg'); J= imsubtract(K,I);
subplot(1,3,1),imshow(I);title(‘原图1’); subplot(1,3,2),imshow(K); title(‘相加图’); subplot(1,3,3),imshow(J2); title(‘相减后图’)
3、几何运算 （1）、对任意一幅图像调用imresize函数，分别实现对原图的2倍比例放大，0.8倍比例缩小，非比例放大和缩小(人为定义调整后的图像长宽尺寸)，并分别显示结果，对比分析不同的插值方法效果。
参考程序如下：
A=imread(‘lena.jpg’);%假设原图为256*256 B=imresize(A,2);%2倍放大 C=imresize(A,0.8);%0.8倍缩小
D=imresize(A,[300 300]);%非比例放大到300*300 E=imresize(A,[200 200]);%非比例缩小到200*200
%在以上几条语句的基础上，也可在Imresize后面添上一个参数method，指明插值方式，如’nearest’,’bilinear’,’bicubic’ subplot(2,2,1),imshow(B),title(‘放大2倍’) subplot(2,2,2),imshow(C),title(‘0.8倍缩小’) subplot(2,2,3),imshow(D),title(‘非比例放大’) subplot(2,2,4),imshow(E),title(‘非比例缩小’) （2）、任选一副图像，使用imrotate函数，分别顺时针60度，逆时针90度旋转，选择裁剪性旋转或非裁剪性旋转，在同一个窗口中显示处理后的图像。 参考程序如下：
I=imread(‘hua.jpg’);
J11=imrote(I,-60,’bilinear’);%顺时针旋转60度，非裁剪型旋转 J12=imrote(I,-60,’bilinear’，‘crop’);%顺时针旋转60度，裁剪型旋转 J21=imrote(I,90,’nearest’);%逆时针旋转90度，非裁剪型旋转 J22=imrote(I,90,’nearest’，‘crop’);%逆时针旋转90度，裁剪型旋转
subplot(2,2,1),imshow(J11),title(‘放大2倍’) subplot(2,2,2),imshow(J12),title(‘0.8倍缩小’) subplot(2,2,3),imshow(J21),title(‘非比例放大’) subplot(2,2,4),imshow(J22),title(‘非比例缩小’) 显示结果如下：
（3)、任选一幅图像，调用imcrop对该图片进行剪裁，显示处理后的图像。 I=imread(‘hua.jpg’); J=imcrop;
subplot(2,1,1),imshow(I);title(‘原图’); subplot(2,1,2),imshow(J);title(‘剪裁后图’); （4）、任选一幅图像，编程实现该图像的平移和镜像(选作)。
?x'??10x0??x?
?'???x'?x?x0???y? y?01y即：???0???? ?1??001??1??y'?y?y0??????
参考程序如下：
function [I]=hmove(i,x0,y0);
%编写实现图像平移的函数hmove，平移量为x0,y0,平移前图像矩阵为i， [r,c]=size(i);
I(r+x0,c+y0)=0; %平移后图像矩阵为I for x=1:r;
for y=1:c; x1=x+x0; y1=y+y0;
I(x1,y1)=i(x,y); end; end; 主程序:
RGB=imread(‘hua.jpg’); subplot(2,2,1)
imshow(RGB),title(‘原图’) subplot(2,2,2)
gray1=rgb2gray(RGB);title(‘灰度图’) imagesc(gray1),colormap(gray); subplot(2,2,3)
I1=hmove(gray1,100,20);title(‘平移图’) subimage(gray1),axis('image');
subplot(2,2,4),imagesc(I1),colormap(gray),axis([1,700],[1,820]); 镜像
沿纵轴翻转，水平镜像：a(x,y) = -x; b(x,y) = y;
沿横轴翻转，垂直镜像：a(x,y) = x; b(x,y) = -y; 对于实际的图像来说,即将图像的左右或上下翻转， 程序如下：
I0=imread(‘hua.jpg’); I=rgb2gray(I0); [r,c]=size(I); I1=zeros(r,c); %水平镜像
I1(1:r,1:c)=I(1:r,c:1); %垂直镜像
I1(1:r,1:c)=I(r:1,1:c);
（5)、了解用imtransform和maketform函数实现各种空间变换的方法。 演示一下程序：
仿射变换，可以用以下函数来描述：
，其中，A是变形矩阵，b是平移矩阵。
尺度变换：
CLF;I=checkerboard(20,2);
subplot(121);imshow(I);axis on;title('原图') s=1.5;T=[s 0;0 s;0 0]; tf=maketform('affine',T);
I1=imtransform(I,tf,'bicubic','FillValues',0.3); subplot(122);imshow(I1);axis on;title('尺度变换') 伸缩变换 【例】
CLF;I=checkerboard(20,2);
subplot(121);imshow(I);axis on;title('原图') t=2;T=[1 0;0 t;0 0];
tf=maketform('affine',T);
I1=imtransform(I,tf,'bicubic','FillValues',0.3); subplot(122);imshow(I1);axis on;title('伸缩变换') 旋转变换：
CLF;I = checkerboard(20,2);
subplot(1,2,1);imshow(I);title('原图') angle=20*pi/180;
sc=cos(angle);ss=sin(angle); T=[sc -ss; ss sc;0 0]; tf=maketform('affine',T);
I1=imtransform(I,tf,'bicubic','FillValues',0.3); subplot(122);imshow(I1);title('旋转变换') 综合变换 【例】
CLF;I = checkerboard(20,2);
subplot(1,2,1);imshow(I);title('原图') s=2;As=[s 0;0 s]; % 尺度 t=2;At=[1 0;0 t]; % 伸缩 u=1.5;Au=[1 u;0 1]; % 扭曲
st=30*pi/180;sc=cos(angle);ss=sin(angle); Ast=[sc -ss; ss sc]; % 旋转 T=[As*At*Au*Ast;3 5]; tf=maketform('affine',T);
I1=imtransform(I,tf,'bicubic','FillValues',0.3); subplot(122);imshow(I1);title('综合')
四、图像邻域操作、区域操作、图像统计及MATLAB实现 1、运行以下程序,观察显示结果： （1）、%用函数mean作滑动处理 CLF
I=imread('tire.tif');
I2=nlfilter(I,[5 5],'mean2'); subplot(121),imshow(I,[]); subplot(122),imshow(I2,[]); （2）%用自定义函数作滑动处理 CLF
I=imread('tire.tif'); f=inline('max(x(:))'); I2=nlfilter(I,[3 3],f); subplot(1,2,1),imshow(I); subplot(1,2,2),imshow(I2); （3）%快速滑动领域操作 CLF
I=imread('tire.tif');
I2=colfilt(I,[5 5],'sliding','mean'); subplot(121),imshow(I,[]); subplot(122),imshow(I2,[]); （4）%用列操作函数实现滑动 CLF
I=imread('ic.tif');
I1=im2col(I,[3 3],'sliding');
I1=uint8([0 -1 0 -1 4 -1 0 -1 0]*double(I1)); I2=col2im(I1,[3,3],size(I),'sliding'); subplot(121),imshow(I,[]); subplot(122),imshow(I2,[]); （5）%图像块操作（分离） CLF
I=imread('tire.tif');
f=inline('mean2(x)*ones(size(x))'); I2=blkproc(I,[8 8],f);
subplot(1,2,1),imshow(I,[]) subplot(1,2,2),imshow(I2,[]) （6）%快速块操作（分离） CLF
I=imread('tire.tif');
f=inline('ones(64,1)* mean(x)'); I2=colfilt(I,[8 8],'distinct',f); subplot(1,2,1),imshow(I,[]) subplot(1,2,2),imshow(I2,[])
（7）%用列操作函数实现块操作（分离） CLF
I=imread('tire.tif');
I1=im2col(I,[8 8],'distinct'); I1=ones(64,1)* mean(I1);
I2=col2im(I1,[8,8],size(I),'distinct'); subplot(121),imshow(I,[]); subplot(122),imshow(I2,[]);
从显示结果中，试分析滑动领域与分离领域操作的区别？
答：滑动领域各邻域之间有重叠，而分离领域操作各领域间没有重叠部分，因此同样做均值处理，滑动邻域操作后，图像相对光滑，而分离领域操作后图像有块状马赛克效应，颗粒感强。
2、练习roipoly, roicolor, roifilter函数的使用（参考相关教学资料） 可演示以下部分程序来体会函数的用法; 【例】多边形选择区域 CLF
I=imread('eight.tif');
c=[222 272 300 270 221 194]; r=[21 21 75 121 121 75]; BW=roipoly(I,c,r);
subplot(121),subimage(I); subplot(122),subimage(BW); 【例】灰度选择法 CLF
I=imread('rice.tif'); BW=roicolor(I,128,255); subplot(121),subimage(I); subplot(122),subimage(BW); 【例】对指定区域进行锐化 CLF
I=imread('eight.tif');
c=[222 272 300 270 221 194]; r=[21 21 75 121 121 75]; BW=roipoly(I,c,r);
h=fspecial('unsharp'); % 滤波函数 J=roifilt2(h,I,BW);
subplot(121),subimage(I); subplot(122),subimage(J); I=imread('eight.tif');
c=[222 272 300 270 221 194]; r=[21 21 75 121 121 75]; J=roifill(I,c,r);
subplot(121),subimage(I); subplot(122),subimage(J);
【思考题】
1、 MatLab软件可以支持哪些图像文件格式？
答：Matlab可支持的图像格式包括BMP,JPEG,TIFF,GIF等等。
2、MATLAB中有哪些类型的图像，在MATLAB中分别用什么类型的数组来存储？ 答：MATLAB图像处理工具箱包括5种类型的图像，如下 二进制图像，使用uint8或双精度类型数组存储
索引图像，数据矩阵使uint8,uint16或双精度类型，颜色映射矩阵每个元素均为[0 1]间的双精度浮点型数据。
灰度图像，uint8,uint16或双精度类型的数组来描述 RGB图像，双精度浮点型，uint8,uint16类型 多帧图像
3、由图像算术运算的运算结果分析，思考图像减法运算在什么场合上发挥优势？ 答：图像减法是一种常用于检测图像变化及运动物体的图像处理方法。在很短时间内，可认为图像背景固定不变，可以直接运用减法运算来检测变化或运动的物体。 4、非比例缩放得到的图片和比例缩放得到的图片有何差异？是否失真？
答：非比例缩放图像对原图有畸变，而按照比例缩放的图像没有畸变。由于x,y方向的缩放比例不同，改变了原始图像像素的相对位置，产生了几何畸变。 5、图像的旋转会导致图像的失真么？
图像的旋转有时会发生失真，尤其旋转不是90度整数倍时。 【实验设备】
装有MATLAB软件的电脑，u盘
【实验报告要求】
1、实验报告写明实验题目，可简写实验原理，实验内容，详细写实验过程和实验结果分析，要求附相关程序及显示结果 2、回答思考题，写在实验报告上
实验3 MATLAB实现数字图像的变换技术（2学时）
【实验内容】 一、傅里叶变换
1、选择任意一幅图像进行平移和旋转，显示原始图像与处理后图像，分别对其进行FFT傅里叶变换，显示变换后结果，并要求使用FFTshift搬移频谱中心，分析原图的傅里叶谱与平移，旋转后傅里叶频谱的对应关系。
参考程序如下（参考程序中是人为构造一个图像矩阵，形成一幅图像，也可选择现成的图像进行平移，旋转，而后求解傅里叶变换）：
【例】图像的平移 CLF
f1 = zeros(256,256);f1(108:148,108:148) = 1;
f2 = zeros(256,256);f2(158:198,58:98) = 1;%与f1相比，就是平移效果
F1 = fft2(f1);F1 = log(1+abs(F1)); F2 = fft2(f2);F2 = log(1+abs(F2)); subplot(221),imshow(f1,[]) subplot(222),imshow(f2,[]) subplot(223),imshow(F1,[]) subplot(224),imshow(F2,[])
平移性：原图像平移，变换图像幅值图没有变化，但相位会旋转。
【例】图像的旋转 CLF
f1 = zeros(256,256);f1(58:198,108:148) = 1; f2 = imrotate(f1,45,'crop');
F1 = log(1+abs(fftshift(fft2(f1)))); F2 = log(1+abs(fftshift(fft2(f2)))); subplot(221),imshow(f1,[]) subplot(222),imshow(f2,[]) subplot(223),imshow(F1,[]) subplot(224),imshow(F2,[]) 图像旋转，谱旋转
注：
2、运行如下程序： bw = imread('text.png)
a = bw(32:45,88:98); %a=imcrop(bw); figure,subplot(1,2,1),imshow(bw); subplot(1,2,2), imshow(a);
C = real(ifft2(fft2(bw) .* fft2(rot90(a,2),256,256))); figure,subplot(1,2,1), imshow(C,[]) ; max(C(:));thresh = 60;
subplot(1,2,2), imshow(C > thresh);% 解释说明FFT变换图像特征识别的原理； 傅立叶变换可用来分析图像的相关性，相关性可确定一幅图像的特征，此意义下，相关性通常称为模板匹配。自己与自己的相关性最强，此时相关函数可达峰值。 而FFT变换可用于图像的相关函数计算，在数字图像处理相关的运算中用于图像匹配，可以对某些模板对应的特征进行定位。算法是：以待定的目标为模板在待识别图像上滑动，对运算接过取适当阈值，确定待定位目标的位置。
二、DCT变换
1、运行如下语句，分析显示结果： %应用MATLAB实现图像的DCT变换。
A=imread('pout.tif'); %读入图像
I=dct2(A); %对图像作DCT变换 subplot(1,2,1)
imshow(A); %显示原图像 subplot(1,2,2)
imshow(log(abs(I)),[0 5]);
显示图像DCT变换后的图像，发现DCT变换后数值较大的主要集中在一部分区域(左上角)，而非所有区域。
2、
运行如下语句，分析显示结果 üT变换重构图像
RGB=imread('autumn.tif'); figure(1),imshow(RGB); I=rgb2gray(RGB); figure(2),imshow(I);
J=dct2(I);figure(3),imshow(log(abs(J)),[]),colormap(jet(64)),colorbar;
J(abs(J)<10)=0; K=idct2(J)/255;
figure(4),imshow(K);
对原图像DCT变换后，保留大于某阈值的DCT系数值，其余为0，利用idct重
构图像，与原图像画质相比，损失很少，这为图像压缩提供了思路。
3、运行如下程序，分析运行结果： D= imread('D:\\lena.jpg'); D1=rgb2gray(D); s=size(D1); D2=dct2(D1); D3=idct2(D2);
P=zeros(s);P1=P;P2=P;P3=P; P1(1:10,1:10)=D2(1:10,1:10); P2(1:30,1:30)=D2(1:30,1:30); P3(1:60,1:60)=D2(1:60,1:60); E1=idct2(P1); E2=idct2(P2); E3=idct2(P3);
subplot(2,3,1); imshow(D1); subplot(2,3,2); imshow(D2); subplot(2,3,3);image(D3)
subplot(2,3,4); image(E1); subplot(2,3,5); image(E2); subplot(2,3,6); image(E3)
通过以上实例，分析DCT变换应用于图像恢复与压缩的原理。
对原图像DCT变换后，分别保留DCT系数矩阵大小分别为10*10，20*20,30*30的值，其余为0，而后利用idct变换进行重构图像，与原图像画质相比，并比较各重构图像之后的效果图，可以发现只要保留最主要的一些系数，无需全部DCT系数来重构，图像就能达到较好的重构效果，这为图像压缩提供了思路。
三、练习并简单了解Radon变换 利用demos演示完成
运行demos,而后查找radon变换,查询相关的使用实例
【思考题】
1、FFT变换具有哪些性质?应用于数字图像处理后，如何体现出来 答：平移性：原图像平移，变换图像幅值图没有变化。
旋转不变性：如果原函数旋转一个角度，则变换函数同样旋转相同的角
比例性质（尺度变换特性）：原图像的缩放，影响到变换图像的缩放。时域中信号被压缩，频域中信号就被拉伸。 2、FFT为什么能用于图像压缩？
答：FFT变换系数刚好表现的是各个频率点上的幅值。在小波变换没有提出时，用来进行压缩编码。考虑到高频反映细节、低频反映景物概貌的特性。往往认为可将高频系数置为0，骗过人眼。这样只需主要考虑保留低频系数，达到了压缩的目的。
3、解释FFTshift函数的功能？
答：傅里叶变换以零点为中心，导致谱图像最亮点在图像的左上角。为符合正常习惯，将谱图像的原零点从左上角移到显示屏的中心。函数fftshift(I)
可将变换后图像频谱中心从矩阵原点移到矩阵中心。即矩阵I一、三象限和二、四象限进行互换。而图像的四个象限是指以矩形图像的纵向以及对横向对称轴所分割的四个区域。
4、图像还有哪些正交变换，为什么正交变换可以用于图像压缩？ 答；图像还有DCT,沃尔什，哈达玛变换，K_L，小波变换等，
信号经过正交变换，分解成相互独立的若干分量，去除分量中的某一部分分量对于其他分量并无影响。于是，在可以容忍的前提下，选择正交分量中非本质的部分去除掉，留下信号的本质分量。整个信号量减小，即被压缩。
5、简单叙述Radon变换用于直线检测的原理？
答：对原始图像进行radon变换，radon变换矩阵中的峰值，即原始图像中的直线。因为存在直线时，该直线的点经过radon变换后沿着某个方向（角度）有多次投影，存在累计最大值。
【实验报告要求】
1、实验报告写明实验题目，可简写实验原理，实验内容，详细写实验过程和实验结果分析，要求附相关程序。
2、回答思考题，写在实验报告上
实验4 MATLAB实现数字图像增强（2学时）
【实验内容】 1、灰度变换
(1)、利用直接灰度变换法对任选的一副图像进行灰度变换（imjust函数），要求将原图像0.5到0.75的灰度级扩展到范围[0 1]，再实现明暗转换(负片图像)，显示原图，直方图，处理后的图像。 参考程序如下：
I=imread('flower.tif'); Y=rgb2gray(I); M=255-Y;%负片
J=imadjust(Y,[0.5,0.75],[0,1]);% 将原图像0.5到0.75的灰度级扩展到范围[0 1]
subplot(3,3,1),imshow(Y); subplot(3,3,2),imshow(J); subplot(3,3,3),imshow(M); subplot(3,3,4),imhist(Y); subplot(3,3,5),imhist(J); subplot(3,3,6),imhist(M);
显示结果：
(2）、一幅灰度图像如图4.4左(也可选用任意一副灰度图像),运用imadjust函数，将0~60灰度级压缩到0~30范围内,压缩比1/2;60~180的灰度级扩大到30~240,比率为190/120;将180~255灰度级压缩到240~255范围内,压缩比为15/75。效果图如图右。
图4.4 线性灰度变换例图
参考程序如下：
I1=imread(‘fruit.jpg’);
I2=imadjust(I1,[0 ,60]/255,[0,30]/255,0.5,[60,180]/255,[30,240]/255,19/12,[180,255]/255,[240,255]/255,15/75);
subplot(2,1,1),imshow(I1,[]),title(‘原图’);
subplot(2,1,2),imshow(I2,[]),title(‘灰度变换后图’) (3）、对数变换
对数变换常用来扩展低值灰度，压缩高值灰度，这样可以使低值灰度的图像细节更容易看清楚。对数变换的表达式为： g(x,y)=log[f(x,y)+1]
运行如下程序，分析显示结果 I=imread('pout.png'); subplot(121),imshow(I); I=double(I);J=log(I+1); subplot(122),imshow(J,[]);
可以看出，原图像的灰度有所变化。 (4）、任选一副灰度图像对其进行幂次变换，选择不同的r(分别取r=0.5,r= 1,r=5),观察分析其显示结果 参考程序如下： clear all; close all;
I0=imread(‘test.jpg’); I0=I0/255;
I1=I0.^0.5;%r=0.5幂次变换 I2=I0.^5; %r=5幂次变换
subplot(1,3,1),imshow(I0,[]),title(‘原图’)
subplot(1,3,2), imshow(I0,[]),title(‘r=0.5幂次变换’); subplot(1,3,2), imshow(I0,[]),title(‘r=5幂次变换’);
从显示结果可以看出来，r<1时，黑色区域被扩展，变得清晰，当r<1,黑色区域被压缩，变得几乎不可见。
(5)、练习教学资料上分段线性灰度变换例子。
可练习的参考程序如下，该程序可实现灰度变换(2)题的操作：
x1=imread('pout.tif');
subplot(221),imshow(X1);title('原图') f0=0;g0=0; f1=60;g1=30; f2=180;g2=240; f3=255;g3=255; % 绘制变换曲线
subplot(222),plot([f0 f1 f2 f3],[g0 g1 g2 g3]);
axis tight,xlabel('f'),ylabel('g'),title('灰度变换曲线') 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:m for J=1:n f=X2(I,J); g(I,J)=0;
if (f>=0)&(f<=f1) g(I,J)=r1*f+b1; elseif (f>f1)&(f<=f2) g(I,J)=r2*f+b2; elseif (f>f2)&(f<=f3) g(I,J)=r3*f+b3; end end end
subplot(223),imshow(mat2gray(g));title('灰度变换图')
显示结果如下：
2、直方图变换（直方图均衡化和直方图规定化）实现图像增强(histeq，imhist函数) （1）、任选一幅灰度图像，显示原图像，绘制其及直方图，然后对图像进行均衡化处理，并显示处理后图像及其直方图，与原图像作比较； 参考程序如下：
I=imread('Lenna.tif'); M=histeq(I);
J0=histeq(I,32);
subplot(2,4,1),imshow(I);title(‘原图’)；
subplot(2,4,2),imhist(I); title(‘原图直方图’)；
subplot(2,4,3),imshow(M); title(‘直方图均衡化图像’)；
subplot(2,4,4),imhist(M); title(‘直方图均衡化后的直方图’)；
subplot(2,4,5),imshow(J0); title(‘32个灰度级直方图均衡化图像’)； subplot(2,4,6),imhist(J0); title(‘32个灰度级直方图均衡化的直方图’)；
[counts,x]=imhist(J); J1=histeq(Q,counts); subplot(2,4,7),imshow(J1); title(‘32个灰度级直方图规定化’)；
subplot(2,4,8),imhist(J1); title(‘32个灰度级直方图规定化的直方图’)； %可将J0%与J1进行比较
显示结果如下：
（2）、对原图做直方图规定化，目标直方图服从正态分布----N(4,8),观察经过规定化处理的图像效果。 clear all; close all;
I0=rgb2gray(imread('lena.jpg')); %构造正态分布直方图
N=32;2个灰度级 i=0:N-1;
Hist0=exp(-(i-1).^2/8); Hist0=Hist0/sum(Hist0); %Hist_cu=cumsum(Hist0); I1=histeq(I0,Hist0);
subplot(3,2,1),imshow(I0,[]),title('原图'); subplot(3,2,2),imhist(I0),title('原图直方图')
subplot(3,2,3),stem([0:31],Hist0),title('正态分布目标直方图');%正态分布直方图
subplot(3,2,5),imshow(I1,[]),title('直方图规定化后的图像'); subplot(3,2,6),imhist(I1),title('直方图规定化后的图像直方图');
图像做过均衡化处理后比原图的对比度增强了，图像的细节更加清楚了，均衡化后的直方图比原直方图分布均匀了；直方图规定化是以直方图均衡化做基础，可以按人的意愿去变换，选择合适的目标直方图来实现规定化，可以突出感兴趣的灰度范围，使用起来更灵活。
3、空间域滤波
（1）、选择任一副图像作为源图像，使用函数imnoise分别添加两种不同的噪声，用ordfilt2实现图像增强，然后设置不同的参数，比较显示结果有什么不同。 参考考程序如下：
I = imread('moon.tif');
J = imnoise(I,'salt & pepper',0.04); G = imnoise(I,'gaussian',0.04);
K1 = ordfilt2(J,1,ones(3,3));%最小滤波 K2 = ordfilt2(J,9,ones(3,3));%最大滤波 K3 = ordfilt2(J,5,ones(3,3));%中值滤波 K4 = ordfilt2(G,1,ones(3,3)); K5 = ordfilt2(G,9,ones(3,3)); K6 = ordfilt2(G,5,ones(3,3));
subplot(3,3,1);imshow(I);title(‘原图’)
subplot(3,3,2);imshow(J); title(‘添加椒盐噪声后’) subplot(3,3,3);imshow(G); title(‘添加高斯噪声后’) subplot(3,3,4);imshow(K1); title(‘最小滤波椒盐噪’) subplot(3,3,5);imshow(K2); title(‘最大滤波椒盐噪’) subplot(3,3,6);imshow(K3); title(‘中值滤波椒盐噪’) subplot(3,3,7);imshow(K4); title(‘最小滤波高斯噪’) subplot(3,3,8);imshow(K5); title(‘最大滤波高斯噪’) subplot(3,3,9);imshow(K6); title(‘中值滤波高斯噪’)
（2）、请采用二维中值滤波函数medfilt2对受椒盐噪声干扰的任意一副图像滤波，窗口分别采用3*3,5*5，显示结果(提示：首先对图像添加椒盐噪声，最后调用mefilter2函数滤波) 参考程序如下：
I = imread('moon.tif');
J = imnoise(I,'salt & pepper',0.04);
K1 = medfilt2(J,[3 3]); K2= medfilt2(J,[5 5]);
subplot(2,2,1);imshow(I);title(‘原图’)
subplot(2,2,2);imshow(J);title(‘添加椒盐噪声后的图像’) subplot(2,2,3);imshow(K1);title(‘3*3中值滤波’) subplot(2,2,4);imshow(K2);title(‘5*5中值滤波’)
(3）、选择合适的MATLAB函数对受高斯噪声干扰的一副图像进行均值滤波（提示：首先对图像添加高斯噪声，然后用fspecial函数构造均值滤波器h，最后调用filter2函数滤波）； 参考程序如下：
I=imread('Lenna.tif');
K=imnoise(I,'gaussian',0.04); h1=fspecial('average',[3,3]); h2=fspecial('average',[5,5]); J=filter2(h1,I); M=filter2(h2,I);
subplot(2,2,1),imshow(I,[]);title(‘原图’)
subplot(2,2,2),imshow(K,[]);title(‘加高斯噪声后’)
subplot(2,2,3),imshow(J,[]);title(‘3*3均值滤波高斯噪声后’) subplot(2,2,4),imshow(M,[]); title(‘5*5均值滤波高斯噪声后’) (4）、任意采用三种不同锐化算子对任一副图像进行锐化处理，显示原图像及处理后的图像效果。(可以用fpecial,也可使用直接写出算子) 参考程序如下：
I=imread('Lenna.tif'); h1=fspecial('unsharp',0);
I1=filter2(h1,I);h2=fspecial('sobel');%h2=[]可直接写sobel算子矩阵 I2=filter2(h2,I);
h3=fspecial('prewitt');%h3=可直接写算子矩阵 I3=filter2(h3,I);
subplot(221),imshow(I,[]); subplot(222),imshow(I1,[]); subplot(223),imshow(I2,[]); subplot(224),imshow(I3,[]);
（5）、选择任意幅图像作为源图像，练习使用fspecial和filter2函数实现各种算子如：拉普拉斯，sobel，prewitt等完成图像锐化，比较各算子的操作效果有何不同。
参考程序如下：注:各算子也可直接写出算子的形式
I=imread('coins.png');
h1=fspecial('unsharp',0); % 二维拉普拉斯滤器h1=[-1 -1 -1;-1 9 -1;-1 -1 -1]; I1=filter2(h1,I); h2=fspecial('sobel'); I2=filter2(h2,I);
h3=fspecial('prewitt'); I3=filter2(h3,I);
subplot(221),imshow(I,[]); subplot(222),imshow(I1,[]);
subplot(223),imshow(I2,[]); subplot(224),imshow(I3,[]);
新疆大学信息科学与工程学院
这三种锐化算子都使边缘和轮廓线模糊的图像变得清晰，边界更加明显
4、频域滤波实现图像增强
分别调试运行以下程序，修改可能的错误，比较分析实验结果，并回答问题 1）、低通滤波， A、高斯低通滤波
I=double(imread(‘lena.jpg’));
H=fspecial(‘gaussian’,7,3);%构造什么滤波器？参数？ FD=double(filter2(H,I));
Figure,subplot(1,2,1),imshow(I); subplot(1,2,2),imshow(FD);
B、%各种频率低通滤波器对gray.bmp图像的增强的MATLAB程序清单： clc;
[I,map]=imread(‘gray.bmp’);
noisy= imnoise(I,’gaussian’,0.01); Imshow(noisy,map); [M N]=size(I); F=fft2(noisy); fftshift(F);
%以下程序段功能是什么 Dcut=100; DO=150; D1=250; for u=1 : M
for v=1 : N
D( u,v)=sqrt(u^2+v^2);
BUTTERH(u,v)=1/(1+(sqrt(2)-1)*(D(u,v)/ Duct)^2); XEPOTH（u,v)=exp(log(1/(sqrt(2))*(D(u,v)/ Duct)^2); if D(u,v) TRAPEH(u,v) =1; elseif D(u,v) <=D1 TRAPEH(u,v)= D(u,v) –D1)/(D0-D1); else TRAPEH(u,v)=0; end end end BUTTERG=BUTTERH.*F; BUTTERfiltered=ifft2(BUTTERG); EXPOTG=EXPOTH.*F; EXPOTfiltered=ifft2(EXPOTG); TRAPEG=TRAPEG.*F; TRAPEfiltered=ifft2(TRAPEG); 31 新疆大学信息科学与工程学院 subplot(2,2,1),imshow(noisy); subplot(2,2,2),imshow(BUTTERfilteered) subplot(2,2,3),imshow(EXPOTfiltered) subplot(2,2,4),imshow(TRAPEfiltered) 2）、频域高通滤波法对图像进行增强 A、频域高通滤波程序 I1=imread('blood1.tif'); [M N]=size(I1); F=fft2(I1);F=fftshift(F); m1=fix(M/2);n1=fix(N/2); for u=1:M for v=1:N D=sqrt((u-m1)^2+(v-n1)^2); if D==0 H(u,v)=0; else H(u,v)=1/(1+0.414*(5/D)^4); end end end F1=H.*F;F1=ifftshift(F1); I2=abs(ifft2(F1)); subplot(121),imshow(I1,[]); subplot(122),imshow(I2,[]); %什么类型高通滤波器？ B、选择各种频域高通滤波器对图像进行增强 clc; [I,map]=imread(‘gray.bmp’); noisy= imnoise(I,’gaussian’,0.01); [M N]=size(I); F=fft2(noisy); Fftshift(F); Dcut=100; DO=250; D1=150; for u=1 : M for v=1 : N D( u,v)=sqrt(u^2+v^2); BUTTERH(u,v)=1/(1+(sqrt(2)-1)*(Dcut /D(u,v))^2); XEPOTH（u,v)=exp(log(1/(sqrt(2))*(Dcut /D(u,v))^2); if D(u,v) TRAPEH(u,v) =0; 32 新疆大学信息科学与工程学院 elseif D(u,v) <=D0 TRAPEH(u,v)= D(u,v) –D1)/(D0-D1); else TRAPEH(u,v)=1; end end end BUTTERG=BUTTERH.*F; BUTTERfiltered=ifft2(BUTTERG); EXPOTG=EXPOTH.*F; EXPOTfiltered=ifft2(EXPOTG); TRAPFG=THPFH.*F; TRAPEfiltered=ifft2(THPFG); subplot(2,2,1),imshow(noisy); subplot(2,2,2),imshow(BUTTERfilteered) subplot(2,2,3),imshow(EXPOTfiltered) subplot(2,2,4),imshow(TRAPEfiltered) 【思考题】 1、该实验可应用到哪些实际问题中？ 答：图像平滑的主要目的是减少图像噪声，减少噪声的方法可以在空间域或频率域进行，空间域主要用中值滤波或均值滤波，而频域则主要用低通滤波技术。 而图像锐化的目的则是使边缘和轮廓线模糊的图像变得清晰，图像锐化可以在空间域中利用对图像微分进行，也可以在频域中运用高通滤波技术处理。 应用：印刷中的细微层次强调。弥补扫描、挂网对图像的钝化 超声探测成象，分辨率低，边缘模糊，通过锐化来改善 图像识别中，分割前的边缘提取 锐化处理恢复过度钝化、暴光不足的图像 图像创艺（只剩下边界的特殊图像） 尖端武器的目标识别、定位 2、利用中值滤波与均值滤波分别对图像处理，两者原理与处理效果有何不同？ 答:均值滤波：这种方法的基本思想是用几个像素灰度的平均来代替一个像素原来的灰度值，实现图像的平滑。中值滤波是用一个有奇数点的滑动窗口，将窗口中心点的值用窗口各点的中值代替。 相比之下： (1) 对大的边缘高度，中值滤波较邻域均值好得多，而对于较小边缘高度，两种滤波只有很少差别。 （2）中值滤波是非线性的。 （3）中值滤波在抑制图像随机脉冲噪声方面甚为有效。且运算速度快，便于实时处理。 （4）中值滤波去除孤立线或点干扰，而保留空间清晰度较平滑滤波为好； 33 新疆大学信息科学与工程学院 但对高斯噪声则不如平滑滤波。 3、图像锐化有哪些哪些算子？请写出相应的算子矩阵 答：图像锐化的算子有如下： Roberts算子模板： Sobel算子模板： Prewitt算子模板： Laplacian算子模板： LOG算子： 4、低通和高通滤波分别对数字图像起到什么效果？ 答：低通滤波对图像起平滑作用，高通滤波对图像起锐化作用 5、低通的理想滤波器，巴特沃斯滤波器，指数滤波器平滑滤波在平滑数字图像有何区别？ 答：既然H(u,v)是理想的矩形特性，那么它的反变换h(x,y)的特性必然会产生无限的振铃特性。经与f(x,y)卷积后则给g(x,y)带来模糊和振铃现象，D0越小这种现象越严重，当然，其平滑效果也就较差。这是理想低通不可克服的弱点。 振铃效果——理想低通滤波器的一种特性 。 Butterworth低通过滤器的分析 在任何经BLPF处理过的图像中都没有明显的振铃效果，这是过滤器在低频和高频之间的平滑过渡的结果 低通滤波是一个以牺牲图像清晰度为代价来减少干扰效果的修饰过程 34 新疆大学信息科学与工程学院 与理想低通滤波器的处理结果相比，经Butterworth滤波器处理过的图像模糊程度会大大减少。因为它的H(u,v)不是陡峭的截止特性，它的尾部会包含有大量的高频成份。 BLPF处理过的图像中都没有振铃效果.这是由于在滤波器的通带和阻带之间有一平滑过渡的缘故。另外，由于图像信号本身的特性，在卷积过程中的折迭误差也可以忽略掉。由此可知，Butterworth 低通滤波器的处理结果比理想滤波器为好。 由于指教低通滤波器有更快的衰减率，所以，经指数低通滤波的图像比Butterworth低通滤波器处理的图像稍模糊一些。由于指数低通滤波器的传递函数也有较平滑的过渡带，所以图像中也没有振铃现象。 由于梯形滤波器的传递函数特性介于理想低通滤波器和具有平滑过渡带滤波器之间，所以其处理效果也介于其两者中间。梯形滤波法的结果有一定的振铃现象。 用低通滤波器进行平滑处理可以使噪声伪轮廓等寄生效应减低到不显眼的程度，但是由于低通滤波器对噪声等寄生成份滤除的同时，对有用高频成份也滤除，因此，这种去噪声的美化处理是以牺牲清晰度为代价而换取的。 【实验报告要求】 1、实验报告写明实验题目，可简写实验原理，实验内容，详细写实验过程和实验结果分析，要求附相关程序。 2、回答思考题，写在实验报告上 35 新疆大学信息科学与工程学院 实验5 MATLAB实现图像的复原（2学时） 【实验内容】 1、任选一幅彩色风景图片作为源图像，用imadd和imnoise给图像添加不同类型的噪声，显示噪声图像。 参考程序如下： a=imread('peppers.png'); b=imnoise(a,'gaussian',0.05);%高斯噪声 c=0.1*randn(size(a));%随机噪声 d=imadd(a,im2uint8(c)); a=a(10+[1:256],222+[1:256],:); subplot(1,3,1),imshow(a);title(‘原图’) subplot(1,3,2),imshow(b); subplot(1,3,3),imshow(d); 2、设置不同的模糊参数实现任一副图像的运动模糊(fspecial,imfilter函数)，显示结果。 a=imread('peppers.png'); a=a(10+[1:256],222+[1:256],:); LEN=50; THETA=39; PSF=fspecial('motion',LEN,THETA); e=imfilter(a,PSF,'circular','conv'); LEN=100; THETA=39; PSF=fspecial('motion',LEN,THETA); f=imfilter(a,PSF,'circular','conv'); LEN=20; THETA=39; PSF=fspecial('motion',LEN,THETA); m=imfilter(a,PSF,'circular','conv'); LEN=50; THETA=10; 36 新疆大学信息科学与工程学院 PSF=fspecial('motion',LEN,THETA); E=imfilter(a,PSF,'circular','conv'); LEN=100; THETA=10; PSF=fspecial('motion',LEN,THETA); F=imfilter(a,PSF,'circular','conv'); LEN=20; THETA=10; PSF=fspecial('motion',LEN,THETA); M=imfilter(a,PSF,'circular','conv'); subplot(2,4,1),imshow(a); subplot(2,4,2),imshow(e); subplot(2,4,3),imshow(f); subplot(2,4,4),imshow(m); subplot(2,4,5),imshow(E); subplot(2,4,6),imshow(F); subplot(2,4,7),imshow(M); 3、对1,2产生的图像分别进行复原，选用逆滤波法，改进逆滤波器法进行图像复原.（需要编程，为选作题），显示处理结果。 4、对1,2产生的图像分别进行复原，选用维纳滤波器进行图像复原，显示处理结果。 参考程序如下： 在1,2产生的图像后可加一段： figure, wnr1=deconvwnr(b,PSF); subplot(2,3,1),imshow(wnr1); wnr2=deconvwnr(d,PSF); subplot(2,3,2),imshow(wnr2); wnr3=deconvwnr(e,PSF); subplot(2,3,3),imshow(wnr3); wnr4=deconvwnr(f,PSF); subplot(2,3,4),imshow(wnr4); wnr5=deconvwnr(m,PSF); subplot(2,3,5),imshow(wnr5); 5、对1,2产生的图像使用维纳滤波图像复原(deconvwnr)、约束最小平方滤波图像复原(deconvreg)2种方法任一种去实现，设置不同参数，比较复原效果。（选作，提示：可参看相关参考书实现deconvreg函数功能） 【思考题】 1、逆滤波器存在什么病态问题，如何改进？ 答：复原由退化函数H退化的图像最直接的方法是直接逆滤波。在方法中，用退化函数除退化图像的傅立叶变换来计算原始图像的傅立叶变换。 37 新疆大学信息科学与工程学院 即使我们知道退化函数，也可能无法准确复原未退化的图像。因为噪声是一个随机函数，其傅氏变换未知。实际应用逆滤波复原方法时存在病态的问题，即在频谱域中对图像信号的那些频谱点上，若：H(u,v)=0,而噪声N(u,v)<>0,则 N(u,v)/H(u,v)>>F(u,v),面目全非，这种恢复会出现病态性。由上述方法进行图像复原时，由于H(u,v)在分母上，当u-v平面上的某引起点或区域H(u, v)很小或等于零，即出现了零点时， 就会导致不稳定解。因此，即使没有噪声，一般也不可能精确地复原f(x, y)。如果考虑噪声项N(x, y)， 则出现零点时， 噪声项将被放大，零点的影响将会更大，对复原的结果起主导地位， 这就是无约束图像复原模型的病态性质。它意味着退化图像中小的噪声干扰在H(u, v)取得很小值的那些频谱上将对恢复图像产生很大的影响。由简单的光学分析知道，在超出光学系统的绕射极限时，H(u, v)将很小或等于零，因此对多数图像直接采用逆滤波复原会遇到上述求解方程的病态性。 改进： 即如果退化为零或非常小的值，则N(u,v)/H(u,v)之比很容易决定复原函数的值。 1)、在H(u,v)=0处，人为仔细设置H-1(u,v)值，使得这些频谱点附近， N(u,v)/H(u,v)对^F(u,v)影响小些 2)、窄带H(u,v),相比噪声具低通特性。 2、维纳滤波适用于什么场合？有何缺陷？ 答：逆滤波较简单，但对噪声有放大作用，而维纳滤波则抑制噪声，维纳滤波即寻找一个滤波器，使得复原后的图像与原始图像的方差最小，即 1)、适用于假设图像和噪声不相关，且h(x,y)有零均值(平稳随机过程)，退化系统线性时不变 2)、信噪比较高时，等同逆滤波 3)、大多数情况下较满意，但信噪比较低时，效果不好，由于以平稳随机过程为基础，线性时不变条件不满足，与实际情况有差距；另外，最小均方误差准则与人眼的视觉准则并不匹配。 【实验报告】 1、对上述操作，要求编程，并附图像的显示结果。 2、实验报告写明实验题目，可简写实验原理，实验内容，详细写实验过程和实验结果分析，要求附相关程序。按格式书写，将结果贴入相应位置 3、回答思考题，写在实验报告上 38 新疆大学信息科学与工程学院 实验6 MATLAB实现图像分割（2学时） 【实验内容】 一、图像分割 1、任选一副图像，对实验图像用edge函数实现roberts算子、sobel、log等算子检测边缘，并对比显示其结果，评价一下各算子对于噪声条件下边界检测的性能。 参考程序如下： I=imread('rice.png'); K=imnoise(I,'salt & pepper',0.04); BW1=edge(I,'roberts'); BW2=edge(I,'sobel'); BW3=edge(I,'log'); BW4=edge(K,'roberts'); BW5=edge(K,'sobel'); BW6=edge(K,'log'); subplot(3,3,1),imshow(I),title('?ê?í???'); subplot(3,3,2),imshow(K),title('?ó?·????éùí?'); subplot(3,3,3),imshow(BW1),title('?í?roberts±??μ?ì2a'); subplot(3,3,4),imshow(BW2),title('?í?sobel±??μ?ì2a'); subplot(3,3,5),imshow(BW3),title('?í?log±??μ?ì2a'); subplot(3,3,6),imshow(BW4),title('?ó??roberts±??μ?ì2a'); subplot(3,3,7),imshow(BW5),title('?ó??sobel±??μ?ì2a'); subplot(3,3,8),imshow(BW6),title('?ó??log±??μ?ì2a'); 参考程序如下: I=imread('eight.tif'); imshow(I) BW1=edge(I,'roberts'); figure ,imshow(BW1),title('用Roberts算子') BW2=edge(I,'sobel'); figure,imshow(BW2),title('用Sobel算子 ') BW3=edge(I,'log'); figure,imshow(BW3),title('用拉普拉斯高斯算子') 得到： 39 新疆大学信息科学与工程学院 比较提取边缘的效果可以看出，sober算子是一种微分算子，对边缘的定位较精确，但是会漏去一些边缘细节。而Laplacian-Gaussian算子是一种二阶边缘检测方法，它通过寻找图象灰度值中二阶过零点来检测边缘并将边缘提取出来，边缘的细节比较丰富。通过比较可以看出Laplacian-Gaussian算子比sober算子边缘更完整，效果更好。 2、任选一副背景和前景区别较大的灰度图像，对该图像进行二值化处理(im2Bw，设置不同阈值实现)，绘制该图像直方图，利用双峰法和迭代法估计设置其阈值，然后利用计算出的阈值对测试图像进行分割处理，结果显示原图像、处理后的分割图像及相应的分割阈值，观察不同阈值的分割效果。 参考程序如下： clear; close all; I0=imread(‘cell.tif’); I1=im2bw(I0，0.3);%二值化处理 I1=im2bw(I0，0.7);%二值化处理 figure， subplot(3,2,1),imshow(I0),title(‘原图’); subplot(3,2,2),imhist(I0),title(‘原图直方图’); subplot(3,2,3),imshow(I1),title(‘二值化图,T=0.3’); subplot(3,2,4),imhist(I1),title(‘二值化图的直方图,T=0.3’); 40 新疆大学信息科学与工程学院 subplot(3,2,3),imshow(I2),title(‘二值化图,T=0.7’); subplot(3,2,4),imhist(I2),title(‘二值化图的直方图,T=0.7’); %简单的阈值分割法,利用函数产生全局阈值 T=graythresh(f); GT=I0>=T; figure, subplot(1,2,1), imshow(I0);title(‘原图’) subplot(1,2,2), imshow(GT);title(‘自动阈值分割效果’); %利用双峰法 i=imread('w01.tif'); subplot(1,2,1); imhist(i); title('原始图像直方图'); thread=130/255;%先观察直方图，再定阈值 subplot(1,2,2); i3=im2bw(i,thread); imshow(i3); title('分割结果'); 41 新疆大学信息科学与工程学院 根据原图像的直方图，发现背景和目标的分割值大约在130左右，并将灰度图像转为二值图像，分割效果比较理想。 %迭代阈值法进行分割 I0=imread(‘eight.tif’); I1=I0; I0=double(I0); T0=(min(I0(:))+max(I0(:)))/2; done=false; i=0; while ~done r1=find(I<=T0); r2=find(I>=T0); Tnew=(mean(I0(r1))+mean(I0(r2)))/2; done=abs(Tnew-T0)<1; T0=Tnew; i=i+1; end T0 I1(r1)=0; I1(r2)=1; subplot(1,2,1),imshow(I0),title(‘原图’); subplot(1,2,2),imshow(I1),title(‘迭代阈值法分割后图’) 3、对一幅典型的二值图像进行腐蚀和膨胀(erode和dilate)，开闭运算，使用方形模板3*3和5*5，显示结果 I=imread('circles.png');%假设原图为二值图像 %若原图为灰度图像，则先需二值化处理，才能腐蚀或膨胀 SE1=eye(5);%或SE1=strel(‘square’,5); BW1=IMerode(I,SE1); BW2=IMdilate(I,SE1); BW3=imopen(I,SE1); BW4=imclose(I,SE1); subplot(2,3,1),imshow(I),title('原图'); subplot(2,3,2),imshow(BW1),title('5*5腐蚀图'); subplot(2,3,3),imshow(BW2),title('5*5膨胀图'); subplot(2,3,4),imshow(BW3),title('5*5开运算图'); subplot(2,3,5),imshow(BW4),title('5*5闭运算图'); figure， SE2=eye(3); 或SE1=strel(‘square’,3); BW1=IMerode(I,SE2); BW2=IMdilate(I,SE2); BW3=imopen(I,SE2); BW4=imclose(I,SE2); subplot(2,3,1),imshow(I),title('原图'); subplot(2,3,2),imshow(BW1),title('3*3腐蚀图'); subplot(2,3,3),imshow(BW2),title('3*3膨胀图'); subplot(2,3,4),imshow(BW3),title('3*3开运算图'); 42 新疆大学信息科学与工程学院 subplot(2,3,5),imshow(BW4),title('3*3闭运算图'); 4、区域分割方法，分别利用区域生长法和区域分裂合并法对图像进行分割处理。 %区域生长法 参考程序如下： I0=imread(‘peppers_gray.bmp’); subplot(1,2,1),imshow(I0),title(‘原图’); seedx=[256,128,480]; seedy=[128,256,384]; hold on; plot(seedx,seedy,’gs’,’linewidth’,1); title(‘原始图像及种子位置’); I=double(I0); markerim=I=I(seedy(1),seedx(1)); for i=2:length(seedx) markerim=markerim|(I==I(seedy(i),seedx(i))); end thresh=[15,10,15]; maskim=zeros(size(I)); for i=1；length(seedx) g=abs(I-I(seedy(i),seedy(i)))<=thresh(i); maskim=maskim|g; end [g,nr]=bwlabel(imreconstruct(markerim,maskim),8); g=matgray(g); subplot(1,2,2),imshow(g);title(‘三个种子点区域生长结果’)； 区域分割法： function quadtree(x) I=imread(x); q=2^nextpow2(max(size(I))); [m n]=size(I); I=padaarray(I,[q-m,q-n],’post’); mindim=2; s=qtdecomp(I,@split,mindim,@predicate); Lmax=full(max(s(:))); g=zeros(size(I)); 43 新疆大学信息科学与工程学院 marker=zeros(size(I)); for k=1:lmax [vals,r,o]=qtgetblk(I,s,k); if ~isempty(vals) for i=1:length(r) xlow=r(i);ylow=c(i); xhigh=xlow+k-1; yhigh=ylow+k-1; region=I(xlow:xhigh,ylow:yhigh)=1; end end end end g=bwlabel(imreconstruct(marker,g)); g=g(1:m,1:n); I=I(1:m,1:n); subplot(1,2,1),imshow(I),title(‘原始图像’); subplot(1,2,2),imshow(g),title(‘分割后图像’) end function v=split(b,mindim,fun) k=size(b,3); v(1:k)=false; for i=1:k quadrgn=b(:,:,i); if size(quadrgn,1)<=mindim v(i)=false continue; end flag=feval(fun,quadrgn); if flag v(i)=true; end end 44 新疆大学信息科学与工程学院 end function flag=predicate(region) sd=std2(region); m=mean2(region); flag=(sd>5)&(m>0)&(m<200); end 5、 图6.5 例图 图6.6 例图 图6.7 例图 1）.运行如下程序实现对图6.5，实现去除图像中的噪声，请分析显示结果。 I=imread('test1.tif'); J=im2bw(I); se = strel('diamond',2); K=imerode(J,se); 45 新疆大学信息科学与工程学院 subplot(131),imshow(I) subplot(132),imshow(J) subplot(133), imshow(K) 2）.运行如下程序，实现将图6.6转化为二值图像，并计算图中鸡块中骨头的比重，请分析显示结果。 I=imread('test2.tif'); J=im2bw(I); total1=bwarea(J); a=size(J); s1=a(1),s2=a(2); s=s1*s2 k=total1/s 46 新疆大学信息科学与工程学院 3）.运行如下程序，实现去除图6.7中的矩形区域外的噪声，并填充矩形区域内部的小孔，请显示结果。 I=imread('test3.tif'); se = strel('diamond',5); J=imerode(I,se); se = strel('diamond',10); K=imdilate(J,se); subplot(131),imshow(I) subplot(132),imshow(J) subplot(133), imshow(K) 二、图像描述(根据课程内容，可选作其中部分内容) 查阅相应参考文献，完成 1、选择一幅彩色图像，编程计算其在RGB和HIS空间下的颜色矩，绘制颜色直方图 参考程序如下： %颜色矩，假设原图为rgb图像 I=imread('a.jpg'); I1=rgb2hsv(I);%如果需要求解RGB空间下的颜色矩，则将此条语句去掉即可 [n,m,q]=size(I1); sum1=n*m; h=I1(:,:,1); s=I1(:,:,2); v=I1(:,:,3); y=zeros(1,9); %求一阶矩（均值） for i=1:sum1 y(1)=y(1)+h(i); y(2)=y(2)+s(i); y(3)=y(3)+v(i); end y(1)=y(1)/sum1; y(2)=y(2)/sum1; y(3)=y(3)/sum1; %求二阶矩（方差） for i=1:sum1 y(4)=y(4)+(h(i)-y(1))^2; y(5)=y(5)+(s(i)-y(2))^2; y(6)=y(6)+(v(i)-y(3))^2; end 47 新疆大学信息科学与工程学院 y(4)=sqrt(y(4)/sum1); y(5)=sqrt(y(5)/sum1); y(6)=sqrt(y(6)/sum1); %求三阶矩 for i=1:sum1 y(7)=y(7)+(h(i)-y(1))^3; y(8)=y(8)+(s(i)-y(2))^3; y(9)=y(9)+(v(i)-y(3))^3; end y(7)=(y(7)/sum1)^(1/3); y(8)=(y(8)/sum1)^(1/3); y(9)=(y(9)/sum1)^(1/3); y %绘制颜色直方图： 参考程序如下： %读入图像 Image=imread(B); %end [M,N,O] = size(Image); %M = 256; %N = 256; % 计算每一幅图像的颜色直方图 [h,s,v] = rgb2hsv(Image); H = h; S = s; V = v; h = h*360; %将hsv空间非等间隔量化： % h量化成16级； % s量化成4级； % v量化成4级； for i = 1:M for j = 1:N if h(i,j)<=15||h(i,j)>345 H(i,j) = 0; end if h(i,j)<=25&&h(i,j)>15 H(i,j) = 1; end if h(i,j)<=45&&h(i,j)>25 H(i,j) = 2; end if h(i,j)<=55&&h(i,j)>45 48 新疆大学信息科学与工程学院 H(i,j) = 3; end if h(i,j)<=80&&h(i,j)>55 H(i,j) = 4; end if h(i,j)<=108&&h(i,j)>80 H(i,j) = 5; end if h(i,j)<=140&&h(i,j)>108 H(i,j) = 6; end if h(i,j)<=165&&h(i,j)>140 H(i,j) = 7; end if h(i,j)<=190&&h(i,j)>165 H(i,j) = 8; end if h(i,j)<=220&&h(i,j)>190 H(i,j) = 9; end if h(i,j)<=255&&h(i,j)>220 H(i,j) = 10; end if h(i,j)<=275&&h(i,j)>255 H(i,j) = 11; end if h(i,j)<=290&&h(i,j)>275 H(i,j) = 12; end if h(i,j)<=316&&h(i,j)>290 H(i,j) = 13; end if h(i,j)<=330&&h(i,j)>316 H(i,j) = 14; end if h(i,j)<=345&&h(i,j)>330 H(i,j) = 15; end end end for i = 1:M for j = 1:N if s(i,j)<=0.15&&s(i,j)>0 S(i,j) = 0; 49 新疆大学信息科学与工程学院 end if s(i,j)<=0.4&&s(i,j)>0.15 S(i,j) = 1; end if s(i,j)<=0.75&&s(i,j)>0.4 S(i,j) = 2; end if s(i,j)<=1&&s(i,j)>0.75 S(i,j) = 3; end end end for i = 1:M for j = 1:N if v(i,j)<=0.15&&v(i,j)>0 V(i,j) = 0; end if v(i,j)<=0.4&&v(i,j)>0.15 V(i,j) = 1; end if v(i,j)<=0.75&&v(i,j)>0.4 V(i,j) = 2; end if v(i,j)<=1&&v(i,j)>0.75 V(i,j) = 3; end end end %将三个颜色分量合成为一维特征向量：L = H*Qs*Qv+S*Qv+v；Qs,Qv分别是S和V的量化级数, L取值范围[0,255] %取Qs = 4; Qv = 4 for i = 1:M for j = 1:N L(i,j) = H(i,j)*16+S(i,j)*4+V(i,j); end end %计算L的直方图 for i = 0:255 Hist(i+1) = size(find(L==i),1); end % Hist = Hist/sum(Hist); T=Hist; stem(hist(:)); 50