技术领域
[0001] 本
发明涉及X射线成像系统与成像方法领域，尤其涉及一种X射线光栅相衬成像滤波迭代重建方法和系统。
背景技术
[0002] 传统X射线成像技术利用X射线穿透物体时，物质对X射线的吸收能
力的差异进行成像。当不同物质之间的
密度差异较大时，该方法可以获得很好的成像衬度，例如人体成像时，骨骼可以很清晰地从周边软组织中分辨出来。但是当两种物质密度接近时，该方法往往成像效果不佳。
相位衬度成像(PCI)的
对比度是基于复折射率
实部的变化，复折射率实部的变化比X射线
能量范围内吸收指数
虚部的变化更为敏感。在PCI中，X射线通过物体的偏
角记录了相位信息。因此，基于PCI的计算机
断层扫描对于弱
吸收材料尤其是软组织具有良好的成像
质量。其中，基于光栅的相位衬度成像(GPCI)是一种最有潜力的临床技术，它可以利用光栅干涉仪重建和估计样品的衰减系数、相位信息和散射信息。而光栅相衬CT(GPC-CT)仍有许多技术局限性，如步进过程、
数据采集时间、
辐射剂量和连续旋转
机架等。重建
算法的改进是临床应用的必要条件。GPC-CT的标准重建方法是一种解析重建算法，即修正
滤波反投影(mFBP)。为了提高图像质量和抑制伪影，解析算法需要完整的投影数据，但这会增加
辐射剂量。
[0003] 近年来，GPCI断层成像因其具有高
分辨率成像的潜力和广阔的临床应用前景而成为一个研究热点[B]。Faris等人首先从理论上提出了光束偏转光学
层析成像，得到了定量的高分辨率成像结果。这是首次断层扫描技术可以应用于相衬成像，而不是传统的吸收成像。Pfeiffer等人利用
X射线管源设计并实现了GPC-CT。近年来，为了减少辐射剂量和体积，Wang等人提出了利用改进的系统矩阵实现差分相位衬度成像(DPCI)迭代重建的微分代数重建技术(DART)算法。Xu等人提出了Kaiser-Bessel窗函数来构造系统矩阵，并采用了两种迭代重建模型来进行DPCI层析成像。
[0004] 现有的GPC-CT迭代重建算法将微分投影数据直接映射到重建图像。这些算法本质上是解决了投影与图像之间的数学矩阵函数。然而，这些GPC-CT的迭代方法需要找出具体的系统矩阵，这与传统的迭代重建算法不同。
[0005]
现有技术所提供的光栅相衬重建方法，为了有快速的重建速度，大多采用了滤波反投影的解析重建技术，这些非迭代算法，所需投影角度数多，样品接收的X射线辐射剂量大，限制了该项技术的发展。
[0006] 现有技术所提供的光栅相衬迭代重建方法，都是由光栅相衬微分投影数据直接重建出光栅相衬重建图像，采用了一种端到端的重建方法，需要给出复杂的重建数学模型。
[0007] 因此，本领域的技术人员致力于开发一种X射线光栅相衬滤波迭代重建算法以及其重建系统，利用常规投影数据与微分相位数据的关系，通过严格的数学推导和物理模型，提出一种新的迭代重建算法，解决了传统重建算法的不可重复性和不可扩展性问题。
发明内容
[0008] 有鉴于现有技术的上述
缺陷，本发明所要解决的技术问题是现有技术所提供的光栅相衬迭代重建方法，是由光栅相衬微分投影数据直接重建出光栅相衬重建图像，采用了一种端到端的重建方法，需要给出复杂的重建数学模型。
[0009] 为实现上述目的，本发明提供了一种X射线光栅相衬成像滤波迭代重建方法，包括以下步骤：
[0011] 所述X射线束穿透被检测物体后，通过光栅干涉仪的同步辐射装置，采集对应的差分相位衬度投影数据；
[0012] 通过相位恢复算法提取所述投影数据对应的折射率信息；
[0013] 对所述投影数据进行滤波处理，将所述差分相位衬度投影数据转换为正积分投影数据；
[0014] 采用迭代算法对经过滤波处理后的相位信息进行图像重建。
[0015] 在本发明的较佳实施方式中，所述投影数据是平行投影、扇形投影或锥形束投影；所述重建是二维重建或
三维重建。
[0016] 在本发明的较佳实施方式中，所述迭代算法包括ART、Siddon-ART和基于全变分正则化算法中的任一种。
[0017] 在本发明的较佳实施方式中，所述差分相位衬度投影数据表示为：
[0018]
[0019] 待重建的目标函数为：
[0020]
[0021] 其中，(s,t)表示旋转θ度的
坐标系，s是投影方向，t是垂直于s的方向；δ表示衍射指数衰减，Δ(t,θ)对应于坐标系(s,t)中线积分的投影函数，Γ(ω,θ)是Δ(t,θ)的
傅立叶变换，F-1表示反傅立叶变换，ω为法线积分的滤波函数；
[0022] 所述同步辐射装置采集到的投影数据用下面的公式表示：
[0023]
[0024] 则，待重建的目标函数表示为：
[0025]
[0026] 其中，P(ω,θ)是p(t,θ)的傅立叶变换，
滤波器函数H(ω)＝(ω)/2πiω，根据式(3)，P(t,θ)表示为：
[0027]
[0028] 通过分析式(2)和式(4)，建立两个不同线积分的傅立叶变换之间的联系，用下面的公式表示：
[0029]
[0030] 其中，滤波器L(ω)＝1/2πiω，用于将所述差分相位衬度投影数据转换为正积分投影数据。
[0031] 在本发明的较佳实施方式中，提供了一种X射线光栅相衬成像滤波迭代重建系统，包括：分布式X射线源，用于向被检测物体发射X射线束；固定光栅模
块，包括具有同步辐射装置的光栅干涉仪，用于采集所述X射线束穿透被检测物体后的差分相位衬度投影数据，所述被检测物体位于所述分布式X射线源与固定光栅模块之间；X射线探测器，位于所述被检测物体和所述光栅模块的后面，用于将所述差分相位衬度投影数据转换为可识别的电
信号；
数据处理模块，用于处理所述
电信号的数据，并从中对所述X射线束的差分相位衬度投影数据进行滤波迭代重建处理，由此重建出被检测物体的图像信息。
[0032] 在本发明的较佳实施方式中，所述投影数据是平行投影、扇形投影或锥形束投影；所述重建是二维重建或三维重建。
[0033] 在本发明的较佳实施方式中，所述迭代重建处理包括ART、Siddon-ART和基于全变分正则化重建处理中的任一种。
[0034] 在本发明的较佳实施方式中，所述数据处理模块包括滤波单元和重建单元，其中，滤波单元，用于将所述差分相位衬度投影数据转换为正积分投影数据，包括以下公式：
[0035] 所述差分相位衬度投影数据表示为：
[0036]
[0037] 待重建的目标函数为：
[0038]
[0039] 其中，(s,t)表示旋转θ度的坐标系，s是投影方向，t是垂直于s的方向；δ表示衍射指数衰减，Δ(t,θ)对应于坐标系(s,t)中线积分的投影函数，Γ(ω,θ)是Δ(t,θ)的傅立叶变换，F-1表示反傅立叶变换，ω为法线积分的滤波函数；
[0040] 所述同步辐射装置采集到的投影数据用下面的公式表示：
[0041]
[0042] 则，待重建的目标函数表示为：
[0043]
[0044] 其中，P(ω,θ)是p(t,θ)的傅立叶变换，滤波器函数H(ω)＝(ω)/2πiω，根据式(3)，P(t,θ)表示为：
[0045]
[0046] 通过分析式(2)和式(4)，建立两个不同线积分的傅立叶变换之间的联系，用下面的公式表示：
[0047]
[0048] 其中，滤波器L(ω)＝1/2πiω；
[0049] 重建单元采用迭代算法，将上述投影数据进行重构，得到重建后的被检测物体图像信息。
[0050] 在本发明的另一较佳实施方式中，提供了一种X射线光栅相衬成像滤波迭代重建装置，包括：
存储器，用于存储程序；处理器，用于通过执行所述存储器存储的程序以实现如上所述的方法。
[0051] 在本发明的另一较佳实施方式中，提供了一种计算机可读存储介质，包括程序，所述程序能够被处理器执行以实现如上所述的方法。
[0052] 与现有技术相比，本发明提供的技术方案具有以下优点：本发明提出的方法不直接计算系统矩阵，而是利用常规投影数据与微分相位数据的关系，通过严格的数学推导和物理模型加以说明。因此，将迭代重建过程分为两个简单的部分，包括滤波操作和常规迭代重建。本发明提出的方法有效、准确地简化了现有技术中的迭代重建过程。
[0053] 以下将结合
附图对本发明的构思、具体结构及产生的技术效果作进一步说明，以充分地了解本发明的目的、特征和效果。
附图说明
[0054] 图1是本发明的一个较佳
实施例中的X射线光栅相衬成像滤波迭代重建方法
流程图；
[0055] 图2是本发明的另一个较佳实施例中采集得到的被测物体差分相位衬度投影数据坐标图；
[0056] 图3是本发明的另一个较佳实施例中的普通滤波反投影中的Ram-Lak滤波器示意图；
[0057] 图4是本发明的另一个较佳实施例中的针对微分相衬投影所用的Hilbert滤波器示意图；
[0058] 图5是本发明的另一个较佳实施例中的针对微分相衬迭代重建设计的滤波器示意图；
[0059] 图6是本发明的另一个较佳实施例中的标准Shepp-Logan头部数字模型图像、本
专利算法重建图像以及两图像在同一
位置的灰度线对比图；
[0060] 图7是本发明的另一个较佳实施例中的瑞士同步辐射
光源采集的正常大鼠肾组织的在不同投影角度下的重建图像对比图；
[0061] 图8是本发明的另一个较佳实施例中的X射线光栅相衬成像滤波迭代重建系统示意图。
具体实施方式
[0062] 以下参考
说明书附图介绍本发明的多个优选实施例，使其技术内容更加清楚和便于理解。本发明可以通过许多不同形式的实施例来得以体现，本发明的保护范围并非仅限于文中提到的实施例。
[0063] 在附图中，结构相同的部件以相同数字标号表示，各处结构或功能相似的组件以相似数字标号表示。附图所示的每一组件的尺寸和厚度是任意示出的，本发明并没有限定每个组件的尺寸和厚度。为了使图示更清晰，附图中有些地方适当夸大了部件的厚度。
[0064] 本实施例中的一些术语：
[0065] GPC-CT：光栅相衬CT；
[0066] GPCI：基于光栅的相位衬度成像；
[0067] DPCI：差分相位衬度成像；
[0068] DART：微分代数重建技术；
[0069] mFBP：修正滤波反投影。
[0070] 如图1所示，一种X射线光栅相衬成像滤波迭代重建方法，包括以下步骤：
[0071] S01：向被检测物体发射X射线束。
[0072] S02：所述X射线束穿透被检测物体后，通过光栅干涉仪的同步辐射装置，采集对应的差分相位衬度投影数据；其中，所述投影数据是平行投影、扇形投影或锥形束投影。
[0073] S03：通过相位恢复算法提取所述投影数据对应的折射率信息。其中，所述的相位恢复算法为现有标准的提取光栅相衬成像相位信号的算法技术，具体请参考文献J.R.Fienup,"Phase retrieval algorithms:a personal tour[Invited]."Applied Optics,vol.52pp.45-56,Dec.2013.
[0074] S04：对所述投影数据进行滤波处理，将所述差分相位衬度投影数据转换为正积分投影数据。
[0075] S05：采用迭代算法对经过滤波处理后的相位信息进行图像重建。其中，所述重建是二维重建或三维重建。所述迭代算法包括ART、Siddon-ART和基于全变分正则化算法中的任一种。
[0076] 在本实施例中的S03和S04步骤中，参见图2-图5，所述差分相位衬度投影数据转换为正积分投影数据，包括以下公式：
[0077] 所述差分相位衬度投影数据表示为：
[0078]
[0079] 待重建的目标函数为：
[0080]
[0081] 其中，(s,t)表示旋转θ度的坐标系，s是投影方向，t是垂直于s的方向；δ表示衍射指数衰减，Δ(t,θ)对应于坐标系(s,t)中线积分的投影函数，Γ(ω,θ)是Δ(t,θ)的傅立叶1
变换，F-表示反傅立叶变换，ω为法线积分的滤波函数；
[0082] 所述同步辐射装置采集到的投影数据用下面的公式表示：
[0083]
[0084] 则，待重建的目标函数表示为：
[0085]
[0086] 其中，P(ω,θ)是p(t,θ)的傅立叶变换，滤波器函数H(ω)＝(ω)/2πiω，根据式(3)，P(t,θ)表示为：
[0087]
[0088] 通过分析式(2)和式(4)，建立两个不同线积分的傅立叶变换之间的联系，用下面的公式表示：
[0089]
[0090] 其中，滤波器L(ω)＝1/2πiω；用于将所述差分相位衬度投影数据转换为正积分投影数据。如图3-图5所示，不同的投影数据需要采用不同的滤波器。
[0091] S05步骤采用迭代算法，将经过滤波处理后的相位信息进行重构，得到重建后的被检测物体图像信息。
[0092] 在本实施中，采用标准的Shepp-Logan头部模型作为数值模拟模型，对本发明提出的一种X射线光栅相衬成像滤波迭代重建方法进行数值模拟实验。本实验中采用了ART迭代重建算法。其中，重建切片的
体模尺寸为512×512。差分投影采用等角度间距的90个投影角度。
像素的灰度表示折射率系数的近似值。得到了模拟GPCI的雷登变换梯度。在这种情况下，省略了光栅相位步进过程，因此模拟重建不需要相位恢复步骤。如图6中(C)图所示，通过构建灰度曲线来比较模型图像和重建图像之间的差异，重建前后的图像灰度和对比度非常相近。
[0093] 在本实施中，还采用了瑞士X射线同步辐射光源光束线放射学实验(TOMCAT)的真实实验数据进行验证：以正常大鼠肾组织为样本，采用五步相位步进法记录折射率信息，并采用相位重聚算法进行相位恢复。通过该验证方法，不仅研究了稀疏角度重建的成像质量，而且展示并对比了不同数据采集角度下的重建结果。如图7所示，改进的FBP重建与本实施例中的重建方法比较表明，本实施例中的重建方法能够满足GPCI层析成像任务，重建结果非常接近。本实施例中的重建方法可以有效地重建GPC-CT数据。通过不同角度，本实施例中的重建方法的结果研究了投影角度数量对重建图像质量影响。其中，75个及以上个投影角度的图像质量没有明显下降，50个投影角度和60个投影角度的重建图像具有相对更明显的噪声像素。
[0094] 与现有技术相比，本实施例提供的重建方法具有以下优点：本实施例中提出的重建方法不直接计算系统矩阵，而是利用常规投影数据与微分相位数据的关系，通过严格的数学推导和物理模型加以说明。因此，将迭代重建过程分为两个简单的部分，包括滤波操作和常规迭代重建。本实施例中提出的重建方法有效、准确地简化了现有技术中的迭代重建过程。
[0095] 在另一个实施例中，如图8所示，还提出了一种X射线光栅相衬成像滤波迭代重建系统，包括：
[0096] 分布式X射线源A01，用于向被检测物体A00发射X射线束。
[0097] 固定光栅模块A02，包括具有同步辐射装置的光栅干涉仪，用于采集所述X射线束穿透被检测物体A00后的差分相位衬度投影数据，所述被检测物体A00位于所述分布式X射线源A01与固定光栅模块A02之间；所述投影数据是平行投影、扇形投影或锥形束投影。
[0098] X射线探测器A03，位于所述被检测物体A00和所述光栅模块A02的后面，用于将所述差分相位衬度投影数据转换为可识别的电信号；
[0099] 数据处理模块A04，用于处理所述电信号的数据，并从中对所述X射线束的差分相位衬度投影数据进行滤波迭代重建处理，由此重建出被检测物体A00的图像信息。所述重建是二维重建或三维重建。迭代重建处理包括ART、Siddon-ART和基于全变分正则化重建处理中的任一种。
[0100] 所述数据处理模块A04包括滤波单元和重建单元，其中，滤波单元，用于将所述差分相位衬度投影数据转换为正积分投影数据，包括以下公式：
[0101] 所述差分相位衬度投影数据表示为：
[0102]
[0103] 待重建的目标函数为：
[0104]
[0105] 其中，(s,t)表示旋转θ度的坐标系，s是投影方向，t是垂直于s的方向；δ表示衍射指数衰减，Δ(t,θ)对应于坐标系(s,t)中线积分的投影函数，Γ(ω,θ)是Δ(t,θ)的傅立叶变换，F-1表示反傅立叶变换，ω为法线积分的滤波函数；
[0106] 所述同步辐射装置采集到的投影数据用下面的公式表示：
[0107]
[0108] 则，待重建的目标函数表示为：
[0109]
[0110] 其中，P(ω,θ)是p(t,θ)的傅立叶变换，滤波器函数H(ω)＝(ω)/2πiω，根据式(3)，P(t,θ)表示为：
[0111]
[0112] 通过分析式(2)和式(4)，建立两个不同线积分的傅立叶变换之间的联系，用下面的公式表示：
[0113]
[0114] 其中，滤波器L(ω)＝1/2πiω。如图3-图5所示，在不同的滤波单元中，不同的投影数据需要不同的滤波器。
[0115] 重建单元采用迭代算法，将上述投影数据进行重构，得到重建后的被检测物体图像信息。
[0116] 与现有技术相比，本实施例提供的重建系统具有以下优点：本实施例中提出的重建系统不用计算系统矩阵，而是利用常规投影数据与微分相位数据的关系，通过严格的数学推导和物理模型加以说明。因此，将现有技术中的迭代重建模块分为两个简单的部分，包括滤波单元和常规迭代重建单元。本实施例中提出的重建系统有效、准确地简化了现有技术中的迭代重建系统。
[0117] 本领域技术人员可以理解，上述实施方式中各种方法的全部或部分功能可以通过
硬件的方式实现，也可以通过
计算机程序的方式实现。当上述实施方式中全部或部分功能通过计算机程序的方式实现时，该程序可以存储于一计算机可读存储介质中，存储介质可以包括：
只读存储器、随机存储器、磁盘、光盘、
硬盘等，通过计算机执行该程序以实现上述功能。例如，将程序存储在设备的存储器中，当通过处理器执行存储器中程序，即可实现上述全部或部分功能。另外，当上述实施方式中全部或部分功能通过计算机程序的方式实现时，该程序也可以存储在
服务器、另一计算机、磁盘、光盘、闪存盘或移动硬盘等存储介质中，通过下载或复制保存到本地设备的存储器中，或对本地设备的系统进行版本更新，当通过处理器执行存储器中的程序时，即可实现上述实施方式中全部或部分功能。
[0118] 以上详细描述了本发明的较佳具体实施例。应当理解，本领域的普通技术无需创造性劳动就可以根据本发明的构思作出诸多
修改和变化。因此，凡本技术领域中技术人员依本发明的构思在现有技术的
基础上通过逻辑分析、推理或者有限的实验可以得到的技术方案，皆应在由
权利要求书所确定的保护范围内。
