CN101317194A - 用于自动4d冠脉建模和运动矢量场估计的方法和设备 - Google Patents

Abstract

一种用于对解剖对象进行计算机辅助自动四维(4D)建模的方法，包括采集表示对象在整个周期的多个静态的三维(3D)模型的集合。在该3D模型的集合上执行4D对应估计以确定该3D模型中的哪些点最可能彼此对应，其中该4D对应估计包括下列中的一项或多项：(i)定义参考相位；(ii)执行脉管导向的对应估计；以及(iii)对4D运动数据进行后处理。该方法还包括用前向传播算法进行自动3D建模。

Description

用于自动4D冠脉建模和运动矢量场估计的方法和设备

本实施例通常涉及根据诊断图像数据对三维解剖对象进行计算机辅助重建，并且更具体地，涉及一种用于自动4D冠脉建模和运动矢量场估计的方法和设备。

在注射造影剂之后，能够用介入X射线系统对冠状动脉进行成像。由于冠脉的运动，从二维(2D)投影的集合生成三维(3D)重建仅能使用属于同一心脏相位的有限数量投影，这导致图像质量很差。因此，已经开发出从两个或两个以上的投影导出冠脉树的3D模型的方法。这些方法有的是基于一个X射线血管造影中的初始2D中心线和采用极线约束在同一心脏相位的其它血管造影中寻找对应的中心线点。结果，该算法对呼吸和其它残余的非周期运动非常敏感。

另一种方法是基于在3D中的前向传播算法。在后面的方法中，用于控制前向传播的速度函数定义为前向的边界体素属于脉管的概率。该概率通过将体素向前投影到同一心脏相位的每个脉管经滤波的投影中并乘以响应值来估算。应当注意到的是，这种算法对于在不同血管造影之间不对应的残余运动较不敏感。然而，这种在3D中的前向传播算法仅是半自动的。

例如，必须手动定义作为前向传播起点的3D种子点。必须手动定义每条脉管的3D终点。从终点到种子点，3D前向传播算法自动寻找关于速度函数的最快连接路径。在3D前向传播算法的一个方面，终点根据所考虑的重建体积大小导出。然而，这是非常不具体的标准，若设置过小则使该算法遗漏脉管分支；否则若该值设置过大，则前向传播超出脉管树体积的边界。很可能在大多数情况下，该准则不存在对整个脉管树避免上述伪影的单一值。需要一种对每条脉管进行优化的更加具体的准则。

此外，关于3D前向传播算法，对于不同脉管和脉管段根据其相关性进行寻找和排列称为“构建”。在3D前向传播算法的工作流程中，用户通过手动选择特定脉管并且手动为每条脉管定义种子点和终点而进行排列，因而手动地实现“构建”。

此外，3D前向传播算法仅提取单个心脏相位的冠脉模型和中心线。为了从来自不同心脏相位的模型或中心线的集合导出四维(4D)运动场，必须给出用于导出3D中心线上对应点的方法。

图1示意性示出了包括在同一心脏相位中借助于X射线荧光检查采集的两个(2)二维(2D)投影1和2的诊断投影数据集。注意，可以使用任何适当类型的心脏相位监测，例如，在采集X射线投影的同时记录心电图(ECG)。以不同投影角度记录的投影1和2中的每个都示出了患者的分支血管3。因此投影图像1和2从不同视角示出了同一血管3。为了采集投影数据集，将造影剂施用给患者，使得血管3在投影中显示为深色。

为了根据3D前向传播方法重建血管3的三维结构，种子点5初始设置在重建体积4内。然后通过按照传播准则在每种情况下定位体积4内属于血管3的相邻点而在体积4内重建血管3。为此，在二维投影1和2内分别属于各自点5的局部区域6和7在每种情况下都单独地服从于数学分析。在对与种子点5相邻的点进行定位之后，对与该点相邻的各点依次重复该过程，直到在体积4内重建了血管3的整个结构。

如果局部区域6和7的数学分析对属于投影数据集的所有或者大多数投影都给出正结果，那么将在每种情况下用每个传播步骤调查的点识别为属于血管(即，分别在该示例的投影1和2中)。通过按照两个投影1和2中记录的投影方向将点5投影到这两个投影的对应平面中来确定局部区域6和7。这是在图1中分别用箭头8和9来指示的。注意，虽然关于同一心脏相位的两个(2)投影对该公知3D前向传播方法进行了描述，但是其并不局限于两个(2)投影。

因此，期望一种用于克服本领域中这些问题的改进的方法和系统。

根据本公开的实施例，一种用于对解剖对象进行计算机辅助自动四维(4D)建模的方法，包括对表示出在整个周期对象的多个静态的三维(3D)模型的集合进行自动采集。在该3D模型的集合上执行4D对应估计以确定该3D模型的哪些点最有可能彼此对应，其中4D对应估计包括下列中的一项或多项：(i)定义参考相位；(ii)执行脉管导向的对应估计；以及(iii)对4D运动数据进行后处理。该方法还能够实现为成像系统以及以计算机程序产品形式。此外，根据本公开的一个实施例的方法还包括能够用前向传播算法进行自动3D建模。

图1示意性示出了包括两个(2)二维(2D)投影图像的诊断投影数据集；

图2是用根据本公开的一个实施例的建模方法获得的全自动提取后向投影到潜在心脏相位的两个投影图像中的3D中心线的示例；

图3是示出用根据本公开的一个实施例的建模方法获得的、沿在两个不同心脏相位提取的脉管的三个正交轴的投影的示例的示意图；以及

图4是根据本公开的另一实施例的成像设备的部分方框图。

在附图中，相同的附图标记表示相同的元件。此外，应注意的是附图可能不是按比例绘制的。

自动3D建模：

根据本公开的一个实施例，一种方法包括使用前向传播方法从门控旋转血管造影术X射线投影的自动3D脉管中心线提取。尤其，该方法包括用于从门控3D旋转X射线投影对冠脉中心线树进行自动提取的非交互算法，即，没有人交互。该方法利用前向传播逼近来选择属于冠状动脉的体素。前向传播速度受到3D脉管概率的控制，该概率定义为：将所考虑的体素前向投影到同一心脏相位的每个脉管经滤波的投影中，拾取2D响应像素值并将其组合。该方法还包括将2D响应值组合成3D脉管概率的不同方式。该方法又包括利用多个单相模型建立组合的多相模型。

还描述了另一方式，该方法包括用于从门控3D旋转X射线投影对冠脉中心线树进行提取的全自动算法。该算法可用于在舒张末期的心脏相位使用优质投影时。通过使用备选版本的前向传播算法，能够显著减少来自在收缩相位中几乎接触的脉管的近路-伪影以及鬼影(ghost)脉管的伪影。所有算法版本都具有有限的运动补偿能力，因而在找到最佳心脏相位之后，能够对具有残余呼吸运动的投影进行中心线提取。此外，也能够将单相模型进行组合，以便确定最佳心脏相位并降低不正确追踪脉管的概率。此外，能够找到在不同单相模型中的对应点，以便用该逼近生成全4D冠脉运动场。

因此，这里讨论的前向传播方法能够在没有人交互的情况下对冠脉血管中心线树进行自动提取。如上所述进一步的，前向传播模型对残余运动，尤其是呼吸所引起的残余运动相对不敏感。根据一个实施例，必须从ECG门控模型的集合中确定表示在最小运动的心脏相位处的冠脉脉管形状的模型。在中心线提取算法中，该算法能够基于前向传播逼近进行全自动的冠脉脉管中心线提取。

如在此所讨论的，自动3D前向传播算法使用门控投影作为输入。门控是根据同时记录的心电图(ECG)信号来执行的。该算法包括多个预备和分析步骤，包括：(i)对门控投影进行预滤波；(ii)查找种子点；(iii)前向传播；(iv)对所有脉管候选：(a)查找终点；(b)回溯；以及(c)剪切和构建；(v)查找“根弧(root arc)”；(vi)连接；(vii)加权；以及(viii)输出和进行连接以供输出。

对门控投影进行预滤波

在第一步骤中，将投影分类到相对于ECG信号的R峰具有相同延迟的组中。门控投影数据集包括与来自每个心脏周期的给定门控点最邻近的投影。该算法的所有下列步骤都在门控投影集上执行。在下一步骤中，使用多尺度脉管滤波器对投影进行滤波，其滤波器宽度从1到7像素。结果是2D响应矩阵R^2D的集合，该响应矩阵提供每个像素是否属于脉管的概率。多尺度脉管滤波器定义为所有尺度的Hessian矩阵的特征值中的最大值。为了避免边界伪影，能够通过半径为大约(0.98*投影宽度)的圆形掩模对脉管经滤波的投影进行剪切。

查找种子点

对于每个体素

能够通过使用锥束向前投影计算每个投影上的对应像素。锥束向前投影能够具有特征在于：若用n表示当前投影，则

和

是检测器平面的法向矢量，

是检测器原点，

是焦点，以定义每个投影的轨线数据。是所考虑的体素而

是其投影。检测器平面的尺寸由w_x和w_y(以mm为单位的宽度和高度)以及P_x和P_y(以像素为单位的宽度和高度)来确定。

在3D中检测器平面上的投影像素计算如下：

$\stackrel{\→}{P_n}=\left(\frac{(\stackrel{\→}{D_n}-\stackrel{\→}{F_n})\·\stackrel{\→}{e_{n,z}}}{(\stackrel{\→}{x^{3D}}-\stackrel{\→}{F_n})\·\stackrel{\→}{e_{n,z}}}\right)\·(\stackrel{\→}{x^{3D}}-\stackrel{\→}{F_n})+\stackrel{\→}{F_n}$ (公式1)

于是，投影上的对应的(x，y)坐标为：

$v_{n,x/y}=(\stackrel{\→}{P_n}-\stackrel{\→}{D_n})\·{\stackrel{\→}{e}}_{n,x/y}\·\frac{P_{x/y}}{w_{x/y}}$ (公式2)

因为系统几何数据对每个投影都是具体的，所以像素坐标v还取决于当前投影n。

假设在不同投影之间没有运动，将要在脉管中定位到体素的概率R^3D能够通过对所有对应像素乘以2D脉管结果值R^2D来获得：

$R^{3D}\left(\stackrel{\→}{x^{3D}}\right)=\underset{n=0}{\overset{N}{\mathrm{\Π}}}{R}^{2D}\left(\stackrel{\→}{v_n}\left(\stackrel{\→}{x^{3D}}\right)\right)$ (公式3)

因此种子点是通过在某个子体积中选择具有最大响应的体素来查找到的。当前，占整个体积大约11％的子体积通过这种方式接受检查，这是因为假设重要脉管(最好是根弧)位于该体积的颅半部分并且位于中心，所以子体积如下的确定：

0.25·x_max≤x＜0.75·x_max

0.25·z_max≤z＜0.75·z_max    (公式4)

0.5·y_max≤y≤0.95·y_max

其中，将y轴标定为沿头足方向。最大y值不应当达到y_max，这是因为脉管经滤波的投影的残余边界伪影可能会影响对适当种子点的寻找。

为了进一步加速，使用所有N个投影并不能完全计算出每个体素的3D响应值。如果，在计算了n个投影的乘积之后，中间值低于当前最高响应值之下，那么不需要计算剩余的N-n个投影，这是因为随着每次额外的乘法运算，中间响应值只能进一步减小。这产生取决于源数据的2到5的额外加速因子。

前向传播

在找到适当种子点之后，可以开始前向传播。对于之前已经检查过的每个体素，将存储特性值，这指示出从种子点开始向该体素前向传播有多“快”。因此，该值称为时间值并且在种子点处设为零。因比这些时间值随任意路径的增加对于可能是好的血管而言应当较低，而对“不好的”血管和伪影而言应当较高(陡峭的)。

在每个迭代步骤中，从前向上当前具有最低时间值的体素开始，计算每个相邻体素的3D脉管响应值，并且将其倒数添加到所考虑的开始体素的时间值上。如果在之前已经考虑过相邻体素，那么不再重复计算该值。因而，在λ₀步骤之后到达体素

的时间值

表示从种子点开始的最佳可能路径的历史记录，这是因为其包含所有先前体素的响应值：

$T\left(\stackrel{\→}{x^{3D}}\left({\mathrm{\λ}}_0\right)\right)=\underset{\mathrm{\λ}=1}{\overset{{\mathrm{\λ}}_0}{\mathrm{\Σ}}}{\left(R^{3D}\left(\stackrel{\→}{x^{3D}}\left(\mathrm{\λ}\right)\right)\right)}^{-1}$ (公式5)

有多种方式来计算每个体素的适当响应值R^3D。算法的整体质量主要取决于这里使用的逼近的质量。因而，已经尝试了不同的逼近，但是其中只有三种证明是可行的。

第一前向传播逼近(FP1)

一种简单且稳定的方式是将每个经滤波的投影上对应像素的所有响应值相乘：

$R^{3D}\left(\stackrel{\→}{x^{3D}}\right)=\underset{n=0}{\overset{N}{\mathrm{\Π}}}{R}^{2D}\left(\stackrel{\→}{v_n}\left(\stackrel{\→}{x^{3D}}\right)\right)$ (公式6)

其中，n覆盖门控投影并且R^2D是当前经滤波的投影上的对应像素值，如上所述，其坐标由

给出。因而，对于较好的响应而言R^3D较高，反之亦然。乘法运算实际上对非常低的R^2D响应而言也没有问题，这是因为即使远离脉管结构，R^3D响应实际上并不达到零。

如果在集合的几乎所有投影上的脉管是相似的并且具有相对高的质量，那么该逼近给出合理的结果。其在追踪弱且细的脉管上存在问题，因此甚至由于较大的脉管变细而在其实际上终结时不能被追踪。前向传播快速朝向“好的”脉管，但是当其变小时，前向过程变得越来越没有偏向性并且趋向于朝着脉管边界传播。因此，使用相对质量差的投影对整个脉管树进行合理地追踪将由于执行了许多次迭代(例如对于512³的分辨率需要大约3-5百万次)而消耗很多计算能力。不过，可能仍然无法完全追踪到脉管的外端。

第二前向传播逼近(FP2)

对于先前章节中描述的追踪细脉管问题的解决方案可以是优选相对于那些明显根本不在脉管上的体素具有低响应的体素。因此第二前向传播逼近试图强调与那些其反向投影像素的响应值更加不同的体素相比在所有投影中具有相对均匀响应的体素。该决定可能是错误的，这是因为甚至“正确”的体素也可能由于运动或者不好的投影/预滤波质量而在一些投影中具有不好的响应值。因为将每个经滤波的投影正规化到1，所以结果能够通过将其提升到低于1的幂(power)而得到强调，并通过将其提升到超过1的强度而得到抑制。为了描述某体素

的2D响应值如何均匀地分布，指数

现在作为正规化方差来计算：

$\mathrm{\η}\left(\stackrel{\→}{x^{3D}}\right)=\frac{\underset{m=1}{\overset{N}{\mathrm{\Σ}}}|R^{2D}\left(\stackrel{\→}{v_m}\left(\stackrel{\→}{x^{3D}}\right)\right)-\stackrel{\‾}{R_n^{2D}}|}{N\·{\left(\stackrel{\‾}{R_n^{2D}}\right)}^2}$ (公式7)

以及

$\stackrel{\‾}{R_n^{2D}}=\frac{\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{R}^{2D}\left(\stackrel{\→}{v_n}\left(\stackrel{\→}{x^{3D}}\right)\right)}{N}$ (公式8)

并且使用如下：

$R^{3D}\left(\stackrel{\→}{x^{3D}}\right)=\underset{n=0}{\overset{N}{\mathrm{\Π}}}{R}^{2D}{\left(\stackrel{\→}{v_n}\left(\stackrel{\→}{x^{3D}}\right)\right)}^{\mathrm{\η}\left(\stackrel{\→}{x^{3D}}\right)}$ (公式9)

该逼近优选弱的脉管但是将降低运动补偿能力。在一些情况下其趋向于不稳定。

第三前向传播逼近(FP3)

第三前向传播逼近用以说明两个投影m和n之间的投影角度差α_m-α_n，以对于取自相似角度的视图的那些优选从垂直视图选取的信息。这将最小化对两个投影中深度信息的误解。因为有两个以上可用的投影，所以所有投影(1...n₀)都要成对考虑并且各自的结果是通过乘法运算组合得到的。每对投影的响应值通过乘以其一致的2D响应值并且用投影差角度的正弦对其进行加权而计算：

$R^{3D}\left(\stackrel{\→}{x^{3D}}\right)=\underset{m=0}{\overset{N-1}{\mathrm{\Π}}}\underset{n=m+1}{\overset{N}{\mathrm{\Π}}}\left|\sin ({\mathrm{\α}}_m-{\mathrm{\α}}_n)\right|\·R^{2D}\left(\stackrel{\→}{v_n}\left(\stackrel{\→}{x^{3D}}\right)\right)\·R^{2D}\left({\stackrel{\→}{v}}_m\left(\stackrel{\→}{x^{3D}}\right)\right)$ (公式10)

正弦通过将指向体积中心M到检测器D的矢量的叉积除以其各自长度来计算获得：

$\sin ({\mathrm{\α}}_m-{\mathrm{\α}}_n)=\frac{|(\stackrel{\→}{D_m}-\stackrel{\→}{M})\×(\stackrel{\→}{D_n}-\stackrel{\→}{M})|}{|\stackrel{\→}{D_m}-\stackrel{\→}{M}|\·|\stackrel{\→}{D_n}-\stackrel{\→}{M}|}$ (公式11)

该第三前向传播逼近在追踪细的脉管时执行的很好并补偿残余运动。此外，第三前向传播逼近与第二前向传播逼近相比可以更稳定。

终止前向传播

取决于体积分辨率和投影的质量，存在合理迭代数目的经验值：

i_0，FP1≈0.03·像素的数目    (公式12)

关于第一前向传播，对于256³个体素，大约500k次迭代是足够的，而512³将需要大约4,000k次迭代以使前向传播到相似区域。然而，后者的迭代数目消耗大约八(8)倍之多的存储器和计算时间。第二和第三FP逼近仅需要大约一半的迭代就能获得相似的结果。

查找脉管段

在找到终点之后，对脉管中心线进行追踪、剪切，并且分别存储其各部分。用相同的方式处理连续脉管。因此对每个脉管候选及和其子脉管分别完成下列三个步骤：(1)查找终点；(2)回溯；以及(3)剪切和构建。

(1)查找终点

在完成前向传播之后，必须为每个脉管找到适当的终点。这通过将整个体积划分为n³个子体积来实现，其中在该阶段n＝50。在每个体积中，选择具有最高时间值的体素。该体素位于脉管外缘，这是因为在每个脉管中心处前向传播快，并且然后朝其边界缓慢加宽(引起高时间值)。

(2)回溯

使用最陡梯度法执行回溯。假定终点，引导回溯朝着关于当前体素具有最大时间值减少的体素。通过跟随在每个步骤处的最大减少，计算返回种子点的最佳路径。从前向传播的表面开始，其直接通向脉管中心并且然后沿着中心线通向种子点。如果通过先前的迭代在之前已经追踪了路径，将不再进行追踪。这通过3D位图来管理，在3D位图中标记追踪到的体素并在其每侧加以两个体素的额外安全区。这样可防止对相似(平行)路径的双重追踪。

(3)剪切和构建

应当注意的是，位于脉管边界处的体素不属于中心线，并且由此需要对这种体素进行剪切。通过递归算法来完成剪切，其中递归算法的任务在于将追踪到的中心线分裂成不同质量的段。在开始回溯的点处的段具有最差的质量并且从而被除去。

递归剪切算法假设每个脉管靠近种子点的质量最好并且向其回溯的起点降低。计算最初四分之一当前脉管体素的平均值，其中该计算得到的值然后在朝着追踪起点扫描时用作阈值。可以偶尔超过该阈值数倍，但是如果这些超过的数目超出公差值(例如，连续十(10)次的最大值)，那么将该特定斑点视为明显的质量裂口(breach)，并且将该脉管分裂成两个部分。这意味着，从质量较好的脉管段中切除质量最差的段，然后将其存储作为独立的脉管。然后以相同的方式处理该第二脉管，因而分离出独立脉管的段，等等。若剩余部分小于最小长度(例如，在十(10)个体素的位数上)，则中止递归算法。位于追踪起点的边界体素根据最小长度准则切除，或者如果其长度超过十(10)个体素，那么通过稍后讨论的加权算法将其评价为不重要的。

查找“根弧”

如在此所述，用于前向传播的种子点并不必须对应于根弧，根弧是冠状动脉树的流入节点。因此，追踪每个脉管回到该“错误”起点。为了估计根弧的实际位置，使用最长的三个单个脉管段的最佳入颅点。根据需要，然后用种子点与新的顶点之间连接的脉管段延伸其它脉管。

连接

到目前为止，脉管彼此之间没有关系。每个脉管的结束是由下列三个原因之一所导致的：i)已经到达了根弧，因而不需要连接；ii)该脉管原来是更长脉管的一部分，并且通过上述剪切和构建算法得到分离；以及iii)存在分叉，分叉意味着有另一脉管交叉，这已经在回溯阶段检测到。到这点为止，仅知道是某条路径是否之前已经追踪到了，但是并不知道哪个脉管使用它。通过选择几何上与每个脉管段的终点最接近的点来确定正确的后继脉管。因为在回溯阶段，所有血管都以升序进行索引，所以仅需要寻找比所考虑的脉管索引低的血管上的点。在连接之后，通过沿连接路径将所有脉管段的长度相加能够容易地计算出每个血管(从终点到根弧)的全长。

加权

在上述的步骤中，已经提取了大量的路径，但是其中只有一些是表示实际存在的脉管，而大多数是由伪影导致的，例如投影质量的不足、残余运动、透视缩短等。因此，必须确定它们中哪些最可能表示真实脉管。用于提取的路径候选的总体显著性的测量S能够包括多种因素：i)脉管段的长度或者总长度；ii)由时间值所确定的质量；iii)3D位置(可能借助于预先定义的模型)；以及iv)形状。根据显著值S，能够对所有路径候选进行分类，这使得能够选择最显著的路径来输出，其中用于输出的路径的最大数目能够由系统用户进行设置。显著值S的计算仍有待改善，这是因为在此的误判可能导致错误(“鬼影”)脉管的输出。在一个实施例中，S计算如下：

$S=|y_{\mathrm{end}}-y_{\mathrm{root}\_\mathrm{arc}}|\·\frac{l_{\mathrm{part}}^2}{T\left(\stackrel{\→}{x^{3D}}\left({\mathrm{\λ}}_{\mathrm{end}}\right)\right)}$ (公式13)

其中，y_end和y_{root_arc}分别是如上所述确定的当前脉管段终点和根弧的y坐标(沿头足旋转轴)。数量l_part是脉管段以体素为单位的长度，而

是脉管段终点的时间值。可以使用例如梯度准则自动估计合理数量的可提取脉管中心线。

输出和进行连接以供输出

当将中心线数据存储到文件中时，由于可以不选择一个或多个段的连接的路径以供输出，因此必须对连接进行检查并且必须重新连接脉管的某些部分。

根据本公开的实施例，改善的前向传播算法将现有公知的半自动3D算法转换为全自动4D算法。该方法解决了上面所讨论的各种问题，并提供如下解决方案：

1.种子点：根据一个实施例，通过对上述在每个血管造影中可观测的3D体积的中央颅侧子体积中的3D脉管响应进行评价，并且选择具有最大3D响应的点来自动定义种子点。在采集对应3D响应的X射线投影的同时可以使用任何适当类型的心脏相位监测，例如，心脏相位监测可以包括心电图(ECG)的记录。最大3D响应点位于脉管树上，但是不必须在主干分叉的流入节点上。一种备选方法是在上述体积表面的颅部上选择具有最大3D响应的点。在后面的实例中，这提供了位于充满造影剂的导管上的种子点，其经由主动脉自颅侧进入。

2.停止前向传播：执行前向传播迭代的数目从下列两者中任一导出：(i)前向传播体积的体素分辨率，或者(ii)通过对3D响应值沿提取的脉管减小进行分析。

3.终点：能够通过一种或者多种不同方法确定脉管的潜在终点。在第一实施例中，将前向传播体积分成大量子体积(例如，50³或者50*50*50)。在每个子体积内，选择具有最近前向到达的点作为回溯算法的起点。回溯算法如下：速度场沿着具有最陡梯度的路径追溯到种子点。在第二实施例中，在前向传播期间，该算法沿最陡梯度追踪路径，并且若检测到3D脉管响应显著减少则停止。在任何情况下，对潜在脉管终点的精确估计不是非常关键的，这是因为在以下构建步骤中，根据血管段的相关性对其进行分析和加权。

4.构建：通过动态构建算法将脉管分成不同段。动态构建算法确定所提取的具有匀质(homogenous)3D脉管响应的中心线的部分(section)。对每个脉管段进行的加权根据下列不同准则执行：(i)长度；(ii)3D脉管响应(对应于质量)；(iii)中心线的形状和位置(或者任选地，基于先验的冠脉模型)。自动选择经由最相关的加权脉管并使其构成3D算法的输出。图2包含向后投影到潜在心脏相位的两个投影(22和24)中的全自动提取的3D中心线的示例(20)，其是使用根据本公开的一个实施例的建模方法获得的。

4D算法：

根据本公开的一个实施例，自动4D冠状建模和运动矢量场估计方法在输入处需要表示贯穿整个心脏周期的所有静态的3D模型的集合，其通过对每个可区分的心脏相位重复上述过程而获得。该方法通过将不同模型的分叉和其它形状特性进行匹配来确定不同模型的对应点。一种其中采用4D信息的可能应用是导出用于门控或者运动补偿3D重建的最佳心脏相位。

根据本公开的实施例的方法提供一种用于冠脉中心线提取和建模的全自动、鲁棒的4D算法。该方法能够处理由于残余运动而导致相同心脏相位的血管造影中的不一致。此外，根据本公开的实施例的方法提供了对现有公知3D前向传播算法的改进，其中改进允许诸如4D运动补偿重建和建模的新应用。

表示贯穿整个心脏周期的所有静态的3D模型的集合能够通过对每个可区分的心脏相位重复3D建模过程而获得。取决于旋转运转期间的最小心博率f_h，min(以每分钟心跳为单位，bpm)和采集帧率f_a(以1/s为单位)，可区分的心脏相位p_N的数目等于：

$\mathrm{PN}=60\frac{s}{\min }\·\frac{f_a}{f_{h,\min }}$

这意味着已经创建了p_N个独立3D模型。该值的范围大约是：对于采集帧率f_a为25fps(帧每秒)和心博率f_h为100bpm(每分钟心跳)时大约是15，到f_a为30fps和心博率f_h为45bpm时的大约40。4D对应估计的任务在于确定模型中的那些点最可能彼此对应，这使得能够估计某部分脉管树在整个心脏周期中的运动。必须考虑到脉管的纵向运动和在3D建模过程中导致的模糊问题，这些问题使得4D对应估计更为困难。对应估计通过执行下列步骤来执行：

1.参考相位(稳定相位)的定义；

2.脉管导向的对应估计；

3.4D运动数据的后处理。

1.参考相位的定义

为了估计稳定的4D对应，需要决定在所述步骤期间提取的多个潜在脉管结构中的哪些在整个心脏周期期间具有最高显著性。在3D算法中，根据血管段假定的显著性对其进行加权，但这对于每个单个3D模型是独立完成的，其导致了所提取的血管在不同心脏相位的波动。因此，必须在对应估计之前定义具有所有期望的被提取脉管的参考相位p_r(稳定相位)。这可以自动或者手动完成。

自动定义：选择表示与35％RR最接近的相位的3D模型或者选择包含三个最长血管的模型，该模型在实践中非常可能是低运动的相位并且因此是良好提取质量的相位；或者选择包含三个最长脉管的模型。注意，RR表示由ECG的两个随后R峰定义的时间间隔，其中该ECG受控于R峰，并且每个R峰表示在心脏收缩之前的电脉冲。

手动定义：根据对所有提取的3D模型进行目视检查(例如，使用如图

手动定义：根据对所有提取的3D模型进行目视检查(例如，使用如图3中所示所有模型具有投影的总体绘图30)，用户能够手动定义最适当的心脏相位并重新开始该算法。图3示出了所提取的脉管在不同心脏相位的两个投影的示例30。表示43.5％RR的心脏相位的上列32示出了三个正确提取脉血管，其将该相位限定为潜在参考相位，而在下列34(5％RR)中示出的脉管的质量较差。

2.脉管导向的对应估计

使用每个模型中的一个稳定点对在参考相位p_r处的每个提取的脉管独立地执行对应估计。当第一次执行该步骤时，主干分叉(“根弧”)用作稳定点，而在后面的迭代中，使用具有可能更高精度的子分叉点。该算法利用了如下事实：在心脏周期期间，脉管的弧长λ不相当地改变(总体上小于2％)。3D坐标为：

$\stackrel{\→}{x^{3D}}=\stackrel{\→}{x^{3D}}\left(\mathrm{\λ}\right)$

任何脉管点用脉管的弧长λ进行参数化，其取决于所考虑的相位数p、所考虑的脉管数v以及沿脉管路径的体素数i：λ＝λ(p，v，i)。如果，在下文中，文本中提及整个脉管，那么省略体素数i。

创建等间距形式(version)的当前所考虑的参考相位脉管λ(p_r，v_r)和当前目标相位脉管λ(p，v)，其保持预定的间距s(当前设置为2mm)，这是因为最初的3D模型的点对点距离根据√3或者更多因素变化，这是由对角体素距离和连接间隙所造成。它们表示从稳定点到脉管末端的整个路径。在进行等距离间隔之前先对脉管点坐标进行低通滤波，以消除源自前向传播的体素表达的量化效应，从而提供一个稳定的弧长准则。低通形式的脉管λ(p，v)表示为λ′(p，v)。对两个脉管进行逐点比较，并且总体相似性准则C计算如下：

$C(v_r,v)=\frac{\underset{i=0}{\overset{i_{\max }-1}{\mathrm{\Σ}}}\∈\·|\stackrel{\→}{x^{3D}}\left({\mathrm{\λ}}^{\′}(p_r,v_r,i)\right)-\stackrel{\→}{x^{3D}}\left({\mathrm{\λ}}^{\′}(p,v,i)\right)|}{i_{\max }^2\·\underset{i=0}{\overset{i_{\max }-1}{\mathrm{\Σ}}}\∈},$

i_max＝min[λ′(p_r，v_r，end)，λ′(p，v，end)]，

∈＝i+1.

较小的相似性准则C指示两个当前脉管之间更好的对应。因此，认为具有最小C的脉管组合是相等同的。对源脉管v_r和目标相位脉管v的每个组合以及每个可能的目标相位p≠p_r重复该过程。对应血管的所有对应坐标最终都存储在具有索引[0...p_N-1](相位)以及[0...i_max-1](对应3D点)的动态阵列A(p，i)中(称为运动场)。

3.4D运动数据的后处理

在对应估计过程中，将每个对应血管表示成从参考点开始(通常是根弧)，这导致脉管树的多个部分被多次表示。这产生高的局部点密度，需要加以变稀疏，以避免奇异性和其它模糊。这种降低通过以下方式来实现：计算属于某相位的点的每个组合之间的欧几里德距离d，并且如果该距离低于阈值就删除所述点中的一个，其中t＝0.5，s＝1mm：

$d(p,i_1,i_2)=|\stackrel{\→}{A}(p,i_1)-\stackrel{\→}{A}(p,i_2)|$

能够对贯穿所有心脏周期的得到的对应“根弧”点进行异常值检查。如果特定相位的根弧到中值(平均)位置的距离超过给定阈值，那么从该模型中排除该心脏相位。能够以类似的方式对所有其它分叉和单个点进行处理。

现在翻到图4，其中示出的成像设备是C型臂X射线设备，其包括C型臂10，该C型臂借助于托架11悬挂于例如天花板(未示出)。X射线源12和X射线图像转换器13可在C型臂10上移动地引导，从而可以在不同投影角度对平躺于C型臂10中心的病床14上的患者15的多个二维投影X射线图像进行记录。X射线源12与X射线图像转换器3的同步移动受控制于控制单元16。在图像记录期间，X射线源12和X射线转换器13同步地围绕患者15行进。将X射线图像转换器13所生成的图像信号传送到受控的图像处理单元17。使用ECG设备18对患者15的心搏进行监测。ECG设备18将控制信号传送到图像处理单元17，使得后者能够存储在心博周期的相同相位的情况下的多个二维投影，以执行冠状动脉的血管造影检查。图像处理单元17包括程序控制，借助于该程序控制，根据3D前向传播方法能够执行用所采集的投影数据集检测的脉管树的三维模型。此外，图像处理单元17包括另一程序控制，借助于该程序控制，根据本公开的实施例能够执行4D建模。然后可以在与图像处理单元17相连接的监控器19上以任何适当方式显示4D建模以及一条或多条重建的血管。

尽管上面仅详细描述了一些示范性实施例，但是本领域的技术人员应当容易理解的是，在本质上不脱离本公开的实施例的新颖教导和优点的情况下，在这些示范性实施例中可以进行许多修改。例如，本公开的实施例可以用于诸如心脏静脉的其它周期性移动结构或者更通常为树状结构。因此，所有这种修改都旨在包括在如在下面的权利要求书中所限定的本公开实施例的范围内。在权利要求书中，模块加功能的条款旨在覆盖作为执行所引用的功能的所述结构，并且不但是结构的等同，而且是等同的结构。

此外，在一项或者多项权利要求中置于括号内的任何附图标记不应解释为对权利要求的限制。词语“包括”和“包含”等并不排除存在任何权利要求或说明书中所列出的那些作为整体之外的元件或步骤。单数标记的元件并不排除复数标记的这种元件，反之亦然。一个或多个实施例可以借助于包括多个独立元件的硬件来实现，和/或借助于适当编程的计算机来实现。在列举了多个模块的装置权利要求中，这些模块中的多个可以实施为一个和相同项的硬件。在互不同的从属权利要求中引用的某些测量的这个事实并不指示这些测量的组合不能用来得到优点。

Claims (35)

1、一种对解剖对象进行计算机辅助建模的方法，包括：

采集所述解剖对象的门控旋转X射线投影；以及

使用前向传播方法从所述门控旋转X射线投影中自动提取三维(3D)脉管中心线，其中，所述前向传播方法包括自动查找在不同单相前向传播中的点。

2、如权利要求1所述的方法，其中，响应于查找在所述不同单相前向传播中的对应点，能够作为所述对应点的函数来生成四维(4D)冠脉运动场。

3、如权利要求1所述的方法，其中，自动提取3D脉管中心线的步骤包括下列的一项或多项：

(i)对所述门控旋转X射线投影进行预滤波，其中，预滤波包括将所述门控投影分类到数据集中，其中，所述门控投影数据集包括与来自每个心脏周期的给定门控点最邻近的投影；

(ii)查找种子点，其中，所述种子点包括在给定子体积中具有最大3D脉管响应的体素；

(iii)执行前向传播，其中，对所述前向传播执行的迭代数目得自：(a)前向传播体积的体素分辨率；或者(b)通过分析沿提取的脉管候选三维(3D)响应的减少；

(iv)对所述提取的脉管候选和对应子脉管执行：(a)查找脉管终点；(b)将脉管中心线沿着具有最陡梯度的路径回溯到所述种子点；以及(c)剪切和构建，其中，所述剪切和构建将所述脉管分成不同段，并且进一步确定所提取的具有匀质3D脉管响应的中心线的部分；

(v)查找根弧，所述根弧对应于冠状动脉树的流入节点；

(vi)将相关脉管段彼此连接，其中，通过选择几何上与给定脉管段的终点最接近的点来确定对应的后继脉管段；以及

(vii)对脉管段进行加权，其中，对每个脉管段进行的加权根据下列不同准则中的一个或多个来执行，包括：(a)脉管段的长度；(b)3D脉管响应；以及(c)所述中心线的形状和位置。

4、如权利要求3所述的方法，其中进一步，所述投影数据集相对于ECG信号的R峰具有相同的延迟。

5、如权利要求3所述的方法，其中，预滤波还包括使用多尺度脉管滤波器对所述门控旋转X射线投影进行滤波，所述多尺度脉管滤波器定义为所有尺度的Hessian矩阵的特征值中的最大值。

6、如权利要求3所述的方法，其中，预滤波还包括用半径为投影数据集宽度的大约百分之九十八(98％)的圆形掩模对所述投影数据集进行剪切。

7、如权利要求1所述的方法，其中，对所述门控旋转X射线投影进行的门控是根据同步记录的心电图(ECG)信号来执行的。

8、如权利要求1所述的方法，还包括：

对所述门控旋转X射线投影进行预滤波，其中，将所述投影分类到相对于ECG信号的R峰具有相同延迟的组中。

9、如权利要求1所述的方法，还包括：

从具有残余呼吸运动的所述门控旋转X射线投影中确定最佳心脏相位；以及

使用所述前向传播方法从所述门控旋转X射线投影中自动提取三维(3D)脉管中心线，进一步作为所述最佳心脏相位的函数。

10、如权利要求1所述的方法，还包括：

使用3D脉管概率控制所述前向传播方法的速度。

11、如权利要求10所述的方法，其中，所述3D脉管概率通过将所考虑的体素向前投影到同一心脏相位的每个脉管滤波的投影中，选择二维(2D)响应像素值，并将所述2D响应像素值组合成所述3D脉管概率来定义。

12、如权利要求1所述的方法，其中，所述前向传播选择属于冠状动脉的体素。

13、如权利要求1所述的方法，其中，所述前向传播模型利用一个以上的单相前向传播以建立组合的多相前向传播。

14、如权利要求1所述的方法，还包括：

查找在不同所述单相前向传播中的对应点；以及

生成作为所述不同单相前向传播中的所述对应点的函数的四维(4D)冠脉运动场。

15、一种成像设备，包括：

用于生成投影数据集的模块，所述集合包括从不同投影方向记录的患者身体部分的多个旋转X射线投影；并且具有用于根据所述投影数据集重建三维对象的计算机模块，其中，所述计算机模块包括计算机控制，所述计算机控制用于根据权利要求1所述的方法对所述对象执行计算机辅助建模。

16、如权利要求15所述的成像设备，还包括ECG控制，其中对旋转X射线投影进行的记录能够按照所述患者的所述心脏周期来控制。

17、一种计算机程序产品，包括：

具有计算机可执行的指令集的计算机可读介质，用于根据如权利要求1所述的方法对对象执行计算机辅助建模。

18、一种对解剖对象进行计算机辅助四维(4D)建模的方法，包括：

采集表示所述对象在整个周期的多个静态的三维(3D)模型的集合；以及

在所述3D模型集合上执行4D对应估计以确定所述3D模型中的哪些点最可能彼此对应，其中，所述4D对应估计包括下列的一项或多项：(i)定义参考相位；(ii)执行脉管导向的对应估计；以及(iii)对4D运动数据进行后处理。

19、如权利要求18所述的方法，其中，采集步骤包括采集表示贯穿整个心脏周期的所有静态的3D模型的集合。

20、如权利要求18所述的方法，其中，所述周期包括心脏周期，并且其中，采集所述3D模型的集合的步骤还包括通过将3D建模过程重复所述心脏周期的可区分心脏相位的数目来进行采集。

21、如权利要求20所述的方法，其中，可区分心脏相位的所述数目取决于旋转运转期间的最小心博率和采集帧率。

22、如权利要求18所述的方法，其中，所述4D对应估计使得能够对整个心脏周期的某部分脉管树的运动进行估计。

23、如权利要求18所述的方法，其中，所述参考相位包括预定义的稳定相位，其在所述脉管导向的对应估计之前被定义。

24、如权利要求18所述的方法，其中，定义所述参考相位的步骤包括自动定义或手动定义中的一种。

25、如权利要求24所述的方法，其中，所述自动定义选择下列中的一项：(i)表示与给定百分比RR最接近的期望相位的3D模型，其中对应于良好提取质量的相位，所述期望相位具有较低的运动；或者(ii)包含三个最长脉管的3D模型。

26、如权利要求24所述的方法，其中，所述手动定义包括：(i)对提取的3D模型进行目视检查；(ii)根据所述目视检查的3D模型手动定义最合适的心脏相位；以及(iii)用所述手动定义的参考相位开始所述4D对应估计。

27、如权利要求18所述的方法，其中，使用每个3D模型中的稳定点对所述参考相位处的每个提取的脉管独立地执行脉管导向的对应估计。

28、如权利要求27所述的方法，其中，对于初始的脉管导向的对应估计，所述稳定点包括主干分叉，而对于脉管导向的对应估计的一个或多个后续迭代，所述稳定点包括子分叉点。

29、如权利要求27所述的方法，其中，所述脉管导向的对应估计：(i)用所述脉管的弧长λ对任意脉管点的3D坐标进行参数化，这取决于所考虑的相位数p、所考虑的脉管数v、以及沿所述脉管路径的体素数i；(ii)创建等间距形式的当前所考虑的参考相位脉管和当前目标脉管，其保持预定的间距；(iii)对脉管点坐标执行低通滤波以提供稳定的弧长准则；(iv)逐点地对两个脉管进行比较；以及(v)计算整体相似性准则，作为所述两个脉管的所述逐点比较的函数。

30、如权利要求29所述的方法，其中，所述脉管导向的对应估计还包括对源脉管和目标相位脉管的每个组合以及除所述参考相位之外的每个可能的目标相位同样地重复步骤(i)-(v)，并且还包括将对应脉管的所有对应坐标存储到以相位和对应3D点作为索引的动态运动场阵列中。

31、如权利要求18所述的方法，其中，对4D运动数据进行后处理的步骤包括对于异常值检查在整个所述心脏周期上的点，并且响应于查找特定相位中的根弧点到高于给定阈值的中值位置的距离，对4D运动数据进行后处理的步骤还包括从4D建模中排除所述心脏相位。

32、如权利要求18所述的方法，其中，所述对4D运动数据进行后处理的步骤包括计算属于某个相位的点的每个组合之间的欧几里德距离d，并且如果所述距离低于阈值就放弃所述点中的一个。

33、一种成像设备，包括：

用于生成投影数据集的模块，所述集合包括从不同投影方向记录的患者身体部分的多个二维投影；并且具有用于根据所述投影数据集重建三维对象的计算机模块，其中，所述计算机模块包括计算机控制，所述计算机控制用于根据权利要求18所述的方法对所述对象执行计算机辅助四维建模和运动补偿重建。

34、如权利要求33所述的成像设备，还包括ECG控制，其中对二维投影进行的记录能够按照所述患者的所述心脏周期来控制。

35、一种计算机程序产品，包括：

具有计算机可执行的指令集的计算机可读介质，用于根据如权利要求18所述的方法对对象执行计算机辅助四维建模和运动补偿重建。

Images