CN1658795A - 用于快速发散波束断层摄影术的方法和装置 - Google Patents

Abstract

提供了一种用于发散波束反向投影的快速方法，用于从预处理的发散波束窦腔造影(22)生成电子图像，该窦腔造影是发散波束投影的集合。该方法能够包括以下步骤：将窦腔造影(22)细分(1010、1012、1014、1016)为多个子窦腔造影(g₁、g₂、g₃、g₄)(1010、1012、1014、1016)；执行所述子窦腔造影(1018、1020、1022、1024)的加权的反向投影，以生成多个相应的子图像(f₁、f₂、f₃、f₄)；以及会聚(1026)所述子图像(f₁、f₂、f₃、f₄)；以创建电子图像(f，26)。从窦腔造影到子窦腔造影的细分能够以递归方式执行。所提供方法的双重目的是提供一种用于再投影电子图像的装置，例如，从该图像生成发散波束窦腔造影。

Description

用于快速发散波束断层摄影术的方法和装置

美国专利Nos.6263096、6282257、6307911、6322035和6351548的内容全部引用于此作为参考。

技术领域

本发明涉及用于发散波束断层摄影术的方法和装置，特别涉及降低了计算的复杂度而没有在视觉上可察觉的衰减或严重的数值不精确的、用于发散波束断层摄影术的方法和装置。

背景技术

断层成像再现是几乎是所有诊断成像形式的基础的一种众所周知的技术，包括x-射线计算断层摄影术(CT)、正电子放射断层摄影术(PET)、单光子发射断层摄影术(SPECT)和用于磁共振成像(MRI)的某些获取方法。它还广泛用在制造业，用于非破坏可靠性评估(NDE)，而最近还被用于机场行李的检查。当前CT扫描器的主要数据获取模式是通常与螺旋形(helical)或螺线形(spiral)扫描相结合的扇形波束几何形状。

为了更快速的图像获取，CT扫描器技术正在移位到使用区域检测器的锥形波束获取。所生产的许多CT扫描器已经具有多行检测器(一般是16行)，从而导致了(低角度)锥形波束获取。用于广角度锥形波束获取的区域检测器将在未来几年中被引入到医学CT扫描器中。当前存在对用于检测隐藏武器和爆炸器件的高速3D断层摄影成像的仍未被满足的需求。在诸如PEY或SPECT的其它形式中，锥形波束获取已经是主要的形式。

在扇形波束和锥形波束CT中的再现问题是从一组沿着分别来自在指定获取轨线或轨道上传播的源的诸如扇形或锥形发散的线的线积分投影中恢复一目标的2D或3D图像。因此，扇形波束和锥形波束几何形状是发散波束结构的实例。

用于断层成像再现的选择的方法是滤波反向投影(FBP)或卷积反向投影(CBP)，它们都是用加权的反向投影步骤。这个步骤是在技术上的计算瓶颈，计算要求标度(scaling)对2D的N×N像素图像为N³，而对3D的N×N×N体元图像则至少为N⁴。因此，图像分辨率从N倍增到2N将导致计算量增加大约8倍(或在3D情况下增加16倍)。虽然计算机已经变的更加快速，能够实时收集更大量的数据(例如，具有多行检测器的心脏成像，干涉成像)的新技术的出现，以及向3D锥形波束几何形状的移位，但是仍然存在对快速再现技术的不断增长的需要。快速再现能够加快图像形成过程，降低专用图像再现计算机(见下)的成本，或同时达成两者。

反向投影的双重(duel)操作是再投影，它是计算电存储图像的投影的处理。这个处理在断层成像再现中也起着基础的作用。反向投影和再投影的组合还能够被用来构造用于螺线锥形波束几何形状中的长国标问题的快速再现算法，其是人类对象的实际3D成像的关键。而且，在各种应用中，最好甚至必须使用迭代再现算法，其中，为了单一图像的再现，反向投影和再投影步骤都被执行多次。加快反向投影和再投影步骤将决定这样的迭代方法的经济可行性。

在讨论快速方法中，区别断层成像数据的两种主要形式是很重要的：(i)平行波束；和(ii)发散波束及其特殊情况。虽然平行波束形式最方便用于全部理论的和计算的操作，但是发散波束是在商业医学和工业CT扫描器中最常见到的形式。尽管最初是2D方法，扇形波束再现是用于螺旋和多层面螺旋3D(测定体积的)再现的当前技术发展水平的方法中的关键组分。可以预料，新的和未来的设计将会移位到锥形波束几何形状。因此，不论是在2D还是在3D中，发散波束几何形状是主要的成像形式，并且在可预见的将来也很可能保持这样。

由于到目标的很大源距离的限制，发散波束几何形状减到平行波束几何形状。从而，发散波束几何形状的所有数据处理方法(包括本发明的方法)也都适用于平行波束几何形状。然而，反之则不然：在实践中对于很小或中等源距离，两种几何形状非常不同，以致当直接应用到发散波束数据时，平行波束处理方法生成低等的或不能令人满意的结果。因此，发散波束几何形状需要不同于平行波束几何形状的处理。

本发明致力于发散波束断层摄影术。一些再现方法依赖于称为数据再装箱(rebinning)的处理，以便将发散波束数据再排列(或插入)成为平行波束数据，然后用平行波束算法进行处理。3D发散波束几何形状的其它方法依赖于将锥形波束数据变换为3D拉登(Radon)变换数据(或其衍生物)。但是，所定义的自然发散波束方法(本发明为快速版本)是直接对发散波束数据进行处理，而没有预先再装箱为平行波束数据，或将数据变换为3D拉登变换域。在自然发散波束方法中的再现处理包括预处理投影(例如，加权和滤波)，接着进行加权的发散波束反向投影操作，以及或许进行发散波束再投影操作，或者一系列这样的操作。这样的方法因此称为发散波束滤波反向投影(DB-FBP)算法。正如平行波束情况一样，它们的计算通常由反向投影和再投影操作支配。

发散波束再投影能够通过以下步骤执行：首先执行平行波束再投影，然后对结果进行数据再装箱，但是这里也最好使用不进行数据再装箱的自然(native)发散波束再投影算法。

数据再装箱的缺点包括限制加速的可能的假象和增加的计算量。数据再装箱还需要在处理能够开始之前对大部分数据的获取和处理，这再一次限制了速度。依赖于到3D拉登域的变换的方法具有类似的缺点。因此，需要能够克服这些缺点的用于发散波束断层摄影术的方法和装置，它们具有高度的灵活性，并与使用传统发散波束反向投影的方法相比提供很大的性能增益。

专用硬件是加速反向投影处理的传统方法。使用定制芯片的特殊类型的计算机、或多个处理器、或两者的组合被设计来执行需要的计算。

这个方案的缺点包括多个处理器或定制硬件的成本，以及随着通用计算机的性能的高速增长，所需的特殊用途的体系结构很快变得陈旧的事实。因此，存在不需要专用硬件就能够对发散波束几何形状进行快速处理的需求，并且这种处理能够容易地在标准串行或并行体系结构上实现，使得它们有更好的成本效率。

发明内容

因而，本发明的一个目的是提供新的和改进的用于发送波束断层摄影术的方法和装置。

另一个目的是提供新的和改进的用于发散波束断层摄影术的方法和装置，其降低了计算的复杂度而没有在视觉上可察觉的衰减或严重的数值不精确。

与本发明的一个方面一致，提出了一种方法，用于根据适合于反向投影的预先处理的发散波束窦腔造影生成电子图像，所述窦腔造影是发散波束投影的集合。该方法能够包括下面步骤：将窦腔造影细划分为多个子窦腔造影；在固定到目标的整体坐标系中执行子窦腔造影的加权的反向投影，以便在整体坐标系中的适当位置处生成多个对应的子图像；和会聚子图像以创建所述电子图像。从窦腔造影到子窦腔造影的细划分能够以递归方式执行，直到每个子窦腔造影都表现为期望尺寸的子图像为止，其中细划分步骤包括期望数量的精确细划分和期望数量的近似细划分。

与本发明的另一个方面一致，提出一种用于再投影电子图像、即，从该图像生成发散波束窦腔造影的方法。该方法能够包括下面步骤：将图像划分为多个子图像；在整体坐标系中计算每个子图像的子窦腔造影；和会聚所述子窦腔造影以创建窦腔造影。图像的细划分能够以递归方式执行，直到每个子图像都具有期望的尺寸、并且子窦腔造影的计算能够在期望数量的递归级次中是近似的，而在剩下的递归级次中是精确的为止。

与传统反向投影和再投影算法比较，这些方法提供了类似于FFT的加速。节约了基本的计算量，而没有可察觉的衰减或严重的数值不精确的损失。

附图说明

结合附图参考下面对发明的实施例的描述，本发明上述和其它特征以及获得它们的方式将变得更加清楚，发明本身也将被最好地理解，其中：

图1是本发明使用的装置的框图；

图2是平面等间距锥形波束几何形状的图形；

图3是平面等间距螺旋锥形波束几何性状的图形；

图4是说明具有圆形源轨道和平面等间距检测器的锥形波束断层摄影术的图形；

图5是说明共线性等间距扇形波束几何形状的图形；

图6是说明具有圆形源轨道和圆柱形检测器的锥形波束断层摄影术的图形；

图7是等角扇形波束几何形状的图形；

图8是用于M×M子图像的共线性等间距扇形波束投影几何形状的图形；

图9是说明将一个图像分割为多个子图像的图形；

图10是说明扇形波束反向投影的精确分解的图形；

图11是由窦腔造影分解算子

执行的处理的流程图；

图12是扇形波束反向投影的近似分解的图形；

图13是二级二重图像分割的图形；

图14是扇形波束反向投影的二级近似分解的图形；

图15是扇形波束反向投影的二级混合(一个精确、一个近似)分解的图形；

图16是用于快速扇形波束分级反向投影的伪码；

图17是说明扇形波束再投影的精确分解的图形；

图18是由投影增大算子

执行的处理的流程图；

图19是扇形波束再投影的近似分解的图形；

图20是扇形波束再投影的二级近似分解的图形；

图21是扇形波束再投影的二级混合(一个精确、一个近似)分解的图形；

图22是用于子图像的等角扇形波束投影几何形状的图形；

图23是用于子图像的平面等间距锥形波束投影几何形状的图形；

图24是将一个3D图像分割为多个子图像的图形；

图25是将3D图像各向异性地分割为子图像的图形；

图26是3D图像的二级非均匀递归分割的图形；

图27(a)-27(d)是比较用于共线性等间距检测器的精确和快速扇形波束反向投影的图；

图28(a)-28(d)是比较用于等角检测器的精确和快速扇形波束反向投影的图；

图29(a)-29(b)是通过图28(a)-28(d)的反向投影的片段取得的曲线图；

图30(a)-30(b)是点扩展函数的比较；

图31(a)-31(b)是通过点扩展函数的片段的曲线图的比较；

图32是3D锥形波束再现的图的集合；

图33是3D锥形波束再现的图的另一个集合；

图34是通过在图32中的图像的片段的曲线图的集合；以及

图35是通过在图33中的图像的片段的曲线图的集合。

具体实施方式

本发明可以应用于各种成像装置，包括CT扫描器。典型的像装置10(图1)包括从诸如头部的目标获得数据、并将相应于发散波束投影14的原始数据发送给投影预处理器16的扫描器12。投影预处理器16对数据施加各种变换、规格化和校正，以及可以被移位变化的加权和滤波。投影预处理器16的输出是发散波束窦腔造影₁18，其被馈送到窦腔造影更新处理器20。窦腔造影更新处理器20可能使用来自发散波束窦腔造影₂34的信息修改输入窦腔造影₁18，例如校正包括波束硬化现象的各种假象，或者作为多级的或迭代的再现处理的一部分。

窦腔造影更新处理器20的输出是被输入到快速反向投影处理器24的发散波束窦腔造影₃22。快速反向投影处理器24一般是任何适当类型的、被编程以执行这里描述的精确和近似反向投影算法的计算机或专用硬件，或它们的任意组合。

快速反向投影处理器24的输出是被输入到图像调整处理器28的电子图像₁26。图像调整处理器28执行电子图像所需的后处理，可能包括假象图像的识别和提取，或者在将来的多级的或迭代的再现过程中处理的图像。

如果希望，图像调整处理器28能够生成将被馈送到快速再投影处理器32的电子图像₂30。快速再投影处理器32一般是任何适当类型的、被编程以执行这里描述的精确和近似再投影算法的计算机或特殊用途的硬件，或它们的任意组合。如果希望，这个处理器能够与反向投影处理器24共享使用相同的计算机和硬件。

快速再投影处理器32的输出是被反馈到窦腔造影更新处理器20的发散波束窦腔造影₂34。所述反向投影/再投影处理能够继续，直到获得适当的结果为止。虽然再投影并不是总需要，但是它在许多情况下是有帮助的。

当电子图像₁26是适当的时，图像调整处理器28生成被馈送到存储/分析/显示设备38的电子图像₃36。意图是希望电子图像₃36能够被存储在计算机存储器中，和/或对例如异常或危险物质进行电子分析，和/或显示，和/或以某种可视的形式打印。

本发明的快速发散波束反向投影算法依赖用于在美国专利No.6287257中描述的平行波束断层摄影术的快速分级反向投影算法(FHBP)。同样的，本发明的快速发散波束再投影算法依赖用于在美国专利Nos.6263096和6351548中描述的平行波束断层摄影术的快速分级再投影算法(FHRP)。反向投影和再投影紧密相关，并且它们的快速算法也是得自同一个原理。在这里将更详细地描述更普遍使用的快速反向投影。

为了理解平行波束FHBP算法，考虑一个附有整体坐标系的N×N图像，所述坐标系的原点位于图像的中心。考虑分解这个N×N大的图像为四个较小的N/2×N/2尺寸的图像，其称为子图像，每一个都位于整体坐标系的不同的象限。拉登变换的两个性质被用在FHBP算法的导出中。(i)蝴蝶结(bow-tie)性质，即对于以整体坐标系的原点为中心的半幅图像，对角变量的投影的频谱支持也减少了一半。这意味着有居中的半幅图像能够从包括一半投影数量的窦腔造影中再现。(ii)移位性质，即被移位图像的平行投影对应于原始图像的被移位的平行投影。

现在假设包括p滤波后投影的窦腔造影能够被用于再现整个图像f_o为了再现多个子图像f_i之一，P投影首先被截断以支持子图像的投影，然后被移位以对应于子图像的居中版本，最后按照角方向分样成P/2投影。根据这些新的P/2投影再现子图像f_i的居中版本，然后该子图像被移位回到其在整体坐标系中的正确位置。对四个子图像的每一个执行这个处理，当把它们会聚在一起时就提供了整个图像f的再现。对四个子图像的再现需要总共4c(P/2)(N/2)²＝cPN²/2次操作，其中c是常数。这减少了原始再现成本的一半。递归地应用这种分解以在每一步骤减少图像尺寸的一半，从而允许使用一半数量的投影对半幅子图像进行格式化，使总的计算成本能够减少到O(N²logN)。

这个分解和会聚方案的双重操作被用在美国专利Nos.6263096和6351548中描述的快速分级平行波束再投影算法，即FHRP。

用于发散波束反向投影和再投影的本发明是基于分级分解图像、然后使用包含较少发散波束投影的发散波束子窦腔造影再现(或再投影)子图像的思想。然而，当试图将平行波束FHBP或FHRP的思想应用到发散波束情形时，就出现困难，即上述移位性质(ii)不再适用。也就是被移位的子图像的发散波束投影并不对应于原始子图像的被移位的发散波束投影。由于平行波束FHBP和FHRP使用这个性质，所以，平行波束方案必须针对发散波束情形进行修改。

本发明的算法对所述分解使用不同的处理，其不包括子图像相对于固定到目标、即全幅图像上的整体坐标系的中心的移位。与平行波束FHBP类似，新的发散波束反向投影算法具有两种类型的分解，一种精确但是缓慢，另一种近似但是快速。新算法的精确发散波束分解根本不必包括任何数据移位，它只包括到对应于所述子图像的子窦腔造影的窦腔造影的截断。这些子窦腔造影然后在整体坐标系中进行反向投影，从而在整体坐标系中的正确位置处生成相应的子图像。近似发散波束分解使用不同于在平行波束FHBP使用的处理，以便减少用于再现较小尺寸的子图像的投影的数量。这种新的处理只包括投影的截断和移位，而不包括子图像的截断和移位。这种新的处理还会包括在FHBP中没有用到的加权。类似的陈述也适用针对发散波束再投影的新的处理。

为了再现半幅子图像，适当截断的(加权的和滤波的)发散波束投影被移位一个由在整体坐标系中子图像中心投影位置所决定的距离，然后在源位置参数中进行二分样，并以相同量移位回去。所述子图像能够根据这个减少的投影组被非常近似地再现。对于在子图像上的精确和近似加权的反向投影，所有涉及图像坐标、包括加权的计算都是在参考实际尺寸图像的原点的绝对或整体图像坐标中完成的。以递归形式组合精确和近似分解，就得到发散波束几何形状的快速分级反向投影算法。就像它们的平行波束FHBP的对应部分，这些算法通过选择精确分解对近似分解的数量以及通过调整各种算法参数在计算量和精确度之间提供平衡。

发散波束成像几何形状

一般3D几何形状的发散波束获取几何形状包括在围绕成像目标的3D空间中的抛物线上移动的发散射线的顶点(通常是源)，以及检测器表面，在其上测量沿着经过成像目标的源射线的线积分。其它具有会聚而不是发散波束几何形状的成像形式，诸如在发射断层摄影术中出现的那些形式(SPECT和PET)，能够被变换到这种几何形状。

获取几何形状一般通过检测器的形状和相关联的射线采样方案，以及通过源抛物线或轨道来分类。最普通的两种检测器的类型是(i)具有等间距检测器的平面检测器表面；以及(ii)具有在方位角中射线采样等间距、以及在高度中线性等间距的圆柱形检测器。在利用表面检测器的3D成像中，发散波束几何形状也被称为锥形波束几何形状。

在发散波束2D获取的一个特殊情况中，检测器表面被减小到检测器的单线或弧线，从而导致所谓的共线性等间距和等角扇形波束几何形状。这些扇形波束几何形状被在当前的商业系统中广泛地开发。

在已经提出的各种源轨道中，最普通的是单圆形轨道和螺线(或螺旋)轨道。然而，诸如两个圆形、或一个圆形和一条线的其它轨道也有实际意义，并且这些各种轨道的任意混杂的版本在实践中也有出现。

本发明的处理、算法和装置可被应用于一般的发散波束几何形状，并且能够被专用于任何早先提到的特定几何形状，或专用于发散波束几何形状的其它实例。类似的，本发明的方法适用于任何数学地等同的数据组，包括但不局限于发散波束几何形状。为了说明的目的，该方法将对下面的情况进行一些详细的解释：

1.单圆形源轨道

(a)共线性等间距扇形波束几何形状

(b)等角扇形波束几何形状

(c)平面锥形波束几何形状

2.具有平面检测器的一般锥形波束几何形状

3.螺线源轨道

将对一般处理如何被专用于这些情况进行说明，从这些说明可以明白，用于一个检测器形状和源轨道的算法可以被很容易地修改成用于不同的检测器形状和源轨道。

具有任意轨道的发散波束成像几何形状

在图2中说明了检测器在平面表面上均匀地隔开的平面等间距锥形波束几何形状。在这种情况下，检测器平面允许对每个源位置具有任意地位置和方位。

发散射线的源S利用由参数θ(θ_min≤θ≤θ_max)规定的轨道上的源位置在围绕3D目标的轨道上传播，该3D目标由相对于固定到该目标的整体坐标系的中心O的源位置向量 $\stackrel{\→}{s}\left(\mathrm{\θ}\right)=[x\left(\mathrm{\θ}\right),y\left(\mathrm{\θ}\right),z\left(\mathrm{\θ}\right){]}^T$ 定义。对于每个θ，定义了附加到检测器的3D坐标系，以致原点O_θ是轨道点

在检测器平面上的正交投影。

正交向量 和

是在从O_θ指向源的检测器平面中选择的，而

是正交于检测器平面的单位向量。在检测器平面中的点 $\stackrel{\→}{t}=t_1{\hat{t}}_1+t_2{\hat{t}}_2$ 是按照它的坐标(t₁，t₂)定义的。在这个系统中，源轨道点被表示为

其中D(θ)是源和检测器之间的距离。(在SPECT中，D(θ)是锥形波束准直仪的焦距。)

用 $\stackrel{\→}{r}={[x,y,z]}^T$ 表示在3D图像(例如，目标)中的位置，而用 $\stackrel{\→}{\mathrm{\τ}}(\mathrm{\θ},\stackrel{\→}{r})={[{\mathrm{\τ}}_1(\mathrm{\θ},\stackrel{\→}{r}),\mathrm{\τ}(\mathrm{\θ},\stackrel{\→}{r})]}^T$ 表示与源射线通过点

的检测器平面的焦点的t₁，t₂位置(见图2)。目标f在源轨道位置θ和检测器位置(t₁，t₂)的发散波束投影(Pf)(θ，t₁，t₂)是沿着用参数(θ，t₁，t₂)表示的源射线的线积分，并由下式给出：

$\left(\mathrm{Pf}\right)(\mathrm{\θ},t_1,t_2)={\∫}_{\stackrel{\→}{r}\∈R^3}f\left(\stackrel{\→}{r}\right)\mathrm{\δ}[(t_1,t_2)-\stackrel{\→}{\mathrm{\τ}}(\mathrm{\θ},\stackrel{\→}{r})]d\stackrel{\→}{r}---\left(1\right)$

其中，

是狄拉克三角函数(Dirac delta function)。同样的获取几何形状能够被用来描述一种情形，其中，源射线向着源会聚而不是从其发散。从而，本发明的处理同样地适用于会聚波束几何形状。

在由沿着源轨道的参数θ_p＝pΔθ，p＝0，...P-1定义的P离散源位置处获得投影，但通常不必具有统一的间距Δθ＝(θ_max-θ_min)/P。(Pf)(θ_p，.)项(对于所有的 值)被称为在源位置θ_p的发散波束投影。通常利用在t₁和t₂轴上的不同间距T₁和T₂在统一的矩形栅格上对检测器平面进行采样，

不同的源位置处的发散波束投影集合将被称为发散波束窦腔造影。同一项也将被用于一组预处理的发散波束投影，或任何其它的数学等效的数据组。在窦腔造影中的投影的集合可以是不完全的，最少只包含一个单一的投影。

现在将描述一般发散波束几何形状的一些特殊情况。

具有平面检测器的锥形波束螺线断层摄影术

在这种情况下的源轨道是如图3所示围绕轴的螺旋或螺线，由 $\stackrel{\→}{s}\left(\mathrm{\θ}\right)=[D\cos \left(\mathrm{\θ}\right),D\sin \left(\mathrm{\θ}\right),\mathrm{h\θ}/2\mathrm{\π}]$ 定义，其中半径D是常数，并且h被称为螺旋的斜度。

参数θ在这种情况下指出在(x，y)平面中相对于x轴的源角。假设检测器平面一般包含z轴，并且围绕着源运动旋转，以致对每个θ来讲，它都垂直于源向量 在(x，y)平面上的投影。在这种情况下，在检测器平面中的原点O_θ落在z轴上。t₂轴与z轴重合，而t₁轴平行于(x，y)平面。(其它普遍的选择能够通过坐标旋转减少到这里描述的情况，例如，其中，

在位置θ平行于轨道的切线，而 正交于该轨道的切线。)

虽然在实践中，检测器平面可以定位不同，但是，利用沿着(t₁，t₂)坐标的统一间距，简单的线性坐标变换将把在实际检测器平面上获得的投影转换为在上面假设的虚拟检测器上获得的投影。

具有单一圆形轨道和平面检测器的锥形波束成像

图4中所示这种情况下的源轨道是螺旋轨道的一个特殊情况，具有零斜度，h＝0，即，其位于在x，y平面中在原点O处进入的半径为D的单圆形 $\stackrel{\→}{s}\left(\mathrm{\θ}\right)=[D\cos \left(\mathrm{\θ}\right),D\sin \left(\mathrm{\θ}\right),0]$ 上。假设检测器平面不失一般性地包含z轴，并且垂直于源到中心的线 SO。

考虑在目标中的点 $\stackrel{\→}{r}={[x,y,z]}^T.$ 对检测角θ，在检测器上的锥形波束投影的位置 由下面公式给出：

${\mathrm{\τ}}_1(\mathrm{\θ},\stackrel{\→}{r})=\frac{D(y\cos \mathrm{\θ}-x\sin \mathrm{\θ})}{D+x\cos \mathrm{\θ}+y\sin \mathrm{\θ}}---\left(2\right)$

${\mathrm{\τ}}_2(\mathrm{\θ},\stackrel{\→}{r})=\frac{D_Z}{D+x\cos \mathrm{\θ}+y\sin \mathrm{\θ}}---\left(3\right)$

这些清楚的表达式对在后面对快速算法的解释是很有用的。

具有圆形轨道和共线性等间距检测器的扇形断层摄影术

如图5所示，在这种情况下，所述目标的单一横截面被成像，并且再现2D图像。所述源轨道仍然是在(x，y)平面中以原点O为中心、并具有由源角θ指出的源位置的半径为D的圆形。然而，该检测器被限制在(x，y)平面中的单线(t轴)上，并且不失一般性地假设通过源的旋转O的中心，并且垂直于源到中心的线 SO。检测器元素在t轴上被均匀地分隔开。

用 $\stackrel{\→}{r}={[x,y]}^T$ 表示在2D图像中的位置，并用

表示与通过点 的源射线的t-轴相交的位置t(见图5)。目标f在源角θ和检测器位置t处的扇形波束投影(Pf)(θ，t)是沿着由参数(θ，t)表示的源射线的线积分，并由下式给出：

$\left(\mathrm{Pf}\right)(\mathrm{\θ},t)={\∫}_{\stackrel{\→}{r}\∈R^2}f\left(\stackrel{\→}{r}\right)\mathrm{\δ}[t-\mathrm{\τ}(\mathrm{\θ}-\stackrel{\→}{r})]d\stackrel{\→}{r}---\left(4\right)$

投影是用统一间距Δθ＝θ_max/P在P离散源角θ_p＝pΔθ，p＝0，...P-1处获得的。(Pf)(θ_p，.)项(对于所有的t值)是在源角θp处的投影，而用于不同p的投影集合则是窦腔造影。在全扫描的情况下，最大源角θ_max等于2π，但是在快速扫描或过扫描的其它情况下也能够采用其它值。部分数据组能够包括在任意角度子集处的投影。

在发散波束图形是一个曲面而不是平面的情况下，经常利用经修改的射线采样方案使用检测器表面。例如，对于螺旋或单一圆形轨道锥形波束成像来讲，检测器表面通常被选择位于以源为中心的圆柱形的表面，并且射线采样在方位角中是等角的，在高度z中是等间距的。图6示出了所述单一圆形源轨道，其中曲面检测器被选择包含一z轴。

在圆柱形和平面检测器形状之间的关系一般通过研究扇形波束几何形状来获得。在图7所示的具有等角射线采样的扇形波束几何形状中，检测器d在以源S为中心的弧线上，在角度上以间距Δt等间距。因此，在这种成像几何形状中，扇形角度t代替了在共线性等间距扇形波束几何形状中出现的检测器位置t。如图7所示，现在用

表示通过点

的射线的扇形角度位置。然后再次用公式(4)给出投影。

试图使本发明也适用于同样是发散波束几何形状的特殊情形的扇形波束几何形状的其它实例。

发散波束滤波反向投影再现和再投影

发散波束数据的自然(直接)倒置、无论是近似地还是精确的、都能够被作为加权的滤波的反向投影来阐明。首先，发散波束窦腔造影被预处理，例如，包括窦腔造影的投影被加权和滤波，以生成包括相应于源位置θp的投影

的改进的发散波束窦腔造影。在扇形波束或近似锥形波束再现的情况下，该滤波可以是简单的不变斜度移位的滤波，或在精确锥形波束再现的情况下，是更加涉及移位变化的滤波。在用于从部分扫描的长目标的螺旋锥形波束数据的再现的精确或准精确方法的情况下，加权和滤波甚至可以更加一般的生成两个或多个投影组，其贡献被相加。为了简短和一般性，将被反向投影并进行了用于从原始发散波束投影数据中获得它的预处理的数据都将被称作“发散波束投影数据”或“发散波束窦腔造影”。对于固定的p，它将被称作“发散波束投影”(在源位置θp)，并用函数g(p，.)表示。

然后，利用传统的发散波束离散-θ反向投影来再现所述3D图像：

$f\left(\stackrel{\→}{r}\right)=\underset{p=0}{\overset{P-1}{\mathrm{\Σ}}}W(\mathrm{p\Δ\θ},\stackrel{\→}{r})g[p,\stackrel{\→}{\mathrm{\τ}}(\mathrm{p\Δ\θ},\stackrel{\→}{r})]\mathrm{\Δ\θ},---\left(5\right)$

其中，

是一个适当的加权函数。这个离散反向投影公式与用于使用对所有θ测量的投影的反向投影的积分表示近似。实际上，因为离散检测器的使用，

也是在 中被采样的。由于

通常不对应可用的采样位置，所以，这就要求在变量

中内插g以实现反向投影。除了τ是标量而 是二维的以外，用于2D扇形波束情形的该反向投影公式具有同样的形式。

3D锥形波束反向投影的计算成本对N×N×N图像来讲是cN³P，其中，常数c依赖于诸如内插复杂度的实施细节。相反，当滤波是不变移位、且使用FFT执行卷积的时候，加权和斜坡滤波的计算成本只是O(PN²logN)。即使更精确形式的滤波也通常能够以类似的成本执行。从而，反向投影的成本支配传统锥形波束再现的成本，当通常P＝O(N)的情况下，其具有O(PN³)或O(N⁴)的成本。情形与2D扇形波束再现中相似，其中滤波和反向投影步骤的复杂度分别为O(N²logN)和O(N³)。

电存储图像f的自然(直接)发散波束再投影是公式(1)的实施。实际上，f只在离散形式下可用，并且只计算在

的

中的采样。例如，f可以被表示为一个系列：

$f\left(\stackrel{\→}{r}\right)=\underset{i,j,k=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\β}}_{\mathrm{ijk}}{\mathrm{\φ}}_{\mathrm{ijk}}\left(\stackrel{\→}{r}\right),---\left(6\right)$

其中，系数β_ijk是f的离散表示，而 是诸如样条函数的张量积的局部积函数或标准像素基函数，即，在第ijk个体元上的指标函数，对于2D情形，这描述在美国专利Nos.6263096和6287257中，而对于3D情形，这描述在美国专利No.6332035中。所述投影能够通过下面公式计算：

$\left(\mathrm{Pf}\right)({\mathrm{\θ}}_p,\stackrel{\→}{t})=\underset{i,j,k=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\β}}_{\mathrm{ijk}}{\widetilde{\mathrm{\φ}}}_{\mathrm{ijk}}({\mathrm{\θ}}_p,\stackrel{\→}{t}),---\left(7\right)$

其中，如公式(1)所定义的， ${\widetilde{\mathrm{\φ}}}_{\mathrm{ijk}}({\mathrm{\θ}}_p,\stackrel{\→}{t})=\left(P{\mathrm{\φ}}_{\mathrm{ijk}}\right)({\mathrm{\θ}}_p,\stackrel{\→}{t})$ 是基函数

的发散波束投影。基函数的这些投影通常能够使用一个解析表达式来计算，或存储在一查询表中。如反向投影情况一样，自然再投影的成本对于3D为O(N⁴)，或对于2D为O(N³)。

在曲线检测器的情况下，其中t₁，t₂是在检测器上的适当坐标，现在用

表示在通过点

的源射线的检测器平面上的(t₁，t₂)位置。仍然由公式(1)给出的该投影能够使用公式(7)计算，其中 ${\widetilde{\mathrm{\φ}}}_{\mathrm{ijk}}({\mathrm{\θ}}_p,\stackrel{\→}{t})=\left(P{\mathrm{\φ}}_{\mathrm{ijk}}\right)({\mathrm{\θ}}_p,\stackrel{\→}{t})$ 表示在曲线的检测器上基函数

的发散波束投影。同样的，再现能够再次表示为经适当修改的投影 的加权反向投影，其中t₁，t₂是在曲线检测器上的适当的坐标。和以前一样，所述反向投影由同样的公式(5)给出，但使用了不同的加权函数 例如，在2D等角扇形波束几何形状的情况下，扇形角t代替了在等间距共线性几何形状的线检测器上的位移t。基于

的离散性的考虑与平面检测器算法相同，就是计算复杂度的考虑。

与已经讨论过的一样，不论是在2D或是在3D的情况下，反向投影步骤(和再投影，如果用到的话)通常都具有最大的计算量，因而这就是本发明的目标。将描述一种适用于所有发散波束几何形状的快速反向投影的新颖的一般处理。为了简化说明，将首先被描述2D扇形波束几何形状的快速反向投影和再投影。

用于共线性等间距扇形波束几何形状的快速自然反向投影和再投影

如图8所示，考虑f的子图像 的反向投影操作。使

是由下面公式定义的图像截断算子：

其中， ${\left|\right|\stackrel{\→}{r}\left|\right|}_{\∞}=\max \{x,y\}.$ 因此， ${f^{\′}=K}_M\left[\stackrel{\→}{\mathrm{\δ}}\right]f$ 是f的子图像，该子图像以 为中心、且尺寸为M×M。使用在源角pΔθ，p＝0...Q-1处的Q投影在子图像f’(位于整体坐标系中其正确位置处)上的整体坐标系中的反向投影由下面公式给出：

$f^{\′}\left(\stackrel{\→}{r}\right)={\mathrm{\β}}_{M,Q}\left[\stackrel{\→}{\mathrm{\δ}}\right]g\left(\stackrel{\→}{r}\right)=K_M\left[\stackrel{\→}{\mathrm{\δ}}\right]\underset{p=0}{\overset{Q-1}{\mathrm{\Σ}}}W(\mathrm{p\Δ\θ},\stackrel{\→}{r})g[p,\mathrm{\τ}(\mathrm{p\Δ\θ},\stackrel{\→}{r})]\mathrm{\Δ\θ},\mathrm{\Δ\θ}={\mathrm{\θ}}_{\max }/Q---\left(9\right)$

项表示相关的反向投影算子。因此，特别是， $f={\mathrm{\β}}_{N,p}\left[\stackrel{\→}{0}\right]g.$

因为反向投影的局限性，只有部分投影g(p，.)将所述反向投影分配给f’。这部分被表示为

其中

是算子，其对于每个θ_p在变量t中将g(p，t)截断到子图像 $f^{\′}=K_M\left[\stackrel{\→}{\mathrm{\δ}}\right]f$ 的支集的投影的支集Ξ＝ AB(见图8)。对所有源角θ_p能够预先计算确定

的截断间隔，以及所关心的子图像尺寸M和位置

或者，能够使用更简单的保守的近似，在螺旋算法中有些增加的计算成本承认需要移位和分样的投影的所需的支持的过度估计。

其符合反向投影的定位，如果 $f={\mathrm{\β}}_{N,P}\left[\stackrel{\→}{0}\right]g,$ 则

$f^{\′}=K_M\left[\stackrel{\→}{\mathrm{\δ}}\right]f={\mathrm{\β}}_{M,P}\left[\stackrel{\→}{\mathrm{\δ}}\right]{\hat{K}}_M\left[\stackrel{\→}{\mathrm{\δ}}\right]g,---\left(10\right)$

即，能够通过近似截断的投影的尺寸(M，P)的反向投影获得f’上的反向投影。

现在考虑在不相重叠的子图像J＝(N/M)²中的一部分图像f，每个子图象的尺寸是M×M，

$f=\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{K}_M\left[{\stackrel{\→}{\mathrm{\δ}}}_j\right]f,---\left(11\right)$

其中，向量 j＝1，...J是子图像在整体坐标系的中心，如图9所示。应用公式(10)，我们获得下面在子图像上的反向投影中的反向投影的精确分解。

$f={\mathrm{\β}}_{N,P}\left[\stackrel{\→}{0}\right]g=\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{\mathrm{\β}}_{M,P}\left[{\stackrel{\→}{\mathrm{\δ}}}_j\right]{\hat{K}}_M\left[{\stackrel{\→}{\mathrm{\δ}}}_j\right]g---\left(12\right)$

这种分解在图10中对N/M＝2进行了说明。发散波束窦腔造影₃22，也在公式(12)中用g表示，并被划分为多个对应刚描述的子图像的子窦腔造影。所述子窦腔造影能够是任何期望的尺寸，和一个像素或一个体元一样小。由截断算子1010、1012、1014和1016生成的子窦腔造影分别单独地反向投影在整体坐标系中的1018、1020、1022和1024。反向投影的窦腔造影生成子图像f₁、f₂、f₃和f₄，它们在1026会聚以生成图1中的电子图像₁26。这过程生成扇形波束反向投影的精确分解。对N/M＝2，分解成四个子窦腔造影和子图像，如图10所示。

再参考图10，线22下的P表示窦腔造影或子窦腔造影中的投影编号。通过自身，与单步反向投影β_N，P相比，这种分解并不提供加速。实际上，虽然单个子图像反向投影β_M，P的计算成本cM²P是比全部图像的单步β_N，P小的(N/M)²＝J倍，但是对需要的J的总成本，这样的子图像反向投影是保持同样的(可能具有薄记的一些小的额外开销)。精确分解仍然能够被用来通过分割在J平行处理器之间的反向投影来提供加速，或允许重复的使用需要少J倍的存储器的单一的处理器。

快速反向投影算法使用具有下面附加性质的分解思想。对于固定图像分辨率(带宽)，再现图像所需要的投影的数量与图像的尺寸成比例。即，如果P投影需要再现N×N图像，则P’＝(M/N)P投影足够再现同样分辨率的M×M图像。

快速反向投影算法因而是基于从减少了数量的投影中再现子图像 $f^{\′}=K_M\left[\stackrel{\→}{\mathrm{\δ}}\right]f$ 的。这种较少过程将通过算子

执行的，其自身由几个现在将要定义的算子组成。

投影移位算子

和 沿着t轴移位每个投影

量被定义如下：

${\hat{M}}_{-}\left[\stackrel{\→}{\mathrm{\δ}}\right]g(\mathrm{\θ},t)=g[\mathrm{\θ},t-\mathrm{\τ}(\mathrm{\θ},\stackrel{\→}{\mathrm{\δ}})]---\left(13\right)$

${\hat{M}}_{+}\left[\stackrel{\→}{\mathrm{\δ}}\right]g(\mathrm{\θ},t)=g[\mathrm{\θ},t+\mathrm{\τ}(\mathrm{\θ},\stackrel{\→}{\mathrm{\δ}})]$

L折叠角分样算子D_↓L通过后面是角子采样的角滤波将一组P投影中图像的数量减少到P/L个投影。滤波可以简单为平均在相邻角处的投影，或者可以更精致，并且被表示为：

$\tilde{g}(p,t)=h_D(p,t)*g(p,t)---\left(14\right)$

其中*表示卷积而h_D是滤波器内核。滤波可以只发生在用于每个固定t的p(即，在离散θ)中，或可以是在θ和t中分离的或不分离的组合滤波。除了每L个滤波的投影

以外，所述角子采样只保留一个。

算子

即截断、移位和角分样操作的合成被定义如下：

$g^{\′}(q,t)=O[L,M,\stackrel{\→}{\mathrm{\δ}}]g(q,t)={\hat{M}}_{-}\left[\stackrel{\→}{\mathrm{\δ}}\right]W_2{D}_{\↓L}{W}_1{\hat{M}}_{+}\left[\stackrel{\→}{\mathrm{\δ}}\right]{\hat{K}}_M\left[\stackrel{\→}{\mathrm{\δ}}\right]g(q,t),q=0,...,P\text{'}-1---\left(15\right)$

其中，P’＝P/L。算子

包括在分样之前和之后利用加权函数W_k(θ，t)，k＝1，2的乘法加权。一个可能的选择是 $W_1(\mathrm{\θ},t)=W[\mathrm{\θ},\mathrm{\τ}(\mathrm{\θ},\stackrel{\→}{\mathrm{\δ}})],$ 和W₂(θ，t)＝1/W₁(θ，t)。一般地说，可以选择加权以最优化算法的精确度，或完全忽略。

图11更详细地表示能够处理子窦腔造影的方式。例如，首先在步骤1110能够进行截断 接着是在步骤1112的移位算子

在步骤1114的加权，接着是在步骤116的分样，步骤1118的另一个加权，和步骤1120的第二移位。这个处理描述了窦腔造影分解算子

利用这些定义，用于在子图像f’上反向投影的精确公式(10)被近似地代替为：

$f^{\′}=K_M\left[\stackrel{\→}{\mathrm{\δ}}\right]f={\mathrm{\β}}_{M,P/L}\left[\stackrel{\→}{\mathrm{\δ}}\right]g^{\′}---\left(16\right)$

$={\mathrm{\β}}_{M,P/L}\left[\stackrel{\→}{\mathrm{\δ}}\right]O[L,M,\stackrel{\→}{\mathrm{\δ}}]g---\left(17\right)$

其中，β_M，P/L是使用由公式(9)定义的P/L投影在M×M子图像上的在整体坐标系中的反向投影。这导致为近似为公式(12)分割的图像的反向投影操作的近似分解，

$f={\mathrm{\β}}_{N,P}\left[\stackrel{\→}{0}\right]g=\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{\mathrm{\β}}_{M,P/L}\left[{\stackrel{\→}{\mathrm{\δ}}}_j\right]O[L,M,\stackrel{\→}{\mathrm{\δ}}]g---\left(18\right)$

这个分解在图12中进行了说明，该图示出了用于L＝N/M＝2并具有减少再投影数量的与扇形波束反向投影相关的近似分解。为了执行近似分解，快速反向投影处理器24被编程以划分发散波束窦腔造影₃，如图12中所指出的，并且包括由用g₁、g₂、g₃和g₄表示的、每个都包括 投影的1210、1212、1214和1216生成的P投影四个子窦腔造影。子窦腔造影g₁、g₂、g₃和g₄是由1218、1220、1222和1224在整体坐标系中反向投影的，以便在整体坐标系中的正确位置上生成子图像f₁、f₂、f₃和f₄。子图像在1226处会聚以生成电子图像26，仍在图12中指示为f。

参考图12，对于在每个J＝(N/M)²子图像上反向投影β_M，P/L的操作数量是cM²P/L，生成用于全部反向投影的总共cN²P/L。这比单步反向投影β_N，P少L倍。基于前面关于利用图像尺寸的投影数量的标度的观测，L能够最大被选择为L＝(N/M)。当包括 的计算成本时，子图像的最有数量能够针对在美国专利No.6263096中所示的单级分解来确定。

在快速反向投影算法的一个优选实施例中，根据公式(11)，在进一步将子图像分割为更小的子图像的每一级处递归地应用在公式(12)和(18)中的分解。在图13中示出了一个图像的二级(N/M＝2)二重分割的例子。还示出了在不同级的一些子图像中心位置。请注意，它们中所有的都是指第零级图像的原点O，即，指整体坐标系。分别在图14和15中示出了二级近似分解的一个对应的例子，和有一个精确分解和一个近似分解的级联组成的二级分解。在图14中，发散波束子窦腔造影₃22包括P投影，在快速反向投影处理器24中通过将窦腔造影划分为每个都包括P/2投影的子窦腔造影1410、1412、1414和1416来进行处理。这个近似细分是通过图11的流程图描述的近似分解算子来执行的。子窦腔造影1410的进一步近似分解仍使用图11的流程图描述的算子在1418、1420、1422和1424生成，但是这次使用对应图13中的更小的子图像的不同的参数。这些包括P/4投影的子窦腔造影的每一个都在整体坐标系中被反向投影在1426、1428、1430和1432，以便在如图13所示的整体坐标系中的正确位置生成多个子图像f_1，1、f_1，2、f_1，3、f_1，4也是在整体坐标系中的正确位置，在1434会聚所述子图像以生成子图像f₁。

由1416表示的子窦腔造影也类似地在1436、1438、1440和1442进行细分。这些子窦腔造影分别在1444、1446、1448和1450处被反向投影，并在1452处会聚以生成子图像f₄。虚线方框1454和1456分别表示f₂和f₃的相应的第二级分解。在图14中，所有的分解都是对于L＝N/M＝2的扇形波束投影的近似分解。

图15说明对于L＝N/M＝2的扇形波束投影的二级混合(一个精确，一个近似)分解。在图15中，包括P投影的窦腔造影22分解为精确子窦腔造影1510、1512、1514和1516，其中的每个也都包括P投影。分别表示为g₁、g₂、g₃和g₄的精确子窦腔造影使用图11中说明的近似分解算子(利用适当的参数)被分解为近似子窦腔造影，其中的每个近似子窦腔造影都包括P/2投影。例如，g₁在1518、1520、1522和1524被分解为子窦腔造影g’_1，1、g’_1，2、g’_1，3和g’_1，4。那些子窦腔造影在整体坐标系中被反向投影在1526、1528、1530和1532，以便在整体坐标系的正确位置处生成子图像f_1，1、f_1，2、f_1，3、f_1，4。那些子图像在1534处会聚以生成子图像f₁。

在图15中表示为g₄的1516处生成的子窦腔造影还在1536、1538、1540和1542进一步被分解为四个子窦腔造影以便生成子窦腔造影g’_4，1、g’_4，2、g’_4，3和g’_4，4，它们在1552处会聚以生成子图像f₄。虚线方框1554和1556分别表示对于g₂和g₃的相应的第二级混合分解。那些过程生成子图像f₂和f₃，并且子图像f₁、f₂、f₃和f₄在1558被会聚以生成图像26，在图15中用f指示。

递归能够被继续直到子图像达到期望的最小尺寸M_min为止，然后使用公式(5)执行在整体坐标系中的反向投影β_Pmin，Mmin。最佳的M_min通常是依赖于实施的。它可以被选择为小到M_min＝1，以致最小的子图像是1个像素或体元的大小。假设N是在分解的所有级的二和N/M＝2的乘方，这里将会有log₂N这么多的乘方。再假设L＝(N/M)＝2和P＝kN，其中k是一个小的整数，最后一级将涉及N₂个单一像素反向投影β_P/N，1，总成本为cPN。假设固定长度的内插在

的实施中被用于分样和移位操作，它能够表示每一级的计算成本将是c₁PN，其中c₁是常数。从而，log₂N级的总成本为O(PNlogN)，并且因为P＝O(N)，递归分解算法的总成本变为O(N²logN)。

实际上，因为离散检测器的使用，投影数据以间隔Δt在t中采样。因此，投影移位算子 将以离散的形式实施，但将需要移位Δt的非整数倍。因为投影在反向投影以前已经被滤波，所以它们是有频带限制的，并且能够使用例如在美国专利No.6287257中描述的内插和再采样的组合来移位。特别是，移位Δt的非整数倍能够分为移位Δt的整数倍和移位Δt的分数倍的级联。整数倍的移位对应简单的改变符号，而分数倍的移位能够通过在t变量中的内插来执行。此外，这个内插能够在公式(14)中的一般的分样滤波步骤中，通过相应于具有滤波器h_D(p，t)的在径向维度中卷积的径向滤波来来执行。

重要的是，利用这里提出的特定投影移位方案，投影的采样间隔保持统一和恒定，以致内插和再采样能够作为离散对离散映射来执行，并向简单的数字滤波器一样有效地实施。在这样的分数t-移位中涉及的误差能够随着内插的长度而指数地衰减，以致即使是简单的线性内插的短内插也通常足够。t-内插精度通过对t中的投影的过采样还能够增加，如果需要，其能够通过投影的数字内插(上采样后的滤波)来有效地执行。

使用近似分解获得的反向投影的精度依赖于处理中的各种参数，包括P、N、M、L，图像中的分辨率和投影的带宽，与目标相交的源射线最大扇形角度，检测器响应，算法中使用的内插器。这些参数能够在发散波束算法中适当的最优化。影响精确度的一个重要的参数是在反向投影中使用的投影的数量和子图像的尺寸之间的比率(P/L)/M。通过增加这个比率，也称为角过采样，近似分解的精确度增加了。

这里有几种方式来增加角过采样，包括：(i)投影g的数字角内插，以便在分解开始之前增加P；(ii)在近似分解的各级处选择L＜(N/M)；以及(iii)使用精确分解步骤，例如，如图15中的例子所示。为了理解方法(iii)，在公式(12)中精确分解的重新调用减少了子图像的尺寸，而不是子窦腔造影中的投影数量。又此，这将导致由因数(N/M)引起的角过采样的增加，并且改善了精确度。

事实上，除可以给出方法(ii)或(iii)中之一的比其它方法有很小计算优势的实施细节以外，由于无论何时都能够保持

位于中心在 的子图像的子图像的中心的恒等式 ${\hat{K}}_{M/2}\left[{\stackrel{\→}{\mathrm{\δ}}}_{\mathrm{ij}}\right]{\hat{K}}_M\left[{\stackrel{\→}{\mathrm{\δ}}}_i\right]={\hat{K}}_{M/2}\left[{\stackrel{\→}{\mathrm{\δ}}}_{\mathrm{ij}}\right],$ 所以，角过采样的方法(iii)等同于方法(ii)。从而，它遵守下述原则，即，具有二重分割的一个精确和一个近似的二级分解等同于分割为具有N/4尺寸的16个子图像的单一近似分解，其中分样L＝2，而不是最大可能的L＝4。因此，两种方法提供2角过采样的同样因数。另一方面，通过选择L＝3，单级近似分解到16子图像能够提供4/3因数的角过采样，因此提供在操作点的选择上的灵活性。一般的，增加角过采样的方法导致算法的计算成本增加常数因数倍。因此，它们提供在计算量和精确度之间的平衡。

算法的优选实施例包括许多精确细分和许多近似细分，或具有小于它们的细分因数的分样因数的近似细分。精确细分的数量和/或分样因数控制过采样因数，并且被选择以获得在计算量和精确度之间的平衡。

最简单的实施例是在每级处使用N/M＝2的常数因数的细分，其中在近似细分级处L＝2，但是其它选择也同样可能。

如所示，该算法要求P被在所有级中的分样因数L的乘积除尽。当这可以够通过分样因数的选择而方便地或有效地实现时，通过角内插和投影的再采样来改变P。

为了进一步说明算法的优选实施例的递归结构，图16示出了M/N＝L＝2的被写成递归函数的算法的伪码的例子。由于伪码是自身解释的，所以将被简要地描述。

在对全窦腔造影G的第一最高级调用呼叫中，函数FAN_FHBP是递归函数操作。在对其自身的后续连续呼叫中，这个函数将该图像分解为子图像并执行窦腔造影分解直到实现最小图像尺寸N_min为止。这些分解中的Q个是精确的，剩余的是近似的。当达到图像最小图像尺寸N_min时，对相应于子图像的子窦腔造影执行在整体坐标系中的精确反向投影。对于大于N_min的子图像，将所述投影截断成对应于所述子图像。这构成了精确窦腔造影细分。当完成精确分解时(当Q＜1)，执行近似分解的附加步骤。在图16中，投影被移位、加权分样，然后再移位。子图像在 $\mathrm{FAN}\_\mathrm{FHBP}=f={\mathrm{\Σ}}_{i=1}^4{f}_i\mathrm{FAN}$ 处会聚。

子图像中心也被进行递归的计算，中心在 处的较大子图像的子图像i的中心

表示为 ${\stackrel{\→}{\mathrm{\δ}}}_i=\stackrel{\→}{\mathrm{\δ}}+{\stackrel{\→}{\mathrm{\ξ}}}_i,$ 其中，

是子图像f_i相对于

的适当的位置矢量。根据定义，全尺寸图像的中心是整体坐标系的原点。因此所有后续子图像的中心被表示在整体坐标系中。

在快速反向投影算法的优选实施例中，通过下面的观测来简化递归分解。投影的操作在t中的移位、投影截断和投影加权换算(利用适当的坐标变换)。从而，一个级的操作 能够与下面级的 组合为一个移位和一个加权操作。同样的，在最后子图像反向投影β_Mmin，Pmin之前出现的移位和加权操作

能够与在反向投影自身中的加权和内插操作。因此，在递归实现中， 只需要一个投影移位步骤，以及每级(最多)一个加权。

本领域的技术人员应当明白，在单一处理算法设计和实施中，定义本发明的过程的各种操作能够以不同于由这里表示的图描述的那些顺序的各种顺序被再配置，或平行执行，而不改变过程的功能或原理。一些这样的再配置可能从实施、操作或成本观点上更优越。例如，用于分样一定数量与一子图像相关的投影的操作可以开始得和对该子图像的两个或多个投影已经被截断并移位的的速度一样快，而不用等待所有的P投影被截断并移位。这样的再配置能够被用来管线再现算法，并用数据获取对其进行交迭。

此外，本发明的方法能够被应用以随着时间获取的投影连续地更新断层摄影的影像。这能够通过使用提出的方法反向投影由在最近获取的投影和以前的一组投影之间的区别组成的窦腔造影来实现。然后反向投影的结果被加到但前图像以生成更新的图像。这样的更新能够以连续的方式来继续。

快速再投影

反向投影和再投影是双重操作，或彼此互相为伴随的操作。从而，快速自然扇形波束算法的实施是紧密相关的。从而，由于可以理解到快速反向投影算法的所有讨论都在再投影算法中具有一个类似物，所以，快速自然扇形波束再投影算法将不那么详细地进行描述。

下面考虑图8所示的f的M×M子图像 $f^{\′}=K_M\left[\stackrel{\→}{\mathrm{\δ}}\right]f.$ 在源角θ_q＝qΔθ，q＝0...Q-1处的扇形波束投影为：

$g(q,t)=P_{M,Q}\left[\stackrel{\→}{\mathrm{\δ}}\right]f^{\′}(q,t)=\left({\mathrm{PK}}_M\right[\stackrel{\→}{\mathrm{\δ}}\left]f\right)(q,t)---\left(19\right)$

我们用 表示中心在位置

的M×M子图像在整体坐标系中的相关的投影算子。特别是， $\left(\mathrm{Pf}\right)=p_{N,P}\left[\stackrel{\→}{0}\right]f$ 是整个图像的扇形波束投影。P_M，Q的计算成本对某些常数c是cM²Q。

考虑在图9中说明的公式(11)的下一个图像部分。因为投影算子的线性，对于源角的同样数量P，分解到较小的M×M图像的J再投影的总和的再投影的下面精确分解能够如下来获得：

$\left(\mathrm{Pf}\right)=\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{P}_{M,P}\left[{\stackrel{\→}{\mathrm{\δ}}}_j\right]K_M\left[{\stackrel{\→}{\mathrm{\δ}}}_j\right]f---\left(20\right)$

这种再投影的精确分解在图17中针对N/M＝2进行说明。在图17中，指示为f的电子图像₂30被截断算子1710、1712、1714和1716分别细分为子图像f₁、f₂、f₃和f₄。子图像被1718、1720、1722和1724再投影在整体坐标系中以生成子窦腔造影g₁、g₂、g₃和g₄。子窦腔造影在1726会聚以生成在图17中指示为g的窦腔造影34。

以反向投影的公式(12)的精确分解的情况下相同的原因，与单步反向投影 $\left(\mathrm{Pf}\right)=P_{N,P}\left[\stackrel{\→}{0}\right]f$ 比较，这种再投影的分解不提供加速，但是有其它可能的应用，特别是与近似分解组合。

快速再分解算法使用具有下面在早先提到的附加性质的分解思想：对于固定图像分辨率(带宽)，描述图像的扇形波束投影的数量；或等同的，在θ变量中的投影的带宽，与图像的尺寸成比例。因此，如果P投影描述一个N×N图像f，则P’＝P(N/M)投影描述子图像f’，并且剩余的投影能够从这个较小的组中获得。

从而快速再投影算法基于计算一组M×M子图像’的P’投影，并使用“投影增加处理”获得更大组的P投影。这种增加过程将通过算子

执行，其自身有若干算子组成。投影移位算子

和

是以前在公式(13)中定义的。L折叠角内插算子I_↑L通过角滤波之后的L折叠角上采样，增加一组Q投影的影像的数量为QL投影。角上采样在每两个投影之间引入了L-1个零值的投影。用g(p，t)指示的上采样的投影，滤波如公式(14)进行描述并被相关联地讨论。结果的角内插算子I_↑L是在快速反向投影算法中出现的角分样算子D_↓L的伴随矩阵。所述移位和角内插操作的组成部分的投影增加算子定义为：

$g^{\′}=(p,t)=v[L,M,\stackrel{\→}{\mathrm{\δ}}]g(q,t)={\hat{M}}_{-}\left[\stackrel{\→}{\mathrm{\δ}}\right]W_2{I}_{\↑L}{W}_1{\hat{M}}_{+}\left[\stackrel{\→}{\mathrm{\δ}}\right]g(q,t),p=0,...,P\text{'}-1---\left(21\right)$

其中P’＝QL。利用与

类似的函数和选择标准，算子

包括用加权函数W_k(θ，t)，k＝1，2在内插之前和之后进行乘法加权。例如，如图18所示，第一移位函数1810后面接着是第一加权函数1812和内插函数1814，第二加权函数1816和第二移位函数1818。如果需要，也能够使用投影增加算子的其它组合。

由于这些定义，子图像f’的近似再投影为：

$g^{\′}=P_{M,P}[\stackrel{\→}{\mathrm{\δ}}]f^{\′}=v[L,M,\stackrel{\→}{\mathrm{\δ}}\left]P_{M,P/L}\right[\stackrel{\→}{\mathrm{\δ}}]f^{\′}---\left(22\right)$

其中，P_M，P/L是在只计算P/L投影的整体坐标系中的再投影。这导致类似于公式(20)的分割的图像的再投影操作的近似分解

$g=P_{N,P}\left[\stackrel{\→}{0}\right]f=\underset{j=1}{\overset{J}{\mathrm{\Σ}}}v[L,M,\stackrel{\→}{\mathrm{\δ}}]P_{M,P/L}\left[{\stackrel{\→}{\mathrm{\δ}}}_j\right]K_M\left[{\stackrel{\→}{\mathrm{\δ}}}_j\right]f---\left(23\right)$

这种近似的分解在图19中对L＝N/M＝2进行说明。在图19中，来自图1的快速再投影处理器32被编程以接收图19中指示为f电子图像₂30，并分别在1910、1912、1914和1916中划分为子图像f₁、f₂、f₃和f₄。整体坐标系中的再投影在1918、1920、1922和1924执行以分别生成子窦腔造影g₁、g₂、g₃和g₄，每个具有P/2投影。在1926、1928、1930和1932增加子窦腔造影，如前面所述，结果的每个都具有P投影的子窦腔造影g₁、g₂、g₃和g₄在1934会聚以生成图19中指示为g窦腔造影34。

如反向投影的近似分解的情况一样，所述计算成本是小于单步再投影P_N，P的L折叠，并且应用相同的讨论。

在快速反向投影算法的优选实施例是根据公式(11)，在进一步将子图像分割为更小子图像每一级处，递归地应用公式(20)和(23)中的分解，并如图13所示。分别在图20和图21中示出了二级近似分解的相应例子，和由一个精确分解和一个近似分解的级联组成的二级分解的例子。在图20中，L＝N/M＝2的扇形波束再投影的二级近似分解开始于电子图像₂30，在图20中指示为f。图像f分别在2010、2012、2014和2016被划分为子图像f₁、f₂、f₃和f₄，而子图像f₁在2018、2020、2022和2024进一步细分为子图像f_1，1、f_1，2、f_1，3、f_1，4。这些子图像在整体坐标系中分别由2026、2028、2030和2032再投影，以生成子图像g”_1，1、g”_1，2、g”_1，3和g”_1，4，其中的每一个都包括P/2投影，它们在2042处会聚以生成包括P/2投影的子窦腔造影g₁’。在2044处的附加增加之后，包括P投影的子窦腔造影g₁在2046处与子窦腔造影g₂、g₃和g₄的会聚。

子图像f₄进一步在2048、2050、2052和2054处被分解，并且子图像类似地分别在2056、2058、2060和2062处进行再投影以生成子窦腔造影g”_4，1、g”_4，2、g”_4，3和g”_4，4。在2064、2066、2068和2070处的增加之后，生成的子窦腔造影g’_4，1、g’_4，2、g’_4，3和g’_4，4在2072处会聚，而生成的子窦腔造影进一步在2074处增加。生成的子窦腔造影g₄在2046处与其它子窦腔造影会聚以生成图1中的子窦腔造影34，也在图20中指示为g。虚线方框2076和2078分别表示相应的f₂和f₃的第二级分解。

图21示出了L＝N/M＝2的扇形波束再投影的二级混合(一个精确，一个近似)分解。在图21中，来自图1的电子图₂30、在图21中被指定为f，分别由2110、2112、2114和2116分解为子图像f₁、f₂、f₃和f₄。子图像f₁被2118、2120、2122和2124细分为子图像f_1，1、f_1，2、f_1，3和f_1，4。子图像在整体坐标系中分别被2126、2128、2130和2132再投影，以生成子图像g’_1，1、g’_1，2、g’_1，3和g’_1，4，其中的每一个都包括P/2投影。包括P投影的子图像g_1，1、g_1，2、g_1，3和g_1，4每个都是在被2134、2136、2138和2140增加后生成的。在2142处的会聚生成包括P投影的子窦腔造影g₁，其在2144处与子窦腔造影g₂、g₃和g₄会聚以生成包括P投影的窦腔造影g。类似地，子图像f₄在2146、2148、2150和2152处被进一步细分以生成f_4，1、f_4，2、f_4，3和f_4，4。由2154、2156、2158和2160的再投影生成子窦腔造影g’₁、g’₂、g’₃和g’₄，其被2162、2164、2166和2168增加以生成子窦腔造影g_4，1、g_4，2、g_4，3和g_4，4。这些子窦腔造影在2170会聚以生成子窦腔造影g₄。虚线方框2172和2174分被表示g₂和g₃的相应的第二级分解。

继续递归直到子图像成为期望的最小尺寸M_min，然后使用公式(7)计算再投影P_Pmin，Mmin。

下面考虑快速再投影算法，类似于对反向投影情形的讨论：

1.在快速再投影算法中的各种参数的选择，包括P、N、M、L和M_min。

2.处理t中的离散。

3.使用角过采样和精确与近似分解步骤的组合来控制精确度和速度之间的平衡。

4.通过组合在连续的级之间的移位和加权的算子 和 导致每一级只有单一移位和(最多)单一加权操作。

5.操作顺序的可能再配置。

6.生成的分级再投影算法的O(N²logN)计算复杂度。

等角扇形波束几何形状的快速自然算法

对于等角检测器的情形，能够应用类似的方案以分级的分解图像为较小的子图像。图22说明在这种情况下子图像f’的投影几何形状。

使用如公式(8)中的K_M同样的定义，子图像反向投影算子 由公式(9)给出，其中

如早前定义的等角几何形状一样(见图6)。g(p，.)的对在f’上的方向投影有贡献的部分仍然是 其中现在

截断扇形角中的g(p，.)为由子图像的 $f^{\′}=K_M\left[\stackrel{\→}{\mathrm{\δ}}\right]f$ 的支持的投影的角t_A和t_B确定的角。应用如在共线性等间距扇形波束几何形状相同的讨论和论点，从而导致公式(10)-(12)以及图9所示的精确分解。

使用下面处理，近似分解再次减少了用于再现子图像f’的投影的数量。令t_OO’是在分别经过f和子图像f’的中心O和O’的两个射线之间的角。在直接模拟在共线性等间距检测器情况下以对应于f’的中心O’的投影的间隔OO″对被截断的投影进行t移位的过程中，以扇形角将用于该子图像的被截断投影t移位t_OO’。以相同角间隔Δt均匀间隔的采样点是以t＝0为中心获得的。因此，再等角情形中角变量t担当出像在等间距几何形状中的t的作用，并且如图7中所定义的。否则，在这里应用前述的相同递归分解过程。特别是，相应于子图像的截断的扇形波束投影已经在t中截断时，他们对源角θ以L分样，并且接着剩余的程序。除了已经列出的改变，各种算子的定义几乎没有变化。近似分解仍用图11描述。

类似的论点和修改应用到快速自然再投影算法。

具有平面等间距检测器的锥形波束几何形状的快速自然反向投影和再投影算法

一般3D发散波束几何形状的处理具有为2D情况描述的相同成分，并具有下面主要不同：

1.图像在3D中沿着所有三个坐标分解。

2.因为投影是二维的，所以单独的投影的截断、移位和内插是二维的。

3.在源位置θ中的投影被沿着轨道分样(或内插)，而不是在圆周上的源角中。

4.因为源轨道通常在3D中不是各向同性的，所以优选的分解可以不是各向同性的，并且可以是空间变化的。

考虑图23所示的f的子图像

的反向投影操作。用

表示由公式(8)定义的图像截断算子，所以 $f^{\′}=K_M\left[\stackrel{\→}{\mathrm{\δ}}\right]f$ 中心的

的尺寸为M×M×M的f子图像。仍应用公式(9)和(10)，并且τ用2D向量

代替，并且截断整个目标的投影的截断算子

为子图像 $f^{\′}=K_M\left[\stackrel{\→}{\mathrm{\δ}}\right]f$ 的M×M×M支持的投影的2D支持Ω(见图23)。

现在考虑J＝(N/M)³非重叠子图像中的图像f的一部分，每个尺寸M×M×M，

$f=\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{K}_M\left[{\stackrel{\→}{\mathrm{\δ}}}_j\right]f---\left(24\right)$

其中，向量 ${\stackrel{\→}{\mathrm{\δ}}}_j,j=1,...J$ 是子图像的中心，如图24所示M/N＝2的情形。

应用公式(12)，提供反向投影的精确分解给到子图像上的反向投影。除了它具有八个而不是四个“信道”以外，其描述类似于图9。除了如2D情形中的比率J、β_N，P和β_M，p的计算成本分别为cN³P和cM³P以外，计算成本的讨论类似于2D的情形。

如2D扇形波束情形，快速3D锥形波束反向投影算法使用具有下面附加性质的分解思想。对于固定图像分辨率(带宽)，再现图像所需的投影的数量与图像的尺寸成比例。即，如果P投影需要再现一个N×N×N图像，则P’＝(M/N)P投影足够以相同分辨率再现一个M×M×M图像。

子图像 $f^{\′}=K_M\left[\stackrel{\→}{\mathrm{\δ}}\right]f$ 仍是通过使用算子

执行减少而从减少数量的投影中再现的。它由如前面的类似算子合成所定义，并具有以下改变：

1.如前面所解释的，现在投影截断算子 截断到2D区域Ω。

2.投影移位算子 通过向量

在t₁、t₂坐标中执行2D移位，而不是在扇形波束几何形状情况下的1D移位。这种移位可以用分离的方式执行，例如，通过首先在t₁中移位然后在t₂中移位。

3.算子D_↓L在源位置中而不是在角中进行分样。当轨道由两个或多个不同的部分的联合组成时，诸如两个垂直的圆周，对于每部分分离地执行源位置中的分样；

4.涉及由公式(4)的分样的 变量中的滤波是二维而不是一维。它也可以用分离的形式在t₁和t₂中执行，虽然

滤波可能能够也可能不能与θ中的滤波分离。

近似分解仍由公式(18)给出，其中现在β_M，P/L是使用如公式(9)定义的P/L投影的在M×M×M子图像上的整体坐标系中的发散波束反向投影。对于L＝N/M＝2，除了具有八个而不是四个“信道”以外，这个分解的描述类似于图10的描述，。计算成本的讨论类似于2D情形，除了β_N，P和β_M，P/L的计算成本分别为cN³P和cM³P/L，并且J＝(N/M)³。J子图像反向投影的总成本仍是小于传统的整幅图像反向投影β_N，P的成本的L倍。

在快速反向投影算法的优选实施例中，根据公式(24)，在进一步将该子图像分解为更小子图像的每一级处，递归地应用公式(12)和(18)的分解。。扇形波束几何形状的过程和各种特性的描述，以及前面所示的改进也应用到一般锥形波束几何形状，其中的改变与精确和近似分解本身列出的类似。特别是，与传统锥形波束反向投影的O(N⁴)相比，使用递归分解算法的锥形波束反向投影的总成本变为

对于这一点，我们已经描述了3D图像分解的最简单的实施例，其中图像在所有三个轴中均等地分解。然而，在锥形波束几何形状的一些实例中，因为源轨道在3D中不是各向同性的，优选的分解也可能不是各向同性的，并且可能是空间变化的。重要的例子是单圆形和螺旋形锥形波束轨道(图4、6和3)。这些情况展示了围绕z-轴的源轨道的圆柱形对称(或对称)，并且在z轴中选择的图像部分可以不用于在x和y轴中的图像。即使当整个轨道不具有对称轴，它有时也可以分解成具有对称部分的单元，如源轨道由两个垂直圆周组成的情形。所述反向投影操作能够被作为两个反向投影的总和来执行，其中的每一个都具有其自己的分离的轴对称轨道。

在这个例子中，不太频繁地分割一个图像的多个部分或将其分割成Z轴上具有大尺寸的子图像就已经足够了。在一个极端的例子中，可以将所述分割限制在(x，y)方向上，如图25所示。图26示出了更加一般的各向异性和空间的非同一分割。这可能通过减少所需对沿着t₂轴的投影的移位、内插或滤波操作导致节省计算。

这样的减少的分割起作用的原因在单一圆周源轨道的情况下通过检查公式(2)和(3)可以看出。公式描述了一个点在检测器上投影位置的θ相关性。考虑一个具有很小的z和/或很小的

的点 即，一个靠近源轨道平面和/或接近z轴的点。注意对于这样的点，

随着θ的变化将比

慢。这能够用来表示其中心这样定位的子图像能够从减少的数量的影像中精确地再现，即使其z-维度没有与其(x，y)维度同样的比例减少。结果，在(x，y)平面中的精细分割已足够保证从减少的数量的投影中的精确地再现。相应的论点也应用到螺旋轨道的情况，如果首先对投影应用在t₂中整体移位hθ/(2π)，以补偿在检测器上沿着z轴的坐标中心随着θ改变的移位(见图3)。

由于这种各向异性和可能不统一的分割，递归分解因此而被修改。如果图25的分割是被递归地使用，相应的分解在只涉及在t₁变量中的移位的每级递归处都将具有四个信道，并且将由图10-15描述。另一方面，如果图26中的分割被用在二级近似分解中，第一级分解将具有八个信道(相应于立方体的粗略分割的八个卦限)，但是然后在下一级中每个分解将只具有对应于在每个分割的粗糙图像卦限中的七个图像的七个信道。

包括先前讨论的发散波束几何形状的特定实例的对一般发散波束几何形状的快速再投影的设计是相应快速反向投影算法的双重算法。如同快速扇形波束再投影具有快速扇形波束反向投影一样，所述快速再投影具有与快速发散波束几何形状反向投影相同的关系。因为所述关系已经详细地描述了，所以在这里不再说明。

实施和实验

新算法已经利用MATLAB^TM和C编程语言实现，并成功地在几种涉及在2D情形中的N×N＝512×512尺寸的谢普-罗庚(Shepp-logan)头部模型(用在断层摄影算法的数字评估中的标准测试)的再现的数字进行了实验。我们使用所谓的未改变的模型，其具有在头骨和脑之间的现实的(和高度的)对比。这表现了对算法的富有挑战的测试，因为即使是由头骨的再现所创建的轻微的假象也会在脑背景上突显出来。

我们首先描述2D扇形波束再现的结果。对于全部等间距的和等角的扇形波束成像几何形状，我们设置源到图像中心的距离为图像尺寸的1.25倍，即，D＝1.25N，并设置扇形角为2α₀ ＝1.17弧度，在扇形上具有1025个检测器。我们使用对通过椭圆的线积分的解析表达式生成具有P＝1536个平均间隔的源角的整个扇形的扇形波束投影数据，对10-椭圆模型θ_p∈[0，2π]。

在我们的仿真中，在每级都分解图像为一半尺寸的4个子图像。对于近似分解，我们使用分样因数L＝2。图像被迭代地分解成单一像素的子图像，即，M_min＝1。为了增加角过采样，首先的Q级分解是没有角分样的精确分解。我们使用Q＝2和Q＝3，其中Q是指“切换参数”。投影阶段算子

的参数(“注脚”)是以增加一些计算量的代价使用简单保守近似随时计算的。

图27和28针对两个扇形波束成像几何形状比较了精确反向投影算法和我们的快速扇形波束反向投影算法。作为参考和基准，图27(a)和28(a)表示原始的模型，图27(b)和28(b)是对等间距和等角检测器的使用精确扇形波束反向投影的再现的图像。图27(c)、28(c)、27(d)和28(d)是使用分别具有切换参数Q＝3和Q＝2的快速算法的再现。快速算法分别对切换Q＝3和Q＝2的加速是30×和60×。使用精确和快速算法的再现表示简直没有可只觉得质量区别，即使是Q＝2。通过图28中的再现的图像的行和列的曲线被表示在图29(a)和29(b)中。图29(a)比较通过图28的图像的第410行的片层，而图29(b)比较通过图28的图像的第356列的片层。在这些曲线中，实线表示原始模型，虚线表示精确扇形波束反向投影，点划线表示具有切换3的快速算法，而点线表示具有切换2的快速算法。

接下来，我们比较对于等间距的检测器几何形状的精确和快速扇形波束反向投影算法的点分布函数(PSF’s)。我们在尺寸256×256的图像上运行精确和快速算法，四个点目标分别位于位置(10，10)、(128，128)、(49，37)和(201，93)。我们使用的参数类似于之前所使用的：源到中心距离D＝1.25N，源角数量P＝2N，具有737等间距检测器的1.20弧度的总扇形角。在所有四个位置的点的响应都很类似，都暗示了主要得移位不变式响应。利用精确快速算法(以Q＝3)获得的在(128，128)点目标的点分布函数(PSF’s)在图30(a)和30(b)中进行比较。图31(a)表示精确算法(实线)和具有Q＝3的快速算法(虚线)的通过第128行的片层。图31(b)表示精确算法(实线)和具有Q＝3的快速算法(虚线)的通过第128列的片层。两种算法的PSF’s几乎相同，确认了快速算法的精确度。

接下来，我们讨论具有在(x，y)平面中的单一圆形源轨道的锥形波束再现的实验的结果。从源到体积中心的距离是侧面长度的1.25倍。锥形波束投影是以512源角取得的，其统一地以[0，2π)为间隔。扇形角(方位角)是1.20弧度，而锥形角是1.5弧度。我们考虑分别在图4和图6中示出的全部平面等间距的和圆柱形的检测器几何形状。在每种情况下，检测器表面由沿着垂直(t₂)轴等间距的、或在t₁坐标中等间距或等角的间隔的375×375检测器组成。

具有平面检测器的单一圆周锥形波束轨道的传统再现算法是著名的Feldkamp、Davis和Kress(FDK)算法，其由对投影进行加权、沿着t₁坐标滤波以及然后执行加权的锥形波束反向投影组成。加权和滤波与具有平面检测器的扇形波束几何形状相同，不同是使用锥形波束而不是扇形波束反向投影。这是一个近似算法，并且它在大锥形角时经受一些固有的假象。但是，它是这种锥形波束几何形状的标准基准。对于圆柱型检测器几何形状，我们使用的近似算法是通过将等角检测器几何形状的扇形波束算法扩展到锥形波束几何形状，通过将加权的锥形波束反向投影合并到圆柱形检测器上。

为了示范的目的，在图25中示出了具有简单分割方案的使用具有前面描述的递归分解的反向投影的相应的快速算法。在我们的仿真中，我们在子图像尺寸减小到4×4×128的时候终止分解过程，这是因为在算法的特别实施中涉及的花销抵消了这样小的体积的反向投影算子的加速。快速算法使用与传统FDK算法同样的滤波和加权。

传统FDK再现和使用本发明处理的快速再现的结果在图32和图33中分别对平面和圆柱形检测器几何形状进行了比较。图32比较单一圆形轨道的3D锥形波束再现和平面等间距的检测器几何形状。顶行是在z＝-0.25的xy片层，中间行是在x＝-0.13的yz片层，而底行是在y＝-0.04的xz片层。左边的列是使用传统FDK算法的再现结果，并且中间和右边的列分别使用具有切换2和1的本发明的快速算法。图33是单一源轨道的3D锥形波束再现和圆柱形检测器几何形状的比较。顶行是在z＝-0.25的xy片层，中间行是在x＝-0.13的yz片层，而底行是在y＝-0.04的xz片层。左边的列是使用传统FDK算法的再现结果，中间和右边的列分别使用具有切换使用具有切换2和1的本发明的快速算法。图34表示在图32所示相应行中贯穿所述图像的片层。顶部的曲线说明了行102，中间的曲线说明行80，而底部的曲线说明列81。实线来自传统FDK再现。利用切换2的快速算法用点划线示出，利用切换1的快速算法用点线示出。图35说明在图33中的相应行中通过图像的片层。顶部曲线说明行102，中间曲线说明行80，而底部曲线说明列81。实线来自传统FDK再现，点划线表示使用具有切换2的快速算法的结果，而点线表示使用具有切换1的快速算法的结果。附加的细节由图34和35提供，其表示沿着通过图32和33的图像的特定线的再现密度的曲线。给出了对全部切换＝1和切换＝2的结果。

使用本发明的具有切换Q＝2的FDK算法的快速版本大约比传统Feldkamp算法快三倍，并且几乎没有任何图像质量上的衰减。即使切换因数Q＝1，由于以7-折叠加速，其图像质量对某些应用也是可以接受的，。我们希望对512×512×512图像量的加速将增加几十倍，以及如图24或26中的完全合格的分割和相应的分解。(事实上，具有接近圆柱形对称的锥形波束几何形状，诸如圆形或螺旋轨道，将对任何z-尺寸H的512×512×H的图像获得几十倍的加速)。

试验结论

在我们2D中的试验中，新的扇形波束反向投影算法为512×512图像提供30×-60×的实际加速，并没有精确度的可知觉损失。在3D中，以单一圆形轨道和平面检测器锥形波束几何形状，并使用最简单的分解形式，对128×128×128的加速是3×-7×。几十倍的加速因数应当能够保证本发明的提出的分解的完全合格的实施，并且对于512×512×512图像。

试验对不同的检测器几何形状示范了本发明的过程的有效性。这些对简单单一圆形轨道几何形状的结果将扩展到更一般的发散波束几何形状，诸如螺旋锥形波束。因为在再投影算法和反向投影算法之间的原理和结果中的类似性，类似的结果应当对本发明的快速发散波束再投影算法也成立。再投影和反向投影算法单独都是有用的；组合起来，它们能够同于各种断层摄影再现的迭代方案，或用在断层摄影术中的假象校正，如例如在美国专利Nos.6263096和6351548中解释的。反向投影和再投影的组合也能够用来构造在螺旋锥形波束几何形状中的长目标问题的快速再现算法。

虽然本发明的原理已经结合特定装置和应用在上面进行了描述，应当理解这种描述只是为了示例，而不是对发明的范围进行限制。

Claims (42)

1.一种用于根据发散波束窦腔造影(22)生成电子图像的方法，所述发散波束窦腔造影适合于反向投影，该方法包括步骤：

将窦腔造影(22)细分(1010、1012、1014、1016)为多个子窦腔造影(g₁、g₂、g₃、g₄)；

在整体坐标系(1018、1020、1022、1024)中执行所述子窦腔造影的加权的反向投影，以便在整体坐标系中的适当位置处生成多个相应的子图像(f₁、f₂、f₃、f₄)；以及

会聚(1026)所述子图像(f₁、f₂、f₃、f₄₎以创建所述电子图像(f，26)。

2.根据权利要求1所述的方法，其中，所述窦腔造影细分是近似的(1210、1212、1214、1216)。

3.根据权利要求1所述的方法，其中，所述窦腔造影以递归方式细分为多个子窦腔造影，直到每个子窦腔造影都表示期望尺寸的子图像为止，这里所述细分步骤包括期望数量的精确细分(1510、1512、1514、1516)和期望数量的近似细分(1518、1520、1522、1524)。

4.根据权利要求3所述的方法，其中，所述精确细分步骤包括没有移位的截断。

5.根据权利要求3所述的方法，其中，所述子窦腔造影对应于在尺寸上与一个像素或一个体元一样小的子图像。

6.根据权利要求1所述的方法，其中，所述会聚步骤是以递归方式执行的。

7.根据权利要求1所述的方法，其中，所述电子图像是断层摄影图像。

8.根据权利要求1所述的方法还包括预处理(16)，其中，数据坐标中的过采样被用来改善电子图像的精确度。

9.根据权利要求1所述的方法，其中，窦腔造影包括预处理的发散波束投影。

10.根据权利要求1所述的方法，其中，所述窦腔造影细分步骤是近似的，并包括窦腔造影的截断(1110)、移位(1112、1120)和分样(1116)，从而使每个所生成的子窦腔造影对应于在整体坐标系中的一子图像。

11.根据权利要求10的方法，其中，所述近似细分步骤还包括加权(1114、1118)。

12.根据权利要求10的方法，其中，所述细分被递归地执行，并且其中所述在递归的后续级中的移位操作被组合以减少计算成本。

13.根据权利要求1所述的方法，其中，细分被在空间上各向异性地和不均匀地执行，从而在给定递归级中的子图像不需要是相同的形状或尺寸。

14.根据权利要求13所述的方法，其中，根据源轨道的对称性质选择分割，以达到在计算成本和图像精确度之间所期望的平衡。

15.根据权利要求1所述的方法，其中，该方法被单独地应用到从多个源轨道中获得的多个窦腔造影，以生成可以被组合去生成最后图像的多个图像。

16.根据权利要求1所述的方法，其中，精确细分步骤被配置成使在子图像的窦腔造影或子窦腔造影中的投影的截断操作在一个或多个投影变成该子图像可用之后开始。

17.根据权利要求1所述的方法，其中，细分步骤被配置成使所述子图像的窦腔造影的细分在一个或多个投影变成为该子图像可用之后开始。

18.根据权利要求1所述的方法，当它们变成可用时，将它们依时间连续地应用到多个发散波束窦腔造影上。

19.根据权利要求1所述的方法，其中，发散波束窦腔造影是在会聚波束几何形状中获得的。

20.用于生成一目标的电子图像(36)的装置，包括：

用于扫描该目标以生成表示该目标图像的数据的装置(12)；

用于操纵所述数据以生成经处理的窦腔造影(22)的装置(16、20)；

用于将所述窦腔造影(22)细分为多个子窦腔造影(1010、1012、1014、1016)的装置(24)；

用于在整体坐标系(1018、1020、1022、1024)中反向投影每个所述子窦腔造影，以便在整体坐标系中的适当位置处生成多个相应的子图像(f₁、f₂、f₃、f₄)的装置(24)；

用于会聚所述子图像以生成电子图像的装置(24)；以及

用于存储和/或显示和/或分析电子图像的装置(38)。

21.根据权利要求20所述的装置，其中，所述窦腔造影细分是近似的(1210、1212、1214、1216)。

22.根据权利要求20所述的装置，其中，所述窦腔造影以递归方式被细分成多个子窦腔造影，直到每个子窦腔造影都对应于一期望尺寸的子图像为止，并且所述用于细分的装置执行期望数量的精确细分(1510、1512、1514、1516)和期望数量的近似细分(1518、1520、1522、1524)。

23.根据权利要求22所述的装置，其中，所述精确细分包括没有移位的截断。

24.根据权利要求22所述的装置，其中，所述子窦腔造影对应于在尺寸上与一个像素或一个体元一样小的子图像。

25.根据权利要求20所述的装置，其中，所述用于会聚的装置以递归方式操作。

26.根据权利要求20所述的装置，其中，所述电子图像是断层摄影图像。

27.根据权利要求20所述的装置，其中，所述用于操纵的装置执行在数据坐标中(16)的过采样以改善电子图像的精确度。

28.根据权利要求20所述的装置，其中，所述窦腔造影细分步骤是近似的，并包括窦腔造影的截断(1110)、移位(1112、1120)和分样(1116)，从而使每个生成的子窦腔造影对应于整体坐标系中的一个子图像。

29.根据权利要求28所述的装置，其中，所述近似细分步骤还可以包括加权(1114、1118)。

30.一种用于从电子图像生成发散波束窦腔造影(34)的方法，包括步骤：

将该图像划分(1710、1712、1714、1716)为多个子图像(c)；

计算(1718、1720、1722、1724)在整体坐标系中每个所述子图像的子窦腔造影(g₁、g₂、g₃、g₄)；以及

会聚所述子窦腔造影(1726)以创建窦腔造影(34)。

31.根据权利要求30所述的方法，其中，所述子窦腔造影的计算是近似的(1926、1928、1930、1932)。

32.根据权利要求30所述的方法，其中，所述细分是以递归方式执行的，直到每个所述子图像具有期望的尺寸为止，并且所述子窦腔造影的计算在期望数量的递归级中是近似的，而在剩余的递归级中是精确的。

33.根据权利要求30所述的过程，其中，所述会聚以递归方式执行。

34.根据权利要求30所述的过程，其中，所述精确会聚是通过添加子窦腔造影而没有移位执行的。

35.根据权利要求30所述的过程，其中，所述子窦腔造影的计算是近似的，并且通过移位(1810、1818)、内插(1814)和加入子窦腔造影来执行，并且每个移位和内插的子窦腔造影都对应于整体坐标系中的一个子图像。

36.根据权利要求35所述的过程，其中，近似计算所述子窦腔造影的步骤还包括加权(1812、1816)。

37.用于创建一目标的图像(38)的装置，包括：

扫描器(12)，其从目标生成数据；

处理器(16)，用于创建该图像的至少一个发散波束投影；

装置(24)，用于根据至少一个投影重现或在迭代再现装置的情况下反向投影一个图像；

用于检测在由所述再现装置生成的图像中的错误的装置(28)；

用于在纠错之后再投影图像、对窦腔造影执行附加校正并将校正后的窦腔造影馈送到所述再现或反向投影装置(24)的装置(32)；以及

用于显示由再现装置在错误被校正以后创建的图像的装置(38)，

其中，所述再投影装置将图像划分为多个子图像(f₁、f₂、f₃、f₄)，计算在整体坐标系中(g₁、g₂、g₃、g₄)所述子图像的子窦腔造影，并会聚所述子窦腔造影(1726)以生成图像窦腔造影(g)。

38.根据权利要求37所述的装置，其中，所述子窦腔造影的计算是近似地。

39.根据权利要求37所述的装置，其中，所述细分以递归方式执行，直到每个所述子图像具有期望的尺寸为止，并且所述子窦腔造影的计算在期望数量的递归的级中是近似的，而在剩余的递归的级中是精确度。

40.根据权利要求37所述的装置，其中，所述会聚是以递归方式执行的。

41.根据权利要求37所述的装置，其中，所述精确会聚是通过添加子窦腔造影而没有移位执行的，。

42.根据权利要求37所述的装置，其中，所述子窦腔造影的计算是近似的，并且是通过移位、内插和加入子窦腔造影执行的，并且每个经移位和内插的子窦腔造影对应于整体坐标系中的一个子图像。

Images