三维图割算法结合随机游走算法的PET‐CT肺肿瘤分割方法与流程

本发明属于生物医学影像处理领域，利用高斯混合模型(Gaussian Mixture Models)优化的图割(Graph Cut)方法结合随机游走(Random walk)算法进行肺肿瘤的精确分割方法。

背景技术：

肺肿瘤大多数起源于支气管粘膜上皮，因此也称支气管肺癌。近 50年来，全世界肺癌的发病率明显增高，据统计，在欧美某些国家和我国大城市中，肺肿瘤的发病率已居男性各种肿瘤的首位。对于肺肿瘤的治疗需要对肿瘤进行定位、大小形状分析，而PET(Positron Emission Computed Tomography)和CT(Computed Tomography)作为两种定量的分子‐结构成像技术，已广泛使用于肺肿瘤的分析与诊断中， PET图像能够提供人体新陈代谢的功能信息，CT图像能够提供人体组织器官的结构信息，两者的结合能够更好地定位诊断肺肿瘤。

然而目前国内外很多现有肺肿瘤方法都是基于单一模态(PET或者CT)，例如：区域生长方法、Level Set方法等，其分割的准确性和可靠性都不是很高，并不能真正给临床应用带来质的改变。

技术实现要素：

本发明所要解决的技术问题在于针对上述背景技术中存在的问题，提供一种利用了随机游走算法的鲁棒性和图割算法的精确性使分割结果更加准确鲁棒的三维图割算法结合随机游走算法的PET‐CT肺肿瘤分割方法。

为解决上述技术问题，本发明采用的技术方案为：三维图割算法结合随机游走算法的PET‐CT肺肿瘤分割方法，包括如下步骤：

1)对原始PET图像进行线性上采样操作，并对PET和CT图像进行仿射配准；并手动标定图像的肿瘤种子点和非肿瘤种子点；

2)结合肿瘤种子点和非肿瘤种子点对PET图像进行随机游走算法分割，得到目标肺肿瘤区域的分割结果；

3)对于步骤2)中的随机游走分割结果作形态学膨胀操作，得到完全包括目标肺肿瘤区域的前景目标区域Ro，除Ro以外的区域作为非肺肿瘤区域的背景区域Rb；

4)将前景区域Ro和背景区域Rb分别建立高斯混合模型，并用最大期望算法估计完善高斯混合模型参数；

5)根据前景及背景的高斯混合模型分别计算能量项并运用图割算法得到最终分割结果。

前述的三维图割算法结合随机游走算法的PET‐CT肺肿瘤分割方法，步骤2)中，随机游走算法基于图的分割方法，在给定一部分种子点后能计算图像中非种子点的像素点首次到达种子点的概率，通过将概率最大的一类像素点归为目标区域从而得到最终的随机游走分割结果；其中随机游走算法中相邻像素点之间的权重表示为： Wij＝exp(-β(gi-gj)²),

其中β为权重参数，gi和gj表示邻接像素点i和j的灰度值；求解非种子点的像素点首次到达种子点的概率等同于求解组合狄利克雷问题，离散狄利克雷积分可以表示为：

其中x表示每个非种子点的像素点到达种子点的概率，L是组合拉普拉斯矩阵，xi为像素点i到达种子点的概率，xj为像素点j到达种子点的概率，eij表示邻接像素点i、j之间的一条边，x^T表示每个非种子点的像素点到达肿瘤种子点和非肿瘤的种子点的各自概率组成的矩阵的转置，wij为邻接点i、j之间组成的边eij的权重，定义如下：

其中Lij是PET图修改过的邻接矩阵，di＝∑w(eij)是顶点的度； w(eij)表示邻接点i、j之间组成的边eij的权重，∑w(eij)是指与当前像素点相邻接的边的权重和；adjacent表示像素点i、j是相邻接的， otherwise表示像素点i、j不相邻接且像素点i、j不是同一像素点；通过结合肿瘤种子点和非肿瘤种子点对PET图像进行随机游走算法分割，得到目标肺肿瘤区域的分割结果。

在步骤2)中，已被标记的肿瘤种子点及非肿瘤种子点和未被标记的非种子点像素，已标记点的集合(即种子点，包括肿瘤种子点和非肿瘤种子点的像素点集合)用VM表示，未标记点的集合(即不是种子点的像素点集合)用VU表示，VM和VU满足VM∩VU＝0，VM∪VU＝V，因此

其中xM和xU分别表示标记像素点和未标记像素点的概率，而xM^T和 xU^T则表示xM和xU的转置；是对组合拉普拉斯矩阵L根据标记点和非标记点分解的一个组合矩阵，其中B和B^T分别表示组合拉普拉斯矩阵L在标记点处的子拉普拉斯矩阵及其转置，LU表示组合拉普拉斯矩阵L在未标记点处的子拉普拉斯矩阵，LM表示组合拉普拉斯矩阵L在已标记点的拉普拉斯矩阵；求解上式对xU的微分，可以得到最终的概率，微分得：

通过结合用户输入的肿瘤种子点和非肿瘤种子点对PET图像进行随机游走算法分割，也就是求解上式，得到各个未标记的像素点分别对于所标记的肿瘤种子点和非肿瘤种子点的概率大小，当某一像素点属于前景的概率大于0.5，则认为其属于前景，反之亦然；通过概率阈值0.5，最终得到目标肺肿瘤区域的分割结果。

步骤3)中，形态学膨胀操作是将目标前景像素点接触的一定相邻背景像素点合并到目标前景当中，使得前景边界向外部扩张。

前述的三维图割算法结合随机游走算法的PET‐CT肺肿瘤分割方法，步骤4)中，混合高斯模型定义如下：

其中k表示高斯混合模型中第k个高斯函数，K表示高斯混合模型中高斯函数总数，β代表每个高斯函数的权重是高斯密度函数，μk和分别为均值和方差：

并采用最大期望算法来估计完善高斯混合模型参数，首先通过 K‐means算法将前景区域和背景区域分别聚类成独立的五个区域，并且得到初始各个类的参数，包括均值、方差以及权重，作为各个高斯分量的初始参数值。再通过EM算法，迭代计算最终参数，使得前景高斯混合模型和背景高斯混合模型参数最优，最终得到PET图像的前景和背景混合高斯模型。

前述的三维图割算法结合随机游走算法的PET‐CT肺肿瘤分割方法，步骤5)中，图割算法通过对能量的最小化得到全局最优解，将 PET和CT图像统一建图，并根据PET和CT所提供的信息不同设置上下文代价能量项Econtext，其中对于PET图像的能量项设置为：

其中三个能量项分别为高斯混合模型能量项，梯度下降能量项和边界代价能量项，u、v为PET 图像中的像素点，GPET是PET图像中邻域点的集合，NPET是PET图像中邻域点的集合，f为每个像素点的标记，α、β、γ分别对应每个能量项的对应权重系数，fu为PET图像中u像素点的标记，fv为PET图像中v像素点的标记；

CT图像的能量项设置为：

相应的其中三项能量项分别为高斯混合模型能量项，形状约束能量项和边界代价能量项。GCT是 CT图像中邻域点的集合，NCT是CT图像中邻域点的集合，f为每个像素点的标记，α、β、γ分别对应每个能量项的对应权重系数，fu为 CT图像中u像素点的标记，fv为CT图像中v像素点的标记；

上下文代价能量项设置为：

Econtext(f)＝p(1-|Nu-Nu'|+q),

其中p是缩放因子，NU和NU'分别是PET和CT区域项的归一化， q是PET和CT判定相反时的惩罚值。

通过上述的PET和CT图的能量项，以及PET和CT判定相反时的惩罚值，运用最大流/最小割算法实现图割算法，最终得到PET和CT 图的肺肿瘤分割结果。

相对于现有技术，本发明具有以下有益效果：

本发明首次提出了运用高斯混合模型来优化图割算法，并且结合随机游走算法来分割肺肿瘤。本方法充分运用了PET图像的功能信息和CT图像的结构信息，使得随机游走算法和图割算法进行互补，最终得到更加精确的肺肿瘤分割结果。

附图说明

图1图割算法PET‐CT联合建图；

图2(a)PET图像(b)CT图像；

图3为本发明方法与传统分割方法的DSC系数对比图；

图4为本发明方法在PET和CT上的部分切片分割结果对比，红色为金标准，蓝色为此发明方法的分割结果；

图5为本发明的算法原理框图。

具体实施方式

下面结合附图对本发明作更进一步的说明。

本发明的肺肿瘤分割方法是基于PET‐CT两种模态的联合分割，意在充分合理运用两种模态的特征信息进行精准定位和分割。如图2 为肺部PET和CT的一个切片图像。具体分割方法如下：

1、对原始PET图像进行线性上采样操作，使得PET图像和CT图像具有相同的像素点，并对PET和CT图像进行仿射配准，最终使得 PET和CT图像的像素点一一对应。并通过人眼判断是否属于肿瘤或非肿瘤区域，手动标定图像的肿瘤种子点和非肿瘤种子点，同时在临床肿瘤专家的监督下确定肺肿瘤的金标准，即图像中确切的肺肿瘤区域，以便最后比较本发明分割方法的分割准确性。

2、随机游走算法(Random Walk)是一种基于图的分割方法，在给定一小部分种子点后能快速计算图像中非种子点的像素点首次到达种子点的概率，通过将概率最大的一类像素点归为目标区域从而得到最终的随机游走分割结果。其中随机游走中相邻像素点之间的权重表示为：

Wij＝exp(-β(gi-gj)²),

其中β为权重参数，gi和gj表示邻接像素点i和j的灰度值。数学上已经证明求解随机游走的概率等同于求解组合狄利克雷问题 (combinatorial Dirichlet problem)，离散狄利克雷积分可以表示为：

其中x表示每个非种子点的像素点到达种子点的概率，L是组合拉普拉斯矩阵，xi为像素点i到达种子点的概率，xj为像素点j到达种子点的概率，eij表示邻接像素点i、j之间的一条边，x^T表示每个非种子点(即没有被定义为肿瘤种子点或者非肿瘤种子点的像素点) 的像素点到达肿瘤种子点和非肿瘤的种子点的各自概率组成的矩阵的转置，wij为邻接点i、j之间组成的边eij的权重；非肿瘤种子点代表不是肿瘤而是正常组织的种子点，非种子点是除了肿瘤种子点和非肿瘤种子点以外的点；定义如下：

其中Lij是PET图的邻接矩阵，di＝∑w(eij)是顶点的度，w(eij)表示邻接点i、j之间组成的边eij的权重，∑w(eij)是指与当前像素点相邻接的边的权重和；adjacent、otherwise分别表示邻接点、不邻接点， adjacent表示像素点i、j是相邻接的，otherwise表示不属于前两种情况的情况，也就是不相邻接且像素点i、j不是同一像素点。

求解组合狄利克雷问题即是求解一个调和函数来满足这个狄利克雷积分的边界值，并且调和函数是最小化狄利克雷积分的函数，由于拉普拉斯矩阵L是半正定的，所以(1)式的最小值在它的极点位置。根据图像中的已标记点和未标记点，也就是已被标记的肿瘤种子点及非肿瘤种子点和未被标记的非种子点像素，已标记点(种子点)的集合用VM表示，未标记点(非种子点)的集合用VU表示，VM和VU满足 VM∩VU＝0，VM∪VU＝V，因此上式可以分解为：

其中xM和xU分别表示标记像素点和未标记像素点的概率，而xM^T和 xU^T则表示xM和xU的转置。是对组合拉普拉斯矩阵L根据种子点和非种子点分解的一个组合矩阵，其中B和B^T分别表示组合拉普拉斯矩阵L在种子点处的子拉普拉斯矩阵及其转置，LU表示组合拉普拉斯矩阵L在未标记点(非种子点)处的子拉普拉斯矩阵，LM表示组合拉普拉斯矩阵L在已标记点(种子点处)的拉普拉斯矩阵。所以求解上式对xU的微分，可以得到最终的概率，微分得：

通过结合用户输入的肿瘤种子点和非肿瘤种子点对PET图像进行随机游走算法分割(Random Walk)，也就是求解上式，得到各个未标记的像素点分别对于所标记的前景(肿瘤)和背景(非肿瘤)种子点的概率大小，本发明中规定像素点属于概率大于0.5的一类，即当某一像素点属于前景的概率大于0.5，则认为其属于前景，反之亦然。通过概率阈值0.5，最终得到目标肺肿瘤区域的分割结果。前景种子点即为肿瘤种子点，背景种子点即为非肿瘤种子点。

3、对于步骤2中的随机游走分割结果作形态学膨胀操作，形态学膨胀是将随机游走分割结果所表示的像素点接触的一定相邻像素点合并到分割结果当中，使得边界向外部扩张的一个过程，在本发明中规定相邻像素点为四个像素点。具体步骤为：用4x4的结构元素来依次覆盖扫描random walk分割结果图的各个分割结果像素点做形态学膨胀操作，并与其覆盖的random walk分割结果图的各个像素点做“与”操作，如果结果都为0，则判定该4x4的结构元素为0，该结构元素不加入分割前景区域；反之则判定为1，该结构元素加入分割前景区域。因此，通过膨胀处理会使分割结果扩大，以至于完全包括了肺肿瘤区域，标定这一扩大以后的区域为前景区域Ro，余下的称为背景区域Rb。

4、混合高斯模型是运用多个高斯密度函数来精确归类目标图像。在本发明中，我们将PET图像中的像素分布视为若干个高斯分布的混合重叠，因此将前景区域Ro和背景区域Rb分别建立高斯混合模型，混合高斯模型定义如下：

其中k表示高斯混合模型中高斯函数的个数，K表示高斯混合模型中高斯函数总数，y表示每个像素点的灰度值，β代表每个高斯函数的权重是高斯密度函数，μk和分别为均值和方差：

并采用最大期望算法(EM)来估计完善高斯混合模型参数，首先通过K‐means算法将前景区域和背景区域分别聚类成独立的五个区域，并且得到初始各个类的参数，包括均值、方差以及权重，作为各个高斯分量的初始参数值。再通过EM算法，迭代计算最终参数，使得前景高斯混合模型和背景高斯混合模型参数最优，最终得到PET 图像的前景和背景混合高斯模型。

5、本发明采用的图割算法(Graph Cut)是基于图论的一种分割方法，通过对其能量的最小化来得到全局最优解，本发明中将PET和 CT图像统一建图(如图1)，其中黄色像素点表示PET图像的各个像素点，红色像素点表示CT图像的像素点；S、T分别代表前景和背景；红线表示各个像素点属于前景或者背景的权重值；绿线表示相邻像素点之间的相似的权重值；蓝线表示PET和CT的上下文代价权重值(图 1只是示例，本发明实际分割的PET和CT图像像素点个数为120*120 个)。Graph cut即是寻找一些线并将其断开，得到两个像素点子集，分别对应图像的前景像素集和背景像素集，就相当于完成图像分割。并且需要使得分割的线的权重值总和最小，也就是最小割，来达到最佳分割，而福特‐福克森定理表明，网路的最大流max flow与最小割 min cut相等，所以由Boykov和Kolmogorov发明的max‐flow/min‐cut (最大流/最小割)算法就可以用来获得S‐T图的最小割。

其中对于PET图像的能量项设置为：

其中三个能量项分别为高斯混合模型能量项，梯度下降能量项和边界代价能量项，u、v 为PET图像中的像素点，GPET是PET图像中邻域点的集合，NPET是PET 图像中邻接边的集合，f为每个像素点的标记(或为前景或为背景，根据比对为前景所累加的E大小与为背景所累加的E大小最终来判定是前景还是背景)，α、β、γ分别对应每个能量项的对应权重系数， fu为PET图像中u像素点的标记，fv为PET图像中v像素点的标记。 CT图像的能量项设置为：

相应的其中三项能量项分别为高斯混合模型能量项，形状约束能量项和边界代价能量项。GCT是CT图像中邻域点的集合，NCT是CT图像中邻接边的集合，f为每个像素点的标记(或为前景或为背景，根据比对为前景所累加的E大小与为背景所累加的E大小最终来判定是前景还是背景)，α、β、γ分别对应每个能量项的对应权重系数，fu为CT图像中u像素点的标记， fv为CT图像中v像素点的标记。上下文代价能量项设置为：

Econtext(f)＝p(1-|Nu-Nu'|+q),

其中p是缩放因子，NU和NU'分别是PET和CT区域项的归一化， q是PET和CT判定相反时的惩罚值。

max-flow/min-cut算法是Graph cut的核心算法，它的解法可分为两类：增广路径法和推进重标记法。在本发明中使用增广路径法解max-flow/min-cut算法，它的基本思想是如果从源点到汇点所经过的各个边的流量没有饱和(这些边就被称为增广路径)，那么计算它们能允许通过的最大流，然后加入流量，使它们饱和，然后寻找下一条非饱和边，直到路径的流量饱和。具体实现方式如下：首先建立两颗搜索树S和T，S为源点，T为汇点。树S中的结点分为主动结点和被动结点，主动结点在树的前端，可以使搜索树不断地生长，扩冲结点，被动结点包围着主动结点，它不能继续生长。搜索树在生长过程中，对结点进行标记，然后不断地扩展，更新被标记的结点，形成新的搜索树，重复以上操作直到没有可增广路径。算法可分为三个阶段：1、生长阶段：搜索树从源点出发，在生长过程中遇到结点就将它包含在树中，直到找到汇点，停止生长。2、扩展阶段：若有非饱和边，则扩展路径，搜索树变成森林。3、收养阶段：收养孤立结点，恢复搜索树。通过上述max-f l ow/mi n-cut算法得到图像的饱和边，由最大流与最小割的等价性可知，从源点到汇点方向的饱和边为割边时所花费的代价函数总和最小，也就是使得分割最精确，由此便可实现graph cut算法。

通过上述给定的PET和CT图的能量项，也就是每个点对应前景和背景的权重值，以及惩罚分割PET和CT图不一致的权重值，运用 max-flow/min-cut算法实现graph cut算法，最终得到PET和CT图的肺肿瘤分割结果。

6、本发明所采用的肿瘤影像数据来自苏州大学附属第一医院影像科，其中选取20个数据作为分析此发明方法分割结果准确性和鲁棒性的比较，通过DSC系数(Dice similarity coefficient)来直观统计分割结果与金标准的重合度，即分割的准确性。结果分析如图3所示。其中所对比的传统分割方法是基于图割的传统方法，但缺少本发明所提出的高斯混合能量项以及随机游走分割结果做初始化的传统方法。其中DSC系数的定义为：

其中U1和U2分别表示此发明分割结果和金标准。图4是本发明方法在PET和CT上的部分切片分割结果和金标准的对比示例图，图 5列出了本发明算法的原理框图。

本发明技术关键点包括：(1)将高斯混合模型融合到图割算法的能量项中，以提高肿瘤分割的准确性。(2)利用随机游走算法的鲁棒性对肿瘤进行预分割得到初步分割结果作为图割算法的种子点。

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。