CN113804709B - 基于拟蒙特卡罗和强制探测的校正ct散射信号的方法 - Google Patents

基于拟蒙特卡罗和强制探测的校正ct散射信号的方法 Download PDF

Info

Publication number
CN113804709B
CN113804709B CN202111107623.5A CN202111107623A CN113804709B CN 113804709 B CN113804709 B CN 113804709B CN 202111107623 A CN202111107623 A CN 202111107623A CN 113804709 B CN113804709 B CN 113804709B
Authority
CN
China
Prior art keywords
photon
scattering
point
detector
probability
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
CN202111107623.5A
Other languages
English (en)
Other versions
CN113804709A (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
Southwest University of Science and Technology
Original Assignee
Tsinghua University
Southwest University of Science and Technology
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, Southwest University of Science and Technology filed Critical Tsinghua University
Priority to CN202111107623.5A priority Critical patent/CN113804709B/zh
Publication of CN113804709A publication Critical patent/CN113804709A/zh
Application granted granted Critical
Publication of CN113804709B publication Critical patent/CN113804709B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N23/00—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00
    • G01N23/02—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material
    • G01N23/04—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and forming images of the material
    • G01N23/046—Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and forming images of the material using tomography, e.g. computed tomography [CT]
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2223/00—Investigating materials by wave or particle radiation
    • G01N2223/03—Investigating materials by wave or particle radiation by transmission
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2223/00—Investigating materials by wave or particle radiation
    • G01N2223/10—Different kinds of radiation or particles
    • G01N2223/101—Different kinds of radiation or particles electromagnetic radiation
    • G01N2223/1016—X-ray
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2223/00—Investigating materials by wave or particle radiation
    • G01N2223/30—Accessories, mechanical or electrical features
    • G01N2223/303—Accessories, mechanical or electrical features calibrating, standardising

Landscapes

  • Health & Medical Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Pulmonology (AREA)
  • Radiology & Medical Imaging (AREA)
  • Theoretical Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Biochemistry (AREA)
  • General Health & Medical Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Immunology (AREA)
  • Pathology (AREA)
  • Analysing Materials By The Use Of Radiation (AREA)

Abstract

本发明实施例公开了一种基于拟蒙特卡罗和强制探测的校正CT散射信号的方法，包括：模拟光子从X射线源S通过被测物体M到达探测器D的散射信号；将每个探测器像素检测到的信号减去模拟得到的散射信号与预设系数的乘积，根据差值确定穿过被测模体M后的X射线信号。通过计算光子理论上发生散射的概率，把QMC方法的优势和强制固定探测(FFD)技术结合起来，能够准确地得到探测器探测到的光子的散射信号，从而实现对探测器探测到的X射线的散射信号的校正。

Description

基于拟蒙特卡罗和强制探测的校正CT散射信号的方法
技术领域
本发明涉及计算机断层成像技术(CT)领域，特别是关于一种基于拟蒙特卡罗和强制探测的校正CT散射信号的方法。
背景技术
计算机断层成像技术(Computed Tomography，CT)以其无损检测、图像不重叠和高分辨率等优势，广泛应用于材料科学、无损检测、医学等诸多领域，尤其在医学诊断，彻底改变了放射诊断学。
CT图像重建是指通过测量X射线穿过被测模体M后的衰减数据，重建出被测模体M的线性衰减系数μtot(x,E)的过程。设I0(l)和I(l)为沿着路径l穿过模体M前和后的X射线信号。在探测X射线衰减数据时，探测器D会同时探测到主信号和散射信号，这会导致测量到的I(l)的增加。Lambert-Beer定律描述了在理想条件下，不考虑光子散射，X射线穿过被测模体M的衰减。这样，由Lambert-Beer定律
I(l)＝I0(l)exp[-∫lμtot(x,E)dx]， (0)
得到的μtot(x,E)会被低估。
散射信号会在测量得到的μtot值中引起非线性误差，并对重建的CT图像产生复杂的影响，可能导致高衰减图像区域之间的条状伪影、均质物体中的杯状伪影以及重建CT图像的对比度分辨率降低等。此外，当在CT中使用高能X射线或多排、平板探测器时，散射信号占总信号(主信号和散射信号的总和)的比重会增加。
因此，为了能够提高线性衰减系数μtot(x,E)的值的准确性，需要对探测器探测到的X射线的散射信号进行校正。
发明内容
本发明的目的在于提供一种基于拟蒙特卡罗和强制探测的校正CT散射信号的方法，来克服或至少减轻现有技术的上述缺陷中的至少一个。
为实现上述目的，本发明提供一种基于拟蒙特卡罗和强制探测的校正CT散射信号的方法，包括：
步骤10，模拟光子从X射线源S通过被测物体M到达探测器D的散射信号；包括：
步骤11，通过沃克混叠方法，用预先生成的4n维Sobol'点u＝(u1,u2,...,u4n)的第一个分量u1，从预先给定的能谱分布φ(E0)中抽样光子的初始能量E0，其中n表示光子在被测物体内发生碰撞的次数；
步骤12，记初始入射光到探测器D的立体角为Ω0，用4n-维Sobol'点u的第二个和第三个分量u2和u3在Ω0中均匀抽样单位入射方向
Figure BDA0003273083400000021
得到初始流强I0(u2,u3)以及初始入射X射线沿单位入射方向
Figure BDA0003273083400000022
运动与被测物体M的交点A0，得到光子以权重W0＝I0(u2,u3)沿
Figure BDA0003273083400000023
方向从A0开始在被测物体M中运动；
步骤13，在确定第i-1阶相互作用点Ai-1,i＝2,…,n的位置后，抽样光子在Ai-1的散射类型和散射角方向；包括：使用所述4n维Sobol’点u的第4(i-1)+1个分量u4(i-1)+1，通过重要性抽样方法，抽样散射类型，并使用u的第4(i-1)+2和4(i-1)+3个分量u4(i-1)+2和u4(i-1)+3，通过RITA算法抽样散射单位方向
Figure BDA0003273083400000024
光子从当前位置Ai-1,i＝1,…,n开始，沿着单位散射方向
Figure BDA0003273083400000025
逃出被测物体M的概率为：
Figure BDA0003273083400000026
其中
Figure BDA0003273083400000027
是光子初始单位入射方向，
Figure BDA0003273083400000028
是光子在Ai-1发生碰撞后的散射方向，δi-1＝0,1，
Figure BDA0003273083400000029
和
Figure BDA00032730834000000210
分别表示光子在Ai-1发生康普顿和瑞利散射后的剩余能量，ci-1是光子从Ai-1运动到Ci-1的长度，Ci-1是光子沿
Figure BDA0003273083400000031
方向运动到与被测物体M边界上的交点，
Figure BDA0003273083400000032
μtot表示线性衰减系数；主信号，即光子从光源S出发，沿单位入射方向
Figure BDA0003273083400000033
逃出被测物体M并到达探测器D的概率，是
Figure BDA0003273083400000034
光子以权重
Figure BDA0003273083400000035
在被测物体M中进行第i次相互作用；
步骤14，基于下式确定下一个相互作用点Ai(i＝1,2,...,n)：
Figure BDA0003273083400000036
其中si,i＝1,…,n表示位置Ai-1到下一个相互作用点Ai的路径长度，
Figure BDA0003273083400000037
u4(i-1)+4是4n维Sobol'点的第(4(i-1)+4)个分量；
步骤15，确定光子在Ai,i＝1,2,...,n出发生第i次相互作用后，沿单位强制方向
Figure BDA0003273083400000038
到达固定探测器像素Dj,j＝1,...,m×m的概率：
Figure BDA0003273083400000039
其中，
Figure BDA00032730834000000310
和
Figure BDA00032730834000000311
是光子在相互作用点Ai发生康普顿散射和瑞利散射的概率，
Figure BDA00032730834000000312
是探测器像素Dj与Ai相对应的立体角，
Figure BDA00032730834000000313
Figure BDA00032730834000000314
是探测器像素Dj的面积，αi表示
Figure BDA00032730834000000315
与探测器像素Dj法线方向的夹角,
Figure BDA00032730834000000316
是光子沿单位强制方向
Figure BDA00032730834000000317
运动与被测物体M交点，
Figure BDA00032730834000000318
和
Figure BDA00032730834000000319
是康普顿散射和瑞利散射极角的概率密度，且是归一化的线性衰减系数
Figure BDA00032730834000000320
Figure BDA00032730834000000321
步骤16，获得光子沿路径li,j:S→A0→A1→…→Ai→Bi,j→Dj,i∈{1,...,n}的概率为：
P(li,j)＝WiPi(Ai→Dj)，
光子沿路径S→A0→A1的概率与预先生成的4n维Sobol’点u的分量u1、u2、u3和u4相关，故：
P(l1,j):＝P1(A1→Dj)g0(u1,...,u4)，
其中，光子运动路径S→A0→A1→…→Ai是马尔科夫链，光子运动路径S→A0→A1由u1、u2、u3和u4抽样所得，光子由Ai-1运动到下一阶相互作用点Ai,i＝2,…,n由u4(i-1)+1、u4(i-1)+2、u4(i-1)+3和u4(i-1)+4抽样所得，且相互作用点Ai只与上一阶相互作用点Ai-1有关，因此得到：
Figure BDA0003273083400000041
步骤17，获得光子在被测物体M内发生n次碰撞后，固定探测器像素Dj,j＝1,...,m×m上累计接收到的散射总信号为
Figure BDA0003273083400000042
所述方法还包括步骤20，将每个探测器像素检测到的信号减去模拟得到的散射信号与预设系数的乘积，根据差值确定穿过被测模体M后的X射线信号。
在一种实施方式中，所述线性衰减系数μtot(x,E)表示光子穿过单位路径与被测物体发生相互作用的概率，是关于空间坐标x和能量E的函数，且
Figure BDA0003273083400000043
其中
Figure BDA0003273083400000044
是康普顿散射系数，
Figure BDA0003273083400000045
是瑞利散射系数，μa(x,E)是光电吸收系数。
本发明由于采取以上技术方案，其具有以下优点：
通过计算光子理论上发生散射的概率，把QMC方法的优势和强制固定探测(FFD)技术结合起来，能够准确地得到探测器探测到的光子的散射信号，从而实现对探测器探测到的X射线的散射信号的校正。
附图说明
图1为本发明实施例提供的基于拟蒙特卡罗和强制探测的校正CT散射信号的方法的流程示意图。
图2为本发明实施例提供的光子从X射线源S射出经过被侧物体M到达探测器D的过程示意图。
图3示出本发明实施例提供的光子的能谱分布示意图。
图4示出均匀铝模体和人头模体的几何示意图。
图5示出对均匀铝模体，gMCFRD、gQMCFRD、gMCFFD和gQMCFFD四种方法计算的散射信号。
图6示出对人头模体，gMCFRD、gQMCFRD、gMCFFD和gQMCFFD四种方法计算的散射信号。
具体实施方式
在附图中，使用相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面结合附图对本发明的实施例进行详细说明。
在本发明的描述中，术语“中心”、“纵向”、“横向”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明保护范围的限制。
在不冲突的情况下，本发明各实施例及各实施方式中的技术特征可以相互组合，并不局限于该技术特征所在的实施例或实施方式中。
下面结合附图以及具体实施例对本发明做进一步的说明，需要指出的是，下面仅以一种最优化的技术方案对本发明的技术方案以及设计原理进行详细阐述，但本发明的保护范围并不仅限于此。
本文涉及下列术语，为便于理解，对其含义说明如下。本领域技术人员应当理解，下列术语也可能有其它名称，但在不脱离其含义的情形下，其它任何名称都应当被认为与本文所列术语一致。
本发明实施例提供一种基于拟蒙特卡罗和强制探测的校正CT散射信号的方法，图1示出该方法的流程示意图，包括：
步骤10，模拟光子从X射线源S通过被测物体M到达探测器D的散射信号。
图2是光子从X射线源S射出经过被侧物体M到达探测器D的过程示意图。其中，该散射信号的模拟可以是针对信号，也可以是针对信号的其它参数，本文对此不做限制。
该图示仅为了便于理解本发明，不用于其它限定。
上述步骤10具体包括下述步骤11-17。
步骤11，通过沃克混叠方法，用预先生成的4n维Sobol'点u＝(u1,u2,...,u4n)的第一个分量u1，从预先给定的能谱分布φ(E0)中抽样光子的初始能量E0，其中n表示预设的光子在被测物体内发生碰撞的次数。
其中，能谱分布φ(E0)如图3所示，该图示仅为示例，不同X射线源不同光子的能谱分布可以不同。
步骤12，记初始入射光到探测器D的立体角为Ω0，用4n-维Sobol'点u的第二个和第三个分量u2和u3在Ω0中均匀抽样单位入射方向
Figure BDA0003273083400000061
得到初始流强I0(u2,u3)以及初始入射X射线沿单位入射方向
Figure BDA0003273083400000062
运动与被测物体M的交点A0，得到光子以权重W0＝I0(u2,u3)沿
Figure BDA0003273083400000063
方向从A0开始在被测物体M中运动。
本步骤中，设Ω0是初始入射光到探测器D的立体角，
Figure BDA0003273083400000064
是初始光流强，
Figure BDA0003273083400000065
是初始入射光的单位入射方向。
步骤13，在确定第i-1阶相互作用点Ai-1,i＝2,…,n的位置后，抽样光子在Ai-1的散射类型和散射角方向；包括：使用所述4n维Sobol’点u的第4(i-1)+1个分量u4(i-1)+1，通过重要性抽样方法，抽样散射类型，并使用u的第4(i-1)+2和4(i-1)+3个分量u4(i-1)+2和u4(i-1)+3，通过RITA算法抽样散射单位方向
Figure BDA0003273083400000071
光子从当前位置Ai-1开始，沿着单位散射方向
Figure BDA0003273083400000072
逃出被测物体M的概率为：
Figure BDA0003273083400000073
其中，
Figure BDA0003273083400000074
是光子的初始单位入射方向，
Figure BDA0003273083400000075
是光子在Ai-1发生碰撞后的散射方向，δi-1＝0,1，
Figure BDA0003273083400000076
和
Figure BDA0003273083400000077
分别表示光子在Ai-1发生康普顿和瑞利散射后的剩余能量，ci-1是光子从Ai-1运动到Ci-1的长度，Ci-1是光子沿
Figure BDA0003273083400000078
方向运动到与被测物体M边界上的交点，
Figure BDA0003273083400000079
μtot表示线性衰减系数；主信号，即光子从光源S出发，沿单位入射方向
Figure BDA00032730834000000710
逃出被测物体M并到达探测器D的概率，是
Figure BDA00032730834000000711
光子以权重
Figure BDA00032730834000000712
在被测物体M中进行第i次相互作用。
本步骤中，用线性衰减系数μtot(x,E)来反映光子与被测物体的相互作用。线性衰减系数μtot(x,E)表示光子穿过单位路径与被测物体发生相互作用的概率，是关于空间坐标x和能量E的函数，且
Figure BDA00032730834000000713
其中
Figure BDA00032730834000000714
是康普顿散射系数，
Figure BDA00032730834000000715
是瑞利散射系数，μa(x,E)是光电吸收系数。
记光子从其当前位置Ai-1到下一个相互作用点Ai的路径长度为si,i＝1,…,n，其遵循以下分布
Figure BDA00032730834000000716
其中
Figure BDA00032730834000000717
Figure BDA00032730834000000718
是光子在相互作用点Ai-1发生碰撞后的单位散射方向，δi-1＝0,1，
Figure BDA00032730834000000719
和
Figure BDA00032730834000000720
分别是光子在Ai-1发生康普顿和瑞利散射后的剩余能量；t是积分变量。
步骤14，基于下式确定下一个相互作用点Ai(i＝1,2,…,n)：
Figure BDA0003273083400000081
其中u4(i-1)+4是4n维Sobol'点的第(4(i-1)+4)个分量。
容易理解，在确定下一个相互作用点后，可以再根据步骤13的方式确定光子从该下一个相互作用点碰撞后逃出被测物体M的概率。
步骤15，确定光子在Ai,i＝1,2,…,n出发生第i次相互作用后，沿单位强制方向
Figure BDA0003273083400000082
到达固定探测器像素Dj,j＝1,…,m×m的概率：
Figure BDA0003273083400000083
其中，
Figure BDA0003273083400000084
和
Figure BDA0003273083400000085
是光子在相互作用点Ai发生康普顿散射和瑞利散射的概率，
Figure BDA0003273083400000086
是探测器像素Dj与Ai相对应的立体角，
Figure BDA0003273083400000087
Figure BDA0003273083400000088
是探测器像素Dj的面积，αi表示
Figure BDA0003273083400000089
与探测器像素Dj法线方向的夹角,
Figure BDA00032730834000000810
是光子沿单位强制方向
Figure BDA00032730834000000811
运动与被测物体M交点，
Figure BDA00032730834000000812
和
Figure BDA00032730834000000813
是康普顿散射和瑞利散射极角的概率密度，且是归一化的线性衰减系数
Figure BDA00032730834000000814
Figure BDA00032730834000000815
本步骤中，光子在Ai,i＝1,2,…,n处发生第i次散射后，理论上将以一定的概率从被测物体M逃出并到达探测器D。为了提高效率，只考虑理论上到达固定探测器像素的光子，这种技术称为强制固定检测(fixed forced detection,FFD)。
如图2所示，光子在Ai,i＝1,2,…,n出发生第i次相互作用后，光子沿单位强制方向
Figure BDA0003273083400000091
以一定概率到达固定探测器像素Dj,j＝1,...,m×m，该概率可以通过上式(5)得到。
步骤16，获得光子沿路径li,j:S→A0→A1→…→Ai→Bi,j→Dj,i∈{1,...,n}的概率为：
P(li,j)＝WiPi(Ai→Dj)， (8)
光子沿路径S→A0→A1的概率与预先生成的4n维Sobol’点u的分量u1、u2、u3和u4相关，故：
P(l1,j):＝P1(A1→Dj)g0(u1,...,u4)， (9)
其中，光子运动路径S→A0→A1→…→Ai是马尔科夫链，光子运动路径S→A0→A1由u1、u2、u3和u4抽样所得，光子由Ai-1运动到下一阶相互作用点Ai,i＝2,...,n由u4(i-1)+1、u4(i-1)+2、u4(i-1)+3和u4(i-1)+4抽样所得，且相互作用点Ai只与上一阶相互作用点Ai-1有关，因此得到：
Figure BDA0003273083400000092
步骤17，获得光子在被测物体M内发生n次碰撞后，固定探测器像素Dj,j＝1,...,m×m上累计接收到的散射总信号为
Figure BDA0003273083400000093
因此，fn,j是与u＝(u1,u2,…,u4n)有关的函数，在gQMCFFD(基于GPU的拟蒙特卡罗(quasi-Monte Carlo，QMC)方法与强制固定探测(fixed forceddetection，FFD)技术的结合)方法下，散射估计量的问题构造维数为d＝4n。
Figure BDA0003273083400000094
的拟蒙特卡罗估计量如下
Figure BDA0003273083400000101
其中
Figure BDA0003273083400000102
是取值[0,1]d上的Sobol’点。
步骤20，将每个探测器像素检测到的信号减去模拟得到的散射信号与预设系数的乘积，根据差值确定穿过被测模体M后的X射线信号。
校正投影数据中的散射信号，并重建CT图像的过程如下：1，读取投影数据；FDK算法重建得到体数据图像2，使用发明的基于GPU的拟蒙特卡罗方法和强制探测器技术相结合的散射模拟方法，快速模拟出散射信号；3，把每个探测器像素接收到的信号减去估计到的散射信号；4，光束硬化校正；5，使用校正散射污染的投影数据，根据FDK算法重建CT图像。6，重复步骤2-5，直到图像收敛。
在计算光子散射信号所需要的数据(例如几何参数、X射线光谱和物理数据)初始化仿真。所有的这些数据传输到GPU内存中。具体来说，由于GPU的缓存功能，体素化的模体数据和衰减系数被存储在GPU的纹理存储器中以实现快速访问。DCS数据存储在GPU全局内存中，其余变量存储在GPU常量内存中。所有图像计数器都设置为零。之后，以批处理方式进行仿真，进而使用批处理之间的结果来估计不确定性。
本发明实施例提供的上述方法可以简化为如下步骤：
1.赋n具体值，预先生成4n维Sobol’点u＝(u1,u2,…,u4n)；
2.用u1，通过沃克混叠(Walker's aliasing)方法从预先给定能谱抽样初始能量E0；用u2和u3抽样初始入射X射线的单位方向
Figure BDA0003273083400000103
并获得初始流强I0(u2,u3)以及计算初始入射X射线沿
Figure BDA0003273083400000104
方向运动与被测物体M的交点A0。令i＝1，通过(4)式，用u4抽样第一阶相互作用点A1；
3.如果i≥n，跳到步骤5；否则令i＝i+1，用u4(i-1)+1去抽样光子在相互作用点Ai-1,i＝2,…,n的散射类型，用u4(i-1)+2和u4(i-1)+3，通过RITA算法抽样散射方向
Figure BDA0003273083400000105
计算
Figure BDA0003273083400000106
和剩余能量
Figure BDA0003273083400000107
4.用u4(i-1)+4通过(4)式抽样相互作用点Ai；
5.令
Figure BDA0003273083400000108
计算P(li,j)＝WiPi(Ai→Dj)。如果i＜n，跳到步骤3，否则，跳到步骤6；
6.返回P(li,j),i＝1,...,n,j＝1,...,m×m。
为了说明本发明实施例提供的校正方法的效果，下面提供实验结果。
分别计算X射线穿过均匀铝模体和人头模体后到达探测器D的散射信号。实验所用的铝模体体积为160×28×160mm3。在计算中化成160×28×160个体素处理，每个体素大小为1×1×1mm3。所用的人头模体由软组织、骨头和空气组成，体积为135×200×215mm3。在计算过程中化成270×400×430个体素处理，每个体素大小为0.5×0.5×0.5mm3。所有均匀铝模体和人头模体的几何示意图见图4，其中a为均匀铝模体，b为人头模体。所用探测器由512×512个探测器像素组成，每个探测器的大小为0.8×0.8mm2。所用能谱为16-120keV，X射线源S到被测物体M的距离和M到探测器D的距离都是500mm。
通过计算相对差异(RD)和有效提高因子(EIF)来说明gQMCFFD算法的有效性和稳健性。
RD＝||r-t||2/||r||2， (13)
其中r是由gMCFRD算法计算探测器D接收到的散射信号，t可代表由gQMCFFD算法、gMCFFD算法或gQMCFRD算法计算探测器D接收到的散射信号。与gQMCFFD算法比较，gMCFFD算法用随机序列代替Sobol’序列进行抽样。gMCFRD算法基于GPU的MC方法和FRD技术相结合、gQMCFRD算法是基于GPU的QMC方法和FRD技术相结合,详见参考文献[14]。
Figure BDA0003273083400000111
其中
Figure BDA0003273083400000112
T表示运行时间，σ是计算结果的标准差，在实际计算中计算运行100次后的标准差。FOMgMCFRD表示gMCFRD算法的时间和方差成绩的倒数。A可取gQMCFFD、gMCFFD或gQMCFRD。EIF大于1意味着算法A相对于gQMCFRD更有效，如果EIF小于1，则相反。模拟中，gMCFRD、gQMCFRD、gMCFFD、gQMCFFD分别使用229、229、215、215个光子。硬件方面我们使用配备4352个处理器和11.0GB全局内存的GeForce RTX 2080Ti显卡。在模拟中，使用了Geant4(GEometryANdTracking)的PENELOPE物理数据库。
下述表1和表2分别显示了对均匀铝模体和人头模体，gMCFRD运行229个光子，gQMCFRD运行229个光子，gMCFFD运行215光子，gQMCFFD运行215光子的运行时间T，标准差σ(运行一百次)，以及将gQMCFRD、gMCFFD、gQMCFFD的计算结果与gMCFRD的计算结果进行比较的相对差异(RD)和EIF。从表1和表2可以看出，对铝模体gQMCFFD的有效因子提高46倍，对人头模体gQMCFFD的有效因子提高4倍，而RD小于1.5％。
表1示出关于均匀铝模体，gQMCFRD、gMCFFD、gQMCFFD方法计算的散射信号与gMCFRD比较，其中，N表示模拟次数，T表示运行时间，RD表示其余三种方法的计算结果与gMCFRD方法的计算结果的相对差异，FOM等于运行时间和方差乘积的倒数，EIF为有效提高因子。
表2示出关于人头模体，gQMCFRD、gMCFFD、gQMCFFD方法计算的散射信号与gMCFRD比较，其中，N表示模拟次数，T表示运行时间，RD表示其余三种方法的计算结果与gMCFRD方法的计算结果的相对差异，FOM等于运行时间和方差乘积的倒数，EIF为有效提高因子。
方法 log<sub>2</sub>(N) T(分钟) RD σ FOM EIF
gMCFRD 29 1.51 \ 1.14×10<sup>-2</sup> 5096 \
gQMCFRD 29 1.83 1.22％ 2.38×10<sup>-3</sup> 96470 18.93
gMCFFD 15 0.44 1.29％ 1.14×10<sup>-2</sup> 17488 3.43
gQMCFFD 15 0.48 1.17％ 2.98×10<sup>-3</sup> 234599 46.04
表1
方法 log<sub>2</sub>(N) T(分钟) RD σ FOM EIF
gMCFRD 29 15.17 \ 1.17×10<sup>-2</sup> 482 \
gQMCFRD 29 18.73 1.38％ 6.29×10<sup>-3</sup> 1362 2.83
gMCFFD 15 8.55 1.68％ 1.35×10<sup>-2</sup> 642 1.33
gQMCFFD 15 8.63 1.33％ 7.74×10<sup>-3</sup> 1934 4.01
表2
图5和图6分别展示了对均匀铝模体和人头模体，gMCFRD、gQMCFRD、gMCFFD和gQMCFFD四种方法计算的散射信号。图像视觉上的高度一致性说明了gQMCFFD方法的准确性。图5(a)-(d)分别示出gMCFRD、gQMCFRD、gMCFFD和gQMCFFD对Al模体计算的总散射信号结果，(e)对应上述4种方法，探测器测器中心且平行于x轴位置接收到的散射信号曲线图；(f)示出主信号。图6(a)-(d)分别示出gMCFRD、gQMCFRD、gMCFFD和gQMCFFD对人头模体计算的总散射信号结果，(e)对应上述4种方法，探测器测器中心且平行于x轴位置接收到的散射信号曲线图，(f)示出主信号。
通过采用本发明实施例提供的方法，计算光子理论上发生散射的概率，用公式(4)抽样下一阶相互作用点，这不同于已有的基于MC的散射校正方法。gQMCFFD算法使用Sobol’序列抽样，并在光子与被测物体发生碰撞后，计算其到达每个固定探测器像素的理论上的概率，把QMC方法的优势和强制固定探测(FFD)技术结合起来，能够准确地得到探测器探测到的光子的散射信号。
最后需要指出的是：以上实施例仅用以说明本发明的技术方案，而非对其限制。本领域的普通技术人员应当理解：可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (2)

1.一种基于拟蒙特卡罗和强制探测的校正CT散射信号的方法，其特征在于，包括：
步骤10，模拟光子从X射线源S通过被测物体M到达探测器D的散射信号；包括：
步骤11，通过沃克混叠方法，用预先生成的4n维Sobol点u＝(u1，u2，...，u4n)的第一个分量u1，从预先给定的能谱分布φ(E0)中抽样光子的初始能量E0，其中n表示光子在被测物体内发生碰撞的次数；
步骤12，记初始入射光到探测器D的立体角为Ω0，用4n-维Sobol点u的第二个和第三个分量u2和u3在Ω0中均匀抽样单位入射方向
Figure FDA0003665620490000011
得到初始流强I0(u2，u3)以及初始入射X射线沿单位入射方向
Figure FDA0003665620490000012
运动与被测物体M的交点A0，得到光子以权重W0＝I0(u2，u3)沿
Figure FDA0003665620490000013
方向从A0开始在被测物体M中运动；
步骤13，在确定第i-1阶相互作用点Ai-1，i＝2，…，n的位置后，抽样光子在Ai-1的散射类型和散射角方向；包括：使用所述4n维Sobol点u的第4(i-1)+1个分量u4(i-1)+1，通过重要性抽样方法，抽样散射类型，并使用u的第4(i-1)+2和4(i-1)+3个分量u4(i-1)+2和u4(i-1)+3，通过RITA算法抽样散射单位方向
Figure FDA0003665620490000014
光子从当前位置Ai-1，i＝1，…，n开始，沿着单位散射方向
Figure FDA0003665620490000015
逃出被测物体M的概率为：
Figure FDA0003665620490000016
其中
Figure FDA0003665620490000017
是光子初始单位入射方向，
Figure FDA0003665620490000018
是光子在Ai-1发生碰撞后的散射方向，δi-1＝0，1，
Figure FDA0003665620490000019
和
Figure FDA00036656204900000110
分别表示光子在Ai-1发生康普顿和瑞利散射后的剩余能量，ci-1是光子从Ai-1运动到Ci-1的长度，Ci-1是光子沿
Figure FDA00036656204900000111
方向运动到与被测物体M边界上的交点，
Figure FDA0003665620490000021
μtot表示线性衰减系数；主信号，即光子从光源S出发，沿单位入射方向
Figure FDA0003665620490000022
逃出被测物体M并到达探测器D的概率，是
Figure FDA0003665620490000023
光子以权重
Figure FDA0003665620490000024
在被测物体M中进行第i次相互作用；
步骤14，基于下式确定下一个相互作用点Ai(i＝1，2，..．，n)：
Figure FDA0003665620490000025
其中si，i＝1，..．，n表示位置Ai-1到下一个相互作用点Ai的路径长度，
Figure FDA0003665620490000026
u4(i-1)+4是4n维Sobol点的第(4(i-1)+4)个分量；
步骤15，确定光子在Ai，i＝1，2，...，n出发生第i次相互作用后，沿单位强制方向
Figure FDA0003665620490000027
到达固定探测器像素Dj，j＝1，...，m×m的概率：
Figure FDA0003665620490000028
其中，
Figure FDA0003665620490000029
和
Figure FDA00036656204900000210
是光子在相互作用点Ai发生康普顿散射和瑞利散射的概率，
Figure FDA00036656204900000211
是探测器像素Dj与Ai相对应的立体角，
Figure FDA00036656204900000212
Figure FDA00036656204900000213
是探测器像素Dj的面积，αi表示
Figure FDA00036656204900000214
与探测器像素Dj法线方向的夹角，
Figure FDA00036656204900000215
是光子沿单位强制方向
Figure FDA00036656204900000216
运动与被测物体M交点，
Figure FDA00036656204900000217
和
Figure FDA00036656204900000218
是康普顿散射和瑞利散射极角的概率密度，且是归一化的线性衰减系数
Figure FDA00036656204900000219
Figure FDA00036656204900000220
步骤16，获得光子沿路径li，j：S→A0→A1→...→Ai→Bi，j→Dj，i∈{1，...，n}的概率为：
P(li，j)＝WiPi(Ai→Dj)，
光子沿路径S→A0→A1的概率与预先生成的4n维Sobol点u的分量u1、u2、u3和u4相关，故：
P(l1，j)；＝P1(A1→Dj)g0(u1，...，u4)，
其中，光子运动路径S→A0→A1→...→Ai是马尔科夫链，光子运动路径S→A0→A1由u1、u2、u3和u4抽样所得，光子由Ai-1运动到下一阶相互作用点Ai，i＝2，..．，n由u4(i-1)+1、u4(i-1)+2、u4(i-1)+3和u4(i-1)+4抽样所得，且相互作用点Ai只与上一阶相互作用点Ai-1有关，因此得到：
Figure FDA0003665620490000031
步骤17，获得光子在被测物体M内发生n次碰撞后，固定探测器像素Dj，j＝1，...，m×m上累计接收到的散射总信号为
Figure FDA0003665620490000032
所述方法还包括步骤20，将每个探测器像素检测到的信号减去模拟得到的散射信号与预设系数的乘积，根据差值确定穿过被测模体M后的X射线信号。
2.如权利要求1所述的基于拟蒙特卡罗和强制探测的校正CT散射信号的方法，其特征在于，所述线性衰减系数μtot(x，E)表示光子穿过单位路径与被测物体发生相互作用的概率，是关于空间坐标x和能量E的函数，且
Figure FDA0003665620490000033
其中
Figure FDA0003665620490000034
是康普顿散射系数，
Figure FDA0003665620490000035
是瑞利散射系数，μa(x，E)是光电吸收系数。
CN202111107623.5A 2021-09-22 2021-09-22 基于拟蒙特卡罗和强制探测的校正ct散射信号的方法 Active CN113804709B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111107623.5A CN113804709B (zh) 2021-09-22 2021-09-22 基于拟蒙特卡罗和强制探测的校正ct散射信号的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111107623.5A CN113804709B (zh) 2021-09-22 2021-09-22 基于拟蒙特卡罗和强制探测的校正ct散射信号的方法

Publications (2)

Publication Number Publication Date
CN113804709A CN113804709A (zh) 2021-12-17
CN113804709B true CN113804709B (zh) 2022-08-12

Family

ID=78939950

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111107623.5A Active CN113804709B (zh) 2021-09-22 2021-09-22 基于拟蒙特卡罗和强制探测的校正ct散射信号的方法

Country Status (1)

Country Link
CN (1) CN113804709B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117796827A (zh) * 2022-09-26 2024-04-02 同方威视技术股份有限公司 用于成像设备的标定方法、装置、成像设备

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107202805A (zh) * 2017-05-31 2017-09-26 中国人民解放军信息工程大学 基于卷积核的锥束ct散射伪影校正方法
CN110702706A (zh) * 2019-09-20 2020-01-17 天津大学 一种能谱ct系统输出数据的模拟方法

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP3409742B2 (ja) * 1999-03-19 2003-05-26 日本軽金属株式会社 モンテカルロシミュレーションを用いたｅｐｍａ分析法
GB201115419D0 (en) * 2011-09-07 2011-10-19 Univ Leuven Kath Non-invasive in-situ radiation dosimetry
CN108618796A (zh) * 2018-02-09 2018-10-09 南方医科大学 一种基于任务驱动的蒙特卡洛散射光子模拟方法
CN110009707B (zh) * 2019-04-09 2023-01-03 上海联影医疗科技股份有限公司 图像散射校正方法、装置、计算机设备和存储介质
CN112763519B (zh) * 2020-12-22 2022-03-01 清华大学 一种用拟蒙特卡罗方法模拟光子散射的方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107202805A (zh) * 2017-05-31 2017-09-26 中国人民解放军信息工程大学 基于卷积核的锥束ct散射伪影校正方法
CN110702706A (zh) * 2019-09-20 2020-01-17 天津大学 一种能谱ct系统输出数据的模拟方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Quasi-Monte Carlo method for calculating X-ray scatter in CT;Lin, GY 等;《OPTICS EXPRESS》;20210421;全文 *
X射线横向移动成像法识别土壤中物体的蒙特卡罗模拟;朱鹏飞等;《核电子学与探测技术》;20111220(第12期);全文 *

Also Published As

Publication number Publication date
CN113804709A (zh) 2021-12-17

Similar Documents

Publication Publication Date Title
CN108292428B (zh) 图像重建的系统和方法
US8488902B2 (en) Detection values correction apparatus
Kyriakou et al. Combining deterministic and Monte Carlo calculations for fast estimation of scatter intensities in CT
Lasio et al. Statistical reconstruction for x-ray computed tomography using energy-integrating detectors
US8553831B2 (en) Imaging system and method using primary and scattered radiations
US10255696B2 (en) System and method for image reconstruction
US7889834B2 (en) Method for preparing reconstructed CT image data records and CT system
CN107802280B (zh) 校正曲线生成方法、投影图像的校正方法、系统及存储介质
EP2702563B1 (en) Multi-energy imaging
EP1700137B1 (en) Correction of artifacts caused by the heel effect
US11475611B2 (en) System and method for image reconstruction
EP2652709B1 (en) Imaging system for imaging a region of interest
JP2017528735A (ja) 材料特性評価のための断層撮影再構成
US10593070B2 (en) Model-based scatter correction for computed tomography
JP7341879B2 (ja) 医用画像処理装置、ｘ線コンピュータ断層撮影装置及びプログラム
Colijn et al. Experimental validation of a rapid Monte Carlo based micro-CT simulator
CN113804709B (zh) 基于拟蒙特卡罗和强制探测的校正ct散射信号的方法
US7379527B2 (en) Methods and apparatus for CT calibration
La Rivière Approximate analytic reconstruction in x-ray fluorescence computed tomography
US20130343510A1 (en) Method of reconstructing image for polychromatic x-ray tomography
Baer et al. Scatter correction methods in dimensional CT
US11353411B2 (en) Methods and systems for multi-material decomposition
Danyk et al. PHYSICAL BASES FOR DETERMINATION OF SCATTERING KERNELS FROM INCOMPLETE DATA IN GRID-LESS X-RAY IMAGING.
Vijayalakshmi et al. Comparison of algebraic reconstruction methods in computed tomography
Michielsen et al. Iodine quantification in limited angle tomography

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