CN102034246B - 基于拓扑约束和匈牙利算法的高密度细胞追踪方法 - Google Patents

Abstract

本发明提供的是一种基于拓扑约束和匈牙利算法的高密度细胞追踪方法。(1)采用水平集算法与局部灰度阈值法相结合的图像分割方法对细胞图像序列进行分割，并对每一帧分割后的细胞初始标号；(2)对第k帧中的任意待匹配细胞根据距离限制在k+1帧中建立追踪的搜索区域，将区域中的细胞列为候选细胞；(3)建立系数矩阵Q。若k+1帧中细胞j是第k帧中细胞i的候选细胞，则根据拓扑约束进行数据关联计算其相似度Q_ij，否则将其相似度赋一较大值。(4)利用匈牙利算法对系数矩阵进行变换，找出独立零元素，其行列代表的细胞相匹配；(5)找出矩阵变换后没有独立零元素的行列，对其所对应的细胞分别进行考虑；(6)将k加1，跳到第2步，重复进行，直到图像序列的最后一帧。本发明能实现较高效率的细胞追踪。

Description

基于拓扑约束和匈牙利算法的高密度细胞追踪方法

技术领域

本发明涉及的是一种细胞运动分析方法，具体地说是一种采用拓扑约束和匈牙利算法相结合的方法实现高密度细胞追踪。

背景技术

细胞运动分析通常是对目标细胞进行追踪，得到细胞的运动速度、位移、轨迹等信息，在细胞行为、药物和疾病的研究中是至关重要的。我们对一个目标细胞进行详细和准确的追踪，特别是对大量细胞进行追踪时，要尽量避免其他细胞的干扰，这是细胞追踪的难点，也是近些年细胞追踪的研究热点。

本发明提到的追踪方法是基于拓扑约束和匈牙利算法相结合的追踪方法，是一种基于分割的追踪方法。基于图论的细胞追踪方法是张乐林等人2007年在Proceedings of the 7thIEEE International Conference on BIBE发表的论文“Graph Theory Application in Cellnuleus Segmentation，Tracking and Identification”中提出的细胞追踪方法。该方法利用图模型描述细胞之间的拓扑关系，将分割图像中追踪细胞特征问题转变为相似结构图之间顶点匹配问题。

图论的方法已经证实能够有效的追踪高密度细胞图像，且对于细胞的形变不敏感。但该方法中几个关键参数的人工设定限制了其可适性，并且对于含有稀疏细胞的图像或图像中的部分区域处理能力不强，尤其是当相邻细胞发生分裂时易造成细胞的错误追踪。

发明内容

本发明的目的在于提供一种能够实现较高效率的高密度细胞追踪方法：基于拓扑约束和匈牙利算法的高密度细胞追踪方法。

本发明的具体实现步骤如下：

1、一种基于拓扑约束和匈牙利算法的高密度细胞追踪方法，其特征是：

(1)、细胞分割：采用水平集算法与局部灰度阈值法相结合的序列细胞图像分割方法进行分割，然后将图像每一帧分割后的细胞初始标号，将第k帧中标号为i的细胞记录为R(k，i)；

(2)、候选细胞的选取：对第k帧中的任意待匹配细胞R(k，i)根据距离限制在k+1帧中建立追踪的搜索区域，将区域中的细胞列为候选细胞；

(3)、系数矩阵的构建：

①若第k+1帧中细胞R(k+1，j)是第k帧待匹配细胞R(k，i)的候选细胞，则根据下式进行数据关联计算其相似程度Q_ij，

$Q_{\mathrm{ij}}=Q_{\mathrm{Dist}}(R(k,i),R(k+1,j))\·Q_{\mathrm{SSN}}(R(k,i),R(k+1,j))\·Q_{\mathrm{area}}(R(k,i),R(k+1,j))$

$\left\{\begin{array}{c}{Q}_{\mathrm{Dist}}(R(k,i),R(k+1,j))=\exp \left(\frac{\mathrm{Dist}(R_{\mathrm{center}}(k,i),R_{\mathrm{center}}(k+1,j))}{{\mathrm{\δ}}_D}\right)\\ Q_{\mathrm{area}}(R(k,i),R(k+1,j))=\exp \left(\frac{\mathrm{abs}(R_{\mathrm{area}}(k,i)-R_{\mathrm{area}}(k+1,j))}{{\mathrm{\δ}}_S}\right)\\ Q_{\mathrm{SSN}}(R(k,i),R(k+1,j))=\underset{k=1}{\overset{M}{\mathrm{\Σ}}}\exp \left(\frac{\mathrm{Angle}(D_{R(k,i)}\left({\mathrm{\α}}_k\right),D_{R(k+1,j)}\left({\mathrm{\β}}_k\right))}{{\mathrm{\δ}}_{\mathrm{\θ}}}\right)\end{array}\right.$

其中Q_Dist、Q_SSN、Q_area分别表示的是候选细胞与待匹配细胞之间的空间相似性、邻域相似性、细胞大小的相似性。R_center(·)、R_area(·)分别表示的是细胞顶点的中心点坐标和细胞的面积，Dist(·)表示的是两中心坐标之间的距离，Angle(·)是指候选细胞和待匹配细胞各自邻域细胞所构成的邻域角度分布集合中相对应的角度值的差别，这里邻域角度分布集合中所有元素以角度值减小的方向排列。常数δ_D、δ_S和δ_θ分别用于调整Q_ij对于距离、面积和角度的灵敏度。

②若不是候选细胞，将其相似程度Q_ij赋一较大的值M。

则由Q_ij组成的矩阵Q＝(Q_ij)_m*n就是所要构建的系数矩阵，其中m，n分别为第k帧和第k+1帧所分离出的细胞数；

(4)、细胞全局匹配：利用匈牙利算法对目标矩阵进行变换，使得变换后的矩阵中出现独立的零元素，该零元素行列所代表的细胞相匹配，则将行所代表的细胞标号赋值给所代表的细胞；

(5)、未匹配细胞分析：找出矩阵变换后没有独立零元素的行列，对其所对应的细胞分别进行考虑；

(6)、将k加1，跳到第2步，重复进行，直到图像序列的最后一帧。

2、未匹配细胞的具体分析方法为：

①丢失的细胞：找出没有独立0元素的行i，表示该行所代表的细胞R(k，i)在第k+1帧中没有找到与之相匹配的细胞。若该细胞不是位于图像边缘的细胞，首先找出该细胞的候选细胞，如果所有的候选细胞都已经与别的待匹配目标细胞相匹配，则考虑这些已匹配细胞对之间的关系，根据下面三个判断条件来判断是否是由于欠分割导致该细胞的丢失。如果满足这些条件则认为是由于该细胞与其邻域细胞在下一帧中欠分割所造成的丢失：

R(k+1，j).area＞1.7*R(k，l).area

abs(fw_angle(R(k，i)，R(k，l))-R(k+1，j).Orientation)＜10

R(k+1，j)MajorAxisLength＞1.7*R(k，l)MajorAxisLength

其中细胞R(k+1，j)是细胞R(k，i)的候选细胞，细胞R(k+1，j)与细胞R(k，1)相匹配；fw_angle(R(k，i)，R(k，1))表示的是细胞R(k，i)和R(k，1)的质心连线与x轴的夹角；R(k+1，j).Orientation表示的是细胞R(k+1，j)的离心率；R(k，j).MajorAxisLength表示的是细胞R(k+1，j)的长轴；

若由于细胞R(k，i)与R(k，1)在k+1帧中的欠分割导致细胞R(k，i)的丢失，则接着判断细胞R(k，i)与R(k，1)是否是惰性细胞，若是惰性细胞，通过复制的方法将R(k，i)与R(k，1)两细胞替代细胞R(k+1，j)。否则可能是由于该细胞图像边缘与背景的对比度的非常弱，以致没有从背景中识别分割出来，从而导致该细胞的丢失。若细胞R(k，i)是位于图像边界的细胞，则认为细胞移出边界，将该消失细胞的信息存储，以防细胞再一次出现；

②新出现的细胞：找出没有独立0元素的列j，则表示该列所代表的细胞R(k+1，j)在第k帧中没有找到与之相匹配的细胞，为新出现的细胞，然后判断其新出现的原因。

首先根据细胞生物分裂前后的特性判断该细胞是否是细胞分裂成的子细胞。接着查看该新出现的细胞是否是位于边缘的细胞，若是，则认为细胞移入图像边界，然后将该细胞信息与前几帧丢失的细胞信息进行对比。这里主要是对质心位置信息进行对比，若是前几帧消失的细胞，则将消失细胞的ID号给该新出现的细胞，否则认为是新出现的细胞，给其赋予一个新的ID号并存储。

附图说明

图1是细胞图像序列第一帧的细胞图像原图；

图2是细胞图像序列第一帧的细胞分割图像；

图3是细胞的邻域分布图；

图4是本发明的具体实现流程图；

图5是利用拓扑约束方法对50帧细胞图像进行追踪的追踪轨迹图；

图6是利用本发明方法对50帧细胞图像进行追踪的追踪轨迹图；

图7是利用本发明方法对100帧细胞图像进行追踪的追踪轨迹图。

具体实施方式

本发明采用拓扑约束和匈牙利算法相结合进行细胞追踪，其具体实现过程步骤如下：

1、细胞分割：由于本文采用的细胞追踪方法是一种基于分割的追踪方法，所以首先要对细胞图像序列进行分割。本文所追踪的细胞序列，是由光学显微镜成像的未经染色的神经干细胞图像序列。该图像序列的特点是目标与背景的对比度较弱，多数图像帧中存在着细胞粘连、团簇。针对该序列的特点，本文采用了水平集算法与局部灰度阈值法相结合的序列细胞图像分割方法，细胞图像序列第一帧的细胞图像原图及分割图分别如图1、图2所示。然后将图像每一帧分割后的细胞初始标号。为了方便描述我们将第k帧中标号为i的细胞记录为R(k，i)。

2、候选细胞的选取：对第k帧中的任意细胞R(k，i)来说，先根据距离限制在k+1帧中建立追踪的搜索区域，将区域中的细胞列为候选细胞，本发明中经试验该距离设置为细胞平均直径的2倍。

3、目标矩阵的构建：对细胞R(k，i)来说，若k+1帧中细胞R(k+1，j)是它的候选细胞，则根据(1)式进行数据关联计算其相似程度Q_ij。若不是候选细胞，将其相似程度Q_ij赋一较大的值M，在本发明试验中M设置为100。则由Q_ij组成的矩阵Q＝(Q_ij)_m*n就是所要构建的目标矩阵，其中m，n分别为第k帧和第k+1帧所分离出的细胞数。

4、细胞全局匹配：利用匈牙利算法对目标矩阵进行变换，使得变换后的矩阵中出现独立的零元素，该零元素行列所代表的细胞相匹配，则将行所代表的细胞标号赋值给所代表的细胞。

5、未匹配细胞分析：找出矩阵变换后没有独立零元素的行列，对其所对应的细胞分别进行考虑。

①首先找出没有独立0元素的行i，表示该行所代表的细胞R(k，i)在第k+1帧中没有找到与之相匹配的细胞。若细胞R(k，i)不是位于图像边缘的细胞，首先找出该细胞的候选细胞，如果所有的候选细胞都已经与别的目标细胞相匹配，则根据进一步计算这些已匹配细胞对之间的关系，来判断是否是由于欠分割导致细胞R(k，i)的丢失。在本实验仿真中，如果满足下列条件则认为是由于欠分割造成的细胞R(k，i)的丢失。

R(k+1，j).area＞1.7*R(k，l).area

abs(fw_angle(R(k，i)，R(k，l))-R(k+1，j).Orientation)＜10

R(k+1，j)MajorAxisLength＞1.7*R(k，l)MajorAxisLength

其中细胞R(k+1，j)是细胞R(k，i)的候选细胞，细胞R(k+1，j)与细胞R(k，1)相匹配。fw_angle(R(k，i)，R(k，1))表示的是细胞R(k，i)和R(k，1)的质心连线与x轴的夹角。R(k+1，j).Orientation表示的是细胞R(k+1，j)的离心率。R(k，j).MajorAxisLength表示的是细胞R(k+1，j)的长轴。

若由于R(k，i)与R(k，1)在k+1帧中的欠分割导致细胞R(k，i)的丢失，则判断细胞R(k，i)与R(k，1)是否是惰性细胞，若是惰性细胞，我们可以通过复制的方法将R(k，i)与R(k，1)两细胞替代细胞R(k+1，j)。如图5所示本发明实验图像序列第20帧和第21帧的追踪结果。

若判断出R(k，i)不是和其领域细胞在下一帧中由于欠分割造成的丢失，则可能是由于该细胞图像边缘与背景的对比度的非常弱，以致没有从背景中识别分割出来，从而导致该细胞的丢失。若细胞R(k，i)是位于边缘的细胞，则认为细胞移出图像边缘。我们将该消失细胞的信息存储，以防细胞再一次出现。

②然后找出没有独立0元素的列j，则表示该行所代表的细胞R(k+1，j)在第k帧中没有找到与之相匹配的细胞，为新出现的细胞，然后判断其新出现的原因。在细胞的整个追踪过程中，搜索细胞分裂并存贮细胞分裂前后的有关信息是一个非常重要的方面。我们利用细胞分裂前后的特征来分析在追踪过程中新出现的细胞，是否是由分裂产生的。但是，也有可能由于过分割将一个细胞错误地分割成两个细胞，并且也满足细胞分裂前后的特征，这样我们就有可能错误地认为是由于细胞分裂造成的新细胞的出现。首先根据细胞分裂前后的特征找出“母细胞”R(k，i)，方法与上面介绍欠分割时的方法刚好相反。在本实验仿真中，判断后几帧中“母细胞”R(k，i)区域的细胞特征，若只有一个细胞，并且该细胞面积，长轴和“母细胞”R(k，i)的面积，长轴基本上相差不大，则认为是细胞R(k，i)在k+1帧中过分割造成R(k+1，j)的出现。否则认为是细胞R(k，i)在k+1帧中发生了细胞分裂，R(k+1，j)是其子细胞，对此将细胞分裂信息进行存储。

若新出现的细胞R(k+1，j)是位于边缘的细胞，则认为细胞移入图像边缘，首先将该细胞信息与前几帧丢失的细胞信息进行对比，主要是对质心位置进行比较，以此来判断是否是前几帧消失的细胞。若是前几帧消失的细胞，则将消失细胞的ID号给该新出现的细胞，否则认为是新出现的细胞，给其赋予一个新的ID号并存储。

6、将k加1，跳到第2步，重复进行，直到图像序列的最后一帧。

实验结果表明：基于拓扑约束和匈牙利算法相结合进行细胞追踪的效果要好于单独基于拓扑约束算法进行细胞追踪的追踪效果。首先分别用这两种方法对50帧细胞图像进行追踪，其追踪轨迹图分别如图5、图6所示，以此验证本发明在细胞追踪方面的有效性和可行性；然后利用本发明对100帧该细胞图像序列进行追踪，其追踪轨迹图如图7所示，并对实验结果数据进行统计，如表1所示。

表1拓扑约束法与本文方法的实验结果对比

由表1可以看出，基于拓扑约束和匈牙利算法相结合进行细胞追踪的最终准确率比基于拓扑约束算法进行细胞追踪的最终准确率提高了5.27％。实验结果表明，本发明进行细胞追踪有更好的追踪效果，准确性更高。

Claims (1)

1.一种基于拓扑约束和匈牙利算法的高密度细胞追踪方法，其特征是：

(1)、细胞分割：采用水平集算法与局部灰度阈值法相结合的序列细胞图像分割方法进行分割，然后将图像每一帧分割后的细胞初始标号，将第k帧中标号为i的细胞记录为R(k，i)；

(2)、候选细胞的选取：对第k帧中的任意待匹配细胞R(k，i)根据距离限制在k+1帧中建立追踪的搜索区域，将区域中的细胞列为候选细胞；

(3)、系数矩阵的构建：

①若第k+1帧中细胞R(k+1，j)是第k帧待匹配细胞R(k，i)的候选细胞，则根据下式进行数据关联计算其相似程度Q_ij，

其中Q_Dist、Q_SSN、Q_area分别表示的是候选细胞与待匹配细胞之间的空间相似性、邻域相似性、细胞大小的相似性；R_center(·)、R_area(·)分别表示的是细胞顶点的中心点坐标和细胞的面积，Dist(·)表示的是两中心坐标之间的距离，Angle(·)是指候选细胞和待匹配细胞各自邻域细胞所构成的邻域角度分布集合中相对应的角度值的差别，邻域角度分布集合中所有元素以角度值减小的方向排列，常数δ_D、δ_S和δ_θ分别用于调整Q_ij对于距离、面积和角度的灵敏度；

②若第k+1帧中细胞R(k+1，j)不是第k帧待匹配细胞R(k，i)的候选细胞，将其相似程度Q_ij赋一较大的值M；

则由Q_ij组成的矩阵Q＝(Q_ij)_m*n就是所要构建的系数矩阵，其中m，n分别为第k帧和第k+1帧所分离出的细胞数；

(4)、细胞全局匹配：利用匈牙利算法对系数矩阵进行变换，使得变换后的矩阵中出现独立的零元素，该零元素行列所代表的细胞相匹配，则将行所代表的细胞标号赋值给所代表的细胞；

(5)、未匹配细胞分析：找出矩阵变换后没有独立零元素的行列，对其所对应的细胞分别进行考虑；

(6)、将k加1，跳到第2步，重复进行，直到图像序列的最后一帧；

所述的未匹配细胞的具体分析方法为：

①丢失的细胞：找出没有独立零元素的行i，表示该行所代表的细胞R(k，i)在第k+1帧中没有找到与之相匹配的细胞，若该细胞不是位于图像边缘的细胞，首先找出该细胞的候选细胞，如果所有的候选细胞都已经与别的待匹配细胞相匹配，则考虑这些已匹配细胞对之间的关系，根据下面三个判断条件来判断是否是由于欠分割导致该细胞的丢失，如果满足这些条件则认为是由于该细胞与其邻域细胞在下一帧中欠分割所造成的丢失：

R(k+1，j).area＞1.7*R(k，l).area

abs(fw_angle(R(k，i)，R(k，l))-R(k+1，j).Orientation)＜10

R(k+1，j).MajorAxisLength＞1.7*R(k，l).MajorAxisLength

其中细胞R(k+1，j)是细胞R(k，i)的候选细胞，细胞R(k+1，j)与细胞R(k，l)相匹配；fw_angle(R(k，i)，R(k，l))表示的是细胞R(k，i)和R(k，l)的质心连线与x轴的夹角；R(k+1，j).Orientation表示的是细胞R(k+1，j)的离心率；R(k+1，j).MajorAxisLength表示的是细胞R(k+1，j)的长轴；

若由于细胞R(k，i)与R(k，l)在k+1帧中的欠分割导致细胞R(k，i)的丢失，则接着判断细胞R(k，i)与R(k，l)是否是惰性细胞，若是惰性细胞，通过复制的方法将R(k，i)与R(k，l)两细胞替代细胞R(k+1，j)，否则是由于该细胞图像边缘与背景的对比度非常弱，以致没有从背景中识别分割出来，从而导致该细胞的丢失，若细胞R(k，i)是位于图像边界的细胞，则认为细胞移出边界，将该丢失的细胞的信息存储，以防细胞再一次出现；

②新出现的细胞：找出没有独立0元素的列j，则表示该列所代表的细胞R(k+1，j)在第k帧中没有找到与之相匹配的细胞，为新出现的细胞，然后判断其新出现的原因；

首先根据细胞生物分裂前后的特性判断该细胞是否是细胞分裂成的子细胞，接着查看该新出现的细胞是否是位于边缘的细胞，若是，则认为细胞移入图像边界，然后将该细胞信息与前几帧丢失的细胞信息进行对比，这里主要是对质心位置信息进行对比，若是前几帧丢失的细胞，则将丢失的细胞的ID号给该新出现的细胞，否则认为是新出现的细胞，给其赋予一个新的ID号并存储。 

Images