一种基于双能同轴相位CT的材料分解方法与流程

一种基于双能同轴相位ct的材料分解方法
技术领域
1.本发明涉及一种ct(英文全称为“computer tomography”，中文“电子计算机断层摄影”)的材料分解方法，特别是关于一种基于双能同轴相位ct的材料分解方法。


背景技术：

2.传统ct成像技术被广泛应用在很多领域，例如医学诊断、工业检测等等。由于该成像原理是基于不同物质之间衰减的差异性，因此对于弱吸收的样品(low-z样品，如脑，乳腺，肺等)很难进行可视化。幸运的是，对物体的相移特性敏感的相位衬度成像技术已经出现。low-z材料的相移差异比它们的吸收差数值上大三个数量级左右，这对研究预临床样品的高空间分辨率非常有利。相位衬度成像有干涉法、衍射增强法、光栅微分法和同轴法。其中同轴法具有成像实验设备简单，与吸收ct实验光路类似，易于实现的特点而受到关注。
3.材料分解是ct成像中一个重要的应用。传统的双能ct材料分解方法，比如基于图像的方法，也有基于投影的方法，还有迭代方法。但是，这些种成像技术对于弱吸收样品来说，它们的衰减是微乎其微的，因此很难对该类样品有效成像。


技术实现要素：

4.本发明的目的在于提供一种基于双能同轴相位ct的材料分解方法来克服或至少减轻现有技术的上述缺陷中的至少一个。
5.为实现上述目的，本发明提供一种基于双能同轴相位ct的材料分解方法，其包括：
6.步骤1，初始化第一基材ct图像fm和第二基材ct图像gm，预设迭代终止条件；
7.步骤2，通过下式描述的旋转角度和第i能量ei下的原始强度数据与第一、二基材ct图像f、g的交涉，计算和ei下的强度残差ii；
[0008][0009]
式中，i为能量的索引，为傅里叶变换操作，m为和ei下第m轮迭代的强度，μ
1，i
、δ
1，i
分别为第一基材在ei下的线性衰减系数、相移因子，μ
2，i
、δ
2，i
分别为第二基材在ei下的线性衰减系数、相移因子，z为样品到探测器的距离，τ2为频率域形式，为梯度算子，为下的投影算子；
[0010]
步骤3，令和通过下式计算f
err
和g
err
；
[0011][0012]
式中，i1、i2分别为和e1、e2下的强度残差，μ
1，1
、μ
1，2
分别为第一基材在e1、e2下的线性衰减系数，μ
2，1
、μ
2，2
分别为第二基材在e1、e2下的线性衰减系数，δ
1，1
、δ
1，2
分别为第一基材在e1、e2下的相移因子，δ
2，1
、δ
2，2
分别为第二基材在e1、e2下的相移因子；
[0013]
步骤4，更新f的第m 1轮迭代的值f
m 1
和g的第m 1轮迭代的值g
m 1
；
[0014]
步骤5，若未达到迭代终止条件，则令m＝m 1，并返回步骤2；
[0015]
步骤6，返回fm和gm。
[0016]
进一步地，步骤4采用下式迭代更新f
m 1
和g
m 1
；
[0017][0018]
式中，分别为第m轮迭代f、g的第j个像素值，分别为第m 1轮迭代f、g的第j个像素值，ε为松弛因子，表示f、g的第j个像素值在旋转角度下x射线沿着探测器的第u个探测器单元下的投影路径，j＝1,2,
…
j，u＝1，2，
…
u，u为探测器的探测器单元总数，j为像素的总数，
[0019]
进一步地，步骤2的的获取方法包括：
[0020]
使用下式描述的菲涅尔衍射，模拟贴合实际物理机理下的正向传播，得到
[0021][0022]
式中，h
z,i
为ei下的菲涅尔传播因子，为和ei下的透射振幅函数，为卷积算子，||为中括号。
[0023]
进一步地，h
z,i
的获取公式为：
[0024][0025]
式中，i表示虚数单位，k为ei下的波数，x为空间域坐标，λ为ei下的波长。
[0026]
进一步地，当ii为二维数组时，h
z,i
的获取公式为：
[0027]
[0028]
式中，i表示虚数单位，k为ei下的波数，(x，y)为空间域坐标，λ为ei下的波长。
[0029]
本发明由于采取以上技术方案，其具有以下优点：
[0030]
本发明提供的一步法，即在基于双能量同轴相位衬度ct(omd-ppct)中直接从原始数据进行材料分解，得到基材ct图像。
附图说明
[0031]
图1为本发明实施例的数值实验中使用的模体，其中：(a)、(b)、(c)对应使用的模体分别为数值体模、teflon材料、pm(英文全称为“polyme metha”)材料。
[0032]
图2为本发明实施例提供的在无噪声情况下的材料分解示意图。
[0033]
图3为本发明实施例提供的在有噪声情况下的材料分解示意图。
[0034]
图4为本发明实施例提供的使用的模体分解结果剖面图，其中：(a)和(a’)分别是无噪声和有噪声情况下的pm材料结果；(b)和(b’)分别是无噪声和有噪声情况下的teflon材料结果。
具体实施方式
[0035]
下面结合附图和实施例对本发明进行详细的描述。
[0036]
本发明实施例提供的基于双能同轴相位ct的材料分解方法包括：
[0037]
步骤1，初始化第一基材ct图像fm和第二基材ct图像gm，预设迭代终止条件。其中，第一基材ct图像fm、第二基材ct图像gm、m的初始化值均可为0，当然，也可以根据实际情况设置成其他数值。迭代终止条件例如可以是最大迭代次数。
[0038]
步骤2，在一个实施例中，强度传输方程(tie，英文全称为“transport of intensity equation”)还可以描述x射线波在z中的传播，如下式(1)所示：
[0039][0040]
式中，为梯度算子，为旋转角度下距离样品z处的强度数据，为z＝0时(x，y)所在的坐标系的原点通常设置在样品中心，间距则根据射线间隔确定。
[0041]
假设在中没有大的强度梯度，为了简化该数学模型，强度传输方程在一个实施例也可以描述成下式(2)：
[0042][0043]
为了进一步简化该数学模型，假设当时，强度传输方程在另一个实施例还可以描述成下式(3)：
[0044][0045]
当然，强度传输方程还可以式(3)的基础之上，继续变形，比如：在又一实施例中描
述成下式(4)：
[0046][0047]
式中，μ(x，y)为在空间域点(x，y)的衰减函数，δ(x，y)为在空间域点(x，y)的相移函数。
[0048]
对于两个或更多的样品基本材料，衰减和相移可以由一组基函数定义本实施例中，对于双材料分解，如下式(5)：
[0049][0050]
式中，f(x，y)、g(x，y)分别为第一、二基材离散化ct图像，μ1、δ1分别为第一基材的线性衰减系数、相移因子，μ2、δ2分别为第二基材的线性衰减系数、相移因子。
[0051]
在一个实施例中，加入基材信息，强度传输方程可以描述成下式(6)：
[0052][0053]
在本实施例中，i＝1，2，e1，e2为不同能量下的强度。对上式(6)进行快速傅里叶变换，得到下式(7)：
[0054][0055]
本实施例中，为傅里叶变换操作，τ2为频率域形式。
[0056]
令f＝(f1，f2，
…fj
)
τ
和g＝(g1，g2，
…gj
)
τ
代表f(x，y)和g(x，y)的离散图像，在本实施例中，fj和gj为f(x，y)和g(x，y)在第j像素的像素值。j为总像素数目，τ为转置算子。为在旋转角度为下的投影算子，代表fj和gj在旋转角度为x射线沿着u下的投影路径，u为探测器单元总个数。本发明实施例获得在m轮迭代下的强度表示为式(8)：
[0057][0058]
得到旋转角度和第i能量ei下的强度残差ii，如下式(9)所示：
[0059][0060]
式中，μ
1，i
、δ
1，i
分别为第一基材在第i能量ei下的线性衰减系数、相移因子，μ
2，i
、δ
2，i
分别为第二基材在第i能量ei下的线性衰减系数、相移因子，z为样品到探测器的距离，为傅里叶变换操作，τ2为频率域形式，为梯度算子，为在旋转角度下的投影算子。
[0061]
使用式(9)描述的原始数据与第一、二基材ct图像f、g的交涉，计算旋转角度和第i能量ei下的强度残差ii。
[0062]
步骤3，令和根据式(9)有下
式(10)：
[0063][0064]
再通过式(11)计算f
err
和g
err
；
[0065][0066]
步骤4，更新第一基材ct图像f在第m 1轮迭代中的值f
m 1
和第二基材ct图像g在第m 1轮迭代中的值g
m 1
；
[0067]
步骤5，若未达到迭代终止条件，则令m＝m 1，并返回步骤2；
[0068]
步骤6，返回fm和gm。
[0069]
本发明实施例通过一步法，直接从原始数据中获得基材ct图像，而在这个过程中不需要进行相位恢复。
[0070]
在一个实施例中，步骤4采用式(12)迭代更新f
m 1
和g
m 1
；
[0071][0072]
式中，表示f、g的第j个像素值在旋转角度下x射线沿着探测器的第u个探测器单元下的投影路径，为第m轮迭代在第一基材ct图像f的第j个像素值，为第m轮迭代在第二基材ct图像g的第j个像素值，ε为松弛因子，根据模拟实验确定最佳数值，一般其具体范围为0～1，u为探测器的探测器单元总数，j为像素的总数。
[0073]
本实施例通过式(12)的反馈和校正，可以有效地抑制噪声。
[0074]
在一个实施例中，步骤2的原始数据的获取方法包括：
[0075]
使用式(13)描述的菲涅尔衍射，模拟更贴合实际物理机理下的正向传播，得到第m轮迭代在旋转角度和能量i下的强度数据
[0076][0077]
式中，h
z,i
为在第i能量ei下的菲涅尔传播因子，为在旋转角度和第i能量ei下的透射振幅函数，为卷积算子，||为中括号。
[0078]
本实施例符合低能下的物理机理，由于菲涅耳传播函数是类高斯函数，式(13)采
用卷积运算，可以使得在演化过程中对波前起到扩展和平滑作用。
[0079]
在上述实施例中，当iz为一维数组时，h
z,i
的获取公式为式(14)：
[0080][0081]
式中，i表示虚数单位，x为空间域坐标，k为第i能量ei量下的波数，λ为第i能量ei量下的波长，z为样品到探测器的距离(odd)，那么在探测器处的强度为iz。
[0082]
在上述实施例中，当ii为二维数组时，h
z,i
的获取公式为式(15)：
[0083][0084]
式中，i表示虚数单位，(x，y)为空间域坐标，k为第i能量ei量下的波数，λ为第i能量ei量下的波长。
[0085]
作为比较，本发明选择了传统image-basedct方法和使用基于光谱传播的相位成像(md-spbi)的材料分解。
[0086]
本发明设计了两种低z基材的样品。如图1所示，该样品由三个同心均匀材料的球形层组成，不同密度的材料缠绕在芯部，即teflon和polyme metha(pm)。在图1中的(c)中，本发明假设pm(1)是pm的标准密度，pm(2)＝1.05*pm(1)，pm(3)＝1.1*pm(4)。样品的直径尺寸为1.7mm。在模拟中，平行光束设置用于获取180度等间距的360个投影。测试了106光子发射通量对应的无噪声和泊松噪声。对于材料分解，选择teflon和pm作为多材料对象中的基材，重建图像的大小为512*512。
[0087]
基材分解的结果如图2和图3所示，本发明还放大了矩形的内容。
[0088]
图4显示了图3中橙色线条的轮廓。从轮廓和放大区域，发现现有的image-based和spbi两种方法对于pm-2具有非常弱的空间分辨率，并且所提出的算法比其他方法具有更好的材料分解和噪声抑制。这些结果为本发明的omd-ppct中精确材料分解的能力提供了重要证据。
[0089]
最后需要指出的是：以上实施例仅用以说明本发明的技术方案，而非对其限制。本领域的普通技术人员应当理解：可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。