一种多细胞相互作用跟踪系统

摘要

本发明提供了一种多细胞相互作用跟踪系统，原始图片输入后，通过三个模块：多细胞自动分割技术模块、细胞交互的多模型粒子滤波跟踪技术模块、通信与身份管理模块的分析，实现多细胞的跟踪；采用本发明提供的方法能够成功分割多细胞图像，特别包含重叠或粘连细胞的图像；对多细胞碰撞、分裂等情形，能实现自动跟踪，并给出所有细胞动力学特征参数（轨迹、速度、角速度），是实现多细胞的自动跟踪。

说明书

技术领域

本发明涉及一种多细胞相互作用跟踪系统，特别是一种多细胞相互作用跟踪系统。

背景技术

多细胞运动分析对于生物学、药理学和病理学方面的研究都有着十分重要的学术价值。 目前细胞跟踪的方法主要分为三类，分别是基于模型调整或变化的跟踪方法、基于检测的方 法和基于贝叶斯概率估计技术。基于模型调整或变化的跟踪方法，该方法主要应用于细胞的 轮廓信息难以提取等场合，该方法将给定的细胞初始外形轮廓，通过参数优化方法使之变形 来拟合图像中的目标对象，细胞的运动轨迹就从不同帧里自动提取出，通常所用的主动轮廓 法（即Snake方法）、水平集方法（Level-Set）和均值偏移法（Mean-Shift）都属于此类跟 踪方法，这些方法都有一个共同特点，即用上一帧所得到的模型参数作为当前帧图像分析的 起始条件或先验信息。基于检测技术的跟踪方法，该方法分为分割与跟踪两步，分割的任务 就是在每帧中提取细胞信息（形状、重心、面积等），跟踪的任务就是将多帧中所提取的细胞 进行关联，找出各细胞的运动轨迹、计算细胞运动速度等。基于贝叶斯概率估计方法利用先 验的一些信息和测量信息来估计细胞状态的后验概率密度函数。由于细胞跟踪问题是一个非 线性、非高斯参数估计问题，因此，所对应的Bayesian滤波没有解析解，然而可用数值近 似的方法进行求解，如粒子滤波器（Particle Filter，有时称Sequential Monte-Carlo）。 Smal（参考文献）将粒子滤波器与交互式多模型IMM和Mean-Shift聚类方法相结合来对细 胞进行运动分析，所得到跟踪性能明显好于人工跟踪方法。在此基础上，Smal（参考文献） 又提出了一种基于粒子滤波器的二维或三维多细胞跟踪方法，在噪声密度比较大的情况下， 所得到的跟踪性能优于确定性的跟踪方法。粒子滤波器也可与水平集或主动轮廓方法相结合， 虽然所得到的跟踪性能要优于确定性跟踪方法，但是它们只能解决单细胞跟踪问题，对于同 时跟踪多个细胞却显得力不从心。Cui（参考文献）提出了一种Monte Carlo跟踪方法，目 标的状态被描述成白细胞的位置，目标的动力学模型是通过学习方法得到的，假设目标的运 动是随机行为，利用以前的细胞运动知识，在预测值附近生成样本，测量信息的获取是基于 局部的图像强度，所得到的跟踪性能要优于Snake跟踪方法。Nhat H.Nguyen（参考文献） 等人提出了一种自动跟踪多个碰撞细胞的方法，首先对每个碰撞状态采用两个单独的运动模 型，这两个模型分别表示碰撞内运动状态和碰撞外运动状态，其次利用状态过渡间的碰撞假 设来处理碰撞中或碰撞外的突然变化，这种方法降低了碰撞细胞跟踪的失误率。

细胞运动分析面临着众多难题，如细胞碰撞、细胞分裂、多细胞动力学特性差异、细胞 重叠、细胞数量变化、低对照比细胞图像序列等等。现有的多细胞跟踪技术还很少能处理这 些情况，跟踪的精确性和稳定性都不足，对于多细胞的动力学特性的差异，其跟踪性能未加 以研究等等。本发明旨在解决当细胞动力学特性存在差异、细胞数目时变、细胞碰撞、分裂 等跟踪难题，通过混合检测方法获得细胞特征信息，利用细胞交互的多模型粒子滤波跟踪方 法进行跟踪，根据细胞特征差异进行通信与身份管理，从而实现多细胞自动跟踪。

发明内容

本发明要解决的技术问题是提供一种鲁棒的、精确的、面向多细胞碰撞、分裂情形且能 估计多细胞动力学特性的多细胞相互作用跟踪系统，本发明利用混合检测方法分割包含重叠 或粘连细胞的图像后,利用细胞交互的多模型粒子滤波跟踪方法进行多细胞运动分析，为了减 小关联的不确定性，利用细胞特征差异进行通信与身份管理，从而实现细胞碰撞、细胞分裂、 细胞数量时变的多细胞自动跟踪，给出所有细胞动力学特征参数（轨迹、速度、角速度等）。

原始图片输入后，通过该系统三个模块：多细胞自动分割技术模块、细胞交互的多模型 粒子滤波跟踪技术模块、通信与身份管理模块的分析，最终实现多细胞的跟踪；多细胞自动 分割技术针对包含有重叠细胞，粘连细胞的图像序列，利用高帽滤波和低帽滤波对图像进行 预处理后，提出了一种有效的混合检测方法，此方法包含阈值处理、孔洞填充、噪声去除、 图像膨胀和形状与边界限制等几个步骤，从而获得细胞质心位置、面积特征参数；细胞交互 的多模型粒子滤波跟踪技术对多细胞运动特性进行分析，建立合适的细胞运动模型，并估计 细胞运动的角速度；通信与身份管理模块把细胞的运动特征和拓扑特征结合起来，即利用距 离和面积两个特征差异来设计通信法则，并进行细胞身份管理，大大提高了细胞关联的正确 率。

本发明的具体实施方式为：

一、通过多细胞自动分割模块检测分割粘连或重叠的细胞，获取细胞的特征参数如质心 位置、面积等。其实现是通过一种有效的混合检测方法，其具体步骤为：

1）图像增强：利用高帽和低帽滤波变换的组合来增加图像的对比度，其公式如下：

H_k＝A_k-B_k+I_k

其中A_k为高帽滤波图像，B_k为低帽滤波图像，I_k为原始图像，H_k是增强后的图像。

2）细胞混合检测方法：由阈值处理、孔洞填充、噪声去除、图像膨胀和形状与边界限制 五个步骤组成。

K均值聚类的方法对图像没有特别要求，且抗干扰能力较强，故选择K均值聚类的方法 进行阈值处理，其步骤如下：

（1）设计阈值初值；

（2）利用阈值把图像分割为目标和背景两个部分，产生两个子集，像素值大于阈值的 认为是目标（“1”），小于阈值的认为是背景（”0”）；

G₁＝{f(m,n):f(m,n)＞T}(目标像素)

G₂＝{f(m,n):f(m,n)≤T}(背景像素)

（3）计算每个集的平均强度；

m₁＝G₁的平均强度

m₂＝G₂的平均强度

（4）根据平均强度m₁和m₂，重新计算阈值；

T′＝(m₁+m₂)/2

（5）转到第二步直到整个图像被全部分割。

经过阈值处理后的细胞图像存在边界破损、孔洞，故需边界修补和孔洞填充。对 二值图像中进行边界修补时,如果在各个方向（水平、垂直、对角线）上三个像素的和 为2则中间一个像素就由“0”变为“1”，图像中的孔洞可通过标记图像和掩膜图像填 充，再利用中值滤波器滤除高频噪声，最后对图像进行膨胀并且设计宽高比（0.5-5） 限制去除非目标干扰提取细胞的位置和面积等特征信息。

二、通过细胞交互的多模型粒子滤波跟踪模块对存在碰撞、分裂、细胞数目变化等情况 的多细胞进行跟踪。多模型粒子滤波一般由4步构成，混合和交互、粒子滤波、模型概率更 新和状态更新，其步骤如表1所示。

表1 一次循环的多模型粒子滤波过程

本发明提出的目标交互的跟踪方法与传统IMMPF相比，区别在于两个核心问题，一是细 胞角速度ω_k的准确估计，二是确定合适的细胞运动模型。

1）角速度ω_k的估计：在细胞运动过程中ω_k是一个变量，由下面公式进行估计： ω_k,m＝ω_k-1,m+δ_k-1,m，其中ω_k-1,m，是前一时刻角速度，δ_k-1,m是高斯分布协方差为Δ_k-1,m。

2）细胞运动模型：观测细胞运动经常有三种情况发生：一是细胞未与其他细胞发生相 互作用，二是一个细胞在不同帧中与其他细胞发生碰撞，三是碰撞融合的细胞分开或一个细 胞分裂成多个细胞。根据这三种情况定义了三个细胞运动模型分别是增广恒速模型（ACV）、 针对细胞碰撞情况的增广变速率模型（ACT1）和针对细胞分裂情况的增广变速率模型（ACT2）。

$\mathrm{Mode}1\left(\mathrm{ACV}\right):F_k^{\left(1\right)}=\left(\begin{array}{cccccc}1& 0& T& 0& 0& 0\\ 0& 1& 0& T& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right)$

$\mathrm{Mode}2\left(\mathrm{ACT}1\right):F_k^{\left(2\right)}\stackrel{\Delta }{=}\left(\begin{array}{cccccc}1& 0& \frac{\sin {\omega }_{k,1}T}{{\omega }_{k,1}}& \frac{\cos {\omega }_{k,1}T}{{\omega }_{k,1}}& f_{11,k}& 0\\ 0& 1& \frac{1-{\cos \omega }_{k,1}T}{{\omega }_{k,1}}& \frac{{\sin \omega }_{k,1}T}{{\omega }_{k,1}}& f_{21,k}& 0\\ 0& 0& {\cos \omega }_{k,1}T& -{\sin \omega }_{k,1}T& f_{31,k}& 0\\ 0& 0& \sin {\omega }_{k,1}T& \cos {\omega }_{k,1}T& f_{41,k}& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& \alpha \end{array}\right)$

其中$M_{1,k}\stackrel{\Delta }{=}{[f_{11,k},f_{21,k},f_{31,k},f_{41,k}]}^T$用下面公式计算

$f_{11,k}=\frac{{\stackrel{·}{x}}_k({\omega }_{k,1}T\cos {\omega }_{k,1}T-\sin {\omega }_{k,1}T)}{{\omega }_{k,1}^2}-\frac{{\stackrel{·}{y}}_k({\omega }_{k,1}T\sin {\omega }_{k,1}T-1+\cos {\omega }_{k,1}T)}{{\omega }_{k,1}^2}$

$f_{21,k}=\frac{{\stackrel{·}{x}}_k({\omega }_{k,1}T\sin {\omega }_{k,1}T+\cos {\omega }_{k,1}T-1)}{{\omega }_{k,1}^2}-\frac{{\stackrel{·}{y}}_k({\omega }_{k,1}T\cos {\omega }_{k,1}T-\sin {\omega }_{k,1}T)}{{\omega }_{k,1}^2}$

$f_{31,k}=-T({\stackrel{·}{x}}_k\sin {\omega }_{k,1}T+{\stackrel{·}{y}}_k\cos {\omega }_{k,1}T)$

$f_{41,k}=T({\stackrel{·}{x}}_k\cos {\omega }_{k,1}T-{\stackrel{·}{y}}_k\sin {\omega }_{k,1}T)$

Mode3(ACT2):

$F_k^{\left(3\right)}\stackrel{\Delta }{=}\left(\begin{array}{cccccc}1& 0& \frac{\sin {\omega }_{k,2}T}{{\omega }_{k,2}}& \frac{\cos {\omega }_{k,2}T}{{\omega }_{k,2}}& f_{12,k}& 0\\ 0& 1& \frac{1-\cos {\omega }_{k,2}T}{{\omega }_{k,2}}& \frac{\sin {\omega }_{k,2}T}{{\omega }_{k,2}}& f_{22,k}& 0\\ 0& 0& \cos {\omega }_{k,2}T& -\sin {\omega }_{k,2}T& f_{32,k}& 0\\ 0& 0& \sin {\omega }_{k,2}T& \cos {\omega }_{k,2}T& f_{42,k}& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& \beta \end{array}\right)$

其中$M_{2,k}\stackrel{\Delta }{=}{[f_{12,k},f_{22,k},f_{32,k},f_{42,k}]}^T$的计算与$M_{1,k}\stackrel{\Delta }{=}{[f_{11,k},f_{21,k},f_{31,k},f_{41,k}]}^T$上面相同，只是把

ω_k,1用ω_k,2代替。

三、通过细胞通信与身份管理模块同时跟踪和识别多个细胞。通信的作用是利用测量进 行细胞关联，身份管理的作用是识别并对感兴趣的细胞进行标记。下表给出了该模块的 细节。

表2 细胞通信与身份管理模块

采用本发明提供的系统，具有以下特点：

1）成功分割多细胞图像，特别包含重叠或粘连细胞的图像。实验中分割二个图像序列具 有不同的细胞特征的图像共300帧，精确度可达97.03%；

2）对多细胞碰撞、分裂等情形，能实现自动跟踪，并给出所有细胞动力学特征参数（轨 迹、速度、角速度）；

3）能实现对进入、离开图像的多细胞自动跟踪。

附图说明

下面结合附图对本发明的具体实施方式作进一步详细的说明，其中：

图1多细胞自动分割方法框图

图2细胞交互的多模型粒子滤波结构图

(a)单个细胞的跟踪模块（为了简单删除了细胞上标）

(b)模型概率更新

(c)状态更新

图3本系统的结构图

图4细胞跟踪的原始图像序列1

图5细胞跟踪的原始图像序列2

图6细胞分割结果

检测从左到右依次为原始图像、增强图像、我们所提出方法的分割结果和分水岭方法的分割 结果

图7-10为序列1的结果分析，主要跟踪碰撞而且细胞数目变化的情况

图7细胞碰撞、数目在不同帧中变化的多细胞跟踪结果

a)我们提出的方法的跟踪结果

b)每个细胞在x和y方向上的位置估计

图8不同方法对瞬时速度的估计

a)x方向上细胞1的瞬时速度

b)y方向上细胞1的瞬时速度

图9不同方法对所有细胞平均速度的估计

a)x方向上所有细胞的平均速度

b)y方向上所有细胞的平均速度

图10用我们所提出的方法进行角速度估计

a)细胞1

b)细胞4

图11-13为序列2的结果分析，主要跟踪细胞分裂情况

图11序列2细胞分裂的跟踪结果

a）本方法得到的跟踪结果

b)用我们所提出的方法两个细胞在x方向上的位置估计

c)用我们所提出的方法两个细胞在y方向上的位置估计

图12用我们所提出的方法细胞1-1在x和y方向上的瞬时速度

a)x方向

b)y方向

图13用我们所提出的方法对两个细胞角速度的估计

a)细胞1-1

b)细胞1-2

具体实施方式

图3为本系统的结构图。如图3所示，原始图片输入后，通过该系统三个模块：可靠的 多细胞自动分割技术模块、细胞交互的多模型粒子滤波跟踪技术模块、通信与身份管理模块 的分析，最终实现多细胞的跟踪。可靠的多细胞自动分割技术针对包含有重叠细胞，粘连细 胞的图像序列，利用高帽滤波和低帽滤波对图像进行预处理后，提出了一种有效的混合检测 方法，此方法包含阈值处理、孔洞填充、噪声去除、图像膨胀和形状与边界限制等几个步骤， 从而获得细胞质心位置、面积等特征参数。细胞交互的多模型粒子滤波跟踪技术对多细胞运 动特性进行分析，此方法的核心是建立合适的细胞运动模型，并估计细胞运动的角速度。由 于细胞运动的非平滑性，近邻法不能取得满意的关联效果，因此需建立基于细胞特征差异的 通信与身份管理模块，在此模块中把细胞的运动特征和拓扑特征结合起来，即利用距离和面 积两个特征差异来设计通信法则，并进行细胞身份管理，大大提高了细胞关联的正确率。

通过多细胞自动分割模块检测并分割粘连或重叠的细胞，获取细胞的特征参数如质心位 置、面积等。其实现是通过一种有效的混合检测方法，其具体步骤为：

1）图像增强：利用高帽和低帽滤波变换的组合来增加图像的对比度，其公式如下：

H_k＝A_k-B_k+I_k

其中A_k为高帽滤波图像，B_k为低帽滤波图像，I_k为原始图像，H_k是增强后的图像。

2）细胞混合检测方法：由阈值处理、孔洞填充、噪声去除、图像膨胀和形状与边界限制五 个步骤组成。

用我们所提出的方法检测两个序列300帧图像，精确率可达97.03%。

从图6看出我们所提出的方法分割效果良好，而分水岭方法分割存在过分割现象。

通过细胞交互的多模型粒子滤波跟踪模块对存在碰撞、分裂、细胞数目变化等情况的多细胞 进行跟踪。多模型粒子滤波一般由4步构成，混合和交互、粒子滤波、模型概率更新和状态 更新。此方法成功应用需解决两个核心问题，一是细胞角速度ω_k的准确估计，二是确定合适 的细胞运动模型。

1）角速度ω_k的估计：在细胞运动过程中ω_k是一个变量，由下面公式进行估计： ω_k,m＝ω_k-1,m+δ_k-1,m，其中ω_k-1,m，是前一时刻角速度，δ_k-1,m是高斯分布协方差为Δ_k-1,m。

2）细胞运动模型：观测细胞运动经常有三种情况发生：一是细胞未与其他细胞发生相 互作用，二是一个细胞在不同帧中与其他细胞发生碰撞，三是融合的细胞分开或一个细胞分 裂成多个细胞。根据这三种情况定义了三个细胞运动模型分别是增广恒速模型（ACV）、针对 细胞碰撞情况的增广变速率模型（ACT1）和针对细胞分裂情况的增广变速率模型（ACT2）。

$\mathrm{Mode}1\left(\mathrm{ACV}\right):F_k^{\left(1\right)}=\left(\begin{array}{cccccc}1& 0& T& 0& 0& 0\\ 0& 1& 0& T& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\end{array}\right)$

$\mathrm{Mode}2\left(\mathrm{ACT}1\right):F_k^{\left(2\right)}\stackrel{\Delta }{=}\left(\begin{array}{cccccc}1& 0& \frac{\sin {\omega }_{k,1}T}{{\omega }_{k,1}}& \frac{\cos {\omega }_{k,1}T}{{\omega }_{k,1}}& f_{11,k}& 0\\ 0& 1& \frac{1-\cos {\omega }_{k,1}T}{{\omega }_{k,1}}& \frac{\sin {\omega }_{k,1}T}{{\omega }_{k,1}}& f_{21,k}& 0\\ 0& 0& \cos {\omega }_{k,1}T& -\sin {\omega }_{k,1}T& f_{31,k}& 0\\ 0& 0& \sin {\omega }_{k,1}T& \cos {\omega }_{k,1}T& f_{41,k}& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& \alpha \end{array}\right)$

其中$M_{1,k}\stackrel{\Delta }{=}{[f_{11,k},f_{21,k},f_{31,k},f_{41,k}]}^T$用下面公式计算

$f_{11,k}=\frac{{\stackrel{·}{x}}_k({\omega }_{k,1}T\cos {\omega }_{k,1}T-\sin {\omega }_{k,1}T)}{{\omega }_{k,1}^2}-\frac{{\stackrel{·}{y}}_k({\omega }_{k,1}T\sin {\omega }_{k,1}T-1+\cos {\omega }_{k,1}T)}{{\omega }_{k,1}^2}$

$f_{21,k}=\frac{{\stackrel{·}{x}}_k({\omega }_{k,1}T\sin {\omega }_{k,1}T+\cos {\omega }_{k,1}T-1)}{{\omega }_{k,1}^2}-\frac{{\stackrel{·}{y}}_k({\omega }_{k,1}T\cos {\omega }_{k,1}T-\sin {\omega }_{k,1}T)}{{\omega }_{k,1}^2}$

$f_{31,k}=-T({\stackrel{·}{x}}_k\sin {\omega }_{k,1}T+{\stackrel{·}{y}}_k\cos {\omega }_{k,1}T)$

$f_{41,k}=T({\stackrel{·}{x}}_k\cos {\omega }_{k,1}T-{\stackrel{·}{y}}_k\sin {\omega }_{k,1}T)$

Mode3(ACT2):

$F_k^{\left(3\right)}\stackrel{\Delta }{=}\left(\begin{array}{cccccc}1& 0& \frac{\sin {\omega }_{k,2}T}{{\omega }_{k,2}}& \frac{\cos {\omega }_{k,2}T}{{\omega }_{k,2}}& f_{12,k}& 0\\ 0& 1& \frac{1-\cos {\omega }_{k,2}T}{{\omega }_{k,2}}& \frac{\sin {\omega }_{k,2}T}{{\omega }_{k,2}}& f_{22,k}& 0\\ 0& 0& \cos {\omega }_{k,2}T& -\sin {\omega }_{k,2}T& f_{32,k}& 0\\ 0& 0& \sin {\omega }_{k,2}T& \cos {\omega }_{k,2}T& f_{42,k}& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& \beta \end{array}\right)$

其中$M_{2,k}\stackrel{\Delta }{=}{[f_{12,k},f_{22,k},f_{32,k},f_{42,k}]}^T$的计算与$M_{1,k}\stackrel{\Delta }{=}{[f_{11,k},f_{21,k},f_{31,k},f_{41,k}]}^T$上面相同，只是把 ω_k,1用ω_k,2代替。

图3给出了本发明提出的细胞交互的多模型粒子滤波跟踪算法的结构图。

通过细胞通信与身份管理模块同时跟踪和识别多个细胞。通信的作用是利用测量进行细 胞关联，身份管理的作用是识别并对感兴趣的细胞进行标记。下表给出了该模块的细节。

表2 细胞通信与身份管理模块

经过细胞交互的多模型粒子滤波跟踪技术模块、通信与身份管理模块，细胞跟踪结果如 下面例子（例1，例2）所示。

图7-10为序列1（例1）的序列跟踪结果分析，序列1的结果分析，主要跟踪碰撞而且 细胞数目变化的情况。由图7看出细胞2和细胞4在36帧碰撞38帧分离，细胞2和细胞4 在40帧又一次碰撞42帧分离，细胞6向右移动在36帧已经部分离开视觉范围，38帧中完 全消失，细胞7和细胞8在46帧中进入视觉范围，用我们提出的方法所有的细胞都被可靠跟 踪到。图8给出了细胞1在x和y方向上的瞬时速度，图9给出了所有细胞的平均速度，图 10给出了细胞1和细胞4的角速度。

图11-13为序列2（例2）的结果分析，主要跟踪细胞分裂情况。

由图11看出细胞从69帧开始分裂，70帧完成分裂，我们所提出的方法跟踪到了这一结果。图 12给出了用我们所提出的方法细胞1-1在x和y方向上的瞬时速度，图13给出了用我们所提出 的方法对两个细胞转弯角速度的估计。