专利名称：基于Meanshift算法的视频显微图像细胞自动跟踪方法
技术领域：
本发明属于数字图像处理领域，涉及视频显微图像细胞自动跟踪的方法，该方法
可用于细胞运动分析。
背景技术：
研究细胞的运动十分重要，从细胞学建立初，最开始研究的就是细胞的运动，那时 由于实验设备的制约，人眼能够看到的微观级别只是在显微镜下看到整个细胞，而不像今 天，在高倍显微镜下，细胞的基本机构都能观察的到，所以那时细胞学研究更多的是关注细 胞的运动情况，时过境迁，今天科技的发展使得科研的进行方便而有效，虽然更多的科研活 动已经深入到细胞更微观的领域，可是细胞运动研究的重要性并没有随着时代发展而淡 却，反而以另一种方式活跃于细胞研究的领域。 在一种新药的研制过程中，往往需要在显微镜下对细胞进行注药试验数千次，观
察活体细胞在注药之后的反应，当实验成熟之后才会进一步的在生物体内测试，当研究某
种病毒细胞的扩散方式和程度时，同样会采用这种方式，先要在培养皿中观察病毒细胞入
侵正常细胞的过程，当然还有许多需要细胞运动研究的地方，其应用范围在不断的扩大，而
且新的技术也赋予这项研究以新的内容，很多科研组织都已经开发出基于新方法研究细胞
的新平台，可以预见，细胞运动研究的发展会由于其重要作用而越来越迅速。 目前，国内细胞跟踪的实现方法主要还是显微设备辅助下的人工定时检测和记
录，其自动化程度较低， 一方面要求工作人员付出辛苦的劳动，具有较高的劳动强度；另一
方面由于细胞形式的复杂多变性，如细胞分裂、结合、聚集、交叉和视频图像本身的噪声和
干扰等因素，使得人工观测过程更加困难，并且有较高的出错率。 近几年，国外研究机构在细胞自动跟踪方面取得了一些研究成果，借助图像分析 软件研制出一批细胞自动跟踪系统，如近日，英国一个由生物学家和计算机软件开发专家 组成的研究小组研制出一套可以分析显微镜下观察到的一系列静止细胞数字图像的系统， 美国莱斯大学的研究者在一项突破性研究中，设计新型计算机图像分析系统，能够追踪存 活在生命组织模型中的三维个体癌细胞，但在商业化方面它们都还不成熟，这些设备虽然 能够通过图像处理软件和分析软件计算出细胞在图像上的位置，但对于细胞的分裂、细胞 聚合和多细胞运动轨迹的交叉等特殊情况难以处理。
发明内容
本发明的目的在于针对现有细胞跟踪方法的不足，提出一种基于Meanshift算法 的视频显微图像细胞自动跟踪方法，以实现对运动细胞进行自动跟踪，提高对细胞运动中 出现特殊情况时的处理能力。 实现本发明的技术方案是将细胞运动跟踪问题看作多目标跟踪问题，利用 Meanshift算法对预处理后的细胞二值图像进行聚类，将得到的聚类中心作为待跟踪细胞 中心下一时刻的预测位置，以这些预测位置为中心，为每个待跟踪细胞建立矩形搜索匹配
4区域，计算待跟踪细胞与它的每一个可疑匹配细胞之间的代价函数值，将代价函数值最小 的可疑匹配细胞作为待跟踪细胞的最佳匹配细胞，完成对待跟踪细胞的跟踪，具体实现步 骤如下 1)对细胞原始图像进行预处理，利用分水岭方法对预处理后的图像进行二值分 割，并对分割后的二值图像中的细胞区域进行标记，提取每个细胞的中心位置、周长和面积 信息； 2)为每个待跟踪细胞建立匀速运动模型，并对模型中的状态向量及误差协方差矩 阵进行初始化； 3)利用卡尔曼滤波器对当前时刻所有待跟踪细胞的下一时刻运动进行预测，将获 得的预测值作为初始聚类中心，利用Meanshift算法对下一时刻二值图像进行聚类，并将 得到的聚类中心作为待跟踪细胞中心下一时刻的预测位置，以这些预测位置为中心，为每 个待跟踪细胞建立矩形搜索匹配区域，在下一时刻中，将中心点落入这些区域的细胞称为 待跟踪细胞的可疑匹配细胞； 4)计算待跟踪细胞与它的每一个可疑匹配细胞之间的代价函数值为
cosT (i迈，jn) = a D (i迈，jn) + P P (i迈，jn) + y S (i迈，jn) 1《m《M， 1《n《Nm
其中，im表示第m个待跟踪细胞，jn为待跟踪细胞im的第n个可疑匹配细胞，M为 待跟踪细胞总数，Nm为第m个待跟踪细胞的可疑匹配细胞总数，D(im， jn)为待跟踪细胞im
的预测中心与可疑匹配细胞Jn的中心之间距离，P(im，j》为待跟踪细胞im与可疑匹配细胞 J'n的周长差，S(im， jn)为待跟踪细胞im与可疑匹配细胞Jn的面积差，CI为距离权值，|3为
周长差权值，Y为面积差权值，a+|3 + Y = l，本实例选取a = 0. 5， |3 = 0. 2， Y = 0. 3 ;
5)将代价函数值最小的可疑匹配细胞作为待跟踪细胞的最佳匹配细胞，完成对待 跟踪细胞的跟踪； 6)更新步骤2)中的待跟踪细胞匀速运动模型中的状态向量及误差协方差矩阵；
7)循环步骤3)到步骤6)直到最后一帧图像。
本发明与现有的技术相比具有以下优点 1.由于Meanshift算法具有计算速度快和聚类效果好的优点，本发明利用它对细 胞二值图像进行聚类，能够快速得到较为准确的待跟踪细胞预测位置； 2.本发明计算待跟踪细胞与它的每一个可疑匹配细胞之间的代价函数值，将代 价函数值最小的可疑匹配细胞作为待跟踪细胞的最佳匹配细胞，可以在出现细胞的分裂、 细胞聚合和多细胞运动轨迹的交叉等特殊情况时，仍然能够较准确的完成待跟踪细胞的跟 踪； 3.仿真结果表明，本发明的待跟踪细胞跟踪位置误差在很小的范围内，有较好的 跟踪效果。
图1是本发明的流程示意图； 图2是本发明仿真使用的细胞视频序列原图； 图3是本发明仿真实现的图像预处理后得到的细胞二值图像； 图4是本发明仿真实现的部分细胞跟踪轨迹 图5是本发明仿真实现的部分细胞跟踪位置误差图,
具体实施例方式参照图1 ，本发明的具体实现过程如下 步骤l，对细胞原始图像进行二值分割，提取每个细胞的信息； la)采用0tsu方法计算细胞原始图像的阈值，对细胞原始图像进行阈值处理； lb)在阈值处理后的图像上检测细胞轮廓，将轮廓所在像素点灰度值赋为l，将轮
廓以外像素点灰度值赋为O，将得到结果作为分水岭变换的内部控制标记； lc)将内部控制标记依次进行距离变换和分水岭变换，将得到的结果作为分水岭
的外部控制标记； Id)利用内部控制标记和外部控制标记，对输入图像的梯度图像进行修正；
le)将修正后的梯度图像进行分水岭变换，完成对细胞原始图像的二值分割；
If)对分割后的细胞二值图像中的细胞区域进行标记，采用4邻域连通标记算法， 经过图像标记后，细胞部分就成了一个一个相互分离的区域，对各个细胞的中心位置、周长 和面积信息进行提取。 步骤2，为每个待跟踪细胞建立匀速运动模型，并对模型中的状态向量及误差协方 差矩阵进行初始化； 本匀速运动模型建立在一个二维平面上，以水平方向为X轴，以垂直方向为Y轴， 以图像左下角为原点来建立坐标系，模型中的状态向量定义为X二 (x y Xvy》T，观领U向量 定义为Z = (x y)T，其中(x，y)为运动细胞在图像上的位置坐标，(xv，yv)为运动细胞在水 平方向和垂直方向上的速度。
2a)建立运动模型为
Xk+1 = FXk+wk 其中，Xk是状态向量，F为状态转移矩阵，F =
系统噪声，是一个零均值方差为Q的高斯白噪声；
2b)建立测量模型为
Zk = HkXk+vk 其中，Zk是观测向量，H为量测转移矩阵，/Z:
i r o o
0 10 0
o o i r
0 0 0 1
T为帧间时间间隔，Wk是
10 0 0 0 10 0
vk是量测噪声，是-
-水
零均值方差为R的高斯白噪声； 2c)状态向量及误差协方差矩阵进行初始化；状态向量X二(x y xv yv)T中，x， y初始化为待跟踪细胞在第-
和纵坐标，Xv y,初始化为O;
4贞中的横坐标
误差协方差矩阵初始化为i5 =
画0.22 0 0 0 0 0.22 0 0 0 0 0.22 0 0 0 0 0.22 步骤3，利用卡尔曼滤波器对所有待跟踪细胞的下一时刻运动进行预测，将获得的 预测值作为初始点，利用Meanshift算法对下一时刻二值图像进行聚类，为每个待跟踪细 胞建立矩形搜索匹配区域，具体步骤如下 3a)利用卡尔曼滤波器进行预测， 一步预测值;^w和预测误差的协方差阵4H分
别为
/v,FA,尸+rer丄 其中，f t—,表示在k-1时刻对待跟踪细胞的状态估计，,w一表示k-1时刻误差协
方差矩阵，F为状态转移矩阵，r为噪声矩阵，Q为系统噪声Wk的方差； 3b)以卡尔曼滤波器获得的预测值作为初始聚类中心，利用Meanshift算法对下
一时刻二值图像进行聚类，将得到的聚类中心作为待跟踪细胞中心下一时刻的预测位置；Meanshift算法是1975年由Fukunaga和Hostetler提出的一种基于非参数的核
密度估计理论，它将每个数据点转移到密度函数的局部最大值，即在概率空间中求解概率
密度极值的优化迭代算法，具有计算速度快和聚类效果好的优点，其要点是 假设在d维欧式空间Rd中给定n个数据点{xj i =卜n，定义Meanshift矢量为 =
, =1
一K
-—JC
1 其中，h是核窗半径，x为初始点，Xi为待聚类数据点，G(X)为核函数，公式右边第 一部分称为初始点的样本均值。 Mean shift算法的思想就是在每次的当前点计算ms矢量，然后从当前点漂移ms 矢量到另外一个点，重复此过程直到收敛为止，最终收敛的点即为所要求的聚类中心点。
本发明将Meanshift算法应用在对细胞二值图像的聚类中，具体步骤为
首先，给定一个循环跳出的阈值e ，从卡尔曼滤波器所得的一步预测值出发
3bl)计算所有初始聚类中心x的样本均值mh(x); 附力("=.
,=1 A
其中，h是核窗半径，x为初始聚类中心，Xi为待聚类的数据点，核函数G(X)选取 1
高斯核函数，即G(X卜exp(
2
4) 3b2)若I |mh(X)-x| | > e ，把1%00赋给x，返回2bl);否则，停止迭代，mh(x)即
为聚类中心，将得到的聚类中心作为待跟踪细胞中心下一时刻的预测位置；
7
3c)以待跟踪细胞的预测位置为中心，以L = 30为边长，建立搜索匹配区域，边长 是根据经验知识确定的，在下一时刻中，将中心点落入这些区域的细胞称为待跟踪细胞的 可疑匹配细胞。 步骤4，计算待跟踪细胞与每一个可疑匹配细胞之间的代价函数值； 4a)计算待跟踪细胞im与可疑匹配细胞j的周长差P(im， jn)、面积差S(im， jn)以
及待跟踪细胞im预测中心与可疑匹配细胞jn中心之间距离D(im， jn); 尸(，厶)—^ ，其中4，&分别为待跟踪细胞im和可疑匹配细胞jn的周长； 外 ,人)=|\ —^1，其中^,^分别为待跟踪细胞im和可疑匹配细胞jn的面积； D(U")二^/(气-、)2+Cy,m ->^)2 ，其中气,;^分别为待跟踪细胞im的横坐标和 纵坐标，~ ，力 分别为可疑匹配细胞jn的横坐标和纵坐标； 4b)将4a)中计算的待跟踪细胞im和可疑匹配细胞jn的周长差P (im， jn)、面积差 S(im，jn)以及待跟踪细胞预测中心与可疑匹配细胞中心之间距离D(im，Jn)代入如下公式中 计算待跟踪细胞与可疑匹配细胞之间的代价函数值 cosT (i迈，jn) = a D (i迈，jn) + P P (i迈，jn) + y S (i迈，jn) 1《m《M， 1《n《Nm
其中，im表示第m个待跟踪细胞，jn为待跟踪细胞im的第n个可疑匹配细胞，M为 待跟踪细胞总数，Nm为第m个待跟踪细胞的可疑匹配细胞总数，D(im， jn)为待跟踪细胞im
的预测中心与可疑匹配细胞Jn的中心之间距离，P(im，j》为待跟踪细胞im与可疑匹配细胞 J'n的周长差，S(im， jn)为待跟踪细胞im与可疑匹配细胞Jn的面积差，CI为距离权值，|3为
周长差权值，Y为面积差权值，a+|3 + Y = l，本实例选取a = 0. 5， |3 = 0. 2， Y =0.3。
步骤5，将代价函数值最小的可疑匹配细胞作为待跟踪细胞的最佳匹配细胞，完成 对待跟踪细胞的跟踪； 代价函数反映了可疑匹配细胞与待跟踪细胞之间的差异，可疑匹配细胞的代价函 数越小则说明它与待跟踪细胞匹配成功可能性越大，通过比较代价函数值的大小，将代价 函数值最小的可疑匹配细胞作为待跟踪细胞的最佳匹配细胞，完成对待跟踪细胞的跟踪， 如果某个待跟踪细胞找到最佳匹配细胞，利用最佳匹配细胞的中心位置、周长和面积信息 来更新这个待跟踪细胞的信息，在下一帧来到时，继续对它进行跟踪；如果某个待跟踪细胞 没有搜索到最佳匹配细胞，则认为这个待跟踪细胞消失了 ，下一帧来到后，不再对它进行跟
踪o 步骤6，利用更新公式更新状态向量》4|4和误差协方差矩阵^^ ，分别为 < =》,+ -《_,) & =,i^(/^^/^ '/^w 其中，K为滤波增益，^ = 4^^(/^^/T1 ，Am为一步预测值，4^为预测
误差的协方差阵，zk为新获得的量测值，即细胞中心位置，H为量测矩阵，R为量测噪声V的 协方差矩阵。 步骤7，循环步骤3到步骤6，不断的进行待跟踪细胞从当前时刻到下一时刻的跟 踪，直到跟踪到最后一帧图像，停止跟踪。
本发明的效果可以通过以下仿真进一步说明
1.仿真内容 采用蒙特卡罗仿真实验50次，对本发明进行运动细胞跟踪实验，其中图像序列来
自于显微细胞视频序列，如图2所示。图2中的图像序列信息如下 图像每帧大小650像素X 515像素 图像序列数目74帧 图像帧间时间间隔T = l/24s硬件平台为Intel Core2 Duo CPU E6550@2. 33GHZ、2GB RAM
软件平台为MATLAB 7. 0
2.仿真结果及分析 利用分水岭方法对图2中显微细胞图像进行二值分割，得到的二值图像，如图3所 示。 利用图3所示的二值分割后的图像序列中的待跟踪细胞进行跟踪的效果，如图4 所示。 利用50次蒙特卡罗仿真实验的仿真结果求得待跟踪细胞的位置跟踪误差，如图5 所示。 从图4中可以看出，本发明方法的跟踪轨迹非常接近真实运动轨迹，从图5中可以 看出本发明方法的位置跟踪误差在很小的范围内。仿真结果表明，本发明能够对视频显微 图像中的运动细胞进行自动跟踪，并取得较好的跟踪效果。
权利要求
一种基于Meanshift算法的视频显微图像细胞自动跟踪方法，包括如下步骤1)对细胞原始图像进行预处理，利用分水岭方法对预处理后的图像进行二值分割，并对分割后的二值图像中的细胞区域进行标记，提取每个细胞的中心位置、周长和面积信息；2)为每个待跟踪细胞建立匀速运动模型，并对模型中的状态向量及误差协方差矩阵进行初始化；3)利用卡尔曼滤波器对当前时刻所有待跟踪细胞的下一时刻运动进行预测，将获得的预测值作为初始聚类中心，利用Meanshift算法对下一时刻二值图像进行聚类，并将得到的聚类中心作为待跟踪细胞中心下一时刻的预测位置，以这些预测位置为中心，为每个待跟踪细胞建立矩形搜索匹配区域，在下一时刻中，将中心点落入这些区域的细胞称为待跟踪细胞的可疑匹配细胞；4)计算待跟踪细胞与它的每一个可疑匹配细胞之间的代价函数值为cosT(im，jn)＝αD(im，jn)+βP(im，jn)+γS(im，jn) 1≤m≤M，1≤n ≤Nm其中，im表示第m个待跟踪细胞，jn为待跟踪细胞im的第n个可疑匹配细胞，M为待跟踪细胞总数，Nm为第m个待跟踪细胞的可疑匹配细胞总数，D(im，jn)为待跟踪细胞im的预测中心与可疑匹配细胞jn的中心之间距离，P(im，jn)为待跟踪细胞im与可疑匹配细胞jn的周长差，S(im，jn)为待跟踪细胞im与可疑匹配细胞jn的面积差，α为距离权值，β为周长差权值，γ为面积差权值，α+β+γ＝1，本实例选取α＝0.5，β＝0.2，γ＝0.3；5)将代价函数值最小的可疑匹配细胞作为待跟踪细胞的最佳匹配细胞，完成对待跟踪细胞的跟踪；6)更新步骤2)中的待跟踪细胞匀速运动模型中的状态向量及误差协方差矩阵；7)循环步骤3)到步骤6)直到最后一帧图像。
2. 根据权利要求1所述的细胞自动跟踪方法，其中步骤3)所述的利用Meanshift算法 对下一时刻二值图像进行聚类，具体步骤如下2a)计算所有初始聚类中心x的样本均值mh(x):其中，h是核窗半径，x为初始聚类中心，Xi为待聚类的数据点，核函数G(X)选取高斯核函数，即G(Z)"xp(-l^112);2b)若||mh(X)-x|| > e ，把mh(x)赋给x，返回2bl);否则，停止迭代，nih(x)即为聚类 中心。
3. 根据权利要求1所述的细胞自动跟踪方法，其中步骤4)所述的计算待跟踪细胞im 和可疑匹配细胞jn的周长差P(im， jn)，通过如下公式计算其中化，^分别为待跟踪细胞"和可疑匹配细胞jn的周长。
4. 根据权利要求1所述的细胞自动跟踪方法，其中步骤4)所述的计算待跟踪细胞^一^和可疑匹配细胞jn的面积差S(im， jn)，通过如下公式计算<formula>formula see original document page 3</formula>其中《m，A;分别为待跟踪细胞L和可疑匹配细胞jn的面积。
5.根据权利要求1所述的细胞自动跟踪方法，其中步骤4)所述的计算待跟踪细胞im 预测中心与可疑匹配细胞jn中心之间距离D(im， jn)，通过如下公式计算<formula>formula see original document page 3</formula>其中^，J^分别为待跟踪细胞i,n的横坐标和纵坐标，、，力，分别为可疑匹配细胞l的横 坐标和纵坐标。
全文摘要
本发明公开了一种基于Meanshift算法的视频显微图像细胞自动跟踪方法，主要解决现有技术在出现特殊情况时无法对细胞准确跟踪的问题，其步骤是1)对输入图像进行二值分割，提取细胞信息；2)为待跟踪细胞建立匀速运动模型，并初始化；3)利用Meanshift算法对下一时刻二值图像进行聚类，获得最佳预测位置；4)计算待跟踪细胞与可疑匹配细胞之间的代价函数值；5)将代价函数值最小的可疑匹配细胞作为待跟踪细胞的最佳匹配细胞；6)更新步骤2)中状态向量及误差协方差矩阵；7)循环步骤3)到步骤6)直到最后一帧图像。本发明具有确定细胞预测位置准确，在特殊情况下跟踪效果好的优点，可用于医学图像处理。
文档编号G06T7/20GK101777187SQ20101001372
公开日2010年7月14日 申请日期2010年1月15日 优先权日2010年1月15日
发明者刘芳, 吴建设, 尚荣华, 李阳阳, 焦李成, 王爽, 郑堃, 韩红, 马文萍 申请人:西安电子科技大学