CN1867926A - 快速分层x线断层造影方法和装置 - Google Patents

Abstract

通过反投射(100)选定投射以生成中间图像(I₁，m)，并执行数字图像坐标变换(102)和/或对选定中间图像重新采样(图31的186、192、196)，像素图像f从投射(q₁...q_p)创建。数字图像坐标变换(102)被选为说明中间图像的构成投射的视角、并说明其傅立叶特征，从而中间图像可由稀疏样本准确地表示。结果中间图像被聚合成子集(104)，且该过程以递归方式重复直到已处理和聚合了足以形成像素图像f的投射和中间图像。数字图像坐标变换可包括旋转(图18的102)、切变(图10B的120、122)、拉伸、收缩(109)等。重新采样可包括上采样(101、106)、下采样(109)等。通过执行数字图像坐标变换(202)和/或重新采样(204)和/或抽选(图32的204；图33的212)、和/或重新投射最终中间图像(208)，投射(图32的P_θ1...P_θ18)可从像素图像(f)创建。

Description

快速分层X线断层造影方法和装置

相关申请

本申请是序列号为60/501,350的于2003年9月9日提交的临时申请的后续部分，该临时申请通过全部引用结合于此。

技术领域

本发明涉及X线断层造影，尤其涉及用于从投射创建像素图像的方法和装置。

背景技术

断层图像重建是一种众所周知的技术，几乎是所有诊断显像医疗器械的基础，包括：计算机X线断层造影(CT)、正电子断层扫描仪(PET)、单光子断层扫描仪(SPET)、以及用于磁振造影(MRI)的某些获得方法。它还发现在关于非破坏性评价(NDE)、安全扫描的制造业、合成孔径雷达(SAR)、射电天文学、地球物理学和其它区域中的应用。

断层图像数据有若干种主要格式：(i)平行光束，其中沿多组平行线执行线积分；(ii)发散光束，其中沿多组发散为扇形或圆锥形的线执行线积分；(iii)曲线，其中沿多组诸如圆、椭圆或其它闭合或开曲线的曲线执行积分。断层图像重建的一个问题是从其线积分投射集中重建2D或3D图像。断层图像重建的另一个问题是从其面积分投射集(即其在曲面族上的积分)中重建3D图像。例如，3D Radon变换涉及图像在距离原点各个方位和距离的2D平面族上的积分。断层图像重建的一些问题、以及一些重建方法在标准参考书中描述，诸如：F.Natterer的“TheMathematics of Computerized Tomography”(计算机X线断层造影的数学)，1986年John Wiley出版社在Chichester出版；F.Natterer和F.Wubbeling的“MathematicalMethods in Image Reconstruction”(图像重建中的数学方法)，2001年在费城工业与应用数学学会发表；A.C.Kak和M.Soaney的“Principles of ComputerizedTomographic Imaging”(计算机X线断层造影的原理)，1988年IEEE出版社在纽约出版；以及S.R.Deans的“The Radon Transform and Some of Its Applications”(“Radon变换及其一些应用”)，1983年Wiley出版社在纽约出版。

选择用于断层图像重建的方法是过滤反投射(FBP)或回旋反投射(CBP)，使用不加权的(平行光束或Radon变换情形中)或加权的(大多数其它情形中)反投射步骤。该步骤是技术中的计算瓶颈，其中计算需求规模为：对2D的N×N像素图像为N³，对3D的N×N×N体素图像为至少N⁴。因而，从N到2N倍增图像的分辨率导致计算量增加约8倍(或者3D中为16倍)。当计算机已变得快得多时，同时出现了能够实时收集更大量数据的新技术(例如使用多行检测器的心脏成像、介入成像)、以及3D采集几何结构的盛行，对快速重建技术有越来越多的需求。快速重建可加快成像过程、降低专用图像重建计算机的成本，或两者皆可。

反投射的二元操作是再投射(reprojection)，它是计算电子存储图像的投射的过程。该过程也在断层图像重建中起了基本作用。反投射和再投射的组合也可用来构建对螺旋型锥面光束几何结构中长对象问题的快速重建算法，这对人体的实际3D成像是关键的。此外，在各种应用中，使用迭代重建算法是有利甚至必要的，其中对单个图像的重构执行反投射和再投射步骤若干次。加快反投射和再投射步骤将确定这些迭代方法的经济可行性。这些年来已提出了若干种加快重建的方法。例如，Brandt等人的美国专利5,778,038描述了使用多层分解、每个阶段产生覆盖整个视野且分辨率递增的图像的2D平行光束断层造影方法。Nillson等人的美国专利6,151,377揭示了其它分层反投射方法。尽管这些系统都具有价值，但仍然需要产生更准确图像、提供准确度和速度之间的更大灵活性的方法和装置。

因此，本发明的一个目的是提供用于计算机X线断层造影(CT)扫描的新颖和经改进的方法和装置。

另一个目的是提供产生更准确图像、并提供准确度和速度之间的更大灵活性的、用于CT扫描的新颖和经改进的方法和装置。

发明内容

这些目的由本发明实现或超越。通过反投射选定投射来产生中间图像、并执行数字图像坐标变换和/或对中间图像重新采样，像素图像得以根据投射创建。数字图像坐标变换被选择用来说明中间图像的组成投射的视角、及其傅立叶特征，从而中间图像可通过稀疏样本来准确地表示。结果的中间图像被聚合成子集，且该过程以递归方式重复，直到已处理和聚合了足够的投射和中间图像，以形成像素图像。

数字图像坐标变换可包括旋转、切变、拉伸、收缩等。重新采样可包括上采样、下采样等。

通过执行数字图像坐标变换和/或重新采样和/或抽选，并重新投射最后的中间图像，投射可得以从像素图像创建。

附图说明

通过参阅附图和本发明一实施例的以下描述，本发明的上述和其它特征以及获得它们的方式将变得显而易见，且本发明本身将得到最佳的理解，在附图中：

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

图2A、2B和2C是用于本发明一些实施例中的采样图(sampling pattern)的示图；

图3A、3B和3C是用于本发明一些实施例中的其它采样图；

图4是示出一种已知反投射方法的示图；

图5A是示出本发明一实施例的算法的示图；

图5B是示出在图5A的实施例中产生中间图像的方式的示图；

图6是示出用于产生图5A的中间图像的傅立叶特征的示图；

图7A、7B和7C是示出当坐标变换是数字图像旋转时用于图5A中示出的反投射算法的中间图像的傅立叶支持的示图；

图8是示出用于本发明的另一实施例中的算法的示图；

图9是示出图8算法中光谱支持的演变的示图，其中框(1……9)对应于图8中的相应各点；

图10A是描述用于本发明实施例的一种算法的示图；

图10B示出用于图10A的实施例中的图像坐标变换；

图11A是示出切变比例反投射的示图，而图11B是示出分层切变比例反投射的示图；

图12A和12B是示出图像在中间图像的光谱支持上的切变效果的示图；

图13是示出本发明另一个实施例的算法的示图；

图14是找出最优切变因子的算法；

图15示出用于本发明另一实施例的算法；

图16示出用于本发明又一实施例的算法；

图17是示出具有圆形扫描轨迹的通用扇形光束几何结构的示图；

图18示出用于本发明再一实施例的算法；

图19示出图18中所述的算法中中间图像的加权；

图20A、20B和20C示出图18的第二分层结构层的采样点；

图21A和21B是用于图18算法中的采样图的示图；

图22是示出使用图20A-22C中所示方法获得的原始交点的示图；

图23A、23B、23C和23D示出用于图18中使用的上采样的旋转采样点；

图24示出在由图18的算法产生的中间图像的一点上的本地光谱支持；

图25是在图18的算法中使用的不统一采样图的示图；

图26A和26B示出用于再次采样和坐标变换的采样图；

图27A和27B示出可用于本发明的选择性采样方案；

图28包括用于本发明的发散光束的两个示图；

图29A是示出锥面光束的示图，而图29B是示出再次采样的示图；

图30是用于再次采样投射的算法的示图；

图31是用于在本发明一实施例中再次采样的另一种算法；

图32是用于快速分层再投射的算法的示图；

图33是用于快速分层再投射的另一种算法的示图；

图34是示出使用本发明的实验结果的示图；

图35包括用本发明产生的样本图像；

图36A是使用常规算法的重建图像的显示，而图36B显示使用本发明的快速算法获得的结果；以及

图37A和37B是与常规算法相类似的算法的点分布函数的示图，而图37C显示本发明快速算法的点分布函数。

具体实施方式

符号和字体

以下数学符号和字体系统将用来提高明确度。

空间域中的函数由小写字母表示(例如f(x))，而其傅立叶变换由大写字母表示(F(ω))。

两个变量函数的索引取决于方便性可不同地表示。函数f的以下三种表示是等同的：f(x₁，x₂)，f( x)，以及 $f\left(\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]\right).$

连续域和离散域函数分别通过其索引中括号的样式来区分：f(x₁，x₂)是两个连续变量的函数(即， )，而f[m₁，m₂]是f( x)的采样版本，因此是2-D数组。

线性算子及其相应矩阵通过字体样式来区分。假设 是矩阵，则A是其相关联的线性算子。例如，如果A是坐标变换，则g( x)＝(Af)( x)＝f( Ax)。

相同的算子有时在框图内外会不同地表示。在框图外它可表示为A(α)，在框图内它表示为A_a。

硬件概述

本发明在包括CT扫描仪的各种成像装置中都有应用。典型的成像装置10(图11)包括扫描仪12，它从诸如头的对象中采集数据，并向投射预处理器16发送与例如具有不同光束几何形状的线积分投射14相对应的原始数据。该投射预处理器16对数据应用各种转换、归一化、和校正、以及加权和可以是漂移变化的滤波。投射预处理器16的输出是预处理后投射的集合，下文中简称为投射，也称为sinogram₁(窦腔X线造影)18，它可提供给窦腔X线造影更新处理器20。该窦腔X线造影更新处理器20使用来自sinogram₂34的信息，可能更改所输入的sinogram₁18，例如校正包括射线束硬化的各种人工因素、或作为多步骤或迭代重建过程的一部分。

窦腔X线造影更新处理器20的输出是sinogram₃22，它被输入给快速反投射处理器24。该快速反投射处理器24通常是经编程和/或接线用以执行这里所述算法的任何适当类型的计算机或专用硬件，或其任意组合。

该快速反投射处理器24的输出是电子图像₁26，它被输入给图像调整处理器28。该图像调整处理器28执行电子图像的必要后处理，可能包括非源自脑中的电波图像、或用于在多步骤或迭代重建过程中作进一步处理的图像的标识和提取。

在需要时，图像调整处理器28可产生提供给快速再投射处理器32的电子图像₂30。该快速再投射处理器32通常是经编程和/或接线用以执行这里所述算法的任何适当类型的计算机或专用硬件，或其任意组合。需要时该处理器能共享反投射处理器24采用的同一计算机和硬件的运用。

快速再投射处理器32的输出是反馈回窦腔X线造影更新处理器20的sinogram₂34。反投射/再投射过程可继续，直到获得适当结果。尽管再投射并不总是需要的，但在许多情形中都有帮助。

当电子图像₁26适当时，图像调整处理器28产生提供给存储/分析/显示设备38的电子图像₃36。预期电子图像₃36可被存储在计算机存储器中，和/或为例如异常或危险物质对其作电子分析，和/或显示，和/或用一些可查看形式打印。

本发明的反投射和再投射方法的概述

本发明的反投射方法使用各种技术来创建由像素(图片元素)和/或体素(3D图片元素)(下文中统称为像素)组成的图像，以下将用一般方法进行介绍。

本说明书使用常常在多维信号处理中使用的术语和过程，例如1983年Prentice-Hall出版社在Englewood Cliffs出版的D.Dudgeon和R.Mersereau的“Multidimensional Digital Signal Processing”(多维数字信号处理)中所述。本发明说明书中的一些术语用于以下上下文中。术语“采样图”指其位置相对于坐标系统限定的空间点集。采样图的示例如图2A-2C以及3A-3C所示。笛卡尔采样图指由两组相互垂直的平行线的交点形成的点集。术语连续图像指在坐标系统上定义的函数，例如f(x，y)和f(x，y，z)分别是2D和3D函数。数字图像是采样图上连续图像的各个值的阵列。更广泛地，连续图像可由数字阵列表示，这些数字用作参照诸如样条的一些基集的序列扩展中的系数。零阶样条的笛卡尔积产生用于将数字图像显示为持续图像的熟悉的方形像素形式。在下文中，该数字阵列也被称为数字图像。所有存储在数字计算设备中的图像必须是数字的。为简短起见，数字图像和连续图像两者在下文中都经常被简称为图像，其意思可从上下文中推断。使用该术语，像素图像是与格子采样图相对应的数字图像，即均匀分隔的周期性图案、通常并非必然是笛卡尔。

如果一种采样图产生比连续图像少的样本总量，则它被称为比另一种采样图稀疏。通常较为稀疏的采样图将具有较低的采样密度。过采样指使用比以期望准确度表示连续图像所必需的样本更多的样本。相应的数字图像将称为被过采样。给定与一采样图的连续图像相对应的数字图像，用期望准确度产生与不同采样图上的相同连续图像相对应的新数字图像的过程称为数字图像重新采样。上采样和下采样分别是对较密集或较稀疏采样图重新采样的特定情形。此外，上采样或下采样1倍不涉及数字图像中的变化，被视为是上采样或下采样的一种形式。数字图像附加指在相对于基的扩展情形中，相对于同一采样图或同一基定义的数字图像的逐点附加。低维数字过滤指多维阵列仅在其多个维的子集上的数字过滤，例如2D矩形数字图像的每个列的单独1D过滤。

连续图像的坐标变换指诸如旋转、切变、拉伸、或收缩的运算。要定义一数字坐标变换，考虑通过坐标变换关联的两个连续图像、以及表示相对于公共采样图的连续图像的相应数字图像。从另一个数字图像产生一个数字图像的过程称为数字图像坐标变换。这可通过数字过滤来实现，即通过对数字图像的离散指数运算。各特定示例包括(但不限于)数字图像旋转、数字图像切变、数字图像拉伸、数字图像收缩、以及这些运算的组合。用于执行数字图像坐标变换的各种方法是众所周知的，例如，如IEEE Transactions Image Processing 1995年第4卷1371-1381页的M.Unser、P.Thevenaz、和L.Yaroslavsky的“Convolution-based interpolation for fast，high-quality rotation ofimages”(用于图像的快速、高质量旋转的基于回旋内插)中所述。

一些数字图像坐标变换在图2A-2C和3A-3C中示出。图2A示出连续图像(矩形)的轮廓、以及表示连续图像的数字图像的采样图。较大点上连续图像的值包括在数字图像中。图2-B和2-C分别示出对同一采样图的旋转后连续图像和切变后连续图像，其中较大点示出包括在图2-A数字图像的旋转后/切变后数字版本中的值。

图3A也示出一连续图像以及定义数字图像的采样图。图3B示出在x和y维上拉伸某常数倍后的连续图像。该数字上拉伸后的图像由较大点上拉伸后连续图像的值来定义。图3C示出与图3A相同的连续图像，但其中采样图要更密一定拉伸倍。图3C中的对应的数字图像将与图3B中的数字图像相同。更一般地，对于具有一定规则度的采样图，数字图像拉伸或收缩对数字图像上采样或下采样是等同的。

注意，某些坐标变换的应用不会使数字图像改变，诸如旋转0度、切变参数为零地切变、或拉伸或收缩1倍，因此在过程的描述中可包括在内或略去，而不会改变结果。

连续图像的傅立叶变换使人能通过采样理论确定采样图的属性，从而相应数字图像将连续图像表示到期望准确度，如例如1983年Prentice-Hall出版社在Englewood Cliffs出版的D.Dudgeon和R.Mersereau的“Multidimensional DigitalSignal Processing”(多维数字信号处理)中所述。类似地，数字图像的离散时间傅立叶变换(DTFT)使人能确定什么数字图像坐标变换将产生将变换后的连续图像表示到期望准确度的数字图像。连续图像的傅立叶变换的相关属性、以及数字图像的DTFT将统称为傅立叶特征。

加权反投射运算需要对每个投射加权取决于反投射所产生像素的位置的权重。取决于例如获取投射的源的位置，不同的权重可用于不同的投射，如1988年IEEE出版社在纽约出版的A.C.Kak和M.Slaney的“Principles of ComputerizedTomographic Imaging”(计算机X线断层造影的成像原理)中所述。作为一特定情形，加权倍数可等于1，它对应于无加权但的确是加权倍数。除非特别指出，经加权和未经加权的反投射将统称为反投射。

有了这样的背景技术信息，将描述本发明的若干实施例。

反投射处理器28(图1)被编程为执行用来实践本发明的算法。这些算法将详细地进行讨论，但首先将以更普通的术语进行描述。反投射过程中的步骤用标号在框图中示出。

图4示出基于旋转的反投射，作为通过步骤50中以0度反投射各个投射q₁...q_p而形成的旋转后中间图像的和。反投射产生在52进行坐标变换的图像，并在54聚合以生成图像 该结构单独地等同于常规反投射，并在运算计数中不提供缩减。然而，它用作引入本发明的一些快速分层反投射方法的进阶。

图5A示出用于从多个投射q₁...q_p中创建像素图像

的一种分层反投射方法。每个投射q_m都在100反投射以产生多个中间图像I_1，1...I_1，p这是分层反投射的零层或准备层。在102对选定中间图像执行数字图像坐标变换。变换后中间图像的子集(在此为成对)在104聚合以产生聚合中间图像I_2，I...I_2，p/2。这是分层反投射的第一层。第一层的聚合中间图像用作输入到分层反投射的下一层的新的中间图像。在106将数字图像坐标变换应用于选定中间图像、以及在108聚合选定中间图像以产生新的中间图像的过程继续，直到已处理完所有中间图像、并在116聚合成最终图像

在需要时，可组合跨一些层和这些层内的运算。例如，对于一些投射、以及来自通过在两个或多个(对于图5A中的实施例正好为2)选定投射q_p的选定视角上的反投射112产生的集I_2，1...I_2，p/2的相应中间图像，第0和1层可组合，如图5B所示。或者，一些初始中间图像可通过不涉及显式反投射的等同过程生成，诸如直接傅立叶重建方法，如2001年在费城召开的工业和应用数学学会上发表的F.Natterer和F.Wubbeling的“Mathematical Methods and Image Reconstruction”(数学方法和图像重建)中所述。

各个数字图像坐标变换的参数被选择，用以说明中间图像的构成投射的视角，并说明中间图像的傅立叶特征，从而中间图像的聚合可由稀疏样本准确地表示，如以下将解释的。这些傅立叶特征集中在中间图像的必要光谱支持上，即傅立叶(频率)域中的区域，其中中间图像的傅立叶变换远不同于零。采样理论教导：连续图像的光谱支持确定产生表示基础连续图像的、以及从中能可靠地重建连续图像的数字图像的采样图的特性。特别地，图6示出分层算法中中间图像 的典型楔形光谱支持。

图7(A)和7(B)示出在坐标变换被选为数字图像旋转时，图5(A)所示的二进制分层反投射算法中中间图像的傅立叶支持。图7(A)示出虚拟中间图像 I的傅立叶域支持。虚拟图像

由包括在相应图像I_l，m中的多个投射的反投射组成。图7(B)示出通过选择坐标变换的参数来说明中间图像的构成投射的视角以及中间图像的傅立叶特征，中间图像的纵向带宽(虚线矩形的高度)可最小化，从而允许稀疏采样并降低计算需求。图7(C)示出I_l，m、I_l-1，2m和I_l-1，2m-1的空间域采样方案。采样点在空间域中虚线的交点处。

图8示出本发明的另一个实施例，它执行中间图像的三重聚合，其中P＝2×3^L个投射。为了具体说明，示出P＝18个投射的情形。象图5(A)中一样，投射q_θ1，q_θ2，…，q_θ18在100反投射以产生多个中间图像I_1，1 ^d…I_1，18 ^d。这是分层反投射的第零层。

各投射以3个为一组，且每个组中的选定投射(图8的3个中的2个)在102上进行数字图像的坐标变换。变换后中间图像I₁ ^d的子集(图8中的三元组)在104上聚合，以生成聚合中间图像I_2，1 ^d…I_2，6 ^d。这是分层反投射的第一层(1＝1)。

在分层反投射的第二层上(1＝2)，选定聚合中间图像I_2，m ^d在106进行由沿y坐标拉伸或上采样组成的坐标变换，而在108对选定的聚合中间图像进行另一次坐标变换K_2，m。在此，变换后的中间图像在110还是聚合成三个一组，以生成中间图像I_3，1 ^d，I_3，2 ^d。来自层2的选定中间图像在112进行数字图像坐标变换，并在114聚合以生成下一层的中间图像I_4，1 ^d。取决于投射的数量，该过程重复到在116生成图像

所必需的程度。在112最后一层上由K标注的数字图像坐标变换是一数字图像旋转，且在前面各层102和108上的坐标变换也可选为数字图像旋转。在此，各个数字图像坐标变换的参数被选为说明中间图像的构成投射的视角、和中间图像的傅立叶特征，从而中间图像的聚合可由稀疏样本来准确地表示，如以下将解释的。

图9示出图8的三重反投射算法中光谱支持的进展。编号(1)，...，(9)对应于图8算法框图中的对应点。所示光谱是数字图像I_l，m的离散时间傅立叶变换(DTFT)。

图10(A)描述了本发明另一个实施例的算法，该算法使用三重基于两切变的分层反投射。图10(A)的实施例与图8的实施例相似，且来自图8的标号也在适用时使用。如在图8中，示出P＝2×3^L个投射，且L＝2的情形。图10(a)的实施例与图8的不同之处有两个方面。首先，沿x坐标的附加上采样步骤被包括在101第一层的坐标变换中。其次，在102和108由K标注的数字图像坐标变换的至少一部分由两个数字图像切变运算序列组成，第一个(120)沿y坐标(122)，第二个沿x坐标，如图10(b)所示。

图11A示出切变缩放(shear-scale)反投射，而图11B示出等效的分层的切变缩放反投射。在图11A中，多个投射q₁，…，q₄在140反投射，并在142进行切变缩放坐标变换。最后的中间图像在144聚合。

在图11B中，多个投射q₁，…，q₄在146反投射，在148进行切变缩放坐标变换。并在150聚合成子集。通过聚合产生的中间图像在152再次进行切变缩放变换，且最后的中间图像在154聚合。该过程以递归方式继续，直到所有的投射和中间图像都已得到处理并聚合，用以形成图像

在此，数字切变变换和中间图像采样的参数被选为说明视角和傅立叶特征，以便于降低采样需求和所需计算量。

图12示出图像切变对中间图像的光谱支持的效果。图12A示出当投射应出现在最终图像中时其某个子集的光谱支持。图12(B)示出由相同投射组成的中间图像的光谱支持，其中坐标变换参数被选为使纵向最高角频率最小化。

图13示出本发明又一实施例—三重分层切变缩放反投射的算法。投射q_θ1，…，q_θ18在160反投射并在162进行切变缩放数字图像坐标变换。结果图像的子集在164聚合，且结果中间图像在166和168进行上采样，并在168进行数字切变坐标变换。这些图像的子集在170聚合，且选定结果图像在172进行其它坐标变换。该过程继续直到所有投射和中间图像已处理并在174聚合，生成图像

在172所示的最终坐标变换K_-π/2仅涉及当像素图像

的采样图是笛卡尔时像素的重新排列。

图14是使用图12所示中间图像的傅立叶属性的属性，寻找先前所述实施例的坐标变换参数的一种算法。

图15示出本发明的再一实施例—过采样三重两切变分层反投射的算法。图15的实施例与图10A的实施例相似，且来自图10A的标号在适用处使用。然而，除了图10A的各个步骤之外，图15的实施例包括最后步骤之前一个步骤中的下采样步骤109，它在图15中为第二层。在此，各个数字图像坐标变换的参数被选为说明中间图像的构成投射的视角和中间图像的傅立叶特征，从而中间图像的聚合可由稀疏样本来准确表示。然而，一定程度的过采样用来改进随后处理的准确度，如将作解释的。如果需要，则为了改进计算效率，下采样步骤109可与第二个x轴切变组合，使过程108和109一起被图15(B)中所示的过程所替换，其中第二个x轴切变包括产生切变缩放变换的两切变数字图像坐标变换108(如图10B所示)。

图16示出本发明的又一实施例—过采样三重分层的切变缩放反投射的算法。图16的实施例与图13的实施例相似，且来自图13的标号在适用处使用。然而，图16包括中间层中上采样161的附加步骤，以及在最后层之前一层中的下采样169，其中中间图像被分别上采样和下采样。与图15的实施例相似，在图16的实施例中，当沿x坐标的下采样步骤169↓_xU_l，m在x切变步骤168 S_x ^αl，m之后时，为了计算效率两者可组合成单个数字图像切变缩放。

用于快速分层反投射算法的非笛卡尔采样方案

图8、10和13中所示的分层反投射的不同实施例中的中间图像具有便于有效非笛卡尔采样的特殊光谱支持。特别地，第1层中的基础连续域图像I_l，m在傅立叶空间中占据一楔形，如例如在图6、7、9和12中可见。多维采样理论告诉我们，对于具有这种光谱支持的图像，在笛卡尔栅格上的采样不如在适当非笛卡尔栅格上的采样有效。非笛卡尔采样可通过将基带频谱的副本更紧地封装于傅立叶平面而改进采样效率。对于2D采样的解释参见[？]。特别地，梅花形采样方案减少中间图像的大小，因此使该算法的计算成本降低几乎2倍。对周期性的非笛卡尔采样图的数字图像坐标变换可使用一维移位不变滤波器来有效地执行。类似地，所有先前描述的用于选择数字图像变换的参数的方法也在这种采样图情形中应用。因此，前述所有实施例扩展到使用周期性非笛卡尔采样图的各个实施例。

已经描述的本发明各个实施例的变体可应用于各种3D投射形式的3D反投射，例如正电子断层扫描仪中产生的X线变换；以及在磁振造影中产生的3D Radon变换。

三维(3D)X线变换是3D中沿平行线集在各个方位上的积分的集合。每个3D X线投射是可由线定向的3D角度表征的两维函数。3D X线数据的分层反投射的框图与前述的相似，诸如图8、10、13、15或16。该情形中的中间图像是三维的，而非前述示例中的二维。每个中间图像在各构成投射的3D角度的平均值方向上最稀疏的3D采样图上采样。随着该算法的展开，沿该稀疏采样(缓慢)方向上样本的密度增大，以适应该方向上的渐增带宽，如3D X线投射的傅立叶分析所述。因此在该算法的每个阶段，在聚合图像之前每个图像都必须沿该缓慢方向上采样。3D实施例中可用的额外维数还提供可用于算法中各个步骤的更多坐标变换，诸如3D中的旋转。像2D情形中一样，这些数字图像坐标变换的参数可选为说明构成视角和中间图像的傅立叶特征。这些数字图像坐标变换可被分解成一维运算序列，诸如如前所述的切变和切变缩放。像2D情形中一样，可迫使三维的任何子集中的过采样，以改进图像质量。

3D Radon变换投射是一维函数——沿平行2D平面集的以平面离原点的距离为参数的积分的集合。投射视角是与平面集垂直的矢量的3D方位角的视角。3DRadon投射的分层反投射的框图象图8、10、13、15或16中一样。在第一层中多个投射被Radon反投射到3D图像域上。这些图像是沿两维恒定的，且因此仅需要在垂直于该两个恒定方向的方向上采样。当多个组的3D Radon投射组合时，聚合图像的带宽可取决于这些构成投射的视角在一个或两个维上增加。因此，有必要在坐标变换中间图像并添加到其它中间图像之前，在可能两个维中上采样该中间图像。像2D情形中一样，坐标变换可单独执行、可与上采样运算组合、并且过采样可在中间图像上实施。

在本发明的各个实施例中，数字图像坐标变换、以及下采样或上采样运算可通过一系列较低(一)维线性数字过滤运算来执行。此外，当所使用的采样图具有一些周期性时，这些数字滤波器可以是移位不变的，如将要更详细描述的。为了计算效率，数字滤波器可使用递归滤波器结构、或使用众所周知的FFT来实现。一种确定数字滤波器的较佳值的方法是使用样条内插理论，如IEEE Transactions onPattern Analysis and Machine Intelligence 1991年第13卷第277-285页的M.Unser、A.Aldroubi和M.Eden的“Fast B-Spline transforms for continuous imagerepresentation and interpolation”(用于连续图像表示和内插的快速B样条变换)；IEEE Transactions Signal Processing 1993年第41卷第821-832页的M.Unser、A.Aldroubi和M.Eden的“B-spline signal processing：Part I-theory”(B样条信号处理：部分I-理论)；以及IEEE Transactions Signal Processing 1993年第41卷第834-848页的M.Unser、A.Aldroubi和M.Eden的“B-spline signal processing：Part II-efficientdesign and applications”(B样条信号处理：部分II-有效设计和应用)中所述。

本发明的分层反投射方法可应用于广泛范围的X线断层成像几何结构，诸如具有任意源轨迹的包括扇形光束和锥面光束的发散光束几何结构。具有圆形扫描轨迹的通用扇形光束几何结构在图17中示出。射线源绕半径为R的图像在半径为D的圆形扫描轨迹上移动。源角度β上的扇形光束投射对应于沿以扇形角度γ为参数的射线集的线积分度量。T_ST是沿源射线到图像盘的最近边缘的距离，T_END是沿源射线到图像盘的最远边缘的距离。

图18示出用于本发明一实施例的一种算法——分层的三重基于旋转反投射，可应用于扇形光束加权反投射。图18的实施例与图8的实施例相似，且图8的标号在适用时使用。图18的实施例与图8所示实施例的不同之处在于若干方面。第一，示出P＝4×3^L个投射，且L＝2的情形。第二，在第0层上，初始中间图像I_l，m ^d通过在99的标示为W的加权反投射从扇形光束投射中生成。第三，在最后阶段，聚合4个而不是2个中间图像。第四，在分层结构的大多数较早层中使用的采样图最好选为非笛卡尔，如将要解释的。在此，如果笛卡尔采样图用于倒数第二层中的中间图像I_3，m ^d，则最后的数字图像坐标变换同样仅需要重新排列图像像素。此外，像图5(B)中一样，分层结构的第0层和第1层中的运算可组合，从而中间图像I_2，m ^d可通过其三个构成投射的直接加权扇形光束的反投射、或通过其它手段生成。

在选定数量的层中，通过将附加加权步骤包括在上采样步骤和数字图像旋转(图18中的步骤106和108)的级联前后(如图19所示)，来更改图18的实施例是有利的。中间图像I_2，m分别在步骤106和108前后在180和182通过在空间上更改权重来加权。如将要解释的，该加权可用来减少中间图像的采样需求，从而减少计算量。

中间图像的采样方案影响算法的性能。期望采样方案是使用最少样本但不丢失图像信息的方案。在此，各个数字图像坐标变换和重新采样运算的参数可选为说明中间图像的构成投射的视角、以及中间图像的傅立叶特征，从而中间图像的聚合可通过稀疏样本来准确地表示，如将要解释的。如将要解释的，用于选择这些参数的选择性方法基于特殊射线或曲线集的交点。

图20A-20C示出在源角度以间隔Δ_β均匀隔开的情形中，中间图像通过图18所示三重分层结构算法的各个层次的进度。示出了算法中组成中间图像的多个投射的扇形。在零层上，每个中间图像I_1，m ^d由具有方位为β＝0的扇形的单个投射组成，如图20(A)所示。在递归的第一层之后，每个中间图像由三个投射组成，其中扇形由图20(B)所示。在第二层之后，每个图像由9个反投射扇形组成，如图20(C)所示。

算法中用于对具有如图20(C)所示构成扇形的中间图像I_3，m选择坐标变换和采样/重新采样图的基于交点方法如图21A和21B所示。I_3，m的采样图由位于中央构成扇形的点组成，与图20(A)中所示的一致。如图21(A)和21(B)所示，I_3，m的采样点由中央扇形的一半与中央扇形两侧上的极值构成扇形的交点所确定。

通过应用两个约束条件来更改结果采样图，以改进反投射的准确度并减少计算需求是有利的。首先，沿射线的样本密度受到限制，不超过最终图像的期望采样密度。其次，采样图被迫使在每条射线上包含每一侧的图像盘外的至少一个样本。图22中的曲线显示采样点沿中央扇形的扇形角度为γ_r的特定射线上的位置。采样点位于距离γ’的固定间隔上，其中γ’是图21(a)、(b)和(c)所示极值扇形之一的扇形角度。落在图22中连续曲线上的各个点是使用图21A和21B中所示的方法获得的原始交点。虚曲线上的各个点是使用以上两个约束条件更改的。图23D示出使用该更改交点方法获得的采样图的一个示例。所示扇形是已为之生成该采样图的中间图像的中央扇形。

当在本发明前述实施例的情形中时，将重新采样和数字图像旋转运算分解成一维运算系列是有利的。图18中标为↑U的框表示将对中央扇形采样的中间图像上采样成同一扇形上采样点的较小集。这可通过沿扇形的每条射线的单独上采样来实现。对于图18中标有↑U和K的上采样和旋转的级联，联合分解成1D运算是方便的，如图23(A)-23(D)所示。在4个画面的每一个中，两个虚线圆表示半径为D的图像的边界和半径为R的(更大)圆上的源轨迹。采样点由小圆圈表示。具有图23(A)所示的中央扇形(Fan_a)和采样图的数字中间图像要重新采样、并旋转到图23(D)所示中央扇形(Fan_d)的角度，具有图23(d)所示的最终采样图。这用以下两个步骤来完成：

(i)分别沿着Fan_a的每条射线将图像上采样成图23(b)中所示的采样图，它由Fan_a和Fan_d的交点定义。这些点因此将位于图23(c)所示的Fan_d上；

(ii)分别沿着Fan_d的每条射线将图像上采样成图23(d)中所示的采样图。

这些步骤使用1D上采样运算来实现重新采样和旋转的组合运算。

在对周期性采样图采样的中间图像情形中，为本发明各实施例描述的傅立叶分析技术被扩展成设计分散光束情形中的空间变化采样方案。这些技术足够普遍以应用于两维和三维中任意轨迹上的投射和反投射。对于非周期性直线或曲线系统上的反投射情形，局部光谱支持概念代替光谱支持概念。这在图24中示为由图24左侧示出的单个扇区的反投射所产生的中间图像。如以下将解释的，在所示点上以r和θ为参数的连续中间图像的局部光谱是图24(b)右侧所示傅立叶域中的线段。由多个扇形光束投射的反投射组成的中间图像的点上的局部光谱支持如图25所示。在左侧，示出了点的位置、以及用于中间图像的构成投射的视角范围。右侧的蝶形领结状区域是相应的局部光谱支持。

该局部光谱支持分析用来确定中间图像的局部采样需求。结果在空间上不均匀的采样图由图26A和26B中的虚线弧表示，用于两个典型的中间扇形反投射图像。图像边界由虚线圆表示，且仅有极少的点需要在该边界外取得。该局部傅立叶采样方法可直接应用于在直线、曲线或任意维的平面上寻找任意投射几何结构的采样方案。

确定用于重新采样的采样图、以及用于分层扇形光束反投射的坐标变换的基于傅立叶的方法可进一步地扩展到位于除中间图像的中央扇形之外的直线或曲线上的采样图，如下所述。该有用集的示例如图27(A)和27(B)所示。

普通发散光束算法

所述用于扇形光束反投射的方法直接扩展到其它发散光束几何结构，包括现代诊断放射医学中最重要内容之一的3D锥形光束，如下所述。类似于扇形光束几何结构，其中发散射线的源在围绕成像对象的轨迹(即圆)上移动，在普通发散光束的几何结构中发散射线的源在围绕成像对象的3D空间中的轨迹(并非必须为圆)上移动。检测器表面通过成像对象测量线积分。

用于普通发散光束几何结构的分层反投射算法的一个实施例可再次通过类似于图18的框图描述，但用以下方法更改。首先，在第0层，发散光束投射在99用与名义“零”源位置相对应的适当发散光束的单视圈反投射W₁来零反投射，从而产生初始中间图像。其次，因为源的轨迹并非必须是圆，中间图像的构成发散光束不仅可像在扇形光束几何结构一样旋转，而且还能彼此平移。相应地选择18中的坐标变换K。再次，取决于源轨迹和位置中的对称性的出现，可以有或没有“自由”的坐标变换，诸如替代扇形光束算法中π/2旋转的像素重新排列。像扇形光束情形一样，初始中间图像由该算法分层处理。与扇形光束情形相类似地，在位置和方位上彼此接近的中间图像被聚合，以便于聚合后的中间图像可稀疏地采样。通过研究旋转和相对彼此平移的构成加权反投射发散光束的结构，中间图像的3D采样图得以确定。如图28所示的一种采样图，是样本沿与中间图像的中央构成光束相对应的发散光束的射线分布的。沿每条射线的样本间隔被选择为确保该中间图像的所有构成发散光束都足以沿该射线采样。或者，更普通的方法是使用局部傅立叶方法来寻找采样图，如前所述的扇形光束情形。知道了中间图像如何由其构成投射组成，每个中间图像的局部3D傅立叶结构得以确定。中间图像的每一点上的3D局部采样矩阵函数被发现与局部傅立叶支持的采样需求相匹配，如扇形光束情形中所述。该矩阵函数然后在图像域上(可能数字地)积分，以确定样本的位置。

组合旋转和平移地上采样成新的采样光束的可分开方法类似于扇形光束情形实现。它将3D坐标变换和重新采样运算减少成1D重新采样运算序列。如图29(a)所示，每个发散光束被视为按照方位角分布的一组垂直坡面扇形光束、与在不同仰角上分布的倾斜坡面扇形光束集的交点。可分开坐标变换的步骤如下(如图29(b)中所示)：

1.原始发散光束被重新采样成原始射线与新发散光束的垂直平面的交点集。这些点因此位于最终采样点所共享的平面上。

2.来自扇形光束情形中可分开重新采样的步骤1和2对每个平面单独执行，以重新采样成最终点集。

使用适当的有效采样方案，用于发散光束的快速分层反投射算法可期望实现具有期望准确度的较大加速。如在扇形光束情形中可能使用在连续层中更改的伪光束(例如使用移到距离原点更远的顶点位置)。

对本领域技术人员显而易见的是，本发明诸方法并不限于在此所述的成像几何结构或特定实施例的示例。这些方法同样可应用于具有其它几何结构的反投射的一般问题。

图30示出基于重新采样的反投射，作为通过源位置β_p上的各个投射的一般化反投射所形成的上采样后中间图像的和。它与基于旋转的反投射4相似，但不同之处在于两个方面。第一，在第一步骤中(可能是处理后的)第p个投射在源位置或方位β_p进行加权反投射184，而不是4情形中的零反投射。例如在扇形光束投射的情形中，每个投射是在其源角度方位上反投射的。其次，在188通过累加聚合之前，每个初始中间图像在186进行上采样运算。这是必须的，因为这些P个单投射中间图像中的每一个的采样图被选为有效和稀疏的图，所以通常对每个投射是不同的，且常常是非笛卡尔的。在中间图像聚合之前，它们需要被重新采样成普通和更密集的栅格。该结构本身并不提供运算计数的减少。然而，它用作引入本发明的一些快速分层反投射方法的进阶。

图31示出本发明的另一个实施例——用于从多个投射q₁...q_p中创建像素图像的基于重新采样的分层反投射。图31是图30的二进制分层版本。首先，与非分层情形中一样，每个投射在其各个方位上反投射184成适于该方位的采样图，以产生中间数字图像。这是分层结构的第0层。在分层结构的第一层中，这些中间数字图像在186被上采样成图像选定对共用的一种更密集采样图。该上采样通常将是非笛卡尔采样图，但可通过一维重新采样运算序列执行，如图23或图29(b)所示。选定的结果上采样图像在190成对地聚合，以产生新的中间图像。在第二层中，新的中间图像在192再次上采样并在194聚合，以产生新的中间图像。该过程继续直到所有中间图像和投射都已得到处理，以在最后聚合步骤198之后产生最终图像

像前述实施例一样，运算可在层内或跨层组合。

分层结构中每个阶段上的采样图和图31所示实施例中的上采样运算的参数可通过任一前述方法选择。例如，在扇形光束投射的情形中，对给定中间图像的一种可能采样图将位于两个扇形的交点上：第一个扇形朝向中央构成源位置；第二个扇形朝向极值构成源位置。或者上采样步骤的采样图和参数可基于傅立叶或局部傅立叶属性、以及包括在中间图像中的投射视角来确定。特别地，对于非周期性的采样图，所述用于扇形光束的基于旋转算法的局部傅立叶方法可用来寻找采样图：知道投射如何组合以形成中间图像导致局部光谱支持的确定，该局部光谱支持用来计算在图像域上积分时产生中间图像的采样图的局部采样矩阵函数。

图32是快速分层重新投射的框图。重新投射是从给定电子图像计算X线断层投射的过程。重新投射算法通过将流程图换位的过程应用于反投射算法的任何框图来发现，其中在加权运算中可能有一些改变。在该过程中，运算被其共轭或二元运算所替代。因此重新投射的框图显现为与反投射的相应框图的反向版本类似。图32中所述的重新投射过程是从图8所示反投射算法中获取的一种这样的重新投射实施例。

在第一层中，输入图像f的副本保留在示图的顶部分支中作为顶部中间图像，并且在底部分支图像f在200旋转-π/2，以产生底部中间图像。在第二层中，在示图的上半部分，未旋转的中间图像在202进行三次独立的数字图像坐标变换，其中一部分不会使图像改变，从而产生三个不同的顶部中间图像。类似的过程应用于底部分支，从而产生三个底部中间图像。顶部和底部中间图像中的每一个(一共6个)然后在204进行抽选过程(低通滤波然后是下采样)，以产生新的中间图像。在图32所示实施例的实例中，仅有2·3^L(其中L＝2)个视角，所以第三层是递归分层结构中的最后一层。在第三和最后一层上，中间图像进行单独的坐标变换(206)，其中一部分不会使图像改变，从而产生18个中间图像。不是递归一部分的最后一个步骤是不同的：每个中间图像在208进行零度的重新投射。零度重新投射等同于累加垂直像素行以产生一维信号。这些一维信号(518)是该算法产生的输出投射。

算法中数字图像坐标变换的参数根据中间图像的傅立叶特征知识来选择。这些参数仅相关于相应反投射算法的参数。易于看到，因为重新投射框图是反投射框图的流程图换位，重新投射框图的每个分支在反投射框图中都有相应的分支，且相应分支中的坐标变换是彼此的数学共轭。在与图10双切变反投射算法相对应的该重新投射算法的版本中，第二层的坐标变换(202)是x切变加y切变，而最后一层的坐标变换(206)是x切变加y切变，以及x中的部分抽选(这三个运算可减为切变缩放)。这些切变参数是用于反投射算法的相应参数的负。x中抽选的参数与反投射算法的第一层中x上采样的相同。在本算法的切变缩放版本中(对应于图13所示的切变缩放反投射算法)，第二层的坐标变换(202)是x中的切变(且每个切变的参数是图13中相应参数的负)。最后一层中的坐标变换(206)是切变缩放。在重新投射算法中使用的切变缩放是在反投射中使用的相应切变缩放的数学共轭。抽选因子的参数也与反投射中相应的上采样因子相同。

与反投射算法相似，算法中的中间图像的过采样可通过在算法开始时首先上采样图像并在算法结束时下采样相同因数来实施。此外，这些运算可在层内或跨层组合。

图33是基于抽选的加权重新投射的框图。它是图31的流程图换位。它显示将图像重新投射成18个不同源角度上的投射集。开始时，给定图像沿两条平行路径进行处理。在每条路径中图像进行三次平行重新采样(210)，成为三个不同的更稀疏采样图。使用一维抽选(低通滤波、随后是重新采样成更稀疏栅格)，该重新采样能以可分开方法执行。滤波器的参数根据中间图像的傅立叶特征和期望投射来确定。期望投射的局部傅立叶分析用于投射不排列于平行直线上的情形中。在算法的最终层中，每个中间图像再次进行三次平行重新采样(212)，成为更稀疏的采样图。最后对图像执行加权投射以产生投射p_θ。这涉及要产生一维投射的图像像素的加权和。

实现和实验结果

本发明的较佳实施例用C编程语言实现，并用数字实验在具有384MB RAM的Sun Ultra 5工作站上测试。该测试图像是Shepp-Logan头部模型(head phantom)(用于X线断层造影算法的数字估算中的标准测试图像)。通过改变该算法的参数，可在准确度和计算成本之间找到平衡点。准确度指重建图像的质量。尽管视觉质量并不易于计量，所以我们测量重建图像和原始图像之间的误差，从中可用数字计算Radon变换。所使用误差的测量是相对均方根误差(RRMSE)。从f[m₂，m₁]的X线断层投射重建N×N图像 的RRMSE如下进行计算：

$\mathrm{RRMSE}=\frac{\sqrt{\frac{1}{N^2}{\mathrm{\Σ}}_{m_2=1}^N{\mathrm{\Σ}}_{m_1=1}^N{\left|\hat{f}[m_1,m_2]-f[m_1,m_2]\right|}^2}}{{\max }_{m_1,m_2}f[m_1,m_2]}---\left(1\right)$

对于平行光束数据，测试图像大小为256×256，视角数量为486，且Shepp-Logan滤波器(理想V型滤波器乘以正弦函数)用来滤波各个投射。图30显示在过采样参数γ为0.75～1.0之间的各个值上两种算法的RRMSE误差对运行时间。双切变算法由圆圈表示，而切变缩放算法由方块表示。每种算法使用两种类型的滤波器——称为MOMS 16的三阶(虚线)和五阶(实线)样条滤波器——来运行。对于该算法的每种特性，当γ减小时算法的误差降低而运行时间增加。未与任何其它点相连的曲线点是由相连点表示的算法的非过采样版本。作为比较，使用线性内插的常规算法的运行时间为14秒，且其输出图像的RRMSE误差为0.0486(比在此显示的快速算法差)。

来自用于平行光束数据的算法输出的一些样本图像，在图35中显示。从左到右按列地，它们是来自常规反投射、双切变、和切变缩放算法的输出图像。γ＝0.82的过采样应用于两种快速算法。下面一行图像详细地显示了上面一行相应图像的一部分。

用于扇形光束几何结构的快速分层算法在512×512  2D Shepp-Logan模型上成功测试。所考虑的获得几何结构具有源半径D＝1.06×N＝544，972个源角度，以及1025个等角检测器。快速算法的常规和过采样版本使用各种内插方法来实现。结果的重建图像与使用线性内插的常规扇形光束算法作比较。在所有实验中，投射都用Shepp-Logan滤波器滤波。

图36(A)和36(B)显示重建图像，图36(A)来自常规算法，而图36(B)来自快速算法。快速算法的结果可与常规算法的作比较。快速算法的点分布函数可与常规的作比较。图37(B)显示常规算法的PSF，图37(C)显示使用线性内插和γ＝0.4过采样的快速算法的PSF。图37(c)通过常规算法、非过采样的和γ＝0.4的过采样快速算法的psf的x轴比较各样条。PSF的相似性确认快速算法中相当的图像质量。

所使用的采样方案是基于扇形交点的方法，无需下文中提议的增强。除了先前提及的本方案缺点之外，为了简单实现，算法准确运算所不需要的许多采样点用于实验中测试的实施例中。尽管这些低效性，对于N＝512和D＝1.41N，常规算法中(数据相关)乘法运算与快速(线性)算法的比例为6.4，且加法运算的比例为3.0。注意，在该实验中使用的具有极接近原点的源(D＝2.06R)的几何结构，因为接近源顶点的样本的高密度所以对于该未增强实现特别具有挑战性。在更实用的系统中源更远一些，减少了这些效应。此外，使用前述选择性采样方案将导致高得多的加速因数。

用于本发明的反投射和重新投射算法的详细描述

反投射是用来从已处理投射中创建图像的过程。重新投射是逆过程，用来从给定图像计算投射。两种运算都用于根据投射的图像重建，如图1所示。首先将描述常规反投射，随后描述本发明的反投射和重新投射方法。将首先描述对2D图像的平行光束投射，然后描述更一般的2D和3D几何结构。

用来从其投射估算图像的典型和较佳方法是滤波反投射(FBP)算法。FBP首先涉及用所谓V型或更改后V型滤波器来滤波投射，然后反投射这些过滤后的投射。附加预处理步骤也可在反投射之前应用于该投射。在视角θ_p(p＝1，...，P)上取得的处理后投射是投射索引，取决于对变量t或θ的相关性是否需要显式地示出，将被交替地表示为q(t，p)或q_p或q_θr。为简单起见，处理后的投射在下文中将称为投射。

从一组P投射(在由长度为P的矢量

的分量指定的角度上)的FBP重建

被示为：

$\hat{f}(x_1,x_2)=\left(B_{\stackrel{\→}{\mathrm{\θ}}}q\right)(x_1,x_2)---\left(2\right)$

其中

定义0.1反投射运算符由以下定义：

$\left(B_{\stackrel{\→}{\mathrm{\θ}}}q\right)(x_1,x_2)\stackrel{\mathrm{\Δ}}{=}\underset{p=1}{\overset{P}{\mathrm{\Σ}}}q(x_1\cos {\mathrm{\θ}}_p+x_2{\sin \mathrm{\θ}}_p,p)---\left(3\right)$

本发明的基于O(N²logP)坐标变换的分层反投射算法基于反投射按照坐标变换的分层分解。若干较佳实施例的细节用坐标变换的不同选择来描述，以最一般的坐标变换结束。

特别地，本发明的基于O(N²logP)旋转的分层反投射算法基于反投射按照图像旋转的分解。

将图像

映射成其旋转版本(K(θ)f)

的旋转θ运算符由以下定义：

定义0.2

$\left(\mathrm{\κ}\left(\mathrm{\θ}\right)f\right)\left(\stackrel{\→}{x}\right)\stackrel{\mathrm{\Δ}}{=}f\left(K_{\mathrm{\θ}}\stackrel{\→}{x}\right)$

其中

定义0.3在平面上旋转角度θ的矩阵为 $K_{\mathrm{\θ}}=\left[\begin{array}{cc}\cos \mathrm{\θ}& -\sin \mathrm{\θ}\\ \sin \mathrm{\θ}& \cos \mathrm{\θ}\end{array}\right]$ 1D函数q(t)在单个角度θ上的反投射被定义为 $\left(B_{\mathrm{\θ}}q\right)\left(\stackrel{\→}{x}\right)\stackrel{\mathrm{\Δ}}{=}\mathrm{\πq}(x_1\cos \mathrm{\θ}+x_2\sin \mathrm{\θ}).$ 特别地，零角度反投射(对于θ＝0)为

$\left(B_{0}q\right)\left(\stackrel{\→}{x}\right)=q(x_1\cos 0+x_2\sin 0)=q\left(x_1\right)=q\left(\left[\begin{array}{cc}1& 0\end{array}\right]\stackrel{\→}{x}\right)---\left(4\right)$

定义0.4单个θ_p反投射的中间图像是 $I_{{\mathrm{\θ}}_p}\left(\stackrel{\→}{x}\right)\stackrel{\mathrm{\Δ}}{=}\left(B_{{\mathrm{\θ}}_P}{q}_p\right)\left(\stackrel{\→}{x}\right)$ 如图4可见，方程式(3)中的反投射可重写如下。在此q_p指其视角为θ_p的投射。

$\hat{f}\left(\stackrel{\→}{x}\right)=\underset{p=1}{\overset{P}{\mathrm{\Σ}}}\left(\mathrm{\κ}(-{\mathrm{\θ}}_p)B_0{q}_p\right)\left(\stackrel{\→}{x}\right)---\left(5\right)$

分层反投射的一个实施例基于若干连续旋转的累积结果仍然是旋转的事实。特别地，

         κ(θ₁)κ(θ₂)...κ(θ_N)＝κ(θ₁+θ₂+...θ_N)                (6)

从方程式(6)可得出，对于某整数L的p＝2^L，图4的框图可被重新安排成图5所示的分层树结构。第1层的第m个分支中的中间图像被示为I_l，m。在初始层中，

$I_{1,m}\stackrel{\mathrm{\Δ}}{=}{B}_0{q}_{m}m=\mathrm{1,2},...,P$ (7)

中间图像通过反投射个别投射来产生，并且在后续层中，即，对于l＝1，2，3，...，log₂P

$I_{l+1,m}\stackrel{\mathrm{\Δ}}{=}\mathrm{\κ}\left({\mathrm{\δ}}_{l,2m-1}\right)I_{l,2m-1}+\mathrm{\κ}\left({\mathrm{\δ}}_{l,2m}\right)I_{l,2m}$ (8)

中间图像进行坐标变换(在该实施例中为旋转)并通过累加聚合。对于任意组的投射角θ_i，i＝1，...，P，中间旋转角δ_l，m。可选为确保图4和5中的结构等效，从而图5所示的分层反投射算法产生期望图像

实际上，因为有比视角(约束)多得多的中间旋转角(自由参数)，所以选择δ_l，m时有许多自由度。此外，对实践中使用的数字图像使用数字图像坐标变换，其准确度和计算成本取决于采样图和实现的选择。这些各种自由度或参数在本发明中用于降低分层反投射算法的计算需求，如以下所述。

定义0.5中间图像I_l，m的构成投射的指数集N_l，m是有从图5中输入到中间图像I_l，m的路径的投射的指数集。例如，如图5中可见，N_3，P/4＝{P-3，P-2，P-1，P}。使用方程式(8)或在检查框图之后，便于显示

             N_l，m＝{2^l-1(m-1)+k∶k＝1，2，...，2^l-1}               (9)因而N_l，m列示组成中间图像I_l，m的投射。如果与由图5所示的方法处理相反，这些由集N_l，m指示的相同投射在其相应视角上一起直接反投射，则这将产生虚拟中间图像

定义0.6 ${\tilde{I}}_{l,m}={\mathrm{\Σ}}_{p{\∈N}_{l,m}}{I}_{{\mathrm{\θ}}_p}$

在检查图5中框图之后，可看到中间图像中的投射之间的相对角度在最终的重建图像 中得到保留。换言之，对于算法中的所有中间图像，即l，m，

$I_{l,m}=\mathrm{\κ}\left({\mathrm{\α}}_{l,m}\right){\stackrel{\‾}{I}}_{l,m}$ 对于某α_l，m∈(-π ，π](10)

分层算法的中间旋角δ_l，m完全根据α_l，m如下确定。从N_l，m的定义或方程式(9)中可简便地得出：N_l+1，m＝N_l，2m-1UlN_l，2m，且从而通过定义 I得出：

${\tilde{I}}_{l+1,m}={\stackrel{\‾}{I}}_{l,2m-1}+{\tilde{I}}_{l,2m}---\left(11\right)$

方程式(8)、(10)和(11)隐含

δ_l，2m-1＝α_l+1，m-α_l，2m-1和δ_l，2m＝α_l+1，m-α_l，2m      (12)

中间旋转角(δ_l，m)和因此相对角α_l，m被选择成用以减少中间图像的带宽。较小带宽的图像可由少数样本表示，并因此降低整个算法的计算成本。通过理解傅立叶域中的算法、并跟踪中间图像的光谱支持通过分层结构的进展，这些中间图像的带宽得以确定。

傅立叶域中的X线断层造影投射

数字图像坐标变换的参数被选为说明选定投射的视角θ_p，如前所述。参数还为它们对中间图像的傅立叶特征的效应进行选择。这些考虑也被用来确定选择将哪些中间图像聚合在一起。

解释傅立叶域中的反投射算法的关键是使投射P₀f的一维傅立叶变换(F₁)与图像f的两维傅立叶变换(F₂)相关的投射-片(projection-slice)定理。投射-片定理4称

(F₁P₀f)(ω)＝(F₂f)(ωcosθ，ωsinθ)

图4的反投射算法可在傅立叶域内解释。零反投射的运算符产生其光谱支持受限于水平频率轴ω₁的图像：(B₀q)(x₁，x₂)＝q(x₁)(F₂B₀q)(ω₁，ω₂)＝2πQ(ω₁)δ(ω₂)    (13)

众所周知，在空间域中旋转一图像导致其傅立叶变换也旋转相同角度。

$g\left(\stackrel{\→}{x}\right)=\left(B_{0}q\right)\left(K_{\mathrm{\θ}}\stackrel{\→}{x}\right)\⇔G\left(\stackrel{\→}{\mathrm{\ω}}\right)=\left(F_2{B}_{0}q\right)\left(K_{\mathrm{\θ}}\stackrel{\→}{\mathrm{\ω}}\right)---\left(14\right)$

可得出：方程式(5)中的反投射运算符的功能是将每个投射的光谱分量旋转到重建图像的适当角度，并将它们累加在一起。

假设具有单侧带宽(t变量中)等于π的投射，分层算法中典型的虚拟连续中间图像 I_l，m因此在由构成投射确定的角度范围内具有楔形的光谱支持，如图6所示。特别地，设I_l，m为图5(A)中所示算法的框图中的中间图像，其中N_l，m＝{b，b+l，...，e}，即该图像的构成投射的视角为{θ_p∶p＝b，b+l，...，e}。因为方程式(10)中的关系，I_l，m的光谱支持仅是 I_l，m的具有旋转角α_l，m的旋转版本。可得出，图6中楔形占据的角度的范围为[γ-β，γ+β]＝[θ_b-α_l，m，θ_e-α_l，m]。分别在第一和第二坐标系中I_l，m的带宽Ω_s和Ω_f在楔形的外围矩形上示出，如图6中的虚线所示。因而，α_l，m的选择提供了I_l，m带宽的控制方法。

采样理论指示这种连续中间图像I_l，m如何由笛卡尔采样图上的样本表示。特别地，第一和第二坐标系中具有带宽Ω_s和Ω_f的连续图像可在与坐标轴对齐的笛卡尔图上采样，其中第一和第二坐标系中采样间隔为Δ_s和Δ_f，从而Δ_s＜π/Ω_s和Δ_f＜π/Ω_f。为了使表示连续图像所需的样本数量最小化，应当最小化与外围矩形的面积成比例的数量1/(Δ_sΔ_f)＞π²Ω_sΩ_f。使该区域最小化、并因此使对中间图像I_l，m的采样需求最小化的旋转角为：

${\mathrm{\α}}_{l,m}^{*}=({\mathrm{\θ}}_b+{\mathrm{\θ}}_e)/2---\left(15\right)$

选择用于坐标变换的参数

方程式(15)用于方程式(12)中，以在图5所示的基于旋转的分层反投射实施例中确定用于数字图像坐标旋转的最优旋转参数δ_l，m。指数b和e被确定为分别对应于集N_l，m中的最小值和最大值。有了这种选择，中间图像 $I_{l,m}=K\left({\mathrm{\α}}_{l,m}^{*}\right){\stackrel{\‾}{I}}_{l,m}$ 及其在第一坐标系中的带宽(慢带宽)和第二坐标系中的带宽(快带宽)分别为：

$\begin{array}{cc}{\mathrm{\Ω}}_s\left(I_{l,m}\right)=\mathrm{\π}\sin \frac{{\mathrm{\θ}}_e-{\mathrm{\θ}}_b}{2}& {\mathrm{\Ω}}_f\left(I_{l,m}\right)=\mathrm{\π}\end{array}---\left(16\right)$

为了使带宽进一步最小化，应当最小化差值θ_e-θ_b，这通过选择在每个步骤上聚合从其视角之间具有最小的最大角距的投射产生的中间图像来实现。

除了旋转角和对要聚合哪些中间图像的选择以外，有必要选择算法中用于各个中间图像的数字图像坐标变换的适当采样图。这些使用采样要求进行选择，而采样要求又使用中间图像的光谱支持和带宽来确定。

特别地，对于通过单个投射的反投射形成的初始中间图像，如图5中的I_l，m＝B₀qθ_m，慢带宽Ω_s(I_l，m)＝0，且沿x₂坐标的单个样本充足。慢带宽将更大，且沿x₂坐标的更多样本将为如图3-b通过反投射一个以上投射形成的初始中间图像所需。

两个虚拟中间图像之和 $({\tilde{I}}_{l,m}={\stackrel{\‾}{I}}_{l-\mathrm{1,2}m-1}+{\tilde{I}}_{l-\mathrm{1,2}m})$ 及其相应中间图像如图7(A)-7(C)在空间域和频域中示出。图7(A)和7(B)示出 I和I的傅立叶域支持。图7(C)示出I_l，m、I_l-1.2m、I_l-1，2m-1的空间域采样方案。各采样点在空间域中的水平线和垂直线的交点上。聚合后的中间图像I_l，m具有比其每个分量更大的慢带宽，需要更密集的采样图。相反，在分层算法较早层上的中间图像需要更稀疏的采样图，从而节约计算。

计算成本和节约

计算成本是易于估算的。假设在[0，π)中等同地有P个投射视角，具有Δ_θ＝π/P的角度间隔以及(9)中指定的N_l，m，方程式(16)(使用b＝2^l-1(m-1)+1和e＝2^l-1(m-1)+2^l-1)暗示Ω_s(I_l，m)＝πsin (π(2^l-1-1)/(2P))和Ω_f(I_l，m)＝π。在第1层数字图像I_l，m ^d的大小(点的数量)因此是

$\mathrm{size}\left(I_{l,m}^d\right)\≈(N/{\mathrm{\Δ}}_s)(N/{\mathrm{\Δ}}_f)=({\mathrm{N\Ω}}_s\left(I_{l,m}\right)/\mathrm{\π})({\mathrm{N\Ω}}_f\left(I_{l,m}\right)/\mathrm{\π})$

$=N^2\sin (\mathrm{\π}(2^{l-1}-1)/\left(2P\right))$

$<N^2({\mathrm{\&pi;}2}^{l-1}/\left(2P\right))=O(N^2{2}^l/P).$

假设旋转S个像素的数字图像的成本是O(S)且 $\mathrm{size}\left(I_{l,m}^d\right)=O(N^2{2}^l/P),$ 该算法的复杂度是

$\underset{l=1}{\overset{\log P}{\mathrm{\&Sigma;}}}\frac{P}{2^{l-1}}O\left(\frac{N^2{2}^l}{P}\right)=\underset{l=1}{\overset{\log P}{\mathrm{\&Sigma;}}}O\left(N^2\right)=O\left(N^2\log P\right)---\left(17\right)$

对于典型的P＝O(N)，这是O(N²logN)，对该图像大小成本要比常规反投射的O(N³)在比例上有利得多。

经改进的三重基于旋转的分层反投射

容易看出：采样图像旋转某特定角度——0，±π/2和π——是计算上廉价的运算，因为它们不涉及内插而最多只是现有像素的重新排列。本发明算法的计算效率可通过将这些免运算(free operation)结合其中来改进。二进制算法(图5)可更改成在每个阶段累加三个而非两个图像。每个这种三元组的中央图像旋转0弧度-一免运算。为了使用-π/2的自由旋转，它被包括在最后阶段中：构成图像之一旋转-π/2并累加到另一未旋转图像上。双重和三重阶段的这个特定组合导致一组大小为2×3^L的视角(L为某整数)。尽管算法可适用于任意组投射和任意个投射，该算法的描述和分析可通过假设其数量正好为2×3^L且均匀分布而简化如下：

${\mathrm{\&theta;}}_l=-\frac{\mathrm{\&pi;}}{4}(1-3^{-L})+{\mathrm{\&Delta;}}_0(i-1)$ 其中i＝1，2，...，2·3^L且Δ_θ＝π/(2·3^L)                    (18)

因此，视角可被分成两组，每组3^L个，一组以θ＝0为中心而另一组以θ＝π/2为中心。该三重算法的框图以L＝2的特定情形在图8中示出，即一组P＝2×3²＝18个投射的特定情形。数字中间图像I_l，m ^d是基础的连续中间图像I_l，m的采样版本。在该框图中，标为B₀的框表示零角度反投射运算符。标为K_δ的框表示数字图像旋转运算符。如下所述，用于数字图像和用于数字图像旋转的采样周期和旋转角度被选为使计算成本最低。

选定的聚合中间数字图像可伸展以产生新的中间图像。图8中标有↑_y U_l，m的框表示沿第二个(慢)坐标系上采样，这对算法进展时容纳中间图像的渐增带宽是有必要的。对于(18)和(22)中视角的配置，第1层上中间图像的带宽为Ω_s(I_l，m)＝πsin(Δ_θ(3^l-1-1)/2)和Ω_f(I_l，m)＝π。在算法标为(1)，...(9)的各关键点上的标准化光谱支持如图9所示。随着算法的进展(当1增加时)，慢带宽增大，而所需的慢采样周期Δ_s＜π/Ω_s(I_l，m)＝sin((3^l-1-1)/2)^-1减小。各个上采样块中的上采样因数U_l，m可选择为满足这些要求。

小节0.0.0.1描述涉及两个切变的各个离散图像的可分旋转如何结合到反投射算法中。为了改进反投射后图像的质量，中间图像的系统过采样被引入。小节0.0.0.1描述了经更改包括了该过采样的算法。

用于坐标变换的中间旋转角可如下选择。对于一般L，分层结构包括(L+1)层。在第零层，每个投射q₀被零反投射(B₀)以产生相应图像B₀q_θ。然后各图像被分成三组，并组合成1/3多的图像。在层1上的分组由以下关系式定义。在后续层中，l＝1，...，L。

 I_l+1，m＝κ(δ_l，3m 2)I_l，3m-2+κ(δ_l，3m-1)I_l，3m-1+κ(δ_l，3m)I_l，3m    (19)

和               I_L+2，1＝κ(δ_L+1，1)I_L+1，1+κ(δ_L+1，2)I_L+1，2         (20)

对于(18)中的视角配置，最优中间旋转角是δ_l，3m..1＝0，l＝1，2，...，L、δ_L+1，1＝-π/2和δ_L+1，2＝-π/2。在检查图10之后，容易看出构成图像I_l，m的投射的指数N_l，m如下：N_l，m＝{3^l-1(m-1)+1，3^l-1(m-1)+2，...，3^(l-1)m}l＝1，2，...，L+1       (21)

通过方程式(21)，b＝3^l-1(m-1)+1和e＝3^l-1(m-1)+3^l-1，且使用方程式(18)和(15)， ${\mathrm{\&alpha;}}_{l,m}^{*}=-\frac{\mathrm{\&pi;}}{4}(1-3^{-L})+{\mathrm{\&Delta;}}_{\mathrm{\&theta;}}(3^{l-1}(m-1/2)-1/2)$ 。与二元的相似，可得出对于l＝1，...，L，

对于最后阶段我们使用自由旋转π/2和0弧度。

0.0.0.1坐标变换可基于切变，如下所述。中间图像的数字图像旋转可由图10(B)所示的两个数字图像切变序列代替，其中x和y坐标中的切变定义如下：

定义0.7 x切变和y切变由以下定义

$\left(S_x\left(\mathrm{\&alpha;}\right)f\right)\left(\stackrel{\&RightArrow;}{x}\right)=f\left(S_x^{\mathrm{\&alpha;}}\stackrel{\&RightArrow;}{x}\right)$

$\left(S_y\left(\mathrm{\&alpha;}\right)f\right)\left(\stackrel{\&RightArrow;}{x}\right)=f\left(S_y^{\mathrm{\&alpha;}}\stackrel{\&RightArrow;}{x}\right)$

其中x切变和y切变矩阵分别为 $S_x^{\mathrm{\&alpha;}}=\left[\begin{array}{cc}I& \mathrm{\&alpha;}\\ 0& 1\end{array}\right]$ 和 $S_y^{\mathrm{\&alpha;}}=\left[\begin{array}{cc}1& 0\\ \mathrm{\&alpha;}& 1\end{array}\right].$

该选择性实施例从旋转矩阵的双切变因数分解导出： $K_{\mathrm{\&theta;}}=S_y^{\tan \mathrm{\&theta;}}{S}_x^{-\sin \mathrm{\&theta;}\cos \mathrm{\&theta;}}{S}_c,$ 其中 $S_c=\left[\begin{array}{cc}\cos \mathrm{\&theta;}& 0\\ 0& 1/\cos \mathrm{\&theta;}\end{array}\right].$ 该分解暗指：在数字图像的情形中，使用完美的限带宽内插和带宽充分受限的数字图像作为双切变数字图像坐标变换的输入，该双切变图像仅是图像旋转θ，且有效地在x下采样并在y上采样1/cosθ因数。为了校正采样图中的这些变化，在每次执行双切变运算时都需要非整数的重新采样。相反，通过不同地采样中间图像，仅需要对最终图像解决累计效果。

方便的选择是：通过在x上采样初始数字中间图像并在y上下采样，在一开始分层算法时就校正重新采样。当如图10所示实施例中通过单视圈反投射产生初始中间图像时，y上下采样是自由的，因为它仅涉及初始反投射到具有不同y密度的采样图。因此未在图10中框图明显示出。第一阶段中在x上采样避免了在后面阶段中的伪信号(aliasing)问题。对于初始图像中的每一个，对最终图像的累计有效下采样/上采样因素1/П^K _k＝1cosθ_k在对到最终图像的路径的所有双切变变换上计算，而逆运算应用于所述初始图像。可得出：初始图像I_l，p ^d，p＝1，...，P中的每一个用不同的x和y采样密度进行重新采样。

因此10(A)所示的结果算法包括：首先重新采样所有图像适当因数B₀q_θp，p＝1，...，P(取决于θ_p)，然后执行双切变数字图像坐标变换，如图10(B)所示。该双切变考虑每个中间图像的变化采样图来执行。

切变缩放算法

基于切变缩放的分层反投射算法以反投射的定义为基础，该反投射是被缩放(拉伸/收缩)和切变的单投射图像之和。

定义0.8 x切变缩放运算符 由以下定义

$\left(C\left(\stackrel{\&RightArrow;}{\mathrm{\&sigma;}}\right)f\right)\left(\stackrel{\&RightArrow;}{x}\right)=f\left(C_{\stackrel{\&RightArrow;}{\mathrm{\&sigma;}}}\stackrel{\&RightArrow;}{x}\right)$

其中x切变缩放矩阵 $C_{\stackrel{\&RightArrow;}{\mathrm{\&sigma;}}}=\left[\begin{array}{cc}{\mathrm{\&sigma;}}_1& {\mathrm{\&sigma;}}_2\\ 0& 1\end{array}\right].$

所以，方程式(3)可重写为 $\hat{f}\left(\stackrel{\&RightArrow;}{x}\right)=\underset{p=1}{\overset{P}{\mathrm{\&Sigma;}}}\left[C({\cos \mathrm{\&theta;}}_p,{\sin \mathrm{\&theta;}}_p)B_0{q}_p\right]\left(\stackrel{\&RightArrow;}{x}\right)---\left(23\right)$

以下是描述可便于显示的累积切变和切变缩放的效果的一些结果：

$\underset{i=1}{\overset{n}{\mathrm{\&Pi;}}}{S}_x\left({\mathrm{\&alpha;}}_i\right)=S_x\left(\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i\right)---\left(24\right)$

$C\left({\stackrel{\&RightArrow;}{\mathrm{\&sigma;}}}^{\&prime;}\right)C\left({\stackrel{\&RightArrow;}{\mathrm{\&sigma;}}}^{\&prime;\&prime;}\right)=C({\mathrm{\&sigma;}}_1^{\&prime;}{\mathrm{\&sigma;}}_1^{\&prime;\&prime;},{\mathrm{\&sigma;}}_1^{\&prime;\&prime;}+{\mathrm{\&sigma;}}_1^{\&prime;}{\mathrm{\&sigma;}}_2^{\&prime;\&prime;})---\left(25\right)$

与方程式(24)相似，方程式(25)可归纳地应用，以证实一个n切变缩放序列也是切变缩放。使用该属性，方程式(23)可用分层结构算法结构来重写。图11(A)是方程式(23)的框图表示(P＝4)，而图11(B)是其分层等效体。

两种结构的等同性要求建立给定投射角θ_p和未知切变缩放参数 σ_l，m之间相等的系统，其解总是存在。

使用整数缩放因数而不是非整数缩放因数在计算上是有利的(根据运算或存储空间)。例如，使用x切变缩放 (其中整数切变参数σ₁＝1.0)是有利的，它实际上是纯粹的x切变S_x(σ2)。将算法设计成在层l＝2，3，...，L上仅使用纯粹的x切变并避免x切变缩放是可能的。

与用一系列切变缩放替换方程式(23)中的每个切变缩放相反，使用方程式(24)和(25)每个切变缩放都用单个切变缩放加一系列切变来替换：

$C({\cos \mathrm{\&theta;}}_p,{\sin \mathrm{\&theta;}}_p)=S_x\left({\mathrm{\&alpha;}}_n\right)S_x\left({\mathrm{\&alpha;}}_{n-1}\right)...S_x\left({\mathrm{\&alpha;}}_1\right)C({\cos \mathrm{\&theta;}}_p,\mathrm{\&sigma;})$

$\begin{array}{cc}s.t.& {\mathrm{\&sigma;}+\cos \mathrm{\&theta;}}_p\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&alpha;}}_i={\sin \mathrm{\&theta;}}_p\end{array}---\left(26\right)$

具有自由参数α₁，α₂，...，α_n。

结果分层配置在图13中显示。方便的选择是在层l＝1上执行x切变缩放 并在后续层l＝2，3，...，L上x切变。所以第一层中基础的连续域图像的表达式如下：

$I_{2,m}=\underset{i=0}{\overset{2}{\mathrm{\&Sigma;}}}C\left({\stackrel{\&RightArrow;}{\mathrm{\&sigma;}}}^{3m-i}\right)B_0{q}_{3m-i}---\left(27\right)$

对于l＝2，3，...，L，中间图像I_l+1，m由前一层中的三个中间图像({I_l，3m-1}_l＝0 ²)的x切变版本组成：

         I_l+1，m＝S_x(α_l，3m..2)I_l，3m-2+I_l，3m-1+S_x(α_l，3m)I_l，3m    (28)

并且为了使用0和-π/2弧度的自由旋转，在最后一层中：

                      I_L+2，1＝I_L+1，1+κ(-π/2)I_L+1，2

图13中所示算法中坐标变换的参数(用于各种切变缩放和切变)尚未确定。实际上，仅有第一层的(非整数)缩放因数σ₁根据相关联投射的角度被完全确定。中间切变因素({α_l，m}_{l＝2，3...L})被选为使中间图像的带宽最小，以便于可最稀疏地采样。给定任一组{α_l，m}，在层l＝1上σ₂(切变因数)总是可选为确保由图13的整个分层算法产生的反投射后图像 等于反投射表达式(23)。特别地，分层算法初始层中的切变缩放参数取决于组{α_l，m}，如下所示：

${\mathrm{\&sigma;}}_2^p=\left\{\begin{array}{cc}{\sin \mathrm{\&theta;}}_p-{\cos \mathrm{\&theta;}}_p{\mathrm{\&Sigma;}}_{l=2}^L{\mathrm{\&alpha;}}_{l,\mathrm{\&mu;}(p,l)}& :p\&le;\frac{P}{2}\\ \sin ({\mathrm{\&theta;}}_p-\frac{\mathrm{\&pi;}}{2})-\cos ({\mathrm{\&theta;}}_p-\frac{\mathrm{\&pi;}}{2}){\mathrm{\&Sigma;}}_{l=2}^L{\mathrm{\&alpha;}}_{l,\mathrm{\&mu;}(p,l)}& :p>\frac{P}{2}\end{array}\right.$

其中

描述构中从层1中第p个投到根的路径。

通过检查图13中框图的上半部分p≤P/2，然后从最后一层反向序列切变，可以看到在该分层结构的上半部分，基础连续中间图像I_l，m通过切变与其相关联的虚拟中间图像

相关，即

$I_{l,m}=S_x\left({\mathrm{\&beta;}}_{l,m}\right){\stackrel{\&OverBar;}{I}}_{l,m}---\left(29\right)$

与方程式(12)的推导相似，可示出中间切变因数取决于{β_l，m}，如下所示：

α_l，3m-i＝β_l+1，m-β_l，3m-i    i＝0，1，2.                     (30)

选择切变系数α_l，m的自由度提供选择β参数中的自由度，它可用于使中间投射的带宽最小化。这将降低其采样要求，并改进算法的计算效率。

设中间图像I_l，3m-1的光谱支持由W_l，3m-i表示，且其y带宽(即慢带宽)由Ω_y(W_l，3m-i)表示。然后使Ω_y(W_l，3m-i)最小化的最优β_l，3m-i对i＝0，1，2为

${\mathrm{\&beta;}}_{l,3m-i}^{*}=\underset{\mathrm{\&beta;}\&Element;R}{\arg \min }\left\{{\mathrm{\&Omega;}}_y\left(S_y\left(\mathrm{\&beta;}\right){\stackrel{\&OverBar;}{W}}_{l,3m-i}\right)\right\}---\left(31\right)$

如果各视角均匀分布，则可发现 ${\mathrm{\&beta;}}_{l,3m-1}^{*}\&ap;{\mathrm{\&beta;}}_{l+1,m}^{*}.$ 然后一种有利的选择是 ${\mathrm{\&beta;}}_{l,3m-1}^{*}\stackrel{\mathrm{\&Delta;}}{=}{\mathrm{\&beta;}}_{l+1,m}^{*},$ 结果 ${\mathrm{\&alpha;}}_{l,3m-1}^{*}=0,$ 使所需切变运算消除约1/3。

组{β_l，m}通过从分层结构中最后一层到第一层反向地跟踪算法中中间图像的光谱支持来确定。根据片投射定理，显然 I_l，m的光谱支持是

${\tilde{W}}_{l,m}=\{(\mathrm{\&gamma;}{\cos \mathrm{\&theta;}}_p,{\mathrm{\&gamma;}\sin \mathrm{\&theta;}}_p):p\&Element;N_{l,m}\mathrm{\&gamma;}\&Element;[-\mathrm{\&pi;},\mathrm{\&pi;}\left]\right\}---\left(32\right)$

为了采样要求，光谱支持 W_l，m由其标示为

的端点来表征：

${\tilde{E}}_{l,m}=\{({\cos \mathrm{\&theta;}}_p,{\sin \mathrm{\&theta;}}_p):p\&Element;N_{l,m}\}---\left(33\right)$

并如图12所示。该集 的上下带宽边界分别由

和 表示。所以如果N_l，m＝{b，b+1，...，e}，则 ${\tilde{E}}_{l,m}^{+}=({\cos \mathrm{\&theta;}}_e,{\sin \mathrm{\&theta;}}_e)$ 且 ${\tilde{E}}_{l,m}^{-}=({\cos \mathrm{\&theta;}}_b,{\sin \mathrm{\&theta;}}_b).$ 由于三重组合规则，这些光谱支持端点的上、中、下三组被标示为 ${\tilde{E}}_{l+1,m}^0,{\tilde{E}}_{l+1,m}^1,{\tilde{E}}_{l+1,m}^2,$ 即 ${\tilde{E}}_{l+1,m}^i={\tilde{E}}_{l,3m-i}.$ 并从(30)和有关仿射变换的傅立叶理论得出 $E_{l,m}=S_y\left({\mathrm{\&beta;}}_{l,m}\right){\tilde{E}}_{l,m}.$ 然后最优切变因数为

${\mathrm{\&alpha;}}_{l,3m-i}^{*}=-\frac{{\left(E_{l+1,m}^i\right)}_y^{+}+{\left(E_{l+1,m}^i\right)}_y^{-}}{{\left(E_{l+1,m}^i\right)}_x^{+}+{\left(E_{l+1,m}^i\right)}_x^{-}}i=\mathrm{0,2}---\left(34\right)$

图14中示出的算法FINDALPHAS寻找当在切变缩放反投射算法中从层L到层2向下遍历该分层结构时的最优切变因数α_l，m ^*。该算法的输入是在第(L+1)层开始时图像的光谱支持的端点集：E_L+1，m，m＝1，2。特别地，在均匀分布角度的情形中，E_L+1，1＝E_L+1，2＝{(cosθ_p，sinθ_p)：p＝1，2，...，P/2}。

寻找上采样因数

诸算法和本发明不同实施例的y上采样因数{U_l，m}确定慢方向中中间图像的采样密度。这些上采样因数可选为符合中间图像的采样要求，同时限制为整数值。限制为整数有计算上的优点，并出于同样原因可扩展成具有小值分母的有理数。

选择这些因数的问题在切变缩放算法中比双切变情形略简单些。在前一情形中，查看图13容易看出：中间图像I_l，m的慢采样间隔是 ${\mathrm{\&Pi;}}_{l^{\&prime;}=l}^L{U}_{l^{\&prime;},\mathrm{\&mu;}(l^{\&prime;},l,m)},$ 其中

述了从I_l，m到算法最后节点的路径。采样理论要求慢方向采样间隔如下：

${\mathrm{\&Delta;}}_s^{l,m}=\underset{l^{\&prime;}=l}{\overset{L}{\mathrm{\&Pi;}}}{U}_{l^{\&prime;},\mathrm{\&mu;}(l^{\&prime;},l,m)}\mathrm{\&pi;}/{\mathrm{\&Omega;}}_y\left(W_{l,m}\right)---\left(35\right)$

其中Ω_y(W_l，m)是I_l，m的y带宽。假设一组上采样因数为 $u\stackrel{\mathrm{\&Delta;}}{=}\{U_{l,m}\&Element;Z:l=\mathrm{2,3},...,L;$ m＝1，2，...，6.3^L-1}整个算法的计算成本可示为：

$j\left(u\right)=\underset{l=2}{\overset{L}{\mathrm{\&Sigma;}}}\underset{m=1}{\overset{{6.3}^{L-1}}{\mathrm{\&Sigma;}}}\left(\frac{c_{l,m}}{{\mathrm{\&Pi;}}_{l^{\&prime;}=l}^L{U}_{l^{\&prime;},\mathrm{\&mu;}(l^{\&prime;},l,m)}}\right)+c---\left(36\right)$

其中常数c_l，m表示应用于中间图像I_l，m ^d的数字坐标变换的相对计算成本。这些常数通过所使用数字滤波器的阶、以及坐标变换的特定实现来确定，如将要在坐标变换的有效实现中所述。满足约束方程式(35)的使J(U)最小化的最佳集U可通过动态编程或综合性搜索来解出。对于一给定视角组，上采样因数的最佳集可预先计算并存储在查询表中。

在双切变旋转算法的情形中，问题因为使用双切变旋转导致中间图像的实际分数上下采样的事实而略为复杂化。特别地，算法相邻层中的慢采样周期Δ_l，m相关如下： 其中

所以对慢方向采样间隔的限制变成如下：

${\mathrm{\&Delta;}}_s^{l,m}=\underset{l^{\&prime;}=l}{\overset{L}{\mathrm{\&Sigma;}}}\left(U_{l^{\&prime;},(l^{\&prime;},l,m)}{\mathrm{\&kappa;}}_{l^{\&prime;},(l^{\&prime;},l,m)}\right)\mathrm{\&pi;}/{\mathrm{\&Omega;}}_y\left(W_{l,m}\right)---\left(37\right)$

其中Ω_y(W_l，m)是I_l，m的y带宽。方程式(36)给出的计算成本现在可以方程式(37)中的约束为条件最小化。

算法的过采样版本

过采样被结合在本算法中，以改进重建的准确度。一种这样做的方法是在算法中的所有中间图像上统一：在算法中执行内插时，制作运算的图像至少过采样某预定常数γ。特别地，进行数字坐标变换或重新采样的每个中间图像的Nyquist(尼奎斯特)速率与采样频率(慢方向和快方向上)之比最好应小于γ。该过采样最好在尽量少地更改本发明算法时结合。

慢方向上的过采样

在先前所述的本发明算法中，慢方向上的采样频率由上采样因数U_l，m控制。这证实是在算法过采样版本中保持过采样条件的有用工具，因此导致对本发明算法的结构没有改变。以1/γ倍上采样通过简单地更改慢方向采样间隔上的约束γ倍来实现，即需要 ${\mathrm{\&Delta;}}_s^{l,m}\&le;\mathrm{\&gamma;\&pi;}/{\mathrm{\&Omega;}}_s\left(I_{l,m}\right),l=\mathrm{2,3},...,L.$

快方向上的过采样

在快方向上，计算上廉价的整数上采样并不用来控制过采样，所以算法被更改为包括分数上采样。在双切变分层反投射算法中，这种慢方向上的分数上采样因数{U_l，m∶m＝1，2，...，2.3^L}已包括在l＝1层中。因此我们简单地增大上采样因数U_l，m以在执行了上一数字坐标变换之后结合该过采样并下采样该图像，以将该图像返回给所需采样方案(其中Δ_f＝Δ_s＝1.0)。这种对双切变算法的更改因此仅涉及分数重新采样的一个附加层。该x坐标重新采样与第L层中的数字x切变组合，用于改进计算效率。该算法的框图如图15所示。在切变缩放分层反投射算法中，x上采样运算与x切变缩放(在算法开始时为 )组合，以避免额外层重新采样。容易看出，x上采样仅仅是x切变缩放的特定情形，其中切变因数设定为0。该两种运算的组合仍然是x切变缩放：

C(σ₁，σ₂)C(1/U，O)＝C(σ₁/U，σ₂/U)。

算法结束时的下采样与第L层中的x切变(S_αL，m|_{m＝1，...，6})相组合，以使它们成为有效的数字图像的切变缩放。可示出x切变α加上x下采样U是有效的x切变缩放：C(U，O)S_x(α)＝C(U，1+U_α)。

这些上采样和下采样因数的确切值根据参数γ和中间图像的光谱结构确定。在快方向中，过采样条件是快方向(y坐标)采样间隔满足 ${\mathrm{\&Delta;}}_f^{l,m}\&le;\mathrm{\&gamma;\&pi;}/{\mathrm{\&Omega;}}_f\left(I_{l,m}\right),l=\mathrm{2,3},...,L.$ 所以所必需的附加上采样为

$\mathrm{\&eta;}=\frac{{\max }_{2\&le;l\&le;L}{\mathrm{\&Omega;}}_f\left(I_{l,m}\right)}{\mathrm{\&gamma;\&pi;}}---\left(38\right)$

因此，过采样算法的上采样因数

根据非过采样的双切变算法的U_l，m更改如下：

${\tilde{U}}_{1,m}=U_{1,m},{\tilde{U}}_{L+1,m}=\mathrm{\&eta;}---\left(39\right)$

换言之，第一层中的上采样因数被更改成确保算法中的所有中间图像都根据参数γ来过采样。因此，在最后旋转之后，图像具有快速采样间隔1/η；因此，在第L层，它们必须被下采样η以向单元采样方案返回它们。

在此必须注意，因为输入投射被假设为确切地以Nyquist速率采样，快方向中的过采样条件将不被第一层图像B₀q_0p，p＝1，...，P满足。

该三重过采样的基于切变缩放的算法框图如图16所示。

相似运算的压缩序列

在同一单个坐标上有两个运算序列时，它们可相组合以改进计算成本和重新采样准确度。以下是除了前述的组合x切变与x上采样/下采样、或x切变缩放与x-上采样/下采样的另一个示例。

图10的非过采样的双切变算法包括第L阶段中6个中间图像的4个的y切变加上x切变。我们通过将两个运算序列——该阶段的x下采样加上x切变——压缩成一个，来把图像的最终下采样结合成图15的过采样版本。这使内插级联的长度与非过采样情形中未有改变。

任意基数/分支系数的分层

对于视角组不是形式2*3^L的情形，所有这些算法都可简便地更改。尽管较佳实施例用于包括聚合中间图像的三元组的分层算法的所有分支，或使用±π/2或0的旋转，但可在每个阶段上组合任意数量的中间图像。

假设在算法的特定层上有任意数量(假设为M)的中间图像，它们中的3×M/3可组合成三个组(其中x是小于或等于x的最大整数)。如果剩余数量的中间图像是两个，则它们可聚合成对以在下一层上产生图像。如果仅剩下单个中间图像，则它可不加更改地传递给分层算法的下一层。

分层算法的分支系数(聚合在分层算法节点上的分支数量)可变化，不仅是为了容纳任意数量的视角，而且是为了减少分层算法的深度，从而改进图像质量。在该情形中，不成对或不以三元组、而是以更大的组来聚合图像是有用的。

坐标变换的参数的先前规定可简便地扩展成具有任意分支系数的节点。例如，在基于旋转的算法中，其中 $I_{l+1,m}={\mathrm{\&Sigma;}}_{i=1}^M\mathrm{\&kappa;}\left({\mathrm{\&delta;}}_{l,m_i}\right)I_{l,m_i},$ 程式(10)和(15)中所述的关系仍然成立，因此旋转角度由δ_l，mi＝α_l+1，m-α_l，m，指定。在基于切变缩放算法的情形中，其中 $I_{l+1,m}={\mathrm{\&Sigma;}}_{i=1}^M{S}_x\left({\mathrm{\&alpha;}}_{l,m_i}\right)I_{l,m_i},$ 方程式(29)仍然成立。将概念E_l，m ^l扩展成指正在聚合的M组投射的第i组，可获得等于方程式(34)右侧的α_l，m ^*的方程式。

0.1基于其它图像变换的分层算法

对于任何矩阵A_θ，反投射方程式(3)可写成形式

$\hat{f}\left(\stackrel{\&RightArrow;}{x}\right)=\frac{1}{P}\underset{p=1}{\overset{P}{\mathrm{\&Sigma;}}}\left[B_0{q}_p\right]\left(A_{{\mathrm{\&theta;}}_P}\stackrel{\&RightArrow;}{x}\right)---\left(40\right)$

只要A_θ的第一行为[cosθsinθ]。对A_θ的剩余两个条目选择任意值的自由度允许在设计分层算法中使用的坐标变换时的灵活性。矩阵A_θ可被表示为与θ相关的某些参数向量δ_l的因式： $A_{\mathrm{\&theta;}}={\mathrm{\&Pi;}}_{l=1}^{L}A\left({\mathrm{\&delta;}}_l\right),$ 但由于A_θ两个底部条目的自由度而没有完全确定。该因式分解可用来得出反投射方程式(40)的分层分解，诸如图5所示的相应框图，其中坐标变换步骤由表示定义为 $\left(\mathrm{\&kappa;}\left({\mathrm{\&delta;}}_{l,m}\right)f\right)\left(\stackrel{\&RightArrow;}{x}\right)=f\left(A\left({\mathrm{\&delta;}}_{l,m}\right)\stackrel{\&RightArrow;}{x}\right)$ 的图像坐标变换运算符的κ(δ_l，m)标示。

显然，在此所述的本发明特定实施例是矩阵A_θ及其因式分解、及相关联坐标变换的较普遍选择的特定情形。这些坐标变换对中间图像的傅立叶光谱的效果与已经描述的情形相似地加以分析，因为由矩阵A在空间域中的仿射变换的效果是矩阵A^-T(A的逆转置矩阵)在频率域中的仿射变换。因此类似的考虑用来选择变换中的自由参数，其目标是降低计算要求。因而，用于本发明的分层反投射算法中的一类数字图像坐标变换除了所述特定实施例之外还包括许多其它变换。此外，因为矩阵A(δ_l，m)可被因式分解成三角矩阵的积，所以坐标变换可按需作为单轴坐标变换的级联来执行。

0.2数字图像坐标变换和重新采样的有效和准确实现

本发明的分层反投射和重新投射算法的准确度和速度取决于各个数字图像坐标变换和重新采样运算的实现细节。经改进的准确度一般需要通常增加计算成本的高阶滤波器或内插。高阶滤波或内插的成本可通过将运算分解成较低阶运算来降低。如果所使用的滤波器是移位不变的，则可获得计算和/或存储器需要的其它降低。高阶有限脉冲响应滤波器可通过使用低阶递归滤波器、或通过使用快速傅立叶变换(FFT)，来有效地实现。

特别地，如果对连续图像假设可分的表示基础，则数字图像的数字y切变可通过数字图像阵列中各个竖直行的部分延迟来获得。类似地，数字x切变可表达为某时一行的部分延迟。从而，1D信号的部分延迟可使用移位不变滤波器来完成。3D图像和切变运算的类似分解是众所周知的。

重新采样数字图像通常也可使用常常是移位不变的低维运算来执行。例如，沿一坐标上采样整数倍U的图像可被分解成U个不同的有效计算的部分延迟。这实际上是数字信号处理中众所周知的所谓多相分解。沿一坐标的有理但非整数重新采样可被分解成每个都有效执行的整数下采样和上采样的级联。

更普通的数字图像重新采样也可被分解成低维运算。考虑数字2D图像f从位于标示为CF₁的曲线族上样本S₁的采样图，重新采样成具有位于不同曲线族CF₂上的点S₂的另一个图，从而产生数字图像h。如果两个曲线族足够密度地相交，则参照图23所述的用于从一个扇形重新采样成另一旋转后扇形的方法可用于一般曲线。否则可引入以所需密度与CF₁和CF₂相交的第三曲线族CF₃。数字图像可从CF₁重新采样到其与CF₃的交点，然后到CF₃与CF₂的交点，最后到CF₂，成为所需采样图。

通过例如考虑面而非曲线，该过程归一化为3D，以首先减少对表面的重新采样之一的过程，然后使用2D过程进一步减少对表面的重新采样以对曲线重新采样。

数字坐标变换和重新采样常常可组合，从而改进计算效率。例如，用于图13中切变缩放算法的数字x切变缩放运算可被分解成对各条水平线的重新采样运算。

1发散光束的快速分层反投射算法

1.1扇形光束投射和反投射

考虑在围绕对象的圆上放置的等角分隔检测器的情形。两维图像f(x，y)的源角度为β的扇形光束X线断层造影投射由(R_βf)(γ)标示，并被定义为沿以原点为圆心半径为D的圆上的源位置为中心、参数为γ的扇形射线的线积分组。函数

被假设为在半径为R的盘外为零值。

如图17所示的扇形光束射线V(β，γ，T)定义如下：

$V(\mathrm{\&beta;},\mathrm{\&gamma;},T)=D\left[\begin{array}{c}-\sin \mathrm{\&beta;}\\ \cos \mathrm{\&beta;}\end{array}\right]+T\left[\begin{array}{c}\sin (\mathrm{\&beta;}+\mathrm{\&gamma;})\\ -\cos (\mathrm{\&beta;}+\mathrm{\&gamma;})\end{array}\right]---\left(41\right)$

在源角度为β扇形角度为γ的扇形光束投射为

$\left(R_{\mathrm{\&beta;}}f\right)\left(\mathrm{\&gamma;}\right)={\&Integral;}_{T=0}^{\&infin;}f\left(V(\mathrm{\&beta;},\mathrm{\&gamma;},T)\right)\mathrm{dT}.$ 因为f在半径为R的盘外为0，所以仅需在盘内即T_ST和T_END之间执行积分。

在使用扇形光束几何结构的计算X线断层造影中，投射在离散源角度集{β_p∶p＝1，2，...，P}上可用，且在每个扇形内射线的角度由{γ_J∶j＝1，2，...，J}索引。在等角扇形光束几何结构的情形中，检测器均匀分布在以该源为中心的圆弧上，所以扇形角度是均匀间隔的。在等间隔检测器的情形中，检测器均匀分布在与从源位置到原点的线垂直的线上。在此所述的用于扇形光束几何形状的快速反投射算法假设为等角分布，其中J为奇数，且γ₃＝Δ_γ·(j-(J+1)/2)。但是，诸算法可简便地扩展成其它扇形角几何结构。

根据一组P个扇形光束投射 $\{R_{{\mathrm{\&beta;}}_i}\left(\mathrm{\&gamma;}\right):i=\mathrm{1,2},...,P\}$ 的重建算法，可表达为每个扇形光束投射加上加权反投射(5)的缩放和滤波：

$\hat{f}\left(\stackrel{\&RightArrow;}{x}\right)=\underset{p=1}{\overset{P}{\mathrm{\&Sigma;}}}W\left(T(\stackrel{\&RightArrow;}{x},{\mathrm{\&beta;}}_p)\right)q_{{\mathrm{\&beta;}}_p}\left(\mathrm{\&gamma;}\left(\stackrel{\&RightArrow;}{x},{\mathrm{\&beta;}}_p\right)\right)---\left(42\right)$

其中W(T)是适当的加权函数，

和

是分别是沿射线的源和源位置β的图像点

之间的距离，以及该射线的扇形角度，而q_β(γ)是加权和滤波后的投射

                Q_β(γ)＝R′_β(γ)*g(γ)，R′_β(γ)＝R_β(γ)·D·cosγ其中g(γ)是适当的滤波器。

定义1.1单个源角度β上函数q的加权反投射定义为

$\left(W_{\mathrm{\&beta;}}q\right)\left(\stackrel{\&RightArrow;}{x}\right)\stackrel{\mathrm{\&Delta;}}{=}W\left(T(\stackrel{\&RightArrow;}{x},\mathrm{\&beta;})\right)q\left(\mathrm{\&gamma;}(\stackrel{\&RightArrow;}{x},\mathrm{\&beta;})\right)---\left(43\right)$

并且加权反投射运算符 $W_{\stackrel{\&RightArrow;}{\mathrm{\&beta;}}}(\stackrel{\&RightArrow;}{\mathrm{\&beta;}}\&Element;{\left[\mathrm{0,2}\mathrm{\&pi;}\right)}^p)$ 定义为

$\left(W_{\stackrel{\&RightArrow;}{\mathrm{\&beta;}}}{\left\{q_p\right\}}_{p=1}^P\right)\left(\stackrel{\&RightArrow;}{x}\right)\stackrel{\mathrm{\&Delta;}}{=}\underset{p=1}{\overset{P}{\mathrm{\&Sigma;}}}\left(W_{{\mathrm{\&beta;}}_p}{q}_p\right)\left(\stackrel{\&RightArrow;}{x}\right)---\left(44\right)$

因此 $\hat{f}\left(\stackrel{\&RightArrow;}{x}\right)=\left(W_{\stackrel{\&RightArrow;}{\mathrm{\&beta;}}}{\left\{q_p\right\}}_{p=1}^P\right)\left(\stackrel{\&RightArrow;}{x}\right).$

容易示出

                     W_β＝Κ(-β)W₀                (45)

其中Κ(β)表示β运算符的旋转，因此

$\left(W_{\stackrel{\&RightArrow;}{\mathrm{\&beta;}}}{\left\{q_p\right\}}_{p=1}^P\right)\left(\stackrel{\&RightArrow;}{x}\right)=\underset{p=1}{\overset{P}{\mathrm{\&Sigma;}}}\left(\mathrm{\&kappa;}(-{\mathrm{\&beta;}}_p)W_0{q}_p\right)\left(\stackrel{\&RightArrow;}{x}\right)---\left(46\right)$

因而，像平行光束情形中一样，P个扇形光束投射的加权后投射可被表达为加权后零反投射图像之和，如图4可见，其中B₀＝W₀。实际上，方程式(43)和(44)中定义的反投射运算符的接近相似性更普遍地应用于从一般形式的投射重建两维或更多维的函数，其中函数 和

适当定义。本发明的诸方法扩展到这些具有适当定义的坐标变换K的其它应用。

1.2快速分层反投射

用于扇形光束几何结构的快速分层反投射算法与平行光束情形的相似。它使用了由投射角为B彼此接近的投射所形成的中间图像可稀疏采样的事实，来组合三重分层结构中的扇形光束投射。框图如图18所示。相似的算法可对双重或其它(可能混合)基数结构得出。

支配图18中基础连续图像的组合的方程式如下：

$\begin{array}{c}{I}_{l,m}\stackrel{\mathrm{\&Delta;}}{=}{W}_0{q}_m---\left(47\right)\\ I_{l+1,m}\stackrel{\mathrm{\&Delta;}}{=}\mathrm{\&kappa;}\left({\mathrm{\&delta;}}_{l,3m-2}\right)I_{l,3m-2}+I_{l,3m-1}+\mathrm{\&kappa;}\left({\mathrm{\&delta;}}_{l,3m}\right)I_{l,3m},l=\mathrm{1,2},...,L---\left(48\right)\\ I_{L+\mathrm{1,1}}\stackrel{\mathrm{\&Delta;}}{=}{I}_{L,1}+\mathrm{\&kappa;}(-\mathrm{\&pi;}/2)I_{L,2}+\mathrm{\&kappa;}(-\mathrm{\&pi;})I_{L,3}+\mathrm{\&kappa;}(-3\mathrm{\&pi;}/2)I_{L,4}---\left(49\right)\end{array}$

在随后所述的实施例中，假设源角度在[0，2π)中均匀问隔如下：

${\mathrm{\&beta;}}_3=-\frac{\mathrm{\&pi;}}{4}(1-\frac{1}{9})+{\mathrm{\&Delta;}}_{\mathrm{\&beta;}}(j-1)$ 其中△_β＝π/(2·3^L)

然后中问旋转角使用方程式(22)像平行光束情形一样进行选择，其中Δ_β替换Δ_θ。

具有图18框图的算法的三层中的中间图像的构成扇形如图20所示。

1.2.1扇形光束采样方案

扇形光束算法的一个实施例使用类似于平行光束情形导出的算法中的中间图像的采样方案，这接下来描述。可选方案和改进在后面描述。

单个扇形光束反投射图像(由单个扇形光束的加权反投射产生的图像W₀q)具有可进行稀疏采样的结构。可从方程式(43)导出(W_βq)(V(β，γ，T))＝W(T)q(γ)。因此在方位为β＝0的扇形中由γ标示的射线上的特定T＝T₁上的单个样本足以沿整条射线指定图像的值：当T是从源到射线上点的距离时，它仅与W(T)成比例。

算法中的中间图像是这种单扇形光束的反投射图像的若干旋转后版本之和，每一个都从构成投射中生成。我们将在各构成投射所占角度间隔中央的角度上的扇形指为中央构成扇形。这可对应于实际投射——例如在三重算法的情形中，或对应于虚拟投射——例如在双重算法的情形中。例如，图20(a)中的扇形是图20(b)和20(c)的中央构成扇形。

回想在平行光束算法中，中间图像在与中央成分对齐的笛卡尔采样图上采样。在基于旋转的平行光束情形中，中间图像沿平行的竖直射线采样。样本沿这些平行射线中的每一条的间隔被选为确保图像的其它构成投射足以沿该射线采样(根据Nyquist采样标准)。该必要间隔实际上等于竖直射线的交点的间隔，其中一组旋转后平行射线对应于极值构成投射，即视角中距离中央投射最远的投射。

在平行光束情形的直接模拟中，在此中间图像沿中央构成扇形的射线采样。沿中央扇形的射线的采样点是中央扇形与极值构成扇形的交点。这导致样本未沿每条射线均匀分布。考虑两个扇形 $V(\mathrm{\&beta;},{\mathrm{\&gamma;}}_j,\&CenterDot;){|}_{j=1}^J$ 和 $V({\mathrm{\&beta;}}^{\&prime;},{\mathrm{\&gamma;}}_j^{\&prime;},\&CenterDot;){|}_{j=1}^J.$ 假设等角扇形光束几何结构具有奇数J个检测器，γ_i＝Δ_γ·(j-(J+1)/2)。与β’扇形相交的β扇形的第γ条射线上的点(对应于扇形角γ_r)，由方程式 $V(\mathrm{\&beta;},{\mathrm{\&gamma;}}_r,T)=V({\mathrm{\&beta;}}^{\&prime;},{\mathrm{\&gamma;}}_j^{\&prime;},T^{\&prime;}){|}_{j=1}^J,$ 确定，其解为

$T_{\mathrm{\&beta;},{\mathrm{\&beta;}}^{\&prime;},{\mathrm{\&gamma;}}_r}\left({\mathrm{\&gamma;}}_j^{\&prime;}\right)=\frac{D[\sin (\mathrm{\&beta;}-{\mathrm{\&beta;}}^{\&prime;}-{\mathrm{\&gamma;}}_j^{\&prime;})+\sin \left({\mathrm{\&gamma;}}_j^{\&prime;}\right)]}{\sin (\mathrm{\&beta;}-{\mathrm{\&beta;}}^{\&prime;}+{\mathrm{\&gamma;}}_r-{\mathrm{\&gamma;}}_j^{\&prime;})}$

函数T_{β，β′，γr}(γ′)描述扇形V(β′，·，·)如何沿扇形V(β，·，·)的第r条射线变化。它传送有关沿扇形射线的变化采样速率的信息。斜率越大，则样本越稀疏。任何γ’上的本地采样速率与(dT/dγ′)^-1成正比。典型的T(γ’)由图22中的虚线所示。

对T_{β，β′，γr}(γ′)作两次更改就得到新的采样函数

1.中间图像要在半径R的盘内采样，其中边界为盘外每条射线上至少一个样本(具有T＜T_ST的一个样本和具有T＞T_END的另一个)。定义 ${\mathrm{\&gamma;}}_{\mathrm{END}}^{\&prime;}=T_{\mathrm{\&beta;},{\mathrm{\&beta;}}^{\&prime;},{\mathrm{\&gamma;}}_r}^{-1}\left(T_{\mathrm{END}}\right):$ γ’的值对应于T_END。具有相邻源角度的扇形可因为T＞T_END而缺少交点。为了校正之，对T＞T_END， 的斜率固定为T_{β，β′，γr}(γ′_END)的斜率。

2.最终图像使用单元采样间隔在笛卡尔图上采样。这是中间图像需要采样的最小间隔。获得沿第γ条射线的单元采样速率的点γ′通过解方程式d/(dγ)T_{β，β，γr}(γ)＝＝1/Δ_γ来得出γ。为了保持单元采样间隔的约束， γ′≥γ′_的斜率被设置为1/Δ_γ。

图22示出典型β和γ的T_{β，β′，γr}(γ)和 它还显示是集 $\{{\tilde{T}}_{\mathrm{\&beta;},{\mathrm{\&beta;}}^{\&prime;},{\mathrm{\&gamma;}}_r}\left(j{\mathrm{\&Delta;}}_{\mathrm{\&gamma;}}\right):j\&Element;Z\}$ 的射线上实际采样点的T值。整数j被选择为在图像边界的每一侧上都有至少一个样本的边界。

显然，由在此概述的原理规定的样本点位置根据扇形的几何形状、构成投射集N_l，m的方程式(21)、中间图像的选定旋转角δ_l，m(对于等距离分布视角情形，在方程式(22)中给出)计算，并从 中选择。

用于扇形光束采样的可分旋转和上内插

如在概述中所述，上采样运算、以及与旋转运算组合的上采样可被分解成在计算上有效的一维重新采样运算。

1.2.2基于本地傅立叶结构的采样方案

假设图像 $f\left(\stackrel{\&RightArrow;}{x}\right):R^n\&RightArrow;R,$ 我们想要找出一采样函数 $t\left(\stackrel{\&RightArrow;}{x}\right):R^n\&RightArrow;R^n,$ 使

具有较小的基本带宽，因此可在点集 $\{t\left(\stackrel{\&RightArrow;}{m}\right):\stackrel{\&RightArrow;}{m}{\&Element;Z}^n\}$ 上采样。我们发现该采样函数如下。如将要简述的，知道了 如何由其构成投射组成，我们发现描述函数f应如何在点

上局部采样的矩阵函数 $\mathrm{\&upsi;}\left(\stackrel{\&RightArrow;}{x}\right)=\left[\begin{array}{cc}{\mathrm{\&upsi;}}_{11}\left(\stackrel{\&RightArrow;}{x}\right)& {\mathrm{\&upsi;}}_{12}\left(\stackrel{\&RightArrow;}{x}\right)\\ {\mathrm{\&upsi;}}_{21}\left(\stackrel{\&RightArrow;}{x}\right)& {\mathrm{\&upsi;}}_{22}\left(\stackrel{\&RightArrow;}{x}\right)\end{array}\right].$ 然后我们在图像域上组合该矩阵函数，以得到采样函数

在我们的算法中，我们知道中间图像f是这样的形式：对于某些角 ${\mathrm{\&delta;}}_p\&Element;[{\mathrm{\&beta;}}_{\min },{\mathrm{\&beta;}}_{\max }],f={\mathrm{\&Sigma;}}_p\mathrm{\&kappa;}\left({\mathrm{\&delta;}}_p\right)W_0{q}_p.$ 略去W(T)的加权，并假设投射q_p实际上以Nyquist速率采样，我们知道图像f中每个点上的光谱支持的局部结构。因为每个带宽受限投射q_p的扇形反投射，从单个投射进行扇形反投射的图像在扇形的辐条方向上具有可略去的空间带宽，而在垂直方向上其带宽与距扇形顶点的距离成反比，如图24所示。当一组这些扇形被旋转并累加在一起，该结果图像的局部光谱支持是各个扇形的光谱支持的合并。例如当构成扇形具有β_min和β_max之间的源角度时，如图25所示，图像域中点 上的局部光谱支持是方位在θ_min和θ_max之间具有半径、以及角度θ上径向带宽ΩR/r_θ的弯曲楔形。此处Ω是在图像平面中央经反投射的投射的空间带宽，假设以Nyquist速率采样各投射。在等角情形中，Ω是投射的样本数量除以弧通过经反投射的投射所覆盖图像平面的原点的长度。数学上光谱支持是

$\{({\cos \mathrm{\&theta;\&omega;}}_{\mathrm{\&theta;}},{\sin \mathrm{\&theta;\&omega;}}_{\mathrm{\&theta;}}):\mathrm{\&theta;}\&Element;[{\mathrm{\&theta;}}_{\min },{\mathrm{\&theta;}}_{\max }]$ 和 ${\mathrm{\&omega;}}_{\mathrm{\&theta;}}\&Element;[-{\mathrm{\&Omega;}}_{\mathrm{\&theta;}},{\mathrm{\&Omega;}}_{\mathrm{\&theta;}}]\}$ (50)

其中：

${\mathrm{\&theta;}}_{\max }=\tan \frac{R\sin {\mathrm{\&beta;}}_{\max }-x_1}{R\cos {\mathrm{\&beta;}}_{\max }+x_2}$

${\mathrm{\&theta;}}_{\min }=\tan \frac{R\sin {\mathrm{\&beta;}}_{\min }-x_1}{R\cos {\mathrm{\&beta;}}_{\min }+x_2}$

Ω_θ＝RΩ/r_θ

$r_{\mathrm{\&theta;}}=\sqrt{{(x_1-R\sin \mathrm{\&theta;})}^2+{(x_2-R\cos \mathrm{\&theta;})}^2}$

给定图像中点 上的局部两维傅立叶结构的该认知，该点上的矩阵函数 是这样的采样矩阵：如果该点上的光谱支持在整个图像上都是统一的(诸如平行光束情形)，则会有效采样图像。在感兴趣的中间图像中，这产生两个不同的小带宽和大带宽采样方向。我们将第一个坐标固定为小带宽方向。

我们通过 ${\left(\left[\begin{array}{cc}{s}_1& 0\\ 0& s_2\end{array}\right]\left[\begin{array}{cc}1& {\mathrm{\&alpha;}}_1\\ 0& 1\end{array}\right]\left[\begin{array}{cc}1& 0\\ {\mathrm{\&alpha;}}_2& 1\end{array}\right]\right)}^T$ 来近似该采样矩阵函数 此时 ${\mathrm{\&alpha;}}_2=-\frac{{\mathrm{\&Omega;}}_{\min }{\sin \mathrm{\&theta;}}_{\min }+{\mathrm{\&Omega;}}_{\max }{\sin \mathrm{\&theta;}}_{\max }}{{\mathrm{\&Omega;}}_{\min }{\cos \mathrm{\&theta;}}_{\min }+{\mathrm{\&Omega;}}_{\max }{\cos \mathrm{\&theta;}}_{\max }}$ 和 ${\mathrm{\&alpha;}}_1=\frac{{\mathrm{\&Omega;}}_{\min }{\cos \mathrm{\&theta;}}_{\min }-{\mathrm{\&Omega;}}_{\max }{\cos \mathrm{\&theta;}}_{\max }}{{\mathrm{\&Omega;}}_{\max }{\sin \mathrm{\&theta;}}_{\min }+{\mathrm{\&Omega;}}_{\min }{\sin \mathrm{\&theta;}}_{\min }}$ $S_1={\left(\left[\begin{array}{cc}1& 0\end{array}\right]\left[\begin{array}{cc}1& {\mathrm{\&alpha;}}_1\\ 0& 1\end{array}\right]\left[\begin{array}{cc}1& 0\\ {\mathrm{\&alpha;}}_2& 1\end{array}\right]\left[\begin{array}{c}{\cos \mathrm{\&theta;}}_{\mathrm{mid}}\\ {\sin \mathrm{\&theta;}}_{\mathrm{mid}}\end{array}\right]{\mathrm{\&Omega;}}_{\mathrm{mid}}\right)}^{-1}$ 和 $S_2={\left(\left[\begin{array}{cc}0& 1\end{array}\right]\left[\begin{array}{cc}1& 0\\ {\mathrm{\&alpha;}}_2& 1\end{array}\right]\left[\begin{array}{c}{\cos }_{\max }\\ {\sin }_{\max }\end{array}\right]{\mathrm{\&Omega;}}_{\max }\right)}^{-1}---\left(51\right)$

此处 ${\mathrm{\&Omega;}}_{\min }\stackrel{\mathrm{\&Delta;}}{=}{\mathrm{\&Omega;}}_{{\mathrm{\&theta;}}_{\min }}$ 且 ${\mathrm{\&Omega;}}_{\max }\stackrel{\mathrm{\&Delta;}}{=}{\mathrm{\&Omega;}}_{{\mathrm{\&theta;}}_{\max }}.$ 角度θ_mid指与源角度(β_min+β_max)/2相对应的投射角度，且 ${\mathrm{\&Omega;}}_{\mathrm{mid}}\stackrel{\mathrm{\&Delta;}}{=}{\mathrm{\&Omega;}}_{{\mathrm{\&theta;}}_{\mathrm{mid}}}.$ 通过光谱支持上几何变量选择的这些参数确保对采样矩阵变换的光谱支持受限于[-π，π]×[-π，π]。

整个图像的采样函数通过在整个图像上求该采样矩阵函数的积分来寻找——解微分方程集： $\frac{{\mathrm{dt}}_i}{{\mathrm{dn}}_j}={\mathrm{\&upsi;}}_{\mathrm{ij}}\left(t\left(\stackrel{\&RightArrow;}{n}\right)\right),i,j=\mathrm{1,2}.$ 这可得出用数字表示的解。即使未使用该确切指定的图，局部傅立叶支持分析将估计任何采样图的有效性。

一些中间图像的结果采样图如图26(A)和26(B)所示。对于扇形光束情形，该局部的基于傅立叶方法产生与因前述基于交点方法所导致的相似的采样图。从一个采样图重新采样成另一个的前述可分方法也可用于这些图。该局部傅立叶采样方法可直接应用，以寻找用于直线上、任意维上曲线或平面的任意投射几何结构的采样方案。在平行光束几何结构的情形中，它减少成先前所述的采样图。

1.2.3可选采样方案

样本位于其顶点在源轨迹上的扇形上的采样方案(称为“采样扇形”)具有一些缺点。对每条射线单独选择采样点，导致在两维意义上并非必然最优的方案。尽管算法中的最终图像在均匀矩形网格上用单位间隔采样，但是使用上述方案的中间图像比图像特定区域内(例如接近扇形顶点)所必需地更密集地采样。为了调整这个，上述方案可被更改成使该过采样区域中的采样稀疏。两种这样的可能性在图27中示出。这两种可能性结合这样的事实：图像在第一层中对扇形有效采样，而在最后一层则需要在矩形网格上采样。这些方案尝试结合从对扇形的采样到网格的逐步转换。

在图27(A)中，样本在其顶点位于距离原点比源半径更远的伪扇形上选择。在连续各层中，中间图像从更大范围的源角度的扇形中形成，且伪扇形顶点离原点的距离增大。在顶点在无穷远处的最终层上；即，射线是矩形网格的平行线。在图27(B)中，样本不位于扇形上，而是位于距离源位置越近越不发散的光束上。在连续各层上，光束变得越来越平行但越来越不发散。取两维(或者对于3D图像，三维)视点的这些或其它采样方案将有助于更快速的算法。

1.2.4其它优化

·在具有全轨迹的扇形光束情形中，源角度分布在[0，2π)中。这使得我们能利用不涉及内插而仅涉及像素的重新排列的-π/2、-π和-3π/2。

·对于算法的所有层，采样图中各点的位置可预先计算，并存储在查询表中，或者与数据处理一起在进行中计算。在任一情形中，可从23(A)和23(B)中看出，计算和/或存储可通过利用采样位置作为射线指数的函数的平滑变化来减少。

1.2.5过采样后的快速分层反投射

像平行光束情形中一样，过采样倍数γ＜0可用来改进准确度。在扇形光束情形中实现它的一种方法是使用扇形角度间隔Δ_β′＝Δ_β·γ来确定更密集的采样图。

1.2.6用于短期扫描和可选角度分布的算法

本发明的诸方法也可应用于所谓短期扫描数据获取格式，其中源角度的范围比[0，2π]小。不同之处在于反投射中使用的加权：额外的加权因数(例如Parker加权)被引入以说明短期扫描。类似地，诸算法归一化为源角度的、甚至不允许完整图像重建的多个源角度小组的不均匀分布。此外，算法中运算的顺序可用本领域技术人员所众所周知的各种方法更改(例如流水线技术)。例如，运算可被安排为单个投射一可用处理就开始的顺序。

Claims (41)

1.一种用于从投射(q₁...q_p)创建像素图像

的方法(图5A)，包括以下步骤：

(a)从选定投射(q₁...q_p)产生(100)中间图像(I₁，P)；

(b)对选定中间图像(I₁，P)执行数字图像坐标变换(102)，坐标变换的参数被选为说明从中产生中间图像的投射的视角、并说明所述中间图像的傅立叶特征；

(c)聚合在步骤(b)中生成的变换后中间图像的子集(104)，以生成聚合中间图像(I₂，P/2)；以及

(d)用递归方式重复步骤(b)、(c)，直到已处理和聚合所有投射和中间图像以形成像素图像

其中坐标变换参数被选为中间图像的聚合(104)可由稀疏样本用所需准确度表示。

2.如权利要求1所述的方法，其特征在于，所述聚合(104、108)通过累加数字图像来执行。

3.如权利要求1所述的方法，其特征在于，至少部分中间图像(I_n，m)通过反投射选定投射(q₁...q_p)在步骤(a)生成。

4.如权利要求1所述的方法，其特征在于，至少部分中间图像(I_n，m)的每一个都通过在步骤(a)反投射两个或更多个选定投射(q₁...q_p)而形成。

5.如权利要求1所述的方法，其特征在于，至少部分聚合中间图像(I_n，m)的每一个都通过在步骤(d)聚合三个或更多个选定变换后中间图像而形成。

6.如权利要求1所述的方法，其特征在于，数字图像坐标变换使用数字滤波来执行。

7.如权利要求1所述的方法，其特征在于，选定坐标变换包括数字图像旋转。

8.如权利要求1所述的方法，其特征在于，选定坐标变换包括数字图像切变(图10B 120、122)或切变缩放。

9.如权利要求1所述的方法(图15)，其特征在于，选定坐标变换包括数字图像的上采样(101、106)和/或下采样(109)。

10.如权利要求8所述的方法，其特征在于，所述数字图像切变通过一维线性数字滤波器来执行。

11.如权利要求9所述的方法，其特征在于，所述数字图像上采样和/或下采样通过一维线性数字滤波器来执行。

12.如权利要求10所述的方法，其特征在于，至少部分所述数字滤波器是移位不变的。

13.如权利要求11所述的方法，其特征在于，至少部分所述数字滤波器是移位不变的。

14.如权利要求12或13所述的方法，其特征在于，至少部分所述数字滤波器是递归的。

15.如权利要求12或13所述的方法，其特征在于，至少部分所述数字滤波器使用快速傅立叶变换(FFT)来实现。

16.如权利要求1所述的方法，其特征在于，选定过采样应用于选定中间图像、和/或变换后中间图像、和/或聚合中间图像。

17.如权利要求1所述的方法，其特征在于，使用了非笛卡尔采样图。

18.如权利要求1所述的方法，其特征在于，多个选定坐标变换可在分层算法的一个层内或跨相邻层组合。

19.一种用于沿线、曲线或表面的集合从投射(q₁...q_p)创建像素图像

的方法(图31)，包括以下步骤：

(a)产生(184)中间图像(I₁，m)；

(b)对选定中间图像执行数字图像重新采样(186)，所选定样本的位置被选为说明选定投射的视角、并说明所述中间图像的傅立叶特征；

(c)聚合(190)重新采样后中间图像的选定子集，以生成聚合中间图像(I_z，m)；以及

(d)用递归方式重复步骤(b)、(c)，在递归的每一层上增加中间图像的样本密度，直到已处理和聚合所有投射和中间图像以形成所述像素图像；

其中采样方案被选为重新采样后的中间图像的聚合可由稀疏样本用所需准确度表示。

20.如权利要求19所述的方法，其特征在于，至少部分中间图像通过选定投射的加权反投射(图19 180、182)在步骤(a)生成。

21.如权利要求19所述的方法，其特征在于，至少部分中间图像的每一个都通过在步骤(a)加权反投射两个或更多个选定投射(图19 180、182)而形成。

22.如权利要求19所述的方法，其特征在于，至少部分聚合中间图像的每一个都通过在步骤(d)聚合三个或更多个选定变换后中间图像而形成。

23.如权利要求19所述的方法，其特征在于，所述中间图像具有位于线、曲线或表面的族上的样本。

24.如权利要求19所述的方法，其特征在于，所述数字图像重新采样通过使用位于线、曲线或平面族的交点处的中间采样方案，由一系列低维数字滤波运算来执行。

25.如权利要求19所述的方法，其特征在于，选定过采样度应用于所选定的重新采样后的中间图像、以及聚合后的中间图像。

26.如权利要求19所述的方法，其特征在于，所述聚合通过累加数字图像来执行。

27.如权利要求19所述的方法，其特征在于，所述重新采样和聚合可跨连续层组合。

28.如权利要求30所述的方法，其特征在于，至少一个中间图像在所包括的重新采样步骤(b)前后加权。

29.如权利要求19所述的方法，其特征在于，采样密度的改变通过数字滤波来完成。

30.一种用于从投射(q₁...q_p)创建像素图像

的方法(图18)，包括以下步骤：

(a)产生(99)多个中间图像(I₁，m)，其中至少一个中间图像对应于非笛卡尔和/或非周期性采样图；

(b)对选定中间图像执行数字图像上采样或下采样(106)；

(c)对上采样后/下采样后的中间图像执行数字图像坐标变换；

(d)聚合(110)在步骤(c)中生成的变换后中间图像的子集，以生成聚合中间图像；以及

(e)用递归方式重复步骤(b)、(c)和(d)，直到已处理和聚合所有投射和中间图像以形成像素图像；

其中数字图像坐标变换的至少之一用非笛卡尔和/或非周期性的采样图来执行，且坐标变换参数被选为中间图像的聚合可由稀疏样本用所需准确度表示。

31.如权利要求30所述的方法，其特征在于，所述聚合通过累加数字图像来执行。

32.如权利要求30所述的方法，其特征在于，至少部分中间图像通过反投射选定投射在步骤(a)生成。

33.如权利要求30所述的方法，其特征在于，至少部分中间图像的每一个都通过在步骤(a)反投射两个或更多个选定投射而形成。

34.如权利要求30所述的方法，其特征在于，至少部分聚合中间图像的每一个都通过在步骤(c)聚合三个或更多个选定变换后中间图像而形成。

35.如权利要求30所述的方法，其特征在于，数字图像坐标变换使用数字滤波来执行。

36.如权利要求30所述的方法，其特征在于，选定坐标变换包括数字图像旋转。

37.如权利要求30所述的方法，其特征在于，所述数字图像重新采样的上采样和/或下采样通过一维线性数字滤波器来执行。

38.如权利要求37所述的方法，其特征在于，至少部分所述数字滤波器是移位不变的。

39.如权利要求37所述的方法，其特征在于，至少部分所述数字滤波器是递归的。

40.如权利要求38所述的方法，其特征在于，至少部分所述数字滤波器使用快速傅立叶变换(FFT)来实现。

41.如权利要求30所述的方法，其特征在于，选定过采样应用于选定中间图像、和/或变换后中间图像、和/或聚合中间图像。

Images