基于Gibbs先验的磁共振灰度图像分割方法

381理论计算机科学电子笔记46（2001）网址：http://www.elsevier.nl/locate/entcs/volume46.html12页基于Gibbs先验的Eilat Vardi埃拉特·瓦迪1，2宾夕法尼亚大学生物工程系美国嘉宝Herman3天普大学关闭PA，USA摘要在早期的工作中，提出了一种利用吉布斯先验从投影重建某类二维二值图像的快速随机方法。在本研究中，我们介绍了一个随机分割的磁共振灰度图像的骨小梁吉布斯先验。我们展示了一些结果以及一些后处理，可用于清理有噪声的灰度图像的分割。1介绍分割处理从灰度图像创建二进制图像的问题（其中，在二进制图像中，1表示存在并且0表示不存在我们希望从灰度图像分割的对象文[1]、[5]提出了一种利用Gibbs先验从投影重建二值图像的随机算法，给出了较好的结果。这里报告的工作的主要目标是开发一个分割算法使用类似的方法。1第一作者得到了美国国立卫生研究院（NIH）培训补助金T32 CA 74781的支持。第二作者得到了NIH资助HL 28438以及美国国家科学基金会资助DMS 9612077的支持。 感谢实验室宾夕法尼亚大学视网膜微电路研究所为第一作者慷慨提供空间和计算机设施2电邮地址： eilat@seas.upenn.edu3电子邮件地址：gaborherman@netscape.net2001年由ElsevierScienceB出版。 诉 在CCBY-NC-ND许可下开放访问。瓦尔迪-赫尔曼公司382≤≤××在本文中，图像是一个矩形的二维像素阵列灰度图像是其中每个像素具有0和255之间的值（其中0对应于黑色，255对应于白色）的图像。二值图像是每个像素的值为0或1的图像。在接下来的二进制图像的讨论中，我们可以互换地使用黑色和0，以及白色和1。关于要分割的图像类别的适当的先验信息对于分割过程可能是有用的除了灰度图像中的固有信息之外，所需二值图像的局部属性也可以提供有用的信息。二进制图像ω的局部性质，定义在（一个任意的，但固定的数组包含）H像素（ω（h）是黑色或白色的每个像素索引为整数1h H），可以由吉布斯分布的形式1Q（ω）=exp（γZZHh=1 Ih（ω）），其中Q（ω）是图像ω出现的概率，Z是归一化因子（确保Q是在H像素的固定集合上定义的所有二进制图像的集合上的概率函数），γ是定义吉布斯分布的“峰值”的正参数由于我们取γ>0，Q（ω）随局部势之和的增加而增加在本文中，我们定义了以像素h为中心的团C（h），由h的3 3邻域中的9个像素组成。我们将注意力限制在吉布斯分布上，其中像素h处的局部势仅取决于二值图像ω中C（h）的元素的值。（为了使这个定义和后面的定义适用于任何地方，我们将二进制图像边缘附近的像素视为图像由两层像素的边界扩展在这里描述的实验中，我们使用全为零的边界。）对于（实际;即，而不是扩展的）二值图像，像素h的颜色只影响C（h）中像素的局部电位值这9个局部势只依赖于h的5 × 5邻域中的25个像素的值。因此，由切换像素h（在二值图像ω中）的颜色引起的概率Q（ω）的变化仅取决于ω中的那25个像素的值。我们假设，我们给出了一组典型的（给定的应用领域）灰度图像和其相应的黄金标准分割表示为二进制图像。我们将这些图像称为训练集。局部电位的正确定义对图像分割的成功起着重要的作用。 定义应反映特定应用领域典型图像的特征。局部势的定义方式是鼓励某些局部配置，例如均匀的白色或黑色像素簇。有很多可能的方法瓦尔迪-赫尔曼公司383×−≤≤≤≤≤≤≤ ≤ ≤≤定义局部潜能。下面是一个简单的方法。我们将512个可能的33二进制集团配置的集合划分为6个集团配置类型：0）全黑配置，1）全白配置，2）白色像素形成边缘的集团配置，3）白色像素形成凸角的集团配置，4）白色像素形成凹角的集团配置，以及5）任何其他集团配置。集团配置类型的工作定义如下。只有当在集团的八个外部像素中有k个连续的白色像素和8k个连续的黑色像素时，配置才是第一个五种类型之一如果中心像素是黑色的，则配置是黑色区域、凸角、边缘或凹角，当且仅当k = 0，1 k2，k= 3，或4k 5。类似地，如果中心像素是白色的，则配置是凸角、边缘、凹角或白色区域，分别当且仅当3k 4、k= 5、6k 7或k=8。 否则，集团配置是其他类型设th（ω）为图像ω（0≤th（ω）6）。<（有关局部势和集团配置的更多信息，类型，参见[3]。） 我们通过以下方式从二进制训练图像中创建先验信息。我们创建一个大小为6的表r我们将二进制训练图像中的配置类型t的内部团的数量除以类型t的不同团配置的数量存储在其第t<(By内部集团，我们指的是完全包括在实际图像中的那些集团，而没有前面讨论的边界表值r（t）给出了某个类型t的团出现的平均次数。局部势Ih（ω）定义为ln（r（th（ω））+1）。这些局部电位为我们提供了关于二进制训练图像中不同团配置类型出现频率的信息所得分布的有用性取决于训练图像集合的大小以及它们在应用领域的代表性正如我们在第3节中进一步讨论的那样，我们调整吉布斯分布Q（ω）定义中参数γ的值，使得吉布斯分布随机图像中白色像素的预期百分比将是二进制训练图像中白色像素的平均百分比。我们还使用训练图像来关联二值图像和灰度图像的特征设g（h）为C（h）中灰度值的平均值。对于每个团配置类型t（0t 6），我们计算给定二值图像集中所有类型t的像素h上g（h）的相应值的平均值和标准差。我们分别用µt和σt表示每种类型t的平均值在本文中，我们把μh（ω）写成μth（ω），把σh（ω）写成σth（ω）。分割算法应该找到一个二进制图像，它不仅是统计上与灰度图像一致，但这也是二值图像的吉布斯分布的典型特征因此，如果q（g|ω）是概率瓦尔迪-赫尔曼公司384≤≤≤≤HHg（包含g（h）值的图像，1h H），给定二值图像ω，则贝叶斯M（ω）={Q（ω）q（g|ω）}β，其中Q（ω）是如上且β= 1。 假设上一段中的g（h）（1h H）独立且正态分布，我们有：H1−（g（h）−µh（ω））2q（g|ω）=h=12016 - 05 -25 01：01：02（2 σ2（ω））。由于我们的任务是找到使M（ω）最大化的ωε，并且这个ωε与β无关（只要β>0），出于算法原因，我们使用β >0时M（ω）的更一般形式，而不仅仅是β= 1。退火时间表是一种指定初始逆温度、递增时间表和停止标准的时间表这里，参数β是温度的倒数。它不是保持恒定的，而是在分割算法中递增，如退火时间表所述我们使用迭代Metropolis算法[4]来估计使M（ω）最大化的二值图像ωε我们从一个完全白色的图像开始算法。在迭代步骤中，通过随机选取单个像素hJ并改变其颜色，将当前图像ω1变为ω2设p定义为p=M（ω2）M（ω1）= exp（β{γ（Ih（ω2）−Ih（ω1））h∈C（hJ）σ h（ω 2）（g（h）− µh（ω 2））2（g（h）−µh（ω1））2−（ln（σ）+−）}）。（ω）2σ2（ω）2σ2（ω）一小时二小时一小时如果p大于或等于1，则图像ω1被ω2替换，如果p小于1，则图像ω1以概率p被ω2替换这相当于用概率min（p，1）的ω2替换图像ω1该程序在不同的β值下运行120，000次循环。在每个循环中，算法执行H次迭代。程序输出序列中具有最高M（ω）值的图像ω2减少运行时间我们设计了几种方法来减少分割程序的运行时间创建查找表：设ω为二值图像，hJ为像素。对应于改变像素hJ的颜色的p的值取决于像素的局部电势瓦尔迪-赫尔曼公司3852σ2（ω2）−H哈哈哈。−--×啊啊啊JHJHJHJ其集团包含HJ。因此，p取决于（但不限于）ω中hJ的5×5邻域中的25个二进制像素值。考虑两个图像ω1和ω2，其中ω2通过改变随机选择的像素hJ的颜色从ω1获得。注意，p的指数中的一部分项可以写为：（g（h）−µh（ω2））2H（g（h）−µh（ω1））22σ2（ω1）{σ2（ω1）−σ2（ω2）}=g2（h）·h h2σ2（ω1）σ2（ω2）H H{σ2（ω2）µh（ω1） −σ2（ω1）µh（ω2）}H+g（h）·Hσ2（ω1）σ2（ω2）H H{σ2（ω1）µ2（ω2）−σ2（ω2）µ2（ω1）}2σ2（ω1）σ2（ω2）我们定义H Hσh（ω2）BhJ（ω1）=（γ{Ih（ω2）−Ih（ω1）}−ln（σ）（ω）H1h∈C（h）2 2 2 2{σh（ω1）µh（ω2）−σh（ω2）µh（ω1）}）。2σ2（ω1）σ2（ω2）H H设B是BhJ（ω1）的值，BJ是BhJ（ω2）的值。 很容易证明B=BJ。此外，值对B、BJ不取决于hJ的颜色。如果我们给出24个像素的颜色，5 5邻域（不包括hJ的值），我们可以同时计算B和BJ。然而，由于等式B=BJ，我们只需要存储这两个值中的一个（我们将在下一段中讨论如何做到这一点）。设ωb为图像ω，其中hJ的颜色改变为黑色（如果需要），ωw为图像ω，其中hJ的颜色改变为白色（如果需要）。让ΣT0（ω）=[α（γ{I（ω）−I（ω）}−ln（σh（ωw））HJH W HBh∈C（hJ）σh（ωb）{σ2（ωb）µ2（ωw）−σ2（ωw）µ2（ωb）}-2σ2（ω）σ2（ω））]，hbhw其中α是一个大的正整数，使得T0（ω）可以存储为所有团配置的4字节整数（它用于提高值的准确性），[]表示舍入到最接近的整数。 对于所有24元邻域配置，我们将T0（ω）的值存储在查找表中。查找表有224个4字节整数类型的条目，因此它需要只 有26兆 字 节。 由 于 值T0（ω ）不取决于像素hJ的颜色，请注意，在Metropolis算法的迭代步骤中，我们将总是有T0（ω）=T0（ω）。hJ1hJ2−瓦尔迪-赫尔曼公司386−2h1 1h21−122βββββ≤≤M（ω）HWHH我们再定义两个变量，σ2（ω）µ（ω）σ2（ω）µ（ω）T（ω）=[αh h]，hσ2（ω1）σ2（ω2）H Hσ2（ω）σ2（ω）T（ω）=[α hh]。h2σ2（ω1）σ2（ω2）H H这些值仅取决于th（ω1）、th（ω2）的值（分别为图像ω1、ω2中像素hJ的我们存储T1（ω）和T2（ω）在两个附加的查找表中。这些桌子相对来说H H小;它们有6× 6个条目，每个条目都是4字节整数类型。注意，T0（ω）、T1（ω）和T2（ω）的值不依赖于g（h）的值，1≤h≤H。创建随机数文件：我们创建以下随机数表：一个表用于行拾取，一个表用于列拾取，以及一个表用于退火中的每个β值包含整数类型的明细表αlnx，其中x是随机数均匀分布在（0， 1）范围这些文件的大小将在第4节中进一步讨论。算术运算：g（h）和g2（h）的值因此，在分割程序开始时计算一次要分割的灰度级图像的整数值[g（h）]和[g（h）]2，并存储以备将来使用。我们将上面介绍的表值以及[g（h）]和[g（h）]2的值保存为整数，以便分割程序中的所有算法都是精确的。这样，不会出现舍入误差，也不需要重新校准。减少计算：我们感兴趣的是找到使函数M（ω）最大化的二进制图像ω，而不是M（ω）的实际值令ωW为二进制白色图像（对于所有1h H，ωW（h）= 1）。注意，搜索使函数M（ω）最大化的二进制图像ω等价于搜索使函数M（ω）=M（ω）最大化的二进制图像ω，这等价于搜索使函数1lnM（ω）最大化的二进制图像ω·，即，ω=ω=ω·。因此，我们开始最大化函数1lnM（ω），它不依赖于β的值。分割算法用图像ωW初始化，其中1lnM（ωW）= 0。通过最大化1lnM（ω）而不是M（ω），我们不必计算M（ωW）的初始值瓦尔迪-赫尔曼公司387β∈JHβββββM（ω1）M（ω1）H3执行假设我们有一个训练集，实现的第一步是计算局部势，如第1节所述第二步是为给定的训练集找到合适的γ 设f是 二进制训练图像中前景像素（值为1的像素）的百分比。(Inf是给定训练图像中的平均骨体积分数。 我们运行上面描述的Metropolis算法的变体，以从分布Q（ω）产生随机二进制图像。我们跟踪这些随机图像中前景像素的百分比我们寻找这些百分比接近f的γ值。第三步是计算0≤t 6时的µt和σt值。这是这样做的：对于一个给定的团配置类型t0，找到所有二进制训练图像中的所有像素，其团配置是内部的，输入t0，从灰度训练图像中计算g（h）的相应值，并计算g（h）值的均值和标准差最后，我们讨论了如何使用第2节中描述的表和4字节随机形式的整数整数αlnx. 对于任何图像ω和像素hJ，我们定义两个Σ122ShJ（ω）={[g（h）]Th（ω）+[g（h）]h C（hJ）Th（ω）}，<$T0（ω）−ShJ（ω），如果ω（hJ）= 0，vhJ（ω）=0J−ThJ（ω）− S hJ（ω），如果ω（h）= 1。使用前面给出的定义，很容易检查在迭代步骤中关于Metropolis算法vhJ（ω1）n=αlnM（ω2）=αlnM（ω2） （等于一些舍入误差），而与HJ的颜色无关。因此，无论何时决定改变hJ的颜色，αlnM<$（ω2）的结果值可以通过将整数vhJ（ω1）与已知值αlnM< $（ω1）相加来估计。(For大的α，vhJ（ω）值是αlnMj（ω）值的差的一个很好的近似经过多次迭代，我们无法确定累积值与αlnM（ω）的实际值有多接近，这要求对所描述过程的有用性进行实验验证在这里。）此外，这是一个正确的方式来决定是否ω1应该在Metropolis算法中由yω2表示的算法是检查是否αβlnxβ ≤vhJ（ω1）。 实际上，由于vhJ（ω）是整数，这等价于lnx≤αvhJ（ω1），等价于x≤exp（βvJ（ω））n=M（ω2）=p，αM（ω1）概率是p因为x在区间上是均匀分布的β瓦尔迪-赫尔曼公司388Σβββββ×≤×（0， 1）。在运行分割程序之前，我们创建第2在分割程序开始时，我们计算并存储[g（h）]和[g（h）]2的值。 我们从白色图像ωW开始算法。在每次迭代中，我们(i) 读取随机行和随机列，这对应于随机选择像素h，J，(ii) 查找相应的表值T0（ω）并计算ShJ（ω1），(iii) 计算vhJ（ω1）的值，hj1(iv) 读取一个随机整数αlnx得双曲正弦值.Σ Σ(v) 如果 αlnx≤vhJ（ω1），则改变ω1中hJ的颜色，并将vhJ（ω1）添加到αlnM< $（ω1），以获得αlnM<$（ω2）的值β β该算法执行120，000个循环，其中每个循环是如上所述的H次迭代我们从迄今为止具有最高αlnM（ω）值的图像开始10，000个周期中的每一个（除了从白色图像开始的前10，000个输出是具有获得αlnM（ω）（这是其中当α>0时，得到1lnM（ω4实验作为初步实验，我们用高分辨率扫描了一张128 × 128的骨小梁核磁共振（NMR）灰度图像（取自[2]）;见图1中的左图 使用Adobe Photoshop版本5中的工具。5（Adobe System Inc.，San Jose，CA）和图像处理工具包版本3。0（Reindeer Games Inc.，Raleigh，NC），我们创建了请参见图1中的右图在这个实验中，训练集只包括一个灰度图像和它对应的二值分割图像在这里介绍。从我们的训练集中，我们首先计算局部势的值（9。三十七，七。三十二，四。12，3。二十六，三。31，0。11），对于0 t6。 百分比 二值训练图像中的前景像素（具有值1的像素）的值为f = 17。百分之六我们搜索γ的值（到两个有效数字），对于该值，来自分布Q（ω）的随机图像中的前景像素的平均数尽可能接近f。 在这个实验中，这是由值γ = 0来填充的。083.我们需要大约20亿行/列/整数，每一个在120000个周期的一个分割的128128骨小梁图像。我们声称一个包含5000万个随机行/列/整数的文件足够大，如果瓦尔迪-赫尔曼公司389×β−−β≤ββFig. 1. 左边的图像是从印刷品上扫描得到的骨小梁的NMR灰度图像（128 × 128像素）。 (In核磁共振图像， 骨头没有信号，因此是黑色的。）右边的图像是用于读取的随机起始位置16。每1000个周期从每个文件中选择400万个数字我们创建了以下随机数表，每个表包含5000万个随机数：一个表包含1字节无符号字符用于行拾取，一个表包含1字节无符号字符用于列拾取，一个表包含4字节无符号字符用于列拾取。铅字格αlnx对于退火时间表中的每个β值，如所描述的以上在这个实验中，我们选择α=4，500，000。 我们使用以下简单的退火时间表：从β = 0开始。1，将β的值增加0。每10，000次循环1次，总共120，000次循环。（x的最小值是1。在我们的情况下，RAND MAX =2311，RAND−MAXβ值= 0。1，αlnx约为-1，000， 000， 000，可以是存储为4字节整数。）我们发现µt的值（112。27，187。45，164。六十四，一百四十五。07，174。93，166。第80页）和σ t（32. 07，5。六十，十六。06，20。68，12。85，17。60），分别为0t 6，从灰度图像和其对应的二值图像。我们创建第2节中描述的三个查找表我们的分割算法的输出示例见图2中的左图。 一共有8个。6pt像素在两个分割数据之间存在差异这里显示的图像和二进制训练图像。的价值1lnM（ω）（其不取决于β的值）为-75，737，而对于分割的二进制图像，它是-73，124请注意，此处提供的1lnM（ω）比二进制训练图像的要高。这并不确定--由于二值化训练图像有23个。6%的前景像素在二进制训练图像，相比之下，17。6%的前景像素在我们的分割输出很明显，我们算法的输出包含许多小的（连通的）瓦尔迪-赫尔曼公司390图二.左边的图像是我们的分割算法的输出图像的示例。右边的图像是通过使用阈值对左边的图像进行后处理而获得的：threshold-black=30和threshold-white=60（见正文）。不属于骨小梁的组成部分（骨小梁主要是一种连接结构）。这些成分似乎来自我们试图分割的噪声灰度图像为了删除这些组件，需要进行一些后处理5后处理从我们的分割算法输出的图像包含需要被移除的小的白色和黑色连接分量（颜色需要被切换）。设threshold-white为白色像素的最大8连通分量中将被移除的像素数（颜色将变为黑色）。设threshold-black为黑色像素的最大4连通分量中将被移除的像素数（颜色将变为白色）。假设连通分量的大小是分量中的像素数我们做了以下实验。我们找到并删除所有大小小于或等于threshold-white的白色像素的8-连通分量接下来，我们找到并删除大小小于或等于threshold-black的黑色像素的所有4-连通分量。我们得出的结论是，在分割算法期间以上述方式去除成分是不必要的;从输出分割中去除小成分是足够的得出这一结论后，我们的任务是为threshold-white和threshold-black找到好的值。我们发现，在这个特定的应用中，使用threshold-white=60和threshold-black=30给出了良好的结果。使用较大的阈值会导致丢失部分骨骼，而使用较小的阈值会导致二值图像中出现一些不需要的小对于这些阈值，后处理的分割图像和金标准分割之间的不同像素的百分比为4。6，相比之下，8。未后处理的分割图像和金标准图像之间的6%不同像素瓦尔迪-赫尔曼公司391图三.这些是离体拍摄的骨小梁的二维图像，其被处理为类似于体内图像。参见图2。一共有19个。7%的前景像素在后处理分割图像（相比17. 6%的训练图像和23。6%在未后处理的分割图像中）。6讨论在这里报告的初步结果的鼓舞下，我们计划分割骨小梁的体内NMR图像在图3中，我们显示了几个离体拍摄的松质骨的二维图像，这些图像经过处理以类似于体内图像。为了成功分割这些图像，有必要选择适当的团配置类型。很明显，在第4节中描述的实验中使用的配置类型对于这些包含更小细节的图像来说不是一个好的选择。必须找到一个合适的后验信息模型，经过初步的实验，在这种情况下，与本文提出的后验信息不同的后验信息类型似乎更有用最后，一些实验是必要的，以确定什么是一个很好的退火时间表，这些分割。7结论吉布斯先验似乎是一个潜在的有用的工具，在嘈杂的图像中的复杂对象的分割引用[1] 卡瓦略湾M.，G. T.赫尔曼，S。马泰角Salzberg和E. 三平面心血管造影的二值断层扫描，医学影像信息处理，A。Kuba，M.Samal，A.Todd-Pokropek（eds.），柏林斯普林格（1999），29-41.瓦尔迪-赫尔曼公司392[2] 黄氏S. N.，F. W. Wehrli和J.L. Williams，3D NMR图像中基于概率的结构参数作为骨小梁强度的预测因子，Med.Phys.24（1997），1255[3] 廖氏H.是的，和G. T. Herman，用于二元层析成像的Gibbs先验参数的自动估计，理论计算机科学电子笔记46（2001），网址：http://www.elsevier.nl/locate/entcs/volume46.html。[4] 纽 约 大 都 会 ， A. W. Rosenbluth ， M. N. Rosenbluth ， A. H. Teller 和 E.Teller，用快速计算机计算状态方程，J.Chem. Phys.21（1953），1087[5] Vardi，E.，G. T. Herman和T. Y.孔，从有限投影方向加速二值图像的随机重建，线性代数及其应用：离散层析成像特刊（即将出版）。