CN101611411B - 追踪生物材料运动的系统和方法 - Google Patents

Abstract

本发明提供一种自动系统和方法，能使用户追踪对象的运动，例如多个细胞中的至少一个细胞随着时间从一个点到另一个点的运动，以使用户能确定该至少一个细胞如何工作。用户能追踪该至少一个细胞的运动并确定是否该至少一个细胞的结构发生了改变，例如该至少一个细胞正在经历细胞分裂或两个或更多细胞融合到一起。此用户能最优地追踪至少一个细胞来找到随着其在时间和距离上的移动，该至少一个细胞的大致位置。这样，本发明向用户提供一种通过追踪细胞运动来研究细胞功能的方法。

Description

追踪生物材料运动的系统和方法

相关申请的交叉引用

本申请要求于2007年2月9日提交的申请号为60/889,006的美国临时申请的优先权；该临时申请在此完整地引入作为参考。

发明领域

本发明涉及追踪生物材料的运动的系统和方法。

发明背景

通常，虽然已经有成功的基因组测序程序和其他类型的程序，大多数基因、细胞和病毒在人和其他生物体中的作用方式仍是相对未知的。因此，需要能使人学习和理解基因功能的高通量(high-throughput)筛选。高通量筛选使得研究和分析成百上千的基因产物成为可能。这些分析的结果使人可以研究细胞内发生的生物学过程，这对分析细胞是必要的，因为对生物学检测的科学研究而言，细胞追踪(cell tracking)是重要的。

典型地，需要使用计算机来追踪细胞，以对众多细胞中的单个细胞染色标记，并在细胞从一帧(frame)向另一帧运动的过程中追踪它。这一追踪被标记细胞的运动的过程容易出错且耗费时间，这是因为尽管这一个细胞已被标记，但仍然难以在从一帧至另一帧的众多细胞中被识别。此外，用于标记细胞的染料可能会妨碍正常的细胞周期，因此无法获得细胞从一帧向另一帧运动的过程中如何工作的真实图景。

有几种不同的用于计算机追踪细胞的算法，诸如卡尔曼滤波(Kalman filter)和粒子滤波(particle filter)。卡尔曼滤波是一种有效的回归滤波算法，它根据一系列不完整和有噪声的测量来估计动态系统的状态。卡尔曼滤波使用基于二阶统计量的追踪方法以追踪对象的运动，例如细胞从一点到另一点。卡尔曼滤波使用的追踪方法假定动态和测量模型与高斯噪声存在线性关系。然而，对于细胞图像序列，有很多干扰因素，例如背景杂物，因而使生成干净的细胞界线变得困难，而这经常会导致卡尔曼滤波追踪的崩溃。

粒子滤波是序列蒙特卡洛法的一种形式，它使用基于随机采样的复杂的模型估算技术。粒子滤波不需要假定高斯分布，动态和测量模型也可以是非线性的。粒子滤波不受到卡尔曼滤波的缺点的限制。但这样的代价是速度的降低。然而，尽管如此，粒子滤波并不慢。就速度而言，使用粒子滤波替代卡尔曼滤波相当于从极快变为对于大多数应用足够快。

粒子滤波和卡尔曼滤波的根本不同在于它们表示状态的方法。卡尔曼滤波为每个参数存储一个数值，并以协方差矩阵的形式存储每个参数之间的对应差异和相互关系。粒子滤波采用一系列“粒子”，它们各自对应于包含状态向量及其对应权重的组合。被追踪对象的整体行为可以通过统计粒子的集合来得到。单一对象可以通过使用几十个、几百个甚至几千个粒子来追踪。粒子滤波需要细胞具有特定形状以便保持对它的追踪，但细胞并非具有确定形状的物体，因此有时难以追踪细胞。

因此，为了研究细胞功能，需要有一种自动追踪细胞从而能生成描绘细胞运动的真实图像的结果的系统和方法。

发明概述

本发明基于上述技术背景完成，其目标是提供一种追踪细胞运动的系统和方法。

在本发明的一种优选实施方案中，给出了追踪细胞的方法。该方法包括：提供至少一个细胞在前一帧中带有位置的图像；根据前一帧中该至少一个细胞的位置来确定一段时间后该至少一个细胞在当前帧中的新位置；测量该至少一个细胞在当前帧中的新位置；以及根据测定的该至少一个细胞在当前帧中的新位置来更新该至少一个细胞在当前帧中的位置。

在本发明的另一优选实施方案中，给出了用于追踪细胞的方法。该方法包括：接收至少一个细胞在前一帧中的图像；提供与该至少一个细胞在前一帧中的图像相关的多个粒子，其中，多个粒子追踪来自前一帧的该至少一个细胞向当前帧中的至少一个细胞的运动；以及通过多个粒子来确定与带有(with)多个粒子的前一帧中该至少一个细胞相关联的带有多个粒子的当前帧中该至少一个细胞的大致位置。

在另外一种本发明的优选实施方案中，给出了用于追踪细胞的系统。成像系统被设置用于接收至少一个细胞在前一帧中的图像，其中成像系统连接到图像接收设备。图像接收设备被设置用于接收至少一个细胞在前一帧中的图像；提供在前一帧中与该至少一个细胞相关的多个粒子，其中多个粒子追踪来自当前帧的至少一个细胞到前一帧中的该至少一个细胞的运动；以及确定与带有多个粒子的前一帧中该至少一个细胞相关联的带有多个粒子的当前帧中该至少一个细胞的大致位置。

附图简述

结合附图阅读如下描述，将使本发明的这些以及其他优点更加明确，其中：

图1是依据本发明的实施方案的用于追踪细胞运动的系统的结构图；

图2是依据本发明的图1的图像接收设备的原理图；

图3的流程图描绘的是图1的细胞追踪系统如何根据本发明使用；

图4显示的是根据本发明生成众多细胞的视场(Field of View)的计算机屏幕截图；

图5显示的是根据本发明从众多细胞中生成单一细胞的视场的计算机屏幕截图，它显示这个细胞的特征；

图6A显示的是根据本发明得到的与图4中众多细胞的视场相关的特征的计算机屏幕截图；

图6B显示的是根据本发明挑选出的细胞的特征的计算机屏幕截图；

图6C显示的是根据本发明正在选择的输入特征的计算机屏幕截图；

图7描绘的是根据本发明的贴近度(Proximity)算法的流程图；

图8描绘的是根据本发明的粒子滤波算法的流程图；

图9是根据本发明生成双曲柯西密度的方法的图形描绘；

图10是描述根据本发明使用相对[逆向]距离来确定细胞融合或细胞死亡的表格的图形描绘；以及

图11是描述根据本发明使用相对[逆向]距离来确定细胞分裂或细胞生成的表格的图形描绘。

发明详述

本发明的这一优选实施方案结合图形描述，其中同样的(like)组成部分以相同的数字表示。优选实施方案的描述具有例示性，不意味着限定本发明的范围。

图1描绘的是传统荧光显微系统的结构图，其中包括了用于追踪细胞运动的系统。该自动显微系统100包括：可选的常规计算机103、光源105、光学检测器107、分光(dichroic)镜109、物镜111、载物台113、样本115以及样本架119。

物镜111下方是样本架119，可以为典型的微量滴定板(micro titerplate)、显微镜载物片、芯片(chip)、玻璃片、Petri氏培养皿、塑料或硅或其他任何类型的用于在载物台113上支持样本115的典型的支持物。本发明中也可以使用两个或更多个样本架119和/或两个或多个载物台113。

在另一实施方案中，显微系统100可以通过通信连接117以电或无线的方式连接到常规计算机103。通信连接117可以是能在自动显微系统101与计算机103之间传输数据的任意网络，比如局域网络(LAN)、无线局域网络、无线广域网络(WAN)、通用服务总线(USB)、以太网连接、光纤以及其他。显微镜也可以有多个物镜111。

显微系统100可以为影像传输设备，能抓取放在载物台113上的样本115或任意类型的对象的影像的成像设备或成像系统，这是通过使用光学检测器107或典型的显微目镜来完成的。此外，显微系统100也可以是例如由位于New Jersey的Piscataway的GE Healthcare生产的INCELL^TM Analyzer 1000或3000。样本115可以是活的生物组织、生物细胞、非生物样本以及其他。显微系统100可以是典型的共焦显微镜、荧光显微镜、落射荧光显微镜、相差显微镜、微分干涉差显微镜，或任意为本领域中具有一般技能的人员所知的显微镜。在另一实施方案中，显微系统100可以是典型的高通量和高含量的亚细胞成像分析设备，它能快速地检测、分析并提供生物体或其他的图像。显微系统100也可以是自动化细胞和亚细胞成像系统。

光源105可以是灯、激光、多束激光、发光二极管(LED)、多个发光二极管或任意为本领域中的一般技术人员所知的光源，它生成光线106，此光线通过主光路106到达载物台113以照亮样本115。对于激光扫描显微镜，扫描镜109位于样本115上方；该扫描镜109将来自光源105的光在显微镜的视场中扫描，以获取光学检测器107上的样本的图像。对于荧光显微镜，扫描镜109也可以是分光镜109，其向样本反射激发光，将来自样本的荧光导入光学检测器107。

光学检测器107接收来自样本的反射光或荧光，它可以是光倍增管、电荷耦合器件(CCD)、互补金属氧化物半导体(CMOS)图像检测器，或为本领域中一般技术人员所使用的任意光学检测器。光学检测器107通过通讯连接117以电或无线的方式连接到计算机103。在另一实施方案中，光学检测器107可以替换成典型显微镜目镜或目镜，用来与物镜111一起进一步放大中间图像以观察样本的细节。

计算机103可以为图像接收设备103或图像检测设备103。在本发明的另一实施方案中，图像接收设备103可以位于图像传输设备100内。图像接收设备103充当典型的计算机，它能从光学检测器107接收样本115的图像，接下来，图像接收设备103能使用标准图像处理软件程序、算法或公式来显示、保存或处理图像。计算机103也可以是个人数字助理(PDA)、膝上型电脑、笔记本电脑、移动电话、基于硬盘的设备或其他任何能通过通讯连接117接收、发送和存储信息的设备。尽管本发明中使用了一台计算机，但也可使用多台计算机来替代计算机103。

图2描绘的是图1中细胞追踪系统的图像接收设备的原理图。图像或成像接收设备103包括与常规计算机相关的典型元件。图像接收设备103也可以存储于图像传输系统100上。图像接收设备103包括：处理器103a、输入/输出(I/O)控制器103b、大容量存储器103c、存储器(memory)103d、视频适配器103e、连接接口103f和系统总线103g，它们通过电或无线的方式可操作地将前面提到的系统元件连接到处理器103a。系统总线103g也通过电或无线的方式可操作地将典型的计算机系统元件连接到处理器103a。处理器103a也可为处理单元、中央处理单元(CPU)、多个处理单元或并行处理单元。系统总线103g可以是与常规计算机相关的典型的总线。存储器103d包括只读存储器(ROM)和随机存储器(RAM)。ROM包括典型的输入/输出系统，其包括基本例程，在启动过程中协助计算机的元件之间传输信息。

输入/输出控制器103b通过总线103g连接到处理器103a，其中输入/输出控制器103b作为接口，能允许用户通过细胞追踪图形用户界面(GUI)和如键盘和定位设备的输入设备104将命令和信息输入到计算机。使用的典型定位设备有操纵杆、鼠标、游戏摇杆或其他。显示器201通过视频适配器103e以电或无线的方式连接到系统总线103g。显示器201也可以是典型的计算机显示器、等离子电视、液晶显示器(LCD)或任何能显示计算机103生成的字符和/或静态图像的设备。与计算机103的视频适配器103e相邻的是连接接口103f。连接接口103f可以是网络接口，它按照上述方式通过通讯连接117连接到光学检测器107。图像接收设备103也可以包括网络适配器或调制解调器，它们能将图像接收设备103连接到其他计算机。

存储器103d上方是大容量存储器103c，它包括：1.用于对硬盘执行读取和写入的硬盘驱动器组件(未显示)和硬盘驱动器接口(未显示)，2.磁盘驱动器(未显示)和硬盘驱动器接口(未显示)以及3.用于从诸如CD-ROM的移取式光盘或其他光学介质读取或写入的光盘驱动器(未显示)和光盘驱动器接口(未显示)。上述驱动器及其相关的计算机可读取介质为计算机可读的操作指令、数据结构、程序模块和其他用于计算机103的数据提供非易失性存储。上述驱动器也包括本发明的追踪细胞移动的算法、软件或公式的技术效果，这些将于图3、7和8的流程图中描述。

软件具备细胞追踪图形用户界面(GUI)。细胞追踪图形用户界面是专门编写的GUI，它具备和典型GUI相同的一些功能，典型的GUI是被设计用于使计算机用户容易地和计算机103交流的软件程序。细胞追踪GUI包括一张屏幕截图，它展示：1.图6A、6B和6C中讨论的细胞的特点(“测度(measures)”)。

图3描绘的是细胞追踪软件应用程序如何追踪细胞运动的流程图。该描述基于单个细胞，但此描述也可以用于显微系统100追踪的多个细胞。用户将多个细胞插入显微系统100，并通过成像接收系统103查看。

在块(block)301中，用户使用定位设备104来初始化存储于处理器103a上的协议中的追踪细胞的贴近度或粒子滤波算法，用户从而操作细胞追踪GUI，使用图6B中的屏幕截图上的指示器来核对(check off)追踪的细胞运动，以便选择贴近度或粒子滤波算法。接下来，在块303中，显微系统100处的用户观察多个细胞的视场(FOV)，称为“场景(Scene)”，在特定的状态或帧中也被称为图像接收设备103处的“前一帧”，如图4的计算机屏幕截图所示。可选地，细胞追踪图形用户界面可以在计算机屏幕截图中显示细胞的特征。参见图6A，计算机屏幕截图显示了细胞的特征。对显微系统100上的每一个细胞，基于由Piscataway，NJ的GE Healthcare生产的INCELL^TM Analysis Workstation或INCELL^TM Developer中使用的标准图像分析软件，每个细胞的特征都是已知的。

细胞的特征包括细胞核的密度(intensity)，或细胞核或细胞的核区的致密性等。这些特征意味着核的中心是必要的。大量的细胞特征可以被显示，这些均为本领域的一般技术人员所熟悉，在此不一一列举。每个特征被赋予权重(weight)，如图6C的计算机屏幕截图中所示，该权重描述每个特征对追踪过程的影响程度，它以一个大于0且小于等于1的数表示。

在此，用户也可以为与细胞相关的所有特征赋予特定的权重。例如，对于Nuc cg X的特征，用户可以指定1.00的权重，而对于Nuc cg Y，用户可以指定1.00的权重。Nuc cg X和Nuc cg Y中，X和Y分别是细胞核质心的坐标。细胞核质心相当于细胞核的“重量中心”。细胞核细胞用位于如下位置的M个像素表示：x_i1＜＝i＜＝M，y_i1＜＝i＜＝M，Nuc cg＝-x1/M M∑i＝1x_i，而Nucg＝-y＝1/M M∑i＝1y_i。用户也可以选择在块315处当细胞追踪过程完成时他想收到的输出数据，例如图6C的计算机屏幕截图上所示的标记(Label)、追踪ID(Track ID)、事件(Event)、置信度(Confidence)和距离(Distance)。标记是指定给细胞的唯一标记，它不是细胞在给定时间点的表示，而是代表了真实的物理实体。追踪ID是表示父子关系的数字，一般来说，一个细胞的所有图像结果(imaged progeny)都共享其追踪ID。事件与给定时间点的给定细胞相关，例如无、新生、分裂或碰撞。置信度是由算法生成的关于将当前帧中的当前时间点的指定细胞指定给前一个时间点中它的来源(parent)是否正确的置信度。距离是细胞在从前一帧到当前帧中移动的距离。

细胞追踪软件程序包括标准分割(segmentation)和量化程序，用于分割和量化多个细胞以便检测和分析来自多个细胞的图像中的至少一个细胞，并计算该至少一个细胞的特征，即特征向量(feature vectors)、细胞位置、面积、密度、比例、分类等。“比例”可能是恒定特征的一个例子——即使对象改变了尺寸但不改变形状，该数值依然保持不变。分类结果是更多的特征。例如，在其生命周期中，细胞会经过多个阶段，这些阶段可以用分类器来确定。根据生命周期中所处阶段，我们可以进一步使用不同组的特征。例如，在细胞分裂前，细胞的尺寸和亮度均会增长并会伸长。这三个特征，尤其是伸长或“形状(form)系数”，可以帮助辨别下一帧中的同一个细胞。用户能使用定位设备104(图2)来操作计算机103的细胞追踪图形用户界面(GUI)，以从多个细胞中选择至少一个细胞。

对于图5，在多个细胞经历量化和分割流程后，用户选择的在一段时间内被追踪的细胞将在计算机屏幕截图中显示，其中细胞会放大，细胞追踪GUI显示用于在时间点之间浏览的按钮“<<”和“>>”以及在同一帧中的兄弟细胞之间浏览的按钮“<”和“>”。细胞的计算机屏幕截图也包括一个用于追踪细胞的追踪ID 99以及用于与其他细胞区分该细胞的细胞编号98，如图4所示。

在块305中，图像接收系统103中的处理器103a确定正在查看的帧是第一帧还是图像接收系统103的视场(FOV)中的前一帧。如果该第一帧不是第一帧或前一帧，则在块307中进行帧匹配流程。

在块307中，用户能使用图像接收设备103操作细胞追踪GUI，用于预测在一段时间内，在该至少一个细胞从当前或现在的帧中被观察到的位置到块301中的前一帧中的细胞的变化。如图6C中所示，使用块701中的贴近度算法或块801中的粒子滤波算法，用户能追踪当前帧中的多个细胞中的至少一个细胞到前一帧中的多个细胞中的该至少一个细胞的运动。贴近度算法描述的是与前一帧中的至少一个细胞的位置相等的当前帧中多个细胞中的至少一个细胞(最近的细胞)的位置。粒子滤波算法描述的是当前帧中多个细胞中的至少一个细胞的位置，该位置与多个细胞中的该至少一个细胞的前一个位置或前一帧+速度×时间+1/2加速度×时间的平方+随机变化相等。该粒子滤波使用动态模型。贴近度和粒子滤波算法的目标是识别联系当前(当前帧)和过去(之前的任意已处理的帧；通常是前一帧)的事件。

图7描述的是如何使用贴近度算法预测细胞在从前一位置向下一位置移动时将到达何处的流程图。在块701中，细胞追踪软件生成了以字母T代表的新的帧或当前帧，其中包含了与用字母P代表的前一帧相同的多个细胞，帧的时长例如是10⁶ms到5×10⁹ms，时间间隔10⁴ms。这样，对于每个前一帧中被追踪的对象，当前帧中的与前一帧中的对象或细胞足够近的对象或细胞被当作“候选匹配”。当前帧T的特征向量(位置)以字母t表示，而前一帧P的特征向量(位置)以字母p表示。候选匹配是与被考虑为对象的对象最近的对象，或被考虑的对象与其最近匹配之间距离不远于q的匹配。符号p_i的一组候选匹配表示为M_i。当前帧中与前一帧中的最近目标距离也称两个向量间的标准欧几里德距离，其中每个坐标按上所述进行加权。为了获得当前帧和前一帧间的距离，使用了如下公式：

1.m_i＝m(p_i)＝min_t∈T(||wp_i，wt||)，p_i∈P

2.M_p≡{p∈P：||wt_i，wp||≤qm_i}

3.C≡{(p，t，d)，p∈P，t∈M_p，d＝||wt，wp||}

第一个公式中，m_i表示特征向量p_i和来自T(下一帧，用加权距离衡量，它与p_i最近)的特征向量之间的距离。第二个公式表示的是包含T中每个特征向量t的P中每个向量p的候选匹配集合，这样，t和T之间的加权距离就小于等于q乘以最小距离(如公式1中定义的)。第三个公式形式上表示P和T帧之间的映射C，它表示前一帧和当前帧之间的多个细胞的所有匹配的信息，每组三个数据，包括：P的特征向量p，T的特征向量t(提取自p的候选匹配集合Mp)，以及p和t的加权距离。C是包含以前帧中的向量p、被认为是P的候选匹配的t以及对应距离d的所有元组的集合。两条竖杠的括号并不表示数的绝对值。竖杠表示的是向量的欧几里德范数(norm)，即点wp和wt之间的向量长成空间(span)。

块703中，细胞追踪软件中，字母D被定义为取自C的元素的加权距离的全域(元组中的第三个元素)。C的所有符合距离d(m)＞med(D)+10h_out·mad(D)的元素m均作为异常从集合中移除。

常数h_out由用户提供。符号med表示全域的中位数。符号mad表示绝对偏差的中位数。

同样，在块703中，如上所述，事件被识别。例如，根据以来自P和T的与其各自的特征向量p和t相关的特定元素开始或结束的一系列值，事件可以是碰撞或分裂。在块705中，根据其特征向量p和t，不与P或T直接关联的对象为孤立(orphaned)对象或分别产生分离或到达事件的细胞。

同样，在块705中，为每个匹配记录分配了置信度指标。置信度的定义是任意的(arbitrary)，不和匹配过程直接相关，而是与输入条件相关。根据以下公式，分配了多种事件：

4.无事件和碰撞事件：1-d/s_p，其中，s_p是p和与T第二近的匹配之间的距离。

5.分裂：b_t/s_p，其中b_t是T和P中与之最近的向量之间的距离。

6.到达：1-m_q∈P/b_t

7.离开：1-m_q∈P/m_p

对于无事件和碰撞事件的公式，s_p是p和与T第二近的匹配之间的距离。在分裂公式中，b_t是t和P中与之最近的向量之间的距离。对于到达和离开事件的公式，全域D的平均值除以b_t或m_p分别在数值1中扣除。在此，当前帧中的对象或该至少一个细胞即将通过块309进行标记。

在块309中，来自当前帧或T的该至少一个细胞基于前述过程中列举的事件被标记。这一标记步骤将指派标记或名字，它们将在对象的生命期保持不变，并用于标明物理对象的图像。标记将使程序能以能让用户掌握这些联系的方式输出数据。这些标签是以一个场景的第一帧中从1开始的序列的自然数，其中每个场景分别标记。同样，这一标记过程也识别了来自同一细胞或前一帧中的细胞后代的世系树。与对象标签类似的，世系树也以自然数标记。标签根据发生在当前帧和其前身之间的事件被分配给子对象，其规则如下表中所示：

接下来，在块315中，从前一帧到当前帧的输出数据在一段时间内显示。该输出数据如下表中所示。此表格用表格的形式描述了当前帧和前一帧或父子关系。该表格中的信息也可包括事件标签、细胞ID和父子关系树ID(parentage tree ID)。

如果被考虑的帧是场景中的第一帧，需要使用块311中的追踪算法对数据结构进行初始化。

接下来，在块313中有对被追踪进入上述多种事件的对象的标记，这些事件包括碰撞、到达、离开等。在块315中，如上所述，当前帧和前一帧在一段时间内的输出数据被显示。在块317，图像接收设备103确定图像接收设备103是否取得了场景的最后一帧。必要地，用户告知图像接收设备103或查看图像接收设备103上是否有用于检查若干个当前帧的数量在1到100帧，优选地，60到80帧的范围内。例如，100,000-1,000,000ms时间段内的100个当前帧，用于获取对象或视场内多个细胞中的细胞如何在时间段内运动的精确描述。如果确定该场景的最后一帧并未获得，流程回到块303。如果确定获得了场景的最后一帧，且输出数据也已获得，流程则进入块319。在块319，当前帧和前一帧之间的中间数据将被删除，此流程结束。

本发明中使用的粒子滤波算法包括两个独立的层，一个内层和一个外层。内层的每个组件可以被视为职责是追踪单个目标或细胞的次级追踪器。外层可以被视为职责是追踪所有目标或细胞的追踪器。内层包括一个“粒子”云(cloud)，它像一个蜂群一样，追寻一个简单的目标，例如一个对象或者一个细胞。外层包括一组云，其中的每一个均试图追寻一个目标。内层和外层均有通常的预测、观察/测量和更新步骤。外层有通过对其云构成中的每一个进行指派预测来预测每个云的行踪的预测阶段。换句话说，外层自己不进行预测，而内层自己做预测。

外层也有一个包括两个步骤的测量阶段。首先，对每个目标的每个云，云到每个目标的距离被计算出来。外层测量阶段假定了两点。首先，它假定计算机103知道如何测定每个云和每个目标之间的“距离”，这将在后面关于交叉熵(cross entropy)的讨论中阐明。第二点是计算机103将知道使用这些距离做什么，例如确定事件、无事件、融合、分裂、细胞新生和细胞死亡。在测量阶段的末尾，计算机103应该有与对象或细胞相等数量的粒子。直至这一末尾，计算机103需要在以下之间确定和识别：1.细胞分裂、2.细胞融合、3.细胞消失和4.细胞新生。关于对象和云之间的距离表的对象分析如图10和11中所示。

外层测量阶段的下一步包括计算机103构建一个使用步骤1中计算出的距离的“责任矩阵”(Responsibility Matrix，RM)。接着，计算机103针对细胞的分裂、融合、碰撞等分析该RM。接着，计算机103确定与每个靶标的一组粒子云相等的输出。

外层更新阶段指派每个粒子云更新其自身状态。外层测量阶段是相当复杂的，并包括非常多的新的方面。外层操作图像中的全部多个对象，其中对象以独立的分割算法定义，该算法诸如由位于Piscataway，NJ的GE Healthcare生产的INCELL^TM 1000Developer Toolbox所提供的。

与外层类似，内层是一个典型的贝叶斯(Bayersian)追踪滤波器，它包含了独立的预测、观察和更新阶段。粒子滤波器的内层包括用于追踪简单对象的粒子云。内层的预测阶段为：

预测：

x_t←p(x_t|x_t-1)            理论

漂移：x′(t)＝x(t-1)+vx(t-1)

y′(t)＝y(t-1)+vy(t-1)    实际

扩散：x(t)＝x’(t)+ex(t)

y(t)＝y’(t)+ey(t)

上述例子显示的是一个使用了零加速度的简单模型。总的来说，我们可以使用计算x和y的新的或最终位置模型，该模型考虑到三种不同的动态模型：1.最终位置＝初始位置+速度×时间+噪声；2.最终位置＝初始位置+噪声；以及3.最终位置或3.初始位置+速度×时间+.5×加速度×时间×时间+噪声。

预测阶段的理论表明，计算机103应该从密度p(x_t|x_t-1)直接取得样本，脱氧核但这首先是不可能的。在实际中，该密度对计算机103来说无法得到。为了避免与使用卡尔曼滤波器相关的缺陷，谨慎地选择了一种不同的密度，其本身不基于高斯统计。所选的密度是一个长尾(long tailed density)密度。概率密度函数的使用影响到细胞追踪器的稳健性。稳健统计领域背后的原理是将一个或多个异常值能产生的影响最小化。其中的一个方法是将估计基于厚尾密度(thick tailed densities)，其例子之一是图9中的双曲柯西密度。请注意在我们已经追踪大量细胞之后，实施理论上的采样成为可能，该采样由粒子滤波器的分别称为x_t←p(x_t|x_t-1)的预测阶段控制，如M.Isard和A.Blake的Isard和Blake在此文中所述：“Condensation--conditional density propagation forvisual tracking”，International Journal of Computer Vision 29(1)，pp.5-28，1998，http://citeseer.ist.psu.edu/isard98condensation.html

双曲柯西密度的对数概率图显示出与高斯对数概率不同的表现。对于大距离，其表现是线性而非抛物线型的。对于小距离，这两个测量几乎是相同的，但对于大距离，双曲柯西更宽容。这是稳健估计的要素，因为异常值的权重更小。概率分布函数中重的尾部意味着异常值将被给予更小的权重。这一密度对小偏差使用L2范数，对大偏差则使用L1范数。

高斯统计不稳健，这是因为相关距离测量会将异常值平方。净结果是一个非线性的追踪器，对于小距离使用高斯，对大样本则准确地使用非高斯。这应当能提高对异常值的稳健性。特征之间的距离通过我们使用更稳定的“距离”的量度来隐含地处理了。

H-Cauchy：小x＝加平方

大x＝加绝对值

重要性抽样是使用重尾(heavy tailed)密度的另一个好理由。关于卡尔曼滤波器和粒子滤波器的讨论一般因被大量多维积分的存在导致过于复杂。粒子滤波可被认为是数值蒙特卡洛积分的一种使用重要性抽样方法的形式。在重要性抽样中，一个好的有效建议密度能尽可能地接近真实统计。然而，为了数值稳定性起见，重尾的建议密度是可取的。

综上所述，重尾能达到三个目标：

1.对异常值的稳健性(即突发的大的变化)

2.在查找和探寻新状态方面更积极

3.更稳定的数值积分

接下来，粒子滤波器的内层包括一个测量组件。理论是等于p(z_t|x_t-1)的权重。这一实践与密度的概率密度函数在封闭形式下存在的理论相同。在这一阶段，每个粒子获得一个取决于用概率密度函数衡量的其和目标的接近度的权重。

最后，内层包括一个更新组件，在此样本会根据其各自的权重重新采样。生成好的预测的例子被进行多于一次的采样。生成不好的预测的例子将获得小的权重并可能根本不被采样。更新步骤对粒子重新采样和替换。如下所示的Matlab代码所示，这是相当容易实施的：

function sample＝sampleUnif(nsamp，wts)

     csum＝cumsum(wts)；

     rng＝csum(length(csum))；

     sample＝zeros(nsamp，1)；

     for j＝l:nsamp

           choice＝rng*rand(1)；

           l＝find(csum＞＝choice)；

           sample(j)＝l(1)；

     end

在产品中使用的相应的C语言代码更好，这是因为基于CDF是单调增加且预先排序的这一事实，代码中使用了二进制搜索。

重新采样的实施例

输入：(9，0.167)(16，0.333)(7，0.000)(20，0.000)(18，0.500)(23，0.000)

输出：(9，0.167)(16，0.167)(16，0.167)(18，0.167)(18，0.167)(18，0.167)

第一个数是被采样的，第二个数是对应的权重，在输出中，9出现一次，16出现两次，18出现三次。20、7、23没有出现，因为它们的权重为0。请注意权重被重新正则化以使重新采样后的总和保持为1。在这一例中，值9、16、7、20、18、23为标量，但没有理由表明它们不能是向量。一般来说，它们恰恰就是向量。

在图8的块801中，计算机初始化粒子滤波算法，为每个目标分配其自己的追踪器。追踪器的起始状态是基于与细胞图像的初始位置和状态相关联的细胞的可观察属性，例如细胞的位置、图像密度、细胞面积、形状系数，或一般而言，其他可测量的特征向量等。在块803，基于前述的可观测属性，通过推断比如速度等隐藏属性的数值，以其他贝叶斯学习算法相同的方式，粒子滤波算法能预测细胞图像的位置。

细胞追踪软件程序包括典型的蒙特卡洛粒子滤波算法，它使用户能追踪前一帧中在其初始状态或位置细胞图像到当前帧中的新的状态或位置，其中蒙特卡洛程序的多个粒子像蜂群一样在一段决定于用户的时间内从实质上跟随原细胞从其起始位置到新的位置。粒子滤波器可以包括多达一百个或更多个“粒子”，它们形成粒子云，追寻被指定的细胞的新位置的目标。

细胞追踪软件程序将例如双曲柯西密度等式的标准概率密度函数(probability density function，p.d.f.)用于在当前或下一帧中的该至少一个细胞，以便从前一帧中的多个细胞中确定该至少一个细胞在当前帧中的大概位置。块803中，算法也估算该至少一个细胞在当前帧中的位置。概率密度函数的使用影响了细胞追踪器的稳健性。稳健统计背后的原则是将一个或多个异常值的影响最小化。达到这个目的的一个办法是基于厚尾密度进行估计，而双曲柯西密度正是上述厚尾密度的一个例子。

接着，在块805，在单个目标被追踪的粒子滤波算法内层，基于特征向量、细胞密度、细胞方位等的细胞的新位置被查找，这是基于前一帧中细胞图像初始位置到当前帧中细胞图像的新位置的距离或基于相对于特征空间内细胞被观察到的位置基于粒子的对数概率的一般范数(不只是2D空间，而是定义了特征的多维空间)。这一距离可以是典型的欧几里德距离(L₂范数)或Metropolis距离(L₁范数)或特定的基于下述对数概率的一般距离，该距离能在块807，即重采样的更新阶段提供每个重采样需要的粒子的权重。

对于外层，即块305的追踪器管理系统，计算机103需要计算上述的取决于距离的一般测量的责任矩阵。本领域内的技术人员众所周知，似然函数的对数等效于距离。经常讨论的经典例子是关于独立相同分布数据的高斯似然函数。在这一特例中，等效距离量度是欧几里德范数，而最小自乘的标准方法等效于将似然函数最大化(即距离等效于负的对数概率)。

我们将构成云的粒子之间的总距离用负的对数概率总和表示：

$d_{\mathrm{jk}}=-\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\log w_{\mathrm{ijk}}$

其中w_ijk是关于给定粒子在特征空间内的每个细胞的位置的概率密度函数。(请注意，将对数概率求和等效于将概率相乘，这是计算总概率的恰当方式，因为粒子在统计学上是相互独立的)

求和过程实际上是伪装的蒙特卡洛积分。如果粒子云中的粒子取自密度p(x)，我们将得到：

$d_{\mathrm{jk}}=-\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\log w_{\mathrm{ijk}}\&ap;-\&Integral;p\left(x\right)\log w\left(x\right)\mathrm{dx}$

等式的右边是交叉熵。

进一步地，注意交叉熵可以被分解为：

$-\&Integral;p\left(x\right)\log w\left(x\right)\mathrm{dx}=\&Integral;p\left(x\right)\log \frac{p\left(x\right)}{w\left(x\right)}\mathrm{dx}-\&Integral;p\left(x\right)\log p\left(x\right)\mathrm{dx}$

-∫p(x)log p(x)dx是自熵。熵通常被认为是混乱的量度。在这里，它是密度p(x)的统一性的量度。粒子滤波器的重采样步骤(即“更新”步骤)将粒子的权重重设为统一的数值，从而使每个粒子等效。

是相对熵(Kullback Leibler Divergence，KLD)。它通常被用作概率密度之间距离的量度。(严格地说，它并不是距离的真实量度，因为它并不对称)。当p(x)和w(x)相同时，KLD最小。将我们对“距离”的测量最小化将使KLD最小化。例如当目标位于在粒子云的正中时，将产生最小KLD(因此粒子群的分布是依照具有最小偏移的粒子的真实分布)。在此情况下，我们可以作出结论：基于追踪器的粒子云锁定在目标上。

在图10中是五个云和六个目标的选定列最大值交叉熵。云5对应于目标5和6。我们推断可能曾有分裂。对图10，选定行最大值，目标6没有匹配的云。对数据的另一个可能的解释是云5正在追踪目标5，而目标6是全新的。解决这一冲突的办法是使用表格中的对数概率来计算分裂的概率并将它与被追踪的单个细胞的概率进行比较，例如通过使用贝叶斯决策理论或贝叶斯网络理论来将这些概率转化成对应的图形结构。

在图11中，对于选定行最大值(Maxima)，上方表格的上半部显示已知目标和存在的粒子云之间的交叉熵。表格的下半部显示交叉熵的行最大值以及其发生位置。云6没有匹配的目标，这意味着云3和6原来追踪的目标可能已经融合。

必须对云3和6对应的两个目标是否融合，或作为替换，仅仅是云6追踪的目标消失了作出判断。另一个(可能性小得多)解释是云3追踪的目标消失了，而目标5正在被云6追踪。可以通过对联合计算表格中条目求和来获得每种可能的概率，以便得到必需的边际概率或条件概率。正确的推断可以通过标准贝叶斯定理来正式地获得。注意交叉熵是上方每个细胞的概率的对数。对表格的粗略检查似乎可以表明正确的判断是被云6追踪的目标消失了，因为云6对目标5，或对任意其他目标的总体概率仅有很小的贡献。

对于列最大值，表格显示已知目标和存在的粒子云之间的交叉熵。下半部分显示交叉熵的列最大值以及其产生的位置。云6没有匹配的目标，这意味着该目标很可能消失了。在块809，将确定图像的有限序列中的最后一帧是否已产生，如果是这样则过程结束，否则过程回到块803。

在一些情况下，细胞融合并非是因为它们已经从物理上融合了，而是因为由Piscateway，NJ的GE Healthcare生产的叫做INCELLDeveloper toolbox的软件程序无法识别两个细胞，因为它们离得很近。这样的情况下，可能需要选择性地使用追踪器或软件来识别出事实上存在着两个细胞，并将它们分离。

本发明提供一个允许用户追踪对象的运动的自动化的系统和方法，例如追踪多个细胞中的至少一个细胞在一段时间内从一个点向另一点的移动，以便帮助用户确定该至少一个细胞如何工作。用户能追踪该至少一个细胞的运动并确定是否该至少一个细胞的结构发生了变化，例如该至少一个细胞正在进行细胞分裂，或两个或多个细胞发生了融合。这一用户能选择性地追踪该至少一个细胞来找到它在时间和距离上移动时的该至少一个细胞的大致位置。因此，本发明向用户提供了一种通过追踪细胞运动来研究细胞功能的方法。

尽管上面使用特定的实施方案描述了本发明，对本领域的技术人员显而易见地，在不脱离其如下面权利要求中所述的实质和范围的情况下，本发明可以有多种修改和变化。

Claims (1)

1.一种用于从前一帧到当前帧的多个细胞中追踪单个细胞的方法，所述方法包括：

确定在前一帧中单个细胞的位置和标记；

确定在当前帧中单个细胞的位置；

基于已经确定的在各帧中单个细胞的位置，确定联系当前帧中单个细胞与前一帧中单个细胞的事件，所述事件包括：

无事件、到达事件、离开事件、分裂事件以及碰撞事件；

基于所述联系事件以及前一帧中细胞的标记对当前帧中单个细胞进行标记；

其中所述确定事件的步骤包括采用粒子滤波算法，所述粒子滤波算法包括内层和外层；

其中所述内层的每个组件包括一个“粒子”云，配置用于追踪单个细胞在前一帧中的位置至在当前帧中的新位置；

所述外层包括对于每一内层的粒子云，测定当前帧中每个细胞的距离，并且由所述距离确定事件。

Images