CN106446572A - 基于边界元模型和局部区域修正的肺部呼吸运动获取方法 - Google Patents

Abstract

本发明提供了一种基于边界元模型和局部区域修正的肺部呼吸运动获取方法，主要包括以下步骤：一、事先获取病人的胸部CT图像；二、在序列图像上对肺部区域进行图像分割；三、将不同呼吸时刻的胸部序列图像配准；四、获取肺组织呼吸运动的边界条件；五、结合边界元理论计算肺部表面边界元模型中的各参数并存储；六、在不同呼吸时刻的序列图像中选取约三十组特征点；七、利用线性插值方法对肺内部组织运动进行运动估计；八、根据特征点的运动信息更新肺内部区域运动信息。与现有技术相比，本发明能够体现呼吸运动的非均一性以及各向异性，具有计算速度快、实现方便、应用灵活等优点。

Description

基于边界元模型和局部区域修正的肺部呼吸运动获取方法

技术领域

本发明涉及肺部呼吸运动信息获取领域，特别涉及一种基于边界元模型和局部区域修正的肺部呼吸运动信息获取方法。

背景技术

精确治疗是当今医学的重要发展方向，是高速发展的计算技术、影像学技术与医学交叉发展的产物，可以为病人提供有效并且附加伤害最小的治疗，最大程度地降低风险，改善病人预后。肺癌的诊断率也随之有了一定的提高。目前，大多数的穿刺手术是基于传统三维CT图像引导下进行的，也是诊断及鉴别肺癌的一种新兴的有效方法，其可以实现微创下直接获取病变标本，并且将现代医学影像学技术与病理学相结合，做出组织病理学诊断，具有诊断率高，创伤小等特点。但由于人体解剖结构的动态特性，例如呼吸运动、胃肠道蠕动的变化等会导致病人器官和胸部、腹部等部位的肿瘤位置在治疗过程中发生运动。而三维CT图像是在呼吸周期中瞬时扫描所得，是静态的，基本不包含运动信息。为此穿刺活检往往需要重复多次扫描引导穿刺针进入并观察有无并发症，辐射剂量大，一定程度上限制了应用。因此，建立一种安全，高效，精确可靠并且实时的肺部肿瘤运动估计技术具有重要意义。

经过对现有文献的检索发现，现有的肺部呼吸运动信息获取方法大体可分为两大类：基于图像灰度的数学模型的运动估计和基于生物力学模型的运动估计。前者通过使两个或多幅图像之间进行配准使得图像间达到最大的相似度，从而得到一个肺部运动的运动场。其中，光流场方法在配准领域有着广泛的应用，假定序列图像中相邻图像间的时间间隔很短，同时相邻图像间的差异也很小，通过构建一个包含时间和空间梯度的信息函数达到配准的目的。这类方法需要对患者多次进行扫描CT。不仅对患者产生巨大的辐射，而且得到的运动场没有考虑到生理学的信息。造成了配准结果缺乏物理信息也有可能不满足弹性力学规律。基于生物力学模型的运动估计方法大多基于有限元的理论，虽然这种方法较好的解决了软组织各向异性及非均一性的问题，并且能够较好的描述肺部形变信息。但总体来说基于有限元的方法运算速度较慢，不能满足临床对于实时性的要求。

发明内容

本发明要解决的关键技术问题有：

1、利用不同时刻CT图像信息，结合配准方法得到的边界条件，建立边界元模型的肺部运动模型，此运动模型要具有人与人之间的各异性；

2、针对边界元模型只是对肺部表面进行仿真的不足，本发明将结合肺内特征点的运动信息，提供基于边界元的内部运动信息估计和区域局部修正方法；

3、对于临床对于运动估计实时性的要求，本发明将通过预计算和预存储技术，使得边界元模型在临床仿真过程中只需计算个别几个特征点的局部区域误差修正，从而实现肺部运动的实时仿真。

为了解决上述技术问题，本发明中披露了一种基于边界元模型和局部区域修正的肺部呼吸运动获取方法，本发明的技术方案是这样实施的：

基于边界元模型和局部区域修正的肺部呼吸运动获取方法，包括依次执行的以下步骤：

S1：获取同一人体的两组不同呼吸时刻的胸部CT图像；

S2：处理从步骤S1中获取的两组胸部CT图像，对其进行肺部组织的图像分割，得出肺部表面模型以及肺部血管树模型；

S3：将所述两组胸部CT图像进行配准，处理后得出肺部组织在不同呼吸时刻对应的序列图像；

S4：根据配准后的序列图像信息，提取出肺部组织呼吸运动的边界条件，从而获取肺部表面的运动信息；

S5：利用边界元理论和边界元条件迭代求解肺部表面模型中的各参数并存储，从而获取肺部表面随呼吸运动的边界元运动模型，简称肺部表面模型；

S6：在步骤S2中提取出的肺部血管树模型上设定10个以上的特征点，以这些特征点作为肺内部组织运动局域修正的依据；

S7：肺内部组织运动信息获取：以肺内部组织任一点作为目标点，由该点向肺部表面作六个方向的射线，获得射线与肺部表面模型的六个交点的运动信息，结合这六个交点的运动信息，通过线性插值法算出目标点的运动信息；

S8：目标点运动信息修正：确定肺部表面模型上与目标点距离最近的一个点，称该点为表面点；确定与目标点距离最近的两个特征点；以表面点信息以及该两个特征点信息对目标点的运动信息进行修正，最终得到整个肺部的运动模型。

优选地，步骤S6中，所述特征点位于血管交叉处。

优选地，步骤S3中配准的公式为：

式中，表示为每一个像素所在位置坐标的下舍入数值；同理得到系数 为控制点n_i×n_j×n_k的位移量，B样条的基函数表示为：

其中，0≤t＜1。

优选地，步骤S5中边界元三维弹性体的位移与应力基本解可表示为：

其中，υ为泊松系数，σ_ij为应变，G＝E/2(1+υ)，r代表P、Q两点间的距离，n表示表面的向量，δ为Dirac函数，以Q′为圆心半径为ε＝ε→0的球面。

优选地，步骤S5中边界元的离散边界积分方程可表示为：

其中，

边界Γ分为n各子边界Γ_j。

优选地，步骤S7中线性插值法公式为：

其中，目标点设为X，其坐标为(x，y，z)，P_xl，P_xr，P_ya，P_yp，P_zs，P_zi分别为该点向六个方向所作射线与肺部表面模型的交点，d(R_x)、d(R_y)、d(R_z)为X随呼吸运动的位移量，d(P_xl)、d(P_xr)、d(P_ya)、d(P_yp)、d(P_zs)、d(P_zi)为六个交点随呼吸运动的位移量。

优选地，步骤S8中目标点运动信息修正公式表示为：

其中，d(R)为线性插值的结果，由d(R_x)、d(R_y)、d(R_z)可得出；E_s为与X最近的表面点线性插值结果与边界元计算得到结果的数值差；E_a和E_b为两个特征点实际位移距离与插值得到的结果直接的数值差；D_s为X与最近表面点的欧氏距离；D_a和D_b为X与两个特征点的欧氏距离。

实施本发明的有益效果是：

1、采用配准不同时刻CT图像的方法获取模型的边界条件，从而得到了patient-specific的边界元模型，该运动模型解决了肺部运动的各向异性及非均一性问题，并且能够较好的描述肺部形变信息；

2、边界元模型有效地减少了计算量，提高了运动仿真的实时性；

3、针对边界元模型不能对内部组织的运动进行表达的缺点，本发明基于特征点的局部修正方法实现了对肺内部组织运动的精确计算方法，并且满足临床实时性的要求。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一种实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为一个实施例中，本发明的流程示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

在一个具体实施中，采用CPU为Xeon E5405 2.0GHZ，显卡为NVIDIA Quadro VX200，内存为4.0GB的计算机配置，实现方式为C++，Visualization Toolkit，

一种基于边界元模型和局部区域修正的肺部呼吸运动获取方法，包括依次执行的以下步骤：

S1：获取同一人体的两组不同呼吸时刻的胸部CT图像；

S2：处理从步骤S1中获取的两组胸部CT图像，对其进行肺部组织的图像分割，得出肺部表面模型以及肺部血管树模型；

S3：将所述两组胸部CT图像进行配准，处理后得出肺部组织在不同呼吸时刻对应的序列图像；

S4：根据配准后的序列图像信息，提取出肺部组织呼吸运动的边界条件，从而获取肺部表面的运动信息；

S5：利用边界元理论和边界元条件迭代求解肺部表面模型中的各参数并存储，从而获取肺部表面随呼吸运动的边界元运动模型，简称肺部表面模型；

S6：在步骤S2中提取出的肺部血管树模型上设定10个以上的特征点，以这些特征点作为肺内部组织运动局域修正的依据；

S7：肺内部组织运动信息获取：以肺内部组织任一点作为目标点，由该点向肺部表面作六个方向的射线，获得射线与肺部表面模型的六个交点的运动信息，结合这六个交点的运动信息，通过线性插值法算出目标点的运动信息；

S8：目标点运动信息修正：确定肺部表面模型上与目标点距离最近的一个点，称该点为表面点；确定与目标点距离最近的两个特征点；以表面点信息以及该两个特征点信息对目标点的运动信息进行修正，最终得到整个肺部的运动模型。

边界元法(boundary element method)，又称边界元理论，是一种继有限元法之后发展起来的一种新数值方法，与有限元法在连续体域内划分单元的基本思想不同，边界元法是只在定义域的边界上划分单元，用满足控制方程的函数去逼近边界条件。它与基于偏微分方程的区域解法相比，由于降低了问题的维数，而显著降低了自由度数，边界的离散也比区域的离散方便得多，可用较简单的单元准确地模拟边界形状，最终得到阶数较低的线性代数方程组。

本发明首次将边界元模型应用于肺部呼吸运动(信息)的获取(其流程参照图1)，只需对人体拍两组CT图像，通过图像配准后提取出肺部组织呼吸运动的边界条件，便可求解出肺部表面模型；而目前国外能精确得出肺部组织运动模型的方法，如光流场方法，一般需要在短时间内对人体进行连续十次的CT拍照，得到十组CT图像，随后通过有限元法进行求解，得出整个肺部组织的运动模型，但有限元属于区域法，其剖分涉及到整个区域，而边界元只需对边界离散，因此边界元法相比于有限元法，可以降低求解问题的维数。由于原理上的根本差别，本发明对于目标点信息的获取速度将大大优于现有技术；一般来说，本发明可以通过几秒钟获取到目标点的信息，而现有技术则需要花几个小时。

应用边界元理论的计算，需要获取计算对象的边界条件，而肺部组织呼吸运动的边界条件可以通过呼吸序列图像中肋骨的内部边界来确定。

本发明特别适用于直径30mm的局部区域的肺部呼吸运动信息的获取，由于肿瘤的直径一般在10～15mm之间，因此本发明能为穿刺活检手术提供极高精度的位置信息。

本发明中，特征点信息的采集一般通过探针引导微型传感器进入肺部来实现，由此获得特征点的真实坐标；当然，本领域技术人员也可以采取其他常规手段实现该目的，在此不作特别限制。特征点要求最少选定10个以上，可以确保获得充分的空间特征，进而保证精度；在更高精度要求的手术中，为保证肺内部组织运动得到精确的修正，根据实验经验可以随机设定约三十个特征点。

步骤S7中所述“六个方向”一般指正上、正下、正左、正右、正前、正后面方向；采用六个方向的射线，目的是为了确定目标点相对于肺部表面模型的三维位置关系。当然本领域技术人员也可以根据其他方法来进行目标点与肺部表面模型相对位置的判断，在此不作限制。

在一个优选实施例中，血管树的末端由于受CT清晰度以及断层扫描的影响，其位置不能保证精确，因而步骤S6中所述特征点位于血管交叉处，位置唯一确定。

在一个优选实施例中，步骤S3中配准的公式为：

式中，表示为每一个像素所在位置坐标的下舍入数值；同理得到系数 为控制点n_i×n_j×n_k的位移量，B样条的基函数表示为：

其中，0≤t＜1。

图像配准也可以采用其他现有技术所公开的方法进行，上述公式仅属于本发明的实施例之一，其本身并不对图像配准所可采取的其他方法或者公式进行限制或约束。

在一个优选实施例中，步骤S5中边界元三维弹性体的位移与应力基本解可表示为：

其中，υ为泊松系数，σ_ij为应变，G＝E/2(1+υ)，r代表P、Q两点间的距离，n表示表面的向量，δ为Dirac函数，以Q′为圆心半径为ε＝ε→0的球面。

优选地，步骤S5中边界元的离散边界积分方程可表示为：

其中，

边界Γ分为n各子边界Γ_j。

在一个优选实施例中，步骤S7中线性插值法公式为：

其中，目标点设为X，其坐标为(x，y，z)，P_xl，P_xr，P_ya，P_yp，P_zs，P_zi分别为该点向六个方向所作射线与肺部表面模型的交点，d(R_x)、d(R_y)、d(R_z)为X随呼吸运动的位移量，d(P_xl)、d(P_xr)、d(P_ya)、d(P_yp)、d(P_zs)、d(P_zi)为六个交点随呼吸运动的位移量。

在一个优选实施例中，步骤S8中目标点运动信息修正公式表示为：

其中，d(R)为线性插值的结果，由d(R_x)、d(R_y)、d(R_z)可得出；E_s为与X最近的表面点线性插值结果与边界元计算得到结果的数值差；E_a和E_b为两个特征点实际位移距离与插值得到的结果直接的数值差；D_s为X与最近表面点的欧氏距离；D_a和D_b为X与两个特征点的欧氏距离。

上述所有的公式或方程仅属于本发明实施中的一部分例子，本领域技术人员也可以采用其他现有的公式或方程进行计算，但均应落入本发明保护的范围之内。

上述列举的各种实施例，在不矛盾的前提下，可以相互组合实施，本领域技术人员可结合附图和上文对实施例的解释，作为对不同实施例中的技术特征进行组合的依据。

需要指出的是，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.基于边界元模型和局部区域修正的肺部呼吸运动获取方法，其特征在于：包括依次执行的以下步骤：

S1：获取同一人体的两组不同呼吸时刻的胸部CT图像；

S2：处理从步骤S1中获取的两组胸部CT图像，对其进行肺部组织的图像分割，得出肺部表面模型以及肺部血管树模型；

S3：将所述两组胸部CT图像进行配准，处理后得出肺部组织在不同呼吸时刻对应的序列图像；

S4：根据配准后的序列图像信息，提取出肺部组织呼吸运动的边界条件，从而获取肺部表面的运动信息；

S5：利用边界元理论和边界元条件迭代求解肺部表面模型中的各参数并存储，从而获取肺部表面随呼吸运动的边界元运动模型，简称肺部表面模型；

S6：在步骤S2中提取出的肺部血管树模型上设定10个以上的特征点，以这些特征点作为肺内部组织运动局域修正的依据；

S7：肺内部组织运动信息获取：以肺内部组织任一点作为目标点，由该点向肺部表面作六个方向的射线，获得射线与肺部表面模型的六个交点的运动信息，结合这六个交点的运动信息，通过线性插值法算出目标点的运动信息；

S8：目标点运动信息修正：确定肺部表面模型上与目标点距离最近的一个点，称该点为表面点；确定与目标点距离最近的两个特征点；以表面点信息以及该两个特征点信息对目标点的运动信息进行修正，最终得到整个肺部的运动模型。

2.根据权利要求1所述的肺部呼吸运动获取方法，其特征在于：步骤S6中，所述特征点位于血管交叉处。

3.根据权利要求2所述的肺部呼吸运动获取方法，其特征在于：步骤S3中配准的公式为：

式中，表示为每一个像素所在位置坐标的下舍入数值；同理得到系数 为控制点n_i×n_j×n_k的位移量，B样条的基函数表示为：

$B_0\left(t\right)=\frac{{(1-t)}^3}{6},$

$B_1\left(t\right)=\frac{3t^3-6t^2+4}{6},$

$B_2\left(t\right)=\frac{-3t^3+3t^2+3t+1}{6},$

$B_3\left(t\right)=\frac{t^3}{6};$

其中，0≤t＜1。

4.根据权利要求3所述的肺部呼吸运动获取方法，其特征在于：步骤S5中边界元三维弹性体的位移与应力基本解可表示为：

$u_{lk}^{*}(P,Q)=\frac{1}{16\mathrm{\π}(1-\mathrm{\υ})Gr}\[(3-4\mathrm{\υ}){\mathrm{\δ}}_{ij}+r_{,l}{r}_{,k}\],$

$p_{lk}^{*}(P,Q^{\′})=-\frac{1}{8\mathrm{\π}(1-\mathrm{\υ})r^2}\{\frac{\∂r}{\∂n}\[\left(1-2\mathrm{\υ}\right){\mathrm{\δ}}_{lk}+3r_{,l}{r}_{,k}\]-(1-2\mathrm{\υ})(r_{,l}{r}_{,k}-r_{,k}{r}_{,l})\},$

其中，υ为泊松系数，σ_ij为应变，G＝E/2(1+υ)，r代表P、Q两点间的距离，n表示表面的向量，δ为Dirac函数，以Q′为圆心半径为ε＝ε→0的球面。

5.根据权利要求4所述的肺部呼吸运动获取方法，其特征在于：步骤S5中边界元的离散边界积分方程可表示为：

${\mathrm{Cu}}_i+{\mathrm{\Σ}}_{j=1}^n(\∫{\mathrm{\Γ}}_j{P}^{*}(P_i^{\′},Q^{\′}\left)d\mathrm{\Γ}\right(Q^{\′}\left)\right)u_j={\mathrm{\Σ}}_{j=1}^n(\∫{\mathrm{\Γ}}_j{U}^{*}(P_i^{\′},Q^{\′}\left)d\mathrm{\Γ}\right(Q^{\′}\left)\right)p_j,$

其中，

边界Γ分为n各子边界Γ_j。

6.根据权利要求5所述的肺部呼吸运动获取方法，其特征在于：步骤S7中线性插值法公式为：

$d\left(R_x\right)\≈\frac{x_{P_{xr}}-x}{x_{P_{xr}}-x_{P_{xl}}}d\left(P_{xl}\right)+\frac{x-x_{P_{xl}}}{x_{P_{xr}}-x_{P_{xl}}}d\left(P_{xr}\right),$

$d\left(R_y\right)\≈\frac{y_{P_{ya}}-y}{y_{P_{ya}}-y_{P_{yp}}}d\left(P_{ya}\right)+\frac{y-y_{P_{yp}}}{y_{P_{ya}}-y_{P_{yp}}}d\left(P_{yp}\right),$

$d\left(R_z\right)\≈\frac{z_{zs}-z}{z_{P_{zs}}-z_{P_{zi}}}d\left(P_{zs}\right)+\frac{z-z_{P_{zi}}}{z_{P_{zs}}-z_{P_{zi}}}d\left(P_{zi}\right),$

其中，目标点设为X，其坐标为(x，y，z)，P_xl，P_xr，P_ya，P_yp，P_zs，P_zi分别为该点向六个方向所作射线与肺部表面模型的交点，d(R_x)、d(R_y)、d(R_z)为X随呼吸运动的位移量，d(P_xl)、d(P_xr)、d(P_ya)、d(P_yp)、d(P_zs)、d(P_zi)为六个交点随呼吸运动的位移量。

7.根据权利要求6所述的肺部呼吸运动获取方法，其特征在于：步骤S8中目标点运动信息修正公式表示为：

$d{\left(R\right)}_c=\frac{\frac{1}{D_s}}{\frac{1}{D_s}+\frac{1}{D_a}+\frac{1}{D_b}}\·E_s+\frac{\frac{1}{D_a}}{\frac{1}{D_s}+\frac{1}{D_a}+\frac{1}{D_b}}\·E_a+\frac{\frac{1}{D_b}}{\frac{1}{D_s}+\frac{1}{D_a}+\frac{1}{D_b}}\·E_b+d\left(R\right),$

其中，d(R)为线性插值的结果，由d(R_x)、d(R_y)、d(R_z)可得出；E_s为与X最近的表面点线性插值结果与边界元计算得到结果的数值差；E_a和E_b为两个特征点实际位移距离与插值得到的结果直接的数值差；D_s为X与最近表面点的欧氏距离；D_a和D_b为X与两个特征点的欧氏距离。