一种基于正则化图割的大脑组织分割方法与流程

本发明属于数字图像处理技术领域，涉及大脑磁共振图像处理方法，更为具体地说，是涉及一种基于正则化图割的脑组织分割方法。

背景技术：

磁共振成像(mri)已被广泛用于在临床应用和神经科学研究中检查人脑的解剖结构。与其他医学成像方式相比，mri具有高空间分辨率和良好软组织对比的优势，可以对不同类型的组织进行精细的区分。精确分割这些组织对于一些应用至关重要。

为了方便对mri图像进行处理，引入超体素的概念。超体素(supervoxel)或超像素(superpixel)是近年来快速发展的一种图像预处理技术，它是指图像中局部的、具有一致性的、能够保持一定图像局部结构特征的子区域，相比于传统处理方法中的基本单元——像素，超体素更有利于局部特征的提取与结构信息的表达，并且大幅度降低后续处理的计算复杂度，在计算机视觉领域尤其是图像分割中得到了广泛应用。

当面对mri中噪声、偏置场等因素时，精准的分割变得困难。现在的主流方法主要有三类，水平集法、分类法和基于图谱的方法。当应用在自然图像中的水平集法被拓展到大脑组织分割时，其对用户的初始化设置和参数设置敏感。体素分割法能够实现优良的性能，但是，在mri中，面对数以万计的超体素时，就会遇到很高的计算复杂度问题，这限制了其在mri中的应用。

技术实现要素：

为解决上述问题，本发明公开了一种基于正则化图割的脑组织分割方法，首先基于强度距离和空间相似度，设计新的体素间相似度计算方法，从而对体素进行聚类，把大脑mri图像分割为一系列均匀并且较好地贴合图像边缘的超体素；随后明通过把大脑不同组织的先验概率融入到图割框架中，设计一个能量计算公式，计算每个超体素在分配不同标签时各个部分的能量值，从而使用图割方法对超体素分割，把mri图像分割成不同的组织。

为了达到上述目的，本发明提供如下技术方案：

一种基于正则化图割的脑组织分割方法，包括如下步骤：

步骤1：计算体素点到一定区域内的各个种子点的相似度，聚类到相似度最高的种子点，然后重新计算聚类中心的强度和中心位置，通过数次的迭代，使得具有同一性质，相似度高的体素聚到一类，形成超体素；具体包括如下子步骤：

步骤1-1，初始化种子点，假设图像一共有n个体素，期望生成的超体素的个数为m，则每个超体素的大小为：

v＝n/m

其中m也表示初始的种子点的数量；

步骤1-2，每个超体素的边长为：

上式中，l即为每个超体素的步长；

步骤1-3，计算每个种子点周围l×l×l区域内体素的灰度值的平均值，作为种子点的灰度值i，计算公式如下：

其中(xi,yi,zi)代表种子点i的坐标，ix,y,z是(x,y,z)坐标下的像素值；

步骤1-4，计算体素与种子点间的相似度，对于每个体素，在其2l×2l×2l的区域内，计算其与各个种子体素的距离，计算公式如下：

dint＝||g*ini-ic||

d＝dint+γdspa

其中g代表标准的高斯内核，*代表卷积操作；iint代表以体素ni为中心，立体图像块的强度(灰度)矩阵；ic代表种子点的强度；dspa代表体素和种子之间的空间相似度；dint代表以体素ni为中心的图像块与种子的强度距离；γ代表正则化系数，代表体素强度和空间相似度之间的相关度；

步骤1-5，通过步骤1-4中的公式，计算区域内各个体素点到种子点的相似度，并把体素点归到相似度最高的种子点；

步骤1-6，对于新生成的聚类中心，计算其范围内体素的平均强度和空间位置中心，作为新中心的强度和位置；

步骤1-7，重复上述步骤1-5，1-6，当新的中心点和上一次中心点的距离小于一定的阈值时，停止迭代，完成超体素的生成；

步骤2，融入组织的先验概率，加上从第一步的超体素中提取出的特征，使用图割框架对超体素聚类，具体包括如下子步骤：

步骤2-1，生成三种组织的先验概率分布图；

步骤2-2，通过下式计算每种标签的能量公式

e(c)＝μ(λrs(cs)+(1-λ)*atlass(cs))+∑bs,t

其中，s为超体素，cs为标签，μ和λ代表公式总数据项、概率图项、成对项之间的相关程度；

rs(cs):描述超体素s被分到标签cs时的好坏程度，计算方式如下：

rs(cs)＝-lnpa(c|s)

其中pa(c|s)可通过每种组织分布的高斯模型来得到；

atlas(cs)：描述每个超体素的组织概率，计算公式如下：

atlass(cs)＝-lnpa(cs)

bs,t：根据定义计算不同标签的超体素之间的权重

bs,t＝1/(1+ρ)

其中特征f代表超体素的平均强度，σ代表超体素间的方差，ρ为洛伦兹误差；

步骤3，根据步骤2中计算得到的能量值，使用图割来对其进行处理，得到能量公式的最小值，从而完成超体素的精准分割。

进一步的，γ取值为0.2。

进一步的，所述步骤2-1中生成先验概率分布图的过程为：使用国际脑图测绘组织概率图集，通过模板图像到给定图像的图谱配准工具elastix来计算得到先验概率图，表示每个位置体素属于每种组织的先验概率。

与现有技术相比，本发明具有如下优点和有益效果：

本发明首先把大脑mri图像分割为一系列均匀并且较好地贴合图像边缘的超体素，并利用已有的脑组织先验知识与超体素相结合，使用图割来对脑部mri图像分割、聚类，能够从最初的脑部mri分割出三种脑组织，分割结果中各个组织间边界贴合度高。与已有的mri图像分割方法相比，本发明分割效果更好，边界贴合度更高，效率更高，处理速度更快，可以较好地抑制噪声的影响。

附图说明

图1为本发明提供的基于正则化图割的大脑组织分割方法步骤流程图。

图2为传统k-means方法和本发明提出的方法中改进的k-means方法的计算范围示意图，其中(a)为传统k-means方法，(b)为改进的k-means方法。

图3为大脑mri图像中三种组织的概率图谱，其中(a)为wm白质图像，(b)为gm灰质图像，(c)为cfs脑脊液图像。

具体实施方式

以下将结合具体实施例对本发明提供的技术方案进行详细说明，应理解下述具体实施方式仅用于说明本发明而不用于限制本发明的范围。

本发明提供的一种基于正则化图割的大脑组织分割方法，其过程主要分为两部分：首先，使用基于强度距离和空间相似度的超体素生成方法，将体素划分为规则的，降低噪声影响，能较好贴合图像边缘地区的超体素；然后，结合脑组织的先验知识，使用图割，将各个超体素切割为特定的脑组织。本发明的具体流程如图1所示，包括如下步骤：

步骤1：超体素的生成：

本发明使用一种基于强度距离和空间相似度的超体素生成方法，设计一种新的相似度测量方式，计算体素点到一定区域内的各个种子点的相似度，聚类到相似度最高的种子点，然后重新计算聚类中心的强度和中心位置，通过数次的迭代，使得具有同一性质，相似度高的体素聚到一类，形成超体素，具体包括以下步骤：

步骤1-1，初始化种子点，假设图像一共有n个体素，期望生成的超体素的个数为m，则每个超体素的大小为：

v＝n/m

其中m也表示初始的种子点(即中心点)的数量；

步骤1-2，每个超体素的边长为

上式中，l即为每一个超体素的边长；

步骤1-3，为了降低噪声的影响，计算每个种子点周围l×l×l区域内体素的灰度值的平均值，作为种子点的灰度值i，计算公式如下：

其中(xi,yi,zi)代表种子点i的坐标，ix,y,z是(x,y,z)坐标下的像素值。

步骤1-4，对于每个体素，在其2l×2l×2l的区域内，计算其与所有种子点的相似程度，计算公式如下：

dint＝||g*ini-ic||

d＝dint+γdspa

其中g代表标准的高斯内核，其适用于一定的噪声等级，通过绝对中值差来评估，*代表卷积操作；

ini代表以体素ni为中心，立体图像块的强度(灰度)矩阵；

ic代表种子点的强度；

dspa代表体素点和种子之间的空间相似度，即体素与种子点之间的空间欧拉距离；

dint代表以体素ni为中心的图像块与种子的强度距离；

d为每个体素与种子点之间的距离，距离越近，相似度越高；

γ代表正则化参数，代表空间距离的权重，代表体素强度和空间相似度之间的相关度，经验值为0.2。

步骤1-5，采用步骤1-4中的公式，计算区域内各个体素点到种子点的相似度，并把体素点归到相似度最高的种子点，即生成聚类中心；

步骤1-6，对于新生成的聚类中心，计算其范围内体素的平均强度和空间位置中心，作为新中心的强度和位置；

步骤1-7，重复上述步骤1-5，1-6，当新的中心点和上一次中心点的距离小于一定的阈值时，停止迭代，完成超体素的生成。

传统的k-means算法是计算整幅图像中所有点到各种子点的相似度,如图2(a)所示，而本算法，对于一个体素，计算到其周围一块局部区域内各个种子点的相似度，如图2(b)所示，极大地减少了计算量，从而提高了算法的效率。由于每个超体素的初始大小为l×l×l，因此将局部区域设置为以体素为中心，大小为2l×2l×2l。

步骤2，基于图割的超体素分类：

本发明提供一种基于图割的超体素分类方法，相比于以往的方法，根据从超体素中提取出的特征来分类，本方法加入了组织概率图的先验知识，从而克服了大脑磁共振图中噪声的影响。本步骤中融入组织的先验概率，加上从第一步的超体素中提取出的特征，使用图割框架对超体素聚类，从而实现脑组织的分割。

步骤2-1，生成三种组织的概率分布图，使用国际脑图测绘组织概率图集，通过模板图像到给定图像的图谱配准工具elastix来计算得到先验概率图，表示每个位置超体素属于每种组织的先验概率，脑脊液图像(cfs)、灰质图像(gm)、白质图像(wm)概率图谱。

步骤2-2，计算每种标签的能量公式

e(c)＝μ(λrs(cs)+(1-λ)*atlass(cs))+∑bs,t

公式中包含三项：数据项、概率图项、成对项，μ和λ代表三项之间的相关程度。计算超体素s分配不同的标签cs时的能量公式，从而使用图割框架来切割，使得上式取得最小值，从而完成超体素的分割。

rs(cs):描述超体素s被分到标签cs时的好坏程度，计算方式如下：

rs(cs)＝-lnpa(c|s)

其中pa(c|s)可通过每种组织分布的高斯模型来得到。

atlas(cs)：描述每个超体素的组织概率

atlass(cs)＝-lnpa(cs)

使用步骤2-1中的三种组织的概率分布图，在每种组织的概率图谱中，计算每个超体素范围内各个体素的平均概率，作为该超体素属于该种组织的概率pa(cs),再通过如上的公式计算得到超体素的组织概率；

rs(cs)和atlas(cs)所使用的概率模型，使用不同的概率模板；

bs,t：根据定义计算不同标签的超体素之间的权重，其公式如下：

bs,t＝1/(1+ρ)

其中，s,t代表不同的超体素，特征fs,ft代表超体素s,t的平均强度，σ代表超体素间的方差；计算不同标签的超体素之间的洛伦兹误差ρ(s,t),接着做第二个公式的计算，得到两个不同标签的超体素间的相似度项；对所有不同标签的超体素计算此项，并求和，即得到上面能量公式中的第三项。

步骤3，根据步骤2中计算得到的能量值，使用图割来对步骤1中生成的超体素进行处理，得到能量公式的最小值，从而完成超体素的精准分割。

使用两大数据集brainweb18和internetbrainsegmentationrepository18，对本方法和其他方法比较。

表1详细地列出了四种方法在两大数据集上分割三种组织的结果以及用时，表中显示，本方法在三种组织上的dice系数(dicesimilaritycoefficient)值要高于其他三种方法，而且用时最少。

表1

本发明方案所公开的技术手段不仅限于上述实施方式所公开的技术手段，还包括由以上技术特征任意组合所组成的技术方案。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也视为本发明的保护范围。