CN102324092B - 一种数字图像中颗粒状对象的自动分割方法 - Google Patents

Abstract

本发明公开了一种数字图像中颗粒状对象的自动分割方法，属于数字图像处理技术领域。本发明针对数字图像，尤其是显微图像中颗粒状对象的灰度、结构分布与几何形状等特点，先应用自动阈值法将目标与背景分离；然后计算其梯度向量场，在梯度向量场中搜索关键点，理想的关键点在8邻域均有相应的梯度向量分布，其梯度值为零，获取的关键点作为每一个颗粒状对象的中心；接着定义一个新的基于灰度与空间位置的有效能量函数用以计算方向梯度，将其代替传统的灰度梯度；最后应用主动轮廓模型搜寻颗粒状对象的边界。本发明能准确有效地分割聚集型的颗粒状对象，尤其是生物医学显微图像中存在的大量粘连或重叠的微小颗粒，为图像分析与识别提供帮助。

Description

一种数字图像中颗粒状对象的自动分割方法

技术领域

本发明属于数字图像分析与处理技术领域，特别是数字图像中颗粒状对象的自动分割。

背景技术

数字图像，尤其是生物医学显微图像中，许多对象的结构或微结构呈颗粒状，如淀粉类食物土豆、小麦的淀粉颗粒、嗜热菌、球形球菌、红细胞、小淋巴细胞等的显微图像，其特点为：(1)颗粒数量较多，呈聚集分布状态；(2)颗粒间距离很小，存在大量粘连及重叠现象；(3)形状规则，大多为圆形或近似圆形。在进行生物学、医学基础研究时，常常需要对这些颗粒的结构与分布特性进行分析和处理，如细胞/颗粒计数、病毒特征提取与识别、异常判断、生长过程监测与动态评估等。细胞/颗粒状对象的分割，是进行上述分析的基础与重要环节，然而，由于有限的光照条件，以及大量颗粒状对象发生粘连或重叠等原因，使得它们的自动分割变得十分困难。

图像分割一直是图像处理领域引起人们重点关注的热点与难点问题之一。众多图像分割方法中，比较有代表性的是基于灰度与轮廓几何形状的分割方法，主动轮廓模型(active contourmodel，ACM or snake)，最早由Kass等人于1988年提出，(Kass M，Witkin A，Terzopoulos D.Snakes：Active contour models，Int.J.Comput.Vision.321-331，(1988))。它建立一个包含基于图像灰度外能与轮廓内能的函数，通过轮廓曲线内力与灰度梯度外力的驱动，使初始轮廓产生形变，最终收敛到真实的目标边界，从而实现目标的分割。随后，出现了许多改进的方法，如Xu等人(Xu C，Prince J L. Snakes，shapes，and gradient vector flow.IEEE Trans ImageProcessing 7(3)：359-369(1998)；Xu C.，Prince J L. Generalized gradient vector flow externalforces for active contours，Signal Processing 71(2)：131-139(1998))提出的梯度向量流(gradientvector flow，GVF)，扩大了传统Snake的捕捉范围，对噪声具更强的健壮性。

近年来，另一类应用很广泛的方法，水平集(Level set)，也被用于细胞或颗粒状对象的显微图像的分割(Solorzano De C，Malladi R，Lelièvre S，et al.Segmentation of nuclei and cells usingmembrane related protein markers，Journal ofMicroscopy 201：404-415(2001)；Vese L. Chan T.Amultiphase level set framework for image segmentation using the Mumford and Shah model，International Journal of Computer Vision 50(3)：271-293，(2002)；Yan P，Zhou X，Shah M，et al.Automatic segmentation of high throughput RNAi fluorescent cellular images，IEEE Transaction onInformation Technology in Biomedicine 12(1)：109-117(2008))。基于几何形变模型的分割，无论是传统的Snake，还是改进的GVF，及气球(Ballon)力模型，或B样条Snake等方法，它们分割单一的目标时，经过合适的初始化后，能快速准确地搜寻到目标边界；水平集算法，虽然能将二维轮廓的搜索问题变换到三维空间中表面拓扑结构的演化，使其更适合分割多目标，但它存在构建速率函数困难、控制参数较多、时间开销大等不足。更为重要的是，这两类方法，对于图像中存在大量粘连与重叠的对象，其性能并不理想，它们无法检测到正确的对象边界，实施准确有效的分割。

发明内容

本发明目的在于针对现有技术存在的缺点与不足，提供一种数字图像中颗粒状对象的自动分割方法。

本发明的目的通过如下技术方案实现：

一种数字图像中颗粒状对象的自动分割方法，包括如下步骤：

第一步分离颗粒状对象与背景：先将图像转化为二值图像，再采用自动阈值法将待分割的对象与背景进行分离；

第二步搜索颗粒状对象的中心：计算图像的梯度向量场；求颗粒中心点；

第三步计算方向梯度：先构建一个基于灰度与空间位置的能量函数，计算每个像素对应的能量最大的中心点，作为此像素同质性最强的中心；计算每一个像素与其对应的中心的方向梯度；

第四步初始化圆形轮廓：从每一颗粒的中心出发，分别沿上、下、左、右四个方向，搜索第一个非零值的像素，将中心与该像素的距离，作为初始化圆的半径；

第五步分割颗粒状对象：利用GVF Snake形变模型实施分割，其外力项的梯度采用第三步中的方向梯度。

第三步中定义一个基于灰度与空间位置的有效能量函数：

E_effect＝k₁E_pos+k₂E_homo，

其中 $E_{\mathrm{pos}}=e^{{-\mathrm{\γ}}_1}\left|\right|P-C\left|\right|,$ $E_{\mathrm{homo}}=e^{{-\mathrm{\γ}}_2|f(i,j)-{\mathrm{\μ}}_k|};$

式中，k₁与k₂为权值参数；

γ₁与γ₁为大于0的常数；

P为像素；C为颗粒中心；

μ_k为颗粒C_k的灰度均值。

第三步中的方向梯度定义为：图像中某一像素点P(p_x，p_y)指向某一中心C(c_x，c_y)的方向梯度向量为d(x，y)＝(d_x，d_y)，其中 $d_x=\frac{c_x-p_x}{\left|\right|P-C\left|\right|},$ $d_y=\frac{c_y-p_y}{\left|\right|P-C\left|\right|}.$

与现有技术相比，本发明具有如下有益效果：

(1)自动定位颗粒中心点，确定颗粒数量。

(2)通过建立基于灰度与空间位置的同质函数，每一像素均能唯一对应与其同质性最高的中心点，即每一中心点均有一个同质局部区域，或称力场范围。

(3)以方向梯度代替传统的灰度梯度，一方面，由于只需要计算一个方向的梯度，因此可以减少计算量，另一方面，可提高分割结果的准确性。

(4)分割过程不需要任何人工参与，分割过程快速，结果准确。

(5)对聚集性较强颗粒的分割具有显著特异性，能有效分割粘连或重叠的颗粒。

附图说明

图1为本发明一种数字图像中颗粒状对象的自动分割方法的流程图。

图2为本发明8邻域方向及向量之间相似性。

图3为本发明实施例土豆淀粉颗粒显微图像。

图4为本发明实施例图像的直方图。

图5为本发明实施例图像的颗粒与背景分离后二值图。

图6为本发明实施例图像的GVF场与颗粒中心。

图7为本发明实施例指向各颗粒中心的向量图。

图8为本发明实施例求出的方向梯度图。

图9为本发明实施例的初始圆形轮廓图。

图10为本发明实施例土豆淀粉颗粒显微图像的分割结果。

图11～13为本发明实施例的三个研究对象。

图14～16为本发明实施例的三个研究对象对应的分割结果。

具体实施方式

下面结合实施例和附图对本发明的技术方案作进一步的描述，但本发明的实施方式不限于此。

实例采用土豆淀粉颗粒的显微图像，由DP71 Olympus光学显微镜采集。

如图1所示，本发明包括以下步骤：

(1)输入待分析的土豆淀粉颗粒的显微图像；

(2)分离颗粒状对象与背景；

从图4淀粉颗粒显微图像及其直方图可知，图像主要有两个相隔较大的主峰，较低的峰值区间是灰度值较小的淀粉颗粒，较高的峰值区间为高亮度背景。

采用全局自适应阈值法将对象与背景分离：首先给定一个初始阈值(如图像灰度平均值)，然后通过迭代算法，逐渐使阈值逼近最优值，最后求出准确的阈值，对图像实施分割。具体过程如下：

应用该阈值将图像分为两类，分别计算两类对象的平均密度，求两个平均密度的均值并作为新的阈值，再对图像进行分类，比较连续两次的阈值之差，或判断迭代次数是否达到最大值，决定是否结束处理过程。其过程可描述为：

a)选择初始估计值T₀，给定很小的终止值t，及最大迭代次数N_max；

b)用T₀作为阈值将图像分为C₁与C₂两类；

c)对C₁与C₂中所有像素计算平均密度

${\mathrm{\μ}}_1=\frac{1}{N_{C_1}}\underset{i\∈c_1}{\mathrm{\Σ}}f\left(i\right),$ ${\mathrm{\μ}}_2=\frac{1}{N_{C_2}}\underset{j\∈c_2}{\mathrm{\Σ}}f\left(j\right)---\left(1\right)$

其中f(·)为图像密度，与

分别为C₁与C₂的像素数量。

d)计算新阈值

$T_1=\frac{{\mathrm{\μ}}_1+{\mathrm{\μ}}_2}{2}---\left(2\right)$

e)计算连续二次阈值之差ΔT＝|T₁-T₀|，如果ΔT＜t，或迭代次数等于N_max，则结束；否则，将新阈值T₁赋给T₀，重复步骤b)-e)。

图5为经自动阈值法处理后，颗粒与背景分离后的结果。

(3)搜索颗粒状对象的中心

当得到颗粒与背景被分离后的二值图像，首先计算梯度向量场，再求颗粒中心。

(3.1)计算梯度向量场(GVF)

基于梯度向量场v(x，y)＝(u(x，y)，v(x，y))的能量函数定义为：

$E_{\mathrm{gvf}}\left(v\right)=\∫\∫\mathrm{\μ}(u_x^2+u_y^2+v_x^2+v_y^2)+{|\▿f|}^2{|v-\▿f|}^2\mathrm{dxdy}---\left(3\right)$

其中，f(x，y)为图像边缘图，u_x与u_y分别为u(x，y)和v(x，y)在x与y方向的偏导，||表示向量的模。

求(3)式能量函数的最小值，可应用变分法求解以下两个Euler方程：

${\mathrm{\μ}\▿u}^2-(u-f_x)(f_x^2+f_y^2)=0---\left(4\right)$

${\mathrm{\μ}\▿v}^2-(v-f_y)(f_x^2+f_y^2)=0---\left(5\right)$

其中

为Laplacian算子，f_x与f_y分别为f在x与y方向的偏导。

GVF场通过以下迭代法求得：

$u_t(x,y,t)={\mathrm{\μ}\▿}^2{u}_t(x,y,t)-(u(x,y,t)-f_x(x,y))(f_x{(x,y)}^2+f_y{(x,y)}^2)---\left(6\right)$

$v_t(x,y,t)={\mathrm{\μ}\▿}^2{v}_t(x,y,t)-(v(x,y,t)-f_y(x,y))(f_x{(x,y)}^2+f_y{(x,y)}^2)---\left(7\right)$

其中t为迭代次数。

(3.2)求颗粒中心

由于颗粒大多为圆形或近似圆形，理想的中心点具有两个特征：(a)GVF向量指向8邻域方向，从中心点出发，向外的8邻域方向均有向量分布；(b)梯度值为0。但实际图像中，由于光照不均匀或噪声的影响，中心点的GVF向量并不是完全指向8邻域方向，而是与8邻域方向相近，因此，如果两个向量的相似性足够大，则认为二者是一致的。向量的相似性测度定义如下：

其中S＝(S_x，S_y)表示考察的中心点，

是像素S的8邻域像素， $V_i=(V_x^i,V_y^i)(i=\mathrm{1,2}...,8)$ 为N_i的灰度梯，·表示向量点积，

为从S指向N_i的方向向量：

上述向量的相似性测度，实际上是两个向量夹角的余弦值，取值范围为[0，1]，其值越大，则向量越相似。本发明实例中，如果两个向量的相似值大于或等于0.8，则认为它们是一致的，在8邻域方向中，若存在6个以上方向的GVF向量，则该像素为中心点。

8邻域方向及向量之间的相似性如图2所示，图6为求得的颗粒中心点。

(4)计算方向梯度

(4.1)计算有效能量

大多数基于形变的轮廓搜索方法，很难分割粘连或重叠的颗粒，主要原因在于大量颗粒的聚集，力场较为复杂，其中灰度梯度较大的颗粒，对轮廓具有较强的吸引力，往往会驱动周围灰度较小颗粒的轮廓向其运动，从而收敛到错误的边界。

为了解决此问题，本发明定义以下能量函数：

E_effect＝k₁E_pos+k₂E_homo    (10)

其中k₁与k₂为权系数。首项E_pos，定义如下：

$E_{\mathrm{pos}}=e^{-{\mathrm{\γ}}_1\left|\right|P-C\left|\right|}---\left(11\right)$

其中P为像素，C为颗粒中心，γ₁为大于0的常数。

第二项则为：

$E_{\mathrm{homo}}=e^{-{\mathrm{\γ}}_2|f(i,j)-{\mathrm{\μ}}_k|}---\left(12\right)$

γ₂为大于0的常数，μ_k为颗粒C_k的灰度均值。

(10)式的能量函数，本发明中称为有效能量，E_pos与空间位置有关，反映了像素与颗粒中心的距离关系，其值是像素与中心之间距离的负指数函数，即二者之间距离越小，其值越大，随着距离的增大，其值呈负指数递减。第二项E_homo则反映了像素与颗粒的灰度差异，它是像素与颗粒灰度均值之差的负指数函数，即二者的灰度相差越小，其值越大，反之亦然。该能量函数本质上反映了当前像素与颗粒的同质性，如果当前像素的灰度与颗粒灰度均值相近，距离又较小，则其有效能量值较大，反之则较小。每一像素相对每一颗粒及其中心，都有一个有效能量值，取能量最大者，作为当前像素同质性最强的中心，以便计算方向梯度。

(4.2)计算方向梯度

由于步骤(4.1)中，每一像素均有唯一的同质性最强的颗粒相对应，因此，可从该像素P(p_x，p_y)出发，得到指向对应颗粒中心C(c_x，c_y)的向量D(d_x，d_y)为：

$d_x=\frac{c_x-p_x}{\left|\right|P-C\left|\right|}---\left(13\right)$

$d_y=\frac{c_y-p_y}{\left|\right|P-C\left|\right|}---\left(14\right)$

本发明以此方向的灰度梯度代替传统的灰度梯度，可以减弱灰度梯度值较大的颗粒对附近梯度值较小颗粒的影响，使轮廓收敛到准确的颗粒边界。实际图像中，向量D(d_x，d_y)也存在与8邻域方向不完全一致的情况，此时，与求颗粒中心时的方法类似，计算两向量的夹角余弦相似性，如果二者的相似性值大于或等于0.8，则认为二者是一致的。

图7，图8分别是指向中心点的方向矢量，及对应的方向GVF场。图7表明，每一像素均指向与其同质性最强的中心点，即每一中心点有一个力场范围。

(5)初始化圆形轮廓

由于颗粒为圆形或近似圆形，因此，轮廓初始化为一定半径的圆较为合适。对于背景与颗粒分离后的二值图像，可从每一颗粒的中心出发，分别沿上、下、左、右四个方向，搜索第一个非0值的像素，将中心与该像素的距离，作为初始化圆的半径，由此，每一颗粒的初始轮廓均为一内切圆。图9为初始化圆形轮廓。

(6)分割颗粒状对象

每一颗粒经过初始化后，GVF Snake形变模型的外力项中的梯度采用方向梯度，在内外力的驱动下，初始轮廓会自动收敛到颗粒的真实边界。

图10为图3的最终分割结果，图11～13显示了更多的淀粉颗粒显微图像，图14～16依次为图11～13对应的分割结果。

图3～10说明了实例的处理与分割结果。图10为本发明提出的分割方法主要步骤的结果；图3为待分割的淀粉颗粒显微图像；图4为其直方图；图5则为步骤(1)中应用自动阈值法将颗粒与背景分离后得到的二值图像；图6是步骤(3)中的GVF场与求得的颗粒中心；图7为步骤(4)中所述的经计算有效能量后求得的指向各颗粒中心的向量；图8是步骤(4.2)求出的方向梯度；图8为步骤(5)中的初始轮廓；图10是最终分割后的结果。

本例中处理数据集合所设定的参数如下：

步骤(2)中的阈值初始估计值T₀＝图像平均密度；

步骤(2)中的终止值t＝0.05；

步骤(2)中的最大迭代次数N_max＝300；

步骤(3.2)中向量相似性阈值为0.8，即当两向量夹角余弦值大于或等于0.8，则认为二者是一致的；

步骤(4.1)的公式(10)中，k₁＝0.5，k₂＝0.1；

步骤(4.1)的公式(11)中，γ₁＝0.75；

步骤(4.1)的公式(12)中，γ₂＝1.0。

本发明的实现不限于上述实例所公开的范围，可采用不同于上述实例的形式实现上述技术方案。

Claims (2)

1.一种数字图像中颗粒状对象的自动分割方法，其特征在于包括如下步骤：

第一步分离颗粒状对象与背景：先将图像转化为二值图像，再采用自动阈值法将待分割的对象与背景进行分离；

第二步搜索颗粒状对象的中心：计算图像的梯度向量场；求颗粒中心点；

第三步计算方向梯度：先构建一个基于灰度与空间位置的能量函数，计算每个像素对应的能量最大的中心点，作为此像素同质性最强的中心；计算每一个像素与其对应的中心的方向梯度；

第四步初始化圆形轮廓：从每一颗粒的中心出发，分别沿上、下、左、右四个方向，搜索第一个非零值的像素，将中心与该像素的距离，作为初始化圆的半径；

第五步分割颗粒状对象：利用GVF Snake形变模型实施分割，其外力项的梯度采用第三步中的方向梯度；

所述第三步中，定义一个基于灰度与空间位置的有效能量函数：

E_effect=k₁E_pos+k₂E_homo，

其中 $E_{\mathrm{pos}}=e^{{\mathrm{\γ}}_1}\left|\right|P-C\left|\right|,$ $E_{\mathrm{homo}}=e^{-{\mathrm{\γ}}_2|f(i,j)-{\mathrm{\μ}}_k|};$

式中，k₁与k₂为权值参数；

γ₁与γ₂均为大于0的常数；

P为像素；C为颗粒中心；

μ_k为颗粒C_k的灰度均值；

f(x,y)为图像边缘图。

2.根据权利要求1所述一种数字图像中颗粒状对象的自动分割方法，其特征在于第三步中的方向梯度定义为：图像中某一像素点P(p_x,p_y)指向某一中心C(c_x,c_y)的方向梯度向量为d(x,y)=(d_x,d_y)，其中 $d_x=\frac{c_x-p_x}{\left|\right|P-C\left|\right|},$ $d_y=\frac{c_y-p_y}{\left|\right|P-C\left|\right|}.$

Images