CN103761745A - 一种肺部运动模型估计方法及系统 - Google Patents

Abstract

本发明适用于图像处理领域，提供了一种肺部运动模型估计方法及系统，所述方法包括：提取参考相位的肺部三维CT图像中控制点；将所述控制点在上一时刻的跟踪位置作为预测位置，跟踪所述控制点在连续离散相位三维CT图像的位置；对所述控制点在连续离散相位的跟踪轨迹进行最小二乘曲线拟合，确定所述控制点在任意连续时刻的运动位置；根据所述控制点在任意连续时刻的运动位置构建当前时刻肺部的三维形变模型，得到参考相位肺部的任意体素在任意时间的运动位置估计。本发明采用了基于标志点的图像配准算法，该方法计算量小，计算速度快，可以有效估计肺部运动模型。

Description

一种肺部运动模型估计方法及系统

技术领域

本发明属于图像处理领域，尤其涉及一种肺部运动模型估计方法及系统。 

背景技术

目前三维医学成像技术的应用领域十分广泛。X射线计算机断层摄影(Computed tomography，CT)成像技术是利用特殊的X射线设备来获取人体的组织结构图像技术，利用三维CT成像技术，可以根据病人的生理情况和病理特征重建三维医学图像。由于人体的器官会随着呼吸周期而循环往复地运动，三维成像技术无法处理连续时间序列的器官运动模型。 

最近兴起的四维医学成像技术，在三维数据的基础上增加了时间轴，解决了三维成像技术在时间轴上的局限性。四维CT成像是连续周期的序列三维CT图像集，可以用于构造肺部呼吸运动模型。 

肺部呼吸运动模型在精确放射性治疗上具有重要的意义。在肺部肿瘤的放射治疗过程中，如3D适形放射治疗（3DCRT）、调强放射治疗（IMRT）、立体定向放射治疗（SBRT），可以使放射剂量集中在肿瘤位置，同时将周围的敏感组织进行隔离。但由于呼吸运动产生肿瘤位移，传统方法需要在放疗过程中扩大放射的范围，这就会对肿瘤周围的正常组织造成伤害。精确放射治疗是通过与靶区形状一致的精确辐射剂量分布来减少照射体积，通过对器官运动估计，指引照射线准确作用在肿瘤位置。准确构建肺部运动模型，是肺部疾病放射性治疗努力解决的关键问题，也是图像引导放疗技术（IGRT）亟待解决的问题。 

肺部运动模型估计可以通过四维图像配准来实现，通过图像配准构建肺部连续形变模型。四维医学图像配准技术不同于三维医学图像配准技术，它是在三维医学图像的基础上增加了时间维度，需要对连续时间序列的三维体数据进 行综合性的配准，而不是像三维医学图像配准技术那样只是在一对体数据间进行各种配准算法。目前，已经有一些学者针对四维图像数据的特点研究了一些适合不同领域的算法。Perperidis等提出了一个时空非线性配准算法，他们在两个心脏序列之间确定形变来建立物体形状和运动的概率集。GJ Kein等研究出一个含有12个参数的全局放射运动模型，对舒张末期的心脏PET序列中不同的呼吸门进行配准。D.Sarrut、D.Yang、D.Sarrut、Peyrat、J.M.等也提出了一些基于或类似于光流法的图像配准方法。其他的配准方法还有Yan,Brock提出的利用FEM(有限元建模)的生物机械模型法。 

目前根据研究情况，可将肺部呼吸运动模型的四维医学图像配准方法分为两类：基于图像特征的方法和基于图像强度的方法。 

基于图像特征的配准算法利用匹配的标志点对集和插值模型来确定图像之间的位移量，这里通常使用自动或手动的方法来确定图像之间对应的标志点对，同时利用插值模型来确定感兴趣区域的所有体素的形变位移。这类方法是基于控制点的，避免了求解复杂的方程。然而，人工选取数量庞大的对应标志点对于临床应用来说是不切实际的，因此自动地选取标志点是研究的热门领域。 

基于图像强度的方法是根据一个数学模型来描述图像间的位移。Barron和Beauchemin利用光流方法计算三维肺部数据间的形变；为了获得呼吸运动的信息，Castillo,Li等人使用压缩流方程来模拟四维CT肺组织图像的运动位移；Lucas和Kanade提出了一个适合于较大形变位移的光流配准算法。 

综上所述，目前肺部呼吸运动模型使用的四维图像配准方法通常具有计算量很大，收敛时间很长等缺点，不能够快速地实现肺部运动模型的估计。 

发明内容

本发明所要解决的技术问题在于提供一种肺部运动模型估计方法及系统，旨在解决四维CT成像技术中肺部运动模型的估计问题。 

本发明是这样实现的，一种肺部运动模型估计方法，所述方法包括： 

提取参考相位的肺部三维CT图像中控制点； 

将所述控制点在上一时刻的跟踪位置作为预测位置，跟踪所述控制点在连续离散相位三维CT图像的位置； 

对所述控制点在连续离散相位的跟踪轨迹进行最小二乘曲线拟合，确定所述控制点在任意连续时刻的运动位置； 

根据所述控制点在任意连续时刻的运动位置构建当前时刻肺部的三维形变模型，得到参考相位肺部的任意体素在任意时间的运动位置估计。 

本发明还提供了一种肺部运动模型估计系统，所述系统包括：检测单元、跟踪单元、拟合单元以及构建单元；其中， 

所述检测单元，用于提取参考相位的肺部三维CT图像中控制点； 

所述跟踪单元，用于将所述控制点在上一时刻的跟踪位置作为预测位置，跟踪所述控制点在连续离散相位三维CT图像的位置； 

所述拟合单元，用于对所述控制点在连续离散相位的跟踪轨迹进行最小二乘曲线拟合，确定所述控制点在任意连续时刻的运动位置； 

所述构建单元，用于根据所述控制点在任意连续时刻的运动位置构建当前时刻肺部的三维形变模型，得到参考相位肺部的任意体素在任意时间的运动位置估计。 

本发明与现有技术相比，有益效果在于：通过快速确定控制点，对所述控制点进行快速跟踪，得到控制点在连续相位肺部三维CT图像中的运动位置，确定了构建肺部运动模型所需要的控制点运动轨迹；利用最小二乘法对跟踪的控制点离散运动轨迹进行曲线拟合，既可以得到控制点在任意连续时刻的运动位置，同时也可以修正跟踪误差，提高参考相位图像中的控制点与其他相位图像中的控制点之间的对应精度；利用任意连续时刻的控制点位置构建任意连续时刻的肺部形变模型，将时间维信息引入肺部形变过程，构建了包含时间信息的平滑的肺部运动模型；能够解决精确放射治疗中的呼吸运动估计问题，在自由呼吸的状态下，精确控制肿瘤放射治疗规划中的放射剂量，有效提高放射治 疗的效果，减少对病人正常组织的损伤，该方法计算量小，计算速度快。 

附图说明

图1是本发明实施例提供的肺部运动模型估计方法的实现流程图； 

图2是本发明实施例提供的肺部运动模型估计系统的结构示意图。 

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。 

图1示出了本发明实施例提供的肺部运动模型估计方法的实现流程，如图1所示，所述详述如下： 

步骤S101，提取参考相位的肺部三维CT图像中的控制点。 

在三维CT图像中，I_x,I_y,I_z分别代表三维CT图像中当前体素

沿x,y,z方向的梯度，计算该体素的自相关矩阵C： 

根据 $C\left(\stackrel{\→}{x}\right)=\left[\begin{array}{ccc}{I}_x^2& I_x{I}_y& I_x{I}_z\\ I_x{I}_y& I_y^2& I_y{I}_z\\ I_x{I}_z& I_y{I}_z& I_z^2\end{array}\right]$ 计算所肺部三维CT图像中体素

的自相关矩阵C，其中，I_x,I_y,I_z分别是三维CT图像在(x,y,z)处沿x方向、y方向和z方向的一阶导；以所述体素

为中心和原点，利用高斯平滑函数

对局部邻域内的各体素的自相关矩阵C的各项分别进行平滑，(x,y,z)是高斯平滑窗口内各体素坐标，

是高斯函数的方差。Harris角点（也称为Harris算子）可以定义为下式的局部区域最大值：R＝Det(C)-k·Trace(C)，其中，Det(C)是矩阵C的行列式值，Trace(C)是矩阵C的迹，k为调整系数。当体素

的Harris算子R大于预先设定的阈值T时,则该体素为控制点。 

步骤S102，将所述控制点在上一时刻的跟踪位置作为预测位置，利用mean shift迭代方法跟踪该控制点在连续离散相位三维CT图像的位置。 

以参考相位三维CT图像提取的控制点为跟踪目标，分别跟踪其在其他离散相位三维CT图像中对应的位置。根据呼吸阶段的特性，不管是从最大吸气阶段到最大呼气阶段，还是从最大呼气阶段到最大吸气阶段都是一个循序渐进的过程，也就是每个特征点在各离散相位的三维CT图像中的位置具有稳定的运动趋势，根据这种运动趋势，可以预测参考相位三维CT图像中各控制点的初始位置，即将前一时刻的跟踪结果集作为后一时刻的预测点集即预测位置，可以体现四维数据的连续性。 

假设y₀是当前相位三维CT图像的控制点，假设y^*是参考相位肺部三维CT图像的控制点，目标模式的概率密度方程为： 

其中δ为Kronecher delta函数，

是点间欧几里德距离。{x_i}_i=1...n为以y^*为中心的三维目标模型的体素坐标，x_i=(x,y,z)，函数b:R³→{1...m}将位置x_i处的体素进行量化。搜索区域的直径h=[h_x,h_y,h_z]。C_h是使 

满足的归一化常量，n是三维CT图像局部邻域内体素的个数。对于以y₀为中心的目标候选模式的概率密度函数计算方法相同， 

$p_u\left(y_0\right)=C_p\underset{i=1}{\overset{n}{\mathrm{\Σ}}}k\left({\left|\right|(y_0-x_i)/h\left|\right|}^2\right)\mathrm{\δ}[b\left(x_i\right)-u]$

其中，C_p是使

满足的归一化常量，n是三维CT图像局部邻域内体素的个数。利用Bhattacharyya系数计算目标模式和目标候选模式之间的相似程度，控制点y₀的新位置

可以通过以下公式计算得到： 

通过mean shift迭代直至收敛。如果单纯地以Bhattacharyya系数作为相似性测度存在判断偏差，为了进一步提高mean shift的跟踪精度，在mean shift迭代中增加了判断体数据间相似性的过程。即算法根据Bhattacharyya系数进行mean shift迭代，对于每个控制点，保存其在迭代过程中所有的跟踪位置。每个控制点迭代完成后，计算预测点所处局部邻域和所有保存的跟踪位置所处的局部邻域与以y^*为中心参考体的局部邻域的均方误差；选择均方误差较小的跟踪位置作为当前控制点的最终跟踪位置。 

步骤S103，对所述控制点在连续离散相位的跟踪轨迹进行最小二乘曲线拟合，确定所述控制点在任意连续时刻的运动位置。 

为了获取任意连续时刻的三维数据控制点集，需要对连续时刻的离散控制点作基于最小二乘法的曲线拟合。假设某个控制点在离散时刻的位置为(x_s,y_s,z_s)，s=0,1,...,m，m为相位个数。 

将三维数据分别沿x,y,z三个方向作基于最小二乘法的多项式拟合，即分别沿x,y,z三个方向拟合二次多项式： 

$p_x\left(t\right)=\underset{k=0}{\overset{2}{\mathrm{\Σ}}}{a}_k{t}^k,$ $p_y\left(t\right)=\underset{k=0}{\overset{2}{\mathrm{\Σ}}}{b}_k{t}^k,$ $p_z\left(t\right)=\underset{k=0}{\overset{2}{\mathrm{\Σ}}}{c}_k{t}^k$

其中a₀,a₁,a₂，b₀,b₁,b₂和c₀,c₁,c₂分别为拟合二次多项式系数。以拟合多项式系数a₀,a₁,a₂的求解为例，构造线性方程组， 

$\left[\begin{array}{ccc}m+1& \underset{t=0}{\overset{m}{\mathrm{\Σ}}}t& \underset{t=0}{\overset{m}{\mathrm{\Σ}}}{t}^2\\ \underset{t=0}{\overset{m}{\mathrm{\Σ}}}t& \underset{t=0}{\overset{m}{\mathrm{\Σ}}}{t}^2& \underset{t=0}{\overset{m}{\mathrm{\Σ}}}{t}^{2+1}\\ \underset{t=0}{\overset{m}{\mathrm{\Σ}}}{t}^2& \underset{t=0}{\overset{m}{\mathrm{\Σ}}}{t}^{2+1}& \underset{t=0}{\overset{m}{\mathrm{\Σ}}}{t}^4\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\end{array}\right]=\left[\begin{array}{c}\underset{t=0}{\overset{m}{\mathrm{\Σ}}}{x}_t\\ \underset{t=0}{\overset{m}{\mathrm{\Σ}}}{\mathrm{tx}}_t\\ \underset{t=0}{\overset{m}{\mathrm{\Σ}}}{t}^2{x}_t\end{array}\right]$

通过该方程组解出a₀,a₁,a₂。拟合多项式系数b₀,b₁,b₂和c₀,c₁,c₂的求解类似。 

根据基于最小二乘法的多项式曲线拟合方法，可以将连续时刻的离散控制点拟合成平滑的曲线，获取到任意时刻的控制点集(p_x(t),p_y(t),p_z(t)),t∈[0,m]。 

步骤S104，根据所述控制点在任意连续时刻的运动位置构建当前时刻肺部的三维形变模型，得到参考相位肺部的任意体素在任意时间的运动位置估计； 

针对任意时刻的控制点集，将参考相位三维CT图像形变到当前时刻。假设参考相位三维CT图像中的控制点坐标为(x_i(0),y_i(0),z_i(0))，i=1,2,L,N。对任意时刻t，对应的控制点坐标为(x_i(t),y_i(t),z_i(t))，p_x(t)=x_i(t)，p_y(t)=y_i(t)，p_z(t)=z_i(t)，该时刻的形变函数为f_t(x,y,z)， 

$f_t(x,y,z)=a_{1,t}+a_{x,t}x+a_{y,t}y+a_{z,t}z+\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{w}_{i,t}U\left(\right||(x_i\left(t\right),y_i\left(t\right),z_i\left(t\right))-(x,y,z)|\left|\right)$

其中，U(r_i)为薄板样条的基函数U(r_i)=r_i，径向距离r_i定义为：  $r_i=\sqrt{{(x_i\left(t\right)-x)}^2+{(y_i\left(t\right)-y)}^2+{(z_i\left(t\right)-z)}^2}.$

形变参数a_t=(a_1,t,a_x,t,a_y,t,a_z,t)和w_t=(w_1,t,w_2,t,...,w_N,t)由下述公式得到： 

$\left\{\begin{array}{c}K{w}_t+P{a}_t=v\\ P^T{w}_t=O\end{array}\right.$

其中， $v=\left[\begin{array}{c}{x}_1\left(t\right),y_1\left(t\right),z_1\left(t\right)\\ M\\ x_N\left(t\right),y_N\left(t\right),z_N\left(t\right)\end{array}\right],$ K是N×N矩阵：k_ij=U(r_ij)，  $r_{\mathrm{ij}}=\sqrt{{(x_i\left(t\right)-x_j\left(t\right))}^2+{(y_i\left(t\right)-y_j\left(t\right))}^2+{(z_i\left(t\right)-z_j\left(t\right))}^2},$ P是N×4矩阵：  $P=\left[\begin{array}{cccc}1& x_1\left(0\right)& y_1\left(0\right)& z_1\left(0\right)\\ 1& x_2\left(0\right)& y_2\left(0\right)& z_2\left(0\right)\\ L& L& L& L\\ 1& x_N\left(0\right)& y_N\left(0\right)& z_N\left(0\right)\end{array}\right],$ O为4×4的零矩阵。 

图2示出了本发明实施例提供的肺部运动模型估计系统的结构示意，为了便于说明，仅示出了本发明实施例相关的部分。其中，该肺部运动模型估计系统可以是软件单元，硬件单元或者软硬结合的单元。在本实施例中，该肺部运动模型估计系统包括：检测单元21、跟踪单元22、拟合单元23、构建单元24； 

其中， 

所述检测单元21，用于提取参考相位的肺部三维CT图像中控制点； 

所述跟踪单元22，用于将所述控制点在上一时刻的跟踪位置作为预测位置，跟踪所述控制点在连续离散相位三维CT图像的位置； 

所述拟合单元23，用于对所述控制点在连续离散相位的跟踪轨迹进行最小二乘曲线拟合，确定所述控制点在任意连续时刻的运动位置； 

所述构建单元24，用于根据所述控制点在任意连续时刻的运动位置构建当前时刻肺部的三维形变模型，得到参考相位肺部的任意体素在任意时间的运动位置估计。 

进一步地，所述检测单元21，具体用于计算所述肺部三维CT图像中任一体素的自相关矩阵C；利用R=Det(C)-k·Trace(C)计算所述体素的Harris角点，当所述体素的Harris角点大于预先设定的阈值T时，则提取所述体素为控制点；其中，Det(C)是所述自相关矩阵C的行列式值，Trace(C)是所述自相关矩阵C的迹，k为预先设置的调整系数。 

进一步地，所述检测单元21，具体计算所述肺部三维CT图像中任一体素的自相关矩阵C包括： 

根据 $C\left(\stackrel{\→}{x}\right)=\left[\begin{array}{ccc}{I}_x^2& I_x{I}_y& I_x{I}_z\\ I_x{I}_y& I_y^2& I_y{I}_z\\ I_x{I}_z& I_y{I}_z& I_z^2\end{array}\right]$ 计算所肺部三维CT图像中体素

的自相关矩C，其中，I_x,I_y,I_z分别是三维CT图像在(x,y,z)处沿x方向、y方向和z方向的一阶导；以所述体素为中心和原点，利用高斯平滑函数

对局部邻域内的各体素的自相关矩阵C的各项分别进行平滑，(x,y,z)是高斯平滑窗口内各体素坐标，

是高斯函数的方差。 

进一步地，所述跟踪单元22，具体用于计算控制点在当前相位的新位置，根据Bhattacharyya系数对控制点在当前相位进行mean shift迭代，保存每个控 制点在迭代过程中所有跟踪位置；计算预测位置和跟踪位置在当前体中局部邻域与控制点在参考体中的局部邻域的均方误差；选择均方误差较小的跟踪位置作为当前控制点在当前相位的最终跟踪位置；在各连续离散相位重复上述过程，确定所述控制点在各连续离散相位的最终跟踪位置。 

进一步地，所述跟踪单元22，具体用于计算控制点在当前相位的新位置

为： 

$p_u\left(y_0\right)=C_p\underset{i=1}{\overset{n}{\mathrm{\Σ}}}k\left({\left|\right|(y_0-x_i)/h\left|\right|}^2\right)\mathrm{\δ}[b\left(x_i\right)-u];$ 其中，y₀是当前相位肺部三维CT图像的控制点，{x_i}_i=1...n是以y₀为中心的当前相位三维CT图像中局部邻域的体素坐标，x_i=(x,y,z)。h是局部邻域的直径大小， $k\left(\right|\left|x\right|\left|\right)=\frac{3}{4}(1-|\left|x\right|\left|\right),\left|\right|x\left|\right|\≤1,$ ||||是点间欧几里德距离。函数

将位置x_i处的体素强度f(x_i)进行量化，是下取整函数，δ为Kronecher delta函数，y^*是当前参考相位肺部三维CT图像的控制点，其中{x_i}_i=1...n是以y^*为中心的参考相位三维CT图像中局部邻域的体素坐标；C_h是使

满足的归一化常量，C_p是使

满足的归一化常量，n是三维CT图像局部邻域内体素的个数。 

进一步地，所述拟合单元23，具体用于估计参考相位肺部的任意空间位置在任意连续时间的运动位置。控制点在离散时刻的位置为(x_s,y_s,z_s)，s=0,1,...,m，m为相位个数，分别沿x,y,z三个方向拟合二次多项式：

其中a₀,a₁,a₂，b₀,b₁,b₂和c₀,c₁,c₂分别为拟合二次多项式系数，p_x(t)，p_y(t)，p_z(t)分别是x,y,z三个方向拟合的二次多项式； 

通过构造线性方程组，计算得到拟合二次多项式系数a₀,a₁,a₂，b₀,b₁,b₂和c₀,c₁,c₂； 

根据基于最小二乘法的多项式曲线拟合方法，将连续时刻的离散控制点拟合成平滑的曲线，获取到控制点集(p_x(t),p_y(t),p_z(t)),t∈[0,m]。 

进一步地，所述构建单元24，具体用于将参考相位肺部三维图像在任意时刻的控制点形变到当前时刻。假设参考相位三维CT图像中的控制点坐标为(x_i(0),y_i(0),z_i(0))，i=1,2,L,N。对任意时刻t∈[0,m]，对应的控制点坐标为(x_i(t),y_i(t),z_i(t))，p_x(t)=x_i(t)，p_y(t)=y_i(t)，p_z(t)=z_i(t)，该时刻的形变函数为f_t(x,y,z)： 

$f_t(x,y,z)=a_{1,t}+a_{x,t}x+a_{y,t}y+a_{z,t}z+\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{w}_{i,t}U\left(\right||(x_i\left(t\right),y_i\left(t\right),z_i\left(t\right))-(x,y,z)|\left|\right)$

其中，U(r_i)为薄板样条的基函数U(r_i)=r_i，径向距离r_i定义为： 

(x,y,z)为参考相位中任意体素的坐标。 

形变参数a_t=(a_1,t,a_x,t,a_y,t,a_z,t)和w_t=(w_1,t,w_2,t,...,w_N,t)由下述公式得到： 

$\left\{\begin{array}{c}K{w}_t+P{a}_t=v\\ P^T{w}_t=O\end{array}\right.$

其中， $v=\left[\begin{array}{c}{x}_1\left(t\right),y_1\left(t\right),z_1\left(t\right)\\ M\\ x_N\left(t\right),y_N\left(t\right),z_N\left(t\right)\end{array}\right],$ K是N×N矩阵：k_ij=U(r_ij)，  $r_{\mathrm{ij}}=\sqrt{{(x_i\left(t\right)-x_j\left(t\right))}^2+{(y_i\left(t\right)-y_j\left(t\right))}^2+{(z_i\left(t\right)-z_j\left(t\right))}^2},$ P是N×4矩阵：  $P=\left[\begin{array}{cccc}1& x_1\left(0\right)& y_1\left(0\right)& z_1\left(0\right)\\ 1& x_2\left(0\right)& y_2\left(0\right)& z_2\left(0\right)\\ L& L& L& L\\ 1& x_N\left(0\right)& y_N\left(0\right)& z_N\left(0\right)\end{array}\right],$ O为4×4的零矩阵。 

本发明实施例中，利用Harris角点检测方法快速确定控制点，通过利用mean shift迭代方法对所述控制点进行快速跟踪，得到了控制点在连续相位肺部三维CT图像中的运动位置，确定了构建肺部运动模型所需要的控制点运动轨迹；利用最小二乘法对跟踪的控制点离散运动轨迹进行曲线拟合，既可以得到控制点在任意连续时刻的运动位置，同时也可以修正跟踪误差，提高参考相位图像中的控制点与其他相位图像中的控制点之间的对应精度；利用任意连续时刻的控制点位置构建任意连续时刻的肺部形变模型，将时间维信息引入肺部形变过程，构建了包含时间信息的平滑的肺部运动模型；通过快速肺部运动模型估计，可以用于解决精确放射治疗中的呼吸运动估计问题，在病人自由呼吸的状态下，精确控制肿瘤放射治疗规划中的放射剂量，有效提高放射治疗的效果，减少对病人正常组织的损伤。 

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

Claims (14)

1.一种肺部运动模型估计方法，其特征在于，所述方法包括： 

提取参考相位的肺部三维CT图像中控制点； 

将所述控制点在上一时刻的跟踪位置作为预测位置，跟踪所述控制点在连续离散相位三维CT图像的位置； 

对所述控制点在连续离散相位的跟踪轨迹进行最小二乘曲线拟合，确定所述控制点在任意连续时刻的运动位置； 

根据所述控制点在任意连续时刻的运动位置构建当前时刻肺部的三维形变模型，得到参考相位肺部的任意体素在任意时间的运动位置估计。 

2.如权利要求1所述的方法，其特征在于，所述提取参考相位的肺部三维CT图像中控制点包括： 

计算所述肺部三维CT图像中任一体素的自相关矩阵C； 

利用R=Det(C)-k·Trace(C)计算所述体素的Harris角点，当所述体素的Harris角点大于预先设定的阈值T时，则提取所述体素为控制点； 

其中，Det(C)是所述自相关矩阵C的行列式值，Trace(C)是所述自相关矩阵C的迹，k为预先设置的调整系数。 

3.如权利要求2所述的方法，其特征在于，所述计算所述肺部三维CT图像中任一体素的自相关矩阵C包括： 

根据

计算所肺部三维CT图像中体素

的自相关矩阵C，其中，I_x,I_y,I_z分别是三维CT图像在(x,y,z)处沿x方向、y方向和z方向的一阶导；以所述体素为中心和原点，利用高斯平滑函数

对局部邻域内的各体素的自相关矩阵C的各项分别进行平滑，(x,y,z)是高斯平滑窗口内各体素坐标，

是高斯函数的方差。 

4.如权利要求1所述的方法，其特征在于，所述将所述控制点在上一相位的跟踪位置作为预测位置，跟踪所述控制点在连续离散相位三维CT图像的位置包括： 

计算控制点在当前相位的新位置，根据Bhattacharyya系数对控制点在当前相位进行mean shift迭代，保存每个控制点在迭代过程中所有跟踪位置； 

计算预测位置和跟踪位置在当前体中局部邻域与控制点在参考体中局部邻域的均方误差； 

选择均方误差最小的跟踪位置作为当前控制点在当前相位的最终跟踪位置； 

在各连续离散相位重复上述过程，确定所述控制点在各连续离散相位的最终跟踪位置。 

5.根据权利要求4所述的方法，其特征在于，所述计算控制点在当前相位的新位置为： 

其中，y₀是当前相位肺部三维CT图像的控制点，

为控制点y₀的新位置，{x_i}_i=1...n是以y₀为中心的当前相位三维CT图像中局部邻域的体素坐标，x_i=(x,y,z)；h是局部邻域的直径大小，是点间欧几里德距离；函数

将位置x_i处的体素强度f(x_i)进行量化，

是下取整函数，δ为Kronecher delta函数，y^*是参考相位肺部三维CT图像的控制点，{x_i}_i=1...n是以y^*为中心的参考相位三维CT图像中局部邻域的体素坐标，C_h是使满 足的归一化常量，C_p是使

满足的归一化常量，n是三维CT图像局部邻域内体素的个数。 

6.如权利要求1所述的方法，其特征在于，所述对所述控制点在连续离散相位的跟踪轨迹进行最小二乘曲线拟合，确定所述控制点在任意时刻的运动位置为： 

控制点在离散时刻的位置为(x_s,y_s,z_s)，s=0,1,...,m，m为相位个数，分别沿x,y,z三个方向拟合二次多项式：

其中a₀,a₁,a₂，b₀,b₁,b₂和c₀,c₁,c₂分别为拟合二次多项式系数，p_x(t)，p_y(t)，p_z(t)分别是x,y,z三个方向拟合的二次多项式； 

通过构造线性方程组，计算得到拟合二次多项式系数a₀,a₁,a₂，b₀,b₁,b₂和c₀,c₁,c₂； 

根据基于最小二乘法的多项式曲线拟合方法，将连续时刻的离散控制点拟合成平滑的曲线，获取到控制点集(p_x(t),p_y(t),p_z(t)),t∈[0,m]。 

7.如权利要求1所述的方法，其特征在于，所述根据所述控制点在任意连续时刻的运动位置构建当前时刻肺部的三维形变模型，得到参考相位肺部的任意体素在任意时间的运动位置估计为： 

假设参考相位的肺部三维CT图像中的控制点坐标为(x_i(0),y_i(0),z_i(0))，i=1,2,L,N，对任意时刻t∈[0,m]，对应的控制点坐标为(x_i(t),y_i(t),z_i(t))，p_x(t)=x_i(t)，p_y(t)=y_i(t)，p_z(t)=z_i(t)，参考相位的三维CT图像在时刻t的形变函数为 

其中，U(r)为薄板样条的基函数U(r_i)=r_i，径向距离r_i定义为： 

(x,y,z)为参考相位中任意体素的坐标； 

形变参数a_t=(a_1,t,a_x,t,a_y,t,a_z,t)和w_t=(w_1,t,w_2,t,...,w_N,t)由

得到；其 中，

K是N×N矩阵：k_ij=U(r_ij)， 

P是N×4矩阵： O为4×4的零矩阵。 

8.一种肺部运动模型估计系统，其特征在于，所述系统包括：检测单元、跟踪单元、拟合单元以及构建单元；其中， 

所述检测单元，用于提取参考相位的肺部三维CT图像中控制点； 

所述跟踪单元，用于将所述控制点在上一时刻的跟踪位置作为预测位置，跟踪所述控制点在连续离散相位三维CT图像的位置； 

所述拟合单元，用于对所述控制点在连续离散相位的跟踪轨迹进行最小二乘曲线拟合，确定所述控制点在任意连续时刻的运动位置； 

所述构建单元，用于根据所述控制点在任意连续时刻的运动位置构建当前时刻肺部的三维形变模型，得到参考相位肺部的任意体素在任意时间的运动位置估计。 

9.如权利要求8所述的系统，其特征在于，所述检测单元，具体用于计算所述肺部三维CT图像中任一体素的自相关矩阵C；利用R=Det(C)-k·Trace(C)计算所述体素的Harris角点，当所述体素的Harris角点大于预先设定的阈值T时，则提取所述体素为控制点；其中，Det(C)是所述自相关矩阵C的行列式值，Trace(C)是所述自相关矩阵C的迹，k为预先设置的调整系数。 

10.如权利要求9所述的系统，其特征在于，所述检测单元，具体计算所述肺部三维CT图像中任一体素的自相关矩阵C包括： 

根据

计算所肺部三维CT图像中体素的自相 关矩阵C，其中，I_x,I_y,I_z分别是三维CT图像在(x,y,z)处沿x方向、y方向和z方向的一阶导；以所述体素

为中心和原点，利用高斯平滑函数对局部邻域内的各体素的自相关矩阵C的各项分别进行平滑，(x,y,z)是高斯平滑窗口内各体素坐标，

是高斯函数的方差。 

11.如权利要求8所述的系统，其特征在于，所述跟踪单元，具体用于计算控制点在当前相位的新位置，根据Bhattacharyya系数对控制点在当前相位进行mean shift迭代，保存每个控制点在迭代过程中所有跟踪位置；计算预测位置和跟踪位置在当前体中局部邻域与控制点在参考体中局部邻域的均方误差；选择均方误差最小的跟踪位置作为当前控制点在当前相位的最终跟踪位置；在各连续离散相位重复上述过程，确定所述控制点在各连续离散相位的最终跟踪位置。 

12.根据权利要求11所述的系统，其特征在于，所述跟踪单元，具体用于计算控制点在当前相位的新位置为： 

其中，y₀是当前相位肺部三维CT图像的控制点，{x_i}_i=1...n是以y₀为中心的当前相位三维CT图像中局部邻域的体素坐标，x_i=(x,y,z)；h是局部邻域的直径大小，

||||是点间欧几里德距离；函数

将位置x_i处的体素强度f(x_i)进行量化，

是下取整函数，δ为Kronecher delta函数，y^*是参考相位肺部三维CT图像的控制点，其中{x_i}_i=1...n是以y^*为中心的参考相 位三维CT图像中局部邻域的体素坐标，C_h是使

满足的归一化常量，C_p是使

满足的归一化常量，n是三维CT图像局部邻域内体素的个数。 

13.如权利要求8所述的系统，其特征在于，所述拟合单元，具体用于对所述控制点在连续离散相位的跟踪轨迹进行最小二乘曲线拟合，确定所述控制点在任意时刻的运动位置为： 

控制点在离散时刻的位置为(x_s,y_s,z_s)，s=0,1,...,m，m为相位个数，分别沿x,y,z三个方向拟合二次多项式：

其中a₀,a₁,a₂，b₀,b₁,b₂和c₀,c₁,c₂分别为拟合二次多项式系数，p_x(t)，p_y(t)，p_z(t)分别是x,y,z三个方向拟合的二次多项式； 

通过构造线性方程组，计算得到拟合二次多项式系数a₀,a₁,a₂，b₀,b₁,b₂和c₀,c₁,c₂； 

根据基于最小二乘法的多项式曲线拟合方法，将连续时刻的离散控制点拟合成平滑的曲线，获取到控制点集(p_x(t),p_y(t),p_z(t)),t∈[0,m]。 

14.如权利要求8所述的系统，其特征在于，所述构建单元，具体用于根据所述控制点在任意连续时刻的运动位置构建当前时刻肺部的三维形变模型，得到参考相位肺部的任意体素在任意时间的运动位置估计为： 

假设参考相位的肺部三维CT图像中的控制点坐标为(x_i(0),y_i(0),z_i(0))，i=1,2,L,N，对任意时刻t∈[0,m]，对应的控制点坐标为(x_i(t),y_i(t),z_i(t))，p_x(t)=x_i(t)，p_y(t)=y_i(t)，p_z(t)=z_i(t)，参考相位的三维CT图像在时刻t的形变函数为 

其中，U(r)为薄板样条的基函数U(r_i)=r_i，径向距离r_i定义为：

(x,y,z)为参考相位中任意体素的坐标； 

形变参数a_t=(a_1,t,a_x,t,a_y,t,a_z,t)和w_t=(w_1,t,w_2,t,...,w_N,t)由得到；其中，

K是N×N矩阵：k_ij=U(r_ij)， 

P是N×4矩阵： O为4×4的零矩阵。 

Images