CN103052959A - 具有优化递归的高效计算层析摄影术 - Google Patents

Abstract

一种系统，包括：一个或多个发射器，将激发能量发射到被观察对象；一个或多个检测器，响应于发射到被观察对象中的激发能量，产生投影空间数据，投影空间数据编码了一个或多个检测器接收的能量；控制器，控制一个或多个发射器发射激发能量，并且控制一个或多个接收器产生投影空间数据；以及具有至少一个处理器的图像重构器，接收投影空间数据，并且通过以下操作处理投影空间数据：计算第一量，所述第一量表征投影空间数据与预测的投影数据之间的差，其中投影空间数据与预测的投影数据之间的差与投影增益相关联；将编码了所述差的第一数据记录在数据存储设备中；基于来自在前迭代的记录的第一数据，计算修正的第一量；使用修正的第一量来计算更新值；并且使用更新值重构表示被观察对象的对象空间图像。

Description

具有优化递归的高效计算层析摄影术

相关申请的交叉参考

本PCT专利申请要求2010年7月16日递交的美国专利申请No.12/838,205的优先权，其全部内容通过引用合并于此。

技术领域

本发明涉及图像重构领域。

背景技术

层析摄影术，即，基于n维空间中的区域的m维投影来计算(估计)该区域中的密度(表示为像素；通常0＜m＜n)，通常分为两个主要类别：滤波的反向投影(back-projection FBP)以及迭代重构(IR)，迭代重构(IR)也被称作代数重构技术(ART)的变型。FBP和传统的IR技术有特定缺点，尤其是在一个投影通常包含大概从2000个像素(一维)一直到2000×2000个像素(二维)时通常需要计算10⁶-10¹⁰个图像元素(即，体素)的现代设置中。

与FBP相关联的主要缺点包括需要多个投影来实现有限的定量精度。投影的数目典型地数以百计，但是投影并未有效加以使用。例如，计算的密度估计可以为负值，而已知物理强度不能为负。有限的定量精度有时可以通过后续迭代细化来提高。与FBP相关联的其他缺点包括不能根据体素位置改变数据权重，如在Wood，S.L，Morf，M.，A fastimplementation of a minimum variance estimator for computerizedtomography image reconstruction，IEEE，Trans，on Biomed.Eng.，Vol.BME-28，No.2，pp.56-68，1981中所讨论的，并且不能高效地考虑各种约束。

由于FPB方法通常依赖于快速傅里叶变换，因此FBF方法的强项包括其计算速度、中心切片理论、以及可能使用适合并预计算的滤波器。

FBP和传统IF技术的备选技术使用矩阵运算，并且可以适用于体积或画面元素(体素或像素)的数量数以千计的“小型”问题。然而，对于典型的层析摄影重构设置，通过这些基于矩阵的技术在可预见未来中不能处理计算负担。这是因为计算机运算次数与N的幂成比例增加。例如，在3D情况下，运算次数成比例增加快过N⁵(其中，N是表示感兴趣区域的重构立方体的边的像素数目)，如Wood和2001年12月18日的美国专利No.6,332,035中所讨论的。

IR技术的优点包括能够考虑约束，具体地能确保密度估计为非负。确保非负的密度估计可以引起显著的图像对比度提高以及更好的定量结果。

与传统IR技术相关联的缺点包括需要对逆运算以及正向投影运算进行重复求解，以便获得用于后续迭代的校正项。这些运算的许多递归配对所需的伴生物通常使IR技术较慢，并且依赖于特定的不稳定实现方式。

因为这些传统IR技术使用各种(迭代)优化非负，因此这些技术朝向解的收敛速率依赖于特定技术对赫赛(Hessian)矩阵的执行效果，赫赛矩阵是要优化的目标函数(即相对于图像元素的性能准则)的第二导数矩阵。最流行的技术可以是IR的变型。然而由于用于图像重构的赫赛矩阵的尺寸较大，因此赫赛矩阵的结构(以及其逆矩阵的结构)通常被忽略或粗略近似。此外，由于赫赛矩阵的特征值的广泛分布，因此当前优化技术超出一定次数的迭代(典型地，数以十计到数以百计)倾向于示出没有改进，并且仅处理较大特征值。

这些算法的多网格变型可以有帮助，但是最终仍会失败，这是因为赫赛矩阵的尺寸涉及细化网格。多网格分辨率在这里是指随着迭代的执行所使用的逐步细化的分辨率。

现代层析摄影术设置在一个投影包含大概从2000个像素(一维)一直到65,000×65,000个像素(二维)时通常需要计算10⁶-10¹⁴个图像元素(即，体素)。诸如线性化和次问题维度的缩减等近似技术可以帮助使针对高效率计算层析摄影术(HECT)的计算可行。然而，根据现代层析摄影术设置的图像重构例如由于典型的非线性估计而复杂。

发明内容

本发明可以提供一种系统，包括：一个或多个发射器，将激发能量发射到被观察对象；一个或多个检测器，响应于发射到被观察对象中激发能量，产生投影空间数据，投影空间数据编码了一个或多个检测器接收的能量；控制器，控制一个或多个发射器发射激发能量，并且控制一个或多个接收器产生投影空间数据；以及具有至少一个处理器的图像重构器，接收投影空间数据，并且通过以下操作处理投影空间数据：计算第一量，所述第一量表征投影空间数据与预测的投影数据之间的差，其中投影空间数据与预测的投影数据之间的差与投影增益相关联；将编码了所述差的第一数据记录在数据存储设备中；计算与至少一个脉冲响应相对应的第二量，每个脉冲响应对应于单位强度的参考体素；使用修正的第一量和第二量的逆函数来计算更新值，其中更新值与体素增益相关联；并且使用更新值重构表示被观察对象的对象空间图像。

本发明可以提供一种工作站，包括：一个或多个处理器；以及一个或多个数据存储设备，存储要由一个或多个处理器执行的软件，所述软件包括：接收输入数据的软件代码，输入数据表示从对象空间投影到投影空间的对象；计算第一量的软件代码，第一量表征输入数据与预测的投影数据之间的差，其中输入数据与预测的投影数据之间的差与投影增益相关联；计算第二量的软件代码，第二量对应于至少一个脉冲响应，每个脉冲响应对应于对象空间中位置处单位强度的参考体素；使用第一量和第二量的逆函数计算更新值的软件代码；在一个或多个数据存储设备或另外一个或多个数据存储设备中记录编码了更新值的第二数据的软件代码；基于来自在前迭代的记录的第二数据修正更新值的软件代码；以及使用修正的更新值重构表示被观察对象的对象空间图像的软件代码，其中更新值与体素增益相关联。

本发明可以提供一种由执行一个或多个数据存储设备上存储的软件代码的一个或多个处理器来实现的方法，该方法包括：从以下中的至少一个接收输入数据：投影图像获取系统、一个或多个数据存储设备、另外一个或多个数据存储设备、一个或多个处理器执行的软件仿真、由另外一个或多个处理器执行的软件仿真，输入数据表示从对象空间投影到投影空间的对象；计算第一量，第一量表征分辨率网格上的输入数据与分辨率网格上的预测的投影数据之间的差，其中，预测的投影数据与投影增益相关联；在一个或多个数据存储设备或另外一个或多个数据存储设备中记录编码了所述差的第一数据；基于来自在前迭代的记录的第一数据，计算修正第一量；计算与至少一个脉冲响应相对应的第二量，每个脉冲响应对应于对象空间中位置处单位强度的参考体素；使用第一量和修正的第一量之一以及第二量的逆函数计算更新值，其中更新值与体素增益相关联；在一个或多个数据存储设备或另外一个或多个数据存储设备中记录编码了更新值的第二数据；基于来自在前迭代的记录的第二数据，修正更新值；使用更新值重构表示被观察对象的图像；并且向输出设备输出重构图像，其中输出设备是以下中的至少一个：一个或多个数据存储设备、另外一个或多个数据存储设备、显示设备、和打印设备；其中，根据一个或多个数据存储设备上存储的软件代码，由一个或多个处理器执行接收输入数据、计算第一量、记录第一数据、计算修正的第一量、计算第二量、记录第二数据、修正更新值、重构图像和输出重构图像。

附图说明

现在结合附图描述示例性实施例。

图1示出了根据本发明示例性实施例的系统图。

图2示出了根据本发明示例性实施例的成像系统。

图3示出了根据本发明示例性实施例的流程图。

图4A和4B分别示出了基于传统滤波的反向投影(FBP)重构和根据本发明示例性实施例的方法的重构的细胞图像。

图5示出了仿真人体模型(phantom)的正面样本投影。

图6A和6B分别示出了基于仿真迭代重构技术(SIRT)以及根据本发明示例性实施例的方法的重构误差。

图7A示出了碘对比染料注入到跳动液力冠状人体模型中期间的一个X射线投影。

图7B示出了基于根据本发明示例性实施例的方法而获得的人体模式重构图像。

图8A示出了根据本发明示例性实施例的扩展流程图。

图8B示出了根据本发明示例性实施例的另一扩展流程图。

图9A至9D分别示出了10、50、130和170度轴倾角处的输入投影数据。

图10A至10D分别示出了10度轴倾角处针对迭代1、2、5和6的预测投影矩阵。

图11A至11D分别针对迭代1、2、5和6示出了原始输入投影数据与预测投影矩阵之间的残差。

图12A至12D分别针对迭代1、2、5和6示出了10度轴倾角处生成的投影增益。

具体实施方式

以下详细讨论本发明的示例性实施例。在描述示例性实施例时，为了清楚起见采用特定术语。然而，本发明并不意在受限于所选的特定术语。本领域技术人员应认识到，在不背离本发明广义构思的前提下可以采样其他等同部件并且开发其他方法。本文所引用的所有参考文献通过引用合并，也可以单独合并。

图1示出了根据本发明一些实施例的系统图。可以使用图像重构器105将投影的像素数据104变换成要在输出设备107上可视化的重构图像106，投影的像素数据104包含编码了对象的内部结构的信息。投影的像素数据104可以来自实验获取101、或者计算机上的仿真103、或者包含记录的投影像素数据104的数据存储设备102，记录的投影像素数据104来自在前实验或仿真。实验获取101、数据存储设备102、和仿真103可以直接或经由网络(例如，局域网(LAN)、广域网(WAN)、互联网等)远程地将投影的像素数据104提供给图像重构器105。尽管这里仅示出了投影的像素数据104，但是通常数据可以是将任何形式的激发能量传递到被观察对象中的结果，因此除了包括来自计算层析摄影术的数据以外，还包括来自以下的数据：例如，电子显微术、磁共振(MR)、正电子发射层析摄影术(PET)、单正电子发射计算层析摄影术(SPECT)、超声、荧光、多光子显微术(MPM)、光相干层析摄影术(OCT)等。数据存储设备1 02例如可以是计算机存储器、CD-ROM、DVD、磁光(MO)盘、硬盘、软盘、压缩盘、闪速驱动器等。

图2示出了根据本发明一些实施例成像系统。该成像系统可以包括：控制器205；接口206(例如，图形用户界面、键盘、操纵杆等)，与研究者211(例如，人类用户、计算机等)信号通信，并且与控制器205同步；发射器201，响应于来自控制器205的控制信号，产生和发射激发能量202；以及检测器207，配置为产生投影的像素数据104。投影的像素数据104可以按照数字格式存储在存储设备中。成像系统还可以包括：平移或旋转台204，配置为在其上容纳对象203，并且可操作于相对于发射器201和检测器207平移或旋转；以及图像重构器105、209以电方式或者经由可选数据传送器208(例如，互联网、其他网络、或数据总线)耦合至检测器207和控制器205。例如，数据总线可以是电线的子系统，在计算机内部的计算机部件之间或者在计算机之间传送数据。尽管，相连部件之间的数据流是单向的，但是相连部件之间的通信也可以是双向的。

激发能量202可以辐射能量的形式，例如，X射线能量、电磁(EM)能量、光能量、红外(IR)能量、粒子能量(例如，电子、中子、原子束)、振动能量(例如超声波)等。激发能量202可以照射到对象203上，对象203例如可以是人体模型、患者、标本或其任何组合。激发能量202可以由对应能量类型的发射器201发射。激发能量202可以通过对象203传播，并且其一部分可以由适当的检测器207接收。检测器207可以将接收到的能量转换成可测量的电信号，还可以将测量的电信号转换成数字格式的投影像素数据104。

控制器205可以是将发射器201和检测器207相连的电路，并且可以向发射器201和检测器207发送控制信号，以使激发能量202的发射与检测器207的操作同步。电路可以是模拟、数字或混合形式。控制器205还可以是具有一个或多个处理器以及一个或多个存储器的计算机，以向发射器202和检测器207发送控制信号。

图像重构器105、209可以响应于控制器205，并且接收投影的像素数据104，基于根据本发明一些实施例的方法重构对象203的图像，以高计算效率产生对象的高保真图像。图像重构器105、209可以是具有一个或多个数据存储设备的计算机，一个或多个数据存储设备存储根据本发明示例性实施例操作计算机的软件代码。计算机可以包括一个或多个处理器以及一个或多个存储器，来读取和执行一个或多个数据存储设备上存储的软件代码。在另一示例性实施例中，图像重构器105、209可以包括一个或多个根据本发明示例性实施例可操作的程序存储和执行设备。图像重构器105、209可以产生重构图像106。

输出设备107、210可以接收重构图像106。输出设备107、210可以是可视化设备或数据存储设备。可视化设备例如可以是显示设备或打印设备。示例显示设备可以包括例如阴极射线管(CRT)、发光二极管(LED)显示器、液晶显示器(LCD)、数字光投影(DLP)监视器、真空荧光显示器(VFD)、表面传导电子发射显示器(SED)、场发射显示器(FED)、硅上液晶(LCOS)显示器等。示例打印设备可以包括，例如基于墨粉的打印机、液态喷墨打印机、固态墨打印机、染料升华打印机、和无墨打印机(例如，热感打印机和紫外线(UV)打印机)等。打印设备可以按照3维(3D)进行打印。

成像系统还可以包括研究者211。研究者211可以是人类或编程计算机，来从输出设备210或图像重构器209接收重构图像106，且然后应用算法(例如，与编程例程、人工智能、机器学习等)来提取关于对象203的诊断信息、或者针对发射器201、检测器207、台204、图像重构器209等的细化调谐控制参数。在一些实施例中，接口206或输出设备210不是必要的。在一些实施例中，研究者211、控制器205和图像重构器105、209可以驻留在同一计算机或分离的计算机上。

本发明的一些实施例可以提供包括一个或多个处理器的工作站，一个或多个处理器配置为以类似于图像重构器105、209的方式重构图像。工作站可以接收来自成像系统、数据存储设备或计算机中的至少一个的输入数据。输入数据可以经由数据总线、电缆、有线网络、无线网络等接收。工作站还可以包括输出设备107、210来接收重构图像。输出设备可以是数据存储设备、显示设备、打印设备等。如下讨论示例数据存储设备、显示设备和打印设备。

图3示出了根据本发明示例性实施例的方法的流程图。

在框301，可以选择适当的极限(微调)函数f(x)，例如：

-a＜f(x)＜b    方程(1)

其中，-a＜0＜b，并且对于x＝0的一些邻域，f(x)根据以下方程近似于x：

|f(x)-x|＜ε，对于适当小的ε＞0    方程(2)

函数f可以微调过度值，特别是在早期的迭代下。例如，一个这样的函数是atan(x)。极限函数可以放松，并且不存在严格的截止值a和b。另一示例极限函数可以是h(x)＝c·x+(1-c)·atan(x)，0＜c＜1。Box，G.E.P.，Jenkins，G.M.，Time Series Analysis：Forecasting and Control，SanFrancisco，CA：Holden-Day，1976包含可以用于函数f的非线性数据变换的选择。

在框302中，m维(m-D)像素的输入投影数据104可以针对每个给定投影方向置于改变分辨率的一系列分辨率网格(例如，每个后续分辨率网格具有比在前分辨率网格更少的像素，而同时保持图像范围)上。例如，一系列分辨率网格中的分辨率网格可以具有一系列分辨率网络中在前分辨率网格的分辨率的一半。最粗尺寸例如可以是单个m-D像素。

应当注意，通常1＜m＜n，并且m对于所有投影可以相等。然而通常，m可以随着每个投影而变化，并且m＞＝0(例如，对于特定噪声估计问题)。

同样在框302中，n维(n-D)体素的初始集合可以置于n-D空间中，使得它们的投影覆盖一系列分辨率网格上的可用m-D投影像素。例如，针对这些初始体素的适合值集合可以是与平均投影密度相对应的值；并且对于每个体素j，可以设置v_Lj＝log(v_j)。除了对数以外还可以使用单调函数。例如，可以将初始单n-D体素密度v_i设置为等于平均投影密度。

在框303中，对于本网格分辨率，n-D体素v和m-D像素p均可选地可以利用适合的核来平滑，其中m-D核近似为n-D核的投影。该运算可以对应于具有噪声观察的卡尔曼滤波器的预测符算法(没有更新)(其中v根据v_L获得)。例如，示例函数可以是v＝exp(v_L)。可以采用包括不不连续性的其他函数来解决附加密度约束。示例性平滑核可以具有sigma宽度的高斯形状，例如，在0.4和2个体素间隔之间。然而，稍后参照框307讨论可能的修改。

根据需要，可以经由针对特定值的映射来约束体素v。例如，在一些区域，先验已知体素密度为零。同样，对于每个小体素密度为了避免数值下溢，v_Lj可以被限制到极小的可忽略值。示例性极限值可以表示最大体素密度值

的大约为10^-15的密度。

在框304中，可以基于n-D空间中的n-D体素v集合来计算预测的投影m-D像素p_p。预测的投影m-D像素p_p可以在来自一系列分辨率网格的分辨率网格上。如果网格分辨率在在前分辨率网格上有所增加，则可以经由从更粗/在前分辨率网格的插值获得预测的投影和体素估计v。

在框305中，可以计算对分辨率网格上预测的投影数据与输入投影数据104之间的差加以表征的第一量。

如果输入投影数据104的投影像素值p_i大于零，则第一量例如可以是

u_i＝f(log(L_i))    方程(3)

其中，L是输入投影数据104的测量的p_i和预测投影数据的预测像素值p_pi之间的误差比，例如如以下方程(4)表示的：

L_i＝p_i/p_pi    方程(4)

例如，其他适合的关系也可以用于L_i，其中，L_i是随着p_i而增加并且随着p_pi的增加而减小的函数。这些函数可以捕获输入投影数据104和预测投影数据之间的误差测量。

以上log函数可以是误差比的对数，但是可以是输入与预测投影值的更一般函数(例如，发射电子显微术中丢失楔形(wedge)问题中使用的函数)。

如果输入投影数据104的投影像素值p_i不大于零，则第一量例如可以是

u_i＝-a    方程(5)

其中，a指示界限。

在框306中，可以计算与至少一个脉冲响应R_p相对应的量，所述至少一个脉冲响应R_p针对n-D空间中单位强度的参考体素。可以针对所有p个投影和反向投影，根据当前分辨率网格为每个期望参考n-D体素计算脉冲响应R_p。脉冲响应R_p可以依赖于参考n-D体素的位置。可以根据优选操作序列，针对每个m-D投影或者针对n-D体素计算脉冲响应R_p及其逆，优选操作序列继而可以依赖于期望精度。这类似于FBP。可以选择参考n-D体素的有限子集用于评价脉冲响应_Rp。在一些应用中，可以在重构空间的中心处选择最少一个n-D参考体素。

可以按照若干方式为对应的参考n-D体素计算脉冲响应，并且这些脉冲响应的逆可以用于与第一量的后续组合，如稍后关于框307所讨论的。

在示例性实施例中，n-D重构/对象空间的区域中特定位置处的脉冲响应R_p可以是考虑单位强度的参考体素时的线性系统。

当将特定位置处该单个参考n-D体素的中心投影到所有p个投影方向时，可以创建p个单个像素投影(不必要针对所有相同维度n)。当这些p个投影是针对重构的反向投影时，可以创建原始体素的模糊版本。模糊版本与原始体素线性相关，但是可以在整个n-D空间扩展。该模糊函数R_p可以是针对该特定位置的重构系统的脉冲响应(或者点扩展函数PSF)。

对于有限重构体积，该脉冲响应R_p可以依赖于n-D体素的位置。例如，当投影数目较少并且网格足够细化，脉冲响应R_p可以类似于蜘蛛(例如，示出了“蜘蛛腿”状扩展，或者从蜘蛛中心位置条状发散)。类似地，可以将邻近体素形成为“蜘蛛腿”状延伸，但是可以具有不同长度，这是因为有限的重构区域可以切断边界处的延伸。此外，在几何形状(例如，使用p锥面波束投影的锥面波束)中，“蜘蛛腿”的方向和相互角度可以随着位置而改变，因为“蜘蛛腿”可以指向p锥面的固定尖端。

然而为了近似的目的，这些脉冲响应函数的慢空间变化可以用于共享n-D空间中体素的足够小邻域上相同脉冲响应函数。“蜘蛛腿”的相似度的几何考虑可以确定共享相同近似脉冲响应的良好匹配体素区域的尺寸。例如，“蜘蛛腿”的尖端处相对于它们内容几个(例如，数十个)百分点的整体失配是完全可接受的，尤其在后续迭代可以校正这些误差的情况下。相反，FBP方法通常不考虑这种截断效果，而Wood中示出的矩阵公式化可以考虑。

同样在框306中，然后可以将针对给定体素的脉冲响应R_p与来自噪声偏差的共享相组合来获得R，R是第二量。可以测量输入投影数据以获得噪声偏差。该组合可以增加“蜘蛛”中心处的值。可以根据优选操作序列，为每个m-D投影或n-D体素计算R及其逆，优选操作序列继而可以依赖于期望精度。较小增加(例如，当假定噪声级较低时)可以显著增加对R的逆的精确计算的要求(例如对应于卡尔曼/维纳滤波器增益)。然而传统上，对多个分辨率网格的计算可以以适度/较低级别保持残差投影误差中的信号部分，因此支持有效率的R^-1计算。信号部分可以是重构未考虑的残差投影误差中的结构信息。

在框307中，然后可以对计算的第二量进行逆运算，以获得P，例如表示为P＝R^-1(代替使用逆赫赛)。该步骤可以在傅里叶域中执行。本发明的示例性实施例可以在傅里叶域中将逆函数表征为卡尔曼或维纳增益，如下所述。

例如，针对测量的输入投影数据的观察模型可以是Y＝HX+H，其中X可以表示对象密度，H可以表示系统变换函数，Y可以表示输入投影数据104。

在矩阵表示中，X可以根据可用数据Y来估计。数据Y可以包括来自测量噪声的贡献。X可以具有以下统计属性：

X₀＝E(X)                  方程(6)

V_X＝E[(X-X₀)·(X-X₀)^t]    方程(7)

并且投影噪声W统计可以表示为：

E[W]＝0                   方程(8)

V＝E[W·W^t]，并且      方程(9)

E[W·X^t]＝0            方程(10)

例如，如在Sage(Sage，A.P.，J.L.Melsa，Estimation Theory withApplications to Communications and Control，New York：McGraw-Hill，1971，page 262)中讨论的最小偏差(卡尔曼)增益K可以由如下方程给出：

K＝V_XH^t[HV_XH^t+V]^-1＝V_XH_tR_Y ^-1    方程(11)

或者由以下方程给出：

K＝{H^t[I+V(HV_XH^t)^-1]}^-1         方程(12)

所提供的V_XH^t是不可逆的。在所建议的计算的早期低分辨率迭代中满足该条件，其中，输入投影数据中的像素数目可以超过估计的图像体素的数目。在层析摄影术的上下文中，R_Y ^-1对于逆运算而言太大，但是方程(12)建议有用的H修改来获得实际K。通过方程(11)，提供满秩V可以从矩阵R_Y ^-1中去除可能的零特征值。HV_XH^t和V之间的关系可以确定图像X的估计对X的先验统一的依赖程度。

使用复频谱表示，对象的功率谱密度可以是S_XX，观察传递函数可以是H，投影噪声功率谱密度可以为S_VV，并且维纳(例如，参见Wiener，Norbert，Extrapolation，Interpolation，and Smoothing ofStationary Time Series.New York：Wiley，1949)增益可以表示为：

W＝H*S_XX/(H*S_XXH+S_VV)      方程(13)

＝[H(U+S_VV/(H*S_XXH))]^-1    方程(14)

对于层析摄影术应用而言频谱公式化在计算上是可行的。在层析摄影重构中，对象密度的不确定性(例如，S_XX或V_X可以相对较大)可以主导投影噪声(例如，S_VV或V可以相对较小)，并且找到维纳增益W(或卡尔曼增益K)的主要任务可以是找到对方程(13)(或者同样对方程(11))有用的近似。

因此，可以测量输入投影数据104以获得如上表示的V和V_X。例如可以基于测量V_X的来计算HV_XH^t和H*S_XXH。来自框306中提及的噪声偏差的贡献可以指由于例如V、HV_XH^t、或H*S_XXH的贡献。尽管可以直接计算针对卡尔曼增益的方程(12)和针对维纳滤波器增益的方程(14)，但是其他修改也是可能的。例如，对于测量噪声W的较小互相关(例如，V是对角线主导)以及投影对象密度的适中互相关(例如，HV_XH^t是对角线主导，特别是在迭代过程期间)，可以通过针对各个单独投影(或投影组)计算逐块计算增益来简化方程(12)或(14)的复合表达式。例如，也可以使用基于逐体素(n-D投影)或逐像素(m-D近似处理)修改的局部脉冲响应的逆来简化逐块处理。例如，局部脉冲响应的修改可以由方程(12)中所见的噪声与信号结构的元素V(HV_XH^t)^-1或者类似地通过方程(14)中所见的S_VV/(H*S_XXH)来给出。可以在体素到体素的基础上应用修改。可以通过增加“蜘蛛”中心区域的值来近似该修改(例如，在脉冲响应函数的峰值附近，对应于方程(12)中的V(HV_XH^t)^-1矩阵或方程(14)中的S_VV/(H*S_XXH)矩阵的(近)对角元素)。在W是白噪声并且投影的对象不确定性是估计均方根(RMS)信噪比为r的白噪声的情况下，例如可以通过因子(1+r²)来增加蜘蛛中心。在这种修改之后，例如可以根据方程(12)和(14)通过数值逆运算来计算卡尔曼增益和维纳增益。

应当注意，如上所述，R_p与在方程(11)和(13)中表示的上述矩阵(或传递函数)H相关。可以通过将用卡尔曼增益K或维纳增益W表示的单个步骤拆分成体素操作的两个步骤来改变处理序列。例如，两个步骤可以是首先反向投影投影数据，然后利用参考体素K的n-D滤波器

来滤波；一些周围邻近体素K可以使用生成的滤波器输出。两个步骤还可以是首先例如利用p个m-D滤波器P_iK(0＜＝i＜p)对投影数据进行预滤波，其中可以从被信噪比修改的n-D脉冲响应(PSF)R_PK的投影的逆获得每个P_iK，然后对滤波器输出进行反向投影；同样，一些周围邻近体素K(或其投影)可以使用生成的滤波器输出(更新)。因此，可以将方程(11)和(13)中示出的极大矩阵的逆矩阵分解成数目更小的组R_K的中等组的逆，其中，每个R_K可以在对象空间中仅缓慢变化，允许重复使用。此外，可以通过傅里叶技术非常有效率地针对参考体素计算每个R_K。

总体上，用于计算增益P的近似的两个方面可以有利于高统计效率。

首先，对于低频而言，卡尔曼/维纳滤波器增益(现在由增益P代替)可以通过在脉冲响应函数的原点处添加小常数来估计，脉冲响应函数可以在足够宽的频率范围上考虑信噪比。可以通过在原点处添加“窄”高斯钟而不是仅添加常数来修改脉冲响应函数中心附加的小邻域，以考虑例如“粉红”测量噪声。

其次，对于较高频率而言，与增益P相对应的滤波器可以如高频率信噪比所表征的被逐渐减弱/平滑。为了提高数值精度，可以对当前分辨率网格上图像密度更新和/或估计图像密度执行该减弱/平滑函数。在对一系列分辨率网格的迭代计算期间的每个阶段处，估计图像密度的低通滤波可以帮助逐渐减少高频图像分量，降低计算负担。

在框308中，可以使用第一量和P来计算更新值。例如，P可以与第一量相组合来获得针对给定体素(或像素)图像元素的未处理更新值i。例如，P可以与第一量卷积以获得未处理的更新值。随着P在n-D重构空间中缓慢变化，通过对给定体素的适合但不足够小邻域应用P作为近似，来节省计算时间。然后可以将类似于上述方程(1)和(2)中讨论的f构造的极限函数g应用于未处理更新值i，以获得修改的更新值i_g＝g(i)。这种修改可以防止过度校正。在本发明的示例性实施例的开发期间，许多实验不能直接使用未处理更新值。本发明人任何设想和开发极限函数及其恰到好处的使用，以限制过度更新值。修改的更新可以是用于后续处理的更新值。

在框304-308中，简化的附加考虑可以有助于锥面波束几何重构(例如，如在血管造影术中所使用的)，其中，投影的数目p可以较小(p＜＜N，其中，N可以表示沿着穿过重构体积的典型投影路径的体素的数目)。在这种情况下，各个单独“蜘蛛腿”可以是截然不同的，尤其在多个重构网格的高分辨率最终阶段。针对该高分辨率n-D空间中的许多小区域计算脉冲响应R_p是耗时的。作为使用n-D脉冲响应的近似，可以使用维度m的p个低维度投影的脉冲响应(0＜＝i＜p)。然后它们对应修改的逆或卡尔曼/维纳增益可以迭代地应用于投影的各个单独误差图像。这些m-D脉冲响应可以是上述n-D脉冲响应到相应m-D投影空间的投影(可以是蜘蛛状)。它们的属性如上所述，即，在投影空间中缓慢变化。然而，用于计算p个投影的脉冲响应、这些计算的脉冲响应的逆、以及后续卷积(目前应用于m-D投影)的计算负荷远小于在n-D空间中计算它们，尤其是在p＜＜N时。计算的脉冲响应的修改方式可以类似于R。也可以对计算的逆脉冲响应进行带通滤波，以降低层析摄影重构中的噪声。这种用于层析摄影逆运算的技术可以类似于滤波的反向投影，在滤波的反向投影中在重构之前对投影本身进行滤波。该技术不限于锥面波束几何形状，但是也可以应用于其他几何形状，例如，平行波束几何形状。

在框309中，可以将修改的更新值i_g与n-D体素集合中对应图像密度的值相加，以修正n-D体素集合v。相加可以由与上述方程(1)和(2)类似的极限函数来限制。可以针对整个n-D体素集合的一部分执行相加。为了精确，可以对体素值v进行缩放，使得它们的投影最佳地匹配未处理图像数据。

在框310中，如果分辨率网格是最细化的，并且如果计算的预测与输入投影数据之间的差(例如，残差)足够小(在统计上)，则然后可以向输出设备输出表示密度估计v、v_L或二者(并且，根据需要，估计的预测和误差)的修正的n-D体素集合，如框311所示。可以在输出之前应用平滑。否则，可以将分辨率网格变为尺寸大于当前分辨率网格的分辨率网格(或者可能地，可以降低n-D平滑核的重要性/范围)，并且在框312中，流可以前进至框303。

在框311中终止时，可以获得对象106的高保真重构图像。该重构图像106可以非侵入且非破坏性地揭示对象203的内部信息，在诊断医学成像和非破坏性评估方面具有广泛应用。重构图像106可以在可视化设备上可视，或者存储在计算机存储器或存储介质中，如上针对输出设备107所讨论的。

根据Trachtenberg，S.，Dorward，L.M.，Speransky，V.V.，Jaffe，H.，Andrews，S.B.，Leapman，R.D.，″Structure of the Cytoskeleton ofSpiroplasma melliferum BC3 and Its Interactions with the CellMembrane，″J.Mol.Biol.，378，pp.776-787，2008的投影像素数据104可以用于展示经由本发明示例性实施例可得到的改进。图4A示出了基于对投影像素数据104的滤波反向投影(FBP)重构的重构图像106。图4B示出了根据本发明示例性实施例的重构图像106。图4A和图4B之间的比较说明了：通过各种密度(例如，细胞间与细胞外空间、细胞膜、核糖体、脊髓液细节等)的强烈对比的测量结果，图4B中的图像清晰度更好。用户更期望图4B中的图像。两种重构在近似+/-70度的范围上均使用近似166个双倾斜发射电子显微术(TEM)投影。

图5示出了仿真的人体模型的正面样本投影。仿真的人体模型在尺寸为256×256×64个体素的3D盒体中的随机位置处包含大约100个尺寸为6×6×6个体素的立方体。体素具有0.85个单位周围密度中所嵌入的5个单位的密度。添加了白噪声的总共49个投影(单倾斜，+/-72度)可以由计算机仿真103来获得，以产生投影像素数据104。

根据统计模型验证中公知的黄金标准，通过检验残差可以最严格地评价重构质量，残差即是真实图像(例如，图5中仿真的人体模型，被称作先验)与重构图像106之间的差异。高保真度重构可以获得随机并且包含极少来自于3D对象结构(例如，图5中仿真的人体模型中的立方体)的贡献的残差。图6A示出了同时迭代重构技术(SIRT)的残差，而图6B示出了基于根据本发明示例性实施例的方法获得那些残差。图6A中相当醒目的残差方形图案(在15次迭代之后)和图6B中出现的显著白噪声(在一个高分辨率迭代之后)还说明本发明在图像保真度和计算效率方面的优势。针对该示例可以优化SIRT的迭代次数，因为针对SIRT更多或更少次的迭代导致图像质量劣化。

本发明还可以应用于投影数目较少的情况，例如，血管造影术中的锥面波束几何形状。图7A示出了在塑料管制成的跳动液力冠状人体模型中注入碘对比染料期间(例如，一般在心导管血管造影术中采用)的一个X射线投影。由于注射的瞬态属性，所使用的投影像素数据104来自于下面人体模型的仅6数目字减影投影。图7B示出了经由本发明示例性实施例的重构图像106，精确地表示跳动液力人体模型的几何形状。可以通过滤波反向投影(FBP)或同时迭代重构技术(SIRT)或者甚至最大熵(MENT)方法来获得重构图像，这是因为投影数目通常较少，和/或利用当前计算资源迭代次数变得不实际。

如上所述，本发明的一些示例性实施例组合了计算/处理技术，并且获得对感兴趣区域中期望的n-D体素的快速且精确的计算/估计，而同时允许利用任意m-D投影(和方向)。应当注意，m可以是固定的，并且依赖于投影索引。技术组合的效率可以应用于广泛的投影数目范围。可用投影(测量的)信息的高效使用可以使该方法在投影数目非常大的情况以及投影数目有限(较少或不完整)(例如，在冠状血管造影术)的情况下有效率。

本发明非线性逆运算技术可以组合：(i)通过将n-D图像估计分解到重构域内的两个非线性相关域中而将非线性卡尔曼滤波器的特定部分线性化(近似贝叶斯估计)；(ii)使用减小(与“台阶”可比)的图像平滑程度进一步细化分分辨率的多网格计算，与“阶梯情况”可比；以及(iii)根据需要，通过使用逆赫赛估计图像噪声所修改的线性系统的逆脉冲响应(或传递函数)的高效计算，获得用于对象空间中选定区域的近似卡尔曼/维纳滤波器增益。

具体地，在Sage，Jarisch W.R.& Detwiler J.S.，″StatisticalModeling of Fetal Heart Rate Variability，″IEEE Trans.On Biomed.Eng.，Vol.BME-27，No.10，pp.582-589，Oct.1980中已经讨论了非线性卡尔曼滤波器，并且，在Jarisch W.，Determination of Material Properties byLimited Scan X-ray Tomography，Technical Report USAF AFWAL-TR-81-4050，NTIS，HC A07/MF A01，Defense Technical Information Center，Cameron Station，Alexandria，Virginia 22314，1981中将卡尔曼滤波器应用于图像重构。可以代替逆运算的线性部分使用备选方法(例如，FBP、FFT和维纳滤波器)。重构区域内两个非线性相关域中的任一个可以表示卡尔曼系统状态。

通过上述组合，可以利用每次迭代实现有效率的误差减少，类似于高效地使用优化的赫赛。所实现的效率可以接近四次收敛的效率。该方法可以绕过其他方法中遇到的高维度和不良收敛速率的严重问题，例如，在Kolehmainen，V.，Brandt，S.S.，Structure-From-MotionWithout Correspondence From Tomographic Projections by BayesianInversion Theory，IEEE Trans.Medical Imaging，ITMID4，Vol.26，No.2，pp.238-248，Feb.2007中的那些其他方法。

备选地，本发明的示例性实施例可以视为IR(空间可变/不稳定/分段的)与FBP的组合，在本文被称作S-FBP。根据近似的需要，存在与S-FBP相关联的一个或多个区段。通常，FBP在重构的场上不稳定。然而，场属性的缓慢变化可以通过将重构场拆分成交叠区域/区段来实现近似，这些交叠区域/区段在这些区域之间具有加权过渡。S-FBP的速度、效率和灵活性可以通过若干实现方式的适合选择和组合来实现，如下所述。

例如，S-FBP部分可以应用于非线性变换数据，而不是未处理投影数据。

例如，多网格计算可以支持线性化，并且使每次迭代的运算次数保持极低。例如，对最细化分辨率网格的最终计算之前总运算次数可以低于S-FBP对最细化分辨率网格的运算次数。

例如，体素密度可以在两个非线性相关域(经由映射连接)中表示：一个域用于经由S-FBP(可以是n维空间中的变量，依赖于期望的近似技术)更新；并且另一个域用于计算投影预测。S-FBP的滤波器权重可以依赖于体素位置，但是与相同权重一起使用的成组体素可以提高计算效率。

例如，可以根据线性化系统脉冲响应(或传递函数)的逆运算计算S-FBP部分的滤波器权重，来使用任意投影维度和方向。线性化系统脉冲响应(或传递函数)可以与根据对输入投影数据的测量而确定的噪声偏差相组合。

例如，适当选择的极限(微调)函数可以防止计算期间的数值奇异点和不稳定性。

例如，在每次迭代时恰到好处地应用于体素和像素值的适当选择的平滑函数可以有利于数值计算的稳定性和精度。体素和像素值的这种平滑在功能上比得上产生更新值之前的卡尔曼预测符算法。平滑的迭代应用可以类似于卡尔曼滤波器中的状态预测符算法。

由于在最细化分辨率网格上计算之前的迭代的精度是足够的，因此在整个最细化分辨率网格上的迭代在计算上是不必要的，因此可以免除误差校正项的计算。在这些情况下，重构速度可以变得比得上S-FBP的重构速度。这是因为在最细化分辨率网格上的计算可以实现比得上FBP运算次数(例如，在因子＞1内)，而针对较粗分辨率网格上的所有在前较低分辨率迭代的组合运算次数可以少于对最细化分辨率网格的计算的运算次数，尤其对于较高维度图像重构。

对于各种原因，如上所述，非线性更新技术可以不执行利用单次迭代的理想校正。取而代之，在连续迭代中，存在更新太小并且残差在若干次迭代上持续的区域，或者存在更新太大，并且残差在若干次迭代上持续呈现有限次循环(例如，在正负号之间振荡)的区域，即，对于具有足够高增益的反馈回路是公知效果。有利地，来自在前迭代的残差可以用于实现更高效的更新。

具体地，可以根据在前校正/残差的历史r_h来评价每个图像点处的增益(对于投影空间中给定像素的投影增益，或者对象空间中体素的体素增益)。例如，本发明的校正/残差可以表示为r₀(两个维度中给定位置ij处，或者三个维度中给定位置ijk)，并且针对在前迭代的残差表示为r₁。位置处的增益因子(例如，当m＝2时在位置ij处)通常在r₀和r₁具有相同符号时增大；类似的，增益在r₀和r₁具有相反符号时减小(并且不能被其他因子驱动)。能够影响增益修改的选择的其他因子可以包括例如，与校正r_h相关联的不确定性、来自于其他投影的估计驱动力的不确定性(例如，当应用于来自其他投影的残差的增益因子主导重构的校正时)、在前增益的有限精度、投影图像的噪声级别、和/或稳定解的估计收敛程度。

因此，可以根据在前迭代中的历史图案来自适应地调节局部投影和体素增益。图8A示出了根据本发明示例性实施例的扩展流程图。除了附加框以外，图8A的扩展流程图类似于图3。

在框801中，例如可以将编码了预测投影数据与输入投影数据之间的差的第一数据记录(作为投影残差)在数据残差设备中，例如残差设备102或其他设备。该差可以从框305获得。

同样在框801中，可以利用增益矩阵G_t对迭代t下的第一量r_t(来自框305)进行加权，以基于记录的来自在前迭代的第一数据获得修正的第一量。

可以在给定低(例如，4×4或更低)网格分辨率下针对所有像素i选择与单个固定值相关联的适合初始增益矩阵G₀(G_0i＝G₀＞0)。然后通过针对第t次迭代将矩阵g_t中的增益改变因子g_t(例如，通过在对应位置处的逐元素相乘)应用于增益矩阵G_t的对应区域，来增大或减小增益矩阵G_t的初始固定元素值G_t。增益改变因子g_t可以依赖于观察到的更新/递归的稳定性/振荡，如下所述。

谨慎选择的初始值G₀可以足够小，以安全地避免估计的对象密度的振荡，并且可以足够大，以实现估计的对象密度高效校正。一些数值尝试可以获得给定重构设置的适合初始值。如下所述，对最优值的选择在这点处可以不重要，因为可以在后续迭代中寻求最优值。

当前两次迭代之后获得了两个或更多连续的投影残差集合，针对每个单独投影像素的增益矩阵G_t的调节开始。

调节增益矩阵G_t的目标可以是避免密度估计的振荡(例如，交替地)，太高增益的特性。处理备选投影残差的可能示例调节(或等同地，如与框803相关联地振荡对象密度)在方程15中示出：

$g_t\≈\frac{r_s}{r_s-r_t}$ 方程(15)

其中，g_t是估计的增益改变因子，并且r_s和r_t表示给定像素位置处来自连续迭代s和t(s＝t-1)的残差。

方程15是增益改变因子g_t的一个示例。然而，方程15不能直接使用，因为残差r_s和r_t中的噪声可以典型地引起从增益改变g_t的序列获得的增益矩阵G_t的过度波动。

为了防止增益矩阵G_t的元素具有例如除以r_s-r_t的差值的近零值而产生的过多值，可以实现若干预防措施，例如：

(i)例如通过估计投影残差的局部噪声功率(例如，由sig_ref表示的参考区域的标准偏差所测量的)来获得残差的噪声的局部测量。

(ii)当r_s和r_t符号相反时，但是|r_s|和|r_t|均在以下sig_ref时，或者当r_s和r_t符号相反并且|r_s-r_t|＜2·sig_ref时，可以将g_t设置为1。

(iii)在不应用(ii)的情况下，当方程15的分母中r_s或r_t的绝对值降到sig_ref之下时，r_s或r_t可以由sign(r_x)·sig_ref代替，根据需要使|1-g_t|最小化，其中，x∈{r，t}。

(iv)当r_t的幅度持续大于r_s，而r_t和r_s具有相同符号，在特定位置处会出现控制的潜在丢失。局部增益改变因子g_t的增大可以缓和控制的丢失。在该位置处对局部增益改变(增大)的可能选择可以是：

$g_t=\sqrt{\frac{{\mathrm{sig}}_{\mathrm{ref}}+f\·|r_s-r_t|}{|r_q-r_s|+{\mathrm{sig}}_{\mathrm{ref}}}}$ 方程(16)

其中，s＝t-1，q＝s-1，并且

$f=\frac{{\mathrm{sig}}_{\mathrm{ref}}^2}{{\mathrm{sig}}_{\mathrm{ref}}^2+0.1\·r_s\·r_i}$ 方程(17)

(v)通过设置最大和最小极限值来限制局部增益改变因子g_t的范围：

g_min＜g_t＜g_max          方程(18)

其中，g_max＝C·n_pro     方程(19)

g_min＝1/g_max            方程(20)

并且，其中n_pro是使用的投影数目。C的适合值典型地可以在1以下，例如，0.25。可以以类似方式约束G_t的元素值，例如，对于投影的所有位置i和j，max{G_ti}/min{G_tj}＜0.5*n_pro，其中，n_pro是可用投影的数目。

(vi)可以通过在小邻域上对局部增益改变因子进行平均(平滑)来减小局部增益改变因子g_t的噪声。对于2-D投影数据，3×3平均阵列可以倾向于显著提高迭代重构过程的收敛。

在框801中，当记录的第一数据接续在前迭代具有相同符号时，如上所述的修正的第一量可以倾向于增大投影增益，并且当记录的第一数据接续在前迭代具有相反符号时修正的第一量可以倾向于减小投影增益。

因此，描述了依赖于来自在前迭代的投影残差的自适应增益矩阵。自适应增益矩阵还可以依赖于像素在投影空间中的位置。

在框802中，可以使用P、来自方框307的逆函数、以及来自框801的修正第一量或来自框305的第一量，按照与关于框308讨论的类似方式计算更新值。

在框803中，可以将编码了更新值的第二数据记录在数据存储设备中，例如，存储设备102或其他设备。可以基于来自在前迭代的记录的第二数据修正更新值。可以根据在前迭代，基于更新值序列u_t以及可能的投影残差的子集，按照与以上关于方框801描述的类似方式来计算用于自适应调谐更新值的增益矩阵G_t。与方框801相比，对应更新值应用增益可以通过更直接地减小估计误差，自适应地估计对象密度估计。在方框803中，更新值可以与体素增益相关联。当记录的第二数据接续在前迭代具有相同符号时，修正的更新值可以倾向于增大体素增益，并且当记录的第二数据接续在前迭代具有相反符号时修正的第二量可以倾向于减小体素增益。

在框804中，来自框803的修正的更新值或者来自框803的更新值可以与n-D体素集合v中的对应对象密度的值相加，以修正n-D体素集合v。通过类似于上述方程(1)和(2)的极限函数来限制相加。此外，可以将估计的对象密度约束到特定先验已知值，例如，零/可忽略的密度。此外，可以针对整个n-D体素集合的一部分来执行相加。为了精确，可以对体素值v进行缩放，使得它们的投影最佳地匹配未处理投影数据。

因此，还描述了依赖于来自在前迭代的对象估计中的误差的自适应增益矩阵。自适应增益矩阵可以考虑重构的几何因子以及数据相关因子。例如，几何因子可以是已知的空区域(零密度)，对于这些空区域不需要增益，并且将估计的对象密度重设为零。

图9A至9D分别示出了10、50、130和170度轴倾角处的示例输入投影数据。可以从图8A的框302获得输入投影数据。对于来自图8A中框303到框312的每个迭代循环，可以针对每个投影交叠获得预测投影矩阵。图9A至图9D中示出的输入投影数据来自于发射电子显微术/扫描发射电子显微术(TEM/STEM)，并且通过Karlsruhe技术学院(KIT)以一度投影角度递增地提供。

图10A至10D分别示出了根据图8A的框304获得的针对迭代1、2、5和6的预测投影矩阵。大约140个投影的集合用于重构。预测的投影数据对应于选定的10度单轴倾角。

图11A至11D分别针对迭代1、2、5和6示出了10度倾角下输入投影数据与预测投影数据之间的示例残差(在最高分辨率网格处)。

图12A至12D分别针对迭代1、2、5和6示出了在10度轴倾角处从框801中获得的生成的投影增益矩阵G_t。具有持续残差值的区域可以利用增加的投影增益(例如由一些中心珠增加的亮度所指示)来处理，而呈现有限次循环的其他区域(即，在连续在前迭代上正负符号交替的振荡)利用减小的增益来处理。显示为对数刻度。每幅图像左侧上参考刻度表示对于每个灰度级变化10％。典型的动态范围增益(即，像素或体素的高低比)在因子10与100之间，暗示着对应的计算加速。

上述用于修正第一量和更新值的两个自适应增益矩阵在来自图8A的框303至框310的反馈回路中，以计算对象密度估计。传统上，可以将选择“良好”校正增益因子应用于迭代的误差测量(即，残差)，能够与整形矩阵相结合作为例如在Pan，T.-S.，Yagle，A.E.，Accelerationand Filtering in the Generalized Landweber Iteration Using a VariableShaping Matrix，IEEE Trans.Med.Imaging，ITMID-4，Vol.12，pp.278-292，June 1993中讨论的用于有效率计算的重要设计参数。过于小的校正增益因子会导致向着解的收敛过慢，而过大的校正增益因子会导致迭代的限制次循环或振荡，并且可能导致解的发散。因此对“良好”校正增益因子的恰当好处的选择由此显著提高了迭代算法的性能。

如本发明人所发现的，利用适中的高校正增益因子的残差分析揭示了有限次的循环/振荡典型地可以受限于重构图像(或对应的投影图像)的残差的特定区域。此外，可以在对应的重构图像区域中找到投影的有限次循环/振荡，并且暗示着重构对象密度的特定部分经历有限次循环或振荡。基于本发明人的这种新颖观察，可以在投影空间和/或对象空间中使用缓和区域增益的空间变化，使得观察到有限次循环/振荡。可以减小投影增益或体素增益，并且在收敛较慢的情况下，可以减小投影增益或体素增益。

此外，增益矩阵的空间结构可以用于优化增益值。如图9A至9D以及12A至12D所见，增益值的图案倾向于遵照输入投影数据的图案。为了利用这些系统局部关系，并且减少噪声增益矩阵中不确定性的量，可以使用智能平滑来获得更精确的增益估计。例如，当使用平面投影(维度m＝2)时，可以通过利用回归模型在投影空间或对象空间中拟合元素周围适合的(通常较小的)邻域，来平滑估计的增益矩阵g_t。用于拟合(fitting)覆盖感兴趣区域的若干像素g_ti的可能回归模型可以是：

$g_{\mathrm{ti}}={\mathrm{\Σ}}_{k=0}^K{\mathrm{\α}}_k{p}_{\mathrm{ti}}^k+{\mathrm{\ϵ}}_{\mathrm{ti}}$ 方程(18)

其中，t表示迭代t，i可以表示来自感兴趣区域的小元素集合，α_K表示K个回归系数的集合中的元素，p可以表示投影像素值，ε表示不确定性。例如，K值为2可以足以用于实现足够速度的算法。此外，估计g_ti(例如，最小平方估计)可以用于获得局部估计g_ti ^reg。更平滑的g_ti ^reg可以用于修正框801和图11A至11D中示出的m维度第一量。可以对与对象密度估计的n维度更新(在框803中示出)相关联的增益改变因子应用类似方法。根据该方法，可以增大符号持续更新的区域中的增益，并且减小有限次循环/正负号交替的那些区域中增益。如框303中所示，可以针对m维投影p_i...j或计算的n维对象密度p^obj _i...k或参考值的其他组合来实现平滑过程。

最后，上述过程直接实现会集中需要存储器，这是因为需要保留来自框801中框802中的在前迭代的第一数据和/或第二数据。为了降低改进估计的存储器需要，可以使用缩减的精确实现方式。例如，残差和校正可以仅是未处理图像数据或估计的对象密度的幅度的一部分。因此，与原始输入投影数据或对象密度表示相比，以较低精度按照压缩格式数字地免除第一数据或第二数据以及来自框305的第一量或来自框306的第二量。此外，中等精度(例如，几个比特)的增益因子估计通常是足够的。总之，精度降低的数据表示可以将存储器需要保持在中等水平，并且在计算效率上几乎没有不利影响。

图8B示出了根据本发明另一示例性实施例的另一扩展流程图。除了框306和307以外，图8B的扩展流程图类似于图8A，并且图8A的802被图8B的框805代替。在框805中，可以经由适合的层析摄影重构算法(例如，使用公知的中央分片理论(CST)的算法、滤波反向投影算法等)，基于来自框802的修正的第一量计算更新值。对于高质量重构，尤其在使用锥面波束几何形状时，针对感兴趣区域的典型体素的局部投影几何形状(例如，由交叉投影波束的方向给出)可以自计算第一量时应用，并且用作至所选层析摄影重构算法的输入。在锥面波束重构的情况下，期望选择典型体素集合。对于针对典型体素集合而执行的第一量的全体积层析摄影重构中的每一个，可以通过所述交叉投影波束的特性(例如如上关于框306描述的)变得显著的距离划分典型体素的间隔因此，针对围绕典型体素的各个单独区域的输出可以通过适当地加权它们的输出来接合。针对与特定重构区域相对应的输出的权重可以主导所有其他权重，例如，在参考体素处逼近1。然后接合的重构输出可以由增益矩阵G_t’来修正，并且用于修正n-D体素集合。尽管在框805的输入和输出处可以分别指定不同的增益矩阵G_t和G_t’，但是出于多种原因单个增益矩阵是有利的。例如，使用多个矩阵可以增加通过框3 12回到框303的反馈回路的计算不稳定的风险。

关于图3、8A和8B所示的方法讨论的本发明的示例性实施例可以作为数据存储设备(例如，CD-ROM、DVD、磁光(MO)盘、硬盘、软盘、压缩盘、闪速盘等)上存储的软件代码来实现。存储的软件代码是由具有一个或多个处理器、一个或多个处理设备(例如，随机存取存储器(RAM)设备、动态RAM(DRAM)设备、闪速存储设备以及静态RAM(SRAM)设备等)的计算机可读并且可执行的，来执行以上关于图3、8A和8B讨论的示例性方法。

本发明的示例性实施例可以提供一个或多个程序存储和执行设备(例如专用集成电路(ASIC)和现场可编程逻辑门阵列(FPGA)、复杂可编程逻辑器件(CPLD)、专用指令集处理器(ASIP)等)，以存储和执行以上关于图3、8A和8B讨论的示例性方法。

本文描述的示例和实施例是非限制示例。

关于示例性实施例描述了本发明，并且根据上述本发明对于本领域技术人员显而易见的是，在不背离本发明的更广泛方面可以进行改变和修改，因此权利要求中限定的本发明应覆盖落在本发明真实精神内的所有这样的改变和修改。

Claims (20)

1.一种系统，包括：

一个或多个发射器，将激发能量发射到被观察对象中；

一个或多个检测器，响应于发射到被观察对象中的激发能量，产生投影空间数据，投影空间数据编码了所述一个或多个检测器接收的能量；

控制器，控制所述一个或多个发射器发射激发能量，并且控制所述一个或多个接收器产生投影空间数据；以及

具有至少一个处理器的图像重构器，接收投影空间数据，并且通过以下操作处理投影空间数据：

计算第一量，所述第一量表征投影空间数据与预测的投影数据之间的差，其中投影空间数据与预测的投影数据之间的差与投影增益相关联；

将编码了所述差的第一数据记录在数据存储设备中；

基于来自在前迭代的所记录的第一数据，计算修正的第一量；

使用修正的第一量来计算更新值；并且

使用更新值重构表示被观察对象的对象空间图像。

2.根据权利要求1所述的系统，其中，激发能量是以下中的至少一个：电磁(EM)能量、X射线能量、粒子束、红外能量、光能量、和振动能量。

3.根据权利要求1所述的系统，其中，上述图像重构器还通过以下操作来处理投影空间数据：

在数据存储设备中记录编码了所述更新值的第二数据；

基于来自在前迭代的所记录的第二数据，修正所述更新值，其中使用修正的更新值来重构对象空间图像。

4.根据权利要求1所述的系统，还包括：

输出设备，接收对被观察对象加以表示的对象空间图像，其中，输出设备是数据存储设备、显示设备和打印设备中的至少一个。

5.根据权利要求1所述的系统，还包括：

研究者计算机，接收对象空间图像，并且被编程为执行从对象空间图像提取诊断信息，或者针对所述一个或多个发射器、所述一个或多个检测和图像重构器中的至少一个细化调谐参数。

6.一种工作站，包括：

一个或多个处理器；以及

一个或多个数据存储设备，存储要由所述一个或多个处理器执行的软件，所述软件包括：

接收输入数据的软件代码，输入数据表示从对象空间投影到投影空间的对象；

计算第一量的软件代码，第一量表征输入数据与预测的投影数据之间的差，其中输入数据与预测的投影数据之间的差与投影增益相关联；

使用第一量计算更新值的软件代码；

在所述一个或多个数据存储设备或另外一个或多个数据存储设备中记录编码了所述更新值的第二数据的软件代码；

基于来自在前迭代的所记录的第二数据修正所述更新值的软件代码；以及

使用修正的更新值重构表示被观察对象的对象空间图像的软件代码，其中修正的更新值与体素增益相关联。

7.根据权利要求6所述的工作站，其中，软件还包括：

在所述一个或多个数据存储设备或另外一个或多个数据存储设备中记录编码了所述差的第一数据的软件代码；以及

基于来自在前迭代的所记录的第一数据修正第一量的软件代码。

8.根据权利要求6所述的工作站，其中，

利用第一压缩技术存储所述第一量或所述第二量，利用第二压缩技术存储所述第一数据或所述第二数据，并且所述第二压缩技术与第一压缩技术相比具有可比的精度。

9.根据权利要求6所述的工作站，其中，输入数据是从以下中的至少一个中接收的：投影图像获取系统、所述一个或多个数据存储设备、另外一个或多个数据存储设备、所述一个或多个处理器执行的软件仿真、和另外一个或多个处理器执行的软件仿真。

10.根据权利要求6所述的工作站，其中，输入数据是经由数据总线、电缆、有线网络或无线网络接收的。

11.一种由执行一个或多个数据存储设备上存储的软件代码的一个或多个处理器来实现的方法，该方法包括：

从以下中的至少一个接收输入数据：投影图像获取系统、一个或多个数据存储设备、另外一个或多个数据存储设备、所述一个或多个处理器执行的软件仿真、另外一个或多个处理器执行的软件仿真，输入数据表示从对象空间投影到投影空间的对象；

计算第一量，第一量表征分辨率网格上的输入数据与分辨率网格上的预测的投影数据之间的差，其中，输入数据与预测的投影数据之间的差与投影增益相关联；

在一个或多个数据存储设备或另外一个或多个数据存储设备中记录编码了所述差的第一数据；

基于来自在前迭代的所记录的第一数据，计算修正的第一量；

使用修正的第一量计算更新值，其中更新值与体素增益相关联；

在所述一个或多个数据存储设备或另外一个或多个数据存储设备中记录编码了所述更新值的第二数据；

基于来自在前迭代的所记录的第二数据，修正所述更新值；

使用修正的更新值重构表示被观察对象的图像；并且

向输出设备输出重构图像，其中输出设备是以下中的至少一个：所述一个或多个数据存储设备、另外一个或多个数据存储设备、显示设备、和打印设备；

其中，根据所述一个或多个数据存储设备上存储的软件代码，由所述一个或多个处理器执行接收输入数据、计算第一量、记录第一数据、计算修正的第一量、计算第二量、重构图像和输出重构图像。

12.根据权利要求11所述的方法，其中，

当所述记录的第一数据接续在前迭代具有相同符号时，所述修正的第一量倾向于增大投影增益；

当所述记录的第一数据接续在前迭代具有相反符号时，所述修正的第一量倾向于减小投影增益；

当所述记录的第二数据接续在前迭代具有相同符号时，所述修正的更新值倾向于增大体素增益；并且

当所述记录的第二数据接续在前迭代具有相反符号时，所述修正的更新值倾向于减小体素增益。

13.根据权利要求12所述的方法，其中，利用第一压缩技术存储所述第一量或所述第二量，利用第二压缩技术存储所述第一数据或所述第二数据，并且所述第二压缩技术与第一压缩技术相比具有可比的精度。

14.根据权利要求13所述的方法，还包括：

基于空间模型，针对所述更新值平滑所述第一量或所述修正的第一量，其中，所述空间模型合并了以下中的至少一个：预测的投影数据、表示对象的重构图像、和编码了投影空间或对象空间的噪声等级或噪声容限的信息。

15.根据权利要求11所述的方法，还包括：

计算与至少一个脉冲响应相对应的第二量，每个脉冲响应对应于对象空间中位置处的单位强度的参考体素，其中所述计算更新值还包括使用第二量的逆函数。

16.根据权利要求15所述的方法，其中，针对覆盖对象空间中单位强度的参考体素的位置的周围区域，计算逆函数。

17.根据权利要求15所述的方法，其中，重构图像包括：

将第一量与第二量的逆函数卷积，以获得更新值。

18.根据权利要求11所述的方法，其中，使用层析摄影重构算法来计算更新值。

19.根据权利要求11所述的方法，还包括：

使用更新值修正对象空间体素集合。

20.根据权利要求19所述的方法，还包括：

在修正对象空间体素集合之后，变化到尺寸比所述分辨率网格大的不同分辨率网格。

Images