CN113506355B - 散射校正方法、装置、成像系统及计算机可读存储介质 - Google Patents

散射校正方法、装置、成像系统及计算机可读存储介质 Download PDF

Info

Publication number
CN113506355B
CN113506355B CN202111060795.1A CN202111060795A CN113506355B CN 113506355 B CN113506355 B CN 113506355B CN 202111060795 A CN202111060795 A CN 202111060795A CN 113506355 B CN113506355 B CN 113506355B
Authority
CN
China
Prior art keywords
events
event data
coincidence
energy
scatter
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
CN202111060795.1A
Other languages
English (en)
Other versions
CN113506355A (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.)
Raycan Technology Co Ltd
Original Assignee
Raycan Technology Co Ltd
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 Raycan Technology Co Ltd filed Critical Raycan Technology Co Ltd
Priority to CN202111060795.1A priority Critical patent/CN113506355B/zh
Priority to PCT/CN2021/124039 priority patent/WO2023035360A1/zh
Priority to PCT/CN2021/124040 priority patent/WO2023035361A1/zh
Publication of CN113506355A publication Critical patent/CN113506355A/zh
Application granted granted Critical
Publication of CN113506355B publication Critical patent/CN113506355B/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

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Nuclear Medicine (AREA)

Abstract

本申请提出了一种散射校正方法、装置、成像系统及计算机可读存储介质，该方法包括利用校正参数构建探测器探测到的响应线上的符合事件内的初始单事件数据与校正单事件数据之间的对应关系；构建所述对应关系的目标函数；计算所述目标函数取最大值时所述校正单事件数据的取值，并且将所述校正单事件数据的取值确定为所述符合事件内的与初始单事件数据对应的未散射单事件数；以及利用符合事件数和所述未散射单事件数计算所述符合事件内的散射符合事件数。根据本申请的一些实施例，对所采集的符合事件进行校正时，并不需要活度图像和衰减图像，实现简单，准确性高。

Description

散射校正方法、装置、成像系统及计算机可读存储介质
技术领域
本申请涉及数据处理领域，具体而言，涉及一种散射校正方法、装置、成像系统及计算机可读存储介质。
背景技术
正电子发射断层成像（Positron Emission Tomography，简称PET）技术是当前全球尖端的分子影像技术之一，其通过对生物体内的标记有放射性核素的化合物进行成像，能够无创、定量、动态地评估生物体内各个功能器官的代谢水平、生化反应和功能活动，具有高灵敏度和准确性。
PET的工作原理为将发射正电子的放射性核素标记到能够参与活体组织血流或代谢过程的化合物上，并将标有放射性核素的化合物注射到生物体内，放射性核素在生物体内发射出的正电子移动大约0.1mm后与生物体内的负电子结合，从而导致发生电子对的湮灭事件，产生能量相等、方向相反的两个伽马光子。由于这两个伽马光子的飞行方向不同，所以探测器探测到这两个伽马光子的时间也不同。如果探测器中的位于响应线（Line ofResponse，简称LOR）上的两个闪烁晶体分别在规定的符合时间窗（例如，0~15纳秒）内探测到两个伽马光子，则探测到这两个伽马光子的事件可以称为符合事件。
符合事件一般可以包括真符合事件、散射符合事件和随机符合事件。其中，真符合事件是指同一个湮灭事件产生的两个伽马光子到达探测器中的位于响应线上的两个闪烁晶体的时间差在符合时间窗内的事件。随机符合事件是一种假符合事件，在随机符合事件中，探测到的两个伽马光子来自不同的湮灭事件，但在符合时间窗内被误认为是“同时”发生的2个伽马光子。散射符合事件是指以下事件：针对探测到的同一个湮灭事件产生的2个伽马光子，其中一个伽马光子由于在飞行过程中发生康普顿散射和/或瑞利散射等物理效应而改变飞行方向。
在这三种符合事件中，由于针对随机符合事件和散射符合事件所采集到的数据可能是错误的，这可能会影响PET成像的分辨率、对比度以及定位精度。因此，对所采集的符合事件进行校正显得至关重要。
目前，较为常用的散射校正方法主要包括多能窗技术、卷积/反卷积技术、以及基于模拟的技术等。其中，在这些技术中，最准确、最为广泛应用的是基于模拟的技术。
基于模拟的技术一般包括单散射模拟方法、双散射模拟方法以及蒙特卡洛模拟方法。其中，单散射模拟方法是针对每条响应线，通过在所输入的活度图像和衰减图像中选取一个散射点来获得单散射（即，一对伽马光子总共发生一次散射）的光子运动路径，通过计算该条响应线所对应的所有光子运动路径上产生的单散射事件来得到该条响应线上的单散射事件，最后，利用仅包含散射事件的尾部数据，通过尾部拟合技术（Tail Fitting，简称TF）对所得到的单散射事件进行拟合以得到散射符合事件的拉伸因子，并将该拉伸因子作用于所有数据以实现散射校正。双散射模拟方法通常与单散射模拟方法结合使用，其与单散射模拟方法的区别主要在于双散射模拟方法是通过两个散射点来确定光子运动路径。蒙特卡洛模拟方法主要通过对基于所输入的活度图像和衰减图像确定的每个伽马光子对的运动进行模拟来确定总散射分布。
在实现本申请过程中，发明人发现现有技术至少存在如下问题：
（1）单散射模拟方法仅考虑单散射的情况，然而在实际中会出现多散射的情况，这可能会导致散射校正结果的准确性较低，从而导致重建图像的质量较低。双散射模拟方法和蒙特卡洛模拟方法虽然考虑多散射的情况，但计算过程复杂并且计算量很大。
（2）现有技术很难获得准确的活度图像和衰减图像。
（3）由于视野范围之外的活度图像很难获得，因此，现有技术通常没有考虑外部辐射（即，位于视野范围之外的伽马光子对产生的散射事件），这可能会导致校正结果的准确性较低，从而影响重建图像的质量。
发明内容
本申请提供了一种散射校正方法、装置、成像系统及计算机可读存储介质，以解决现有技术中存在的至少一个问题。
根据本申请的一方面，提出了一种散射校正方法，包括利用校正参数构建探测器探测到的响应线上的符合事件内的初始单事件数据与校正单事件数据之间的对应关系，其中，所述初始单事件数据包括在所述探测器探测到的伽马光子的各个实测能量满足对应的预设条件下的单事件数；构建所述对应关系的目标函数，计算所述目标函数取最大值时所述校正单事件数据的取值，并且将所述校正单事件数据的取值确定为所述符合事件内的与所述初始单事件数据对应的未散射单事件数；以及利用符合事件数和所述未散射单事件数计算所述符合事件内的散射符合事件数。
可选地，所述初始单事件数据通过以下方式来确定：
将所述符合事件拆分成单事件，统计每个所述单事件所对应的所述实测能量，并统计同一所述实测能量所对应的单事件数以获得所述初始单事件数据；
将所述符合事件拆分成单事件，根据所述单事件所在的探头位置绘制至少一个能谱，并且获取各个所述能谱中的所述初始单事件数据；或者
将所述符合事件拆分成单事件，根据所述单事件所在的探头位置绘制至少一个能谱，根据伽马光子到达所述探头的时间差将每个所述至少一个能谱划分成多个子能谱，并且获取各个所述子能谱中的所述初始单事件数据。
可选地，针对所获取的每个所述初始单事件数据，按照公式（1）构建所述对应关系：
Figure 478183DEST_PATH_IMAGE001
公式（1），
其中，
Figure 488864DEST_PATH_IMAGE002
，
其中，X为校正单事件数据，并且，Xj为在伽马光子的真实能量Et满足E j≤Et＜Ej+1下的单事件数；Y为初始单事件数据，并且
Figure 456820DEST_PATH_IMAGE003
，Yi为在实测能量Ed满足 E i≤Ed＜Ei+1下的单事件数；
Figure 615269DEST_PATH_IMAGE004
为Y的期望值，并且
Figure 326873DEST_PATH_IMAGE005
；P为校正参数，P中的元素Pij为真实能量Et满足 E j≤Et＜Ej+1的伽马光子在实测后的实测能量Ed同时满足 E i≤Ed＜Ei+1的概率；Ej和Ej+1分别为预设能窗{Ej}中的第j个能量和第j+1个能量；Ei和Ei+1分别为实测能窗{Ei}中的第i个能量和第i+1个能量；m、n、i和j均为正整数。
可选地，当m=n且{Ei}={Ej}时，通过公式（2）计算所述校正参数P的第i行第j列的元素：
Figure 457640DEST_PATH_IMAGE006
公式（2），
其中，
Figure 596497DEST_PATH_IMAGE007
为在忽略散射符合事件下对真实能量为Et的伽马光子进行扫描所得到的初始单事件数据；所述真实能量Et满足 E i≤Et＜Ei+1；i、j、s、m为正整数；
将计算步骤
Figure 914346DEST_PATH_IMAGE008
和矩阵运算φ=Pθ替换为φ=convolution(P i , θ)；将计算步骤
Figure 429641DEST_PATH_IMAGE009
和矩阵运算θ=P T φ替换为θ=convolution(P i ,φ)，其中，
Figure 227964DEST_PATH_IMAGE010
，m=n，θ={θ j }，j=1,2，…，n，φ={φ i }，i=1,2，…，m，θ、φ被替换为任意具有相同维度的向量，包括校正单事件数据X或初始单事件数据Y。
可选地，所述伽马光子的真实能量Et包括包括202keV、307keV、511keV、622keV或1.27MeV。
可选地，所述构建所述对应关系的目标函数的步骤包括：
利用极大似然法构建所述目标函数，所述目标函数通过公式（3）确定：
Figure 537723DEST_PATH_IMAGE011
公式（3），
其中，β为超参数；R(X)为惩罚项。
可选地，所述惩罚项R(X)按照公式（4）进行计算：
Figure 342868DEST_PATH_IMAGE012
公式（4），
其中，M为常数。
可选地，计算所述目标函数取最大值时所述校正单事件数据的取值的步骤包括通过公式（5）或公式（6）计算所述校正单事件数据的取值：
Figure 396274DEST_PATH_IMAGE013
公式（5），
或
Figure 236054DEST_PATH_IMAGE014
公式（6）
其中，
Figure 779031DEST_PATH_IMAGE015
，其中，
Figure 337051DEST_PATH_IMAGE016
和
Figure 928570DEST_PATH_IMAGE017
为校正单事件数据的取值；
Figure 888436DEST_PATH_IMAGE018
为第k次迭代的中间变量；k为迭代次数。
可选地，利用公式（7）-公式（12）计算所述符合事件内的所述散射符合事件数：
SC≈SC+ 公式（7）
SC≈SC- 公式（8）
SC≈1/[(1/ SC+)+(1/ SC-)] 公式（9）
Figure 539997DEST_PATH_IMAGE019
公式（10）
SC≈0.5(SC++SC-) 公式（11）
或
Figure 319734DEST_PATH_IMAGE020
公式（12）
其中，
Figure 530922DEST_PATH_IMAGE021
，
Figure 345294DEST_PATH_IMAGE022
；S tot 为总单事件数，并且为符合事件数的2倍；S unsc 为未散射单事件数的总和；SC为散射符合事件数。
根据本申请的另一个方面，还提供了一种散射校正方法，包括：
对探测器探测到的多条响应线进行降采样以获得降采样响应线；
利用校正参数构建所述降采样响应线上的符合事件内的初始单事件数据与校正单事件数据之间的对应关系，其中，所述初始单事件数据包括在所述探测器探测到的伽马光子的各个实测能量在对应的预设条件下的单事件数；
构建所述对应关系的目标函数，计算所述目标函数取最大值时所述校正单事件数据的取值，并且将所述校正单事件数据的取值确定为所述符合事件内的与所述初始单事件数据对应的未散射单事件数；以及
利用降采样后的符合事件数和所述未散射单事件数计算所述符合事件内的散射符合事件数。
可选地，所述方法还包括：
按照预设方法对所述降采样响应线进行上采样以获得上采样响应线；
利用所述降采样响应线上的所述散射符合事件数计算所述上采样响应线上的散射符合事件数。
可选地，所述初始单事件数据通过以下方式来确定：
将所述符合事件拆分成单事件，统计每个所述单事件所对应的所述实测能量，并统计同一所述实测能量所对应的单事件数以获得所述初始单事件数据；
将所述符合事件拆分成单事件，根据所述单事件所在的探头位置绘制至少一个能谱，并且获取各个所述能谱中的所述初始单事件数据；或者
将所述符合事件拆分成单事件，根据所述单事件所在的探头位置绘制至少一个能谱，根据伽马光子到达所述探头的时间差将每个所述至少一个能谱划分成多个子能谱，并且获取各个所述子能谱中的所述初始单事件数据。
可选地，针对所获取的每个所述初始单事件数据，按照公式（13）构建所述对应关系：
Figure 167757DEST_PATH_IMAGE001
公式（13），
其中，
Figure 434790DEST_PATH_IMAGE002
，
其中，X为校正单事件数据，并且X={Xj}，j=1, 2，…，n，Xj为在伽马光子的真实能量Et满足 E j≤Et＜Ej+1下的单事件数；Y为初始单事件数据，并且 Y={YI}，i=1, 2，…，m，Yi为在实测能量Ed满足 E i≤Ed＜Ei+1下的单事件数；
Figure 633690DEST_PATH_IMAGE004
为Y的期望值，并且
Figure 364886DEST_PATH_IMAGE023
；P为校正参数，P中的元素Pij为真实能量Et满足 E j≤Et＜Ej+1的伽马光子在实测后的实测能量Ed同时满足 E i≤Ed＜Ei+1的概率；Ej和Ej+1分别为预设能窗{Ej}中的第j个能量和第j+1个能量；Ei和Ei+1分别为实测能窗{Ei}中的第i个能量和第i+1个能量；m、n、i和j均为正整数。
可选地，当m=n且{E i }={E j }时，通过公式（14）计算所述校正参数P的第i行第j列的元素：
Figure 358250DEST_PATH_IMAGE006
公式（14），
其中，
Figure 112579DEST_PATH_IMAGE007
为在忽略散射符合事件下对真实能量为Et的伽马光子进行扫描所得到的初始单事件数据；所述真实能量Et满足 E i≤Et＜Ei+1；i、j、s、m为正整数；
将计算步骤
Figure 115170DEST_PATH_IMAGE008
和矩阵运算φ=Pθ替换为φ=convolution(P i , θ)；将计算步骤
Figure 638555DEST_PATH_IMAGE009
和矩阵运算θ=P T φ替换为θ=convolution(P i ,φ)，其中，
Figure 615870DEST_PATH_IMAGE010
，m=n，θ={θ j }，j=1,2，…，n，φ={φ i }，i=1,2，…，m，θ、φ被替换为任意具有相同维度的向量，包括校正单事件数据X或初始单事件数据Y。
可选地，所述伽马光子的真实能量Et包括202keV、307keV、511keV、622keV或1.27MeV。
可选地，所述构建所述对应关系的目标函数的步骤包括：
利用极大似然法构建所述对应关系的目标函数，所述目标函数通过公式（15）确定：
Figure 857495DEST_PATH_IMAGE011
公式（15），
其中，β为超参数；R(X)为惩罚项。
可选地，所述惩罚项R(X)按照公式（16）进行计算：
Figure 398198DEST_PATH_IMAGE012
公式（16），
其中，M为常数。
可选地，所述计算所述目标函数取最大值时所述校正单事件数据的取值的步骤包括通过公式（17）或者公式（18）计算所述校正单事件数据的取值：
Figure 776090DEST_PATH_IMAGE013
公式（17），
或
Figure 376835DEST_PATH_IMAGE024
公式（18）
其中，
Figure 840178DEST_PATH_IMAGE015
，
Figure 246888DEST_PATH_IMAGE016
和
Figure 479287DEST_PATH_IMAGE017
均为校正单事件数据的取值；
Figure 250934DEST_PATH_IMAGE018
为第k次迭代的中间变量；k为迭代次数。
可选地，利用公式（7）-公式（12）计算所述符合事件内的所述散射符合事件数。
可选地，所述预设方法包括：双线性插值方法、4D线性插值方法或5D线性插值方法。
可选地，当所述预设方法为4D线性插值法时，按照公式（19）、公式（20）来计算所述上采样响应线上的所述散射符合事件数：
Figure 201572DEST_PATH_IMAGE025
公式（19）
Figure 959444DEST_PATH_IMAGE026
公式（20）
其中，表示上采样响应线上的散射符合事件数；SC K,L 表示降采样响应线上的散射符合事件数；I和J表示构成上采样响应线的晶体的编号；K和L表示构成降采样响应线的晶体的编号；W tot 为归一化因子；N为一条降采样响应线中包含的上采样响应线数；TI和TJ分别表示上采样晶体I和J对应的降采样中心晶体、降采样轴向相邻晶体、降采样径向相邻晶体以及降采样相对晶体的集合；
Figure 311927DEST_PATH_IMAGE027
表示对于上采样晶体I而言降采样晶体K的权重；
Figure 254476DEST_PATH_IMAGE028
表示对于上采样晶体J而言降采样晶体L的权重。
可选地，所述利用所述降采样响应线上的所述散射符合事件数计算所述上采样响应线上的散射符合事件数的步骤包括：
计算每条所述降采样响应线上的所述散射符合事件数与符合事件数的比例；
根据每条所述降采样响应线所包含的所述上采样响应线上的符合事件数和所述比例计算所述上采样响应线上的所述散射符合事件数。
可选地，所述方法还包括：
从所述降采样响应线所包含的上采样响应线上的符合事件数中剔除所述散射符合事件数。
可选地，所述方法还包括：
在对所述初始单事件数据进行校正之前，对所述符合事件进行随机校正。
根据本申请的另一个方面，还提供了一种散射校正装置，包括：
构建单元，其被配置为利用校正参数构建探测器探测到的响应线上的符合事件内的初始单事件数据与校正单事件数据之间的对应关系，其中，所述初始单事件数据包括在所述探测器探测到的伽马光子的各个实测能量满足对应的预设条件下的单事件数；
处理单元，其被配置为构建所述对应关系的目标函数，计算所述目标函数取最大值时所述校正单事件数据的取值，并且将所述校正单事件数据的取值确定为所述符合事件内的与所述初始单事件数据对应的未散射单事件数；以及
第一计算单元，其被配置为利用符合事件数和所述未散射单事件数计算所述符合事件内的散射符合事件数。
根据本申请的另一个方面，还提供了一种散射校正装置，包括：
存储器，其上存储有计算机程序；
处理器，其与所述存储器联接，并且当所述计算机程序被所述处理器执行时，实现上述方法。
根据本申请的另一个方面，还提供了一种成像系统，包括：
上述散射校正装置；
探测器，与所述散射校正装置相连接。
可选地，所述探测器包括PET探测器、PET-CT探测器、CT探测器、或者MR探测器。
根据本申请的另一个方面，还提供了一种计算机可读存储介质，其上存储有计算机程序指令，所述程序指令被执行时实现上述方法。
根据本申请的一些实施例，通过利用校正参数构建响应线上的符合事件内的初始单事件数据与校正单事件数据之间的对应关系，构建对应关系的目标函数，计算目标函数取最大值时所述校正单事件数据的取值，并且将校正单事件数据的取值确定为符合事件内的与初始单事件数据对应的未散射单事件数，并且利用未散射单事件数计算符合事件内的散射符合事件数，从而实现对符合事件的散射校正，而并不需要活度图像和衰减图像，使得符合事件的散射校正更易于实现。
根据本申请的一些实施例，考虑了多种散射情况，例如多散射和外部散射，使得散射校正结果的准确性较高，从而可以提高后续重建图像的质量。和现有技术中的双散射模拟方法和蒙特卡洛模拟方法相比，根据本申请的散射校正方法具有计算过程简单并且计算量较小。
附图说明
为了更清楚地说明本申请实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍。
图1示出了根据本申请实施例的一种散射校正方法的流程图。
图2示出了根据本申请实施例绘制的能谱示意图。
图3示出了根据本申请实施例的另一种散射校正方法的流程图。
图4示出了根据本申请实施例的一种合并多条响应线的示意图。
图5示出了根据本申请实施例的另一种合并多条响应线的示意图。
图6示出了根据本申请实施例的另一种散射校正方法的流程图。
图7示出了根据本申请实施例的又一种散射校正方法的流程图。
图8示出了根据本申请实施例的又一种散射校正方法的流程图。
图9示出了一个包含三个冷区的圆柱形假体示意图。
图10示出了针对图9中的圆柱形假体利用OSL-EDR方法计算出的未散射单事件数的误差图。
图11示出了针对图9中的圆柱形假体利用OT-EDR方法计算出的未散射单事件数的误差图。
图12示出了针对图9中的圆柱形假体计算出降采样后的散射符合事件数的误差图。
图13示出了根据本申请实施例的一种散射校正装置的框图。
图14示出了根据本申请实施例的另一种散射校正装置的框图。
图15示出了根据本申请实施例的另一种散射校正装置的框图。
图16示出了根据本申请实施例的又一种散射校正装置的框图。
图17示出了根据本申请实施例的一种成像系统的框图。
具体实施方式
现在将参考附图更全面地描述示例实施例。然而，示例实施例能够以多种形式实施，且不应被理解为限于在此阐述的实施例；相反，提供这些实施例使得本申请将全面和完整，并将示例实施例的构思全面地传达给本领域的技术人员。在图中相同的附图标记表示相同或类似的部分，因而将省略对它们的重复描述。
所描述的特征、结构或特性可以以任何合适的方式结合在一个或更多实施例中。在下面的描述中，提供许多具体细节从而给出对本公开的实施例的充分理解。然而，本领域技术人员将意识到，可以实践本公开的技术方案而没有这些特定细节中的一个或更多，或者可以采用其它的方式、组元、材料、装置或操作等。在这些情况下，将不详细示出或描述公知结构、方法、装置、实现、材料或者操作。
附图中所示的流程图仅是示例性说明，不是必须包括所有的内容和操作/步骤，也不是必须按所描述的顺序执行。例如，有的操作/步骤还可以分解，而有的操作/步骤可以合并或部分合并，因此实际执行的顺序有可能根据实际情况改变。
本申请的说明书和权利要求书及上述附图中的术语“第一”、“第二”等是用于区别不同对象，而不是用于描述特定顺序。此外，术语“包括”和“具有”以及它们任何变形，意图在于覆盖不排他的包含。例如包含了一系列步骤或单元的过程、方法、系统、产品或设备没有限定于已列出的步骤或单元，而是可选地还包括没有列出的步骤或单元，或可选地还包括对于这些过程、方法、产品或设备固有的其他步骤或单元。
需要说明的是，除了另有说明，否则本文中不带下标的大写字母可以表示向量或 矩阵，而带下标的大写字母可以表示向量或矩阵中的元素，并且大、小写字母可以表示不同 含义。另外，本文所述的求和符号默认表示对整个取值范围求和，如Σ可以表示为
Figure 505460DEST_PATH_IMAGE029
。
下面将参照附图，对根据本申请的具体实施例进行详细说明。
本申请提供的散射校正方法可以适用于各种成像系统，例如，PET成像系统、电子计算机断层扫描（Computed Tomography，简称CT）成像系统、PET-CT成像系统、磁共振（Magnetic Resonance，简称MR）成像系统以及各种探测器采集数据后进行图像重建的应用环境中。探测器可以包括多个探头，探测到符合事件的两个探头可构成一个探头对，并且每个探头对上可以形成一条或多条响应线。
图1示出了根据本申请实施例的一种散射校正方法的流程图。下面参照图1，对该散射校正方法进行详细说明。
在步骤S101，利用校正参数构建探测器探测到的响应线上的符合事件内的初始单事件数据与校正单事件数据之间的对应关系。
校正参数可以是指于校正探测器对光子能量（即，伽马光子的真实能量）产生的模糊效应的参数。其中，模糊效应通常是指针对一个具有特定真实能量（例如，511keV）的伽马光子而实际探测到多个实测能量的现象。例如，针对真实能量为E的伽马光子，但实际探测到的光子能量不一定为E，而是可能会出现E1、E2、…、Em+1等多个实测能量。
初始单事件数据可以包括在探测器探测到的伽马光子的各个实测能量满足对应的预设条件下的单事件数。例如，初始单事件数据可以包括在实测能量Et满足E i≤Et＜Ei+1下的单事件数Yi，i为正整数。当然预设条件也可以是满足实际需求的其它条件，在此并不限制。初始单事件数据所包括的单事件数与实测能量之间的关系可以用离散或连续函数等来表示，或者用统计表的形式来表示。而且，可以通过以下方式之一来获得初始单事件数据：
根据一个实施例，将每条响应线上的所有符合事件都拆分成单事件，然后统计每个单事件所对应的伽马光子的实测能量，并且统计同一实测能量所对应的单事件数，从而获得初始单事件数据。
根据另一个实施例，将每条响应线上的所有符合事件都拆分成单事件，按照单事件所在的探头位置绘制一个或多个能谱，如图2所示，所绘制的能谱指示了该探头上的单事件数与实测能量之间的对应关系，从而可以从各个能谱中获取对应的初始单事件数据。针对存在多个能谱的情况，总共可以获取多个初始单事件数据。另外，针对每个探头分别绘制对应的能谱，并且能谱的数量与每条响应线所覆盖的探头数相同。绘制能谱的具体过程可以参照现有技术，在此不再赘叙。
根据另一个实施例，将每条响应线上的所有符合事件都拆分成单事件，按照单事件所在的探头位置绘制一个或多个能谱。然后，根据伽马光子到达探头的时间差将至少一个能谱划分为多个子能谱，并且从各个子能谱中获取初始单事件数据，从而可以获得多个初始单事件数据。例如，针对飞行时间（Time of Flight，简称TOF）散射，在获得能谱后，根据伽马光子到达探头的时间差将每个能谱划分成多个子能谱，每个能谱所对应的子能谱的数量与时间差的段数相同。每一个子能谱均可以指示某个时间间隔内的单事件数与光子能量之间的对应关系。例如，以时间差200ps为例，当时间段数为10时，则每个能谱可以划分成10个子能谱，并且每个子能谱所对应的时间间隔为200ps。这种通过伽马光子到达探头的时间差将每个能谱划分成多个子能谱来获取每条响应线上的符合事件内的初始单事件数据的方式，可以使得后续计算出的散射符合事件数更加准确。
在一个实施例中，针对所获取的每个初始单事件数据，都可以按照公式（1）来构建初始单事件数据与校正单事件数据之间的对应关系：
Figure 191656DEST_PATH_IMAGE001
公式（1）
其中，
Figure 398646DEST_PATH_IMAGE002
，
其中，X为校正单事件数据，并且X={Xj}，j=1, 2，…，n，Xj为在伽马光子的真实能量Et满足 E j≤Et＜Ej+1下的单事件数；Y为初始单事件数据，并且 Y={YI}，i=1, 2，…，m，Yi为在实测能量Ed满足 E i≤Ed＜Ei+1下的单事件数；
Figure 512096DEST_PATH_IMAGE030
为Y的期望值，并且
Figure 437326DEST_PATH_IMAGE023
；P为校正参数，P中的元素Pij为真实能量Et满足 E j≤Et＜Ej+1的伽马光子在实测后的实测能量Ed同时满足 E i≤Ed＜Ei+1的概率；Ej和Ej+1分别为预设能窗{Ej}中的第j个能量和第j+1个能量；Ei和Ei+1分别为实测能窗{Ei}中的第i个能量和第i+1个能量；m、n、i和j均为正整数，并且i小于或等于m，j小于或等于n。
除了上述的矩阵之外，校正参数也可以采用其他形式，在此并没有限制。另外，上述校正参数可以通过利用探测器对具有不同真实能量的多个伽马光子进行探测以获得这些伽马光子的能量分布，并且对该能量分布进行归一化来获得。
在另一个实施例中，当m=n且{E i }={E j }时，可以通过公式（2）计算校正参数P的第i行第j列的元素：
Figure 723951DEST_PATH_IMAGE006
公式（2）
其中，
Figure 51027DEST_PATH_IMAGE007
为在忽略散射符合事件下对真实能量为Et的伽马光子进行扫描所得到的初始单事件数据；所述真实能量Et满足 E j≤Et＜Ej+1，Ej和Ej+1分别为预设能窗{Ej}中的第j个能量和第j+1个能量，通常n=m，{Ei}={Ej}，并且这两个预设能窗可以是根据实际需求来人为设置的，E1为预设能窗的低阈值，En、Em为预设能窗的高阈值；i、j、s、m为正整数。而且，伽马光子的真实能量Et可以为202keV、307keV、511keV、622keV或1.27MeV等。
例如，可以仅对E=511 keV的伽马光子进行扫描得到整个矩阵P。首先，可以根据公式（2）获得矩阵P在第i行的数值，从而可以根据第i行的数值得到整个矩阵，所得到的矩阵P的形式如下：
Figure 69799DEST_PATH_IMAGE031
公式（A）
从上式（A）可见，对于矩阵P在第i-a行或第i+a行的数值就是将矩阵P在第i行的数值左移或右移了a个单位，a为正整数。此时，若记
Figure 482326DEST_PATH_IMAGE010
，即矩阵P在第i行的数值逆序构成的向量。
通过在步骤S101利用校正参数构建探测器探测到的每条响应线上的符合事件内的初始单事件数据与校正单事件数据之间的对应关系，可以在后续计算过程中消除探测器对光子能量的模糊效应的影响。
在步骤S103，构建对应关系的目标函数，计算目标函数取最大值时校正单事件数据的取值，并且将校正单事件数据的取值确定为符合事件内的与初始单事件数据对应的未散射单事件数。
针对每个初始单事件数据，在构建其与校正单事件数据之间的对应关系之后，可以构建该对应关系的目标函数，计算目标函数的最大值，计算该目标函数取最大值时校正单事件数据的取值，并且将校正单事件数据的取值确定为计算符合事件内的与该初始单事件数据对应的未散射单事件数。其中，“与初始单事件数据对应的未散射单事件数”可以是指：针对一个初始单事件数据，可以计算出一个未散射单事件数。
根据一个实施例，针对所构建的每个对应关系，采用带惩罚的极大似然法构建如公式（3）所示的目标函数。但本申请不限于此，任一能够进行参数估计的方法均适用构建目标函数。
Figure 510325DEST_PATH_IMAGE011
公式（3）
其中，
Figure 502027DEST_PATH_IMAGE032
是超参数，其具体数值可提前通过实验来确定；R(X)为惩罚项。
接着，可以对公式（3）所示的目标函数进行求解以获得其最大值，计算此时校正单事件数据的取值，并且将该取值确定为该条响应线上的符合事件内的与该初始单事件数据对应的未散射单事件数。关于计算目标函数的最大值的过程可以参照现有技术中的对应描述，在此不再赘叙。
根据一个实施例，为了使所得到的未散射单事件数更加准确，在对上述目标函数进行求解之前，可以根据先验知识，即光子计数在小于放射源的能量处较少，在大于放射源的能量处为0，几乎都集中在放射源的能量处，将惩罚项R(X)定义为如公式（4）所示。
Figure 957279DEST_PATH_IMAGE012
公式（4）
对R(X)求偏导，得到公式（B）。
Figure 591523DEST_PATH_IMAGE033
公式（B），
其中，M为常数，一般大于100。
另外，可以利用以下两种方法来计算目标函数取最大值时校正单事件数据的取值。
根据一个实施例，利用基于迟一步法（One Step Later，简称OSL）的原理推导出的算法OSL-EDR（Eliminate Detector Response，简称EDR），如公式（5）所示，进行迭代计算公式（3）所示的目标函数取最大值时的校正单事件数据的取值，也即未散射单事件数。
Figure 157633DEST_PATH_IMAGE013
公式（5），
其中，
Figure 459302DEST_PATH_IMAGE016
和
Figure 85455DEST_PATH_IMAGE017
分别表示第k次和第k+1次迭代得到的在伽马光子的真实能量Et满足 E j≤Et＜Ej+1下的单事件数，即，所要计算的校正单事件数据的取值；k为迭代次数，其取值范围一般为200~500；
Figure 269312DEST_PATH_IMAGE034
可根据公式（D）来计算；b为1~n之间的正整数。
根据另一个实施例，利用基于优化迁移（Optimization Transfer，简称OT）的原理推导出的算法OT-EDR，如公式（6）所示，进行迭代计算公式（3）所示的目标函数取最大值时的校正单事件数据的取值，也即未散射单事件数。
Figure 639113DEST_PATH_IMAGE024
公式（6）
其中，
Figure 795288DEST_PATH_IMAGE015
，其为第k次迭代的中间变量。
在步骤S105，利用符合事件数和未散射单事件数计算符合事件内的散射符合事件数。
在计算出一条或多条响应线上的符合事件内的所有未散射单事件数之后，可以利用该条响应线上的所有未散射单事件数和符合事件数（即，符合事件的数量）计算每条响应线上的符合事件内的散射符合事件数。
公式（7）~公式（12）给出了六种计算每条响应线上的散射符合事件数的方法。
SC≈SC+ 公式（7）
SC≈SC- 公式（8）
SC≈1/[(1/ SC+)+(1/ SC-)] 公式（9）
Figure 592343DEST_PATH_IMAGE019
公式（10）
SC≈0.5(SC++SC-) 公式（11）
或
Figure 201179DEST_PATH_IMAGE020
公式（12）
其中，SC为散射符合事件数；
Figure 922141DEST_PATH_IMAGE021
，
Figure 932822DEST_PATH_IMAGE022
；S tot 为总单事件数，并且为符合事件数的2倍，如公式（C）所示；S unsc 为未散射单事件数的总和，如公式（D）所示。
Figure 900778DEST_PATH_IMAGE035
公式（C）
Figure 996910DEST_PATH_IMAGE036
公式（D）
其中，在公式（C）和公式（D）中，S0表示湮灭产生的两个伽马光子均未发生散射的符合事件数；S1表示湮灭产生的两个伽马光子中只有一个伽马光子发生了一次或多次散射的符合事件数；S2表示湮灭产生的两个伽马光子都发生了一次或多次散射的符合事件数。
需要说明的是，当针对每条响应线只绘制出一个能谱或子能谱时，未散射单事件数的总和则是指计算出的一个未散射单事件数本身，而当针对每条响应线绘制出多个能谱或多个子能谱时，未散射单事件数的总和则是指计算出的多个未散射单事件数的加和。
根据另一个实施例，也可以利用公式（7）-（12）计算所有响应线上的总散射符合事件数，此时，SC表示总散射符合事件数，Stot表示多条响应线上的总单事件数，并且其为总符合事件数的2倍，Sunsc表示总未散射单事件数。
通过图1所示的实施例，通过确定响应线上的符合事件内的初始单事件数据，利用校正参数构建初始单事件数据与校正单事件数据之间的对应关系，构建对应关系的目标函数，计算目标函数取最大值时校正单事件数据的取值，并且将校正单事件数据的取值确定为符合事件内的与初始单事件数据对应的散射符合事件数，从而实现对符合事件的校正，而并不需要活度图像和衰减图像，这使得符合事件的散射校正更易于实现。另外，图1所示的实施例考虑了所有的散射情况，包括多散射和外部散射，使得所得到的校正结果的准确性较高，从而可以提高后续重建图像的质量。相对于现有技术中的双散射模拟方法和蒙特卡洛模拟方法，图1所示的散射校正方法的计算过程相对简单并且计算量较小。
图3示出了根据本申请实施例的另一种散射校正方法的流程图。下面参照图3，对该散射校正方法进行详细说明。该方法包括：
为了减少图1所示的散射校正方法的计算量，在从探测器获得多条响应线上的符合事件后，执行如图3所示的步骤301，对多条响应线进行降采样以获得降采样响应线。
具体地，可以按照响应线所在的位置或响应线之间的夹角对所有响应线进行分组，然后将每个分组中的所有响应线分别合并成一条响应线。例如，将每个探头对之间的所有响应线合并成一条响应线，如图4所示，或者也可以将相邻的多个（例如，轴向相邻的2个）探头对之间的所有响应线合并成一条响应线，如图5所示；还可以将夹角小于或等于预设角度（例如，15度）的所有响应线合并成一条响应线。需要说明的是，关于相邻的多个探头对的具体数量可以根据实际需求来选取，在此并不限制。通过对所有响应线进行降采样，可以便于绘制能谱，并且减少数据计算量。
在步骤S303，利用校正参数构建降采样响应线上的符合事件内的初始单事件数据与校正单事件数据之间的对应关系。
在步骤S305，构建对应关系的目标函数，计算目标函数取最大值时校正单事件数据的取值，并且将校正单事件数据的取值确定为降采样后的符合事件内的与初始单事件数据对应的未散射单事件数。
在步骤S307，利用降采样后的符合事件数和未散射单事件数计算符合事件内的散射符合事件数。
步骤S303~S307与图1中的步骤S101-S105相对应，二者的区别仅在于处理的是降采样后的响应线上的相关数据。因此，关于步骤S303~S307的具体执行过程，可以参照上述实施例中对图1中的步骤S101-S105的相关描述，在此就不再赘述。
在本申请的另一实施例中，如图6所示，该方法还包括：
在步骤S309，利用预设方法对降采样响应线进行上采样以获得上采样响应线，并且利用降采样响应线上的散射符合事件数计算上采样响应线上的散射符合事件数。
该预设方法可以包括线性插值方法，例如，双线性插值方法或4D（4D表示响应线的轴向和径向的平移、轴向和径向的旋转）线性插值方法，也可以包括其它上采样方法，在此并不限制。另外，针对TOF散射，上述线性插值方法可以包括5D（即，在4D的基础上还包括时间差）线性插值方法。
关于对降采样响应线进行上采样的具体方法，可以参照现有技术中的相应描述，在此不再赘叙。
根据一个实施例，针对利用4D线性差值法对降采样响应线进行上采样的方式，可以按照公式（19）-（20）计算每条降采样响应线中包含的上采样响应线上的散射符合事件数。
Figure 708514DEST_PATH_IMAGE025
公式（19）
Figure 636019DEST_PATH_IMAGE026
公式（20）
其中，
Figure 774876DEST_PATH_IMAGE037
表示一条降采样响应线包含的N条上采样响应线上的散射符合事件数；SC K,L 表示一条降采样响应线上的散射符合事件数；I和J表示构成上采样响应线的“小”晶体的编号；K和L表示构成降采样响应线的“大”晶体的编号；Wtot为归一化因子；N为一条降采样响应线中包含的上采样响应线数；TI和TJ分别表示上采样晶体I和J对应的降采样中心晶体、降采样轴向相邻晶体、降采样径向相邻晶体以及降采样相对晶体的集合；
Figure 358305DEST_PATH_IMAGE027
表示对于上采样晶体I而言降采样晶体K的权重，
Figure 873600DEST_PATH_IMAGE028
表示对于上采样晶体J而言降采样晶体L的权重，并且
Figure 593294DEST_PATH_IMAGE027
和
Figure 903052DEST_PATH_IMAGE028
可以通过对降采样晶体的轴向长度、径向长度以及上采样晶体的中心距离中心晶体的径向和轴向距离进行计算得到。其中，上采样晶体为探头中的原始晶体，降采样晶体为降采样后的一个或多个晶体的集合，中心晶体为上采样晶体所在的降采样晶体，轴向相邻晶体为与上采样晶体轴向最邻近的晶体，径向相邻晶体为与上采样晶体径向最邻近的晶体，相对晶体为与上采样晶体对角方向最邻近的晶体。
根据一个实施例，也可以通过以下方式来计算上采样响应线上的散射符合事件数：首先计算每条降采样响应线上的散射符合事件数与符合事件数的比例，然后根据每条降采样响应线所包含的上采样响应线上的符合事件数和该比例计算上采样响应线上的散射符合事件数，如公式（E）所示。
Figure 786826DEST_PATH_IMAGE038
公式（E）
其中，
Figure 840233DEST_PATH_IMAGE039
表示第l条降采样响应线包含的所有上采样响应线上的散射符合事件数；
Figure 680013DEST_PATH_IMAGE040
表示第l条降采样响应线包含的所有上采样响应线上的符合事件数，也即，降采样前的对应响应线上的符合事件数；
Figure 160673DEST_PATH_IMAGE041
表示第l条降采样响应线上的散射符合事件数与符合事件数的比例，l为正整数。
在本申请的另一实施例中，可以将所有上采样响应线上的散射符合事件数相加以得到总散射符合事件数。
根据图6所示的实施例，通过对降采样响应线进行上采样，可以将降采样响应线恢复至原始状态。通过利用上述公式（19）-（20）或公式（E）计算上采样响应线上的散射符合事件数，可以使散射符合事件数的计算结果更加准确，从而可以提高散射校正结果的准确性。
图7示出了根据本申请实施例的另一种散射校正方法的流程图。如图7所示，为了更好地利用散射校正结果，在图3中的步骤S309之后，执行步骤S711。
在步骤S711，从降采样响应线所包含的上采样响应线上的符合事件数中剔除散射符合事件数。
在一个实施例中，可以针对每条降采样响应线，从其所包含的上采样响应线上的符合事件数中剔除对应的散射符合事件数。
在另一个实施例中，可以在计算出所有上采样响应线上的总散射符合事件数之后，从这些上采样响应线上的总符合事件数中剔除总散射符合事件数。
需要说明的是，图7中示出的步骤S701-S709分别与图3中的步骤S301-S309相同，在此不再赘叙。另外，上述图7仅是示意，也可以在执行图1中的步骤S105之后执行上述步骤S711。
图8示出了根据本申请实施例的又一种散射校正方法的流程图。如图8所示，为了使散射校正结果更加准确，在执行图1中的步骤S101之前，或者在执行图3中的步骤S301之前或之后，执行步骤S801。
例如，在步骤S301之前执行步骤S801，则在步骤S801，在对初始单事件数据进行校正之前，对符合事件进行随机校正，以消除随机符合事件的影响。
根据一个实施例，首先采用延迟窗的方式采集探测数据并且对所采集的探测数据进行符合处理以获得随机符合事件，然后利用可用的随机校正方法从符合事件中去除随机符合事件。
根据另一个实施例，参照前述对符合事件进行散射校正的类似方式，从符合事件中去除随机符合事件。
需要说明的是，图8中示出的步骤S803-S811分别对应于图3中的步骤S301-S309，在此不再赘叙。
下面以具体应用实例来说明根据本申请实施例所提供的散射校正方法所具有的有益效果。
图9示出了一个包含三个冷区的圆柱形假体示意图。灰色区域称为背景区域，充满放射源，且由水填充；白色区域没有放射源，且分别通过空气、水、聚四氟乙烯（PTFE）填充。这三种材料密度不同，伽马光子在其中发生康普顿散射的概率也不同。
图10和图11分别示出了针对图9中所示的假体的探测数据，利用前述OSL-EDR和OT-EDR两种方法计算出的一条降采样线上的未散射单事件数的误差图，图中误差条表示误差的标准差。其中，横坐标为一条降采样响应线上的未散射单事件数，纵坐标表示误差，且误差等于计算出的未散射单事件数和实际的未散射单事件数的差值与实际的未散射单事件数的比值。
通过图10和图11可以看出，OSL-EDR和OT-EDR这两种方法的计算误差都很小，均值都在5%以内，且误差的标准差较小。由此可见，利用这两种方法计算出的未散射单事件数的准确性较高，从而使得后续计算出的散射符合事件数的准确性也较高。
图12示出了利用公式（7）~（8）和（10）所示的方法（分别对应于图中的SSTSC-1、SSTSC-2和SSTSC-3）计算降采样后的散射符合事件数的误差图，其中，横坐标为一条降采样响应线内的散射符合事件数，纵坐标为降采样响应线集合的误差均值。从图12可以看出上述三种方法的计算误差都较小，从而可以保证散射校正结果的准确性。
图13示出了根据本申请实施例的一种散射校正方法装置的框图。下面参照图13，对该散射校正装置进行详细说明。
如图13所示，该散射校正装置包括构建单元1301、处理单元1303和计算单元1305。
构建单元1301用于利用校正参数构建探测器探测到的响应线上的符合事件内的初始单事件数据与校正单事件数据之间的对应关系，其中，初始单事件数据包括在探测器探测到的伽马光子的各个实测能量满足对应的预设条件下的单事件数。处理单元1303用于构建对应关系的目标函数，计算目标函数取最大值时校正单事件数据的取值，并且将校正单事件数据的取值确定为符合事件内的与初始单事件数据对应的未散射单事件数。第一计算单元1305用于利用符合事件数和未散射单事件数计算符合事件内的散射符合事件数。
在另一实施例中，如图14所示，该散射校正装置还可以包括降采样单元1300，其可以用于在校正单元对初始单事件数据进行校正之前对多条响应线进行降采样以获得降采样响应线。相应地，处理单元1303还可以用于计算降采后的未散射单事件数，并且计算单元1305还可以用于利用降采样后的未散射单事件数和符合事件数计算降采样响应线上的符合事件内的散射符合事件数。
在另一实施例中，如图14所示，该散射校正装置也还可以包括上采样单元1307和第二计算单元1309。其中，上采样单元1307可以用于按照预设方法对降采样响应线进行上采样以获得上采样响应线；第二计算单元1309可以用于利用降采样响应线上的散射符合事件数计算上采样响应线上的散射符合事件数。另外，第一计算单元1305和第二计算单元1309可以集成于同一个模块内。
在另一实施例中，如图15所示，该散射校正装置也还可以包括随机校正单元1200，其可以用于在对初始单事件数据进行校正之前，对符合事件进行随机校正。该随机校正单元1200可以设置在降采样单元1300与构建单元1301之间，也可以设置在降采样单元1300之前。
关于上述各个单元的详细描述，可以参照上述方法实施例中的对应描述，在此不再赘叙。
根据一些实施例，通过利用本申请实施例提供的散射校正装置，可以在不需要活度图和衰减图的情况下就可以实现对符合事件的散射校正，并且散射校正结果的准确性较高，从而可以提高后续重建图像的质量。
图16示出了根据本申请实施例的另一种散射校正装置。图16示出的散射校正装置仅仅是一个示例，不应对本申请实施例的功能和使用范围带来任何限制。
如图16所示，该散射校正装置以通用计算设备的形式表现。该散射校正装置的组件可以包括但不限于：至少一个处理器210、至少一个存储器220、连接不同系统组件（包括存储器220和处理器210）的总线230、显示单元240等。其中，存储器220存储有程序代码，程序代码可以被处理器210执行，使得处理器210执行本说明书描述的根据本申请各种示例性实施方式的方法。例如，处理器210可以执行如图1中所示的方法。
存储器220可以包括易失性存储单元形式的可读介质，例如随机存取存储单元（RAM）2201和/或高速缓存存储单元2202，还可以进一步包括只读存储单元（ROM）2203。
存储器220还可以包括具有一组（至少一个）程序模块2205的程序/实用工具2204，这样的程序模块2205包括但不限于：操作系统、一个或者多个应用程序、其它程序模块以及程序数据，这些示例中的每一个或某种组合中可能包括网络环境的实现。
总线230可以为表示几类总线结构中的一种或多种，包括存储单元总线或者存储单元控制器、外围总线、图形加速端口、处理单元或者使用多种总线结构中的任意总线结构的局域总线。
散射校正装置也可以与一个或多个外部设备300（例如键盘、指向设备、蓝牙设备等）通信，还可与一个或者多个使得用户能与该散射校正装置交互的设备通信，和/或与使得该散射校正装置能与一个或多个其它计算设备进行通信的任何设备（例如路由器、调制解调器等等）通信。这种通信可以通过输入/输出（I/O）接口250进行。而且，散射校正装置还可以通过网络适配器260与一个或者多个网络（例如局域网（LAN），广域网（WAN）和/或公共网络，例如因特网）通信。网络适配器260可以通过总线230与散射校正装置的其它模块通信。应当明白，尽管图中未示出，可以结合散射校正装置使用其它硬件和/或软件模块，包括但不限于：微代码、设备驱动器、冗余处理单元、外部磁盘驱动阵列、RAID系统、磁带驱动器以及数据备份存储系统等。
图17示出了根据本申请实施例的一种成像系统的框图，下面参照图17，对根据本申请实施例的一种成像系统进行详细说明。
如图17所示，根据本申请实施例的一种成像系统包括探测器1701和散射校正装置1703。其中，散射校正装置1703可以包括如图13、图14、图15或图16所示的散射校正装置；探测器1701可以包括PET探测器、PET-CT探测器、CT探测器或MR探测器等。关于这些探测器的详细描述可以参照现有技术，在此不再赘叙。
通过以上的实施方式的描述，本领域的技术人员易于理解，这里描述的示例实施方式可以通过软件实现，也可以通过软件结合必要的硬件的方式来实现。根据本申请实施方式的技术方案可以以软件产品的形式体现出来，该软件产品可以存储在一个计算机可读存储介质（可以是CD-ROM、U盘、移动硬盘等）中或网络上，包括若干计算机程序指令以使得一台计算设备（可以是个人计算机、服务器、或者网络设备等）执行根据本申请实施方式的上述方法。
软件产品可以采用一个或多个可读介质的任意组合。可读介质可以是可读信号介质或者可读存储介质。可读存储介质例如可以为但不限于电、磁、光、电磁、红外线、或半导体的系统、装置或器件，或者任意以上的组合。可读存储介质的更具体的例子（非穷举的列表）包括：具有一个或多个导线的电连接、便携式盘、硬盘、随机存取存储器（RAM）、只读存储器（ROM）、可擦式可编程只读存储器（EPROM或闪存）、光纤、便携式紧凑盘只读存储器(CD-ROM)、光存储器件、磁存储器件、或者上述的任意合适的组合。
计算机可读存储介质可以包括在基带中或者作为载波一部分传播的数据信号，其中承载了可读程序代码。这种传播的数据信号可以采用多种形式，包括但不限于电磁信号、光信号或上述的任意合适的组合。可读存储介质还可以是可读存储介质以外的任何可读介质，该可读介质可以发送、传播或者传输用于由指令执行系统、装置或者器件使用或者与其结合使用的程序。可读存储介质上包含的程序代码可以用任何适当的介质传输，包括但不限于无线、有线、光缆、RF等等，或者上述的任意合适的组合。
可以以一种或多种程序设计语言的任意组合来编写用于执行本申请操作的程序代码，程序设计语言包括面向对象的程序设计语言—诸如Java、C++等，还包括常规的过程式程序设计语言—诸如C语言或类似的程序设计语言。程序代码可以完全地在用户计算设备上执行、部分地在用户设备上执行、作为一个独立的软件包执行、部分在用户计算设备上部分在远程计算设备上执行、或者完全在远程计算设备或服务器上执行。在涉及远程计算设备的情形中，远程计算设备可以通过任意种类的网络，包括局域网（LAN）或广域网（WAN），连接到用户计算设备，或者，可以连接到外部计算设备（例如利用因特网服务提供商来通过因特网连接）。
上述计算机可读介质承载有一个或者多个程序指令，当上述一个或者多个程序指令被一个该设备执行时，使得该计算机可读介质实现前述功能。
本领域技术人员可以理解上述各模块可以按照实施例的描述分布于装置中，也可以进行相应变化唯一不同于本实施例的一个或多个装置中。上述实施例的多个模块可以合并为一个模块，也可以进一步将一个模块拆分成多个子模块。
根据本申请的一些示例实施例，不需要活度图和衰减图的情况下就可以实现符合事件的散射校正，并且散射校正结果的准确性较高，从而可以提高后续重建图像的质量。
虽然本申请提供了如上述实施例或流程图所述的方法操作步骤，但基于常规或者无需创造性的劳动在所述方法中可以包括更多或者更少的操作步骤。在逻辑性上不存在必要因果关系的步骤中，这些步骤的执行顺序不限于本申请实施例提供的执行顺序。
本说明书中的各个实施例均采用递进的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其它实施例的不同之处。
以上对本申请实施例进行了详细介绍，本文中应用了具体个例对本申请的原理及实施方式进行了阐述，以上实施例的说明仅用于帮助理解本申请的方法及其核心思想。同时，本领域技术人员依据本申请的思想，基于本申请的具体实施方式及应用范围上做出的改变或变形之处，都属于本申请保护的范围。综上所述，本说明书内容不应理解为对本申请的限制。

Claims (29)

1.一种散射校正方法，其特征在于，包括：
利用校正参数构建探测器探测到的响应线上的符合事件内的初始单事件数据与校正单事件数据之间的对应关系，其中，所述初始单事件数据包括探测到的伽马光子的各个实测能量满足预设条件下的单事件数；
构建所述对应关系的目标函数，计算所述目标函数取最大值时所述校正单事件数据的取值，并且将所述校正单事件数据的取值确定为所述符合事件内的与所述初始单事件数据对应的未散射单事件数；以及
利用符合事件数和所述未散射单事件数计算所述符合事件内的散射符合事件数。
2.根据权利要求1所述的散射校正方法，其特征在于，所述初始单事件数据通过以下方式来确定：
将所述符合事件拆分成单事件，统计每个所述单事件所对应的所述实测能量，并统计同一所述实测能量所对应的单事件数以获得所述初始单事件数据；
将所述符合事件拆分成单事件，根据所述单事件所在的探头位置绘制至少一个能谱，并且获取各个所述能谱中的所述初始单事件数据；或者
将所述符合事件拆分成单事件，根据所述单事件所在的探头位置绘制至少一个能谱，根据伽马光子到达所述探头的时间差将每个所述至少一个能谱划分成多个子能谱，并且获取各个所述子能谱中的所述初始单事件数据。
3.根据权利要求2所述的散射校正方法，其特征在于，针对所获取的每个所述初始单事件数据，按照公式(1)构建所述对应关系：
Figure FDA0003315378250000021
其中，
Figure FDA0003315378250000022
X为校正单事件数据，并且X＝{Xj},j＝1,2,K,n，Xj为在伽马光子的真实能量Et满足Ej≤Et＜Ej+1下的单事件数；Y为初始单事件数据，并且Y＝{Yi},i＝1,2,K,m，Yi为在实测能量Ed满足Ei≤Ed＜Ei+1下的单事件数；
Figure FDA0003315378250000023
为Y的期望值，并且
Figure FDA0003315378250000024
P为校正参数，P中的元素Pij为真实能量Et满足Ej≤Et＜Ej+1的伽马光子在实测后的实测能量Ed同时满足Ei≤Ed＜Ei+1的概率；Ej和Ej+1分别为预设能窗{Ej}中的第j个能量和第j+1个能量；Ei和Ei+1分别为实测能窗{Ei}中的第i个能量和第i+1个能量；m、n、i和j均为正整数。
4.根据权利要求3所述的散射校正方法，其特征在于，当m＝n且{Ei}＝{Ej}时，通过公式(2)计算所述校正参数P的第i行第j列的元素：
Figure FDA0003315378250000025
其中，YE＝{Ys E},s＝1,2,K,m为在忽略散射符合事件下对真实能量为Et的伽马光子进行扫描所得到的初始单事件数据；所述真实能量Et满足Ej≤Et＜Ej+1；i、j、s、m为正整数。
5.根据权利要求3或4所述的散射校正方法，其特征在于，所述伽马光子的真实能量Et包括202keV、307keV、511keV、622keV或1.27MeV。
6.根据权利要求3所述的散射校正方法，其特征在于，所述构建所述对应关系的目标函数的步骤包括：
利用极大似然法构建所述目标函数，所述目标函数通过公式(3)确定：
Φ(X)＝∑i(-∑jPijXj+Yiln(∑jPijXj))-βR(X) 公式(3)，
其中，β为超参数；R(X)为惩罚项。
7.根据权利要求6所述的散射校正方法，其特征在于，所述惩罚项R(X)按照公式(4)进行计算：
Figure FDA0003315378250000031
其中，M为常数。
8.根据权利要求7所述的散射校正方法，其特征在于，计算所述目标函数取最大值时所述校正单事件数据的取值的步骤包括通过公式(5)或公式(6)计算所述校正单事件数据的取值：
Figure FDA0003315378250000032
或
Figure FDA0003315378250000033
Figure FDA0003315378250000038
其中，
Figure FDA0003315378250000034
Figure FDA0003315378250000035
和
Figure FDA0003315378250000036
为校正单事件数据的取值；
Figure FDA0003315378250000037
为第k次迭代的中间变量；k为迭代次数。
9.根据权利要求2所述的散射校正方法，其特征在于，利用公式(7)-公式(12)计算所述符合事件内的所述散射符合事件数：
SC≈SC+ 公式(7)
SC≈SC- 公式(8)
SC≈1/[(1/SC+)+(1/SC-)] 公式(9)
Figure FDA0003315378250000041
SC≈0.5(SC++SC-) 公式(11)
或
Figure FDA0003315378250000042
其中，SC+＝Stot-Sunsc，
Figure FDA0003315378250000043
Stot为总单事件数，并且为符合事件数的2倍；Sunsc为未散射单事件数的总和；SC为散射符合事件数。
10.一种散射校正方法，其特征在于，包括：
对探测器探测到的多条响应线进行降采样以获得降采样响应线；
利用校正参数构建所述降采样响应线上的符合事件内的初始单事件数据与校正单事件数据之间的对应关系，其中，所述初始单事件数据包括在所述探测器探测到的伽马光子的各个实测能量在对应的预设条件下的单事件数；
构建所述对应关系的目标函数，计算所述目标函数取最大值时所述校正单事件数据的取值，并且将所述校正单事件数据的取值确定为所述符合事件内的与所述初始单事件数据对应的未散射单事件数；以及
利用降采样后的符合事件数和所述未散射单事件数计算所述符合事件内的散射符合事件数。
11.根据权利要求10所述的散射校正方法，其特征在于，所述方法还包括：
按照预设方法对所述降采样响应线进行上采样以获得上采样响应线；
利用所述降采样响应线上的所述散射符合事件数计算所述上采样响应线上的散射符合事件数。
12.根据权利要求11所述的散射校正方法，其特征在于，所述初始单事件数据通过以下方式来确定：
将所述符合事件拆分成单事件，统计每个所述单事件所对应的所述实测能量，并统计同一所述实测能量所对应的单事件数以获得所述初始单事件数据；
将所述符合事件拆分成单事件，根据所述单事件所在的探头位置绘制至少一个能谱，并且获取各个所述能谱中的所述初始单事件数据；或者
将所述符合事件拆分成单事件，根据所述单事件所在的探头位置绘制至少一个能谱，根据伽马光子到达所述探头的时间差将每个所述至少一个能谱划分成多个子能谱，并且获取各个所述子能谱中的所述初始单事件数据。
13.根据权利要求12所述的散射校正方法，其特征在于，针对所获取的每个所述初始单事件数据，按照公式(13)构建所述对应关系：
Figure FDA0003315378250000051
其中，
Figure FDA0003315378250000061
其中，X为校正单事件数据，并且X＝{Xj}，j＝1,2，…，n，Xj为在伽马光子的真实能量Et满足Ej≤Et＜Ej+1下的单事件数；Y为初始单事件数据，并且Y＝{YI}，i＝1,2，…，m，Yi为在实测能量Ed满足Ei≤Ed＜Ei+1下的单事件数；
Figure FDA0003315378250000062
为Y的期望值，并且
Figure FDA0003315378250000063
P为校正参数，P中的元素Pij为真实能量Et满足Ej≤Et＜Ej+1的伽马光子在实测后的实测能量Ed同时满足Ei≤Ed＜Ei+1的概率；Ej和Ej+1分别为预设能窗{Ej}中的第j个能量和第j+1个能量；Ei和Ei+1分别为实测能窗{Ei}中的第i个能量和第i+1个能量；m、n、i和j均为正整数。
14.根据权利要求13所述的散射校正方法，其特征在于，当m＝n且{Ei}＝{Ej}时，通过公式(14)计算所述校正参数P的第i行第j列的元素：
Figure FDA0003315378250000064
其中，YE＝{Ys E},s＝1,2,K,m为在忽略散射符合事件下对真实能量为Et的伽马光子进行扫描所得到的初始单事件数据；所述真实能量Et满足Ei≤Et＜Ei+1；i、j、s、m为正整数。
15.根据权利要求13或14所述的散射校正方法，其特征在于，所述伽马光子的真实能量Et包括202keV、307keV、511keV、622keV或1.27MeV。
16.根据权利要求13所述的散射校正方法，其特征在于，所述构建所述对应关系的目标函数的步骤包括：
利用极大似然法构建所述对应关系的目标函数，所述目标函数通过公式(15)确定：
Φ(X)＝∑i(-∑jPijXj+Yiln(∑jPijXj))-βR(X) 公式(15)，
其中，β为超参数；R(X)为惩罚项。
17.根据权利要求16所述的散射校正方法，其特征在于，所述惩罚项R(X)按照公式(16)进行计算：
Figure FDA0003315378250000071
其中，M为常数。
18.根据权利要求17所述的散射校正方法，其特征在于，所述计算所述目标函数取最大值时所述校正单事件数据的取值的步骤包括通过公式(17)或者公式(18)计算所述校正单事件数据的取值：
Figure FDA0003315378250000072
或，
Figure FDA0003315378250000073
其中，
Figure FDA0003315378250000074
其中，
Figure FDA0003315378250000075
和
Figure FDA0003315378250000076
均为校正单事件数据的取值；
Figure FDA0003315378250000077
为第k次迭代的中间变量；k为迭代次数。
19.根据权利要求11所述的散射校正方法，其特征在于，利用公式(7)-公式(12)计算所述符合事件内的所述散射符合事件数：
SC≈SC+ 公式(7)
SC≈SC- 公式(8)
SC≈1/[(1/SC+)+(1/SC-)] 公式(9)
Figure FDA0003315378250000081
SC≈0.5(SC++SC-) 公式(11)
或
Figure FDA0003315378250000082
其中，SC+＝Stot-Sunsc，
Figure FDA0003315378250000083
Stot为总单事件数，并且为符合事件数的2倍；Sunsc为未散射单事件数的总和；SC为散射符合事件数。
20.根据权利要求11所述的散射校正方法，其特征在于，所述预设方法包括：双线性插值方法、4D线性插值方法或5D线性插值方法。
21.根据权利要求20所述的散射校正方法，其特征在于，当所述预设方法为4D线性插值法时，按照公式(19)、公式(20)来计算所述上采样响应线上的所述散射符合事件数：
Figure FDA0003315378250000084
Figure FDA0003315378250000085
其中，
Figure FDA0003315378250000086
表示上采样响应线上的散射符合事件数；SCK,L表示降采样响应线上的散射符合事件数；I和J表示构成上采样响应线的晶体的编号；K和L表示构成降采样响应线的晶体的编号；Wtot为归一化因子；N为一条降采样响应线中包含的上采样响应线数；TI和TJ分别表示上采样晶体I和J对应的降采样中心晶体、降采样轴向相邻晶体、降采样径向相邻晶体以及降采样相对晶体的集合；
Figure FDA0003315378250000091
表示对于上采样晶体I而言降采样晶体K的权重；
Figure FDA0003315378250000092
表示对于上采样晶体J而言降采样晶体L的权重。
22.根据权利要求11所述的散射校正方法，其特征在于，所述利用所述降采样响应线上的所述散射符合事件数计算所述上采样响应线上的散射符合事件数的步骤包括：
计算每条所述降采样响应线上的所述散射符合事件数与符合事件数的比例；
根据每条所述降采样响应线所包含的所述上采样响应线上的符合事件数和所述比例计算所述上采样响应线上的所述散射符合事件数。
23.根据权利要求11所述的散射校正方法，其特征在于，所述方法还包括：
从所述降采样响应线所包含的上采样响应线上的符合事件数中剔除所述散射符合事件数。
24.根据权利要求11所述的散射校正方法，其特征在于，所述方法还包括：
在对所述初始单事件数据进行校正之前，对所述符合事件进行随机校正。
25.一种散射校正装置，其特征在于，包括：
构建单元，其被配置为利用校正参数构建探测器探测到的响应线上的符合事件内的初始单事件数据与校正单事件数据之间的对应关系，其中，所述初始单事件数据包括在所述探测器探测到的伽马光子的各个实测能量满足对应的预设条件下的单事件数；
处理单元，其被配置为构建所述对应关系的目标函数，计算所述目标函数取最大值时所述校正单事件数据的取值，并且将所述校正单事件数据的取值确定为所述符合事件内的与所述初始单事件数据对应的未散射单事件数；以及
第一计算单元，其被配置为利用符合事件数和所述未散射单事件数计算所述符合事件内的散射符合事件数。
26.一种散射校正装置，其特征在于，包括：
存储器，其上存储有程序代码；
处理器，其与所述存储器联接，并且当所述程序代码被所述处理器执行时，实现权利要求1至24中任一项所述的方法。
27.一种成像系统，其特征在于，包括：
如权利要求25或26所述的散射校正装置；
探测器，与所述散射校正装置相连接。
28.根据权利要求27所述的成像系统，其特征在于，所述探测器包括PET探测器、PET-CT探测器、CT探测器、或者MR探测器。
29.一种计算机可读存储介质，其特征在于，其上存储有程序指令，所述程序指令被执行时实现权利要求1-24中任一项所述的方法。
CN202111060795.1A 2021-09-10 2021-09-10 散射校正方法、装置、成像系统及计算机可读存储介质 Active CN113506355B (zh)

Priority Applications (3)

Application Number Priority Date Filing Date Title
CN202111060795.1A CN113506355B (zh) 2021-09-10 2021-09-10 散射校正方法、装置、成像系统及计算机可读存储介质
PCT/CN2021/124039 WO2023035360A1 (zh) 2021-09-10 2021-10-15 散射校正方法、装置、成像系统及计算机可读存储介质
PCT/CN2021/124040 WO2023035361A1 (zh) 2021-09-10 2021-10-15 图像重建方法、装置、系统及计算机可读存储介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111060795.1A CN113506355B (zh) 2021-09-10 2021-09-10 散射校正方法、装置、成像系统及计算机可读存储介质

Publications (2)

Publication Number Publication Date
CN113506355A CN113506355A (zh) 2021-10-15
CN113506355B true CN113506355B (zh) 2021-12-03

Family

ID=78017142

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111060795.1A Active CN113506355B (zh) 2021-09-10 2021-09-10 散射校正方法、装置、成像系统及计算机可读存储介质

Country Status (2)

Country Link
CN (1) CN113506355B (zh)
WO (2) WO2023035361A1 (zh)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101300601A (zh) * 2005-11-01 2008-11-05 皇家飞利浦电子股份有限公司 使用部分事件数据进行pet图像重建的方法和系统
CN110749916A (zh) * 2019-10-25 2020-02-04 上海联影医疗科技有限公司 获取pet探测器晶体时间延迟量的方法、装置和计算机设备
CN112529977A (zh) * 2020-12-04 2021-03-19 江苏赛诺格兰医疗科技有限公司 一种pet图像重建的方法和系统

Family Cites Families (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8204172B1 (en) * 2010-03-17 2012-06-19 General Electric Company System and method of prior image constrained image reconstruction using short scan image data and objective function minimization
WO2014167522A1 (en) * 2013-04-11 2014-10-16 Koninklijke Philips N.V. Method for modeling and accounting for cascade gammas in images
CN107978002B (zh) * 2016-10-25 2021-01-05 上海东软医疗科技有限公司 一种pet图像重建方法和装置
CN106491153B (zh) * 2016-12-29 2017-10-27 上海联影医疗科技有限公司 一种pet散射校正方法、pet成像方法及pet成像系统
CN107137101A (zh) * 2017-04-24 2017-09-08 沈阳东软医疗系统有限公司 一种时间校准方法和装置
US10410383B2 (en) * 2017-08-26 2019-09-10 Uih America, Inc. System and method for image data processing in positron emission tomography
CN110426729B (zh) * 2019-03-27 2021-02-26 湖北锐世数字医学影像科技有限公司 单事件校正方法、图像重建方法、装置及计算机存储介质
CN110327067B (zh) * 2019-06-10 2023-05-30 沈阳智核医疗科技有限公司 图像重建方法、装置、终端设备及pet系统
US11096634B2 (en) * 2019-08-26 2021-08-24 Siemens Medical Solutions Usa, Inc. Scatter correction based on energy response
CN110680370B (zh) * 2019-09-30 2023-02-28 东软医疗系统股份有限公司 图像重建方法、装置、控制台设备及pet系统
CN112998732B (zh) * 2021-02-08 2023-07-18 上海联影医疗科技股份有限公司 Pet数据校正方法、装置、计算机设备以及pet图像重建方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101300601A (zh) * 2005-11-01 2008-11-05 皇家飞利浦电子股份有限公司 使用部分事件数据进行pet图像重建的方法和系统
CN110749916A (zh) * 2019-10-25 2020-02-04 上海联影医疗科技有限公司 获取pet探测器晶体时间延迟量的方法、装置和计算机设备
CN112529977A (zh) * 2020-12-04 2021-03-19 江苏赛诺格兰医疗科技有限公司 一种pet图像重建的方法和系统

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于TOF PET改善低剂量CT性能的重建算法初探;程李等;《医疗卫生设备》;20190228;第40卷(第2期);第6-10页 *

Also Published As

Publication number Publication date
WO2023035361A1 (zh) 2023-03-16
CN113506355A (zh) 2021-10-15
WO2023035360A1 (zh) 2023-03-16

Similar Documents

Publication Publication Date Title
Taguchi et al. Spatio‐energetic cross‐talk in photon counting detectors: Numerical detector model (Pc TK) and workflow for CT image quality assessment
US7840053B2 (en) System and methods for tomography image reconstruction
US7888651B2 (en) Method and system for using tissue-scattered coincidence photons for imaging
Brasse et al. Correction methods for random coincidences in fully 3D whole-body PET: impact on data and image quality
Lewitt et al. Overview of methods for image reconstruction from projections in emission computed tomography
US8767908B2 (en) Exact and approximate rebinning of time-of-flight PET positron emission tomography data
CN101278318B (zh) 使用替代图像进行pet图像重建的方法和系统
WO2016161844A1 (zh) 能谱ct成像系统及数据采集和重建能谱ct图像的方法
Gillam et al. Sensitivity recovery for the AX-PET prototype using inter-crystal scattering events
Alvarez Near optimal energy selective x‐ray imaging system performance with simple detectors
Kadrmas et al. Application of reconstruction-based scatter compensation to thallium-201 SPECT: implementations for reduced reconstructed image noise
Rahman et al. Fisher information analysis of list-mode SPECT emission data for joint estimation of activity and attenuation distribution
Feng et al. Influence of Doppler broadening model accuracy in Compton camera list-mode MLEM reconstruction
Ghaly et al. Optimization of energy window and evaluation of scatter compensation methods in myocardial perfusion SPECT using the ideal observer with and without model mismatch and an anthropomorphic model observer
Ortuño et al. 3D-OSEM iterative image reconstruction for high-resolution PET using precalculated system matrix
Cajgfinger et al. Fixed forced detection for fast SPECT Monte-Carlo simulation
US20060065838A1 (en) Re-binning method for nuclear medicine imaging devices
Frandes et al. Image Reconstruction Techniques for Compton Scattering Based Imaging: An Overview [Compton Based Image Reconstruction Approaches]
Bouwens et al. Image-correction techniques in SPECT
CN113506355B (zh) 散射校正方法、装置、成像系统及计算机可读存储介质
Álvarez-Gómez et al. Fast energy dependent scatter correction for list-mode PET data
Cheng et al. Maximum likelihood activity and attenuation estimation using both emission and transmission data with application to utilization of Lu‐176 background radiation in TOF PET
Barrett Image reconstruction and the solution of inverse problems in medical imaging
Ghaly et al. Collimator optimization in myocardial perfusion SPECT using the ideal observer and realistic background variability for lesion detection and joint detection and localization tasks
Sun et al. Evaluation of image quality improvements when adding patient outline constraints into a generalized scatter PET reconstruction algorithm

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