CN103886615A - 一种x射线血管造影图像中多运动参数的分离估计方法 - Google Patents

Abstract

本发明公开了一种X射线血管造影图像中多运动参数的分离估计方法;首先根据所述的造影图像序列确定心脏运动信号周期以及平移运动变化帧的序列，然后通过对所述的造影图像序列中的血管结构特征点跟踪得到点的运动序列。再通过多变量参数寻优以及傅里叶频域滤波对所述的运动序列进行处理，根据所述的平移运动变化帧的序列，心脏运动信号的周期，呼吸运动信号周期的范围和高频运动信号周期的范围，可以分离出最佳的平移运动曲线、心脏运动曲线、呼吸运动曲线和高频运动曲线。本发明由于考虑到X射线血管造影图像是反映了多种运动信号的综合，提出的多运动参数的分离估计，能够使分离出来的各个运动更加精确，为医生的诊断提供准确有效的帮助。

Description

一种X射线血管造影图像中多运动参数的分离估计方法

技术领域

本发明属于造影图像处理技术领域，更具体地，涉及一种X射线血管造影图像中多运动参数的分离估计方法。 

背景技术

胸廓有节律的扩大和缩小，从而完成吸气与呼气，这就是呼吸运动。呼吸运动会造成心脏在三维空间中整体的平移，在X射线造影系统中，由于呼吸运动的影响冠状动脉血管在造影面上会发生二维运动。而且在造影过程中，医生为了能够使造影序列覆盖整个冠脉系统，可能会移动导管床，导致一个成像序列的不同帧之间整体存在二维平移。 

因此，冠脉造影图像一方面记录有心脏自身运动在二维平面上的投影，同时也叠加有人体的呼吸运动引起的冠状动脉在造影面上的二维平移运动和病人的二维平移运动。那么，要得到更接近真实情况下的二维血管造影图，则需将人体的心脏运动，呼吸运动和平移运动等分离开来。 

分离平移运动一般方法是进行帧间图像配准。许多文献报道的医学图像配准通常是分外部和内部。外部配准方法是通过在造影前设置一些明显的标记，在造影过程中跟踪，但是这种标记法通常是具有侵入性的；内部配准方法是在造影完成后又分为标记法和分割法等。标记法是在图中选取一些解剖结构点来进行配准，然而并不是每幅图都具有这样的解剖结构点，而且也需要对人体解剖结构比较了解。而分割法可以通过分割出来的解剖结构线来配准，可以适用于刚体模型也可适用于形变模型，但是它只能提取出解剖结构线的运动，却不能单独分离出整体的刚性平移运动。然而造影图中提取出来的血管运动不仅包含了病人平移运动，还有血管自身的生理运动和呼吸运动。虽然有些造影图中包含了不会有生理运动的脊椎等骨 骼可以作为内部标记点来通过内部标记法配准，但不是所有的造影图都存在这样的特征可以进行配准。当造影图不存在这样的特征时，平移运动是很难被提取出来的。而且这样的平移运动的提取在医学图像领域上少有提及。 

文献报导在提取人体呼吸运动的时候都要通过预先设置标记点，然后对它们进行序列跟踪。根据呼吸运动的特点，人在进行呼吸的时候体内的一些器官如心脏，膈肌等会随着肺的运动一起同步进行三维空间运动。所以跟踪不会随着心脏一起运动的结构特征点得到其运动曲线，可近似认作呼吸运动。其中在X射线造影中用的比较多的是，直接手动跟踪造影图中的非心脏结构特征点或者在造影同时就记录这些点的运动。很显然，这两种方案都是有缺陷的。前者主要在于它的适用性很差，因为我们并不能保证每一帧造影图中都存在符合这种条件的标记点，而且找这种点也需要对人体解剖结构比较了解。所以，在造影图中不存在以上特征点时，呼吸运动是很难被提取出来的。而后者的实现则需要大量实验控制，对一般的临床应用不合适。此外，还有一种方法是在双臂X射线造影条件下实现的，其分离心脏运动与呼吸运动的思想是，取同一时刻不同投影角度的两幅造影图，对其中相对应的冠脉血管进行三维重建，获得该时刻的血管三维空间分布，那么，对一个呼吸周期中的所有造影图对进行匹配和重建后，得到的则是一组三维结构序列，而它们间的空间位移矢量便是呼吸运动。相对来说，通过该方法能比较可靠的呼吸运动结果。但是由于双臂x射线造影条件的约束，不能广泛的应用在医学实践中。 

在专利一种从X射线造影图中提取呼吸运动参数的方法200910273528.5中，虽然也是使用的傅里叶频域分离，但是它只能通过提取心脏运动后再将剩余曲线视为呼吸运动，且不能分离出平移运动和高频成分，这样得到的呼吸运动不够准确，此方法的鲁棒性不高。 

发明内容

针对现有技术的以上缺陷或改进需求，本发明提供了一种X射线血管造影图像中多运动参数的分离估计方法，其目的在于基于频域滤波和多变量寻优，能够精确分离各个运动参数，简单快速、无伤害。 

本发明提供了一种X射线血管造影图像中多运动参数的分离估计方法，包括下述步骤： 

（1）对血管结构特征点在X射线血管造影图像序列中进行跟踪，得到点的跟踪曲线s_i(n)，i=1,…，I，n=1,...,N；I为特征点个数，N为X射线血管造影图像序列中图像的总帧数； 

（2）根据平移运动变化帧的序列{N_t,t=1,....,T}和平移运动的斜率角序列{α_t,t=1,....,T-1}得到模拟平移曲线d_α(n)；其中T为运动方向变化的次数； 

（3）根据所述X射线血管造影图像序列确定心脏运动信号周期N_c；根据所述跟踪曲线s_i(n)和所述模拟平移曲线d_α(n)得到去除了平移运动后的多运动的综合运动曲线根据心脏运动周期N_c对所述综合运动曲线

进行傅立叶频域滤波处理，得到心脏运动曲线

（4）根据所述综合运动曲线

和所述心脏运动曲线

得到不包含平移运动信号和心脏运动信号的剩余曲线

在呼吸运动信号周期[3N_c,10N_c]的范围内，根据每一个呼吸周期对所述剩余曲线进行傅立叶频域滤波，得到相应的呼吸运动曲线；以最接近所述剩余曲线

的拟合曲线

为最优准则，得到相对于当前模拟平移曲线的最佳呼吸运动曲线

和最佳呼吸运动信号周期

（5）根据所述剩余曲线和所述拟合曲线

得到的曲线 

判断所述曲线

的振幅是否小于3个像素，若是，则没有高频运动信号；若否，则在高频运动信号周期[1/7N_c,5/7N_c]的范围内，根据每一个高频运动周期对所述剩余曲线

进行傅立叶频域滤波，得到 相应的高频运动曲线；以最接近所述曲线

为最优准则，得到相对于当前模拟平移曲线的最佳高频运动曲线和最佳高频运动信号周期

（6）根据模拟平移曲线d_α(n)、呼吸运动曲线

心脏运动曲线

和高频运动曲线获得综合运动估计曲线  ${\hat{s}}_{\mathrm{i\α}}\left(n\right)=d_{\mathrm{\α}}\left(n\right)+{\hat{r}}_{\mathrm{i\α}}\left(n\right)+{\hat{c}}_{\mathrm{i\α}}\left(n\right)+{\hat{h}}_{\mathrm{i\α}}\left(n\right);$ 以所述综合运动估计曲线

与所述跟踪曲线s_i(n)最接近为最优准则，得到最佳的平移运动曲线，心脏运动曲线，呼吸运动曲线和高频运动曲线。 

总体而言，通过本发明所构思的以上技术方案与现有技术相比，由于考虑到X射线血管造影图像是反映了多种运动信号如平移运动信号、心脏运动信号、呼吸运动信号和高频运动信号等综合，提出了多运动参数的分离估计，能够使分离出来的各个运动更加精确，为医生的诊断提供准确有效的帮助。 

附图说明

图1是本发明实施例提供的X射线血管造影图像中多运动参数的分离估计方法的实现流程框图； 

图2（a）是傅立叶变换提取中造影图及骨骼点的标记,对应的投影角度为（-26.5°，-20.9°）； 

图2（b）是傅立叶变换提取中造影图及骨骼点的标记,对应的投影角度为（42.3°，26.8°）； 

图3（a）是图2（a）的骨架图及特征点标记。 

图3（b）是图2（b）的骨架图及特征点标记。 

图4（a）是对图3（a）中所有标记点在造影图序列中跟踪得到的原始曲线，其中，实线为各个特征点的横向坐标（X轴坐标）的变化，虚线为各个特征点的纵向坐标（Y轴坐标）的变化。 

图4（b）是对图3（b）中所有标记点在造影图序列中跟踪得到的原始 曲线，其中，实线为各个特征点的横向坐标（X轴坐标）的变化，虚线为各个特征点的纵向坐标（Y轴坐标）的变化。 

图4（c）是对图3（a）中原始曲线分解得到的心脏运动曲线，其中，实线为各个特征点的横向坐标（X轴坐标）的变化，虚线为各个特征点的纵向坐标（Y轴坐标）的变化。 

图4（d）是对图3（b）中原始曲线分解得到的心脏运动曲线，其中，实线为各个特征点的横向坐标（X轴坐标）的变化，虚线为各个特征点的纵向坐标（Y轴坐标）的变化。 

图4（e）是对图3（a）中原始曲线分解得到的呼吸运动信号曲线，其中，实线为各个特征点的横向坐标（X轴坐标）的变化，虚线为各个特征点的纵向坐标（Y轴坐标）的变化。 

图4（f）是对图3（b）中原始曲线分解得到的呼吸运动信号曲线，其中，实线为各个特征点的横向坐标（X轴坐标）的变化，虚线为各个特征点的纵向坐标（Y轴坐标）的变化。 

图4（g）是对图3（a）中原始曲线分解得到的高频成分运动信号曲线，其中，实线为各个特征点的横向坐标（X轴坐标）的变化，虚线为各个特征点的纵向坐标（Y轴坐标）的变化。 

图4（h）是对图3（b）中原始曲线分解得到的高频成分运动信号曲线，其中，实线为各个特征点的横向坐标（X轴坐标）的变化，虚线为各个特征点的纵向坐标（Y轴坐标）的变化。 

图4（i）是对图3（a）中原始曲线分解得到的X轴平移运动信号曲线与跟踪图2（a）中的骨骼点得到的X轴曲线的对比，其中，实线为分解得到的X轴平移运动信号曲线，虚线为跟踪骨骼点得到的X轴曲线。 

图4（j）是对图3（b）中原始曲线分解得到的X轴平移运动信号曲线与跟踪图2（b）中的骨骼点得到的X轴曲线的对比，其中，实线为分解得到的X轴平移运动信号曲线，虚线为跟踪骨骼点得到的X轴曲线。 

图4（k）是对图3（a）中原始曲线分解得到的Y轴平移运动信号曲线与跟踪图2（a）中的骨骼点得到的Y轴曲线的对比，其中，实线为分解得到的Y轴平移运动信号曲线，虚线为跟踪骨骼点得到的Y轴曲线。 

图4（l）是对图3（b）中原始曲线分解得到的Y轴平移运动信号曲线与跟踪图2（b）中的骨骼点得到的Y轴曲线的对比，其中，实线为分解得到的Y轴平移运动信号曲线，虚线为跟踪骨骼点得到的Y轴曲线。 

图4（m）是对图3（a）中原始曲线分解得到各个运动得到剩余噪声曲线，其中，实线为各个特征点的横向坐标（X轴坐标）的变化，虚线为各个特征点的纵向坐标（Y轴坐标）的变化。 

图4（n）是对图3（b）中原始曲线分解得到各个运动曲线得到剩余噪声曲线，其中，实线为各个特征点的横向坐标（X轴坐标）的变化，虚线为各个特征点的纵向坐标（Y轴坐标）的变化； 

图4中横坐标“Image Frame”表示“图像帧数”；纵坐标'displacement/mm'表示“运动位移/mm”。 

具体实施方式

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

本发明具体涉及一种血管结构特征点的选取和傅里叶级数展开相结合的数据分离方法。该方法可以解决人体的平移运动信号，心脏运动信号和呼吸运动信号等多种运动的分离难题。 

本发明的目的在于提供一种X射线血管造影图像中多运动参数的分离估计方法。本方法基于频域滤波和多变量寻优，能够精确分离各个运动参数，简单快速、无伤害，而且无论单臂还是双臂造影都可以广泛应用。 

X射线血管造影图像不仅记录了心脏自身运动在冠状动脉上的投影，同时还包含了人体的呼吸运动引起的冠状动脉在造影面上的二维运动和病人的二维平移运动。 

其步骤包括： 

傅里叶级数展开和多变量寻优的方法分离平移信号，心脏信号与呼吸运动信号。具体算法如下： 

Step1：对血管结构特征点在X射线血管造影图像序列中进行跟踪，得到点的跟踪曲线s_i(n)，i=1,…，I，n=1,...,N。其中I为特征点个数，N为X射线血管造影图像序列中图像的总帧数。其中，血管结构特征点需要能够真实综合反映出血管整体的运动信息，因此，根据生理学和解剖学等知识的指导，所选的特征点为血管的分叉点。血管结构特征点的得到是首先对X射线血管造影图像进行预处理，通过分割和细化得到血管骨架，然后通过血管骨架选取血管结构特征点。本发明采用的跟踪方法是通过X射线血管造影图像骨架序列不同帧的每个分叉点建立对应关系得到每个特征点在每一帧的坐标。 

Step2：设平移运动变化帧的序列{N_t,t=1,....,T}，平移运动的斜率角序列{α_t,t=1,....,T-1}，其中T为运动方向变化的次数；根据下式可以得到模拟平移曲线d_α(n)。 

$d_{\mathrm{\α}}\left(n\right)=\left\{\begin{array}{cc}0& n\∈(0,N_1]\\ \tan \left({\mathrm{\α}}_1\right)(n-N_1)& n\∈(N_1,N_2]\\ d_{\mathrm{\α}}\left(N_1\right)+\tan \left({\mathrm{\α}}_2\right)(n-N_2)& n\∈(N_2,N_3]\\ ...& \\ d_{\mathrm{\α}}\left(N_{T-2}\right)+\tan \left({\mathrm{\α}}_{T-1}\right)(n-N_{T-1})& n\∈(N_{T-1},N_T]\\ d_{\mathrm{\α}}\left(N_T\right)& n\∈(N_T,N]\end{array}\right.$

通过造影图确定X轴平移运动变化帧的序列

其中Tx为X轴移动方向变换的次数，以斜率角

范围为(-π/2,π/2)可以 通过上述方法得到X轴模拟平移曲线

同样，通过造影图确定Y轴移动变化帧的序列

其中T_y为Y轴移动方向变换的次数，以斜率角

范围为(-π/2,π/2)可以通过上述方法得到Y轴模拟平移曲线

使 $d_{\mathrm{\α}}\left(n\right)=({d_{\mathrm{\α}}}_x\left(n\right),{d_{\mathrm{\α}}}_y\left(n\right))$ 表示病人的平移运动。 

Step3：根据造影图序列确定心脏运动信号周期N_c。对去除了平移运动后的多运动的综合信号

即

根据心脏运动周期Nc进行傅立叶频域滤波处理，得到心脏运动信号

Step4：对不包含平移运动信号和心脏运动信号的剩余分量

以呼吸运动信号周期N_αr=m_cN_c,m_c∈[3,10]的范围，对每一个呼吸周期进行傅立叶频域滤波，以最接近

的拟合曲线

为最优准则搜索得到相对于当前平移曲线的最佳呼吸运动信号曲线和最佳呼吸运动信号周期

Step5：对不包含平移运动信号和心脏运动信号的剩余分量

减去其拟合曲线

后的曲线

的振幅若小于3个像素，则可以认为没有高频成分运动信号；否则，以高频成分运动信号周期N_αh=m_hN_c,m_h∈[1/7,5/7]的范围，对每一个高频成分运动周期进行傅立叶频域滤波，以最接近

为最优准则搜索得到相对于当前平移曲线的最佳高频成分运动信号曲线

和最佳高频成分运动信号周期

Step6：综合运动估计曲线为 ${\hat{s}}_{\mathrm{i\α}}\left(n\right)=d_{\mathrm{\α}}\left(n\right)+{\hat{r}}_{\mathrm{i\α}}\left(n\right)+{\hat{c}}_{\mathrm{i\α}}\left(n\right)+{\hat{h}}_{\mathrm{i\α}}\left(n\right)$ i∈[1,I]，以 与s_i(n)最接近为最优准则，得到最佳的曲线平移运动信号曲线，心脏运动曲线，呼吸运动信号曲线和高频成分运动信号。 

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

以下结合附图和实例对本发明进一步说明： 

本发明提出了一种直接从X射线造影序列图像中提取平移运动信号，呼吸运动信号和心脏运动信号等多运动参数的方法，通过在冠脉血管上选取一些特征点并在序列中进行跟踪得到它们随时间的运动曲线，然后通过傅里叶级数展开，频域滤波和多变量寻优的方法分离得到单独的平移信号，心脏信号与呼吸运动信号等信号。具有很好的临床适用性。 

临床上测得的正常的心脏运动为60～100次/分钟，其周期为0.6-1.0s，而将频率高于140次/分钟的周期运动则诊断为病理性异常运动，周期小于0.42s。如果能够用冠脉造影图分析出运动频率高于正常心脏运动，而且该高频信号的幅度超过2-3像素时，便可以合理的认为该高频信号为某种有用信号（或诊断病变信号，本文称为高频成分）。这是因为，在对冠脉造影图像进行处理的过程中，因分割或骨架提取的不精确性会导致一定的误差，但是这种误差一般不会超过3个像素。所以，当我们分析出存在这种高频运动而且幅度大于3个像素时，就可以为医生提供一种新的有用诊断信号。而呼吸运动信号的频率相比这两种运动则要小的多，呼吸运动信号的周期也就长得多，一般为3-6s，安静的时候可能更长。 

根据造影图中各运动成分（包括平移信号、心脏运动信号、呼吸运动信号及高频成分运动信号等）的这些特点，本文通过傅立叶级数展开结合频域滤波的方法对这些成分进行分离。 

下面首先介绍一下离散傅立叶级数变换。 

1）离散傅立叶级数变换基础 

由傅立叶级数展开的原理，连续时间周期信号可以用傅立叶级数来表达，与此相应，周期序列也可用离散傅立叶级数来表示。若周期序列x_p(n)的周期为N_p，那么 

x_p(n)=x_p(n+rN_p)（r为任意整数）     （1） 

x_p(n)可以通过公式（2）和（3）进行离散的傅立叶级数变换 

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

$x_p\left(n\right)=\frac{1}{N_p}\underset{k=0}{\overset{N_p-1}{\mathrm{\Σ}}}{X}_p\left(k\right)\·e^{j\left(\frac{2\mathrm{\π}}{N_p}\right)\mathrm{nk}}---\left(3\right)$

对于式（3）作如下解释：式中

是周期序列的基频分量，

就是k次谐波分量，各次谐波的系数为X_p(k)；全部谐波分量只有N_p个是独立的，因为

因此，级数取和的项数是k=0,...,N_p-1，共N_p个独立谐波分量。而式（2）正是由x_p(n)决定系数X_p(k)的求和公式。 

2）分离算法 

假设造影图图像序列中冠脉血管上第i个特征点P_i(x,y)的X轴坐标的运动曲线为x_i(n)（n为造影帧帧数，代表离散化的时间，序列总长度为N），Y轴坐标的运动曲线为y_i(n)，令s_i(n)=(x_i(n),y_i(n))，可将s_i(n)表示为: 

s_i(n)=d_i(n)+r_i(n)+c_i(n)+h_i(n)+t_i(n)i∈[1,I]     （4） 

其中I表示特征点的个数，d(n)=(x_d(n),y_d(n))表示随时间变化的病人的刚性平移运动，r_i(n)=(x_ri(n),y_ri(n))表示呼吸运动引起的运动，c_i(n)=(x_ci(n),y_ci(n))表示心脏的运动引起的血管上各点的运动，h_i(n)=(x_hi(n),y_hi(n))表示病人可能具有的高频运动引起的运动，t_i(n)=(x_ti(n),y_ti(n))表示其他噪声分量（如因为血管分割、细化等原因造成的）。为方便表示，后面都用s_i(n)表示提取的血管特征点沿X轴，Y轴的坐标变化曲线，d(n),r_i(n),c_i(n),h_i(n)和t_i(n)分别表示 由平移运动信号、呼吸运动信号、心脏运动信号、高频成分运动信号以及噪声分量沿X轴、Y轴的坐标变化曲线。 

假设心脏运动信号周期为N_c，呼吸运动信号周期为N_r，高频运动信号周期为N_h，则c(n)、r(n)和h(n)根据式（3）可变为： 

$\begin{array}{c}c\left(n\right)=\frac{1}{N_c}\underset{k=0}{\overset{N_c-1}{\mathrm{\Σ}}}C\left(k\right)\·e^{j\left(\frac{2\mathrm{\π}}{N_c}\right)\mathrm{nk}}\\ r\left(n\right)=\frac{1}{N_r}\underset{k=0}{\overset{N_r-1}{\mathrm{\Σ}}}R\left(k\right)\·e^{j\left(\frac{2\mathrm{\π}}{N_r}\right)\mathrm{nk}}\\ h\left(n\right)=\frac{1}{N_h}\underset{k=0}{\overset{N_h-1}{\mathrm{\Σ}}}H\left(k\right)\·e^{j\left(\frac{2\mathrm{\π}}{N_h}\right)\mathrm{nk}}\end{array}---\left(5\right)$

其中，

是心脏运动信号的谐波分量，是呼吸运动信号的谐波分量，

是高频成分运动信号的谐波分量,C(k)、R(k)和H(k)分别表示c(n)、r(n)和h(n)的各次谐波系数。要使c(n)、r(n)和h(n)完全分离，则其谐波分量不能发生重合（假设此时平移运动d(n)已消除），即 

$\∀k_1\left(n_1\right)=1,...,N_c-1,k_2\left(n_2\right)=1,...,N_r-1,k_3\left(n_3\right)=1,...,N_h-1,$

（6） 

$\∃\left(\frac{2\mathrm{\π}}{N_c}\right)n_1{k}_1\≠\left(\frac{2\mathrm{\π}}{N_r}\right)n_2{k}_2\≠\left(\frac{2\mathrm{\π}}{N_h}\right)n_3{k}_3,$ （n₁,n₂,n₃为时域自变量，k₁,k₂,k₃为频域自变量） 

要满足式（6），必须使N_c、N_r和N_h两两互质。同时要分离c(n)、r(n)和h(n)，还必须满足s(n)的序列长度大于或等于s(n)的一个周期，即大于N_c·N_r·N_h。 

然而上述要求很难达到，就算N_c、N_r和N_h两两互质，第二个条件也几乎无法满足。因为一般情况下，造影剂在人体内停留的时间不会太长，一个造影图序列只会持续6s左右，不会超过10s，假设帧速为80ms/frame（实际可能更小），则N_c≈10,N_r>30，N_c·N_r·N_h>N_c·N_r>300>10/0.08，很难获得 一个完整周期的序列。而且病人的平移动信号曲线，呼吸运动信号周期确定的要求较高。 

虽然无法直接对心脏动信号、呼吸动信号、高频成分动信号和平移动信号曲线进行理想的分离，但是根据它们各自运动的特点，可以通过傅立叶级数展开，频域滤波和多变量寻优的方法来进行近似的分离。 

分离算法基于如下假设： 

（1）假设周期为N_a,当序列长度为周期的整数倍时且不是其他运动周期的整数倍时，该运动只在频率

上存在分量，这一点已经通过式（3）证明；而其他运动频率分量在处平滑。由于无论是呼吸运动信号，心脏运动信号还是高频成分运动信号都具有平滑性，且互相为对方的基频分量的可能性较小，完全有理由这样假设。 

（2）由于平移运动信号是随时间变化的非周期性的刚体运动，所以平移运动信号估计假设为折线运动，如下式： 

$d_{\mathrm{\α}}\left(n\right)=\left\{\begin{array}{cc}0& n\∈(0,N_1]\\ \tan \left({\mathrm{\α}}_1\right)(n-N_1)& n\∈(N_1,N_2]\\ d_{\mathrm{\α}}\left(N_1\right)+\tan \left({\mathrm{\α}}_2\right)(n-N_2)& n\∈(N_2,N_3]\\ ...& \\ d_{\mathrm{\α}}\left(N_{T-2}\right)+\tan \left({\mathrm{\α}}_{T-1}\right)(n-N_{T-1})& n\∈(N_{T-1},N_T]\\ d_{\mathrm{\α}}\left(N_T\right)& n\∈(N_T,N]\end{array}\right.---\left(7\right)$

其中，{N_t,t=1,....,T}为平移运动信号变化帧，{α_t,t=1,....,T-1}为平移运动信号的斜率角，其中T为运动方向变化的次数。 

如图1所示，为整体流程图。 

具体算法如下： 

Step1：对血管结构特征点在X射线血管造影图像序列中进行跟踪，得到点的跟踪曲线s_i(n)，i=1,…，I，n=1,...,N。其中I为特征点个数，N为X射线血管造影图像序列中图像的总帧数。其中，血管结构特征点需要能够 真实综合反映出血管整体的运动信息，因此，根据生理学和解剖学等知识的指导，所选的特征点为血管的分叉点。如图2所示，在投影角分别为（-26.5°，-20.9°）和（42.3°，26.8°）的一对造影图中，分别存在着5个编号后的特征点（如图3（a）、图3（b）中的标记点）。血管结构特征点的得到是首先对X射线血管造影图像进行预处理，通过分割和细化得到血管骨架，然后通过血管骨架选取血管结构特征点。本发明采用的跟踪方法是通过X射线血管造影图像骨架序列不同帧的每个分叉点建立对应关系得到每个特征点在每一帧的坐标。 

在理想情况下,可将s_i(n)表示为: 

s_i(n)≈d(n)+r_i(n)+c_i(n)+h_i(n)i∈[1,I]     （8） 

Step2:平移运动信号曲线估计； 

若存在病人平移，则通过造影图确定X轴移动变化帧的序列 其中T_x为X轴移动方向变换的次数，以及Y轴移动变化帧的序列

其中T_y为Y轴移动方向变换的次数。并以斜率角  $\{{{\mathrm{\α}}_t}_x,t_x=1,...,T_x-1\},$ $\{{{\mathrm{\α}}_t}_y,t_y=1,...,T_y-1\}$ 范围为(-π/2,π/2)分别得到X轴模拟平移曲线和Y轴模拟平移曲线

如公式（7）。使

表示病人的平移运动信号。 

去除了平移运动后的多运动的综合信号为

Step3:心脏运动曲线估计 

根据造影图序列确定心脏运动信号周期N_c，一般情况下，N_c≈0.8/f，其中，N_c表示在一个心动周期内的造影图帧数，0.8表示一个心动周期的时间，单位为秒，f表示帧频，即每秒钟造影图被拍摄的数量。通过对

选取长度为n_sc的序列

n_sc为N_c的整数倍。将

分别分解为X方向的运动曲线x_ic(n)与Y方向的运动曲线y_ic(n)，再分别对它们进行离散傅里叶展开， 得到各次谐波系数X_ic(k)和Y_ic(k)，k=0,...,n_sc-1; 

不包含平移运动信号和心脏运动信号的剩余运动信号曲线的频域响应为 

${{\hat{S}}_{\mathrm{i\α}}^2}^x\left(k\right)=\left\{\begin{array}{cc}{X}_i\left(k\right),& (k\≠\frac{n_{\mathrm{sc}}}{N_c}l,l=1,...,N_c-1)\\ (X_i\left(k\right)+X_i(k+1))/2& (k=\frac{n_{\mathrm{sc}}}{N_c}l,l=1,...,N_c-1)\end{array}\right.---\left(9\right)$

心脏运动的频域响应为： $C_{\mathrm{i\α}}^x\left(k\right)=X_i\left(k\right)-{{\hat{S}}_{\mathrm{i\α}}^2}^x\left(k\right)(k=1,...,n_{\mathrm{sc}}-1)$

对

进行离散傅里叶逆变换得到X轴心脏运动曲线

同理可得 

即心脏运动曲线

对不包含平移运动信号和心脏运动信号的剩余运动曲线进行傅里叶逆变换可得

Step4:呼吸运动信号曲线估计 

对呼吸运动信号周期N_αr=m_cN_c,m_c∈[3,10]，当原始运动信号序列长于呼吸运动信号周期N_αr，则通过对进行与心脏运动信号分离相同的方法进行傅立叶频域滤波得到的呼吸运动信号曲线

当呼吸运动信号周期长于原始运动信号序列，则不包含呼吸运动信号的一个周期，则用以下步骤分离呼吸运动信号： 

a.从序列

第j点开始取长度为N_αr/2的序列 

对其进行傅里叶离散展开得到各次谐波系数 ${\hat{S}}_{{\mathrm{iN}}_{\mathrm{\αr}}j}^{21}\left(k\right),k=0,...,N_{\mathrm{\αr}}-1.$

b.对于长度为呼吸周期信号一半的曲线将认为其包含了周期为N_αr/2的运动曲线。该运动只在频率

上存在分量。通过傅立叶频域滤波得到该运动曲线 $r_{\mathrm{i\α}{N}_{\mathrm{\αr}}j\left(n\right),j=1,...,n_{\mathrm{sc}}-N_{\mathrm{\αr}}/2+1.}$

c.将

取均值得到近似呼吸运动信号曲线 

然后，对

进行曲线拟合得到呼吸运动信号估计曲线

（因为呼吸运动信号的周期最长，同时比高频成分运动信号的值要大得多，所以该曲线是呼吸运动信号的大致形状）。通过下式（10）为最优准则得到最佳呼吸运动信号周期

以及最佳呼吸运动信号曲线

$e_{\mathrm{\αr}}=\underset{N_{\mathrm{\αr}}\∈[3N_c,10N_c]}{\min }\left\{\underset{i=1...I}{\max }\underset{n=1}{\overset{N}{\mathrm{\Σ}}}\right[r_{\mathrm{i\α}{N}_{\mathrm{\αr}}}\left(n\right)-{\hat{r}}_{\mathrm{i\α}}^{\′}\left(n\right){]}^2\}---\left(10\right)$

Step5:高频成分运动信号曲线估计 

对不包含平移运动信号和心脏运动信号的剩余分量

减去其拟合曲线

后的曲线即

若

振幅不超过3个像素，则认为不具有高频成分（因为在对冠脉造影图进行处理的过程中，因分割或骨架提取的不精确性导致的误差一般不会超过2-3个像素）。否则，对于高频成分运动信号曲线周期N_αh=m_hN_c,m_h∈[1/7,5/7]，通过对与心脏运动信号分离相同的方法进行傅立叶频域滤波得到高频成分运动信号曲线 

通过下式（11）为最优准则得到最佳高频成分运动信号周期

以及最佳高频成分运动信号曲线

$e_{\mathrm{\αh}}=\underset{N_{\mathrm{\αh}}\∈[1/7N_c,5/7N_c]}{\min }\left\{\underset{i=1...I}{\max }\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{[h_{\mathrm{i\α}{N}_{\mathrm{\αh}}}\left(n\right)-{\hat{h}}_{\mathrm{i\α}}^{\′}\left(n\right)]}^2\right\}---\left(11\right)$

Step6:运动估计曲线为 ${\hat{S}}_{\mathrm{i\α}}\left(n\right)=d_{\mathrm{\α}}\left(n\right)+{\hat{r}}_{\mathrm{i\α}}\left(n\right)+{\hat{c}}_{\mathrm{i\α}}\left(n\right)+{\hat{h}}_{\mathrm{i\α}}\left(n\right)$ i∈[1,I]，通过（12）式为最优准则，得到最佳平移运动信号曲线以及相应的心脏运动曲线 

呼吸运动信号曲线

和高频成分运动信号曲线

$e=\underset{\mathrm{\α}\∈(-\mathrm{\π}/2,\mathrm{\π}/2)}{\min }\left\{\underset{i=1...I}{\max }\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{[s_i\left(n\right)-{\hat{s}}_{\mathrm{i\α}}\left(n\right)]}^2\right\}---\left(12\right)$

通过对造影图序列的观察确定心脏运动周期为11frames，图4是对图3中标记点的运动进行运动分离的结果。 

从图4中可以看到，特征点的原始运动曲线受呼吸运动影响较大，比较凌乱，而且不同特征点的运动幅度不同，这是因为不同特征点位于心脏的不同部位，而心脏不同部位的运动情况是不一样的。分离后的心脏运动曲线显示了良好的周期性，将所有标记点的呼吸运动曲线在一起进行比较，可以发现同一个造影面内所有特征点的呼吸运动曲线相差不大（如图4（e）和图4（f））。图4（g）和图4（h）为分离得到的高频成分大约是心脏运动的0.54倍，周期约为6frames，且其幅度基本都超过3个像素，若是因为分割或骨架提取等导致的误差应该是在3个像素范围内，如图4（i），图4（j），图4（k）和图4（l），因此可以假定它是某种有规律的心脏异常运动（如心脏震颤等）。图4（i），图4（j），图4（k）和图4（l）为分离的平移曲线与手动跟踪的平移曲线的对比，可以看出两者除了因为手动跟踪引起的误差外，基本一致。图4（m）和图4（n）为分离各个运动的剩余曲线，其幅度都少于3个像素，可以看作是因为分割，骨架提取等算法误差引起的，运动分量基本完全分离。 

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

Claims (1)

1.一种X射线血管造影图像中多运动参数的分离估计方法，其特征在于，包括下述步骤:

(1)对血管结构特征点在X射线血管造影图像序列中进行跟踪，得到点的跟踪曲线s_i(n)，i＝１，…，Ｉ，n＝１，...，N;Ｉ为特征点个数，Ｎ为X射线血管造影图像序列中图像的总帧数;

(2)根据平移运动变化帧的序列｛N_t，t＝1，....，T}和平移运动的斜率角序列{α_t，t＝1，....，T-1}得到模拟平移曲线d_α(n);其中Ｔ为运动方向变化的次数;

(3)根据所述X射线血管造影图像序列确定心脏运动信号周期Ｎ_c;根据所述跟踪曲线s_i(n)和所述模拟平移曲线d_α(n)得到去除了平移运动后的多运动的综合运动曲线

根据心脏运动周期N_c对所述综合运动曲线

进行傅立叶频域滤波处理，得到心脏运动曲线

(4)根据所述综合运动曲线

和所述心脏运动曲线

得到不包含平移运动信号和心脏运动信号的剩余曲线在呼吸运动信号周期[3Ｎ_c，10N_c]的范围内，根据每一个呼吸周期对所述剩余曲线

进行傅立叶频域滤波处理，得到相应的呼吸运动曲线;以最接近所述剩余曲线

的拟合曲线为最优准则，得到相对于当前模拟平移曲线的最佳呼吸运动曲线

和最佳呼吸运动信号周期

(5)根据所述剩余曲线

和所述拟合曲线得到的曲线

判断所述曲线

的振幅是否小于3个像素，若是，则没有高频运动信号;若否，则在高频运动信号周期[1/7N_c，5/7N_c]的范围内，根据每一个高频运动周期对所述剩余曲线

进行傅立叶频域滤波处理，得到相应的高频运动曲线;以最接近所述曲线

为最优准则，得到相对于当前模拟平移曲线的最佳高频运动曲线

和最佳高频运动信号周期

(6)根据模拟平移曲线d_α(n)、呼吸运动曲线

心脏运动曲线

和高频运动曲线

获得综合运动估计曲线 ${\hat{s}}_{\mathrm{i\α}}\left(n\right)=d_{\mathrm{\α}}\left(n\right)+{\hat{r}}_{\mathrm{i\α}}\left(n\right)+{\hat{c}}_{\mathrm{i\α}}\left(n\right)+{\hat{h}}_{\mathrm{i\α}}\left(n\right);$ 以所述综合运动估计曲线

与所述跟踪曲线s_i(n)最接近为最优准则，得到最佳的平移运动曲线，心脏运动曲线，呼吸运动曲线和高频运动曲线。

Images