CN101719278A

CN101719278A - 基于khm算法的视频显微图像细胞自动跟踪方法

Info

Publication number: CN101719278A
Application number: CN200910254423A
Authority: CN
Inventors: 韩红; 王爽; 郑堃; 刘芳; 李阳阳; 马文萍; 尚荣华; 吴建设
Original assignee: Xidian University
Current assignee: Xidian University
Priority date: 2009-12-21
Filing date: 2009-12-21
Publication date: 2010-06-02
Anticipated expiration: 2029-12-21
Also published as: CN101719278B

Abstract

本发明公开了一种基于KHM算法的视频显微图像细胞自动跟踪方法，主要解决现有跟踪方法自动化程度较低以及出现特殊情况时无法准确跟踪的不足。其步骤包括：1)利用分水岭方法完成细胞图像的二值分割，提取目标细胞；2)为每个目标细胞建立匀速运动学动力模型，并初始化；3)利用KHM算法对下一时刻二值图像进行聚类，获得所有细胞中心位置；4)通过目标细胞的关联，确定目标细胞在下一时刻的位置，完成运动细胞的跟踪；5)更新步骤2)中目标细胞模型的状态向量及其误差协方差矩阵；6)循环步骤3)到步骤5)直到跟踪到最后一帧图像。本发明具有确定细胞中心位置准确；在特殊情况下跟踪效果好的优点，可用于医学图像处理。

Description

基于KHM算法的视频显微图像细胞自动跟踪方法

技术领域

本发明属于数字图像处理领域，涉及一种视频显微图像细胞自动跟踪的方法，该方法可用于医学图像处理。

背景技术

细胞动态学是众多生物医学处理中重要的一个方面，它直接引导着人类健康的走向。细胞运动分析是细胞动态学研究中重要的组成部分之一。分析细胞运动，通常是对目标细胞进行追踪，进而对细胞的运动速度、位移、轨迹、形状等属性进行定性或定量的分析，它是进行细胞活性、细胞迁徙和细胞趋向性等细胞生物学和生物制药研究的有效方法和必要手段，在生物学、药理学和病理学方面都具有十分重要的研究意义和实用价值。

目前，国内细胞跟踪的实现方法主要还是显微设备辅助下的人工定时检测和记录，其自动化程度较低，一方面要求工作人员付出辛苦的劳动，具有较高的劳动强度；另一方面由于细胞形式的复杂多变性，如细胞分裂、结合、聚集、交叉和视频图像本身的噪声和干扰等因素，使得人工观测过程更加困难，并且有较高的出错率。

近几年，国外研究机构在细胞自动跟踪方面取得了一些研究成果，借助图像分析软件研制出一批细胞自动跟踪系统，如University of Virginia、University ofAberdeen和European Molecular Biology Laboratory等，但在商业化方面还不成熟，这些设备虽然能够通过图像处理软件和分析软件计算出细胞在图像上的位置，但对于细胞的分裂与结合、细胞聚合和多细胞运动轨迹的交叉等特殊情况难以处理。

发明内容

本发明的目的在于针对现有细胞跟踪方法的不足，提出一种基于KHM算法的视频显微图像细胞自动跟踪方法，以实现对运动细胞进行自动跟踪，提高对细胞运动中出现特殊情况时的处理能力。

实现本发明的技术方案是：将细胞运动跟踪问题看作多目标跟踪问题，利用KHM算法对预处理后的细胞二值图像进行聚类，获得所有细胞的位置，以最近邻方法为准则进行目标细胞关联，当关联过程中出现一对多，多对一的特殊情况时，计算不同关联结果的差异代价函数，以差异代价函数最小的关联结果作为最终关联结果，完成运动细胞的自动跟踪，提高对细胞运动中出现特殊情况时的处理能力。具体实现步骤如下：

1)对细胞运动图像进行预处理，利用分水岭方法完成目标细胞与背景以及目标细胞之间的二值分割，并从获得的目标二值图像中提取目标细胞；

2)根据目标细胞运动特点，为每个目标细胞建立匀速运动学动力模型，并对模型中的状态向量及误差协方差矩阵进行初始化；

3)利用卡尔曼滤波对当前时刻所有目标细胞在下一时刻的位置进行预测，以这些预测值作为初始点，利用KHM算法对下一时刻二值图像进行聚类，获得图像中所有细胞中心位置；

4)以最近邻方法为准则进行目标细胞关联，当关联过程中出现一对多，多对一的特殊情况时，计算不同关联结果的差异代价函数，以差异代价函数最小的关联结果作为最终关联结果，确定每个目标细胞在下一时刻的位置，完成运动细胞当前时刻到下一时刻的跟踪；

5)更新步骤2)中的目标细胞模型的状态向量及误差协方差矩阵。

6)循环步骤3)到步骤5)直到跟踪到最后一帧图像。

本发明与现有的技术相比具有以下优点：

1.由于KHM算法具有对初始值不敏感的特性和聚类性能良好的优点，本发明利用它对细胞二值图像进行聚类，能够较为准确的得到细胞中心位置；

2.本发明由于引入差异代价函数，可以在出现细胞的分裂与结合、细胞聚合和多细胞运动轨迹的密集交叉等特殊情况时，能够准确的完成目标细胞的跟踪；

3.仿真结果表明，本发明的目标细胞跟踪位置误差在很小的范围内，有较好的跟踪效果。

附图说明

图1是本发明的流程示意图；

图2是本发明使用的细胞视频序列原图；

图3是本发明实现的图像预处理后得到的细胞二值图像；

图4是本发明实现的部分细胞跟踪效果示意图；

图5是本发明实现的部分细胞跟踪位置误差示意图。

具体实施方式

参照图1，本发明的具体实现过程如下：

步骤1，对细胞运动图像进行预处理，利用分水岭方法完成细胞目标与背景以及细胞目标之间的二值分割，并从获得的二值图像中提取目标细胞。

1a)采用Otsu方法计算图像阈值，对输入图像进行阈值处理，得到初始的目标与背景二值分割结果；

1b)将检测到的目标轮廓作为分水岭的内部控制标记。获取内部控制标记的方法是，对于阈值化方法检测到的初始目标轮廓，将轮廓所在像素点灰度值赋为1，将轮廓以外像素点灰度值赋为0，得到的二值图像，将该二值图像作为分水岭变换的内部控制标记；

1c)将内部控制标记依次进行距离变换和分水岭变换，并将得到的分水线结果作为分水岭的外部控制标记；利用内部控制标记和外部控制标记，对输入图像的梯度图像进行修正；

1d)将修正后的梯度图像进行分水岭变换，检测图像中的目标轮廓，提取目标细胞。

步骤2，根据目标细胞运动特点，为每个目标细胞建立匀速运动学动力模型，并对模型中的状态向量及其误差协方差矩阵进行初始化；

本模型是建立在一个二维平面上，以X轴表示水平方向，以Y轴表示垂直方向，以图像左下角为原点来建立坐标系。由于相邻两幅图像的时间间隔较短，近似认为目标在相邻两幅图像之间是匀速运动，状态向量定义为X＝(x y x_vy_v)^T，量测向量定义为Z＝(x y)^T，x，y分别为目标中心横坐标和纵坐标，x_v，y_v分别为目标在x轴方向和y轴方向的速度，状态转移矩阵和量测矩阵分别定义为：

F = [\begin{matrix} 1 & T & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & T \\ 0 & 0 & 0 & 1 \end{matrix}]

和

H = [\begin{matrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{matrix}]

其中，T为帧间时间间隔。

2a)建立目标运动模型为：

X_i(k+1)＝F_iX_i(k)+Γ_iW_i(k)

其中，X_i(k)是第i个目标的状态向量，F_i和Γ_i分别是第i个目标的状态转移矩阵和噪声矩阵；W_i(k)是与第i个目标关联的系统噪声；

2b)建立测量模型为：

Z_i(k)＝H_iX_i(k)+V_i(k)

其中，Z_i(k)是第i个目标的量测向量，H_i是第i个目标的量测矩阵，V_i(k)是与第i个目标关联的量测噪声。

2c)状态向量及其误差协方差矩阵进行初始化；

状态向量X＝(x y x_v y_v)^T中，x，y初始化为目标细胞在第一帧中的横坐标和纵坐标，x_v y_v初始化为0；

协方差矩阵初始化为

步骤3，利用卡尔曼滤波对当前时刻所有目标细胞在下一时刻的位置进行预测，以这些预测值作为初始点，利用KHM算法对下一时刻二值图像进行聚类，获得图像中所有细胞中心位置；

3a)利用卡尔曼滤波器进行状态估计和预测，一步预测值和预测误差的协方差阵

分别为：

{\hat{X}}_{k | k - 1} = F {\hat{X}}_{k - 1 | k - 1}

{\hat{P}}_{k | k - 1} = F {\hat{P}}_{k - 1 | k - 1} F^{T} + ΓQ Γ^{T}

其中，表示在k-1时刻对目标的状态估计，

表示k-1时刻误差协方差矩阵，F为状态转移矩阵，Γ为噪声矩阵，Q为系统噪声W的方差；

3b)以预测位置为初始聚类中心，利用KHM方法对细胞二值图像进行聚类，获得细胞中心位置；

KHM方法是1999年由Zhang提出的一种基于中心迭代的聚类算法，用所有数据点到每个聚类中心的调和平均值的和作为目标函数，具有对初始值不敏感的特性和聚类性能良好的优点，其要点是：

假设是一组数据元素X＝{X₁，X₂，…X_N}，其中X_i＝{x_i1，x_i2，...，x_im}表示具有m个属性的数据对象，聚类个数为k，KHM算法的目标函数表示为：

KHM (X, C) = Σ_{i = 1}^{N} \frac{k}{Σ_{j = 1}^{k} \frac{1}{d^{2} (X_{i}, C_{l})}}

这里C_l＝{c_l1，c_l2，...，c_lm}是聚类中心，d(X_i，C_l)是距离测度，简写为d_i，l，采用欧式

距离，聚类中心的更新公式如下：

C_{l} = \frac{Σ_{i = 1}^{N} \frac{1}{Σ_{j = 1}^{K} {(\frac{d_{i, l}^{2}}{d_{i, j}^{2}})}^{2}} x_{i}}{Σ_{i = 1}^{N} \frac{1}{Σ_{j = 1}^{K} {(\frac{d_{i, l}^{2}}{d_{i, j}^{2}})}^{2}}}

通过对聚类中心的不断更新，使得目标函数不断减小直到稳定，获得最终的聚类结果。

相对于其它基于中心点迭代的聚类算法，KHM算法用数据点与所有聚类中心的距离的调和平均替代了数据点与聚类中心的最小距离，从而引入了聚类中心对数据点的条件概率和每次迭代过程中数据点的动态加权，实际上起到了将“硬”聚类“软化”的作用。

本发明将KHM算法应用在对细胞二值图像的聚类中，具体步骤为：

3b1)将卡尔曼滤波器所得的一步量测预测值作为初始聚类中心，将细胞二值图像中非零像素点作为待聚类的数据点；

3b2)计算每个聚类中心对二值图像中非零像素点的条件概率P_KHM，以及迭代权函数w：

w (x_{i}) = \frac{Σ_{j = 1}^{K} {| | x_{i} - c_{j} | |}^{- p - 2}}{{(Σ_{j = 1}^{K} {| | x_{i} - c_{j} | |}^{- p})}^{2}}

P_{KHM} (c_{j} / x_{i}) = \frac{{| | x_{i} - c_{j} | |}^{- p - 2}}{Σ_{j = 1}^{K} {| | x_{i} - c_{j} | |}^{- p - 2}}

该权函数w表示每次迭代过程中，数据点对聚类中心的影响程度，该条件概率P_KHM表示数据点归属于聚类中心所代表的类的可能性，其中c_j是聚类中心，x_i是二值图像中非零像素点；

3b3)利用条件概率P_KHM，迭代权函数w，计算新的聚类中心：

c_{j} = \frac{Σ_{i = 1}^{n} P_{KHM} (c_{j} / x_{i}) w (x_{i}) x_{i}}{Σ_{i = 1}^{n} P_{KHM} (c_{j} / x_{i}) w (x_{i})};

3b4)重复步骤3b2)和3b3)，直到收敛或达到最大迭代次数，停止迭代，得到聚类中心，即细胞中心位置。

步骤4，进行目标细胞的关联，确定每个目标细胞在下一时刻的位置，完成运动细胞当前时刻到下一时刻的跟踪；

4a)计算目标细胞一步量测预测值与聚类中心之间的马氏距离，利用最近邻关联准则，用离目标细胞一步量测预测位置最近的聚类中心作为目标细胞的更新位置，完成目标细胞的关联；

4b)若前后两帧图像出现细胞增加，细胞消失，细胞分裂等特殊情况时，前后两帧图像中细胞数目会发生变化，出现一对多，多对一的特殊情况，通过计算不同关联结果的差异代价函数，以差异代价函数最小的关联结果作为最终关联结果，确定目标细胞在当前时刻的位置，差异代价函数可以表示为：

cosT_i，j＝αD(i，j)+βP(i，j)+γS(i，j)

其中，D(i，j)为关联细胞与目标细胞之间马氏距离，P(i，j)为关联细胞与目标细胞之间的周长差，S(i，j)为关联细胞与目标细胞之间面积差，差异代价函数越小，则这两个细胞关联的可能性就越大，α为距离权值，β为周长差权值，γ为面积差权值，三个权值分别取值为：α＝0.5，β＝0.2，γ＝0.3；

步骤5，更新步骤2中的目标细胞模型的状态向量及误差协方差矩阵。

利用更新公式更新状态向量

和误差协方差矩阵

分别为：

{\hat{X}}_{k | k} = {\hat{X}}_{k | k - 1} + K (z_{k} - H {\hat{X}}_{k | k - 1})

{\hat{P}}_{k | k} = {\hat{P}}_{k | k - 1} - {\hat{P}}_{k | k - 1} H^{T} {(H {\hat{P}}_{k | k - 1} H^{T} + R)}^{- 1} H {\hat{P}}_{k | k - 1}

其中，K为滤波增益，

为一步预测值，

为预测误差的协方差阵，z_k为新获得的量测值，即细胞中心位置，H为量测矩阵，R为量测噪声V的协方差矩阵。

步骤6，循环步骤3到步骤5，不断的进行目标细胞从当前时刻到下一时刻的跟踪，直到跟踪目标细胞到最后一帧图像，停止跟踪。

本发明的效果可以通过以下仿真进一步说明

1.仿真内容：

采用蒙特卡罗仿真实验50次，对本发明进行运动细胞跟踪实验，其中图像序列来自于显微细胞视频序列，如图2所示。图2中的图像序列信息如下：

图像每帧大小：650像素X515像素

图像序列数目：74帧

图像帧间时间间隔：T＝1/24s

硬件平台为：Intel Core2 Duo CPU E6550 @ 2.33GHZ、2GB RAM

软件平台为：MATLAB 7.0

2.仿真结果及分析：

利用分水岭方法对图2中显微细胞图像进行二值分割，得到的二值图像，如图3所示。

利用本发明方法对目标细胞跟踪的效果，如图4所示，。

利用50次蒙特卡罗仿真实验的仿真结果求得目标细胞的位置跟踪误差，如图5所示。

从图4中可以看出，本发明方法的跟踪轨迹非常接近目标的真实运动轨迹，从图5中可以看出本发明方法的位置跟踪误差在很小的范围内。仿真结果表明，本发明能够对视频显微图像中的运动细胞进行自动跟踪，并取得较好的跟踪效果。

Claims

1.一种基于KHM算法的视频显微图像细胞自动跟踪方法，包括如下步骤：

6)循环步骤3)到步骤5)直到跟踪到最后一帧图像。

2.根据权利要求1所述的细胞自动跟踪方法，其中步骤3)所述的用KHM算法对下一时刻二值图像进行聚类，按如下步骤进行：

2a)将卡尔曼滤波器所得的一步量测预测值作为初始聚类中心，将细胞二值图像中非零像素点作为待聚类的数据点；

2b)计算每个聚类中心对二值图像中非零像素点的条件概率P_KHM，以及迭代权函数w：

w (x_{i}) = \frac{Σ_{j = 1}^{K} {| | x_{i} - c_{j} | |}^{- p - 2}}{{(Σ_{j = 1}^{K} {| | x_{i} - c_{j} | |}^{- p})}^{2}}

P_{KHM} (c_{j} / x_{i}) = \frac{{| | x_{i} - c_{j} | |}^{- p - 2}}{Σ_{j = 1}^{K} {| | x_{i} - c_{j} | |}^{- p - 2}}

2c)利用条件概率P_KHM，迭代权函数w，计算新的聚类中心：

c_{j} = \frac{Σ_{i = 1}^{n} P_{KHM} (c_{j} / x_{i}) w (x_{i}) x_{i}}{Σ_{i = 1}^{n} P_{KHM} (c_{j} / x_{i}) w (x_{i})}

2d)重复步骤2b)和2c)，直到收敛或达到最大迭代次数，停止迭代，得到聚类中心，即细胞中心位置。

3.根据权利要求1所述的细胞自动跟踪方法，其中步骤4)所述的以最近邻方法为准则进行目标细胞关联，按如下步骤进行：

3a)计算目标细胞一步量测预测值与聚类中心之间的马氏距离，利用最近邻关联准则，用离目标细胞一步量测预测位置最近的聚类中心作为目标细胞的更新位置，完成目标细胞的关联；

3b)若前后两帧图像出现细胞增加，细胞消失，细胞分裂等情况时，前后两帧图像中细胞数目会发生变化，出现一对多，多对一的特殊情况，通过计算不同关联结果的差异代价函数，以差异代价函数最小的关联结果作为最终关联结果，确定目标细胞在当前时刻的位置，差异代价函数可以表示为：

cosT_i，j＝αD(i，j)+βP(i，j)+γS(i，j)

其中，D(i，j)为关联细胞与目标细胞之间马氏距离，P(i，j)为关联细胞与目标细胞之间的周长差，S(i，j)为关联细胞与目标细胞之间面积差，差异代价函数越小，则这两个细胞关联的可能性就越大，α为距离权值，β为周长差权值，γ为面积差权值，三个权值分别取值为：α＝0.5，β＝0.2，γ＝0.3。