CN101224110A - 三维心肌形变应变计算方法 - Google Patents

Abstract

本发明公开了一种三维心肌形变应变计算方法。本发明在完成加标记的心脏核磁共振影像序列分割处理的基础上，将心肌按时间空间坐标划分多个局部(子)区间，使用前馈神经网络(BPNN)、多项式或支持向量机(SVM)拟合局部位移场，再通过牛顿迭代方法对建立的局部连续位移场进行计算，求解心肌质点的运动参数和使用非线性插补技术求解心肌应变参数。本发明物理意义明确，算法简单有效，模型适合于并行计算，能够计算任意心肌质点的前向与后向运动，应变计算结果可作为三维心脏超声应变测试结果的评判标准。

Description

三维心肌形变应变计算方法

技术领域

本发明属于加标记的心脏核磁共振成像(tagged MRI)的心肌形变和应变分析技术，特别是一种基于局部插值的3D(三维)心肌应变计算方法。

背景技术

心脏运动是非刚体的复杂3D运动，在周期性的搏动中发生形变(包括位移，剪切，旋转，收缩或者扩张等)。通过计算心肌质点3D的运动，可获得心肌各点的3D应变量。3D心肌应变数据是心脏疾病的临床诊断和心脏运动的研究的重参数。因此，分析3D心脏的应变的时空分布具有重要的研究价值，是目前国内外生物医学领域研究的热点之一。

当前，非介入性(noinvasive)的心脏医学成像主要方法有心脏超声成像和核磁共振成像(MRI)。加标记的心脏核磁共振成像技术(Tagged MRI)是目前心脏影像分析研究中采用最多的成像方式之一，它将磁饱和模式在射频脉冲的作用下加到心肌上形成标记面(tag面)，使心肌内的tag面随着心肌一起运动，通过观测标记面与影像平面交线(tag线)的运动研究心肌内部的复杂运动，如图1所示。

但tagged MRI是在2D成像平面上获取的离散的心肌质点位移场信息。不能直接得到心肌的3D位移。因此必须通过一些特殊方法，利用所获得的稀疏的2D位移信息来重建心肌的3D形变和计算心肌的应变。目前使用的方法有：(1)生物力学方法。将心肌假设为不可压缩的弹性固体物质，满足材料力学中的压力——应力关系，服从Hooke定理；再用有限元方法对心肌运动进行分析和计算。(P.Shi.Image Analysis of3D CardiacMotion Using Physical and Geometrical Models[D]Ph.D.Dissertation，Yale University，1996)(2)随机模型方法。通过构造随机过程或随机场来得到度量模型，再使用估计方法求解模型。(L.Yan，T.S.Denny.Unsupervised Estimation of left Ventricular displacementfrom MR tagged Images using Markov Random field edge Priors[C].IEEE，1998.)(3)可变形模型方法。应用几何方法对心肌形状进行描述，用力学模型来跟踪心脏的运动和形变，用逼近理论来解决模型求解问题。是当前使用较多的一种方法。(4)B样条模型方法。借助于B样条可以表示单方向的tag线位移和tag面的变形，描绘3D立体时空连续统。(A.Amini.R.W.Curwen，John C.Gore.Snakes and splines for tracking no-rigid heartmotion[C].An European Conference on Computer vision，University of Cambridge，UK，April 1996：251-261.)(5)光流方法。光流(optical flow)方法可以从相邻的图像帧中分析物体真实的2D运动情况。可用该方法对加标记线的MRI图像进行了分析。但此方法在MRI图像的处理应用上还有待研究上。

但上述方法都存在着算法复杂、耗时多、计算精度低、(心肌连续性和等容性)前提假设有争议等问题。基于局部拟合的心肌质点位移方法提出了使用前馈神经网络进行拟合心肌位移的方法，但没有解决计算心肌应变的问题(朱近王平安夏德深：基于BPNN方法的心肌形变计算，计算机研究与发展2005 Vol.42 No.12 pp：2143-2149)。

发明内容

本发明的目的在于提供一种由离散2D心肌运动场数据重建心肌质点3D连续运动模型及计算3D心肌形变、应变的方法。

实现本发明目的的技术解决方案为：一种三维心肌形变应变计算方法，在完成加标记的心脏核磁共振影像序列分割处理的基础上，首先进行数据后处理，对心肌按时间、空间坐标划分多个局部区间，拟合局部位移场，然后对位移场进行计算，再通过牛顿迭代法求解心肌质点的运动和心肌形变参数，最后使用非线性局部插补技术计算心肌应变参数。

本发明与现有技术相比，其显著优点：(1)物理意义明确。计算中仅使用一个心动周期的标记线交点信息建立心肌质点运动的计算模型，模型中没引入心肌是各向均匀的和不可压缩等假设，因而是可行的。(2)算法简单有效。不需要建立和求解复杂的偏微分方程或高维B样条曲面。(3)模型适合于并行计算。由于每个函数的拟合与计算仅与局部区间数据相关，在硬件环境允许时，可以很容易地实现并行算法，以提高计算速度。(4)能够计算任意心肌质点的前向与后向运动。(5)应变计算结果可作为三维心脏超声应变测试结果的评判标准。

下面结合附图对本发明作进一步详细描述。

附图说明

图1是现有技术的成像平面与标记面关系图。

图2是现有技术的PBNN结构图。

具体实施方式

本发明三维心肌形变应变计算方法，在完成加标记的心脏核磁共振影像序列分割处理的基础上，首先进行数据后处理，对心肌按时间、空间坐标划分多个局部区间，拟合局部位移场，然后对位移场进行计算，再通过牛顿迭代法求解心肌质点的运动和心肌形变参数，最后使用非线性局部插补技术计算心肌应变参数。该方法的操作过程如下：

(1)原始数据整理

原始数据来自对加标记的心脏核磁共振影像(Tagged MRI)序列的分割结果。短轴方向的tag线为网格型，间隔为ds；长轴方向的tag线为平行线，间隔为dl。n个采样时间点T＝{t₀，t₁，...t_n-1}，其中t₀为舒张末期(参考)时刻。短轴影像面间隔为dz；长轴影像面间隔为dθ。在对Tagged MRI序列进行了分割处理后获得了各层图像tag线的网格交点和tag线与心肌内外轮廓的交点，并对所有交点进行了编码。整理后的分割结果见表1和表2。令p_i表示在时间t₀，编码为i的交点(x(t₀)，y(t₀)，z(t₀)，表1中t₀列的m_s个点{p₁，p₂，...p_ms}构成了短轴网格节点集合Sp；同样表2中t₀列的m_l个点{p₁，p₂，...p_ml}构成长轴网格节点集合Lp。

表1整理后的短轴分割结果

交点编码(i)	t0(Sp)			t1		...	tn-1
									x(t0)	y(t0)	z(t0)	Δx(t1)	Δy(t1)	...	Δx(tn-1)	Δy(tn-1)
	1
2
	...
ms

其中2-D位移场分量计算：Δx(t_i)＝x(t_i)-x(t₀)，Δy(t_i)＝y(t_i)-y(t₀)。

表2整理后的长轴分割结果

交点编码(i)	t0(Lp)			t1	...	tn-1
							θ(t0)	r(t0)	z(t0)	Δz(t1)	...	Δz(tn-1)
	1
2
	...
ml

表2中θ(t₀)＝arctg(y(t₀)/x(t₀))， $r\left(t_0\right)=\sqrt{x^2\left(t_0\right)+y^2\left(t_0\right)},$ Δz(t_i)＝z(t_i)-z(t₀)。

已建立的t₀时刻描述心肌外壁和内壁的数学模型：

$\left\{\begin{array}{c}\mathrm{\θ}=\mathrm{arctg}(y\left(t_0\right)/x\left(t_0\right))\\ R=f_{\mathrm{out}}(\mathrm{\θ}\left(t_0\right),z\left(t_0\right))\\ r=f_{\mathrm{in}}(\mathrm{\θ}\left(t_0\right),z\left(t_0\right))\end{array}\right.---\left(I\right)$

式中R，r分别表示在t₀时刻，柱坐标中心肌外壁和内壁与z轴的距离。通过公式(I)可以判断一个空间点p(x(t₀)，y(t₀)，z(t₀))在心肌内的条件为：

$f_{\mathrm{in}}(\mathrm{\θ}\left(t_0\right),z\left(t_0\right))\≤\sqrt{x^2\left(t_0\right)+y^2\left(t_0\right)}\≤f_{\mathrm{out}}(\mathrm{\θ}\left(t_0\right),z\left(t_0\right))---\left(\mathrm{II}\right)$

(2)局部连续位移场方程组的建立

从表1和2中可得到采样时刻t_i位于离散网格节点上心肌质点的位移分量。由这些稀疏的网格节点建立描述心肌质点连续位移场的方程：

$\left\{\begin{array}{c}\mathrm{\Δx}\left(t\right)=f_x(x\left(t_0\right),y\left(t_0\right),z\left(t_0\right),t)\\ \mathrm{\Δy}\left(t\right)=f_y(x\left(t_0\right),y\left(t_0\right),z\left(t_0\right),t)\\ \mathrm{\Δz}\left(t\right)=f_z(x\left(t_0\right),y\left(t_0\right),z\left(t_0\right),t)\end{array}\right.---\left(\mathrm{III}\right)$

真实心肌质点的运动模式是未知的，而无法有效确定方程(III)的形式。为此，我们将心肌按时间空间坐标分别建立以(p_i，p_j)和(p_i’，t_j)为区域中心的n*ms和n*ml个局部(子)区间。对各子区间可以使用前馈神经网络(BPNN)、多项式或支持向量机(SVM)方式建立拟合函数。对S_p和L_p中节点i在时间t_j(j＝0，1，...n-1)上，实现对应区域上的拟合方法如下：

a.局部(子)区间的划分

取表1中i行的x_i0＝x(t₀)，y_i0＝y(t₀)，z_i0＝z(t₀)和t_j为(4D)区域中心(i＝1，2，...m_s；j＝0，1，...n-1)，寻找所有满足条件：|x_k0-x_i0|≤d∧|y_k0-y_i0|≤d∧|z_k0-z_i0|≤z_ds的k行，作为短轴子区间i；其中d＞0为短轴影像面上初始tag线的间距的1.5倍，z_ds＞0为短轴影像面间距的1.5倍。同理，对长轴取表2中的p_i’行的θ_i0＝θ(t₀)，r_i0＝r(t₀)，z_i0＝z(t₀)和t_j为区域中心，寻找所有满足条件：|θ_k0-θ_i0|≤θ_d∧ |z_k0-z_i0|≤z_dt的p_k’行，其中θ_d＞0为长轴影像面的间夹角的1.5倍，作为长子区间i’；z_dt＞0为长轴初始tag面的间距的1.5倍。

b.训练集和教师集的建立

在短轴子区间中将(x_k0，y_k0，z_k0，t_j-1)，(x_k0，y_k0，z_k0，t_j)和(x_k0，y_k0，z_k0，t_j+1)加入短轴训练输入矢量集PT_s中；取同一行中3个相邻时间点的2D位移(Δx_k(t_j-1)，Δy_k(t_j-1))，(Δx_k(t_j)，Δy_k(t_j))和(Δx_k(t_j+1)，Δy_k(t_j+1)加入短轴的教师矢量集S_s中。同样，将(θ_k0，r_k0，z_k0，t_j-1)，(θ_k0，r_k0，z_k0，t_j)和(θ_k0，y_k0，z_k0，t_j+1)加入长轴训练集PT_l中；取同一行中3个时间点的1D位移分量Δz_k(t_j-1)，Δz_k(t_j)和Δz_k(t_j+1)加入长轴教师矢量集S_l中。

c.局部拟合函数的建立

根据所选用的拟合函数形式，用对应训练输入集、教师集数据和均方误差最小原则计算各拟合函数f的具体参数(系数)，即建立了局部连续位移场方程组。

(3)位移场分量的计算

方程组(III)的建立解决了心肌位移场的计算问题。由此基础上可求解在t₀时刻位于左心室上的任意心肌质点p(x，y，z)，在t时刻位移分量(Δx(t)，Δy(t)和Δz(t))，计算的过程如下：

a.心肌质点的检测

对点p(x，y，z)，判断其在t₀时刻是否位于心肌内的计算为： $f_{\mathrm{in}}(\mathrm{\θ}\left(t_0\right),z\left(t_0\right))\≤\sqrt{x^2\left(t_0\right)+y^2\left(t_0\right)}\≤f_{\mathrm{out}}(\mathrm{\θ}\left(t_0\right),z\left(t_0\right)).$ 式中的f_in和f_out由公式(II)得到。

b.计算区域中心

令：t_j∈T且min{|t-t_j||j＝0，1，...n-1}，p_i∈S_p且min{|p-p_i||i＝1，2，...m_s}；其中|p-p_i|表示p与p_i两点的欧氏距离。此时(p_i，t_j)即为心肌质点p在t时刻，短轴方向的位移场计算的区域中心。同理设：p_k∈L_p且min{|p-p_k||k＝1，2，...m_l}。(p_k，t_j)为心肌质点p在t时刻，长轴方向的位移场计算的区域中心。

c.计算位移场分量

在获取的短轴与长轴位移场计算区域中心后，选择对应的局部拟合函数(f_x，f_y和f_z)，再由(III)式计算出p点的3个位移场分量Δx(t)，Δy(t)和Δz(t)。

(4)3D心肌质点位移计算

为计算点p的真实位移(ΔX，ΔY，ΔZ)，我们建立了牛顿迭代公式：

Δz⁰＝Δx⁰＝Δy⁰＝0

$\left\{\begin{array}{c}\mathrm{\Δ}{z}^{(i+1)}=f_z(x\left(t_0\right)+\mathrm{\Δ}{x}^{\left(i\right)},y\left(t_0\right)+\mathrm{\Δ}{y}^{\left(i\right)},z\left(t_0\right),t)\\ \mathrm{\Δ}{x}^{(i+1)}=f_x(x\left(t_0\right),y\left(t_0\right),z\left(t_0\right)+\mathrm{\Δ}{z}^{(i+1)},t)\\ \mathrm{\Δ}{y}^{(i+1)}=f_y(x\left(t_0\right),y\left(t_0\right),z\left(t_0\right)+\mathrm{\Δ}{z}^{(i+1)},t)\end{array}\right.---\left(\mathrm{IV}\right)$

式中时间t的取值应在心动周期内，(x(t₀)，y(t₀)，z(t₀))是心肌质点p在参考时刻的坐标。迭代结束条件为：|Δx⁽ⁱ⁺¹⁾-Δx⁽ⁱ⁾|＜ε，|Δy⁽ⁱ⁺¹⁾-Δy⁽ⁱ⁾|＜ε和|Δz⁽ⁱ⁺¹⁾-Δz⁽ⁱ⁾|＜ε，在实际计算时取ε＝0.01mm，一般在2-5次迭代后即可收敛。取ΔX＝Δx⁽ⁱ⁺¹⁾，ΔY＝Δy⁽ⁱ⁺¹⁾和ΔZ＝Δz⁽ⁱ⁺¹⁾。通过对心肌质点运动的逐点计算，即可完成对整个心脏的形变的计算。(IV)式可以简单表示为：

(ΔX，ΔY，Δz)＝f_xyz(x(t₀)，y(t₀)，z(t₀)，t)    (V)

心肌质点p在t时刻的坐标p(t)＝p(t₀)+(ΔX(t)，ΔY(t)，ΔZ(t))。

(5)基于非线性局部插补的心肌应变计算方法

应变是弹性力学中的概念，能够包括形变的全部信息。因此被广泛应用于描述弹性塑性体的形变。针对心脏应变分析的方法主要采用基于Green应变方法。

设在t时刻心肌质点p(t)具有坐标(y₁，y₂，y₃)，如经过常数时间Δt，p(t+Δt)点移动到(x₁，x₂，x₃)位置。数学关系为：x_i＝x_i(y₁，y₂，y₃)i＝1，2，3，y_i＝y_i(x₁，x₂，x₃)i＝1，2，3，心肌质点p的位移为：x_i＝y_i+u_i  i＝1，2，3，u_i称为心肌质点p的相对位移。p点的Green应变计算公式：

$e_{\mathrm{ij}}=\frac{1}{2}(\∂x_i/\∂y_j+\frac{\∂x_i}{\∂y_i}+\underset{k=1}{\overset{3}{\mathrm{\Σ}}}\∂x_k/\∂y_i\·\∂x_k/\∂y_j),i,j=\mathrm{1,2},3---\left(\mathrm{VI}\right)$

有6个独立应变分量{e₁₁，e₁₂，e₁₃，e₂₂，e₂₃，e₃₃}，其中e₂₁＝e₁₂，e₃₂＝e₂₃。

a.心肌质点的初始坐标计算

首先需要计算心肌质点p(t)在参考时刻(舒张末期)T＝0的初始坐标值。即由p(t)反向计算p(0)。我们设计了迭代方法计算：

p⁽⁰⁾∈LV

$\left\{\begin{array}{c}(\mathrm{\Δ}{x}^{\left(i\right)},\mathrm{\Δ}{y}^{\left(i\right)},\mathrm{\Δ}{z}^{\left(i\right)})=f_{\mathrm{xyz}}(x{\left(0\right)}^{\left(i\right)},y{\left(0\right)}^{\left(i\right)},z{\left(0\right)}^{\left(i\right)},t)\\ x{\left(0\right)}^{(i+1)}=x\left(t\right)-k_x\left(i\right)\mathrm{\Δ}{x}^{\left(i\right)}\\ y{\left(0\right)}^{(i+1)}=y\left(t\right)-k_y\left(i\right)\mathrm{\Δ}{y}^{\left(i\right)}\\ z{\left(0\right)}^{(i+1)}=z\left(t\right)-k_z\left(i\right)\mathrm{\Δ}{z}^{\left(i\right)}\end{array}\right.---\left(\mathrm{VII}\right)$

式中p₀∈LV表示初始点p₀位于心肌内部。p(t)＝(x(t)，y(t)，z(t)＝(y₁，y₂，y₃)为已知量，迭代结束条件为：|x(0)⁽ⁱ⁺¹⁾-x(0)⁽ⁱ⁾|＜ε，|y(0)⁽ⁱ⁺¹⁾-y(0)⁽ⁱ⁾|＜ε和|z(0)⁽ⁱ⁺¹⁾-z(0)⁽ⁱ⁾|＜ε。ε是预定一个小正数，一般取ε＝0.1mm。f_xyz是计算(正向)心肌质点位移的公式(V)。k_x，k_y，k_z是修正因子，与迭代次数i有关，取k_x＝k_y＝k_z＝1/(i+1)。结果为：p(0)＝(x(0)⁽ⁱ⁺¹⁾，y(0)⁽ⁱ⁺¹⁾，z(0)⁽ⁱ⁺¹⁾)。(VII)式可以简单表示为：

$(x\left(0\right),y\left(0\right),z\left(0\right))=f_{\mathrm{xyz}}^{-1}(x\left(t\right),y\left(t\right),z\left(t\right),t)---\left(\mathrm{VIII}\right)$

b.心肌质点的相对位移计算

已知心肌质点p(t)＝(y₁，y₂，y₃)，计算该点经过Δt时间后p(t+Δt)＝(x₁，x₂，x₃)的相对位移(u₁，u₂，u₃)方法：

$\left\{\begin{array}{c}(x\left(0\right),y\left(0\right),z\left(0\right))=f_{\mathrm{xyz}}^{-1}(y_1,y_2,y_3,t)\\ (x_1,x_2,x_3)=f_{\mathrm{xyz}}(x\left(0\right),y\left(0\right),z\left(0\right),t+\mathrm{\Δt})\\ (u_1,u_2,u_3)=(x_1,x_2,x_3)-f_{\mathrm{xyz}}(x\left(0\right),y\left(0\right),(z\left(0\right),t))\end{array}\right.---\left(\mathrm{IX}\right)$

c.网格节点的心肌应变计算

设P(t)＝(y₁，y₂，y₃)是在t时刻位于网格节点上的心肌质点，选取同一时刻相邻影像层上相邻tag面的网格节点P’(t)＝(y₁’，y₂’，y₃’)。在t+Δt时刻P点的位置为P(t+Δt)＝(x₁，x₂，x₃)＝(y₁+u₁，y₂+u₂，y₃+u₃)，P’(t+Δt)点的坐标为(x₁’+u₁’，x₂’+u₂’，x₃’+u₃’)，可由公式(IX)计算得到。心肌变形后，两质点的位移差为：Δu＝(Δu₁，Δu₂，Δu₃)＝(u₁-u₁，u₂-u₂，u₃-u₃)。

在x方向(j＝1)，Δx＝y₁’-y₁，其x方向(i＝1)位移的偏微分为x_i/y_j＝x₁/y₁≈Δu₁/Δx。可用同方法分别计算y(i，j＝2)和z(i，j＝3)方向的情况。再通过公式(VI)可计算出各网格节点P的6个应变分量{E_ij}(i，j＝1，2，3)。

d.用局部插值方法计算非网格的心肌质点的应变

为了减少计算量，对其它非网格节点p的应变分量{e_ij}(i，j＝1，2，3)的计算采用通过对已求得的相邻网格节点{P}的应变分量{E}插值获得：

首先将感兴趣心肌区域的网格节点的应变{E_ij}(i，j＝1，2，3)全部计算出来。设非网格节点的心肌质点p坐标为(x，y，z)，建立对应的网格节点子集S_p＝{P(y₁，y₂，y₃)||x-y₁|＜d∧|y-y₂|＜d∧|z-y₃|＜d’，P∈网格节点}，d取tag线间隔的1.5倍，d’取短轴影像面间隔的1.5倍，S_p中节点数量为n_p。设网格节点P₁的坐标为(x₁，y₁，z₁)，应变为{E_ij}₁(1＝1，2，...n_p，i，j＝1，2，3)，心肌质点p的应变{e_ij}_p用非线性插补方法计算：

$\left\{\begin{array}{c}{d}_l^2={(x_l-x)}^2+{(y_l-y)}^2+\mathrm{\ϵ}\\ w_l=d_l^2/\underset{i=1}{\overset{n_k}{\mathrm{\Σ}}}{d}_i^2,l=\mathrm{1,2},\mathrm{\Λ},n_p\\ {\left\{e_{\mathrm{ij}}\right\}}_p=\underset{l=1}{\overset{n_l}{\mathrm{\Σ}}}{\left\{E_{\mathrm{ij}}\right\}}_l\·w_l,i,j=\mathrm{1,2,3}\end{array}\right.---\left(X\right)$

式中的ε用于平滑噪声影响避免出现锯齿状分布，取ε＝0.01。用公式(X)对感兴趣的心肌质点进行计算，即可实现所需的应变计算。

实施例：本发明的三维心肌形变应变计算方法过程为：

1、原始数据整理

采样时刻为T＝{11，15，...49，...60}，共有n＝15个时间点，其中t₀＝11为心脏舒张期末；x和y方向的tag面数目各取7，z方向的tag面数为9；短轴影像面数取10(ds＝0.6132cm)，长轴影像面数取8(dθ＝22.5°＝0.3927弧度)。

表1短轴分割结果

编码(i)	t0(Sp)			t1		...	tn-1
										x(11)	y(11)	z(11)	Δx(15)	Δy(15)	...	Δx(60)	Δy(60)
1	-1.7345	-1.3297	-1.5107	0.2811	0.0116		-0.0902	-0.0094
									2	-1.7345	-1.1563	-1.5107	0.2590	0.0322		-0.085	-0.0076
...
									ms	0.9954	0.5782	4.0073	-0.1465	0.4410		0.9887	-0.0067

表2长轴分割结果

编码(i)	t0(Lp)			t1	...	tn-1
							θ(11)	r(11)	z(11)	Δz(15)	...	Δz(60)
	1	-3.1416	2.1959	-1.4493	-0.0378								0.0903
2							-3.1416	1.3430	-1.4493	0.0017		0.0637
	...
Ml							-0.3927	0.5287	3.9460	0.1675		-0.0472

2、建立以(p₁，t₀)为区域中心的局部拟合函数

采用前馈神经网络(BPNN)模型，区域中心(-1.7345，-1.3297，-1.5107，15)，取阈值：d＝0.5782*1.5＝0.8673；

建立短轴训练集：PT_s＝{(-1.7345，-1.3297，-1.5107，11)，(-1.7345，-1.3297，-1.5107，15)，(-1.7345，-1.1563，-1.5107，11)，...}；教师集：S_s＝{(0，0)，(0.2811，0.0116)，(0，0)，...}；

建立长轴训练集：PT₁＝{(-3.1416，2.1959，-1.4493，11)，(-3.1416，2.1959，-1.4493，15)，(-3.1416，1.3430，-1.4493，11)，...}；教师集：S₁＝{0，-0.0378，0，...}。

采用2层结构BP神经网络进行拟合，见附图。

所有BPNN的输入节点为4(x，y，z，t)长轴输出层节点为1(Δz)，短轴输出层节点为2(Δx，Δy)。因心肌几何形状的不规则性，对不同的区域中心所建立的训练集(教师集)的大小也不一样。PT_s和S_s中元素个数为18～54，PT_l和S_l的中元素为16～36。中间层采用双曲正切S型(tansig)传递函数，输出层为线性(purelin)传递函数可以满足位移场的拟合要求。短轴PBNN中间层神经元数目在2～8之间，长轴PBNN中间层神经元数目在2～5之间随训练集的大小而变化。用对BPNN进行训练，即建立了局部拟合函数(III)。

3、位移场分量的计算

取t₀＝11时刻位于心肌内的心肌质点p(x，y，z)＝(-1.70，-1.30，-1.50)，时间t＝16，得到计算区域中心(-1.7345，-1.3297，-1.5107，15)，选择对应的局部BPNN(f_xy和f_z)，再由(IV)式计算出p点的相对于舒张末期t₀＝11时刻3个位移场分量Δx(16)，Δy(16)和Δz(16)。

4、3D心肌质点位移计算

为计算心肌质点p＝(-1.70，-1.30，-1.50)在时间t＝16的真实位移，将坐标值和时间值代入公式(V)：即(Δx，Δy，Δz)＝f_xyz(-1.70，-1.30，-1.50，16)可计算出p点的3D真实位移(Δx，Δy，Δz)。

在完成加标记的心脏核磁共振影像序列分割处理的基础上，将心肌按时间空间坐标划分多个局部(子)区间，使用前馈神经网络(BPNN)、多项式或支持向量机(SVM)拟合局部位移场，再通过牛顿迭代方法对建立的局部连续位移场进行计算，求解心肌质点的运动和心肌形变参数。

Claims (6)

1.一种三维心肌形变应变计算方法，其特征在于：在完成加标记的心脏核磁共振影像序列分割处理的基础上，首先进行数据后处理，对心肌按时间、空间坐标划分多个局部区间，拟合局部位移场，然后对位移场进行计算，再通过牛顿迭代法求解心肌质点的运动和心肌形变参数，最后使用非线性局部插补技术计算心肌应变参数。

2.根据权利要求1所述的三维心肌形变应变计算方法，其特征在于：含有短轴方向和长轴方向的影像序列作分割数据后处理，即来自对加标记的心脏核磁共振影像Tagged MRI序列的分割结果，其短轴方向的tag线为网格型，间隔为ds；其长轴方向的tag线为平行线，间隔为dl；影像序列的n个采样时间点T＝{t₀，t₁，...t_n-1)，其中t₀为舒张末期作为参考时刻，短轴影像面间隔为dz；长轴影像面间隔为dθ，在对Tagged MRI序列进行分割获得各层图像tag线的网格交点和tag线与心肌内外轮廓的交点，经分割处理后得到所有交点的编码i、交点参考时刻坐标p_i和t_j时刻的位移分量Δx(t_j)及Δy(t_j)的短轴分割处理结果，其中编码i＝1，2，...m_s；p_i为x_i0＝x(t₀)，y_i0＝y(t₀)，z_i0＝z(t₀)；Δx(t_j)＝x(t_j)-x(t₀)，Δy(t_j)＝y(t_j)-y(t₀)，j＝1，2...n-1；以及包括所有交点的编码i、交点参考时刻极坐标q_i和t_j时刻的位移Δz(t_j)的长轴分割结果，其中编码i＝1，2，...m_l·q_i为θ_i0＝θ(t₀)，r_i0＝r(t₀)，z_i0＝z(t₀)；θ(t₀)＝arctg(y(t₀)/x(t₀))， $r\left(t_0\right)=\sqrt{x^2\left(t_0\right)+y^2\left(t_0\right)},$ Δz(t_j)＝z(t_j)-z(t₀)，j＝1，2...n-1。

3.根据权利要求1所述的三维心肌形变应变计算方法，其特征在于：首先，进行局部区间的划分，即：

(1)将心肌按时间空间坐标分别建立以(P_i，t_j)和(q_i，t_j)为区域中心的n*m_s和n*m_l个局部区间，p_i表示在时间t₀，编码为i的短轴交点(x(t₀)，y(t₀)，z(t₀))，q_i表示在时间t₀，编码为i的长轴交点(θ(t₀)，r(t₀)，z(t₀))；

(2)以(p_i，t_j)和(q_i，t_j)为区域中心来确定局部区间，取短轴分割处理结果编码为i的交点p_i为x_i0＝x(t₀)，y_i0＝y(t₀)，z_i0＝z(t₀)和t_j为4D区域中心(i＝1，2，...m_s；j＝0，1，...n-1)，寻找所有满足条件为|x_k0-x_i0|≤d∧|y_k0-y_i0|≤d∧|z_k0-z_i0|≤z_ds的k行，作为短轴局部区间I；其中d＞0为短轴影像面上初始tag线的间距的1.1～1.9倍，z_ds＞0为短轴影像面间距的1.1～1.9倍；同理，对长轴分割处理结果编码为i的交点q_i的θ_i0＝θ(t₀)，r_i0＝r(t₀)，z_i0＝z(t₀)和t_j为4D区域中心，寻找所有满足条件：|θ_k0-θ_i0|≤θ_d∧|z_k0-z_i0|≤z_dt的q_k交点，其中θ_d＞0为长轴影像面的间夹角的1.1～1.9倍，作为长轴局部区间I’；z_dt＞0为长轴初始tag面的间距的1.1～1.9倍；

然后，建立训练集和教师集，即：

在短轴局部区间I中将(x_k0，y_k0，z_k0，t_j-1)，(x_k0，y_k0，z_k0，t_j)和(x_k0，y_k0，z_k0，t_j+1)加入短轴训练矢量集PT_s中，取3个相邻时间点的2D位移(Δx_k(t_j-1)，Δy_k(t_j-1))，(Δx_k(t_j)，Δy_k(t_j)和(Δx_k(t_j+1)，Δy_k(t_j+1))加入短轴的教师矢量集S_s中；同样，在长轴局部区间I’中将(θ_k0，r_k0，z_k0，t_j-1)，(θ_k0，r_k0，z_k0，t_j)和(θ_k0，y_k0，z_k0，t_j+1)加入长轴训练集PT_l中，取3个相邻时间点的1D位移分量Δz_k(t_j-1)，Δz_k(t_j)和Δz_k(t_j+1)加入长轴教师矢量集S_l中；

最后，局部拟合函数的建立，即：

选用前馈神经网络BPNN、多项式或支持向量机SVM的拟合函数形式，用对应训练输入集、教师集数据和均方误差最小原则计算各拟合函数f的具体参数，即建立局部连续位移场方程组，

$\left\{\begin{array}{c}\mathrm{\Δx}\left(t\right)=f_x(x\left(t_0\right),y\left(t_0\right),z\left(t_0\right),t)\\ \mathrm{\Δy}\left(t\right)=f_y(x\left(t_0\right),y\left(t_0\right),z\left(t_0\right),t)\\ \mathrm{\Δz}\left(t\right)=f_z(x\left(t_0\right),y\left(t_0\right),z\left(t_0\right),t)\end{array}\right.$

x(t₀)，y(t₀)，z(t₀)是心肌质p在参考时刻t₀的坐标点；Δx(t_j)，Δy(t_j)和Δz(t_j)是p在时刻t的位移量；f_x，f_y，和f_z是拟合函数，具体形式依所选用的前馈神经网络BPNN、多项式或支持向量机SVM而定。

4.根据权利要求1所述的三维心肌形变应变计算方法，其特征在于：位移场的计算，在局部连续位移场方程组的基础上求解在t₀时刻位于左心室上的任意心肌质点p(x，y，z)，计算在t时刻位移分量(ΔX(t)，ΔY(t)和ΔZ(t)，其过程如下：

(1)心肌质点的检测，对质点p(x，y，z)，判断其在t₀时刻是否位于心肌内，即p(x，y，z)是否在心肌内壁f_in和心肌外壁f_out之间，判断式为

$f_{\mathrm{in}}(\mathrm{\θ}\left(t_0\right),z\left(t_0\right))\≤\sqrt{x^2\left(t_0\right)+y^2\left(t_0\right)}\≤f_{\mathrm{out}}(\mathrm{\θ}\left(t_0\right),z\left(t_0\right));$

(2)选择区域中心，令t_j∈T且min{|t-t_j||j＝0，1，...n-1}，p_i∈S_p且min{|p-p_i||i＝1，2，...m_s}；其中|p-p_i|表示p与p_i两点的欧氏距离，此时(p_i，t_j)即为心肌质点p在t时刻，短轴方向的位移场计算的区域中心；同理设p_k∈L_p且min{p-p_k||k＝1，2，...m_l}，(p_k，t_j)为心肌质点p在t时刻，长轴方向的位移场计算的区域中心；

(3)计算位移场，在获取的短轴与长轴位移场区域中心后，选择对应的局部拟合函数(f_x，f_y和f_z)，再由局部连续位移场方程组计算出p点位移场的3个分量Δx(t)，Δy(t)和Δz(t)。

5.根据权利要求1所述的三维心肌形变应变计算方法，其特征在于：3D心肌质点位移计算，即为计算点p的位移分量ΔX(t)，Δ(t)和ΔZ(t)，建立牛顿迭代公式：

Δz⁰＝Δx⁰＝Δy⁰＝0

$\left\{\begin{array}{c}{\mathrm{\Δz}}^{(i+1)}=f_z(x\left(t_0\right)+\mathrm{\Δ}{x}^{\left(i\right)},y\left(t_0\right)+\mathrm{\Δ}{y}^{\left(i\right)},z\left(t_0\right),t)\\ \mathrm{\Δ}{x}^{(i+1)}=f_x(x\left(t_0\right),y\left(t_0\right),z\left(t_0\right)+\mathrm{\Δ}{z}^{(i+1)},t)\\ \mathrm{\Δ}{y}^{(i+1)}=f_y(x\left(t_0\right),y\left(t_0\right),z\left(t_0\right)+\mathrm{\Δ}{z}^{(i+1)},t)\end{array}\right.$

式中时间t的取值应在心动周期内，(x(t₀)，y(t₀)，z(t₀))是心肌质点p在参考时刻的坐标，迭代结束条件为：|Δx⁽ⁱ⁺¹⁾-Δx⁽ⁱ⁾|＜ε，|Δy⁽ⁱ⁺¹⁾-Δy⁽ⁱ⁾|＜ε和|Δz⁽ⁱ⁺¹⁾-Δz⁽ⁱ⁾|＜ε，在实际计算时取ε＝0.05～0.1mm，取质点位移分量为ΔX＝Δx⁽ⁱ⁺¹⁾，ΔY＝Δy⁽ⁱ⁺¹⁾和ΔZ＝Δz⁽ⁱ⁺¹⁾，通过对心肌质点运动的逐点计算，即完成对整个心脏的形变的计算。上式简化表示为：(ΔX，ΔY，ΔZ)＝f_xyz(x(t₀)，y(t₀)，z(t₀)，t)

6.根据权利要求1所述的三维心肌形变应变计算方法，其特征在于应变计算的步骤为：

(1)心肌质点的初始坐标计算，即在已知T＝t时刻心肌质点p(t)p(t)＝(x(t)，y(t)，z(t))＝(y₁，y₂，y₃)，计算该点在参考时刻(舒张末期)T＝0的初始坐标值：p⁽⁰⁾∈LV

$\left\{\begin{array}{c}(\mathrm{\Δ}{x}^{\left(i\right)},\mathrm{\Δ}{y}^{\left(i\right)},\mathrm{\Δ}{z}^{\left(i\right)})=f_{\mathrm{xyz}}(x{\left(0\right)}^{\left(i\right)},y{\left(0\right)}^{\left(i\right)},z{\left(0\right)}^{\left(i\right)},t)\\ x{\left(0\right)}^{(i+1)}=x\left(t\right)-k_x\left(i\right)\mathrm{\Δ}{x}^{\left(i\right)}\\ y{\left(0\right)}^{(i+1)}=y\left(t\right)-k_y\left(i\right)\mathrm{\Δ}{y}^{\left(i\right)}\\ z{{\left(0\right)}^{(i+1)}=z\left(t\right)-k_z\left(i\right)\mathrm{\Δ}{z}^{\left(i\right)}}^{}\end{array}\right.$

p₀∈LV表示初始点p₀位于心肌内部，迭代结束条件为：|x(0)⁽ⁱ⁺¹⁾-x(0)⁽ⁱ⁾|＜ε，|y(0)⁽ⁱ⁺¹⁾-y(0)⁽ⁱ⁾|＜ε和|z(0)⁽ⁱ⁺¹⁾-z(0)⁽ⁱ⁾|＜ε，ε取值范围：ε＝0.1～0.5mm；f_xyz是计算(正向)心肌质点位移的公式，k_x，k_y，k_z是修正因子：k_x＝k_y＝k_z＝1/(i+1)；结果为：p(0)＝(x(0)⁽ⁱ⁺¹⁾，y(0)⁽ⁱ⁺¹⁾，z(0)⁽ⁱ⁺¹⁾)；以上计算简单表示为：

$(x\left(0\right),y\left(0\right),z\left(0\right))=f_{\mathrm{xyz}}^{-1}(x\left(t\right),y\left(t\right),z\left(t\right),t);$

(2)心肌质点的相对位移计算，即在已知心肌质点p(t)＝(y₁，y₂，y₃)，计算该点经过Δt时间后p(t+Δt)＝(x₁，x₂，x₃)的相对位移(u₁，u₂，u₃)方法：

$\left\{\begin{array}{c}(x\left(0\right),y\left(0\right),z\left(0\right))=f_{\mathrm{xyz}}^{-1}(y_1,y_2,y_3,t)\\ (x_1,x_2,x_3)=f_{\mathrm{xyz}}(x\left(0\right),y\left(0\right),z\left(0\right),t+\mathrm{\Δt})\\ (u_1,u_2,u_3)=(x_1,x_2,x_3)-f_{\mathrm{xyz}}(x\left(0\right),y\left(0\right),z\left(0\right),t)\end{array}\right.;$

(3)网格节点的心肌应变计算，即设P(t)＝(y₁，y₂，y₃)是在t时刻位于网格节点上的心肌质点，选取同一时刻相邻影像层上相邻tag面的网格节点P’(t)＝(y₁’，y₂’，y₃’)；在t+Δt时刻P点的位置为P(t+Δt)＝(x₁，x₂，x₃)＝(y₁+u₁，y₂+u₂，y₃+u₃)，P’(t+Δt)点的坐标为(x₁’+u₁’，x₂’+u₂’，x₃’+u₃’)；心肌变形后，两质点的位移差为：Δu＝(Δu₁，Δu₂，Δu₃)＝(u₁’-u₁，u₂’-u₂，u₃’-u₃)；

在x方向(j＝1)，Δx＝y₁’-y₁，其x方向(i＝1)位移的偏微分为x_i/y_j＝x₁/y₁≈Δu₁/Δx；用同步骤(3)的方法分别计算y(i，j＝2)和z(i，j＝3)方向的情况，再通过公式：

$E_{\mathrm{ij}}=\frac{1}{2}(\∂x_i/\∂y_i+\frac{\∂x_i}{\∂y_i}+\underset{k=1}{\overset{3}{\mathrm{\Σ}}}\∂x_k/\∂y_i\·\∂x_k/\∂y_i),i,j=\mathrm{1,2,3}$

计算出各网格节点P的6个应变分量{E_ij}(i，j＝1，2，3)；

(4)用局部插值方法计算非网格的心肌质点的应变，即对其它非网格节点p的应变分量{e_ij}(i，j＝1，2，3)的计算采用通过对已求得的相邻网格节点{P}的应变分量{E}插值获得；

设非网格节点的心肌质点p坐标为(x，y，z)，建立对应的网格节点子集S_p＝{P(y₁，y₂，y₃)||x-y₁|＜d∧|y-y₂|＜d∧|z-y₃|＜d’，P∈网格节点}，d取tag线间隔的1.1～2倍，d’取短轴影像面间隔的1.1～2倍，S_p中节点数量为n_p；设网格节点P₁的坐标为(x₁，y₁，z₁)，应变为{E_ij}₁(1＝1，2，...n_p，i，j＝1，2，3)，心肌质点p的应变{e_ij}_p用非线性插补方法计算：

$\left\{\begin{array}{c}{d}_l^2={(x_l-x)}^2+{(y_l-y)}^2+\mathrm{\ϵ}\\ w_l=d_l^2/\underset{i=1}{\overset{n_k}{\mathrm{\Σ}}}{d}_i^{2,}l=\mathrm{1,2},\mathrm{\Λ},n_p\\ {\left\{e_{\mathrm{ij}}\right\}}_p=\underset{l=1}{\overset{n_k}{\mathrm{\Σ}}}{\left\{E_{\mathrm{ij}}\right\}}_l\·w_l,i,j=\mathrm{1,2,3}\end{array}\right.$

式中的ε用于平滑噪声影响避免出现锯齿状分布，取ε＝0.001～0.1。