第7章 MATLAB数字图像处理

在数学建模中有时要用MATLAB进行数字图像处理，这介绍几个最基本的方法及命令。MATLAB中的绝大多数的运算都是通过矩阵这一形式进行的。这一特点也就决定了MATLAB在处理数字图像上的独特优势。理论上讲，图像是一种二维的连续函数，然而在计算机上对图像进行数字处理的时候，首先必须对其在空间和亮度上进行数字化，这就是图像的采样和量化的过程。二维图像进行均匀采样，就可以得到一幅离散化成M×N样本的数字图像，该数字图像是一个整数阵列，因而用矩阵来描述该数字图像是最直观最简便的了。而MATLAB的长处就是处理矩阵运算，因此用MATLAB处理数字图像非常的方便。

MATLAB支持五种图像类型，即索引图像、灰度图像、二值图像、RGB图像和多帧图像阵列；支持BMP、GIF、HDF、JPEG、PCX、PNG、TIFF、XWD、CUR、ICO等图像文件格式的读，写和显示。MATLAB对图像的处理功能主要集中在它的图像处理工具箱（Image Processing Toolbox）中。图像处理工具箱是由一系列支持图像处理操作的函数组成，可以进行诸如几何操作、线性滤波和滤波器设计、图像变换、图像分析与图像增强、二值图像操作以及形态学处理等图像处理操作。

7.1 图像的读取

将图1存入某一文件夹，比如MATLAB中，取名floqwers.tif.

图7.1 照片：花儿

MATLAB中从图像文件中读取数据用函数imread()，这个函数的作用就是将图像文件的数据读入矩阵中.用imshow()显示。 例1：图像数据及图像信息的读取

[A,M]=imread('flowers','tif'); %图像数据的读取，将图像数据放入矩阵A中，颜色数据放入矩阵M中.

subplot(1,2,1);imshow(A,M);title('原图像'); M(:,1)=0;%将颜色数据矩阵的一列置零

subplot(1,2,2);imshow(A,M);title('改变颜色后的图像')

原图像改变颜色后的图像

图7.2 图像数据的读取

MATLAB还提供了将数据写入一个文件的函数imwrite以及不同类型文件相互转换的函数，在此就不一一赘述了，可以参考MATLAB 的帮助文件。

7.2 灰度直方图及直方图均衡化

灰度直方图用于显示图像的灰度值分布情况，是数字图像处理中最简单和最实用的工具。MATLAB中提供了专门绘制直方图的函数imhist()。用它可以很简单的绘制出一幅图像的灰度直方图。 例2：直方图的显示

subplot(1,2,1);imshow('flowers.tif');title('原图像') %显示原图像 A=imread('flowers.tif','tif');

imwrite(rgb2gray(A),'PicSampleGray.tif'); A=rgb2gray(A); %将彩色图片灰度化并保存 subplot(1,2,2);imhist(A),title('对应直方图')

对应直方图8000原图像700060005000400030002000100000100200

图7.3 直方图的显示

在图像处理中，点运算是简单而又重要的一种技术，其中最常用的一种应用就是直方图的均衡化（见例3）。 %例3：直方图均衡化

subplot(2,3,1);imshow('flowers.tif');title('原图像') %显示原图像

A=imread('flowers.tif','tif');

imwrite(rgb2gray(A),'PicSampleGray.tif'); A=rgb2gray(A); %将彩色图片灰度化并保存

subplot(2,3,2);imhist(A),title('对应直方图')

%从得到的直方图可以看出，图像的对比度很低，灰度级集中在70-160范围内，如果只取 %这个范围内的灰度，并扩展到[0,255]，则会明显增强图像对比度 J=imadjust(A,[70/255 160/255],[]);

subplot(2,3,3);imshow(J),title('经灰度级调整后的图') subplot(2,3,4);imhist(J),title('灰度级调整后的直方图')

%MATLAB还提供了histeq函数（自动直方图均衡化） K=histeq(A);

subplot(2,3,5);imshow(K),title('经直方图均衡化后的图') subplot(2,3,6);imhist(K),title('直方图均衡化后的直方图')

对应直方图原图像8000600040002000004灰度级调整后的直方图x 10经灰度级调整后的图100200直方图均衡化后的直方图32100100200经直方图均衡化后的图10000500000100200

图7.4 直方图均衡化

7.3图像的代数运算

代数运算是指对两幅输入图像进行点对点的加、减、乘和除计算而得到输出图像的运算。对于相加和相乘的情形，可能不止有两幅图像参加运算。图像相加的一个重要应用是对同一场景的多幅图像求平均值。这点被经常用来有效地降低加性（additive）随机噪声的影响。

%例4：图象加噪声再通过多次相加求平均的方法祛除噪声

[I,M]=imread('flowers.tif','tif'); J=imnoise(I,'salt & pepper',0.005); subplot(1,3,1),imshow(I,M),title('原图像'); subplot(1,3,2),imshow(J,M),title('加噪声后图像');

K=zeros(size(J)); for i=1:100

J=imnoise(I,'salt & pepper',0.005); J1=im2double(J); K=K+J1; end

K=K/100;%ÇóͼÏñµÄƽ¾ù

subplot(1,3,3);imshow(K),title('相加求平均后的图像');

原图像加噪声后图像相加求平均后的图像

图7.5 多次平均去噪声

7.4 图像滤波处理

在数字图像处理中，常常会遇到图像中混杂有许多的噪声。因此，在进行图像处理中，有时要先进行祛除噪声的工作。最常用的祛除噪声的方法是用滤波器进行滤波处理。MATLAB的图像处理工具箱里也设计了许多的滤波器。如均值滤波器、中值滤波器、维纳滤波器等。用户可以很方便的运用一些函数完成数字滤波工作。。

%例5：用滤波器祛除图象噪声（分别用均值滤波，中值滤波，及维纳滤波器祛除加入高斯噪声的图象）

I=imread('flowers.tif','tif'); %将彩色图片灰度化并保存

imwrite(rgb2gray(I),'PicSampleGray.tif'); I=rgb2gray(I);

J=imnoise(I,'gaussian',0,0.002);%加入高斯噪声 %进行均值滤波

h=fspecial('average',3); I2=uint8(round(filter2(h,I))); %进行中值滤波 I3=medfilt2(J,[3,3]); %进行维纳滤波

I4=wiener2(J,[3,3]);%进行一次维纳滤波 I5=wiener2(I4,[3,3]);%进行二次维纳滤波 subplot(2,3,1),imshow(I),title('原图象') subplot(2,3,2),imshow(J),title('加噪声图象') subplot(2,3,3),imshow(I2),title('均值滤波后图象') subplot(2,3,4),imshow(I3),title('中值滤波后图象') subplot(2,3,5),imshow(I4),title('维纳滤波后图象') subplot(2,3,6),imshow(I5),title('两次维纳滤波后图象')

原图象加噪声图象均值滤波后图象中值滤波后图象维纳滤波后图象两次维纳滤波后图象

图7.6 滤波器祛除图像噪声 注：彩色图片要灰度化

7.5、傅立叶变换

傅立叶变换是线性系统分析的一个有力的工具。它在图像处理，特别是在图像增强、复原和压缩中，扮演着非常重要的作用。实际中一般采用一种叫做快速傅立叶变换（FFT）的方法，MATLAB中的fft2指令用于得到二维FFT的结果，ifft2指令用于得到二维FFT逆变换的结果。（见例6）

%例6：近似冲击函数的二维快速傅立叶变换（FFT） x=1:49;y=1:49; [X,Y]=meshgrid(x,y); A=zeros(49,49); A(24:26,24:26)=1; B=fft2(A);

subplot(2,2,1),imshow(A),xlabel('空域图象'); subplot(2,2,2),imshow(B),xlabel('时域图象'); subplot(2,2,3),mesh(X,Y,A),xlabel('空域'),grid on; subplot(2,2,4),mesh(X,Y,abs(B)),xlabel('时域'),grid on;

图6 近似冲击函数的傅立叶变换

7.6 图像压缩

在图像的变换和压缩中，常常用到离散余弦变换（DCT）。DCT具有能使图像的最重要的信息集中在DCT的几个系数上的性能。正是基于此，DCT通常应用于图像的压缩。（见例7） %例7：DCT变换用于图象的压缩实例 I=imread('flowers.tif'); % 输入图像

imwrite(rgb2gray(I),'PicSampleGray.tif'); I=rgb2gray(I); I = im2double(I); % 数据类型转换 T = dctmtx(8); % 计算二维离散DCT矩阵 dct = @(x)T * x * T'; % 设置函数句柄 B = blkproc(I,[8 8],dct); % 图像块处理 mask = [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];

B2 = blkproc(B,[8 8],@(x)mask.* x); % 图像块处理 invdct = @(x)T' * x * T; % 设置函数句柄 I2 = blkproc(B2,[8 8],invdct); % 图像块处理 subplot(1,2,1);imshow(I);title('原始图像') subplot(1,2,2);imshow(I2);title('压缩重构图像')

原始图像压缩重构图像

图7.8 DCT变换用于图像压缩

可以看出，尽管由于85%的DCT系数被抛弃而使恢复后的图像质量有所降低，图像内容仍能清晰可辨，达到了图像压缩的目的。

7.7图像对象的面积计算

有时候，需要就是图像对象的面积，bwarea的功能就是就是二进制图像对象的面积。 例8 利用bwarea计算图像对象的面积.

图 7.9 circles.tif 图7.9. circles.tif I=imread('circles.tif'); imshow(I);

imwrite(rgb2gray(I),'PicSampleGray.tif'); I=rgb2gray(I); bwarea(I)

ans =

470

注：该面积和二值图像中对象的像素数目不一定相等。

7.8 删除小面积对象

函数bwareaopen删除小面积对象,如 I2 = bwareaopen(BW,P,conn)

删除二值图像BW中面积小于P的对象，默认情况下使用8邻域(conn)。 例9：利用bwareaopen删除二值图像中面积小于P的对象 originalBW = imread('text.png'); subplot(1,2,1);imshow(originalBW) bwAreaOpenBW = bwareaopen(originalBW,50); subplot(1,2,2);imshow(bwAreaOpenBW)

图7.10 text.png

例7.10. 利用bwareaopen删除非二值图像中面积小于P的对象 I=imread('circles2.png');

imwrite(rgb2gray(I),'PicSampleGray.tif'); I=rgb2gray(I);%灰度处理 %二值处理 PS=I;

for i=0:255

I(find(PS<150))=0; I(find(PS>=150))=250; end

subplot(1,2,1);imshow(I) bwAreaOpenBW = bwareaopen(I,50); subplot(1,2,2);imshow(bwAreaOpenBW)

图7.11 circles2.png

7.9 图像的边缘检测

MATLAB的图像处理工具箱中提供的edge函数可以实现检测边缘的功能，其语法格式如下：BW = edge(I,'sobel') BW = edge(I,'sobel',direction) BW = edge(I,'roberts') BW = edge(I,'log')

这里BW = edge(I,'sobel')采用Sobel算子进行边缘检测。BW = edge(I,'sobel',direction)可以指定算子方向，即：

direction=’horizontal’，为水平方向； direction=’vertical’，为垂直方向； direction=’both’，为水平和垂直两个方向。

BW = edge(I,'roberts')和BW = edge(I,'log')分别为用Roberts算子和拉普拉斯高斯算子进行边缘检测。

例：用三种算子进行边缘检测。 I=imread('circles2.png ');

imwrite(rgb2gray(I),'PicSampleGray.tif'); I=rgb2gray(I);%灰度处理 subplot(2,2,1);imshow(I) BW1=edge(I,'roberts');

subplot(2,2,2);,imshow(BW1),title('用Roberts算子') BW2=edge(I,'sobel');

subplot(2,2,3);,imshow(BW2),title('用Sobel算子 ') BW3=edge(I,'log');

subplot(2,2,4);,imshow(BW3),title('用拉普拉斯高斯算子')

用Roberts算子用Sobel算子 用拉普拉斯高斯算子

图12 边缘检测

7.10 三维重建

在对生物组织、器官等的形态的研究中，需要知道组织、器官的三维形态。其中的一种途径是通过样本切片来重建其三维图象。其大体的步骤为：用切片机对研究的样本进行平行切片；将切片拍照并采样得到平行切片的数字图象；运用计算机建立三维重建的算法，根据二维数字图象的信息重现原样本的三维形态。

我们现在研究的对象是一种血管样本。现在已经得到了该血管的100张平行切片的图象，每张图象的象素点为512×512个，格式为bmp。每幅图转换后的大小为512×512个像素，重建血管的三维形态。（2001 年大学生数学建模比赛中第一题，可在如下网址上

下载：http://www.mcm.edu.cn/html_cn/node/2f00df27925ba6e9735be9df4b8872f0.html） 整体三维重建步骤为：

第1步 、数据采集： MATLAB 提供了图像数据采集的命令 imread ，对采集到的 100 幅 512x512 象素的平行切片图像（图像文件名依次为 0 . bmp 、 l . bmp 、 … 、 99 . bmp ，格式均为 BMP ) 构造出三维体数据集 D ，其中 D 是 512x512x100 的三维数据集。把 100 张图像文件放到 MATLAB6的搜索路径中，建立 m 文件，运行如下语句即可得到数据集 D 。 for i=1:100

n=strcat(num2str(i-1), '.bmp ') ; %生成第 n 幅图像的文件名 [x,m]=imread(n) ; %读取第 n 幅图像的数据 D(:,:,i)=x ; %利用三维数组构造数据矩阵D end ;

第2步 、数据处理：由上面步骤构造出的矩阵数据量大，占用内存多，给处理带来一定困难。根据实际情况，在不影响结果的情况下，对数据进行预处理，以减少所处理的数据量。 [ x y z D]＝reducevolume ( D,[a b c]) ; D = smooth3 ( D ，’filter’)；其中 a 、 b 、 c 三个参数为在 x 、 y 、 z 轴上的数据抽取比例， filter 可以为’ gaussian ’或者’ box ' （缺省）表示采用滤子对数据平滑处理，书写如下代码：

[ x y z D] = reducevolume ( D,[ 4 4 1]) ; ％在 x 、 y 、 z 轴上按比例 [ 4 4 1]抽取数据， D 变为 128xl28xl00 的三维数组

D = smooth3 ( D , ' gaussian ' ) ; ％使用高斯滤子对数据 D 平滑处理；

第3步、体数据在标准透明表面的绘制：利用 isosurface 函数计算体数据集在显示平面累计投影。

fv = isosurface ( x , y , z , D , isovalue ) ；其中 x , y , z , D 在上面已经得到， isovalue 根据实际情况自定，接着运行如下代码：

fv = isosurface ( x , y , z , D ,’ noshare’ ) ;

第4步 、构造结果图像碎片：使用 patch 函数来对图像子区域进行分类，可以定义结果图像的颜色、光线等信息。

pl = patch ( fv ，’FaceColor ‘, yourscolor1 , ‘EdgeColor ‘, yourscolor2) ；其中 yourscolor1 和 yourscolor2 为自定义的所期望显示的颜色，如 ' green ’、’ blue ' ，也可为’ none ' ; 第5步 、计算碎片的法线方向：利用 isnormals 函数计算第四步图形子区域顶部法线方向。

isonomials ( x , y , z , D , p )

p = patch ( fv, ’FaceColor’ ,’ red ‘, ‘EdgeColor’ , ‘none’ ) ;

第6步 、计算三维图像的几何边界：利用 isocaPs 函数计算三维图像的几何边界。 fvc = isocaps ( x , y , z , D , isovalue ) ;

第7步 、利用第 6 步所得结果替换 fv ，重复第 4 步；

第8步 、图像显示：利用 View 、 daspect 、 colormap 、 camlight 、 lishting 等函数显示图像。其中 view 函数定义观察视角， daspect 定义 x 轴、 y 轴和 Z 轴的显示比例， colormap 定义图像的颜色集， camligh 定义照相机光线的位置， lighting 定义显示图像的光线阴影。

血管局部三维重建

（1） 、局部数据的提取：利用 subvolume 函数对整体三维重建中得到的数据进行局部提取。 [ x y z D ] ＝ suhvolume ( D ， [ xl yl 21 x2 y2 z2]) ; 其中 [ xl yl 21 x2 y2 z2] 为数据集中所要提取的局部信息的范围.

（2） 、重复整体重建的第 3 、 4 、 5 、 6 、 7 、 8 步。

clear;

for i=1:100

n=strcat(num2str(i-1), '.bmp ') ; %生成第 n 幅图像的文件名 [x,m]=imread(n) ; %读取第 n 幅图像的数据 D(:,:,i)=x ; %利用三维数组构造数据矩阵D end ; close all

[x y z D]=reducevolume(D,[4 4 1]); D=smooth3(D,'gaussian');

fv = isosurface ( x , y , z , D, 'noshare ');

pl = patch(fv ,'FaceColor','red' ,'EdgeColor','none') ; p = patch ( fv, 'FaceColor' ,' red ', 'EdgeColor' , 'none' ) ; view(3); axis tight; daspect([1 1 .4]); colormap(gray(4000)); isonormals(x,y,z,D,p);

lighting phong material shiny

light('color','w','position',[-1 0.5 1],'style','local') axis equal box

图7.9 血管三维重建