一种基于图谱字典学习的图像分割方法与流程

本发明属于图像分割技术领域，具体涉及一种基于图谱字典学习的图像分割方法。

背景技术：

正电子发射断层扫描技术(positronemissiontomography，pet)很适合用于测量动物的代谢强度，但是基于起分辨率受限等原因，因此必须对测量的数据进行适当地处理，字典学习算法在近几年收到许多学者的关注，并不断发展成熟，成为了医学图像数据处理的有力工具。

从原理上来看，pet扫描器获得的原始数据是一系列由探测器获得，由正子与电子湮灭产生的一对光子的并发事件。每个并发事件背后，有一个正电子逸出，从而引发一个湮灭事件，在空间中同时射出背向的两个光子并被捕捉到。pet成像的分辨率主要受限于以下四个因素：(1)正电子在湮灭和发射γ射线之间有一段小距离；(2)两束γ射线之间并不是沿着理想情况下的180°的角度发射；(3)γ射线在可能晶体中任意厚度处发生反应；(4)相机的失真。其中(2)是分辨率损失的主要来源，目前pet图像的分辨率普遍在4～7mm。信噪比是医学图像中的一个重要参数，这取决于成像系统的分辨率，因为其分辨率的限制，为了在较低信噪比的情况下分割出有用信息，近年来人们先后提出了多种分割算法，其中字典学习的分割算法因其快速和高精度而倍受学者们关注。

基于字典学习图像分割算法(如mod)，面临着对运算速度要求高精确度低等问题，一个原因是这些字典学习的图像分割方法只基于图像像素或体素的强度(灰度)，没有考虑到其位置所能代表的信息，通过使用基于图集的方法可以增加先验判断来把像素的位置信息也考虑进去，把体素的位置信息体现到图集中，并在每次更新字典时也更新图集。

技术实现要素：

鉴于上述，本发明提供了一种基于图谱字典学习的图像分割方法，其通过位置信息更新先验统计概率，在分割过程中保留了图片像素的位置信息，结合像素强度信息，可以更加有效地分割图片。

一种基于图谱字典学习的图像分割方法，包括如下步骤：

(1)获取一定数量已精确分割的图像，利用em算法(expectationmaximizationalgorithm，期望最大化算法)对这些图像重新进行分割并根据图谱信息以列向量的形式建立矩阵a1，同样以列向量的形式根据这些图像精确分割后的结果建立矩阵a2；

(2)对于待分割的n幅图像，利用em算法对这些图像进行分割并根据图谱信息以列向量的形式建立矩阵b，n为大于0的自然数；

(3)对矩阵a1按列进行降采样得到矩阵d1，根据矩阵a1的降采样结果对矩阵a2中对应的列进行取样得到矩阵d2；

(4)利用图谱字典d和稀疏系数矩阵x来表示矩阵b，建立如下目标函数：

其中：λ为权重系数，稀疏系数矩阵x的维度为m×n，m为步骤(1)中获取到已精确分割的图像数量；

(5)使矩阵d1作为图谱字典d的初值，对所述目标函数进行最优化求解得到稀疏系数矩阵x，使d2×x并经后期聚类处理即可得到这n幅图像的分割结果。

所述步骤(1)和步骤(2)中利用em算法对图像重新进行分割，该em算法中先验概率的计算表达如下：

其中：对于任一类别r以及任一像素i，pr(i)表示图像中像素i属于类别r的先验概率，z为像素i周围像素的集合，j为集合z中的任一像素，s为集合z中的像素数量，labelr(j)为像素j的分类标签，若像素j属于类别r则labelr(j)＝1，否则labelr(j)＝0。

所述矩阵a1由m列图谱向量组成，每一列即对应一幅图像，矩阵a1的图谱向量由一条图像灰度向量以及k条分类灰度向量依次级联组成，k为设定图像分割后的类别数量；所述图像灰度向量中各元素值依次对应图像中各像素的灰度值，第k条分类灰度向量中各元素值依次对应图像中各像素的类别灰度值，即其中某一像素属于第k类则该像素的类别灰度值为255，否则该像素的类别灰度值为0，k为自然数且1≤k≤k。

所述矩阵a2由m列图谱向量组成，每一列即对应一幅图像，矩阵a2的图谱向量中各元素值依次对应图像中各像素的类别灰度值，所述类别灰度值对应有k种取值，k为设定图像分割后的类别数量，若其中某一像素属于第k类则该像素的类别灰度值为255(k-1)/(k-1)，k为自然数且1≤k≤k。

所述矩阵b由n列图谱向量组成，每一列即对应一幅待分割的图像，矩阵b的图谱向量由一条图像灰度向量以及k条分类灰度向量依次级联组成，k为设定图像分割后的类别数量；所述图像灰度向量中各元素值依次对应图像中各像素的灰度值，第k条分类灰度向量中各元素值依次对应图像中各像素的类别灰度值，即其中某一像素属于第k类则该像素的类别灰度值为255，否则该像素的类别灰度值为0，k为自然数且1≤k≤k。

所述步骤(5)中采样k-svd算法对目标函数进行最优化求解。

本发明基于图谱的字典学习方法进行图像分割，相对于传统的字典学习方法，本方法可以通过对图谱进行一定的处理来达到需要的目标，图谱的获得采用基于马尔可夫场考虑的em算法，通过字典学习获得图谱与所需分割图像的关系，目标函数与传统的字典学习方法一样，但是字典的初始化来自于图谱的取样，字典的每一列都是一幅图谱，采用的是k-svd的方法来进行字典的更新，可以提高运行速度以及避免字典陷入局部最优。

附图说明

图1为本发明图像分割方法的流程示意图。

图2为结构模板的示意图像。

图3(a)加入了高斯噪声的结构模板示意图像。

图3(b)为加入了瑞利噪声的结构模板示意图像。

图3(c)为加入了边缘模糊的结构模板示意图像。

图4(a)为采用本发明方法对图3(a)分割的结果示意图像。

图4(b)为采用本发明方法对图3(b)分割的结果示意图像。

图4(c)为采用本发明方法对图3(c)分割的结果示意图像。

具体实施方式

为了更为具体地描述本发明，下面结合附图及具体实施方式对本发明的技术方案进行详细说明。

如图1所示，本发明基于图谱字典学习的图像分割方法，包括如下步骤：

(1)根据已有的精确分割图像图片像素点强度以及其位置关系建立图谱矩阵a1及其对应的精确分割矩阵a2。

图谱是一个大小为n+1的元组{a0～an}，其中a0是扫描得到的原图像，a1～an-1是分类的先验值，分别代表n-1种分类，其实际的意义是一幅根据先验分割得到的仅包含此类像素点的图像；an为已有的专业分割图像，每类的像素有不同的标签，是a0的精确标记。

首先，实验中的每个a0(原图像)为一幅可以表示为一个一维向量，向量的长度d为图片像素点的数目，向量元素的大小表示图谱像素点的灰度值，则记此向量为y＝{y1,y2,…yd}对于y的每个元素，即原图的像素，都有一个对应它的分类，将此分类向量记为z＝{z1,z2,…zd}，zi一共有n-1种取值：e1～en-1，分别代表一种分类。

对每个点都有一个概率函数：

f(zi＝ek)＝p(k)

假设同一类的组织中像素强度是服从高斯分布的，则有下式：

f(yi|zi＝ek)＝gσ(yi-uk)

其中：σ和u分别为高斯分布参数，前者为标准差，后者为均值。

为了获得a1～an-1即各组织的分割，本实施方式采用基于马尔可夫从(mrf)的em算法，更新每一类高斯分布的两个参数以及每个体素的分类情况。em算法分为两步，估计每一类点按照高斯分布的参数和估计所有点所属的分类；然后依据其中一步的结果修正另一步的结果，多次循环直到参数变化极小。在估计体素所属某一类的时候，需要用到先验概率乘以每一类高斯分布得到，取概率最大的哪一类，即：

p(zi＝ek)＝pprior(zi＝ek)×f(yi|zi＝ek)

其中的pprior为先验概率，根据先验知识，可以认为一个体素很大概率上与其相邻体素属于同一类，本实施方式为了让体素的分类不仅仅基于体素强度，还基于周围体素所属的分类，采用了马尔可夫场来获得先验概率，将一个点的周围点所属于的类的作为先验概率的来源，具体表达如下：

其中：pprior(i)为属于第i个类的先验概率，s为其周围的点数量，j为其标号，labelj(i)的取值为0或1，若像素点j属于第i类，则取值为1，否则取值为0。

下一步是初始化像素点的标签，通过聚类的方法实现，仅通过像素强度将像素打上初始化标签，此标签会在后续的步骤中更新。

首先对每一类，求属于该类的体素强度的均值作为此类高斯分布的均值参数，求得它们的方差作为此类高斯分布的方差参数，由此可以得到n-1个类中体素的高斯分布函数，表达式如下：

f(yi|zi＝ek)＝gsigma(yi-uk)

根据此函数以及对应的先验概率相乘，得到所有体素属于每类的各自的后验概率，然后选择后验概率最大的一类作为此体素所属的类，代替原来用聚类得到的分类，然后重复上述步骤，直到循环次数达到预定的上限，或者每类的高斯分布所有的参数的变化率都小于预定的阈值。

对于已有的a1～an-1，将其写成向量，表达式如下：

f(i)＝[a1(i)^t,a2(i)^t,…an-1(i)^t]^t

将m个这样的列向量组合即可以形成一个二维矩阵，即字典，其形式为nd×m，记为a1：

a1＝[f(1),f(2)…,f(m)]

同样地定义a2为d×m的矩阵，其为m个an形成的字典。

a2＝[an(1),an(2)…an(m)]

(2)建立相应的待分割图像的图谱向量b(j)。

与{a1～an}类似，我们定义一个{s1～sn}的元组，其中s1代表我们需要分割的图像，也是以向量的形式表示，s2～sn则分别表示n-1个类，首先将a1做放射变化配准到s1，然后再将已有的a2～an做相同的变化，得到s2～sn。将每个需要分割的图像中的s1～sn类似地定义一个向量，表达式如下：

b(j)＝[s1(j)^t,s2(j)^t,…,sn(j)^t]

(3)稀疏表达目标函数。

从a1中抽样若干列形成字典d1，取a2相应地列形成字典d2，并用字典d1和系数向量x稀疏表示b(j)，稀疏表达目标函数的表达式如下：

(4)固定步骤(3)得到的系数向量x，更新字典d1，使误差减少，字典更新的目标函数表达式如下：

其中：表示最优的更新字典。

交替重复步骤(3)和(4)，直到误差达到预期的阈值；通过omp算法得到最优的系数求解，交替循环时定义新的原子dk+1，表达式如下：

其中：γk为无法被前k个原子表示的向量。

通过k-svd方法进行字典更新，对于字典，每次只更新一列，目标函数如下：

其中：ek为全局误差，是仅不考虑第k列时的误差，才用奇异值分解求解。

(5)得到最终的分割结果。

将得到的最优系数x作用域字典d2，可以得到初步的分解结果，经过后期聚类处理，即可得到最终的分割结果，表达式如下：

b(j)＝d2×x(j)

其中：b(j)为第j幅图的分割，b(j)中的像素点并非完全只有n-1类，按照从0～255的强度将像素点分成n-1类，即可得到原图的n-1分割。

以下我们通过对模拟的图像数据进行实验从而验证本发明系统分割结果的准确性，具体方法为在结构图片上加噪声，图2为实验所用结构模板示意图，将不同的区域分为三个区域。实验运行环境为：8g内存，2.30ghz，64位操作系统，cpu为inteli5-53000u；在本次试验中，高斯噪声、瑞利噪声和边缘模糊三种不同的条件下进行实验，实验生成数据总量为784幅结构图以及其各自的三种测试图，实验过程为随机抽取其中50幅作为已知分割结果的图像，另外再随机抽取50幅作为测试图像。图3(a)～图3(c)分别为加入高斯噪声、瑞利噪声以及边缘模糊的图像示例，分割结果如图4所示，其中图4(a)为高斯噪声图像的分割结果示例，图4(b)为瑞利噪声图像的分割结果示例，图4(c)为边缘模糊图像的分割结果示例。

此验证实验的评价参数采用dice值，表征了两幅图像的相似度，其定义如下式：

其中：此处的x与y分别代表两幅图，交集符号表示两者强度值相同的像素数，绝对值符号表示集合中元素的个数，我们将50幅图分别求其dice系数，并求平均。

经过实验验证，对高斯噪声处理后dice系数为0.8927，瑞利噪声处理后dice系数为0.8732，对边缘模糊处理后dice系数为0.9234。从分割结果数据和分割示例图中可以看出，本发明对三种类型噪声的图像分割都取得了较好的效果，其中对瑞利噪声图像的分割性能稍差，对边缘模糊的分割性能较好。

上述对实施例的描述是为便于本技术领域的普通技术人员能理解和应用本发明。熟悉本领域技术的人员显然可以容易地对上述实施例做出各种修改，并把在此说明的一般原理应用到其他实施例中而不必经过创造性的劳动。因此，本发明不限于上述实施例，本领域技术人员根据本发明的揭示，对于本发明做出的改进和修改都应该在本发明的保护范围之内。