图像连通域分析及相关算法研究——很详细的原理以及简单的代码
1、连通域的基本概念
连通域一般是指图像中具有相同像素值且位置相邻的像素点组成的图像区域。连通域分析是指将图像中的各个连通区域找出并标记。连通区域分析在图像分析处理的众多应用领域非常常用。连通区域分析处理的对象一般是一张二值化后的图像。像素邻域关系一般有四邻域和八邻域两种，下面是以四邻域做介绍。
2、图像灰度化
（1）最大类间方差法原理
上面说到连通域分析处理的一般是二值化后的图像，那么灰度图像怎么二值化呢？我们可以用到otsu算法，即最大类间方差法。下面简单介绍一下最大类间方差法。
最大类间方差法最重要的是这个公式：
可以看到原理非常简单，类间方差和方差的原理差不多，方差反应所有数据的偏差情况，而类间方差反应的是数据分成两部分之后两部分的偏差情况。举个简单例子来说：
假设有5个数据1、2、3、4、5，那么它们的平均值为3，那么方差就是1/5{(1-3)2+(2-3)2+(3-3)2+(4-3)2+(5-3)2}。那么类方差呢？假设我们设置阈值T为3，小于等于3的作为A部分，大于3的作为B部分，那么A部分的均值为2，概率为3/5；B部分的均值为4.5，概率为2/5；类方差就是3/5(2-3)2+2/5(4.5-3)2。
（2）最大类间方差法实现
matlab实现最大类间方差法：
close;clear;clc;
[filename,pathname]=uigetfile({'*.jpg';'*.bmp';'*.gif'},'选择原图片');
img = imread([pathname,filename]);
img = rgb2gray(img);
level=graythresh(img);
img =double(img);
[M,N]=size(img); %得到图像行列像素
number_all=M*N; %总像素值
hui_all=0.0; %预设图像总灰度值为0
ICV_t=0.0; %预设最大方差为0
k=0;
%得到图像总灰度值
for i=1:M
for j=1:N
hui_all=hui_all+img(i,j);
end
end
all_ave=hui_all/number_all; %图像灰度值的总平均值
%t为某个阈值，把原图像分为A部分（每个像素值>=t）与B部分（每个像素值<t）
for t=0:255 %不断试探最优t值
hui_A=0.0; %不断重置A部分总灰度值
hui_B=0.0; %不断重置B部分总灰度值
number_A=0.0; %不断重置A部分总像素
number_B=0.0; %不断重置B部分总像素
for i=1:M %遍历原图像每个像素的灰度值
for j=1:N
if (img(i,j)>=t) %分割出灰度值>=t的像素
number_A=number_A+1; %得到A部分总像素
hui_A=hui_A+img(i,j); %得到A部分总灰度值
elseif (img(i,j)<t) %分割出灰度值《t的像素
number_B=number_B+1; %得到B部分总像素
hui_B=hui_B+img(i,j); %得到B部分总灰度值
end
end
end
PA=number_A/number_all; %得到A部分像素总数与图像总像素的比列
PB=number_B/number_all; %得到B部分像素总数与图像总像素的比列
A_ave=hui_A/number_A; %得到A部分总灰度值与A部分总像素的比例
B_ave=hui_B/number_B; %得到B部分总灰度值与B部分总像素的比例
ICV=PA*((A_ave-all_ave)^2)+PB*((B_ave-all_ave)^2); %Otsu算法
if (ICV>ICV_t) %不断判断，得到最大方差
ICV_t=ICV;
k=t; %得到最大方差的最优阈值
end
end
k %显示阈值
（3）优缺点分析
优点：（1）原理简单，速度较快。
（2）数学上可以证明取最佳阈值时，两部分之间差别最大。
缺点：（1）受噪声影响大（解决方法二维otsu算法）
（2）类间方差准则函数可能呈现双峰或多峰，此时效果不好（解决方法局部阈值法）
可以看到一维最大类间方差方法实现的效果不理想，有很多噪声点。
（4）二维最大类间方差法
二维最大类间方差法相比一维otsu算法抗噪声能力更强，基本原理可以参考这篇博客，也可以找相关论文进行学习，这里用matlab实现了快速二维最大类间方差法。
matlab实现的二维最大类间方差法
clear;
clc;
[filename,pathname]=uigetfile({'*.jpg';'*.bmp';'*.gif'},'选择原图片');
I = imread([pathname,filename]);
Igray=rgb2gray(I);
imageBinary=Igray;
% Igray32=uint32(Igray); %强制类型转换
[Height,Width]=size(Igray);
Igray_ex=[zeros(1,Width+2);zeros(Height,1),Igray,zeros(Height,1);zeros(1,Width+2)];%扩展图像
Igray_ex32=uint32(Igray_ex); %强制类型转换
%%定义变量
xsumm=0; %行像素和
xsumn=0; %
summ=0;
sumn=0;
flag=1; %while循环标志位
graygrade=zeros(256,256);
possiblity=zeros(256,256);
%%求初始阈值
for m = 2:(Height) %每行循环,从第二行第二列开始
xsumm=0; %行像素和清零
xsumn=0; %
for n = 2:(Width) %每列循环
xsumm=xsumm+Igray_ex32(m,n);
temp=Igray_ex32(m-1,n-1)+Igray_ex32(m-1,n)+Igray_ex32(m-1,n+1)...
+Igray_ex32(m,n-1)+Igray_ex32(m,n)+Igray_ex32(m,n+1)...
+Igray_ex32(m+1,n-1)+Igray_ex32(m+1,n)+Igray_ex32(m+1,n+1);
xsumn=xsumn+temp/9;
temp=0;
end
xsumm=xsumm/(Width);
xsumn=xsumn/(Width);
summ=summ+xsumm;
sumn=sumn+xsumn;
end
%%求初始阈值
t0=summ/(Height); %灰度域值邻域平均值
s0=sumn/(Height); %灰度域值
t=t0;
s=s0;
%%求最终阈值
while(flag~=0)
for m=1:256
for n=1:256
graygrade(m,n)=0;
possiblity(m,n)=0;
end
end
%%目标类区域和背景类区域清零
w0a=0;
w1a=0;
w0b=0;
w1b=0;
u0m=0;
u1m=0;
u0n=0;
u1n=0;
%%
for m = 2:(Height) %每行循环,从第二行第二列开始
for n = 2:(Width) %每列循环
f=Igray_ex32(m,n);%像素灰度值
temp=Igray_ex32(m-1,n-1)+Igray_ex32(m-1,n)+Igray_ex32(m-1,n+1)...
+Igray_ex32(m,n-1)+Igray_ex32(m,n)+Igray_ex32(m,n+1)...
+Igray_ex32(m+1,n-1)+Igray_ex32(m+1,n)+Igray_ex32(m+1,n+1);
g=temp/9; % 像素邻域灰度平均值
if(f~=0 && g~=0) %防止下标为0越界
graygrade(f,g)=graygrade(f,g)+1;
end
end
end
%%
for m = 1:256 %每行循环,从第二行第二列开始
for n = 1:256 %每列循环
possiblity(m,n)=graygrade(m,n)/((Width)*(Height));
if(m<=t)
w0a=possiblity(m,n)+w0a;
else
w1a=possiblity(m,n)+w1a;
end
if(n<=s)
w0b=possiblity(m,n)+w0b;
else
w1b=possiblity(m,n)+w1b;
end
end
end
%%
for m = 1:256 %每行循环,从第二行第二列开始
for n = 1:256 %每列循环
if(m<=t)
u0m=m*possiblity(m,n)+u0m;
else
u1m=m*possiblity(m,n)+u1m;
end
if(n<=s)
u0n=n*possiblity(m,n)+u0n;
else
u1n=n*graygrade(m,n)+u1n;
end
end
end
%%
u0m=u0m/w0a;
u0n=u0n/w0b;
u1m=u1m/w1a;
u1n=u1n/(w1b*(Width)*(Height));
tfinal = (u0m+u1m)/2;
sfinal = (u0n+u1n)/2;
if(tfinal>256 || sfinal>256 || tfinal<0 || sfinal<0)
t = t+1;
s = s+1;
else
if (abs(tfinal-t)>1 || abs(sfinal-s)>1)
t = tfinal;
s=sfinal;
else
flag=0;
end
end
end
%%强制类型转换
t=uint8(t);
s=uint8(s);
tfinal=uint8(tfinal);
sfinal=uint8(sfinal);
%%二值分割
for m=1:Height
for n=1:Width
if imageBinary(m,n) >tfinal
imageBinary(m,n) = 255;
else
imageBinary(m,n) = 0;
end
end
end
%%画图
level=graythresh(Igray);
BW=im2bw(Igray,level);
level=uint8(level*256);
subplot(2,2,1),imshow(Igray);title(' 原图')
subplot(2,2,2);imshow(imageBinary);title(' Ostu二值图')
subplot(2,2,4);imshow(BW);title(' 自带函数二值图')
可以看到二维otsu实现的效果相比一维otsu以及matlab自带的阈值分割函数graythresh来说效果更好。但是我们发现二维的otsu算法实现的阈值分割二值图像仍然有很多噪声，我们可以通过开运算进行消除，开运算就是图像先进行腐蚀运算在进行膨胀运算实现，具体原理这里不详细介绍，可以通过matlab自带的函数imopen进行实现。效果图如下：
可以看到开运算能够除去孤立的小点，毛刺，而总的位置和形状不变。
3、连通域标记算法
将图像二值化后我们就可以进行连通域标记了，这里讲解两种连通域标记算法的实现，一个是Two-Pass算法，也叫两步法，另一个是Seed-Filling算法，也叫种子填充法。
（1）Two-Pass算法
1）Two-Pass算法原理
两遍扫描法，指的是通过扫描两遍图像，就可以将图像中存在的所有连通区域找出并标记。思路：第一遍扫描时赋予每个像素位置一个label，扫描过程中同一个连通区域内的像素集合中可能会被赋予一个或多个不同label，因此需要将这些属于同一个连通区域但具有不同值的label合并，也就是记录它们之间的连通关系；第二遍扫描就是将具有连通关系的标记的像素归为一个连通区域并赋予一个相同的label。
下面给出Two-Pass算法的简单步骤，步骤参考的是这篇博客：
（1）第一次扫描：
访问当前像素B(x,y)，如果B(x,y) == 1：a、如果B(x,y)的领域中像素值都为0，则赋予B(x,y)一个新的label：
label += 1， B(x,y) = label；
b、如果B(x,y)的领域中有像素值 > 1的像素Neighbors：
1）将Neighbors中的最小值赋予给B(x,y):
B(x,y) = min{Neighbors}
2）记录Neighbors中各个值（label）之间的相等关系，即这些值（label）同属同一个连通区域；
（2）第二次扫描：
访问当前像素B(x,y)，如果B(x,y) > 1：a、找到与label = B(x,y)同属相等关系的一个最小label值，赋予给B(x,y)；
完成扫描后，图像中具有相同label值的像素就组成了同一个连通区域。
算法图解如下：
2）Two-Pass算法实现
两步法的matlab代码如下
function [ outImg, labels ] = MyBwLabel( inputImg )
%MyBwLabel 自制的连通区域分析函数
% [ outImg, labels ] = MyBwLabel( inputImg )
% inputImg: 输入的图像，要求是二值化图像，且最大值为1
% outputImg: 输出的图像，不同的连通区域的；label不同
% labels：连通区域的个数
%
% example : a = false(400);
% a(rand(400) > 0.5) = true;
% [b , l] = MyBwLabel(a);
% imshow(b , [])
if ~islogical(inputImg)
error( '========do not support the data type!!============' )
end
labels = 1;
outImg = double( inputImg );
markFlag = false( size(inputImg) );
[height , width] = size( inputImg );
%% first pass
for ii = 1:height
for jj = 1:width
if inputImg(ii,jj) > 0 % 若是前景点，则进行统计处理
neighbors = []; % 记录符合要求的邻域中的前景点，三列依次是行、列、值
if (ii-1 > 0)
if (jj-1 > 0 && inputImg(ii-1,jj-1) > 0)
neighbors = [neighbors ; ii-1 , jj-1 , outImg(ii-1,jj-1)];
end
if inputImg(ii-1,jj) > 0
neighbors = [neighbors ; ii-1 , jj , outImg(ii-1,jj)];
end
elseif (jj-1) > 0 && inputImg(ii,jj-1) > 0
neighbors = [neighbors ; ii , jj-1 , outImg(ii,jj-1)];
end
if isempty(neighbors)
labels = labels + 1;
outImg(ii , jj) = labels;
else
outImg(ii ,jj) = min(neighbors(:,3));
end
end
end
end
%% second pass
[r , c] = find( outImg ~= 0 );
for ii = 1:length( r )
if r(ii)-1 > 0
up = r(ii)-1;
else
up = r(ii);
end
if r(ii)+1 <= height
down = r(ii)+1;
else
down = r(ii);
end
if c(ii)-1 > 0
left = c(ii)-1;
else
left = c(ii);
end
if c(ii)+1 <= width
right = c(ii)+1;
else
right = c(ii);
end
tmpM = outImg(up:down , left:right);
[r1 , c1] = find( tmpM ~= 0 );
if ~isempty(r1)
tmpM = tmpM(:);
tmpM( tmpM == 0 ) = [];
minV = min(tmpM);
tmpM( tmpM == minV ) = [];
for kk = 1:1:length(tmpM)
outImg( outImg == tmpM(kk) ) = minV;
end
end
end
u = unique(outImg);
for ii = 2:1:length(u)
outImg(outImg == u(ii)) = ii-1;
end
labels = length( u ) - 1;
end
（2）Seed-Filling算法
1）Seed-Filling算法原理
种子填充方法来源于计算机图形学，常用于对某个图形进行填充。思路：选取一个前景像素点作为种子，然后根据连通区域的两个基本条件（像素值相同、位置相邻）将与种子相邻的前景像素合并到同一个像素集合中，最后得到的该像素集合则为一个连通区域。
种子填充法的原理参考的是这篇博客
（1）扫描图像，直到当前像素点B(x,y) == 1：
a、将B(x,y)作为种子（像素位置），并赋予其一个label，然后将该种子相邻的所有前景像素都压入栈中；
b、弹出栈顶像素，赋予其相同的label，然后再将与该栈顶像素相邻的所有前景像素都压入栈中；
c、重复b步骤，直到栈为空；
此时，便找到了图像B中的一个连通区域，该区域内的像素值被标记为label；
（2）重复第（1）步，直到扫描结束；
扫描结束后，就可以得到图像B中所有的连通区域；
算法图解如下：
2）Seed-Filling算法实现
种子填充法的matlab代码如下：
clear;
clc;
img=imread('5.4元.jpg');
imgGray=rgb2gray(img);
figure;
imshow(imgGray);
level=graythresh(imgGray);
imBW=imbinarize(imgGray,level);%二值化
figure;
imshow(imBW);
se=strel('disk',6);
imBW = imopen(imBW,se); %开运算
figure('NumberTitle', 'off', 'Name', '二值图像');
imshow(imBW);
outImg=uint8(imBW);
label=1; %标记
flag=0; %看是上下左右哪一个
num=0; %用于计算连通域的个数
[Height,Width]=size(imBW);
arr = [];%模拟堆栈
for m = 2:(Height-1)
for n = 2:(Width-1)
x=m;
y=n;
if(outImg(x,y)==0)
label=label+1;
arr(end+1)=x;%记录坐标
arr(end+1)=y;
num=num+1;
while(~isempty(arr))
x=arr(end-1);
y=arr(end);
arr(end)=[];
arr(end)=[];
outImg(x,y)=label;
if(outImg(x,y-1)==0)
arr(end+1)=x;%记录坐标
arr(end+1)=y-1;
% flag=1;%左
end
if(outImg(x,y+1)==0)
arr(end+1)=x;%记录坐标
arr(end+1)=y+1;
% flag=2;%右
end
if(outImg(x+1,y)==0)
arr(end+1)=x+1;%记录坐标
arr(end+1)=y;
% flag=3;%下
end
if(outImg(x-1,y)==0)
arr(end+1)=x-1;%记录坐标
arr(end+1)=y;
% flag=4;%上
end
end
end
end
end
figure('NumberTitle', 'off', 'Name', 'outImg');
imshow(outImg);
%%着色
R=img(:,:,1); %red
G=img(:,:,2); %green
B=img(:,:,3); %blue
for m=2:label
[O,P]=find(outImg==m);
for q=1:size(O,1)
R(O(q),P(q))=255;
G(O(q),P(q))=0;
B(O(q),P(q))=0;
end
end
img(:,:,1)=R; %red
img(:,:,2)=G; %green
img(:,:,3)=B; %blue
figure('NumberTitle', 'off', 'Name', 'RGB');
imshow(img);
（3）两种算法实现结果
可以看出两种算法和matlab自带的连通域标记函数bwlabel函数标记的结果相同。
3万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
