CN101169868A - 用于分割数字医学图像的方法 - Google Patents

Abstract

本发明公开了一种用于分割数字医学图像的方法。一种基于马尔可夫随机场(MRF)的技术被描述以用于执行由少量或有限数据表征的图像的聚类。所提出的方法是统计分级模型，该模型基于图像像素的统计和上下文信息的描述来标记图像像素。除了评价源自K平均聚类方案的定义的像素统计之外，所述模型通过像素和它们的标记(上下文)之间的空间相关性的描述来扩展分析，因此相对于纯K平均聚类的结果而导致分割输出的不均匀性减小。

Description

用于分割数字医学图像的方法

技术领域

本发明涉及一种诸如乳房摄影照片之类的数字医学图像的上下文分割方法。

背景技术

女性乳房的解剖结构随着年龄而变化。在生殖年龄期间，乳房主要由导管、腺和脂肪组织组成。这散布有提供支撑和附着到胸壁的纤维组织。腺和纤维组织被共同称为纤维腺组织。乳房腺组织被称为乳房实质，并且由20..25个小叶(腺)组成，这些小叶负责生产乳汁，并且由集合在一起以形成较大导管的许多微管(导管)朝着乳头引流。每个乳汁生产小叶包含一簇或一圈细胞。小叶和导管的部分由脂肪围绕以进行保护，并且由纤维组织支撑。随着年龄的增长，导管和腺单元经历萎缩变化，并且逐渐由脂肪组织替代。乳房被将乳房组织附着到胸肌的韧带固定就位。除了乳头和围绕它的乳晕之外，乳房到处被普通皮肤覆盖。

纤维、导管和腺组织在X射线乳房摄影照片上呈现为暗的或“致密的”。另一方面，脂肪具有透明或明亮的外观。术语“乳房摄影密度(MD)”和“乳房摄影模式”被广泛地用来描述呈现在乳房摄影照片上的致密/明亮区域在乳房中的比例。在过去已经提出了乳房摄影实质模式的不同分级方法，例如Nottingham分级(5种模式，例如正常(N)，腺(G2，G1，G0)，发育异常(DS-轻度，DM-中度，DY-重度)，导管显著(P1，P2)和不确定(IND))、Wolfe分级(4种类别)和Tabar-Dean分级(5种模式)。

乳癌的两种主要类别是小叶和导管癌。

-小叶原位癌(LCIS)是在乳房的乳汁生产小叶中包含的细胞的数量、外观和异常行为急剧增加的情况。术语“原位”是指癌症的早期，并且被用来指示异常癌细胞存在但还未扩散超过它们初始发展处的组织的边界。尽管并不认为LCIS是癌症，但是诊断患有LCIS(也称为小叶瘤形成)的妇女在一生中以后有发展成乳癌的更高风险。

-导管原位癌(DCIS)是乳房中早期癌发展的最常见情况。“原位”再次描述了还未移动到它初始发展处的身体的区域之外的癌症。对于DCIS，癌细胞被限制在乳房中的乳汁导管，并且还未扩散到乳房脂肪组织中或身体的任何其他部分(例如淋巴结)。DCIS可以作为微钙斑(称为微钙化)出现在乳房摄影照片上。

LCIS和DCIS都可以发展成侵袭性癌症(浸润性小叶癌或浸润性导管癌)，其中癌症扩散到乳房脂肪组织中或身体的任何其他部分(例如淋巴结)，这被称为转移。

到目前为止，乳房摄影法在乳癌的早期症状的检测中已经成为最常用和最成功的工具，所述早期症状常常可以由微钙化或肿块的存在而得知。然而，由放射科医生执行的视觉分析仍然是非常复杂的任务，并且开发了许多计算机辅助检测/诊断(CAD)系统来支持他们的检测和分级。实际上，CAD系统分别对有经验的和无经验的放射科医生的检测效率的影响例如在C.Balleyguier、K.Kinkel、J.Fermanian、S.Malan、G.Djen、P.Taourel、O.Helenon的“Computer-aided detection(cad)inmammography：Does it help the junior or the senior radiologist？”，EuropeanJournal of Radiology 54(2005)(2005)90-96中进行了调查。在这两种情况下，CAD系统已经证明是一种用于检测的有效的支持工具，尽管其自主性仍有疑问。因此，由于问题的复杂性，自动或半自动系统仍然仅仅扮演放射科医生的信令工具的角色。

在CAD环境中，图像处理的作用之一将是检测对于给定筛检或诊断应用而言需要进一步处理的感兴趣区域(ROI)。一旦检测到ROI，后续任务将涉及所述区域的表征以及将它们分级成几个类别中的一个。

乳房摄影病变的分级系统由美国放射学会提供，并且被称为BI-RADS(乳房成像报告和数据系统)；除了X射线乳房摄影之外，包括关于乳房的超声和磁共振成像(MRI)的部分。描述乳房摄影的发现的特征由描绘所述特征的线形图来说明，之后是几个乳房摄影的例子。

乳房摄影中ROI的例子是(a)钙化、(b)肿瘤和肿块、(c)胸肌、(d)乳房轮廓或皮肤-空气边界。分割是将图像分成其组成部分、对象或ROI的过程。在图像或其组成部分(即乳房摄影病变)的检测、描述、识别或分级可以进行之前，分割是必需步骤。

诸如乳房摄影照片之类的放射图像典型地由三个主要区域组成：

-诊断区域包括对应于患者解剖结构(即乳房)的像素。通常，该成像区域的轮廓可以采用任何形状。

-直接曝光区域是接收未衰减辐射的图像区域。尽管该区域具有仅仅由噪声破坏的恒定强度，但是入射能量(例如X射线源Heel效应)和接受器(例如计算机射线摄影中变化的存储磷光体灵敏度)中的不均匀性会使该模式失真。在欧洲专利申请1256907中公开了一种方法，该方法从诊断图像中回顾地估计这些整体不均匀性，并根据外推的背景信号来使所有图像部分中的响应变平。

-如果有的话，准直区域作为高度衰减的像素出现在图像上。这些区域的形状典型地是直线的，但是圆形的或弯曲的准直形状也可适用。

在放射图像中的这些主要区域之间，可以考虑三个不同的区域过渡类型：诊断/直接曝光、诊断/准直区域、以及直接曝光/准直区域的边界。

分割算法旨在检测并分离构成受分析的对象(一个或多个)的一组像素。这些技术可以根据应用于图像的处理的类型被广义地分类。基于区域的算法根据合适的相似性标准对图像中的像素进行分组。在欧洲专利申请EP887769中公开了一种基于区域的算法，该算法通过根据灰度值直方图的形心聚类对像素进行分组来分割直接曝光区域。基于边缘的算法根据相邻区域的灰度值差来分离图像中的高对比度区域中的图像像素。在欧洲专利申请610605和欧洲专利申请742536中公开了一种基于边缘的算法，该算法检测并描绘在单次或多次曝光图像上的准直区域和诊断区域之间的边界。在基于区域和基于边缘的方法中，模型可以被用来限制所分割的图像区域的外观或形状，以便遵守预定义的光度或几何约束。该范例的例子是所谓的主动外观模型和主动形状模型(AAM和ASM)。

由于针对内乳房区域的分析，所以这些技术都不适用，原因是：(a)它们是整体分割技术，其不适于乳房结构和病变的特定内容并且仅仅产生主要实体，例如乳房轮廓，或者(b)受到特定几何或光度约束，其不适用于乳房结构和病变的乳房摄影外观的大的可变性。所以，在当前的专利申请中，焦点集中于这样的技术，所述技术可靠地分割乳房皮肤轮廓内部的结构和密度。

依赖于分割步骤的检测步骤通常被认为是非常复杂的任务。例如，肿块是比周围组织更致密地成簇在一起的细胞群，并且可以在乳房摄影照片上由相对较小的强度变化来表示。此外，乳房摄影扫描记录了乳房中存在的所有结构，其结构在尺寸、均匀性、位置和医学意义上不同。最后，数字乳房摄影在其真正本性上固有的特征是由于现实世界和数字表示之间的转换而引入的误差，即量化噪声。所有这三个因素增加了乳癌诊断任务的复杂性。

发明内容

本发明的一个目的是提供一种分割数字医学图像的方法，所述方法处理特别阻碍有用分割的数字乳房摄影照片的一般不均匀性。

上述目的由如权利要求1中所述的分割数字医学图像的方法来实现。在从属权利要求中阐述了本发明的优选实施例的特定特征。

当使用像素灰度值作为局部图像特征时本发明中的基本假设是，图像由具有基本上一致的灰度值的区域组成。

在图像不由基本上一致的灰度值的区域组成的情况下优选地使用多维特征，例如在局部导数的滤波器组的基础上计算出的多维特征。

通常，聚类方案通过为每个像素一个分配表示级别的值来将输入集(具有其强度值的图像)分成一系列组，所述组被称为簇。有效聚类的目的是加重存在于图像中的对象之间的边界，而不管尺寸、位置或可见度水平。聚类的最常用和最公知的技术之一是如普通教科书例如R.O.Duda，P.E Hart，D.G.Stork，Pattern classification，John Wiley and Sons，Inc.，2001中所述的K平均。该技术的特性和缺陷在乳房摄影分析的上下文中进行论述。

根据本发明，马尔可夫随机场模式被发展，并且作为扩展被应用于一般K平均方法。显示出了简单相关因子的引入如何显著地改善结果。

在下面描述了K平均方案并且引入MRF算法。举例说明了在乳房摄影数据分析上的主要应用。根据本发明的基于MRF的聚类模型的详细定义与其参数和特征被一起公开。最后，由所述算法获得的关于微钙化和肿块的检测的结果被注释，并且与经典K平均例程的结果进行比较。

根据后面的描述和附图，本发明更多的优点和实施例将变得显而易见。

附图说明

图1是来自MIAS数据库的具有(a)微钙化和(b)指示的肿块的例子。

图2是图1(a)的放大片段，其中(a)发生微钙化，(b)是其直方图均衡形式。

图3说明利用5个级别聚类的K平均。

图4示出利用10个级别并且利用(左)和不利用(右)形心重算的聚类的结果。

图5示出表示微钙化的簇。

图6是利用7个(左)和8个(右)级别聚类的结果；指示的区域是记录的肿块。

图7是由低密度组织(左)和聚类的结果(右)表征的乳房的乳房摄影照片的片段。

图8是由高密度组织(左)和聚类的结果(右)表征的乳房的乳房摄影照片的片段。

图9示出n¹邻域系统和n¹中的团集(clique)类型。

图10示出n²邻域系统和n²中的团集类型。

图11是马尔可夫随机场聚类(MRFC)方法的操作的总框图。

图12是操作的详细框图。

图13示出像素的3×3邻域。

图14说明MRFC方法的多分辨率逼近。

具体实施方式

在当前的专利申请中，一种新的像素聚类模型被公开并且应用于数字乳房摄影照片的分析。所述聚类表示乳房摄影分析的更一般方法中的第一步，并且旨在创建简洁数据集(簇)，以用于一方面是乳房摄影结构(例如胸肌和乳房轮廓)、另一方面是病变(例如肿块和微钙化)的自动检测和分级。病变典型地是在乳癌的早期诊断中分析的第一症状。图像像素由它们的强度(灰度级)或多维特征向量来描述。可以通过应用多尺度滤波器组或纹理算子来获得多维分解。

引入了基于马尔可夫随机场(MRF)的技术，该技术适合用于执行由少量或有限的数据所表征的图像的聚类。所提出的方法是统计分级模型，其基于图像像素的统计和上下文信息的描述来标记图像像素。除了评价源自K平均聚类方案的定义的像素统计之外，所述模型通过像素和它们的标记(上下文)之间的空间相关性的描述来扩展分析，因此相对于纯K平均聚类的结果导致了分割输出的不均匀性的减小。

乳房摄影数据

被称为乳房摄影照片的乳房摄影检查被用作一种筛检工具，用于在无症状的妇女中检测早期乳癌以及在有症状(例如硬块、疼痛或乳头溢液)的妇女中检测和诊断乳房疾病。不能感觉但是可以在乳房摄影照片上看到并且告知可能的癌发展的异常之一是存在微钙化，所述微钙化是乳房组织中钙的累积。微钙化的尺寸从0.1mm到5mm之间变化，并且在乳房摄影照片上钙群呈现为由较暗的不均匀背景围绕的小点。告知可能的癌发展的异常的另一重要组是存在肿块，所述肿块是乳房组织中致密细胞的累积。在乳房摄影照片上细胞群呈现为由较暗的不均匀背景围绕的亮区。

与所有X射线图像一样，乳房摄影图像被各种性质的噪声源破坏，例如(a)由于X射线光子通量的离散特性而造成的采集噪声，(b)由于量化仓(bin)的有限范围而造成的数字化噪声，以及(c)由于投影图像中解剖结构的叠加而造成的解剖学噪声。因此，乳房摄影照片上的图像区域本来是不均匀的，并且受到这些噪声源的组合效应影响。

图1(a)是来自公众可得到的MIAS(乳房图像分析协会)数据库的在尺寸上减小的示例乳房摄影照片。所给出的图像是发生微钙化的情况；然而，异常的尺寸太小以至于不能容易地识别。在图2(a)上更容易识别该钙联系(liaison)，在该图中仅仅给出了相同乳房摄影照片的片段。仍然难以识别实际微钙化和周围组织之间的边界。即使是相同片段的直方图均衡形式的图2(b)也不允许容易地分割所述结构。图1(b)是发生肿决的情况，所述肿块用所划圈来指示。在识别实际肿块和周围组织之间的边界上仍然有困难。

K平均统计数据聚类

数值数据的聚类形成许多分级和建模算法的基础。聚类的目的是从大数据集中识别数据的自然分组以产生系统行为的简洁表示。K平均例程根据关于实体所提供的数据将输入集分成K簇。K平均聚类原则在诸如R.O.Duda，P.E.Hart，D.G.Stork的Pattern classification，John Wileyand Sons，Inc.，2001之类的关于模式识别的基础教科书中进行了介绍。数据例如可以基于像素的灰度值通过特征空间中每个实体一个点来表示，或者基于多尺度滤波器组分解或一组纹理描述符，它可以是多维特征空间。开始，K级形心被定义，并且每个点被分配一个级别，其描述“最接近”所述形心。每次分配之后在属于所述级别的点中重算形心。聚类以迭代方式重复所述分配，直到所有点和它们各自的级形心之间的距离之和被最小化。图3给出了通过应用K平均方案并提供灰度级作为像素的描述来将图像像素分成5个级别的结果。

可以清楚地看到，输入图像的空间不均匀性对所得到的簇集的形状和不一致性具有强烈影响。K平均聚类例程的已知缺陷是，它不使用空间约束并且假设每个簇由恒定的强度来表征。后一种假设很少是正确的，并且也适用于乳房摄影图像。减少问题的不完全解决方案可以通过将像素的坐标与灰度级值一起结合在所用数据中来实现。然而，在这样的情况下，K平均将图像分成K个空间相异簇，这是非常不希望的。也可以定义非常大数量的级别以满足乳房摄影照片聚类的要求。然而，它将大大增加输出的复杂性。

利用马尔可夫随机场建模

MRF理论为在贝叶斯框架内构造的视觉问题提供统计工具，其结合了结构和统计信息。MRF理论在S.Geman，D.Geman的公开文献“Stochastic relaxation，Gibbs distribution，and the Bayesian restoration ofimages”，IEEE Trans.on Pattern Analysis and Machine Intelligence，Vol.6(1984)721-741中进行了介绍。

然而，为了与前述公布中的理论一致，聚类被定义为像素标记，像素被称为位点(site)。

一组位点S＝{1，...，n}表示被标记的原始对象，即像素。向所有位点进行标记分配被定义为S上的随机变量族的实现，并且被表示为{X＝x}＝{X₁＝x₁，...，X_n＝x_n}(为了简单起见，记作x＝{x₁，...，x_n})。描述位点的特征空间也被看作是具有其实现y＝{y₁，...，y_n}的随机观测场。通常，该组观测向量不需要具有与该组位点相同的大小，例如当观测被定义成用于成对位点时。然而，在本发明的上下文中，观测的数量等于位点的数量。对于监督技术，位点的训练集S_tr也被定义。在该情况下的训练集由先前提取并标记(常常是人工的)的基元组成。本发明的技术是非监督技术。所以没有定义用于聚类的训练集。最后，标记集L={1，...m}被表示，其中每个标记对应于一个级别。值得一提的是，输出可以由被分配相同标记的许多空间相异像素簇组成。

MRF中的最佳标记满足最大后验概率(MAP-MRF)准则，该准则要求标记的后验概率P(X=x|Y=y)(在下文被表示为P(x|y))最大化，可以假设其遵循下列形式的吉布斯分布：

$P\left(x\right|y)=Z^{-1}\×e^{-\frac{U\left(x\right|y)}{T}}---\left(1\right)$

其中配分函数Z是归一化常数，U(x|y)是后验能量，T是所谓的系统温度并且在大多数视觉问题中被分配值1。吉布斯概率密度函数通过使用简化的贝叶斯法则P(x|y)∝P(x)p(y|x)被结合在判定函数中，其导致下列形式：

P(x)p(y|x)∝e^-U(x)e^-U(y|x)＝e^{-(U(x)+U(y|x))}    (2)

因此，来自等式1的后验能量现在被表示为先验能量项U(x)和似然能量项u(y|x)之和：

U(x|y)=U(x)+U(y|x)    (3)

吉布斯分布模型定义了位点i的标记，这取决于S中所有其他位点的标记。然而，被称为Hammersley-Clifford定理的重要成果在吉布斯随机场(GRF)和马尔可夫随机场(MRF)之间建立了一一对应等价。MRF模型是一种条件概率模型，其中位点标记的概率仅仅取决于给定邻域内的位点标记。给定从本节开始时的记法，S的邻域系统由N＝{N_i|i∈S)给出，其中N_i是与i相邻但不包括f的所有位点的集合。

$P\left(x_i\right|x_{S-\left(i\right)})=P\left(x_i\right|x_{N_i})---\left(4\right)$

其中

$x_{N_i}=\left\{x_{i^{\′}}\right|i^{\text{\′}}\∈N_i\}---\left(5\right)$

依赖项的该“局部化”对于模型的复杂性的减小具有重要影响。相邻的位点被分组成所谓的团集，取决于模型的复杂性，所述团集可以由任何数量的相邻位点组成。团集的阶指示所包括的位点的数量，并且它的增大通常增加了模型描述的复杂性以及其计算复杂性。

这里定义了一阶的一元团集(其包括一个位点{i})和二阶的二元团集(其包括一对相邻位点{i，i′})。能量可以作为在团集上定义的局部势能之和被估计。给定位点的标记，每个团集的团集势能V(x_i)和V(x_i，x_i′)和似然势能V(y_i|x_i)被估计，其中i和i′标引S中的位点。

先验能量是所有团集势能之和：

$U\left(x\right)=\underset{i\∈S}{\mathrm{\Σ}}V\left(x_i\right)+\underset{i\∈S}{\mathrm{\Σ}}\underset{i^{\′}\∈N_i}{\mathrm{\Σ}}V(x_i,x_{i^{\′}})---\left(6\right)$

而在由标记x_i调节的观测是相互独立的假设下，似然能量是所有似然势能之和：

$U\left(y\right|x)=\underset{i\∈S}{\mathrm{\Σ}}V\left(y_i\right|x_i)---\left(7\right)$

通过最小化如等式(3)中所定义的能量来实现从等式(1)导出的MAP-MRF解：

$\hat{x}=\underset{x}{\arg \min }U\left(x\right|y)---\left(8\right)$

找到能量的全局最小值是一个组合问题，其通常使用以下算法来实现：如J.Besag的“Spatial interaction and the statistical analysis of latticesystems”，J.Royal.Statistical Society Vol.36(1974)p.192-226中所介绍的确定性迭代条件模式算法(ICM)，或者如在S.Geman，D.Geman的“Stochastic relaxation，Gibbs distribution，and the Bayesian restoration ofimages”，IEEE Trans.on Pattern Analysis and Machine Intelligence Vol.6(1984)p.721-741中所介绍的随机模拟退火(SA)算法。定义的组合特性暗示了进一步进行描述的模型的某些特征。

伊辛模型

伊辛模型是MRF中的自动逻辑(auto-logistic)模型的特殊情况，也被称为多级逻辑(MLL)模型。它被表示为在相邻位点上定义的简单约束。当相邻位点i和i′被分配两个不同标记时的情况由团集势能中的小常数来“处罚”，由此加强标记的均匀性；它可以以下列形式被表示：

$V(x_i,x_{i^{\′}})=\left\{\begin{array}{cc}\mathrm{\β},& x_i\≠x_{i^{\′}}\\ -\mathrm{\β},& x_i=x_{i^{\′}}\end{array}---\left(9\right)\right.$

其中β是小正常数。伊辛模型的简单性暗示了容易对标记过程包括某些约束。已知先验信息可以在模型中被编码以便增强所需的特性(一般(非)均匀性、定向(非)均匀性、纹理参数、...)。

操作的一般流程

图11示出操作的一般流程的图，其由下列步骤组成：

-初始化和开始条件的评价。

初始化的过程由计算级形心的位置和它们的γ函数的相应值组成。

-在所有图像像素上迭代

迭代的类型取决于在存在的几种可能性中选择的采样方案。特别重要的两种情况是：

○ICM-迭代条件模式。在该算法中沿着有序的像素列表来访问像素。该列表从开始贯穿至结束并且从结束贯穿至开始。该方案中的标记也被排序，并且该算法依次评价它们中每一个的能量。在访问所有像素之后，评价终止条件。

将ICM算法应用于乳房摄影分析的优点在于，一旦被正确地初始化，它将把较高标记分配给较亮区域，因而为分割区域的进一步分析自动地提供有用信息。

○SA-模拟退火。该算法在图像中随机选择像素。然后还随机选择标记，并且仅仅为随机选择的标记和像素的当前标记评价能量。

模拟退火需要明显更高数量的迭代。SA算法的优点归因于这样的事实，即它可以避免分割陷入局部最小值。缺陷可能是当应用于乳房摄影处理时SA会交换较高和较暗区域之间的标记，由此减少标记的含义。

-重算级参数。

每个标记改变触发了两个级形心和各自γ函数值的重算。要求重算的两个级别的选择由当前标记和新分配像素标记来确定。

-终止条件。

该算法评价在单次迭代期间改变其标记的像素的数量。重复迭代一直到不再有图像像素改变其标记。乳房摄影照片的最终分割由对应于每个标记的图像区域组成。

本发明的详细描述

聚类方案的设计困难之一是级别的数量的选择。一些例程被设计用于聚类成数量尽可能少的级别，并且具有最小的信息损失。然而在本情况中，尽可能少的方法将由于乳房摄影照片的特性而不会产生期望的结果。以非常简单的方式，搜索乳房摄影照片上的肿块的过程可以被描述为搜索很大网格的较亮区域。而且，被搜索对象的出现不明显，并且通常它们的视觉特征非常类似于背景。实际上，强度上的微小差别是第一突出特征，并且众所周知的是肿块形成乳房组织内的稍微更亮的区域。在本情况下，数据由一组具有12位深度的强度图像来表示，所以，数据的特征在于相对较低的分辨率。不过，已经设计了相干的和相对较简单的MRF模型，其克服了所述问题并且执行有效聚类。

令S＝{1，...n}是位点即像素的集合，其中每个像素i由其灰度级值y_i来描述。表示每个像素的标记分配的随机变量的集合由x＝{x₁，...，x_n}表示，并且是所谓的配置。如前一节中所示，由于标记分配过程部分地取决于标记本身，因此配置x扩展了模型描述。最后，我们定义标记集L＝{1，...，m}，其中每个下标表示一个级别，该级别是一维灰度级特征空间中的子区域。

邻域系统和团集

本发明的应用涉及在图像像素的有限N₁×N₂矩形网格上定义为L＝{(m，n)∶1≤m≤N₁，1≤n≤N₂}的离散随机场。当且仅当像素的邻域η_mn仅仅包括作为(i，j)的邻域的(k，l)的像素时，网格L的像素的子集的集合是L上的邻域系统。这里使用邻域系统的分级有序的序列。也被称为一阶邻域系统的最近邻域模型使用每个像素的最近四个邻域，二阶模型使用所有八个邻域。与网格-邻域对(L，η)相关联的团集c被定义为L的子集，使得c由单个像素组成c＝{i}，或由一对相邻位点组成c＝{i，i′}，或由三个相邻位点组成c＝{i，i′，i″}组成，等等。这些邻域系统和它们的相关团集在图9和图10中被描绘。

分割的特性可以与分配给每个团集类型的它们的控制参数一起被控制：单团集位点的参数控制图像中的每个区域类型中的像素的百分比；包含多于一个的位点的团集的参数控制聚类的大小和方向。

初始标记

基于MRF的例程的效率常常取决于起始点。在本申请中，该方案被定义为确定过程；所以，存在收敛到并不表示MAP-MRF解的局部最小值的危险。所以，通过为每个像素分配一个级别来执行初始标记，其灰度级最接近所述级别的形心。所得到的图像是标记的散布集，在该图像上容易看到当前噪声和不均匀性的影响(参见图3)。为了简单起见，我们将一维特征空间分成等尺寸区域，并且将级形心置于所述区域的中心，即形心沿着特征空间均匀分散。初始标记根据预定数量的级别产生抖动图像。在这一点上值得注意的是，初始标记等于K平均算法的第一次迭代的结果。

能量函数

在现有的问题中，能量函数扮演双重角色。一方面，它必须引导例程给像素i分配标记，所述标记关于像素的强度值、因此其距级形心的距离而指向其级别。另一方面，它必须引导例程将相邻像素分组成一致的、在空间上均匀的簇。

关于在2.3节中描述的和在等式(3)中以简单形式给出的基本能量的形式，似然能量被定义为似然势能之和。在本申请中，势能表示像素距级形心的欧几里德距离，从而产生：

$V\left(y_i\right|x_i)=\left|\right|y_i-C_{x_i}\left|\right|=|y_i-C_{x_i}|---\left(10\right)$

其中C_j，是级别j的形心的值。所定义的似然势能导致标记的正确分配，不过，它由于在乳房摄影照片中存在的噪声而导致空间非相干性。如果团集势能未被定义(V(x)=0)，则所给出形式中的标记等效于K平均，并且欧几里德距离被选作距离量度。由该例程所获得的结果被显示在图3中。在等式(9)中定义的团集势能表示能量函数的第二角色。该部分负责从标记中除去空间非相干性。它使用的唯一信息是一对相邻像素的实际标记，并且采用如在S.Z.Li“Markov Random FieldModelling in Computer Vision”，Springer Verlag，1995中给出的改进的类伊辛形式：

$V(x_i,x_{i^{\′}})=\left\{\begin{array}{cc}0,& x_i=x_{i^{\′}}\\ 1,& x_i\≠x_{i^{\′}}\end{array}---\left(11\right)\right.$

势能“阻碍”了例程将不同标记分配给两个相邻像素。

能量函数和归一化参数的一般形式

组合两个所定义类型的势能的问题起因于不同的值范围。似然势能的范围取决于级别的数量、灰度级域的大小和级形心的位置，而团集势能的范围取决于邻域系统的大小。所以，在试图归一化势能的过程中，两个参数α和β被引入以平衡(归一化)能量函数。实际上，利用这两个参数，可以在分配给势能的类型的偏好方面调节模型的性能。通过如下对其扩展来修改来自等式(3)的能量形式：

$U\left(x\right|y)=\mathrm{\α}\underset{i\∈S}{\mathrm{\Σ}}V\left(y_i\right|x_i)+\mathrm{\β}\underset{i\∈S}{\mathrm{\Σ}}\underset{i^{\′}\∈N_i}{\mathrm{\Σ}}V(x_i,x_{i^{\′}})---\left(12\right)$

上面给出的形式假设参数α和β是常数。然而，在后面各节中显示出，它们需要与级形心一起被重算。这两个参数的函数在前面已被描述。它们对算法的影响是在期望性能的方面平衡例程并且导致它收敛。它们在根据两种特殊情况的极端之间影响例程的性能：

·α＝0和β＞0

由于当使用各向同性邻域(例如3×3邻域)时β的值在所有方向上相等，因此标记的所得到的集合将仅仅由一个标记组成。在初始标记中出现最多的标记跨过网格分散。

·α＞0和β＝0

仅仅评价似然势能，并且标记方案简化为经典K平均算法。

值α和β对于能量函数的实际全局最小值具有直接影响。对于这两个参数的人工设置可以迫使模型呈现某些特性(加强均匀性，最小化在像素值和它们各自的级形心之间的距离，...)。可以以这样一种方式来设置参数，使得两个相应势能的“重要性”被平衡，即它们根据势能的范围被近似，以便在试图归一化势能的过程中使二者相等。在一个级别内，即在分配给一个级别的像素的集合中，我们定义在灰度级域中的最大发生距离：

$\mathrm{MaxDist}\left(y_i\right|x_i)=\underset{j,x_j=x_i}{\max }|y_j-C_{x_i}|---\left(13\right)$

如下计算所定义的最大可能距离α：

α(y_i|x_i)＝MaxDist(y_i|x_i)^-1(14)

显然，参数的该定义暗示了其对级形心的位置的依赖性，因此，当形心移动它们的位置时需要重算。以类似的方式计算参数β；仅仅团集势能的域取决于假定的邻域系统。由于我们将邻域系统定义为8个邻域，因此当具有标记x_i的像素i被具有不同于x_i的标记的8个像素围绕时出现最大值。所以，我们计算：

β＝8^-1(15)

8个邻域保证了在再分配标记时的各向同性；可以通过考虑如图9和图10中给出的定向团集来加入方向性。

改进的能量函数

如前一节中所示，参数β是常数，所以，它可以通过下列代入从等式中被除去：

$\mathrm{\γ}\left(y_i\right|x_i)=\frac{\mathrm{\α}\left(y_i\right|x_i)}{\mathrm{\β}}---\left(16\right)$

在等式(12)中，α也是常数，但是可能会显示出，它的值遵循标记分配，而且，对于每个级别需要评价它。原因在于这样的事实，即尽管级形心改变了它们的位置，但是级别的范围也发生变化，这暗示了像素的强度和形心之间的距离分别遵循所述变化。由于α已被定义为级别内距离量度的归一化参数，因此它必须与级形心一起被重算，并且由于新定义的参数γ是所分配标记和描述像素的特征的函数，因此能量函数采用改进的形式：

$U\left(x\right|y)=\underset{i\∈S}{\mathrm{\Σ}}\mathrm{\γ}\left(y_i\right|x_i)V\left(y_i\right|x_i)+\underset{i\∈S}{\mathrm{\Σ}}\underset{i^{\′}\∈N_i}{\mathrm{\Σ}}V(x_i,x_{i^{\′}})---\left(17\right)$

级形心的值和γ′的值同时被重算；而且，遵循K平均的定义，在标记每次被分配给像素之后进行重算。

所述例程根据如在J.Besag的“Spatial interaction and the statisticalanalysis of lattice systems”，J.Royal Statistical Society，Vol.36(2974)p.192-226中所介绍的ICM采样方案沿着图像进行迭代。在每次迭代中，每个像素被访问，并且对于S中可用的每个标记评价其能量。产生最低能量的标记如接下来所解释的那样进行选择，各个级形心和函数γ(...)被重算，并且算法进入下一个像素。最后，每次迭代之后是计数改变其标记的像素，并且当计数结果低于预定义阈值时程序终止。

最低能量标记的选择

通过省略等式(17)中的图像环，由强度值y所表征的单个像素i的标记x的能量作为两个部分之和被估计：

$U\left(x_i\right|y_i)=\mathrm{\γ}\left(y_i\right|x_i)V\left(y_i\right|x_i)+\underset{i^{\′}\∈N_i}{\mathrm{\Σ}}V(x_i,x_{i^{\′}})$

其中V(y_i|x_i)是由函数γ(y_i|x_i)加权的似然势能，并且V(x_i，x_i′)是团集势能。N_i是如图13中所绘的相邻像素的下标的列表，代表性像素i用黑色绘出，并且其邻域用灰色绘出。N_i是这些灰色像素的下标的列表。

似然势能是在像素的强度值与由实际标记x_i指示的形心的值之间的距离的量度。

团集势能是所谓的罚值之和，并且当像素i和i′被分配不同标记时采用值一，否则采用零。这在当前的分割方法中是新颖的上下文术语，原因在于相邻像素的标记影响当前像素的新标记的分配。它负责填充在分割标记的图像中的局部间隙，由此减小由于噪声或乳房组织的解剖学变化而引起的分段外观。它因此克服了经典分割方法(例如K平均算法)的缺陷。

将所有可用标记依次分配给像素i并且估计相应的能量。选择产生能量的最低值的标记，并且将该标记分配给像素。选择“最佳”标记之后，该算法根据访问图像中的像素的顺序进入下一个像素。

操作的详细流程

图12的详图扩展了图11的总图。

值“变化计数”形成算法的终止条件。在迭代期间当没有图像像素改变其标记时，算法终止稳定的标记。在每次迭代之前它被设置为零。

值“倒序”控制迭代的方向。在经典形式中，ICM采样方案从第一个到最后一个、然后从最后一个到第一个地访问有序像素的集合。一旦图像像素已被“来回”访问，迭代就被终止并且终止条件被检查。

关于“变化计数”是否等于零的判定形成了算法的终止条件。

被称为“选择第一像素”的过程具有的含义是：按顺序选择第一像素，即，当“倒序”被分配值零时选择图像的第一像素，而当“倒序”被设置为一时选择图像的最后一个像素。

过程“根据顺序选择下一个像素”也由“倒序”的值控制。在零的情况下，算法移到图像的下一个像素，在一的情况下它进入前一个像素。

模型的性能

模型的性能在被应用于实际图像上时关于两个方面来分析：(a)级形心重算的影响，以及(b)级参数的数量的选择。在下文中使用具有微钙化和肿块的乳房摄影照片的代表性例子来评论这些方面。

形心重算。在每次标记分配之后，形心的重算在K平均方法的定义中具有其来源。如前所示，函数γ由于形心的重定位也需要重算。然而，当形心被设置成静态并且重算未被执行时，最终结果(参见图4)明显发生变化。没有重算的算法在较少数量的迭代中收敛，然而，结果指示了聚类的较低均匀性和所提取像素组的较低平滑度。

形心重算的原因是为了所公开算法的另一特征出现。在一些情况下，可能发生的是在标记的过程中级别中的一个将从图像消失，当形心为静态时这不会发生。实际上，由乳房摄影照片的高维度所产生的灰度级域的低密度和位点的高数量向我们保证，在初始标记之后每个级别将具有像素成员。然而，在进行中，如果算法检测到级别还未获得任何成员，即没有像素被分配给所述级别，则它不能重算形心和函数γ(...)。在这样的情况下，在初始标记之后两个参数被重置为由计算产生的数值。重置这些值由这样的事实证明，即如果一个级别从图像中消失，则算法允许它再次重现。

级别的数量。被人工分配给算法的参数是级别的数量。实验表明，小于10的数量对于聚类来说是满意的，其将不会从最终结果中除去微钙化和肿块。然而，进行进一步的分析以评估描绘这些病变所需的级别的最小数量。

为此，使用10作为级别的数量，具有记录的肿块的图像的集合被聚类。肿块总是被正确地描绘，这允许计算它们的灰度级特性：属于肿块的像素的灰度级范围覆盖由12位所表示的范围的大约1/7.5。所以，确定了8是正确提取肿块边界所需的级别的最小数量。实际上，图6示出分别使用7个和8个级别来聚类的结果。可以看到，聚类成7个级别并没有提取可疑区域。然而，另一方面，8作为级别的数量却成功地提取肿块。

对于微钙化适用类似的推理。同样在这里，小于10个级别的数量对于聚类来说是满意的，如果在围绕图像中期望微钙化出现的位置的感兴趣区域中局部地执行分析，其将不会从最终结果中除去微钙化。图5示出局部聚类成5个级别的结果。所指示的区域是保持它们的形状和大小的所记录的微钙化。然而，当应用于整个图像时，级别的数量必须增加以便保持这样的小区域。而且，级别的数量的选择由乳房摄影照片中呈现的乳房的密度来规定。由于X射线乳房摄影照片的叠加性质，因此不同组织彼此重叠，并且不同组织可以与微钙化重叠。出现在乳房摄影照片中的微钙化和其他结构之间的主要区别在于它们的大小和空间相干性：微钙化通常作为像素的圆形小组出现，其具有比周围背景更高的强度值。

图7和图8给出了存在微钙化的两个乳房摄影照片的两个片段，一个在主要由低密度组织组成的乳房中明显存在钙联系，在另一个中高密度的乳房组织构成难以检测的背景；微钙化的两种存在已由放射科医生的检查确认。当算法在整个图像上运行之后提取片段。第二图像明显显示出乳房摄影照片的更高密度；因此与用于第一乳房摄影照片的15个级别相对，它需要20个级别来提取微钙化。

聚类的较高分辨率(即更多可能的级别)(a)导致在低可检测条件中更好地检测微钙化，以及(b)导致更小微钙化或对象的提取；因此它增加了算法的真阳性率和灵敏度。然而，级别的较低数量避免了提取不期望的对象；因此它导致低假阴性率并且增加了算法的特异性。在分析给出了由低组织密度表征的乳房的乳房摄影照片的情况下，在存在微钙化时它产生更清楚的视图；然而，如所示，它也在致密组织中导致较低的检测率。

结果的分析让我们定义了一个模型，该模型实际上是无参数的。实际上，8作为级别的数量在90％的情况下做得好。然而，这种类型的估计取决于数据集的特性，例如它取决于对于数字化胶片乳房摄影照片所用的扫描器的特性。尽管未根据经验分配，但是级别的数量必须在扫描器数据或特定数据集的基础上进行估计。

算法加速

可以通过使用图像的多尺度分解来显著地加速马尔可夫随机场聚类(MRFC)算法。在该分解中，高斯低通(LP)滤波被应用于原始图像，并且该图像随后在每个维度上由因子2二次采样(由符号↓表示)，从而产生第1级模糊图像。该过程在每一等级i上被迭代许多次，从而产生所谓的高斯锥。在每次迭代i，上标度(由符号↑表示)粗糙细节的形式可以从前一精细细节的形式提取，从而产生带通滤波图像。这些带通滤波图像的序列被称为拉普拉斯锥。通过将带通滤波锥等级b_i转换成b′_i并且重组它们，可以获得增强的图像。在图14的虚线框中概括了这些操作。

从足够粗糙的等级开始(例如图14中的等级g₄)，上面给出的马尔可夫随机场聚类算法被应用于尺寸减小的图像上，从而产生分割的图像s₄。该等级的分割标记的初始化通过K平均聚类算法来获得。由于在该等级上像素的数量被大大减小，因此获得了ICM采样方案的非常快的执行。

该等级的分割的图像然后在每个维度上通过在每个维度上分割标记的像素复制由因子2扩展，并且该图像构成下一个较高分辨率等级(其在图1 4的例子中是等级3)的初始化。使用等级3的像素数据，MRFC算法确定新的分割标记，从而产生分割的图像s₃。由于已经从较粗糙等级正确地获得了大多数标记，因此像素将仅仅局部地改变标记，并且因此该迭代也将在非常少量的时间内到达其终止条件。

在前一段中描述的该过程在下一个较精细等级上被重复，直到达到最高分辨率等级，从而产生原始图像的分割s₀。在图14的虚线框的左边描绘了这些操作的链。所述操作也可以应用于增强的图像，其中对应于肿块的微小区域关于周围背景组织的对比度被增大。该操作链在图14中的虚线框的左边被描绘，并且在该例子中从被处理等级p₄开始。通常，增强图像的MRFC分割s′₀将不正好等于原始图像的分割s₀。然而，在临床环境中，原始粗数据大多数在PACS站不再可用，并且仅仅对比度增强的图像将可用。对比度增强处理的例子是Agfa-Gevaert N.V.的MUSICA^TM算法。

其他分级数据结构可以被使用，例如四叉树，正如在关于计算机视觉的普通教科书例如M.Sonka，V.Hlavac，R.Boyle的“Image Processing，Analysis，and Machine Vision”，2nd edition，Brooks/Cole PublishingCompany，1999中所述的。

可以通过使用快速多分辨率计算方案来缓解级别的合适数量的选择，原因在于优化(例如ICM)是耗时的步骤。

为此，许多级别被输入到加速算法，并且与每个级别数量相关联的输出分割被转发到CAD算法以供进一步处理。

朝著在多维特征向量上操作的MRF的扩展

如先前概述的MRF模型由于根据所考虑的团集的相邻像素的标记分配的相互作用而引入局部空间上下文。较大的空间上下文可以通过应用纹理算子来获得，而不是仅仅使用像素的灰度值；因此，由于关于类似纹理模式而不是仅仅像素的强度的均匀性，所以可以推动朝着类似标记的分组。

在这样的框架中，似然势能现在是像素的多维纹理特征向量和实际标记x_i所指示的多维级形心之间的距离的量度。团集势能仍是所谓的罚值之和，当像素i和i′被分配不同的标记时它采用值一，否则采用零。

滤波器组可以由不同尺度和方向的一组导数滤波器组成；每个滤波器被称为通道，并且提供多维特征向量的一个特征。例如，边缘滤波器对应于一阶导数；条形滤波器对应于二阶导数滤波器。这两个滤波器可以在多个尺度和方向上被计算。由于导数运算不是各向同性的，因此可以通过仅仅使用在所有方向上导致旋转不变性的最大滤波器响应来使负责给定尺度的滤波器各向同性。由带通滤波器组成的滤波器组的例子基于Gabor滤波器，该滤波器是高斯函数和复正弦的乘积。各向同性滤波器包括在不同尺度上计算的高斯滤波器和高斯拉普拉斯(LOG)滤波器。小波变换是可以作为数字滤波器组的树实施的另一个例子。

当使用多维特征时，优选地使用在纯欧几里德距离上的其他距离量度来解释给定级别的特征值中的可变性。

例如，可以使用加权欧几里德距离，所述权重与标准偏差成反比例，或者可以使用马哈拉诺比斯距离，其中权重由方差/协方差矩阵的逆矩阵表示。

对乳房摄影数据的应用

该算法的主要目的是描绘呈现在乳房摄影照片中的解剖学结构。该算法通过同时分析所有图像像素的灰度级值来实现该目的，这允许在像素所属的结构之间进行区分，同时增强被分配相同级别的标记的像素的空间分组。结果的特征在于解剖学结构的一致描绘，而不会使它们的边界变形。

该算法的另一优点是：当应用于原始输入即未经处理的图像时，实现了高性能聚类。实际上除了标记分组本身之外，以团集势能的形式的依赖因子的公式提供了典型使用的预处理技术(即模糊化、直方图锐化等)的替代，所述预处理技术被用来实现更一致或更平滑的输入以供进一步处理。在乳房摄影照片分析中所述特征尤其重要，原因是乳房摄影照片的特征在于高灰度级分辨率的事实；所以任何类型的处理需要一种方法来处理被扫描的乳房摄影胶片或数字检测器的自然噪声。该算法减少了对预处理图像的需要，并且以非常简单的方式克服了图像中存在的噪声的问题。

该算法的又一个优点是一致数据集的提取以适合进一步的处理。该算法的输出是在乳房摄影照片上的标记区域的集合。标记的值与图像的灰度级范围直接相关，即越高的标记是指越亮的区域。所以当涉及进一步处理时，比如肿块检测，根据标记的值来分析区域以提取ROI可能就已足够，而不用大量分析原始图像数据。

最后，该算法可以作为CAD(计算机辅助检测/诊断)系统中的快速预处理步骤而被集成。模型的数学公式允许简单的实施。当采样方案被选择成ICM时，该算法可以作为滑动窗口滤波器被实施。在ICM方案中，从开始直到末尾然后返回地访问所有图像像素。在每个点处，像素及其八个相邻像素被用在计算中以考虑到相邻标记的各向同性影响。所以，使用窗口宽度被设置为三的滑动窗口方案就已足够，并且仅仅窗口函数需要被实施，所述函数将选择中心像素的最低能量标记。

该方案提供了一种优化代码以便实现所需的实时性能的容易的方式。

本发明特别适合于分割存在于乳房摄影照片中的乳房摄影组织类型，例如

-纤维腺盘

-脂肪组织

-皮肤

-胸肌

-病变的感兴趣区域，例如包括微钙化簇的区域或包括肿瘤块的区域。

利用图6作为例子，纤维腺盘由标记的小集合所表示的图像区域的集合来表示。微钙化或肿块当存在时被表示为与小至仅仅一个分割标记相关联的纤维腺盘中的叠加或嵌入区域。在乳房区域周围的脂肪组织由与一个或两个标记相关联的细长区域来表示。胸肌在分割结果中自身呈现为对应于少量标记的图像的角落中的区域，其中最小的密度标记具有典型为矩形的边界。考虑到用于乳房摄影的CAD系统，这些区域的进一步自动分析现在在这些区域的几何和放射测量特性的基础上成为可能。

已经详细描述了本发明的优选实施例，现在对于本领域技术人员而言显而易见的将是，可以在其中进行大量修改而不脱离如由所附权利要求书中限定的本发明的范围。

Claims (22)

1.一种用于分割数字医学图像的方法，包括以下步骤：

a.在K平均聚类算法的基础上计算每个像素的初始分割级标记，并且计算相关的级形心和伽马函数；伽马函数对能量函数中的似然势能和团集势能的重要性进行加权；

b.计算与依次将每个标记分配给所述像素相关联的当前像素的能量量度，所述能量量度是以下两个部分之和：

-像素的特征值与由当前分割标记所指示的形心的值之间的距离的量度；

-被分配标记和至少一个相邻像素的标记的差的量度；

c.为当前像素选择产生最低能量的标记；

d.重算分割级形心和伽马函数；

e.进入下一个像素，直到图像中的所有像素都被访问；

f.迭代步骤b-e，直到图像像素的级标记不发生变化，以便产生图像的最终分割。

2.根据权利要求1所述的方法，其中距离的量度是欧几里德距离。

3.根据权利要求1所述的方法，其中所述距离的量度是加权距离。

4.根据权利要求1所述的方法，其中所述距离的量度是马哈拉诺比斯距离。

5.根据权利要求1所述的方法，其中所述相邻像素属于各向同性邻域。

6.根据权利要求5所述的方法，其中所述邻域被限定为3×3邻域。

7.根据权利要求1所述的方法，其中根据特定局部方向来选择相邻像素。

8.根据权利要求1所述的方法，其中像素由其灰度值表示。

9.根据权利要求1所述的方法，其中像素由多维纹理算子所提供的多维特征来表示。

10.根据权利要求9所述的方法，其中在局部导数的滤波器组的基础上计算所述多维特征。

11.根据权利要求1所述的方法，其中步骤(d)和(e)的顺序被交换。

12.一种用于分割数字医学图像的方法，包括以下步骤：

a)计算数字图像的多分辨率(MR)表示，

b)在所述MR表示中在粗分辨率等级选择一个图像作为原始图像，并且将K平均聚类算法应用于所述原始图像，

c)计算与依次将每个标记分配给所述像素相关联的当前像素的能量量度，所述能量量度是以下两个部分之和：

-像素的特征值与由当前分割标记所指示的形心的值之间的距离的量度；

-被分配标记和至少一个相邻像素的标记的差的量度；

d)为当前像素选择产生最低能量的标记；

e)重算分割级形心和伽马函数；

f)进入下一个像素，直到图像中的所有像素都被访问；

g)迭代步骤b-f，直到图像像素的级标记不发生变化，以便产生图像的最终分割；

j)将标记扩展到多分辨率表示的下一个更精细的分辨率等级；

h)迭代步骤(c)-(j)，直到原始图像被分割。

13.根据权利要求12所述的方法，其中所述多分辨率表示是高斯锥表示。

14.根据权利要求12所述的方法，其中所述多分辨率表示是四叉树。

15.根据权利要求12所述的方法，其中通过在所述多分辨率表示中的更高分辨率等级将标记局部复制到图像来执行扩展。

16.根据权利要求12所述的方法，其中步骤(e)和(f)的顺序被交换。

17.根据权利要求1所述的方法，其被应用于乳房摄影图像。

18.根据权利要求1 2所述的方法，其被应用于乳房摄影图像。

19.一种计算机程序产品，适于在计算机上运行时执行权利要求1所述的方法。

20.一种计算机程序产品，适于在计算机上运行时执行权利要求12所述的方法。

21.一种计算机可读介质，包括适于执行权利要求1所述的步骤的计算机可执行程序代码。

22.一种计算机可读介质，包括适于执行权利要求12所述的步骤的计算机可执行程序代码。

Images