CN110806598A - 基于剂量分布的信息提取方法及x射线相衬成像系统 - Google Patents

基于剂量分布的信息提取方法及x射线相衬成像系统 Download PDF

Info

Publication number
CN110806598A
CN110806598A CN201911098795.3A CN201911098795A CN110806598A CN 110806598 A CN110806598 A CN 110806598A CN 201911098795 A CN201911098795 A CN 201911098795A CN 110806598 A CN110806598 A CN 110806598A
Authority
CN
China
Prior art keywords
mma
information extraction
ray
dose
information
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201911098795.3A
Other languages
English (en)
Other versions
CN110806598B (zh
Inventor
张丽
吴承鹏
高河伟
邢宇翔
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Tsinghua University
Original Assignee
Tsinghua University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Tsinghua University filed Critical Tsinghua University
Priority to CN201911098795.3A priority Critical patent/CN110806598B/zh
Publication of CN110806598A publication Critical patent/CN110806598A/zh
Application granted granted Critical
Publication of CN110806598B publication Critical patent/CN110806598B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01T—MEASUREMENT OF NUCLEAR OR X-RADIATION
    • G01T1/00—Measuring X-radiation, gamma radiation, corpuscular radiation, or cosmic radiation
    • G01T1/29—Measurement performed on radiation beams, e.g. position or section of the beam; Measurement of spatial distribution of radiation
    • G01T1/2914—Measurement of spatial distribution of radiation
    • G01T1/2964—Scanners
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B6/00—Apparatus for radiation diagnosis, e.g. combined with radiation therapy equipment
    • A61B6/50—Clinical applications
    • A61B6/502—Clinical applications involving diagnosis of breast, i.e. mammography
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01T—MEASUREMENT OF NUCLEAR OR X-RADIATION
    • G01T1/00—Measuring X-radiation, gamma radiation, corpuscular radiation, or cosmic radiation
    • G01T1/29—Measurement performed on radiation beams, e.g. position or section of the beam; Measurement of spatial distribution of radiation
    • G01T1/2914—Measurement of spatial distribution of radiation
    • G01T1/2992—Radioisotope data or image processing not related to a particular imaging system; Off-line processing of pictures, e.g. rescanners

Abstract

一种基于剂量分布的信息提取方法及X射线相衬成像系统，所述信息提取方法包括：以机械步进方式或等同于机械步进的方式在设置物体和不设置物体的情况下分别进行曝光；X射线探测器在曝光后进行探测，获取背景位移曲线和物体位移曲线；采用一信息提取算法分析背景位移曲线和物体位移曲线，以提取出物体的相衬或暗场信息；其中，机械步进方式或等同于机械步进的方式中，各步的X射线剂量不是均匀的，各步的X射线剂量遵循一特定剂量分布，所述特定剂量分布使得总剂量保持固定的情况下相衬或者暗场信息的方差最小且均值不变。该分布降低了数据的噪声，可以在相同总剂量下得到更低噪声水平的对比度信息或者在相同的噪声水平下降低采集需要的总剂量。

Description

基于剂量分布的信息提取方法及X射线相衬成像系统
技术领域
本公开属于X射线成像技术领域，涉及一种基于剂量分布的信息提取方法及X射线相衬成像系统。
背景技术
X射线相衬成像技术，能够实现微米甚至亚微米级上的局部结构分辨，是对传统X射线衰减成像技术的很好补充。该技术能够同时提取吸收、相衬和暗场三种对比度信息，适用于低原子序数、低密度物质，特别是针对包括乳腺在内的生物软组织结构。
根据成像原理的不同，目前在常规X光源上实现的X射线相衬成像主要包括基于分析晶体的成像(ABI，Analyzer-based Imaging)、基于光栅的成像(GI，Grating-basedImaging)、基于斑点的成像(SI，Speckle-based Imaging)以及边缘照射成像(EI，Edge-Illumination)。为了获得由于射线在物质中发生微小折射和散射等包含的信息，除了基于斑点的成像，以上另外三种X射线相衬成像系统通常都通过步进扫描方式在每个探测器像素上获得一条位移曲线，如在X射线光栅相衬成像技术中通过相位步进的方式来获得相位步进曲线，该曲线近似为余弦曲线，而在边缘照射成像中获得的位移曲线近似为高斯曲线。分别在有物体和没有物体的情况下获得背景位移曲线和物体位移曲线，然后通过信息提取方法从两条曲线的差异中获得物体对应的衰减信息、折射信息和小角散射信息，即分别对应于吸收、相衬和暗场三种图像对比度。
已有的分析方法大部分是通过优化信息提取算法来降低数据的噪声。仍有需要来进一步改善噪声问题。
发明内容
(一)要解决的技术问题
本公开提供了一种基于剂量分布的信息提取方法及X射线相衬成像系统，首次提出通过变化步进扫描方式中每一步的X射线剂量，使得X射线的剂量呈现一定规律的分布，该分布降低了数据的噪声，可以在相同总剂量的情况下得到更低噪声水平的对比度信息或者在获得相同噪声水平下的对比度信息需要更少的采集总剂量。
(二)技术方案
根据本公开的一个方面，提供了一种基于剂量分布的信息提取方法，用于X射线相衬成像系统，所述信息提取方法包括：以机械步进方式或等同于机械步进的方式在设置物体和不设置物体的情况下分别进行曝光；X射线探测器在曝光后进行探测，获取不设置物体的背景位移曲线和设置物体的物体位移曲线；采用一信息提取算法分析背景位移曲线和物体位移曲线，以提取出物体的相衬或暗场信息；其中，所述机械步进方式或所述等同于机械步进的方式中，各步的X射线剂量不是均匀的，各步的X射线剂量遵循一特定剂量分布，所述特定剂量分布使得总剂量保持固定的情况下相衬或者暗场信息的方差最小且均值不变。
所述信息提取算法包括如下算法的一种或几种：D-MMA、GD-MMA、TA-MMA以及DB-MMA。
在本公开的一实施例中，当所述信息提取算法为D-MMA、GD-MMA、TA-MMA或DB-MMA时，所述特定剂量分布与信息提取算法中的核函数相关。
在本公开的一实施例中，当所述信息提取算法为D-MMA、GD-MMA或TA-MMA时，所述特定剂量分布中各步的X射线剂量与信息提取算法中的核函数成正比，特定剂量分布以核函数作为权重。
当所述信息提取算法为DB-MMA时，所述特定剂量分布满足：
T为常数
Figure BDA0002267856810000022
x＝eiRF-1＝[x0，x1，...，xN-1]
Figure BDA0002267856810000031
其中，t(n)表示各步的剂量；n＝1，...，N，N为总步数；ei表示在第i个位置为1其它均为0的三维单位向量；R表示三种对比度信息的矩阵表达式中的雅可比矩阵，上述雅可比矩阵中，各个元素分别为对比度信息对应的核函数；F是一个循环矩阵，该循环矩阵中各个元素为背景位移曲线。
在本公开的一实施例中，所述各步的X射线剂量通过调整每一步的X射线源的电流、曝光时间、X射线源的电压或者X射线源的光强中的一个或多个来控制各步的X射线剂量，从而使得各步的X射线剂量遵循特定剂量分布。
在本公开的一实施例中，提供所述X射线的X射线源为如下光源中的一种：X光源、微焦点X光源、以及同步辐射X光源；
所述X射线探测器包括如下探测器中的一种：能量积分型探测器或者光子计数型探测器。
在本公开的一实施例中，所述X射线相衬成像系统包括如下成像系统中的一种：基于分析晶体的成像系统、基于光栅的成像系统以及边缘照射成像系统；
其中，所述基于光栅的成像系统包括Talbot-Lau类型、几何投影类型和双相位光栅类型中的一种。
根据本公开的另一个方面，提供了一种X射线相衬成像系统，所述X射线相衬成像系统以机械步进方式或等同于机械步进的方式获得背景位移曲线和物体位移曲线；采用一信息提取算法分析背景位移曲线和物体位移曲线，以提取出物体的相衬或暗场信息；其中，所述机械步进方式或所述等同于机械步进的方式中，各步的X射线剂量不是均匀的，各步的X射线剂量遵循一特定剂量分布，所述特定剂量分布使得总剂量保持固定的情况下相衬或者暗场信息的方差最小且均值不变。
在本公开的一实施例中，所述信息提取算法包括如下算法的一种或几种：D-MMA、GD-MMA、TA-MMA以及DB-MMA。
在本公开的一实施例中，当所述信息提取算法为D-MMA、GD-MMA、TA-MMA或DB-MMA时，所述特定剂量分布与信息提取算法中的核函数相关。
在本公开的一实施例中，当所述信息提取算法为D-MMA、GD-MMA或TA-MMA时，所述特定剂量分布中各步的X射线剂量与信息提取算法中的核函数成正比，特定剂量分布以核函数作为权重。
当所述信息提取算法为DB-MMA时，所述特定剂量分布满足：
Figure BDA0002267856810000041
T为常数
Figure BDA0002267856810000042
x＝eiRF-1＝[x0，x1，...，xN-1]
其中，t(n)表示各步的剂量；n＝1，...，N，N为总步数；ei表示在第i个位置为1其它均为0的三维单位向量；R表示三种对比度信息的矩阵表达式中的雅可比矩阵，上述雅可比矩阵中，各个元素分别为对比度信息对应的核函数；F是一个循环矩阵，该循环矩阵中各个元素为背景位移曲线。
(三)有益效果
从上述技术方案可以看出，本公开提供的基于剂量分布的信息提取方法及X射线相衬成像系统，具有以下有益效果：
(1)首次提出通过变化步进扫描方式中每一步的X射线剂量，使得各步X射线的剂量呈现一定规律的分布，所述特定剂量分布使得总剂量保持固定的情况下相衬或者暗场信息的方差最小且均值不变，使得提取出的相衬或暗场信息数值不变但噪声最小，该分布降低了数据的噪声，可以在相同总剂量的情况下得到更低噪声水平的对比度信息；在保证相同噪声水平的对比度信息条件下，降低数据采集过程中的总体剂量水平；
(2)该变化剂量分布的求解方法适用于多种X射线相衬成像系统，例如为ABI、EI和GI，并且GI系统类别不限，即包括Talbot-Lau类型、几何投影类型和双相位光栅类型等各种类型的系统，适用范围广泛。
附图说明
图1为根据本公开一实施例所示的基于剂量分布的信息提取方法的流程图。
图2为根据本公开一实施例所示的(a)X射线相衬成像系统示意图以及(b)两种位移曲线的示意图。
图3为采用如图1所示的信息提取方法获取图像信息的过程示意图。
图4为采用本公开一实施例所示的信息提取方法中的最优曝光时间分布与传统方式对比的结果示意图。
图5为采用D-MMA算法处理模拟数据情况，(a)均匀分布的曝光时间和(b)优化后的加权分布曝光时间采集方式提取的相衬信息剖面线对比图。
图6为根据本公开一实施例所示的采用D-MMA算法处理模拟数据情况，均匀分布的曝光时间和优化后的加权分布曝光时间采集方式提取的相衬信息噪声水平对比图。
具体实施方式
现有的采集背景位移曲线和物体位移曲线这两种位移曲线的方法通常都是基于均匀间隔的机械结构步进，同时每一步的X射线剂量是相同的，即X射线源电压、电流、曝光时间等保持不变。
本公开基于MMA类数据分析方法的特点，提出一种基于剂量分布的信息提取方法及X射线相衬成像系统，首次提出通过变化步进扫描方式中每一步的X射线剂量，使得X射线的剂量呈现一定规律的分布，该分布降低了数据的噪声，可以在相同总剂量的情况下得到更低噪声水平的对比度信息。换言之，可以在保证相同噪声水平的对比度信息条件下，降低数据采集过程中的总体剂量水平，这对于X射线相衬成像的很多实际应用具有非常重要的意义。
本公开的X射线相衬成像系统的剂量分布方法及系统是针对于基于机械步进方式的成像系统，或者其它可以获得类似背景位移曲线和物体位移曲线的成像方式或系统，诸如公开号为“WO2016070771A1”、发明名称为“X射线相衬成像系统与成像方法”中的基于垂直于光栅条纹方向的分布式非相干X射线源依次曝光方式，或者中国专利申请号201910439158.1中的交错式光栅或倾斜式光栅和平行于光栅条纹方向的分布式非相干X射线源依次曝光方式，均能够形成类似机械步进方式的背景位移曲线和物体位移曲线的X射线相衬成像成像系统。本申请的剂量分布方法适用的系统类型包括ABI、EI和GI，并且GI系统类别不限，即包括Talbot-Lau类型、几何投影类型和双相位光栅类型等各种类型的系统。此外物体位置既可以是放置在光源与第一块光栅之间，也可以是放置在两块光栅之间。此外X射线源既可以是常规X光源，也可以是微焦点X光源或者同步辐射X光源。此外X射线探测器既可以是能量积分型探测器也可以是光子计数型探测器。
为使本公开的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本公开进一步详细说明。
第一实施例
在本公开的第一个示例性实施例中，提供了一种X射线相衬成像系统的剂量分布方法。
在其它条件不变的情况下，X射线的剂量正比于X射线源的电流以及曝光时间，也与X射线源的电压有一定的关系，当然也正比于X射线源的光强，从而根据比尔定律正比于探测器接收到的光强。本实施例中，为了比较明确地调节每一步的X射线剂量，通过改变X射线源的电流或者曝光时间来实现X射线的剂量调控。需要明确的是，改变每一步的剂量大小会使得原有的信息提取方法出现问题，进行剂量分布控制需要克服对于信息提取方法造成的不利影响同时还需要辅助该信息提取方法降低数据噪声。
图1为根据本公开一实施例所示的X射线相衬成像系统的剂量分布方法的流程图。
参照图1所示，本公开的基于剂量分布的信息提取方法，包括：
步骤S11：以机械步进方式或等同于机械步进的方式在设置物体和不设置物体的情况下分别进行曝光；其中，所述机械步进方式或所述等同于机械步进的方式中，各步的X射线剂量不是均匀的，各步的X射线剂量遵循一特定剂量分布，所述特定剂量分布使得总剂量保持固定的情况下相衬或者暗场信息的方差最小且均值不变；
在本公开的一实施例中，所述信息提取算法包括如下算法的一种或几种：D-MMA、GD-MMA、TA-MMA以及DB-MMA。
在本公开的一实施例中，当所述信息提取算法为D-MMA、GD-MMA、TA-MMA或DB-MMA时，所述特定剂量分布与信息提取算法中的核函数相关。
在本公开的一实施例中，当所述信息提取算法为D-MMA、GD-MMA或TA-MMA时，所述特定剂量分布中各步的X射线剂量与信息提取算法中的核函数成正比，特定剂量分布以核函数作为权重。
在本公开的一实施例中，当所述信息提取算法为DB-MMA时，所述特定剂量分布满足：
Figure BDA0002267856810000071
T为常数
Figure BDA0002267856810000072
x＝eiRF-1＝[x0，x1，...，xN-1]
Figure BDA0002267856810000073
其中，t(n)表示各步的剂量；n＝1，...，N，N为总步数；ei表示在第i个位置为1其它均为0的三维单位向量；R表示三种对比度信息的矩阵表达式中的雅可比矩阵，上述雅可比矩阵中，各个元素分别为对比度信息对应的核函数；F是一个循环矩阵，该循环矩阵中各个元素为背景位移曲线。
后续以具体实例来介绍上述特定剂量分布的推导过程。
步骤S12：X射线探测器在曝光后进行探测，获取不设置物体的背景位移曲线和设置物体的物体位移曲线；
在本公开的一实施例中，所述X射线成像系统包括如下成像系统中的一种：基于分析晶体的成像系统(ABI)、基于光栅的成像系统(GI)以及边缘照射成像系统(EI)；
其中，所述基于光栅的成像系统(GI)包括Talbot-Lau类型、几何投影类型和双相位光栅类型中的一种。
图2为根据本公开一实施例所示的(a)X射线相衬成像系统示意图以及(b)两种位移曲线的示意图。
本实施例中以Talbot-Lau类型的GI系统作为示例，参照图2中(a)所示，该X射线相衬成像系统主要包括X射线源S、源光栅G0、相位光栅G1、分析光栅G2以及X射线探测器。
提供X射线的X射线源为如下光源中的一种：X光源(常规X光源)、微焦点X光源、以及同步辐射X光源。
源光栅G0、相位光栅G1和分析光栅G2形成固定光栅模块P，源光栅G0和相位光栅G1之间的距离为L，相位光栅G1和分析光栅G2之间的距离为D。此外扫描物体W的位置既可以是位于光源与第一块光栅之间，也可以是放置在其中两块光栅之间。
所述X射线探测器包括如下探测器中的一种：能量积分型探测器或者光子计数型探测器。
参照图2中(b)所示，步进位置为横坐标，X射线探测器探测到的信号强度为纵坐标，不放置物体对应的背景位移曲线用圆圈示意，放置物体对应的物体位移曲线用正方形示意。下面用f(φ)和s(φ)分别表示背景位移曲线和物体位移曲线，其中φ表示步进过程中的位置，本实施例中，步骤S11调整的是步进过程中每一步的X射线剂量，剂量在图2中(b)没有体现。
步骤S13：采用一信息提取算法分析背景位移曲线和物体位移曲线，以提取出物体的相衬或暗场信息；
本实施例中，所述信息提取算法包括如下算法的一种或几种：D-MMA、GD-MMA、TA-MMA以及DB-MMA。其中，D-MMA方法在提出时并不适用于X射线光栅综合成像系统，申请人发现如果对背景位移曲线进行平移变换使得最大值位于中心位置，同时对相应的物体位移曲线平移相同距离，该方法就能够适用于GI系统。在上述基础上，申请人已有专利-中国专利申请号为201910593224.0发现了两种位移曲线均满足一种积分特性，即：位移曲线和任意一个函数的内积只包含其傅里叶级数中的零阶项和一阶项，于是建立了一种广义的D-MMA方法，也就是GD-MMA，拓宽了核函数的范围并且该专利申请提出了一种最优化的核函数形式。在申请人的专利申请号为201910207542.9的专利中提出了另一种解析MMA类信息提取算法，原文档记载为：ASAXS(analytic small angle X-ray scattering，解析式小角散射)方法，这里将其称之为TA-MMA，本申请中的基于剂量分布的信息提取方法在上述三个专利中的信息提取算法均是适用的。
图3为采用如图1所示的信息提取方法获取图像信息的过程示意图。
下面先具体介绍D-MMA、DB-MMA、GD-MMA和TA-MMA的提取信息的方法，然后参照图3以其中一种特定成像系统下(例如以Talbot-Lau类型的GI系统为例)、采用特定信息提取算法(例如以GD-MMA信息提取算法为例)提取出特定信息(例如以相衬信息为例，当然也可以是暗场信息)为例说明基于剂量分布的信息提取方法中确定各步的X射线剂量遵循一特定剂量分布，以及特定剂量分布的表达形式的推导过程。
在本实例中，在Talbot-Lau类型的GI系统下，采用GD-MMA信息提取算法提取出相衬信息，各步的X射线剂量不是均匀的，各步的X射线剂量遵循一特定剂量分布。
(一)MMA类信息提取算法介绍
1、MMA方法
基于小角散射(SAXS，Small Angle X-ray Scattering)的多阶矩分析(MMA，Multi-order Moments Analysis)为基于卷积的数据分析方法。MMA类方法的基本假设为：物体位移曲线s(φ)可以表示为背景位移曲线f(φ)和小角散射分布g(φ)的卷积，表达式如下：
2、DB-MMA法
在DB-MMA中，DB-MMA为基于解卷积的MMA方法(DB-MMA，Deconvolution-BasedMMA)，首先通过Lucy-Richardson迭代的方法，从采集到的每个X射线探测器像素上的背景位移曲线f(φ)和物体位移曲线s(φ)中解卷积获得小角散射分布，这里为了突出重点，括号里面的参量-步进过程中的位置φ进行了省略，全文省略方式相同，解卷积过程中的第k次迭代计算公式为：
Figure BDA0002267856810000101
其中，
Figure BDA0002267856810000102
表示f关于原点的镜像对称，通常选择g0＝s作为初始值。然后可以将传统基于余弦模型的傅里叶成分分析(FCA，Fourier Component Analysis)-FCA方法提取出的吸收(A)、相衬(P)和暗场(D)三种对比度信息表示为小角散射分布的零阶矩M0、归一化一阶矩和归一化二阶中心矩
Figure BDA0002267856810000104
A→M0(g)＝∫g(φ)dφ (3)
Figure BDA0002267856810000105
Figure BDA0002267856810000106
其中，Mn(g)表示g(φ)原始的n阶矩，
Figure BDA0002267856810000107
表示归一化后的n阶矩，
Figure BDA0002267856810000108
表示归一化后的n阶中心矩。
DB-MMA优势在于不需要对位移曲线进行模型假设，并且获取的图像噪声水平较低，而且不会出现相位卷绕的问题，还可以计算出三阶矩及以上的更丰富信息，但是限制因素在于解卷积迭代过程计算时间较长，并且需要手动设置最优的迭代次数，以及迭代可能导致部分结构细节的损失等。
3、D-MMA法
在DB-MMA基础上，D-MMA也从公式(1)中的卷积关系出发，根据卷积的傅里叶变换性质可以得到：
Figure BDA0002267856810000109
其中，和
Figure BDA00022678568100001011
分别是s(φ)、f(φ)和g(φ)的傅里叶变换，ω是傅里叶空间中的变量。然后根据多阶矩的傅里叶变换性质，可以得到：
Figure BDA0002267856810000111
由公式(6)可以得到
Figure BDA0002267856810000112
的零阶、一阶和二阶导数为：
Figure BDA0002267856810000113
Figure BDA0002267856810000114
Figure BDA0002267856810000115
将公式(8)-(10)带入公式(7)中进行化简以及归一化等可以获得与吸收、相衬和暗场相对应的多阶矩表达式为：
Figure BDA0002267856810000116
Figure BDA0002267856810000118
D-MMA相比于DB-MMA最大的优势就在于避免了复杂而又耗时的解卷积迭代过程，并且没有引入更多的模型假设，同样可以获得相同的多阶矩信息。
4、GD-MMA法
特别的，在基于光栅的X射线相衬成像(GI)中，根据波动光学中的菲涅尔衍射等理论以及几何光学中的莫尔偏折法等理论，研究者发现在只考虑±1阶衍射的情况下，背景位移曲线可以看成是余弦函数模型，并且在实际实验中也发现该模型与实际数据吻合得十分好，因此在GI系统中，可以将背景位移曲线表示为：
Figure BDA0002267856810000119
基于以上表达式和三角函数的正交性，申请人的中国专利申请号201910593224.0发现了两种位移曲线均满足一种积分特性，即：位移曲线和任意一个函数的内积只包含其傅里叶级数中的零阶项和一阶项。从上面介绍的D-MMA方法中，可以发现D-MMA主要是计算f(φ)或者s(φ)与另一个核函数h(φ)(比如对于相衬是hP(φ)＝φ，对于暗场是hD(φ)＝φ2)的内积。将位移曲线的这种积分特性应用其中，中国专利申请号201910593224.0建立了一种广义的D-MMA方法，也就是GD-MMA。在GD-MMA中，分别对应于相衬和暗场信息计算过程中的核函数为：
此外，中国专利申请号201910593224.0还找出了一种针对噪声抑制的最优GD-MMA核函数，即：
Figure BDA0002267856810000123
Figure BDA0002267856810000124
5、TA-MMA法
申请人的中国专利申请号201910207542.9提出了另一种解析MMA类信息提取算法TA-MMA，其提取信息方法如下所示：
Figure BDA0002267856810000125
Figure BDA0002267856810000126
Figure BDA0002267856810000127
其中，M0为小角散射分布(散射角分布)的零阶矩；M1为小角散射分布的一阶矩；M2为角散射分布的二阶矩。
其中，SinM1(·)和CosM1(·)定义如下，
SinM1(y)＝∫sin(φ)y(φ)dφ (22)
CosM1(y)＝∫cos(φ)y(φ)dφ (23)
(二)MMA类信息提取算法中各步X射线剂量遵循的剂量分布
MMA方法的前提是采集到了背景位移曲线和物体位移曲线，通常的采集方式是均匀间隔地进行机械步进，并且保证每一步的X射线剂量是相同的(即X射线源电压、电流、曝光时间等保持不变)，本公开的改变各步剂量的方式需要克服对于信息提取方法造成的不利影响同时还需要辅助该信息提取方法降低数据噪声。
根据剂量研究的观点，在其它条件不变的情况下，X射线的剂量正比于X射线源的电流以及曝光时间，也与X射线源的电压有一定的关系，当然也正比于X射线源的光强，从而根据比尔定律正比于探测器接收到的光强。因此，为了比较明确地调节每一步的X射线剂量，可以改变X射线源的电流或者曝光时间。
1、GD-MMA信息提取算法和D-MMA信息提取算法
下面以在Talbot-Lau类型的GI系统下，采用GD-MMA信息提取算法提取出相衬信息进行分析，讨论在数据采集过程总体剂量不变的情况下，改变每一步的剂量来优化最终提取出的相衬信息噪声水平。当然，本实施例中，GD-MMA中的核函数为D-MMA中核函数的拓宽形式，该推导同样适用于D-MMA法。
从实际操作简单考虑，仅以改变每一步的曝光时间来改变剂量，当然，改变曝光时间或者其他因素也是可以的，只要能够调节X射线的剂量即可。假设采集总步数为N，每一步的曝光时间为t(n)，其中n＝1，...，N，那么在GD-MMA方法中，根据公式(12)中的相衬信息表达式，来求解满足如下条件的优化问题：所述特定剂量分布使得总剂量保持固定的情况下(对应公式(25)的约束条件)相衬(本实施例中对应公式(24)中的相衬信息表达式)或者暗场信息的方差最小(对应公式(24))且均值不变(对应公式(26))，由于采集背景位移曲线f(φ)时没有物体，因此不需要考虑这部分对于扫描物体造成的剂量，在公式(12)中可以不考虑第二项。因此可以将上述优化问题表示为如下形式：
Figure BDA0002267856810000131
Figure BDA0002267856810000132
Figure BDA0002267856810000133
其中，Var(·)表示变量的方差；min表示最小值函数；
Figure BDA0002267856810000141
表示约束条件为总剂量保持固定；Mean表示均值；hP(n)表示解析方法中的核函数，在公式(12)中，hP(n)＝φ；s(n)为各步对应的曝光时间t(n)，其中n＝1，...，N下的物体位移曲线表达式。
其中
Figure BDA0002267856810000142
可以与公式(12)中第一项对应。
另一方面，按照现有的数据采集方案，各步对应的剂量是相等的，均等于总剂量除以总步数，即t0(n)＝ΔT＝T/N，据此得到各步等剂量条件下的物体位移曲线s0(n)，并且通常认为s0(n)和背景位移曲线一样也满足余弦模型，参照公式(14)的形式可知，s0(n)可以表示为：
根据剂量正比于曝光时间，同时正比于探测器接收到的光强，可以得到任意曝光时间分布t(n)下的物体位移曲线表达式s(n)，即：
Figure BDA0002267856810000144
对于公式(24)中的分母而言，可以考虑成
Figure BDA0002267856810000146
和
Figure BDA0002267856810000147
两顶相加。由于通常情况下系统对比度且
Figure BDA0002267856810000149
存在一半为负数，可近似认为
Figure BDA00022678568100001410
从而可得：
Figure BDA00022678568100001411
此外，再考虑到探测器的统计噪声，通常认为s(n)中的噪声满足泊松分布，从而可表示为：
Figure BDA00022678568100001412
其中，
Figure BDA00022678568100001413
是一个标准正态分布随机变量。将公式(28)(29)(30)代入，公式(24)中的目标函数可以化简为：
Figure BDA00022678568100001414
为了让约束/限制条件(25)满足，即均值不变，需要将核函数除掉曝光时间分布t(n)，从而得到目标函数如下：
Figure BDA0002267856810000151
从上式(32)可以看出括号内第一项为均值，与t(n)无关保持不变，第二项为方差来源，因此可以将最终的优化问题表示为：
Figure BDA0002267856810000152
Figure BDA0002267856810000153
利用拉格朗日乘子法来进行优化求解，首先得到拉格朗日函数为：
其中，λ为常数。然后分别对t(n)，n＝1，...，N和λ求偏导并令其为0：
Figure BDA0002267856810000156
由第一个方程(36)可得：
从而可以得到：
Figure BDA0002267856810000158
将公式(39)带入第二个方程(37)可以得到最终的最优曝光时间分布表达式为：
Figure BDA0002267856810000159
以上就是针对D-MMA和GD-MMA类信息提取算法的相衬信息进行优化的最优曝光时间分布，或者说最优剂量分布，即为上述介绍的特定剂量分布。
从上述公式(40)可知，该实例中，所述特定剂量分布与信息提取算法中的核函数相关。
在某些应用场景下，s0(n)不易获得，则可以近似认为s0(n)＝常数，从而得出以下近似分布，即：在一实例中，所述特定剂量分布中各步的X射线剂量与信息提取算法中的核函数成正比，特定剂量分布以核函数作为权重。仅通过算法中的核函数即可确定数据采集时的曝光时间：
Figure BDA0002267856810000161
2、TA-MMA信息提取算法
此外，针对TA-MMA算法的剂量分布优化，同样可以利用公式(40)得出的结论，只是其对应于相衬信息和暗场信息的核函数根据公式(20)和(21)获得，例如对于相衬信息的核函数可表示为如下所示：
hP(φ)＝CosM1(f)sin(φ)-SinM1(f)cos(φ) (42)
由上述内容可知，当所述信息提取算法为D-MMA、GD-MMA、或TA-MMA时，所述特定剂量分布与信息提取算法中的核函数相关。
在一实例中，当所述信息提取算法为D-MMA、GD-MMA或TA-MMA时，所述特定剂量分布中各步的X射线剂量与信息提取算法中的核函数成正比，特定剂量分布以核函数作为权重。
3、DB-MMA信息提取算法
此外，针对于DB-MMA算法，虽然其效率较低但是适用的范围更加广泛，因此下面通过噪声传播的方法推导了其剂量分布的最优形式。下面的公式中，采用加粗的字体表示矩阵参量，除非另有定义的除外，例如，
Figure BDA0002267856810000162
为一个数值。
首先将DB-MMA的卷积假设关系公式(1)用矩阵形式可以表示为：
Figure BDA0002267856810000163
其中，s表示物体位移曲线s(φ)的矩阵形式；Fg表示背景位移曲线f(φ)和小角散射分布g(φ)的卷积的矩阵形式；g表示小角散射分布g(φ)的矩阵形式；
s＝[s(0)，s(1)，...，s(N-1)]T，g＝[g(0)，g(1)，...，g(N-1)]T，而
Figure BDA0002267856810000171
是一个循环矩阵，该循环矩阵中各个元素为背景位移曲线。根据该关系可以通过最小二乘噪声传播关系得到g的协方差矩阵Vg：
其中，Vs表示s的协方差矩阵；(矩阵)T表示矩阵的转置，(矩阵)-1表示逆矩阵；
然后，将DB-MMA信息提取公式(3)-(5)同样用矩阵形式表示为：
Figure BDA0002267856810000173
其中，M表示三种对比度信息：小角散射分布的零阶矩M0、归一化一阶矩
Figure BDA0002267856810000174
阳归一化二阶中心矩
Figure BDA0002267856810000175
的矩阵形式；
M＝[M0(g)，M1(g)，M2(g)]T，雅可比矩阵
Figure BDA0002267856810000176
其中，上述雅可比矩阵中，第一行、第二行、第三行中的各个元素分别为对比度信息对应的核函数；根据该关系，同样通过最小二乘噪声传播关系可得到M的协方差矩阵VM：
Figure BDA0002267856810000177
当对三种对比度信息里面的第i种信息进行优化时，可以将三种对比度信息其中一个的方差表示为：
Figure BDA0002267856810000178
Figure BDA0002267856810000179
为矩阵对角线上的元素，其中，ei是在第i个位置为1其它均为0的三维单位向量，i＝1，2，3。然后对该方差优化，可以通过对其求导，令导数为0获得，即：
Figure BDA0002267856810000181
其中，将eiRF-1简记为x，即x＝eiRF-1＝[x0，x1，...，xN-1]。通常情况下，每一步的噪声满足泊松分布，从而与曝光时间成正比，即可得到：
根据公式(48)(49)，可以得到：
Figure BDA0002267856810000183
因此，只要满足公式(50)的剂量分布即为DB-MMA的最优分布，其中，t(n)满足公式(25)。由公式(50)可知，该剂量分布与R有关，即与信息提取算法中的核函数有关，另外，公式(50)中x＝eiRF-1＝[x0，x1，...，xN-1]也说明该剂量分布还与背景位移曲线的数值有关，基于公式(50)得到剂量分布在非均匀化的情况下能够实现提取出的相衬或者暗场信息数值不变但噪声最小。
由上述内容可知，当所述信息提取算法为DB-MMA时，所述特定剂量分布与信息提取算法中的核函数相关。
为了更直观地对比以上得出的最优曝光时间分布和传统分布的区别，
图4为采用本公开一实施例所示的信息提取方法中的最优曝光时间分布与传统方式对比的结果示意图。图4中，D-MMA和GD-MMA两种核函数优化后的曝光时间分布与传统的均匀分布对比，总体剂量保持不变。D-MMA中采用hP(φ)＝φ加权，GD-MMA噪声最优形式中采用hP(φ)＝sin(φ)优化，根据公式(40)得出的曝光时间分布进行对比。
图5为采用D-MMA算法处理模拟数据情况，(a)均匀分布的曝光时间和(b)优化后的加权分布曝光时间采集方式提取的相衬信息剖面线对比图。
以D-MMA中hP(φ)＝φ加权为例，图5中(a)和(b)分布展示出在处理模拟数据的情况下，均匀分布曝光时间和加权分布曝光时间的相衬信息剖面线，对比图5中(a)和(b)，可以明显看出相同剂量下以核函数加权分布的曝光时间产生的噪声水平更低。
为了更近一步地定量对比，在模拟数据N＝19的条件下计算了公式(33)中的目标函数，分别在传统的均匀分布和D-MMA核函数加权优化后的曝光时间值，进而计算出理论上可以降低的噪声比例为：
Figure BDA0002267856810000191
然后通过改变不同的光子数，模拟出在不同泊松噪声水平下两种数据采集方式实际产生的噪声水平，以及加权分布降低的噪声比例，如图6所示，可以看出基本在公式(51)给出的理论值附近，证明了以上推导过程的合理性。
第二实施例
在本公开的第二个示例性实施例中，提供了一种X射线相衬成像系统。本实施例的X射线相衬成像系统中各步的X射线剂量遵循一特定剂量分布。本公开的X射线相衬成像系统也可以采用本公开提及的基于剂量分布的信息提取方法实现三种对比度信息提取。
所述X射线相衬成像系统以机械步进方式或等同于机械步进的方式获得背景位移曲线和物体位移曲线；采用一信息提取算法分析背景位移曲线和物体位移曲线，以提取出物体的相衬或暗场信息；其中，所述机械步进方式或所述等同于机械步进的方式中，各步的X射线剂量不是均匀的，各步的X射线剂量遵循一特定剂量分布，所述特定剂量分布使得总剂量保持固定的情况下相衬或者暗场信息的方差最小且均值不变。
在本公开的一实施例中，所述特定剂量分布与信息提取算法中的核函数相关。例如，信息提取算法为D-MMA、GD-MMA、TA-MMA或者GB-MMA。
在本公开的一实施例中，所述特定剂量分布与信息提取算法中的核函数成正比。信息提取算法为D-MMA、GD-MMA或者TA-MMA。
在本公开的一实施例中，当所述信息提取算法为DB-MMA时，所述特定剂量分布满足：
T为常数
Figure BDA0002267856810000201
x＝eiRF-1＝[x0，x1，...，xN-1]
Figure BDA0002267856810000202
其中，t(n)表示各步的剂量；n＝1，...，N，N为总步数；ei表示在第i个位置为1其它均为0的三维单位向量；R表示三种对比度信息的矩阵表达式中的雅可比矩阵，上述雅可比矩阵中，各个元素分别为对比度信息对应的核函数；F是一个循环矩阵，该循环矩阵中各个元素为背景位移曲线。
综上所述，本公开提供了一种基于剂量分布的信息提取方法及X射线相衬成像系统，首次提出通过变化步进扫描方式中每一步的X射线剂量，使得各步X射线的剂量呈现一定规律的分布，所述特定剂量分布使得总剂量保持固定的情况下相衬或者暗场信息的方差最小且均值不变，该分布降低了数据的噪声，可以在相同总剂量的情况下得到更低噪声水平的对比度信息；并可以在保证相同噪声水平的对比度信息条件下，降低数据采集过程中的总体剂量水平；该变化剂量分布的求解方法适用于多种X射线相衬成像系统，例如为ABI、EI和GI，并且GI系统类别不限，即包括Talbot-Lau类型、几何投影类型和双相位光栅类型等各种类型的系统，适用范围广泛。
单词“包含”或“包括”不排除存在未列在权利要求中的元件或步骤。位于元件之前的单词“一”或“一个”不排除存在多个这样的元件。此外，除非特别描述或必须依序发生的步骤，上述步骤的顺序并无限制于以上所列，且可根据所需设计而变化或重新安排。并且上述实施例可基于设计及可靠度的考虑，彼此混合搭配使用或与其他实施例混合搭配使用，即不同实施例中的技术特征可以自由组合形成更多的实施例。
以上所述的具体实施例，对本公开的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本公开的具体实施例而已，并不用于限制本公开，凡在本公开的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本公开的保护范围之内。

Claims (10)

1.一种基于剂量分布的信息提取方法，用于X射线相衬成像系统，其特征在于，所述信息提取方法包括：
以机械步进方式或等同于机械步进的方式在设置物体和不设置物体的情况下分别进行曝光；
X射线探测器在曝光后进行探测，获取不设置物体的背景位移曲线和设置物体的物体位移曲线；
采用一信息提取算法分析背景位移曲线和物体位移曲线，以提取出物体的相衬或暗场信息；
其中，所述机械步进方式或所述等同于机械步进的方式中，各步的X射线剂量不是均匀的，各步的X射线剂量遵循一特定剂量分布，所述特定剂量分布使得总剂量保持固定的情况下相衬或者暗场信息的方差最小且均值不变。
2.根据权利要求1所述信息提取方法，其特征在于，所述信息提取算法包括如下算法的一种或几种：D-MMA、GD-MMA、TA-MMA以及DB-MMA。
3.根据权利要求2所述的信息提取方法，其特征在于，当所述信息提取算法为D-MMA、GD-MMA、TA-MMA或DB-MMA时，所述特定剂量分布与所述信息提取算法中的核函数相关。
4.根据权利要求3所述的信息提取方法，其特征在于，当所述信息提取算法为D-MMA、GD-MMA或TA-MMA时，所述特定剂量分布中各步的X射线剂量与所述信息提取算法中的核函数成正比，特定剂量分布以核函数作为权重。
5.根据权利要求3所述的信息提取方法，其特征在于，当所述信息提取算法为DB-MMA时，所述特定剂量分布满足：
Figure FDA0002267856800000011
T为常数
x＝eiRF-1＝[x0，x1，...，xN-1]
Figure FDA0002267856800000022
其中，t(n)表示各步的剂量；n＝1，...，N，N为总步数；ei表示在第i个位置为1其它均为0的三维单位向量；R表示三种对比度信息的矩阵表达式中的雅可比矩阵，上述雅可比矩阵中，各个元素分别为对比度信息对应的核函数；F是一个循环矩阵，该循环矩阵中各个元素为背景位移曲线。
6.根据权利要求1所述的信息提取方法，其特征在于，所述各步的X射线剂量通过调整每一步的X射线源的电流、曝光时间、X射线源的电压或者X射线源的光强中的一个或多个来控制各步的X射线剂量，从而使得各步的X射线剂量遵循特定剂量分布。
7.根据权利要求1所述的信息提取方法，其特征在于，
提供所述X射线的X射线源为如下光源中的一种：X光源、微焦点X光源、以及同步辐射X光源；
所述X射线探测器包括如下探测器中的一种：能量积分型探测器或者光子计数型探测器。
8.根据权利要求1所述的信息提取方法，其特征在于，所述X射线相衬成像系统包括如下成像系统中的一种：基于分析晶体的成像系统、基于光栅的成像系统以及边缘照射成像系统；
其中，所述基于光栅的成像系统包括Talbot-Lau类型、几何投影类型和双相位光栅类型中的一种。
9.一种X射线相衬成像系统，其特征在于，以机械步进方式或等同于机械步进的方式获得背景位移曲线和物体位移曲线；采用一信息提取算法分析背景位移曲线和物体位移曲线，以提取出物体的相衬或暗场信息；其中，所述机械步进方式或所述等同于机械步进的方式中，各步的X射线剂量不是均匀的，各步的X射线剂量遵循一特定剂量分布，所述特定剂量分布使得总剂量保持固定的情况下相衬或者暗场信息的方差最小且均值不变。
10.根据权利要求9所述的X射线相衬成像系统，其特征在于，所述信息提取算法包括如下算法的一种或几种：D-MMA、GD-MMA、TA-MMA以及DB-MMA；
可选的，当所述信息提取算法为D-MMA、GD-MMA、TA-MMA或DB-MMA时，所述特定剂量分布与所述信息提取算法中的核函数相关；
进一步可选的，当所述信息提取算法为D-MMA、GD-MMA或TA-MMA时，所述特定剂量分布中各步的X射线剂量与所述信息提取算法中的核函数成正比，特定剂量分布以核函数作为权重；
进一步可选的，当所述信息提取算法为DB-MMA时，所述特定剂量分布满足：
Figure FDA0002267856800000031
T为常数
x＝eiRF-1＝[x0，x1，...，xN-1]
Figure FDA0002267856800000033
其中，t(n)表示各步的剂量；n＝1，...，N，N为总步数；ei表示在第i个位置为1其它均为0的三维单位向量；R表示三种对比度信息的矩阵表达式中的雅可比矩阵，上述雅可比矩阵中，各个元素分别为对比度信息对应的核函数；F是一个循环矩阵，该循环矩阵中各个元素为背景位移曲线。
CN201911098795.3A 2019-11-11 2019-11-11 基于剂量分布的信息提取方法及x射线相衬成像系统 Active CN110806598B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911098795.3A CN110806598B (zh) 2019-11-11 2019-11-11 基于剂量分布的信息提取方法及x射线相衬成像系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911098795.3A CN110806598B (zh) 2019-11-11 2019-11-11 基于剂量分布的信息提取方法及x射线相衬成像系统

Publications (2)

Publication Number Publication Date
CN110806598A true CN110806598A (zh) 2020-02-18
CN110806598B CN110806598B (zh) 2021-06-25

Family

ID=69502140

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911098795.3A Active CN110806598B (zh) 2019-11-11 2019-11-11 基于剂量分布的信息提取方法及x射线相衬成像系统

Country Status (1)

Country Link
CN (1) CN110806598B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111795980A (zh) * 2020-08-04 2020-10-20 合肥工业大学 一种基于逐像素高斯函数拟合法的x射线边界照明成像方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2007049464A1 (ja) * 2005-10-24 2007-05-03 Konica Minolta Medical & Graphic, Inc. Ｘ線画像撮影システム及びｘ線画像撮影方法
CN101756707A (zh) * 2009-12-31 2010-06-30 苏州和君科技发展有限公司 用Micro-CT成像系统对长目标物体进行扫描重建的方法
CN102579066A (zh) * 2012-02-17 2012-07-18 天津大学 一种x射线同轴相衬成像方法
EP3021110A1 (en) * 2014-11-11 2016-05-18 Paul Scherrer Institut System for obtaining quantitative x-ray images using Hilbert transform on imaged fringes
CN109975334A (zh) * 2019-04-25 2019-07-05 兰州大学 一种单次曝光的x射线二维相衬成像方法
CN110095481A (zh) * 2019-05-24 2019-08-06 清华大学 X射线光栅成像系统与成像方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2007049464A1 (ja) * 2005-10-24 2007-05-03 Konica Minolta Medical & Graphic, Inc. Ｘ線画像撮影システム及びｘ線画像撮影方法
CN101756707A (zh) * 2009-12-31 2010-06-30 苏州和君科技发展有限公司 用Micro-CT成像系统对长目标物体进行扫描重建的方法
CN102579066A (zh) * 2012-02-17 2012-07-18 天津大学 一种x射线同轴相衬成像方法
EP3021110A1 (en) * 2014-11-11 2016-05-18 Paul Scherrer Institut System for obtaining quantitative x-ray images using Hilbert transform on imaged fringes
CN109975334A (zh) * 2019-04-25 2019-07-05 兰州大学 一种单次曝光的x射线二维相衬成像方法
CN110095481A (zh) * 2019-05-24 2019-08-06 清华大学 X射线光栅成像系统与成像方法

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111795980A (zh) * 2020-08-04 2020-10-20 合肥工业大学 一种基于逐像素高斯函数拟合法的x射线边界照明成像方法
CN111795980B (zh) * 2020-08-04 2022-04-26 合肥工业大学 一种基于逐像素高斯函数拟合法的x射线边界照明成像方法

Also Published As

Publication number Publication date
CN110806598B (zh) 2021-06-25

Similar Documents

Publication Publication Date Title
US7187794B2 (en) Noise treatment of low-dose computed tomography projections and images
Wunderlich et al. Image covariance and lesion detectability in direct fan-beam x-ray computed tomography
Lu et al. Noise properties of low-dose CT projections and noise treatment by scale transformations
Mumcuoglu et al. Accurate geometric and physical response modelling for statistical image reconstruction in high resolution PET
US9406154B2 (en) Iterative reconstruction in image formation
Delaney et al. A fast and accurate Fourier algorithm for iterative parallel-beam tomography
Burger et al. EM-TV methods for inverse problems with Poisson noise
US8204173B2 (en) System and method for image reconstruction by using multi-sheet surface rebinning
Xu et al. Statistical projection completion in X-ray CT using consistency conditions
US20140363067A1 (en) Methods and systems for tomographic reconstruction
CN110806598B (zh) 基于剂量分布的信息提取方法及x射线相衬成像系统
US8184876B2 (en) NNLS image reconstruction
Hutton et al. Iterative reconstruction methods
Sotthivirat et al. Image recovery using partitioned-separable paraboloidal surrogate coordinate ascent algorithms
Humphries et al. Comparison of regularized and superiorized methods for tomographic image reconstruction
O'Sullivan A study of least squares and maximum likelihood for image reconstruction in positron emission tomography
Deift et al. Universality for eigenvalue algorithms on sample covariance matrices
US8014580B2 (en) Determining a pixon map for image reconstruction
Nuyts et al. The SNR of positron emission data with Gaussian and non-Gaussian time-of-flight kernels, with application to prompt photon coincidence
US10515467B2 (en) Image reconstruction system, method, and computer program
Hamilton et al. A simulation of QCD radiation in top quark decays
Thompson et al. Inference of non-thermal electron energy distributions from hard X-ray spectra
CN110197486B (zh) X射线光栅综合成像信息提取方法、系统及存储介质
Brune et al. Forward backward EM-TV methods for inverse problems with Poisson noise
Wang et al. Prospective prediction and control of image properties in model-based material decomposition for spectral CT

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant