CN102289814B - 一种心脏核磁共振图像分割方法 - Google Patents

Abstract

本发明涉及一种心脏核磁共振图像分割方法，包括以下步骤：一、对图像进行高斯滤波预处理；二、计算边缘保持法向梯度矢量流的外力场；三、定义心脏左心室内膜初始化轮廓位置；四、在曲线演化过程中添加圆形能量约束，分割出内膜；五、将内膜的最终分割轮廓结果定义为外膜的初始化轮廓位置；六、将原始边缘图中内膜轮廓所包围区域的边缘强度置为0，重新计算外力场；七、在曲线演化过程中添加圆形能量约束，分割出外膜。本发明捕捉范围大、抗噪能力强、且对弱边界泄漏有更好的鲁棒性，能准确地分割左室壁内、外膜。

Description

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

技术领域

本发明涉及一种图像分割方法，特别涉及一种心脏核磁共振图像分割方法，属于图像处理技术领域。

背景技术

心脏MRI(magnetic resonance imaging)能够提供高解析的、较好的软组织对比图像，对心脏的解剖结构和功能进行准确的描述，是当前医学图像分析领域的研究热点之一。由于左心室负责供血，目前的研究主要集中在左心室的功能上。为了充分利用图像中的解剖信息，为临床诊断提供量化、直观的参考，首先必须分割出左室壁的内、外膜。然而，由于心脏的运动和血液的高速流动，图像受噪声干扰，使得心脏MR图像的分割仍是一个值得深入研究的问题。

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

目前，心脏左心室内、外膜分割存在以下难点：首先，心脏MR图像在成像过程中由于受血液流动的影响，会在血池中产生伪影，使得图像灰度不均；其次，目标边界往往受到乳突肌，呼吸等因素的影响而变得模糊不清甚至边界断裂；另外，由于左心室与右心室及周围其他组织如肝脏等灰度非常接近，形成弱边界，这时基于主动轮廓模型的方法分割左心室外膜时往往发生泄露。现有的方法对于这些问题没有提出很好的解决方案。

发明内容

本发明的目的是针对心脏左心室内、外膜分割存在的难点，提出一种高效、鲁棒、准确的图像分割方法，分割左室壁内、外膜。

本发明的方法是基于广义的边缘保持法向梯度矢量流(epnGGVF)模型提出的，其思想如下：

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)是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(2\right)$

当方程式(2)的解收敛时，就得到了待分割目标的轮廓。这时，可以将Snake轮廓的运动过程看成其内、外力的平衡过程，αc_ss(s)-βc_ssss(s)称为Snake模型的内力，称为其外力。外力在Snake模型的演化中起决定性作用，对外力的研究是Snake模型研究的一个重要方面。由于式(2)定义的Snake模型外力

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

来约束动态轮廓线，并将其定义为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模型扩展了边缘映射的梯度向量，同时通过一个各向同性的扩散过程来抑制噪声。GVF的扩散过程依赖于拉普拉斯算子，对图像f(x，y)而言，其拉普拉斯算子为

Δf＝f_xx+f_yy                (5)

其中，f_xx为沿x轴方向的二阶导数，f_yy为沿y轴方向的二阶导数。这是在迪卡尔坐标系下的表达方式，那么在Gauge坐标系下则可写为

Δf＝f_TT+f_NN                (6)

其中，f_TT为沿切线方向的二阶导数，f_NN为沿法线方向的二阶导数。尽管GVF已经取得了巨大成功，但仍有改进的余地。最近提出的NGVF模型是对GVF的一种改进，它从插值的角度出发，把GVF看成一个插值的过程，认为从插值的角度来说，f_NN是最好的插值算子，Δf其次，f_TT是最差的，因此提出了NGVF模型：

$\left\{\begin{array}{c}{u}_t=\mathrm{\μ}{u}_{\mathrm{NN}}-{|\▿f|}^2(u-f_x)\\ v_t=\mathrm{\μ}{v}_{\mathrm{NN}}-{|\▿f|}^2(v-f_y)\end{array}\right.---\left(7\right)$

其中，μ为权重系数，

是边缘图的梯度向量，u_t和v_t分别为u和v对时间的偏导数，u_NN和v_NN分别为u和v沿法线方向的二阶导数，f_x和f_y分别为边缘映射f沿x轴方向和y轴方向的一阶导数。然而，从GVF的扩散过程来分析，GVF并不是一个插值的过程，而是一个逼近的过程，当目标为细长的凹陷区域时，GVF仍难以进入。同时扩散分解的理论告诉我们，沿切线方向的扩散有助于去除噪声并保持边缘。基于这样的考虑，本发明提出了一种新的外部力场：广义的边缘保持法向梯度矢量流(epnGGVF)模型，并采用这一外力模型来分割左室壁。

epnGGVF模型在NGVF的基础上，引入了权重因子和边缘保持项，使得在边缘附近法线方向的扩散减弱，平坦区域的扩散增强，这样有利于保持边缘，去除噪声并扩大捕捉范围。新的能量泛函如下：

$\mathrm{\ϵ}=\∫\∫g\left(\right|\▿f\left|\right){|\▿V|}^2+h\left(\right|\▿f\left|\right)\left(m\right|{J_{V}p|}^2+{|V-\▿f|}^2)---\left(8\right)$

其中，

$\left\{\begin{array}{c}g\left(\right|\▿f\left|\right)=\exp (-{|\▿f|}^2/k^2)\\ h\left(\right|\▿f\left|\right)=1-g\left(\right|\▿f\left|\right)\end{array}\right.---\left(9\right)$

m为权重系数，|J_Vp|²为边缘保持项，

J_V是关于力场V的雅可比矩阵。使用变分原理，得到下列欧拉方程：

$\left\{\begin{array}{c}{u}_t=g\left(\right|\▿f\left|\right)\mathrm{\Δu}-h\left(\right|\▿f\left|\right)(u-f_x-{\mathrm{mp}}_1^2{u}_{\mathrm{xx}}-{\mathrm{mp}}_2^2{u}_{\mathrm{yy}}-2{\mathrm{mp}}_1{p}_2{u}_{\mathrm{xy}})\\ v_t=g\left(\right|\▿f\left|\right)\mathrm{\Δv}-h\left(\right|\▿f\left|\right)(v-f_y-{\mathrm{mp}}_1^2{v}_{\mathrm{xx}}-{\mathrm{mp}}_2^2{v}_{\mathrm{yy}}-2{\mathrm{mp}}_1{p}_2{v}_{\mathrm{xy}})\end{array}\right.---\left(10\right)$

将拉普拉斯算子Δu，Δv用更好的插值算子u_NN，v_NN代替，则epnGGVF模型可通过解下列欧拉方程获得：

$\left\{\begin{array}{c}{u}_t=g\left(\right|\▿f\left|\right)u_{\mathrm{NN}}-h\left(\right|\▿f\left|\right)(u-f_x-{\mathrm{mp}}_1^2{u}_{\mathrm{xx}}-{\mathrm{mp}}_2^2{u}_{\mathrm{yy}}-2{\mathrm{mp}}_1{p}_2{u}_{\mathrm{xy}})\\ v_t=g\left(\right|\▿f\left|\right)v_{\mathrm{NN}}-h\left(\right|\▿f\left|\right)(v-f_y-{\mathrm{mp}}_1^2{v}_{\mathrm{xx}}-{\mathrm{mp}}_2^2{v}_{\mathrm{yy}}-2{\mathrm{mp}}_1{p}_2{v}_{\mathrm{xy}})\end{array}\right.---\left(11\right)$

这里，g(.)是关于

的单调递减函数，相应地，h(.)是关于

的单调递增函数，有利于力场根据图像的结构信息来调整扩散过程并顺利进入细长的凹陷区域。边缘保持项可以沿着边缘方向平滑力场而不穿过力场，有利于防止弱边界泄露，因而具有较好的边缘保持功能。此外，m是一个可为特定应用而调整的常量，它决定着要保持的边缘对比度。当需要保持的边界是弱边缘，这个值很小，反之亦然。

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

一、根据方程

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

表示卷积运算；

二、计算边缘保持法向梯度矢量流的外力场，记为epnGGVF外力场，具体方法为：

1)定义图像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]，作为epnGGVF外力场F_out的初始值；

2)按照epnGGVF外力场的迭代公式，迭代计算epnGGVF外力场，迭代次数由用户设定，迭代计算到epnGGVF外力场V(x，y)＝[u(x，y)，v(x，y)]稳定为止；

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

$\left\{\begin{array}{c}{u}_t=g\left(\right|\▿f\left|\right)u_{\mathrm{NN}}-h\left(\right|\▿f\left|\right)(u-f_x-{\mathrm{mp}}_1^2{u}_{\mathrm{xx}}-{\mathrm{mp}}_2^2{u}_{\mathrm{yy}}-2{\mathrm{mp}}_1{p}_2{u}_{\mathrm{xy}})\\ v_t=g\left(\right|\▿f\left|\right)v_{\mathrm{NN}}-h\left(\right|\▿f\left|\right)(v-f_y-{\mathrm{mp}}_1^2{v}_{\mathrm{xx}}-{\mathrm{mp}}_2^2{v}_{\mathrm{yy}}-2{\mathrm{mp}}_1{p}_2{v}_{\mathrm{xy}})\end{array}\right.---\left(12\right)$

其中，m为权重系数， $p=[{-f}_y/|\▿f|,f_x/|\▿f\left|\right]=[p_1,p_2],$ $g\left(\right|\▿f\left|\right)=\exp (-{|\▿f|}^2/k^2),$ k为调整常量，

u_t和v_t分别为u和v对时间的偏导数，u_NN和v_NN分别为u和v沿法线方向的二阶导数，u_xx和v_xx分别为u和v沿x轴方向的二阶导数，u_yy和v_yy分别为u和v沿y轴方向的二阶导数，u_xy和v_xy分别为u和v沿x轴y轴方向的二阶混合偏导；

三、定义心脏左心室内膜初始化轮廓位置，初始轮廓任意选取位于内膜范围内的一个圆；

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

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

$E_{\mathrm{circle}}=\frac{\mathrm{\λ}}{2}{\∫}_0^1{\left((R\left(s\right)-\stackrel{\‾}{R})\right)}^2\mathrm{ds}---\left(13\right)$

其中， $R\left(s\right)=\sqrt{{(x\left(s\right)-x_c)}^2+{(y\left(s\right)-y_c)}^2},$ $x_c={\∫}_0^{1}x\left(r\right)\mathrm{dr},$ $y_c={\∫}_0^{1}y\left(r\right)\mathrm{dr},$ $\stackrel{\‾}{R}={\∫}_0^{1}R\left(s\right)\mathrm{ds},$ λ为权重因子，将(x_c，y_c)看作Snake轮廓线的质心；

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

$\left\{\begin{array}{c}\mathrm{\λ}(x\left(s\right)-x_c-\stackrel{\‾}{R}\cos \left(2\mathrm{\πs}\right))=0\\ \mathrm{\λ}(y\left(s\right)-y_c-\stackrel{\‾}{R}\sin \left(2\mathrm{\πs}\right))=0\end{array}\right.---\left(14\right)$

写成离散形式为

$\left\{\begin{array}{c}\mathrm{\λ}(x_i-x_c-\stackrel{\‾}{R}\cos (2\mathrm{\πi}/n))=0\\ \mathrm{\λ}(y_i-y_c-\stackrel{\‾}{R}\sin (2\mathrm{\πi}/n))=0\end{array}\right.---\left(15\right)$

其中， $\stackrel{\‾}{R}=\frac{1}{n}\underset{i=0}{\overset{n-1}{\mathrm{\Σ}}}{R}_i,$ $R_i=\sqrt{{(x_i-x_c)}^2+{(y_i-y_c)}^2},$ $x_c=\frac{1}{n}\underset{i=0}{\overset{n-1}{\mathrm{\Σ}}}{x}_i,$ $y_c=\frac{1}{n}\underset{i=0}{\overset{n-1}{\mathrm{\Σ}}}{y}_i,$ i＝0，…，n-1；从而得到圆形能量约束场：

$F_{\mathrm{circle}}=(x_i,y_i)=[x_c+\stackrel{\‾}{R}\cos (2\mathrm{\πi}/n),y_c+\stackrel{\‾}{R}\sin (2\mathrm{\πi}/n)]---\left(16\right)$

2)构造添加了圆形能量约束的曲线迭代公式：

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

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

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

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

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

七、根据步骤六计算出的epnGGVF外力场，在初始化轮廓确定的情况下分割出外膜，在曲线演化过程中需添加圆形能量约束，具体的左心室外膜分割过程等同步骤四的左心室内膜分割过程。

有益效果

本发明提出了广义的边缘保持法向梯度矢量流外力模型epnGGVF，该外力场除了采用了权重因子外，还加入了新的边缘保持项，具有捕捉范围大、抗噪能力强、且对弱边界泄漏有更好的鲁棒性。在左室壁内膜的分割而言，考虑到左室壁的近似为圆形的特点，采用了圆形约束的能量项，这种形状约束有利于克服由于图像灰度不均、乳突肌等而导致的局部极小。对于左室壁外膜的分割，利用内膜的分割结果初始化，即通过重新组合梯度分量来构造的外力场.这种外力场能有效克服原始梯度矢量流的不足，使得室壁外膜边缘很弱时也能得到保持。实验结果表明，该方法能准确地分割左室壁内、外膜。

具体实施方式

下面详细说明本发明的优选实施方式。

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

一、根据方程

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

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

二、计算epnGGVF外力场：

1)定义图像I(x，y)的边缘映射f(x，y)，使其在图像边缘附近取值较大，而在均匀区域内部取值较小；边缘图的梯度向量

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

2)按照epnGGVF外力场的迭代公式，迭代计算epnGGVF外力场，迭代次数由用户设定，迭代计算到epnGGVF外力场V(x，y)＝[u(x，y)，v(x，y)]稳定为止；

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

$\left\{\begin{array}{c}{u}_t=g\left(\right|\▿f\left|\right)u_{\mathrm{NN}}-h\left(\right|\▿f\left|\right)(u-f_x-{\mathrm{mp}}_1^2{u}_{\mathrm{xx}}-{\mathrm{mp}}_2^2{u}_{\mathrm{yy}}-2{\mathrm{mp}}_1{p}_2{u}_{\mathrm{xy}})\\ v_t=g\left(\right|\▿f\left|\right)v_{\mathrm{NN}}-h\left(\right|\▿f\left|\right)(u-f_y-{\mathrm{mp}}_1^2{v}_{\mathrm{xx}}-{\mathrm{mp}}_2^2{v}_{\mathrm{yy}}-2{\mathrm{mp}}_1{p}_2{v}_{\mathrm{xy}})\end{array}\right.---\left(18\right)$

其中，m为权重系数， $p=[{-f}_y/|\▿f|,f_x/|\▿f\left|\right]=[p_1,p_2],$ $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).$

三、定义心脏左心室内膜初始化轮廓位置，初始轮廓可随意选取位于内膜范围内的一个圆。

四、在初始化轮廓确定的情况下分割出内膜：

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

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

$E_{\mathrm{circle}}=\frac{\mathrm{\λ}}{2}{\∫}_0^1{\left((R\left(s\right)-\stackrel{\‾}{R})\right)}^2\mathrm{ds}---\left(19\right)$

其中， $R\left(s\right)=\sqrt{{(x\left(s\right)-x_c)}^2+{(y\left(s\right)-y_c)}^2},$ $x_c={\∫}_0^{1}x\left(r\right)\mathrm{dr},$ $y_c={\∫}_0^{1}y\left(r\right)\mathrm{dr},$ $\stackrel{\‾}{R}={\∫}_0^{1}R\left(s\right)\mathrm{ds},$ λ为权重因子，(x_c，y_c)可以认为是Snake轮廓线的质心，它和

是随着Snake曲线演化而动态变化的。E_circle这一能量就度量了Snake轮廓上的点与圆心为(x_c，y_c)，半径为

的圆之间的差异。当Snake轮廓不受外力作用时，该能量项将使Snake轮廓保持为圆形。在分割左心室内膜过程中当Snake曲线演化到伪影和乳突肌时，由于受到圆形约束的限制，曲线能绕过伪影和乳突肌向着我们需要的目标特征继续演化。根据变分法原理，公式(17)对应的欧拉方程为

$\left\{\begin{array}{c}\mathrm{\λ}(x\left(s\right)-x_c-\stackrel{\‾}{R}\cos \left(2\mathrm{\πs}\right))=0\\ \mathrm{\λ}(y\left(s\right)-y_c-\stackrel{\‾}{R}\sin \left(2\mathrm{\πs}\right))=0\end{array}\right.---\left(20\right)$

写成离散形式为

$\left\{\begin{array}{c}\mathrm{\λ}(x_i-x_c-\stackrel{\‾}{R}\cos (2\mathrm{\πi}/n))=0\\ \mathrm{\λ}(y_i-y_c-\stackrel{\‾}{R}\sin (2\mathrm{\πi}/n))=0\end{array}\right.---\left(21\right)$

其中， $\stackrel{\‾}{R}=\frac{1}{n}\underset{i=0}{\overset{n-1}{\mathrm{\Σ}}}{R}_i,$ $R_i=\sqrt{{(x_i-x_c)}^2+{(y_i-y_c)}^2},$ $x_c=\frac{1}{n}\underset{i=0}{\overset{n-1}{\mathrm{\Σ}}}{x}_i,$ $y_c=\frac{1}{n}\underset{i=0}{\overset{n-1}{\mathrm{\Σ}}}{y}_i,$ i＝0，…，n-1。将公式(19)放在时间演化框架下求解并采用半隐式格式，得到

$\left\{\begin{array}{c}\frac{x_i^{t+1}-x_i^t}{\mathrm{\Δt}}+\mathrm{\λ}{x}_i^{t+1}-\mathrm{\λ}(x_c^t+{\stackrel{\‾}{R}}^t\cos (2\mathrm{\πi}/n))=0\\ \frac{y_i^{t+1}-y_i^t}{\mathrm{\Δt}}+\mathrm{\λ}{y}_i^{t+1}-\mathrm{\λ}(y_c^t+{\stackrel{\‾}{R}}^t\sin (2\mathrm{\πi}/n))=0\end{array}\right.---\left(22\right)$

2)构造添加了圆形能量约束的曲线迭代公式：

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

c(s)＝λ₁F_int+λ₂F_out+λ₃F_circle      (23)

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

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

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

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

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

七、根据步骤六计算出的epnGGVF外力场，在初始化轮廓确定的情况下分割出外膜，在曲线演化过程中需添加圆形能量约束，具体的左心室外膜分割过程等同步骤四的左心室内膜分割过程。

通过在一套心脏MR图像上使用本实施方式中所述方法验证上述分割策略，并与手工分割的结果进行定量比较。这里所用的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.3，μ＝0.15，k＝0.1，m＝0.01.计算环境为Matlab7.1，CPU 3.39G，RAM 1.0G，Windows XP Professional。

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

$\mathrm{mad}(S,M)=0.4(\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(24\right)$

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

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

Claims (1)

1.一种心脏核磁共振图像分割方法，包括以下步骤：

一、根据方程

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

表示卷积运算；

二、计算边缘保持法向梯度矢量流的外力场，记为epnGGVF外力场，具体方法为：

1)定义图像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]，作为epnGGVF外力场的初始值；

2)按照epnGGVF外力场的迭代公式，迭代计算epnGGVF外力场，迭代次数由用户设定，迭代计算到epnGGVF外力场V(x，y)=[u(x，y)，v(x，y)]稳定为止；

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

$\left\{\begin{array}{c}{u}_t=g\left(\right|\▿f\left|\right)u_{\mathrm{NN}}-h\left(\right|\▿f\left|\right)(u-f_x-{\mathrm{mp}}_1^2{u}_{\mathrm{xx}}-{\mathrm{mp}}_2^2{u}_{\mathrm{yy}}-{2\mathrm{mp}}_1{p}_2{u}_{\mathrm{xy}})\\ v_t=g\left(\right|\▿f\left|\right)v_{\mathrm{NN}}-h\left(\right|\▿f\left|\right)(v-f_y-{\mathrm{mp}}_1^2{v}_{\mathrm{xx}}-{\mathrm{mp}}_2^2{v}_{\mathrm{yy}}-{2\mathrm{mp}}_1{p}_2{v}_{\mathrm{xy}})\end{array}\right.---\left(12\right)$

其中，m为权重系数， $p=[-f_y/|\▿f|,f_x/|\▿f\left|\right]=[p_1,p_2],$ $g\left(\right|\▿f\left|\right)=\exp (-{|\▿f|}^2/k^2),$ k为调整常量，

u_t和v_t分别为u和v对时间的偏导数，u_NN和v_NN分别为u和v沿法线方向的二阶导数，u_xx和v_xx分别为u和v沿x轴方向的二阶导数，u_yy和v_yy分别为u和v沿y轴方向的二阶导数，u_xy和v_xy分别为u和v沿x轴y轴方向的二阶混合偏导；

三、定义心脏左心室内膜初始化轮廓位置，初始轮廓任意选取位于内膜范围内的一个圆；

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

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

$E_{\mathrm{circle}}=\frac{\mathrm{\λ}}{2}{\∫}_0^1{\left((R\left(s\right)-\stackrel{\‾}{R})\right)}^2\mathrm{ds}---\left(13\right)$

其中， $R\left(s\right)=\sqrt{{(x\left(s\right)-x_c)}^2+{(y\left(s\right)-y_c)}^2},$ $x_c={\∫}_0^{1}x\left(r\right)\mathrm{dr},$ $y_c={\∫}_0^{1}y\left(r\right)\mathrm{dr},$ $\stackrel{\‾}{R}={\∫}_0^{1}R\left(s\right)\mathrm{ds},$ λ为权重因子，将(x_c，y_c)看作Snake轮廓线的质心；

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

$\left\{\begin{array}{c}\mathrm{\λ}(x\left(s\right)-x_c-\stackrel{\‾}{R}\cos \left(2\mathrm{\πs}\right))=0\\ \mathrm{\λ}(y\left(s\right)-y_c-\stackrel{\‾}{R}\sin \left(2\mathrm{\πs}\right))=0\end{array}\right.---\left(14\right)$

写成离散形式为

$\left\{\begin{array}{c}\mathrm{\λ}(x_i-x_c-\stackrel{\‾}{R}\cos (2\mathrm{\πi}/n))=0\\ \mathrm{\λ}(y_i-y_c-\stackrel{\‾}{R}\sin (2\mathrm{\πi}/n))=0\end{array}\right.---\left(15\right)$

其中， $\stackrel{\‾}{R}=\frac{1}{n}\underset{i=0}{\overset{n-1}{\mathrm{\Σ}}}{R}_i,$ $R_i=\sqrt{{(x_i-x_c)}^2+{(y_i-y_c)}^2},$ $x_c=\frac{1}{n}\underset{i=0}{\overset{n-1}{\mathrm{\Σ}}}{x}_i,$ $y_c=\frac{1}{n}\underset{i=0}{\overset{n-1}{\mathrm{\Σ}}}{y}_i,$ i＝0，…，n-1；从而得到圆形能量约束场：

$F_{\mathrm{circle}}=(x_i,y_i)=[x_c+\stackrel{\‾}{R}\cos (2\mathrm{\πi}/n),y_c+\stackrel{\‾}{R}\sin (2\mathrm{\πi}/n)]---\left(16\right)$

2)构造添加了圆形能量约束的曲线迭代公式：

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

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

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

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

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

七、根据步骤六计算出的epnGGVF外力场，在初始化轮廓确定的情况下分割出外膜，在曲线演化过程中需添加圆形能量约束，具体的左心室外膜分割过程等同步骤四的左心室内膜分割过程。