CN110009707A - 图像散射校正方法、装置、计算机设备和存储介质 - Google Patents
图像散射校正方法、装置、计算机设备和存储介质 Download PDFInfo
- Publication number
- CN110009707A CN110009707A CN201910281959.XA CN201910281959A CN110009707A CN 110009707 A CN110009707 A CN 110009707A CN 201910281959 A CN201910281959 A CN 201910281959A CN 110009707 A CN110009707 A CN 110009707A
- Authority
- CN
- China
- Prior art keywords
- photon
- detector
- string
- image
- scattering
- 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
Links
- 239000006185 dispersion Substances 0.000 title claims abstract description 20
- 238000003860 storage Methods 0.000 title claims description 13
- 238000000342 Monte Carlo simulation Methods 0.000 claims abstract description 90
- 230000005855 radiation Effects 0.000 claims abstract description 61
- 238000005070 sampling Methods 0.000 claims abstract description 60
- 238000000034 method Methods 0.000 claims abstract description 47
- 238000012937 correction Methods 0.000 claims abstract description 40
- 238000005259 measurement Methods 0.000 claims abstract description 34
- 238000001514 detection method Methods 0.000 claims description 90
- 238000004590 computer program Methods 0.000 claims description 43
- 238000004088 simulation Methods 0.000 claims description 28
- 238000012545 processing Methods 0.000 claims description 25
- 230000002123 temporal effect Effects 0.000 claims description 23
- 238000004422 calculation algorithm Methods 0.000 abstract description 17
- 238000004364 calculation method Methods 0.000 abstract description 15
- 238000009826 distribution Methods 0.000 description 19
- 239000000126 substance Substances 0.000 description 16
- 230000005251 gamma ray Effects 0.000 description 12
- 238000010586 diagram Methods 0.000 description 11
- 230000008569 process Effects 0.000 description 6
- 230000003993 interaction Effects 0.000 description 5
- 150000002500 ions Chemical class 0.000 description 5
- 230000005012 migration Effects 0.000 description 4
- 238000013508 migration Methods 0.000 description 4
- 238000002591 computed tomography Methods 0.000 description 3
- 239000013078 crystal Substances 0.000 description 2
- 238000002059 diagnostic imaging Methods 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 238000011156 evaluation Methods 0.000 description 2
- 238000005457 optimization Methods 0.000 description 2
- 238000005086 pumping Methods 0.000 description 2
- 230000009467 reduction Effects 0.000 description 2
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 2
- WQZGKKKJIJFFOK-GASJEMHNSA-N Glucose Natural products OC[C@H]1OC(O)[C@H](O)[C@@H](O)[C@@H]1O WQZGKKKJIJFFOK-GASJEMHNSA-N 0.000 description 1
- 240000007594 Oryza sativa Species 0.000 description 1
- 235000007164 Oryza sativa Nutrition 0.000 description 1
- 230000002776 aggregation Effects 0.000 description 1
- 238000004220 aggregation Methods 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000006399 behavior Effects 0.000 description 1
- 230000008901 benefit Effects 0.000 description 1
- 230000008859 change Effects 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 238000003745 diagnosis Methods 0.000 description 1
- 235000014113 dietary fatty acids Nutrition 0.000 description 1
- 230000002708 enhancing effect Effects 0.000 description 1
- 229930195729 fatty acid Natural products 0.000 description 1
- 239000000194 fatty acid Substances 0.000 description 1
- 150000004665 fatty acids Chemical class 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 238000010304 firing Methods 0.000 description 1
- 239000008103 glucose Substances 0.000 description 1
- 238000009499 grossing Methods 0.000 description 1
- 230000036541 health Effects 0.000 description 1
- 238000003702 image correction Methods 0.000 description 1
- 238000003709 image segmentation Methods 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000001727 in vivo Methods 0.000 description 1
- 239000004973 liquid crystal related substance Substances 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 230000004060 metabolic process Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- VIKNJXKGJWUCNN-XGXHKTLJSA-N norethisterone Chemical compound O=C1CC[C@@H]2[C@H]3CC[C@](C)([C@](CC4)(O)C#C)[C@@H]4[C@@H]3CCC2=C1 VIKNJXKGJWUCNN-XGXHKTLJSA-N 0.000 description 1
- 238000009206 nuclear medicine Methods 0.000 description 1
- 102000039446 nucleic acids Human genes 0.000 description 1
- 108020004707 nucleic acids Proteins 0.000 description 1
- 150000007523 nucleic acids Chemical class 0.000 description 1
- 102000004169 proteins and genes Human genes 0.000 description 1
- 108090000623 proteins and genes Proteins 0.000 description 1
- 235000009566 rice Nutrition 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 238000000528 statistical test Methods 0.000 description 1
- 230000001360 synchronised effect Effects 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
- 230000007704 transition Effects 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Classifications
-
- G06T5/70—
-
- 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
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10072—Tomographic images
- G06T2207/10104—Positron emission tomography [PET]
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Measurement Of Radiation (AREA)
- Nuclear Medicine (AREA)
Abstract
本申请涉及一种图像散射校正方法，方法包括：获取图像原始数据，所述图像原始数据包括衰减图像以及测量弦图；根据所述测量弦图获得放射图像；根据衰减图像以及放射图像进行蒙特卡罗模拟，得到散射弦图；对散射弦图进行去噪声处理，得到低噪声散射弦图；根据图像原始数据对低噪声散射弦图进行统计缩放，得到校正散射弦图。通过对蒙特卡罗模拟中各项随机抽样算法进行优化，能够提高蒙特卡罗模拟的计算效率，减小计算量。通过对散射弦图进行去噪声处理，降低对蒙特卡罗模拟算法的计算量需求，进一步的节省了计算资源以及计算时间，并能够更好的满足临床PET扫描对扫描时间的限制要求。
Description
技术领域
本申请涉及医学成像技术领域，特别是涉及一种图像散射校正方法、装置、计算机设备和存储介质。
背景技术
蒙特卡罗模拟方法，又被称为随机抽样或统计试验方法。蒙特卡罗模拟方法的基本思想就是当所要求解的问题是某种事件出现的概率，或者是某个随机变量的期望值时，它们可以通过某种“试验”的方法，得到这种事件出现的频率，或者这个随机变数的平均值，并用它们作为问题的解。蒙特卡罗方法通过抓住事物运动的几何数量和几何特征，利用数学方法来加以模拟，即进行一种数字模拟实验。
目前蒙特卡罗模拟方法已经被应用到医学成像的技术领域，在PET(正电子发射型计算机断层显像)系统的散射校正中，由于蒙特卡罗模拟方法对PET系统的还原度高并且校正精准，因此经常被作为各种校正方法的参考标准。
但是由于对PET系统进行蒙特卡罗模拟所需要的计算量大，需要消耗大量的计算资源以及时间，难以满足临床PET扫描对扫描时间的限制要求。
发明内容
基于此，有必要针对上述技术问题，提供一种能够减少计算量并能够节省计算资源以及计算时间的图像散射校正方法、装置、计算机设备和存储介质。
一种图像散射校正方法，所述方法包括：获取图像原始数据，所述图像原始数据包括衰减图像以及测量弦图；根据所述测量弦图获得放射图像；根据所述衰减图像以及放射图像进行蒙特卡罗模拟，得到散射弦图；对所述散射弦图进行去噪声处理，得到低噪声散射弦图；根据图像原始数据对所述低噪声散射弦图进行统计缩放，得到校正散射弦图。
在其中一个实施例中，根据图像原始数据对所述低噪声散射弦图进行统计缩放，得到校正散射弦图之后，还包括：根据所述校正散射弦图以及测量弦图进行图像重建，得到校正放射图像；根据预设迭代次数，多次获取校正散射弦图并对所述校正放射图像进行迭代校正，得到最终放射图像。
在其中一个实施例中，所述根据所述衰减图像以及放射图像进行蒙特卡罗模拟，得到散射弦图包括：根据所述衰减图像以及放射图像确定蒙特卡罗模拟边界以及初始的光子信息，所述图像原始数据包括：衰减图像以及测量弦图，光子信息包括：光子位置、光子方向以及光子能量；根据所述蒙特卡罗模拟边界以及初始的光子信息进行光子运动模拟，得到光子到达探测器的探测信息；根据所述探测信息，得到散射弦图。
在其中一个实施例中，所述根据所述衰减图像以及放射图像确定蒙特卡罗模拟边界以及初始的光子信息包括：根据所述衰减图像设定蒙特卡罗边界；根据所述衰减图像以及放射图像，确定初始的光子信息。
在其中一个实施例中，所述根据所述蒙特卡罗模拟边界以及初始的光子信息进行光子运动模拟，得到光子到达探测器的探测信息包括：对光子迁移量进行随机抽样；根据抽样的光子迁移量、初始的光子信息以及蒙特卡罗模拟边界进行光子运动模拟，得到光子到达探测器的探测信息。
在其中一个实施例中，所述根据抽样的光子迁移量、初始的光子信息以及蒙特卡罗模拟边界进行光子运动模拟，得到光子到达探测器的探测信息包括：以预设步长移动相应光子；当光子移动距离未到达相应的光子迁移量，且光子移动位置达到蒙特卡罗模拟边界时，则将光子沿当前光子方向移动至探测器，根据光子信息得到光子到达探测器的探测信息，所述探测信息包括：散射标记、光子到达探测器的位置、光子到达探测器的方向、光子到达探测器的能量以及光子达到探测器的输运时间。
在其中一个实施例中，所述以预设步长移动相应光子之后还包括：当光子移动距离到达相应的光子迁移量，则根据康普散射角度抽样得到光子散射角度；根据所述光子散射角度，生成散射标记；根据光子散射角度以及光子迁移量更新光子信息；当更新后的光子信息所记录的光子能量大于第一预设阈值，则继续进行光子运动模拟，直到光子移动位置到达蒙特卡罗模拟边界。
在其中一个实施例中，所述以预设步长移动相应光子之后还包括：当光子移动距离未到达相应的光子迁移量，且光子移动位置未达到蒙特卡罗模拟边界时，则继续以预设步长移动相应光子，直至相应光子移动距离到达光子迁移量或相应光子移动位置达到蒙特卡罗模拟边界。
在其中一个实施例中，所述根据所述探测信息，得到散射弦图包括：获取至少一个光子对的探测信息；根据至少一个光子对的所述探测信息，得到有效光子对的探测信息；根据有效光子对的探测信息，得到散射弦图。
在其中一个实施例中，所述根据至少一个光子对的所述探测信息，得到有效光子对的探测信息包括：获取相应探测器的能量分辨率；当光子对中的两个光子都到达探测器，则对光子到达探测器的能量根据所述能量分辨率进行高斯随机处理；当高斯随机处理后光子对中的两个光子到达探测器的能量都高于第二预设阈值，则将相应光子对作为有效光子对；获取有效光子对的探测信息。
在其中一个实施例中，所述根据至少一个光子对的所述探测信息，得到有效光子对的探测信息包括：获取相应探测器的能量分辨率；当光子对中的两个光子都到达探测器，则对光子到达探测器的能量根据所述能量分辨率进行高斯随机处理；根据高斯随机处理后光子对中的两个光子到达探测器的能量、能量分辨率以及第二预设阈值，计算光子对被探测器接收的概率；根据所述光子对被探测器接收的概率，获取有效光子对的探测信息。
在其中一个实施例中，所述获取有效光子对的探测信息还包括：获取预设飞行时间重建标识；当所述预设飞行时间重建标识为进行飞行时间重建，则对有效光子对进行飞行时间重建；获取飞行时间重建后的有效光子对探测信息。
在其中一个实施例中，所述对有效光子对进行飞行时间重建包括：获取有效光子对中两个光子的光子达到探测器的输运时间；根据有效光子对中两个光子的光子达到探测器的输运时间进行时间分辨率高斯抽样；对时间分辨率高斯抽样的有效光子对，划分数据子集。
在其中一个实施例中，所述根据有效光子对中两个光子的光子达到探测器的输运时间进行时间分辨率高斯抽样包括：获取相应探测器的时间分辨率；根据有效光子对中两个光子的光子达到探测器的输运时间，计算有效光子对中两个光子到达探测器的时间差；对所述时间差根据所述探测器的时间分辨率进行时间分辨率高斯抽样。
一种图像散射校正装置，所述装置包括：获取模块，用于获取图像原始数据；图像重建模块，用于根据所述测量弦图获得放射图像；蒙特卡罗模拟模块，用于根据所述衰减图像以及放射图像进行蒙特卡罗模拟，得到散射弦图；去噪声处理模块，用于对所述散射弦图进行去噪声处理，得到低噪声散射弦图；统计缩放处理模块，用于根据图像原始数据对所述低噪声散射弦图进行统计缩放，得到校正散射弦图。
一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，所述处理器执行所述计算机程序时实现上述任一种所述方法的步骤。
一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现上述任一种所述的方法的步骤。
上述图像散射校正方法、装置、计算机设备和存储介质，通过获取图像原始数据，并根据原始数据中的测量弦图获得放射图像，根据衰减图像以及放射图像进行蒙特卡罗模拟得到散射弦图，再对得到的散射弦图进行去噪处理，得到低噪声散射弦图，利用图像原始数据对低噪声散射弦图进行统计缩放，得到校正散射弦图。通过对蒙特卡罗模拟中各项随机抽样算法进行优化，能够提高蒙特卡罗模拟的计算效率，减小计算量。通过对散射弦图进行去噪声处理，降低对蒙特卡罗模拟算法的计算量需求，进一步的节省了计算资源以及计算时间，并能够更好的满足临床PET扫描对扫描时间的限制要求。
附图说明
图1为一个实施例中图像散射校正方法的流程示意图；
图2为一个实施例中利用蒙特卡罗模拟得到散射弦图的方法流程示意图；
图3为一个实施例中获取光子到达探测器的探测信息方法流程示意图；
图4为一个实施例中光子迁移距离抽样结果的图示；
图5为一个实施例中300～600keV内衰减系数能量依赖关系的二阶多项式拟合曲线；
图6为一个实施例中不同能量伽马射线的康普散射角度抽样分布与理论分布的比较；
图7为一个实施例中根据探测信息计算散射弦图的方法流程示意图；
图8为一个实施例中为本发明计算所得弦图与GATE仿真弦图对比图；
图9为一个实施例中图像散射校正装置的结构框图；
图10为一个实施例中计算机设备的内部结构图。
附图标记：获取模块100、图像重建模块200、蒙特卡罗模拟模块300、去噪声处理模块400、统计缩放处理模块500。
具体实施方式
为了使本申请的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本申请进行进一步详细说明。应当理解，此处描述的具体实施例仅仅用以解释本申请，并不用于限定本申请。
正电子发射型计算机断层显像(Positron Emission Computed Tomography，PET)，是核医学领域比较先进的临床检查影像技术。是将某种物质，一般是生物生命代谢中必须的物质，如：葡萄糖、蛋白质、核酸、脂肪酸，标记上短寿命的放射性核素(如18F，11C等)，注入人体后，放射性核素在衰变过程中释放出正电子，一个正电子在行进十分之几毫米到几毫米后遇到一个电子后发生湮灭，从而产生方向相反的一对能量为511KeV的光子。这对光子，通过高度灵敏的照相机捕捉，并经计算机进行散射和随机信息的校正。经过对不同的正电子进行相同的分析处理，我们可以得到在生物体内聚集情况的三维图像，从而达到诊断的目的。
在一个实施例中，如图1所示，提供了一种图像散射校正方法，包括以下步骤：
步骤S102，获取图像原始数据。
具体地，图像原始数据包括：衰减图像以及测量弦图。根据衰减图像能够获取被测物体的几何信息。其中几何信息具体是指被测物体在PET系统中所处的位置以及占据的空间大小。更具体的，衰减图像可以直接由CT(电子计算机断层扫描设备)图像计算得到，也可以根据MR(磁共振成像设备)以图像分割算法计算得到。衰减图像中各像素点的值体现了伽马光子与被测物体发生相互作用的概率。测量弦图可以由PET图像的符合数据转化得到，测量弦图中各像素点的值为PET系统符合线上的符合计数。其中，每一条符合线表示可能发生符合关系的一对探测器晶体。
步骤S104，根据所述测量弦图获得放射图像。
步骤S106，根据衰减图像以及放射图像进行蒙特卡罗模拟，得到散射弦图。
具体地，散射弦图体现了散射事件在测量弦图中的强度分布。根据衰减图像以及放射图像设定蒙特卡罗模拟中湮灭光子的初始光子位置、初始光子方向以及初始光子能量。在根据湮灭光子的初始光子位置、初始光子方向以及初始光子能量进行模拟，最终得到散射弦图。
步骤S108，对散射弦图进行去噪声处理，得到低噪声散射弦图。
具体地，去噪声的目的是为了减小散射弦图上的统计噪声。也就是需要对蒙特卡罗模拟得到的散射弦图进行平滑处理。其中，去噪声处理可以采用多种卷积平滑处理或插值法进行平滑处理。选取去噪声处理方法的依据是散射事件在空间上是否为连续分布。优选的，可以采用中值滤波算法进行去噪声处理，以散射弦图像中各像素点附近一定范围数据的平均值代替元像素的值。去噪声包含了在正常散射弦图维度进行去噪声处理，也包括将散射弦图转换为michelogram后的去噪声处理。
步骤S110，根据图像原始数据对低噪声散射弦图进行统计缩放，得到校正散射弦图。
具体地，统计缩放的目的是为了使蒙特卡罗模拟所得到的散射弦图与图像原始数据中的测量弦图保持统计量的一致。更具体的，将蒙特卡罗模拟得到的散射弦图中真事件的总数与散射事件的总数求和，将测量弦图中真事件的总数与散射事件的总数求和。并将散射弦图的事件总和与测量弦图的事件总和的比值作为缩放系数。也可以将蒙特卡罗模拟得到的散射弦图与图像原始数据中的测量弦图做线性最小二乘拟合，得到缩放系数。还可以将散射弦图和测量弦图拆分为不同的片段分别处理，得到各个片段的缩放系数。根据缩放系数，得到校正散射弦图。
在其中一个实施例中，得到校正散射弦图之后，根据校正散射弦图以及测量弦图进行图像重建，得到校正放射图像。
在其中一个实施例中，根据校正散射弦图以及测量弦图进行图像重建，得到校正放射图像之后还包括：根据预设迭代次数，多次获取校正散射弦图并对放射图像进行迭代校正，得到最终放射图像。
上述图像散射校正方法、装置、计算机设备和存储介质，通过获取图像原始数据，并根据原始数据中的测量弦图获得放射图像，根据衰减图像以及放射图像进行蒙特卡罗模拟得到散射弦图，再对得到的散射弦图进行去噪处理，得到低噪声散射弦图，利用图像原始数据对低噪声散射弦图进行统计缩放，得到校正散射弦图。通过对蒙特卡罗模拟中各项随机抽样算法进行优化，能够提高蒙特卡罗模拟的计算效率，减小计算量。通过对散射弦图进行去噪声处理，降低对蒙特卡罗模拟算法的计算量需求，进一步的节省了计算资源以及计算时间，并能够更好的满足临床PET扫描对扫描时间的限制要求。
在一个实施例中，如图2所示，提供了一种利用蒙特卡罗模拟得到散射弦图的方法，包括以下步骤：
步骤S202：根据衰减图像以及放射图像确定蒙特卡罗模拟边界以及初始的光子信息。
具体地，光子信息包括：光子位置、光子方向以及光子能量。更具体的，根据衰减图像设定蒙特卡罗边界；根据衰减图像以及放射图像，确定初始的光子信息。其中，蒙特卡罗边界是指在该边界外不存在被测物体，只有空气或PET系统探测器。具体的，根据衰减图像各像素点的像素值，能够判断PET系统中哪些空间不存在被探测物体只有空气。在实际的操作过程中，蒙特卡罗边界一般被设定为方形或椭圆形。通过设定蒙特卡罗边界能够使蒙特卡罗模拟空间覆盖全部散射体，同时能够排除无散射体的视野，在光子运动模拟的过程中省去对边界外光子运动的模拟。现有技术在确定初始的光子位置时，需要根据放射图定义光子在每一个位置产生的概率，在对光子进行随机抽样。在本实施例中，初始的光子位置不经过随机抽样，而是采用遍历的方式，根据统计量依次遍历所有有效体素。每个有效体素内含有的湮灭光子对数为：
其中，有效体素是指根据衰减图像或放射图像，在该体素内确定存在被测量物体，并且被测量物体内带有放射性。其判断标准为对衰减图像或放射图像各自设定极小值阈值，在衰减图像中低于此极小值阈值的体素内部为空气，在放射图像中低于此极小值阈值的体素为噪声。当体素在衰减图像和放射图像中均大于极小值阈值时，当前体素为有效体素。优选的，将衰减图像和放射图像中像素点均不为零的体素作为有效体素。例如，放射图像中某体素的有效强度为0.1，放射图像总有效强度为100，蒙特卡罗计算统计量为10000，在在该体素位置存在10对湮灭光子对。其中，体素有效强度是指在放射图像中相应位置的像素值大小。总有效强度是指所有有效体素强度的总和。计算统计量为预先设定的值。蒙特卡罗计算统计量设定的越高，最终计算的结果越精确。该计算统计量可以设定为固定值，要求其满足所有被测量物体的计算精度需求；该计算统计量也可以根据被测量物体的大小进行动态设定，被测量物体越大，设定的计算统计量越大；被测量物体越小，设定的计算统计量越小。被测量物体的大小可以由衰减图像上的有效体素数量决定。初始的光子方向为：对一对湮灭光子对中的一个光子以各向同性的方向分布进行随机抽样，另一个光子方向与其相反。初始的光子能量一般设定为511keV。
步骤S204：根据蒙特卡罗模拟边界以及初始的光子信息进行光子运动模拟，得到光子到达探测器的探测信息。
具体地，对光子迁移量进行随机抽样得到抽样的光子迁移量，根据抽样的光子迁移量、初始的光子信息以及蒙特卡罗模拟边界进行光子运动模拟，得到光子到达探测器的探测信息。更具体地，以预设步长移动相应光子，当光子移动距离未到达相应的光子迁移量，且光子移动位置达到蒙特卡罗模拟边界时，则将光子沿当前光子方向移动至探测器，根据光子信息得到光子到达探测器的探测信息。光子运动模拟为蒙特卡罗模拟过程中计算量最大的部分，在本实施例中，对光子运动模拟中涉及的各项物理过程及抽样算法进行了加速优化。
步骤S206：根据探测信息，得到散射弦图。
具体地，获取至少一个光子对的探测信息，根据至少一个光子对的探测信息，得到有效光子对的探测信息，根据有效光子对的探测信息，得到散射弦图。
在一个实施例中，如图3所示，提供了一种获取光子到达探测器的探测信息的方法，包括以下步骤：
步骤S302：对光子迁移量进行随机抽样。
具体地，对光子迁移量进行随机抽样得到抽样的光子迁移量。PET系统针对的是复杂的人体系统，衰减图像和放射图像中的每一个像素点都有其特定的值，这使得光子在不同体素中发生与物质相互作用的概率不同。如果对每一个体素分别计算光子与物质的相互作用概率，再进行随机抽样，会增加大量的计算量，降低蒙特卡罗模拟的效率。在本实施例中，首先定义衰减系数与光子迁移距离的乘积为光子迁移量，然后直接对光子迁移量进行抽样。更具体的，根据光子在物质中的强度随着迁移距离l的分布密度P(l)满足指数分布：
P(1)＝μ·exp(-μl)
其中，P(l)表示迁移距离的分布密度，l表示光子迁移距离，μ表示物质的衰减系数。迁移距离的物理意义为光子从发射到与物质发生相互作用的距离，其值根据迁移量与衰减系数的比值确定。根据逆变换法则可以得到：
μl＝-ln(ξ)
其中，ξ表示在(0，1]内均匀分布的随机数。因为蒙特卡罗模拟中不同体素的衰减系数不同，设置当满足下述公式时，光子在当前体素中发生了与物质的相互作用。
∑μidl≥-ln(ξ)
其中，μi表示在不同位置读取到的衰减图上对应位置的衰减系数，dl表示设定的光子迁移步长。-ln(ξ)表示抽样的光子迁移量。
当判断光子发生相互作用时，对光子的位置沿当前移动方向的反方向回退一定距离Δl，其计算公式如下。
上式中，μn表示光子迁徙时调用的最后一个位置的衰减系数。对光子位置的回退，可以保证光子迁移距离的连续性，不会因步长的设置而离散化。
如图4，图4为一个实施例中光子迁移距离抽样结果的图示。以400mm内的光子迁移距离抽样为例。定义前200mm的衰减系数等同于水，后200mm的衰减系数为水的一半。以步长为20mm进行抽样计算，得到图4所示结果。
更具体的，衰减系数由衰减图像特定位置像素值及光子能量计算所得。在PET系统中衰减图像中给出的是针对511keV伽马射线在物质中的衰减系数。但衰减系数具有一定的能量依赖性，将衰减图像中的像素值直接应用至不同能量伽马射线会带来计算误差。因此，就需要根据511keV伽马射线在物质中的衰减系数计算各种能量的伽马射线在物质中的衰减系数。具体的，光子在物质中的衰减系数和光子与物质的相互作用截面有关。由于PET系统针对的光子能量一般在300～511keV，近似认为所有能量的光子与物质发生的相互作用均为康普顿散射，忽略其它物理过程。则有：
其中，σ表示光子与物质中原子相互作用的总截面；σc表示康普顿散射截面，NA表示阿伏伽德罗常数；ρ表示物质密度；MA表示摩尔质量；μ表示衰减系数。由上述公式可知：衰减系数可近似认为与康普顿散射截面成正比μ∝σc。从而可以通过不同能量光子的康普顿散射截面之比间接获得不同能量光子的衰减系数。根据康普顿散射积分截面，可知不同能量对应的衰减系数μ(E′)为：
其中，σc(E0)表示伽马射线在能量为E0＝511keV的康普顿散射截面；σc(E′)表示伽马射线在所求能量的康普顿散射截面；μ(E′)表示所求能量的衰减系数；μ(E0)为衰减图像给出的511keV伽马射线的衰减系数；α＝Eγ/mec2为伽马射线能量与电子静止质量之比。
在其中一个实施例中，由于上述公式中，σc(E′)/σc(E0)的计算过于复杂，直接计算会增加一定计算量。在PET系统内的散射模拟计算，只需要考虑300～511keV以内的伽马射线。因此，可以利用多项式拟合的方式，重新描述300～600keV能量区内物质衰减系数的能量依赖性。更具体的，也就是只需要简单的二项式拟合，即可保证300～600keV内对物质衰减系数的高精度拟合，大大减少了计算量。如图5所示，图5为一个实施例中300～600keV内衰减系数能量依赖关系的二阶多项式拟合曲线。多项式拟合足以对一定能量区域内σc(E′)/σc(E0)的能量依赖关系提供足够精确的描述，以减少PET模拟中的计算量。利用二阶多项式拟合结果，衰减系数μ(E′)的计算公式简化为：
μ(E′)＝μ(E0)×(1.7589-1.0994×α+0.3393×α2)
其中，α＝Eγ/mec2为伽马射线能量与电子静止质量之比；μ(E0)为衰减图像给出的511keV伽马射线的衰减系数。
步骤S304：以预设步长移动相应光子。
具体地，根据抽样得到的光子迁移量，以预先设定的步长移动相应的光子，并记录光子移动预设步长与所在位置光子能量的对应衰减系数的乘积。因为衰减系数由衰减图像特定位置像素值及光子能量计算所得，因此在每一次以预设步长移动光子后，需要重新计算光子能量及对应衰减图上的位置的衰减系数。
步骤S306：当光子移动距离未到达相应的光子迁移量，且光子移动位置达到蒙特卡罗模拟边界时，则将光子沿当前光子方向移动至探测器，根据光子信息得到光子到达探测器的探测信息。
具体地，在以预设步长移动相应光子之后，首先判断光子移动距离是否到达相应的光子迁移量，如果光子移动距离未到达相应的光子迁移量；再判断光子移动位置是否到达蒙特卡罗模拟边界，如果光子移动位置到达蒙特卡罗模拟边界，将光子沿当前光子方向移动至探测器，并获取光子到达探测器的探测信息。其中，探测信息包括：散射标记、光子到达探测器的位置、光子到达探测器的方向、光子到达探测器的能量以及光子达到探测器的输运时间。
在其中一个实施例中，判断光子移动距离是否到达相应的光子迁移量，当光子移动距离到达相应的光子迁移量，则根据康普散射角度抽样得到光子散射角度。更具体的，康普顿散射角θc的随机抽样，需依据Klein-Nishina公式给出的微分截面分布。其中，Klein-Nishina公式是用于计算电磁场对静止电子的康普顿散射微分散射截面与初、末态光子能量的公式。康普顿散射在散射角为θc上的微分截面dσc，e/dΩ为：
其中，α＝Eγ/mec2为伽马射线能量与电子静止质量之比，r0为经典电子半径。由于Klein-Nishina公式过于复杂，需要通过拒绝采样法进行随机抽样，但拒绝采样法需要合适的初始分布以提高抽样效率。
在其中一个实施例中还提供了一种能量相关的初始抽样分布，能够极大地提高抽样效率。具体的抽样分布为：
其中，α＝Eγ/mec2为伽马射线能量与电子静止质量之比。对康普顿散射角的θc的随机抽样为：
首先产生在[0，1]内均匀分布的随机数ξ，令θc＝2*arcsin(ξ)，则使得θc满足coS(θc/2)形式的随机分布。
再计算当前抽样随机角度的接收概率Pacc：
其中，公式中dσc，e/dΩ为Klein-Nishina公式给出的当前角度伽马射线的康普顿散射微分截面。
最后产生在[0，1]内均匀分布的随机数u，若u≤Pacc，则将当前θc为康普顿散射角，否则重新计算，直至产生满足要求的散射角。本实施例康普散射角度的抽样算法，可以保证对能量在300～511keV范围内伽马射线的康普顿散射角抽样准确度，满足对PET系统蒙特卡罗模拟的需求，并且抽样速度相对均匀的初始抽样分布能提高40％的蒙特卡罗模拟效率。如图6所示，图6为一个实施例中不同能量伽马射线的康普散射角度抽样分布与理论分布的比较。图中实线为理论分布，直方图为抽样分布。当100keV时，抽样结果才会明显偏离理论值。而对于PET系统，由于LLD(Low-Level Discriminator低电平甄别器)的存在，低能量的事件不会被系统记录，因此本算法不会影响蒙特卡罗模拟PET系统的精确性。
以预设步长移动相应光子，判断光子移动距离是否到达相应的光子迁移量，当光子移动距离到达相应的光子迁移量，则根据康普散射角度抽样得到光子散射角度。根据光子散射角度，生成散射标记。散射标记表示光子在运动模拟的过程中是否发生散射。当存在散射角度时，则散射标记为有散射标记；当不存在散射角度时，则散射标记为没有有散射标记。根据光子散射角度以及光子迁移量更新光子信息。根据光子的散射角度以及移动步长，更新当前光子的光子位置、光子方向以及光子能量，并在当前光子信息中添加散射标记。更具体的，当判断光子发生相互作用时，将光子的位置回退一定距离Δl。
其中，μn表示光子运动时调用的最后一个位置的衰减系数。对光子位置的回退，可以保证光子迁移距离的连续性，不会因步长的设置而离散化。
根据光子散射角度以及光子迁移量更新光子信息，当更新后的光子信息所记录的光子能量大于第一预设阈值，则继续进行光子运动模拟，直到光子移动位置到达蒙特卡罗模拟边界。更具体的，通过第一预设阈值对光子能量进行限定，当发生散射后的光子能量小于第一预设阈值，则停止对当前光子的运动的模拟。第一预设阈值的设定可以省去对不需要的低能量光子的模拟追踪，降低计算量。此处设定的第一预设阈值不等于PET系统本身的低能量阈值。在考虑系统存在一定的能量分辨后，此处设定的第一预设阈值ELLD应满足：
其中，为PET系统低能量阈值，E0为511keV，ΔE0为PET系统探测器对511keV伽马射线的能量分辨。
在其中一个实施例中，以预设步长移动相应光子，判断光子移动距离是否到达相应的光子迁移量，如果光子移动距离未到达相应的光子迁移量；再判断光子移动位置是否到达蒙特卡罗模拟边界，如果光子移动位置未到达蒙特卡罗模拟边界，则继续以预设步长移动相应光子，直至相应光子移动距离到达光子迁移量或相应光子移动位置达到蒙特卡罗模拟边界。也就是再次以预设步长移动相应光子，直至相应光子移动距离到达光子迁移量或相应光子移动位置达到蒙特卡罗模拟边界。
在一个实施例中，如图7所示，提供了一种探测信息计算散射弦图的方法，包括以下步骤：
步骤S402：获取至少一个光子对的探测信息。
步骤S404：根据至少一个光子对的探测信息，得到有效光子对的探测信息。
具体的，获取相应探测器的能量分辨率；当光子对中的两个光子都到达探测器，则对光子到达探测器的能量根据能量分辨率进行高斯随机处理；当高斯随机处理后光子对中的两个光子到达探测器的能量都高于第二预设阈值，则将相应光子对作为有效光子对；获取有效光子对的探测信息。更具体的，首先判断光子对中的两个光子是否都到达探测器位置，此处需结合实际探测其的位置分布，若光子所在位置无探测晶体存在，如探测器的模块缝隙，则不进行后续计算。其次需要对光子能量进行高斯随机处理，模拟实际探测器的能量分辨。若随机处理后光子能量高于第二预设阈值，则判断这一对光子对为有效事件。第二预设阈值为PET系统低能量阈值。
在其中一个实施例中，获取相应探测器的能量分辨率；当光子对中的两个光子都到达探测器，则对光子到达探测器的能量根据能量分辨率进行高斯随机处理；根据高斯随机处理后光子对中的两个光子到达探测器的能量、能量分辨率以及第二预设阈值，计算光子对被探测器接收的概率；根据光子对被探测器接收的概率，获取有效光子对的探测信息。
在其中一个实施例中，获取有效光子对的探测信息还包括：获取预设飞行时间重建标识；当预设飞行时间重建标识为进行飞行时间重建，则对有效光子对进行飞行时间重建；获取飞行时间重建后的有效光子对探测信息。其中，预设飞行时间重建标识为预先设定，当此次蒙特卡罗模拟需要进行飞行时间重建时，则预先将预设飞行时间重建标识设定为进行飞行时间重建；当此次蒙特卡罗模拟不需要进行飞行时间重建时，则预先将预设飞行时间重建标识设定为不进行飞行时间重建。
在其中一个实施例中，对有效光子对进行飞行时间重建包括：获取有效光子对中两个光子的光子达到探测器的输运时间；根据有效光子对中两个光子的光子达到探测器的输运时间进行时间分辨率高斯抽样；对时间分辨率高斯抽样的有效光子对，划分数据子集。更具体的，进行时间分辨率高斯抽样包括：获取相应探测器的时间分辨率；根据有效光子对中两个光子的光子达到探测器的输运时间，计算有效光子对中两个光子到达探测器的时间差；对时间差根据探测器的时间分辨率进行时间分辨率高斯抽样。首先根据光子对中两个光子的光子达到探测器的输运时间，计算两光子达到探测器的时间差；根据两光子分别抵达的探测器编号，更改时间差方向。探测器对光子探测时，记录光子的时间信息有一定误差，此误差分布满足高斯分布。因此，对模拟所得光子理想时间差进行随机抽样，即可得到探测器可能记录的光子到达探测器的时间差。根据时间分辨率高斯抽样的有效光子对的到达探测器的时间差，划分数据子集。例如，有效光子到达探测器的时间差范围为-1500ps到+1500ps，则将该范围分为15个子集，时间差为[-1300,-1100]ps的有效光子对记录为第1个子集；(-1300,-1100]ps的有效光子对记录为第2个子集，依次类推，共15分为15个子集。通过对有效光子对进行飞行时间重建能够提高重建图像的质量。
步骤S406：根据有效光子对的探测信息，得到散射弦图。
上述图像散射校正方法，对PET系统蒙特卡罗模拟的算法进行优化及模型简化，提高了PET系统蒙特卡罗校正算法的效率、降低了计算统计量。使得优化后的蒙特卡罗模拟算法能达到PET系统临床扫描的时间要求。并且通过对散射弦图进行去噪声处理，消除了散射弦图上的涨落噪声。由于散射弦图的形状分布为平滑过渡形，本方法在降低蒙特卡罗计算统计的同时，进一步的能够保证足够好的计算精度。
如图8所示，图8为一个实施例中本发明计算所得弦图与GATE仿真弦图对比图。此图为本发明计算结果(实线)与GATE仿真计算结果(虚线)的对比图。GATE仿真程序为开源蒙卡仿真程序，采用业内公认通用蒙卡仿真算法，具有极高精确度。本发明算法的计算速度比GATE仿真算法有数百倍的提升。如图所示，经过上述实施例中的算法计算得到的弦图与GATE仿真弦图吻合的很好。且尽管蒙特卡罗统计计算量的降低增加了弦图上的噪声，但平滑后的散射弦图与真实数据相符。
应该理解的是，虽然图1、2、3、7的流程图中的各个步骤按照箭头的指示依次显示，但是这些步骤并不是必然按照箭头指示的顺序依次执行。除非本文中有明确的说明，这些步骤的执行并没有严格的顺序限制，这些步骤可以以其它的顺序执行。而且，图1、2、3、7中的至少一部分步骤可以包括多个子步骤或者多个阶段，这些子步骤或者阶段并不必然是在同一时刻执行完成，而是可以在不同的时刻执行，这些子步骤或者阶段的执行顺序也不必然是依次进行，而是可以与其它步骤或者其它步骤的子步骤或者阶段的至少一部分轮流或者交替地执行。
在一个实施例中，如图9所示，提供了一种图像散射校正装置，包括：获取模块100、图像重建模块200、蒙特卡罗模拟模块300、去噪声处理模块400以及统计缩放处理模块500，其中：
获取模块100，用于获取图像原始数据，所述图像原始数据包括衰减图像以及测量弦图。
图像重建模块200，用于根据所述测量弦图获得放射图像；
蒙特卡罗模拟模块300，用于根据衰减图像以及放射图像进行蒙特卡罗模拟，得到散射弦图。
去噪声处理模块400，用于对散射弦图进行去噪声处理，得到低噪声散射弦图。
统计缩放处理模块500，用于根据图像原始数据对低噪声散射弦图进行统计缩放，得到校正散射弦图。
关于图像散射校正装置的具体限定可以参见上文中对于图像散射校正方法的限定，在此不再赘述。上述图像散射校正装置中的各个模块可全部或部分通过软件、硬件及其组合来实现。上述各模块可以硬件形式内嵌于或独立于计算机设备中的处理器中，也可以以软件形式存储于计算机设备中的存储器中，以便于处理器调用执行以上各个模块对应的操作。
在一个实施例中，提供了一种计算机设备，该计算机设备可以是终端，其内部结构图可以如图10所示。该计算机设备包括通过系统总线连接的处理器、存储器、网络接口、显示屏和输入装置。其中，该计算机设备的处理器用于提供计算和控制能力。该计算机设备的存储器包括非易失性存储介质、内存储器。该非易失性存储介质存储有操作系统和计算机程序。该内存储器为非易失性存储介质中的操作系统和计算机程序的运行提供环境。该计算机设备的网络接口用于与外部的终端通过网络连接通信。该计算机程序被处理器执行时以实现一种图像散射校正方法。该计算机设备的显示屏可以是液晶显示屏或者电子墨水显示屏，该计算机设备的输入装置可以是显示屏上覆盖的触摸层，也可以是计算机设备外壳上设置的按键、轨迹球或触控板，还可以是外接的键盘、触控板或鼠标等。
本领域技术人员可以理解，图10中示出的结构，仅仅是与本申请方案相关的部分结构的框图，并不构成对本申请方案所应用于其上的计算机设备的限定，具体的计算机设备可以包括比图中所示更多或更少的部件，或者组合某些部件，或者具有不同的部件布置。
在一个实施例中，提供了一种计算机设备，包括存储器和处理器，存储器中存储有计算机程序，该处理器执行计算机程序时实现以下步骤：
获取图像原始数据，所述图像原始数据包括衰减图像以及测量弦图；根据所述测量弦图获得放射图像；根据衰减图像以及放射图像进行蒙特卡罗模拟，得到散射弦图；对散射弦图进行去噪声处理，得到低噪声散射弦图；根据图像原始数据对低噪声散射弦图进行统计缩放，得到校正散射弦图。
在一个实施例中，处理器执行计算机程序时还实现以下步骤：
根据校正散射弦图以及测量弦图进行图像重建，得到校正放射图像；根据预设迭代次数，多次获取校正散射弦图并对所述校正放射图像进行迭代校正，得到最终放射图像。
在一个实施例中，处理器执行计算机程序时还实现以下步骤：
根据衰减图像以及放射图像确定蒙特卡罗模拟边界以及初始的光子信息，图像原始数据包括：衰减图像以及测量弦图，光子信息包括：光子位置、光子方向以及光子能量；根据蒙特卡罗模拟边界以及初始的光子信息进行光子运动模拟，得到光子到达探测器的探测信息；根据探测信息，得到散射弦图。
在一个实施例中，处理器执行计算机程序时还实现以下步骤：
根据衰减图像设定蒙特卡罗边界；根据衰减图像以及放射图像，确定初始的光子信息。
在一个实施例中，处理器执行计算机程序时还实现以下步骤：
对光子迁移量进行随机抽样；根据抽样的光子迁移量、初始的光子信息以及蒙特卡罗模拟边界进行光子运动模拟，得到光子到达探测器的探测信息。
在一个实施例中，处理器执行计算机程序时还实现以下步骤：
以预设步长移动相应光子；当光子移动距离未到达相应的光子迁移量，且光子移动位置达到蒙特卡罗模拟边界时，则将光子沿当前光子方向移动至探测器，根据光子信息得到光子到达探测器的探测信息，探测信息包括：散射标记、光子到达探测器的位置、光子到达探测器的方向、光子到达探测器的能量以及光子达到探测器的输运时间。
在一个实施例中，处理器执行计算机程序时还实现以下步骤：
当光子移动距离到达相应的光子迁移量，则根据康普散射角度抽样得到光子散射角度；根据光子散射角度，生成散射标记；根据光子散射角度以及光子迁移量更新光子信息；当更新后的光子信息所记录的光子能量大于第一预设阈值，则继续进行光子运动模拟，直到光子移动位置到达蒙特卡罗模拟边界。
在一个实施例中，处理器执行计算机程序时还实现以下步骤：
当光子移动距离未到达相应的光子迁移量，且光子移动位置未达到蒙特卡罗模拟边界时，则继续以预设步长移动相应光子，直至相应光子移动距离到达光子迁移量或相应光子移动位置达到蒙特卡罗模拟边界。
在一个实施例中，处理器执行计算机程序时还实现以下步骤：
获取至少一个光子对的探测信息；根据至少一个光子对的探测信息，得到有效光子对的探测信息；根据有效光子对的探测信息，得到散射弦图。
在一个实施例中，处理器执行计算机程序时还实现以下步骤：
获取相应探测器的能量分辨率；当光子对中的两个光子都到达探测器，则对光子到达探测器的能量根据能量分辨率进行高斯随机处理；当高斯随机处理后光子对中的两个光子到达探测器的能量都高于第二预设阈值，则将相应光子对作为有效光子对；获取有效光子对的探测信息。
在一个实施例中，处理器执行计算机程序时还实现以下步骤：
获取相应探测器的能量分辨率；当光子对中的两个光子都到达探测器，则对光子到达探测器的能量根据能量分辨率进行高斯随机处理；根据高斯随机处理后光子对中的两个光子到达探测器的能量、能量分辨率以及第二预设阈值，计算光子对被探测器接收的概率；根据光子对被探测器接收的概率，获取有效光子对的探测信息。
在一个实施例中，处理器执行计算机程序时还实现以下步骤：
获取预设飞行时间重建标识；当预设飞行时间重建标识为进行飞行时间重建，则对有效光子对进行飞行时间重建；获取飞行时间重建后的有效光子对探测信息。
在一个实施例中，处理器执行计算机程序时还实现以下步骤：
获取有效光子对中两个光子的光子达到探测器的输运时间；根据有效光子对中两个光子的光子达到探测器的输运时间进行时间分辨率高斯抽样；对时间分辨率高斯抽样的有效光子对，划分数据子集。
在一个实施例中，处理器执行计算机程序时还实现以下步骤：
获取相应探测器的时间分辨率；根据有效光子对中两个光子的光子达到探测器的输运时间，计算有效光子对中两个光子到达探测器的时间差；对时间差根据探测器的时间分辨率进行时间分辨率高斯抽样。
在一个实施例中，提供了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现以下步骤：
获取图像原始数据，所述图像原始数据包括衰减图像以及测量弦图；根据所述测量弦图获得放射图像；根据衰减图像以及放射图像进行蒙特卡罗模拟，得到散射弦图；对散射弦图进行去噪声处理，得到低噪声散射弦图；根据图像原始数据对低噪声散射弦图进行统计缩放，得到校正散射弦图。
在一个实施例中，计算机程序被处理器执行时还实现以下步骤：
根据校正散射弦图以及测量弦图进行图像重建，得到校正放射图像；根据预设迭代次数，多次获取校正散射弦图并对所述校正放射图像进行迭代校正，得到最终放射图像。
在一个实施例中，计算机程序被处理器执行时还实现以下步骤：
根据衰减图像以及放射图像确定蒙特卡罗模拟边界以及初始的光子信息，图像原始数据包括：衰减图像以及测量弦图，光子信息包括：光子位置、光子方向以及光子能量；根据蒙特卡罗模拟边界以及初始的光子信息进行光子运动模拟，得到光子到达探测器的探测信息；根据探测信息，得到散射弦图。
在一个实施例中，计算机程序被处理器执行时还实现以下步骤：
根据衰减图像设定蒙特卡罗边界；根据衰减图像以及放射图像，确定初始的光子信息。
在一个实施例中，计算机程序被处理器执行时还实现以下步骤：
对光子迁移量进行随机抽样；根据抽样的光子迁移量、初始的光子信息以及蒙特卡罗模拟边界进行光子运动模拟，得到光子到达探测器的探测信息。
在一个实施例中，计算机程序被处理器执行时还实现以下步骤：
以预设步长移动相应光子；当光子移动距离未到达相应的光子迁移量，且光子移动位置达到蒙特卡罗模拟边界时，则将光子沿当前光子方向移动至探测器，根据光子信息得到光子到达探测器的探测信息，探测信息包括：散射标记、光子到达探测器的位置、光子到达探测器的方向、光子到达探测器的能量以及光子达到探测器的输运时间。
在一个实施例中，计算机程序被处理器执行时还实现以下步骤：
当光子移动距离到达相应的光子迁移量，则根据康普散射角度抽样得到光子散射角度；根据光子散射角度，生成散射标记；根据光子散射角度以及光子迁移量更新光子信息；当更新后的光子信息所记录的光子能量大于第一预设阈值，则继续进行光子运动模拟，直到光子移动位置到达蒙特卡罗模拟边界。
在一个实施例中，计算机程序被处理器执行时还实现以下步骤：
当光子移动距离未到达相应的光子迁移量，且光子移动位置未达到蒙特卡罗模拟边界时，则继续以预设步长移动相应光子，直至相应光子移动距离到达光子迁移量或相应光子移动位置达到蒙特卡罗模拟边界。
在一个实施例中，计算机程序被处理器执行时还实现以下步骤：
获取至少一个光子对的探测信息；根据至少一个光子对的探测信息，得到有效光子对的探测信息；根据有效光子对的探测信息，得到散射弦图。
在一个实施例中，计算机程序被处理器执行时还实现以下步骤：
获取相应探测器的能量分辨率；当光子对中的两个光子都到达探测器，则对光子到达探测器的能量根据能量分辨率进行高斯随机处理；当高斯随机处理后光子对中的两个光子到达探测器的能量都高于第二预设阈值，则将相应光子对作为有效光子对；获取有效光子对的探测信息。
在一个实施例中，计算机程序被处理器执行时还实现以下步骤：
获取相应探测器的能量分辨率；当光子对中的两个光子都到达探测器，则对光子到达探测器的能量根据能量分辨率进行高斯随机处理；根据高斯随机处理后光子对中的两个光子到达探测器的能量、能量分辨率以及第二预设阈值，计算光子对被探测器接收的概率；根据光子对被探测器接收的概率，获取有效光子对的探测信息。
在一个实施例中，计算机程序被处理器执行时还实现以下步骤：
获取预设飞行时间重建标识；当预设飞行时间重建标识为进行飞行时间重建，则对有效光子对进行飞行时间重建；获取飞行时间重建后的有效光子对探测信息。
在一个实施例中，计算机程序被处理器执行时还实现以下步骤：
获取有效光子对中两个光子的光子达到探测器的输运时间；根据有效光子对中两个光子的光子达到探测器的输运时间进行时间分辨率高斯抽样；对时间分辨率高斯抽样的有效光子对，划分数据子集。
在一个实施例中，计算机程序被处理器执行时还实现以下步骤：
获取相应探测器的时间分辨率；根据有效光子对中两个光子的光子达到探测器的输运时间，计算有效光子对中两个光子到达探测器的时间差；对时间差根据探测器的时间分辨率进行时间分辨率高斯抽样。
本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，的计算机程序可存储于一非易失性计算机可读取存储介质中，该计算机程序在执行时，可包括如上述各方法的实施例的流程。其中，本申请所提供的各实施例中所使用的对存储器、存储、数据库或其它介质的任何引用，均可包括非易失性和/或易失性存储器。非易失性存储器可包括只读存储器(ROM)、可编程ROM(PROM)、电可编程ROM(EPROM)、电可擦除可编程ROM(EEPROM)或闪存。易失性存储器可包括随机存取存储器(RAM)或者外部高速缓冲存储器。作为说明而非局限，RAM以多种形式可得，诸如静态RAM(SRAM)、动态RAM(DRAM)、同步DRAM(SDRAM)、双数据率SDRAM(DDRSDRAM)、增强型SDRAM(ESDRAM)、同步链路(Synchlink)DRAM(SLDRAM)、存储器总线(Rambus)直接RAM(RDRAM)、直接存储器总线动态RAM(DRDRAM)、以及存储器总线动态RAM(RDRAM)等。
以上实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。
以上所述实施例仅表达了本申请的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本申请构思的前提下，还可以做出若干变形和改进，这些都属于本申请的保护范围。因此，本申请专利的保护范围应以所附权利要求为准。
Claims (17)
1.一种图像散射校正方法，其特征在于，所述方法包括：
获取图像原始数据，所述图像原始数据包括衰减图像以及测量弦图；
根据所述测量弦图获得放射图像；
根据所述衰减图像以及放射图像进行蒙特卡罗模拟，得到散射弦图；
对所述散射弦图进行去噪声处理，得到低噪声散射弦图；
根据图像原始数据对所述低噪声散射弦图进行统计缩放，得到校正散射弦图。
2.根据权利要求1所述的方法，其特征在于，根据图像原始数据对所述低噪声散射弦图进行统计缩放，得到校正散射弦图之后，还包括：
根据所述校正散射弦图以及测量弦图进行图像重建，得到校正放射图像；
根据预设迭代次数，多次获取校正散射弦图并对所述校正放射图像进行迭代校正，得到最终放射图像。
3.根据权利要求1所述的方法，其特征在于，所述根据所述衰减图像以及放射图像进行蒙特卡罗模拟，得到散射弦图包括：
根据所述衰减图像以及放射图像确定蒙特卡罗模拟边界以及初始的光子信息，所述图像原始数据包括：衰减图像以及测量弦图，光子信息包括：光子位置、光子方向以及光子能量；
根据所述蒙特卡罗模拟边界以及初始的光子信息进行光子运动模拟，得到光子到达探测器的探测信息；
根据所述探测信息，得到散射弦图。
4.根据权利要求3所述的方法，其特征在于，所述根据所述衰减图像以及放射图像确定蒙特卡罗模拟边界以及初始的光子信息包括：
根据所述衰减图像设定蒙特卡罗边界；
根据所述衰减图像以及放射图像，确定初始的光子信息。
5.根据权利要求3所述的方法，其特征在于，所述根据所述蒙特卡罗模拟边界以及初始的光子信息进行光子运动模拟，得到光子到达探测器的探测信息包括：
对光子迁移量进行随机抽样；
根据抽样的光子迁移量、初始的光子信息以及蒙特卡罗模拟边界进行光子运动模拟，得到光子到达探测器的探测信息。
6.根据权利要求5所述的方法，其特征在于，所述根据抽样的光子迁移量、初始的光子信息以及蒙特卡罗模拟边界进行光子运动模拟，得到光子到达探测器的探测信息包括：
以预设步长移动相应光子；
当光子移动距离未到达相应的光子迁移量，且光子移动位置达到蒙特卡罗模拟边界时，则将光子沿当前光子方向移动至探测器，根据光子信息得到光子到达探测器的探测信息，所述探测信息包括：散射标记、光子到达探测器的位置、光子到达探测器的方向、光子到达探测器的能量以及光子达到探测器的输运时间。
7.根据权利要求6所述的方法，其特征在于，所述以预设步长移动相应光子之后还包括：
当光子移动距离到达相应的光子迁移量，则根据康普散射角度抽样得到光子散射角度；
根据所述光子散射角度，生成散射标记；
根据光子散射角度以及光子迁移量更新光子信息；
当更新后的光子信息所记录的光子能量大于第一预设阈值，则继续进行光子运动模拟，直到光子移动位置到达蒙特卡罗模拟边界。
8.根据权利要求6所述的方法，其特征在于，所述以预设步长移动相应光子之后还包括：
当光子移动距离未到达相应的光子迁移量，且光子移动位置未达到蒙特卡罗模拟边界时，则继续以预设步长移动相应光子，直至相应光子移动距离到达光子迁移量或相应光子移动位置达到蒙特卡罗模拟边界。
9.根据权利要求7所述的方法，其特征在于，所述根据所述探测信息，得到散射弦图包括：
获取至少一个光子对的探测信息；
根据至少一个光子对的所述探测信息，得到有效光子对的探测信息；
根据有效光子对的探测信息，得到散射弦图。
10.根据权利要求9所述的方法，其特征在于，所述根据至少一个光子对的所述探测信息，得到有效光子对的探测信息包括：
获取相应探测器的能量分辨率；
当光子对中的两个光子都到达探测器，则对光子到达探测器的能量根据所述能量分辨率进行高斯随机处理；
当高斯随机处理后光子对中的两个光子到达探测器的能量都高于第二预设阈值，则将相应光子对作为有效光子对；
获取有效光子对的探测信息。
11.根据权利要求9所述的方法，其特征在于，所述根据至少一个光子对的所述探测信息，得到有效光子对的探测信息包括：
获取相应探测器的能量分辨率；
当光子对中的两个光子都到达探测器，则对光子到达探测器的能量根据所述能量分辨率进行高斯随机处理；
根据高斯随机处理后光子对中的两个光子到达探测器的能量、能量分辨率以及第二预设阈值，计算光子对被探测器接收的概率；
根据所述光子对被探测器接收的概率，获取有效光子对的探测信息。
12.根据权利要求10-11中任一项所述的方法，其特征在于，所述获取有效光子对的探测信息还包括：
获取预设飞行时间重建标识；
当所述预设飞行时间重建标识为进行飞行时间重建，则对有效光子对进行飞行时间重建；
获取飞行时间重建后的有效光子对探测信息。
13.根据权利要求12所述的方法，其特征在于，所述对有效光子对进行飞行时间重建包括：
获取有效光子对中两个光子的光子达到探测器的输运时间；
根据有效光子对中两个光子的光子达到探测器的输运时间进行时间分辨率高斯抽样；
对时间分辨率高斯抽样的有效光子对，划分数据子集。
14.根据权利要求13所述的方法，其特征在于，所述根据有效光子对中两个光子的光子达到探测器的输运时间进行时间分辨率高斯抽样包括：
获取相应探测器的时间分辨率；
根据有效光子对中两个光子的光子达到探测器的输运时间，计算有效光子对中两个光子到达探测器的时间差；
对所述时间差根据所述探测器的时间分辨率进行时间分辨率高斯抽样。
15.一种图像散射校正装置，其特征在于，所述装置包括：
获取模块，用于获取图像原始数据，所述图像原始数据包括衰减图像以及测量弦图；
图像重建模块，用于根据所述测量弦图获得放射图像；
蒙特卡罗模拟模块，用于根据所述衰减图像以及放射图像进行蒙特卡罗模拟，得到散射弦图；
去噪声处理模块，用于对所述散射弦图进行去噪声处理，得到低噪声散射弦图；
统计缩放处理模块，用于根据图像原始数据对所述低噪声散射弦图进行统计缩放，得到校正散射弦图。
16.一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，其特征在于，所述处理器执行所述计算机程序时实现权利要求1至14中任一项所述方法的步骤。
17.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现权利要求1至14中任一项所述的方法的步骤。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910281959.XA CN110009707B (zh) | 2019-04-09 | 2019-04-09 | 图像散射校正方法、装置、计算机设备和存储介质 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910281959.XA CN110009707B (zh) | 2019-04-09 | 2019-04-09 | 图像散射校正方法、装置、计算机设备和存储介质 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110009707A true CN110009707A (zh) | 2019-07-12 |
| CN110009707B CN110009707B (zh) | 2023-01-03 |
Family
ID=67170640
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910281959.XA Active CN110009707B (zh) | 2019-04-09 | 2019-04-09 | 图像散射校正方法、装置、计算机设备和存储介质 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110009707B (zh) |
Cited By (10)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110687585A (zh) * | 2019-09-23 | 2020-01-14 | 上海联影医疗科技有限公司 | 获取晶体效率的方法、装置、计算机设备和存储介质 |
| CN110706175A (zh) * | 2019-09-27 | 2020-01-17 | 上海联影医疗科技有限公司 | Pet校正系数的生成方法、系统、可读存储介质和设备 |
| CN111067560A (zh) * | 2019-12-25 | 2020-04-28 | 东软医疗系统股份有限公司 | 散射校正方法、装置、可读存储介质和电子设备 |
| CN111588399A (zh) * | 2020-05-27 | 2020-08-28 | 上海联影医疗科技有限公司 | 医学成像设备状态监控的方法、设备和计算机设备 |
| CN111627084A (zh) * | 2020-05-29 | 2020-09-04 | 上海联影医疗科技有限公司 | 基于pet扫描的图像校正方法、装置和系统 |
| CN112816510A (zh) * | 2021-03-03 | 2021-05-18 | 明峰医疗系统股份有限公司 | Ct扫描设备的探测信号处理方法、系统及计算机可读存储介质 |
| WO2021191151A1 (en) * | 2020-03-24 | 2021-09-30 | Aarhus Universitet | Device and method for pet image analysis and reconstruction |
| CN113804709A (zh) * | 2021-09-22 | 2021-12-17 | 清华大学 | 基于拟蒙特卡罗和强制探测的校正ct散射信号的方法 |
| CN115429299A (zh) * | 2022-09-26 | 2022-12-06 | 明峰医疗系统股份有限公司 | 一种基于定位片的散射校正的方法、系统、设备及存储介质 |
| CN116563413A (zh) * | 2023-07-04 | 2023-08-08 | 中国科学院深圳先进技术研究院 | 校正参数确定方法、图像重建方法、装置、设备及介质 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2013175352A1 (en) * | 2012-05-21 | 2013-11-28 | Koninklijke Philips N.V. | Fast scatter estimation in pet reconstruction. |
| CN108618796A (zh) * | 2018-02-09 | 2018-10-09 | 南方医科大学 | 一种基于任务驱动的蒙特卡洛散射光子模拟方法 |
-
2019
- 2019-04-09 CN CN201910281959.XA patent/CN110009707B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2013175352A1 (en) * | 2012-05-21 | 2013-11-28 | Koninklijke Philips N.V. | Fast scatter estimation in pet reconstruction. |
| CN104335247A (zh) * | 2012-05-21 | 2015-02-04 | 皇家飞利浦有限公司 | 在pet重建中的快速散射估计 |
| CN108618796A (zh) * | 2018-02-09 | 2018-10-09 | 南方医科大学 | 一种基于任务驱动的蒙特卡洛散射光子模拟方法 |
Cited By (17)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110687585A (zh) * | 2019-09-23 | 2020-01-14 | 上海联影医疗科技有限公司 | 获取晶体效率的方法、装置、计算机设备和存储介质 |
| CN110687585B (zh) * | 2019-09-23 | 2021-05-11 | 上海联影医疗科技股份有限公司 | 获取晶体效率的方法、装置、计算机设备和存储介质 |
| CN110706175B (zh) * | 2019-09-27 | 2022-07-29 | 上海联影医疗科技股份有限公司 | Pet校正系数的生成方法、系统、可读存储介质和设备 |
| CN110706175A (zh) * | 2019-09-27 | 2020-01-17 | 上海联影医疗科技有限公司 | Pet校正系数的生成方法、系统、可读存储介质和设备 |
| CN111067560A (zh) * | 2019-12-25 | 2020-04-28 | 东软医疗系统股份有限公司 | 散射校正方法、装置、可读存储介质和电子设备 |
| CN111067560B (zh) * | 2019-12-25 | 2023-05-02 | 沈阳智核医疗科技有限公司 | 散射校正方法、装置、可读存储介质和电子设备 |
| WO2021191151A1 (en) * | 2020-03-24 | 2021-09-30 | Aarhus Universitet | Device and method for pet image analysis and reconstruction |
| CN111588399A (zh) * | 2020-05-27 | 2020-08-28 | 上海联影医疗科技有限公司 | 医学成像设备状态监控的方法、设备和计算机设备 |
| CN111588399B (zh) * | 2020-05-27 | 2023-05-16 | 上海联影医疗科技股份有限公司 | 医学成像设备状态监控的方法、设备和计算机设备 |
| CN111627084A (zh) * | 2020-05-29 | 2020-09-04 | 上海联影医疗科技有限公司 | 基于pet扫描的图像校正方法、装置和系统 |
| CN111627084B (zh) * | 2020-05-29 | 2023-11-14 | 上海联影医疗科技股份有限公司 | 基于pet扫描的图像校正方法、装置和系统 |
| CN112816510B (zh) * | 2021-03-03 | 2022-06-21 | 明峰医疗系统股份有限公司 | Ct扫描设备的探测信号处理方法、系统及计算机可读存储介质 |
| CN112816510A (zh) * | 2021-03-03 | 2021-05-18 | 明峰医疗系统股份有限公司 | Ct扫描设备的探测信号处理方法、系统及计算机可读存储介质 |
| CN113804709A (zh) * | 2021-09-22 | 2021-12-17 | 清华大学 | 基于拟蒙特卡罗和强制探测的校正ct散射信号的方法 |
| CN115429299A (zh) * | 2022-09-26 | 2022-12-06 | 明峰医疗系统股份有限公司 | 一种基于定位片的散射校正的方法、系统、设备及存储介质 |
| CN116563413A (zh) * | 2023-07-04 | 2023-08-08 | 中国科学院深圳先进技术研究院 | 校正参数确定方法、图像重建方法、装置、设备及介质 |
| CN116563413B (zh) * | 2023-07-04 | 2024-04-09 | 中国科学院深圳先进技术研究院 | 校正参数确定方法、图像重建方法、装置、设备及介质 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110009707B (zh) | 2023-01-03 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110009707A (zh) | 图像散射校正方法、装置、计算机设备和存储介质 | |
| CN110097611B (zh) | 图像重建方法、装置、设备及存储介质 | |
| Iriarte et al. | System models for PET statistical iterative reconstruction: A review | |
| CN111080584B (zh) | 医学图像的质控方法、计算机设备和可读存储介质 | |
| CN108615250B (zh) | 图像重建方法、装置、系统和计算机可读存储介质 | |
| CN110584698B (zh) | 探测器质量控制效验方法、装置、计算机设备和存储介质 | |
| CN110025329B (zh) | 符合计数弦图生成方法、装置、计算机设备和存储介质 | |
| US10210635B2 (en) | Reconstruction quality assessment with local non-uniformity in nuclear imaging | |
| US11860320B2 (en) | Gamma camera dead time determination in real time using long lived radioisotopes | |
| CN111588399B (zh) | 医学成像设备状态监控的方法、设备和计算机设备 | |
| US10126439B2 (en) | Reconstruction with multiple photopeaks in quantitative single photon emission computed tomography | |
| CN110047116B (zh) | Pet图像校正方法、装置、计算机设备和存储介质 | |
| CN108209958A (zh) | 一种归一化校正因子的确定、获取方法及医学成像方法 | |
| CN110136076A (zh) | 医学扫描成像方法、装置、存储介质及计算机设备 | |
| CN110687585B (zh) | 获取晶体效率的方法、装置、计算机设备和存储介质 | |
| EP1631844B1 (en) | Generating detector efficiency estimates for a oet scanner | |
| CN110400361B (zh) | 子集划分及图像重建的方法、装置和计算机设备 | |
| CN110215223A (zh) | 散射校正方法、系统、可读存储介质和设备 | |
| Hebert et al. | The GEM MAP algorithm with 3-D SPECT system response | |
| CN108873047A (zh) | 检测放射源活度的方法、系统、计算机设备和存储介质 | |
| Xiao et al. | A new method for spatial structure detection of complex inner cavities based on 3D γ-photon imaging | |
| EP3104196B1 (en) | Gamma camera dead time compensation using a companion radioisotope | |
| Xiao et al. | A study on scattering correction for γ-photon 3D imaging test method | |
| CN107961028A (zh) | 一种归一化校正因子获取、确定方法及医学成像方法 | |
| EP3835828A1 (en) | Ai-based scintillator response modelling for increased spatial resolution in pet |
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 | ||
| CB02 | Change of applicant information |
Address after: 201807 Shanghai City, north of the city of Jiading District Road No. 2258 Applicant after: Shanghai Lianying Medical Technology Co.,Ltd. Address before: 201807 Shanghai City, north of the city of Jiading District Road No. 2258 Applicant before: SHANGHAI UNITED IMAGING HEALTHCARE Co.,Ltd. |
|
| CB02 | Change of applicant information | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |