CN101968886B - 基于质心跟踪框架的粒子滤波与均值漂移的细胞跟踪方法 - Google Patents

Abstract

本发明公开了基于质心跟踪框架的粒子滤波与均值漂移的细胞跟踪方法，主要解决现有细胞跟踪方法正确率低的问题。其过程为：对视频图像进行二值化分割，提取每个细胞的中心位置；对细胞进行质心跟踪，记录细胞的跟踪轨迹；将轨迹的起始和终止坐标分别记录起始坐标集合和终止坐标集合中，选择待跟踪细胞；采用粒子滤波对待跟踪细胞的轨迹进行一步预测，得到在下一帧图像中的预测坐标点；根据预测的坐标点，适时采用均值漂移方法，选择待跟踪细胞的后续跟踪轨迹；循环预测和选择步骤直到最后一帧图像，完成对所有细胞的跟踪。本发明较其它几种传统的跟踪方法，在跟踪效果及正确率方面都有提高，可用于医学显微视频图像中对运动细胞的分析。

Description

基于质心跟踪框架的粒子滤波与均值漂移的细胞跟踪方法

技术领域

本发明属于数字图像处理技术领域，涉及视频显微图像中运动细胞的自动跟踪，可用于医学显微视频图像中对运动细胞的分析。

背景技术

细胞图像处理作为生物学中的一个研究方向，由于其对细胞研究尤其是癌细胞的研究具有很重要的指导意义，已经成为生物学中的一个新的研究热点。显微技术的发展给观察和研究细胞周期性活动提供了一个很重要的工具。不过，传统的数据人工分析方法对处理这类细胞显微图像数据的作用非常有限，费时费力而且不准确。如今，用于视频显微图像下的细胞自动跟踪系统，作为生物学研究中的一个新兴研究方向，对细胞研究具有重要的指导意义。

细胞跟踪的目标是高速、自动的提供细胞的行为，包括移动、分裂、出现、消失等。广泛的计算机视觉领域的研究已造成强大和灵活的可视化算法。目前，存在的细胞跟踪方法，大致可以分为三类：基于数据关联、单独帧独立分割的方法；基于轮廓变化模型的方法；随机滤波和概率统计的方法。其中：

1.基于数据关联、单独帧独立分割的方法，是细胞跟踪当中比较有效的方法。Al-Kofahi等人在各种匹配假说中使用了线性规则，但是他们的方法不能捕获聚集的细胞以及进入或者离开图像视角的细胞。Dehauwer等人采用了基于欧氏距离匹配的方法，在帧间进行关联。Padfield等人通过跟踪单一细胞训练出的特征，在时空上关联细胞。不过，各种相关的分割匹配技术在关联进入或移出的细胞方面不是很有效，在细胞的分裂跟聚集情况上也都需要特别处理。

2.基于轮廓变化模型的方法，如水平集算法也是常用于细胞跟踪的方法。由于水平集方法可以很容易的在拓扑结构下传递变化信息，为此Yang和Padfield等人采用该方法，将进化水平集用于时空上，能及时有效的关联检测目标。其他相关方法，如Dufour等人提出的水平集方法，采用由先前图像得到的轮廓作为初始，来实时传播变化信息。这种基于轮廓变化模型的方法，虽然能非常有效的处理拓扑结构的变化，但是，在细胞的快速移动或者出现、消失时，都要求重新进行参数的初始化，而大量参数的及时获取和初始化是比较复杂和困难的，无疑会影响该方法的效果。

3.用以跟踪的随机滤波和概率统计的方法，一般情况下都依赖于强有力的模型假设。例如，均值漂移算法可以通过一个基本的外观模型来跟踪目标，这种方法的效果取决于目标定位的精确性，以及模型的实时转移倾向；又如卡尔曼滤波算法，如果该方法采用的状态转移模型可以较好的模拟跟踪目标的运动状态，那么这种滤波方法将非常有效。Kachouie等人提出了一种基于概率准则的最大后验概率，用以细胞追踪，但是由于该算法需要大量的假设，所以该系统能跟踪的细胞数目受到一定限制。事实上，由于细胞自身的非线性运动，不能找到理想的模型去较好的模拟所有细胞的运动，以及视频图像中待跟踪的细胞数目相对较多，增加了跟踪的复杂度，传统的随机滤波和概率统计的方法用于细胞跟踪都不是非常理想。

综上，现有的细胞跟踪方法有如下一些缺点：(1)有些方法不能对细胞运动变化做出准确的判断，如细胞的出现或消失，即使做出判断，对于复杂的运动也不能达到理想的跟踪效果；(2)有些方法对跟踪的目标数目有一定的限制，不能跟踪视频图像中的所有细胞，仅能跟踪一定数目的细胞；(3)有些方法不具有很好的鲁棒性，仅对某些视频图像中具有特定特征或者运动规律的细胞有较为理想的跟踪效果。

发明内容

本发明的目的在于克服上述已有技术的不足，提出一种基于质心跟踪框架的粒子滤波与均值漂移的细胞跟踪方法，以对每个细胞的运动变化做出准确的判断，提供精确的细胞运动轨迹，并且对视频图像中的细胞及跟踪的细胞数目没有限制。

实现本发明目的的技术思路是将粒子滤波与均值漂移跟踪方法引入到基于质心跟踪的框架中，得到跟踪细胞的准确位置。其具体实现步骤包括如下：

(1)对视频的每帧图像运用Otsu方法进行二值化分割，并对二值图像中的细胞区域进行标记，提取每个细胞的中心位置；

(2)根据每个细胞的中心位置，对二值图像中出现的细胞进行质心跟踪，记录得到的细胞跟踪轨迹；

(3)根据每个细胞的跟踪轨迹，判断二值图像是否为细胞跟踪轨迹的起始图像和终止图像，并将起始图像和终止图像中细胞轨迹的起始和终止坐标分别记录于起始坐标集合和终止坐标集合中；

(4)根据轨迹的终止坐标集合和轨迹起始坐标集合，先对视频的相邻两帧二值图像进行判断：如果当前二值图像为细胞跟踪轨迹的终止图像，且下一帧二值图像为细胞跟踪轨迹的起始图像，则把当前二值图像中轨迹终止坐标代表的细胞作为待跟踪细胞；再采用粒子滤波方法，对该待跟踪细胞的轨迹进行一步预测，得到该细胞在下一帧二值图像中的预测坐标点；

(5)根据预测的坐标点，对预测结果进行判断：如果该坐标点在下一帧二值图像中某个细胞区域的内部，且此二值图像是该细胞轨迹的初始图像，则把该细胞和待跟踪细胞确定为同一细胞，并将该细胞的跟踪轨迹作为待跟踪细胞的后续跟踪轨迹；否则，运用均值漂移方法对预测的坐标进行更新，得到新的预测坐标点，再根据新的坐标点进行上述判断，若依旧不满足判断准则，停止对此待跟踪细胞进行跟踪；

(6)循环步骤(4)到步骤(5)直到最后一帧二值图像，完成对所有细胞的跟踪。本发明与现有技术相比具有如下优点：

1)本发明由于在跟踪中采用了粒子滤波与均值漂移方法，对细胞的非线性运动有较好的跟踪有效，而传统的方法仅能有效的跟踪细胞的平稳线性运动；

2)本发明由于将质心跟踪作为整体的跟踪框架，当细胞的运动状态发生变化，如细胞的出现与消失，也能做出准确的判断；

3)本发明由于将质心跟踪作为跟踪第一步，粒子滤波只需要对待跟踪细胞坐标做一步预测，从而解决了常见的粒子退化及运算效率的问题，同时以质心跟踪作为跟踪的基础框架，在粒子滤波预测时，排除了其他细胞对跟踪的干扰，提高了跟踪的正确率。

附图说明

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

图2是本发明仿真使用的2种细胞显微视频图像；

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

图4是本发明仿真实现的某个细胞跟踪轨迹图。

具体实施方式

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

步骤1，对每帧图像进行二值化分割，提取每个细胞的中心位置。

(1a)采用传统的Otsu方法计算原始图像的阈值，对图像进行二值化分割，同时，由于图像噪声的影响，去除二值图像中面积小于20的区域，得到最终二值图像；

(1b)采用4邻域连通标记算法，对最终二值图像中的细胞区域进行标记，提取每个细胞区域的中心，即得到每个细胞的中心位置。

步骤2，根据每个细胞的中心位置，对二值图像中出现的细胞进行质心跟踪，记录得到的细胞跟踪轨迹。

2a)将需要进行质心跟踪的细胞确定为跟踪细胞A；

2b)如果跟踪细胞A的坐标在下一帧二值图像某个细胞B的区域内部，将细胞A和B确定为同一个细胞，并将细胞B的坐标作为跟踪细胞A在下一帧二值图像中的坐标；

2c)将细胞B的坐标加入到跟踪细胞A的跟踪轨迹中，实现对跟踪细胞A的跟踪轨迹更新；

2d)循环步骤2b)到步骤2c)，直到跟踪细胞A的坐标不再在细胞B的区域内部，完成跟踪细胞A的质心跟踪。

步骤3，根据每个细胞的跟踪轨迹，判断二值图像是否为细胞跟踪轨迹的起始图像和终止图像，并将起始图像和终止图像中细胞轨迹的起始和终止坐标分别记录于起始坐标集合和终止坐标集合中。

3a)如果细胞跟踪轨迹的起始坐标在二值图像上，将此二值图像作为细胞跟踪轨迹的起始图像；如果细胞跟踪轨迹的终止坐标在二值图像上，将此二值图像作为细胞跟踪轨迹的终止图像；

3b)将起始图像和终止图像中细胞轨迹的起始和终止坐标分别记录于起始坐标集合和终止坐标集合中；

3c)循环步骤3a)到步骤3b)，直到最后一帧二值图像。

步骤4，选择待跟踪细胞，采用粒子滤波方法，对待跟踪细胞的轨迹进行一步预测，得到该细胞在下一帧二值图像中的预测坐标点。

细胞的运动有时会超出质心跟踪的范围，导致质心跟踪中断，致使本来完整的细胞运动被跟踪成几个平稳的轨迹，这些跟踪轨迹由质心跟踪得来。因此需要将属于同一细胞的跟踪轨迹确定，将其按顺序连接，以形成对细胞的完整跟踪。

根据轨迹的终止坐标集合和轨迹起始坐标集合，先对视频的相邻两帧二值图像进行判断：如果当前二值图像为某些细胞跟踪轨迹的终止图像，且下一帧二值图像为某些细胞跟踪轨迹的起始图像，说明在当前二值图像中，跟踪轨迹终止的细胞，有可能在下一帧图像中消失，也有可能因为细胞的运动变化较大，导致质心跟踪中断，需要与下一帧图像中起始的跟踪轨迹进行匹配，因此，把当前二值图像中轨迹终止坐标代表的细胞作为待跟踪细胞，采用粒子滤波算法对该细胞进行一步预测，具体操作步骤如下；

4a)通过如下方程，初始化粒子集合

$X_k^n=X_k+\frac{3}{2}\sqrt{2}\·\mathrm{randn},(n=1...N)$

其中，

为当前图像中第n个粒子的坐标，X_k为待跟踪细胞的轨迹终止点坐标，randn为随机参数，随机波动范围由图像中细胞运动剧烈程度确定，k表示当前图像在视频中的帧号，N为粒子数，取为30；

4b)对初始化后的粒子集合

采用以下两种状态转移方程进行预测，得到新的粒子集合

当待跟踪细胞的轨迹只有当前终止图像中的坐标X_k时，采用的状态转移方程为：

$X_{k+1}^n=X_k^n+V$

其中，

为第n个粒子的新坐标，V为协方差矩阵，取为

当待跟踪细胞的轨迹包含至少两帧图像的坐标X_k和X_k-1时，采用的状态转移方程为：

$X_{k+1}^n=2\·X_k^n-X_{k-1}+Q$

其中，X_k-1为上一帧图像中细胞坐标，

为第n个粒子在下一帧图像中的新坐标，Q为协方差矩阵，取为

4c)以坐标X_k和为中心分别建立窗宽为H的矩形窗口，设定H为30，计算两个矩形窗的直方图估计概率密度分布，将两个直方图估计概率密度分布的Bhattacharyya距离D_n作为粒子的权值，其公式如下：

$D_n={(1-\underset{u=1}{\overset{B}{\mathrm{\Σ}}}\sqrt{p_n^{\left(u\right)}{q}_n^{\left(u\right)}})}^{1/2}$

其中，B为图像的灰度级数，分别为两个矩形窗的直方图估计概率密度分布；

4d)利用粒子权值，预测待跟踪细胞在下一帧图像中的坐标X_k+1，其公式如下：

$X_{k+1}=\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{D}_n{X}_{k+1}^n.$

步骤5，根据预测的坐标点，对预测结果进行判断，将满足条件的跟踪轨迹连接，以实现对待跟踪细胞的跟踪轨迹更新。

本步骤的具体实现分为三种情况：

5a)如果粒子滤波预测的坐标点在下一帧二值图像中某个细胞区域的内部，且此二值图像是该细胞轨迹的初始图像，根据质心跟踪原理，确定该细胞和待跟踪细胞为同一细胞，并确定该细胞的跟踪轨迹为待跟踪细胞的后续跟踪轨迹，将该细胞的跟踪轨迹连接到待跟踪细胞的跟踪轨迹，实现对待跟踪细胞的跟踪轨迹更新；

5b)如果粒子滤波预测的坐标点不在下一帧二值图像中某个细胞区域的内部，运用均值漂移方法对粒子滤波预测的坐标进行更新，得到新的预测坐标点：

(5b1)将预测的坐标X_k+1作为最初的初始样本x；

(5b2)以初始样本x为中心建立窗宽半径为h的矩形窗口，h取值为15；

(5b3)计算坐标样本均值m_h(x)：

$m_h\left(x\right)=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}G\left(\frac{x_i-x}{h}\right)w\left(x_i\right)x_i}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}G\left(\frac{x_i-x}{h}\right)w\left(x_i\right)}$

其中，h是矩形窗半径，x为初始坐标，x_i为矩形窗中各点的坐标，w(x_i)为坐标点x_i在图像中的灰度值，G(X)选取高斯核函数，其表达式为 $G\left(\frac{x_i-x}{h}\right)=\exp (-\frac{1}{2}{\left|\right|\frac{x_i-x}{h}\left|\right|}^2);$

(5b4)计算坐标样本均值m_h(x)和x的差值，如果‖m_h(x)-x‖大于容许误差ε，将m_h(x)的值赋给x，ε取值为0.1；

(5b5)循环步骤(5b2)到(5b4)，直到‖m_h(x)-x‖小于容许误差ε，所得x即为更新的预测坐标点。

5c)根据均值漂移得到的新的坐标点，进行步骤5a)中的操作，若新的坐标点依旧不在下一帧二值图像中某个细胞区域的内部，停止对此待跟踪细胞进行跟踪，判定该细胞在下一帧图像中消失。

步骤6，循环步骤(4)到步骤(5)直到最后一帧二值图像，确定属于同一细胞的跟踪轨迹，将其按顺序连接，形成完整的跟踪轨迹，以完成对所有细胞的跟踪。

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

1.仿真内容

应用本发明跟踪方法、基于质心跟踪框架的卡尔曼滤波和均值漂移的跟踪方法与基于卡尔曼滤波和均值漂移的跟踪方法，分别对两种显微视频图像中的细胞进行跟踪，如图2所示，其中图2中的图像序列信息如下：

图像每帧大小：图2(a)515像素*650像素图2(b)515像素*650像素

图像序列数目：图2(a)30帧图2(b)200帧

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

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

软件平台为：MATLAB 7.5

仿真内容1.将本发明跟踪方法应用于两种显微视频图像进行跟踪；

仿真内容2.将基于质心跟踪框架的卡尔曼滤波和均值漂移的跟踪方法应用于两种显微视频图像进行跟踪；

仿真内容3.将基于卡尔曼滤波和均值漂移的跟踪方法应用于两种显微视频图像进行跟踪。

2.仿真结果与分析

利用Otsu方法对图2(a)和图2(b)两种显微细胞图像进行二值化分割，得到的二值图像，如图3(a)和图3(b)所示。

利用如图3(b)所示的二值分割后的图像序列，对某个细胞进行跟踪的效果，如图4所示。

从图4可以看出，细胞的运动有非常明显的随机性，这是细胞跟踪的难点，而本发明方法对细胞无规律的运动能够做出准确的判断，即使细胞运动发生突变，也有较好的跟踪效果。

表1列出了仿真内容(1)、(2)、(3)仿真所得的跟踪对比结果。

表1：不同跟踪方法跟踪结果对比

从表1中可以看出，本发明方法能够对视频显微图像中的所有运动细胞进行自动跟踪，并取得较理想的跟踪正确率；相比于其他方法，本发明方法对简单或复杂的细胞视频显微图像都有比较稳健的跟踪效果，鲁棒性较高。

Claims (3)

1.一种基于质心跟踪框架的粒子滤波与均值漂移的细胞跟踪方法，包括如下步骤：

(1)对视频的每帧图像运用Otsu方法进行二值化分割，并对二值图像中的细胞区域进行标记，提取每个细胞的中心位置；

(2)根据每个细胞的中心位置，对二值图像中出现的细胞进行质心跟踪，记录得到的细胞跟踪轨迹；

(3)根据每个细胞的跟踪轨迹，判断二值图像是否为细胞跟踪轨迹的起始图像和终止图像，并将起始图像和终止图像中细胞轨迹的起始和终止坐标分别记录于起始坐标集合和终止坐标集合中；

(4)根据轨迹的终止坐标集合和轨迹起始坐标集合，先对视频的相邻两帧二值图像进行判断：如果当前二值图像为细胞跟踪轨迹的终止图像，且下一帧二值图像为细胞跟踪轨迹的起始图像，则把当前二值图像中轨迹终止坐标代表的细胞作为待跟踪细胞；

(5)再采用粒子滤波方法，对该待跟踪细胞的轨迹进行一步预测，得到该细胞在下一帧二值图像中的预测坐标点：

5a)通过如下方程，初始化粒子集合{

n＝1...N}

$X_k^n=X_k+\frac{3}{2}\sqrt{2}\·\mathrm{randn},,$ n＝1...N，

其中，

为当前图像中第n个粒子的坐标，X_k为待跟踪细胞的轨迹终止坐标，randn为随机参数，N为粒子数，k表示当前图像在视频中的帧号；

5b)对初始化后的粒子集合{

n＝1...N}采用以下两种状态转移方程进行预测，得到新的粒子集合{n＝1...N}：

当待跟踪细胞的轨迹只有当前终止图像中的坐标X_k时，采用的状态转移方程为：

$X_{k+1}^n=X_k^n+V$

其中，

为第n个粒子的新坐标，V为状态协方差矩阵；

当待跟踪细胞的轨迹包含至少两帧图像的坐标X_k和X_k-1时，采用的状态转移方程为：

$X_{k+1}^n=2\·X_k^n-X_{k-1}+Q$

其中，X_k-1为上一帧图像中细胞坐标，

为第n个粒子在下一帧图像中的新坐标，Q为状态协方差矩阵；

5c)以坐标X_k和

为中心分别建立窗宽为H的矩形窗口，计算两个矩形窗的直方图估计概率密度分布，将两个直方图估计概率密度分布的Bhattacharyya距离D_n作为粒子的权值，其公式如下：

$D_n={(1-\underset{u=1}{\overset{B}{\mathrm{\Σ}}}\sqrt{p_n^{\left(u\right)}{q}_n^{\left(u\right)}})}^{1/2}$

其中，B为图像的灰度级数，

分别为两个矩形窗的直方图估计概率密度分布；

5d)利用粒子权值，预测待跟踪细胞在下一帧图像中的坐标X_k+1，其公式如下：

$X_{k+1}=\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{D}_n{X}_{k+1}^n;$

(6)根据预测的坐标点，对预测结果进行判断：如果该坐标点在下一帧二值图像中某个细胞区域的内部，且此二值图像是该细胞轨迹的初始图像，则把该细胞和待跟踪细胞确定为同一细胞，并将该细胞的跟踪轨迹作为待跟踪细胞的后续跟踪轨迹；否则，运用均值漂移方法对预测的坐标进行更新，得到新的预测坐标点，再根据新的坐标点进行上述判断，若依旧不满足判断准则，停止对此待跟踪细胞进行跟踪，

所述运用均值漂移方法对预测的坐标进行更新，具体步骤如下：

6a)将预测的坐标X_k+1作为最初的初始样本x；

6b)以初始样本x为中心建立窗宽半径为h的矩形窗口；

6c)计算坐标样本均值m_h(x)：

$m_h\left(x\right)=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}G\left(\frac{x_i-x}{h}\right)w\left(x_i\right)x_i}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}G\left(\frac{x_i-x}{h}\right)w\left(x_i\right)}$

其中，h是矩形窗半径，x为初始坐标，x_i为矩形窗中各点的坐标，w(x_i)为坐标点x_i在图像中的灰度值，G(X)选取高斯核函数，其表达式为 $G\left(\frac{x_i-x}{h}\right)=\exp (-\frac{1}{2}{\left|\right|\frac{x_i-x}{h}\left|\right|}^2);$

6d)计算坐标样本均值m_h(x)和x的差值，如果||m_h(x)-x||大于容许误差ε，将m_h(x)的值赋给x；

6e)循环步骤6b)到6d)，直到||m_h(x)-x||小于容许误差ε，所得x即为更新的预测坐标点；

(7)循环步骤(4)到步骤(6)直到最后一帧二值图像，完成对所有细胞的跟踪。

2.根据权利要求1所述的细胞跟踪方法，其中步骤(2)所述的对二值图像中出现的细胞进行质心跟踪，按如下步骤进行：

2a)将需要进行质心跟踪的细胞确定为跟踪细胞A；

2b)如果跟踪细胞A的坐标在下一帧二值图像某个细胞B的区域内部，将细胞A和B确定为同一个细胞，并将细胞B的坐标作为跟踪细胞A在下一帧二值图像中的坐标；

2c)更新跟踪细胞A的跟踪轨迹，即将细胞B的坐标加入到跟踪细胞A的跟踪轨迹中；

2d)循环步骤2b)到步骤2c)，直到跟踪细胞A的坐标不再满足2b)中的条件，完成跟踪细胞A的质心跟踪。

3.根据权利要求1所述的细胞跟踪方法，其中步骤(3)所述的判断二值图像是否为细胞跟踪轨迹的起始图像和终止图像，按如下准则判断：

如果细胞跟踪轨迹的起始坐标在二值图像上，将此二值图像作为细胞跟踪轨迹的起始图像；

如果细胞跟踪轨迹的终止坐标在二值图像上，将此二值图像作为细胞跟踪轨迹的终止图像。

Images