专利名称：一种视频显微图像细胞自动跟踪方法
技术领域：
本发明属于细胞生物学与生物制药领域，涉及一种视频显微图像细胞自动跟踪的方法。

背景技术：
细胞跟踪是指在观测时间内对特定细胞的运动轨迹、运动速度、颜色、形状等属性进行定性或定量的分析，它是进行细胞活性、细胞迁徙和细胞趋向性等细胞生物学和生物制药研究的有效方法和必备手段，在生物学、药理学和病理学方面都具有十分重要的研究意义和实用价值。
目前国内细胞跟踪的实现方法主要还是显微设备辅助下的人工定时检测、记录，其自动化程度较低。一方面要求工作人员付出辛苦的劳动，具有较高的劳动强度；另一方面由于细胞运动形式的复杂多变性(如细胞分裂、结合、聚集、交叉)和视频图像本身的噪声和干扰等因素，使得人工观测过程更加困难，并且有较高的出错率。
近几年，国外研究机构在细胞自动跟踪方面取得了一些研究成果，借助图像分析软件研制出一批细胞自动跟踪系统(如University of Virginia、University of Aberdeen和European Molecular Biology Laboratory等)，但在商业化方面还不成熟。这些设备虽然能够通过图像处理和分析软件计算出细胞在图像上的位置，但对于细胞的分裂与结合、细胞聚集和细胞群的分离、多细胞运动轨迹的交叉等复杂现象难以处理。


发明内容
本发明的目的在于针对现有细胞跟方法的不足，为细胞生物学和生物制药等领域研究提供一种自动化程度高、处理能力强的视频显微图像细胞跟踪的方法。
本发明包括以下步骤1.通过相差显微镜获取实时细胞运动图象，将获取的细胞运动视频显微图像逐帧进行增强处理；2.从增强处理后的细胞运动图像中提取目标细胞；3.建立细胞运动动力学模型；4、跟踪运动的目标细胞。
步骤1中视频显微图像的增强处理采用广义模糊增强处理方法，具体步骤如下 1)初始化，输入待处理图像f，设定迭代次数的初始值为r＝1； 2)采用中值滤波器进行平滑滤波，以图像f的每个像素f(i，j)作为窗口(能够覆盖一定图像的单元)的中心，将窗口所覆盖图象的像素灰度平均值作为像素f(i，j)的新的灰度值fij； 3)将新的灰度值fij进行第r次迭代，确定中值滤波后图像的模糊特征平面{μij(r)}； 4)对μij(r)作如下非线性变换，变换结果记为μij′(r)； 5)对μij′(r)作逆变换，得到新的灰度图像{fij′}； 6)灰度变换，设定适当的fmine，fmaxe，计算出广义模糊增强后的图像灰度{fije}； 7)比较第r次和第r-1次增强图像的图像质量评价指标，如果σw(r)小于σw(r-1)，则令r+1r，并返回4)对μij(r)进行迭代计算；否则输出第r-1次增强的图像。
以上算法的计算复杂度可分析如下对于具有m×n个像素的图像，广义模糊增强算法的计算复杂度为O(m×n)。
步骤2中提取目标细胞的方法采用主动轮廓线模型方法(Snake方法)，具体实现包括以下步骤 1)构造能量模型 设定曲线v(s)＝[x(s)y(s)]，s∈
，定义其上的总的能量表示为 其中Eint(v(s))＝(α(s)|vs(s)|2+β(s)|vss(s)|2)(2) Eimage(v(s))＝wlineEline(v(s))+wedgeEedge(v(s))+wtermEterm(v(s))(3) Econ(v(s))＝k(x1-x2)2(4) Eint(v(s))为内部能量，表达了驱使曲线更平滑的力，其中一阶项表达了使相邻点距离更小的拉力，二阶项表达了抵御弯曲的刚性力；α(s)和β(s)表示各自的权重；Eimage(v(s))为图像能量，是从图像得到的引导蛇朝向低灰度或高灰度位置前进的线能量Eline＝I(x，y)、边缘能量用Eedge＝-|▽I(v(s))|和图像中线的终止点和拐角对轮廓线走向的影响的能量Eterm的三个能量项的加权和；wline、wedge和wterm代表图像能量各分量的权重；Econ(v(s))表示吸引轮廓线到图象位置的弹性力，x1和x2分别表示轮廓线和图像位置的指定点。若将外部能量定义为 Eext(v(s))＝Eimage(v(s))+Econ(v(s)) (5) 则总的能量为 2)利用变分法对总能量进行极小化，使轮廓线满足 3)通过曲线包围区域的中心位置确定细胞目标在图像上的位置，该位置作为当前帧该细胞图像的位置测量。
步骤3中建立细胞运动动力学模型包括建立目标运动模型和测量模型，其中目标运动模型为 x(k+1)＝F(k)x(k)+G(k)u(k)+v(k)(8) 测量模型为 z(k)＝H(k)x(k)+w(k)(9) x(k)表示目标细胞的在k时刻的运动状态(位置、速度)，z(k)表示k时刻的图像测量，F(k)、G(k)和H(k)分别表示k时刻的状态转移矩阵、控制输入矩阵和测量矩阵，v(k)和w(k)分别描述了随机系统噪声和测量噪声。
步骤4中运动目标细胞的跟踪方法采用递推Bayesian滤波方法更新每个目标，获得每个目标的当前状态和估计精度，对于多细胞跟踪和细胞分裂、细胞聚集通过数据关联进行处理。
递推Bayesian滤波方法和联合概率数据关联方法均采用成熟的现有技术，如递推Bayesian滤波方法可采用α-β滤波、Kalman滤波、PF滤波等方法，数据关联采用联合概率数据关联方法(由Bar Shalom Y提出)。
本发明根据视频显微图像具有较强噪声和扰动、图像对比度差的特点，采用中值滤波、模糊增强、灰度变换和图像质量等评估环节构成的广义模糊增强方法对原始的图像数据逐帧处理，能较大程度提高图像处理效果。提出的广义模糊增强算法较通常的图像增强方法在处理弱对比度、强噪声图像方面，具有更好的性能。
本发明从获取的原始视频显微图像入手，运用图像处理和分析方法改善图像质量、提取细胞轮廓和确定细胞位置，从而得到时序观测数据。在分析细胞运动特性和图像中噪声与干扰特征的基础上，采用动态系统建模和随机建模方法，建立了细胞运动系统模型，通过递推Bayesian滤波和数据关联技术，对细胞的运动轨迹进行自动跟踪，同时给出了细胞复杂运动的跟踪结果。本发明具有完整的系统性和很强的实用性。

具体实施例方式 视频显微图像细胞自动跟踪的方法包括以下步骤1.通过相差显微镜获取实时细胞运动图象，将获取的细胞运动视频显微图像逐帧进行增强处理；2.从增强处理后的细胞运动图像中提取目标细胞；3.建立细胞运动动力学模型；4、跟踪运动的目标细胞。
步骤1中视频显微图像的增强处理采用广义模糊增强处理方法，具体步骤如下 (1)为了减弱图像中的颗粒噪声，一般情况下，在空间域内可以用邻域平均来减少噪声，在频率域，因为噪声频谱多在高频段，因此可以采用各种形式的低通滤波方法来减少噪声。为了尽量避免图像边缘模糊又能去除脉冲噪声和所谓“椒盐”噪声(Salt-and-pepper noise)，采用中值滤波器进行平滑滤波； (2)采用如下变换方法进行图像增强 其中Fc是模糊特性参数，fij、fmax、fmin分别表示图象中象素(i，j)的灰度值、图象最大灰度和最小灰度值，μij表示象素(i，j)模糊隶属度。运用迭代运算对图像模糊特征平面{μij}进行增强处理，获得新的模糊特征平面{μij′}，并在此基础上进行如下逆变换，得到模糊增强后图像的灰度值 (3)对模糊增强后的图像进行如下灰度变换 式中fije为经过灰度变换t(·)后的图像灰度值，fmine，fmaxe分别为设定的灰度变换后的图像灰度的最大、最小值，fmin′，fmax′分别为模糊增强图像灰度的最大、最小值，且 (4)利用图像灰度对其灰度直方图的标准差加权，得到如下图像质量评价指标 其中σw是图像灰度直方图的加权标准差，pj为第j灰度等级的象素数量在图像总象素N中所占百分比，

是pj的平均值，Δf为图像灰度范围。
(5)重复(2)-(4)步，直到图象质量评价指标σw不再减小为止。
步骤2中提取目标细胞的方法采用主动轮廓线模型方法(Snake方法)，该方法是一种自上而下定位图像特征的机制，首先设定一个初始的轮廓线(“蛇”)，然后通过作用在“蛇点”上的约束力推动轮廓线向图像特征方向前进，最终锁定目标结构是通过极小化动态轮廓线总体能量的积分度量首先的，具体实现包括以下步骤 (1)能量模型的构造 令曲线v(s)＝[x(s)y(s)]，s∈
，则定义在其上的总的能量可以表示为 其中Eint(v(s))＝(a(s)|v(s)|2+β(s)|vss(s)|2) (2) Eimage(v(s))＝wlineEline(v(s))+wedgeEedge(v(s))+wtermEterm(v(s))(3) Econ(v(s))＝k(x1-x2)2 (4) 内部能量Eint(v(s))表达了驱使曲线更平滑的力，其中一阶项表达了使相邻点距离更小的拉力，二阶项表达了抵御弯曲的刚性力，α(s)和β(s)代表各自的权重。图像能量Eimage(v(s))是从图像得到的三个能量项的加权和引导蛇朝向低灰度或高灰度位置前进的线能量，Eline＝I(x，y)；边缘能量用Eedge＝-|▽I(v(s))|表示，从而吸引轮廓线到有着高梯度模值的图像边缘点上；Eterm表示图像中线的终止点和拐角对轮廓线走向的影响。wline、wedge和wterm代表图像能量各分量的权重。Econ(v(s))表达了吸引轮廓线到指定某个位置的弹性力，x1和x2分别表示轮廓线和图像位置的指定点。若记 Eext(v(s))＝Eimage(v(s))+Econ(v(s))(5) 则有 (2)基于变分法的能量极小化计算 轮廓线的最终位置可以通过变分方法获得。将(6)式中的积分项用F(s，vs，vss)代替，则推导的曲线方程应该满足如下欧拉-拉个朗日方程 (3)在主动轮廓线模型方法处理的结果上，可以通过计算曲线包围区域的中心位置来确定细胞目标在图像上的位置。该位置可以作为当前帧图像某细胞的位置测量。
由于细胞运动具有一定的随机性，因此对其运动描述步骤3中建立细胞运动动力学模型采用随机建模的方法，包括建立目标运动模型和测量模型，其中目标运动模型为 x(k+1)＝F(k)x(k)+G(k)u(k)+v(k)(8) 测量模型为 z(k)＝H(k)x(k)+w(k) (9) x(k)表示目标细胞的在k时刻的运动状态(位置、速度)，z(k)表示k时刻的图像测量，F(k)、G(k)和H(k)分别表示k时刻的状态转移矩阵、控制输入矩阵和测量矩阵，v(k)和w(k)分别描述了随机系统噪声和测量噪声。
步骤4中运动目标细胞的跟踪方法采用递推Bayesian滤波方法更新每个目标，获得每个目标的当前状态和估计精度，对于多细胞跟踪和细胞分裂、细胞聚集通过数据关联进行处理。
同一帧图像上往往有多个细胞，这些细胞的运动是无规律的，并且还可能存在细胞的分裂和合并等生物现象。由于图像采集设备性能和实际环境的干扰，图像中的噪声和杂波等级较高。这些因素使得基于图像的细胞跟踪必须解决一个关键技术问题数据关联，即正确判断多个测量信息与多条被跟踪目标航迹的关联性。一旦消除了测量源的不确定性，就可以将多目标跟踪问题转化为多个单目标跟踪问题求解。对于杂波可以通过跟踪门技术和随机建模方法进行处理，细胞的分裂和合并现象可等价为目标航迹的生成和合并问题。
JPDA(联合概率数据关联)是一种解决数据关联问题的次优算法。当跟踪测量空间内接近的多个目标时，测量可能同时落在几个目标的跟踪门内，将这些跟踪门有交集的目标形成“聚”(cluster)，设聚内目标数为nt，落入这些目标跟踪门的测量数为nm，将聚用一个二元关联逻辑矩阵表示。Ωij＝1，表示第i个测量可能源于第j个目标(j＝0表示该测量为杂波)；反之Ωij＝0表示第i个测量不可能源于第j个目标。满足如下三个约束条件的一种测量和目标之间的可能配对事件称之为一个可行事件x(feasible event) 每个目标最多产生一个测量； 每个测量最多来源于一个目标； 落入某个目标的跟踪门之内的候选测量或者源于该目标，或者源于杂波，或者源于其他目标。
将可行事件用二元关联逻辑矩阵Ф表示，由可行事件的约束条件可知，Ф的每行元素之和等于1，每列元素之和等于1或者0(第0列除外)。可行事件可以看作是在测量集和待定目标集形成的所有数学组合中按照三个约束条件遴选出的部分组合。求出每个可行事件的后验概率，并将所有的Фij＝1可行事件的后验概率相加，即第i个测量属于第j个目标的后验概率，用它求出该测量更新该目标时的权重。记表示可行事件x中第i个测量的指标函数，表示可行事件x中第j个目标的指标函数,

表示可行事件x中杂波个数。在源自目标的测量满足正态分布，杂波满足均匀分布，杂波个数满足泊松分布的前提下，得可行事件x发生的后验概率为
其中第一个连乘表示测量属于实际目标的正态分布概率，第二个连乘表示所有目标被检测的概率，第三个连乘表示目标都没有被检测到的概率，c是归一化因子。第i个测量源自第j个目标的后验概率为 第j个目标没有产生任何测量的后验概率为 将所有有效测量加权得融合测量 按照递推Bayesian滤波方法(如α-β滤波、Kalman滤波、PF滤波等方法)更新每个目标，获得每个目标的当前状态(估计量)和估计精度(目标状态协方差矩阵)。
权利要求
1.一种视频显微图像细胞自动跟踪方法，其特征在于该方法包括以下步骤(1)通过相差显微镜获取实时细胞运动图象，将获取的细胞运动视频显微图像逐帧进行增强处理；(2)从增强处理后的细胞运动图像中提取目标细胞；(3)建立细胞运动动力学模型；(4)跟踪运动的目标细胞；
步骤(1)中视频显微图像的增强处理采用广义模糊增强处理方法，具体包括以下步骤
①初始化，输入待处理图像f，设定迭代次数的初始值为r＝1；
②采用中值滤波器进行平滑滤波，以图像f的每个像素f(i，j)作为窗口的中心，将窗口所覆盖图象的像素灰度平均值作为f(i，j)的新的灰度值fij；
③将新的灰度值fij进行第r次迭代，确定中值滤波后图像的模糊特征平面{μij(r)}；
④对μij(r)作如下非线性变换，变换结果记为μij′(r)；
⑤对μij′(r)作逆变换，得到新的灰度图像{fij′}；
⑥对模糊增强后的图像进行如下灰度变换
式中fije为经过灰度变换t(·)后的图像灰度值，fmine，fmaxe分别为设定的灰度变换后的图像灰度的最大、最小值，fmin′，fmax′分别为模糊增强图像灰度的最大、最小值，且
⑦比较第r次和第r-1次增强图像的图像质量评价指标，如果σw(r)小于σw(r-1)，则令r+1r，并返回④对μij(r)进行迭代计算；否则输出第r-1次增强的图像；
步骤(2)中提取目标细胞的方法采用主动轮廓线模型方法，具体实现包括以下步骤
①构造能量模型
设定曲线v(s)＝[x(s) y(s)]，s∈
，定义其上的总的能量表示为
Etotal(v(s))＝∫s(Eint(v(s))+Eimage(v(s))+Econ(v(s)))ds(3)
其中
Eint(v(s))＝(α(s)|vs(s)|2+β(s)|vss(s)|2) (4)
Eimage(v(s))＝wlineEline(v(s))+wedgeEedge(v(s))+wtermEterm(v(s))(5)
Econ(v(s))＝k(x1-x2)2 (6)
Eint(v(s))为内部能量，表达了驱使曲线更平滑的力，其中一阶项表达了使相邻点距离更小的拉力，二阶项表达了抵御弯曲的刚性力；α(s)和β(s)表示各自的权重；Eimage(v(s))为图像能量，是从图像得到的引导蛇朝向低灰度或高灰度位置前进的线能量Eline＝I(x，y)、边缘能量用Eedge＝-|I(v(s))|和图像中线的终止点和拐角对轮廓线走向的影响的能量Eterm的三个能量项的加权和；wline、wedge和wterm代表图像能量各分量的权重；Econ(v(s))表示吸引轮廓线到图象位置的弹性力，x1和x2分别表示轮廓线和图像位置的指定点；
②利用变分法对总能量进行极小化，使轮廓线满足
③通过曲线包围区域的中心位置确定细胞目标在图像上的位置，该位置作为当前帧该细胞图像的位置测量；
步骤(3)中建立细胞运动动力学模型包括建立目标运动模型和测量模型，其中目标运动模型为
x(k+1)＝F(k)x(k)+G(k)u(k)+v(k)(8)
测量模型为
z(k)＝H(k)x(k)+w(k)(9)
x(k)表示目标细胞的在k时刻的运动状态，z(k)表示k时刻的图像测量，F(k)、G(k)和H(k)分别表示k时刻的状态转移矩阵、控制输入矩阵和测量矩阵，v(k)和w(k)分别描述了随机系统噪声和测量噪声；
步骤(4)中运动目标细胞的跟踪方法采用递推Bayesian滤波方法更新每个目标，获得每个目标的当前状态和估计精度，对于多细胞跟踪和细胞分裂、细胞聚集通过数据关联进行处理。
全文摘要
本发明涉及一种视频显微图像细胞自动跟踪的方法。现有的细胞跟踪方法自动化程度较低，不能适应细胞运动形式的复杂多变性和多细胞运动跟踪。本发明的步骤包括将获取的细胞运动视频显微图像逐帧进行增强处理；从增强处理后的细胞运动图像中提取目标细胞；建立细胞运动动力学模型；跟踪运动的目标细胞。本发明在分析细胞运动特性和图像中噪声与干扰特征的基础上，采用动态系统建模和随机建模方法，通过递推Bayesian滤波和数据关联技术，对细胞的运动轨迹进行跟踪，自动化程度高、处理能力强。
文档编号G01N21/84GK101144784SQ20071007107
公开日2008年3月19日 申请日期2007年9月4日 优先权日2007年9月4日
发明者彭冬亮, 林岳松, 金朝阳, 薛安克, 陈华杰, 朱胜利, 郭云飞 申请人:杭州电子科技大学