CN112884853B - 一种基于相位ct的材料分解与识别方法 - Google Patents

一种基于相位ct的材料分解与识别方法 Download PDF

Info

Publication number
CN112884853B
CN112884853B CN202110019682.0A CN202110019682A CN112884853B CN 112884853 B CN112884853 B CN 112884853B CN 202110019682 A CN202110019682 A CN 202110019682A CN 112884853 B CN112884853 B CN 112884853B
Authority
CN
China
Prior art keywords
image
phase
substrate
projection
absorption
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.)
Active
Application number
CN202110019682.0A
Other languages
English (en)
Other versions
CN112884853A (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.)
Capital Normal University
Original Assignee
Capital Normal 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 Capital Normal University filed Critical Capital Normal University
Priority to CN202110019682.0A priority Critical patent/CN112884853B/zh
Publication of CN112884853A publication Critical patent/CN112884853A/zh
Application granted granted Critical
Publication of CN112884853B publication Critical patent/CN112884853B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T11/00—2D [Two Dimensional] image generation
    • G06T11/003—Reconstruction from projections, e.g. tomography
    • G06T11/005—Specific pre-processing for tomographic reconstruction, e.g. calibration, source positioning, rebinning, scatter correction, retrospective gating
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/70—Denoising; Smoothing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10072—Tomographic images
    • G06T2207/10081—Computed x-ray tomography [CT]

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Analysing Materials By The Use Of Radiation (AREA)

Abstract

本发明公开了一种基于相位CT的材料分解与识别方法，所述方法包括：步骤1，令m＝0,初始化第一基材图像f(m)和第二基材图像g(m)；步骤2，根据式(6)计算离散化吸收投影数据
Figure DDA0002888216130000011
的第m轮迭代中的估计值
Figure DDA0002888216130000012
和对应的残差
Figure DDA0002888216130000013
以及离散化微分相位投影数据
Figure DDA0002888216130000014
在第m轮迭代中的估计值
Figure DDA0002888216130000015
和对应的残差
Figure DDA0002888216130000016
步骤3，根据式(5)迭代求解第一基材图像f(m)在第m+1轮迭代中的值f(m+1)和第二基材图像g(m)在第m+1轮迭代中的值g(m+1)；步骤4，对f(m+1)，g(m+1)进行正则化操作；步骤5，若停止条件不满足，令m＝m+1，并返回步骤2。本发明方法联合吸收投影信息及微分相位投影信息进行基材分解，得到的基材图像噪声更小，后续计算的等效原子序数图像的噪声也比常规的更小。

Description

一种基于相位CT的材料分解与识别方法
技术领域
本发明涉及材料识别技术领域，特别是关于一种基于相位CT的材料分解与识别方法。
背景技术
X射线成像技术被广泛应用在很多领域，例如生物医学、工业和安全检查等等。传统X射线吸收CT是基于不同物质或者组织之间吸收辐射的差异性进行成像，对于轻元素构成的材料(如肌肉、血管和乳腺等)成像不够理想。而X射线相位衬度CT是基于物质对穿过的X射线相位的改变而成像，通常其对轻物质的成像比吸收CT有较大的优势。X射线相衬成像方法里常用的是光栅相衬方法和衍射增强法，这两种方法都能同时提供吸收和相位信息。利用吸收信息和相位信息进行的物质分解和识别是相位CT的一个重要应用，在低能范围内，其效果比传统的利用双能CT进行物质分解和识别的方法更好。
目前，常用的基于相位CT的材料分解和识别技术主要是：1)基于投影域的分解技术；2)基于重建图像域的分解技术。由于采集到的数据的相位信息是相位投影的微分，需要经过积分运算才能得到相位投影。这个积分过程容易引起整体的低频噪声，使得计算的相位投影不准确，进而导致在投影域的分解误差扩大，所以在相位CT的材料分解重建中，通常以图像域的分解为主。另外，这两种方法将重建和分解的步骤分离，得到的分解图像的质量并不是太好，会影响到后续的定量分析。因此，为提高分解图像的质量，发展同时利用吸收和相位信息的重建分解算法是有意义的。
基于重建图像域的分解技术方法是基于吸收图像和相位图像可通过两种基材图像线性组合表示，在分别重建出吸收图像和相位图像后，选择所需要的两个基材，查找当前等效能量下两基材的吸收和相位信息作为线性组合的系数，然后求解得到两种基材的分布函数图像。
但是，在图像域中进行一次性的分解，存在噪声放大问题。分解的结果会受到重建图像的影响，而重建图像的质量主要由初始的投影数据的噪声水平和重建算法的选取决定，没有联合吸收及相位的信息进行重建和分解，得到的结果欠佳。
发明内容
本发明的目的在于提供一种基于相位CT的材料分解与识别方法来克服或至少减轻现有技术的上述缺陷中的至少一个。
为实现上述目的，本发明提供一种基于相位CT的材料分解与识别方法，所述方法包括：
步骤1，令m＝0,初始化第一基材图像f(m)和第二基材图像g(m)；
步骤2，根据式(6)计算离散化吸收投影数据
Figure GDA0003631647890000021
的第m轮迭代中的估计值
Figure GDA0003631647890000022
和对应的残差
Figure GDA0003631647890000023
以及离散化微分相位投影数据
Figure GDA0003631647890000024
在第m轮迭代中的估计值
Figure GDA0003631647890000025
和对应的残差
Figure GDA0003631647890000026
步骤3，根据式(5)迭代求解第一基材图像f(m)在第m+1轮迭代中的值f(m+1)和第二基材图像g(m)在第m+1轮迭代中的值g(m+1)；
步骤4，对f(m+1)，g(m+1)进行正则化操作；
步骤5，若停止条件不满足，令m＝m+1，并返回步骤2；
Figure GDA0003631647890000027
Figure GDA0003631647890000028
其中：
Figure GDA0003631647890000031
表示第一基材的fj和第二基材的gj在投影角度
Figure GDA0003631647890000032
下对第u条射线路径的贡献；
Figure GDA0003631647890000033
表示第一基材的fj在第m+1轮迭代中的值；
Figure GDA0003631647890000034
表示第一基材的fj在第m轮迭代中的值；
Figure GDA0003631647890000035
表示第二基材的gj在第m+1轮迭代中的值；
Figure GDA0003631647890000036
表示第二基材的gj在第m轮迭代中的值；
λ表示松弛因子；
Figure GDA0003631647890000037
表示离散化吸收投影数据
Figure GDA0003631647890000038
的第m轮迭代中的估计值；
Figure GDA0003631647890000039
表示离散化微分相位投影数据
Figure GDA00036316478900000310
在第m轮迭代中的估计值；
Figure GDA00036316478900000311
表示第m轮迭代中的吸收投影残差
Figure GDA00036316478900000312
的第u个分量；
Figure GDA00036316478900000313
表示第m轮迭代中的微分相位投影残差
Figure GDA00036316478900000314
的第u个分量。
进一步地，采用式(7)表示的第二优化模型求解
Figure GDA00036316478900000315
的最优值
Figure GDA00036316478900000316
Figure GDA00036316478900000317
式中，
Figure GDA00036316478900000318
表示
Figure GDA00036316478900000319
和当前相位投影残差
Figure GDA00036316478900000320
的二范数距离最小时的相位投影残差
Figure GDA00036316478900000321
进一步地，离散化吸收投影数据
Figure GDA00036316478900000322
和离散化微分相位投影数据
Figure GDA00036316478900000323
表示式(3)：
Figure GDA00036316478900000324
其中：
Figure GDA00036316478900000325
表示投影角度
Figure GDA00036316478900000326
下不同探测器接收到的离散化吸收投影数据，其中的元素
Figure GDA00036316478900000327
表示第u个探测器像素采集的吸收投影值，U表示探测器的总像素个数，τ表示向量的转置；
Figure GDA00036316478900000328
表示投影角度
Figure GDA00036316478900000329
下的投影矩阵，其中的J表示每幅图像的总像素个数；
f＝(f1,f2...,fJ)τ表示第一基材的分布函数图像f(x,y)的离散化图像，其中的元素fj为f(x,y)在被测样品上第j个像素上的采样值；
g＝(g1,g2...,gJ)τ表示第二基材的分布函数图像g(x,y)的离散化图像，其中的元素gj为g(x,y)在被测样品上第j个像素上的采样值；
Figure GDA0003631647890000041
表示投影角度
Figure GDA0003631647890000042
下不同探测器的离散化微分相位投影数据，其中的元素
Figure GDA0003631647890000043
表示第u个探测器像素采集的微分相位投影值。
进一步地，通过求解如下式(4)提供的第一优化模型，得到最后的分解图像：
Figure GDA0003631647890000044
其中：
f*表示第一优化模型取得最小时第一基材的分布函数图像f；
g*表示第一优化模型取得最小时第二基材的分布函数图像g；
R(·)是正则化算子。
进一步地，所述方法还包括：
步骤6，根据基材图像f(m)和g(m)，利用公式(1)得到吸收图像和相位图像；
Figure GDA0003631647890000045
其中：
μ(x,y)表示被测样品的点(x,y)处的吸收图像；
δ(x,y)表示被测样品的点(x,y)处的相位图像；
μ1为第一基材的线性吸收图像；
μ2为第二基材的线性吸收图像；
δ1为第一基材的相位图像；
δ2为第二基材的相位图像；
f(x,y)为第一基材的分布函数图像；
g(x,y)为第二基材的分布函数图像；
步骤7，利用式(9)求解等效原子序数：
Figure GDA0003631647890000051
其中：
Cp、CE、CZ是通过扫描与被测物体相近的模体的数据进行标定求解的参数常量；
Figure GDA0003631647890000052
r0是经典原子半径，h是普朗克常数，c是光速常数；
CKN(E)为Klein-Nishina截面；
μ(E,x,y)表示被测样品的点(x,y)处对能量E的射线的吸收图像；
δ(E,x,y)表示被测样品的点(x,y)处对能量E的射线的相位图像；
ρe(x,y)表示被测样品的点(x,y)处的电子密度；
Z(x,y)表示被测样品的点(x,y)处的等效原子序数。
由于本发明根据分解基材图像的需要，建立了直接由微分相位和吸收投影数据重建基材图像的、能抑制噪声的优化成像模型，结合图像迭代重建方法，提出了一种利用微分相位和吸收的物质分解方法，该方法建立在优化模型的基础上，采用的迭代算法可直接从吸收和微分相位的投影中重建基材成分，进一步地计算可得到等效原子序数，具有抑制噪声和分解准确的优点。
附图说明
图1为Image-based方法和本发明方法分解的各种基材图像结果对比示意图。
图2为使用双能采集和本发明方法计算的等效原子序数图像对比示意图。
具体实施方式
下面结合附图和实施例对本发明进行详细的描述。
根据材质分解方法，在单能的情况下，被测样品的线性吸收图像可以用两种基材的线性吸收图像近似线性表示。同理，被测样品的折射率实部衰减率也可以用两种基材的折射率实部衰减率近似线性表示。于是，可以列出下列方程组(1)：
Figure GDA0003631647890000061
其中：
μ(x,y)表示被测样品的点(x,y)处的吸收图像。
δ(x,y)表示被测样品的点(x,y)处的折射率实部衰减率系数，即文中的相位图像。
μ1为第一基材的线性吸收图像。
μ2为第二基材的线性吸收图像。
δ1为第一基材的相位图像。
δ2为第二基材的相位图像。
f(x,y)为第一基材的分布函数图像。
g(x,y)为第二基材的分布函数图像。
(x,y)表示样品坐标系中任一点的坐标，样品坐标系通常选取样品的几何中心为原点，X轴和Y轴相互垂直，其方向可根据实际需求进行选取。
基于光栅干涉仪的X射线微分相位成像通过光栅步进的采集方式，可以分离出每个投影角度
Figure GDA0003631647890000062
下样品对X射线的吸收和微分相移信息表示为式(2)：
Figure GDA0003631647890000063
其中：
Figure GDA0003631647890000064
为
Figure GDA0003631647890000065
坐标下接收到的吸收投影数据。
Figure GDA0003631647890000066
为
Figure GDA0003631647890000067
坐标下微分相位投影数据。
Figure GDA0003631647890000068
表示测量坐标系中任一点的坐标。
将式(2)可以离散化为式(3)：
Figure GDA0003631647890000071
其中：
Figure GDA0003631647890000072
表示投影角度
Figure GDA0003631647890000073
下不同探测器接收到的离散化吸收投影数据，其中的元素
Figure GDA0003631647890000074
表示第u个探测器像素采集的吸收投影值，U表示探测器的总像素个数，τ表示向量的转置。
Figure GDA0003631647890000075
表示投影角度
Figure GDA0003631647890000076
下的投影矩阵，其中的元素
Figure GDA0003631647890000077
表示fj和gj在投影角度
Figure GDA0003631647890000078
下对第u条射线路径的贡献，J表示每幅图像的总像素个数。
f＝(f1,f2...,fJ)τ表示f(x,y)的离散化图像，其中的元素fj为f(x,y)在被测样品上第j个像素上的采样值。
g＝(g1,g2...,gJ)τ表示g(x,y)的离散化图像，其中的元素gj为g(x,y)在被测样品上第j个像素上的采样值。
Figure GDA0003631647890000079
表示投影角度
Figure GDA00036316478900000710
下不同探测器的离散化微分相位投影数据，其中的元素
Figure GDA00036316478900000711
表示第u个探测器像素采集的微分相位投影值；D是微分算子
Figure GDA00036316478900000712
离散后的矩阵形式。
本发明实施例解决的是根据离散化吸收投影数据
Figure GDA00036316478900000713
和离散化微分相位投影数据
Figure GDA00036316478900000714
重建两种基材的分布函数图像f和g。
本发明实施例是通过求解如下式(4)提供的第一优化模型，得到最后的分解图像：
Figure GDA00036316478900000715
其中：
f*表示第一优化模型取得最小时第一基材的分布函数图像；
g*表示第一优化模型取得最小时第二基材的分布函数图像；
R(·)是正则化算子，如TV、L0或图像滤波操作；
||·||2表示二范数的平方。
也可以采用其它的模型得到最后分解的图像，比如：采用其它范数代替上述第一优化模型的二范数的平方，并通过在吸收项的范数前面加一个参数，控制吸收项和相位项的比重，使得由f和g计算的吸收项和微分相位项尽可能接近探测器探测到的吸收值和微分相位值(保真)。
在式(4)的保真项中，应用如下式(5)表示的迭代步骤：
Figure GDA0003631647890000081
其中：
Figure GDA0003631647890000082
Figure GDA0003631647890000083
表示第一基材的fj在第m+1轮迭代中的值。
Figure GDA0003631647890000084
表示第一基材的fj在第m轮迭代中的值。
Figure GDA0003631647890000085
表示第二基材的gj在第m+1轮迭代中的值。
Figure GDA0003631647890000086
表示第二基材的gj在第m轮迭代中的值。
λ表示松弛因子，其取值在(0，1)之间，可抑制噪声的影响。
Figure GDA0003631647890000087
表示J个
Figure GDA0003631647890000088
的总和。
Figure GDA0003631647890000089
表示U个
Figure GDA00036316478900000810
的总和。
Figure GDA00036316478900000811
表示离散化吸收投影数据
Figure GDA00036316478900000812
的第m轮迭代中的估计值。
Figure GDA00036316478900000813
表示离散化微分相位投影数据
Figure GDA00036316478900000814
在第m轮迭代中的估计值。
Figure GDA00036316478900000815
表示第m轮迭代中的吸收投影残差
Figure GDA00036316478900000816
的第u个分量。
Figure GDA0003631647890000091
表示第m轮迭代中的微分相位投影残差
Figure GDA0003631647890000092
的第u个分量。
计算相位投影残差
Figure GDA0003631647890000093
需要利用微分算子的逆，即D-1，而该算子对噪声有放大作用。为了克服噪声对重建图像的影响，可以采用式(7)表示的第二优化模型求解：
Figure GDA0003631647890000094
式中，
Figure GDA0003631647890000095
表示
Figure GDA0003631647890000096
和当前相位投影残差
Figure GDA0003631647890000097
的二范数距离最小时的相位投影残差
Figure GDA0003631647890000098
在另一个实施例中，可以采用式(8)表示的第三优化模型求解：
Figure GDA0003631647890000099
式中，
Figure GDA00036316478900000910
表示D(e)和当前相位投影残差
Figure GDA00036316478900000911
距离的二范数与e的正则化最小时的相位投影残差e。
在又一个实施例中，还可以采用距离范数和对变量e进行正则化约束，来确定相位投影残差
Figure GDA00036316478900000912
实际迭代中采用的是固定步数的梯度下降法来进行求解。当然，也可以直接求解或其它求解方式。
对于第一模型的正则化项，本实施例使用的策略是在对所有投影数据使用式(5)更新完f(m+1)，g(m+1)后，对其进行一次正则化操作。正则化算子主要选取TV算子，从而可以较好地抑制噪声的影响。
综上，本发明实施例的迭代步骤如下：
步骤1，初始化：f(0)＝0，g(0)＝0，m＝0。
步骤2，根据式(6)计算吸收投影的估计值
Figure GDA00036316478900000913
和对应的残差
Figure GDA00036316478900000914
步骤3，根据式(6)计算相位投影的估计值
Figure GDA00036316478900000915
由式(7)计算其残差
Figure GDA00036316478900000916
步骤4，根据式(5)迭代求解f(m+1)，g(m+1)。
步骤5，对f(m+1)，g(m+1)进行正则化操作。
步骤6，若停止条件不满足，令m＝m+1，并返回步骤2。其中，停止条件可以保真项的二范数是否小于某个值，或者设置固定的迭代步数。
步骤7，返回f(m)和g(m)。
相位CT及本发明方法的一个应用是对物质进行定量分析，利用上述方法得到基材图像f(m)和g(m)，代入式(1)计算吸收图像μ(x,y)和相位图像δ(x,y)，即式(9)中的μ(E,x,y)和δ(E,x,y)，再通过比如查表的方式获得等效能量E。接着，利用以下公式(9)可求解得到等效原子序数：
Figure GDA0003631647890000101
其中：
Cp、CE、CZ是通过扫描与被测物体相近的模体的数据进行标定求解的参数常量，比如：根据被测物体，选取基材相近的已知模体，那么，Z(x,y)、μ(E,x,y)和δ(E,x,y)均为已知量，再查表得到等效能量E，最后可以拟合出Cp、CE、CZ动具体数值。
Figure GDA0003631647890000102
r0是经典原子半径，h是普朗克常数，c是光速常数，各具体数值可通过查表获得。
CKN(E)为Klein-Nishina截面，通常使用公式(10)计算：
Figure GDA0003631647890000103
其中a＝E/511keV是相对电子质量，r0是经典原子半径。
μ(E,x,y)表示被测样品的点(x,y)处对能量E的射线的吸收图像。
δ(E,x,y)表示被测样品的点(x,y)处对能量E的射线的相位图像。
ρe(x,y)表示被测样品的点(x,y)处的电子密度。
Z(x,y)表示被测样品的点(x,y)处的等效原子序数。
此流程得到的吸收图像和相位图像是同时利用了吸收信息和相位信息，而并非原来单纯的由吸收信息得到吸收图像、相位信息得到相位图像。联合两种信息可减少原来单一信息直接计算吸收图像和相位图像的噪声影响，计算的等效原子序数图像的噪声也会更小。
下面是一组实验数据：在上海光源BL13W实验站上，使用基于光栅干涉仪的X射线微分相位成像装置上采集。光栅干涉仪由周期为2.396um的π/2相位光栅和周期为2.4um的吸收光栅构成，两块光栅之间的距离为46.38mm。实验所用的光子能量为20keV，采用有效像素阵列为2048×600的sCMOS探测器，像素尺寸为6.5μm。采集的数据是使用8步步进方法，即对相位光栅在一个周期内移动8步分别采集数据，完成一个角度的采集，而在180度内等间隔采集540个角度的数据。最后根据步进方法提取吸收和微分相位公式计算吸收投影数据和微分相位投影数据。
如图1所示，a是实物照片，模体包含四种成分：低密度聚乙烯(LDPE)、聚甲基丙烯酸甲酯(PMMA)、聚四氟乙烯(PTFE)和水。这几种成分被放在一个试管中的外直径10.7mm的聚乙烯塑料试管，PTFE，LDPE和PMMA的直径分别是2.0mm，4.0mm和5.6mm。
在图1中的b的吸收图像断层重建图中，圆环代表聚乙烯塑料容器，其中放置三种成分不同的小球，即圆环中间的三个灰度不同的圆面积代表三种成分的小球，左边最暗的是LDPE小球,下边最亮的是PTFE小球,而右上方的是PMMA小球，其余部分是水。
在图1中的c的折射率实部衰减率断层重建图中，LDPE和水无法区分开。这表明，利用相位信息重建的断层图像并不一定比利用吸收信息重建的断层图像的衬度高。
图1中的第二第三行是两种方法分解的PTFE基和水基结果，图1中的e和f分别是ImageBase和PAMD-SART方法分解的PTFE基图像，图1中的h和i分别是ImageBase和PAMD-SART方法分解的水基图像，图1中的d和g分别是两种方法下PTFE基和水基结果的剖线对比图。
定量分析的结果如图2所示，与相位CT对比的是使用双能采集计算的结果，在图上可看到相位CT计算的等效原子序数图像是要优于双能的。
图2，分解结果对比图，其中的a20keV下计算的原子序数断层图像，其中的b是12keV的虚拟吸收图像断层图像，其中的c是20keV的虚拟吸收图像断层图像，其中的d是双能计算的原子序数断层图像，其中的e是12keV的真实吸收图像断层图像，其中的f是20keV的真实吸收图像断层图像。
基于投影域的分解技术，基于重建图像域的分解技术和DEXA都可以实现材料分解和识别。但前两种方法属于两步法，存在噪声放大的问题。DEXA当被测样品由弱吸收的低原子序数物质构成时,其物质区分能力往往难以满足要求。
最后需要指出的是：以上实施例仅用以说明本发明的技术方案，而非对其限制。本领域的普通技术人员应当理解：可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (5)

1.一种基于相位CT的材料分解与识别方法，其特征在于，包括：
步骤1，令m＝0,初始化第一基材图像为f(m)和第二基材图像为g(m)；
步骤2，根据式(6)计算离散化吸收投影数据
Figure FDA0003631647880000011
的第m轮迭代中的估计值
Figure FDA0003631647880000012
和对应的残差
Figure FDA0003631647880000013
以及离散化微分相位投影数据
Figure FDA0003631647880000014
在第m轮迭代中的估计值
Figure FDA0003631647880000015
和对应的残差
Figure FDA0003631647880000016
步骤3，根据式(5)迭代求解第一基材图像f在第m+1轮迭代中的值f(m+1)和第二基材图像g在第m+1轮迭代中的值g(m+1)；
步骤4，对f(m+1)，g(m+1)进行正则化操作；
步骤5，若停止条件不满足，令m＝m+1，并返回步骤2；
Figure FDA0003631647880000017
Figure FDA0003631647880000018
其中：
Figure FDA0003631647880000019
表示第一基材的fj和第二基材的gj在投影角度
Figure FDA00036316478800000110
下对第u条射线路径的贡献；
Figure FDA00036316478800000111
表示第一基材的fj在第m+1轮迭代中的值；
Figure FDA00036316478800000112
表示第一基材的fj在第m轮迭代中的值；
Figure FDA0003631647880000021
表示第二基材的gj在第m+1轮迭代中的值；
Figure FDA0003631647880000022
表示第二基材的gj在第m轮迭代中的值；
λ表示松弛因子；
Figure FDA0003631647880000023
表示离散化吸收投影数据
Figure FDA0003631647880000024
的第m轮迭代中的估计值；
Figure FDA0003631647880000025
表示离散化微分相位投影数据
Figure FDA0003631647880000026
在第m轮迭代中的估计值；
Figure FDA0003631647880000027
表示第m轮迭代中的吸收投影残差
Figure FDA0003631647880000028
的第u个分量；
Figure FDA0003631647880000029
表示第m轮迭代中的微分相位投影残差
Figure FDA00036316478800000210
的第u个分量。
2.如权利要求1所述的基于相位CT的材料分解与识别方法，其特征在于，采用式(7)表示的第二优化模型求解
Figure FDA00036316478800000211
的最优值
Figure FDA00036316478800000212
Figure FDA00036316478800000213
式中，
Figure FDA00036316478800000214
表示
Figure FDA00036316478800000215
和当前相位投影残差
Figure FDA00036316478800000216
的二范数距离最小时的相位投影残差
Figure FDA00036316478800000217
3.如权利要求1或2所述的基于相位CT的材料分解与识别方法，其特征在于，离散化吸收投影数据
Figure FDA00036316478800000218
和离散化微分相位投影数据
Figure FDA00036316478800000219
表示式(3)：
Figure FDA00036316478800000220
其中：
Figure FDA00036316478800000221
表示投影角度
Figure FDA00036316478800000222
下不同探测器接收到的离散化吸收投影数据，其中的元素
Figure FDA00036316478800000223
表示第u个探测器像素采集的吸收投影值，U表示探测器的总像素个数，τ表示向量的转置；
Figure FDA00036316478800000224
表示投影角度
Figure FDA00036316478800000225
下的投影矩阵，其中的J表示每幅图像的总像素个数；
f＝(f1,f2...,fJ)τ表示第一基材的分布函数图像f(x,y)的离散化图像，其中的元素fj为f(x,y)在被测样品上第j个像素上的采样值；
g＝(g1,g2...,gJ)τ表示第二基材的分布函数图像g(x,y)的离散化图像，其中的元素gj为g(x,y)在被测样品上第j个像素上的采样值；
Figure FDA0003631647880000031
表示投影角度
Figure FDA0003631647880000032
下不同探测器的离散化微分相位投影数据，其中的元素
Figure FDA0003631647880000033
表示第u个探测器像素采集的微分相位投影值。
4.如权利要求3所述的基于相位CT的材料分解与识别方法，其特征在于，通过求解如下式(4)提供的第一优化模型，得到最后的分解图像：
Figure FDA0003631647880000034
其中：
f*表示第一优化模型取得最小时第一基材的分布函数图像f；
g*表示第一优化模型取得最小时第二基材的分布函数图像g；
R(·)是正则化算子。
5.如权利要求3所述的基于相位CT的材料分解与识别方法，其特征在于，还包括：
步骤6，根据基材图像f(m)和g(m)，利用公式(1)得到吸收图像和相位图像；
Figure FDA0003631647880000035
其中：
μ(x,y)表示被测样品的点(x,y)处的吸收图像；
δ(x,y)表示被测样品的点(x,y)处的相位图像；
μ1为第一基材的线性吸收图像；
μ2为第二基材的线性吸收图像；
δ1为第一基材的相位图像；
δ2为第二基材的相位图像；
f(x,y)为第一基材的分布函数图像；
g(x,y)为第二基材的分布函数图像；
步骤7，利用式(9)求解等效原子序数：
Figure FDA0003631647880000041
其中：
Cp、CE、CZ是通过扫描与被测物体相近的模体的数据进行标定求解的参数常量；
Figure FDA0003631647880000042
r0是经典原子半径，h是普朗克常数，c是光速常数；
CKN(E)为Klein-Nishina截面；
μ(E,x,y)表示被测样品的点(x,y)处对能量E的射线的吸收图像；
δ(E,x,y)表示被测样品的点(x,y)处对能量E的射线的相位图像；
ρe(x,y)表示被测样品的点(x,y)处的电子密度；
Z(x,y)表示被测样品的点(x,y)处的等效原子序数。
CN202110019682.0A 2021-01-07 2021-01-07 一种基于相位ct的材料分解与识别方法 Active CN112884853B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110019682.0A CN112884853B (zh) 2021-01-07 2021-01-07 一种基于相位ct的材料分解与识别方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110019682.0A CN112884853B (zh) 2021-01-07 2021-01-07 一种基于相位ct的材料分解与识别方法

Publications (2)

Publication Number Publication Date
CN112884853A CN112884853A (zh) 2021-06-01
CN112884853B true CN112884853B (zh) 2022-07-08

Family

ID=76047058

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110019682.0A Active CN112884853B (zh) 2021-01-07 2021-01-07 一种基于相位ct的材料分解与识别方法

Country Status (1)

Country Link
CN (1) CN112884853B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117007621B (zh) * 2023-03-13 2024-03-29 北京光影智测科技有限公司 基于微焦点光源的双能同轴相位ct材料分解方法及装置

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112179926A (zh) * 2020-09-24 2021-01-05 首都师范大学 一种基于同轴ct的相位-吸收反演和材料定量成像方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10111638B2 (en) * 2016-05-24 2018-10-30 Toshiba Medical Systems Corporation Apparatus and method for registration and reprojection-based material decomposition for spectrally resolved computed tomography

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112179926A (zh) * 2020-09-24 2021-01-05 首都师范大学 一种基于同轴ct的相位-吸收反演和材料定量成像方法

Also Published As

Publication number Publication date
CN112884853A (zh) 2021-06-01

Similar Documents

Publication Publication Date Title
US10274439B2 (en) System and method for spectral x-ray imaging
Touch et al. A neural network-based method for spectral distortion correction in photon counting x-ray CT
Longo et al. Towards breast tomography with synchrotron radiation at Elettra: first images
Cai et al. A full‐spectral Bayesian reconstruction approach based on the material decomposition model applied in dual‐energy computed tomography
US9808216B2 (en) Material decomposition of multi-spectral x-ray projections using neural networks
US11116470B2 (en) Beam hardening correction in x-ray dark-field imaging
Mohammadi et al. Quantitative evaluation of a single-distance phase-retrieval method applied on in-line phase-contrast images of a mouse lung
CN104939848B (zh) 单色图像的生成
Mechlem et al. Spectral angiography material decomposition using an empirical forward model and a dictionary-based regularization
CN112884853B (zh) 一种基于相位ct的材料分解与识别方法
Holbrook et al. Deep learning based spectral distortion correction and decomposition for photon counting CT using calibration provided by an energy integrated detector
US20190069865A1 (en) Apparatus for multi material decomposition
Chen et al. Convolutional neural network based attenuation correction for 123I-FP-CIT SPECT with focused striatum imaging
He et al. Energy-discriminative performance of a spectral micro-CT system
He et al. Preliminary experimental results from a MARS Micro-CT system
Pivot et al. Scatter correction for spectral CT using a primary modulator mask
Xie et al. Scatter correction for cone-beam computed tomography using self-adaptive scatter kernel superposition
CN116183647A (zh) 一种物质识别方法
Chukalina et al. A hardware and software system for tomographic research: reconstruction via regularization
Choi et al. A unified statistical framework for material decomposition using multienergy photon counting x‐ray detectors
Martz et al. CT dual-energy decomposition into x-ray signatures Rho-e and Ze
Deng et al. A method for material decomposition and quantification with grating based phase CT
Dong et al. An improved physics model for multi-material identification in photon counting CT
Chang et al. Spectrum estimation based on a parametric physical model for CT
Pham et al. A BVMF-B algorithm for nonconvex nonlinear regularized decomposition of spectral x-ray projection images

Legal Events

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