《数字图像处理与机器视觉》案例四基于分水岭算法的粘连物体的分割与计数

一、引言

分水岭算法（Watershed Algorithm），是一种基于拓扑理论的数学形态学的分割方法，其基本思想是把图像看作是测地学上的拓扑地貌，图像中每一点像素的灰度值表示该点的海拔高度，每一个局部极小值及其影响区域称为集水盆，而集水盆的边界则形成分水岭。分水岭的概念和形成可以通过模拟浸入过程来说明。在每一个局部极小值表面，刺穿一个小孔，然后把整个模型慢慢浸入水中，随着浸入的加深，每一个局部极小值的影响域慢慢向外扩展，在两个集水盆汇合处构筑大坝，即形成分水岭。分水岭算法可以应用于分割粘连对象和对象计数。在该算法中，空间上相邻并且灰度值相近的像素被划分为一个区域。一般来讲，直接应用分水岭分割算法的效果往往并不好，如果在图像中对前景对象和背景对象进行标注区别，再应用分水岭算法会取得较好的分割效果。下面以具有粘连大米的分割计数为例，开发了基于MATLAB的分水岭算法大米颗粒的分割计数，算法主要步骤如下图所示。

二、程序代码

clear all;
close all;
clc;
%%步骤1 读图并进行背景处理
I0=imread('.\rice2.jpg');
I_r=I0(:,:,1); % 提取R通道图像
I_g=I0(:,:,2);
I_b=I0(:,:,3);
I_rb=I_r-I_b;
figure,imhist(I_rb);%具有明显的双峰特性
I_bw=1-im2bw(I_rb,45/255);             %使用指定阈值分割
%I_bw=1-im2bw(I_rb,graythresh(I_rb));   %使用OTSU方法进行阈值分割并反色
%I_bw=1-imbinarize(I_rb,'adaptive');
figure,imshow(I_bw);
Obj1_r=uint8(I_bw).*I_r;
Obj1_g=uint8(I_bw).*I_g;
Obj1_b=uint8(I_bw).*I_b;
rgb=cat(3,Obj1_r,Obj1_g,Obj1_b);
figure,imshow(rgb);
%%步骤2 利用梯度实现灰度图像的分割
I = rgb2gray(rgb);
hy = fspecial('sobel');hx = hy';    %使用sobel算子进行边缘检测
Iy = imfilter(double(I), hy, 'replicate');
Ix = imfilter(double(I), hx, 'replicate');
gradmag = sqrt(Ix.^2 + Iy.^2);     %求模
%%步骤3 标记前景对象
se = strel('diamond',3);                  %构建结构元素
Io = imopen(I, se);                         %形态学开操作
Ie = imerode(I, se);                        %腐蚀
Iobr = imreconstruct(Ie, I);             %重建
Ioc = imclose(Io, se);                     %形态学闭操作
Iobrd = imdilate(Iobr, se);               %膨胀
Iobrcbr = imreconstruct(imcomplement(Iobrd), imcomplement(Iobr));
Iobrcbr = imcomplement(Iobrcbr); %求反
%%步骤4 获取局部最大值
fgm = imregionalmax(Iobrcbr);          %获得局部最大值
It1 = rgb(:, :, 1); It2 = rgb(:, :, 2);
It3 = rgb(:, :, 3);
It1(fgm) = 255; It2(fgm) = 0; It3(fgm) = 0;
I2 = cat(3, It1, It2, It3);
%%步骤5
se2 = strel(ones(3,3));
fgm2 = imclose(fgm, se2);                    %闭操作
fgm3 = imerode(fgm2, se2);                 %腐蚀
fgm4 = bwareaopen(fgm3, 20);            %删除
It1 = rgb(:, :, 1);
It2 = rgb(:, :, 2);
It3 = rgb(:, :, 3);
It1(fgm4) = 255;                                     %前景设置为255
It2(fgm4) = 0;
It3(fgm4) = 0;
I3 = cat(3, It1, It2, It3);
%%步骤6 标记背景
bw = im2bw(Iobrcbr, graythresh(Iobrcbr));          %二值转换
D = bwdist(bw);                                                   %计算距离
DL = watershed(D);                                             %分水岭变换
bgm = DL == 0;
gradmag2 = imimposemin(gradmag, bgm | fgm4);     %置最小值
L = watershed(gradmag2);
It1 = rgb(:, :, 1);
It2 = rgb(:, :, 2);
It3 = rgb(:, :, 3);
fgm5 = imdilate(L == 0, ones(9, 2)) | bgm | fgm4;   %前景及边界处置255
It1(fgm5) = 255; It2(fgm5) = 0; It3(fgm5) = 0;
I4 = cat(3, It1, It2, It3);
figure,imshow(I4);
%%步骤7 去粘连、填孔、分割
figure('units', 'normalized', 'position', [0 0 1 1]);
P=rgb2gray(I4);                                 % 彩色图像转灰度
figure,imshow(P);
figure,imhist(P);              %具有明显的双峰特性
P1=im2bw(P,200/255);                      % 灰度图像二值化
SE1=strel('square',3);
P2=imopen(P1,SE1);                        % 以SE为单位对P1进行开运算
SE2=strel('square',2);                       % 构造半径为2圆形元素
P3=imerode(P2,SE2);             % 形态学腐蚀运算,部分目標物有粘连现象,去除粘连
P4 = imfill(P3,'hole');
SE3=strel('square',2);
P5=imopen(P4,SE3);                        % 以SE3作开运算
SE4=strel('disk',3);
%SE4=strel('disk',5);
P6=imerode(P5,SE4);                       %形态学腐蚀运算,进一步去除粘连
P7 = bwareaopen (P6, 100);
figure,imshow(P7);
%%步骤8 计算连通数，输出结果
[L,N]=bwlabel(P7,8);            % 计算连通数,N即为目标个数
figure,imshow(I0);
hold on
for k=1:N
[r,c] = find(L == k);
rbar = mean(r);
cbar = mean(c);
plot(cbar,rbar,'marker','+','markeredgecolor','b','markersize',6); % 对话框显示目标物个数,通过坐标定点,标记点为“+”
end
title(['大米总数为',num2str(N)],'Fontsize',14)
figure,
subplot(1,3,1),imshow(I0),title('存在粘连的大米彩色图像');
subplot(1,3,2),imshow(I4),title('分水岭算法分割图像');
subplot(1,3,3),imshow(P7),title('数学形态学运算后的图像');