基于改进区域生长算法的图像分割方法及实现

一、题目与要求

1.1课程设计题目
采用灰度图像（该图像由高分二号卫星拍摄，是武汉市江夏区矿区的部分灰度图像）分别实现几种图像分割方法：
（1） 经典区域生长方法；
（2） 改进的区域生长方法；
（3） 分水岭方法。
1.2题目具体要求
（1）首先，基于经典区域生长方法，实现图像分割；
（2）然后，基于经典区域生长方法，自己查阅文献，设计改进的区域生长算法，实现图像分割；
（3）目视比较分析2种方法的分割效果，概述改进方法分割效果的优点和不足之处；
（4）实现分水岭分割。

二、需求分析

2.1问题分析
图像分割(image segmentation)是按照某些特性（如灰度值、频谱、纹理等）将图像划分成一些区域，在这些区域内其特性是相同的，或者说是均匀的，两个相邻区域彼此特性则是不同的，其间存在着边缘或边界。图像分割从本质上来说是将图象按照区域内的一致性和区域内的不一致性进行分类的过程，可视为一种特殊的图像识别问题。
图像分割主要分为以下几种方法：
（1）特征空间聚类：包括阈值化、尺度空间聚类、均值迁移等；
（2）边缘检测：包括梯度法、Snake模型等；
（3）区域跟踪：包括区域生长法、分裂合并法等；
（4）全局优化：包括Bayes方法、最小描述长度、最小割算法等；
（5）混合方法：包括形态学方法、表面拟合法、边缘-区域法等；
（6）基于识别的方法。
本次研究区域生长（region growing）方法。所谓区域生长法就是从满足检测准则的点（或者已知点）开始，在各个方向上把具有相似性质的像素集合起来构成一个区域。即在每个分割区域找一个种子像素作为生长的起点，再将种子像素周围邻域中与种子像素有相同或相似性质的像素（根据某种事先确定的生长或相似准则判定）合并到种子像素所在区域中。将这些新像素当做新的种子像素继续进行上面的过程，直到再没有满足条件的像素可被包括进来，这样的区域就生长成了。
综上，要解决该问题，通过区域生长方法实现图像分割，即要通过编程确定种子点、生长准则和停止条件。
2.2运行环境
我是在matlab（R2018a版）下完成的。
MATLAB是用于数值计算的软件产品，主要包括以下几个部分：MATLAB程序设计语言、MATLAB工作环境、MATLAB工具箱、MATLAB的应用程序接口等。MATLAB语言既是一种编程环境,又是一种程序设计语言。它的规则更接近于数学表示,因此集成性更好,应用更为简便。MATLAB超强的矩阵计算能力和图像处理工具箱相结合,使得MATLAB在数字图形、图像的处理领域获得广泛的应用。主要图像处理应用领域有：亮度变换与空间滤波、频域处理、图像复原、图像小波处理、图像压缩、形态学图像处理、图像分割、图像表示与描述、对象识别等。
基于此，matlab是能方便完成目标的不错选择。
2.3预处理
观察原图像，可发现图像上较均匀分散有噪声，影响观察及后续操作，故先对图像进行预处理，使用中值滤波去噪。中值滤波能达到噪声被消除，而边缘还完好保留的效果。

A=imread('jx_3.jpg');
imshow(A);
title('原图像');
B=rgb2gray(A);
C=medfilt2(B,[3,3]);      %利用中值滤波消除噪声
figure,imshow(C);
title('去噪后图像');

三、经典区域生长方法

3.1理论基础
区域生长算法是从单个像素出发，按照一定的生长准则，将具有相似性质的像素合并起来构成同一区域。具体做法如下：
(1) 选出图像中要分割的目标物体内的一点或者说生长起始点;
(2) 在种子区域的基础上不断将其周围的像素点以一定的规则加入其中;
(3) 最终将代表该物体的所有像素点结合成一个区域。
在所有图像中，同一个子区域内部像素会有相似的特性，不同子区域之间像素特性差异一般较大。区域生长算法从种子点开始生长，直到所有符合条件的像素点加入为止。区域生长算法原理很简单，其中最为核心的是生长准则，一般可以表示为一个灰度区间。如果邻域中像素点的灰度值处于这个区间中则被接纳，否则就被拒绝。
整幅图像中所有的像素点分为两类：己检测过的像素点和未被检测过的像素点。首先选择目标区域中一个或多个点(种子区域)加入到已生长区域，并以此作为生长起点，计算已生长区域中所有像素点灰度值平均值和标准差，接着判断目前区域的邻域中是否有符合生长准则的像素点，如果存在就将其划分入已生长区域，这样就完成一次迭代。第一次迭代完成后，重复上述步骤，直到没有满足条件的邻域像素点划分入己生长区域为止，区域生长算法结束。
区域生长算法程序流程图如下：

3.2算法设计
3.2.1鼠标选取单种子点

I=double(C);              %转换为灰度值是0-1的双精度
[M,N]=size(I);
J=zeros(M,N);             %定义输出图像矩阵，初始为零矩阵
msgbox('请用鼠标单击选择一种子点并按回车键', '操作','modal');
[y,x]=getpts;             %通过鼠标获得区域生长起始点
xr=round(x); 
yr=round(y);
seed=I(xr,yr);
J(xr,yr)=1;               %将J中与所取点相对应位置的点设置为白点
graysum=seed;             %储存符合区域生长条件的点的灰度值的总和
pointsum=1;               %储存符合区域生长条件的点的总个数
count=1;                  %每次判断一点周围八点符合条件的新点的数目
threshold=13;             %阈值
while count>0             %判断是否有新的符合生长条件的点，若没有，则结束
 s=0;                     %记录判断一点周围八点时，符合条件的新点的灰度值之和
 count=0;
 for i=1:M
     for j=1:N
       if J(i,j)==1
        if (i-1)>0 && (i+1)<(M+1) && (j-1)>0 && (j+1)<(N+1)       %判断此点是否为图像边界上的点
         for u= -1:1                                              %判断点周围八点是否符合阈值条件
          for v= -1:1                                             %u,v为偏移量
            if  J(i+u,j+v)==0 && abs(I(i+u,j+v)-seed)<=threshold  %判断是否未存在于输出矩阵J，并且为符合阈值条件的点
               J(i+u,j+v)=1;                                      %将符合以上两条件的点在J中与之位置对应的点设置为白点
               count=count+1;   
               s=s+I(i+u,j+v);
            end
          end  
         end
        end
       end
     end
 end
 pointsum=pointsum+count;
 graysum=graysum+s;
 seed=graysum/pointsum;                                     %计算新的灰度平均值
end
figure,imshow(J);
title('分割后图像');

3.2.2鼠标选取多种子点

I=double(C);              %转换为灰度值是0-1的双精度
[M,N]=size(I);
J=zeros(M,N);             %定义输出图像矩阵，初始为零矩阵
msgbox('请用鼠标单击选择多个种子点并按回车键', '操作','modal');
[y,x]=getpts;             %通过鼠标获得区域生长起始点(getpts函数本身就可以获得多个点，并返回多个点的坐标)
xr=round(x); 
yr=round(y);
threshold=13;             %阈值
T=length(xr);
for k=1:T                 %循环，依次以当选取的各点为种子点进行区域分割
 seed=I(xr(k),yr(k));
 J(xr(k),yr(k))=1;         %将J中与所取点相对应位置的点设置为白点
 graysum=seed;             %储存符合区域生长条件的点的灰度值的总和
 pointsum=1;               %储存符合区域生长条件的点的总个数
 count=1;                  %每次判断一点周围八点符合条件的新点的数目   
 while count>0             %判断是否有新的符合生长条件的点，若没有，则结束
  s=0;                     %记录判断一点周围八点时，符合条件的新点的灰度值之和
  count=0;
  for i=1:M
      for j=1:N
        if J(i,j)==1
         if (i-1)>0 && (i+1)<(M+1) && (j-1)>0 && (j+1)<(N+1)       %判断此点是否为图像边界上的点
          for u= -1:1                                              %判断点周围八点是否符合阈值条件
           for v= -1:1                                             %u,v为偏移量
             if  J(i+u,j+v)==0 && abs(I(i+u,j+v)-seed)<=threshold  %判断是否未存在于输出矩阵J，并且为符合阈值条件的点
                J(i+u,j+v)=1;                                       %将符合以上两条件的点在J中与之位置对应的点设置为白点
                count=count+1;   
                s=s+I(i+u,j+v);
             end
           end  
          end
         end
        end
      end
  end
  pointsum=pointsum+count;
  graysum=graysum+s;
  seed=graysum/pointsum;                                     %计算新的灰度平均值
 end
end
figure,imshow(J);
title('分割后图像');

3.3调试分析
最为经典的区域生长方法是单个种子点的，运行该程序，用鼠标选择种子点。

由以上两个测试例可以看出，区域分割效果较好，但结果与种子点的选取位置有关，同时不连通的区域无法被一起分割出来。
以上两测试例中，我们给定阈值13，下面将测试不同阈值对区域生长效果的影响。

由上面三个同一种子点、不同阈值的测试例（分别取10、20、30）可以看出，阈值的选取对图像分割效果有巨大影响。阈值选取过小会造成分割不完全，得到的图像不能将预期部分完整展现，而阈值选取过大会造成分割过度，得到的图像中包括不期望部分。故在手动选取阈值时，需多次尝试，得到最符合预期的分割情况。
在单个种子点的图形分割中，不连通的部分无法分割出，这时就需用多个种子点。（以下测试均在阈值恒定为13的情况下进行）

上述对多种子区域生长的测试例，选取不同种子点，得到的分割结果存在差异。同时，可以看出，种子点位置选取恰当，即可分割出不连通区域。从而，选取多个种子点，可实现将原图像中预期的多部分同时分割出来的效果。

四、区域生长改进方法

4.1改进思路
从上部分对经典区域生长方法的测试中，可以看出该方法存在的一个较大问题，是通过鼠标手动选取种子点，较为繁琐，且起始点的位置选择随机，完全依靠操作人员判断，图像复现性差。理想的图像分割是输入图像后，无需人手动选择，可程序自动确定种子点，得到分割效果较好的一个结果。基于此，对经典方法进行改进。
观察原图像，可发现预期分割出来的部分灰度低、亮度大，对应的图像像素值就要普遍较大。针对图像的这一特征，我们就可以用一个 3x3矩阵去遍历整个图片，然后计算矩阵内像素点的平均值，取平均值最大的区域中间值坐标点为生长起始点进行区域生长。便可实现程序自动选择较理想种子点的目的，称此方法为局部平均值判定法。
4.2算法设计

I=double(C);              %转换为灰度值是0-1的双精度
[M,N]=size(I);
J=zeros(M,N);             %定义输出图像矩阵，初始为零矩阵
for k=1:6					 %一次循环即为一个种子点
x=0;
for i=2:3:(M-mod(M,3)-1)                   %用 3x3 矩阵将整幅图像遍历 
    for j=2:3:(N-mod(N,3)-1) 
        Q=0; 
        for u=-1:1:1 
            for v=-1:1:1                   %计算每一个小矩阵的均值 
                Q=Q+I(i+u,j+v)./9;          
            end 
        end
        if Q>x&&J(i,j)==0                  %将均值最大且未存在于输出矩阵J的矩阵中心点作为生长起始点 
            x=Q;  
            xr=i;                              
            yr=j;                              
        end 
    end 
end 
seed=I(xr,yr);
J(xr,yr)=1;               %将J中与所取点相对应位置的点设置为白点
graysum=seed;             %储存符合区域生长条件的点的灰度值的总和
pointsum=1;               %储存符合区域生长条件的点的总个数
count=1;                  %每次判断一点周围八点符合条件的新点的数目
threshold=13;             %阈值
while count>0             %判断是否有新的符合生长条件的点，若没有，则结束
 s=0;                     %记录判断一点周围八点时，符合条件的新点的灰度值之和
 count=0;
 for i=1:M
     for j=1:N
       if J(i,j)==1
        if (i-1)>0 && (i+1)<(M+1) && (j-1)>0 && (j+1)<(N+1)       %判断此点是否为图像边界上的点
         for u= -1:1                                              %判断点周围八点是否符合阈值条件
          for v= -1:1                                             %u,v为偏移量
            if  J(i+u,j+v)==0 && abs(I(i+u,j+v)-seed)<=threshold  %判断是否未存在于输出矩阵J，并且为符合阈值条件的点
               J(i+u,j+v)=1;                                      %将符合以上两条件的点在J中与之位置对应的点设置为白点
               count=count+1;   
               s=s+I(i+u,j+v);
            end
          end  
         end
        end
       end
     end
 end
 pointsum=pointsum+count;
 graysum=graysum+s;
 seed=graysum/pointsum;                                     %计算新的灰度平均值
end
end
figure,imshow(J);
title('分割后图像');

4.3调试分析
在阈值给定为13的情况下进行测试，调整循环次数（一次循环即为一个种子点）。
1次循环：

2次循环：

3次循环：

4次循环：

5次循环：

6次循环：

4.4优点与不足
由上述测试例可以看出，该改进的区域生长方法较好地实现了区域分割，且能同时分割出不连通的几部分。执行一次循环，即为程序自动确定一种子点并进行分割，通过控制循环次数，即可控制种子点数。程序自动确定种子点，使得阈值和循环次数确定的情况下，图像分割效果一定，结果可复现。从上述1次循环递增至6次循环，可以看出循环次数较少时，可能会分割不完全，得到的图像不能将预期部分完整展现，而循环次数较多会造成分割过度，得到的图像中包括不期望部分。故需多次尝试，确定理想的循环次数，得到最符合预期的分割情况。此改进后的区域生长方法仍存在一较大不足，无法实现阈值自适应，还需人工给定阈值。下一步可从此处着手，寻找解决方案，实现程序自动确定阈值，获得更好的效果。

五、分水岭方法

5.1理论基础
传统分水岭算法（watershed）的基本思想是将图像看做为地理学的拓扑地貌，将图像想象成一个立体的地形表面。图像中的每一个像素的灰度值作为该点的海拔高度，每一个局部极小值及它影响的区域被称为集水盆，则集水盆的边界就形成了分水岭。分水岭的概念和形成能够通过模拟浸入的过程来说明，假设在每个区域极小值的位置刺穿一个小孔，从孔中向该区域不断注水，则地势较低的区域会首先被淹没，在浸入加深的过程中，每一个局部极小值的影响域都处在慢慢扩展过程中，当水淹没到一定高度时两个集水盆会相连合并为同一个区域，此时应在两个集水盆汇合处构筑“大坝”来阻止区域合并，最终形成分水岭，这些“大坝”的边界对应于分水岭的分割线，同时也是由分水岭算法提取出的连续边界轮廓。要实现分水岭方法，可分为以下几步进行：读入图像，求取图像的边界，应用分割算法。

5.2算法设计

I = C;
hy = fspecial('sobel');   %sobel算子
hx = hy';
Iy = imfilter(double(I), hy, 'replicate');   %滤波求y方向边缘
Ix = imfilter(double(I), hx, 'replicate');   %滤波求x方向边缘
gradmag = sqrt(Ix.^2 + Iy.^2);               %求摸
figure,imshow(gradmag,[]),                   %显示梯度
title('Gradient magnitude (gradmag)')
L = watershed(gradmag);                      %直接应用分水岭算法
Lrgb = label2rgb(L);                         %转化为彩色图像
figure,imshow(Lrgb),                         %显示分割后的图像
title('Watershed transform of gradient magnitude (Lrgb)')

5.3调试分析

通过测试可以看出，传统的分水岭方法对该图像分割效果较差，过分割现象明显。传统分水岭算法虽对微弱的边缘具有良好的响应，但对图像中噪声，物体表面细微的灰度变化比较敏感，噪声会恶化图像梯度进而造成分割结果轮廓偏移，通常的分割结果会出现严重的“过分割”现象。因此要想通过分水岭方法得到更好的图像分割效果，必须对其算法进行改进。
PS：暂未进行对传统分水岭方法的改进~

参考文献
【1】《数字图像处理与分析基础——MATLAB和VC++实现》 孙明主编 电子工业出版社
【2】李政国.基于区域生长法的高空间分辨率遥感图像分割与实现[D]．南宁：广西大学，2008．
【3】曹彪．基于区域生长的OTC 图像分割算法研究[D]．北京：北京理工大学，2014．
【4】许宪法.基于显著性检测模型的图像分割算法研究[D]．兰州：兰州理工大学，2018．
【5】闫东阳，明东萍.基于自动多种子区域生长的遥感影像面向对象分割方法[J].工程科学学报，2017.
【6】施兢业.基于改进区域生长法的电力设备红外图像分割[J].光学技术,2017.
【7】李云红,张秋铭,周小计,贾凯莉.基于形态学及区域合并的分水岭图像分割算法[J].计算机工程与应用,2019.