CN101056582A - 用于检查感兴趣对象的计算机断层摄影装置和方法 - Google Patents

Abstract

一种用于检查感兴趣对象(107)的计算机断层摄影装置(100)，该计算机断层摄影装置(100)包括探测元件(123)，适于以能量分辨方式探测从感兴趣对象(107)相干散射的电磁辐射，以及确定单元(118)，适于根据对从探测元件(123)接收的探测信号进行的统计分析，确定关于感兴趣对象(107)的结构信息。

Description

用于检查感兴趣对象的计算机断层摄影装置和方法

技术领域

本发明涉及X射线成像领域。特别地，本发明涉及计算机断层摄影装置，涉及使用计算机断层摄影装置检查感兴趣对象的方法，涉及计算机可读介质和程序单元。

背景技术

在过去的几年中，X射线行李安检已经从完全依赖于操作者的干涉的简单X射线成像***，发展到可以自动辨认某些类型的物质并且在存在危险物质时触发警报的更完善的自动***。安检***已经使用了X射线辐射源，用于发射X射线，X射线透射穿过所检查的包裹或者从其散射回探测器。

计算机断层摄影(CT)是一种使用数字处理而从围绕旋转的单一轴获取的一系列二维X射线图像中产生对象内部的三维图像的方法。可以通过应用适当的算法来完成CT图像的重建。

基于相干散射的X射线光子或量子的成像技术是所谓的“相干散射计算机断层摄影”(CSCT)。CSCT是产生感兴趣对象的散射特性(尤其是小角度)的图像的技术。这些取决于对象的分子结构，使可以产生每个部分的特定材料图。小角度散射的主要部分是相干散射。由于相干散射光谱取决于散射样本的原子排列，相干散射计算机断层摄影(CSCT)是用于穿过二维对象截面对行李或生物组织的分子结构的空间变化成像的灵敏技术。

所需的是拥有一种高度有效的图像重建方案。

发明内容

本发明提供了一种计算机断层摄影装置、一种使用计算机断层摄影装置检查感兴趣对象的方法，一种具有根据独立权利要求的特征的计算机可读介质和一种程序单元。

根据本发明，用于检查感兴趣对象的计算机断层摄影装置包括适于以能量分辨方式(energy-resolving manner)探测从感兴趣对象相干散射的X射线的探测元件，以及适于根据对从探测元件接收的探测信号进行的统计分析而确定关于感兴趣对象的结构信息的确定单元。

此外，本发明提供了一种用计算机断层摄影装置检查感兴趣对象的方法，该方法包括这些步骤，以能量分辨方式探测从感兴趣对象相干散射的X射线，并且根据对探测的信号进行的统计分析来确定关于感兴趣对象的结构信息。

而且，本发明提供了一种计算机可读介质，其中存储了用于使用计算机断层摄影装置检查感兴趣对象的计算机程序，其当由处理器执行时，适于执行上述方法步骤。

除此之外，提供了用于检查感兴趣对象的程序单元，其当由处理器执行时，适于执行上述方法步骤。

通过计算机程序、即通过软件，或通过使用一个或多个特定电子最优电路、即以硬件形式，或者以混合形式、即借助于软件部件和硬件部件，实现根据本发明对感兴趣对象的检查。

本发明的特征尤其具有这样的优点，即，与纯确定性分析相比，使用统计分析实现了从能量分辨探测信号重建三维图像，所述分析的对象是从多个探测元件接收的探测信号。因而，本发明的重建方案考虑光子统计，并且优选允许使用每个光子各自的位置和能量信息进行重建。随着能量分辨探测器实现在本发明的CSCT***中，可以使用每个光子的信息实现可能的最高图像质量。因而，本发明引入了一种对相干散射计算机断层摄影的统计评估模型，从而可以计算概率分布，且其根据在先执行的统计分析来确定，其中感兴趣对象的三维图像确实存在的概率高或最高。因而，与纯确定性模型相比，本发明考虑了获自相干散射电磁辐射、尤其是X射线的能量分辨统计分析的随机信息。

优选用于本发明的统计分析的统计模型是所谓的极大似然模型。根据“极大似然”的统计分析包括具有关于灵敏度的显著优点的重建运算法则。优选根据用于CSCT的获取的对数似然函数，可以分别地考虑每个独立光子的能量信息和位置信息。此外，可以在分析中引入能量相关衰减图，以及任何轨线和任何探测器形状。本发明提供了极大似然函数，其考虑了每个独立光子的贡献，因而提高了信噪比。

通过使用大量可获取的信息，即每个散射光子各自的信息，可以从测量结果获得最大信息，并且因而计算机断层摄影分析的精确度显著增加。

本发明提供了用于CSCT的对数似然函数，其可以通过导出的迭代更新步骤被最大化。对数似然函数考虑光子统计并且分别考虑每个光子的信息。另外，本发明的方法对于任何轨线、探测器形状、束几何形状(例如，扇形束、锥形束等)有效，并且其支持能量相关衰减图。

由于对数似然对于每个q-平面各自有效，对于每个q-平面的独立计算是可能的，其促进了并行处理。

在统计学中，“极大似然”是点估计的方法，其用作无法观察的总体参数、最大化似然函数的参数空间的成员。变量h表示将要估计的可观察总体参数。X表示观察到的随机变量。观察到的结果X＝x(或者随机变量X的概率密度函数的x处的值作为x保持固定的h的函数)的概率是似然函数L_X＝x(h)：

L_X＝x(h)＝P(X＝x|h)    (1)

最大化L(h)的h值是h的极大似然估计值。对似然取对数产生了所谓的对数似然项。

本发明提供了照顾到光子统计和考虑每个光子各自信息的对数似然函数。另外，本发明的方法对于任何轨线和束几何形状(例如，扇形束、锥形束)有效。由于对数似然对于每个q-平面各自有效，对于每个q-平面的独立计算是可能的，其促进了并行处理。与全面对数似然最大化相比，这还增加了重建的速度。另外，对数似然支持能量相关衰减图。

本发明可以有利地实施在相干散射计算机断层摄影(CSCT)中，用于单层探测器和用于多层探测器，尤其其中至少一部分探测器是能量分辨的。

因而，本发明涉及CSCT/CT设备，其中通过考虑到至少一些独立光子统计数字(位置和能量)的测量值而重建图像，并且其中使用测量值和极大似然模型重建图像。

根据本发明，可以考虑并且可以由探测元件测量光子的任何能量值。因而，本发明可以-但不局限于-仅使用两个不同的能量值，但也可以允许至少三个不同的能量值，或者可以允许所探测的光子的连续(或准连续)能量谱。

参考从属权利要求，还将在下文中描述本发明的优选实施例。

接着，将描述本发明的计算机断层摄影装置的优选实施例。这些实施例也可以应用于使用计算机断层摄影装置检查感兴趣对象的方法、应用于计算机可读介质或应用于程序单元。

在计算机断层摄影装置中，确定单元可以适于根据探测信号的位置和能量的分析而确定关于感兴趣对象的结构信息。通过实现能量分辨探测元件，通过考虑-在统计分析模型的框架中-在由多个探测元件形成的阵列探测器处的探测信号的位置，以及所探测的信号的能量，所有可用的信息可以用于重建感兴趣对象的三维图像。

此外，确定单元可以适于根据对电磁辐射的每个量子各自的探测信号的分析而确定关于感兴趣对象的结构信息。对探测元件探测得到的各个光子中的每个的考虑，允许增加本发明统计模型的结果的可靠性。

该确定单元可以适于根据探测信号的极大似然分析而确定关于感兴趣对象的结构信息。虽然任何备选统计模型可以应用于本发明的计算机断层摄影装置，但是将实现的优选统计模型是极大似然，因为该结构尤其适于使用每个独立光子的能量和位置信息，因而增加了感兴趣对象的重建图像的效率和精确性。

确定单元可以适于根据探测信号的对数似然分析，即通过分析极大似然函数的对数，来确定关于感兴趣对象的结构信息。

确定单元可以适于在考虑了吸收系数的实验和/或理论预定或预知值的情况下，确定关于感兴趣对象的结构信息。吸收系数的试验估计可以包括执行透射断层摄影的参考扫描。用于吸收系数的可靠值也可以从相当合理的近似值(例如，参见2003年Med.Phys.30(9)，van Stevendaal，U.等人的“A reconstruction algorithm for coherentscatter computed tomography based on filtered back-projection”)导出。

确定单元可以适于在考虑了能量相关衰减图的情况下确定关于感兴趣对象的结构信息。通过考虑穿过感兴趣对象的电磁辐射衰减的能量相关性，本发明分析结果的可靠性增加。

计算机断层摄影装置的确定单元可以适于在考虑了探测元件的接受度的能量相关性的情况下确定关于感兴趣对象的结构信息。考虑到测量来自感兴趣对象的特定体素(像素的三维类似物)的特定光子的接受度或几何灵敏度，由于其改善了所用的模型，允许重建具有高精确度的感兴趣对象的图像。

而且，确定单元可以适于根据探测信号的迭代分析来确定关于感兴趣对象的结构信息。这种迭代过程-尤其与极大似然方法结合-是获取所需信息的直接方案。

确定单元可以还根据单独为动量传递的不同值、尤其对动量传递的每个平面q执行的探测信号的分析，来确定关于感兴趣对象的结构信息。从光子到散射中心(例如，原子)的这种动量传递可以在散射过程期间发生。特别地，对数似然模型对于每个q-平面有效，从而对于每个q-平面的独立计算可以促进并行处理。

该计算机断层摄影装置的探测元件可以适于探测预定数目的电磁辐射的量子，或者对预定时间间隔探测电磁辐射的量子。通过将分析限制到将探测的给定数目的光子或者作为选择地限制到给定辐射时间，来自不同测量/扫描的分析结果更可比，从而可以获得高度可复制和有意义的结果。

本发明的计算机断层摄影装置可以适于作为相干散射计算机断层摄影装置(CSCT)，即该计算机断层摄影装置可以根据上述的CSCT技术配置和操作。

该计算机断层摄影装置还可以包括适于向感兴趣对象发射电磁辐射的电磁辐射源，并且可以包括布置在电磁辐射源和探测元件之间的准直器，该准直器适于使电磁辐射源发射的电磁辐射束准直，以形成扇形束、锥形束或者任何随其他所需几何形状的束。作为扇形束结构的备选，可以使用锥形束结构。

计算机断层摄影装置的探测元件可以形成单层探测器阵列，或者可选的多层探测器阵列。单层探测器阵列具有结构简单和快速评估探测信号的优点。然而，可以实现多层探测器阵列以达到探测信号的特别高的分辨率，以及经探测信号的较大数量。

根据本发明的计算机断层摄影装置可以配置为以下组中的一个，该组包括行李安检装置、医疗应用装置、材料测试装置和材料科学分析装置。本发明应用的优选领域是行李安检，由于本发明的已定义功能允许对行李内容物的安全和可靠分析，允许探测可疑内容物，甚至允许确定在这种行李中的物质类型。本发明创建了一种高质量自动***，其可以自动辨认某些类型的物质，并且如果需要，在存在危险物质时触发警报。这种安检***可能已经使用了本发明的计算机断层摄影装置，其具有X射线辐射源，用于发射透射穿过受检包裹或者从受检包裹散射到探测器的X射线，这允许以能量分辨方式探测相干散射辐射。

上面定义的方面以及本发明的其他方面从将在下文中描述的实施例的范例中是显而易见的，并且参考实施例的这些范例而说明。

下文将参考实施例的范例更详细的描述本发明，而本发明不局限于所述实施例的范例。

附图说明

图1示出了根据本发明典型实施例的计算机断层摄影装置；

图2和图3示出了根据本发明的计算机断层摄影装置的示意图，基于其说明了本发明的统计分析的典型实施例；

图4示出了将在本发明的计算机断层摄影装置中实现的数据处理设备的典型实施例。

附图中的图示是示意性的。在不同的附图中，相似或相同的元件具有相同的参考标记。

具体实施方式

图1示出了根据本发明的CSCT(相干散射计算机断层摄影)扫描器***的典型实施例。参考该典型实施例，将描述本发明应用在行李安检中以探测行李中的诸如***物的危险物质。然而，应当注意到，本发明不局限于该应用，也可以应用于医学成像领域，或者诸如材料测试的其他工业应用中。

图1中描述的计算机断层摄影装置100是锥形束CT扫描器。然而，本发明也可以实现为具有扇形束几何形状。图1中描述的CT扫描器包括构台(gantry)101，其可围绕旋转轴102旋转。构台101借助于马达103驱动。参考标记104指示诸如X射线源的辐射源，其根据本发明的一方面，发射多色辐射。

参考标记105指示孔径***，其将从辐射源发射的辐射束成形为锥形辐射束106。导向锥形束106使其穿过布置在构台101中心、即在CSCT扫描器的检查区域中的感兴趣对象107，并且撞击在探测器108上。如可以从图1中获取的，探测器108布置在相对于辐射源104的构台101上，从而探测器108的表面由锥形束106覆盖。图1中描述的探测器108包括多个探测器元件123，每一个能够以能量分辨方式探测已经由感兴趣对象107相干散射的X射线。

在扫描感兴趣对象107期间，辐射源104、孔径***105和探测器108沿着构台101在箭头116所指示的方向上旋转。为了使构台101与辐射源104、孔径***105和探测器108一起旋转，马达103连接到马达控制单元117，其连接到计算或确定单元118。

在图1中，感兴趣对象107是放置在传送带119上的行李。在扫描感兴趣对象107期间，在构台101围绕行李107旋转的同时，传送带119沿着平行于构台101的旋转轴102的方向移动感兴趣对象107。这样，沿着螺旋扫描路径扫描感兴趣对象107。在扫描期间，传送带119也可以停止，以便由此测量信号层(slice)。作为提供传送带119的代替，例如在感兴趣对象107是患者的医疗应用中，使用移动台。然而，应当注意到，在所有所述情况下，也可以执行圆扫描，其中在平行于旋转轴102的方向上不存在位移，而仅围绕旋转轴102旋转构台101。

此外，应当强调的是，作为图1中所示锥形束结构的备选，本发明可以由扇形束结构实现。为了产生主扇形束，孔径***105可以配置为狭缝准直器。

探测器108连接到确定单元118。确定单元118接收探测结果，即来自探测器108的探测器元件123的示值读数，并且根据这些示值读数确定扫描结果。而且，确定单元118与马达控制单元117通信，以便用马达103和120调整构台101和传送带119的移动。

确定单元118可以适于使用统计学方法从探测器108的示值读数重建图像。由计算单元118产生的重建图像可以经由接口122输出到显示器(图1中未示出)。

确定单元118可以由数据处理器实现，从而处理来自探测器108的探测器元件123的示值读数。

此外，如可以从图1中获取的，确定单元118可以连接到扩音器121，例如在行李107中探测到可疑物质的情况下自动输出警报。

用于检查感兴趣对象107的计算机断层摄影装置100，包括具有以类似矩阵形式布置的多个探测元件123的探测器108，每个适于以能量分辨方式探测从感兴趣对象107相干散射的X射线。此外，计算机断层摄影装置100包括确定单元118，适于根据对从探测元件123接收的探测信号的统计分析确定关于感兴趣对象107的结构信息。特别地，确定单元118根据对探测器阵列108上的探测信号的位置和能量的分析，确定关于感兴趣对象107的结构信息。优选对X射线源104产生的、并被相干散射到感兴趣对象107上的电磁辐射的每个光子，单独执行该分析。如将在下文中更详细描述的，确定单元118根据对探测信号的及大似然分析确定关于感兴趣对象107的结构信息。在该评估的框架中，关于感兴趣对象107的结构信息由确定单元118在考虑了吸收系数的试验预知值的基础上确定。为了确定关于感兴趣对象107的结构信息，确定单元118考虑能量相关衰减图。而且，考虑探测元件123接受度(acceptance)的能量相关性而重建该结构信息。由确定单元118执行的分析包括迭代。

计算机断层摄影装置100包括X射线源104，适于向感兴趣对象107发射X射线。提供在电磁辐射源104和探测元件123之间的准直器105适于准直自电磁辐射源104发射的电磁辐射束以形成锥形束。作为选择，图1中未示出，可以使用狭缝准直器作为准直器105的代替以产生扇形束。探测元件123形成多层探测器阵列108。计算机断层摄影装置100配置为行李安检装置。

在下文中，参考图2和图3，将更详细地描述确定单元118执行的统计分析。

参考图2，将说明方案200，其图示了相干散射CT的几何形状的划痕(scratch)和光子的可能飞行路径。方案200示出了X射线源104，其发射主束106，其在感兴趣对象107的体素k处散射，其中散射角度由角度θ定义。相干散射束202撞击一个探测元件123以在该特定探测元件上产生能量分辨探测信号。换句话说，光子在X射线源104发射，在体素k处在某角度θ下散射，并且在探测器像素i中测量光子的能量E_i。

在下文中，根据图2中图示的几何方案，将描述基于精确独立光子的导数的极大似然函数。

在探测器元件123中测量光子或量子i具有光子能量E_i。下列考虑的目标是计算该光子在具有坐标x_j的K个体素j的任何一个处相干散射的概率密度。所有吸收系数μ_j从公知的透射断层摄影扫描中精确已知，或者从相当合理的近似值近似已知，参见上述的van Stevendaal等人的工作。

在假设光子已经在具有坐标x_k的体素k处相干散射的假设下，由下列等式给出从源104到体素x_k的衰减A_att，plane ⁱ：

$A_{\mathrm{att},\mathrm{plane}}^i\left(x_k\right)=e^{-{\mathrm{\Σ}}_{j=1}^{k-1}{l}_{\mathrm{ij}}{\mathrm{\μ}}_j}---\left(2\right)$

其中l_ij是体素栅格的基函数，同时从体素x_k到探测器元件123的衰减A_att，scat ⁱ可以通过下列等式计算：

$A_{\mathrm{att},\mathrm{plane}}^i\left(x_k\right)=e^{-{\mathrm{\Σ}}_{j=k}^k{l}_{\mathrm{ij}}{\mathrm{\μ}}_j}---\left(3\right)$

光子i的总衰减A_att ⁱ(x_k)是

$A_{\mathrm{att}}^i\left(x_k\right)=A_{\mathrm{att},\mathrm{plane}}^i\left(x_k\right)*A_{\mathrm{att},\mathrm{scat}}^i\left(x_k\right)---\left(4\right)$

衰减考虑了穿过感兴趣对象107的电磁辐射的吸收。

由于光子能量已知并且在相干散射过程期间保持不改变，可以使用能量相关衰减图。

具有能量E_i的光子i到达探测器108的概率密度p(E_i|x_k)由下列等式给出：

$p\left(E_i\right|x_k)=p\left(A_{\mathrm{att},k}^i\right)\·P\left({\mathrm{\θ}}_i\left(x_k\right)\right)=\frac{A_{\mathrm{att}}^i\left(x_k\right)}{{\mathrm{\Σ}}_{m=1}^K{A}_{\mathrm{att}}^i\left(x_m\right)}\·\frac{1+\mathrm{co}{s}^2{\mathrm{\θ}}_i\left(x_k\right)}{1.5\mathrm{\π}}---\left(5\right)$

函数P(θ_i(x_k))定义了某一角度θ_i(x_k)下光子相干散射的概率(参见图2)，并且由Thompson公式给出，参见等式(5)。Thompson公式描述了光子在特定角度下散射的概率。

测量来自体素k的光子i的几何灵敏度(或接受度)s_k，i包括球面度(角度的三维模型)，其从探测器元件123可见。例如：在具有平行探测器的扇形束几何形状中(参见图3)，该灵敏度由下列等式给出：

$s_{k,i{,E}_i}=\frac{\mathrm{\Δ}{L}_{k,i}}{2{\mathrm{\πR}}_i}\·a_{E_i}---\left(6\right)$

该接受度包括例如由于光子被准直器的叶片吸收而阻止光子探测的可能性。

图3示出了扇形束结构301中来自具有主辐射的源104的X射线束106的方案300，其中主X射线束106在体素k处散射，从而相干散射束202撞击探测器108的探测器元件123。因而，图3示出了具有平行探测器108的扇形束301几何形状的几何灵敏度。在体素k处相干散射光子，并且其在像素i中的探测和其无散射的位置(参见虚线)之间的距离定义了距离R_i。探测器像素123的大小和准直定义了具有半径R_i的圆的长度ΔL_k，i。

在等式(6)中，ΔL_k，i是探测器元件123观察到的半径R_i描述的圆的长度，所述探测器元件123记录光子i。这考虑了探测器108上的所有准直。接受度a_Ei照顾到探测器像素123的能量相关接受度。

具有对固定k落入q_r，k的库(bin)r中的能量Ei的光子在体素x_k处散射和在探测器i测量的概率P(x_k|q_r，k(E_i，x_k))，由下列等式给出：

$P\left(x_k\right|q_{r,k}(E_i,x_k))=\frac{q_{r,k}(E_i,x_k)s_{k,i,E_i}}{{\mathrm{\Σ}}_{i=1,i\∈A}^N{\mathrm{\Σ}}_{m=1}^K{q}_{r,m}(E_i,x_m)s_{m,i,E_i}}---\left(7\right)$

其中q(E_i，x_k)＝2E_isin(Θ_i(x_k)/2)/hc且A_r，k＝{i|q_r，k＜q(E_i，x_k)＜q_r+1，k}，具有N_r，k个元素。为了简化，定义为 ${\mathrm{\Σ}}_{i=1}^{N_{r,k}}:={\mathrm{\Σ}}_{i=1,i\∈A_{r,k}}^{N_{r,k}}.$ 具有(q_r，1，...，q_r，k)的q的分布，通过寻找一组极大似然的q_r，k而确定：

$L(E_1,\·\·\·,E_n|q)=\underset{i=1}{\overset{N}{\mathrm{\Π}}}p\left(E_i\right|q)---\left(8\right)$

其由对测量的所有N个光子的概率密度p(E_i|q)给出。该概率p(E_i|q)由下列等式给出：

$p\left(E_i\right|q)=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}(p\left(E_i\right|x_k)\·P\left(x_k\right|q_{r,k}(E_i,x_k)))$

$=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}(\frac{A_{\mathrm{att}}^i\left(x_k\right)}{{\mathrm{\Σ}}_m{A}_{\mathrm{att}}^i\left(x_m\right)}\·\frac{1+\mathrm{co}{s}^2{\mathrm{\θ}}_i\left(x_k\right)}{1.5\mathrm{\π}}\·\frac{q_{r,k}(E_i,x_k)s_{k,i,E_i}}{{\mathrm{\Σ}}_{i=1}^{N_r}{\mathrm{\Σ}}_m{q}_{r,m}(E_i,x_m)s_{m,i,E_i}})$ (9)

对数似然函数是

$L(E_1,\·\·\·,E_N|q)=\underset{i=1,i\∈B}{\overset{N}{\mathrm{\Σ}}}\ln p\left(E_{i_r}\right|q)=\underset{i=1}{\overset{N_r}{\mathrm{\Σ}}}\ln \left(\underset{k=1}{\overset{K}{\mathrm{\Σ}}}p\left(E_i\right|x_k)q_{r,k}(E_i,x_k)s_{k,i,E_i}\right)$

$-N_r\ln \left(\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\underset{i=1}{\overset{N_{r,m}}{\mathrm{\Σ}}}{q}_{r,k}(E_i,x_k)s_{k,i,E_i}\right)$ (10)

其将被最大化，其中 $B={\∪}_{k=1}^K{A}_{r,k},$ 具有N_r个元素且 ${\mathrm{\Σ}}_{i=1}^{N_r}:={\mathrm{\Σ}}_{i=1,i\∈B}^N.$

等式(10)的对数燃函数可以-如果需要或者合适-由罚函数R(q)扩大，其可以是一阶邻域和成对像素差值(pairwise pixel difference)罚函数，如

$R\left(q\right)=\underset{i=1}{\overset{N_r}{\mathrm{\Σ}}}\frac{1}{2}\underset{k\∈x_i}{\mathrm{\Σ}}\mathrm{\Ψ}(q_j-q_k)---\left(10a\right)$

其中N_j是第j个像素的邻域(例如，左、右、上、下)，而ψ称为势函数(potential function)。然后，对数似然函数读取

$L(E_1,\·\·\·,E_N|q)=\underset{i=1,i\∈B}{\overset{N}{\mathrm{\Σ}}}\ln p\left(E_{i_r}\right|q)=\underset{i=1}{\overset{N_r}{\mathrm{\Σ}}}\ln \left(\underset{k=1}{\overset{K}{\mathrm{\Σ}}}p\left(E_i\right|x_k)q_{r,k}(E_i,x_k)s_{k,i,E_i}\right)$

$-N_r\ln \left(\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\underset{i=1}{\overset{N_{r,m}}{\mathrm{\Σ}}}{q}_{r,k}(E_i,x_k)s_{k,i,E_i}\right)+\mathrm{\βR}\left(q\right)$ (10b)

其中参数β设置成本函数(cost function)R(q)的强度。对于β＝0，重新获得等式(10)(关于函数R的详细信息，参见医学成像手册，第2卷，第1-70页，2000，Eds.：M.Sonka，J.M.Fizpatrick)。

在测量期间，测量时间T可以固定，而总计数N成为实验期间测得的附加随机变量，其可以表示为“当前时刻”。作为选择，计数N的总数量可以固定，该情况下测量时间T成为随机变量(表示为“当前计数”)。在“当前时刻”情况下，对数似然函数改变成：

$L(E_1,\·\·\·,E_N|q)=\underset{i=1}{\overset{N_r}{\mathrm{\Σ}}}\ln p\left(E_{i_r}\right|q)+P\left(N\right|T,q)+\mathrm{\βR}\left(q\right)---\left(11\right)$

而在“当前计数”情况下，改变成：

$L(E_1,\·\·\·,E_N|q)=\underset{i=1}{\overset{N_r}{\mathrm{\Σ}}}\ln p\left(E_{i_r}\right|q)+P\left(T\right|N,q)+\mathrm{\βR}\left(q\right)---\left(12\right)$

其中P(N|T，q)和P(T|N，q)分别描述“当前时刻”和“当前计数”的概率。

下文中，将描述表示的对数似然函数的可能迭代求解。

没有成本函数R(q)的迭代解的导数(例如，对于情况β＝0)倾向于PET(“正电子发射断层摄影”)重建的相似导数，参见Parra，L，Barett，HH“List-Mode Likelihood：EM Algorithm and Image Quality EstimationDemonstrated on 2-D PET”，IEEE Trans.Med.Img.，第17卷，第2号，第228-235页，1998。

等式(9)可以认为是密度混合模型，其中P(E_i|x_k)相应于混合模型的密度，而P(x_k|q(x_k，E_i))相应于混合系数。对于混合模型，该数据由无法观察的变量(unobserved variable)z_rk，i扩大(参见Dempster，A，Laird，N和Rubin，D“Maximum likelihood from incomplete data via the EMalgorithm”，J.Roy.Stat.Soc.第B39卷，第1-38页，1977)。对每个库r，由下列等式定义无法观察的变量z_rk，i：

$z_{\mathrm{rk},i}=\left\{\begin{array}{cc}1,& \mathrm{if\; event\; i\; orginated\; in\; bin\; r},k\\ 0,& \mathrm{otherwise}\end{array}\right.---\left(13\right)$

列向量z_rk仅包含一个非零项。z_rk从 $P\left(x_{z_{r,j}}\right|q)={\mathrm{\Σ}}_k{z}_{\mathrm{rk},i}P\left(x_k\right|q)$ 中单独提取。对于“当前计数”的对数似然完整数据是(E₁，z_r1，...，E_N，z_rN，T)，并且似然函数可以写成：

$L(E_1,z_r,\·\·\·,E_N,z_{\mathrm{rN}},T|N,q\left)\right)=\underset{i=1}{\overset{N_r}{\mathrm{\Σ}}}\ln p(E_i,x_{z_{r,i}}|q)+\ln \left(T\right|N,q)$

$=\underset{i=1}{\overset{N_r}{\mathrm{\Σ}}}\ln P\left(x_{z_{r,i}}\right|q)+\underset{i=1}{\overset{N_r}{\mathrm{\Σ}}}\ln p\left(E_i\right|x_{z_{r,i}},q)+\ln p\left(T\right|N,q)$

$=\underset{i=1}{\overset{N_r}{\mathrm{\Σ}}}\underset{k}{\mathrm{\Σ}}{z}_{\mathrm{rk},i}(\ln P\left(x_k\right|q)+\ln p\left(E_i\right|x_k))+\ln p\left(T\right|N,q)$ (14)

等式(14)遵循等式(13)的定义和等式(9)的混合。测量时间T是N_r个间隔时间(interarrival time)的总和，并且遵循由等式(15)给出的Erlang密度p(T|N_r，q)：

$p\left(T\right|N_r,q)=\frac{{\left({\mathrm{\λ}}_{r}T\right)}^{N_r}}{T}{e}^{-\mathrm{\λT}}/(N_r-1)!---\left(15\right)$

其中 ${\mathrm{\λ}}_r={\mathrm{\Σ}}_{i=1}^{N_r}{\mathrm{\Σ}}_{k=1}^K{q}_{r,k}{s}_{k,j,E_i}$ 是泊松分布(Poisson distribution)的平均比率。

在E-步骤中，给出实际测量数据E₁，...，E_N，T和固定的q^(t)，计算L(...|N，q)的估计值，即，在隐变量上限制L(...|N，q)处于极限。简而言之，这可以由下列等式表示：

Q(q|q^(t))＝E[L(E₁，z_r1，...，E_N，z_rN，T|N，q)|(E₁，...，E_N，T|N_r，q^(t))]    (16)

由于等式(14)在z_rk，i上为线性，这相当于由其预期值 z_rk，i代替z_rk，i，其是：

${\stackrel{\‾}{z}}_{\mathrm{rk},i}\left(q^{\left(t\right)}\right)=P\left(x_k\right|E_i,q^{\left(t\right)})=\frac{p(E_i,x_k|q^{\left(t\right)})}{p\left(E_i\right|q^{\left(t\right)})}=\frac{P\left(x_k\right|q^{\left(t\right)})p\left(E_i\right|x_k)}{{\mathrm{\Σ}}_{i=1}^{K}P\left(x_i\right|q^{\left(t\right)})p\left(E_i\right|x_i)}---\left(17\right)$

在此，p(E_i，x_k|q^(t))表示了x_k和E_i的联合概率(密度)。在M-步骤中，需要计算L预期值的最大值，其通过使导数等于零(vanishing derivates)求解而获得

$\frac{\∂Q\left(q\right|q^{\left(t\right)})}{\∂q_{r,l}}=\frac{\∂}{\∂q_{r,l}}\underset{i=1}{\overset{N_r}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\stackrel{\‾}{z}}_{\mathrm{rk},i}\left(q^{\left(t\right)}\right)\ln \frac{q_{r,k}(E_i,x_k)s_{k,i,E_i}}{\underset{m=1}{\overset{K}{\mathrm{\Σ}}}\underset{i=1}{\overset{N_{r,m}}{\mathrm{\Σ}}}{q}_{r,m}(E_i,x_m)s_{m,i,E_i}}$

$+\frac{\∂}{\∂q_{r,l}}\{N_r\ln \left(T\underset{m=1}{\overset{K}{\mathrm{\Σ}}}\underset{i=1}{\overset{N_{r,m}}{\mathrm{\Σ}}}{q}_{r,m}(E_i,x_m)s_{m,i,E_i}\right)-T\underset{m=1}{\overset{K}{\mathrm{\Σ}}}\underset{i=1}{\overset{N_{r,m}}{\mathrm{\Σ}}}{q}_{r,m}(E_i,x_m)s_{m,i,E_i}\}$

$=\underset{i_r=1}{\overset{N_r}{\mathrm{\Σ}}}{z}_{\mathrm{rl},i}\left(q^{\left(t\right)}\right)\frac{1}{q_{r,l}(E_i,x_k)}\frac{\underset{m=1}{\overset{K}{\mathrm{\Σ}}}\underset{i=1}{\overset{N_{r,m}}{\mathrm{\Σ}}}{q}_{r,m}(E_i,x_m)s_{m,i,E_i}-q_{r,l}(E_i,x_l)\underset{i=1}{\overset{N_r}{\mathrm{\Σ}}}{s}_{l,i,E_i}}{\underset{m=1}{\overset{K}{\mathrm{\Σ}}}\underset{i=1}{\overset{N_{r,m}}{\mathrm{\Σ}}}{q}_{r,m}(E_i,x_m)s_{m,i,E_i}}$

$+N_r\frac{\underset{i=1}{\overset{N_{r,i}}{\mathrm{\Σ}}}{s}_{l,i,E_i}}{\underset{m=1}{\overset{K}{\mathrm{\Σ}}}\underset{i=1}{\overset{N_{r,m}}{\mathrm{\Σ}}}{q}_{r,m}(E_i,x_m)s_{m,i,E_i}}-T\underset{i=1}{\overset{N_{r,l}}{\mathrm{\Σ}}}{s}_{l,i,E_i}$

$+N_r\frac{{\mathrm{\Σ}}_{i=1}^{N_{r,i}}{s}_{l,i,E_i}}{{\mathrm{\Σ}}_{m=1}^K{\mathrm{\Σ}}_{i=1}^{N_{r,m}}{q}_{r,m}(E_i,x_m)s_{m,i,E_i}}-T\underset{i=1}{\overset{N_{r,l}}{\mathrm{\Σ}}}{s}_{l,i,E_i}$

$=\underset{i=1}{\overset{N_r}{\mathrm{\Σ}}}\frac{z_{\mathrm{rl},i}\left(q^{\left(t\right)}\right)}{q_{r,k}(E_i,x_k)}-T\underset{i=1}{\overset{N_{r,l}}{\mathrm{\Σ}}}{s}_{l,i,E_i}=0$

(18)

这给出了最大值 $q_{r,k}^{(t+1)}={\left(T{\mathrm{\Σ}}_{i=1}^{N_{r,l}}{s}_{l,m,E_m}\right)}^{-1}{\mathrm{\Σ}}_{i=1}^{N_r}{\stackrel{\‾}{z}}_{\mathrm{rk},1}\left(q^{\left(t\right)}\right).$ 与等式(7)和(16)一起，以定点迭代得出：

$q_{r,k}^{(n+1)}=\frac{1}{T{\mathrm{\Σ}}_{i=1}^{N_{r,k}}{s}_{k,i,E_i}}\underset{i=1}{\overset{N_r}{\mathrm{\Σ}}}\frac{q_{r,k}^{\left(n\right)}{s}_{k,i,E_i}{A}_{\mathrm{att},k}^{\left(i\right)}(1+{\cos }^2\left({\mathrm{\θ}}_i\left(x_k\right)\right))}{{\mathrm{\Σ}}_{m=1}^K{q}_{r,m}^{\left(n\right)}{s}_{m,i,E_i}{A}_{\mathrm{att},m}^{\left(i\right)}(1+{\cos }^2\left({\mathrm{\θ}}_i\left(x_m\right)\right))}$

$=\frac{1}{T{\mathrm{\Σ}}_{i=1}^{N_{r,k}}{s}_{k,i,E_i}}\underset{i=1}{\overset{N_{r,k}}{\mathrm{\Σ}}}\frac{q_{r,k}^{\left(n\right)}{s}_{k,i,E_i}{A}_{\mathrm{att},k}^{\left(i\right)}(1+{\cos }^2\left({\mathrm{\θ}}_i\left(x_k\right)\right))}{{\mathrm{\Σ}}_{m=1}^K{q}_{r,m}^{\left(n\right)}{s}_{m,i,E_i}{A}_{\mathrm{att},m}^{\left(i\right)}(1+{\cos }^2\left({\mathrm{\θ}}_i\left(x_m\right)\right))}$ (19)

对总和从

到∑_i＝1 ^Nr，k的限制，是由于仅来自q_r，k中的光子对总和有贡献的事实。

参考等式(20)、(21)，示出了等式(19)中导数的逐步计算：

$\frac{\∂}{\∂q_{r,l}}\underset{i=1}{\overset{N_r}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\stackrel{\‾}{z}}_{\mathrm{rk},i}\left(q^{\left(t\right)}\right)\ln \frac{q_{r,k}(E_i,x_k)s_{k,i,E_i}}{{\mathrm{\Σ}}_{m=1}^K{\mathrm{\Σ}}_{i=1}^{N_{r,m}}{q}_{r,m}(E_i,x_m)s_{m,i,E_i}}=$

$\underset{i=1}{\overset{N_r}{\mathrm{\Σ}}}\underset{k=1,k=1}{\overset{K}{\mathrm{\Σ}}}{\stackrel{\‾}{z}}_{\mathrm{rk},i}\left(q^{\left(t\right)}\right)\frac{1}{\frac{q_{r,k}(E_i,x_k)s_{k,i,E_i}}{{\mathrm{\Σ}}_{m=1}^K{\mathrm{\Σ}}_{i=1}^{N_{r,m}}{q}_{r,m}(E_i,x_m)s_{m,i,E_i}}}\frac{s_{k,i,E_i}{\mathrm{\Σ}}_{m=1}^K{\mathrm{\Σ}}_{i=1}^{N_{r,m}}{q}_{r,m}(E_i,x_m)s_{m,i,E_i}-q_{r,k}(E_i,x_k)s_{k,i,E_i}{s}_{l,i,E_i}}{{\left({\mathrm{\Σ}}_{m=1}^K{\mathrm{\Σ}}_{i=1}^{N_{r,m}}{q}_{r,m}(E_i,x_m)s_{m,i,E_i}\right)}^2}=$

$\underset{i=1}{\overset{N_r}{\mathrm{\Σ}}}{z}_{\mathrm{rl},i}\left(q^{\left(t\right)}\right)\frac{1}{q_{r,k}(E_i,x_k)}\frac{{\mathrm{\Σ}}_{m=1}^K{\mathrm{\Σ}}_{i=1}^{N_{r,m}}{q}_{r,m}(E_i,x_m)s_{m,i,E_i}-q_{r,l}(E_i,x_i){\mathrm{\Σ}}_{i=1}^{N_{r,m}}{s}_{l,i,E_i}}{{\mathrm{\Σ}}_{m=1}^K{\mathrm{\Σ}}_{i=1}^{N_{r,m}}{q}_{r,m}(E_i,x_m)s_{m,i,E_i}}$

(20)

$\frac{\∂}{\∂q_{r,l}}\{N_r\ln \left(T\underset{m=1}{\overset{K}{\mathrm{\Σ}}}\underset{i=1}{\overset{N_{r,m}}{\mathrm{\Σ}}}{q}_{r,m}(E_i,x_m)s_{m,i,E_i}\right)-T\underset{m=1}{\overset{K}{\mathrm{\Σ}}}\underset{i=1}{\overset{N_{r,m}}{\mathrm{\Σ}}}{q}_{r,m}(E_i,x_m)s_{m,i,E_i}\}=$

$N_r\frac{1}{T\underset{m=1}{\overset{K}{\mathrm{\Σ}}}\underset{i=1}{\overset{N_{r,m}}{\mathrm{\Σ}}}{q}_{r,m}(E_i,x_m)s_{m,i,E_i}}T\underset{i=1}{\overset{N_{r,m}}{\mathrm{\Σ}}}{q}_{r,l}(E_i,x_l)s_{l,i{,E}_i}-T\underset{i=1}{\overset{N_{r,m}}{\mathrm{\Σ}}}{s}_{l,i,E_i}=$

$N_r\frac{\underset{i=1}{\overset{N_{r,m}}{\mathrm{\Σ}}}{s}_{l,i,E_i}}{\underset{m=1}{\overset{K}{\mathrm{\Σ}}}\underset{i=1}{\overset{N_{r,m}}{\mathrm{\Σ}}}{q}_{r,m}(E_i,x_m)s_{m,i,E_i}}-T\underset{i=1}{\overset{N_{r,m}}{\mathrm{\Σ}}}{s}_{l,i,E_i}$

(21)

总而言之，参考使用能量分辨探测器的相干散射计算机断层摄影(CSCT)，由等式(10b)描述了对每个散射平面(q-平面，q(E_i，x_k)＝2E_isin(Θ_i(x_k)/2)/hc)的任意几何形状和轨迹的罚对数似然函数。在等式(10)和(10b)中，k和r分别是q-平面的体素数量和库。

根据本发明，根据等式(10)的函数应当最大化，以便寻找具有最大似然性的图像。由等式(19)给出了对每个q_r，k的优选迭代更新步骤。

图4示出了根据本发明的数据处理设备400的典型实施例，用于执行根据本发明的方法的典型实施例。图4中描述的数据处理设备400包括中央处理单元(CPU)或者连接到存储器402的图像处理器401，存储器402用于存储描述诸如患者或行李的感兴趣对象的图像。数据处理器401可以连接到多个输入/输出网络或诊断设备，例如MR设备或CT设备。而且，数据处理器401可以连接到显示设备403，例如计算机监视器，用于显示在数据处理器401中计算或适用的信息或图像。操作者或用户可以经由键盘404和/或其它在图4中未示出的输出设备与数据处理器401交互。而且，经由总线***405，还可以将图像处理和控制处理器401连接到例如运动监视器，其监视感兴趣对象的运动。如果，例如，对患者的肺成像，运动传感器可以是呼气传感器。如果对心脏成像，运动传感器可以是心电图(ECG)。

本发明可以有利地应用的典型技术领域，包括行李安检、医学应用、材料测试和材料科学。可以实现改进的图像质量和减少的计算量与低劳动强度的组合。而且，本发明可以应用在心脏扫描领域中以探测心脏疾病。

应当注意到，术语“包括”不排除其它元件和步骤，并且“一”、“一个”不排除多个。而且，关于不同实施例描述的元件可以组合。

还应当注意到，权利要求中的参考标记不解释为限制权利要求的范围。

Claims (19)

1、一种用于检查感兴趣对象(107)的计算机断层摄影装置(100)，该计算机断层摄影装置(100)包括

探测元件(123)，适于以能量分辨方式探测从感兴趣对象(107)相干散射的电磁辐射；

确定单元(118)，适于根据对从所述探测元件(123)接收的探测信号进行的统计分析，确定关于所述感兴趣对象(107)的结构信息。

2、根据权利要求1所述的计算机断层摄影装置(100)，其中，所述确定单元(118)适于根据对探测信号的位置和能量进行的分析，确定关于所述感兴趣对象(107)的结构信息。

3、根据权利要求1所述的计算机断层摄影装置(100)，其中，所述确定单元(118)适于根据对电磁辐射的每个量子的探测信号单独进行的分析，确定关于所述感兴趣对象(107)的结构信息。

4、根据权利要求1所述的计算机断层摄影装置(100)，其中，所述确定单元(118)适于根据对探测信号进行的极大似然分析，确定关于所述感兴趣对象(107)的结构信息。

5、根据权利要求1所述的计算机断层摄影装置(100)，其中，所述确定单元(118)适于根据对探测信号进行的对数似然分析，确定关于所述感兴趣对象(107)的结构信息。

6、根据权利要求1所述的计算机断层摄影装置(100)，其中，所述确定单元(118)适于在考虑了吸收系数的实验和/或理论预定值的情况下，确定关于所述感兴趣对象(107)的结构信息。

7、根据权利要求1所述的计算机断层摄影装置(100)，其中，所述确定单元(118)适于在考虑了电磁辐射源(104)和所述探测元件(123)之间的电磁辐射的衰减的能量相关性的情况下，确定关于所述感兴趣对象(107)的结构信息。

8、根据权利要求1所述的计算机断层摄影装置(100)，其中，所述确定单元(118)适于在考虑了所述探测元件(123)的接受度的能量相关性的情况下，确定关于所述感兴趣对象(107)的结构信息。

9、根据权利要求1所述的计算机断层摄影装置(100)，其中，所述确定单元(118)适于根据对探测信号进行的迭代分析，确定关于所述感兴趣对象(107)的结构信息。

10、根据权利要求1所述的计算机断层摄影装置(100)，其中，所述确定单元(118)适于根据对动量传递的不同值单独执行的探测信号的分析，确定关于所述感兴趣对象(107)的结构信息。

11、根据权利要求1所述的计算机断层摄影装置(100)，其中，所述探测元件(123)适于探测预定数量的电磁辐射的量子，或者探测预定时间间隔内的电磁辐射的量子，用于统计分析。

12、根据权利要求1所述的计算机断层摄影装置(100)，适于为相干散射计算机断层摄影装置。

13、根据权利要求1所述的计算机断层摄影装置(100)，包括电磁辐射源(104)，其适于向感兴趣对象(107)发射电磁辐射，并且所述计算机断层摄影装置包括准直器(105)，其布置在所述电磁辐射源(104)和所述探测元件(123)之间，所述准直器(105)适于准直由所述电磁辐射源(104)发射的电磁辐射束，使其形成扇形束或锥形束。

14、根据权利要求1所述的计算机断层摄影装置(100)，其中，所述探测元件(123)形成单层探测器阵列。

15、根据权利要求1所述的计算机断层摄影装置(100)，其中，所述探测元件(123)形成多层探测器阵列(108)。

16、根据权利要求1所述的计算机断层摄影装置(100)，配置为以下组中的一个，该组包括行李安检装置、医疗应用装置、材料测试装置和材料科学分析装置。

17、一种使用计算机断层摄影装置(100)检查感兴趣对象(107)的方法，该方法包括下列步骤：

以能量分辨方式探测从感兴趣对象(107)相干散射的电磁辐射；

根据对所探测的信号进行的统计分析，确定关于所述感兴趣对象(107)的结构信息。

18、一种计算机可读介质(402)，其中存储了使用计算机断层摄影装置(100)检查感兴趣对象(107)的计算机程序，其当由处理器(401)执行时适于执行下列步骤：

以能量分辨方式探测从感兴趣对象(107)相干散射的电磁辐射；

根据对所探测的信号进行的统计分析，确定关于所述感兴趣对象(107)的结构信息。

19、一种检查感兴趣对象(107)的程序单元，其当由处理器(401)执行时，适于执行下列步骤：

以能量分辨方式探测从感兴趣对象(107)相干散射的电磁辐射；

根据对所探测的信号进行的统计分析，确定关于所述感兴趣对象(107)的结构信息。

Images