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

摘要

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

说明书

技术领域

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

背景技术

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

近年来，对心脏MR(magnetic resonance)图像的分割国内外都有广泛研究，这 些方法大致可以分为基于形态学的方法，基于模糊聚类的方法，基于模板的方 法以及基于主动轮廓模型的方法等.主动轮廓模型能够将有关目标形状的先验知 识和来自图像的知识融入一个统一的过程中，是当前图像分割领域的热点方法， 也是心脏MR图像分割中的主流方法，在国内外都有广泛研究。Hong等人采用 基于Lagrange动力学的B样条Snake模型来提取左室壁内膜，Makowski等人采 用气球Snake模型来分割左室壁内膜，设计了专门的方法来解决轮廓缠绕问题。 Nguyen等人对传统Snake，GVF(gradient vector flow)Snake和气球Snake模型分 割左室壁内膜的结果作了比较，并与手工勾勒的轮廓进行对比验证，其中，GVF  Snake模型性能最好。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}}=\int \frac{1}{2}(\alpha {\left|c_s\right|}^2+\beta {\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{\alpha c}}_{\mathrm{ss}}\left(s\right)-{\mathrm{\beta 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)]，其满足下列能量泛函的最小值：

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

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

$\left(\begin{array}{c}\mathrm{\mu \Delta u}-{|▿f|}^2(u-f_x)=0\\ \mathrm{\mu \Delta 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{\Delta u}-h\left(\right|▿f\left|\right)(u-f_x)=0\\ g\left(\right|▿f\left|\right)\mathrm{\Delta v}-h\left(\right|▿f\left|\right)(v-f_y)=0\end{array}\right)---\left(5\right)$

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

$\mathrm{\Delta f}=\frac{{\partial }^{2}f}{{\partial x}^2}+\frac{{\partial }^{2}f}{{\partial 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{\Delta f}=f\otimes 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\otimes G_{24})-h\left(\right|▿\left|\right)(u-f_x)\\ v_t=g\left(\right|▿f\left|\right)(v\otimes 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\otimes {\mathrm{NS}}_{24}-u\otimes {\mathrm{RM}}_{24})-h\left(\right|▿f\left|\right)(u-f_x)\\ v_t=g\left(\right|▿f\left|\right)(v\otimes {\mathrm{NS}}_{24}-v\otimes {\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\otimes {\mathrm{NS}}_{24}-u\otimes {\mathrm{RM}}_{24})-h\left(\right|▿f\left|\right)(u-f_x)\\ v_t=g\left(\right|▿f\left|\right)(v\otimes {\mathrm{NS}}_{24}-v\otimes {\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}\times {\int }_0^1{((x\left(s\right)-x_c)\cos \left(\theta \right)+(y\left(s\right)-y_c)\sin \left(\theta \right)-r_1\cos (2\mathrm{\pi s}-\theta ))}^2\mathrm{ds}$

(14)

$+\frac{1}{2}\times {\int }_0^1{(-(x\left(s\right)-x_c)\sin \left(\theta \right)+(y\left(s\right)-y_c)\cos \left(\theta \right)-r_2\sin (2\mathrm{\pi s}-\theta ))}^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{\pi i}/n-\theta )\cos \left(\theta \right){+r}_2\sin (2\mathrm{\pi i}/n-\theta )\sin \left(\theta \right)=0\\ y_i-y_c-r_1\cos (2\mathrm{\pi i}/n-\theta )\sin \left(\theta \right)-r_2\sin (2\mathrm{\pi i}/n-\theta )\cos \left(\theta \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\otimes {\mathrm{NS}}_{24}-u\otimes {\mathrm{RM}}_{24})-h\left(\right|▿f\left|\right)(u-f_x)\\ v_t=g\left(\right|▿f\left|\right)(v\otimes {\mathrm{NS}}_{24}-v\otimes {\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}\times {\int }_0^1{((x\left(s\right)-x_c)\cos \left(\theta \right)+(y\left(s\right)-y_c)\sin \left(\theta \right)-r_1\cos (2\mathrm{\pi s}-\theta ))}^2\mathrm{ds}$

(19)

$+\frac{1}{2}\times {\int }_0^1{(-(x\left(s\right)-x_c)\sin \left(\theta \right)+(y\left(s\right)-y_c)\cos \left(\theta \right)-r_2\sin (2\mathrm{\pi s}-\theta ))}^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{\pi i}/n-\theta )\cos \left(\theta \right){+r}_2\sin (2\mathrm{\pi i}/n-\theta )\sin \left(\theta \right)=0\\ y_i-y_c-r_1\cos (2\mathrm{\pi i}/n-\theta )\sin \left(\theta \right)-r_2\sin (2\mathrm{\pi i}/n-\theta )\cos \left(\theta \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。

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

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

其中，S＝{s1，...，sn}，M＝{m1，...，mk}分别表示Snake轮廓和手工轮 廓上的点，对于整套图像而言，左室壁内膜的平均MAD值 为0.45像素，基本上与手工分割结果相同；外膜的平均MAD值为1.26像素，与手 工分割结果非常接近。