从X射线锥形束数据中重建三维图像数据组

摘要

本发明涉及一种利用成像系统(1)产生对象(20)的三维图像数据组的方法，所述方法包括：当所述源(4)沿着由通过计数器参数(λ)连续地编号的一系列源点所描述的、大体上平面的轨迹(40)围绕待成像的对象移动时，采集(100)来自于探测器(10)的锥形束数据的一系列二维阵列；通过执行以下步骤从锥形束数据中重建三维图像：A)在固定的射线方向(α)上关于源轨迹的计数器参数(λ)对锥形束数据微分；B)利用类似Hilbert的滤波器对导数进行滤波以产生经滤波的锥形束数据；C)在步骤A)之前，或者在步骤B)之后，分别将采集的锥形束数据或者经滤波的锥形束数据与冗余加权函数(w)相乘；并且D)对锥形束数据进行反投影以计算三维图像数据组。

说明书

技术领域

本发明一般地涉及一种能够产生对象的三维(3D)图像数据组的方法和成像系统。成像系统包括以锥形的光束(一般地称为锥形束)发射X射线的X射线源和包括用于在一个时间采集对象的一个投影图像(一般地称为锥形束投影)的探测器元件的二维(2D)阵列的探测器。更具体地，本发明涉及一种用于从通过这样的成像系统所采集的锥形束数据中重建三维图像数据组的新方法。

背景技术

在扇形束计算机断层造影(CT)系统中，X射线源投影扇形束，该扇形束被准直以位于笛卡尔坐标系的称为成像平面的xy平面内。X射线束通过被成像的对象(诸如医学患者)并且击中放射线探测器的一维阵列。每个探测器产生一个电信号，该信号是X射线束被对象衰减的度量。在常规的CT系统中的源和探测器阵列在成像平面内在机架上围绕对象180-360°的角度旋转，其中在X射线源和探测器的一个旋转期间获取一组在不同角度方向拍摄的照片。从获取的衰减数据中，可以重建通过对象的xy平面的2D图像。在扇形束CT中执行图像重建的通常方法例如是滤波的反投影。

当X射线束是锥形的并且探测器相应地是探测器元件的2D阵列时，图像重建在数学上变得更复杂。然而，这些几何形状通常被用在介入的或血管造影的X射线系统中，例如C形臂系统，在该C形臂系统中将源和探测器安装于被构造为围绕患者旋转的C形臂的末端上。另一方面，非常期望能够通过围绕患者旋转C形臂一次来获得在患者内部体积的三维图像(在其间获取一组X射线投影的、例如大约180-360°的一次旋转被称为3D扫描)。因此，从圆形的锥形束数据进行3D图像重建已经成为近几十年来活跃的研究领域。在L.A.Feldkamp，L.C.Davis and J.W.Cress，“Practical cone beam algorithm”J.Opt.Soc.Am.A1，612-619(1984)中公开了用于基于圆形的源轨迹的实际的解决方案。Feldkamp方法是一种分析的重建方法。不幸的是，对于每个新的采集几何形状，必须调整分析的重建方法。换言之，Feldkamp方法不能被用于非圆形源轨迹。

近年来，可以使用具有其它采集几何形状的X射线系统，即，其中X射线源和探测器被安装在伸缩的臂上的系统。然而，这些系统没有被用来重建3D图像，因为在这样的X射线系统中，X射线源在一个3D扫描期间将沿着基于平面多边形轨迹移动。不幸的是，没有用于利用非圆形轨迹获得的锥形束投影的可用的重建算法。

发明内容

因此，本发明要解决的技术问题是，提供一种从沿着基于平面多边形的X射线源轨迹获得的锥形束投影来重建3D图像数据组的有吸引力的分析的方法。

该技术问题是通过按照本发明的方法、成像系统和计算机程序产品来解决的。

本发明方法允许从沿着可变半径的平面源轨迹获得的X射线投影中实现类似CT的重建。对于一个给定的几何形状，该类似Feldkamp的重建算法产生在扫描的平面中的精确的结果、在3D视野(FOV)的其余部分内的一般近似的结果以及该FOV的形状的精确的恢复。

探测器元件的二维阵列优选地布置在一个平面中，由此得到平的平板探测器。

在采集期间，X射线源沿着围绕待成像的对象、优选是物体或患者的大体上平面的轨迹移动。“大体上平面的”意味着由于成像系统中结构上的可变性导致的与平面的偏离是允许的，只要重建方法仍产生可接收的结果。在实践中，轨迹可以位于离轨迹的中央平面大约±15mm、优选在大约±10mm或±8mm之内。

然而，本发明方法的主要优点在于，轨迹不一定要是圆形的，或实际上不需要具有在平面内的特别形状。因此，是可变半径的轨迹，仅仅通过由在权利要求中被称为计数器参数的标量λ连续地编号的、空间中的一系列点(源点)来描述。特别地，轨迹可以是非凸的(none-convex)，即，其可以向内和向外弯曲。凸的轨迹不能关于等角点(isocenter)向外弯曲。换言之，轨迹可以具有非凸形状，其中可以通过一条线多于两次地穿过轨迹。在凸的轨迹中λ可以被认为是源极角，但这并不是必须的。优选地，在源和探测器之间的距离沿着轨迹变化，其中，该变化优选为显著，即，不是具有某些容差的基本为圆形的轨迹，而是距离优选变化大约至少±5mm、优选±10mm并且最优选至少±15mm。换言之，当源围绕对象移动、即沿着轨迹移动时，扫描半径R、即在源和扫描的等角点之间的距离是可变的。同样，该变化优选是显著的，即，至少±5mm、优选±10mm，其中非圆形扫描半径是优选的。按照一种优选实施方式，轨迹的至少一部分具有多边形的形状，即，其由有限数量(优选2-30条)的直线组成。

当源沿着轨迹移动时，通过探测器采集锥形数据的一系列二维(2D)阵列。这些2D阵列的每个是对象的投影图像并且以下也将被称为锥形束(CB)投影。在一次扫描期间(即源扫过轨迹一次)所采集的一系列CB投影被称为锥形束数据。如对于计算机断层造影来说通常的，轨迹优选张开至少一个180°的围绕对象的角度扫描间隔(angular scan interval)加上锥形束角度，例如大约200°-300°。

成像系统可以是任何X射线系统，但是特别可以是一个其中X射线源和探测器被安装于伸缩的臂和/或轨道和/或机器人臂或其组合上的X射线系统。

以下将详细并且以数学表达描述重建算法。然而，要注意，使用在权利要求1中定义的步骤的方法的任何其它实现也将落入本专利申请的范围。

在一个扫描期间X射线源沿着其移动的轨迹可以在数学上表达为

a(λ)＝R(λ)(cosλ，sinλ，0)其中λ∈Λ                    (1)

标量λ是对轨迹上的、在其上采集锥形束投影的点的计数器参数，并且如上所述，可以被认为是源极角。Λ可以被认为是角度扫描间隔并且R(λ)描述在源和对象坐标系的原点(0，0，0)之间的距离。该几何结构在图2中示出。因此，参数λ被用来描述沿着源轨迹的位置。探测器10优选是平板探测器并且u和v分别是在水平和垂直方向上测量的探测器坐标。由此，探测器平行于e_u(λ)＝(-sinλ，cosλ，0)和e_v(λ)＝(0，0，1)并且垂直于e_w(λ)＝(cosλ，sinλ0)。探测器-源距离在扫描期间允许非常平滑地变化并且表示为D(λ)。

尽管在图2中，探测器的平面垂直于与对象坐标系的原点和a(λ)二者相交的线，对于该算法来说这并不是必须的。探测器还可以关于中央的射线倾斜。当然，优选不应该倾斜太多，因为视野只能覆盖对象的被投影到探测器上的区域。

对象20在图2中用圆表示并且具有对象密度函数f(x)。

由探测器10采集的锥形束投影用函数g(λ，u，v)表示，使得值g(λ，u，v)相应于沿着连接在a(λ)处的源位置到相应的探测器上的点(u，v)的线的对象密度积分。该线的单位方向在图2中用矢量α表示。

锥形束投影在数学上可以表达为

$g(\lambda ,u,v)={\int }_o^{\infty }\mathrm{dtf}(\underset{‾}{\alpha }\left(\lambda \right)+t\underset{‾}{\alpha }(\lambda ,u,v))---\left(2\right)$

其中，u∈[-u_m，u，_m]和v∈[-v_m，v_m]是沿着e_u(λ)和e_v(λ)测量的坐标，并且其中

$\underset{‾}{\alpha }(\lambda ,u,v)=\frac{u{\underset{‾}{e}}_u\left(\lambda \right)+v{\underset{‾}{e}}_v\left(\lambda \right)-D\left(\lambda \right){\underset{‾}{e}}_w\left(\lambda \right)}{\sqrt{u^2+v^2+D{\left(\lambda \right)}^2}}---\left(3\right)$

在步骤A，计算锥形束数据的导数。优选地关于源轨迹在固定的射线方向α上的计数器参数λ对数据微分，例如使用在文章“A new scheme for view-dependent data differentiation in fan-beam and cone-beam computed tomography”by F.Noo，S.Hoppe，F.Dennerlein，G.Lauritsch and J.Hornegger，Phys.Med.Biol.52(2007)5393-5414中描述的微分方案。在一种实施方式中，关于所有三个空间方向或角度，即视角λ、射线角度u和锥角v(垂直分量)对锥形束数据微分。由于u和v不是角度而是探测器坐标，所以其相应于这样的角度。换言之，在该实施方式中，进行全微分，也包括在扇形束数据中不重要、但是仅在锥形束数据中的垂直分量。这具有如下效果：微分(参见公式5)具有三项。

数学上，如果都满足：|u|≤u_m和|v|≤v_m，则微分可以被表达为

$g_D(\lambda ,u,v)=\frac{\partial }{\partial \mu }g(\mu ,\underset{‾}{\alpha }(\lambda ,u,v)){|}_{\mu =\lambda }.---\left(4\right)$

然而，如果|u|＞u_m或者|v|＞v_m，则设置g_F(λ，u，v)＝○，其中符号○被用来定义在探测器边界外部的数据。选择○具有特征(i)○+x＝○and(ii)○·x＝○，使得其通过计算传播并且在重建中将任何点x设置为○，其至少对于一个λ∈Λ投影在探测器边界外。通过这样做，算法自动地恢复FOV，该FOV具有在可变的半径的几何形状中的一般不是不重要的形状。在实际的实现中，○可以被置为浮点值：not a number(NaNq)，其满足上面描述的要求。

等式(4)可以被展开为

$g_D(\lambda ,u,v)=(\frac{\partial }{\partial \lambda }-\frac{\mathrm{uv}}{D{\left(\lambda \right)}^2}\frac{\partial }{\partial v}+\frac{D{\left(\lambda \right)}^2+u^2}{D\left(\lambda \right)}\frac{\partial }{\partial u}---\left(5\right)$

$+\frac{u(\partial /\partial \lambda )D\left(\lambda \right)}{D\left(\lambda \right)}\frac{\partial }{\partial u})g(\lambda ,\underset{‾}{\alpha }(\lambda ,u,v))$

与在圆形几何形状中的相应等式相比，该等式包含一个附加的项，即，在(5)的第二行的项，其被称为NEW TERM(新项)并且以下将被讨论。

在步骤A之后，在步骤B中利用类似Hilbert的滤波器对锥形束数据的导数进行滤波以产生经滤波的锥形束数据。该经滤波的锥形束数据例如可以是

$g_F(\lambda ,u,v)={\int }_{-u_m}^{u_m}\mathrm{du}\prime h_H(u-u\prime )\frac{D\left(\lambda \right)g_D(\lambda ,u\prime ,v)}{\sqrt{D{\left(\lambda \right)}^2+{u\prime }^2}}---\left(6\right)$

其中，h_H(u)是Hilbert变换的核。

因此，锥形束数据的导数优选仅在行方向(即在u方向)上与核函数卷积(convoluted)，得到类似Hilbert的变换。核h_H(u)优选是1/πu，如在经典的Hilbert变换中那样，但是还可以使用具有不同的变迹法(apodization)或具有不同的带宽限制的修改的Hilbert滤波器。因此，使用术语类似Hilbert的滤波器。

在步骤A之后的某点上，微分结果被乘以长度校正权重、也称为余弦权重。在上面的公式(6)中，该长度校正权重表达为

在下一步骤C中，数据被乘以冗余的加权函数。冗余处理是必要的，以考虑如下事实：通过对象的某些射线路径可能在数据中被包括两次，即，如果该特定的路径已经从相对侧被辐射的话。

按照一种优选实施方式，冗余等待函数被用于特定的轨迹。这允许使用任何可能的轨迹，也包括非凸轨迹。

由此，优选地对于每个轨迹重新计算等待函数。对于非凸轨迹，冗余等待不是不重要的，并且不能求助于已知函数、诸如Parker等待函数。仅对于凸轨迹，λ也具有极角的含义。

在实践中，锥形束数据包含少数冗余数据，因为只有在轨迹的平面中的射线可以是冗余的。然而，在一种实施方式中，算法假定，这样的伪冗余(pseudo-redundant)射线实际上包含与近似相同的信息。

对于圆形轨迹，已知使用所谓的Parker加权函数，其对仅被测量一次的每个射线积分给出1的贡献，并且平衡冗余数据样本的贡献，使用三角函数的特性，从而相应于对于重建考虑了两次的数据样本的权重总计为单数(unity)。最初的Parker加权函数在D.L.Parker“Optimal short scan convolution reconstruction for fan-beam CT”Med.Phys.9(2)：254-257，1982中描述。

对于本发明，优选类似于在Feldkamp方法中使用的过程，在一个给定的经滤波的锥形束投影内的位于相同的列中(即具有相同的u坐标)的所有的点接收相同的权重，表示为w(λ，u)。

按照一种实施方式，用于某一列的权重如下计算：识别在轨迹上所有补充的源点，即，源轨迹与一个垂直于轨迹的平面并且平行于在两个点上(即在a(λ)处和通过探测器列u)的X射线路径的平面相交处的点。然后使用辅助函数c(λ)并且计算冗余权重为

w(λ，u)＝c(λ)/(c(λ)+c_c(λ，u))，                            (7)

其中，c_c(λ，u)是在被确定为是补充的源点处所分析的辅助函数的总和。辅助函数c(λ)原则上可以是由用户预先确定的任意函数，仅有的限制是：对于所有相关的λ，c≥0。c(λ)例如可以是常数，或者可以在边沿处平滑下降。在上面提到的公式中，确定的是，所有的权重相加为1。优选地，在上面的公式中，可以获得类似Parker的冗余加权函数。

优选地对于每个不同的轨迹重新计算函数w(λ，u)，但是可以总是使用相同的辅助函数。

可以在步骤B之后执行冗余处理步骤C。在这种情况下，对经滤波的锥形束数据应用加权函数。替换地，可以在步骤A之前执行冗余处理步骤C。在这种情况下，加权函数乘以最初(采集的)锥形束数据。

在步骤A-C之后，不取决于其顺序地执行最后的重建步骤D、反投影。在步骤D中，经滤波和冗余加权的锥形束数据被反投影到图像体积以计算3D图像数据组其给出对对象密度f(x)的良好近似。优选地，以下公式被用于CB数据，该CB数据是在u上非删简的：

$\hat{f}\left(\underset{‾}{x}\right)=\frac{1}{2\pi }\underset{\Lambda }{\int }\mathrm{d\lambda }\frac{w(\lambda ,u^{*})}{R\left(\lambda \right)-\underset{‾}{x}·{\underset{‾}{e}}_w\left(\lambda \right)}{g}_F(\lambda ,u^{*}{,v}^{*}).---\left(8\right)$

该公式描述了对于所有的λ∈Λ，经滤波的CB数据的加权的反投影g_F(λ，u，v)，其中w(λ，u)是类似Parker冗余加权函数并且u^*和v^*是对于给定的λ，投影x的探测器坐标。

优选地，使用反投影加权，其类似于Feldkamp方法的反投影加权，但是仅取决于在体素和源点之间的距离，而不是平方距离。优选地，在3D图像数据组中累加反投影值，以连续地计算对象密度函数，其中例如一个接一个地使用经处理的锥形束投影。

由于重建算法可以处理非圆形轨迹，因此优选地，源轨迹是非圆形的。在实践中，这不是仅包括在C形臂系统中可能发生的与圆形路径的偏差(例如在圆形轨迹的±2mm内)，而是意味真的非圆形几何形状。例如，在源和探测器之间的距离可以沿着轨迹显著地改变，即，＞2mm，优选＞5mm或＞10mm。优选地，变化是平滑的。优选地，源轨迹的至少部分不能由圆的片段近似。

本发明方法优选被应用于其中源轨迹的至少部分可以通过一个或几个直线来近似的几何形状。这样的轨迹例如可以利用这样的X射线系统实现：在该X射线系统中源和/或探测器被安装用于平移运动，例如，在轨道和/或在伸缩臂上。替换地，轨迹是这样的：这样的伸缩臂的尖端按照直线移动，优选在其伸缩方向上移动。

按照一种优选实施方式，轨迹可以通过多边形的至少一部分来近似，例如由3至10条直线组成的多边形，其张开一个180°至360°、优选200°至300°的极角Λ。

探测器优选地沿着与X射线源相对的相应路径移动，从而在沿着源轨迹的每个点处，以尽可能小的删简(truncation)，可以采集来自对象的锥形束投影。然而，重建的方法甚至能处理删简的数据，例如，当数据外推方法被用来估计在探测器边界外的投影数据。一种可能的外推方法是使得数据值在探测器边界外平滑地下降。

替换地，探测器还可以在采集一系列锥形束投影期间位于相同位置。

本发明可用于在视野内部比较整个对象的情况、或者待检查的对象是未完全包含在视野内部的大对象的一部分的情况。在后一情况下，例如通过锥形束投影数据g(λ，u，v)的平滑下降，可以外推在视野外部的区域。

本发明还涉及一种能够产生待成像的对象的3D图像数据组的成像系统，包括X射线源和探测器，它们优选被构造为沿着上面描述的轨迹移动。

成像系统还包括采集系统，当源沿着轨迹围绕对象移动时，采集系统用于采集一系列二维阵列锥形束数据。因此，该采集例如可以包括对于X射线源和探测器的安装，允许其围绕对象移动，例如在轨道上或伸缩臂上。此外，采集系统可以包括用于控制源和探测器的运动的控制单元。

此外，成像系统包括图像重建系统，其被用于从锥形束数据按照上面描述的方法重建3D图像。该图像重建系统可以在计算机、例如PC或包括CPU和工作存储器(例如RAM)的工作站上实现。

替换地，成像系统可以是基于机器人的C形臂系统或任意其它X射线系统，其中对于不同的投影角度，在源和探测器之间的距离可以改变。

本发明的一个应用是医学成像。然而，其还可以应用于工业领域，特别是(例如金属或其它材料的工件的)非损伤性测试和其它应用。

本发明最后涉及一种计算机程序产品，包含计算机可读的软件代码部分，当程序在计算机上运行时，软件代码部分被用来使得计算机执行以下步骤：当以锥形束发射光子射线的X射线源沿着基本为平面的轨迹围绕待成像的对象移动时，该轨迹由通过计数器参数在空间上连续地编号的一系列点所描述，访问由包括探测器元件的二维阵列的探测器所采集的一系列二维阵列锥形束数据；通过上面描述的方法从锥形束数据中重建3D图像。优选地，计算机程序产品可以直接加载到数字计算机的内部存储器中。本发明还涉及一种计算机程序本身，并且涉及一种计算机可用的介质或数字可读的介质，在其上存储了这样的计算机程序。

附图说明

以下参考附图通过优选实施方式进一步解释本发明。附图中，

图1是在其上可以实现本发明的成像系统的示意图；

图2是可变半径轨迹和对象的顶视图；

图3是本发明方法的一种实施方式的流程图；

图4是重建方法的第一实施方式的流程图；

图5是重建方法的第二实施方式的流程图；

图6-8是按照[0，100]Houndsfield单位(HU)在z＝0mm，z＝80mm和x＝5mm处通过头部模体重建的片；

图9-11是新项对在HU-刻度上的宽为100的窗中在z＝0mm，z＝80mm和x＝5mm的切片处的最后重建的贡献；

图12示出了在按照短扫描直角几何形状的扫描期间在X射线源和探测器中心的示例性轨迹。

具体实施方式

图1是在其上可以实现本发明的成像系统的示意图。患者22可以被置于患者卧榻6上。待成像的对象可以是患者内部的一个区域20。安装X射线源4使得可以通过伸缩臂30和轨道32在几个方向上平移运动，沿着轨道32，X射线源可以在滑块33处滑动。此外，X射线源安装在结合点34上，从而其可以倾斜以使得其发射的射线的锥形束将穿过对象20并且然后击中探测器10。

探测器10优选是平板探测器并且本身可倾斜地通过结合点38安装。结合点38在伸缩臂36上可以被平移地移动。

应当注意，关于成像系统的结构可以有无数变形，特别是关于X射线源4和探测器10的安装。例如，源4可以被安装在房间的天花板上的可能具有伸缩臂的轨道上。此外，如下的配置也是可能的：只有源沿着伸缩臂移动，而探测器围绕对象倾斜。优选地，源和/或探测器分别在至少两个方向上是可移动的或者具有至少两个、优选三个或四个自由度。

可以液压地或者通过电机移动伸缩臂30、36。

图2是可变半径的几何的示意图，其中，轨迹a(λ)位于纸平面中，探测器10与其垂直。待成像的对象20围绕坐标系的原点定位。

图3是显示该方法的实施方式的主要步骤的流程图。在步骤100中，在一次扫描中采集一系列锥形束投影，其中X射线源沿着平面可变半径轨迹围绕对象运动。在一次扫描中，典型地采集大约50-300个、优选150-250个锥形束投影。

在本身包括步骤A-D的步骤102中，将在步骤100中采集的锥形束数据重建为3D图像数据组。可选地，在步骤104中，这样的图像可以被显示在显示屏上或者被打印在纸上，从而图像数据组可以被观察并且可以被进一步处理。如现有技术中公知的，这样的3D图像数据组可以被用来在该过程(session)的后面阶段期间，将3D数据与在相同的成像系统上采集的实时2D投影相关联，以便例如控制介入。

图4示出了重建方法的第一实施例，该方法包括在步骤A中对锥形束数据微分、在步骤B中优选利用Hilbert滤波器滤波并且与长度校正权重相乘。在步骤C中，将步骤B的结果与冗余加权函数相乘。最后，在步骤D中，处理后的锥形束数据被反投影到3D体积以产生3D图像数据组。

按照图5的第二实施例，还可以首先执行步骤C，即冗余处理，随后进行在步骤A中的微分、在步骤B中的滤波和步骤D中的反投影。

例子

在如图6所示的近似直角短扫描几何形状中测试了本发明的方法。在40处标出了在扫描期间X射线源的轨迹，并且在42处标出了探测器的轨迹。每个粗略沿着直角轨迹移动，其中边角被拉平。

在该几何形状中，轨迹半径和源-探测器距离被参数化为

R(λ)＝R_ml(λ)和D(λ)＝D_ml(λ)

其中，R_m＝750mm和D_m＝1200mm分别是最小半径和最小源-探测器距离，并且其中

$l\left(\lambda \right)=\min ({\left|\cos \lambda \right|}^{-1};{\left|\sin \lambda \right|}^{-1};\sqrt{2}-ϵ)$

是调整该几何形状的比例函数。在该定义中，标量ε被用来避免在轨迹中的90°边角并且在我们的分析中被置为0.07，得到图6的轨迹。

使用恒定增量Δλ＝0.5°并假定探测器尺寸U_m＝20mm和U_m＝15cm以及0.5mm大小的方形像素，在区间λ＝[0°，220°]上模拟FORBILD头部模体的CB投影。使用平滑的类似Parker冗余加权函数w(λ，u)执行重建。

在图7-12中展示了结果。这些图提供通过重建的结果的两个轴向和一个冠状切片，其中在视野外部的区域被显示为黑色(左手侧)或灰色(右手侧)。如特别在图7中可以看出的，在z＝0mm处的切片，FORBILD头部模体的良好重建是可能的。图10-12示出了与上面提到的NEW TERM相关的贡献，其出现以便校正在重建结果中的梯度。

本发明提供一种用于来自特别是具有平滑改变的半径的平面源轨迹的平板锥形束重建的、实用的类似Feldkamp重建方法。连同重建一起，该方法还恢复在所考虑的几何形状中可以采取复杂形状的视野。