CN105719294B - 一种乳腺癌病理学图像有丝分裂核自动分割方法 - Google Patents

Abstract

本发明公开了一种乳腺癌病理学图像有丝分裂核自动分割方法，该方法首先根据BR图的不同尺度LoG滤波的最大响应初略估计出前景像素和背景像素。在此基础上对当前的前景像素和背景像素分别训练一个高斯混合模型GMM以描述前景像素和背景像素的颜色分布；并建立s‑t图G＝{V,E}，通过计算图G的最大流实现了对前景和背景的新的估计，然后重复GMM训练和Graph‑Cut分割，直到收敛或满足规定的迭代次数。最后对分割的结果进行形态学的开操作，完成最终的分割。该方法在分割时不仅考虑了单个像素颜色特征属于前景或背景的概率，还考虑了相邻像素间的关系，保证了分割结果的平滑性，提高了精度。

Description

一种乳腺癌病理学图像有丝分裂核自动分割方法

技术领域

本发明属于图像处理技术领域，特别涉及一种乳腺癌病理学图像有丝分裂核自动分割方法。

背景技术

据国家卫生部统计，乳腺癌在我国女性患癌中占据了很大的比例，成为了危害女性的生命安全的最大杀手之一，其发病年龄下从20岁到70、80岁，且发病率随着年龄的增长而呈上升势态。目前我国大多数乳腺癌患者集中在45～55这个年龄段，发病率占全身各种恶性肿瘤的7-10％，在女性中仅次于子宫癌的发病率。而目前关于乳腺癌的病因尚未完全清楚，其早期发现和早期诊断则显得至关重要。

在乳腺癌严重程度评估中，有丝分裂计数是一个很重要的指标。但是，现有的有丝分裂计数全靠人工进行，主要是依靠经验丰富的专家在显微镜下进行人工观察，非常繁琐，效率低，而且不同的专家往往有不同的结果，可重复性低。因此，采用图像处理技术设计合适的算法，利用计算机对显微镜图像进行分析，自动完成有丝分裂核的分割与识别，实现有丝分裂的自动检测与计数，对当前乳腺癌的诊断和筛查均具有极其重要的意义。

经过标准H&E染色后的切片，发生了有丝分裂的细胞核和未发生有丝分裂的细胞核在外观上看起来相差不大，且都拥有不规则的外表，有丝分裂核表现出来的多样性使得对其进行自动检测与计数的研究极具有挑战性，已成为当前的研究热点。

现有方法的标准流程是首先采用图像分割方法提取大量的候选区域，然后提取候选区域的各种形状、纹理等特征并训练分类器完成有丝分裂的自动识别。大量的诸如阈值分割、分水岭分割、形态学方法、基于LoG的blob检测等图像分割方法被用来分割细胞核，然而这些算法过于简单，在分割出有丝分裂核(正样本)的同时往往会分割出大量的非有丝分裂部分(负样本)，其平均正负样本的比例达到1：100，而正负样本分布不平衡给后续的分类器训练带来了很大的困难。

发明内容

本发明提供了一种乳腺癌病理学图像有丝分裂核自动分割方法，其目的在于克服上述现有技术中的不足，通过构建前景和背景的颜色概率模型，同时考虑像素之间的空间关系，采用Graph Cuts的方法完成前景、背景的分割的迭代方法，获得准确的分割结果，且大大减少了后续分类器训练的难度。

一种乳腺癌病理学图像有丝分裂核自动分割方法，包括以下步骤：

步骤1：对待分割的原始图像进行高斯平滑处理，获得平滑图像；

步骤2：从平滑图像中提取r(x,y),g(x,y),b(x,y)通道，按以下公式计算获得BR图像；

步骤3：利用不同尺度的LoG滤波器对BR图像进行卷积滤波，以每个像素点的最大响应作为滤波图像；

步骤4：采用Otsu算法计算滤波图像的二值化阈值T，利用二值化阈值T对滤波图像进行前景和背景分割，获得标记图M(x,y)；

所述标记图用于区分背景区域和前景区域；

步骤5：以标记图中像素值为0的像素对应在待分割的原始图中的像素作为背景像素，待分割的原始图中的其余像素为前景像素，利用背景像素集与前景像素集分别训练前景颜色模型p(z；θ_f)和背景颜色模型p(z；θ_b)；

步骤6：根据待分割的原始图像I(x,y)、标记图M(x,y)、前景颜色高斯模型和背景颜色高斯模型构建s-t图G＝{V,E}，采用最大流的方法计算s-t图G的最小割，并用最小割对s-t图G进行分割，更新标记图M(x,y)；

步骤7：返回步骤5，重新训练前景和背景颜色模型，更新s-t图G中边的权重，计算最大流对s-t图G进行分割，进而更新M(x,y)直到满足迭代中止条件，进入步骤8；

步骤8：对最新的标记图中每个像素与1进行“与”位运算，得到二值分割结果图。

对所述步骤4中的标记图M(x,y)按以下公式进行计算：

所述前景颜色模型p(z；θ_f)和背景颜色模型p(z；θ_b)按以下过程获得：

首先，采用K均值算法分别对前景像素和背景像素进行聚类，获得前景颜色模型的聚类结果数量为K_f，背景颜色模型的聚类结果数量为K_b；

其次，用聚类结果来初始化训练前景高斯模型和背景高斯模型：

前景高斯模型：

背景高斯模型：

其中，表示第k个高斯分布模型，π_k、μ_k和∑_k分别为该高斯分布模型的先验概率、均值和协方差矩阵，N表示所有用于训练模型的像素的个数，N_k表示属于第k类像素的个数，z^T表示矢量z的转置；

所述构建的s-t图G＝{V,E}，其中，V为节点集合，包含待分割的原始图像中的所有像素点以及一个前景标识点和背景标识点，共W×H+2个节点；

E为边的集合，包含相邻边和关系边；

其中，W和H分别表示待分割的原始图像的宽度和高度；

其中，所述相邻边是指每个像素点与其四邻域像素的连接边，共2×[4×W×H-3×(W+H)+2]条边；所述关系边是指每个像素分别到前景标识点和背景标识点的连接边，共2×W×H条边；

相邻边的权重按公式计算获得；

其中，z_i,z_j∈R³分别表示相邻像素i和j的颜色RGB矢量，参数γ＝40，参数β则由下式计算

关系边的权重按以下公式计算获得：

ω_b(x,y)＝-log(p(z；θ_f)),

其中，ω_f(x,y)和ω_b(x,y)分别表示像素(x,y)到前景标识点和背景标识点的权重，参数λ＝9γ。

所述相邻边用于描述分割的光滑性，所述关系边用于描述每个像素到前景标识点和背景标识点的关系；

所述步骤6中用最小割对图G进行分割，更新标记图M(x,y)是指将像素分割成前景点和背景点两类，然后按照以下规则更新标记图：

M(x,y)＝2,if节点(x,y)属于背景点。

采用9×9大小的结构元素对步骤8获得的分割结果进行形态学开操作，得到最终的分割结果图。

所述步骤2中不同尺度的LoG滤波器的尺度取值依次为2⁰，2¹，2²。

有益效果

本发明提供了一种乳腺癌病理学图像有丝分裂核自动分割方法，首先在对图像进行高斯平滑后计算其BR图以突出细胞与背景的对比；再采用一组不同尺度的LoG滤波器对BR图进行滤波，然后计算在每个位置上的各个尺度上的最大响应作为滤波结果；并自适应的计算一个较高的阈值对滤波后的图像进行二值化操作，估计出可能性大的前景像素，其他像素则置为背景像素。在此基础上一方面对当前的前景像素和背景像素的RGB矢量分别训练一个高斯混合模型GMM以描述前景像素和背景像素的颜色分布；另一方面，以RGB图像建立一张s-t图G＝{V,E}，其中节点集合V中包含所有像素点和前景标识点和背景标识点，边的集合E中不仅包含描述每个像素与其四邻域像素之间关系的部分，还包含每个像素到源节点和宿节点关系的部分，并分别计算各条边的权重。采用增广路径的方法计算s-t图G＝{V,E}的最大流(max flow)，从而间接求出该图的最小割(minCut)，完成对整幅图像的分割，实现了对前景区域和背景区域的新的估计，然后重复训练GMM和基于Graph-Cut的分割，直到收敛或满足规定的迭代次数。最后对分割的结果进行形态学的开操作，以过滤小的前景区域，完成最终的有丝分类核的分割。该方法在分割是不仅考虑了单个像素颜色特征属于前景或背景的概率，还考虑了相邻像素间的关系，保证了分割结果的平滑性，在保证高的召回率(成功找到有丝分裂核)的同时能有效降低非有丝分裂核的个数，提高了精度。

附图说明

图1是本发明的流程图；

图2是原始的RGB图像I(x,y)；

图3是根据I(x,y)计算得到的blue ratio图像BR(x,y)；

图4是根据图像BR(x,y)计算得到的最大LoG响应图像L(x,y)；

图5是第一次迭代的结果；

图6是第二次迭代的结果；

图7是第三次迭代的结果；

图8是第四次迭代的结果；

图9是第五次迭代的结果；

图10是经形态学开操作后的结果；

图11是直接对图像L(x,y)进行二值化的分割结果；

图12是对应原始图像I(x,y)的人工分割结果，作为分割的金标准。

具体实施方式

下面将结合附图和实施例对本发明做进一步的说明。

如图1所示，本发明的方法主要包括以下步骤：对原始图像I(x,y)进行预处理，计算其BR图像；用一组不同尺度的LoG滤波器对BR图像进行滤波，计算最大响应图像L(x,y)；根据otsu算法求最大响应图像的全局阈值，进而计算出标记图M(x,y)；根据标记图M(x,y)收集图像I(x,y)上像素，分别训练前景高斯混合颜色模型p(x；θ_f)和背景高斯混合颜色模型p(x；θ_b)；根据图像I(x,y)、标记图M(x,y)和高斯混合模型构建s-t图G＝{V,E}；计算s-t图G的最大流对图G进行分割，并更新标记图M(x,y)；重新训练前景和背景颜色模型，更新s-t图G中边的权重，计算最大流对s-t图G进行分割，进而更新M(x,y)直到满足迭代中止条件；对M(x,y)的每个像素与1进行“与”位操作，并对位操作后的结果进行形态学开操作，得到最终的分割结果作为候选区域。

步骤1：采用高斯滤波器对原始的彩色图像I(x,y)＝[r(x,y)；g(x,y)；b(x,y)](如图2所示)的每个分量进行高斯平滑

即r(x,y):＝h(x,y)*r(x,y),g(x,y):＝h(x,y)*g(x,y),b(x,y):＝h(x,y)*b(x,y)，

这里*表示卷积操作，实际中采用7×7空间模板的方式近似高斯函数，其中σ＝1.0；

分别提取平滑后图像的r,g,b通道，由如下公式计算BR图(如图3所示)；

步骤2：采用一组具有不同尺度的LoG滤波器h_σ(x,y对BR图进行卷积滤波，L_σ(x,y)＝h_σ(x,y)*BR(x,y)；

其中，LoG滤波器h_σ(x,y)的数学表达式为

在滤波时同样采用模板进行近似，模板大小为round(6×σ+1)×round(6×σ+1)，这里round()表示取整运算。然后计算在每个位置(x,y)的最大响应作为最终的滤波结果，即L(x,y)＝max_σL_σ(x,y)，结果如图3所示。

步骤3：采用0-255中所有的数作为阈值对L(x,y)进行二值化分割，并计算每次分割结果像素的类内方差，取使得类内方差最小的值作为全局最优阈值T。然后根据下式计算图像的标记图

步骤4：从图像I(x,y)中收集对应在标记图M(x,y)中为0或2的所有像素作为背景像素，其余像素作为前景像素。分别对这些前景像素和背景像素进行K均值聚类，用得到聚类结果去训练对应的高斯混合模型。具体而言包括：

(1)随机选择K个像素的颜色作为均值的初始估计μ₁,…,μ_K；

(2)计算每个像素颜色矢量z到当前这k个均值的距离d_i＝||z-μ_i||，然后根据最近邻分类准则将其归入第c类，其中c＝arg min_i d_i；

(3)根据分类的结果重新计算每类的均值μ₁,…,μ_K；

(4)重复上述步骤直到计算得到的均值不再发生变化为止。

上述K均值聚类将所有像素分成K类可以用来训练高斯混合模型 具体而言可根据如下公式完成高斯混合模型的参数的估计：

这里N表示所有用于训练模型的像素的个数，N_k表示属于第k类像素的个数，z^T表示矢量z的转置。

对前景像素和背景像素按上述步骤分别训练，可分别得到前景高斯混合模型 和背景高斯混合模型对应的模型参数分别用θ_f和θ_b表示，K_f＝3,K_b＝5。

步骤5：对图像I(x,y)建立一张s-t图G＝{V,E}，其中节点集合V中包含所有像素点和一个源节点和宿节点，共W×H+2个节点，这里W和H分别表示图像的宽度和高度。边的集合E中包含两大类：一类是描述分割的光滑性，由每个像素与其四邻域的像素之间的表示，共2×[4×W×H-3×(W+H)+2]条边；另一类是描述每个像素到源节点和宿节点的关系，共2×W×H条边。同时计算每条边的权重，其中第一类边的权重由下式计算

这里z_i，z_j分别表示相邻像素i和j的颜色矢量，参数γ＝40，参数β则由下式计算

这里C为所有相邻像素对(i,j)的集合。

第二类边的权重由下式计算

ω_b(x,y)＝-log(p(z；θ_f)),

这里ω_f(x,y)和ω_b(x,y)分别表示像素(x,y)到源节点(前景)和宿节点(背景)的权重，参数λ＝9γ。

步骤6：计算s-t图G＝{V,E}的最小割，将图中的所有像素节点分成两类(即属于源节点或宿节点)。根据最大流/最小割定理，图的最小割等价于图的最大流。将s-t图G中各条边的权重视为该边的初始容量，可通过如下迭代算法计算s-t图G的最大流：

(1)每次迭代搜索任意一条从源节点出发到宿节点的增广路径，其流量(flow)恒为正；

(2)路径上所有的边的容量减去对应的流，继续迭代，直至无法找到这样的一条路径为止，此时余量容量为0的那些边将节点分成两部分，即属于源节点或宿节点。

完成对像素节点的分类后，可根据下式重新更新标记图M(x,y)：

M(x,y)＝2,if节点(x,y)属于宿节点。

步骤7：根据更新后的标记图M(x,y)重复步骤4-步骤6，直到收敛或迭代次数超过预定义的次数。

步骤8：对标记图M(x,y)的每个像素与1进行“与”位运算，得到二值分割结果图；图5、图6、图7、图8、图9分别显示经过1次、2次、3次、4次、5次迭代的分割结果。

步骤9：对二值分割结果进行形态学开操作，去除小的区域。最终的分割结果如图10所示。

图11给出了当前主流的采用基于LoG响应图的进行分割的对比结果，图12则给出了人工分割的结果图像(金标准)，可以看出本发明所提出的方法在保证不丢失真实有丝分裂细胞(即召回率＝100％)的同时，大大减少了候选区域的数量，提高了精度。

Claims (5)

1.一种乳腺癌病理学图像有丝分裂核自动分割方法，其特征在于，包括以下步骤：

步骤1：对待分割的原始图像进行高斯平滑处理，获得平滑图像；

步骤2：从平滑图像中提取r(x,y),g(x,y),b(x,y)通道，按以下公式计算获得BR图像；

步骤3：利用不同尺度的LoG滤波器对BR图像进行卷积滤波，以每个像素点的最大响应作为滤波图像L(x,y)；

步骤4：采用Otsu算法计算滤波图像的二值化阈值T，利用二值化阈值T对滤波图像进行前景和背景分割，获得标记图M(x,y)；

步骤5：以标记图中值为0或2的像素对应在待分割的原始图中的像素作为背景像素，待分割的原始图中的其余像素为前景像素，利用背景像素集与前景像素集分别训练前景颜色模型p(z；θ_f)和背景颜色模型p(z；θ_b)；

步骤6：根据待分割的原始图像I(x,y)、标记图M(x,y)、前景颜色高斯模型和背景颜色高斯模型构建s-t图G＝{V,E}，采用最大流的方法计算s-t图G的最小割，并用最小割对s-t图G进行分割，更新标记图M(x,y)；

步骤7：返回步骤5，重新训练前景和背景颜色模型，更新s-t图G中边的权重，计算最大流对s-t图G进行分割，进而更新M(x,y)直到满足迭代中止条件，进入步骤8；

步骤8：对最新的标记图中每个像素与1进行“与”位运算，得到二值分割结果图；

对所述步骤4中的标记图M(x,y)按以下公式进行计算：

所述构建的s-t图G＝{V,E}，其中，V为节点集合，包含待分割的原始图像中的所有像素点以及一个前景标识点和一个背景标识点，共W×H+2个节点；

E为边的集合，包含相邻边和关系边；

其中，W和H分别表示待分割的原始图像的宽度和高度；

其中，所述相邻边是指每个像素点与其四邻域像素的连接边，共2×[4×W×H-3×(W+H)+2]条边；所述关系边是指每个像素分别到前景标识点和背景标识点的连接边，共2×W×H条边；

相邻边的权重按公式计算获得；

其中，z_i,z_j∈R³分别表示相邻像素i和j的颜色RGB矢量，参数γ＝40，参数β则由下式计算：

关系边的权重按以下公式计算获得：

ω_b(x,y)＝-log(p(z；θ_f)),

其中，ω_f(x,y)和ω_b(x,y)分别表示像素(x,y)到前景标识点和背景标识点的权重，参数λ＝9γ。

2.根据权利要求1所述的方法，其特征在于，所述前景颜色模型p(z；θ_f)和背景颜色模型p(z；θ_b)按以下过程获得：

首先，采用K均值算法分别对前景像素和背景像素进行聚类，获得前景颜色模型的聚类结果数量为K_f，背景颜色模型的聚类结果数量为K_b；

其次，用聚类结果来初始化训练前景高斯模型和背景高斯模型：

前景高斯模型：

背景高斯模型：

其中，表示第k个高斯分布模型，π_k、μ_k和∑_k分别为该高斯分布模型的先验概率、均值和协方差矩阵，N表示所有用于训练模型的像素的个数，N_k表示属于第k类像素的个数，z^T表示矢量z的转置；

3.根据权利要求2所述的方法，其特征在于，所述步骤6中用最小割对图G进行分割，更新标记图M(x,y)是指将像素分割成前景点和背景点两类，然后按照以下规则更新标记图：

M(x,y)＝2,if节点(x,y)属于背景点。

4.根据权利要求1-3任一项所述的方法，其特征在于，采用9×9大小的结构元素对步骤8获得的分割结果进行形态学开操作，得到最终的分割结果图。

5.根据权利要求1-3任一项所述的方法，其特征在于，所述步骤2中不同尺度的LoG滤波器的尺度取值依次为2⁰，2¹，2²。