一种基于水平集的前列腺磁共振图像分割方法与流程

本发明涉及一种基于水平集的前列腺磁共振图像分割方法，属于磁共振图像分割技术领域。

背景技术：

随着人口增长和生活习惯的改变，前列腺癌的发病率和死亡率近年来呈明显上升趋势。临床经验表明，前列腺癌若能尽早被发现，得到及时治疗，有很高的存活率，因此，对于前列腺癌诊断和治疗的相关研究具有重要意义。医学影像作为前列腺癌诊断和治疗的重要手段之一，发挥着越来越重要的作用。磁共振图像因其具有对软组织分辨率高，可多参数成像，能对任意断层进行扫描的特点，被认为是目前前列腺癌诊断和辅助治疗的最佳医学影像。

图像分割是基于图像早期诊断和治疗的基础，是首要解决的关键问题。目前，前列腺磁共振图像的分割研究大多集中在前列腺外轮廓分割，其分割方法主要有图论、形变模型、特定理论以及混合的图像分割方法等四大类。对基于磁共振图像的前列腺内外轮廓全分割研究近年来才开始展开。2011年，法国的Makni等人采用C-均值的方法最早实现了基于磁共振图像的前列腺内外轮廓全分割。2012年，荷兰的Litjens等人通过模式识别的方法利用解剖、灰度值和纹理特征对前列腺内部体素进行分类，实现了内外轮廓全分割。这两种方法都是先通过手动的方式实现外轮廓分割，并将其作为内部分割的初始化，使前列腺的全分割费时费力。2013年，美国的Toth等人利用多重耦合水平集活动轮廓模型的方法进行了内外轮廓全分割，该方法的分割效果更多依赖于所分割图像的质量，而且分割一张图像需要训练大量图集，消耗时间较长。2014年，加拿大的Qiu等人利用优化连续最大流模型的方法进行了前列腺的内外全分割，改善了分割效果并提高了分割效率。由于前列腺内部区域磁共振图像存在噪声，灰度显示不均匀，区域边界模糊不清，使基于磁共振图像的前列腺内外轮廓全分割研究一直是个挑战。

技术实现要素：

针对上述问题，本发明要解决的技术问题是提供一种基于水平集的前列腺磁共振图像分割方法，在构建统一水平集能量函数的基础上，首先基于纵向弛豫时间图像，将前列腺从其周围组织中分割出来，即实现前列腺的外轮廓分割，其次基于横向弛豫时间图像，将外轮廓作为约束条件来实现前列腺的内部区域分割，进而实现前列腺的内外轮廓全分割。

上述目的主要通过以下方案实现：

本发明的一种基于水平集的前列腺磁共振图像分割方法，其特征在于：所述方法的具体实现过程为：

步骤一：定义水平集演化方程

在Ω域内定义一个水平集函数能量函数ε(φ)定义为：

ε(φ)＝μR_p(φ)+αε_drive(φ) (1)

其中，R_p(φ)是水平集的距离调整项，ε_drive(φ)是轮廓驱动能量项，μ＞0，α＜0，都为常数；

水平集的距离调整项R_p(φ)定义为：

其中，p是能量密度函数，

能量密度函数构造为：

能量密度函数p(s)具有两个极值点，分别是s＝0和s＝1，其一阶导数和二阶导数为：

式(2)中函数R_p(φ)的加托导数为：

其中，函数d_p定义为：

轮廓驱动能量项ε_drive(φ)定义为：

其中，g是边界约束函数，H是单位阶跃函数，通常将单位阶跃函数H近似地用函数H_ε来代替，且定义为：

H_ε的导数δ_ε为：

轮廓驱动能量函数ε_drive(φ)的加托导数为：

求解梯度流方程的稳态解，

其中，是函数ε(φ)的加托导数；

将式(6)和式(11)代入(12)中，可以得到能量函数ε(φ)的梯度流表达式为：

式(13)所示的偏微分方程就是基于式(1)的前列腺内外轮廓分割的水平集演化方程；

瞬态偏导数可以近似采用正向有限差分方程进行求解，时变函数φ(x,y,t)的离散形式用来表示，则水平集演化方程可以离散为如下所示的有限差分方程：

步骤二：外轮廓分割

读取原始的纵向弛豫时间图像，选择外分割初始化方法-变形椭圆法：

基本椭圆参数方程如式(15)所示：

其中，a_x是x方向的半轴长，a_y是y方向的半轴长；

沿着y轴通过转换基本的椭圆方程获得变形椭圆的参数方程ψ(x_d,y_d)，如式(16)所示：

其中，

将已确定的前列腺变形椭圆所确定的区域设为S_e，则初始水平集函数为：

其中，c₀为正常数；

在式(16)和式(17)中，(x_c,y_c)是变形椭圆的中心坐标，t_y∈[-1,1]是描述椭圆上部沿着y轴方向线性变尖的参数，b_y∈[-1,0)∪(0,1]是描述椭圆下部沿着y轴方向内凹弯曲的参数，调整式(16)和式(17)相应的参数，使得可变形椭圆最大限度逼近前列腺的外轮廓形状；

然后，确定外轮廓边界约束函数：

在纵向弛豫时间图像中，假定I为前列腺图像，定义图像I的边界指示器为：

其中，G_σ是方差为σ的高斯核，将式(19)作为前列腺外轮廓分割的边界约束函数，并给定参数值。

最后对水平集演化方程(14)进行迭代求解，实现前列腺的外轮廓分割；

步骤三：内部区域分割

读取原始横向弛豫时间图像，选择内分割初始化方法-多线段拟合法：

在中央腺内依次选取N个点，使这N个点首尾相连形成一封闭区域，设为S_N，则初始水平集函数为：

其中，c₀为正常数；

然后，确定内轮廓边界约束函数：

采用全向边界梯度作为边界指示器来描述前列腺中央腺图像的边界特征，假定I为前列腺图像，I_i,j为I的某一元素，设定为中心元素，其相邻的8元素分别为I_i-1,j-1，I_i-1,j，I_i-1,j+1，I_i,j-1，I_i,j+1，I_i+1,j-1，I_i+1,j，I_i+1j+1，为求取这8元素与中心元素的差值，定义如下对应的8个卷积模板，

中心元素与相邻8元素的差值计算为：

Dif_lu＝conv2(I,Temp_lu,'same') (29)

Dif_u＝conv2(I,Temp_u,'same') (30)

Dif_ru＝conv2(I,Temp_ru,'same') (31)

Dif_l＝conv2(I,Temp_l,'same') (32)

Dif_r＝conv2(I,Temp_r,'same') 33)

Dif_ld＝conv2(I,Temp_ld,'same') (34)

Dif_d＝conv2(I,Temp_d,'same') (35)

Dif_rd＝conv2(I,Temp_rd,'same') (36)

conv2是卷积运算符，图像I的全向边界梯度函数定义为：

Grad_I＝[Grad_I_x Grad_I_y Grad_I_xy- Grad_I_xy+] (37)

其中，各项分别定义为：

图像I的全向边界梯度模定义为：

|Grad_I|＝sqrt(Grad_I_x²+Grad_I_y²+Grad_I_xy-²+Grad_I_xy+²) (42)

式(12)中的边界约束函数为：

式(43)称为前列腺内轮廓分割的边界约束函数，并给定参数值；

最后对水平集演化方程(14)进行迭代求解，获得前列腺内部中央腺的轮廓；将第二步得到的外轮廓与第三步所得到的中央腺轮廓进行区域相减，便得到前列腺外周带区域，进而实现了前列腺的全面分割。

本发明的有益效果为：

1、提出了基于纵向弛豫时间图像与横向弛豫时间图像相结合的前列腺两步分割法，综合了纵向弛豫时间图像内外信号对比鲜明，可实现外轮廓分割和横向弛豫时间图像内部结构显示清晰，外周带与中央腺信号形成明显对比，可实现内区域的分割的优点，克服了纵向弛豫时间图像难以区分内部结构和横向弛豫时间图像内部清晰的多区域信号灰度值将干扰外轮廓的提取与分割的缺点。

2、所建立的能量函数中融入了距离调整项，在演化过程中可以不断的进行调整，引起的周围扩散效应可以维持期望的形状与期望轮廓附近的距离，在不必重新初始化的同时避免了通用水平集方法由于不断的初始化引起的数值错误。

3、采用了初始轮廓接近内外轮廓的外分割初始化-变形椭圆法和内分割初始化-多线段拟合法来分别初始化水平集函数，可以提高分割效果。

附图说明

为了易于说明，本发明由下述的具体实施及附图作以详细描述。

图1为本发明方法的流程图；

图2为本发明外分割初始化-变形椭圆法中不同a_x参数的轮廓示意图；

图3为本发明外分割初始化-变形椭圆法中不同b_y参数的轮廓示意图；

图4为本发明的人机交互界面。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚明了，下面通过附图中示出的具体实施例来描述本发明。但是应该理解，这些描述只是示例性的，而并非要限制本发明的范围。此外，在以下说明中，省略了对公知结构和技术的描述，以避免不必要地混淆本发明的概念。

如图1、图2、图3、图4所示，本具体实施方式采用以下技术方案：一种基于水平集的前列腺磁共振图像分割方法，其特征在于：所述方法的具体实现过程为：

步骤一：定义水平集演化方程

在Ω域内定义一个水平集函数能量函数ε(φ)定义为：

ε(φ)＝μR_p(φ)+αε_drive(φ) (1)

其中，R_p(φ)是水平集的距离调整项，ε_drive(φ)是轮廓驱动能量项，驱动水平集函数曲线向前列腺轮廓边界运动，μ＞0，α＜0，都为常数；

水平集的距离调整项R_p(φ)定义为：

其中，p是能量密度函数，

为了避免边界效应，能量密度函数构造为：

能量密度函数p(s)具有两个极值点，分别是s＝0和s＝1，其一阶导数和二阶导数为：

式(2)中函数R_p(φ)的加托导数为：

其中函数d_p定义为：

轮廓驱动能量项ε_drive(φ)定义为：

其中，g是边界约束函数，H是单位阶跃函数，通常将单位阶跃函数H近似地用函数H_ε来代替，且定义为：

H_ε的导数δ_ε为：

轮廓驱动能量函数ε_drive(φ)的加托导数为：

为了求取能量函数ε(φ)的最小值，常规方法就是求解梯度流方程的稳态解：

其中，是函数ε(φ)的加托导数；

将式(6)和式(11)代入(12)中，可以得到能量函数ε(φ)的梯度流表达式为：

式(13)所示的偏微分方程就是基于式(1)的前列腺内外轮廓分割的水平集演化方程；

瞬态偏导数可以近似采用正向有限差分方程进行求解，时变函数φ(x,y,t)的离散形式用来表示，这样水平集演化方程可以离散为如下所示的有限差分方程：

步骤二：外轮廓分割

读取原始的纵向弛豫时间图像，选择外分割初始化方法-变形椭圆法：

基本椭圆参数方程如式(15)所示：

其中，a_x是x方向的半轴长，a_y是y方向的半轴长；

由于前列腺横断轴位各层的外轮廓形状是沿着y轴改变的，因此沿着y轴通过转换基本的椭圆方程，可以获得变形椭圆的参数方程ψ(x_d,y_d)，如式(16)所示：

其中，

将已确定的前列腺变形椭圆所确定的区域设为S_e，则初始水平集函数为：

其中，c₀为正常数。

在式(16)和式(17)中，(x_c,y_c)是变形椭圆的中心坐标，t_y∈[-1,1]是描述椭圆上部沿着y轴方向线性变尖的参数，b_y∈[-1,0)∪(0,1]是描述椭圆下部沿着y轴方向内凹弯曲的参数。通过调整式(19)和式(20)相应的参数，使得可变形椭圆最大限度逼近前列腺的外轮廓形状；

然后，确定外轮廓边界约束函数：

在纵向弛豫时间图像中，假定I为前列腺图像，定义图像I的边界指示器为：

其中，G_σ是方差为σ的高斯核，式中卷积用来平滑前列腺图像，减小噪声的影响，将式(19)作为前列腺外轮廓分割的边界约束函数，并给定参数值。

最后对水平集演化方程(14)进行迭代求解，实现前列腺的外轮廓分割；

步骤三：内部区域分割

读取原始横向弛豫时间图像，选择内分割初始化方法-多线段拟合法：

在中央腺内依次选取N个点，使这N个点首尾相连形成一封闭区域，设为S_N，则初始水平集函数为：

其中，c₀为正常数；

然后，确定内轮廓边界约束函数：

采用全向边界梯度作为边界指示器来描述前列腺中央腺图像的边界特征，

假定I为前列腺图像，I_i,j为I的某一元素，设定为中心元素。其相邻的8元素分别为I_i-1,j-1，I_i-1,j，I_i-1,j+1，I_i,j-1，I_i,j+1，I_i+1,j-1，I_i+1,j，I_i+1,j+1。为求取这8元素与中心元素的差值，定义如下对应的8个卷积模板，

中心元素与相邻8元素的差值计算为：

Dif_lu＝conv2(I,Temp_lu,'same') (29)

Dif_u＝conv2(I,Temp_u,'same') (30)

Dif_ru＝conv2(I,Temp_ru,'same') (31)

Dif_l＝conv2(I,Temp_l,'same') (32)

Dif_r＝conv2(I,Temp_r,'same') 33)

Dif_ld＝conv2(I,Temp_ld,'same') (34)

Dif_d＝conv2(I,Temp_d,'same') (35)

Dif_rd＝conv2(I,Temp_rd,'same') (36)

conv2是卷积运算符，图像I的全向边界梯度函数定义为：

Grad_I＝[Grad_I_x Grad_I_y Grad_I_xy- Grad_I_xy+] (37)

其中，各项分别定义为：

图像I的全向边界梯度模定义为：

|Grad_I|＝sqrt(Grad_I_x²+Grad_I_y²+Grad_I_xy-²+Grad_I_xy+²) (42)

式(12)中的边界约束函数为：

式(43)称为前列腺内轮廓分割的边界约束函数，并给定参数值；

最后对水平集演化方程(14)进行迭代求解，获得前列腺内部中央腺的轮廓；将第二步得到的外轮廓与第三步所得到的中央腺轮廓进行区域相减，便得到前列腺外周带区域，进而实现了前列腺的全面分割。

以上显示和描述了本发明的基本原理和主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。