CN104517301A - 多参数模型指导的迭代提取血管造影图像运动参数的方法 - Google Patents

Abstract

本发明公开了一种多参数模型指导的迭代提取血管造影图像运动参数的方法，包括：从造影图序列的医学图像中自动选取I个血管结构特征点，并分别对这些特征点在整个造影图序列中进行自动跟踪，以得到每个特征点的跟踪序列，利用离散傅里叶变换对每个特征点的跟踪序列进行处理，以生成离散傅里叶变换结果S_i(k)，初始化迭代参数j＝0，并获取离散傅里叶变换结果S_i(k)中各频点的幅度和频率范围，在各频点的幅度和频率范围中对该频点的跟踪序列进行傅里叶变换，以得到傅里叶变换结果，将获得的傅里叶变换结果进行逆傅里叶变换，并获取该频点的估计最小均方误差。本发明能够解决现有方法中存在的不能自动提取平移运动、以及不能准确分离呼吸和心脏等运动的技术问题。

Description

多参数模型指导的迭代提取血管造影图像运动参数的方法

技术领域

本发明属于数字信号处理与医学成像的交叉领域，更具体地，涉及一种多参数模型指导的迭代提取血管造影图像运动参数的方法。

背景技术

在X射线造影***中，由于呼吸运动的影响，冠状动脉血管在造影面上的图像会发生二维运动。在造影过程中，医生为了能够使造影序列覆盖整个冠脉***，可能会移动导管床，会导致一个成像序列的不同帧之间整体存在二维平移。因此，冠脉造影图像一方面记录有心脏自身运动在二维平面上的投影，同时也叠加有人体的呼吸运动引起的冠状动脉在造影面上的二维运动和病人的二维平移运动。那么，要从动态二维血管造影图重建三维心血管树，则需将人体的心脏运动，呼吸运动和平移运动等分离开来。

在造影过程中病人平移将直接影响运动分离的结果，而将平移运动分离实际上是进行帧间图像配准，现有的医学图像配准分外部和内部。外部配准方法是在造影前设置一些标记，在造影过程中跟踪它们，但是这种标记法通常是具有侵入性的；内部配准方法分为标记法和分割法等，标记法是在图中选取一些解剖结构点来进行配准。然而并不是每幅图都具有这样的解剖结构点，而且也需要熟悉人体解剖结构。分割法通过分割出来的解剖结构线来配准，既适用于刚体模型也适用于形变模型，但是它只能提取出解剖结构线的运动，不能单独分离出整体的刚性平移运动。然而造影图中表观的血管运动不仅包含了病人平移，还有血管自身的生理运动和呼吸运动。

也有文献中报导了将心脏运动和呼吸运动分离的方法，其中一种方法是预先在体内的一些器官如心脏、膈肌等设置标记点，然后对它们进行跟踪，把跟踪这些结构特征点所得到的运动曲线近似为呼吸运动。无论是直接手动跟踪造影图中的非心脏结构特征点，还是在造影同时就记录这些点的运动，都是有缺陷的。前者主要在于它的适用性很差；后者的实现需要大量实验控制，对一般的临床应用不合适。

另外一种方法是在双臂x射线造影条件下获得呼吸和心脏运动参数，该方法必须先重构出冠脉的三维运动序列而后分解得到呼吸和心脏运动，处理过程十分复杂。在单臂三维重建中，参与重建的两幅不同角度的造影图像呼吸运动时相不同，所以文献的方法不适用于单臂造影***。特别地，实际临床应用中，由于造影剂对人体的毒害性，造影时间通常是很短的，我们得到的造影图像的长度有限，从而导致造影图像中的结构特征点的运动是有限长的，这种有限长的运动中包含有周期分量、非周期分量和周期不完整的分量，因此，常规的离散傅里叶变换不能分离。

在申请号为201310750294.5的发明申请“单臂X射线血管造影图像多运动参数分解估计方法”中，虽然也是使用的傅里叶频域分离，但是它不能自动的将各运动参数通过循环迭代的方法自动提取，这样得到的各运动参数不够准确，此方法的鲁棒性不高。申请号为201310752751.4的发明申请“一种X射线血管造影图像中多运动参数的分解估计方法”中采用的是经验模态分解(Empirical Mode Decomposed，简称EMD)的方法估计各参数，同样得到的各运动参数不够准确，鲁棒性不高。

本发明提出了一种自动不断循环优化迭代从X射线造影序列图像中提取平移运动，呼吸运动和心脏运动等多运动参数的方法，通过在冠脉血管上选取一些特征点并在序列中进行跟踪得到它们随时间的运动曲线，然后在多运动参数模型约束下，以离散傅里叶正-反变换和估计的重建信号与原始信号在频域中各频点的局部均方误差最小为准则，通过循环迭代使原始信号与重建信号的差值信号均方最小的优化方法对这些运动分量进行优化处理，得到二维心脏、震颤、呼吸和平移等最优运动曲线。具有很好的临床适用性。

发明内容

针对现有技术的以上缺陷或改进需求，本发明提供了一种多参数模型指导的迭代提取血管造影图像运动参数的方法，其目的在于，解决现有方法中存在的不能自动提取平移运动、以及不能准确分离呼吸和心脏等运动的技术问题。

为实现上述目的，按照本发明的一个方面，提供了一种多参数模型指导的迭代提取血管造影图像运动参数的方法，包括以下步骤：

(1)从造影图序列的医学图像中自动选取I个血管结构特征点，并分别对这些特征点在整个造影图序列中进行自动跟踪，以得到每个特征点的跟踪序列{s_i(n),i＝1,…,I}，其中n表示医学图像在造影图序列中的帧号：

(2)利用离散傅里叶变换对步骤(1)获得的每个特征点的跟踪序列{s_i(n),i＝1,…,I}进行处理，以生成离散傅里叶变换结果S_i(k)：

(3)初始化迭代参数j＝0，并获取步骤(2)中获得的离散傅里叶变换结果S_i(k)中各频点的幅度和频率范围；

(4)在各频点的幅度和频率范围中对该频点的跟踪序列进行傅里叶变换，以得到傅里叶变换结果；

(5)将步骤(4)获得的傅里叶变换结果进行逆傅里叶变换，并获取该频点的估计最小均方误差；

(6)判断估计最小均方误差是否大于预设的阈值，如果大于，则转入步骤(7)，否则过程结束；

(7)利用多参数迭代优化算法对各频点的频谱进行处理，以得到第j+1次迭代后的时域信号；

(8)通过平移模型对剩余信号进行处理，以得到第j+1次迭代后的平移信号；

(9)将第j+1次迭代后的时域信号和第j+1次迭代后的平移信号叠加后得到第j+1次迭代后的估计混合信号，并利用以下公式计算第j+1次迭代后的最小均方误差：

(10)判断第j+1次迭代后的最小均方误差是否大于步骤(6)中的阈值，如果大于则返回步骤(7)，否则过程结束。

优选地，步骤(1)是采用以下公式：

s_i(n)＝L(n)+r_i(n)+c_i(n)+h_i(n)+t_i(n),i∈[1,I]

其中，L(n)表示平移运动；r_i(n)表示第i个特征点的呼吸运动；c_i(n)表示第i个特征点的心脏运动；h_i(n)表示第i个特征点的震颤运动；t_i(n)表示第i个特征点的其他运动。

优选地，步骤(2)是采用以下公式：

S_i(k)＝L(k)+R_i(k)+C_i(k)+H_i(k)

其中k表示频点，L(k)、C(k)、R(k)和H(k)分别表示L(n)、c(n)、r(n)和h(n)的各次谐波系数。

优选地，步骤(5)中该频点的估计最小均方误差是采用以下公式计算：

${\widehat{\mathrm{\ϵ}}}_i^j=\min \left(\frac{1}{N}\sqrt{\underset{n}{\mathrm{\Σ}}{(s_i\left(n\right)-{\hat{s}}_i^j\left(n\right))}^2}\right)$

其中 ${\hat{s}}_i^j\left(n\right)=L^j\left(n\right)+r_i^j\left(n\right)+c_i^j\left(n\right)+h_i^j\left(n\right).$

优选地，步骤(7)包括以下子步骤：

(7-1)保持L^j(k),不变，采用以下公式计算其在频率范围中的频点k_ic附近的值L^j(k_ic)、

$X_p\left(k\right)=\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}{x}_p\left(n\right)\·e^{-j\left(\frac{2\mathrm{\π}}{N}\right)\mathrm{nk}},$

然后利用公式 $C_i^{j+1}\left(k_{\mathrm{ic}}\right)=C_i^0\left(k_{\mathrm{ic}}\right)-L^j\left(k_{\mathrm{ic}}\right)-R_i^j\left(k_{\mathrm{ic}}\right)-H_i^j\left(k_{\mathrm{ic}}\right)$ 计算得到第j+1次迭代的心脏信号频谱，在该频点的频率范围内求离散傅里叶变换后，利用以下公式得到第j+1次的最优心脏时域信号

${\hat{e}}_i^j=\min \left(\underset{k=\mathrm{\ω}-M}{\overset{\mathrm{\ω}+M}{\mathrm{\Σ}}}{(S_i\left(k\right)-{\hat{S}}_i^j\left(k\right))}^2\right)$

其中，ω为信号经离散傅里叶变换后所得频点；M表示窗口的大小，一般设置为3； ${\hat{S}}_i^j\left(k\right)=L^j\left(k\right)+R_i^j\left(k\right)+C_i^{j+1}\left(k\right)+H_i^j\left(k\right);$

(7-2)保持L^j(k),不变，采用以下公式计算其在频率范围中的频点k_ih附近的值L^j(k_ih)、

$X_p\left(k\right)=\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}{x}_p\left(n\right)\·e^{-j\left(\frac{2\mathrm{\π}}{N}\right)\mathrm{nk}}$

然后利用公式 $H_i^{j+1}\left(k_{\mathrm{ih}}\right)=H_i^0\left(k_{\mathrm{ih}}\right)-L^j\left(k_{\mathrm{ih}}\right)-R_i^j\left(k_{\mathrm{ih}}\right)-C_i^j\left(k_{\mathrm{ih}}\right)$ 计算得到第j+1次迭代的高频信号频谱，在该频点的频率范围内求离散傅里叶变换后，利用以下公式得到第j+1次的最优心脏时域信号

${\hat{e}}_i^j=\min \left(\underset{k=\mathrm{\ω}-M}{\overset{\mathrm{\ω}+M}{\mathrm{\Σ}}}{(S_i\left(k\right)-{\hat{S}}_i^j\left(k\right))}^2\right)$

其中 ${\hat{S}}_i^j\left(k\right)=L^j\left(k\right)+R_i^j\left(k\right)+C_i^{j+1}\left(k\right)+H_i^{j+1}\left(k\right);$

(7-3)保持L^j(k),不变，采用以下公式计算其在频率范围中的频点k_ir附近的值L^j(k_ir)、

$X_p\left(k\right)=\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}{x}_p\left(n\right)\·e^{-j\left(\frac{2\mathrm{\π}}{N}\right)\mathrm{nk}}$

然后利用公式 $R_i^{j+1}\left(k_{\mathrm{ir}}\right)=R_i^0\left(k_{\mathrm{ir}}\right)-L^j\left(k_{\mathrm{ir}}\right)-C_i^{j+1}\left(k_{\mathrm{ir}}\right)-H_i^{j+1}\left(k_{\mathrm{ir}}\right)$ 计算得到第j+1次迭代的呼吸信号频谱，在该频点的频率范围内求离散傅里叶变换后，利用以下公式得到第j+1次的最优心脏时域信号

${\hat{e}}_i^j=\min \left(\underset{k=\mathrm{\ω}-M}{\overset{\mathrm{\ω}+M}{\mathrm{\Σ}}}{(S_i\left(k\right)-{\hat{S}}_i^j\left(k\right))}^2\right)$

其中 ${\hat{S}}_i^j\left(k\right)=L^j\left(k\right)+R_i^{j+1}\left(k\right)+C_i^{j+1}\left(k\right)+H_i^{j+1}\left(k\right).$

优选地，步骤(9)中最小均方误差是采用以下公式计算：

${\widehat{\mathrm{\ϵ}}}_i^{j+1}=\min \left(\frac{1}{N}\sqrt{\underset{n}{\mathrm{\Σ}}{(s_i\left(n\right)-{\hat{s}}_i^{j+1}\left(n\right))}^2}\right)$

按照本发明的另一方面，提供了一种多参数模型指导的迭代提取血管造影图像运动参数的***，包括：

第一模块，用于从造影图序列的医学图像中自动选取I个血管结构特征点，并分别对这些特征点在整个造影图序列中进行自动跟踪，以得到每个特征点的跟踪序列{s_i(n),i＝1,…,I}，其中n表示医学图像在造影图序列中的帧号：

第二模块，用于利用离散傅里叶变换对第一模块获得的每个特征点的跟踪序列{s_i(n),i＝1,…,I}进行处理，以生成离散傅里叶变换结果S_i(k)：

第三模块，用于初始化迭代参数j＝0，并获取第二模块中获得的离散傅里叶变换结果S_i(k)中各频点的幅度和频率范围；

第四模块，用于在各频点的幅度和频率范围中对该频点的跟踪序列进行傅里叶变换，以得到傅里叶变换结果；

第五模块，用于将第四模块获得的傅里叶变换结果进行逆傅里叶变换，并获取该频点的估计最小均方误差；

第六模块，用于判断估计最小均方误差是否大于预设的阈值，如果大于，则转入第七模块，否则过程结束；

第七模块，用于利用多参数迭代优化算法对各频点的频谱进行处理，以得到第j+1次迭代后的时域信号；

第八模块，用于通过平移模型对剩余信号进行处理，以得到第j+1次迭代后的平移信号；

第九模块，用于将第j+1次迭代后的时域信号和第j+1次迭代后的平移信号叠加后得到第j+1次迭代后的估计混合信号，并利用以下公式计算第j+1次迭代后的最小均方误差：

第十模块，用于判断第j+1次迭代后的最小均方误差是否大于第六模块中的阈值，如果大于则返回第七模块，否则过程结束。

总体而言，通过本发明所构思的以上技术方案与现有技术相比，能够取得下列有益效果：

1、由于采用了步骤(2)，本发明自动获取结构特征点得到运动曲线的方法克服了在心脏的表面设置标记进行跟踪获取特征点的运动。

2、采用步骤(7)的多次循环迭代寻优的方法既可以解决文献中不能自动提取平移运动的缺陷；又解决了其他周期运动(如呼吸运动、心脏运动等)不能分离的问题；

3、由于采用了步骤(7)的多次循环迭代寻优的方法解决了专利中所提取的各运动信号不准确且鲁棒性不好的问题。

附图说明

图1是本发明多参数模型指导的迭代提取血管造影图像运动参数的方法的流程图；

图2为获得的一个角度的原始造影图像；

图3为图2中血管骨架图像选取的特征点；

图4(a)和(b)为图2的特征点A1-A4的X轴和Y轴的原始运动曲线；

图5(a)和(b)为图2的第四个特征点的X轴分量频谱和第四个特征点的Y轴分量频谱；

图6(a)和(b)为图2的各特征点分离出的X轴和Y轴的心脏运动曲线；

图7(a)和(b)为图2的各特征点分离出的X轴和Y轴的呼吸运动曲线；

图8(a)和(b)为图2的各特征点分离出的X轴和Y轴的高频运动曲线；

图9(a)和(b)为图2的各特征点分离出的X轴和Y轴平移曲线与手动跟踪骨骼点得到的平移曲线的对比。

图10(a)和(b)为图2的各特征点分离出的X轴和Y轴分离的残余曲线。

图11为获得的另一个角度的原始造影图像；

图12为图11中血管骨架图像选取的特征点；

图13(a)和(b)为图11的特征点A1-A5的X轴和Y轴的原始运动曲线；

图14(a)和(b)为图11的第四个特征点的X轴分量频谱和第四个特征点的Y轴分量频谱；

图15(a)和(b)为图11的各特征点分离出的X轴和Y轴的心脏运动曲线；

图16(a)和(b)为图11的各特征点分离出的X轴和Y轴的心脏运动曲线；

图17(a)和(b)为图11的各特征点分离出的X轴和Y轴呼吸运动曲线。

图18(a)和(b)为图11的各特征点分离出的X轴和Y轴分离的残余曲线。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

以下首先就本发明的技术术语进行解释和说明：

最小均方误差：即估计信号与原始信号差值均方累加求和最小。本文多迭代参数模型利用离散傅里叶变换的性质，采取最小均方误差的方法不断在时域和频域迭代使得残差最小的原则求得各分离信号的最优。

为解决现有技术中存在的问题，本发明提出了结构特征点的多运动参数模型。在冠脉血管造影图上自动选取血管结构特征点，在造影图序列中对其进行自动跟踪，得到它们随时间变化的位移曲线。通过正-反傅里叶变换，重建信号与原始信号在频域中各频点的局部均方误差最小初值估计为基础，结合平移模型，采用循环迭代使原始信号与重建信号的差值信号最小的全局优化得到最终估计的平移运动曲线、呼吸运动曲线以及心脏运动参数。该方法具有广泛的适用性和灵活性，几乎可适用于所有造影序列图(当然需满足有较清晰的血管分布和包含两个或两个以上心动周期，这点不难做到)，也包括双臂X射线造影图像。

如图1所示，本发明多参数模型指导的迭代提取血管造影图像运动参数的方法包括以下步骤：

(1)从造影图序列的医学图像中自动选取I个血管结构特征点(其中I为自然数)，并分别对这些特征点在整个造影图序列中进行自动跟踪，以得到每个特征点的跟踪序列{s_i(n),i＝1,…,I}，其中n表示医学图像在造影图序列中的帧号：

s_i(n)＝L(n)+r_i(n)+c_i(n)+h_i(n)+t_i(n),i∈[1,I]

其中，L(n)表示平移运动；r_i(n)表示第i个特征点的呼吸运动；c_i(n)表示第i个特征点的心脏运动；h_i(n)表示第i个特征点的震颤运动；t_i(n)表示第i个特征点的其他运动；

(2)利用离散傅里叶变换对步骤(1)获得的每个特征点的跟踪序列{s_i(n),i＝1,…,I}进行处理，以生成离散傅里叶变换结果S_i(k)：

S_i(k)＝L(k)+R_i(k)+C_i(k)+H_i(k)

其中k表示频点，L(k)、C(k)、R(k)和H(k)分别表示L(n)、c(n)、r(n)和h(n)的各次谐波系数；

(3)初始化迭代参数j＝0，并获取步骤(2)中获得的离散傅里叶变换结果S_i(k)中各频点的幅度和频率范围；

(4)在各频点的幅度和频率范围中对该频点的跟踪序列进行傅里叶变换，以得到傅里叶变换结果L^j(k),

(5)将步骤(4)获得的L^j(k),进行逆傅里叶变换，然后利用以下公式获取该频点的估计最小均方误差

${\widehat{\mathrm{\ϵ}}}_i^j=\min \left(\frac{1}{N}\sqrt{\underset{n}{\mathrm{\Σ}}{(s_i\left(n\right)-{\hat{s}}_i^j\left(n\right))}^2}\right)$

其中 ${\hat{s}}_i^j\left(n\right)=L^j\left(n\right)+r_i^j\left(n\right)+c_i^j\left(n\right)+h_i^j\left(n\right).$

(6)判断估计最小均方误差是否大于预设的阈值(其为不大于2的正整数)，如果大于，则转入步骤(7)，否则过程结束；

(7)利用多参数迭代优化算法对各频点的频谱L^j(k),进行处理，以得到第j+1次迭代后的时域信号等；

本步骤包括以下子步骤：

(7-1)保持L^j(k),不变，采用以下公式计算其在频率范围中的频点k_ic附近的值L^j(k_ic)、

$X_p\left(k\right)=\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}{x}_p\left(n\right)\·e^{-j\left(\frac{2\mathrm{\π}}{N}\right)\mathrm{nk}},$

然后利用公式 $C_i^{j+1}\left(k_{\mathrm{ic}}\right)=C_i^0\left(k_{\mathrm{ic}}\right)-L^j\left(k_{\mathrm{ic}}\right)-R_i^j\left(k_{\mathrm{ic}}\right)-H_i^j\left(k_{\mathrm{ic}}\right)$ 计算得到第j+1次迭代的心脏信号频谱，在该频点的频率范围内求离散傅里叶变换后，利用以下公式得到第j+1次的最优心脏时域信号

${\hat{e}}_i^j=\min \left(\underset{k=\mathrm{\ω}-M}{\overset{\mathrm{\ω}+M}{\mathrm{\Σ}}}{(S_i\left(k\right)-{\hat{S}}_i^j\left(k\right))}^2\right)$

其中，ω为信号经离散傅里叶变换后所得频点；M表示窗口的大小，一般设置为3；

(7-2)保持L^j(k),不变，采用以下公式计算其在频率范围中的频点k_ih附近的值L^j(k_ih)、

$X_p\left(k\right)=\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}{x}_p\left(n\right)\·e^{-j\left(\frac{2\mathrm{\π}}{N}\right)\mathrm{nk}}$

然后利用公式计算得到第j+1次迭代的高频信号频谱，在该频点的频率范围内求离散傅里叶变换后，利用以下公式得到第j+1次的最优心脏时域信号

${\hat{e}}_i^j=\min \left(\underset{k=\mathrm{\ω}-M}{\overset{\mathrm{\ω}+M}{\mathrm{\Σ}}}{(S_i\left(k\right)-{\hat{S}}_i^j\left(k\right))}^2\right)$

其中 ${\hat{S}}_i^j\left(k\right)=L^j\left(k\right)+R_i^j\left(k\right)+C_i^{j+1}\left(k\right)+H_i^{j+1}\left(k\right).$

(7-3)保持L^j(k),不变，采用以下公式计算其在频率范围中的频点k_ir附近的值L^j(k_ir)、

$X_p\left(k\right)=\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}{x}_p\left(n\right)\·e^{-j\left(\frac{2\mathrm{\π}}{N}\right)\mathrm{nk}}$

然后利用公式 $R_i^{j+1}\left(k_{\mathrm{ir}}\right)=R_i^0\left(k_{\mathrm{ir}}\right)-L^j\left(k_{\mathrm{ir}}\right)-C_i^{j+1}\left(k_{\mathrm{ir}}\right)-H_i^{j+1}\left(k_{\mathrm{ir}}\right)$ 计算得到第j+1次迭代的呼吸信号频谱，在该频点的频率范围内求离散傅里叶变换后，利用以下公式得到第j+1次的最优心脏时域信号

${\hat{e}}_i^j=\min \left(\underset{k=\mathrm{\ω}-M}{\overset{\mathrm{\ω}+M}{\mathrm{\Σ}}}{(S_i\left(k\right)-{\hat{S}}_i^j\left(k\right))}^2\right)$

其中 ${\hat{S}}_i^j\left(k\right)=L^j\left(k\right)+R_i^{j+1}\left(k\right)+C_i^{j+1}\left(k\right)+H_i^{j+1}\left(k\right).;$

(8)通过平移模型对剩余信号 $v_i^{j+1}\left(n\right)=s_i\left(n\right)-c_i^{j+1}\left(n\right)-h_i^{j+1}\left(n\right)-r_i^{j+1}\left(n\right)$ 进行处理，以得到第j+1次迭代后的平移信号L^j+1(n)；

(9)将和L^j+1(n)叠加后得到第j+1次迭代后的估计混合信号然后利用以下公式计算第j+1次迭代后的最小均方误差 ${\mathrm{\ϵ}}_i^{j+1}:$

${\widehat{\mathrm{\ϵ}}}_i^{j+1}=\min \left(\frac{1}{N}\sqrt{\underset{n}{\mathrm{\Σ}}{(s_i\left(n\right)-{\hat{s}}_i^{j+1}\left(n\right))}^2}\right)$

(10)判断第j+1次迭代后的最小均方误差是否大于步骤(6)中的阈值，如果大于则返回步骤(7)，否则过程结束。

图2是某病人1在角度(50.8°，30.2°)下的造影图像以及自动提取该图像序列各结构分叉特征点(图3)的运动分量实验。其中，该病人的心脏运动周期为10帧，造影图序列的采样率为12.5帧/s。由于造影时间较短，本实验中选取在造影序列图中停留时间较长的结构特征点A1～A4用于实验。

图11是某病人4在角度(30.8°，15.3°)下的造影图像以及自动提取该图像序列各结构分叉特征点的运动分量实验。其中，该病人患有心率不齐疾病，利用本文算法可以提取出该病人在心脏周期范围内的两种心脏运动周期分别为9帧和12帧，造影图序列的采样率为12.5帧/s。

从图4中的图(a)和图(b)以及图13中的图(a)和图(b)中可以看到，特征点的原始运动曲线受呼吸运动影响较大，比较凌乱，而且不同特征点的运动幅度不同，这是因为不同特征点位于心脏的不同部位，而心脏不同部位的运动情况是不一样的。图5的图(a)和图(b)以及图14的图(a)和图(b)表示混合信号经过离散傅里叶后的频谱图，从频谱图中可以看出，各频率点呈现的峰值很明显，峰值点的多少说明了该混合信号包含的各频率信号的种类，特别地，如果在0频率点处呈现的峰值很大，则该混合信号一定存在平移信号。图6的图(a)和图(b)、图15的图(a)图(b)以及图16的图(a)和图(b)为分离后的心脏运动曲线，从图中可以看出各特征点运动曲线显示了良好的周期性。图7的图(a)和图(b)以及图17的图(a)和图(b)表示分离后的呼吸运动曲线，从图中可以看出各特征点的呼吸运动曲线不仅具有良好的周期性，而且波峰所在位置基本一致，这是因为造影时间短，各特征点短时间的相位变化非常小。波峰的大小反映了各特征点在血管表面的分布，若特征点离肺部越近，则该特征点的波峰值越大。图8的图(a)和图(b)为分离后的震颤信号，说明了该病人在治疗的过程中出现了比较明显的抖动，抖动的幅度随着各特征点的分布而不一样。图9的图(a)和图(b)分别为各特征点分离的平移曲线与手动跟踪的平移曲线的对比，可以看出两者除了因为手动跟踪引起的误差外，基本一致。图10的图(a)和图(b)以及图18的图(a)和图(b)分别为各特征点分离各个运动的剩余曲线，其幅度都少于2个像素，可以看作是因为分割，骨架提取等算法误差引起的，运动分量基本完全分离。

不同于图2中的某病人1的是，图11中的某病人4没有出现震颤的高频成分，却出现了在心脏信号频率范围内的两种不同频率成分，这类信号的出现，为医生诊断病人病情提供了极大的参考价值。

总而言之，本发明方法考虑到在实际单臂X射线造影的条件下，并非在所有造影图中都能找到合适的特征点(像肋骨的交点等)，所以采用血管结构特征点的选取和傅立叶频域滤波相以及多变量寻优结合的途径来提取多个运动分量，比单纯的用手动跟踪来得到呼吸运动或者平移运动具有更广泛的适用性和灵活性，几乎可适用于所有造影序列图，同时此方法相较于直接在心脏附近组织设置标识点再通过相关成像手段进行跟踪的方法拥有更高安全性和可操作性，这是因为在体内组织添加的标记物一般是可侵入性的，会对人体自身产生或多或少的损害，并且其标记物的添加、成像、排除、呼吸运动提取整个过程都是繁杂的，为实际操作中带来不可避免的麻烦与误差。

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

Claims (7)

1.一种多参数模型指导的迭代提取血管造影图像运动参数的方法，其特征在于，包括以下步骤：

(1)从造影图序列的医学图像中自动选取I个血管结构特征点，并分别对这些特征点在整个造影图序列中进行自动跟踪，以得到每个特征点的跟踪序列{s_i(n),i＝1,…,I}，其中n表示医学图像在造影图序列中的帧号：

(2)利用离散傅里叶变换对步骤(1)获得的每个特征点的跟踪序列{s_i(n),i＝1,…,I}进行处理，以生成离散傅里叶变换结果S_i(k)：

(3)初始化迭代参数j＝0，并获取步骤(2)中获得的离散傅里叶变换结果S_i(k)中各频点的幅度和频率范围；

(4)在各频点的幅度和频率范围中对该频点的跟踪序列进行傅里叶变换，以得到傅里叶变换结果；

(5)将步骤(4)获得的傅里叶变换结果进行逆傅里叶变换，并获取该频点的估计最小均方误差；

(6)判断估计最小均方误差是否大于预设的阈值，如果大于，则转入步骤(7)，否则过程结束；

(7)利用多参数迭代优化算法对各频点的频谱进行处理，以得到第j+1次迭代后的时域信号；

(8)通过平移模型对剩余信号进行处理，以得到第j+1次迭代后的平移信号；

(9)将第j+1次迭代后的时域信号和第j+1次迭代后的平移信号叠加后得到第j+1次迭代后的估计混合信号，并利用以下公式计算第j+1次迭代后的最小均方误差：

(10)判断第j+1次迭代后的最小均方误差是否大于步骤(6)中的阈值，如果大于则返回步骤(7)，否则过程结束。

2.根据权利要求1所述的方法，其特征在于，步骤(1)是采用以下公式：

s_i(n)＝L(n)+r_i(n)+c_i(n)+h_i(n)+t_i(n),i∈[1,I]

其中，L(n)表示平移运动；r_i(n)表示第i个特征点的呼吸运动；c_i(n)表示第i个特征点的心脏运动；h_i(n)表示第i个特征点的震颤运动；t_i(n)表示第i个特征点的其他运动。

3.根据权利要求2所述的方法，其特征在于，步骤(2)是采用以下公式：

S_i(k)＝L(k)+R_i(k)+C_i(k)+H_i(k)

其中k表示频点，L(k)、C(k)、R(k)和H(k)分别表示L(n)、c(n)、r(n)和h(n)的各次谐波系数。

4.根据权利要求3所述的方法，其特征在于，步骤(5)中该频点的估计最小均方误差是采用以下公式计算：

其中

5.根据权利要求4所述的方法，其特征在于，步骤(7)包括以下子步骤：

(7-1)保持L^j(k),不变，采用以下公式计算其在频率范围中的频点k_ic附近的值L^j(k_ic)、

然后利用公式计算得到第j+1次迭代的心脏信号频谱，在该频点的频率范围内求离散傅里叶变换后，利用以下公式得到第j+1次的最优心脏时域信号

其中，ω为信号经离散傅里叶变换后所得频点；M表示窗口的大小，一般设置为3；

(7-2)保持L^j(k),不变，采用以下公式计算其在频率范围中的频点k_ih附近的值L^j(k_ih)、

然后利用公式计算得到第j+1次迭代的高频信号频谱，在该频点的频率范围内求离散傅里叶变换后，利用以下公式得到第j+1次的最优心脏时域信号

其中

(7-3)保持L^j(k),不变，采用以下公式计算其在频率范围中的频点k_ir附近的值L^j(k_ir)、

然后利用公式计算得到第j+1次迭代的呼吸信号频谱，在该频点的频率范围内求离散傅里叶变换后，利用以下公式得到第j+1次的最优心脏时域信号

其中

6.根据权利要求5所述的方法，其特征在于，步骤(9)中最小均方误差是采用以下公式计算：

。

7.一种多参数模型指导的迭代提取血管造影图像运动参数的***，其特征在于，包括：

第一模块，用于从造影图序列的医学图像中自动选取I个血管结构特征点，并分别对这些特征点在整个造影图序列中进行自动跟踪，以得到每个特征点的跟踪序列{s_i(n),i＝1,…,I}，其中n表示医学图像在造影图序列中的帧号：

第二模块，用于利用离散傅里叶变换对第一模块获得的每个特征点的跟踪序列{s_i(n),i＝1,…,I}进行处理，以生成离散傅里叶变换结果S_i(k)：

第三模块，用于初始化迭代参数j＝0，并获取第二模块中获得的离散傅里叶变换结果S_i(k)中各频点的幅度和频率范围；

第四模块，用于在各频点的幅度和频率范围中对该频点的跟踪序列进行傅里叶变换，以得到傅里叶变换结果；

第五模块，用于将第四模块获得的傅里叶变换结果进行逆傅里叶变换，并获取该频点的估计最小均方误差；

第六模块，用于判断估计最小均方误差是否大于预设的阈值，如果大于，则转入第七模块，否则过程结束；

第七模块，用于利用多参数迭代优化算法对各频点的频谱进行处理，以得到第j+1次迭代后的时域信号；

第八模块，用于通过平移模型对剩余信号进行处理，以得到第j+1次迭代后的平移信号；

第九模块，用于将第j+1次迭代后的时域信号和第j+1次迭代后的平移信号叠加后得到第j+1次迭代后的估计混合信号，并利用以下公式计算第j+1次迭代后的最小均方误差：

第十模块，用于判断第j+1次迭代后的最小均方误差是否大于第六模块 中的阈值，如果大于则返回第七模块，否则过程结束。