CN1735907A - 图像分析 - Google Patents

Abstract

一种用于数字图像的自动分析的方法，尤其用于基于在组织载玻片的数字图像中识别的细管组织和上皮细胞的相对比例来评估乳房组织中癌症的存在和严重程度。该方法包括下列步骤：使用选择像素的邻域中强度的局部平均和局部标准偏差的特性从图像中的一些或所有像素生成特性同现矩阵(PCM)，以及基于PCM的分析通过将选择像素标记为属于特定类型来分割图像。这样，表示图像中的上皮细胞的相对较暗和基本上有纹理的区域能够与更亮和更均匀的背景区域相区分。其他步骤包括根据强度、形状和尺寸标准识别图像中表示输送管细胞的像素组，把这些像素扩大一个尺寸到标记为上皮细胞的周围组以对应整个细管结构，并根据在这样扩大后输送管像素的数量和输送管与上皮像素的总数量的比值计算一个尺度。该方法的其他应用可以包括包含特定类型的晶体形成的矿物标本的分析。

Description

图像分析

技术领域

本发明涉及数字图像的自动分析。本发明尤其涉及在组织或细胞标本的数字图像中的不同细胞分组的自动识别，并尤其用于根据在组织切片的数字图像中识别的细管结构(tubule formation)和上皮细胞的相对比例来评估乳房组织中癌症的存在和严重程度，并且在此主要描述的本发明正是关于这一方面。然而，本发明也可以应用于呈现适合于以类似方式识别的图像成分的各种其他种类结构的分析中，例如可以应用于包含特定类型的晶体形成的矿物标本的分析。

背景技术

每年有数以千计的妇女不必要地死于乳腺癌，如果足够早地检测到该癌症，则理论上生存的可能性很高。如果在一个样本中未察觉癌组织的存在，则到下一次进行检查时，癌症可能已经发展并大大降低了存活的机会。因此，在样本中检测癌组织的重要性无论如何强调也不为过。

典型的国家乳房筛选程序使用乳房X线照相术进行触摸不到的机体损害的早期检测。一旦检测到代表乳腺癌的机体损害，则取组织样本并由受过训练的组织病理学家进行检查以建立诊断和预后。这是一个耗费时间、劳动密集和代价昂贵的过程。执行这种检查的资格不容易获得并需要频繁地复查。检查本身需要由眼睛分析彩色图像，以观察者之间和观察者自身分析中值得考虑的变化为特征的高度主观的过程，即不同的组织病理学家和相同的组织病理学家在不同的时间对相同的样本可能出现观察的变化。例如，研究显示检查同样十个样本的两个不同的组织病理学家对这些样本中的三个会给出不同的鉴定，即有30％的误差。这一问题由一些样本的复杂度而加重，尤其在没有一个确定的结论的边缘情况下更是这样。如果没有足够可用的受过训练的人员，则这对于完成该分析的压力造成影响，从而可能会造成错误的评估和诊断的延迟。

这些问题意味着实际上在乳腺癌筛选的程度和效果上有局限，因此一些妇女不能被正确识别为具有该疾病，并且在一些情况下，这样的失败可能造成过早的死亡。相反，其他人被错误地诊断为有乳腺癌，并因此可能不必要地遭受损伤性的治疗。

发明内容

因此，本发明的目的是提供一种图像分析方法，它能够被实现为一种鲁棒(robust)、客观和划算的工具以辅助乳腺癌的诊断和预后，尽管如前所述本发明也可以应用于其他领域。

为了帮助理解这一目的，参照附图1，附图1是典型地出现在乳房组织的组织载玻片(slide)中的各种对象的简化表示。存在包括诸如在由上皮层2包围的用1表示的输送管(duct)的细管结构。输送管显示为各种形状的小的、明亮的区域，同时上皮细胞基本上显示为更有纹理且更暗。诸如用3标识的脂肪细胞显示具有与输送管1的同样的强度，但是通常基本上更大。与输送管1和脂肪细胞3同样强度的伸长区域也可能存在，比如用4表示的区域，并且这是组织中的撕裂(tear)或由收缩引起的裂缝的特征。载玻片的其余部分包括“背景”组织5，它一般显示为比输送管1、脂肪细胞3和撕裂/裂缝4更暗，但比上皮细胞2更亮并且纹理更均匀。健康的组织应当包含大量的细管结构，该细管结构包括通常具有两个上皮细胞边界的输送管。在癌组织中，细管往往断裂而上皮细胞增生，所以任何给定样本中这些结构之间的面积比能够用作癌症存在和严重程度的指示。更特别的是，组织病理学家按照惯例作出尺度M的主观评估，由下式给出：

$M=\frac{T}{D+E}...\left(1\right)$

其中T是由细管结构(输送管加两个上皮细胞的边界)覆盖的载玻片中的表面积，D是由输送管覆盖的表面积，以及E是由所有上皮细胞(包括T中的那些)覆盖的表面积，并且使用阈值将该尺度的值的他们的评估与癌症等级相联系，典型地如下所示：

尺度值          癌症等级

≥75％            1级

≥10％，＜75％    2级

＜10％            3级

表1.用于癌症严重程度的组织病理学阈值其中1级最轻，而3级最严重。

如果相同或类似尺度的客观评估通过一个图像分析的自动方法而获得，则尤其需要区分在一个包括上皮细胞的图像中的那些对象，因此一方面，本发明在于一种用于对包括像素阵列的数字图像进行自动分析的方法，该方法包括下列步骤：利用在所选择像素的邻域(neighbourhood)中强度的局部平均和局部标准偏差的特性从所述的一些或全部像素中生成特性同现矩阵(property co-occurrencematrix)(PCM)；以及作为分析所述PCM的结果通过将所选择的像素标记为属于特定分类来分割图像。

本发明也在于包括执行前述方法的装置的用于数字图像的自动分析的设备，以及在于一种包括在其上具有适于引起计算机执行上述方法的计算机程序代码装置的计算机可读介质的计算机程序产品和一种包括这样做的指令的计算机程序。

特性同现矩阵(PCM)例如在Electronics and CommunicationEngineering Journal的1993年第5卷第2期第71-83页(Co-occurrence Matrices for Image Analysis，JF Haddon和JF Boyce)进行了描述，并且是例如在IEEE Trans.Syst.，Man，Cybern.的1973年Vol SMC-3的第610-621页(Texture Features for ImageClassification，RM Haralick，K Shanmugan和I Dinstein)中描述的标准灰度级同现矩阵的扩展或概括。它们是多维的直方图，其中每个元素是所选择的特性同时出现的频率。通过使用图像像素的邻域中局部平均和局部标准偏差的特性生成一个PCM，因此这样的PCM的分析能够区分构成比如说相对较低局部平均和相对较高局部标准偏差(比如在本发明的这一方面的优选实施方式中表示上皮细胞的暗的、有纹理的区域)区域的像素和构成比如说相对较高局部平均和相对较低局部标准偏差(比如在本发明的这一方面的优选实施方式中表示“背景”组织的更亮、更均匀的区域)区域的像素，或者区分在本发明的其他应用中的那些特性的其他组合的区域的像素。

现在将参照附图借助于例子并在用于在乳房的潜在癌的组织载玻片的数字图像中的细管结构的基础上对癌症分级的自动系统的上下文中更详细地描述本发明的这些和其他方面。

附图说明

在附图中：

图1是乳房组织的组织载玻片中典型对象的简化表示，该组织可以按照本发明的优选实施例进行分析；

图2是在用于获得和分析数字化图像的优选实施例中设备的框图；

图3示出优选实施例中主要算法部分的布局和处理流程；

图4是一幅图像中像素的示意图；

图5示出用于在优选实施例中确定掩蔽(mask)点的过程；

图6(a)示出了直方图中两个峰值之间的一个谷，而图6(b)中没有谷；

图7是示出用于PCM轴的不同标记和边缘分布的图；

图8(a)示出原始直方图，而图8(b)示出在两侧上延伸的这一零点，图8(c)-(f)示出通过在2、4、8和16个二进制数(bin)上进行平均的逐渐增大标度的同一个直方图；对于每一种分辨率示出所有峰值和谷的位置，这是与特定峰值或谷相关联的轨迹长度(tracklength)；

图9示出一个PCM和使用EM算法拟合的两个高斯分布的位置和参数；所述分布的半径以两个标准偏差画出；

图10示出如何计算从PCM中的一个地址到拟合的分布之间的距离，从而能够确定一个像素的分割值；

图11示出与像素X相邻的8个像素(a，b，c，d，e，f，g，h)；

图12示出在优选方法的过程中在图像分割之后的几个斑点(blob)的例子；以及

图13是优选方法的过程中形象化的扩大(dilation)结果。

具体实施方式

通用系统结构

图2示出用于执行本发明的进行图像捕获和处理的基于典型计算机的系统。各切片从乳房组织样本中切割，放在载玻片上并且按照传统技术染色。一个病理学家在显微镜21中细看载玻片，选择按照将执行的分析看起来最有希望的区域，并利用数码相机22对它们进行拍摄。来自照相机22的图像被下载到个人计算机(PC)23，在那里图像按照以下描述进行存储和处理。在本发明人使用的系统中，显微镜提供10X光学放大，并且数字图像是横向1476像素乘以纵向1160像素。可以使用其他放大倍率和数字化尺寸而不损害以下更详细描述的算法，假如诸如细胞大小、扩大中的最大桥接间隙和形状标准之类的一些系统参数相应进行调整的话。显微镜和照相机可以由其他适当的设备如高分辨率平板扫描仪等来替换。自动设备能够替换显微镜/数字化设备组合。PC也可以用适当处理能力的任何通用计算机或专用硬件来替换。这里描述的技术能够应用于数字成像，而不管数据是怎样获得的。

处理概述

图3示出在本发明的优选实施例中的处理部分的概略图，并且这些将在随后的部分中更详细地单独论述。但是，简单地说，该处理如以下所述进行。

在初始数字化之后，如果需要，在第一步骤31将校正图像的彩色平衡和晕光(vignetting)。在步骤32，识别掩蔽点以基于强度值的局部标准偏差和局部平均从在步骤33生成的PCM中排除特定像素(它们既不是上皮细胞的部分，也不是“背景”组织)。在步骤34分析PCM分布以区分分别属于上皮细胞和“背景”组织的像素，并且在步骤35通过将像素标记为上皮或背景(来自步骤34)或掩蔽(来自步骤32)来分割图像。具有相同标记的相邻像素被分组为斑点，并进行过滤以在步骤36清洁图像(或者在掩蔽像素将在步骤32被分组并过滤的情况下)。因此在步骤37识别和再次标记具有代表输送管的形状和尺寸特征的标记为掩蔽的像素的斑点。接着在步骤38将所得到的输送管斑点扩大一个对应于两个上皮细胞的量来与标记为上皮的斑点相邻，从而它们现在覆盖一个对应于细管结构面积的一个面积。在步骤39，基于扩大后的输送管像素与输送管和上皮像素的总数的比率计算一个尺度，并且该尺度在步骤40被转换为癌症严重程度的指示和/或在步骤41用于训练分类器。

数字化和校准

图2所示的显微镜21、照相机22和计算机23的组合提供如图4所示物理尺寸为L₁×L₂并且分辨率为N×M像素单元的数字图像。因此，每个像素的物理尺寸为Δx×Δy，或者L₁/M×L₂/M。每个像素将具有在图像中的一个地址(x，y)，它通常用向量x表示。每个像素具有灰度级强度I或与它相关联的颜色元组。在优选实施例中和本发明开发中使用的设备中，每个像素具有与它相关联的红、绿和蓝强度(I_r，I_g，I_b，)并且是正方形。

在步骤31中，利用载玻片的清晰部分来拍摄一个校准图像。与显微镜相关联的照明的强度增强，直到在一个或多个波段(红、绿或蓝)中的几个像素的最大强度处于饱和或接近饱和。应当使尽可能少的像素饱和，但是一般大约在1％。任何透镜系统均引起诸如图像中的晕光之类的变化，并且这影响整个图像的强度的均匀性。如果这些变化严重的话，则它们在应用任何图像处理算法之前可能需要进行校正。在使用的设备中，晕光效应引起图像的中间和角落之间高达10％的强度变化。这可能影响算法的功效并且最好进行校正。

在优选实施例中，图像晕光通过将像素x的彩色分量缩放F倍进行粗略校正，其中

$F_k\left(x\right)=S\frac{I_k\left(x\right)C_k\left(x\right)}{\underset{x}{\max }{C}_k\left(x\right)}...\left(2\right)$

其中S是定标因子。在优选实施例中，S＝0.95。下标k指波段；红、绿、蓝或灰度或适当的成分。C_k是校准图像的第k个波段，并且函数max表示索引(index)x的最大值。在以上例子中，该处理打算使图像变得一致以便大大小于5％。但是，将进行分析的区域能够进一步被限制，如果想要可以通过在图像上强制加入一个椭圆的掩蔽，从而椭圆的轴是图像尺寸的95％。这进一步有助于确保图像强度在将被分析的区域中充分一致。

彩色平衡也应当一致并且在不同的数字化之间可复制。如果使用红、绿和蓝波段的组合，比如在一个灰度级图像中，则这可能成为关键。彩色平衡校正能够通过假定校准图像中的最大灰度级强度对应于峰值白色并且强制红、绿和蓝分量相等且为最大强度而获得。因此，任何像素的红、绿和蓝分量接着能够被线性定标。也可能使用认为应当是峰值白色的几个像素的平均：

$I_{\mathrm{gr}}\left(x\right)=\frac{I_r\left(x\right)+I_g\left(x\right)+I_b\left(x\right)}{3}...\left(3\right)$

$I_{r^{\′}}\left(x\right)=I_r\left(x\right)\frac{S_c}{I_{\mathrm{rw}}}$

$I_{g^{\′}}\left(x\right)=I_g\left(x\right)\frac{S_c}{I_{\mathrm{gw}}}...\left(4\right)$

$I_{b^{\′}}\left(x\right)=I_b\left(x\right)\frac{S_c}{I_{\mathrm{bw}}}$

其中，I_gr是从红、绿和蓝分量形成的灰度级图像。

I_rw、I_gw、I_bw是对应于应当是峰值白色的像素(或对应于所述像素的那些分量的平均值)的红、绿和蓝分量。

S_c是确定实际峰值白色强度的定标因子；这一般为255。

该处理假定峰值黑色对应于(0，0，0)并且不需要进行校正。对等式(4)的较小改变将使得这一点考虑进去。通过对将进行分析的数字化图像运用彩色和晕光校正，于是产生的图像具有足够的均匀性以适合于以下的分析。如果没有显著的强度变化就能够获得一幅图像，则不需要校准。

以下的处理步骤能够在单个波段图像上操作，该图像可以是红、绿或蓝分量或者它们的某一组合，比如从彩色分量的平均形成的灰度级图像。在实践中，已经发现当使用传统的染色时，红波段包含最多用于在图像的相关部分之间区分的信息。

生成特性同现矩阵

特性同现矩阵(PCM)是这一算法的基础。如前面表示的，它们是对于由Haralick等人定义的标准灰度级同现矩阵的扩展或概括。PCM是多维直方图，其中每个元素是特定特性在一幅图像中同现的频率。形式上，一个PCM S可以定义为：

$S(l_1,l_2,......,l_K)=\underset{x}{\∑}\mathrm{\δ}(l_1;P_1\left(x\right))\mathrm{\δ}(l_2;P_2\left(x\right))......\mathrm{\δ}(l_K;P_K\left(x\right))$

$=\underset{x}{\∑}\underset{k=1,K}{\mathrm{\Π}}\mathrm{\δ}(l_k;P_k\left(x\right))...\left(5\right)$

其中P_k是像素x的第k个特性

δ是Kronecker delta函数，从而

在优选实施例中，在步骤33生成的PCM是2维的，并且所用的两个特性是分别在邻域Δ_m和Δ_σ上计算的强度的局部平均和局部标准偏差。这些邻域是基于上皮细胞的大小。邻域越大，边界的定位越差，但是局部平均和标准偏差的估计越精确。在优选实施例中，已经实现一种折衷，其中Δ_m和Δ_σ的典型值分别是7和13。典型的上皮细胞对于用于这种情况的横向放大约为13个像素。

一个掩蔽成分M已经引入到同现矩阵的形成中，从而

$S(l_1,l_2,......,l_K)=\underset{x}{\mathrm{\Σ}}\mathrm{\δ}(1;M\left(x\right))\underset{k=1,K}{\mathrm{\Π}}\mathrm{\δ}(l_k;P_k\left(x\right))...\left(6\right)$

其中

用在优选实施例中的特性是以标准方式计算的局部平均和局部标准偏差，但是掩蔽了将被排除的像素，即

$P_1\left(x\right)=M\left(x\right)\frac{\underset{v={\mathrm{\Δ}}_m}{\text{\Σ}}M\left(v\right)I\left(v\right)}{\underset{v={\mathrm{\Δ}}_m}{\mathrm{\Σ}}M\left(v\right)}...\left(8\right)$

$P_2\left(x\right)=M\left(x\right)\sqrt{\frac{\underset{v={\mathrm{\Δ}}_{\mathrm{\σ}}}{\mathrm{\Σ}}M\left(v\right)\underset{v={\mathrm{\Δ}}_{\mathrm{\σ}}}{\mathrm{\Σ}}M\left(v\right)I{\left(v\right)}^2-{\left(\underset{v={\mathrm{\Δ}}_{\mathrm{\σ}}}{\mathrm{\Σ}}M\left(v\right)I\left(v\right)\right)}^2}{\underset{v={\mathrm{\Δ}}_{\mathrm{\σ}}}{\mathrm{\Σ}}M\left(v\right)(\underset{v={\mathrm{\Δ}}_{\mathrm{\σ}}}{\mathrm{\Σ}}M\left(v\right)-1)}}...\left(9\right)$

如果像素x将被省略，则P₁和P₂的值是零，但是已经从PCM的形成中被排除。如果等式(8)或(9)的结果是不确定的(比如除以0)，则结果从PCM的形成中被排除。这没有明确地在等式(7)中规定。

掩蔽点的识别

在步骤33中PCM能够从将进行分析的图像中的所有像素来形成。但是，来自已知不感兴趣的区域的像素将引起矩阵复杂度增加，并且降低感兴趣矩阵的部分之间的差异。因此，在优选实施例中，使用如等式(7)定义的掩蔽操作来将特定像素从PCM的形成中排除，并在步骤32中进行识别。将排除的像素是以下那些：

●不是组织的一部分，即在样本之外或撕裂的部分或由收缩引起的间隙。

●是脂肪组织的一部分。

●是输送管的一部分。

●掩蔽在没有用来计算局部平均的像素，或有1个或更少用来计算局部标准偏差的像素的位置。

可以使用识别这些像素的任何算法。但是，在优选实施例中使用图5所示的算法。将被省略的像素的特征为接近饱和(与载玻片上的其他数据相比)，并且具有很低的局部标准偏差。如果有大量像素将被省略，则在由图像中的所有像素形成的灰度级直方图H_g中可能有一个谷，例如图6(a)中所示的61：

$H_g\left(i\right)=\underset{x}{\∑}\mathrm{\δ}(i;I\left(x\right))...\left(10\right)$

如果例如在图6(b)中没有设置阈值的适当的谷，则PCM从局部平均和局部标准偏差形成。这与在接下来的步骤33中形成的PCM不同，并比在接下来的步骤33中形成的PCM分辨率更高。因为局部平均从邻域中计算，并因此具有非整数值，所以有可能增加矩阵的分辨率。例如，如果使用9个像素的邻域，则分辨率能够增加9倍。实际上，增加一般不会这么大。一个局部平均直方图H_m通过把与用于低的局部标准偏差的标准偏差轴平行的PCM相加形成(参见图7)：

$H_m\left(i\right)=\underset{j}{\∑}S(i,j)...\left(11\right)$

其中PCM中的i索引是局部平均，以及j索引是局部标准偏差。如直方图分析所确定的，直方图H_m被限制在主峰值和次峰值的较低点之间，这将在下文中论述。

在设置掩蔽像素的优选实施例中过程的框图在图5中示出，并在下文中说明。关键目的是找到灰度级直方图中或局部平均直方图中的一个适当的阈值，从而输送管细胞、脂肪和组织撕裂能够从接下来步骤33中PCM的形成和分析中被省略。如果不能找到一个简单的阈值，则使用一个基于更复杂的一组标准的阈值。具有一个阈值并不总是可能的或理想的，在这种情况下，不使用阈值操作，并且从接下来的PCM分析中不排除像素。

步骤         操作/测试                        动作

51           形成直方图H_g

52           识别H_g中峰和谷的位置

             v＝主谷的位置

             p₁，p₂＝主峰和次峰的位置

             H_g(p₁)＞H_g(p₂)

53           直方图中有一个谷？            否：转到步骤54

                                           是：

53a          测试谷。它在主峰之上？        否：转到步骤54

             v＞p₁？是：

53b          谷太低吗？                    否：使用谷作为阈值

             v＜200？                      来设置掩蔽点，转到

                                           步骤59

                                           是：转到步骤54

54           使用局部平均和标准偏差建立

             PCM

55           从PCM的一部分建立局部平均直

             方图H_m

56           识别H_m中峰和谷的位置

             v＝主谷的位置

             p₁，p₂＝主峰和次峰的位置

             H_m(p₁)＞H_m(p₂)

57         直方图中有一个谷吗？             否：转到步骤58设置默

                                            认阈值

                                            是：

57a        测试谷。它在主峰之上？           否：搜索主峰之上直

           v＞p₁？                         方图中的另一个谷。

                                            如果找不到谷则转到

                                            步骤58

                                            是：

57b        主峰的幅度大于次峰的X％？        (次是：使用谷作为阈

           峰在正被检测的谷之上)。          值来识别包括在掩

           H_m(p₁)＞XH_m(p₂)             蔽中的点。转到步

           为了降低对噪声的灵敏度，峰的幅   骤59。

           度通过包括矩阵中相邻元素的一

           个分量进行评估。这一般将通过用

           [121]或更大的算子进行卷积

           来完成。

           X一般是150％

                                            否：

57c        低于谷的至少H_m的X部分           是：使用谷作为阈

           (fraction)正在被检测吗？         值来识别包括在掩

           X一般是75％                      蔽中的点。转到步

                                            骤9。

$\frac{\underset{j=0}{\overset{v=1}{\mathrm{\Σ}}}{H}_m\left(j\right)}{\underset{j}{\mathrm{\Σ}}{H}_m\left(j\right)}>X$

                                            否：

          谷小于主峰和次峰幅度的            是：使用谷作为阈

          X部分平均吗？X一般是80            值来识别包括在掩

          ％                                蔽中的点。转到步

                                            骤59

                                            否：转到步骤58

58        如果所有的测试都失败了，则没有掩

          蔽点要被设置，或者数量非常小。使

          用H_m中的一个备用阈值来设置从主

          峰到图像中的最大强度之间的距离

          的X部分。这可能意味着没有掩蔽点

          要被设置。X一般为0.75

59        结束

可以使用任一适当的算法来确定直方图H_g或H_m中峰和谷的位置和相对重要性。但是，在优选实施例中使用以下算法。

用确定峰和谷位置和涉及峰的“重要性”的尺度(轨迹长度)的多分辨率方法来确定直方图中的峰和谷。一个小的孤立的峰被认为与大的峰一样重要。

直方图进行零延伸从而它从-2^n/2到2^3n/2，其中确定n从而2ⁿ是零延伸之前完全包含直方图的最小二进制数。这是通过下述实现的：首先延伸直方图直到其二进制数等于2ⁿ并设置新的二进制数为零，并且接着在两侧进一步延伸直方图50％从而它具有2ⁿ⁺¹个二进制数并设置新的二进制数为零。接着通过平均二进制数组形成一组n个多分辨率直方图，从而连续直方图具有2²，2³，2⁴，.....2^n-1，2ⁿ，2ⁿ⁺¹个二进制数，以及通过在每个分辨率用一个[1 -1]边缘算子卷积并识别符号的变化来确定峰和谷的位置。直方图是平坦的时候，则如果需要谷被适当调整。接着通过每个分辨率跟踪峰和谷的位置，峰和谷所在的位置分辨率越高，则轨迹长度越大并且感觉的峰或谷的重要性越大。这样的一个例子在图8中示出，其中主(或最重要的)峰不是最高的。

该分析使得像素将被掩蔽并从用于把图像分割为上皮和背景像素的步骤33的PCM的形成中排除。但是，非常小的孤立的像素组不应当被包括在掩蔽中。因此，识别出的像素用与以下参照步骤36更全面进行说明的相同的过程进行分组和过滤，对于斑点具有以下的标准。只有符合这些标准的像素组才被包括在掩蔽中。

标准                 掩蔽

            如果符合所有的标准就接受：

X尺寸       ≥3

Y尺寸             ≥3

面积              ≥9

宽高比            ≥0.25

密度              ≥0.25

PCM中分布的分析

必须确定步骤33中生成的PCM中的两个主要的2维高斯分布N(m_1x，s_1x，m_1y，s_1y)和N(m_2x，s_2x，m_2y，s_2y)的位置和范围。这种分布的一个例子在图9中示出。在优选实施例中，这通过使用例如在Journal of Royal Statistical Society B，39：1-38，1977，3(Maximum Likelihood from Incomplete Data via the EMAlgorithm，AP Dempster，NM Laird和DB Rubin)中描述的期望最大化(EM)算法执行方案在步骤34中实现。标准EM算法能够利用两个步骤的迭代应用来实现模型与符合数据的隐藏参数混合，这两个步骤的迭代应用估计用于混合模型的新参数，并且接着最大化该适合度。简单地说：

初始化：	估计初始混合模型参数，例如，分量高斯分布的平均和标准偏差。初始模型参数能够以任何适当的方式确定/估计。
		步骤1：	假定模型参数正确，则得到属于混合模型的给定分量即分布1或2的每个数据点的概率。重新估计数据分量和模型分量之间的关联的概率。这引起确定混合的“加权”数据组。
步骤2：	重新估计模型参数并且从步骤1迭代直到收敛或其他终止标准，比如混合模型中的误差小于某一分数。

在优选实施例中，加入标准算法上的变化以使得将两个分布的相关属性强加上。在这种情况下，模型分量的高斯幅度由它们的x和y标准偏差的积来定标，即拟合的模型成为 $\frac{N(m_{1x},s_{1x},m_{1y},s_{1y})}{s_{1x}{s}_{1y}}$ 和

$\frac{N(m_{2x},s_{2x},m_{2y},s_{2y})}{s_{2x}{s}_{2y}}.$

图像分割

在步骤35，图像被分割为三种类型：

●背景，它通常是亮的，并具有一个小的局部变化。

●上皮细胞，包括输送管周围的那些细胞，它们通常较暗，并具有较大局部变化。

●掩蔽像素。

这是一个硬分割任务T，其对于背景和上皮细胞是基于一个像素的局部特性形成的PCM中的地址和从该地址到PCM中最近的(规一化的)分布(由EM算法确定)的距离：

$T\left(x\right)=\mathrm{\δ}(1;M\left(x\right))\underset{k=\mathrm{1,2}}{\mathrm{k\δ}}\left(\begin{array}{c}\underset{k}{\min }\sqrt{{\left(\frac{\underset{{\mathrm{\Δ}}_m}{\mathrm{mean}}\left(I\left(x\right)\right)-m_{\mathrm{kx}}}{R_k{s}_{\mathrm{kx}}}\right)}^2+{\left(\frac{\underset{{\mathrm{\Δ}}_{\mathrm{\σ}}}{\mathrm{sd}}\left(I\left(x\right)\right)-m_{\mathrm{ky}}}{R_k{s}_{\mathrm{ky}}}\right)}^2};\\ \sqrt{{\left(\frac{\underset{{\mathrm{\Δ}}_m}{\mathrm{mean}}\left(I\left(x\right)\right)-m_{\mathrm{kx}}}{R_k{s}_{\mathrm{kx}}}\right)}^2+{\left(\frac{\underset{{\mathrm{\Δ}}_{\mathrm{\σ}}}{\mathrm{sd}}\left(I\left(x\right)\right)-m_{\mathrm{ky}}}{R_k{s}_{\mathrm{ky}}}\right)}^2}\end{array}\right)...\left(12\right)$

其中

和 (I(x))是局部邻域Δ_m和Δ_σ分别在由向量x索引的图像I中的局部平均和标准偏差。R_k是一个改变两个分布的相对重要性的定标因子，一般R₁＝1且R₂＝0.8。

参照图10这会很清楚。感兴趣的地址是(x，y)，并且该算法试图找到哪个拟合分布最接近：这对应于I₁和I₂的最小值。这些距离能够根据用适当的标准偏差并用定标因子R₁和R₂规一化的a，b，d和e确定。

这导致具有标记为掩蔽的、上皮或背景的像素的分割T具有以下值：

清洁(cleanup)图像

在步骤36，将分割中已经标记为上皮或背景的并且与用一个8像素邻域(见图11)相同标记的像素相邻的像素分组为斑点并按照以下步骤过滤。

清洁分割以消除小洞(与其周围的标记不同的斑点)、小的孤立的斑点、行和任何其他简单的赝象。使用的量度是斑点的X和Y尺寸、面积(按照像素数)、密度和宽高比。这些参数的值将涉及上皮细胞的放大倍数和尺寸。在具有X10大倍数的优选实施例中，如果需要则用于接受一组像素的标准能够改变，但是典型的最小和最大值列在下表中。如果拒绝一个斑点，那么用其周围来代替它。

标准           上皮细胞-洞               上皮细胞-孤立的组

          如果符合任一标准就拒绝       如果符合所有标准就接受：

          (即合并为细胞)：

X尺寸     ≤30                         ≥35，≤10000

Y尺寸     ≤30                         ≥35，≤10000

面积      ≤900                        ≥3000

宽高比    ≤0.15                       ≥0.1

密度      ≤0.1                        ≥0.1

密度被定义为在最小闭合矩形的面积上的斑点中的像素数。宽高比是斑点的短轴和长轴的比。

借助于例子，图12示出分组成斑点的上皮(较暗的阴影)和背景(较亮的阴影)像素。上皮斑点121和122具有分别包括背景斑点123和124的洞，其中大洞123将被保留而小洞124将被拒绝并且再次标记为上皮。小的上皮斑点125、126和127将由于尺寸而被拒绝，而长的薄上皮斑点128将由于宽高比被拒绝，并且所有的都将被再标记为背景。

注意在该步骤36中，掩蔽像素(在步骤32中分组和过滤)被忽略，所以例如由于输送管或脂肪细胞(在该阶段标记为掩蔽)的存在引起上皮斑点中任何尺寸的洞被保留，但是为了便于说明在图12中没有示出。

识别输送管像素

在步骤37，标记为掩蔽的但满足简单形状和尺寸标准的像素或斑点组被再次标记为输送管像素。优选实施例中的典型值将是：

标准                输送管

             如果所有标准都满足就接受：

X尺寸        ≥3，≤100

Y尺寸        ≥3，≤100

面积         ≥9，≤5000

宽高比       ≥0.25

密度         ≥0.25

这些值不是特别关键的，但还是被选择，从而拒绝长的组织撕裂、收缩裂缝和大多数脂肪细胞。

在这个阶段，像素被标记为掩蔽、上皮、背景或输送管，其值如下：

扩大细管种子点

步骤38的目的是以一个等于两个一般上皮细胞直径的量将步骤37中识别的输送管像素的斑点扩大成周围的上皮斑点，从而扩大的输送管斑点于是等于原始的组织样本中对应的细管结构的假定尺寸。

输送管细胞应当由上皮细胞包围。由于载玻片准备的方式，所以这可能不总是输送管偶尔由背景细胞部分包围的情况。已经标记为输送管但其边界没有足够的相邻上皮细胞的任何斑点将被再次分类为掩蔽，而不是输送管。在优选执行方式中，至少边界的15％必须具有通过计数边界像素确定的上皮细胞。

而且，在标记为输送管的像素和标记为上皮的像素之间可能有非常小的间隙。这是由用在分割中的邻域引起的。邻域越大，该间隙就会越大。在优选执行方式中，通过重复应用以下的扩大技术直到整个扩大等于两个上皮细胞的直径来允许该扩大穿过比如说两个像素的间隙。

有很多方法能够用于扩大，最值得注意的是形态学。但是，在优选实施例中，图像I用一个2维高斯核N(0，s)卷积以生成一个新的扩大图像D。高斯核具有一个零平均和一个标准偏差s，从而高斯具有在半径1内的值1，其大于将被桥接的间隙。更详细地参考以下算法：

1.生成一个新图像D，其中受到扩大的分割中的所有像素被设置为1，并且所有其他像素被设置为0。

2.图像D用一个具有标准偏差的高斯核卷积，从而高斯的值在扩大的理想量时为1，即如果将被穿过的最大间隙是2个像素，那么理想的扩大量是3个像素，并且来自中心的3个像素的高斯核的值将是1。卷积的结果被下舍入，即为值0或1。

3.扩大图像中是1的以及将被扩大到原始图像中的分类的所有像素都被设置为扩大分类，否则它们保留在原来的状态。这意味着标记为输送管的像素斑点能够穿过另一个标记的间隙扩大到标记为上皮的斑点。能够穿过的最大间隙是1，小于扩大量。

重复该算法的应用使得输送管斑点穿过较小赝像扩大到上皮斑点而没有不受控制的行为。

一个扩大的单个步骤用下式确定

                D(x)＝δ(T(x)；输送管)

     D₁(x)＝δ(int(D(x)N(0，s))；δ(T(x)；输送管或上皮))

                         ·                            (15)

                         ·

     D_n(x)＝δ(int(D_n-1(x)N(0，s))；δ(T(x)；输送管或上皮))

                 D(x)←D_n(x)

其中是用于数字卷积的算符。开始，将被扩大的图像D(x)包括已经标记为输送管的所有像素。接着这用适当范围的2维高斯核卷积并且结果通过截断转换为整数。现在标记为输送管并且先前标记为输送管或上皮的像素被保留为输送管，所有其他像素被设置为0。该过程重复进行直到实现理想的扩大水平，即等于两个典型的上皮细胞，在优选实施例中这将是13次迭代。

该扩大的形象化结果在图13中示出。在图13(a)中的扩大之前，在上皮斑点132和相邻背景斑点133内有一个输送管斑点131。斑点132表示原始载玻片中大量的上皮细胞，其在一些方向上延伸到包含用斑点131表示的输送管的细管结构以外。图13(b)示出扩大后同样的斑点，其在一些方向上被原始上皮斑点的外边界停止。

评估癌症的严重程度

在步骤39中，癌症的严重程度将基于的尺度M被计算为标记为输送管的扩大像素数与输送管和上皮像素的总数的比：

将看到这有效地对应于在当用眼睛对载玻片分级时病理学家使用的等式(1)中给出的尺度。

如前面表明的，临床医生一般使用10％和75％的判定边界来根据该尺度对癌症严重程度分级。但是，对于步骤40中的自动系统这并不是必然适合的，因为：

●在极端情况下，即真的75％条件下观察的人类视觉系统将倾向于感知得更高，而真的10％将倾向于感知得更低。

●在自动系统中计数的像素可能不需要与将由人类观察者包括的那些像素相同。但是，对一个像素是被排除还是包括在自动系统中的判断比人类观察者更一致。

由于这些原因，将被使用的实际判断边界最好由步骤41中的训练确定。这样阈值能够在尺度M上选择以确定等级1和2之间以及等级2和3之间的边界。这些阈值应当基于典型的和完全训练的数据。这可能在每个实验室基础上意味着很好。

Claims (34)

1、一种用于自动分析包括像素阵列的数字图像的方法，包括下列步骤：

使用所选择像素的邻域中强度的局部平均和局部标准偏差的特性来从所述像素的一些或全部生成特性同现矩阵(PCM)；以及

作为对所述PCM的分析结果，通过将选择的像素标记为属于特定分类来分割图像。

2、根据权利要求1所述的方法，其中使用期望最大化(EM)算法的执行来将各个高斯分布拟合成PCM内的两个主要分布以确定分布参数。

3、根据权利要求1或权利要求2所述的方法，其中像素按照PCM内它们最接近的分布进行标记并包括以下步骤：

将各个标记分配给PCM内的单独分布；

确定各个像素所属的PCM内的点和每个被标记的分布的中心之间的规一化距离；以及

给各个像素分配这样的规一化距离最短的分布的分布标记。

4、根据权利要求3所述的方法，其中定标因子被引入到规一化中以把标记向特定分布偏移。

5、根据前述任一权利要求所述的方法，其中由特定像素的局部特性引起这些像素被从所述PCM的形成中排除。

6、根据权利要求5所述的方法，包括在权利要求1中规定的步骤之前的下列步骤：

从图像像素的一些或全部中形成灰度级直方图；

作为分析所述直方图的结果来建立一个阈值；以及

从所述PCM的形成中排除超过所述阈值的那些像素。

7、根据权利要求6所述的方法，其中所述阈值被建立为在直方图的主峰之上最重要的谷。

8、根据权利要求5所述的方法，包括在权利要求1规定的步骤之前的下列步骤：

使用各个像素的邻域中强度的局部平均和局部标准偏差特性从图像像素中生成特性同现矩阵(PCM)，并具有比第一次提到的PCM更高的分辨率；

通过沿着用于局部标准偏差的小范围的固定的局部平均相加形成局部平均的直方图；

作为分析所述直方图的结果来建立一个阈值；以及

从第一次提到的PCM形成中排除超过所述阈值的那些像素。

9、根据权利要求8所述的方法，其中所述阈值被建立为在直方图的主峰之上的最重要的谷。

10、根据权利要求7或权利要求9所述的方法，包括下列步骤：

延伸各个直方图直到其二进制数等于2ⁿ，其中n是最低值，从而2ⁿ个二进制数完全包含直方图，并设置新的二进制数为零；

在两侧进一步延伸各个直方图50％，从而它具有2ⁿ⁺¹个二进制数，并设置新的二进制数为零；

通过平均二进制数的组形成一组n个多分辨率直方图，从而连续的直方图具有4，8...2ⁿ⁺¹个二进制数；

在每个分辨率用一个边界算子对直方图的值卷积；

通过所述卷积结果的符号变化识别每个分辨率中所有峰和谷的位置；

把每个分辨率中每个峰和谷与随后的分辨率中对应的峰或谷(如果有的话)相关联；以及

分别将主峰和最重要谷识别为以所述分辨率的最高数量存在。

11、根据权利要求5到10中任一所述的方法，其中也通过将从第一次提到的PCM的形成中排除的像素标记为属于与第一次提到的分类不同的特定分类来分割图像。

12、根据前述任一权利要求所述的方法，还包括将标记为属于任何所述分类中相同的一个分类的相邻像素分组为斑点的步骤。

13、根据权利要求12所述的方法，还包括计算涉及所述各个斑点的统计量并按照所述统计量对相同部分进行过滤的步骤。

14、根据权利要求13所述的方法，其中所述统计量包括各个斑点的阵列的一个或多个轴方向上的尺寸、面积、宽高比和密度中的一个或多个。

15、根据权利要求13或权利要求14所述的方法，其中所述过滤包括把在选择斑点中的像素再次标记为属于相应周围斑点中像素的分类。

16、根据权利要求13到15中任一所述的方法，其中所述过滤包括把在选择的斑点中的像素再次标记为属于与任何前述分类不同的新分类。

17、根据权利要求16所述的方法，还包括通过由再次标记为属于所述新分类(Cn)的像素组成的特定量的斑点扩大为由标记为属于第一次提到的特定分类(Co)中选择的一个分类的像素组成的相邻斑点的步骤。

18、根据权利要求17所述的方法，其中所述扩大包括下列步骤：

通过在新图像中分配原始图像中分类Cn的像素以1的值并分配所有其他像素为0来生成一个新图像；

用一个具有零平均和设置等于所述特定量的标准偏差的二维高斯核卷积该新图像，从而在与平均有1的标准偏差时高斯值是1；

截断产生的图像，从而它只包含值1和0；以及

如果在截断后的图像中一个像素具有值1并且在原始图像中被标记为类型Cn或Co，则给它在原始图像中分配类型Cn。

19、根据权利要求17所述的方法，其中所述扩大包括重复执行下列步骤：

通过在新图像中分配原始图像中分类Cn的像素以1的值并分配所有其他像素为0来生成一个新图像；

用一个具有零平均和设置等于预定值(L)的标准偏差的二维高斯核卷积该新图像，从而在与平均有1的标准偏差时高斯值是1；

截断产生的图像，从而它只包含值1和0；以及

如果在截断后的图像中一个像素具有值1并且在原始图像中被标记为类型Cn或Co，则给它在原始图像中分配类型Cn；

从而实现所述扩大的规定的量，尽管在所述由标记为属于类型Cn的像素组成的所述斑点和标记为属于类型Co的像素组成的所述斑点之间，存在标记为既不属于类型Co也不属于类型Cn的不超过L-1个像素的间隙。

20、根据任一前述权利要求所述的方法，还包括计算尺度的步骤，该尺度作为被标记为属于选择的所述类型的像素的数量的函数。

21、根据权利要求20所述的方法，当附加到权利要求17到19中任一个时，其中所述尺度是所述扩大之后标记为类型Cn的像素数与标记为类型Cn的像素数和标记为类型Co的像素数的总和的比。

22、根据前述任一权利要求所述的方法，用于自动分析组织或细胞标本的数字图像。

23、根据权利要求22所述的方法，其中图像是乳房组织的一部分。

24、根据权利要求22或权利要求23所述的方法，其中第一次提到的分割的结果是将选择的像素标记为属于上皮细胞的类型。

25、根据权利要求22到24中任一权利要求在被附加到权利要求16时或者任何其他权利要求在被附加到权利要求16时所述的方法，其中所述新类型被识别为输送管细胞的类型。

26、根据权利要求22到24中任一权利要求在被附加到权利要求17时或者任何其他权利要求在被附加到权利要求17时所述的方法，其中类型Cn被识别为输送管细胞的类型，并且类型Co被识别为上皮细胞的类型。

27、根据权利要求26所述的方法，其中所述扩大是在对应于特定数量的上皮细胞的距离上。

28、根据权利要求22到27中任一所述的方法，当附加到权利要求20或权利要求21时还包括把所述尺度转换为癌症等级的指示的步骤。

29、根据权利要求1到21中任一所述的方法，用于自动分析矿物标本的数字图像。

30、一种用于自动分析包括像素阵列的数字图像的方法，包括下列步骤：

通过将各个像素标记为属于两个或多个类型中的一个类型来分割图像；

将相同类型的相邻像素分组为斑点；

计算涉及各个所述斑点的统计量；

将选择的所述斑点中的像素再次标记为属于不同的所述类型；

将一个所述类型的选择斑点扩大一个特定量到另一个所述类型的斑点；以及

计算一个尺度，该尺度把由扩大的斑点覆盖的总面积和由选择的一个类型或多个类型的斑点覆盖的总面积相关联。

31、一种用于自动分析包括像素阵列的乳房组织的组织标本的数字图像的方法，包括下列步骤：

将像素标记为分别表示上皮细胞和输送管细胞；

将标记为表示输送管细胞的像素组扩大一个与上皮细胞的尺寸相关的特定量到标记为表示上皮细胞的相邻像素组；

计算在这样的扩大之后标记为表示输送管细胞的像素的总数以及标记为表示输送管细胞或上皮细胞的像素的总数；

根据前一步骤的计算来计算一个尺度；以及

把所述尺度转换为癌症等级的指示。

32、用于自动分析数字图像的设备，该设备包括适于执行根据任一前述权利要求所述的方法的装置。

33、一种包括计算机可该介质的计算机程序产品，在该计算机可读介质上具有适于引起计算机执行根据权利要求1到31中任一所述的方法的计算机程序代码装置。

34、一种计算机程序，包括引起计算机执行根据权利要求1到31中任一所述的方法的指令。

Images