CN106952285A - 基于先验统计运动模型及自适应配准的肺部运动估计方法 - Google Patents

Abstract

本发明涉及一种基于先验统计运动模型及自适应配准的肺部运动估计方法，包括：S1：获取多组病人呼吸始末两个相位的胸部CT图像；S2：对其进行肺部组织的图像分割；S3：将每组病人两个不同呼吸相位的胸部CT图像进行配准；S4：提取肺实质的运动信息，并以此建立肺部呼吸运动样本库，每组样本的EI时刻CT图像都与第一组样本的EI时刻CT图像做配准，获取样本间的转换矩阵；S5：采集一组新增肺部运动模型EI时刻的全局肺区CT及一组EE时刻的肺部病灶区域CT；S6：在新增的两种CT中选取十组的特征点对；S7：稀疏先验运动模型的生成；S8：病灶区域的精确运动信息获取。与现有技术相比，本发明具有计算速度快、实现方便、应用灵活等优点。

Description

基于先验统计运动模型及自适应配准的肺部运动估计方法

技术领域

本发明涉及肺部呼吸运动信息获取领域，尤其是涉及一种基于先验统计运动模型及自适应配准的肺部运动估计方法。

背景技术

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

经过对现有文献的检索发现，现有的肺部运动模型按其建立方法主要分为三大类：基于图像灰度的数学模型，生物力学模型和统计学模型。基于图像的方法主要是运用图像配准的技术，通过使两幅或多幅图像之间达到最大相似度，从而得到一个肺部运动的流场。假定序列图像中相邻图像间的时间间隔很短，同时相邻图像间的差异也很小，通过构建一个包含时间和空间梯度的信息函数达到配准的目的。这类方法在建立图像与运动数学关系时，需要多个时刻肺部的CT数据，甚至需要采集患者的4DCT才可以完成上述目的。这势必对患者产生了大量的辐射量。并且由于图像配准要处理很大的计算量，因此，它的实时性就会收到制约。基于生物力学模型是把两个时刻肺部的模型从CT数据中重建出，然后将其中一个时刻模型上的顶点和三角面移动到另一个时刻模型上对应的位置实现呼吸模拟。基于生物力学模型的方法通过建立能量方程的形式代替相似度作为评判肺部呼吸运动的依据，然而复杂的组织间挤压情况，软组织的各向异性和非均一性都使得这类方法在获取生物属性参数及边界条件上存在一定困难，无法快速生成各异性的力学模型。

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

1、利用不同时刻CT图像信息，结合配准方法得到多组样本的运动信息，以此建立肺部呼吸运动模型的数据库。

2、针对当前运动统计模型无法体现个体差异或保留局部呼吸运动细节的不足，本研究将对统计模型的建立进行探索。

3、针对个体化统计学模型生后存在的估计误差，结合二次局部CT图像信息对感兴趣区域的运动进行进一步修正，以提高统计学模型在肺部运动获取中的精度。

发明内容

本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种基于先验统计运动模型及自适应配准的肺部运动估计方法。

本发明的目的可以通过以下技术方案来实现：

一种基于先验统计运动模型及自适应配准的肺部运动估计方法，包括依次执行的以下步骤：

S1：事先获取多组病人呼吸始末两个相位的胸部CT图像，其中两个相位分别为吸气末EI、呼气末EE；

S2：处理从步骤S1中获取的胸部CT图像，对其进行肺部组织的图像分割，得到肺实质；

S3：将每组病人两个不同呼吸相位的胸部CT图像进行配准；

S4：根据配准后的序列图像信息，提取肺实质的运动信息，并以此建立肺部呼吸运动样本库，每组样本的EI时刻CT图像都与第一组样本的EI时刻CT图像做配准，获取样本间的转换矩阵；

S5：采集一组新增肺部运动模型EI时刻的全局肺区CT及一组EE时刻的肺部病灶区域CT；

S6：在新增的两种CT中选取十组的特征点对，获取这些特征点对的运动信息作为步骤S7的输入；

S7：稀疏先验运动模型的生成：将步骤S5中EI时刻的CT与运动样本库中的第一组CT数据进行配准，得到两者之间的转换矩阵，结合步骤S6中的特征点对运动信息获取其在运动样本库中各个样本中的各异性表达，计算获得步骤S6提供的运动信息在运动样本库中的稀疏线性表达，并以此生成肺实质的运动信息；

S8：病灶区域的精确运动信息获取：统计运动模型提供的先验运动信息作为依据，结合步骤S5中的两组图像利用配准算法得到病灶区域的精确运动信息。

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

式中，Ω为图像区域；I_p,EE和I_p,EI为样本库中每个病例在呼气末和吸气末的CT图像；x_i为图像像素的位置；d为像素的运动向量；|Δd|为运动向量的梯度，α为正则化系数。

优选地，所述的步骤S4中样本间的对应关系获取描述为：

采用S3中的配准方法对每组样本的EI时刻CT图像都与第一组样本的EI时刻CT图像进行配准，获取样本间的转换矩阵得到运动样本库中每组样本I_p,EI对应第一组样本EI时刻CT图像I_1,EI的关系为：

优选地，所述的步骤S7中稀疏先验运动模型的生成表示为：

其中y_j为步骤S6中特征点对的运动信息；M为这些特征点在样本库每组样本中对应位置的运动信息，其中k为特征点对数量，n为样本数量；x为一组样本线性表达的系数；e为输入y_j的误差；λ₁和λ₂控制x和e的稀疏量。

优选地，所述的步骤S8中病灶区域的精确运动信息获取表示为：

其中，Ω为图像区域；为B样条配准中控制点影响系数；I_j,EE和I_j,EI为当前病例在呼气末和吸气末的CT图像；g_prior为步骤S7提供的先验统计运动信息；g(X)为B样条控制点的运动量，g_n+1(X)和g_n(X)为第n+1次和第n次配准时控制点的运动量；λ为控制点运动量正则化系数。

优选地，所述的步骤S8中B样条配准中控制点影响系数的表达公式为：

其中，R_i,j,k为第ijk个控制点在EI时刻图像上所影响的区域；为其在EE时刻图像上所影响的区域，为ξ_i,j,k的集合。

与现有技术相比，本发明具有以下优点：

1、采用配准多组病例不同时刻CT图像的方法获取肺部呼吸运动信息，从而建立运动数据库。该运动数据库能够较好的描述不同呼吸模式下肺部运动信息，解决了肺部运动的各向异性及非均一性问题所引起的问题；

2、稀疏先验运动模型有效地反映了新增病例的肺部局部呼吸运动细节；

3、针对稀疏先验运动模型不能提供足够精确的运动信息，本发明基于先验统计运动模型及自适应配准的肺部运动估计方法实现了对肺实质病灶组织运动的精确计算方法。

附图说明

图1为本发明的流程图。

具体实施方式

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

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

一种基于先验统计运动模型及自适应配准的肺部运动估计方法，包括依次执行的以下步骤：

S1：事先获取多组病人呼吸始末两个相位(吸气末EI、呼气末EE)的胸部CT图像；

S2：处理从步骤S1中获取的胸部CT图像，对其进行肺部组织的图像分割，得到肺实质；

S3：将每组病人两个不同呼吸相位的胸部CT图像进行配准；

S4：根据配准后的序列图像信息，提取肺实质的运动信息，并以此建立肺部呼吸运动样本库。每组样本的EI时刻CT图像都与第一组样本的EI时刻CT图像做配准，获取样本间的转换矩阵；

S5：采集一组新增肺部运动模型EI时刻的全局肺区CT及一组EE时刻的肺部病灶区域CT；

S6：在新增的两种CT中选取约十组的特征点对，获取这些特征点对的运动信息作为步骤S7的输入；

S7：稀疏先验运动模型的生成：将步骤S5中EI时刻的CT与运动样本库中的第一组CT数据进行配准，得到两者之间的转换矩阵。结合步骤S6中的特征点对运动信息获取其在运动样本库中各个样本中的各异性表达。计算获得步骤S6提供的运动信息在运动样本库中的稀疏线性表达，并以此生成肺实质的运动信息；

S8：病灶区域的精确运动信息获取：统计运动模型提供的先验运动信息作为依据，结合第五步中的两组图像利用配准算法得到病灶区域的精确运动信息。

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

式中，Ω为图像区域；I_p,EE和I_p,EI为样本库中每个病例在呼气末和吸气末的CT图像；x_i为图像像素的位置；d为像素的运动向量；|Δd|为运动向量的梯度，α为正则化系数。

在一个优选实施例中，步骤S4中样本间的对应关系获取可描述为：

采用S3中的配准方法对每组样本的EI时刻CT图像都与第一组样本的EI时刻CT图像进行配准，获取样本间的转换矩阵得到运动样本库中每组样本对应第一组样本EI时刻CT图像的关系为：

在一个优选实施例中，步骤S7中稀疏先验运动模型的生成可表示为：

其中,y_j为步骤S6中特征点对的运动信息；M为这些特征点在样本库每组样本中对应位置的运动信息，其中k为特征点对数量，n为样本数量；x为一组样本线性表达的系数；e为输入y_j的误差；λ₁和λ₂控制x和e的稀疏量。

在一个优选实施例中，步骤S8中病灶区域的精确运动信息获取可表示为：

其中，Ω为图像区域；为B样条配准中控制点影响系数；I_j,EE和I_j,EI为当前病例在呼气末和吸气末的CT图像；g_prior为步骤S7提供的先验统计运动信息；g(X)为B样条控制点的运动量，g_n+1(X)和g_n(X)为第n+1次和第n次配准时控制点的运动量；λ为控制点运动量正则化系数。

在一个优选实施例中，步骤S8中病灶区域的精确运动信息获取可表示为：步骤S8中B样条配准中控制点影响系数的表达公式为：

其中，R_i，j，k为第ijk个控制点在EI时刻图像上所影响的区域；为其在EE时刻图像上所影响的区域。为ξ_i,j,k的集合。

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

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

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到各种等效的修改或替换，这些修改或替换都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以权利要求的保护范围为准。

Claims (6)

1.一种基于先验统计运动模型及自适应配准的肺部运动估计方法，其特征在于，包括依次执行的以下步骤：

S1：事先获取多组病人呼吸始末两个相位的胸部CT图像，其中两个相位分别为吸气末EI、呼气末EE；

S2：处理从步骤S1中获取的胸部CT图像，对其进行肺部组织的图像分割，得到肺实质；

S3：将每组病人两个不同呼吸相位的胸部CT图像进行配准；

S4：根据配准后的序列图像信息，提取肺实质的运动信息，并以此建立肺部呼吸运动样本库，每组样本的EI时刻CT图像都与第一组样本的EI时刻CT图像做配准，获取样本间的转换矩阵；

S5：采集一组新增肺部运动模型EI时刻的全局肺区CT及一组EE时刻的肺部病灶区域CT；

S6：在新增的两种CT中选取十组的特征点对，获取这些特征点对的运动信息作为步骤S7的输入；

S7：稀疏先验运动模型的生成：将步骤S5中EI时刻的CT与运动样本库中的第一组CT数据进行配准，得到两者之间的转换矩阵，结合步骤S6中的特征点对运动信息获取其在运动样本库中各个样本中的各异性表达，计算获得步骤S6提供的运动信息在运动样本库中的稀疏线性表达，并以此生成肺实质的运动信息；

S8：病灶区域的精确运动信息获取：统计运动模型提供的先验运动信息作为依据，结合步骤S5中的两组图像利用配准算法得到病灶区域的精确运动信息。

2.根据权利要求1所述的一种基于先验统计运动模型及自适应配准的肺部运动估计方法，其特征在于，所述的步骤S3中配准所需的能量公式E(d)为：

$E\left(d\right)=\underset{x_i\∈\mathrm{\Ω}}{\Σ}{\[I_{p,EE}\left(x_i\right)-I_{p,EI}(x_i+d)\]}^2+\mathrm{\α}|\mathrm{\Δ}d{|}^2,$

式中，Ω为图像区域；I_p,EE和I_p,EI为样本库中每个病例在呼气末和吸气末的CT图像；x_i为图像像素的位置；d为像素的运动向量；|Δd|为运动向量的梯度，α为正则化系数。

3.根据权利要求1所述的一种基于先验统计运动模型及自适应配准的肺部运动估计方法，其特征在于，所述的步骤S4中样本间的对应关系获取描述为：

采用S3中的配准方法对每组样本的EI时刻CT图像都与第一组样本的EI时刻CT图像进行配准，获取样本间的转换矩阵得到运动样本库中每组样本I_p,EI对应第一组样本EI时刻CT图像I_1,EI的关系为：

4.根据权利要求1所述的一种基于先验统计运动模型及自适应配准的肺部运动估计方法，其特征在于，所述的步骤S7中稀疏先验运动模型的生成表示为：

$\arg \underset{xe}{min}\left|\right|y_i-Mx-e|{|}_2^2+{\mathrm{\λ}}_1|\left|x\right|{|}_1+{\mathrm{\λ}}_2\left|\right|e|{|}_1,$

其中y_j为步骤S6中特征点对的运动信息；M为这些特征点在样本库每组样本中对应位置的运动信息，其中k为特征点对数量，n为样本数量；x为一组样本线性表达的系数；e为输入y_j的误差；λ₁和λ₂控制x和e的稀疏量。

5.根据权利要求1所述的一种基于先验统计运动模型及自适应配准的肺部运动估计方法，其特征在于，所述的步骤S8中病灶区域的精确运动信息获取表示为：

$\begin{array}{c}E=\underset{X\∈\mathrm{\Ω}}{\mathrm{\Σ}}\left(1-{\mathrm{\ξ}}_{I_c}\right){\[I_{j,EE}\left(X\right)-I_{j,EI}\left(X+g_{prior}\left(X\right)+g\left(X\right)\right)\]}^2+\underset{X\∈\mathrm{\Ω}}{\mathrm{\Σ}}{\mathrm{\ξ}}_{I_c}g{\left(X\right)}^2\\ +\mathrm{\λ}\underset{X\∈\mathrm{\Ω}}{\mathrm{\Σ}}{\left(g_{n+1}\left(X\right)-g_n\left(X\right)\right)}^2,\end{array}$

其中，Ω为图像区域；为B样条配准中控制点影响系数；I_j,EE和I_j,EI为当前病例在呼气末和吸气末的CT图像；g_prior为步骤S7提供的先验统计运动信息；g(X)为B样条控制点的运动量，g_n+1(X)和g_n(X)为第n+1次和第n次配准时控制点的运动量；λ为控制点运动量正则化系数。

6.根据权利要求5所述的一种基于先验统计运动模型及自适应配准的肺部运动估计方法，其特征在于，所述的步骤S8中B样条配准中控制点影响系数的表达公式为：

${\mathrm{\ξ}}_{i,j,k}=\frac{R_{i,j,k}-R_{{\mathrm{EE}}_{i,j,k}}}{R_{i,j,k}},$

其中，R_i,j,k为第ijk个控制点在EI时刻图像上所影响的区域；为其在EE时刻图像上所影响的区域，为ξ_i,j,k的集合。