参考文献
第8章 基于分水岭分割的粘连颗粒图像分析技术研究
粘连颗粒图像无损处理与分析在医疗卫生、工农业生产等领域均有广泛的应用。它是将物质颗粒从图像中分割出来，统计其数目并提取各单个物质颗粒的特征参数(如面积、周长、直径、中心矩、颜色等)，对其品质进行分析。早先这些工作主要由人工观测来完成，工序繁杂、效率低，而且带有一定的主观性，误差大，重现性不好。采用图像处理与分析的方法能够将操作人员从这一繁重的工作中解脱出来，并大大提高计数与分析的精度，从而得到广泛的应用，成为近年来国内外的一大研究热点。
粘连颗粒图像分析的重点在颗粒图像的预处理、图像分割和颗粒计数、特征参数测量。首先对图像进行各种预处理，先将图像转换为灰度图，除去图像中不需要的彩色颜色信息；利用中值滤波算法对图像进行平滑处理，除去各类噪声的干扰；利用最大类间方差法选取合适的阈值，对图像进行二值化，使颗粒从背景中分割出来；另外设计了一种数学形态学滤波器，滤去了图像中颗粒不规则的边缘，使得互相接触部分减少，填充颗粒中的孔洞。接着提出了一种改进的分水岭分割算法，通过距离变换将二值图像转换为新的灰度图，消除了原始灰度图像中大量存在的“伪”局部极小值点，对灰度图进行分水岭分割算法，使得互相粘连的颗粒很好地分割开来，在后期对少量过分割产生的小面积区域，采用面积阈值法滤去面积小于阈值的区域，使得分割的效果更好。大量的实验结果表明该算法可以适用于不同区域的粘连颗粒图像，具有良好的适应性与鲁棒性，处理速度较快。最后提出了两种颗粒自动计数的方法，并利用颗粒的大小、形状、颜色等特征参数来描绘颗粒的特征，提取出各项特征参数的值。我们设计的图像分割与分析算法可以快速准确地将粘连颗粒物体分割开来，颗粒计数准确度达到99%以上，特征参数提取结果误差小。整个算法处理速度快而准确，具有良好的实时性与适应性，为颗粒的识别与归类打下了良好的基础。
8.1 颗粒图像分析基础与预处理
8.1.1 图像处理的基本概念
数字图像处理就是利用计算机技术对现实生活中获取的数字图像进行处理，提取各目标的特征加以分析，最后通过执行与人体视觉相关的识别方式，对不同目标进行相应的自动操作。视觉是人类感觉中最高级的，因此数字图像是计算机获取外界信息的最佳途径。在颗粒图像处理过程中，我们可以通过摄像机、相机、扫描仪来获取颗粒图像，但这些图像会受到环境因素的影响，如光照强度不均匀、元器件特性不稳定、其他颗粒物体干扰等使得采集的图像存在噪声，影响到后续处理。因此，需要对采集的原始图像进行一定的预处理来消除噪声、平滑图像，使得后续处理能更加准确。在这一章，首先对粘连颗粒图像分割处理过程要涉及的一些基本概念作简单的介绍，然后，对图像预处理的原理和具体过程作介绍。
8.1.2 数字图像与像素的表示
一幅图像可以被定义为一个二维函数f(x，y)，x 和y 是平面坐标，而f(x，y)是指在坐标点(x，y)处图像在该点的亮度。图像中(x，y)坐标和f(x，y)亮度是连续的，要将这样一幅图像进行数字化，就需要将坐标值通过取样方法来数字化，将亮度值通过量化方法来数字化。当亮度和坐标值都是有限且离散的量时，就可以称该图像为数字图像。
一副大小为M×N 的数字图像可以用如下的二维数组表示：
该数组中的每一个元素都称为像素或像元，像素的值表示该像素的颜色度。
8.1.3 图像类型
计算机上显示的图像有二值图像、灰度图像、伪彩色图像及真彩色图像等不同格式类型。下面就本节图像处理算法所涉及的三种类型作简单介绍。
(1)彩色图像(Color Image)
彩色图像的表示方法有许多种，最常用的是以RGB(Red/Green/Blue，红色、绿色、蓝色)三基色来表示的RGB彩色图像。RGB颜色空间中含有红色、绿色、蓝色三基色组合成的256×256×256种不同颜色，其模型如图8.1所示。
图8.1 RGB彩色空间
RGB彩色图像都是一个彩色像素的M×N×3三维数组。每个像素点都含有颜色空间中红色、绿色、蓝色三个基色的分量信息，每个分量的颜色信息从亮到暗可以用8个二进制位来表示，取值范围为[0，255]。彩色图像的特点是包含的信息比较多，如物体的色彩、纹理、形状、轮廓等信息。但需要更多的磁盘空间和系统内存，处理速度慢。通常在不需要太多颜色信息的时候会将彩色图像先简化为灰度图像后再处理。
(2)灰度图像(Gray Scale Image)
灰度图像是数字图像的基本形式，只表达图像的黑白亮度信息，包含白色、黑色和灰色部分。每个像素用8个二进制位来表示，取值范围为[0，255]，可表示256个灰度级。考虑到人眼可以分辨20～60级的灰度级别，因此这种表示方法保证人眼的分辨能力的同时符合计算机的数据类型。灰度图像可以表示为一个M×N 的一维数组，具有容易处理并且处理速度快的特点。
(3)二值图像(Binary Image)
二值图像是每个像素的亮度值取值仅为0或1的图像，0表示黑，1表示白，因此也称为1-bit图像。每个像素点仅用2个二进制位来表示，因此二值图像只能反映出图像中物体的基本形状。在图像处理中，二值图像由于其简单、处理速度快等特点而经常被运用到。
本设计中最开始采集到的图像是RGB彩色图像，由于后期处理中灰度图像已经满足图像分析所需要的信息要求，为提高图像处理效率，第一步需要将彩色图像利用转换公式先转换为灰度图像，使得图像只含有8位灰度信息，方便后续处理。具体的转化方法有如下两种：
方法1：由于RGB彩色图像是由R、G、B颜色分量组合而成，而三种颜色分量都以8个位表示，可以直接将分量数据转化为对应的灰度图像。其结果如图8.2所示，玉米颗粒的颜色可看作红色、绿色和蓝色三基色叠加而成。
图8.2 彩色图像转化为3个分量的颜色
(a)玉米图像；(b)R分量灰度图像；
(c)G分量灰度图像；(d)B分量灰度图像
方法2：提取出图像中的每个像素，利用矩阵的分布提取其R、G、B颜色分量，然后利用公式计算灰度图像中对应像素点的灰度值：
实验结果如图8.3所示。
图8.3 灰度图
比较两种转化方法，方法1的图8.2(b)、图8.2(c)和图8.2(d)3幅图分别对应三个分量，观察结果可以看出R分量转化图中黑色背景和灰色的目标(玉米粒)反差最大，效果最好；而B分量图像中目标接近黑色，与黑色背景最相似，该灰度图区分效果最差。造成该结果的原因是黄色的玉米粒在R分量上的灰度级比较高，使得转化为灰度图后亮度比较大。方法2则是结合三个分量的颜色分布计算得到灰度值，结果图中可以很好地将黑色背景与目标(玉米粒)区分开来。
分析以上结果可知，若颗粒图像中，目标颗粒的颜色与背景相差不是很大并且比较复杂，而在R、G、B中某个分量的灰度值比较高，而其他分量的灰度值比较低，则可以利用方法1来处理，灰度级别最高的分量转化成的灰度图要比利用方法2得到的结果好。在实际应用中，使用方法1的前提是处理对象是固定的目标与背景，可以通过实验获得大量的先验知识，从而确定哪个分量的结果图可以作为最佳的转化灰度图。
方法2虽然在转化效果上要比方法1略差，但方法2可以适应不同类型的颗粒图像，即若图像中目标与背景的颜色发生变化，利用方法2也可以很好地转化出效果良好的灰度图。由于本节中的主要实验目标对象是黄色的玉米粒，与黑色的背景差异较大，方法2得到的结果图完全可以满足后续的处理需求。考虑到本设计的整个图像处理方法可以适用于不同类型的颗粒图像分析，即适用范围更广，因此采用方法2将彩色图像转化为灰度图像。
在本设计后续处理中，还会将图像利用阈值法转化为二值图像，具体的处理方法在8.2节作介绍。
8.1.4 图像分割的定义
图像分割是将一幅图像按照不同的特征细分为不同的区域或对象，并从复杂背景中提取出感兴趣的目标区域。细分的程度取决于实际应用中待解决的问题，只要满足需求就可以停止分割。非普通的图像分割是图像处理中最困难的问题之一，分割的准确性是进行图像分析的关键，直接影响着后续的特征提取的结果，从而影响更高层次对象操作的成功与否。图像分割可以定义如下：
假设集合A 代表图像中的所有像素，图像分割相当于将集合A 分成n个满足下面5个条件的子集合(区域)A1，A2，A3，…，An：
(1)
(2)集合A 的所有子集合中，若i≠j则有Ai∩Aj=φ；
(3)任意的i=1，2，3，…，n，P(Ai)=TURE；
(4)若i≠j，则P(Ai∪Aj)=FALSE；
(5)集合A 的所有子集合中的像素都是连通的区域。
以上定义的条件是从像素级来定义的，因此该定义适用于不同目标的图像分割。条件(1)表示图像中所有的像素都属于某一个子区域，所有子区域里的像素总和为整幅图像的像素。条件(2)表示一个像素只能属于各个子区域，不同子区域没有重复的像素。图像分割是根据应用需求按照不同的特征将像素划分为不同区域，因此条件(3)表示同一区域中的像素具有相同的特征，条件(4)表示不同区域的像素具有不同的特征。条件(5)表示分割出来的目标区域是一个连通的区域，只有连通区域在高层次的图像处理中才能代表有意义的目标对象。
本节研究的图像分割主要是针对颗粒物体出现粘连的情况，首先需要将图像中颗粒物体所在的目标区域从背景环境中分割出来，再将互相粘连的颗粒物体分割成单独的目标区域，为单独统计每一颗粒的特征参数打下良好的基础。
8.1.5 图像平滑处理
图像平滑是指消除或减少噪声和物体本身存在的颜色突变部分，改善图像的质量，使得后续图像分析有更好的结果。噪声产生的原因包括图像采集时光照的影响、图像保存或传输时引入的信道噪声、环境中一些颗粒物体的干扰。物体本身的颜色突变部分是指物体由于受损或质量问题存在的一些裂缝、空洞或斑点等，这些都会使图像处理产生干扰。消除这些干扰的常用方法有均值滤波法、中值滤波法。图像中的噪声往往和有用信息交织在一起，这两种滤波方法在使用过程中都会使图像清晰度降低。
(1)均值滤波法
均值滤波法是一种线性滤波方法，它的原理是给各个目标像素一个模板，该模板由以目标像素为中心的周围相邻像素组成，利用模板中的所有像素的平均值来代替目标像素值。
假设f(x，y)是图像中坐标为(x，y)的像素点的颜色值，经过一个M×M 的模板均值滤波后，该像素点的值变为g(x，y)，变换公式为：
该方法的模板其实是一个低通滤波器，滤掉图像中的一些高频的噪声。但是图像中的细节信息主要分布在高频区域，因此图像滤波往往会导致图像的细节变得模糊起来。滤波的能力与模板大小选取有关，模板越大则去噪能力越强，但图像会更模糊，模板越小则滤波能力越小。因此滤波方法需要对消除噪声干扰和保留图像细节两方面进行平衡。
(2)中值滤波法
中值滤波法是一种常用的非线性的平滑滤波处理方法，其原理与均值滤波法相似，不同点在于输出的像素是由邻域像素的中间值而不是均值得到的，适合于消除孤立的噪声，较好地保护图像的细节。假设f(x，y)是图像中坐标为(x，y)的像素点的颜色值，经过一个M×M 的模板均值滤波后，该像素点的值变为g(x，y)，变换公式为：
其中Med{}代表取像素的中间值。中值滤波法的具体步骤为：
①将模板在途中遍历，模板中心与图像中的像素点重合，读取模板下各对应像素的值。
②对读取的各个像素值进行排序，计算其中间值。
③将中间值赋给图像中该像素点作为此像素的颜色值。
④按上述操作遍历所有的像素点。
在实际应用中数字图像中存在的噪声可以分为椒盐噪声与高斯噪声，其中椒盐噪声是由图像传感器、传输信道、解码处理等产生的黑白相间的亮暗点噪声，而高斯噪声是指由于磁场、电场发生无序的变化，经过放大电路产生的各种干扰噪声，其概率密度函数服从正态分布。图8.4为利用5×5的模板对两类噪声干扰图像分别进行均值滤波与中值滤波平滑预处理的结果。
分析以上实验结果可知，均值滤波与中值滤波都可以滤去图像中的噪声，起到平滑图像的作用。对于椒盐噪声干扰的图像，由于其噪声点相对比较少，利用中值滤波处理图像，在消除噪声点的同时还能保持原始图像的边缘部分保持清晰；对于高斯噪声的图像，由于噪声点比较多，利用均值滤波可以比中值滤波更好地消除大部分的噪声影响，所得到的图像模糊程度与滤波时采用的模板有关，模板越大，平滑效果越好，图像也会变得更加模糊。
对于本节研究粘连颗粒图像处理所采用的玉米颗粒图像，图像本身会在传输编码过程中存在一些微小的噪声干扰，目标的边缘位置会出现一些容易与背景混淆的像素点，灰色目标内部也可能会因为细节部分而出现一些突变的点，黑色背景像素点会由于光照过强而反光的原因被误判为目标点。同时为了避免在采集图像的过程中空气中灰尘等小颗粒物体的干扰，需要对图像进行滤波方法的平滑处理，以使后续的分割效果更好。由于图像中出现的噪声点总数不多，可以归为椒盐噪声，因此采用中值滤波会取得更好的平滑效果。图像平滑结果见图8.5。
图8.4 滤波效果图
(a)原始图像；(b)椒盐噪声干扰图像；(c)高斯噪声干扰图像；
(d)椒盐噪声图像均值滤波；(e)高斯噪声图像均值滤波；
(f)椒盐噪声图像中值滤波；(g)高斯噪声图像中值滤波
图像平滑可以使得图像中微小的干扰消失，如图8.5中大的圆圈内的背景像素点受光照影响灰度接近目标，经平滑处理后亮的像素点减少。小的圆圈内的玉米颗粒边缘则存在一些干扰的像素点，颗粒内部的灰度值也不均匀，这些干扰都会对颗粒的特征分析以及颗粒计算造成一定的影响。经过平滑后图像目标颗粒的灰度值变得较为统一，不会出现极小区域的突变点。平滑处理使得后续的图像二值化处理区分目标与背景像素点变得更加准确。
图8.5 平滑效果图
(a)玉米灰度图像；(b)平滑结果
8.1.6 图像二值化
图像二值化是指将灰度图像中的像素点的灰度值设为255或0，使得所有像素点只有白色与黑色两种颜色。二值图像化依然可以反映出图像的整体特征与局部特征，在数字图像处理中有十分重要的地位，因为图像中像素值只有0与255两个值，在计算机中只需要1个位的0值和1值表示，不需要涉及多级的灰度，使得二值图像变得十分简单，处理速度快，并且数据量小，占用内存也小。在实际应用中，一般采用阈值法将灰度图像变换为二值图像。阈值法将图像中灰度特征大于阈值的像素集合划分为目标物体，小于阈值的像素集合被划分为没有意义的背景像素。如最简单的划分方式是选择特定的阈值T，按如下公式将灰度图像f(x，y)转化为二值图像g(x，y)：
常用的阈值法有很多，如平均值阈值法、最大类间方差法、迭代方法选取阈值、最大熵阈值法。阈值T 的选取方法会直接影响像素的判定结果，所以利用阈值法将灰度图像变换为二值图像关键是选取合适的阈值。以下对常用的阈值算法作简单介绍。
(1)平均值阈值法
平均值阈值法是指将图像中所有像素的灰度值的平均值作为阈值来划分图像，该方法思想简单，运算速度快，适合于目标物体与背景颜色反差较大的图像。其中一幅M×N 的图像f(x，y)的阈值T 的计算公式如下：
(2)最大类间方差法
最大类间方差法原理假设用阈值T 将图像划分为目标和背景两个类，计算出两个类的方差，当方差最大时，表示目标和背景之间差异最大，此时的T 值为最佳阈值。因为T 的取值范围为[0，255]，因此需要经过256次的方差计算，求出类间方差最大值对应的T 值。该方法是一种自动的非参数监督的阈值选取方法，耗时少，适用于灰度分布直方图中背景与目标没有明显的双峰，类间方差呈单峰分布的图像，划分效果比较理想。
假设一幅图像其灰度级范围为i=0，1，2，…，L-1，当灰度级为k 的像素个数为nk时，则图像的总像素N 为：
灰度级为i的像素出现的概率为：
假设阈值t将图像像素分为背景和目标两个类，即A0={0，1，2，…，t}，A1={t+1，t+2，…，L-1}，两个类的像素出现的概率分别是：
两个类的均值是：
两个类总的平均值是：
最后求出两个类的方差：
分别对t=0，1，2，…，L-1计算两个类的方差，当方差最大时，对应的t就是最佳阈值T，即T 满足如下公式：
(3)迭代法选取阈值
迭代法是指运用逐次逼近的思想来选取最佳阈值，该方法能区分出图像的目标和背景的主要区域所在，但对图像的细微处没有很好的区分度，不适合目标和背景区域面积差异比较大的图像。该方法的具体步骤如下：
①首先设置一个初始阈值T0，可以将图像中所有像素的平均灰度值作为T0的值，即在一幅M×N 的图像f(x，y)中令：
②根据阈值Tk将图像中灰度级i≥Tk的像素划分为目标，灰度级i＜Tk的像素划分为背景，分别利用公式求出背景区域和目标区域的平均灰度值A0和A1：
③求出新的阈值Tk+1=(A0+A1)/2。
④若Tk+1≠Tk，转到②，继续迭代计算；若Tk+1=Tk则表明阈值不再发生变化，此时的阈值Tk+1为最佳阈值T。
(4)最大熵阈值法
最大熵阈值法是指结合信息论中的信息熵的概念与图像阈值化技术，使得选择的阈值T 分割图像目标和背景区域两部分灰度统计所得的信息量最大。该方法对不同目标大小和信噪比的图像能产生较好的分割效果，由于其中含有对数运算，因此算法效率比较低，对一般的图像而言，其分割效果和最大类间方差法相比要差一些。
假设阈值为t，图像中像素灰度i=0，1，2，…，L-1出现的概率为pi，其中灰度级i＜t的像素划分到背景区域A，灰度级i≥t的像素划分到目标B，则背景区域A 的概率灰度分布为：
目标区域B 的概率灰度分布为：
以上两公式中：
由熵的定义可知，图像中背景区域A 与目标区域B 的熵为：
由背景与目标区域的熵HA(t)和HB(t)可得整幅图像的信息熵φ(t)为：
对所有的灰度值t=0，1，2，…，L-1求信息熵φ(t)，当φ(t)最大时，对应的灰度值t即为最佳阈值：
8.1.7 粘连颗粒图像二值化的结果与分析
在颗粒图像二值化中可以用最大类间方差法、迭代法以及依据像素灰度直方图分布等方法来判断图像的分割阈值。不同的方法各具优缺点，其中依据灰度直方图分布判断阈值的方法效果最好，其结果如图8.6所示。
依据灰度直方图分布判断阈值需要人为地选择最佳阈值，如在直方图中，背景像素点几乎在灰度值40附近，而目标区域集中在灰度值120附近，可以选择灰度值80作为最佳阈值，将图像分为白色的目标像素点与黑色的背景像素点。该方法可以在大量实验的基础上充分利用预先统计的先验知识来求得最佳的阈值，但是当图像的目标或背景发生变化时，该方法的阈值由于已经设定的原因无法适应新的图像类型。因此实际应用中通常也会采用最大类间方差法与迭代法来求阈值，其对于同一幅图像的实验结果如图8.7所示。
图8.6 直方图选取阈值过程图
(a)细胞灰度图；(b)二值化结果图；(c)灰度分布直方图
图8.7 其他方法图像二值化效果图
(a)迭代法求阈值的二值化结果；(b)最大类间方差法二值化结果
通过计算得出迭代法求得的最佳阈值为99.8，最大类间方差法求得的最佳阈值为87.3。从二值化结果中可以看出，利用迭代法时有一部分目标内部像素点由于接近黑色而被划分为背景像素，该方法优点是计算量小、速度快。利用最大类间方差法能更好地将目标与背景像素区分开来，虽然目标内部含有少量的误判像素点，但其结果已经可以满足下一步的图像处理需求。
在进行二值化处理时，由于最大类间方差法得到的分割结果良好，而且适用性强，可以处理不同类型的颗粒图片，综合考虑图像二值化的处理效果与算法的鲁棒性等因素，采用最大类间方差法对玉米颗粒图像进行二值化处理，其结果见图8.8。
图8.8 玉米颗粒图像二值化结果
(a)玉米颗粒图像；(b)二值图像
观察结果可知图像二值化处理可以将原始图像先转化为灰度图像，再转化为黑白两色表示的二值图，去除了图像中复杂的灰度颜色、纹理、细节等信息，每个像素点只需用1位来表示，同时保留了图像中目标的形态、大小等信息。由于后续的图像处理是在二值图像上对目标颗粒的形状进行优化，只需要颗粒的形状信息，因此二值化处理后简单的二值图像为后续的图像快速处理打下了良好的基础。
8.2 图像数学形态学处理与距离变换
8.2.1 图像数学形态学处理
数学形态学诞生于1964年，当时法国G.Matheron和J.Serra在从事铁矿石的定量岩石学分析及预测开采过程中，第一次引入了数学形态学表达式，建立了颗粒的分析方法。图像数学形态学处理是以几何学为基础，着重研究图像的集合结构。其基本思想是利用一定形态的结构元素去度量和提取图像中目标的对应形态，来实现图像的分析和识别。如图8.9给出的二值图像A 与圆形结构元素S，结构元素从左往右、从上往下平移，其中一个位置可以完整地放入结构元素，另一个位置却无法放入结构元素。标记图像内适合放入结构元素的位置，就可以得到与目标结构有关的信息，当然结果受结构元素的大小和形状影响。
图8.9 数学形态学处理示意图
数学形态学的描述语言是集合论，可以运用集合的思想来解决图像处理应用中碰到的问题。它在实际应用中主要是利用不同的数学形态学基本运算(腐蚀、膨胀、开运算、闭运算)的合理组合，对图像进行处理，简化图像或改善图像的质量。与频域或空域的其他图像处理方法相比，数学形态学具有一些明显的优点，如：
(1)在图像恢复过程中，数学形态学的形状滤波器可参考已知的目标形状特征，不仅可以很好地滤掉噪声，图像的初始信息也可以被较完整恢复出来。
(2)可以并行处理，速度比较快。
(3)提取边缘信息比利用微分算法的边缘提取算法要好，对噪声干扰不敏感，而且提取的边缘比较平滑完整。数学形态学方法提取出来的图像目标的骨架也具有连续、断点少的优点。
在研究颗粒图像分割的过程中，利用数学形态学中的腐蚀与膨胀运算，可以很好地达到二值图像简化与改善图像的效果，为后续的分割打下良好的基础。
8.2.2 二值图像的腐蚀运算
二值图像是值只有0和1的图像，利用数学形态学运算处理图像的时候，运算只对值为1的区域进行处理。设一幅图像目标区域为集合A，结构元素(腐蚀因子)为S，令结构元素S 在图像上移动，则对于每个当前位置x，S+x 有三种状态：
(1)S+x⊆A；(2)S+x⊆A ；(3)S+x∩A 和S+x∩A都不为空，具体如图8.10所示。
图8.10 结构元素与目标对象关系图
状态1表示S+x 与A 相关最大，状态2表示S+x 与A 不相关，状态3表示S+x 与A 部分相关。定义图像中所有最大相关的点的集合为S 对A的腐蚀，记为AΘS，集合的表示方式为：
图8.11中A 是被处理的目标区域，S 是结构元素。由定义可知，目标区域经过腐蚀后剩余的点x 都满足S+x⊆A，剩余的阴影部分都在A 的范围之内，并且比A 要小，目标区域A 像被溶解掉一层，因此称这过程为腐蚀。腐蚀的直接作用是消除物体边沿的点。
图8.11 腐蚀示意图
在实际应用当中，腐蚀的结果和结构元素(腐蚀因子)的大小和形状选取有相当大的关系。如图8.12为三种不同大小和形状的结构元素数组表示方式：
图8.12 不同大小和形状的结构元素图
(a)3×3方形结构元素；(b)5×5菱形结构元素；(c)5×5圆形结构元素
利用3×3方形结构元素去腐蚀目标区域时，可以使目标区域沿边界减少一个像素，因为结构元素大小太小的缘故，腐蚀的效果并不明显，可以经过多次腐蚀来增加效果；若利用大一些的5×5圆形结构元素去腐蚀目标区域，一次腐蚀目标会沿边界减少两个像素。利用5×5菱形结构元素虽然一次也能沿边界减少两个像素，由于两个结构元素的形状不一样的原因，如图8.12中方框内表示圆形结构元素比菱形结构元素在四周要多两个值为1的点，在腐蚀过程中就会少腐蚀像素，腐蚀结果就会使目标区域边界呈圆滑性，而菱形结构元素则腐蚀性更好，腐蚀结果边界会呈现一些尖角。当结构元素的大小继续增大时，图像中一些小目标区域长宽小于结构元素，结果会被直接腐蚀掉。在颗粒图像处理中合理利用这些性质，可以有效地简化图像，使粘连颗粒分开一些，滤掉干扰小物体。
8.2.3 二值图像的膨胀运算
膨胀运算与腐蚀运算相反，它的作用是把目标区域沿边界周围的点归并到目标区域，使目标区域面积增大。在一幅图像中，移动结构元素S，对于图像的每个位置x，若S+x 与A 有交集，则所有x 位置的像素集合为S 对A 的膨胀，记为A ⊕S，集合的表示方式为：
图8.13中A 是被处理的目标区域，S 是结构元素。由定义可知，目标区域经过膨胀后的点x 都满足S+x ∪A ≠φ，阴影部分比原先的目标区域A 要大，就像被扩大一层，因此称这过程为膨胀。
图8.13 膨胀示意图
在实际应用中图像经过膨胀运算后，两个距离比较近的目标物体可能会连通在一起，因此可以利用膨胀运算来填补目标物体不光滑的凹陷部分。另外，根据目标物体中的孔洞直径，选择恰当大小的结构元素对图像进行膨胀，也具有一定的修补空洞的作用，但会导致目标物体增大。
8.2.4 开运算与闭运算
数学形态学处理中腐蚀和膨胀运算不是逆运算，不同组合会产生不同的效果。其中开运算是利用相同的结构元素对图像先腐蚀后膨胀的运算，它可以很好地保持图像中目标区域的基本外形，同时除去边界毛刺，使目标区域的凸角变得圆滑一些。利用结构元素S 对图像A 作开运算，可以用符号A。S 表示，其定义为：
闭运算就是指利用相同的结构元素对图像先膨胀后腐蚀的运算，它可以填平小的空洞，弥合小的缝隙，同时使得目标区域的凹角变得更圆滑一些。可以用符号A·S 表示闭运算，其定义为：
8.3 粘连颗粒图像的数学形态学处理
8.3.1 粘连图像数学形态学处理的结果与分析
由于图像中物体表面反光不均匀，本身的边缘不圆滑，外界光强不足，其他杂质混入等原因，即使图像经过平滑预处理，二值化处理后的黑白图像中仍然存在着一些缺陷，如图8.14所示，有些目标物体内部有像素点颜色接近背景而被误判为背景像素，使得物体内部出现孔洞[图8.14(a)颗粒图像用灰色表示，透明部分代表空洞]；有些物体的边缘会出现一些不圆滑的突出部分；小面积的杂质被判定为目标物体。这些缺陷都会严重影响后续的分水岭分割的结果，使得物体计数发生错误，同时在特征提取过程中，也会产生错误的面积、周长、质心位置等特征参数。
图8.14 平滑处理后二值图像
(a)原始图像；(b)二值图
针对上述缺陷，可以利用数学形态学对二值图像进行孔洞填充、开运算、多次腐蚀等一系列的滤波处理，以改善二值图像，使其能更准确地反映真实的图像信息，具体的方法如下：
(1)孔洞填充
孔洞填充的思想是首先在各个封闭区域中，选择一个点作为种子点并填充一新的颜色。再搜索该点的4邻域像素点，判断其颜色是否与种子点颜色相同，或者与边界点颜色相同。前者表明已判断过，后者表明为不需处理的点。否则将点作为孔洞点压入堆栈并填充种子点的颜色。重复进行以上处理直到所有点都判断过，最好再将所有孔洞点填充为目标的颜色。具体的处理流程图如图8.15所示。
图8.15 孔洞填充流程图
对图8.14(b)中的二值图像进行填充，其结果如图8.16所示。
图8.16 二值图像进行填充结果
由结果可知，孔洞填充可以将物体内部的黑色像素都填充为目标像素。但当三个或三个以上的物体互相粘连时，也会出现“伪”孔洞现象，如图8.17所示，此时对图像进行孔洞填充会使不同物体完全粘连在一起成为一个物体，由于粘连处的缝隙被完全填充，而失去凹性特征，因此后期分水岭分割也不能将其分割开来。多次实验过程表明，对于会出现多目标粘连而造成“伪”孔洞现象的图像，解决该问题只能通过增加光照强度，或者选用与目标颜色差异很大的背景色，使得在图像二值化过程中物体不会因为光强不足而出现孔洞。对于很多不会出现“伪”孔洞现象的图像，则加入孔洞填充这一图像处理过程，使得结果更精确。
图8.17 “伪”孔洞图像填充结果图
(a)“伪”孔洞图像；(b)填充效果图
(2)图像开运算处理
当图像中存在着一些面积微小的杂质物体并且目标物体边缘不圆滑时，会使后续图像分割发生错误，可利用数学形态学中的开运算对图像进行优化。选用5×5的圆形结构元素对图像作开运算，即首先对图像作腐蚀运算，消除图像中的杂质物体，同时腐蚀掉物体边缘突出部分；再对图像作相同结构元素的膨胀运算，使目标物体尽量恢复成原形状。对图作开运算得到的结果如图8.18所示，从结果图中可以看出，开运算在保证目标物体基本形状不变的情况下，可以在一定程度上消除物体边缘的突出部分，剔除杂质物体，使图像更加简洁，并且利于分割。
(3)图像腐蚀运算处理
当图像物体与物体过度粘连时，其粘连处的凹形特征不明显，如图8.18橙色圈内的粘连物体，在分割过程里中间接触部分会多分割出一个扁的长方形，分析阶段也会产生错误的结果。可以利用腐蚀因子对图像进行腐蚀运算，减少粘连处的像素点，使得凹形特征更加明显，有利于后续的分水岭分割。对一幅图像进行多次腐蚀，腐蚀的效果与腐蚀次数和结构元素的形状选取有关。如图8.19分别是利用1×1的圆形结构元素对图像作2次腐蚀运算与4次腐蚀运算。若物体之间的粘连部分较少，则两个物体会被直接分割开来；若粘连部分较多，则多次腐蚀使得粘连部分减少。由于腐蚀运算是针对所有的目标像素点，因此在腐蚀过程中，目标物体的大小与形状也会发生变化。
图8.18 开运算结果图
图8.19 不同次数的腐蚀结果图
(a)1×1的圆形结构元素2次腐蚀结果图；(b)1×1的圆形结构元素4次腐蚀结果图
腐蚀运算的结构因素形状也会对结果产生影响，如图8.20中分别是利用5×5的菱形结构与5×5的圆形结构进行腐蚀，可以看出同样大小的菱形结构元素可以更好地保持物体的圆滑性，而圆形结构元素能腐蚀更多的像素点。容易知道与腐蚀相反的膨胀运算也具有该类性质。
图8.20 不同形状结构因素的腐蚀结果图
(a)5×5菱形结构因素腐蚀结果图；(b)5×5圆形结构因素腐蚀结果图
由以上分析可知，腐蚀运算不仅使粘连物体趋向于分开，也会使物体本身的形状发生变化，这对后期的特征参数提取会造成严重的影响。因此需要选择合适的腐蚀结构元素，在保证物体基本形状不发生变化的情况下，使得物体粘连处的像素点尽量减少。
8.3.2 数学形态学滤波器的设计
结合上节中分析的各种数学形态学运算的优点与缺陷，本节提出一种新的数学形态学滤波器来对粘连颗粒图像进行优化。其示意图如图8.21所示。
图8.21 一种新的数学形态学滤波器示意图
由于数学形态学处理中的算法速度都比较快，因此在该滤波器中对二值图像进行两次开运算，可以更好地消除小的杂质物体，使物体边缘那些突出的不规则部分变得圆滑。再用1×1圆形结构元素进行闭运算，可以消除像素直径不大于2μm 小孔洞，并且可以保持“伪”孔洞的形状不变。考虑到圆形结构因素腐蚀能力强，菱形结构能更好地保持物体形状的特性，用4×4圆形结构因素对图像进行腐蚀，用3×3菱形结构因素在腐蚀后的图像上进行膨胀，并重复2次相同的处理。该方法比单独利用2×2的菱形或圆形结构元素直接对图像进行一次腐蚀的效果要好，更能保持物体的基本形状。图8.22是利用该滤波器对本节所采用的粘连玉米颗粒二值图像进行优化处理的结果。
图8.22 数学形态学滤波后的二值图像
(a)原二值图像；(b)滤波后的二值图像
从图像中可以看出数学形态学滤波器可以滤掉图像中的一些小颗粒物与物体边缘不圆滑的部分，同时减少颗粒物体粘连处互相接触的点。该优化结果让类圆性颗粒物体粘连处的凹形特征更加明显，可以在一定程度上抑制分水岭分割过程中的“伪”极小值，有效避免将互相粘连处单独分割成一个细长的长方形。
8.3.3 图像距离变换
(1)图像距离变换的定义
图像距离变换是指将二值图像中目标区域像素(一般值为1)的值转化为该点到达最近的背景点的距离，离区域边界越近的边界像素点距离值越小，离区域核心越近的像素点距离值越大，转化的结果是将二值图像转化为表示距离值的灰度图。1966年Rosenfeld等人提出距离变换的概念，虽然得到的距离只是近视的欧式距离，但为后来的各类快速距离变换奠定了坚实的基础，在那之后学者们提出了各类改进的距离变换算法来提高算法速度和准确性，使得距离变换在目标细化，图像合并、分割，目标重建等图像分析和识别领域得到了广泛的应用。尤其是在颗粒图像处理中，距离变换成功地使二值颗粒重建成灰度颗粒，方便了后续的分割处理。在实际应用中，距离变换可以描述为：一幅N×N 的二值图像A=[aij]中目标区域的像素点集合为M={(x，y)|axy=1}，背景像素点集合为B={(x，y)|axy=0}，距离变换就是利用公式求得目标区域像素到背景像素之间的最短距离，得到距离变换后的图像D[dij]。
在实际的应用中，欧氏距离变换是一种精确的L2范数距离变换，具有非线性。它可以被应用于不同领域的图像处理，但由于在变换过程中计算某个像素到本图像所有背景点的最短距离是一种全局操作，对每个点进行的平方和开放次数为图像的总点数。若图像的大小为100×100dpi，不作简化直接进行距离变换，需要的计算次数最大可达到(100×100)2=108级。由于欧氏距离的计算量大、速度慢、复杂度高，因此一般采用8邻域距离D8(棋盘距离)或4邻域距离D4(街区距离)等近似的欧氏距离来代替欧氏距离。
像素(x，y)和(i，j)之间的D4距离(街区距离)定义为：
所有距像素点(x，y)的D4距离为不大于d 的像素可以组成一个中心点在(x，y)的菱形。如与点(x，y)的D4距离为小于2或等于2的像素组成的菱形如图8.23所示。
(https://www.daowen.com)
图8.23 街区距离图
像素(x，y)和(i，j)之间的D8距离(棋盘距离)定义为：
所有距像素点(x，y)的D8距离为不大于d 的像素可以组成一个中心点在(x，y)的方形。如距点(x，y)的D8距离小于或等于2的像素组成的方形如图8.24所示。
图8.24 棋盘距离图
(2)图像距离变换的实现
图像距离变换的实现函数有很多，直接进行变换计算量过大、速度慢，因此在具体实现过程中会采用并行的方法来提高效率，或者将全局操作分解为局部的操作来简化算法，如用串行扫描的方法。串行扫描的方法是在扫描的过程中利用上一次计算的结果来计算当前像素的值，扫描过程由正向扫描和反向扫描两次扫描组成。具体的扫描方式如图8.25所示。
图8.25 扫描方式
(a)正向扫描；(b)反向扫描
在实际的二维图像中，采用一个类似卷积核的模板以类似卷积的操作在图像中进行扫描来实现图像的距离变换，其过程与一维图像的距离变换相似。首先对图像进行从左到右、从上到下的正向扫描，然后再对图像进行从右到左、从下往上的反向扫描。在扫描过程中模板覆盖的每个像素点要将其值与对应的像素值相加，从而形成一个两项和的集合。对于每个当前像素点即模板中心点对应的像素点，其值为两项和集合中的最小值。图8.26是二维图像距离变换的模板，其中3×3的模板扫描计算的结果是一灰度级图像，其像素灰度级为像素点到最近背景的欧氏距离的3倍，与真实欧氏距离的偏差不超过8%。5×5的模板扫描结果其像素灰度级为像素点到最近背景的欧氏距离的5倍，与真实欧氏距离的偏差不超过2%。由于距离变换的精度与模板的大小有关，模板越大，精度也会越高。因此对于一些精度要求较高的应用领域，很多学者都利用提高模板大小的方法来提高精度，现在已有13×13的模板。另一方面，由于模板越大，计算量也会随之增大，图像处理的效率就会变低；而随着模板的变大，精度提高的效果也会越来越低，因此有必要平衡模板的大小和实际应用中对精度的需求。对于后期的分水岭算法，我们选用5×5的模板就可以满足精度的需求。
图8.26 距离变换3×3的模板与5×5的模板
8.3.4 粘连颗粒图像的距离变换结果与分析
原始灰度图像中包含了许多纹理、细节信息，使得对灰度信息敏感的分水岭分割无法理想地将互相粘连的颗粒分割开来。距离变化将优化后的二值颗粒图像重新转化为灰度图，新产生的灰度图中的灰度值分布均匀，除去了物体的细节纹理等灰度信息，像素点的灰度值反映该点到背景像素的最近距离。
在实际应用中，棋盘距离具有计算简单、距离最小、忽略微小差距的优点，使所得的灰度图的灰度分布更加均匀。可以利用3×3的模板对图像进行两次扫描来实现棋盘距离变换，将颗粒二值图像转变为灰度图像，其结果如图8.27所示。图8.27(c)为二维矩阵表示的距离图，可以看出物体中间的像素点离背景最远的距离为6。将距离图转化为对应的0～255的灰度图如图8.27(b)所示，物体中间距离值最大的点，对应灰度值为255的亮点，而边缘的点则为灰度值接近0的暗点。
图8.27 二值图像距离变换
(a)原图像；(b)距离灰度图；(c)距离图像二维矩阵
对本设计所采用的玉米颗粒图像进行距离变换，其结果如图8.28所示。从图中可以看出，颗粒的灰度分布均匀，由边缘向中间逐级变亮，中间为最亮点，边沿最暗。粘连部分像素点的灰度值比较低。
虽然距离变换可以使灰度图像分布均匀，利于分水岭分割，但由于物体本身形状不具有类圆性或者物体互相粘连部分太多，导致分水岭分割将单独的颗粒物体一分为二，或两个粘连物体分为三部分。在后续的分水岭分割过程中会出现局部极小值的概念，一般局部极小值对应物体的中心点，其数目与物体数目相同。在这里可以通过对灰度图像按灰度级加轮廓的方法来更清楚地观察图像，查看其局部极小值，结果如图8.29所示。其中图8.29(b)中的一个物体由于形状不具类圆性而出现了两个局部极小值，会被从中间分割为两个物体。8.2节中的数学形态学优化可以有效地减少该类误差的出现。
图8.28 玉米颗粒图像距离变换
(a)玉米图像；(b)距离灰度图
图8.29 距离灰度图彩色表示
(a)距离灰度图彩色化；(b)两个局部极值点图
8.4 颗粒图像分水岭分割方法
在粘连颗粒图像处理过程中，通常颗粒物体会粘连在一起，妨碍单独颗粒特征参数的提取。这时需要一种有效的方法将互相粘连的颗粒分割开来。图像分割的结果好坏直接影响着后续图像分析的准确性，但是颗粒图像的复杂性和图像分割结果的不确定性，使得要得到一个良好的分割结果变得十分困难，因此颗粒图像分割是目前图像处理领域研究的热点和难点之一。分水岭分割方法借助地形学概念进行图像分割，分割结果为连通封闭的、单像素宽度的、位置准确的目标边界，对确定图像中目标所在的区域位置非常有利。在分水岭分割过程中，由于无法判断哪些点为合理的局部最小值点，因分水岭算法存在对噪声敏感，很容易产生过分割的问题，经常会将一个单独的目标错分割为几个目标。
本节首先对图像进行数学形态化处理，在简化图像目标和增强颗粒粘连的凹点特性的基础上，将二值图像进行距离变换后变为层次分明的灰度图像，再对图像进行分水岭分割，使得粘连颗粒分离开来。针对分水岭造成的过分割问题提出了改进的算法，即利用面积阈值法来判断分割的区域是否合理，如果不合理则进行区域合并来减少过分割的区域，有效地改善了分割的结果，为后续的图像分析、特征提取打下良好的基础。
8.4.1 分水岭变换原理及算法研究
分水岭算法是一种基于区域的图像分割方法，20世纪70年代Digabel和Lantuejoul第一次将分水岭的概念引入图像处理领域，用于二值图像的分析与处理。但是由于其计算量大、运算速度慢的原因，极大地限制了应用范围。直到1991年Vincent等人提出了基于“排序”与“浸没”的快速分水岭算法，算法能够满足实时性的需求，此后分水岭算法被广泛应用于灰度图像的分割。
对于一幅二维的灰度图像f(x，y)，像素点坐标与对应像素点的灰度值可以组成一个包含山峰和山谷的三维地形，在这地形中，对应于山谷位置的点一般称为局部极小值点(区域)，地形图中其他位置的点与极小区域的点有如下两种关系：
①若将一水滴置于地形图上的某一点，该水滴一定会流入特定的一个局部极小值点。②该水滴会等概率地流入一个以上的局部极小值点。
对于某一个局部极小值点，满足条件①的点的集合称为该极小值点的积水盆，满足条件②的点的集合称为分水岭，它将不同的积水盆相隔开来。
图8.30 分水岭算法示意图
如图8.30，分水岭算法的基本思想是模拟水流淹没地形的过程：①将灰度图像视为一个高低起伏的地形图，原图上较平坦的区域梯度值较小，构成盆地，原图上的边界区域梯度值较大，构成分割盆地的山脊；②水从盆地内最低的地方渗入，随着水位不断涨高，有的盆地将被连通，为了防止洼地之间被连通，就在分割两者的山脊上筑起水坝，水位越涨越高，水坝也越筑越高；③当水坝达到最高山脊高度时，算法结束，此时修筑的水坝就成了分水岭，每一个孤立的积水盆地对应一个分割区域。处理过程如图8.31所示。
8.4.2 分水岭算法的数学描述
根据分水岭算法就是通过水不断溢出淹没地貌的思路，对分水岭算法进行算法的数学描述。假设一幅需要分割的灰度图像f(x，y)，梯度图像为g(x，y)，图像中局部极小值点的位置用M1，M2，M3，…，MR表示。min和max表示梯度图中灰度的最小值和最大值。C(Mi)为极小值点Mi对应的积水盆点的集合，溢流过程是以单灰度值增加的，令n 为当前溢流的深度，T[n]表示g(u，v)＜n像素点的集合，即：
图8.31 分水岭算法处理过程图
(a)原始图像；(b)地形俯视图；(c)开始浸没；(d)积水盘的水逐渐汇聚；
(e)不断汇聚(已有分水岭)；(f)分水岭形成图
水流从灰度最小值处开始溢出，随着每次淹没过程中溢流深度从n 变为n+1，需要统计处在平面g(u，v)=n 下面的像素点集合T[n]。而对于每个局部极小值点Mi所在区域，处在平面g(u，v)=n 下面的像素点集合记为Cn(Mi)，计算公式为：
此公式说明可以用“与”操作将对应于Mi所在区域的那些小于n 的像素点提取出来。图8.32给出了上面给出的各个概念的直观解释，由图可以看出Cn(Mi)可以由C(Mi)和T[n]求并集得出，同理Cn(Mi+1)可以由C(Mi+1)和T[n]求并集得出。
若用C[n]表示溢流深度是n 时所有满足像素值小于n 的像素集合，C[n]计算方法如下：
则C[max+1]表示所有区域像素点的并集，计算方法如下：
图8.32 计算Cn(Mi+1) 示意图
在水平面不断增高的过程中，集合Cn(Mi)与T[n]中的像素会一直属于原集合，而随着n的增加，这两个集合中的像素个数是单调函数。容易知道C[n-1]是C[n]的子集，可得C[n]是T[n]的子集，三者之间的关系可以用如下公式表示：
由上面公式可知，每个C[n-1]中的连通组元都包含在T[n]的连通主元中。
分水岭算法开始进行初始化时，令C[min+1]=T[min+1]，然后利用迭代运算计算出溢流深度n 增加时的结果。假设在步骤n 时已经计算出C[n-1]的结果，要从该结果结算出C[n]，可以令S 代表T[n]中连通区域的集合，集合中每个连通区域s∈S[n]和C[n-1]存在三种关系：
关系1：s∩C[n-1]的结果为空集。
关系2：s∩C[n-1]的结果含有C[n-1]中一个连通区域。
关系3：s∩C[n-1]的结果含有C[n-1]中一个以上的连通区域。
对应三种关系的示意图如图8.33所示。
图8.33 计算C［n］ 的示意图
C[n]的计算过程和建立分水岭的时间都与上述关系中的某种关系成立有关。关系1表明随着n的增加，水平面遇到了新的局部极小值点，因此只要把新的连通区域s加入到C[n-1]就可以得到C[n]；关系2表明连通区域s只属于某个极小值区域，各局部极小值溢出的水并没有流进其他积水盆，此时只要把s加入到C[n-1]就可以得到C[n]；关系3表明水平面已经淹没了一个以上的积水盆，因此s中包含了C[n-1]中一个以上的连通区域，此时需要在s中建立分水岭，防止不同的区域会被连通起来。
8.4.3 分水岭算法及其代码实现
分水岭分割的主要优点是对图像的变化高度敏感，可以提取出图像中目标区域的封闭面，或者连通区域的边界。在实际应用中，如果直接使用分水岭算法对图像进行分割，会受到图像中噪声与目标本身的细节纹理的影响，导致图像出现过分割(over-segmentation)。从算法过程方面来说，是由于梯度噪声、目标本身的细节纹理和量化误差等的影响，均匀的目标内部会产生很多局部的“山峰”和“山谷”，分水岭分割的结果就会在目标内部出现很多“伪”边界，造成过度分割，使得正确的边界被大量不相关的边界淹没，如图8.34所示。
图8.34 分水岭分割结果图
(a)原始图像；(b)直接分水岭分割；(c)平滑处理后的分水岭分割
从图中的分割结果可知，直接使用分水岭分割图像出现的过分割现象非常严重，分割结果无法满足实际需要。针对这个问题学者们提出了很多改进方法，主要从三个方向对分水岭算法提出改进：
(1)使用滤波方法对图像进行平滑预处理，消除图像中的细节与噪声干扰来减少过分割，这种方法容易使图像变得模糊，而且效果不太明显。
(2)在处理过程中，通过控制标记种子点数量的方法来消除过度分割，这种方法能有效地防止过分割的现象，但对标记的提取没有固定的方法，一般需要事先具有足够关于图像目标的先验知识。
(3)在后期处理中，通过对分水岭分割的区域进行合并的方法来减少过分割现象，这种方法可以将过分割区域有效地合并在一起，但由于分割所得的区域数目过大，合并过程的计算量往往很大，寻找最优的合并准则也是难点。
8.4.4 分水岭分割的常用方法
1.基于距离变换的分水岭分割
在分水岭变换算法中经常配合使用的工具是距离变换，二值图像的距离变换反映了每个像素到最近背景点的距离。在实际的应用中，灰度图像中的目标物体由于存在纹理细节，使得灰度值发生突变，造成分割的结果中存在过分割现象。若将灰度图像先变为二值图像，再进行距离变换而产生新的灰度图像，新的灰度图像是一副灰度值均匀变化的图像，因此将消除掉物体本身存在的纹理细节引起的灰度突变，使得原灰度图像中存在的大量伪局部极小值点消失，再进行分水岭分割后，能取得比原先好很多的分割效果。
距离变换所得到的灰度图像中，灰度值与图像中目标的形状有关。若目标的形状不规则，即目标内部的像素点到背景的距离存在不止一个局部极小值点，则图像分割的结果会把一个完整的目标分割成几部分。所以若图像中感兴趣的目标形状是比较规则的，使用距离变换的分水岭才能取得较好的分割效果。
2.基于梯度的分水岭分割
在实际使用分水岭算法分割图像时，通常将梯度图像作为分水岭变换的对象。梯度图像中像素点的灰度值反映了原图像中该点与其领域点的灰度值变化程度。若直接对梯度图像进行分水岭变换，结果图像中会得到很多不感兴趣的分水岭，影响正确分割。一种改善的方法是先对梯度图像进行平滑操作，使得梯度图像中的伪局部极小值点减少，再对图像进行分水岭变换，从而改善分割的效果。
对原图像与梯度图像分别进行平滑操作比对原图像进行多次平滑操作的效果要好，而且该方法鲁棒性好，可以适用于不同应用领域图像的分水岭分割。虽然该方法能改善分割的效果，但是分割的最终结果依然存在过分割，对精度要求高的领域，该方法依然无法满足图像处理的需求。
3.基于控制标记符的分水岭分割
分水岭变换中由于存在噪声与细节纹理会导致过分割，以至于图像分割的结果无法满足实际应用需求。一种解决该问题的方法是将先验知识加入到分割过程当中，限制允许出现的区域数目。该方法具体实现基于标记符的概念，标记符是属于一幅图像的一个连接分量，可以对目标区域计算一个内部标记符集合，对大量的局部极小值逐个做标记。对背景部分也计算一个外部标记符集合。对于计算标记符的方法有很多，如线性滤波、数学形态学处理，当像素符合制定的标准时则将它归入标识符的集合。在这个过程中，符合规则的局部极小值都会归入到集合，此时可以观察标识符集合内的像素的特征，利用强制最小技术来修改图像的灰度值。强制最小技术是指保留小于某一阈值的局部极小值，大于阈值的其他像素点按需求将像素值增大到某个特定值，其结果可以通过调整阈值大小按实际需求减少一定数量的局部极小值。从而减少过分割的区域，改善图像分割效果。
该方法在制定标记符的选择范围时可以充分利用到先验知识，即从统计的目标形状、纹理、颜色特征等方面来制定合理的规则。通过对阈值的调整，可以很有效地控制分割区域的数量，因此该方法可以通过修正阈值来达到一个更好的分割效果，因此经常被用于各种图像的分割。
8.4.5 一种改进的分水岭分割方法
1.整个改进分水岭分割的算法流程
常用的分水岭分割方法直接对图像进行分割后，图像中大量纹理细节的影响，使得灰度图像中有许多伪局部极小值，分割结果中有大量的过分割现象，严重影响后续的物体特征参数提取。对于粘连的颗粒图像，粘连处的颜色介乎目标物体颜色与背景颜色之间，若在分割之前进行图像平滑处理，虽然减少了过分割现象，也会使互相粘连的不同物体因为粘连处颜色信息的丢失而混为一个物体；在分割过程中由于颗粒本身灰度不均匀不可避免地出现过分割现象，而各个过分割区域的特征都相似，虽然根据颜色特征相似规则将过分割区域进行合并从而减少过分割区域，但是合并也会将两个原本分开的物体合并为一个物体，因此无法有效分割粘连颗粒。
本节针对常用方法存在的问题，提出了一种结合数学形态学处理与距离变化的分水岭分割方法。该算法利用距离变换将物体的形状信息转化为均匀变化的颜色灰度信息，由于颗粒物具有类圆性，局部极小值只会出现在每个物体的中心处，因此在一定程度上限制了伪极小值的出现，从而控制分水岭分割出现的区域数。由于距离变换的分水岭只能解决一部分的过分割问题，因此需要在分水岭变换的前期与后期对图像进行优化，使得分割效果更加理想。前期的优化包含灰度图像平滑处理、二值图像形态学处理，这些处理主要消除噪声、杂质、反光不均匀和物体内部的纹理细节、突变点，以及边缘突出部分等因素对分割造成的干扰。而后期的优化则是根据分割开来的区域面积的大小，选取合适的阈值，将面积小于阈值的区域滤去，因为面积小于阈值的区域一般是过分割产生的物体边缘小部分以及两个物体互相接触处的细长的长方形和杂质物体。具体的整个改进分水岭分割的算法流程如图8.35所示。
图8.35 整个改进分水岭分割的算法流程图
在前面的8.2节和8.3节中已经介绍了图像的预处理、二值化和距离变化等处理，图8.36(a)为对优化后的二值图像，图8.36(b)为对图像进行分水岭变换得到的分水岭结果图。由于分水岭结果与原图像都是二值图像，可以将分水岭结果叠加到原二值图像中，得到结果图8.36(c)，可以看出分水岭将粘连的颗粒分割开来，为了更好地显示效果，将图像中不同的区域映射成不同颜色的RGB图像，结果如图8.36(d)所示，可以看出粘连的颗粒都具有不同的颜色，表明分水岭分割可以有效地分割粘连颗粒图像。
图8.36 玉米图像分水岭分割过程图
(a)玉米二值图像；(b)形成的分水岭；
(c)分水岭叠加到二值图像的分割图；(d)不同区域标注不同颜色图
2.分水岭分割算法目标颗粒中的应用
为了验证分水岭分割算法的适应性与鲁棒性，以玉米颗粒为主要处理对象，首先将该算法应用于相同背景与相同目标，对于不同的粘连状况的玉米颗粒图像，图8.37为比较典型的处理结果。从结果图中可以看出，尽管图像中玉米粒粘连严重，分水岭分割仍然可以有效地将两个以及两个以上互相粘连的颗粒分割开来。但是当图像中的玉米粒形状不是类圆形时，玉米颗粒的边缘不规则部分会被分为单独的一个颗粒，而且该部分具有面积小的特点。如图8.37中的黑色圈中的部分，有一单个的颗粒被分为两部分。这是由于几个物体互相粘连与物体本身形状不规则的影响，在距离变换过程中，处在中间的不规则目标物体内部的距离图中出现了两个伪局部伪极小值，导致分水岭变换出现分割错误的现象。从总体上来说，该分水岭分割算法对于具有相同特征(即具有相同的光照、拍摄位置、目标颗粒等)、不同的粘连状况的颗粒图像，都可以有效地进行分割，且分割效果比较理想。对于极少量过分割产生小区域，在后面会采用面积阈值法进行优化，使其不对特征分析与计数造成影响。
图8.37 不同分布的玉米图像分割结果
(a)玉米图像；(b)分割效果图；(c)玉米图像；(d)分割效果图
3.分水岭分割算法在不同领域的颗粒图像中的应用
在实际的自动化处理中会碰到粘连颗粒图像的问题，图8.38为将该算法应用在不同领域分割结果的图片。图8.38(a)为农作物中米粒质量分析的图像分割；图8.38(c)、图8.38(e)为生物医学领域中细胞特征分析的图像分割；图8.38(g)为工业上棒材自动计数过程中需要的捆绑棒材的图像分割。从各分割结果图中可以看出，只要目标物体具有类圆性，而且与背景的颜色具有较大的灰度差异，改进的分水岭算法都能很好地将粘连物体分割为单独的物体。在分割过程中，由于不同类型的图像其颗粒大小不同，因此在数学形态学处理中需要选择不同大小的结构因素对图像进行闭运算与腐蚀运算，使得形态学滤波器能最佳地滤去图像中的杂质、不圆滑的边缘，减少两个粘连物体的接触面积。
另外，对于不同类型的图像，分割效果也有一些差异。如在大米图像中，虽然米粒的形状接近于长方形，但是没有出现互相粘连以致重叠的现象，因此分割效果也比较好。在圆形细胞图像中，虽然可以很好地分割开来，但对于过度粘连的细胞，出现了重叠的状况，在分割后在接触部分会失去一部分的形状特征，这对计数与颜色特征的提取无影响，但对统计形状特征有一定影响。在钢材图像中，钢材截面具有较好的类圆性，接触部分与背景颜色反差较明显，因此分割效果比较理想。
4.分割后期的面积阈值优化方法
对图像进行分水岭分割后，图像中有少量不具类圆性的粘连目标颗粒，其不规则的边缘部分被分割成单独的小面积物体，如图8.39(a)中小圆圈中所标出的部分。为了不影响颗粒计数，可以利用面积阈值法将这部分区域除去。通常用一个区域中的像素点的总数来作为该区域的面积S，通过对同一类型的大量图片进行面积分布统计，可以找出一个合适的面积阈值T，将面积小于T 的区域中的所有像素点的值设为0(除去该区域)，用公式表示为：
图8.38 不同领域的图像分割效果图
(a)大米图像；(b)分割效果图；(c)血细胞图像；(d)分割效果图；
(e)血细胞图像；(f)分割效果图；(g)棒材图像；(h)分割效果图
图8.39为玉米颗粒图像中各个区域的面积统计直方图，从图8.39(c)中可以看出，正常的玉米粒面积都大于400，而错误分割的区域面积小于200，可以设定阈值T=300。由于同一类型的图片中颗粒大小相似，因此该阈值可以处理所有的玉米颗粒图像。利用该阈值对图像进行处理，为了效果更加明显，只显示图像中错误分割的小区域，如图8.39(b)所示。从图中可以看出，面积阈值法可以很好地找出错误的分割部分。在实际应用中，分水岭分割后期利用该方法将这些小区域从图中滤去，可以达到更好的分割效果。
图8.39 面积阈值法示意图
(a)分割后的图像；(b)过分割产生的小区域图；(c)面积分布直方图
8.5 图像颗粒的自动计数与特征参数提取
8.5.1 颗粒自动计数
实现图像中颗粒的自动计数是本设计的主要目的之一。通常颗粒的粘连程度与杂质物体的存在都会影响到计数的准确性，本节在研究过程中设计了两种方法来进行计数。一种是利用面积法来判断颗粒数，另一种是通过统计改进分水岭分割后的区域个数来实现计数功能。这两种方法都能比较好地实现计数功能，但面积法会破坏颗粒本身的形状信息，而区域计数的方法则保持了颗粒的基本信息，在此基础上可以很好地提取颗粒的特征参数。
(1)面积法计数
面积法计数的思想是通过对原图像的二值图像进行数学形态学处理，选择合适的结构因素，对颗粒物体不断地进行腐蚀，使得互相粘连的颗粒接触部分越来越少，直至分割开来，而颗粒物体不会被完全腐蚀掉。最终腐蚀结果是图像中大部分的颗粒被独立地分割出来，只有少数的颗粒还粘连在一起。然后对物体面积进行统计，设定多个阈值T1、T2、T3、T4。根据面积在某个阈值范围内，判断该物体是单独的颗粒，还是由几个粘连颗粒组成的物体。具体的判断公式如下所示：
公式中nums为总颗粒数，当物体面积小于T1时，将该物体判为单独的颗粒，总数加1；当物体面积大于等于T1小于T2时，则判为2个颗粒粘连而成的物体，其他阈值范围也相类似。对本节设计的采用玉米颗粒图像进行数学形态学处理后的结果如图8.40所示。通过大量的实验表明，对图像进行不断腐蚀后粘连在一起的颗粒数不会超过5个，因此在公式(8.40)中，将面积大于T4的物体都视为5个颗粒粘连而成。图8.40(c)为面积分布直方图，根据直方图上面积的分布规律，可以分别设定阈值为800、1600、2200、2800。利用这些阈值对不同分布、不同数目的玉米颗粒图像进行计数，其结果误差都不超过2%。该方法过程简单，由于腐蚀运算可以将大部分的粘连颗粒分割开来，因此省略了分割过程，提高了效率，也使计数结果不受物体粘连程度的影响。
图8.40 面积法计数示意图
(a)玉米图像；(b)多次形态学处理图；(c)面积分布直方图
(2)基于分水岭分割区域的计数
利用改进分水岭算法对图像进行分割后，图像中的颗粒物体都被分割成不同的区域。对不同区域进行顺序标号，使得不同区域中像素点的值等于其区域号。对玉米颗粒图像中各区域进行标号如图8.41所示。其中最大的区域号即为图像中颗粒的总数。该方法是在改进分水岭算法分割的基础上，对颗粒进行计数。在计数过程中，不会破坏颗粒的形状信息，而且计数具有较高的准确率，也可以使用于其他应用领域不同类型的图片。同时计数过程中对各区域进行标号，方便了后面特征参数提取，具体的区域标号方法在面积特征参数提取中作介绍。
图8.41 对不同区域标号图
利用上述两种方法对不同数目、不同分布的玉米颗粒图像进行计数，其计数结果统计如表8.1所示。从表8.1中可以看出，在相同拍摄条件下，随着颗粒数目的增加，颗粒粘连情况越严重，使得两种计数方法的准确率都有所降低。
表8.1 不同数量与分布状况的玉米图像计数结果
8.5.2 颗粒参数提取
在实际应用中一般可以通过颗粒物体的大小、形状和颜色特征等来对颗粒进行描述。颗粒的大小是最容易区分不同颗粒的一个特征，在数字图像处理中可以计算出颗粒的粒径、面积、周长、长轴和短轴等参数来表示大小特征。而颗粒具有各种各样的形状，如定性的分析可分为球状、多角状、针状、枝状等不同形的分水岭，具有封闭性的特点，因此可以通过求出完整区域的颜色特征来区分不同的颗粒状，定量的分析则可根据颗粒的形状因子、偏心率、圆形度等来表示形状特征。
1.颗粒大小特征参数提取
(1)周长特征参数提取
图像颗粒物体的周长是一个简单的边界全局特征，它是指包围目标区域的轮廓的边界长度。目标区域的周长特征值在识别复杂形状目标时特别有用。对于图像中的目标区域R，其边界像素p 满足以下两个条件：①像素p属于目标区域R 的集合。②像素p 的邻域像素中含有不属于R 的像素。若像素只满足第一个条件，则此像素为目标区域内部的像素，只满足第二个条件的像素点为背景像素点。在测量周长的过程中经常包含有许多90°的转弯，使得周长值比实际值要大。周长的表示方法有多种，其测量出的结果也有差异。周长的测量方法有：
1)利用隙码求目标区域的周长。可以将图像中的每个像素都看成一个小正方形的区域，如图8.42所示。目标像素与背景像素的粘连处为小正方形的一条边界，可用隙码表示。目标区域的周长就是计算区域所有隙码的长度和。
2)利用链码求目标区域的周长。链码是图像目标区域边界点的一种编码表示方式，它利用一组规定长度和方向的直线来表示目标的边界。边界的起点需要用绝对坐标表示，其他点可以用后续的方向链码所代表的偏移量来表示。链码表示利用一个方向值来代替一点的二维坐标，因此可以减少数据量。图8.43为4方向链码和8方向链码的示意图。其中4方向的链码长度为1，8方向链码中水平和竖直方向的链码长度为1，对角线方向的链码长度为
。
图8.42 边界编码图形
图8.43 方向与链码的对应关系
(a)4方向链码；(b)8方向链码
对于图8.43所示的目标区域边界，若将点O(0，0)当做起始点，则其4方向链码表示：(0，0)1 1 1 2 3 2 3 2 3 0 0；8方向链码表示为：(0，0)2 2 2 4 5 5 6 0 0。
目标边界的周长即为水平和竖直方向上的链码个数加上
乘以对角线上链码的个数。由于8方向链码具有更高的准确性，通常使用它来计算周长。对于某个目标区域，在计算过程中，所有水平与竖直方向的链码为偶数，其总数记为Ne。所有对角线方向的链码为奇数，其总数记为No，目标周长P 可用如下公式计算：
3)将目标区域的边界像素的总数作为目标的周长，该方法计算简单，具有统计意义。
对图8.43(a)所示的由0、1、2、3连接所构成的目标区域，利用三种不同的方法所得的结果为：缝码计算结果为16，利用链码表示时周长计算结果为
，利用统计像素个数作周长时结果最小为10。
上述三种方法计算出的周长结果都是具有统计学上的特征，与实际的周长有较大的差异，其中方法3)最简单，可以利用该方法来提取玉米颗粒图像中所有玉米粒的周长特征。如对图8.40(a)进行提取，其周长分布如图8.44所示。从结果可以看出，周长在140左右的颗粒数目最多，周长整体分布趋势为中间颗粒数多，越向两边颗粒数越少，呈正态分布。
图8.44 周长分布直方图
(2)面积特征参数提取
面积是目标区域大小特征的一个度量，它只与区域的边界包围的范围有关，与内部像素值无关。实际应用中通常用目标区域中的像素个数总和来表示面积。对于一幅M×N 的二值图像f(x，y)，若二值图像灰度中值含有一个感兴趣的目标，则面积A 为像素值为1的点的个数和，即：
对于含有k个目标区域的图像，首先需要将各个分离的目标区域进行标号，以方便各个目标区域的面积计算。区域标号分给每个目标区域唯一的序号，同一区域内的所有像素的灰度值都记为该区域的区域号，其中目标区域的个数与最大的区域序号相同，背景像素的灰度值依然为0。线性扫描的区域标号方法如下：
首先对一幅二值图像进行从上到下、从左到右的扫描，直到扫描到灰度值为1的像素点。此时由于扫描顺序的原因，该像素点p 的8邻域中左、左上、上、右上四个点一定为已经处理过的像素点，当前像素点的标号与四个像素点的值相关，分三种状况：
1)若四个邻域像素点的灰度值都为0，当前像素点赋予新的标号。
2)若四个邻域像素中有像素点已经被标注，且所有像素的标号一致，则当前像素点标上该标号。
3)若四个领域像素点中的已经标注的点标号不同，则按照左、左上、上、右上的顺序寻找不为0的像素标号，将该标号赋予当前像素点。
用上述的方法对所有的像素进行一次扫描后，所有的像素都被标号，如图8.45(b)所示。由图8.45(b)可知第一次扫描完成后，存在着同一目标区域像素点有几个标号的问题。解决该问题的方法是在第一遍扫描过程中如果发生一个区域多标号登记的情况，当扫描完成后对登记表中的标号进行整理。在第二次扫描过程中对同一区域的不同标号子区域进行匹配，使得一个目标区域只有一个标号。
图8.45 线性扫描法区域标号过程
(a) 二值图像矩阵示意图；(b) 一遍扫描结束时的标号结果；(c) 第二遍扫描后的标号图
完成对不同目标区域的标号之后，不同区域像素的灰度值变为该区域的标号，统计不同灰度级的像素个数，就是各个目标区域的面积。
对玉米颗粒图像进行面积特征参数提取，其结果如图8.46所示。从结果可以看出颗粒的面积大小分布在950左右。
图8.46 面积分布直方图
(3)粒径参数提取
粒径是对物体大小特征的一个度量，在圆形物体中，粒径等同于直径，对于其他不规则的物体，无法直接计算其粒径，因此通常使用等效的方法来测量物体的粒径。在实际应用中，由于测量的面积与周长比较真实地反映出目标物体的形状特征，可以与物体有相同面积的圆的直径来定义物体的等效粒径：
由于在实际检测中面积特征参数比较真实地反映物体的特征，周长特征参数由于方法不同而结果不一样，因此最后得出的等效粒径存在误差。
图8.47 粒径分布直方图
对玉米颗粒图像进行粒径特征参数提取，其结果如图8.47所示。从结果可以看出颗粒的面积大小分布在950左右。
(4)长轴和宽轴特征参数提取
物体的长轴与宽轴是反映物体外形的重要参数之一。在计算过程中，由于物体的不规则性，无法定义其长宽，可以通过外接矩形的长度与宽度来间接地表示物体的长宽。外接矩形是完全包围住物体的矩形。只要计算出最左与最右边界点的x 轴坐标之差就可以算出外界矩形的长度，最上与最下边界点的y 轴坐标之差就可以算出外界矩形的宽度。对于图像中任意放置的物体，由于其摆放方向的原因，其外接矩形可能并不是最佳的，即不是最小外接矩形。可以将物体沿主轴在90°角范围里旋转，每旋转一次记录一次外接矩形的长度和宽度，当矩形面积最小时则判定为最佳外接矩形，计算出的长度和宽度最接近物体本身长宽值。
对玉米颗粒图像进行长轴特征参数提取，其结果如图8.48所示。从结果可以看出颗粒的长轴大小分布在39左右。颗粒的宽轴特征参数提取方法也相类似。
图8.48 长轴分布直方图
2.颗粒形状特征参数提取
(1)偏心率特征参数提取
偏心率是物体外接椭圆两焦点间的距离与长轴长度的比值，它是反映物体紧凑性(目标区域关于中心点的对称程度)的一个重要参数。由于直接计算偏心率比较复杂，在实际应用中通常是通过简化的方法来计算偏心率，如直接计算目标区域的长轴p 与短轴q的比值作为近似的偏心率：
图8.49 偏心率分布直方图
计算出的偏心率越大，物体外接椭圆越扁，偏心率越小则物体的形状越接近于圆。该方法求出的偏心率不受目标物体的旋转、平移等变换的影响，但是由于是近似计算的原因，结果不是非常精确。
对玉米颗粒图像进行偏心率特征参数提取，其结果如图8.49所示。从结果可以看出颗粒的偏心率主要分布在0.5左右。
(2)形状因子特征参数提取
形状因子是反映物体与圆形接近程度的一个特征参数，圆形的形状因子最小为1，其他形状的物体形状因子都大于1。当目标物体与圆形相差越大时形状因子也越大。利用目标区域的周长P 和面积A 可以定义形状因子F，公式如下：
形状因子对物体的旋转、平移等操作不敏感，但是由于在数字图像处理过程中，提取周长与面积的特征参数只是反映物体特征的一个数字化的度量，与实际值成一定的比例，但是也存在着差异。如图8.50给出了不同形状的物体却拥有相同的周长与面积，其计算的形状因子不能准确地反映物体的形状。因此为了更好地识别目标物体，对物体进行分类，最好是结合物体的几种特征参数一起进行判断。
图8.50 形状因子相同的不同形状(A=5，P=12，F1=F2=F3)
对玉米颗粒图像进行偏心率特征参数提取，其结果如图8.51所示。从结果可以看出颗粒的形状因子主要分布在1.4左右。
图8.51 形状因子分布直方图
(3)颜色特征参数的提取
颗粒物体的颜色特征是区分不同颗粒的重要信息，它不受物体的形状与大小的影响。颜色信息可以分为彩色信息和灰度信息。其中彩色信息可以用RGB颜色空间三个分量来表示，而灰度信息则由像素点的灰色程度来表示，两者的取值范围都为[0，255]。在图像中由于颗粒由不同像素组成，求某个颗粒的颜色信息需要统计该颗粒区域中每个像素点的颜色信息。如利用区域的平均灰度值来表示该颗粒的灰度信息，公式为：
式中Sn为标号图像中区域n中像素点的集合，h(i，j)为原始图像中像素点的颜色信息，nums(n)为区域n 的像素点的总个数。利用上述公式计算图中的玉米颗粒的R分量上的颜色信息，结果如图8.52所示。从结果中可以看出玉米颗粒在R分量上的灰度值主要集中在175附近，由于图像中只有玉米颗粒，各个物体的颜色信息相似，分布在[140，200]范围内。
随着计算机技术不断提高，基于颗粒图像处理的自动化处理技术在不同领域的应用越来越多。为了解决粘连颗粒图像分割与分析过程中的一些问题，本节以粘连玉米为实验对象，对整个颗粒图像处理过程做了深入的研究，总结出如下结论：
(1)提出了一种新的图像预处理的方法。将彩色图像转化为灰度图，除去图像中不需要的彩色颜色信息。利用中值滤波算法对图像进行平滑，该算法除去了图像中各类噪声的干扰，提高了图像的质量。再利用最大类间方差法选取合适的阈值，将灰度图像转化为二值图像，使得颗粒对象从背景中分割出来，除去了颗粒的细节纹理信息。同时设计了一种数学形态学滤波器，滤去了图像中颗粒不规则的边缘，使得互相接触部分减少，填充颗粒中的孔洞。该滤波器具有处理速度快的特点，在很大程度上优化了图像。
图8.52 颗粒灰度值分布直方图
(2)在分割互相粘连颗粒的过程中，针对传统分水岭分割算法中经常出现过分割现象，提出了一种改进的分水岭算法。该算法通过距离变换将二值图像转换为新的灰度图，消除了原始灰度图像中大量存在的“伪”局部极小值点，对灰度图进行分水岭分割算法，使得互相粘连的颗粒很好地分割开来，在后期对少量过分割产生的小面积区域，采用面积阈值法滤去面积小于阈值的区域，使得分割的效果更好。改进分水岭算法可以适用于不同领域的粘连颗粒图像，具有良好的适应性与鲁棒性，处理速度较快。
(3)提出了两种颗粒自动计数的方法，一种是利用数学形态学腐蚀运算对图像进行处理，再根据面积分布直方图判断不同面积区域含有的颗粒数；另一种是对改进分水岭算法分割后的单独颗粒进行标号，最大标号数作为颗粒总数。另外，提出了利用颗粒的大小、形状、颜色等特征参数来描绘颗粒的特征，提取出各项特征参数的值。该方法得到的特征参数值误差较小、效率高，为颗粒的识别与归类打下了良好的基础。
参考文献
[1]Kotsialos A，Papageorgiou M，Diakaki C，Pavlis Y，Middelham F.Traffic flow modeling of large-scale motorway using the macroscopic modeling tool METANET[J].IEEE Transactions on Intelligent Transportation Systems，2002，3(4)：282-292.
[2]Yu Huang，Huang TS.A region-based method for model-free object tracking[C].Proceedings 16th International Conference on Pattern Recognition，2002，1：592-595.
[3]Schmaltz Christian.Region-based pose tracking[C].3rd Iberian Conference on Pattern Recognition and Image Analysis，2007，4478：56-53.
[4]Chiou-Ting Hsu，Ming-Shen Hsieh.Region tracking for non-rigid video objects in a nonparametric MAP framework[J].Signal Processing：Image Communication，2006，21(3)：235-251.
[5]柴毅，阳小燕.基于运动区域检测的运动目标跟踪算法[J].计算机应用研究，2010，27(12)：4787-4789.
[6]Comaniciu D，Meer P.Mean Shift Analysis and Applications[C].Proceedings of the Seventh IEEE International Conference on Computer Vision，1999：1197-1203.
[7]彭宁嵩，杨杰.Mean Shift跟踪算法中核函数窗宽的自动选取[J].软件学报，2005，16(9)：1542-4545.
[8]王杰，王加银.Mean Shift算法的收敛性讨论[J].北京师范大学学报，2008，44(5)：472-475.
[9]Li Z，Tang.Q.L，Sang N.Improved mean shift algorithm for occlusion pedestrian tracking[J].Electronics Letters，2008，44(10)：622-623.
[10]聂国梁，卢正鼎.流数据密度估计[J].计算机科学，2006，33(12)：94-96.
