一种基于随机游走和图割的活动轮廓模型的MR图像三维交互分割方法与流程

本发明涉及图像分割
技术领域：
，特别是对磁共振成像的人脑图像进行区域分割的，基于随机游走和图割的活动轮廓模型的MR图像三维交互分割方法。
背景技术：
：磁共振成像(Magneticresonanceimaging，MRI)具有优良的软组织分辨力，可对人体各部位多角度、多平面成像，多种成像方法与脉冲序列技术的应用有助于对垂体瘤的定位、定性诊断。通常临床上医学影像的肿瘤分割会由医生在二维影像上逐序列手工描绘垂体瘤边界来完成。根据实体瘤疗效评价(RECIST)标准，放射医生会识别出二维影像上脑部的高亮度/低亮度组织，测量垂体瘤病灶一、两个方向的直径数值，但由于垂体瘤的形状多变，单一的直径并不能表示垂体瘤的体积。另外由于脑部扫描序列图像的数据量庞大，人工分割耗时费力。而且人工分割和影像信息分析的结果受医生的专业知识和操作熟练度影响，往往因人而异、主观性较大，很容易产生不可预测的结果。多模态垂体瘤磁共振(Magneticresonance，MR)图像处理与分析，综合运用图像处理、图像分析、模式识别、机器学习等现代计算机技术，可以大大改进临床垂体瘤医学影像处理与分析的现状，提供更高效、准确、客观的分析结果。垂体瘤属于脑肿瘤中的一种类型，临床实践中发现多模态脑垂体瘤的医学影像具有复杂的特性。其MR影像与其他脑肿瘤MR影像特点有相似之处，主要表现有：垂体瘤的结构复杂；其空间占位、形状和大小不确定；病变组织、周围软组织的灰度分布不均匀，密度分布有混叠，边界比较模糊；不同个体的组织结构存在较大差异；又由于MRI技术原理和图像采集过程中各种因素的影响，图像存在噪声、伪影等。因此要实现对垂体瘤的正确分割非常困难。脑肿瘤的分割任务包括正常脑组织(脑白质、脑灰质和脑脊液)、病变组织(水肿、肿瘤及其内部的坏死和囊变等)。从20世纪70年代开始，涌现出各种类型的脑肿瘤分割方法，主要可分为基于区域的分割方法、基于边缘的分割方法，以及结合区域和边界的分割方法，还有神经网络、模糊聚类、马尔科夫随机场模型、基于图论和基于图谱的其他分割方法。由于脑部医学图像的复杂性，各类分割方法仍只能就具体问题具体解决，尚没有一种通用的分割方法形成。目前的研究趋势呈现以下几个特点：(1)融合各自优点的多种分割算法结合。如结合模糊聚类和形变模型的图像分割方法、基于区域生长和神经网络的分割方法；还有近年来发展出的，基于图论的分割算法，如图割和主动形状模型联合分割算法、图割和形变模型联合分割算法等。(2)半自动的分割技术仍然是研究的重点方向。全自动分割算法可以大大提高脑肿瘤的分割效率，但在临床中还需允许医生根据具体病例进行特殊分割操作。交互式的半自动分割方法在自动分割结果的基础，使得医生可以评估分割结果、人工干预分割，弥补自动分割算法在进行特殊分割任务时的不足，达到精准分割、定量分析。如应用随机游走的脑肿瘤和肺肿瘤分割算法；也发展出一些交互式分割软件，如3DSlicer等等。近几年出现的自动、半自动交互式的脑肿瘤分割算法，可以大致分为基于机器学习的两个大类：生成算法和判别算法。生成算法从一系列先验图像(已知组织分类的图像)中学习健康脑组织类型的灰度空间分布，使用学习到的分布函数去“预测”健康组织、找出边界分割肿瘤。判别算法直接学习肿瘤和正常脑组织之间的灰度差异，用学习到的判别函数分类未标记的像素/体素。学习步骤既可以使用一些手工标记的病人训练图像进行离线执行，也可以通过用户在当前病人图像上标记种子区域在线完成。生成算法对于正常脑组织分割有优势，但对于任意形状的肿瘤分割则有难度。判别算法的优点是融合了正常脑组织和肿瘤两方面的灰度分布，有利于肿瘤的分割；但它仅使用局部图像特征，对多模态影像数据中不同组织的图像灰度要求相对一致，这是很难达到的，尤其对于MR图像来说。将两种方法结合，融合背景的脑结构解剖知识和病灶的灰度特性，可能会改善分割的正确性。近年来还出现了基于形变模型的脑肿瘤多模态影像分割方法。该类方法从图像中获取约束信息、目标位置、大小和形状等先验知识，既利用了图像本身的底层信息又结合了人体解剖结构的上层信息，适合于脑肿瘤多模态影像的三维数据分割。形变模型中基于水平集理论的几何形变模型，由于其可以解决曲线演化时拓扑结构变化问题，以及其他优于参数形变模型的特点，得到了广泛重视和发展。单个水平集函数只能将图像分割为目标区域和背景区域两个部分即图像两相分割，为满足图像中多个目标区域分割的需要，产生了多相水平集(MultiphaseLevelSet)方法。经过近二十年的发展，多相水平集思想逐渐应用到Mumford-Shah模型、Chan-Vese模型上。Li等人提出的LBF模型引入了高斯核函数，增强了水平集方法对边缘的捕获能力，在某些多相图像上仅使用一个水平集函数，解决了多相水平集算法需采用多个水平集函数所带来的计算复杂性问题；并由Peng等人发展为归一化的LBF模型。但是目前来看水平集方法的求解仍存在着方法复杂、计算量大的问题。以J.Egger为代表的研究团队提出了采用基于图论的图割和图搜索方法，设计出以模板为基础的半自动的垂体瘤分割方法。但从他们发表的研究论文来看，主要采用的病例均是形状规则、近似椭圆形状的垂体结构，分割难度相对不大。在图割算法和形变模型的能量函数中添加形状模型也是分割形状相对固定、变化不大的脑肿瘤的常见分割方法。但大部分的垂体瘤形状没有规则性，且可能会有对其他脑组织的浸润情况出现，所以，模板或形状模型在不规则形状的垂体瘤分割过程中所起作用有限。技术实现要素：本发明要解决的技术问题为：提出一种基于随机游走(Randomwalk，RW)和图割(Graphcut)的活动轮廓模型(Randomwalksandgraphcutsbasedactivecontourmodel，RGCACM)的MR图像三维交互分割方法，实现垂体瘤MR影像的三维分割，提高垂体瘤图像分割的准确性。本发明采取的技术方案具体为：一种基于随机游走和图割的活动轮廓模型的MR图像三维交互分割方法，包括以下步骤：S1，获取初始边界曲面，包括步骤：S11，选取随机游走算法所需的种子点，以该种子点作为三维中心点，从脑部MR三维数据中截取包含垂体瘤的局部三维MR图像数据；S12，利用三维随机游走算法对截取的局部三维MR图像数据进行处理，得到各体素点到达种子点的概率三维图；S13，利用最大期望值算法选取概率阈值，以分割出疑似目标区域，将疑似目标区域内的体素作为活动轮廓模型算法的前景点，疑似目标区域的边界曲面即初始边界曲面；S2，基于初始边界曲面，建立混合活动轮廓模型：对标准活动轮廓模型的能量函数中边界能量项采用标准测地线模型，区域能量项采用局部高斯模型描述的灰度最大后验概率，进而得到混合活动轮廓模型的能量函数；S3，模型离散化，包括：S31，给每个体素点定义一个二元变量，赋值1和0分别对应代表前景点和背景点，以将区域能量项函数离散化；S32，利用割测度将边界能量项函数即测地线模型离散化；S33，得到离散化后的能量函数；S4，构建图：对于截取的局部三维MR图像数据，将其中的每个体素点作为图的节点，采用每个体素点的6邻域进行图的构建；对初始边界曲面内和初始边界曲面外的体素点分别赋予初始值；并根据离散后的能量函数给节点间连接的边，节点与源点连接的边，和节点与汇点连接的边分别赋予对应的权值；S5，基于步骤S4构建的图进行图割计算：采用最大流和最小割算法得到分割结果，提取分割结果中的边界曲面即为分割轮廓；S6，以当前分割轮廓作为初始边界曲面，重复迭代步骤S2至步骤S5，直至分割结果收敛，输出最终分割结果。进一步的，本发明还包括步骤S7，利用三维中值滤波对步骤S5得到的演化曲面进行处理，得到滤波后分割结果；步骤S6以当前滤波后的分割轮廓作为初始边界曲面。优选的，步骤S11中，种子点的选取步骤为，在矢状位的二维显示图像中，选取中心位置切片，在切片中垂体瘤内近似中心位置选取一点，即作为种子点。所述近似中心位置即中心位置附近，此选取步骤为手动选取，故判断所选种子点是否近似中心位置依赖选取人来进行，以尽量靠近中心位置为基准，一般误差范围在垂体瘤直径的30％以内。优选的，步骤S11中，所述包含垂体瘤的局部三维MR图像数据的尺寸为71×71×41个体素。步骤S13中，根据概率图的灰度值分布来确定，取值范围为0～1。优选的，步骤S2中，给定图像域I(x):Ω→R为二维图像；闭合曲线C将图像分成两个独立区域，L(C)表示C的长度；设曲线C内区域为Ω1，曲线C外区域外为Ω2，对于图像域Ω的各像素p，考虑其半径为ρ的圆形邻域，定义为Ox＝{y:|x-y|≤ρ}，x为圆形邻域的中心像素点，y为圆形邻域中除x外的任意像素点；则区域Ox∩Ωi内的像素点灰度概率密度函数表示为：式中ui(x)和σi(x)分别为区域Ox∩Ωi内像素点x的局部灰度均值和方差，I(y)为像素点y的灰度值，边界能量项EEdge采用标准测地线模型，区域能量项EReigon采用局部高斯模型描述的灰度最大后验概率形式的能量函数E表示为：式中β是一个任意正常数，函数w是一个权重函数。进一步的，本发明步骤S2中，权重函数w采用截断高斯核形式，表示为：式中α为使得∫w(x,y)＝1成立的常数。即，本发明对邻域Ox中接近中心点x的像素y的灰度值I(y)赋予大的权重值。本发明步骤S3对模型进行离散化是为了与图割算法的能量函数相融合。优选的，步骤S31中，给每个体素点p定义一个二元变量xp，以对能量函数的区域能量项进行离散化：Object表示前景点即分割目标，Background表示背景点；区域能量项函数离散化后表示为：p、q指代图像中不同的体素，xp、xq指代不同体素的二元值(0或1)，N(p)指体素p的邻域。us(p)、ut(p)和σs(p)、σt(p)分别为ui(x)和σi(x)的离散形式，表示如下：步骤S32中，边界能量项函数离散化后表示为：ωpq为非负边权值，表达为三维形式即：ΔΦpq＝Δφpq·Δψpq＝π2/4，对应6邻域中单位圆的角度间隔，δ为单位格的尺寸，可取值δ＝1，|epq|为边界epq的单位长度，本发明中取值|epq|＝1，D是常数黎曼测度。优选的，本发明中ωpq＝π/4。基于离散后的区域项和边界项，即可得到离散后的混合活动模型能量函数。优选的，为避免图割算法的标准能量函数中区域和边界项的平衡问题，并且区域项会影响图割算法趋于分割出小区域。离散后的能量函数采用乘法形式，即，将区域能量项作为边界能量项的权重项，离散后的能量函数表示为：E＝EEdge(p,q)·E'Region(p)(10)如果邻域体素p,q同时位于分割目标内或外时，有|Es(p)+Et(q)|＞0和|Es(q)+Et(p)|＞0成立；如果p,q位于目标边界的两边，则有|Es(p)+Et(q)|≈0或|Es(q)+Et(p)|≈0成立。所以只有当邻域体素p,q分别位于目标边界的两边时，即p,q分别属于目标和背景区域时，公式(10)才能得到最小值。优选的，步骤S4中，定义体素点之间相互连接的边n-link权值为EEdge(p,q)·E'Region(p)；如果|Es(p)|＞|Et(p)|并且xp＝1，则体素点与汇点T(指代分割目标)相连接的边t-link的权值设置为无穷；如果|Es(p)|＜|Et(p)|并且xp＝0，则体素点与源点S(指代背景)相连接的边t-link的权值设置为无穷。最后即可用最大流/最小割算法求解分割结果，即对图的能量函数求解，得到图中各节点的标号(0或1)，此为现有技术。有益效果本发明方法是基于图割的活动轮廓模型(Graphcutsbasedactivecontourmodel，GCACM)算法的改进和拓展，利用活动轮廓模型解决图割算法的分割边缘阶梯化缺点，同时也利用图割的最大流/最小割算法解决了活动轮廓模型水平集方法的计算复杂性问题，以及容易陷于局部点的缺点。为解决图割算法中的边界种子点选择问题，采用随机游走算法大概确定目标区域，并利用垂体瘤的相对固定解剖位置(鼻子后方、双眼正中、视交叉下方等)，为随机游走算法手动交互提供种子点。结合位置先验知识和图割、活动轮廓模型实现垂体瘤的三维分割，可提高在例如垂体瘤的定位、与周围组织的病理关联以及体积测量等定量分析中的利用有效性。即本发明首次提供了一种具有可行性和有效性的半自动垂体瘤MR图像分割算法，结合了基于图论和基于活动轮廓模型的两类典型图像分割算法，可以实现规则和不规则垂体瘤的分割及体积估计，尤其对浸润型垂体瘤分割也有良好的准确性。交互式操作简单易行，只需用户鼠标选取两个点即可自动完成剩余的分割任务。算法中采用局部高斯分布描述图像区域灰度，可以有效解决MR图像垂体瘤病变区域灰度不均匀的问题。附图说明图1所示为本发明方法流程示意图；图2所示为本发明分割过程和结果示意图，其中a-用户选取一个前景点(黄色点)和一个背景点(绿色点)叠加在原始矢状位二维图像上，红色方框表示截取的感兴趣区；b-叠加在概率图上的初始轮廓；c-叠加在原始图像上的分割轮廓；d-分割结果的三维渲染模型；图3所示为步骤S2中像素的圆形邻域示意图，虚线圆形表示像素x的邻域。具体实施方式以下结合附图和具体实施例进一步描述。本发明基于随机游走和图割的活动轮廓模型的MR图像三维交互分割方法，包括以下步骤：S1，获取初始边界曲面，包括步骤：S11，选取随机游走算法所需的种子点，以该种子点作为三维中心点，从脑部MR三维数据中截取包含垂体瘤的局部三维MR图像数据；S12，利用三维随机游走算法对截取的局部三维MR图像数据进行处理，得到各体素点到达种子点的概率三维图；S13，利用最大期望值算法选取概率阈值，以分割出疑似目标区域，将疑似目标区域内的体素作为活动轮廓模型算法的前景点，疑似目标区域的边界曲面即初始边界曲面；S2，基于初始边界曲面，建立混合活动轮廓模型：对标准活动轮廓模型的能量函数中边界能量项采用标准测地线模型，区域能量项采用局部高斯模型描述的灰度最大后验概率，进而得到混合活动轮廓模型的能量函数；S3，模型离散化，包括：S31，给每个体素点定义一个二元变量，赋值1和0分别对应代表前景点和背景点，以将区域能量项函数离散化；S32，利用割测度将边界能量项函数即测地线模型离散化；S33，得到离散化后的能量函数；S4，构建图：对于截取的局部三维MR图像数据，将其中的每个体素点作为图的节点，采用每个体素点的6邻域进行图的构建；对初始边界曲面内和初始边界曲面外的体素点分别赋予初始值；并根据离散后的能量函数给节点间连接的边，节点与源点连接的边，和节点与汇点连接的边分别赋予对应的权值；S5，基于步骤S4构建的图进行图割计算：采用最大流和最小割算法得到分割结果，提取分割结果中的边界曲面即为分割轮廓；S6，以当前分割轮廓作为初始边界曲面，重复迭代步骤S2至步骤S5，直至分割结果收敛，输出最终分割结果。进一步的，本发明还可包括步骤S7，在S5与S6之间实施，即利用三维中值滤波对步骤S5得到的演化曲面进行处理，得到滤波后分割结果；步骤S6以当前滤波后的分割轮廓作为初始边界曲面。实施例1参考图1所示，本实施例方法的步骤简述为：图像分割种子点选取及初始化边界曲面获取、改进的GCACM算法的迭代分割、分割结果的三维中值滤波后处理。也即下述的初始化步骤、分割步骤和后处理步骤，具体描述如下。1、初始化步骤。主要为用户交互式手动选取Randomwalk算法所需的种子点，以及获取分割算法所需的初始边界曲面。为减少数据计算复杂度，并利用用户的医学背景知识，在原脑部MR三维数据中截取包含垂体瘤的数据立方体。具体做法是在矢状位的二维显示图像中，选取中心位置切片，再在垂体瘤内近似中心位置鼠标选取一点。以该点为中心截取数据立方体(本实验选用71×71×41大小的尺寸)，同时该点也作为Randomwalk算法的前景点，然后在该切片图像内背景区域任意选取一点充当Randomwalk算法的背景点。对截取的数据立方体进行三维Randomwalk算法处理，得到各体素到达前景种子点的概率三维图。利用EM算法选取概率阈值，初步分割出疑似目标的区域，区域内的体素将会被用作GCACM的前景点，区域的边界曲面即是初始边界曲面。获取初始边界曲面的过程如图2示意，以近似垂体瘤中心位置的一幅矢状位二维图像为例，首先a-用户选取一个前景点(黄色点)和一个背景点(绿色点)叠加在原始矢状位二维图像上，红色方框表示截取的感兴趣区；b图中为叠加在概率图上的初始轮廓，c图中为叠加在原始图像上的分割轮廓，d中为分割结果的三维渲染模型。2、分割步骤，包括混合活动轮廓模型的设计、模型离散化、图的构建及求解。(1)混合活动轮廓模型的建立参考图3所示，给定图像域I(x):Ω→R为二维图像。闭合曲线C将图像分成两个独立区域，L(C)表示C的长度。设曲线C内区域为Ω1，曲线C外区域外为Ω2。对于图像域Ω的每个像素p，考虑其半径为ρ的圆形邻域，定义为Ox＝{y:x-y≤ρ}。区域Ox∩Ωi内的像素点灰度概率密度函数可表示如下：式中ui(x)和σi(x)分别为区域Ox∩Ωi内像素点x的局部灰度均值和方差。能量函数结合了边界能量项和区域能量项，其中边界项EEdge采用标准测地线模型，区域项EReigon采用灰度局部高斯模型并且表示为最大后验概率(Maximumaposterioriprobability，MAP)形式。能量函数可写成如下形式：式中β是一个任意正常数。函数w是一个权重函数，对邻域Ox中接近中心点x的像素y灰度值I(y)赋予大的权重值。本算法中采用截断高斯核形式的w，表示如下：其中α是常数，使得∫w(x,y)＝1成立。(2)模型离散化为了与图割算法的能量函数融合，给图中每个体素点p定义一个二元变量xp，赋值1和0分别代表前景点(分割目标)和背景点(背景)，将区域项能量函数离散化，如式(4)所示。再利用Boykov等人提出的割测度(Cutmetric)，将测地线模型离散化表示，则能量函数的边界项表示为：式中ωpq是非负边权值，可表达为三维形式如下：上式中ΔΦpq＝Δφpq·Δψpq＝π2/4对应6邻域中单位圆的角度间隔，δ为单位格的尺寸，δ＝1，|epq|为边界epq的单位长度，|epq|＝1，D是常数黎曼测度，本发明中取ωpq＝π/4。能量函数的区域项离散形式可表达如下：式中p、q指代图像中不同的体素，xp、xq指代不同体素的二元值(0或1)，N(p)指体素p的邻域；us(p)、ut(p)和σs(p)、σt(p)分别为ui(x)和σi(x)的离散形式，如下所示：为避免图割算法的标准能量函数中区域和边界项的平衡问题，并且区域项会影响图割算法趋于分割出小区域。本算法使用乘法形式的能量函数，即使用区域项充当边界项的权重项，如公式(10)、(11)所示：E＝EEdge(p,q)·E'Region(p)(10)如果邻域像素p,q同时位于分割目标(即垂体瘤区域)内或外时，有|Es(p)+Et(q)|＞0和|Es(q)+Et(p)|＞0成立；如果p,q位于目标边界的两边，则有|Es(p)+Et(q)|≈0或|Es(q)+Et(p)|≈0成立。所以只有当邻域像素p,q位于目标边界的两边时，公式(10)才能得到最小值。公式(10)取得最小值代表图的能量为最小值，即能量函数取得最优结果，而能量达到最优时即分割完成，此时分属目标和背景区域的体素点被赋予了不同的标记值。(3)图的构建及求解本发明用MR数据体中的每个体素充当图的节点，采用6邻域构建三维图。图中体素之间相互连接的边n-link权值设为EEdge(p,q)·E'Region(p)；如果|Es(p)|＞|Et(p)|并且xp＝1，则体素与汇点T(对应目标)相连接的边t-link权值设置为无穷；如果|Es(p)|＜|Et(p)|并且xp＝0，则体素与源点S(对应背景)相连接的边t-link权值设置为无穷。最后用最大流/最小割算法求解分割结果，此处为现有技术，不做赘述。3、后处理对分割步骤所得的演化曲面再用三维中值滤波进行后处理，然后根据分割结果进行迭代分割直至结果收敛，得到最终的分割结果本发明的迭代分割流程简述如下：(1)采用Randomwalk算法获取初始边界曲面Sur；初始赋值Sur内的体素p对应的二元变量xp为1，Sur外的体素p对应的xp为0；(2)用公式(8)、(9)计算局部灰度的均值和方差us(p)、ut(p)和σs(p)、σt(p)；(3)用公式(10)表示的能量函数构建图；(4)用最大流/最小割算法求解分割结果，再根据分割结果更新体素p对应的二元变量xp值，实际上，求解得到分割结果的同时也即更新了体素p对应的二元变量值；(5)用三维中值滤波器对分割结果进行光滑后处理；(6)以当前分割轮廓取代初始轮廓，即S2中的初始边界曲面，重复步骤(2)-(5)直到分割结果收敛，求得的分割结果不再变化即为收敛。4、实验结果利用本发明在23位垂体瘤病人的MRT1W数据上进行分割实验，结果示意图如图2中c、d所示。表1给出23个病例分割结果Dice系数值(Dicesimilaritycoefficient，DSC)的最小值、最大值以及均值和标准方差值。DSC值定义为实验结果VE和金标准VG的相对重叠体积比，如公式(12)所示：此外为表现出本发明方法分割效果的优越性，将本算法和另外两种已发表的分割算法(加载于软件3DSlice中的基于区域生长的GrowCut算法和一般GCACM算法)进行分割结果DSC值的比较。比较结果参阅表1，可以看到本算法的分割结果DSC值均值为88.36％，优于GrowCut算法的81.61％和一般GCACM算法的70.88％。说明本算法在垂体瘤分割应用中具有较高的有效性和准确性。表1本算法和两种已发表分割算法的分割结果DSC值ResultsProposedmethodGrowCutGCACMmin74.86％67.62％49.45％max93.35％89.07％78.55％mean±SD88.36％±5.61％81.61％±7.37％70.88％±8.08％以上所述仅是本发明的优选实施方式，应当指出，对于本
技术领域：
的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。当前第1页1 2 3