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

技术特征：

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

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

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

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

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

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

$R_p\left(\mathrm{\φ}\right)\stackrel{\mathrm{\Δ}}{=}{\∫}_{\mathrm{\Ω}}p\left(\right|\▿\mathrm{\φ}\left|\right)dx---\left(2\right)$

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

能量密度函数构造为：

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

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

$\frac{\∂R_p}{\∂\mathrm{\φ}}=-div\left(d_p\right(|\▿\mathrm{\φ}|)\▿\mathrm{\φ})---\left(6\right)$

其中函数d_p定义为：

$d_p\stackrel{\mathrm{\Δ}}{=}\frac{p^{\′}\left(s\right)}{s}---\left(7\right)$

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

${\mathrm{\ϵ}}_{drive}\left(\mathrm{\φ}\right)\stackrel{\mathrm{\Δ}}{=}{\∫}_{\mathrm{\Ω}}gH(-\mathrm{\φ})dx---\left(8\right)$

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

$H_{\mathrm{\ϵ}}=\left\{\begin{array}{cc}\frac{1}{2}\[1+\frac{x}{\mathrm{\ϵ}}+\frac{1}{\mathrm{\π}}cos\left(\frac{\mathrm{\π}x}{\mathrm{\ϵ}}\right)\],& \left|x\right|\≤\mathrm{\ϵ}\\ 0,& x\>\mathrm{\ϵ}\\ 0,& x\<-\mathrm{\ϵ}\end{array}\right.---\left(9\right)$

H_ε的导数δ_ε为：

${\mathrm{\δ}}_{\mathrm{\ϵ}}=\left\{\begin{array}{cc}\frac{1}{2\mathrm{\ϵ}}\[1-sin\left(\frac{\mathrm{\π}x}{\mathrm{\ϵ}}\right)\],& \left|x\right|\≤\mathrm{\ϵ}\\ 0,& \left|x\right|\>\mathrm{\ϵ}\end{array}\right.---\left(10\right)$

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

$\frac{\∂{\mathrm{\ϵ}}_{drive}}{\∂\mathrm{\φ}}={\mathrm{g\δ}}_{\mathrm{\ϵ}}\left(\mathrm{\φ}\right)---\left(11\right)$

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

$\frac{\∂\mathrm{\φ}}{\∂t}=-\frac{\∂\mathrm{\ϵ}}{\∂\mathrm{\φ}}---\left(12\right)$

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

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

$\frac{\∂\mathrm{\φ}}{\∂t}=\mathrm{\μ}div\left(d_p\right(|\▿\mathrm{\φ}|)\▿\mathrm{\φ})+{\mathrm{\αg\δ}}_{\mathrm{\ϵ}}\left(\mathrm{\φ}\right)---\left(13\right)$

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

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

${\mathrm{\φ}}_{i,j}^{k+1}={\mathrm{\φ}}_{i,j}^k+\mathrm{\Δ}tL\left({\mathrm{\φ}}_{i,j}^k\right),k=0,1,2...---\left(14\right)$

步骤二：外轮廓分割

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

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

$\left\{\begin{array}{c}x=a_{x}cos\mathrm{\α}\\ y=a_{y}sin\mathrm{\α}\end{array}\right.,\mathrm{\α}\∈\[-\mathrm{\π},\mathrm{\π}\]---\left(15\right)$

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

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

$\left\{\begin{array}{c}{x}_d=(a_y/b_y-a_y\sin \mathrm{\α})sin\mathrm{\β}+x_c\\ y_d=a_y/b_y-\left(a_y/b_y-a_y\sin \mathrm{\α}\right)cos\mathrm{\β}+y_c\end{array}\right.,\mathrm{\α}\∈\[-\mathrm{\π},\mathrm{\π}\]---\left(16\right)$

其中，

$\mathrm{\β}=\frac{\[t_{y}sin\mathrm{\α}+1\]a_{x}cos\mathrm{\α}}{a_y/b_y-a_y\sin \mathrm{\α}}---\left(17\right)$

将已确定的前列腺变形椭圆所确定的区域设为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\stackrel{\mathrm{\Δ}}{=}\frac{1}{1+|\▿G_{\mathrm{\σ}}*I{|}^2}---\left(19\right)$

其中，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个卷积模板，

$Temp\_lu=\left[\begin{array}{ccc}0& 0& 0\\ 0& 1& 0\\ 0& 0& -1\end{array}\right]---\left(21\right)$

$Temp\_u=\left[\begin{array}{ccc}0& 0& 0\\ 0& 1& 0\\ 0& -1& 0\end{array}\right]---\left(22\right)$

$Temp\_ru=\left[\begin{array}{ccc}0& 0& 0\\ 0& 1& 0\\ -1& 0& 0\end{array}\right]---\left(23\right)$

$Temp\_l=\left[\begin{array}{ccc}0& 0& 0\\ 0& 1& -1\\ 0& 0& 0\end{array}\right]---\left(24\right)$

$Temp\_r=\left[\begin{array}{ccc}0& 0& 0\\ -1& 1& 0\\ 0& 0& 0\end{array}\right]---\left(25\right)$

$Temp\_ld=\left[\begin{array}{ccc}0& 0& -1\\ 0& 1& 0\\ 0& 0& 0\end{array}\right]---\left(26\right)$

$Temp\_d=\left[\begin{array}{ccc}0& -1& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right]---\left(27\right)$

$Temp\_rd=\left[\begin{array}{ccc}-1& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right]---\left(28\right)$

中心元素与相邻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)

其中，各项分别定义为：

$Grad\_I_x=\frac{abs(Dif\_l)+abs(Dif\_r)}{2}---\left(38\right)$

$Grad\_I_y=\frac{abs(Dif\_u)+abs(Dif\_d)}{2}---\left(39\right)$

$Grad\_I_{xy-}=\frac{abs(Dif\_lu)+abs(Dif\_rd)}{2}---\left(40\right)$

$Grad\_I_{xy+}=\frac{abs(Dif\_ru)+abs(Dif\_ld)}{2}---\left(41\right)$

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

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

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

$g=\frac{1}{1+|Grad\_I{|}^2}---\left(43\right)$

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

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