【MATLAB Image Processing Toolbox 入门教程二】快速入门之“亮度校正”和“目标识别”

校正亮度不均匀问题并分析前景对象

一般在分析图像特征之前，需要对图像进行某种预处理，以便分析。此示例演示了对于背景亮度不均的图片（MATLAB自带的大米图片），先将其转为二值图像，以便于识别每个米粒，接着对对象进行分析，例如计算米粒面积、米粒的数量。

图像预处理

本次示例我们采用的图片是MATLAB自带的大米图片（‘rice.png’），读入图片并显示：

srcImage=imread('rice.png');
imshow(srcImage);

可以看出，图片中间背景比周围更亮。为方便分析，对图像预处理，使得背景亮度更加均匀。
首先，使用形态学开运算删除所有前景，此图像中米粒即为前景，开运算的作用是删除图像中无法包含结构元素的对象。首先我们需要使用strel函数创建结构元素（strel=structure+element）。

结构元素包含平面结构元素和非平面结构元素，非平面结构元素用offsetstrel函数构建，此处我们先详述平面结构元素。平面结构元素是一个二维或多维的二值领域，其中真像素（1）用在形态学运算中，假像素（0）不用于运算。结构元素的中心称为原点，对应于被处理图像中正在处理的像素。strel函数创建的结构元素可用于二值图像和灰度图像的处理，在腐蚀和膨胀的操作中起到重要作用。下图就是一个半径为3的碟形结构元素。

strel函数有下面几种使用形式：

1 SE=strel(nhood); %创建一个指定领域的平面结构元素，参数nhood为指定为特定维度的数值数组，nhood的中心或原点是他的中心元素，计算方式为：floor((size(nhood)+1)/2)；floor(x)：将x的每个元素四舍五入到小于或等于该元素的最接近整数。

2 SE=strel('diamond',r); %创建一个菱形的结构元素，r是结构元素原点到菱形的距离。

3 SE=strel('disk',r,n); %创建一个盘形的结构元素，r是盘的半径，n为结构元素为了近似形状所用的周期线结构元素数量。盘形结构元素较为常用，且一般情况下只需要指定盘形的半径。

4 SE=strel('octagon',r); %创建一个八角形结构元素，r为结构元素原点到八角形边的距离，注意r必须为3的非负倍数。

5 SE=strel('line',len,deg); %创建相对领域中心对称的线性结构元素，len是线性元素的长度，deg为线性元素的角度，沿水平轴逆时针测量。

6 SE=strel('rectangle',[m n]); %创建一个矩形结构元素，矩形尺寸为[m n]。

7 SE=strel('square',w); %创建一个正方形结构元素，边长为w个像素。

8 SE=strel('cube',w); %创建一个三维立方结构元素，方块的边长为w个像素。

9 SE=strel('cuboid',[m n p]); %创建一个尺寸为[m n p]的三维立方结构，m行，n列，p个平面。

10 SE=strel('sphere',r); %创建一个半径为r的三维球形结构元素。

既然我们需要用到开运算删除图片中的所有米粒，就需要创建比米粒大的结构元素。采用盘形的结构元素，现在我们要做的就是确定盘形半径。依旧采用imdistline函数，继续输入：

d=imdistline;

可以确定，盘形结构元素半径设置为15即可，只要比米粒大就行。确定了结构元素之后，我们先对形态学开运算函数imopen进行初步了解：

J=imopen(I,SE)
对灰度图像或二值图像I执行形态学开运算，返回图像J。开运算是一种先腐蚀后膨胀的操作，腐蚀操作和膨胀操作使用的都是同一个结构元素SE。

继续输入：

delete(d);
se=strel('disk',15); %创建盘形结构元素后，使用imopen函数进行形态学开运算
background=imopen(srcImage,se);
imshow(background);

这样我们就得到了米粒图片的背景图，用原图减去背景图，就可得到除了米粒之外其余的背景全是黑色的图片，继续输入：

srcImage1=srcImage-background;
imshow(srcImage1);

注意，上面的两个步骤：先进行形态学开运算再从原始图像中减去它，在图像处理中成为顶帽，我们可以用MATLAB中的imtophat函数直接实现它：

J=imtophat(I,SE) 对灰度或二值图像I进行顶帽滤波，并返回滤波后的图像。首先进行生态学开运算，再从原图像中减去它。

当前生成的图像具有均匀的背景，但是米粒的颜色较暗，整个图像的对比度较小，对此图像查看灰度直方图：

imhist(srcImage1);

我们可以使用直方图均衡化增强对比度：

srcImage2=histeq(srcImage1);
imshow(srcImage2);

可以看出，直方图均衡化后图像背景不均匀，更难以进行下一步分析，显然这不是个好方法。好在MATLAB提供了另外一个函数：imadjust。

J=imadjust(I) 将灰度图像I中的像素强度进行重新映射到图像J，默认情况下，imadjust将I的所有像素值中最低的1%和最高的1%进行饱和处理。

让我们使用imadjust重新处理srcImage1，并观察其灰度直方图。
继续输入：

srcImage3=imadjust(srcImage1);
imshow(srcImage3);
figure;
imhist(srcImage3);

可以看出，图像中的米粒比之前更加清晰。接下来我们使用imbinarize函数将灰度图像转为二值图像，并使用bwareaopen函数去除背景噪声。我们先对这两个函数进行了解。

BW=imbinarize(I) 将灰度图像I中所有高于阈值的像素值替换为1，低于阈值的替换为0。默认情况下，imbinarize函数使用的是Otus方法。
BW2=bwareaopen(BW,P) 从二值图像BW中删除少于P个元素的所有连通分量，生成另一个二值图像BW2。

了解这些之后，我们继续输入：

bw=imbinarize(srcImage3); %有些版本的MATLAB中灰度转二值图像使用的是im2bw函数，所以报错的同学可以试试这个函数。
bw1=bwareaopen(bw,50);
imshow(bw1);

至此，我们可以得到清晰的米粒与背景颜色分开的图像，可以进行下一步——对象识别。

对象识别

使用函数bwconncomp查找二值图像bw1中所有的连通分量，对于bwconncomp函数，使用方法为：

CC=bwconncomp(BW,conn) 返回二值图像BW中的连通分量CC，conn指定连通分量所需要的连通性。

对于二维连通，conn取值为4或8，含义如下：

查找bw1中的连通分量，继续输入：

cc=bwconncomp(bw1,4)

返回的cc为：

一共检测出95个连通分量，即95个米粒。若想要查看第49个米粒，可以输入以下代码：

grain_49=false(size(bw1)); %创建一个和二值图像bw1尺寸相同的全0数组
grain_49(cc.PixelIdxList{49})=true; %将第49个米粒所在位置处的元素设置为真（1）
imshow(grain_49);

为了将检测出的所有米粒都显示出来，使用labelmatrix函数根据bwconncomp的输出cc创建标签矩阵。继续输入：

label=labelmatrix(cc);
imshow(label);

可以看出，米粒从左至右逐渐明显，这是为什么呢？查看label数组我们发现，在label中，第n颗米粒的灰度值全为n，所以上图中最明显的米粒灰度值也只不过为95。
使用label2rgb函数将标签图像转化为rgb彩色图像：

rgb_label=label2rgb(label); %label2rgb函数根据标签矩阵中对象的数量确定要分配给每个对象的颜色。
rgb_label1=label2rgb(label,'spring','c','shuffle'); %使用“spring”颜色贴图，将背景像素设置为青色，并随机化将颜色指定给标签。
imshow(rgb_label);
figure;
imshow(rgb_label1);

计算对象的面积统计量

计算面积统计量依旧基于刚才检测的连通分量cc，使用regionprops函数：

stats=regionprops(CC,properties) 返回CC中每一个连通分量的属性，properties指定测量的类型，此处选择basic即可

输入：

graindata=regionprops(cc,'basic');

可以看出，返回的graindata是包含了95个连通分量的area（面积）、centroid（质心）和boundingbox（边界线）的数组，我们将area部分单独分离出来：

grain_areas=[graindata.Area];

%找到最小面积的米粒并将其显示出来
[min_area,num]=min(grain_areas);
grain=false(size(bw));
grain(cc.PixelIdxList{num})=true;
imshow(grain);

这样，我们就得到了最小面积米粒的位置图像。

代码综合部分

下面代码演示了读入图像，预处理操作：背景分离、提高对比度、转化二值图像、去除背景噪声，通过检测连通分量得到识别的目标并将其可视化，最后对目标的面积分量进行统计。废话不多说，上代码：

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%本实例中imbinarize函数在旧版本中为im2bw,histogram函数在旧版本中为hist，如果程序报错替换相应函数即可
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%载入图片
srcImage=imread('rice.png');

%背景分离，得到均匀背景图像
se=strel('disk',15);
background=imopen(srcImage,se);
srcImage1=srcImage-background;

%提高图像对比度
srcImage2=imadjust(srcImage1);

%灰度图像转化为二值图像，并去除背景噪声
bw=imbinarize(srcImage2);
bw1=bwareaopen(bw,50);

%检测二值图像中的连通分量，即米粒
cc=bwconncomp(bw1,4);
num=cc.NumObjects; %num即为图像中检测到的米粒的数量

%创建标签矩阵，并将其转化为可视化伪彩色图像
label=labelmatrix(cc);
rgb_label=label2rgb(label,'spring','c','shuffle');
imshow(rgb_label);

%基于标签分量label创建米粒面积直方图
graindata=regionprops(label,'basic');
grain_areas=[graindata.Area];
figure;
histogram(grain_areas);
title('米粒面积直方图');

本次内容就这些，下篇见：）