CN102509292A - 一种心脏核磁共振图像的快速分割方法 - Google Patents

Abstract

本发明涉及一种心脏核磁共振图像的快速分割方法，包括以下步骤：一、对取得的心脏核磁共振图像进行高斯滤波预处理；二、在预处理后的图像上计算基于扩展邻域和噪声平滑的广义梯度矢量流的外力场；三、定义心脏左心室内膜初始化轮廓位置；四、对心脏左心室内膜进行分割；五、将心脏左心室内膜的最终分割轮廓结果定义为心脏左心室外膜的初始化轮廓位置；六、将原始边缘图中内膜轮廓所包围区域的边缘强度置为0，重新计算基于扩展邻域和噪声平滑的广义梯度矢量流的外力场；七、对心脏左心室外膜进行分割。本发明基于卷积运算，考虑了椭圆形状能量约束，具有运算速度快、捕捉范围大、抗噪能力强的优点，且在弱边界保护和深度凹陷区域的分割上性能卓越，能准确地分割心脏左心室内、外膜。

Description

一种心脏核磁共振图像的快速分割方法

技术领域

本发明涉及一种图像分割方法，特别涉及一种心脏核磁共振图像分割方法，属于医学图像分析领域。

背景技术

心脏MRI(magnetic resonance imaging)能够提供高分辨率、高品质的图像，对心脏的解剖结构和功能进行准确的描述，是当前医学图像分析领域的研究热点之一，也是心脏疾病诊断的重要辅助手段，对心血管疾病的早期无创诊断和准确预后评估具有重要意义.为了充分利用图像中的解剖信息，为临床诊断提供量化、直观的参考，首先必须分割出左室壁的内、外膜.然而，由于心脏的运动和血液的高速流动，图像受噪声干扰，使得心脏MR图像的分割仍是一个值得深入研究的问题.

近年来，对心脏MR(magnetic resonance)图像的分割国内外都有广泛研究，这些方法大致可以分为基于形态学的方法，基于模糊聚类的方法，基于模板的方法以及基于主动轮廓模型的方法等.主动轮廓模型能够将有关目标形状的先验知识和来自图像的知识融入一个统一的过程中，是当前图像分割领域的热点方法，也是心脏MR图像分割中的主流方法，在国内外都有广泛研究。Hong等人采用基于Lagrange动力学的B样条Snake模型来提取左室壁内膜，Makowski等人采用气球Snake模型来分割左室壁内膜，设计了专门的方法来解决轮廓缠绕问题。Nguyen等人对传统Snake，GVF(gradient vector flow)Snake和气球Snake模型分割左室壁内膜的结果作了比较，并与手工勾勒的轮廓进行对比验证，其中，GVFSnake模型性能最好。Jolly等人首先采用极大鉴别分析方法来找到左室壁内膜的大致轮廓，并用Snake模型提取左室壁内膜。Nachtomy等人提出了一种基于阈值的方法提取左心室内、外膜，但由于阈值的局限性，结果并不令人满意。Pednekar等人针对图像的模糊特点，提出了一种基于模糊分析的左室壁内、外膜分割方法。国内对心脏图像的分割也有相关研究，周寿军等人用梯度矢量流(gradient vector flow，GVF)Snake模型分割左心室时，引入广义模糊集合理论，提出了广义模糊梯度矢量流。秦安等人将广义模糊梯度矢量流与几何主动轮廓相结合来分割左心室内膜，然后采用一种区域灰度均值和距离约束的外力来分割左心室外膜。周则明等人将简化Snake模型用于心脏图像分割，并用贪婪算法求解能量泛函的局部极小点。王元全和贾云得提出了二种基于Snake模型的分割策略，引入了形状约束，提出了退化最小曲面梯度矢量流(dmsGVF)和卷积虚拟静电场(CONVEF)外力模型。

心脏左心室分割的困难主要来自如下三个方面：首先图像灰度不均.这种灰度不均可能是成像过程中射频脉冲的干扰或者磁场强度不均，也可能是血液高速运动冲撞心肌壁造成的；其次，乳突肌的干扰.一般来说，乳突肌与心肌相连的部分被认为是心肌的一部分，而漂浮在血池中的部分则被认为不是心肌的一部分；另外，由于左室壁与右室壁及周围其他组织(如肝脏)等灰度非常接近，形成弱边界，这时基于主动轮廓模型的方法分割左心室外膜时往往发生泄露。现有的方法对于这些问题没有提出很好的解决方案。

发明内容

本发明的目的是针对心脏左心室内、外膜分割存在的图像灰度不均、乳突肌的干扰和弱边界的难点，提出一种快速、高效、鲁棒、准确的图像分割方法，分割左室壁内、外膜。

本发明的方法是基于扩展邻域和噪声平滑的广义梯度矢量流(ENGGVF)模型提出的，其原理如下：

Snake模型是一种自顶向下的图像分析方法，具有传统方法无法比拟的优点。Snake模型用曲线c(s)＝(x(s)，y(s))(s∈[0，1])来定义，这是以归一化弧长s作为参数的曲线表达形式。它通过极小化如下的能量泛函来确定目标轮廓：

$E_{\mathrm{snake}}=\∫\frac{1}{2}(\mathrm{\α}{\left|c_s\right|}^2+\mathrm{\β}{\left|c_{\mathrm{ss}}\right|}^2+E_{\mathrm{ext}}\left(c\left(s\right)\right))\mathrm{ds}---\left(1\right)$

其中，α和β为弹性和刚性系数，控制着弹性和刚性能量的大小；一阶导数项c_s刻画了曲线的连续性，是曲线的弹性能量；二阶导数项c_ss刻画了曲线的光滑性，是曲线的刚性能量；这两个导数项构成Snake模型的内部能量；E_ext(c(s))是Snake模型的外部能量；根据变分法原理，能量泛函式(1)的最小化可以通过求解如下Euler方程得到：

$c\left(s\right)={\mathrm{\αc}}_{\mathrm{ss}}\left(s\right)-{\mathrm{\βc}}_{\mathrm{ssss}}\left(s\right)-{\▿E}_{\mathrm{ext}}\left(c\left(s\right)\right)---\left(2\right)$

其中，c_ss(s)为曲线c(s)关于s的二阶导数，c_ssss(s)为曲线c(s)关于s的四阶导数。当方程式(2)的解收敛时，就得到了待分割目标的轮廓。这时，可以将Snake轮廓的运动过程看成其内、外力的平衡过程，αc_ss(s)-βc_ssss(s)称为Snake模型的内力，

称为其外力。外力在Snake模型的演化中起决定性作用，对外力的研究是Snake模型研究的一个重要方面。由于式(2)定义的Snake模型外力

是基于图像梯度的，因此其捕捉范围小，不能进入深度凹陷区域，初始化敏感。针对这些问题，Xu和Jerry提出了用梯度矢量流场(GVF)作为新的外力条件代替式(2)中的

来约束动态轮廓线，并将其定义为V(x，y)＝[u(x，y)，v(x，y)]，其满足下列能量泛函的最小值：

$\mathrm{\ϵ}=\∫\∫\mathrm{\μ}{|\▿V|}^2+{|\▿f|}^2{|V-\▿f|}^2\mathrm{dxdy}---\left(3\right)$

μ为权重系数，f是边缘图，可以由其他边缘检测算子得到或者用图像梯度来近似。使用变分原理，GVF场可以通过解下列欧拉方程获得：

$\left\{\begin{array}{c}\mathrm{\μ\Δu}-{|\▿f|}^2(u-f_x)=0\\ \mathrm{\μ\Δv}-{|\▿f|}^2(v-f_y)=0\end{array}\right.---\left(4\right)$

其中Δ为拉普拉斯算子。作为主动轮廓模型中最成功的外部力场之一，GVF模型扩展了边缘映射的梯度向量，同时通过一个各向同性的扩散过程来抑制噪声。将μ和

用

k为调整常量和

代替，得到GGVF模型，它在深度凹陷区域有更好的收敛效果，GGVF场可以通过解下列欧拉方程获得：

$\left\{\begin{array}{c}g\left(\right|\▿f\left|\right)\mathrm{\Δu}-h\left(\right|\▿f\left|\right)(u-f_x)=0\\ g\left(\right|\▿f\left|\right)\mathrm{\Δv}-h\left(\right|\▿f\left|\right)(v-f_y)=0\end{array}\right.---\left(5\right)$

在4邻域中，拉普拉斯算子的计算可以用下列等式来近似：

$\mathrm{\Δf}=\frac{{\∂}^{2}f}{{\∂x}^2}+\frac{{\∂}^{2}f}{{\∂y}^2}=f(i,j+1)-2f(i,j)+f(i,j-1)+f(i+1,j)-2f(i,j)+f(i-1,j)---\left(6\right)$

$=f(i,j+1)+f(i,j-1)+f(i+1,j)+f(i-1,j)-4f(i,j)$

在图像中，对式(6)的计算的可以借助于掩模和卷积来实现，式(6)对应的掩模运算为：

$\mathrm{\Δf}=f\⊗G_4---\left(7\right)$

其中，

表示卷积，G₄表示4邻域的拉普拉斯算子掩模。

$G_4=\left[\begin{array}{ccc}0& 1& 0\\ 1& -4& 1\\ 0& 1& 0\end{array}\right]---\left(8\right)$

通过将拉普拉斯算子的计算扩展到了更大的邻域，更多的图像信息就会得到利用，基于这一点，我们提出了基于扩展邻域的广义梯度矢量流。这里，在梯度矢量流的拉普拉斯算子的计算中，我们采用24邻域的掩模来取代原来的4邻域的掩模运算。该外力场可以通过解下列欧拉方程获得：

$\left\{\begin{array}{c}{u}_t=g\left(\right|\▿f\left|\right)(u\⊗G_{24})-h\left(\right|\▿\left|\right)(u-f_x)\\ v_t=g\left(\right|\▿f\left|\right)(v\⊗G_{24})-h\left(\right|\▿f\left|\right)(v-f_y)\end{array}\right.---\left(9\right)$

其中，

表示卷积，G₂₄表示24邻域的拉普拉斯算子掩模。

$G_{24}=\left[\begin{array}{ccccc}1& 1& 1& 1& 1\\ 1& 1& 1& 1& 1\\ 1& 1& -24& 1& 1\\ 1& 1& 1& 1& 1\\ 1& 1& 1& 1& 1\end{array}\right]---\left(10\right)$

我们可以将拉普拉斯算子的掩膜分成两部分，一部分是中值滤波的掩膜(RM)，一部分是全通滤波的掩膜(AP)。拉普拉斯算子对噪声敏感的原因在于原始的全通滤波掩膜是全通的，对噪声没有抑制作用。我们把原始的全通滤波掩膜用具有边缘保持和噪声平滑的噪声滤波掩膜(NS)代替，能够获得更好的效果。24邻域的拉普拉斯算子掩模(G₂₄)可分解为

G₂₄＝NS₂₄-RM₂₄                   (11)

其中，NS₂₄和RM₂₄分别为噪声滤波掩膜和中值滤波的掩膜。

${\mathrm{NS}}_{24}=\left[\begin{array}{ccccc}0& 0& 1/12& 0& 0\\ 0& 0& 1/12& 0& 0\\ 1/12& 1/12& 1/3& 1/12& 1/12\\ 0& 0& 1/12& 0& 0\\ 0& 0& 1/12& 0& 0\end{array}\right],$ ${\mathrm{RM}}_{24}=\left[\begin{array}{ccccc}1/24& 1/24& 1/24& 1/24& 1/24\\ 1/24& 1/24& 1/24& 1/24& 1/24\\ 1/24& 1/24& 0& 1/24& 1/24\\ 1/24& 1/24& 1/24& 1/24& 1/24\\ 1/24& 1/24& 1/24& 1/24& 1/24\end{array}\right]$

基于这样的考虑，本发明提出了一种新的外部力场：基于扩展邻域和噪声平滑的广义梯度矢量流(ENGGVF)模型，并采用这一外力模型来分割左室壁。相对于原始的梯度矢量流来说，基于扩展邻域和噪声平滑的广义梯度矢量流模型除了采用了权重因子和扩展邻域卷积运算外，还将拉普拉斯算子模板中加入噪声平滑模板，具有运算速度快、捕捉范围大、抗噪能力强，在弱边界保护和深度凹陷区域的分割上性能卓越。

基于扩展邻域和噪声平滑的广义梯度矢量流(ENGGVF)模型可通过解下列欧拉方程获得：

$\left\{\begin{array}{c}{u}_t=g\left(\right|\▿f\left|\right)(u\⊗{\mathrm{NS}}_{24}-u\⊗{\mathrm{RM}}_{24})-h\left(\right|\▿f\left|\right)(u-f_x)\\ v_t=g\left(\right|\▿f\left|\right)(v\⊗{\mathrm{NS}}_{24}-v\⊗{\mathrm{RM}}_{24})-h\left(\right|\▿f\left|\right)(v-f_y)\end{array}\right.---\left(12\right)$

其中，

表示卷积， $g\left(\right|\▿f\left|\right)=\exp (-|\▿f|/k),$ k为调整常量， $h\left(\right|\▿f\left|\right)=1-g\left(\right|\▿f\left|\right),$ ${\mathrm{NS}}_{24}=\left[\begin{array}{ccccc}0& 0& 1/12& 0& 0\\ 0& 0& 1/12& 0& 0\\ 1/12& 1/12& 1/3& 1/12& 1/12\\ 0& 0& 1/12& 0& 0\\ 0& 0& 1/12& 0& 0\end{array}\right],$ ${\mathrm{RM}}_{24}=\left[\begin{array}{ccccc}1/24& 1/24& 1/24& 1/24& 1/24\\ 1/24& 1/24& 1/24& 1/24& 1/24\\ 1/24& 1/24& 0& 1/24& 1/24\\ 1/24& 1/24& 1/24& 1/24& 1/24\\ 1/24& 1/24& 1/24& 1/24& 1/24\end{array}\right]$ 为卷积模板，该外力场基于卷积运算，计算速度快。

基于以上思想，本发明提供了一种心脏核磁共振图像分割方法，包括以下步骤：

一、对取得的心脏核磁共振图像进行高斯滤波预处理

根据方程

对图像进行高斯滤波预处理，式中I₀为输入的原始图像结构信息，G_σ为标准差为σ的二维高斯函数，

表示卷积运算；

二、在预处理后的图像上计算基于扩展邻域和噪声平滑的广义梯度矢量流的外力场，记为ENGGVF外力场，具体方法为：

1)定义ENGGVF外力场F_out的初始值

定义图像I(x，y)的边缘映射为f(x，y)，设f_x和f_y分别为边缘映射f沿x轴方向和y轴方向的一阶导数，边缘图的梯度向量构成了一个向量场V(x，y)＝[u(x，y)，v(x，y)]＝[f_x，f_y]，作为外力场F_out的初始值；

2)根据外力场F_out的初始值，计算ENGGVF外力场

ENGGVF外力场的迭代公式如下：

$\left\{\begin{array}{c}{u}_t=g\left(\right|\▿f\left|\right)(u\⊗{\mathrm{NS}}_{24}-u\⊗{\mathrm{RM}}_{24})-h\left(\right|\▿f\left|\right)(u-f_x)\\ v_t=g\left(\right|\▿f\left|\right)(v\⊗{\mathrm{NS}}_{24}-v\⊗{\mathrm{RM}}_{24})-h\left(\right|\▿f\left|\right)(v-f_y)\end{array}\right.---\left(13\right)$

其中，

表示卷积， $g\left(\right|\▿f\left|\right)=\exp (-{|\▿f|}^2/k^2),$ k为调整常量， $h\left(\right|\▿f\left|\right)=1-g\left(\right|\▿f\left|\right),$ ${\mathrm{NS}}_{24}=\left[\begin{array}{ccccc}0& 0& 1/12& 0& 0\\ 0& 0& 1/12& 0& 0\\ 1/12& 1/12& 1/3& 1/12& 1/12\\ 0& 0& 1/12& 0& 0\\ 0& 0& 1/12& 0& 0\end{array}\right],$ ${\mathrm{RM}}_{24}=\left[\begin{array}{ccccc}1/24& 1/24& 1/24& 1/24& 1/24\\ 1/24& 1/24& 1/24& 1/24& 1/24\\ 1/24& 1/24& 0& 1/24& 1/24\\ 1/24& 1/24& 1/24& 1/24& 1/24\\ 1/24& 1/24& 1/24& 1/24& 1/24\end{array}\right]$ 为卷积模板，该外力场基于卷积运算，计算速度快；迭代次数由用户设定，迭代计算到ENGGVF外力场V(x，y)＝[u(x，y)，v(x，y)]稳定为止；

三、在预处理后的图像上定义心脏左心室内膜初始化轮廓位置初始轮廓任意选取位于心脏左心室内膜范围内的一个圆；

四、对心脏左心室内膜进行分割

根据步骤二计算出的ENGGVF外力场，在初始化轮廓确定的情况下分割出内膜，在曲线演化过程中需添加椭圆形状能量约束，具体的分割过程如下：

1)构造椭圆形状能量约束场

定义轮廓线为曲线c(s)＝(x(s)，y(s))，其中s∈[0，1]，构造椭圆形状能量约束场，引入椭圆形状能量约束项：

$E_{\mathrm{ellipse}}=\frac{1}{2}\×{\∫}_0^1{((x\left(s\right)-x_c)\cos \left(\mathrm{\θ}\right)+(y\left(s\right)-y_c)\sin \left(\mathrm{\θ}\right)-r_1\cos (2\mathrm{\πs}-\mathrm{\θ}))}^2\mathrm{ds}$

(14)

$+\frac{1}{2}\×{\∫}_0^1{(-(x\left(s\right)-x_c)\sin \left(\mathrm{\θ}\right)+(y\left(s\right)-y_c)\cos \left(\mathrm{\θ}\right)-r_2\sin (2\mathrm{\πs}-\mathrm{\θ}))}^2\mathrm{ds}$

其中，(x_c，y_c)为椭圆中心，θ为椭圆偏转角，r₁，r₂分别为椭圆的两半径，[c_x，c_y，θ，r₁，r₂]可以通过最小二乘法拟合；

根据变分法原理，得到下列欧拉方程

$\left\{\begin{array}{c}{x}_i-x_c-r_1\cos (2\mathrm{\πi}/n-\mathrm{\θ})\cos \left(\mathrm{\θ}\right){+r}_2\sin (2\mathrm{\πi}/n-\mathrm{\θ})\sin \left(\mathrm{\θ}\right)=0\\ y_i-y_c-r_1\cos (2\mathrm{\πi}/n-\mathrm{\θ})\sin \left(\mathrm{\θ}\right)-r_2\sin (2\mathrm{\πi}/n-\mathrm{\θ})\cos \left(\mathrm{\θ}\right)=0\end{array}\right.---\left(15\right)$

其中，i＝0，1，…，n-1；从而得到椭圆形状能量约束场为：

F_ellipse＝(x_i，y_i)＝[x_c+r₁cos(2πi/n-θ)cos(θ)-r₂sin(2πi/n-θ)sin(θ)，

                                                                  (16)

y_c+r₁cos(2πi/n-θ)sin(θ)+r₂sin(2πi/n-θ)cos(θ)]

2)构造曲线迭代公式：

c(s)＝λ₁F_int+λ₂F_out+λ₃F_ellipse    (17)

其中，λ₁，λ₂和λ₃分别为内力场、外力场和圆形约束能量场的权重系数；内力场F_int＝αc_ss(s)-βc_ssss(s)，其中，α和β为弹性和刚性系数，c_ss(s)为曲线c(s)关于s的二阶导数，c_ssss(s)为曲线c(s)关于s的四阶导数；而外力场为步骤二计算出的ENGGVF外力场；椭圆形状能量约束场采用公式(16)；

3)将初始化轮廓作为曲线迭代公式的初始值，迭代计算上述曲线迭代公式(17)，得到一个稳定的解，即曲线收敛到左心室内膜的轮廓；

五、将内膜的最终分割轮廓结果定义为外膜的初始化轮廓位置；

六、将原始边缘图中内膜轮廓所包围区域的边缘强度置为0，这就抹平了左室壁内膜边缘及部分噪声，再采用这一改动的边缘图来重新计算ENGGVF外力场；

七、对心脏左心室外膜进行分割

在步骤六计算出的ENGGVF外力场作用下，根据步骤五定义的外膜初始化轮廓，对心脏左心室外膜进行分割，在曲线演化过程中需添加椭圆形状能量约束，具体的左心室外膜分割过程等同步骤四的左心室内膜分割过程。

有益效果

本发明提出了基于扩展邻域和噪声平滑的广义梯度矢量流ENGGVF，该外力场除了采用了权重因子和扩展邻域卷积运算外，还将拉普拉斯算子模板中加入噪声平滑模板，具有运算速度快、捕捉范围大、抗噪能力强，在弱边界保护和深度凹陷区域的分割上性能卓越。在左室壁内膜的分割而言，考虑到左室壁的近似为椭圆的特点，采用了椭圆形状约束的能量项，这种形状约束有利于克服由于图像灰度不均、乳突肌等而导致的局部极小。对于左室壁外膜的分割，利用内膜的分割结果初始化，即通过重新组合梯度分量来构造的外力场.这种外力场能有效克服原始梯度矢量流的不足，使得室壁外膜边缘很弱时也能得到保持。实验结果表明，该方法能准确地分割左室壁内、外膜。

附图说明

图1为本发明实施例某心脏核磁共振图像，其中1为左心室内膜，2为左心室外膜；

图2为图1预处理后的图像；

图3为该实施例图像左心室内膜的初始轮廓；

图4为该实施例图像左心室内膜的分割；

图5为形状约束对左心室内膜分割的影响，其中左图考虑了形状约束，右图没有考虑形状约束；

图6为形状约束对左心室外膜分割的影响，其中左图考虑了形状约束，右图没有考虑形状约束；

图7为该实施例图像左心室外膜的初始轮廓；

图8为该实施例图像左心室外膜的分割；

图9为一个心动周期内的21幅心脏核磁共振图像的左心室内外膜的分割结果。

具体实施方式

下面结合附图和具体实施方式对本发明作详细说明。

本实施方式具体实现了本发明提出的心脏核磁共振图像分割方法，包括以下步骤：

一、对取得的心脏核磁共振图像(如图1所示)进行高斯滤波预处理根据方程

对取得的心脏核磁共振图像进行高斯滤波预处理，式中I₀为输入的原始图像结构信息，G_σ为标准差为σ的二维高斯函数，

表示卷积运算。通过高斯滤波预处理，可以有效的滤去图像中的噪声，以便能够更好的实现心脏左心室内、外膜的分割。

取得的一幅心脏核磁共振图像见附图1，预处理后的图像见附图2。

二、在预处理后的图像上计算基于扩展邻域和噪声平滑的广义梯度矢量流的外力场，记为ENGGVF外力场，具体方法为：

1)定义ENGGVF外力场F_out的初始值

定义图像I(x，y)的边缘映射为f(x，y)，设f_x和f_y分别为边缘映射f沿x轴方向和y轴方向的一阶导数，边缘图的梯度向量构成了一个向量场V(x，y)＝[u(x，y)，v(x，y)]＝[f_x，f_y]，作为外力场F_out的初始值；

2)根据外力场F_out的初始值，计算ENGGVF外力场

ENGGVF外力场的迭代公式如下：

$\left\{\begin{array}{c}{u}_t=g\left(\right|\▿f\left|\right)(u\⊗{\mathrm{NS}}_{24}-u\⊗{\mathrm{RM}}_{24})-h\left(\right|\▿f\left|\right)(u-f_x)\\ v_t=g\left(\right|\▿f\left|\right)(v\⊗{\mathrm{NS}}_{24}-v\⊗{\mathrm{RM}}_{24})-h\left(\right|\▿f\left|\right)(v-f_y)\end{array}\right.---\left(18\right)$

其中，

表示卷积， $g\left(\right|\▿f\left|\right)=\exp (-{|\▿f|}^2/k^2),$ k为调整常量， $h\left(\right|\▿f\left|\right)=1-g\left(\right|\▿f\left|\right),$ ${\mathrm{NS}}_{24}=\left[\begin{array}{ccccc}0& 0& 1/12& 0& 0\\ 0& 0& 1/12& 0& 0\\ 1/12& 1/12& 1/3& 1/12& 1/12\\ 0& 0& 1/12& 0& 0\\ 0& 0& 1/12& 0& 0\end{array}\right],$ ${\mathrm{RM}}_{24}=\left[\begin{array}{ccccc}1/24& 1/24& 1/24& 1/24& 1/24\\ 1/24& 1/24& 1/24& 1/24& 1/24\\ 1/24& 1/24& 0& 1/24& 1/24\\ 1/24& 1/24& 1/24& 1/24& 1/24\\ 1/24& 1/24& 1/24& 1/24& 1/24\end{array}\right]$ 为卷积模板，该外力场基于卷积运算，计算速度快，而迭代次数由用户设定，迭代计算到ENGGVF外力场V(x，y)＝[u(x，y)，v(x，y)]稳定为止。

三、在预处理后的图像上定义心脏左心室内膜初始化轮廓位置

初始轮廓任意选取位于心脏左心室内膜范围内的一个圆，见附图3；

四、对心脏左心室内膜进行分割

根据步骤二计算出的ENGGVF外力场，在初始化轮廓确定的情况下分割出内膜，在曲线演化过程中需添加椭圆形状能量约束，具体的分割过程如下：

1)构造椭圆形状能量约束场：

为了克服血液的高速运动冲撞心肌壁造成的伪影(artifact)等引起的图像灰度不均，以及乳突肌干扰等对心脏MR图像的影响，我们既需要考虑曲线的光滑性，也需要考虑目标的整体形状。整体形状是一种全局性的约束，有利于克服图像中的噪声。但Snake模型内能只能约束曲线的连续性和光滑性等局部性质，且由于缺乏关于目标形状的全局信息而不能有效地刻画目标的形状。考虑到左室壁内、外膜的形状特点，本发明引入椭圆形状能量约束项，使得Snake轮廓在演化过程中其全局形状得到保持。该能量项如下：

$E_{\mathrm{ellipse}}=\frac{1}{2}\×{\∫}_0^1{((x\left(s\right)-x_c)\cos \left(\mathrm{\θ}\right)+(y\left(s\right)-y_c)\sin \left(\mathrm{\θ}\right)-r_1\cos (2\mathrm{\πs}-\mathrm{\θ}))}^2\mathrm{ds}$

(19)

$+\frac{1}{2}\×{\∫}_0^1{(-(x\left(s\right)-x_c)\sin \left(\mathrm{\θ}\right)+(y\left(s\right)-y_c)\cos \left(\mathrm{\θ}\right)-r_2\sin (2\mathrm{\πs}-\mathrm{\θ}))}^2\mathrm{ds}$

其中，(x_c，y_c)为椭圆中心，θ为椭圆偏转角，r₁，r₂分别为椭圆的两半径，[c_x，c_y，θ，r₁，r₂]可以通过最小二乘法拟合。它们是随着Snake曲线演化而动态变化的。E_ellipse这一能量就度量了Snake轮廓上的点与中心为(x_c，y_c)，偏转角为θ，两半径为r₁，r₂的椭圆之间的差异。当Snake轮廓不受外力作用时，该能量项将使Snake轮廓保持为椭圆。在分割左心室内膜过程中当Snake曲线演化到伪影和乳突肌时，由于受到椭圆形状能量约束的限制，曲线能绕过伪影和乳突肌向着我们需要的目标特征继续演化。根据变分法原理，公式(19)对应的欧拉方程为

$\left\{\begin{array}{c}{x}_i-x_c-r_1\cos (2\mathrm{\πi}/n-\mathrm{\θ})\cos \left(\mathrm{\θ}\right){+r}_2\sin (2\mathrm{\πi}/n-\mathrm{\θ})\sin \left(\mathrm{\θ}\right)=0\\ y_i-y_c-r_1\cos (2\mathrm{\πi}/n-\mathrm{\θ})\sin \left(\mathrm{\θ}\right)-r_2\sin (2\mathrm{\πi}/n-\mathrm{\θ})\cos \left(\mathrm{\θ}\right)=0\end{array}\right.---\left(20\right)$

其中，i＝0，1，…，n-1；从而得到椭圆形状能量约束场：

F_ellipse＝(x_i，y_i)＝[x_c+r₁cos(2πi/n-θ)cos(θ)-r₂sin(2πi/n-θ)sin(θ)，

                                                           (21)

y_c+r₁cos(2πi/n-θ)sin(θ)+r₂sin(2πi/n-θ)cos(θ)]

2)构造曲线迭代公式：

轮廓线用曲线c(s)＝(x(s)，y(s))(s∈[0，1])来定义，构造添加了椭圆形状能量约束的曲线迭代公式：

c(s)＝λ₁F_int+λ₂F_out+λ₃F_ellipse                 (22)

其中，λ₁，λ₂和λ₃分别为内力场、外力场和椭圆形状能量约束场的权重系数，内力场F_int＝αc_ss(s)-βc_ssss(s)，其中，α和β为弹性和刚性系数，c_ss(s)为曲线c(s)关于s的二阶导数，c_ssss(s)为曲线c(s)关于s的四阶导数；而外力场为步骤二计算出的ENGGVF外力场；椭圆形状能量约束场采用公式(21)。

在分割内膜时，如果不采用椭圆形状能量约束，Snake轮廓易受到乳突肌和伪影的干扰而陷入局部极小。采用形状约束后Snake轮廓能绕过伪影和乳突肌收敛到我们需要的目标边界，得到较好的分割结果。另外，在分割外膜时由于左心室与右心室及周围其他组织如肝脏等灰度非常接近易形成弱边界，并且外力场不够完美，若不采用形状约束，由于图像的梯度力在低对比度区域和弱边界处太小，会出现变形曲线泄露的现象。在采用全局形状约束后，可以阻止Snake轮廓从低对比度区域或弱边界区域泄露。

形状约束对左心室内膜分割的影响见图5，对左心室外膜分割的影响见图6。

3)将初始化轮廓作为曲线迭代公式的初始值，迭代计算上述曲线迭代公式(22)，得到一个稳定的解，即曲线收敛到左心室内膜的轮廓。

心脏左心室内膜的分割结果见附图3。

五、采用内膜的分割结果初始化，将内膜的最终分割轮廓结果定义为外膜的初始化轮廓位置，见附图7。

六、将原始边缘图中内膜轮廓所包围区域的边缘强度置为0，这就抹平了左室壁内膜边缘及部分噪声，再采用这一改动的边缘图来重新计算ENGGVF外力场。

七、对心脏左心室外膜进行分割

在步骤六计算出的ENGGVF外力场作用下，根据步骤五定义的外膜初始化轮廓，对心脏左心室外膜进行分割，在曲线演化过程中需添加椭圆形状能量约束，具体的左心室外膜分割过程等同步骤四的左心室内膜分割过程；

心脏左心室外膜的分割结果见附图8。

通过在一个心动周期内的一套心脏核磁共振图像上，使用本实施方式中所述方法，见附图9，验证上述分割策略，并与手工分割的结果进行定量比较。这里所用的MR图像由SIEMENS 1.5T临床系统产生，成像参数如下：原始图像尺寸192×156，切片厚度8mm，重复时间(TR)＝29.16，回波时间(TE)＝1.08，分辨率1.82×1.82，回转角(flip angle)＝50，视野(FOV)＝81.25。实验中使用的参数为α＝0.1，β＝0，λ₁＝0.1，λ₂＝0.1，λ₃＝0.3，k＝0.1.计算环境为Matlab7.1，CPU 3.39G，RAM 1.0G，Windows XP Professional。

我们对分割结果与手工分割结果进行比较，采用平均绝对距离(meanabsolute distance，简称MAD)度量二者之间的差异。设Snake轮廓为S，手工分割结果为M，则

$\mathrm{mad}(S,M)=0.5(\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}d(s_i,M)+\frac{1}{k}\underset{j=1}{\overset{k}{\mathrm{\Σ}}}d(m_j,S))---\left(23\right)$

其中，S＝{s1，...，sn}，M＝{m1，...，mk}分别表示Snake轮廓和手工轮廓上的点，

对于整套图像而言，左室壁内膜的平均MAD值为0.45像素，基本上与手工分割结果相同；外膜的平均MAD值为1.26像素，与手工分割结果非常接近。

Claims (1)

1.一种心脏核磁共振图像的快速分割方法，其特征在于包括以下步骤：

一、对取得的心脏核磁共振图像进行高斯滤波预处理

根据方程对取得的心脏核磁共振图像进行高斯滤波预处理，式中I₀为输入的原始图像结构信息，G_σ为标准差为σ的二维高斯函数，

表示卷积运算；

二、在预处理后的图像上计算基于扩展邻域和噪声平滑的广义梯度矢量流的外力场，记为ENGGVF外力场，具体方法为：

1)定义ENGGVF外力场F_out的初始值

定义图像I(x，y)的边缘映射为f(x，y)，设f_x和f_y分别为边缘映射f沿x轴方向和y轴方向的一阶导数，边缘图的梯度向量

构成了一个向量场V(x，y)＝[u(x，y)，v(x，y)]＝[f_x，f_y]，作为外力场F_out的初始值；

2)根据外力场F_out的初始值，计算ENGGVF外力场

ENGGVF外力场的迭代公式如下：

$\left\{\begin{array}{c}{u}_t=g\left(\right|\▿f\left|\right)(u\⊗{\mathrm{NS}}_{24}-u\⊗{\mathrm{RM}}_{24})-h\left(\right|\▿f\left|\right)(u-f_x)\\ v_t=g\left(\right|\▿f\left|\right)(v\⊗{\mathrm{NS}}_{24}-v\⊗{\mathrm{RM}}_{24})-h\left(\right|\▿f\left|\right)(v-f_y)\end{array}\right.---\left(1\right)$

其中，

表示卷积， $g\left(\right|\▿f\left|\right)=\exp (-{|\▿f|}^2/k^2),$ k为调整常量， $h\left(\right|\▿f\left|\right)=1-g\left(\right|\▿f\left|\right),$ ${\mathrm{NS}}_{24}=\left[\begin{array}{ccccc}0& 0& 1/12& 0& 0\\ 0& 0& 1/12& 0& 0\\ 1/12& 1/12& 1/3& 1/12& 1/12\\ 0& 0& 1/12& 0& 0\\ 0& 0& 1/12& 0& 0\end{array}\right],$ ${\mathrm{RM}}_{24}=\left[\begin{array}{ccccc}1/24& 1/24& 1/24& 1/24& 1/24\\ 1/24& 1/24& 1/24& 1/24& 1/24\\ 1/24& 1/24& 0& 1/24& 1/24\\ 1/24& 1/24& 1/24& 1/24& 1/24\\ 1/24& 1/24& 1/24& 1/24& 1/24\end{array}\right]$ 为卷积模板，迭代次数由用户设定，迭代计算到ENGGVF外力场V(x，y)＝[u(x，y)，v(x，y)]稳定为止；

三、在预处理后的图像上定义心脏左心室内膜初始化轮廓位置

初始轮廓任意选取位于心脏左心室内膜范围内的一个圆；

四、对心脏左心室内膜进行分割

根据步骤二计算出的ENGGVF外力场，在初始化轮廓确定的情况下分割出内膜，在曲线演化过程中需添加椭圆形状能量约束，具体的分割过程如下：

1)构造椭圆形状能量约束场

定义轮廓线为曲线c(s)＝(x(s)，y(s))，其中s∈[0，1]，构造椭圆形状能量约束场，引入椭圆形状能量约束项：

$E_{\mathrm{ellipse}}=\frac{1}{2}\×{\∫}_0^1{((x\left(s\right)-x_c)\cos \left(\mathrm{\θ}\right)+(y\left(s\right)-y_c)\sin \left(\mathrm{\θ}\right)-r_1\cos (2\mathrm{\πs}-\mathrm{\θ}))}^2\mathrm{ds}$

(2)

$+\frac{1}{2}\×{\∫}_0^1{(-(x\left(s\right)-x_c)\sin \left(\mathrm{\θ}\right)+(y\left(s\right)-y_c)\cos \left(\mathrm{\θ}\right)-r_2\sin (2\mathrm{\πs}-\mathrm{\θ}))}^2\mathrm{ds}$

其中，(x_c，y_c)为椭圆中心，θ为椭圆偏转角，r₁，r₂分别为椭圆的两半径，[c_x，c_y，θ，r₁，r₂]可以通过最小二乘法拟合；

根据变分法原理，得到下列离散形式的欧拉方程

$\left\{\begin{array}{c}{x}_i-x_c-r_1\cos (2\mathrm{\πi}/n-\mathrm{\θ})\cos \left(\mathrm{\θ}\right){+r}_2\sin (2\mathrm{\πi}/n-\mathrm{\θ})\sin \left(\mathrm{\θ}\right)=0\\ y_i-y_c-r_1\cos (2\mathrm{\πi}/n-\mathrm{\θ})\sin \left(\mathrm{\θ}\right)-r_2\sin (2\mathrm{\πi}/n-\mathrm{\θ})\cos \left(\mathrm{\θ}\right)=0\end{array}\right.---\left(3\right)$

其中，i＝0，1，…，n-1；从而得到椭圆形状能量约束场为：

F_ellipse＝(x_i，y_i)＝[x_c+r₁cos(2πi/n-θ)cos(θ)-r₂sin(2πi/n-θ)sin(θ)，

                                                    (4)

y_c+r₁cos(2πi/n-θ)sin(θ)+r₂sin(2πi/n-θ)cos(θ)]

2)构造曲线迭代公式为：

c(s)＝λ₁F_int+λ₂F_out+λ₃F_ellipse                            (5)

其中，λ₁，λ₂和λ₃分别为内力场、ENGGVF外力场和椭圆形状能量约束场的权重系数；内力场F_int＝αc_ss(s)-βc_ssss(s)，其中，α和β为弹性和刚性系数，c_ss(s)为曲线c(s)关于s的二阶导数，c_ssss(s)为曲线c(S)关于s的四阶导数；而外力场为步骤二计算出的ENGGVF外力场；椭圆形状约束能量场采用公式(4)；

3)将初始化轮廓作为曲线迭代公式的初始值，迭代计算上述曲线迭代公式(5)，得到一个稳定的解，即曲线收敛到左心室内膜的轮廓；

五、将心脏左心室内膜的最终分割轮廓结果定义为心脏左心室外膜的初始化轮廓位置；

六、将原始边缘图中内膜轮廓所包围区域的边缘强度置为0，这就抹平了左室壁内膜边缘及部分噪声，再采用这一改动的边缘图来重新计算ENGGVF外力场；

七、对心脏左心室外膜进行分割

在步骤六计算出的ENGGVF外力场作用下，根据步骤五定义的外膜初始化轮廓，对心脏左心室外膜进行分割，在曲线演化过程中需添加椭圆形状能量约束，具体的左心室外膜分割过程等同步骤四的左心室内膜分割过程。

Images