CN105580054A - 电子密度图像的联合重建 - Google Patents

Abstract

一种用于重建对象PAT中的电子密度的图像的方法和相关系统(IMA)。X射线辐射成像器被用于将所述对象PAT暴露于辐射以接收投影数据。重建方法将来自两个通道的投影数据进行组合，即将基于康普顿散射的衰减数据p_C和相移数据进行组合。

Description

电子密度图像的联合重建

技术领域

本发明涉及用于生成图像的方法、图像处理系统、计算机程序单元和计算机可读介质。

背景技术

传统X射线成像关注于作为提供在各种类型的人类组织之间的对比的物理量的衰减。最近，如双能量成像或相位对比成像的其他成像技术提供与如被成像组织的电子密度、折射率递减或散射横截面的不同物理量有关的信息。似乎在基于衰减的成像或基于折射率成像中的任一个中图像对比(即，在图像元素之间的相对差异)能够更为显著，并且似乎不能够做出关于这两个成像技术中的一个是更好的综述。

发明内容

因此，可能存在对用于图像生成的方法和相关装置或系统的需要以鲁棒性地改善在广泛的应用中的图像质量。

本发明的目的通过独立权利要求的主题来解决，其中，在从属权利要求中包含另外的实施例。应当注意，本发明的以下描述的方面同样适用于图像处理系统、计算机程序单元和计算机可读介质。

根据本发明的第一方面，提供了一种用于生成对象中的电子密度分布的图像的方法，其包括以下步骤：

接收：

i)针对至少一个能量值，如由辐射敏感探测器探测的相位对比正弦图其中，涉及在辐射通过所述对象之后由所述辐射经历的相移的相移梯度；以及

ii)如由所述辐射敏感探测器探测的或者不同辐射敏感探测器探测的康普顿衰减正弦图p_C，其中，p_C实质上记录到对由所述辐射在所述的通过所述对象中经历的总体衰减的康普顿散射贡献；

基于两个正弦图和p_C两者，联合重建所述电子密度的横截面图像。

在一个实施例中，联合重建基于(统计)迭代重建算法，即从针对电子密度图像的初始估计到针对电子密度图像的一系列后续“经更新的”估计迭代地进行的算法，并且本文提出将康普顿正弦图p_C和微分相位对比正弦图提供到所述(单一或共同)迭代重建算法中。这允许将来自两种不同模态或“通道”(在一方面谱或双能量(其被用于导出康普顿正弦图)以及在另一方面相位对比，其被用于导出相位对比正弦图)的与电子密度有关的投影信息进行组合，以由此更好地管理否则可能损坏图像质量的不利噪声行为。

根据一个实施例，“联合重建”尤其包括在每个迭代步骤中(或者至少在多个迭代步骤中)生成针对两个正弦图(p_C和)的估计，并将其与实际测量的正弦图p_C和进行比较。根据关于针对p_C和的所估计的正弦图两者的“失配”导出针对电子密度图像估计的每次更新。再换言之，正弦图p_C和两者被用作“参考数据”以测量迭代是否收敛于针对电子密度的估计，当被向前投影时，其足够程度地接近实际测量的正弦图p_C和

因为在提出的方法中，来自两种模式的数据被组合成单一重建算法，有用的平衡可以被打破以由此更好地处理或对付两个相反的噪声效应：(i)在双能量材料分解中在更高的空间频率中的主导图像噪声，以及(ii)由相位正弦图的微分属性引起的在微分相位对比成像通道中在更低的空间频率中的主导图像噪声。另外，在提出的方法中，由于关于缺失数据(出于相位缠绕、相干损失等的原因)将来自两个通道(p_C和)的数据合并到(联合)迭代重建中是更鲁棒的，这是因为其他通道仍然可以保持足够的信息来实现具有令人满意的保真度的重建。本文中不一定必需由相同的探测器采集两个正弦图p_C和并且本文设想由专用探测器模块采集每个通道的技术，比如其中的各自的曝光由所使用的X射线源的焦点移动来控制。

根据一个实施例，联合重建中的迭代由成本函数来驱动，所述成本函数包括对如在两个正弦图中记录的噪声进行建模的噪声项。正确地考虑在两个正弦图中的单个噪声水平和/或行为。不需要对图像域中的噪声进行建模，因为相反考虑针对所测量的数据的噪声估计，其更易于实施，因为噪声在图像域中被认为是白噪声。换句话说，重建基本是从投影空间(其中定义正弦图)到图像空间(即，由在扫描器的在成像期间对象驻留在其中并且将要被重建的电子密度图像将要被定义在其中的膛中的点定义的3D空间)的变换。考虑像这样的转换的重建算法，所提出的方法之后能够被认为将成本函数中的噪声项用于模型，或捕捉已经在投影空间中但不在图像空间中的噪声行为。

成本函数被配置为考虑两个信号通道p_C和换言之，成本函数是针对两个通道的共同成本函数。比如，这通过在成本函数中并且针对重建使用具有针对康普顿通道的专用部件(其是线积分的传统离散实施方式)和针对微分相位对比通道的另一专用部件(其是微分线积分的离散实施方式)的(“组合”)系统矩阵(即，向前投影运算符的表示)来实现。成本函数的噪声项被用于对在实际测量的正弦图与在迭代重建的过程期间生成的电子密度估计的投影之间的偏差进行加权。

在一个实施例中，噪声项中的噪声模型通过协方差矩阵来参数化。在一个范例性实施例中，协方差矩阵是对角形式的。在该实施例中，协方差矩阵忽略在两个通道之间的任何可能的协方差。这对良好近似是有效的，但是所提出的发明也能够被实施具有含有至少一个非零非对角元素或者尤其是全协方差矩阵的协方差矩阵。

根据一个实施例，重建合并共轭梯度方法来解决电子密度图像，尤其是使所述成本函数最小化。

根据一个实施例，使用多能量扫描器，尤其是双能量扫描器或具有多层的扫描器，尤其是双层探测器，或者具有光子计数探测器的谱成像扫描器来探测p_C正弦图和/或正弦图。

根据一个实施例，扫描器包括用于采集相位对比正弦图的微分相位对比成像装备。在一个实施例中，所述装备包括干涉仪，所述干涉仪包括多个光栅结构。

在一个实施例中，成本函数基于最大似然方法，并且迭代重建的更新函数通过被应用到成本函数的数值最小化算法来导出。

对象可以是人类或动物患者或其相关部分(感兴趣区域)。所述方法和相关系统也可以在其中对象是诸如在安全检查中的行李件或在非破坏性材料测试中的有机或无机样本的非医学情境中使用。

附图说明

现在将参考以下附图来描述本发明的范例性实施例，其中：

图1示出了成像布置；

图2示出了如在根据图1的布置中使用的相位对比成像装备；

图3示出了生成图像的方法的流程图。

具体实施方式

参考图1，示出了根据一个实施例的成像系统。广泛地，所述成像布置包括基于X射线的CT扫描器IMA和用于控制所述扫描器IMA的操作的工作站或操作控制台CON。

通用计算系统可以充当操作员控制台CON，并且包括输出设备(例如显示器M)和输入设备(例如键盘、鼠标等)。驻留在控制台CON上的软件允许操作员控制系统IMA的操作，例如允许操作员通过选择预生成的成像协议直接地或间接地选择成像参数。显示单元(例如监控器/屏幕)M通信耦合到控制台CON以辅助输入控制信息或观察扫描器的操作状态或观察由扫描器供应的图像或通过图像处理如此供应的图像来获得图像的操作状态。也具有图像处理模块IP，在一个实施例中其通信耦合到控制台CON以实现图像处理。广泛地，图像处理器IP包括输入IN和输出OUT接口，并且提供图像重建功能。下面在图3中将进一步更详细地解释图像处理模块IP的操作。

通过可视化器件VS来图像绘制经重建的图像以供显示在屏幕M上，可视化器件VS被配置为接收图像数据和可视化信息来驱动控制台CON的视频卡以由此实现在屏幕M上的显示。可视化信息可以包括颜色或灰度面板，在所述颜色或灰度面板上通过映射器来映射图像值以实现以在期望范围中的期望颜色编码或灰度编码的显示。

扫描器IMA包括固定机架和由固定机架可旋转地支撑的旋转机架G。旋转机架104绕检查区域关于纵轴或z轴旋转。检查区域被形成作为旋转机架G中的开口或膛。

诸如卧榻的对象支撑物B将对象PAT或目标支撑在检查区域中，并且能够被用于在扫描之前、期间和/或之后相对于x轴、y轴和/或z轴定位对象或目标。

诸如X射线管(“管”)的辐射源XR由旋转机架104支撑，并且随旋转机架G绕检查区域旋转，并且经由焦点放射贯穿检查区域的辐射。在一个实施例中，该源XR能够在至少两个不同电压处操作以由此提供双能量模式功能。

辐射射束以膛的等中心为中心，并且定义(扫描器的)大体圆形视场FoV以重建贯穿重建平面，其大体垂直于射束的中心射线并且其延伸通过等中心。

(在适当的外壳中的)辐射敏感探测器阵列D被定位为跨检查区域与辐射源XR相对。探测器阵列D包括一行或多行探测器像素px，其探测贯穿检查区域的辐射并生成指示探测到的辐射的电流或电压信号。可旋转机架的旋转(以及由此至少X射线源的焦点的旋转，在第四代扫描器中仅仅管在机架中旋转，其中探测器阵列被固定地布置以将机架内部的全部排成行以使得探测器的至少部分能够接收针对任何旋转角度的辐射)由一个或多个控制器和/或一个或多个驱动系统(例如，电机、耦合等)来实现。

在当扫描器处于使用中时的成像运行期间，在对象或目标的扫描期间焦点以给定角度频率围绕膛B以及因此围绕在其中的对象或目标PAT在路径上(大体圆弧或全圆周旋转)的预定角度范围上旋转。针对每个旋转角度探测器探测由焦点以所述角度通过对象或目标之后发射的辐射。当在通过对象或目标期间被认为是波而不是射线时，辐射经历衰减和辐射的一定相移。衰减一般与对象或目标的局部密度成比例，并且相移与材料的局部折射率有关。每个探测器像素(相对的焦点)接收一定量的光子能量，并且如以上简单提到的通过发出对应的电信号来做出响应。然后通过数据采集系统DAS(未示出)来将所述信号的收集转化成相应数字值。所述值对关于由X射线射束经历的衰减量和相移的信息进行编码，并且如下面将更详细解释的能够通过使用干涉仪设置来提取所述信息。衰减和相移两者能够由辐射淹没的组织PAT的局部复数折射率n(x,y,z)(即是复数)来描述，复数折射率的实部描述相移并且虚部描述衰减。然后数字值的收集被整合为数字值阵列，其形成针对给定采集时间和投影角度θ的图像原始数据的1D或2DX射线投影图像。探测器由此输出原始投影图像的收集，针对每个旋转角度θ一个或多个，以由此定义原始图像数据正弦图，即在其中对于每个投影方向/旋转角度θ存在如由探测器D在相应角度θ处探测到的相关联的原始图像数据的数据结构。在一个实施例中，探测器是光子计数类型的。在备选实施例中，探测器是能量积分类型的。

当探测器D是光子计数类型时，成像器的DAS包括能量分辨电路，其可以包括脉冲整形器，所述脉冲整形器被配置为处理由探测器阵列D输出的电子信号，并且生成指示探测到的光子的能量的脉冲，例如电压或其他脉冲。能量鉴别器鉴别能量或脉冲(其是针对探测到的X射线光子的能量的度量)。在一个实施例中，能量鉴别器包括至少一个比较器，其被配置为将脉冲的幅度与对应于感兴趣能量的至少一个能量阈值进行比较。比较器产生指示探测到的光子的能量是否高于或低于阈值的输出信号。计数器增加针对每个阈值的计数值。比如，当针对特定阈值的比较器的输出指示脉冲的幅度超过对应的阈值时，增加针对该阈值的计数值。能量分箱器(binner)库将所计数的脉冲分配到能量箱(bin)，所述能量箱对应于不同的能量范围。例如，可以针对两个阈值之间的能量范围定义箱，其中得到针对较低阈值但不是针对较高阈值的计数的光子将被分配给该箱。换言之，当使用光子计数探测器时，存在由DAS拾取的正弦图的收集，针对每个能量箱或水平一个正弦图。

在成像系统IMA包括(一行或多行)仅仅能量积分探测器像素px的情况下，能够省略能量分辨电路，并且通过使用不同X射线输入谱采集数据来获得能量信息。然而，在优选实施例中，成像器具有光子计数类型的探测器或双层类型的探测器，或者根据任何其他的具有提供能量分辨能力的适当技术的探测器。

如图1进一步示出的，成像扫描器CT包括相位对比成像装备PCI，其允许从记录的投影图像提取相移梯度。如简单提到的，这是因为在由X射线在其通过物质PAT中遭受的衰减上，也存在由介入物质PAT引起的相移

参考图2简单地描述了相位对比成像装备PCI的基础部件，实质上是基于光栅的干涉仪。

在图2中，描绘具有相邻布置的源光栅元件202的X射线源XR。穿透源光栅元件202的X射线辐射114可以被认为是至少部分地空间上相干的。X射线辐射XR包括单个波前，其中的波前210a被描绘为在目标108的穿透之前的波前，而波前210b被描绘为在目标PAT的穿透之后的具有图示的相移的波前。

被布置在目标PAT之后并且与探测器元件104和分析器光栅G₂206两者隔开距离d的是具有间距p的相位光栅204。具有单个探测器像素元素116的探测器D随后地探测X射线辐射XR的干涉图样，其由相位光栅G₁204和分析器光栅G₂206施加。

致动器元件208被示意性地图示为适于相对于具有源光栅202、相位光栅G₁204和X射线探测器D的X射线管XR侧向地移位分析器光栅G₂206(具有间距q)。然而，致动器元件208可以被布置在针对位移x的光栅元件202、204、206中的任何处。所述致动器被用于实现被认为是“相位步进”的功能。在相位步进中，将被用于生成至少部分地空间相干的X射线辐射的源光栅、相位光栅和分析器光栅中的一个相对于其他光栅和X射线探测器元件侧向地移位其各自的光栅间距的分数，例如，组成相位步进状态的各自的光栅元件的光栅间距的四分之一、六分之一或八分之一。图像采集和侧向移位被重复例如四次、六次或八次以用于采集多个相位对比投影，从而一起组成引起能够被分析以还原衰减和相移的变化强度I的信号的相位步进系列。

X射线射束XB的波前的相移可以被认为是直接与沿着X射线辐射的路径的电子密度的积分有关。基于光栅的微分相位对比成像可以允许采用例如具有ΔE/E～10％的相对宽带的X射线源。换言之，多色X射线辐射而不是实质上具有仅仅单一波长的单色X射线辐射可以被采用于相位对比成像。

应当理解，图1以高度示意性方式示出了如插设在X射线源与探测器D之间的PCI的布置。在一个实施例中，PCI光栅G1、G2被集成在探测器D的外壳中，具体地在一个实施例中被安装在探测器D上，并且源光栅元件202被集成在X射线源XR的出口部分。

感兴趣的观察是以下两者，i)诸如双能量成像(或者更一般地，基于光子计数的谱成像技术)的能量分辨成像和ii)相位对比成像，提供与样本PAT的电子密度有关的每个信息。在双能量或光子计数成像中，将衰减分解成光电效应和康普顿效应实际上提供样本PAT的电子密度图。同样地，在相位对比成像中，由样本PAT引起的折射率递减与在X射线波长方案中的电子密度严格地成比例，并且如以上关于图2所解释的，能够通过在基于光栅的相位对比成像设置中的相位步进获得。换言之，基于康普顿效应的衰减和在相位对比成像中的相移记录提供针对电子密度的分布的信息的两个单独通道，并且本文中提出了图像处理器IP，其利用该明显的“冗余”来以使得解决由双重的不良噪声行为提出的挑战的方式将来自两个通道或两种模态的信息进行组合。首先，图像噪声已经被观察以在至少双能量材料分解中增加，其在临床实践中可以解释该方法的相对低摄取。第二，与其中噪声在更高频率中增加的衰减成像的情况相反，当将计算机断层摄影重建算法应用到微分相位对比成像正弦图时，噪声功率已经被观察以朝更低频率增加。

如现在将更详细解释的，如本文提出的图像处理器IP在针对两种模态的“逐联合”重建算法中通过将来自两种模态(谱/双能量或否则能量分辨和相位对比)的与电子密度有关的信息进行组合来解决这种噪声行为挑战。基于由康普顿衰减正弦图和一个或多个相移正弦图提供的投影影像，如本文提出的重建算法允许重建以生成针对每个轴向位置z的样本PAT中的电子密度的横截面图像(“切片(图像)”)。然后能够由适当的绘制器模块将横轴切片图像的收集组合成指示对象PAT的密度的体积图像数据。

在下文中，更详细地解释用于对物理交互作用进行建模的途径，所述物理交互作用引起在探测器D处的记录影像并且其构成由图像处理器IP运行的提出的算法的基础。

初始地，针对每个投影角度θ，在相位步进期间，由探测器像素px拾取辐射信号强度I，并且所述信号在所述探测器像素/通道px处振荡，其中相位步进x扫描近似如下：

其中，A涉及射线的总体线性衰减，涉及波前201b的局部相移梯度，V(“可见性”)是小角度散射，q是光栅间距，并且x是在相位步进期间光栅G1和G2的相对移位。

振荡信号I由适当的采样和保持电路捕获为曲线，并且能够在“相位复原步骤”中通过在所捕获的信号I曲线上运行快速傅里叶变换(FFT)从所述曲线中提取值A、V、然后，傅里叶分解得到A作为DC分量，其中频率分量ν＝1/q的幅度是针对可见性V的度量，并且是频率分量ν＝1/q的相位。

在本文中应当理解，按照(1)的近似一般针对每个像素px是不同的，因为每个像素px一般将看到针对每个x的强度的不同波动，因此值A、V、将需要针对每个像素px单独地进行拟合。波前的相位梯度的投影与通过关系的强度曲线I的度量相位值有关，并且能够从其中还原，其中，x指示相位步进方向，并且d是在光栅G₁与G₂之间的距离，并且λ是近似一般多色辐射的能量谱的设计能量E_d的波长。例如参见TWeitkamp等人的论文(“OpticsExpress”，第13卷，第16号，第6296-6304页(2005))中的方程(1)。

在另一实施例中，不通过其中光栅中的一个被移动的相位步进途径进行相位复原，而是通过傅里叶方法来进行相位复原。在由AMomose等人的文章“Four-dimensionalX-rayphasetomographywithTalbotinterferometryandwhitesynchrotronradiation：dynamicobservationofalivingworm”(OPTICSEXPRESS，第19卷，第9号，第8423页，2011年4月)中更详细地描述了该傅里叶方法。

涉及来自在不同网格位置x处的强度测量的光电效应p_p和康普顿散射p_C对线性衰减系数和相位φ(E_d)的贡献的物理模型的数学表达是：

其中

∫μ(E,l)dl＝p_p(E/E_d)^-3+p_Cf_KN(E)(3)

其中：

E是能量范围或值；

l是通过对象PAT的组织中的路径长度；

I是在成像运行中在网格位置x处记录的强度，然而，I₀是在“空气”或空白扫描中在像素处记录的“参考”强度；

V是可见性，即干涉图样的调制深度。

E_d是设计(或更一般地说参考)能量。这是在图2的干涉仪响应以最佳可见性V处的能量。换言之，设计能量是干涉图样的观察到的调制深度在其最大值处的能量。该调制深度在空白扫描期间被探查并且由在空白扫描期间采集的数据的拟合值V来量化。

x＝网格位置；

p_p、p_C＝光电效应和康普顿散射对射束的衰减的各自的贡献；

＝由目标PAT引起的干涉图样并且经由如按照图2的干涉仪布置在探测器处记录的相移，波的相移梯度ΔΦ可从其中还原，如以上关于Weitkamp论文所提到的；

f_KN＝描述取决于在即将到来的能量范围中的康普顿散射横截面的能量相关性的Klein-Nishina公式。

模型方程(2)、(3)能够被理解为方程(1)的改写，其中，现在在如在探测器D处登记的经衰减的X射线射束的能量谱上对针对每个能量的相应网格位置x处看到强度I(.)求积分。

方程(3)实质上定义总体线性衰减到单独的光电贡献和康普顿散射贡献中的谱分解。也能够看到，由成像器IMA在至少两个能量水平处的图像采集将允许解决各自的衰减贡献。这是“双能量途径”并且在RobertAlvarez的文章“Energy-selectivereconstructionsinXraycomputerizedTomography”(Phys.Med.Biol.，1976，第21卷，第5号，第733页)中进行描述。

在下文中，参考图3中的流程图，描述了由图像处理器IP实施的方法的步骤。方法允许根据投影数据来重建横截面电子密度图像，即康普顿衰减正弦图p_C和相移正弦图

在步骤S305处，成像器IMA用于通过使用成像器IMA的PCI装备来采集针对至少两个能量水平E1、E2的目标PAT的原始图像数据(或者目标PAT的至少部分的原始图像数据)以执行如以上早先关于图2概述的相位步进。

原始图像数据表示在采集期间辐射波210b在其通过对象PAT中经历。的总体变化，尤其是总体折射(即，辐射波前的相移和通过散射和吸收的衰减)。在一个实施例中，通过在双能量模式中操作成像器的管XR来采集针对每个能量水平E1、E2的各自的谱上不同的图像。即，在不同电压(kV)处采集原始图像数据。在其他实施例中，代替双能量，也设想其他谱采集模式，例如使用双层探测器，或者探测器在具有对应于至少两个能量水平E1、E2的箱的光子计数模式中操作。

在步骤310a、b处，在针对每个能量水平E1、E2的各自的原始图像数据上执行相位复原以获得：

以衰减正弦图的形式的衰减的投影，以及

以各自的微分相位正弦图的形式的梯度相移的投影。微分相位正弦图将在每个条目处对各自的投影(即，沿着各自的射线的线积分)编码成相位梯度，即辐射波前在通过对象PAT组织中经历的相位的变化速率。

相位复原可以包括不仅单独地针对每个投影方向的相位步进(如在由Weitkamp描述的经典相位步进中进行的)而且跨针对多个投影角度的多个条目的相位步进。换言之，相位步进是“交错的”，即不局限于仅仅逐行或逐列地构建强度信息，而是比如由IZanette在“Interlacedphasesteppinginphase-contrastXraytomography”(APPLIEDPHYSICSLETTERS98，094101，(2011))中描述的可以跨投影角度对比光栅位置平面构建。交错方法的延伸是“滑动窗口”方法，本文也设想“滑动窗口”方法并且由IZanette在“Trimodallow-doseXraytomography”(PNAS，第19卷，第26号，第10199-10204页，2012年6月26日)中进行描述。在滑动窗口方法中，不仅仅存在跨不同投影角度的相位步进，而且在投影角度对比光栅位置平面中的特定点被重新访问以改进计算效率。

2个或更多个相位对比正弦图(即，每个能量通道E1、E2一个相位对比正弦图)能够被折叠成一个单一“谱”正弦图(各自的能量相关条目简单地依次彼此邻接以由此扩大正弦图向量)，但是这可以不是如此必需的，并且在一些实施例中，相反针对特定能量选择单一正弦图，并且该选择的正弦图正是在下面的步骤中处理的正弦图。在下文中，术语相位正弦图因此包括具有针对单一能量的条目或根据2个或更多个能量的条目的正弦图。

在步骤S315处，谱分解方案被应用到(步骤310a、b的)衰减正弦图数据，所述衰减正弦图数据表示总体衰减贡献以由此根据方程(3)将所述贡献分解成康普顿散射贡献和光电贡献。因为在步骤S305处的采集得到2个总体衰减正弦图(针对每个能量E1、E2一个)，所以现在能够对方程(3)进行求解以得到康普顿贡献。。本文不再考虑并且可以丢弃光电贡献，但是康普顿散射贡献被保留并且被整合成(单一)“康普顿正弦图”p_C，其在每个条目处对由实质上仅仅康普顿散射引起的衰减的相应投影(即，沿着相应射线的线积分)进行编码。换言之，在步骤S315处供应的输出是康普顿正弦图p_C。因为在电子密度与线性衰减之间的先前提到的比例性仅仅严格应用到单独由康普顿散射引起的衰减的那部分，所以该分解允许增加下面计算的准确度。

在步骤S320a、b处，各自的能量相关的微分相位正弦图被重新缩放以获得电子密度的相应投影。在步骤S325处，相似的重新缩放被应用到康普顿正弦图p_C的条目。现在同样地，经如此重新缩放的康普顿正弦图p_C表示电子密度的投影，但是现在根据不同数据集来导出，即康普顿正弦图p_C。如先前针对正弦图预缩放所使用的，下文经重新缩放的正弦图p_C将由相同的符号p_C、指代。能够以已知方式通过使用具有已知电子密度(例如水)的体模的扫描的校准过程来实现重新缩放。

在步骤S330处，然后在图像处理器IP的输入端IN处接收相应康普顿正弦图p_C和相移正弦图输入端IN可以是多通道(针对每个能量水平一个通道，和/或分别针对p_C和一个通道)，一个或输入端IN可以被布置为通过其接收到所有正弦图的单一通道输入端。本文设想，比如在DAS中由预处理器执行或由集成在图像处理器IP中的预处理器执行先前步骤S310-S325。另外可以在并行流水线中执行或相反可以在串行流水线中执行步骤对S310a、b和/或步骤S320a、b、S325。

在步骤S335处，如现在将更详细解释的，然后通过使用正弦图来联合重建电子密度图像。

在优选实施例中，本文设想重建基于迭代重建方案。广泛地，在迭代图像重建中，感兴趣图像(其是在图像空间/域中的目标)不是如在经典的(滤波)反投影(FBP)算法中直接被重建，而是相反算法在图像与投影域之间的许多迭代中循环，在经典的(滤波)反投影(FBP)中正弦图(其是投影空间/域的目标)被直接转变成图像域中的相应图像元素。更具体地，在本情境中，迭代算法以针对电子密度图像的当前估计开始。然后通过向前投影器的数值实施方式来对该所估计的电子密度图像进行向前投影，所述向前投影器将所估计的电子密度图像从图像域带到投影域，以由此获得当前估计的电子密度图像的“合成”或估计的正弦图。在每个迭代步骤中，然后将所估计的正弦图与实际测量的正弦图p_C、进行比较，并且然后通过成本函数C来评估在所估计的正弦图与实际测量的正弦图p_C、之间的偏差(比如，通过获取逐像素差)。成本函数将“成本”(即，实数)赋给当前偏差。然后目的是找到使成本最小化的针对电子密度x的估计。选择更新函数，在一个实施例中其由迭代方案给出，所述迭代方案迭代地输出已知收敛到成本函数的最小元的估计的序列。在一个实施例中，使用最大似然(ML)途径，并且成本函数能够被表示为似然函数的负对数。然后针对最小成本的要求对应于找到负对数似然的最小元。数值算法被应用以解决所述最小元，其继而得到更新函数。如果估计返回改进成本少于预设成本阈值的成本的更新或如果达到迭代的最大数量，则估计被接受作为针对电子密度的最终估计。针对电子密度图像的初始估计可以是非常简单的一个，例如空白图像或粗略地对被成像的预期结构进行编码的图像。

从概率性/统计性角度看，正弦图能够被认为是遭受(白)噪声的随机事件的样本，每个事件是在相应探测器像素处的探测(在所测量的正弦图中的条目)。如果潜在的概率分布被假设是高斯的(或针对光子计数实施例的Poison)，则每个探测器像素能够与针对所述局部分布的均值和方差的高斯参数(μ_G，σ²)相关联。在ML中，然后能够设置与“解释”给出特定(但又未知的)图像x(即，电子密度的空间分布，其继而在图像域中分别与μ(x,y,z)的康普顿部分和μ(x,y,z)的实部有关)的事件b(即，正弦图采样p_C和)的似然性有关的似然函数L(x|b)。在最大似然(ML)途径中，然后可以询问哪个电子密度x将使似然函数最大化的问题。假设针对测量结果的高斯分布并且取L的对数，然后得到能够与正则项组合的“数据项”以得到具有以下结构的针对迭代重建的成本函数：

C＝||W^-1(Ax-b)||²+βR(x)

(x)

其中，平方表达表示数据项，并且βR(，)表示正则项。

更具体地，成本函数C包括以下元素，其中：

A是系统矩阵，即，A是向前投影算子的离散形式；

x是要被重建的未知图像；

b是所测量的数据，即两个正弦图(p_C、)的收集；

W是所测量的数据的协方差矩阵的平方根(由于协方差矩阵是正定性的，所以其被很好地定义)，其量化在投影域中的噪声；

||.||是适当的范数(例如欧几里得距离)以量化图像x与实际测量的数据b的误匹配，

R(.)是粗糙惩罚(例如被应用到当前电子密度的局部差异以量化在图像x中的粗糙度或噪声含量的Huber惩罚)以迫使针对电子密度图像x的估计的特定平滑度；

β是正则化参数，以迫使针对电子密度的估计的特定“形状”。一般地，beta是允许对惩罚项R的响应进行“微调”的正数。比如，选择大beta是强调粗糙度惩罚响应，然而选择小β是缓和惩罚响应。

如早前所提到的，如本文提到的算法使用图像x的投影，即电子密度的空间分布的投影与微分相位对比信号相关而且与康普顿正弦图p_C＝b_C(在谱分解之后)线性相关的事实。因此本文提出，通过定义以上成本函数C一次或同时考虑两个通道，即，系统矩阵A具有针对康普顿通道的分量，其是线积分的传统离散实施方式，并且其具有针对微分相位对比通道的另一分量，其是微分相位线积分的离散实施方式。换言之，如本文提出的，被反馈到成本函数C中用于评估的实际测量的正弦图数据b包括康普顿线积分以及微分相位对比线积分两者。再换言之，在每个迭代步骤中通过同一个成本函数来评估康普顿散射引起的入射波XB的衰减数据以及引发的相位变化两者。迭代重建在这样的意义上“联合”进行，即，在每个迭代中，生成两个正弦图(b_C和b_p)的估计，并且将其与所测量的数据b进行比较，并且然后根据两个正弦图中的“失配”来导出针对图像x的成本。

根据一个实施例，比如由T等人在“IterativeReconstructionforDifferentialPhaseContrastImagingusingsphericallysymmetricbasisfunctions”(MedPhys38(8)，2011，第4542页)一文中描述了向前投影算子A的微分相位线积分分量A_p的离散实施方式，同时理解，本文同样设想其他实施方式。系统矩阵A的每个条目a_ij分量是向前投影算子的离散化实施方式，其量化图像元素j对具有指数i的所测量的像素具有什么影响。针对该投影操作，将要被投影的图像x被认为由图像元素组成。更具体地，图像被表示为一族基函数的线性组合，例如常见棱镜形状的“体素”。然后相应正弦图被“解释”为在投影域中的那些基函数的投影(“足迹”)的线性组合。根据一个实施例，如本文使用的针对系统矩阵分量A_p的实施方式的基函数不是由体素而是由一族Kaiser-Bessel函数(Bessel函数的特殊变型，也称为“blob”)形成。如在由等人的先前提到的论文中观察到的，使用blob作为基函数来实施系统矩阵A的至少A_p部分允许确保许多计算优点。比如，如在正弦图数据b_p中捕获的关于微分数据的投影算子涉及关于相应正弦图数据点的极坐标参数化的距离坐标的偏微分。在迭代重建的向前投影步骤中使用blob作为基函数允许避免所述微分操作的数值实施方式，其引起病态行为，即数值实施方式关于离散误差是高度敏感的。当投影到康普顿分量时，也能够应用基于blob的实施方式以使向前算子离散化。

本文应当理解，为简洁起见以矩阵形式书写的依据表达式(4)的成本函数C。比如，术语b、x，尽管一般表示为矩阵，但是在本文中被认为是“向量”，其中每个条目是多维的以实现依据(4)的特别简洁的符号。应用矩阵乘法，成本函数能够被扩展成项的总和，所述项的总和包括表示i)在图像x的向前投影与所测量的相位对比正弦图(如果相位对比正弦图是谱的，则存在针对两个或更多个能量值E1、E2中的每个的相应偏差贡献项)之间的相应偏差或误匹配贡献，和ii)在图像x的向前投影与所测量的康普顿正弦图之间的相应偏差或误匹配贡献。将所有误匹配贡献整合成单一函数C，从而允许驱动或引导迭代以便得到针对密度图像的估计，其投影近似针对p_C、的相应正弦图估计两者。

如在表达式(4)中能够看到的，成本函数C惩罚(即，引发更高成本)在实际测量的正弦图b与如通过将系统矩阵A应用到图像x合成的所估计的正弦图Ax之间的大偏差(如下面更详细解释的针对普遍的噪声进行加权)。除了该偏差之外，同样地在成本函数中考虑在探测器处的噪声行为，并且将所述噪声行为封装在矩阵表达式W^-1中。噪声项中的噪声矩阵对在实际测量的正弦图b与所估计的正弦图Ax之间的偏差进行加权，以由此更好地使成本函数的返回与用于测量数据b的特定探测器的特性一致。噪声由探测器量子噪声和由相位步进引起的抖动组成。如以上关于ML简单提到的，已经发现可以假设噪声是白噪声并且尤其是高斯白噪声。从概率性角度看，正弦图能够被认为是随机事件的样本，每个事件是在相应探测器单元处的探测(在所测量的正弦图中的条目)。然后能够将每个探测器单元与针对均值和方差的高斯参数(μ，σ²)相关联。

换言之并且在一般情况下，矩阵W^-1(W＝协方差矩阵Σ的平方根)表示在每个探测器像素中的噪声和在探测器像素之间的噪声相关性，并且成本函数C被配置为返回由噪声矩阵W^-1的条目(方差并且在更一般的情况中测量的协方差)加权的针对大偏差||Ax-b||的高成本。这允许更好地对在探测器水平处的期望噪声行为进行建模，并且在投影域中而不是在图像域中并入所述知识。如上所述，W^-1的对角元素基于针对高斯方差σ^s的估计。的条目被拆分成康普顿分量和微分相位分量，其中，Σ_C是包括康普顿方差的对角矩阵，并且Σ_p是包括微分相位方差的对角矩阵。

如比如在由DNGauntt等人在“X-raytubepotential，filtrationanddetectorconsiderationsindualenergychestradiography”(MedPhys21(2)，1994年2月，205页)一文在方程6a-c中描述的，能够通过高斯误差传播途径根据所测量的康普顿正弦图来估计针对康普顿分量的方差估计。

针对微分相位正弦图数据分量，如由Weber等人在“NoiseinX-raygrating-basedphase-contrastimaging”(MedicalPhysics，第38卷，第7号，2011年7月)一文在方程(26)、(29–在行2、列2处的矩阵条目)中定义的，能够通过2/[(均值强度(针对像素x，y的正弦图条目)*n*V²]根据相位对比正弦图来估计方差。V是可见性，并且n是Talbot顺序，其中n＝1、3或5或7。在一些实施例中，也设想n＝9。

在一个实施例中，协方差矩阵忽略在康普顿通道与微分相位通道之间的任何可能的协方差。然后矩阵W^-1是对角矩阵，其中每个对角条目是针对相应探测器像素/信号通道的相应估计的方差σ_i。这对充分近似是有效，但是所提出的方法能够同样利用具有至少一个非对角非零元素的协方差矩阵来实施，或更具体地，也能够利用全协方差矩阵来实施。这可以允许对如下情况进行建模：其中在噪声行为中存在已知的不对称性或相关性，例如在探测器D的特定像素区域中的主要串扰或余辉效应等。用于假设对角形式的协方差矩阵的一个原理在于依据步骤S310a、b的相位复原处理被应用到以逐像素样式的原始探测器原始数据，因此能够假设在不同像素中的噪声在充分近似中是不相关的。因此，协方差矩阵∑通常将甚至在更复杂的情况下具有相对简单的结构

$\mathrm{\Σ}=\left(\begin{array}{cc}{\mathrm{\Σ}}_C& {\mathrm{\Σ}}_{CP}\\ {\mathrm{\Σ}}_{CP}& {\mathrm{\Σ}}_P\end{array}\right)$

其中所有三个子矩阵∑是对角的。

由A、b和W的这些定义限定的成本函数实质上是未知x的函数，并且能够通过如非线性共轭梯度方法的标准技术来最小化。在该实施例中，针对x的更新函数通过标准共轭梯度方法来提供，其中，在“共轭方向”上更新x。在一个实施例中，共轭方向由Fletcher-Reeves公式定义，但是也设想针对步进尺寸的其他变型。

在备选实施例中，在每个迭代中用于更新电子密度图像估计x的更新函数基于如以上提到的T等人的论文中的SPS(可分离抛物面替代)方案。比如通过诸如经典牛顿-拉弗逊方法的牛顿类型算法来进行使似然函数(或对数似然函数)C最小化的搜索。在SPS中，“抛物面替代”是从上面与负对数似然相切的二次函数。凸面考虑得到位于成本函数以上的可分离的二次函数。可分离的替代函数的最小化比对数似然函数的最小化更易处理，并且确保单调收敛于对数似然的最小元。然而，这仅仅是一个实施例，并且本文也设想其他数值方法和/或成本函数结构。具体地，以上在成本函数的最小化方面已经表达了迭代重建，但是再次地，这是范例性实施例。有时，一些迭代重建方案可能要求不同的优化，比如为了使目标函数最大化，并且本文同样设想这样的方案，尤其出于在代数重排之后在最小化问题方面能够正常地表达任何最大化问题的原因。

在任一实施例中，在固定数量(在一个实施例中，其是用户可调节的)的迭代之后终止，或如果成本函数在迭代期间不再改变，例如，如果归因于更新的成本减少小于用户可定义的阈值ε，则终止迭代。

简言之，所提出的方法提供：

所有数据的单个噪声水平被正确地考虑并且被指定在探测器像素水平处，即在投影域中。

能够忽略在相移正弦图中的缺失的数据(出于相位缠绕、相干损失等的原因)，因为迭代重建被认为是关于此相当鲁棒的，并且因为其他通道，即康普顿散射正弦图仍然包括足够的投影信息来弥补在微分相位正弦图中缺失的数据。

即，在区别于诸如在WO2013/027138中描述的一个的其他途径中，由于代替地采取在所测量的数据上的噪声估计(并且在该域中噪声被认为是白噪声)，所以不要求在图像域中的SNR(信噪比)的参数化。

以上描述的方法和相关图像处理器可以被应用到所有微分相位对比CT成像系统，优选地被应用到具有诸如能量选择或光子计数探测器或其他的谱分辨探测器的成像系统。

在任选的步骤(在图3中未示出)处，然后在监控器M上绘制和显示电子密度横截面图像或根据多个电子密度横截面形成的体积。

为了影响可视化的各种模式，如本文提出的图像处理器IP包括可视化模块VS，其驱动系统视频卡以实现在监控器M上对所计算的电子密度的显示。

图像处理器IP和其部件可以被布置为在分布式架构中布置的并且在适当通信网络中连接的单独的模块。

在一个实施例中，图像处理器IP被布置为专用FPGA或作为硬连接独立芯片。

在备选实施例中，图像处理器IP和其部件驻留在工作站CON中，作为软件例程在工作站CON上运行。图像处理器IP和其部件可以在适当的科学计算平台(诸如)中被编程并且可以被转化为被维持在库中的并且当由工作站CON的操作系统调用时被链接的C++或C例程。

尽管在上文中已经描述了基于光栅的设置以从所测量的投影数据提取相移，但是同样地可以使用其他PCI(相位对比成像)原理，例如基于分析器的成像或基于传播的成像。

在本发明的另一范例性实施例中，提供一种计算机程序或计算机程序单元，其特征在于适于在适当系统上执行根据前述实施例中的任一项所述的方法的方法步骤。

因此，计算机程序单元可以被存储在计算机单元中，其也可以是本发明的实施例的一部分。该计算单元可以适于执行上述方法的步骤或诱发对上述方法的步骤的执行。而且，其可以适于操作上述装置的部件。计算单元能够适于自动操作和/或执行用户的顺序。计算机程序可以被加载到数据处理器的工作存储器中。数据处理器可以因此被装备为执行本发明的方法。

本发明的这个范例性实施例涵盖从最开始使用本发明的计算机程序和借助于更新将现有程序转变为使用本发明的程序的计算机程序。

更进一步地，所述计算机程序单元能够提供实现上述方法的范例性实施例的流程的所有必需步骤。

根据本发明的另一范例性实施例，提出了一种计算机可读介质，例如CD-ROM，其中，所述计算机可读介质具有存储在所述计算机可读介质上的计算机程序单元，其中，所述计算机程序单元由前面部分描述。

计算机程序可以存储和/或分布在与其他硬件一起提供或作为其他硬件的部分提供的诸如光学存储介质或固态介质的适当的介质上，但是计算机程序也可以以其他的形式分布，例如经由因特网或其他有线或无线的远程通信系统。

然而，所述计算机程序也可以存在于诸如万维网的网络上并能够从这样的网络中下载到数据处理器的工作存储器中。根据本发明的另一范例性实施例，提供了一种用于使得计算机程序单元能够被下载的介质，其中，所述计算机程序单元被布置为执行根据本发明的之前描述的实施例之一所述的方法。

必须指出，本发明的实施例参考不同主题加以描述。具体而言，一些实施例参考方法类型的权利要求加以描述，而其他实施例参考设备类型的权利要求加以描述。然而，本领域技术人员将从以上和下面的描述中了解到，除非另行指出，除了属于一种类型的主题的特征的任何组合之外，涉及不同主题的特征之间的任何组合也被认为被本申请公开。然而，所有特征能够被组合以提供超过特征的简单相加的协同效应。

尽管已经在附图和前面的描述中详细说明和描述了本发明，但是这样的说明和描述被认为是说明性或范例性的而非限制性的。本发明不限于所公开的实施例。通过研究附图、说明书和从属权利要求，本领域的技术人员在实践所主张的本发明时能够理解和实现所公开的实施例的其他变型。

在权利要求中，词语“包括”不排除其他元件或步骤，并且，词语“一”或“一个”并不排除多个。单个处理器或其他单元可以实现权利要求书中记载的若干项目的功能。在互不相同的从属权利要求中记载了特定措施并不指示不能有利地使用这些措施的组合。权利要求中的任何附图标记不得被解释为对范围的限制。

Claims (14)

1.一种用于生成对象(PAT)中或所述对象的部分中的电子密度的图像的方法，包括以下步骤：

接收(S330)：

i)针对至少一个能量值根据由扫描器(IMA)的辐射敏感探测器探测的数据导出的相位对比正弦图其中，涉及在辐射通过所述对象之后由所述辐射经历的相移的相移梯度；以及ii)根据由所述辐射敏感探测器探测的数据或者根据由不同辐射敏感探测器探测的数据导出的康普顿衰减正弦图p_C，其中，p_C实质上记录到对由在通过所述辐射在所述的通过所述对象中经历的总体衰减的康普顿散射贡献；

基于两个正弦图和p_C两者联合重建(S335)所述电子密度的横截面图像。

2.根据权利要求1所述的方法，其中，所述联合重建迭代地进行。

3.根据权利要求2所述的方法，其中，所述联合重建中的迭代由成本函数来驱动，所述成本函数包括对如在所述正弦图中记录的噪声进行建模的噪声项。

4.根据权利要求3所述的方法，其中，所述噪声项是协方差矩阵。

5.根据权利要求4所述的方法，其中，所述协方差矩阵具有对角形式。

6.根据权利要求4所述的方法，其中，所述协方差矩阵具有至少一个非零非对角元素。

7.根据权利要求3-6中的任一项所述的方法，其中，所述成本函数根据包括所述两个正弦图的对数似然函数来导出。

8.根据前述权利要求中的任一项所述的方法，其中，所述重建步骤基于共轭梯度方法。

9.根据权利要求1-7中的任一项所述的方法，其中，所述重建步骤基于牛顿-拉弗逊方法。

10.根据前述权利要求中的任一项所述的方法，其中，用于采集p_C正弦图和/或正弦图的所述扫描器(IMA)是多能量扫描器，尤其是双能量扫描器，或者其中，所述扫描器具有多层探测器，尤其是双层探测器，或者其中，所述扫描器是具有光子计数探测器的谱成像扫描器。

11.根据权利要求1-10中的任一项所述的方法，其中，所述扫描器(IMA)包括用于采集所述相位对比正弦图的微分相位对比成像装备(PCI)。

12.一种图像处理系统，包括：

一个或多个成像装置(IMA)，其被配置为供应对象的相位对比正弦图和康普顿衰减正弦图p_C；

处理单元(IP)，其被配置为通过根据权利要求1-11中的任一项所述的方法处理相位对比图像和所述康普顿衰减正弦图p_C来重建所述对象中的电子密度的图像；

显示单元(M)，其用于显示如由所述处理单元(IP)重建的所述电子密度的图像。

13.一种用于控制根据权利要求12所述的装置的计算机程序单元，所述计算机程序单元当由处理单元运行时适于执行根据权利要求1-11中的任一项所述的方法的步骤。

14.一种计算机可读介质，在其上存储有根据权利要求13所述的程序单元。