CN116609819A - 一种带散射校正的x射线能谱估计方法 - Google Patents

一种带散射校正的x射线能谱估计方法 Download PDF

Info

Publication number
CN116609819A
CN116609819A CN202310578280.3A CN202310578280A CN116609819A CN 116609819 A CN116609819 A CN 116609819A CN 202310578280 A CN202310578280 A CN 202310578280A CN 116609819 A CN116609819 A CN 116609819A
Authority
CN
China
Prior art keywords
energy spectrum
ray
phi
scattering
model
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.)
Pending
Application number
CN202310578280.3A
Other languages
English (en)
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.)
Southwest University of Science and Technology
Hunan First Normal University
Original Assignee
Southwest University of Science and Technology
Hunan First 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 Southwest University of Science and Technology, Hunan First Normal University filed Critical Southwest University of Science and Technology
Priority to CN202310578280.3A priority Critical patent/CN116609819A/zh
Publication of CN116609819A publication Critical patent/CN116609819A/zh
Pending legal-status Critical Current

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01TMEASUREMENT OF NUCLEAR OR X-RADIATION
    • G01T1/00Measuring X-radiation, gamma radiation, corpuscular radiation, or cosmic radiation
    • G01T1/36Measuring spectral distribution of X-rays or of nuclear radiation spectrometry
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01TMEASUREMENT OF NUCLEAR OR X-RADIATION
    • G01T7/00Details of radiation-measuring instruments
    • G01T7/005Details of radiation-measuring instruments calibration techniques

Landscapes

  • Physics & Mathematics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • High Energy & Nuclear Physics (AREA)
  • Molecular Biology (AREA)
  • Analysing Materials By The Use Of Radiation (AREA)

Abstract

本发明公开了一种带散射校正的X射线能谱估计方法,包括:对样品进行扫描得到原始CT投影数据;构造样品模型,对样品模型进行散射信号模拟,估计出探测器外壳带来的散射信号;对原始CT投影数据进行散射校正,得到散射校正后的投影数据;计算样品模型各个角度下的模型厚度;根据散射校正后的投影数据和各个角度下的模型厚度进行能谱估计。采用本发明提供的带散射校正的X射线能谱估计方案,更能准确快速估计出X射线的探测能谱。

Description

一种带散射校正的X射线能谱估计方法
技术领域
本发明涉及X射线技术领域,特别是关于一种带散射校正的X射线能谱估计方法。
背景技术
在大多数临床CT扫描仪中,X射线源往往是多色的。探测器探测得到的投影数据,除去散射信号影响,还受X射线源能谱分布和探测器响应的影响。
现有技术中常使用如下模型来刻画X射线穿过被测物体的衰减:
其中,I0(l)和I(l)为沿着路径l穿过被测物体前和后的X射线信号,ΩE是能量取值空间,φd(E):=φ(E)R(E)为探测器探测到的能谱,E表示能量值,φ(E)是X射线能谱,R(E)是探测器响应函数,μtot(x,E)线性衰减系数,Scl是与路径l对应的探测器像素接收到的散射信号强度。
由公式(1)可知,准确重建出CT图像(即准确估计出被测物质的线性衰减系数),除了考虑散射信号的影响外,还需要事先估计出准确的X射线能谱。
预先估计出准确的X射线的能谱分布,并应用该能谱进行投影数据校正,能够减少CT中的光束硬化伪影。关于X射线的能谱估计的方法,可以分为下面四类。一是直接使用光子计数探测器测量X射线的能谱,但探测器的能量响应无法准确知道,会导致估计不准。二是基于模型的方法,该方法使用经验或半经验物理模型生成能谱。三是将预估计的能谱表示预先生成能谱的不同权重的线性组合。四是从测量到的投影数据估计CT中X射线的能谱。第四类方法通常分为两个步骤:(i)测量X射线穿过不同厚度的已知材料(如均匀铝)的传输数据;(ii)求解表示X射线衰减过程的线性方程来重建能谱。在2005年,Sidky等提出了使用期望最大化(expectation-maximization,EM)算法来求解关于能谱估计的线性积分方程。
当遇到“数据缺失”或“不完整数据”时,求解使似然函数值极大化的参数值是很复杂的。Rubin等在1997年提出了适用于观测值为不完整数据时,迭代计算最大似然估计的一般方法--EM算法。EM算法的概要如下:
(i)设和/>是两个样本空间,且/>到/>存在映射。记/>为观测到的不完整数据,/>为完整数据,其概率密度为f(x|θ),θ∈Θ为待估的参数,Θ为参数取值空间。记Y=T(X),Y的概率密度
其中EM算法的目的是求出使ln(g(y|θ))取最大值的参数θ。
(ii)求似然函数的期望(E步),即
i=0,1,2,...为迭代次数。
(iii)求对数似然函数期望的最大值点/>(M步),即
(iv)取定初始估计后,重复E步和M步,得到一个EM估计序列/>且满足
Sidky等在2005年提出用EM从传输数据中估计能谱,并用下述式(6)刻画CT成像系统里的传输数据
其中ΩE是能谱φ(E)的取值空间,φd(E)=φ(E)R(E)是探测器能谱,R(E)是探测器响应函数。在一定离散条件下,方程(6)可离散为如下形式
其中I是抽样的能谱数,J是X-射线路径的索引集,{wi}是φd(E)在选定基组/>的扩展系数
根据Sidky等的研究,用EM算法迭代求解式(7),可得
其中k是迭代次数,初始值是在给定初始能谱后通过式(9)得到。
实验结果表明对无散射的投影数据,EM算法是一种准确且稳健的能谱估计方法。但是,散射信号会影X射线能谱估计的准确性,,因此亟待提供一种考虑散射的X射线能谱估计方案。
发明内容
本发明的目的在于提供一种带散射校正的X射线能谱估计方法,来克服或至少减轻现有技术的上述缺陷中的至少一个。
为实现上述目的,本发明提供一种带散射校正的X射线能谱估计方法,包括:
步骤1,对样品进行扫描得到原始CT投影数据;
步骤2,构造样品模型,对样品模型进行散射信号模拟,估计出探测器外壳带来的散射信号;
步骤3,对原始CT投影数据进行散射校正,得到散射校正后的投影数据
其中,
表示校正探测器外壳和被测物体带来的散射信号后的传输数据;I*(lj)和/>分别表示校正探测器外壳带来的散射信号后的穿过样品前和后的X射线信号的强度;Scj是与路径lj对应的探测器像素接收到的散射信号强度;ΩE是预设能量取值空间;φd(E):=φ(E)R(E)为探测器探测到的能谱,φ(E)是X射线能谱,R(E)是探测器响应函数,μtot(x,E)线性衰减系数;lj表示X射线路径,j表示第j条X射线路径,J表示X射线路径的指标集;
步骤4,计算样品模型各个角度下的模型厚度;
步骤5,根据散射校正后的投影数据和各个角度下的模型厚度进行能谱估计;包括:
对式(11)做离散化处理,得到式(12)和式(13):
其中,I是抽样的能谱数;i表示第i个能谱;为预设能谱基;wi为能谱φd(E)的展开系数,如下式:
利用下式(15)迭代求解wi
式(15)中,k为迭代次数,给定初始能谱后,初始值/>可通过式(14)确定,迭代求解得出wi
根据wi和式(14)得到φd(E),再根据φd(E):=φ(E)R(E)得到X射线能谱φ(E)。
优选的,该方法还包括:
步骤6,判断估计出的能谱是否满足预设条件;若满足,则以估计出的能谱作为X射线的探测能谱;否则,返回步骤2。
优选的,预设条件为:估计出的X射线探测能谱φd(E)与参考能谱的均方根误差小于阈值。
优选的,步骤4包括:
记样品模型为M,令x表示样品模型中的点,计算出δ(x)在每条X射线路径lj的积分值,作为X射线路径lj对应角度的模型厚度。
本发明由于采取以上技术方案,其具有以下优点:
采用本发明提供的带散射校正的X射线能谱估计方案,同时考虑了探测器带来的散射信号和样品本身带来的散射信号对投影数据的影响,使用带散射校正的EM算法,更能准确快速估计出X射线的探测能谱。
附图说明
图1为本发明提供的一种带散射校正的X射线能谱估计方法的流程示意图。
图2(a)为采用本发明提供的方法的示例中的用模拟数据生成的均匀铝块的几何示意图。
图2(b)为采用本发明提供的方法的示例中的给定参考能谱和初始能谱。
图3为采用本发明提供的方法的示例中给出的模拟数据的能谱估计的实验结果。
图4为采用本发明提供的方法的示例中给出的模拟数据估计的探测能谱的多色投影值与均匀铝块交线长的对应关系图。
图5为采用本发明提供的方法的示例中给出的实际数据的均匀Al块示意图。
图6(a)为采用本发明提供的方法的示例中估计的探测能谱示意图。
图6(b)-(d)为采用本发明提供的方法的示例中探测能谱的多色投影曲线示意图。
具体实施方式
在附图中,使用相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面结合附图对本发明的实施例进行详细说明。
在本发明的描述中,术语“中心”、“纵向”、“横向”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系,仅是为了便于描述本发明和简化描述,而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作,因此不能理解为对本发明保护范围的限制。
在不冲突的情况下,本发明各实施例及各实施方式中的技术特征可以相互组合,并不局限于该技术特征所在的实施例或实施方式中。
下面结合附图以及具体实施例对本发明做进一步的说明,需要指出的是,下面仅以一种最优化的技术方案对本发明的技术方案以及设计原理进行详细阐述,但本发明的保护范围并不仅限于此。
本文涉及下列术语,为便于理解,对其含义说明如下。本领域技术人员应当理解,下列术语也可能有其它名称,但在不脱离其含义的情形下,其它任何名称都应当被认为与本文所列术语一致。
本发明提供一种带散射校正的X-射线能谱估计方法(EM-SC),适用于含散射信号污染的投影数据。如图1所示,该方法包括:
步骤1,对样品进行扫描得到原始CT投影数据。
样品(或称模体)为铝圆柱或铝长方体等均匀材质物体,在常用条件,扫描样品的原始CT投影数据给定初始能谱信息及几何参数。
步骤2,构造样品模型,对样品模型进行散射信号模拟,估计出探测器外壳带来的散射信号。
在得到CT投影数据后可以进行CT重建。可以利用FDK算法进行CT重建,也可以利用其它算法,本文对此不做限定。可以通过FDK重建算法对CT投影数据重建出来的结果进行阈值分割得到样品模型。
对样品模型进行散射信号模拟可以估计出探测器外壳带来的散射信号。例如可以利用gQMCFFD方法对样品模型进行散射信号模拟,并用卷积方法估计探测器外壳带来的散射信号。利用gQMCFFD对模型进行散射信号模拟的具体方案在专利号为202111107623.5的发明专利中有详细介绍,本文不再赘述。用卷积方法估计探测器外壳带来的散射信号也可以利用现有技术实现。例如,J.Star-Lack,M.Sun,A.Kaestner等人在Physics of MedicalImaging上发表的文章Efficient scatter correction using asymmetric kernels[J],具体介绍了用卷积方法估计探测器外壳带来的散射信号的方法,本文不再赘述。
步骤3,对CT投影数据进行散射校正,得到散射校正后的投影数据
其中,
表示校正探测器外壳和被测物体带来的散射信号后的传输数据;I*(lj)和/>分别表示校正探测器外壳带来的散射信号后的X射线信号;Scj是与路径lj对应的探测器像素接收到的散射信号强度;ΩE是预设能量取值空间;φd(E):=φ(E)R(E)为探测器探测到的能谱,φ(E)是X射线能谱,R(E)是探测器响应函数,μtot(x,E)线性衰减系数;lj表示X射线路径,j表示第j条X射线路径,J表示X-射线路径的指标集。R(E)为预设函数,与具体探测器有关,取值范围为[0,1]。
步骤4,计算样品模型各个角度下的模型厚度。
该步骤包括:记样品模型为M,令x表示样品模型中的点,计算出δ(x)在每条X射线路径lj的积分值,作为每条X射线路径lj对应角度下的模型厚度。
步骤5,根据散射校正后的投影数据和各个角度下的模型厚度进行能谱估计;包括:
对式(11)做离散化处理,得到式(12)和式(13):
其中,I是抽样的能谱数;i表示第i个能谱;为预设能谱基;wi为能谱φd(E)的展开系数,如下式:
利用下式(15)迭代求解wi
式(15)中,k为迭代次数,给定初始能谱后,初始值/>可通过式(14)确定,迭代求解得出/>
根据和式(14)得到第k次迭代后的φd(E),再根据φd(E):=φ(E)R(E)得到X射线能谱φ(E)。
该方法还可以包括:
步骤6,判断估计出的能谱是否满足预设条件;若满足,则以估计出的能谱作为最终能谱;否则,返回步骤2。
其中,预设条件为:估计出的X射线探测能谱φd(E)与参考能谱的均方根误差小于阈值。
本发明还提供一种带散射校正的X射线能谱估计装置,用于执行上述任意实施例或示例中的方法。该装置可以作为一个部件安装在CT扫描仪中,也可以作为独立装置与CT扫描仪连接。
在一种示例中,该装置包括输入接口、处理器、和输出接口,处理器用于执行上述任意实施例或示例中的方法中的各步骤,输入接口用于接收处理器进行X射线能谱估计所需要的外部信息,输出接口用于输出处理器估计出的带散射校正的X射线能谱。
采用本发明提供的带散射校正的X射线能谱估计方案,更能准确快速估计出X射线的探测能谱。
实验结果
本小节给出用带散射校正的EM算法估计出能谱的数值实验,来说明带散射校正的EM算法(EM-SC)的准确性和稳健性。这里先分别给出模拟数据和真实数据的实验结果。对模拟数据,我们先采集X射线沿不同路径穿过体积为80×80×80mm3的均匀铝块的衰减数据。具体操作如下:平移铝块(45度对着X-射线源和探测器方向)进行平移扫描。均匀铝块的几何示意图如图2(a)所示。采集到的衰减数据是带有散射信号污染和噪声污染的。接着用gQMCFFD算法估计X-射线沿不同路径穿过均匀铝块的散射信号,校正投影数据中散射信号。图2(b)表示给定参考能谱和初始能谱,参考能谱也称假定的真实能谱。
图3给出了模拟数据的能谱估计的实验结果。图3(a)是EM算法迭代128次估计的能谱、带散射校正的EM算法(EM-SC)迭代128次估计的能谱与参考能谱的比较。从图3(a)我们可以看出带散射校正后的EM算法估计出的能谱与真实能谱高度吻合,且带散射校正后的EM算法迭代128次估计的能谱与真实能谱的均方根误差(RMSE)为0.001337,EM方法迭代128次估计的能谱与真实能谱的RMSE为0.00535,这说明了带散射校正的EM算法的准确性。图3(b)是EM算法估计的能谱和带散射校正的EM算法估计的能谱与参考能谱比较的RMSE以及其随迭代次数变化的走势。从图3(b)可知EM方法估计的能谱与真实能谱的RMSE随着迭代次数的增加而增长,这可能是迭代过程中引入了噪声所引起的。这说明了对带有散射信号污染的数据,EM算法稳健性较弱。带散射校正的EM算法与参考能谱的RMSE先随着迭代次数的增加而减少,在迭代128次时达到最小,并随着迭代次数的增加保持稳定,这说明了带散射校正的EM算法的稳健性。
为了进一步说明带散射校正的EM算法的准确性,我们给出了模拟数据估计的探测能谱的多色投影值与均匀铝块交线长的对应关系图,如图4所示。图4(a)和(b)分别是未散射校正和散射校正的多色投影值与均匀铝块交线长的对应关系,其中点组成的关系图由测量数据得到,虚线由探测到的探测能谱反解得到。从图4(a)中,我们可以发现散射信号会导致多色投影值的非线性减少,而且随着路线交线长的增加越明显。这说明先校正散射信号的必要性,也解释了为什么带散射校正的EM方法比EM方法出理带散射的数据更准确。从图4(b)中,我们发现用带散射校正的EM算法估计出的能谱计算出的多色投影值与均匀铝块交线长的对应关系曲线和测量到多色投影值与均匀铝块交线长的对应关系曲线,高度拟合,这进一步说明了带散射校正的EM算法的准确性。
为估计实际数据的探测能谱,本文使用40×40×80mm3的均匀Al块,其示意图如图5所示。用带散射校正的EM算法(128迭代)和EM(128迭代)估计的探测能谱如图6(a)所示。相对应的探测能谱的多色投影曲线如图6(b)-(d)所示。带散射校正的EM算法(EM-SC)估计的结果比EM算法更准确,但同时考虑了被测物体本身和探测器外壳带来的散射的EM算法(EM-SC)比只考虑了被测物体自身带来散射的EM算法(EM-SC0)的估计结果要好一些。
最后需要指出的是:以上实施例仅用以说明本发明的技术方案,而非对其限制。本领域的普通技术人员应当理解:可以对前述各实施例所记载的技术方案进行修改,或者对其中部分技术特征进行等同替换;这些修改或者替换,并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (4)

1.一种带散射校正的X射线能谱估计方法,其特征在于,包括:
步骤1,对样品进行扫描得到原始CT投影数据;
步骤2,构造样品模型,对样品模型进行散射信号模拟,估计出探测器外壳带来的散射信号;
步骤3,对原始CT投影数据进行散射校正,得到散射校正后的投影数据
其中,
表示校正探测器外壳和被测物体带来的散射信号后的传输数据;I*(lj)和/>分别表示校正探测器外壳带来的散射信号后的穿过样品前和后的X射线信号的强度;Scj是与路径lj对应的探测器像素接收到的散射信号强度;ΩE是预设能量取值空间;φd(E):=φ(E)R(E)为探测器探测到的能谱,φ(E)是X射线能谱,R(E)是探测器响应函数,μtot(x,E)线性衰减系数;lj表示X射线路径,j表示第j条X射线路径,J表示X-射线路径的指标集;
步骤4,计算样品模型各个角度下的模型厚度;
步骤5,根据散射校正后的投影数据和各个角度下的模型厚度进行能谱估计;包括:
对式(11)做离散化处理,得到式(12)和式(13):
其中,I是抽样的能谱数;i表示第i个能谱;为预设能谱基;wi为能谱φd(E)的展开系数,如下式:
利用下式(15)迭代求解wi
式(15)中,k为迭代次数,给定初始能谱后,初始值/>可通过式(14)确定,迭代求解得出wi
根据wi和式(14)得到φd(E),再根据φd(E):=φ(E)R(E)得到X射线能谱φ(E)。
2.根据权利要求1所述的带散射校正的X射线能谱估计方法,其特征在于,还包括:
步骤6,判断估计出的能谱是否满足预设条件;若满足,则以估计出的能谱作为X射线的探测能谱;否则,返回步骤2。
3.根据权利要求2所述的带散射校正的X射线能谱估计方法,其特征在于,预设条件为:估计出的X射线探测能谱φd(E)与参考能谱的均方根误差小于阈值。
4.根据权利要求1所述的带散射校正的X射线能谱估计方法,其特征在于,步骤4包括:
记样品模型为M,令x表示样品模型中的点,计算出δ(x)在每条X射线路径lj的积分值,作为X射线路径lj对应角度的模型厚度。
CN202310578280.3A 2023-05-22 2023-05-22 一种带散射校正的x射线能谱估计方法 Pending CN116609819A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310578280.3A CN116609819A (zh) 2023-05-22 2023-05-22 一种带散射校正的x射线能谱估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310578280.3A CN116609819A (zh) 2023-05-22 2023-05-22 一种带散射校正的x射线能谱估计方法

Publications (1)

Publication Number Publication Date
CN116609819A true CN116609819A (zh) 2023-08-18

Family

ID=87679540

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310578280.3A Pending CN116609819A (zh) 2023-05-22 2023-05-22 一种带散射校正的x射线能谱估计方法

Country Status (1)

Country Link
CN (1) CN116609819A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117270024A (zh) * 2023-11-20 2023-12-22 合肥综合性国家科学中心人工智能研究院(安徽省人工智能实验室) 能谱响应函数的校正方法、装置、计算机设备及存储介质

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117270024A (zh) * 2023-11-20 2023-12-22 合肥综合性国家科学中心人工智能研究院(安徽省人工智能实验室) 能谱响应函数的校正方法、装置、计算机设备及存储介质
CN117270024B (zh) * 2023-11-20 2024-02-20 合肥综合性国家科学中心人工智能研究院(安徽省人工智能实验室) 能谱响应函数的校正方法、装置、计算机设备及存储介质

Similar Documents

Publication Publication Date Title
US20220035961A1 (en) System and method for artifact reduction of computed tomography reconstruction leveraging artificial intelligence and a priori known model for the object of interest
US11633146B2 (en) Automated co-registration of prostate MRI data
US7203267B2 (en) System and method for boundary estimation using CT metrology
Villarraga-Gómez et al. Empirical approaches to uncertainty analysis of X-ray computed tomography measurements: A review with examples
JP2008502395A (ja) 投影放射線撮影、特に乳房撮影における散乱放射線補正装置および方法
JP7408664B2 (ja) X線撮像システム
US11116470B2 (en) Beam hardening correction in x-ray dark-field imaging
CN116609819A (zh) 一种带散射校正的x射线能谱估计方法
Khan et al. Particle-filter-based multisensor fusion for solving low-frequency electromagnetic NDE inverse problems
Mohan et al. SABER: a systems approach to blur estimation and reduction in x-ray imaging
Liu et al. Nonlinear ultrasonic transmissive tomography for low-contrast biphasic medium imaging using continuous-wave excitation
Niinimaki et al. Multiresolution parameter choice method for total variation regularized tomography
CN109916933B (zh) 基于卷积神经网络的x射线计算机断层成像能谱估计方法
US20220042932A1 (en) Estimating Background Radiation from Unknown Sources
US20080205731A1 (en) Noise Model Selection for Emission Tomography
Duncan MARS imaging for geology
US11302042B2 (en) Method for correcting nonlinearities of image data of at least one radiograph, and computed tomography device
Alsaffar et al. Multi-material blind beam hardening correction in near real-time based on non-linearity adjustment of projections
CN113391341B (zh) 一种考虑散射光子影响的x射线能谱估计方法
Wang et al. On the edge detection of an image by numerical differentiations for gray function
EP3920144A1 (en) Methods and systems for multi-material decomposition
Yin et al. Parametric boundary reconstruction algorithm for industrial CT metrology application
EP2922028A1 (en) Exploiting residual projection data and reference data for processing data produced by imaging modality
US20230154068A1 (en) Device and method for pet image analysis and reconstruction
Li et al. Uncertainty Quantification of Density Reconstruction Using MCMC Method in High-Energy X-ray Radiography

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