高效计算机断层摄影技术的制作方法

专利名称：高效计算机断层摄影技术的制作方法
技术领域：
本发明涉及图像重建领域。
背景技术：
基于η维空间中的区域的m维投影计算(估计)该区域的密度(以像素表示；通 常0<m<n)的断层摄影技术分成两大类滤波反投影(FBP)和迭代重建(IR)，其中，迭代 重建(IR)也被认为是代数重建技术(ART)的变型。FBP和常规的顶技术两者都具有特定 的缺陷，特别是在现代设置中，当一个投影无论在任意处通常包含2000个像素(一维)直 到2000 X 2000个像素(二维)时，通常需要计算IO6 101°个图像元素(即体素)。与FBP相关的主要缺陷包括需要大量的投影以实现有限的量化精度。投影的数量 通常数以百计，但是这些投影没有如其应该的那样被有效地使用。例如，计算的密度估值可 能是负值，然而已知物理强度不可能为负值。有时通过后续的迭代细化可以提高有限的量 化精度。如 Wood，S. L.禾口 Morf，M.在 A fast implementation of a minimum variance estimator for computerized tomography image reconstruction (IEEE, Trans.on Biomed. Eng. ,Vol. BME_28，No. 2,pp. 56-68，1981)中(下文称伍德)所述，与 FBP 相关的其 它缺陷包括不能根据体素位置改变数据权重以及不能有效地考虑各种限制。由于FBP方法通常依赖于快速傅里叶变换、中心切片定理以及可能使用合适且预 先计算的滤波器，因此，FBP方法的优点包括其计算速度。FBP和常规的顶技术的替代技术使用矩阵运算且可能适用于体积或图像元素(体 素或像素)为几千个的“小”问题。然而，对于示例性的断层摄影重建设置，通过这些基于 矩阵的技术在可预见的将来不能解决计算负担。这是由于计算机运算的数量按比例为N的 幂。例如，如伍德和2001年12月18日的美国专利6332035所述，在三维的情况中，运算计 数可以按比例快于N5 (其中N是表示感兴趣区域的重建立方体的一侧的像素数)。IR技术的好处包括其解释限制能力，尤其是保证密度估值不为负值的能力。保证 非负的密度估值可以导致显著的图像对比度增加和更好的量化结果。与常规的顶技术相关的缺陷包括需要重复求解反演运算以及前向投影运算以得 到后续迭代的校正项。所伴随的对于这些运算的许多递归对的需求通常使得顶技术变慢， 并且顶技术依据特定的实施方式变得不稳定。由于这些常规顶技术使用各种(迭代)优化方法，因此，其到解答的收敛速度取 决于特定的技术占作为待优化的目标函数的二阶导数矩阵(即，相对于图像元素的性能标 准)的海赛矩阵(Hessian matrix)的比率。最常用的技术是顶技术的变型。然而，由于用 于图像重建的海赛矩阵的大尺寸，因而，其结构(以及其反演结构)通常被忽略或被拙劣地 近似。此外，由于海赛矩阵的特征值的宽分布，目前的优化技术在超出一定数量的迭代(通 常数以几十到几百计)后往往显示不出改进且可能仅处理大的特征值。这些算法的多网格变型可能有帮助，但由于涉及精细网格的海赛矩阵的尺寸最终 还是失败了。多网格分辨率在这里是指执行迭代时使用逐步精细的分辨率。

发明内容
本发明可以提供一种系统，其包括一个或多个发射器，用于将激发能射入观察的 对象；一个或多个检测器，用于响应于射入观察的对象的激发能来生成对于由所述一个或 多个检测器所接收的能量进行编码的投影空间数据；控制器，用于控制所述一个或多个发 射器发送所述激发能，并控制一个或多个接收器生成所述投影空间数据；以及图像重建器， 用于接收所述投影空间数据并如下处理所述投影空间数据计算表征所述投影空间数据与 预测投影数据之间的差的第一量；计算对应于至少一个脉冲响应的第二量，每个脉冲响应 对应于单位强度的参考体素；使用所述第二量的逆函数和所述第一量计算更新值；以及使 用所述更新值重建表示所述观察的对象的对象空间图像。本发明可以提供一种工作站，其包括一个或多个处理器；以及一个或多个数据 存储装置，用于存储由所述一个或多个处理器执行的软件，所述软件包括用于接收表示从 对象空间投影到投影空间的对象的输入数据的软件代码；用于计算表征所述输入数据与预 测投影数据之间的差的第一量的软件代码；用于计算对应于至少一个脉冲响应的第二量的 软件代码，每个脉冲响应对应于所述对象空间中的位置处的单位强度的参考体素；用于使 用所述第二量的逆函数和所述第一量计算更新值的软件代码；用于使用所述更新值重建表 示所述观察的对象的对象空间图像的软件代码。本发明可以提供一种方法，其由一个或多个处理器执行存储在一个或多个数据存 储装置上的软件代码来实施，所述方法包括从投影图像获取系统、所述一个或多个数据存 储装置、另外一个或多个数据存储装置、由所述一个或多个处理器执行的软件仿真或由另 外一个或多个处理器执行的软件仿真中的至少一个，接收表示从对象空间投影到投影空间 的对象的输入数据；计算表征分辨率网格上的所述输入数据与所述分辨率网格上的预测投 影数据之间的差的第一量；计算对应于至少一个脉冲响应的第二量，每个脉冲响应对应于 所述对象空间中的位置处的单位强度的参考体素；使用所述第二量的逆函数和所述第一量 计算更新值；使用所述更新值重建表示所述观察的对象的图像；以及将重建的图像输出到 输出装置，其中，所述输出装置是所述一个或多个数据存储装置、另外一个或多个数据存储 装置、显示装置或打印装置中的至少一个；其中，根据存储在所述一个或多个数据存储装置 上的所述软件代码，所述一个或多个处理器执行接收所述输入数据、计算所述第一量、计 算所述第二量、计算所述更新值、重建所述图像、和输出重建的图像。


现在将结合附图描述示例性的实施例。图1示出根据本发明的示例性的实施例的系统图。图2示出根据本发明的示例性的实施例的成像系统。图3示出根据本发明的示例性的实施例的流程图。图4A和4B分别示出基于常规的滤波反投影(FBP)和根据本发明的示例性的实施 例的方法重建的细胞图像。图5示出仿真体模的正面样本投影。图6A和6B分别示出基于同步迭代重建技术(SIRT)和根据本发明的示例性的实施例的方法的重建误差。图7A示出在将碘造影剂注入搏动的液压心脏体模期间的一个X射线投影。图7B示出基于根据本发明的示例性的实施例的方法获得的重建的体模图像。
具体实施例方式下面详细论述本发明的示例性的实施例。在描述示例性的实施例时，为清楚起见 使用了特定的术语。然而，本发明不意在限制于如此选择的特定术语上。相关技术领域的 技术人员将认识到在不偏离本发明广义概念的前提下可以使用其它等同组件以及开发其 它方法。并入了此处引用的所有参考文献，如同各自并入每个参考文献。图1示出根据本发明的一些实施例的系统图。使用图像重建器105可以将包含对 对象的内部结构进行编码的信息的投影像素数据104变换成重建的图像106，以在输出装 置107上显现。投影像素数据104可以来自实验所得101、或计算机上的仿真103、或包含 所记录的来自先前实验或仿真的投影像素数据104的数据存储装置102。实验所得101、数 据存储装置102和仿真103可以直接或经由例如局域网(LAN)、广域网(WAN)、因特网等的 网络远程地将投影像素数据104提供给图像重建器105。尽管此处仅示出了投影像素数据 104，通常，该数据可以是将任何形式的激发能传递到观察的对象中的结果，因此除来自电 脑断层摄影技术(CT)的数据以外，可以包括来自例如电子显微镜、磁共振(MR)、正电子发 射断层摄影技术(PET)、单正电子发射计算机断层摄影技术(SPECT)、超声波、荧光、多光子 显微镜(MPM)、光学相干断层摄影技术(OCT)等的数据。数据存储装置102可以是，例如，计 算机存储器、CD-ROM、DVD、磁光(MO)盘、硬盘、软盘、zip盘、闪存驱动器等。图2示出根据本发明的一些实施例的成像系统。成像系统可以包括控制器205 ； 用于与调查者211(如人类用户、计算机等)进行信号通信的接口 206(例如图形用户界面、 键盘、操纵杆等)，且接口 206与控制器205同步；发射器201，其用于响应于来自控制器205 的控制信号生成并发射激发能202 ；以及检测器207，其用于生成投影像素数据104。投影 像素数据104可以以数字格式存储在存储装置中。该成像系统还可以包括移动或转动台 204，其用于容纳其上的对象203且相对于发射器201和检测器207可操地移动或转动；以 及图像重建器105、209，其与检测器207和控制器205电耦合，或经由如因特网、其它网络、 或数据总线的可选数据传输器208与检测器207和控制器205耦合。数据总线例如可以是 在计算机内部的计算机组件之间或在计算机之间传输数据的电线子系统。尽管所连接的组 件之间的数据流是单向的，但所连接的组件之间的通信可以是双向的。例如，激发能202可以是诸如X射线能、电磁(EM)能、光能、红外(IR)能、粒子能 (如电子、中子、原子束)的辐射能、诸如超声波的振动能等的形式。激发能202可以被照射 到例如可以是体模、病人、样本、或其任意组合的对象203上。可以通过相应能量类型的发 射器201发射激发能202。激发能202可以穿过对象203传播且部分由适当的检测器207 接收。检测器207可以将所接收的能量转换成可测量的电信号，从而可以将所测量的电信 号进一步转换成数字格式的投影像素数据104。控制器205可以是连接发射器201和检测器207的电路并且可以向发射器201和 检测器207发送控制信号，以使得激发能202的发射与检测器207的操作同步。该电路可 以是模拟的、数字的、或混合信号的。控制器205也可以是具有一个或多个处理器和一个或多个存储器的计算机，以将控制信号发送到发射器201和检测器207。图像重建器105、209可以响应于控制器205并接受投影像素数据104，以通过根据 本发明的一些实施例的方法重建对象203的图像，从而以高计算效率生成高保真的对象图 像。图像重建器105、209可以是具有一个或多个存储有用以根据本发明的示例性的实施例 来操作计算机的软件代码的数据存储装置。该计算机可以包括一个或多个处理器和一个或 多个存储器以读取并执行存储在一个或多个数据存储装置上的软件代码。在另一示例性的 实施例中，图像重建器105、209可以包括根据本发明的示例性的实施例可操作的一个或多 个程序存储和执行装置。图像重建器105、209可以生成重建的图像106。输出装置107、210可以接收重建的图像106。输出装置107、210可以是可视化装 置或数据存储装置。可视化装置例如可以是显示装置或打印装置。示例显示装置可以包 括例如阴极射线管(CRT)、发光二极管(LED)显示器、液晶显示器(LCD)、数字光投影(DLP) 监视器、真空荧光显示器(VFDs)、表面传导电子发射显示器(SED)、场发射显示器(FEDs)， 硅基液晶(LC0Q显示器等。示例打印装置可以包括例如基于碳粉的打印机、液体喷墨打印 机、固体墨打印机、染料升华打印机以及诸如热敏打印机和紫外线(UV)打印机的无墨打印 机等。打印装置可以以诸如三维造型装置的三维(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，并且根据方程O)，对于χ = 0的一些邻近处，函数f (χ)近似为 χ,|f (χ)-χ I < ε，对于适当小的ε >0。方程O)函数f可以去掉过多的值，尤其在初期迭代中。例如，一种这类函数是atanOO。 可以放宽限制函数并且可以没有严格的切除值a和b。另一示例限制函数可以是h(x)= c ·χ+(1-c) atan (χ))，0 < c < 1 ο Box, G. Ε· P.禾口 Jenkins, G. Μ·的 Time Series Analysis Forecasting and Control (San Francisco, CA :Holden-Day, 1976)nTffiTg^ f W 非线性数据变换的选择。在框302中，对于每个给定投影方向，可以将m维(m_D)像素的输入投影数据104 放置在分辨率改变(例如，在保持图像范围的同时每个后续分辨率网格比先前的分辨率网格具有更少的像素)的一系列分辨率网格上。例如，该系列中的分辨率网格可以具有前一 系列分辨率网格的分辨率的一半。例如，最粗的尺寸可以是单个m维像素。应当注意通常l<m<n且m对于所有的投影可以是相等的。然而，通常，m可 以随每个投影变化且m >= 0(例如，对于某些噪声估计问题)。在框302中，还可以在η维(n_D)空间中放置η维体素的初始集合，使得它们的投 影覆盖所述系列分辨率网格上的可用的m维投影像素。例如，这些初始体素的值的适当集 合可以是对应于平均投影密度的值；并且对于每个体素j，可以设置Aj = Iog(Vj)。可以使 用除对数函数以外的单调函数。例如，可以将初始的单个η维体素密度Vi设置为等于平均 投影密度。在框303中，当前网格分辨率的η维体素ν和m维像素ρ都可以可选地利用适当的 核进行平滑，其中m维核近似为η维核的投影。这个运算可以对应于具有噪声观察的卡尔 曼滤波器(其中由A得到ν)的预测算法(不更新)。例如，示例函数可以是ν = exp (vL)。 可以使用包括间断点的其它函数来解释附加的密度限制。示例性的平滑核可以具有sigma 宽度的高斯形状，例如，在0. 4和2个体素间隔之间。然而，如稍后关于框307所述，可以修 正。如果需要，可以通过映射将体素ν限制到特定的值。例如，在一些区域中，推理可 知体素密度为零。此外，为了避免极小体素密度的数值下溢，可以将 限制到最小可忽略 的值。示例性的限制值可以表示最大体素密度值Vkmax的约10_15的密度。在框304中，可以基于η维空间中的η维体素ν的集合计算预测的投影m维像素 Pp。预测的投影m维像素Pp可以在来自系列分辨率网格的分辨率网格上。如果从先前的分 辨率网格已经增加了网格分辨率，则可以通过从较粗的/先前的分辨率网格插值获得预测 的投影和体素估值V。在框305中，可以计算表征分辨率网格上预测的投影数据与输入投影数据104之
间的差的第一量。如果输入投影数据104的投影像素值Pi大于零，则第一量可以是例如，Ui = Mlog(Li)) 方程(3)其中Li是输入投影数据104的测量Pi与预测投影数据的预测像素值Ppi之间的误 差率，例如表示为，Li = Pi/Ppi 方程对于Li也可以使用其它适当的关系式，例如，诸如Li随Pi增加且随Ppi的增加而 减小的函数。这些函数可以捕获输入投影数据104与预测投影数据之间的测量误差。以上的log函数可以是误差率的对数，但可以是输入对于预测投影值的更普遍的 函数(例如，在透射电子显微镜中丢失楔问题中使用的函数)。如果输入投影数据104的投影像素值Pi不大于零，则第一量例如可以是，Ui =-a. 方程(5)其中a表示边界。在框306中，可以计算与η维空间中的单位强度的参考体素的至少一个脉冲响应 &相对应的量。依据当前的分辨率网格，可以对所有P投影和反投影针对每个期望的参考 η维体素计算脉冲响应民。脉冲响应艮可能取决于参考η维体素的位置。依据优选的运算序列，可以对于每个m维投影或对于η维体素计算脉冲响应&及其逆，优选的运算序列反 过来可以取决于期望的精度。这类似于FBP。可以选择参考η维体素的有限子集来估计脉 冲响应艮。在一些应用中，可以仅选择重建空间的中心处的一个η维参考体素。可以以几种方式对于相应的参考η维体素计算脉冲响应，并且如以后关于框307 所述，可以将这些脉冲响应的逆用于后续与第一量结合。在示例性的实施例中，在考虑单位强度的参考体素时，η维重建/对象空间的区域 中的特定位置处的脉冲响应&可以是线性系统的脉冲响应&。当特定位置处的该单个参考η维体素的中心被投影到所有P投影方向上时，可以 创建P单像素投影(不需要所有都是同一维度m)。当反投影这些P投影以重建时，可以创 建原始体素的模糊版本。该模糊版本关于原始体素可以是线性的，但是可以遍布η维空间。 该模糊函数&可以是对于该特定位置的重建系统的脉冲响应(或点扩散函数PSF)。对于有限的重建体积，该脉冲响应&可以取决于参考η维体素的位置。例如，当投 影的数量小且网格足够精细时，脉冲响应&可以是蜘蛛状(例如，其显示“蜘蛛腿”状的延 伸或从“蜘蛛”中心位置发射出的条纹)。邻近的体素可以类似地形成“蜘蛛腿”状的延伸， 但由于有限的重建区域可能在边界处切断延伸，因此，邻近的体素可以长度不同。此外，在 诸如使用P锥形束投影的锥形束的几何形状中，由于“蜘蛛腿”可以指向P锥的固定顶点， 因此，“蜘蛛腿”的方向和相互角可以随位置而改变。然而，为近似的目的，可以用这些脉冲响应函数的缓慢空间变换以在η维空间中 的足够小的邻近体素上共享同一脉冲响应函数。“蜘蛛腿”的相似度的几何考虑可以确定 共享同一近似脉冲响应的完全匹配体素区域的大小。例如，在“蜘蛛腿”的端点处相对于它 们的内容的几个(例如，十几个)百分比的整体不匹配，尤其当后续迭代可以校正这些误差 时，是完全可以接受的。相反，FBP方法通常不考虑这种截断效应，而如伍德中所示的矩阵 表述则考虑。还是在框306中，之后可以将给定体素的脉冲响应&与来自噪声方差的成分结合 以获得R，于是R将为第二量。可以测量输入投影数据以获得噪声方差。这种结合可以增加 “蜘蛛”中心处的值。依据优选的运算序列，可以对于每个m维投影或对于η维体素计算脉 冲响应&及其逆，优选的运算序列反过来可以取决于期望的精度。诸如当假定的噪声水平 为低时，较小的增量可以显著地增加对于精确计算R的逆(例如，对应于卡尔曼/维纳滤波 器增益)的需要。然而，便利地，多个分辨率网格上的计算可以将残留投影误差中的信号部 分保持在适度的/低的水平，从而支持R-1的有效计算。该信号部分可以是残留投影误差中 的结构信息，其不由根据输入投影数据测量的噪声方差解释。在框307中，之后可以将计算的第二量求逆以获得P，例如表示为P =『乂代替使 用逆海赛)。可以在傅里叶域中执行这个步骤。如下所述，本发明的示例性的实施例可以将 逆函数表征为傅里叶域中的卡尔曼或维纳增益。所测量的输入投影数据Y的观察模型可以是例如Y = HX+W，其中X可以表示对象 密度，H可以表示系统变换函数，并且Y可以表示输入投影数据104。在矩阵符号中，可以根据可用数据Y估计X。数据Y可以包括来自测量噪声的成 分。X可以具有以下统计特性)( = Ε[Χ]，及方程(6)
Vx = E [(X—Xo).(X—Xort方程⑵其中)(。是X的期望值E (即平均值)，Vx是对象密度X的协方差矩阵。可以将投影噪声W表示为E[W] = 0,方程(8)V = E [W · Wt]，及 方程(9)E [W · Xt] = 0方程(10)其中V是投影噪声W的协方差矩阵。例如，如 Sage (Sage，Α. P.，J. L. Melsa, Estimation Theory with Applications to Communications and Control, New York :McGraw_Hill，1971，page 262)中所述，可以 如下给出最小方差(卡尔曼)增益K K = VxHt [HVxH'+V]-1 = VxHtRy-1 方程(11)或者K = {tftl+Vm VX Ht)-1])"1方程(12)假设VxHt可逆，其中I为单位矩阵。在其中输入投影数据的像素数超过估计的图 像体素数的推荐的计算的初期低分辨率迭代中可以满足这个条件。在断层摄影技术的情况 下，Rf1对于求逆而言可能太大，但方程(12)可以建议H的有益变型以获得实际的K。通过 方程(11)，假设满秩的V可以从矩阵R/中去除可能的零特征值。H Vx Ht与V之间的关系 可以确定图像X的估值可以在多大程度上依赖于X的在前统计。使用复谱符号，对象的功率谱密度可以是Sxx，观察传递滤波器可以是H，投影噪声 谱功率密度可以是Svv，并且维纳增益(例如，参见Wiener和Norbert的Extrapolation， Interpolation,and Smoothing of Stationary Time Series(New York :Wiley,1949))可 以表示为W = H*SXX/(H*SXXH+SVV) 方程(13)= [H(U+Svv/(WSxxH)]-1 方程(14)其中*是共轭函数且U是单位传递函数。光谱表述可以计算地可行于断层摄影技术。在断层摄影技术重建中，对象密度的 不确定性(例如，Sxx或Vx可以相对大)可以支配投影噪声(例如，Svv或V可以相对小)， 并且找到维纳增益W(或卡尔曼增益K)的主要任务可以是找到方程(13)的有用近似(或 等价方程(11))。因此，如上所述，可以测量输入投影数据104以获得乂和Vx。例如，可以基于所测量 WVx计算H Vx Ht和H*SXXH。来自框306中提到的噪声方差的成分可以指由于例如V、H Vx Ht或H*SXXH导致的成分。尽管可以直接计算卡尔曼增益的方程(1 和维纳滤波器增益的 方程(14)，但是还可能进一步地简化。例如，对于测量噪声W的小的交叉相关(例如，V可 以是对角占优)和投影对象密度的适度的交叉相关(例如，H VxHt可以是对角占优，尤其是 在迭代过程中)，可以通过对于单个(或成组的)投影按块计算(block-wise)增益来简化 方程(1 或(14)的复合表达。例如，可以通过使用逐一体素(η维处理)或逐一像素(m维 近似处理)的例如修正的局部脉冲响应的逆，来进一步简化按块处理。局部脉冲响应的修 正可以由例如方程(12)中所示的噪声对信号结构V(H VxHt)-1的元素或类似地由方程(14) 中所示的SvvAH^xx H)的元素给出。可以逐一体素应用该修正。可以通过增加“蜘蛛”中心区域(如接近脉冲响应函数的峰，对应于方程(12)中的矩阵V (H Vx HV1或方程(14)中 的Svv/(H^xx H)的(近)对角元)的值来近似该修正。在W是白噪声且投影对象不确定性 以估计的均方根(RMS)噪声信号比r是白的情况下，例如可以将蜘蛛中心增加因子(1+r2)。 经过该修正后，可以通过根据例如方程(1 和(14)的数值求逆来计算卡尔曼增益和维纳 滤波器增益。如上所述，注意与如方程(11)和(13)中表示的上面提到的矩阵(或传递滤波 器)H相关。通过将卡尔曼增益K或维纳增益W中表示的单个步骤分成两个矢量运算步骤 可以改变处理序列。这两个步骤可以是首先反投影例如投影数据，然后用参考体素K的η 维滤波器I3k = R1T1滤波；体素K周围的一些邻近体素可以使用作为结果的滤波器输出。这 两个步骤也可以是首先利用pm维滤波器PiK(0 <= i < ρ)预滤波例如投影数据，其中可以 从由被噪声信号比修正的η维脉冲响应(PSF) Rpk的投影的逆获得每个PiK，然后反投影该滤 波器输出；再者，体素K周围的一些邻近体素(或其投影)可以使用作为结果的滤波器输出 (更新)。为了覆盖整个η维空间，可以使用几个体素位置和它们的邻近位置。因此，可以 将方程(11)和(1 中所示的超大矩阵的逆分解成数量小得多的集合&的适度集合的逆， 其中每个&仅可以在对象空间中缓慢变化，允许重复使用。此外，可以由傅里叶技术非常 有效地计算参考体素的每个IV1。总之，用于计算增益P的两个方面近似可以有助于高的统计效率。首先，对于低频，通过在可以解释足够宽频率范围内的噪声信号比的脉冲响应函 数的原点增加小的常数可以估计卡尔曼/维纳滤波器增益(现在由增益P替代)。通过在 原点处增加“窄”高斯钟(Gaussian bell)而不只是常数，可以修正接近脉冲响应函数中心 的小的近邻以解释例如“粉红”测量噪声。其次，对于较高频率，如高频信噪比所表征，可以渐变/平滑对应于增益P的滤波 器。为了提高数值精度，可以在更新当前分辨率网格上的图像密度和/或估计的图像密度 时执行该递减/平滑函数。在系列分辨率网格上的迭代计算期间的每一阶段，所估计的图 像密度的低通滤波可以帮助递减掉高频图像成分，从而减少计算负担。在框308中，可以使用第一量和P计算更新值。例如，P可以与第一量结合以获得 给定体素(或像素)图像元素的原始更新值i。例如，P可以与第一量卷积以获得原始更 新值。随着P在η维重建空间中缓慢变化，通过将P作为近似值应用于给定体素的适当的、 仍足够小的邻近体素，可以节省计算时间。然后，可以将与上述方程⑴和⑵中的f类似 构建的限制函数g应用于原始更新值i以获得修正的更新值ig = g(i)。该修正可以防止 过度校正。在开发本发明的示例性的实施例时，直接使用原始更新值的许多实验都失败了。 于是，发明人设想并开发了限制函数及其的正确使用，以限制过度更新值。该修正的更新值 可以是后续处理的更新值。在框304-308中，对于投影ρ的数量可能小的(ρ << N，其中N可以表示沿着穿 过重建体积的代表性投影路径的体素的数量)的锥形束几何重建(例如，如血管造影术中 所使用)，额外考虑简化可能有帮助。在这种情况下，尤其是在多重建网格的高分辨末级期 间，单个“蜘蛛腿”可以不同。对于该高分辨η维空间中的许多小区域计算脉冲响应Ι ρ可 能是耗时的。作为使用η维脉冲响应的近似，可以使用维度为Hii (0<=i< ρ)的ρ较低 维投影的脉冲响应。然后，可以将相应的修正的逆或卡尔曼/维纳增益迭代地应用到投影的单个误差图像上。这些m维脉冲响应可以是上述η维脉冲响应在各m维投影空间上的投 影，这些投影可以是蜘蛛状的。它们的属性可以是如上所述，就是说，在投影空间中缓慢变 化。然而，对于计算P投影的脉冲响应，计算负担、这些计算的脉冲响应的逆，以及接下来的 卷积(现在应用到m维投影)可以远小于在η维空间中计算它们，尤其当ρ << N时。可 以以类似于修正R的方式修正计算的脉冲响应。还可以带通滤波计算的脉冲响应以减小断 层摄影重建中的噪声。这种用于断层摄影求逆的技术可以与在重建之前滤波投影自身的滤 波反投影类似。该技术可以不限于锥形束几何形状，而是，例如也可以应用于诸如平行束几 何形状的其它的几何形状。在框309中，可以将修正的更新值ig加到η维体素ν的集合中的相应图像密度的 值上，以修正该η维体素ν的集合。通过与上述方程(1)和( 类似的限制函数可以限制 该相加。可以对于η维体素的整个集合的一部分执行所述加法。为了精确，可以使体素值 ν成比例，使得它们的投影与原始图像数据最佳匹配。在框310中，如果分辨率网格最细且如果所计算的预测与输入投影数据之间的差 (例如，残差)(统计上)足够小，则可以将代表密度估值V、Vl或这两者(以及如果需要， 代表估计的预测和误差)的η维体素的修正集合输出到输出装置，如框311所示。可以在 输出之前应用平滑。否则，分辨率网格可以被改变为具有比当前分辨率网格更大尺寸的分 辨率网格(或者可能可以减小η维平滑核的刚度/范围)，并且在框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 of Spiroplasma melliferum BC3 and Its Interactions with the Cell Membrane, " (J. Mol. Biol. , 378, pp. 776-787, 2008)中 的投影像素数据104被用于证明通过本发明的示例性的实施例可得到的改进。图4A示出 基于投影像素数据104的滤波反投影(FBP)重建的重建图像106。图4B示出根据本发明的 示例性的实施例的重建图像106。正如通过诸如细胞内与细胞外的空间对比、细胞膜、核糖 体、螺原体的细节等的不同密度的更锐利的对比度所测量的，图4A与4B之间的比较说明图 4B中的较好的图像清晰度。图4B中的图像更令用户满意。两种重建都使用在大约+/-70 度范围内的约166双倾斜透射电子显微镜(TEM)投影。图5示出仿真体模的正面样本投影。该仿真体模在大小为256X256X64体素的 3维盒子中的任意位置处包含大小为6X6X6体素的大约100个立方体。这些体素的5个 单位的密度嵌入0. 85个单位的周围密度中。通过计算机仿真103获得具有附加的白噪声 的总共49个投影(单倾斜，+/-72度)，以生成投影像素数据104。根据统计模型验证中的著名的金标准，通过检查残差可以最精密地评估重建质 量，其中，残差是真实图像(例如，图5中的仿真体模，已知为推理的)与重建图像106之间 的差异。高保真的重建可以导致残差是随机的且包含来自3维对象结构(例如，图5中的仿 真体模中的立方体)的非常小的成分。图6A示出同步迭代重建技术(SIRT)的残差，而图6B示出由根据本发明的示例性的实施例的方法获得的那些残差。图6A中的颇为显著的残 留方格图案(经过15次迭代)和图6B中的主导的白噪声外观(经过一次高分辨率迭代)， 进一步说明本发明在图像保真和计算效率上的优势。由于较多或较少的SIRT迭代只会导 致图像质量的下降，因此，对于这个例子优化了 SIRT迭代计数。本发明也可以适用于投影数量小的情况，例如，血管造影术中的锥形束几何形状。 图7A示出在向由塑料管制成的搏动液压心脏体模中注入碘造影剂(例如，如介入性心脏减 影血管造影中通常使用)期间的一个X射线投影。使用的投影像素数据104来自由于注射 的瞬态本性所导致的潜在体模的仅6个数字减影投影。图7B示出通过本发明的示例性的 实施例的重建图像106，其精确地表示了搏动液压体模的几何形状。由于投影的数量通常小 和/或迭代的数量用当前的计算资源变得不实际，因此，通过滤波反投影(FBP)或同步迭代 重建技术(SIRT)或平均最大熵(MENT)方法不能获得重建图像。如上所述，本发明的一些示例性的实施例结合了计算/处理技术且导致感兴趣区 域中的期望的η维体素的快速而精确的计算/估计，同时允许使用任意m维投影(和方向)。 注意m可以不是固定的，且其可以取决于投影指数。技术结合的效率可以适用于宽范围的 投影数量。可用投影(测量)信息的有效使用可以使得当投影数量非常大时和诸如在冠状 血管造影术中投影数量有限(小或不完整)时该方法都有效。创造性的非线性求逆技术可以结合(i)通过将η维图像估计分解成重建区域中 的两个非线性相关的域，将非线性卡尔曼滤波器的某些部分线性化(近似贝叶斯估计)； (ii)通过使用图像平滑度的降低(可比于“步骤”)多网格计算可比于用分辨率的进一步 细分的“楼梯”；及(iii)当需要时，代替使用逆海赛，高效计算由图像噪声的估计修正的线 性系统的逆脉冲响应(或传递函数)，产生对象空间中的选择区域的近似卡尔曼/维纳滤波 器增 O特别地，Sage，JarischW. R.和 Detwiler J. S. ^ "Statistical Modeling 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 by Limited Scan X-ray Tomography,Technical Report USAF AFWAL-TR-81-4050 NTIS, HC A07/MF AOl, Defense Technical Information Center,Cameron Station,Alexandria,Virginia 22314,1981 中已经讨论了卡尔曼滤波器 对于图像重建的应用。可以使用例如诸如FBP、FFT和维纳滤波器的替代的方法代替逆阵中 的线性部分。重建区域中的两个非线性相关域的任一个可以表示卡尔曼系统状态。通过以上结合，与在优化中有效使用海赛类似，利用每个迭代可以实现高效误差 降低。所实现的效率可以接近于二次收敛的效率。该方法可以回避采用其它方法时遇到的 高维度和差收敛率的严重问题，其它方法例如诸如是Kolehmainen，V.和Brandt，S. S.在 Structure-From-Motion Without Correspondence From Tomographic Projections by Bayesian Inversion Theory(IEEE Trans. Medical Imaging, ITMID4, Vol. 26, No. 2, pp. 238-248, Feb. 2007)中采用的那些方法。可替代地，本发明的示例性的实施例可以被看作是顶与(空间可变的/不稳定的 /分段的)FBP的结合，此处以术语S-FBP表示。依据近似的需求，可以有与S-FBP相关的一 个或几个段。通常，FBP在重建的领域上不是稳定的。然而，通过将重建的领域分成用这些区域之间的权重过渡的重叠的区域/段，领域属性的缓慢变化可以允许近似。通过适当选 择并结合如下所述的几个实施方式，可以实现S-FBP的速度、效率和灵活性。例如，可以将S-FBP部分应用于非线性变换数据，而不是原始投影数据。例如，多网格计算可以支持线性化且保持每次迭代的运算计数非常低。例如，最 精细分辨率网格上的最后计算之前的总的运算计数可以小于最精细分辨率网格上的S-FBP 的运算计数。例如，可以将体素密度表现在两个非线性相关的域(经由映射连接)中一个用于 经由S-FBP(其依据期望的近似技术可在η维空间中变化)进行更新，并且一个用于计算投 影预测。S-FBP部分的滤波器权重可以取决于体素位置，但是利用相同的权重将体素集中在 一起来使用可以提高计算效率。例如，通过根据线性化系统脉冲响应(或传递函数)的逆计算S-FBP部分的滤波 器权重，可以使用任意投影维度和方向。可以将线性化系统脉冲响应(或传递函数)与根 据输入投影数据的测量确定的噪声变量结合。例如，适当地选择限制(切除)函数可以防止计算期间数值奇异和不稳定。例如，适当地选择平滑函数，每次迭代时明智地应用于体素和像素值，可以有助于 数值计算的稳定性和精确度。在生成更新值之前，体素和像素值的平滑可以功能性地比得 上卡尔曼预测算法。平滑的迭代应用可以与卡尔曼滤波器中的状态预测算法类似。由于最精细分辨率网格上的计算之前的迭代的精度可能是足够的，通过最精细分 辨率网格的迭代在计算上可以不需要，因此可以节省误差校正项的计算。在这些情况下，重 建速度可以变得比得上S-FBP的重建速度。这可能是由于最精细分辨率网格上的计算可以 实现比得上FBP的运算计数(在因子> 1中)，而对于较粗分辨率网格上的所有之前低分辨 率迭代的结合运算计数可以小于最精细分辨率网格上的计算的结合运算计数，尤其对于较 高维度的图像重建。如关于图3所示的方法所述，例如，本发明的示例性的实施例可以提供为存储在 诸如⑶-ROM、DVD、磁光(MO)盘、硬盘、软盘、压缩盘、闪存驱动器等的数据存储装置上的软 件代码。例如，所存储的软件代码可以是由具有一个或多个处理器、一个或多个存储器装置 的计算机可读和可执行，以执行如上关于图3所述的示例性的方法，其中，存储器装置例如 是诸如随机存取存储器(RAM)装置、动态RAM(DRAM)装置、闪存装置和静态RAM(SRAM)装置寸。本发明的示例性的实施例可以提供一个或多个程序存储和执行装置以存储并 执行如关于图3所述的示例性的方法，例如，专用集成电路(ASICs)、场可编程门阵列 (FPGAs)、复杂可编程逻辑器件(CPLDs)、专用指令集处理器(ASIPs)等。此处说明的例子和实施例是非限制性的例子。关于示例性的实施例详细说明了本发明，且现在对于本领域的技术人员根据前述 以下将变得明显在较宽的方面，进行改变和变型而不偏离本发明，因此，在权利要求中所 限定的意在涵盖落入本发明的真实精神范围内的所有这些改变和变型。
权利要求
1.一种系统，包括一个或多个发射器，用于将激发能射入观察的对象；一个或多个检测器，用于响应于射入观察的对象的激发能来生成对于由所述一个或多 个检测器所接收的能量进行编码的投影空间数据；控制器，用于控制所述一个或多个发射器发送所述激发能，并控制一个或多个接收器 生成所述投影空间数据；以及图像重建器，用于接收所述投影空间数据并如下处理所述投影空间数据计算表征所述投影空间数据与预测投影数据之间的差的第一量；计算对应于至少一个脉冲响应的第二量，每个脉冲响应对应于单位强度的参考体素；使用所述第二量的逆函数和所述第一量计算更新值；以及使用所述更新值重建表示所述观察的对象的对象空间图像。
2.如权利要求1所述的系统，其中，所述激发能是电磁(EM)能、X射线能、粒子束、红外 能、光能或振动能中的至少一个。
3.如权利要求1所述的系统，其中，基于对象空间体素的集合来计算所述预测投影数 据，当将所述对象空间体素的集合投影到投影空间时，所述对象空间体素的集合覆盖多个 不同尺寸的预定分辨率网格。
4.如权利要求3所述的系统，其中，使用所述第二量的逆函数和所述第一量计算更新 值包括将所述第一量与所述第二量的逆函数卷积。
5.如权利要求1所述的系统，其中，使用所述更新值重建表示所述观察的对象的所述 对象空间图像包括利用所述更新值修正所述对象空间体素的集合。
6.如权利要求1所述的系统，还包括输出装置，用于接收表示所述观察的对象的所述对象空间图像，其中所述输出装置是 数据存储装置、显示装置或打印装置中的至少一个。
7.如权利要求1所述的系统，还包括调查者计算机，用于接收所述对象空间图像，并且被编程为从所述对象空间图像提取 诊断信息或对于所述一个或多个发射器、所述一个或多个检测器或所述图像重建器中的至 少一个细调参数。
8.一种工作站，包括 一个或多个处理器；以及一个或多个数据存储装置，用于存储由所述一个或多个处理器执行的软件，所述软件 包括用于接收表示从对象空间投影到投影空间的对象的输入数据的软件代码； 用于计算表征所述输入数据与预测投影数据之间的差的第一量的软件代码； 用于计算对应于至少一个脉冲响应的第二量的软件代码，每个脉冲响应对应于所述对 象空间中的位置处的单位强度的参考体素；用于使用所述第二量的逆函数和所述第一量计算更新值的软件代码；以及 用于使用所述更新值重建表示所述观察的对象的对象空间图像的软件代码。
9.如权利要求8所述的工作站，其中，从投影图像获取系统、所述一个或多个数据存储 装置、另外一个或多个数据存储装置、由所述一个或多个处理器执行的软件仿真或由另外 一个或多个处理器执行的软件仿真中的至少一个接收所述输入数据。
10.如权利要求8所述的工作站，其中，经由数据总线、线缆、有线网络、或无线网络接 收所述输入数据。
11.如权利要求8所述的工作站，还包括输出装置，用于接收图像，其中，所述输出装置是所述一个或多个数据存储装置、另外 一个或多个数据存储装置、显示装置或打印装置中的至少一个。
12.如权利要求8所述的工作站，其中，所述软件还包括用于将所述第一量与所述第二 量的逆函数卷积的软件代码。
13.一种由一个或多个处理器执行存储在一个或多个数据存储装置上的软件代码来实 施的方法，所述方法包括从投影图像获取系统、所述一个或多个数据存储装置、另外一个或多个数据存储装置、 由所述一个或多个处理器执行的软件仿真或由另外一个或多个处理器执行的软件仿真中 的至少一个，接收表示从对象空间投影到投影空间的对象的输入数据；计算表征分辨率网格上的所述输入数据与所述分辨率网格上的预测投影数据之间的 差的第一量；计算对应于至少一个脉冲响应的第二量，每个脉冲响应对应于所述对象空间中的位置 处的单位强度的参考体素；使用所述第二量的逆函数和所述第一量计算更新值；使用所述更新值重建表示所述观察的对象的图像；以及将重建的图像输出到输出装置，其中，所述输出装置是所述一个或多个数据存储装置、 另外一个或多个数据存储装置、显示装置或打印装置中的至少一个；其中，根据存储在所述一个或多个数据存储装置上的所述软件代码，所述一个或多个 处理器执行接收所述输入数据、计算所述第一量、计算所述第二量、计算所述更新值、重建 所述图像、和输出重建的图像。
14.如权利要求13所述的方法，其中，从多个不同尺寸的预定分辨率网格获得所述分 辨率网格。
15.如权利要求14所述的方法，其中，基于对象空间体素的集合计算所述预测投影数 据，当将所述对象空间体素的集合投影到投影空间时，所述对象空间体素的集合覆盖所述 多个不同尺寸的预定分辨率网格。
16.如权利要求15所述的方法，还包括利用平滑核对所述分辨率网格上的所述输入数据或所述对象空间体素的集合中的至 少一个进行平滑。
17.如权利要求15所述的方法，其中，所述对象空间体素的集合被限制在一定值范围内。
18.如权利要求13所述的方法，其中，所述第一量被表征为所述输入数据与所述预测 投影数据之间的测量误差。
19.如权利要求13所述的方法，其中，所述第二量包括来自由测量所述输入数据获得的噪声方差的成分。
20.如权利要求13所述的方法，其中，对于覆盖所述对象空间中的所述单位强度的参 考体素的位置的周围区域计算所述逆函数。
21.如权利要求13所述的方法，其中，重建所述图像包括将所述第一量与所述第二量的逆函数卷积，以获得所述更新值。
22.如权利要求13所述的方法，其中，所述更新值由限制函数进一步限制。
23.如权利要求13所述的方法，还包括 使用所述更新值修正所述对象空间体素的集合。
24.如权利要求23所述的方法，还包括在修正了所述对象空间体素的集合后，改变为具有比所述分辨率网格更大尺寸的不同 分辨率网格。
25.权利要求13中的一个或多个数据存储装置，用于存储由一个或多个处理器执行的 软件以执行如权利要求13所述的方法。
全文摘要
一种系统，包括发射器，用于将激发能射入观察的对象；检测器，用于响应于射入对象的激发能来生成对接收的能量进行编码的投影空间数据；控制器，用于控制发射器并控制接收器生成投影空间数据；图像重建器，用于接收投影空间数据并如下处理所述投影空间数据计算表征投影空间数据与预测投影数据之间的差的第一量；计算对应于至少一个脉冲响应的第二量，每个脉冲响应对应于单位强度的参考体素；使用第二量的逆函数和第一量计算更新值；以及使用所述更新值重建表示观察的对象的对象空间图像。
文档编号G01N23/00GK102132149SQ200980123522
公开日2011年7月20日 申请日期2009年6月29日 优先权日2008年6月27日
发明者沃尔弗拉姆·R·雅里施 申请人:沃尔弗拉姆·R·雅里施