CN103903015A - 一种细胞有丝分裂检测方法 - Google Patents

Abstract

本发明公开了一种细胞有丝分裂检测方法，包括以下步骤：图像矩阵向量化；新矩阵的构造；矩阵低秩分解；稀疏矩阵的每一列构造新矩阵；矩阵求和；低秩分解过程迭代；矩阵与操作。本发明所公开的细胞有丝分裂检测方法，将有丝分裂细胞看作低秩表示中的稀疏部分，非有丝分裂部分看作低秩部分，和已有的基于模式识别的检测方法相比，不需要进行大量的预训练，能够非常快速的得到检测结果，并且在相同的实验条件下，能够获得更高的F-measure和Recall值。

Description

一种细胞有丝分裂检测方法

技术领域

本发明公开了一种细胞有丝分裂检测方法，涉及图像信息处理技术领域。

背景技术

据国家癌症中心和卫生部疾病预防控制局公布的2009年乳腺癌发病数据显示：全国肿瘤登记地区乳腺癌发病率位居女性恶性肿瘤的第1位，中国女性乳腺癌发病率在逐渐上升，目前已从五年前0.017%增加到去年0.052%，上升超过三倍。另据《2012中国肿瘤年报》显示，目前我国肿瘤的发病率为285.91/10万，平均每天每分钟有6人被诊断为恶性肿瘤，乳腺癌在全国恶性肿瘤发病率中排名第6位，其发病率呈逐年上升的趋势，并且发病年龄有年轻化的趋势。乳腺癌已成为当前社会的重大公共卫生问题。自20世纪90年代全球乳腺癌死亡率呈现出下降趋势；究其原因，一是乳腺癌筛查工作的开展，使早期病例的比例增加；二是乳腺癌综合治疗的开展，提高了疗效。乳腺癌已成为疗效最佳的实体肿瘤之一。

图1是病理医生诊断乳腺癌的流程图。在乳腺癌的诊断过程中，病理医生首先从乳腺病变区域中提取组织样本，提取的样本被制作成组织切片。切片通常由H&E(Hematoxylin&Eosin)进行着色。H&E着色是病理学家广泛使用的方法，距今已有一百多年的历史，着色的目的是显示细胞的成份，H着色使得细胞核呈现蓝色，而E着色使得细胞质和结缔组织呈现粉红色。随后病理学家运用高倍显微镜来分析这些活组织切片并根据诺丁汉评分系统（Nottingham gradingsystem）对乳腺癌进行病理分级。

诺丁汉评分系统是由世界卫生组织推荐的乳腺癌国际评分系统。该系统综合评价了H&E着色切片的三种形态学特征：腺小管的构成、细胞核的多形性以及有丝分裂的个数。病理医生根据以上提到的三种形态学特征的良性和恶性程度给出1-3分的打分，最后三个分值的总和（在1-9分范围内）便是该患者的诺丁汉评分系统的分值。医生最后根据该分值的高低确定相应的治疗方案。因此病理医生对三种形态学特征的准确打分非常重要。

现有的基于数字病理切片计算机辅助自动分析方法主要集中在细胞核的检测，也有部分腺小管的检测结果。对乳腺癌评级而言，有丝分裂细胞的个数是一个重要的指标。因为它对肿瘤的扩散和侵略性进行了评估。在临床实验中，病理学家对H&E着色的切片在显微镜放大的情况下人工的计算有丝分裂的个数，通常放大40倍。在40倍视角下的可见区域被称为高能量区域（high power field，HPF）。当前的病理诊断都是基于病理学家的个人经验，对于同一个切片，不同的病理学家由于经验的差别等因素对病理图像的人工分析具有差异性，即使是同一个病理学家在不同的时间段受疲劳程度以及心情好坏等的影响也会产生不同的诊断结果。针对这种情况，迫切需要一种基于图像分析的自动化分析工具（如图1中的虚线矩形框）来对数字病理学切片进行量化评估。这种量化的描述不仅能帮助临床诊断（比如减少不同的观察者或者同一个观察者在不同时间段之间的诊断差异性），而且能够帮助病理学家理解使用某种具体诊断方法的潜在原因。此外，还能帮助研究人员理解疾病产生的生物学机理。

有丝分裂细胞的检测是一个非常具有挑战性的课题。一方面，对于研究者来说缺乏足够多可利用的数据；另一方面，有丝分裂细胞本身比较小并且形状和纹理变化很大，因此很容易跟其他物体、甚至是切片获取过程中引入的噪声混淆。有丝分裂的四个主要阶段为有丝分裂前期、中期、后期和末期。在不同的阶段，细胞核的形状变化非常大。在大多数阶段，有丝分裂的细胞核和非有丝分裂的细胞核非常相似，或者是类似于其它深蓝色的点，没有经过专业训练的人很难区分它们。另一个复杂的问题是，在有丝分裂的末期，细胞核会分裂成两个深蓝色的斑点，即一个有丝分裂会有两个不同的细胞核，但它们还不是完整意义上的单个细胞。一个有丝分裂的末期被认为是一个单一的有丝分裂，不能认为是两个有丝分裂。

目前，只有少量关于有丝分裂自动检测的工作。Belien等人考虑在福尔根（Feulgen）反应着色的乳腺癌区域计算有丝分裂的次数。Schlachter等人对荧光反应（fluorescence）着色的结肠直肠癌进行有丝裂检测。Roullier等人提出对免疫组化着色的乳腺癌切片进行有丝分裂细胞检测，这种着色方法能够突出具体的有丝分裂细胞，方便了检测。最近，不少学者尝试用模式识别的方法来对H&E着色的切片进行处理，以区分有丝分裂细胞和非有丝分裂细胞。Irshad等人利用形态学特征和纹理特征对每一个候选目标构造一个143维的特征集，然后用决策树分类器来对每一个目标进行分类；Tek等人结合常用的形状特征和颜色特征来描述目标，然后用级联的Adaboost对目标自动分类；这些方法对特征的选择非常敏感。针对这种情况，Malon等人提出利用卷积神经网络的方法来自动获取有效的特征表示，并将这些特征和手动设计的细胞核特征结合起来，卷积神经网络的使用从某种程度上来说减少了手动筛选样本的敏感性。以上基于模式识别的方法第一步通常需要将组织病理图像中的有丝分裂细胞和非有丝分裂细胞分割出来，然后再对每个分割块进行特征提取。然而一个合适的分割定义非常难而且费时；另一方面，要训练一个泛化能力较强的分类器，需要大量的训练样本，这对本来就稀少的数据来说具有很大的挑战性。

发明内容

本发明所要解决的技术问题是：针对现有技术的缺陷，提供一种细胞有丝分裂检测方法，该方法不需要进行复杂的特征提取，不需要先对图像进行分割，也不需要进行大量的预训练，只需要简单的灰度特征，使用起来非常快速方便。

本发明为解决上述技术问题采用以下技术方案：

一种细胞有丝分裂检测方法，具体包括以下步骤：

步骤1、用多光谱显微镜对有丝分裂中的细胞进行不同波段的扫描，并用不同的聚焦平面来数字化扫描得到的图片，形成图片矩阵；将扫描得到每幅图片所对应的图片矩阵分别用一个列向量来表示，该列向量的长度与原图片矩阵大小相同；

步骤2、用步骤1产生的列向量构造第一类矩阵，所述第一类矩阵的每一列分别按顺序对应步骤1中产生的每一个列向量；

步骤3、对步骤2得到的第一类矩阵进行低秩分解，分解成大小相同的低秩矩阵和稀疏矩阵；

步骤4、将步骤3得到的稀疏矩阵的每一个列向量都分别转化成相应的第二类矩阵；

步骤5、对步骤4中得到的所有第二类矩阵求和，得到第三类矩阵；

步骤6、循环执行步骤1至步骤5直到完成所有波段的数据处理，得到复数个大小相同的第三类矩阵；

步骤7、对步骤6中得到的复数个第三类矩阵进行与计算，得到最终的检测结果，其中的非零部分即为有丝分裂部分。

作为本发明的进一步优选方案，步骤1中将图片矩阵转化成列向量时，像素点的个数保持不变，大小为M×N的图片应转化成长度为M×N的列向量，其中M、N分别为自然数。

作为本发明的进一步优选方案，步骤3中对矩阵进行低秩分解，所得到的低秩矩阵和稀疏矩阵其大小均与低秩分解前的原第一类矩阵的大小一致。

作为本发明的进一步优选方案，步骤3中所述低秩分解的具体计算方法如下：

D＝[I₁,...,I_n]∈R^m×n表示n个连续的细胞有丝分裂病理图像，其中j∈n，I_j∈R^m表示第j幅图像，每幅图像由包含m个像素的列向量组成；矩阵D分解成低秩矩阵B∈R^m×n和低秩表示中的异常值S∈{0,1}^m×n两部分，B为背景图像，代表非有丝分裂部分，S为前景图像，代表有丝分裂细胞；

其中，ij表示第j幅图像中的第i个像素；

P_S(X)表示矩阵X到由S支持的矩阵线性空间的正交投影：

是与P_S(X)互补的投影，满足

给定一个实数阈值K，则矩阵B应该满足条件：

rank(B)≤K     (3)

其中，K越小，矩阵越相关；

通过马尔科夫随机场来建模：由矩阵D中所有像素点作为顶点υ构造一个图模型G＝(υ,ε)，υ的大小为m×n，ε是边的集合，每条边都是用来连接时间或空间域上的相邻像素，前景图像S应该满足：

$\underset{\mathrm{ij}\∈\mathrm{\υ}}{\mathrm{\Σ}}{u}_{\mathrm{ij}}\left(S_{\mathrm{ij}}\right)+\underset{(\mathrm{ij},\mathrm{kl})\∈\mathrm{\ϵ}}{\mathrm{\Σ}}{\mathrm{\λ}}_{\mathrm{ij},\mathrm{kl}}|S_{\mathrm{ij}}-S_{\mathrm{kl}}|---\left(4\right)$

参数λ_ij＞0是S_ij＝1的惩罚因子，表明前景部分是稀疏的；参数λ_ij,kl表示S_ij和S_kl之间的相关性，表明前景部分是连续的；在所有的位置都令λ_ij＝β，λ_ij,kl＝γ，其中β＞0和γ＞0是正常数；

将（4）式和（5）式合并为：

$\mathrm{\β}\underset{\mathrm{ij}}{\mathrm{\Σ}}{S}_{\mathrm{ij}}+\mathrm{\γ}\underset{(\mathrm{ij},\mathrm{kl})\∈\mathrm{\ϵ}}{\mathrm{\Σ}}|S_{\mathrm{ij}}-S_{\mathrm{kl}}|---\left(6\right)$

在完成对前景图像S和背景图像B建模后，得到信号模型D_o；

当S_ij＝0时，假设D_ij＝B_ij+ξ_ij，其中，ξ_ij为独立同分布的高斯噪声，即D_ij～N(B_ij,σ²)，σ²是高斯噪声的方差；

通过最小二乘法来拟合D_ij，当S_ij＝1时背景完全由前景遮挡，D_ij等于前景的灰度值；

通过最小化以下能量函数来估计B和S：

$\begin{array}{c}\underset{B,S_{\mathrm{ij}}\∈\left\{\mathrm{0,1}\right\}}{\min }\frac{1}{2}\underset{\mathrm{ij}:S_{\mathrm{ij}}=0}{\mathrm{\Σ}}{(D_{\mathrm{ij}}-B_{\mathrm{ij}})}^2+\mathrm{\β}\underset{\mathrm{ij}}{\mathrm{\Σ}}{S}_{\mathrm{ij}}+\mathrm{\γ}\underset{(\mathrm{ij},\mathrm{kl})\∈\mathrm{\ϵ}}{\mathrm{\Σ}}|S_{\mathrm{ij}}-S_{\mathrm{kl}}|,\\ s.t.\mathrm{rank}\left(B\right)\≤K.\end{array}---\left(7\right)$

利用矩阵的三个范数：

表示l₁范数；

表示Frobenius范数；||X||_*表示核范数，它是矩阵奇异值的和；

用核范数来代替矩阵B的秩操作符，将（7）式改写为矩阵形式：

$\underset{B,S_{\mathrm{ij}}\∈\left\{\mathrm{0,1}\right\}}{\min }\frac{1}{2}{\left|\right|P_{S^{\⊥}}(D-B)\left|\right|}_F^2+\mathrm{\α}{\left|\right|B\left|\right|}_{*}+\mathrm{\β}{\left|\right|S\left|\right|}_1+\mathrm{\γ}{\left|\right|\mathrm{Avec}\left(S\right)\left|\right|}_1---\left(8\right)$

其中，A是图G的节点-边的关联矩阵，α是和K相关联的大于零的参数，用来约束背景模型的复杂度。

本发明采用以上技术方案与现有技术相比，具有以下技术效果：

⑴在同样的实验条件下，和现有基于模式识别的检测方法相比，能够获得更高的F-measure和Recall值；

⑵本发明方法不需要进行复杂的特征提取，只需要使用简单的灰度特征；

⑶本发明方法不需要进行大量的预训练，能够对图片直接进行处理，为实际的临床应用提供了可能。

附图说明

图1为有丝分裂细胞的检测流程示意图，其中虚线矩形框为本发明的创新部分。

图2是本发明提出的算法对细胞有丝分裂检测的主要流程示意图。

具体实施方式

面结合附图对本发明的技术方案进行详细说明：

本发明公开的细胞有丝分裂检测方法，如图2所示，包括以下步骤：

步骤1、本发明所用的图片由10波段的多光谱显微镜扫描而成，每个谱波段均为可见光谱。此外，对于每个谱波段，我们用17个不同的聚焦平面来数字化图片，任意两个相邻的平面相差500nm。对每个波段的每幅图片而言，分别用一个列向量来表示，该列向量的长度与原图片矩阵大小相同。

原始图片是以矩阵的形式描述的，以矩阵的每一列为单位，展开成一个大的列向量，向量的长度与原始矩阵的大小相同。

步骤2、用步骤1产生的向量来构造一个新矩阵，矩阵的第一列对应步骤1中第一个列向量，以此按顺序排开；

将步骤1中产生的多个列向量按顺序排成一个大矩阵，即矩阵的第一列为第一个列向量，第二列为第二个列向量，依次类推，从而得到一个新的矩阵。

步骤3、对步骤2得到的新矩阵进行低秩分解，分解成大小相同的低秩矩阵和稀疏矩阵两部分；

本发明中所使用的低秩表示方法具体步骤如下：

本发明中，假设D＝[I₁,...,I_n]∈R^m×n表示n个连续的细胞有丝分裂病理图像。其中，I_j∈R^m是第j幅图像，每幅图像由包含m个像素的列向量组成，ij表示第j幅图像中的第i个像素。假设D中非有丝分裂部分（背景图像）是线性相关的，我们用和D相同大小的矩阵B∈R^m×n表示，则原矩阵D可以分解成低秩矩阵B和低秩表示中的异常值两部分，该异常值（前景）即为有丝分裂细胞。S∈{0,1}^m×n是表示前景支持的二进制矩阵，即：

P_S(X)表示矩阵X到由S支持的矩阵线性空间的正交投影：

是与P_S(X)互补的投影，满足

本发明还会用到矩阵的三个范数，分别是：表示l₁范数，表示Frobenius范数，||X||_*表示核范数，它是矩阵奇异值的和。

在细胞有丝分裂病理图像中，非有丝分裂部分（即背景部分）几乎是不变的。因此，对于图像序列D来说，相邻的两幅图像的背景部分灰度值几乎不变，即线性相关性很大。所以，背景矩阵B应该是一个低秩矩阵。假设给定一个实数阈值K，则矩阵B应该满足条件：

rank(B)≤K     (3)

其中，K越小，矩阵越相关。

由式（2）知，在已知矩阵D的情况下，如果给定前景支持矩阵S，即可求得前景矩阵。因此，对前景部分（即有丝分裂细胞）的建模可以转化成对前景支持矩阵S进行建模。对有丝分裂病理图像来说，有丝分裂细胞非常小，它由一些连续的像素点组成，所以可以通过马尔科夫随机场来建模。由矩阵D中所有像素点作为顶点υ来构造一个图模型G＝(υ,ε)，υ的大小为m×n，ε是边的集合，每条边都是用来连接时间或空间域上的相邻像素。由于有丝分裂细胞非常少，因此连接前景和背景的边应该很少。基于上述描述，前景支持矩阵S应该满足：

$\underset{\mathrm{ij}\∈\mathrm{\υ}}{\mathrm{\Σ}}{u}_{\mathrm{ij}}\left(S_{\mathrm{ij}}\right)+\underset{(\mathrm{ij},\mathrm{kl})\∈\mathrm{\ϵ}}{\mathrm{\Σ}}{\mathrm{\λ}}_{\mathrm{ij},\mathrm{kl}}|S_{\mathrm{ij}}-S_{\mathrm{kl}}|---\left(4\right)$

参数λ_ij＞0是S_ij＝1的惩罚因子，表明前景部分很少，是稀疏的。参数λ_ij,kl主要用来表示S_ij和S_kl之间的相关性，用来惩罚前景和背景的连接，表明前景部分是连续的。本文中为了简单起见，在所有的位置，我们都令λ_ij＝β，λ_ij,kl＝γ，其中β＞0和γ＞0是正常数。这就意味着我们不需要关于目标对象位置的额外的先验。所以，（4）式和（5）式可合并为：

$\mathrm{\β}\underset{\mathrm{ij}}{\mathrm{\Σ}}{S}_{\mathrm{ij}}+\mathrm{\γ}\underset{(\mathrm{ij},\mathrm{kl})\∈\mathrm{\ϵ}}{\mathrm{\Σ}}|S_{\mathrm{ij}}-S_{\mathrm{kl}}|---\left(6\right)$

在完成对前景支持向量S和背景模型B建模后，我们便可以得到信号模型D。当S_ij＝0时，假设D_ij＝B_ij+ξ_ij，其中，ξ_ij为独立同分布的高斯噪声，即D_ij～N(B_ij,σ²)，σ²是高斯噪声的方差。所以，我们可以通过最小二乘法来拟合D_ij。当S_ij＝1时，背景完全由前景遮挡，因此，D_ij等于前景的灰度值。

综合上述模型，我们将通过最小化以下能量函数来估计B和S：

$\begin{array}{c}\underset{B,S_{\mathrm{ij}}\∈\left\{\mathrm{0,1}\right\}}{\min }\frac{1}{2}\underset{\mathrm{ij}:S_{\mathrm{ij}}=0}{\mathrm{\Σ}}{(D_{\mathrm{ij}}-B_{\mathrm{ij}})}^2+\mathrm{\β}\underset{\mathrm{ij}}{\mathrm{\Σ}}{S}_{\mathrm{ij}}+\mathrm{\γ}\underset{(\mathrm{ij},\mathrm{kl})\∈\mathrm{\ϵ}}{\mathrm{\Σ}}|S_{\mathrm{ij}}-S_{\mathrm{kl}}|,\\ s.t.\mathrm{rank}\left(B\right)\≤K.\end{array}---\left(7\right)$

为了便于处理最小化能量函数，我们用核范数来代替矩阵B的秩操作符，核范数已经被证明是秩操作的有效的凸代替，而且还能避免过拟合。则（7）式可以改写为矩阵形式：

$\underset{B,S_{\mathrm{ij}}\∈\left\{\mathrm{0,1}\right\}}{\min }\frac{1}{2}{\left|\right|P_{S^{\⊥}}(D-B)\left|\right|}_F^2+\mathrm{\α}{\left|\right|B\left|\right|}_{*}+\mathrm{\β}{\left|\right|S\left|\right|}_1+\mathrm{\γ}{\left|\right|\mathrm{Avec}\left(S\right)\left|\right|}_1---\left(8\right)$

其中A是图G的节点-边的关联矩阵，α是和K相关联的大于零的参数，用来约束背景模型的复杂度。

步骤4、将步骤3得到的稀疏矩阵的每一个列向量都分别转化成相应的矩阵表示；

具体实施时，若假设原始图片的大小为M×N，则稀疏矩阵的每一列向量长度应为M×N，将第一列向量第1到M-1个元素看作矩阵第一列，第一列向量第M到2M-1个元素看作矩阵第二列，以此类推，便可将稀疏矩阵的第一列转化成一个M×N矩阵来表示；同样的方法对稀疏矩阵的第二列到最后一列都进行操作，可以得到多个不同的矩阵。

步骤5、对步骤4中得到的所有矩阵求和，得到一个新矩阵；

具体实施时，结果矩阵每个元素的值应该等于求和矩阵相同位置各个元素值的和。例如，A、B矩阵求和得到C，则C_ji＝A_ji+B_ji，ij表示第i行第j个元素。

步骤6、循环执行步骤1至步骤5直到所有波段的数据都处理完，得到多个大小相同的新矩阵；步骤1至步骤5是一个反复迭代过程，对每个波段的数据都进行相同的处理，直到所有波段的数据都处理完。

步骤7、对步骤6中得到的多个矩阵进行与操作，得到最终的检测结果，其中，非零部分即为有丝分裂部分；对大小相同的几个矩阵进行与操作时，得到的结果是一个大小相同的0、1矩阵，当且仅当所有矩阵同一个位置的值为非零元素时，结果矩阵在该位置的值才为1，否则为0。

为了便于公众理解本发明技术方案，下面给出一个具体实施例。

本实施例将本发明所提供的技术方案应用在由ICPR2012所举办的有丝分裂检测竞赛所提供的数据中。数据由5个乳腺癌组织切片组成。病理学家从每个切片中选取10个高能量区域，每个高能量区域的大小为512×512μm²（也就是0.262mm²），正好等于直径为0.58mm的显微镜所覆盖的区域。我们通过10波段的多谱显微镜M对这些切片进行扫描，所有波段均为可见光。此外，对于每个波段，通过17个不同的聚焦平面来数字化图片（17层的Z-stack）。由于多谱显微镜的分辨率为0.185μm每像素，因此每个高能量区域数字化后的图片大小为2767×2767。图片通过安置在显微镜顶端的照相机拍摄，而照相机只能生成1360×1360大小的图片，因此我们需要用4个子图来覆盖每个高能量区域，每个子图占总图的1/4。

在本实施例中，将原有的5个切片数据分成训练集和测试集两部分。训练集由3个切片的数据构成（供做对比的5种基于模式识别的方法使用），测试集由2个切片的数据构成。

本实施例的有丝分裂检测过程具体如下：

1、图片矩阵向量化：

对每个波段的的图片而言，将其表示成长度为1360×1360的列向量。每个波段由17幅不同的图片组成，因此可以生成17个长度均为1360×1360的列向量。

2、新矩阵的构造：

用步骤1产生的向量来构造一个新矩阵，矩阵的列数为17，行数为1360×1360。矩阵的第一列对应步骤1中第一个列向量，第二列对应第二个列向量，以此按顺序排开；

3、矩阵低秩分解：

对步骤2得到的新矩阵进行低秩分解，分解成大小相同的低秩矩阵和稀疏矩阵两部分；

4、稀疏矩阵的每一列构造新矩阵；：

步骤3得到的稀疏矩阵的行数为1360×1360，列数为17。将17个长度为1360×1360的列向量都表示成1360×1360矩阵的形式。

5、矩阵求和

对步骤4得到的17个大小相同的矩阵求和，得到一个新的大小为1360×1360的矩阵。

6、低秩分解过程迭代：

循环执行步骤1至步骤5直到10个波段的数据都处理完，得到10个大小均为1360×1360的新矩阵；

7、矩阵与操作：

对步骤6中得到的10个矩阵进行与操作，得到最终的检测结果，其中，非零部分即为有丝分裂部分；

为了验证本发明方法的效果，我们将和已有的基于模式识别的方法进行比较。由于本发明方法只用到了图像的灰度特征，为了比较的公平性，基于模式识别的方法也将只用灰度特征。先用10个波段的灰度特征来训练不同的分类器：决策树（Funcational Tree，FT），神经网络（Multilayer Perceptron，MP），贝叶斯（Bayesian），线性支持向量机（Linear SVM，L-SVM）和非线性支持向量机（NL-SVM）。然后再用训练好的分类器对2个测试切片数据进行测试。实验结果的度量指标主要有3个，分别为：

$\mathrm{Recall}\left(\mathrm{sensitivity}\right)=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}---\left(9\right)$

$\mathrm{Precision}\left(\mathrm{positivepredictivevalue}\right)=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}---\left(10\right)$

$F-\mathrm{measure}=2\×\frac{\mathrm{precision}\×\mathrm{recall}}{\mathrm{precision}\×\mathrm{recall}}---\left(11\right)$

其中，TP表示自动检测结果中，与手动标注结果相同的有丝分裂细胞的个数；FN表示手动标注的有丝分裂细胞中，没有被自动检测到的有丝分裂细胞的个数。最终结果如表1所示，加粗部分为同一种度量方法中精度最高的结果。

从表1可知，本发明提出的方法的检测率在Recall度量指标下的精度达到0.56，高于其它五种方法，比最接近的方法（MP、Bayesian classifier）高了8个百分点，这说明本文方法检测结果中，和手动标注结果相同的有丝分裂细胞个数占所有手动标注细胞个数的比例比较高。在precision的度量指标下，本文方法的精度为0.63，略低于L-SVM Classifier的0.68，与MP Classifier和FT Classifier相同。主要原因有二：一方面说明本文方法自动检测出来的所有结果中，和手动标注结果相同的个数所占比例相比L-SVM Classifier的方法还有待提高；另一方面，我们发现L-SVM Classifier方法的recall指标极低，只有0.25，说明该方法检测结果中，和手动标注结果相同的个数占所有手动标注结果的比例很低，而该方法的precision又较高，只能说明该方法自动检测的结果数量很少。在综合评价precision和recall的F-measure指标中，本文方法的精度是最高的，达到0.59，充分说明了本文方法的有效性。

表1检测结果比较

上面结合附图对本发明的实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下做出各种变化。

Claims (4)

1.一种细胞有丝分裂检测方法，其特征在于，具体包括以下步骤：

步骤1、用多光谱显微镜对有丝分裂中的细胞进行不同波段的扫描，并用不同的聚焦平面来数字化扫描得到的图片，形成图片矩阵；将扫描得到每幅图片所对应的图片矩阵分别用一个列向量来表示，该列向量的长度与原图片矩阵大小相同；

步骤2、用步骤1产生的列向量构造第一类矩阵，所述第一类矩阵的每一列分别按顺序对应步骤1中产生的每一个列向量；

步骤3、对步骤2得到的第一类矩阵进行低秩分解，分解成大小相同的低秩矩阵和稀疏矩阵；

步骤4、将步骤3得到的稀疏矩阵的每一个列向量都分别转化成相应的第二类矩阵；

步骤5、对步骤4中得到的所有第二类矩阵求和，得到第三类矩阵；

步骤6、循环执行步骤1至步骤5直到完成所有波段的数据处理，得到复数个大小相同的第三类矩阵；

步骤7、对步骤6中得到的复数个第三类矩阵进行与计算，得到最终的检测结果，其中的非零部分即为有丝分裂部分。

2.如权利要求1所述一种细胞有丝分裂检测方法，其特征在于：步骤1中将图片矩阵转化成列向量时，像素点的个数保持不变，大小为M×N的图片应转化成长度为M×N的列向量，其中M、N分别为自然数。

3.如权利要求1所述一种细胞有丝分裂检测方法，其特征在于：步骤3中对矩阵进行低秩分解，所得到的低秩矩阵和稀疏矩阵其大小均与低秩分解前的原第一类矩阵的大小一致。

4.如权利要求1所述一种细胞有丝分裂检测方法，其特征在于：步骤3中所述低秩分解的具体计算方法如下：

D＝[I₁,...,I_n]∈R^m×n表示n个连续的细胞有丝分裂病理图像，其中j∈n，I_j∈R^m表示第j幅图像，每幅图像由包含m个像素的列向量组成；矩阵D分解成低秩矩阵B∈R^m×n和低秩表示中的异常值S∈{0,1}^m×n两部分，B为背景图像，代表非有丝分裂部分，S为前景图像，代表有丝分裂细胞；

其中，ij表示第j幅图像中的第i个像素；

P_S(X)表示矩阵X到由S支持的矩阵线性空间的正交投影：

是与P_S(X)互补的投影，满足

给定一个实数阈值K，则矩阵B应该满足条件：

rank(B)≤K     (3)

其中，K越小，矩阵越相关；

通过马尔科夫随机场来建模：由矩阵D中所有像素点作为顶点υ构造一个图模型G＝(υ,ε)，υ的大小为m×n，ε是边的集合，每条边都是用来连接时间或空间域上的相邻像素，前景图像S应该满足：

$\underset{\mathrm{ij}\∈\mathrm{\υ}}{\mathrm{\Σ}}{u}_{\mathrm{ij}}\left(S_{\mathrm{ij}}\right)+\underset{(\mathrm{ij},\mathrm{kl})\∈\mathrm{\ϵ}}{\mathrm{\Σ}}{\mathrm{\λ}}_{\mathrm{ij},\mathrm{kl}}|S_{\mathrm{ij}}-S_{\mathrm{kl}}|---\left(4\right)$

参数λ_ij＞0是S_ij＝1的惩罚因子，表明前景部分是稀疏的；参数λ_ij,kl表示S_ij和S_kl之间的相关性，表明前景部分是连续的；在所有的位置都令λ_ij＝β，λ_ij,kl＝γ，其中β＞0和γ＞0是正常数；

将（4）式和（5）式合并为：

$\mathrm{\β}\underset{\mathrm{ij}}{\mathrm{\Σ}}{S}_{\mathrm{ij}}+\mathrm{\γ}\underset{(\mathrm{ij},\mathrm{kl})\∈\mathrm{\ϵ}}{\mathrm{\Σ}}|S_{\mathrm{ij}}-S_{\mathrm{kl}}|---\left(6\right)$

在完成对前景图像S和背景图像B建模后，得到信号模型D。；

当S_ij＝0时，假设D_ij＝B_ij+ξ_ij，其中，ξ_ij为独立同分布的高斯噪声，即D_ij～N(B_ij,σ²)，σ²是高斯噪声的方差；

通过最小二乘法来拟合D_ij，当S_ij＝1时背景完全由前景遮挡，D_ij等于前景的灰度值；

通过最小化以下能量函数来估计B和S：

$\begin{array}{c}\underset{B,S_{\mathrm{ij}}\∈\left\{\mathrm{0,1}\right\}}{\min }\frac{1}{2}\underset{\mathrm{ij}:S_{\mathrm{ij}}=0}{\mathrm{\Σ}}{(D_{\mathrm{ij}}-B_{\mathrm{ij}})}^2+\mathrm{\β}\underset{\mathrm{ij}}{\mathrm{\Σ}}{S}_{\mathrm{ij}}+\mathrm{\γ}\underset{(\mathrm{ij},\mathrm{kl})\∈\mathrm{\ϵ}}{\mathrm{\Σ}}|S_{\mathrm{ij}}-S_{\mathrm{kl}}|,\\ s.t.\mathrm{rank}\left(B\right)\≤K.\end{array}---\left(7\right)$

结合矩阵的三个范数：

表示l₁范数、表示Frobenius范数、||X||_*表示核范数，它是矩阵奇异值的和；用核范数来代替矩阵B的秩操作符，将（7）式改写为矩阵形式：

$\underset{B,S_{\mathrm{ij}}\∈\left\{\mathrm{0,1}\right\}}{\min }\frac{1}{2}{\left|\right|P_{S^{\⊥}}(D-B)\left|\right|}_F^2+\mathrm{\α}{\left|\right|B\left|\right|}_{*}+\mathrm{\β}{\left|\right|S\left|\right|}_1+\mathrm{\γ}{\left|\right|\mathrm{Avec}\left(S\right)\left|\right|}_1---\left(8\right)$

其中，A是图G的节点-边的关联矩阵，α是和K相关联的大于零的参数，用来约束背景模型的复杂度。

Images