基于局部图动态匹配的植物细胞追踪算法研究

钱伟丽 刘 敏 李洁沁 刘小燕

(湖南大学电气与信息工程学院， 长沙 410082)

引言

图1 植物顶端分生组织的显微图像栈序列。 (a)植物样本；(b)植物顶端分生组织的放大特写；(c)共聚焦显微镜成像系统；(d)通过共聚焦激光扫描显微镜得到的顶端分生组织细胞沿时间方向的图像栈序列(图像可能发生旋转或平移)；(e)基于分水岭算法的分割结果Fig.1 The time series of microscopic plant cell image stack。 (a) A plant sample; (b) Shoot apical meristem; (c) Confocal laser scanning microscope; (d) An example of the cell images taken by confocal laser scanning microscopy at different time instances (images can be translated or rotated); (e) The segmentation results of (d) based on watershed algorithm

生物体的生长发育模式主要包括细胞群体的空间排列规则以及细胞的迁移、扩张和分裂等动态行为规律，因此对细胞群体的行为规律做出精确的理解是建立生物体生长发育模型的关键步骤[1-2]。基于光学显微成像技术，研究活体内的细胞运动模型以及细胞生长规律，进而探索生物体基因的结构功能以及生物体发育过程的科学机理，对于生命科学基础研究具有重大意义。为此，对显微图像数据中细胞群的追踪算法研究具有非常重要的理论和应用价值。近几年来，国内外学者针对细胞追踪问题做了大量富有成效的工作，如基于半局部邻域贝叶斯判别的细胞追踪算法[3]、水平集分割算法[4-5]、拓扑约束法[6-8]、蜂窝划分和反馈纠正的追踪算法[9]、多粒子滤波器和多维分配的追踪算法[10]、概率松弛标记算法[11]、维特比算法[12]、基于细胞纹理特征的追踪算法[13]、基于复合mean-shift算法的细胞追踪方法[14]等。但是，由于植物顶端分生组织图像中细胞之间相互紧挨且形状、灰度值相近这一特殊性(见图1)，这些方法都不能直接用于该图像中细胞的追踪。法国国立计算机及自动化研究院提出了关于植物细胞的多角度图像采集以及图像分析方法[15]，但他们的图像采集方法跟本项目采用的共聚焦激光扫描显微镜的活细胞成像技术[16]有很大不同：他们所采用的激光能量密度比较高，长时间照射会导致细胞死亡，所以不能实现较密集时间间隔的图像数据采集。针对所研究的植物细胞图像，有研究者提出基于条件随机场的细胞追踪算法[17]，其追踪植物细胞的前提是图像经过配准。

早前的研究提出了一种基于局部图匹配的细胞追踪算法[18-19]，该算法对图像中细胞提取面积、相邻细胞距离、方位信息和质心位置特征，然后比较相邻两帧图像中细胞的特征，选择最相似的细胞对作为种子细胞对，再从种子细胞向外扩散，匹配种子细胞的邻细胞。在每次匹配过程中，种子细胞不断更新为新匹配的细胞。该追踪算法在噪声小的图像上能取得较好的追踪结果，但是对于噪声严重的图像在匹配过程中容易出现匹配错误的细胞对，而该对细胞作为新的种子细胞对会造成其邻细胞匹配错误，随着错误的不断累积，将导致最后追踪准确度低。同时，局部图匹配算法通过质心位置与细胞方位信息特征进行细胞匹配，这些特征的使用前提是图像经过配准。但细胞成像过程中由于细胞平面被错位与被旋转带来的噪声，给细胞图像的配准带来很大的难度，所以开发能够对抗旋转与平移噪声的匹配算法至关重要。

因此，笔者提出的局部图动态匹配算法，在以下3个方面改进了之前的算法[18]：

1) 图像出现随机旋转与位移噪声时，细胞与其周围细胞组成的局部图的几何形状以及拓扑结构不变，提取不变特征(即细胞的面积、相邻细胞的夹角与距离)来进行细胞匹配，使本研究提出的算法能对抗成像过程中带入的位移与旋转噪声。

2) 已匹配的细胞作为新增加的种子细胞，在动态扩张的已匹配细胞邻域范围中，每次优先匹配特征距离最小的细胞对，通过动态匹配模式，提高细胞匹配的准确率。

3) 结合细胞分裂与细胞局部图的特点，对分裂的细胞进行检测。

图2给出了本算法的流程。先找到最相似的细胞对作为种子细胞对，然后从种子细胞开始动态匹配种子细胞的邻细胞，最后检测分裂的细胞。

图2 局部图动态匹配算法流程Fig.2 The diagram of the dynamic local graph matching method

1 方法

1.1 细胞图像的分水岭分割

在细胞图像中，每张细胞图像上有近百个细胞，细胞之间结构紧密且形状相似(见图1(d))。因此，对多细胞进行自动追踪时，首先必须将各个细胞从细胞图像中准确地分割出来，明确各个细胞边界，才能精确定位细胞的位置信息以提取相应的其他特征。本研究采用分水岭分割算法对图像进行分割[20]，将图像像素灰度值从小到大进行排序，然后从小到大实现浸没。在分割结束后，图像中的每一个极小值将对应目标图像中的一个独立区域，形成的目标区域数目由图像中局部极小值的数目决定。在早期的研究中，通过引入追踪算法的反馈来实现分割和追踪的共同最优[18]，减少分割错误，目前该算法在细胞图像分割上已经得到较好的结果。分水岭算法对细胞图像分割的结果如图1(e)所示，可以看出细胞之间边界明确，表明分水岭算法对植物细胞图像分割效果较好。

1.2 局部图动态匹配算法

1.2.1细胞之间局部图

由于植物顶端分生组织细胞图像中细胞与细胞之间结构紧密，且形状相似，如果仅仅依靠细胞的面积等特征，很难将所有的细胞加以区分，追踪时也会出现错误。因此，本研究首先用图来描述细胞图像，然后提取细胞之间的局部拓扑结构作为细胞追踪时的一个特征。在图论中，图是由顶点和连接相邻顶点的边构成的。本研究用G表示细胞与其邻细胞之间的局部图，图中细胞用顶点表示，相邻细胞之间通过边连接，中心细胞i与其邻细胞i1～i6通过边连接，组成细胞i的局部图Gt(见图3(a))。

图3 相邻两帧细胞图像中，选择种子细胞时的细胞局部Gt与Gt+1(存在旋转与平移)Fig.3 Local graphs Gt and Gt+1 for seed cell pair selection in two adjacent images (with translation and rotation)

1.2.2选取种子细胞对

由于每幅细胞图像有近百个细胞，且细胞之间结构紧密，形状相似，若将相邻两帧中细胞依次比较，不仅计算量大且极易出现追踪错误。考虑到细胞图像采集的时间间隔较短，连续两帧细胞图像中细胞之间相互位置不会发生剧变，因此可以通过先匹配最相似的细胞对作为种子细胞对，然后从种子细胞出发匹配其邻细胞，缩小搜索范围，提高细胞追踪准确率。

当细胞图像中存在旋转与错位时,追踪时必须利用图像出现旋转与平移时的不变特征。考虑到图像发生旋转与平移时，细胞几何形状、细胞之间的拓扑关系不会发生变化，也就是细胞面积、细胞之间的位置关系和细胞之间的距离不会变化，因此可以提取细胞的面积A、细胞与相邻细胞之间的夹角θ(见图3(a)中细胞i和细胞i5、i6组成的夹角θ1)、细胞与相邻细胞连线的边长L(见图3(a)中细胞i和细胞i4的边长L1)这些特征进行比较，然后结合上述的细胞局部图进行匹配。因此，种子细胞的选取就相当于将连续两帧图像中每个细胞的局部图进行比较(见图3(a)中局部图Gt与(b)中局部图Gt+1)，然后选取最相似的局部图中的中心细胞为种子细胞。给定连续两帧细胞图像It与It+1，种子细胞对的选取步骤如下：

1) 判断细胞的邻细胞个数是否相等。假设图像It中细胞i的邻细胞为{i1,i2，…，im}∈Ni，有m个邻细胞，图像It+1中细胞j的邻细胞为{j1,j2，…，jn}∈Nj，有n个邻细胞。如果m=n，细胞i与j被认定为候选种子细胞对；若m≠n，则认为细胞i与j不是种子细胞对。

2) 通过邻细胞相等条件筛选出候选细胞对后，构建细胞自身面积、邻细胞到中心细胞的距离、相邻两邻细胞与中心细胞的夹角、邻细胞面积的距离函数，分别为f1～f4，有

(1)

式中，Ai与Aj分别表示细胞i与j的面积，Lik与Ljk分别表示细胞i与j和各自的邻细胞ik(k=1，…，m)与jk(k=1，…，m)的距离，θik与θjk分别表示细胞i与j和各自相邻的两个邻细胞的夹角(见图3(a)中细胞i与细胞i5、i6组成的角θ1)，Aik与Ajk分别表示细胞i与j各自的邻细胞ik(k=1，…，m)与jk(k=1，…，m)的面积。

3) 不同距离函数通过参数相加，得总的距离函数，有

D(i,j)=λ1f1(i,j)+λ2f2(i,j)+λ3f3(i,j)+λ4f4(i,j)

(2)

式中，λ1、λ2、λ3、λ4通过本文第1.2.5节部分求解。

假设候选种子细胞对为M，选择距离函数最小的细胞对(u,v)为种子细胞对，有

(i,j=1，…，M)

(3)

1.2.3邻细胞匹配

通过以上步骤找到图像It与It+1中初始种子细胞对(u,v)后，开始匹配种子细胞的邻细胞。由于细胞之间相互紧挨，也就是说种子细胞的邻细胞的位置信息基本确定，因此在进行邻细胞追踪时，可以适当减弱匹配条件，细胞局部图由所有邻细胞组成改为由两个邻细胞且其中一个是种子细胞所组成的三角形局部图结构，如图4中细胞局部图Gt与Gt+1所示。

图4 相邻两帧细胞图像中，邻细胞匹配时细胞局部图Gt与Gt+1(存在旋转与平移)Fig.4 Local graphs Gt and Gt+1 for the neighboring cells matching in two adjacent images (with translation and rotation)

选定种子细胞对(u,v)后，开始匹配邻细胞。种子细胞对(u,v)的邻域分别为{u1,u2,…，um}∈Nu和{v1,v2,…，vm}∈Nv，从{u1,u2,…，um}与{v1,v2,…，vm}中选取一对最相似的邻细胞对。类似找种子细胞对，分别将图像It中细胞u的邻细胞up(p=1，…，m)与图像It+1中细胞v的邻细胞vq(q=1，…，m)的面积、夹角及边长比较，距离函数分别为f1、f2和f3，有

(4)

式中，Aup与Avq分别表示邻细胞up(p=1，…，m)与vq(q=1，…，m)的面积，Lup与Lvq分别表示邻细胞up(p=1，…，m)与vq(q=1，…，m)和各自种子细胞u与v的距离，θup与θvq分别表示邻细胞up(p=1，…，m)与vq(q=1，…，m)和各自相邻细胞及种子细胞的夹角(见图4(a)中细胞u与细胞u1、u2组成的角θ1)。

在细胞追踪时，细胞之间差异越小则细胞是同一细胞的可能性越高。因此，在邻细胞追踪时，选择距离函数最小的邻细胞对作为追踪结果。邻细胞对总的距离函数如下：

(5)

式中，ω1～ω3通过本文第1.2.5节部分求解。

选择距离函数最小的细胞对(us,vs)作为匹配结果，有

(p,q=1，…，m)

(6)

图5 基于局部图动态匹配算法的相邻两帧细胞图像邻域结构动态变化过程Fig.5 Illustration of the dynamic local graph matching in two consecutive images

通过以上步骤追踪到了第二个细胞，接下来需要追踪余下的细胞。在追踪余下细胞时，将已匹配细胞都当作种子细胞，则邻细胞就变成每个种子细胞未匹配邻细胞的并集。在每次追踪过程中，参与比较的细胞数目随着已匹配细胞数目的变化而变化，如图5中紫色细胞的变化所示。同时，由于图像采集间隔较小，图像中细胞的位置发生剧变的可能小。因此，待匹配细胞的邻细胞匹配情况可以作为参考：图像It与It+1中待匹配细胞的邻细胞若已匹配，已匹配的细胞应是同一细胞；已匹配的邻细胞个数越多，则细胞追踪结果越可靠。距离函数可以增加已匹配邻细胞的个数，如下式D2中的f5所示，有

(7)

式中，ω1～ω4通过本文第1.2.5节部分求解。

此时部分种子细胞的邻细胞个数可能不相等，因此假设图像It与It+1中种子细胞的邻细胞分别为M1个和M2个，选取距离函数最小的细胞对(us,vs)为追踪结果，有

(p=1，…，M1；q=1，…，M2)

(8)

在细胞逐步匹配过程中，已匹配的细胞将作为新增加的种子细胞。在动态扩张的已匹配细胞邻域范围中，继续匹配距离函数最小的细胞对，直至距离函数D2大于阈值T1(阈值T1通过本文第1.2.5节部分求解)追踪结束。但已匹配的细胞从邻域中排除，减少计算量，提高追踪效率。细胞追踪过程及局部图的动态变化过程如图5(a)和(b)所示，其中紫色细胞中间的细胞为已匹配细胞，也就是种子细胞，紫色细胞是所有种子细胞的未匹配邻细胞并集，红色箭头代表追踪的顺序。可以看出，整个追踪过程是动态变化的，邻域的不断扩大使追踪结果更可靠，种子点个数的增加避免了因某个种子细胞匹配错误造成邻细胞追踪错误。

1.2.4分裂细胞检测

植物在生长发育过程中会不断进行细胞分裂，随着时间的推移，一个细胞会分裂成两个细胞。如图6(a)和(b)中红色区域所示，在t时刻到t+1时刻之间共有4个细胞发生了分裂。细胞发生分裂后，母细胞近似平均地分为两个子细胞，且细胞分裂时只能沿着t时刻到t+1时刻发生分裂，而不能从t+1时刻到t时刻发生分裂。同时，本研究的植物顶端分生组织细胞拥有相互紧挨的特殊结构，图像采集时间短，细胞与其邻细胞局部结构短时间内不会发生剧变。如图6(c)和(d)所示，母细胞的邻细胞与两个子细胞的邻细胞相对应(标号相同的为同一细胞)。因此，根据细胞分裂与细胞局部图的特点，细胞分裂检测步骤如下：

1) 在追踪后的图像It中寻找未匹配的细胞d，图像It+1中寻找未匹配的相邻细胞d1与细胞d2；

2) 判断图像It中细胞d的邻细胞是否与图像It+1中的两个细胞d1、d2的邻细胞对应；

3) 判断细胞d、细胞d1和细胞d2的面积是否满足下式：

图6 (a)和(b)相邻时刻细胞分裂示意；(c)和(d)分别为(a)和(b)中矩形区域放大图像Fig.6 (a) and (b) are examples of cell division in two consecutive images；(c) and (d) are the enlargements of the rectangular regions in (a) and (b), respectively

(9)

式中，Ad表示图像It中细胞d的面积，Ad1与Ad2分别表示图像It+1中细胞d1、d2的面积，T2、T3为面积阈值(通过本文第1.2.5节求解)。

1.2.5最优参数选择

不同距离函数是通过参数进行相加的，因此参数的选择对追踪有很大影响，本研究通过从粗到精对参数进行优化，具体步骤如下：

1) 从整个细胞图像数据集中选择一组细胞图像数据集作为训练集，在训练集中不断改变某个参数值，同时固定其他参数值进行细胞追踪，得到追踪准确率变化情况，在追踪准确率峰值左右各取一个值构成一个追踪准确率值的区间，将这个区间的参数范围作为最优范围。重复此步骤，获得各个参数的最优范围。

2) 从各个参数的最优范围中均匀采样，获得不同参数值；将各个参数的不同参数值进行组合，得到多个候选参数组合。

3) 在细胞图像数据集上，利用5折交叉验证，从所有参数候选组合中选择最优的一组作为最后的参数值，也就是最优参数值。

2 结果

实验对6组不同植物顶端分生组织(SAM)细胞图像序列A～F进行追踪(D～F是A～C经过配准后的图像序列)，其中，每个序列包含24帧图像，每帧图像包含100～150个细胞。A～C细胞图像序列都是由共聚焦激光扫描显微镜对植物顶端分生组织进行采集，如图1所示。其中每个序列中相邻两帧细胞图像采集时间间隔为3 h，采集周期为72 h。 实验在Windows 7下使用Matlab 2014完成，计算机主频为3.3 GHz，内存为4 GB。

2.1 配准图像序列中的细胞追踪

为了验证本研究提出的局部图动态匹配算法的准确性，实验对已配准的3个图像序列中的细胞进行匹配，之后将本算法匹配结果与人工匹配结果对比，得到追踪准确率，如表1所示。可以看出，与文献[18](表中的M2)、文献[17](表中的M3)中提出的算法的追踪结果相比，本研究提出的算法(表中的M1)追踪准确率提高了4%，分裂细胞检测准确率提高了8%，表明局部图动态匹配算法的追踪效果更好。

表1细胞追踪准确率对比(配准图像序列)

Tab.1Celltrackingaccuracycomparison(registeredimagesequences)

实验数据细胞追踪准确率/%分裂检测准确率/%M1M2 M3 M1M2 M3 D98.4794.7995.5892.0584.7887.11E95.8691.4991.5489.7479.2681.85F93.4288.1689.9984.7271.0372.41

2.2 未配准图像序列中的细胞追踪

图7给出了相邻两帧未配准细胞图像的追踪结果示例，(b)与(d)中相同颜色的细胞表示同一细胞。从(a)与(c)中可以看出，在成像过程中出现了旋转与错位；但从(b)与(d)中可以看出，大部分的细胞都能被追踪，且追踪准确性较高，说明本算法对未配准图像追踪比较可靠。

为了进一步证明本研究提出算法的有效性，用该算法对存在旋转与平移的3组未配准图像序列中的细胞进行追踪，得到的追踪结果如表2、3所示。表2是图像序列A中不同时间间隔的两帧图像追踪结果，表3是对3个图像序列中细胞进行追踪后的结果。由表2、3可以看出，与其他文献中的算法[17-18]相比，在未配准图像序列中得到的追踪准确率提高了30%，达到了90%。表明了本算法基本能对抗成像过程中带入的旋转与平移噪声。

图7 细胞追踪结果。(a)和(c)为源图像；(b)和(d)为追踪结果(图像中存在旋转与平移)Fig.7 Cell tracking results. (a) and (c) are the original images; (b) and (d) are the tracking results

表2不同时间间隔细胞追踪准确率对比(未配准图像)

Tab.2Celltrackingaccuracycomparisonindifferenttimeintervals(unregisteredimages)%

表3细胞追踪准确率对比(未配准图像序列)

Tab.3Celltrackingaccuracycomparison(unregisteredimagesequences)

实验数据细胞追踪准确率/%分裂检测准确率/%M1M2M3M1M2M3A92.9364.0164.5890.5654.7155.42B91.4660.6761.7288.8953.8452.31C89.9055.9954.7684.6151.3950.64

3 讨论

本研究利用局部图动态匹配算法对植物细胞图像中的细胞群进行追踪，根据细胞追踪准确率与分裂检测准确率对提出的算法进行分析。表1～3给出了本算法(M1)与文献[18](M2)、文献[17] (M3)提出的算法在追踪准确率上的比较结果，可以看出，本算法较之其他算法[17-18]在配准图像序列及未配准图像序列上准确率都有提升。

在配准图像序列上，本研究提出的算法最优。这主要是因为在细胞匹配过程中，文献[18]所用的扩散式追踪即以种子细胞为中心，匹配种子细胞所有邻细胞的方式容易造成错误累积；文献[17]利用基于条件随机场的追踪方法，以全局匹配的方式来解决植物细胞匹配的问题，虽然全局匹配方式避免了匹配过程中的错误累积，但是在全局求解时若出现异常值，将会造成追踪错误。本研究在继承文献[18]中局部图匹配的优势上，提出动态匹配细胞，每次匹配时将之前已匹配的细胞都作为种子细胞，然后从这些种子细胞的邻细胞中选择最相似的细胞对作为本次的匹配结果。动态匹配方式扩大了每次匹配时的候选细胞数目，且从中只选择一对最相似的细胞，这种匹配方式在一定程度上降低了错误匹配的风险。

此外，文献[17]与文献[18]中所提的算法在追踪前都需要对图像进行配准，然后再通过提取细胞位置等特征进行追踪与分裂检测，存在很大的局限性。而本研究所提的算法不需要对图像进行配准，直接对未配准图像分析，利用旋转、平移不变特征进行追踪与分裂检测，表2、3中的实验结果证明了本研究所提出算法的有效性。

通过对文献[17]、文献[18]与本研究所提算法的分析与讨论可以看出，算法能否追踪未配准图像与所提取的特征有很大关系，同时细胞匹配顺序是影响追踪准确率的一个重要因素。但与文献[18]中算法在追踪时每次匹配种子细胞所有邻细胞相比，每次只匹配一对细胞的动态匹配方式也增加了算法的运行时间。因此，如何在确保追踪准确率高的情况下提高运行效率，将会是下一步的研究方向。另外，基于共聚焦激光扫描显微成像技术采集的细胞图像序列包括时间上的图像序列和空间上的图像序列，但本算法只利用了细胞在时间序列上的信息。同一细胞在空间上是跨越多层的，因此若结合空间上其他层的细胞图像，则能利用的信息将更加丰富，可能会得到更好的追踪结果，尤其是对噪声区域的细胞追踪，这也将成为后续研究的重要内容。

4 结论

本研究提出了一种基于局部图动态匹配的植物细胞追踪算法。该算法利用细胞面积、相邻细胞之间的夹角与距离作为匹配的基本特征，先选择最相似的细胞对作为种子细胞对，然后从种子细胞出发，每次优先匹配特征距离最小的邻细胞对，同时向已匹配的细胞加入种子细胞，最后结合细胞局部图与细胞分裂特点对分裂细胞进行检测。与其他算法相比，在配准图像序列中，本研究提出的追踪算法的追踪准确率提高了4%，在未配准图像序列中追踪准确率提高了30%，表明该算法对图像中存在的旋转与位移噪声具有一定的鲁棒性。