CN109480892B - 一种图像的生成方法 - Google Patents
一种图像的生成方法 Download PDFInfo
- Publication number
- CN109480892B CN109480892B CN201811636380.2A CN201811636380A CN109480892B CN 109480892 B CN109480892 B CN 109480892B CN 201811636380 A CN201811636380 A CN 201811636380A CN 109480892 B CN109480892 B CN 109480892B
- Authority
- CN
- China
- Prior art keywords
- probability
- positron
- gamma photon
- annihilation
- gamma
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 117
- 239000013078 crystal Substances 0.000 claims abstract description 133
- 238000001514 detection method Methods 0.000 claims abstract description 84
- 230000035515 penetration Effects 0.000 claims abstract description 28
- 238000005259 measurement Methods 0.000 claims abstract description 21
- 238000004088 simulation Methods 0.000 claims abstract description 10
- 230000006870 function Effects 0.000 claims description 78
- 230000004044 response Effects 0.000 claims description 40
- 238000009826 distribution Methods 0.000 claims description 32
- 238000004590 computer program Methods 0.000 claims description 12
- 239000013598 vector Substances 0.000 claims description 8
- 239000000126 substance Substances 0.000 claims description 7
- 238000005070 sampling Methods 0.000 claims description 6
- 238000012545 processing Methods 0.000 claims description 3
- 230000008569 process Effects 0.000 abstract description 33
- 230000000694 effects Effects 0.000 abstract description 17
- 238000005457 optimization Methods 0.000 abstract description 3
- 238000002600 positron emission tomography Methods 0.000 description 17
- 238000000342 Monte Carlo simulation Methods 0.000 description 16
- 238000004364 calculation method Methods 0.000 description 12
- 230000000875 corresponding effect Effects 0.000 description 11
- 238000010586 diagram Methods 0.000 description 9
- 238000005516 engineering process Methods 0.000 description 9
- 239000000523 sample Substances 0.000 description 7
- 238000012986 modification Methods 0.000 description 3
- 230000004048 modification Effects 0.000 description 3
- 239000000700 radioactive tracer Substances 0.000 description 3
- 238000011084 recovery Methods 0.000 description 3
- 238000010521 absorption reaction Methods 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 2
- 238000006243 chemical reaction Methods 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 239000003814 drug Substances 0.000 description 2
- 239000011159 matrix material Substances 0.000 description 2
- 230000000704 physical effect Effects 0.000 description 2
- 230000005855 radiation Effects 0.000 description 2
- 238000009827 uniform distribution Methods 0.000 description 2
- 241000070928 Calligonum comosum Species 0.000 description 1
- 238000009825 accumulation Methods 0.000 description 1
- 101150043733 artA gene Proteins 0.000 description 1
- 230000002238 attenuated effect Effects 0.000 description 1
- 238000003759 clinical diagnosis Methods 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000013480 data collection Methods 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 229940079593 drug Drugs 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 230000004060 metabolic process Effects 0.000 description 1
- 238000009206 nuclear medicine Methods 0.000 description 1
- 230000000149 penetrating effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000012827 research and development Methods 0.000 description 1
Images
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/02—Arrangements for diagnosis sequentially in different planes; Stereoscopic radiation diagnosis
- A61B6/03—Computed tomography [CT]
- A61B6/037—Emission tomography
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/48—Diagnostic techniques
- A61B6/481—Diagnostic techniques involving the use of contrast agents
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B6/00—Apparatus or devices for radiation diagnosis; Apparatus or devices for radiation diagnosis combined with radiation therapy equipment
- A61B6/52—Devices using data or image processing specially adapted for radiation diagnosis
- A61B6/5205—Devices using data or image processing specially adapted for radiation diagnosis involving processing of raw data to produce diagnostic data
-
- 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
Landscapes
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Medical Informatics (AREA)
- Physics & Mathematics (AREA)
- Heart & Thoracic Surgery (AREA)
- Surgery (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Optics & Photonics (AREA)
- Pathology (AREA)
- Radiology & Medical Imaging (AREA)
- Biomedical Technology (AREA)
- Biophysics (AREA)
- Molecular Biology (AREA)
- High Energy & Nuclear Physics (AREA)
- Animal Behavior & Ethology (AREA)
- General Health & Medical Sciences (AREA)
- Public Health (AREA)
- Veterinary Medicine (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Nuclear Medicine (AREA)
- Measurement Of Radiation (AREA)
Abstract
本发明提供了一种基于点扩展函数的图像生成方法。上述图像的生成方法，上述生成方法包括：基于探测晶体空间位置及几何形状信息确定一个探测晶体对A和B探测到每一伽马光子对的概率，上述伽马光子对由在指定位置释放的正电子发生湮灭后产生；基于上述探测晶体对A和B探测到每一伽马光子对的概率确定上述点扩展函数。本发明能够简化晶体衰减和穿透效应的物理过程，仅以少量的仿真或者测量优化物理过程参数，即可获取相应的点扩展函数并基于该点扩展函数生成图像。
Description
技术领域
本发明涉及一种图像的生成方法，尤其涉及一种基于点扩展函数的图像生成方法。
背景技术
正电子发射断层成像(Positron Emission Tomography，PET)技术是临床上广泛应用的一种无创核医学影像诊断技术。该技术以正电子核素为示踪剂，通过重建或创建被扫描对象的图像，了解病灶部位对示踪剂的摄取情况，从而提供活体的新陈代谢等功能信息，在临床诊断、疗效评价、基础医学研究，以及新药研发中发挥着重要作用。
高清图像重建技术是正电子发射断层成像中的关键技术，该技术能够显著提升图像分辨率、恢复系数和信噪比等关键指标。而正电子发射断层成像系统的点扩展函数的正确获得和使用是高清图像重建技术的重要前提。
上述点扩展函数描述从单个正电子出射、发生湮灭至一对伽马光子被正电子发射断层成像探测器探测到的物理过程，具体包括：正电子自由程、一对伽马(γ)光子的非共线性效应、γ光子的飞行轨迹，以及γ光子穿透闪烁晶体过程中的吸收、衰减和散射效应等。
现有的点扩展函数的生成方法包括：解析计算、Monte Carlo仿真方法、实际测量或者其组合。上述解析计算只需获得所有探测晶体的空间位置及几何形状信息即可进行，过程简单，但是一般只能用于较为简单的物理过程，如γ光子的飞行轨迹，以及γ光子穿透闪烁晶体过程中的吸收、衰减效应等。而较复杂的物理效应一般需要通过Monte Carlo仿真方法或者实际测量获得。
由于系统响应线及重建体素的数目巨大，Monte Carlo仿真方法和实际测量都需要巨大的数据量，即需要大量的仿真计算资源或者大量的实际系统机时和复杂的实验装置(如可编程的三维移动装置)，以获得噪音较低的优质数据。而且在完成数据采集后，仍需要进行复杂的数据处理，以获得所需的物理效应模型。
通过上述Monte Carlo仿真方法和上述实际测量所获得的物理模型，仅适用于当前的系统结构。一旦在其他产品型号中稍许调整其系统结构，就需重新进行上述全部仿真和实际测量，效率低下且可重复性低。
因此，现有的正电子发射断层成像技术的点扩展函数生成方法效率低下，而且可重复性低。本领域亟需一种效率高的点扩展函数以用于生成图像，以解决现有技术存在的上述缺陷。
发明内容
以下给出一个或多个方面的简要概述以提供对这些方面的基本理解。此概述不是所有构想到的方面的详尽综览，并且既非旨在指认出所有方面的关键性或决定性要素亦非试图界定任何或所有方面的范围。其唯一的目的是要以简化形式给出一个或多个方面的一些概念以为稍后给出的更加详细的描述之序。
为了克服现有技术存在的上述缺陷，本发明旨在提供一种图像的生成方法，简化晶体衰减和穿透效应的物理过程，以去除现有技术对于大量Monte Carlo仿真方法和实际测量的依赖，从而仅以少量的仿真或者测量优化物理过程参数，即可获取相应的点扩展函数，进而生成图像。
根据本发明的一方面所提供的上述图像的生成方法，用于获取上述图像的点扩展函数，上述生成方法包括：
步骤1：基于所有探测晶体的空间位置及几何形状信息确定一个探测晶体对A和B探测到每一伽马(γ)光子对的概率，上述伽马光子对由在指定位置释放的正电子发生湮灭后产生；
步骤2：基于上述探测晶体对A和B探测到每一伽马光子对的概率确定上述点扩展函数，并基于所述点扩展函数生成图像。
优选地，在本发明所提供的上述图像的生成方法中，任一空间位置出发的正电子发生湮灭后产生的伽马光子被一探测晶体探测到的概率计算公式可以为P＝e-μLa(1-e-μLp)，其中，μ为上述伽马光子在晶体内的衰减系数，La为上述伽马光子在其它晶体内部的总衰减长度，LP为上述伽马光子在目标探测晶体内的穿透长度，上述步骤1可以包括：
步骤11：基于机械图纸获取所有探测晶体对的空间位置及几何形状信息；
步骤12：遍历所有探测晶体对和所有空间位置，其中对每一对探测晶体对的所有与响应线平行的平行线束进行射线追踪处理以获得La和LP，上述射线与上述探测晶体对定义的响应线平行；以及
步骤13：计算探测晶体对A和B探测到上述伽马光子的概率。
优选地，在本发明所提供的上述图像的生成方法中，上述步骤12可以进一步包括：
步骤121：基于各个探测晶体对的形状及位置确定各个探测晶体对的响应线；
步骤122：基于各个探测晶体对对应的响应线旋转坐标系以使得上述响应线平行于上述坐标系的某一坐标轴；以及
步骤123：沿着上述响应线的正向和反向射线追踪入射方向平行于上述响应线的伽马光子的轨迹以计算上述入射方向平行于响应线的伽马光子的总衰减长度La和穿透长度LP。
可选地，在本发明所提供的上述图像的生成方法中，上述点扩展函数可以为
其中，Pp为上述正电子发生湮灭位置的概率、Pn为上述正电子发生湮灭后产生的伽马光子的入射方向发生偏转的概率、PA和PB为上述伽马光子被探测晶体A和B探测到的概率，Ps为上述伽马光子分别在上述探测晶体A和B内部发生散射的概率，编号i为体素中心的编号，以及编号jA和jB为平行线束的编号，上述步骤2还可以包括：
步骤21：获取上述正电子发生湮灭位置的概率、上述正电子发生湮灭后产生的伽马光子的入射方向发生偏转的概率以及上述伽马光子分别在上述探测晶体A和B内部发生散射的概率。
优选地，在本发明所提供的上述图像的生成方法中，上述步骤21可以进一步包括：
步骤211：基于所述正电子发生湮灭位置的概率分布，以体素中心位置进行采样以获得上述正电子发生湮灭位置的概率；
步骤212：基于所述伽马光子发生偏转的概率分布，以过上述体素中心且与上述探测晶体对A和B的响应线平行的平行线束进行采样以获得上述正电子发生湮灭后产生的伽马光子的入射方向发生偏转的概率；以及
步骤213：基于所述探测晶体A和B内部发生散射的概率分布，以过上述体素中心且与上述探测晶体对A和B的响应线平行的平行线束进行采样以获得上述伽马光子在上述探测晶体对A和B内部发生散射的概率。
优选地，在本发明所提供的上述图像的生成方法中，上述正电子发生湮灭位置的概率分布公式为
其中，E为所述正电子的出射能量，
为湮灭密度函数，N(E)为所述正电子能量概率密度；所述伽马光子的入射方向发生偏转的概率分布公式为
其中P(θ)为所述伽马光子在X-Y平面内偏转角度θ的偏转概率，
为所述伽马光子偏离Z轴角度
的偏转概率。
优选地，在本发明所提供的上述图像的生成方法中，以三维高斯函数简化所述正电子发生湮灭位置的概率分布公式，简化后，
其中，
为所述正电子出射位置指向所述正电子湮灭位置的三维空间向量，
为三维高斯曲线的峰值位置，∑1为三维高斯曲线的宽值，
或
以二维高斯函数简化所述正电子发生湮灭位置的概率分布公式，简化后，
其中，
为所述正电子出射位置指向所述正电子湮灭位置的二维空间向量，
为二维高斯曲线的峰值位置，∑1为二维高斯曲线的宽值，
以平均分布概率简化所述偏转概率P(θ)，简化后，
以及以一维高斯函数简化所述偏转概率
简化后，
其中，μ1为期望，σ1 2为方差。
以二维高斯函数简化公式简化所述伽马光子在探测晶体对A和B内部发生散射的概率分布，简化后，
其中，
或
以一维高斯函数简化所述伽马光子在探测晶体对A和B内部发生散射的概率分布，简化后，
其中，μ2为期望，σ2 2为方差。
优选地，在本发明所提供的上述图像的生成方法中，还可以包括：
步骤4：基于蒙特·卡罗(Monte Carlo)仿真方法或实际测量获取的数据优化上述点扩展函数的参数
以及σ2 2中的任意组合。
根据本发明的一方面，本发明还提供了一种计算机装置。
本发明提供的上述计算机装置，包括存储器、处理器，以及存储在存储器上的计算机程序。上述处理器被用于执行存储在上述存储器上的计算机程序时可以实现上述任意一种图像的生成方法的步骤。
根据本发明的一方面，本发明还提供了一种计算机存储介质。
本发明提供的上述计算机存储介质，其上存储有计算机程序。上述计算机程序被执行时可以实现上述任意一种图像的生成方法的步骤。
基于以上描述，本发明可以将复杂的晶体衰减和穿透效应的物理过程进行简化，并将其加入解析计算过程中，以去除现有技术对于大量Monte Carlo仿真方法和实际测量的依赖，从而仅需少量的仿真或者测量优化物理过程参数，即可获取相应的点扩展函数。
更进一步地，本发明还可以将正电子发射断层成像系统的所有探测晶体的空间位置及几何形状信息作为输入，以完成全部的点扩展函数计算。因此，当正电子发射断层成像系统的结构发生变化时，本发明无需大量的Monte Carlo仿真方法和实际测量程序，仅需重复解析计算过程即可获取相应的点扩展函数，从而大大提高了获取点扩展函数的效率，便于生成图像。
附图说明
在结合以下附图阅读本公开的实施例的详细描述之后，能够更好地理解本发明的上述特征和优点。在附图中，各组件不一定是按比例绘制，并且具有类似的相关特性或特征的组件可能具有相同或相近的附图标记。
图1示出了本发明的一个实施例提供的图像的生成方法的流程示意图；
图2示出了本发明的一个实施例提供的确定每一伽马光子对被一探测晶体对探测到的概率的示意图；
图3示出了本发明的一个实施例提供的确定每一伽马光子对被一探测晶体对探测到的概率的流程示意图；
图4示出了本发明的一个实施例提供的伽马光子在S-T-Z立体坐标系中的传播示意图；
图5示出了本发明的一个实施例提供的旋转X-Y坐标轴以获得S-T坐标轴的示意图；
图6示出了本发明的一个实施例提供的获得伽马光子的总衰减长度和穿透长度的流程示意图；
图7示出了本发明的一个实施例提供的简化伽马光子入射方向的示意图；
图8示出了本发明的一个实施例提供的确定点扩展函数的流程示意图；
图9示出了本发明的一个实施例提供的获得正电子发生湮灭位置的概率分布，并以体素中心位置进行采样的示意图；
图10示出了本发明的一个实施例提供的简化伽马光子对非共线性效应的物理过程的示意图；
图11示出了本发明的一个实施例提供的简化伽马光子穿透晶体过程中内部散射效应的物理过程的示意图；
图12示出了本发明的一个实施例提供的利用Monte Carlo仿真方法仿真或者实测数据放置点源并采集数据的示意图；
图13示出了本发明的另一个实施例提供的计算机装置的结构示意图。
为清楚起见，以下给出附图标记的简要说明：
1-3 图像生成方法的步骤；
11-13 确定每一伽马光子对被一探测晶体对探测到的概率的步骤；
121-123 获得伽马光子的总衰减长度和穿透长度的步骤；
211-213 确定点扩展函数的步骤；
A 探测晶体A；
B 探测晶体B；
S S坐标轴；
T T坐标轴；
X X坐标轴；
Y Y坐标轴；
Z Z坐标轴；
La 伽马光子在其它晶体内部的总衰减长度；
LP 伽马光子在目标探测晶体内的穿透长度；
90 计算机装置；
91 存储器；
92 处理器。
具体实施方式
以下由特定的具体实施例说明本发明的实施方式，本领域技术人员可由本说明书所揭示的内容轻易地了解本发明的其他优点及功效。虽然本发明的描述将结合优选实施例一起介绍，但这并不代表此发明的特征仅限于该实施方式。恰恰相反，结合实施方式作发明介绍的目的是为了覆盖基于本发明的权利要求而有可能延伸出的其它选择或改造。为了提供对本发明的深度了解，以下描述中将包含许多具体的细节。本发明也可以不使用这些细节实施。此外，为了避免混乱或模糊本发明的重点，有些具体细节将在描述中被省略。
在本发明的描述中，需要说明的是，除非另有明确的规定和限定，术语“安装”、“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。
另外，在以下的说明中所使用的“上”、“下”、“左”、“右”、“顶”、“底”、“水平”、“垂直”应被理解为该段以及相关附图中所绘示的方位。此相对性的用语仅是为了方便说明之用，其并不代表其所叙述的装置需以特定方位来制造或运作，因此不应理解为对本发明的限制。
能理解的是，虽然在此可使用用语“第一”、“第二”、“第三”等来叙述各种组件、区域、层和/或部分，这些组件、区域、层和/或部分不应被这些用语限定，且这些用语仅是用来区别不同的组件、区域、层和/或部分。因此，以下讨论的第一组件、区域、层和/或部分可在不偏离本发明一些实施例的情况下被称为第二组件、区域、层和/或部分。
正电子发射断层成像(Positron Emission Tomography，PET)技术主要包括以下步骤：
首先将标记有正电子核素的放射性示踪药物注射入被检测体内。以正电子核素发生衰变发射出的正电子，与被检测对象体内的负电子发生湮没反应，从而产生两个方向相反、能量均为511KeV的伽马(γ)光子。
然后，以放置在病灶部位周围的探测单元探测上述γ光子，并对其进行电子学符合探测的处理，以记录下符合要求的γ光子对，其中一次湮灭反应被称为一个符合事件。
之后，以多个符合事件重建或创建被扫描对象的图像。也就是说，重建图像提供了被扫描对象内部放射性核素分布的信息。
为了克服现有技术存在的缺陷，本发明提供一种图像的生成方法的实施例，从而简化晶体衰减和穿透效应的物理过程，以去除现有技术对于大量蒙特·卡罗(MonteCarlo)仿真方法和实际测量的依赖，从而仅以少量的仿真或者测量优化物理过程参数，即可获取相应的点扩展函数。
本领域中的点扩展函数概念与数学领域中的函数概念是完全不同的，数学领域中的函数一般是用于表示自变量与因变量之间的对应法则，而本领域中的点扩展函数实际上是图像的另一种表现方式，获取了一图像的点扩展函数即获取了该幅图像。因此本发明获取点扩展函数的过程即是生成图像的过程。
本实施例提供的上述图像的生成方法，用于获取图像的点扩展函数。如图1所示，上述图像的生成方法可以包括：
步骤1：基于所有探测晶体的空间位置及几何形状信息确定一个探测晶体对探测到每一伽马光子对的概率。
上述探测晶体对可以由如图2所示的晶体A和晶体B组成。晶体A和晶体B泛指目标晶体对。上述伽马光子对由在指定位置释放的正电子发生湮灭后产生。
如图2所示，伽马光子的入射方向沿箭头方向延伸。伽马光子首先在晶体B内发生衰减，其衰减长度为La；然后再穿透晶体A，其穿透深度为Lp。
将晶体对伽马光子的衰减系数定义为μ，并将伽马光子经过上述衰减和穿透效应后被晶体A探测到的概率定义为P，则从任一空间位置出发的正电子，在发生湮灭后产生的伽马光子被探测晶体A探测到的概率计算公式可以为P＝e-μLa(1-e-μLp)，其中，La为上述伽马光子在其它晶体内部的总衰减长度，LP为上述伽马光子在目标探测晶体内的穿透长度。
如图3所示，上述步骤1提供的一个探测晶体对探测到每一伽马光子对的概率，可以具体地通过以下方法来确定：
步骤11：基于机械图纸获取所有探测晶体对的空间位置及几何形状信息，机械图纸为表征所述有探测晶体的空间位置及几何形状信息。
如图4所示，上述所有探测晶体对的形状及位置包括：上述所有探测晶体对在S-T平面内的形状及位置，以及上述所有探测晶体对在T-Z平面内的形状及位置。
本领域的技术人员可以理解，上述S-T坐标轴可以是通过将PET系统初始的X-Y-Z立体坐标系中的X-Y坐标轴，如图5所示地旋转角度
后所获得的，主要用于使上述探测晶体定义的响应线在X-Y平面(即S-T平面)内的投影线与T轴平行，从而便于计算每一伽马光子对被探测晶体A探测到的概率。在其他实施例中，本领域的技术人员即使不采用旋转X-Y坐标轴的方案，也能够实现确定一个探测晶体对探测到每一伽马光子对的概率的目的。
如图3所示，上述步骤1提供的确定一个探测晶体对探测到每一伽马光子对概率的方法还可以包括步骤：
步骤12：遍历所有探测晶体对和所有空间位置，其中对每一对探测晶体对的所有响应线进行射线追踪处理，以获得伽马光子的总衰减长度La和穿透长度LP，上述射线可以与上述探测晶体定义的响应线平行。
如图6所示，上述步骤12提供的伽马光子的总衰减长度La和穿透长度LP，可以通过以下方法来获得：
步骤121：基于各个探测晶体对的形状及位置，确定各个探测晶体对的响应线；
步骤122：基于各个探测晶体对对应的响应线，如图5所示地旋转坐标系以使得上述响应线平行于上述坐标系的一坐标轴(例如：上述S轴或T轴)；以及
步骤123：沿着上述响应线的正向和反向射线，追踪入射方向平行于上述响应线的伽马光子的轨迹，以计算上述入射方向平行于响应线的伽马光子的总衰减长度La和穿透长度LP。
如图7所示，在上述步骤123中，可以假设伽马光子的入射方向全部与响应线平行，从而对其入射方向进行简化。
在分别沿着上述响应线的正向和反向射线，追踪入射方向平行于上述响应线的伽马光子的轨迹的过程中，若伽马光子被其他晶体阻挡，则计算该衰减长度，并对其进行累加；若伽马光子接触到目标晶体，则判定累加结束，并计算获得总衰减长度La。在伽马光子进入目标晶体后，可以计算获得穿透深度Lp。
本领域的技术人员可以理解，上述步骤121-123所提供的获得伽马光子的总衰减长度La和穿透长度LP的方法，只是本实施例基于三维情况提供的一种优选方案，主要用于简化伽马光子的入射方向，从而更简洁、高效地获得伽马光子的总衰减长度La和穿透长度LP。本领域的技术人员可以理解，本案所述的实施方式并不限于三维情况，亦适用于其它情况。
在其他实施例中，本领域的技术人员也可以采用其他方法，以计算伽马光子的总衰减长度La和穿透长度LP。
基于上述步骤，完成了基于机械图纸即可获得探测晶体对A和B探测到伽马光子对的概率，复杂程度远远小于现有技术中通过仿真或实际测量获取探测晶体对A和B探测到伽马光子对的概率的过程，进而大大简化了图像的生成过程。
如图3所示，上述步骤1提供的确定每一伽马光子对概率的方法还可以包括步骤：
步骤13：计算探测晶体对A和B探测到上述伽马光子的概率。
上述步骤13可以根据上述步骤121-123获得的伽马光子的总衰减长度La和穿透长度LP，从而进一步地确定每一伽马光子对被探测晶体对A和B探测到的概率。
本领域的技术人员可以理解，上述步骤11-13所提供的确定每一伽马光子对的概率的方法，只是本实施例提供的一种优选方案，主要用于更简洁、高效地确定每一伽马光子对目标探测晶体探测到的概率。在其他实施例中，本领域的技术人员也可以通过其他点扩展函数的参数，来确定每一伽马光子对被目标探测晶体探测到的概率。
如图1所示，在本实施例提供的上述图像的生成方法中，还可以包括：
步骤2：基于上述探测晶体对A和B探测到每一伽马光子对的概率确定上述点扩展函数，并基于该点扩展函数生成图像。
相应地，本领域的技术人员还可以进一步地通过以下步骤21，以获得更多点扩展函数的具体参数：
步骤21：获取正电子发生湮灭位置的概率、正电子发生湮灭后产生的伽马光子的入射方向发生偏转的概率，以及伽马光子分别在上述探测晶体A和B内部发生散射的概率，从而结合上述探测晶体对A和B探测到每一伽马光子对的概率，以更准确地确定上述点扩展函数。
如图8所示，本实施例提供的上述优选的点扩展函数，可以通过以下方法来确定：
步骤211：基于所述正电子发生湮灭位置的概率分布，以体素中心位置进行采样，以获得上述正电子发生湮灭位置的概率。
如图9所示，根据上述步骤211可以获得正电子发生湮灭位置的概率分布，并以体素中心位置进行采样，其中，星号表示正电子释放位置；圆点表示体素中心位置；虚线表示与响应线平行，且过体素中心位置的平行线束。
上述正电子发生湮灭位置的概率，可以通过概率分布公式
来获得，其中，E为上述正电子的出射能量，
为湮没密度函数，N(E)为上述正电子的出射能量概率密度。
对于某一给定的正电子出射能量E，湮没密度函数
可以用高斯函数或者双高斯函数或者其他函数进行简化。
比如以三维高斯函数简化上述概率分布公式，简化后可获得：
其中，
为所述正电子出射位置指向所述正电子湮灭位置的三维空间向量，
为三维高斯曲线的峰值位置，∑1为三维高斯曲线的宽值，
因此，上述正电子发生湮灭位置的概率，可以进一步通过计算公式
来获得，从而简化上述点扩展函数正电子自由程的物理过程，其中，
x、y、z为向量
在PET系统初始的X-Y-Z立体坐标系中的X、Y、Z坐标轴上的投影值；∑1为协方差矩阵，|∑1|为协方差行列式的值。上述正电子发生湮灭位置的概率的计算公式，以
和∑1为参数。
此处以三维情况为例，但不限于三维情况。例如，在二维情况下，以二维高斯函数简化上述正电子发生湮灭位置的概率分布公式，简化后，
其中，
为所述正电子出射位置指向所述正电子湮灭位置的二维空间向量，
为二维高斯曲线的峰值位置，∑1为二维高斯曲线的宽值，
进一步地，所述正电子发生湮灭位置的概率计算公式为
其中，E为所述正电子的出射能量，
步骤212：基于所述伽马光子发生偏转的概率分布，以过上述体素中心且与上述探测晶体对A和B的响应线平行的平行线束进行采样，以获得上述正电子发生湮灭后产生的伽马光子的入射方向发生偏转的概率。
如图10所示，定义正负电子对湮灭产生的伽马光子对的其中一个伽马光子沿Z轴方向出射，而另一个伽马光子沿
方向出射的概率为
其中，
上述P(θ)可以用均匀分布或者其他分布概率简化，以均匀分布为例，简化后：
而上述
可以用高斯函数或者双高斯函数或者其他函数进行简化，以一维高斯函数简化
简化后：
其中，μ1为期望，σ1 2为方差。可以理解，此处以一维情况为例，但不限于一维情况。
因此，上述正电子发生湮灭后产生的伽马光子的入射方向发生偏转的概率，可以通过计算公式
来获得，从而简化上述点扩展函数伽马光子对非共线性效应的物理过程，其中，μ1和σ1 2为上述正电子发生湮灭后产生的伽马光子的入射方向发生偏转的概率的计算公式的参数。
步骤213：基于所述探测晶体A和B内部发生散射的概率分布，以过上述体素中心且与上述探测晶体对A和B的响应线平行的平行线束进行采样以获得上述伽马光子在上述探测晶体对A和B内部发生散射的概率。
如图11所示，定义伽马光子沿Z轴方向入射时，晶体内部散射发生概率为
上述伽马光子在探测晶体对A和B内部发生散射的概率，可用高斯函数或者双高斯函数或者其他函数进行简化。
比如，以二维高斯函数简化公式简化所述伽马光子在探测晶体对A和B内部发生散射的概率分布，简化后：
其中，
x、y为向量
在PET系统初始的X-Y-Z立体坐标系中的X、Y坐标轴上的投影值；∑2为协方差矩阵，|∑2|为协方差行列式的值，
其中
和∑2为上述晶体内部散射发生概率的计算公式的参数。
此处以二维情况为例，但不限于二维情况。比如，在一维情况下，以一维高斯函数简化所述伽马光子在探测晶体对A和B内部发生散射的概率分布，简化后，
其中，μ2为期望，σ2 2为方差。
本领域的技术人员可以理解，上述步骤211-213提供的确定点扩展函数的方法，只是本实施例提供的一种优选方案，主要用于简化上述点扩展函数从单个正电子湮灭至一对正电子发射断层成像探测器接收到信号的物理过程。在其他实施例中，本领域的技术人员也可以采用其他方法来简化上述点扩展函数。可以理解的是，即使不采用任何上述步骤211-213提供的简化上述点扩展函数的方法，也不会影响基本的简化晶体衰减和穿透效应的物理过程的目的。
根据上述步骤211-213提供的简化上述点扩展函数的物理模型的优选方案，并结合实际的物理过程，可以解析计算出相应的点扩展函数：
其中，Pp为上述正电子发生湮灭位置的概率、Pn为上述正电子发生湮灭后产生的伽马光子的入射方向发生偏转的概率、PA和PB分别为上述伽马光子被探测晶体A和探测晶体B探测到的概率，Ps为上述伽马光子分别在上述探测晶体A和探测晶体B内部发生散射的概率，编号i为体素中心的编号，以及编号jA和jB为平行线束的编号。
本领域的技术人员可以理解，根据上述步骤11-13获得的探测晶体对A和B探测到每一伽马光子对的概率PA和PB、上述步骤211获得的上述正电子发生湮灭位置的概率Pp、上述步骤212获得的上述正电子发生湮灭后产生的伽马光子的入射方向发生偏转的概率Pn，以及上述步骤213获得的上述伽马光子分别在上述探测晶体A和探测晶体B内部发生散射的概率Ps，所确定的点扩展函数
只是本实施例提供的一种优选方案，主要用于进一步地简化上述点扩展函数的物理过程，从而进一步地提升上述点扩展函数的获取效率；并进一步地优化上述点扩展函数，从而进一步提升PET系统的图像分辨率、恢复系数和信噪比等关键指标。
在其他实施例中，本领域的技术人员也可以仅采用其中的一种或任意多种优选的简化方案的组合，以取得相应的效果；即使不采用上述任何优选的简化方案，也不影响实现本发明基本的简化晶体衰减和穿透效应的物理过程的目的。
如图1所示，在本实施例提供的上述图像的生成方法中，还可以包括：
步骤3：基于Monte Carlo仿真方法或实际测量获取的数据优化上述点扩展函数的参数
以及σ2 2中的任意组合。
可以理解，该些基于Monte Carlo仿真方法或实际测量获取的数据优化的参数是基于上述多个物理过程采用的对应简化方法的参数来组合的。
如图12所示，根据点源放置位置和当前的点扩展函数，可以获得PET系统中所有响应线上的计数。将上述所有响应线上的计数与Monte Carlo仿真方法或者实测数据进行比较，并以两者的差距最小作为优化所有参数的依据，即使函数F达到最小：
其中，编号p为所放置的点源的编号，编号k为响应线的编号，Ip,k为根据点源放置位置和当前的点扩展函数计算得到的第k根响应线所收到的第p个点源的计数，Rp,k为Monte Carlo仿真方法或实际采集中的计数。
本领域的技术人员可以理解，上述步骤3提供的基于Monte Carlo仿真方法或实际测量获取的数据，以优化上述点扩展函数的参数的方案，只是本实施例提供的一种优选方案，主要用于进一步地优化上述点扩展函数，从而进一步提升PET系统的图像分辨率、恢复系数和信噪比等关键指标。在其他实施例中，本领域的技术人员即使不进行上述优化，也不影响实现本发明基本的简化晶体衰减和穿透效应的物理过程的目的。
基于以上描述，本实施例提供的上述图像的生成方法，可以将复杂的晶体衰减和穿透效应的物理过程进行简化，并将其加入解析计算过程中，以去除现有技术对于大量Monte Carlo仿真方法和实际测量的依赖，从而仅需少量的仿真或者测量优化物理过程参数，即可获取相应的点扩展函数。
更进一步地，本实施例提供的上述图像的生成方法，还可以将正电子发射断层成像系统的所有探测晶体的空间位置及几何形状信息作为输入，以完成全部的点扩展函数计算。因此，当正电子发射断层成像系统的结构发生变化时，本发明无需大量的Monte Carlo仿真方法和实际测量程序，仅需重复解析计算过程即可获取相应的点扩展函数，从而大大提高了获取点扩展函数的效率。
根据本发明的一方面，本发明还提供了一种计算机装置的实施例。
本实施例提供的上述计算机装置90，如图13所示，包括存储器91和处理器92，以及存储在存储器91上的计算机程序。处理器92被用于执行存储在上述存储器91上的计算机程序时，可以实现上述任意一种图像的生成方法的步骤。
根据本发明的一方面，本发明还提供了一种计算机存储介质的实施例。
本实施例提供的上述计算机存储介质，其上存储有计算机程序。上述计算机程序被执行时可以实现上述任意一种图像的生成方法的步骤。
尽管为使解释简单化将上述方法图示并描述为一系列动作，但是应理解并领会，这些方法不受动作的次序所限，因为根据一个或多个实施例，一些动作可按不同次序发生和/或与来自本文中图示和描述或本文中未图示和描述但本领域技术人员可以理解的其他动作并发地发生。
提供之前的描述是为了使本领域中的任何技术人员均能够实践本文中所描述的各种方面。但是应该理解，本发明的保护范围应当以所附权利要求书为准，而不应被限定于以上所解说实施例的具体结构和组件。本领域技术人员在本发明的精神和范围内，可以对各实施例进行各种变动和修改，这些变动和修改也落在本发明的保护范围之内。
Claims (12)
1.一种图像的生成方法，所述生成方法包括：
步骤1：基于探测晶体的空间位置及几何形状信息确定一探测晶体对A和B探测到每一伽马光子对的概率，所述伽马光子对由在指定位置释放的正电子发生湮灭后产生；
其中，所述步骤1包括：
步骤11：基于机械图纸获取所有探测晶体对的空间位置及几何形状信息；
步骤12：遍历所有探测晶体对和所有空间位置，其中对每一对探测晶体对的所有与响应线平行的平行线束进行射线追踪处理以获得La和LP，所述射线与所述探测晶体对定义的响应线平行；以及
步骤13：基于La和LP计算探测晶体对A和B探测到所述伽马光子的概率，La为所述伽马光子在其它晶体内部的总衰减长度，LP为所述伽马光子在探测晶体A内的穿透长度；
步骤2：基于所述探测晶体对A和B探测到每一伽马光子对的概率确定点扩展函数；
步骤3：基于所述点扩展函数生成图像。
2.如权利要求1所述的生成方法，其特征在于，任一空间位置出发的正电子发生湮灭后产生的伽马光子被探测晶体A探测到的概率计算公式为P＝e-μLa(1-e-μLp)，其中，μ为所述伽马光子在晶体内的衰减系数。
3.如权利要求1所述的生成方法，其特征在于，所述步骤12包括：
步骤121：基于各个探测晶体对的形状及位置确定各个探测晶体对的响应线；
步骤122：基于各个探测晶体对对应的响应线旋转坐标系以使得所述响应线平行于所述坐标系的一坐标轴；以及
步骤123：沿着所述响应线的正向和反向射线追踪入射方向平行于所述响应线的伽马光子的轨迹以计算所述入射方向平行于响应线的伽马光子的总衰减长度La和穿透长度LP。
4.如权利要求1所述的生成方法，其特征在于，所述点扩展函数为
其中，Pp为所述正电子发生湮灭位置的概率、Pn为所述正电子发生湮灭后产生的伽马光子的入射方向发生偏转的概率、PA和PB为所述伽马光子被探测晶体A和B探测到的概率，Ps为所述伽马光子分别在所述探测晶体A和B内部发生散射的概率，编号i为体素中心的编号以及编号jA和jB为平行线束的编号，所述步骤2还包括：
步骤21：获取所述正电子发生湮灭位置的概率、所述正电子发生湮灭后产生的伽马光子的入射方向发生偏转的概率以及所述伽马光子分别在所述探测晶体A和B内部发生散射的概率。
5.如权利要求4所述的生成方法，其特征在于，所述步骤21包括：
步骤211：基于所述正电子发生湮灭位置的概率分布，以体素中心位置进行采样以获得所述正电子发生湮灭位置的概率；
步骤212：基于所述伽马光子发生偏转的概率分布，以过所述体素中心且与所述探测晶体对A和B的响应线平行的平行线束进行采样以获得所述正电子发生湮灭后产生的伽马光子的入射方向发生偏转的概率；以及
步骤213：基于所述探测晶体A和B内部发生散射的概率分布，以过所述体素中心且与所述探测晶体对A和B的响应线平行的平行线束进行采样以获得所述伽马光子在所述探测晶体对A和B内部发生散射的概率。
6.如权利要求5所述的生成方法，其特征在于，所述正电子发生湮灭位置的概率分布公式为
其中，E为所述正电子的出射能量，
为湮灭密度函数，N(E)为所述正电子能量概率密度；所述伽马光子的入射方向发生偏转的概率分布公式为
其中P(θ)为所述伽马光子在X-Y平面内偏转角度θ的偏转概率，
为所述伽马光子偏离Z轴角度
的偏转概率。
7.如权利要求6所述的生成方法，其特征在于，以三维高斯函数简化所述正电子发生湮灭位置的概率分布公式，简化后，
其中，
为所述正电子出射位置指向所述正电子湮灭位置的三维空间向量，
为三维高斯曲线的峰值位置，∑1为三维高斯曲线的宽值，
或
以二维高斯函数简化所述正电子发生湮灭位置的概率分布公式，简化后，
其中，
为所述正电子出射位置指向所述正电子湮灭位置的二维空间向量，
为二维高斯曲线的峰值位置，∑1为二维高斯曲线的宽值，
8.如权利要求6所述的生成方法，其特征在于，以平均分布概率简化所述偏转概率P(θ)，简化后，
以及
以一维高斯函数简化所述偏转概率
简化后，
其中，μ1为期望，σ1 2为方差。
9.如权利要求6所述的生成方法，其特征在于，以二维高斯函数简化公式简化所述伽马光子在探测晶体对A和B内部发生散射的概率分布，简化后，
或
以一维高斯函数简化所述伽马光子在探测晶体对A和B内部发生散射的概率分布，简化后，
其中，μ2为期望，σ2 2为方差。
10.如权利要求7-9中任一项所述的生成方法，其特征在于，还包括：
步骤4：基于蒙特·卡罗仿真方法或实际测量获取的数据优化所述点扩展函数的参数
μ1、σ1 2、
μ2以及σ2 2中的任意组合。
11.一种计算机装置，包括存储器、处理器以及存储在存储器上的计算机程序，其特征在于，所述处理器被用于执行存储在所述存储器上的计算机程序时实现如权利要求1～10中任一项所述的生成方法的步骤。
12.一种计算机存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被执行时实现如权利要求1-10中任一项所述的生成方法的步骤。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811636380.2A CN109480892B (zh) | 2018-12-29 | 2018-12-29 | 一种图像的生成方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811636380.2A CN109480892B (zh) | 2018-12-29 | 2018-12-29 | 一种图像的生成方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109480892A CN109480892A (zh) | 2019-03-19 |
| CN109480892B true CN109480892B (zh) | 2022-11-04 |
Family
ID=65713298
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201811636380.2A Active CN109480892B (zh) | 2018-12-29 | 2018-12-29 | 一种图像的生成方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109480892B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112258506B (zh) * | 2020-11-18 | 2022-11-25 | 上海培云教育科技有限公司 | 基于数值计算的正电子发射型断层成像仿真方法及系统 |
Family Cites Families (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| EP1521982A1 (en) * | 2002-07-17 | 2005-04-13 | European Organization for Nuclear Research | Gamma ray detector for positron emission tomography (pet) and single photon emmission computed tomography (spect) |
| EP1922566B1 (en) * | 2005-08-18 | 2013-02-13 | Koninklijke Philips Electronics N.V. | Positron emission tomography time-of-flight list mode reconstruction with detector response function |
| US8389943B2 (en) * | 2008-09-29 | 2013-03-05 | Siemens Medical Solutions Usa, Inc. | Modeling of the point-spread-function in single-pinhole and multi-pinhole spect reconstruction |
| US8975587B2 (en) * | 2011-03-10 | 2015-03-10 | Shimadzu Corporation | Positron CT apparatus and a reconstruction method |
| US8958622B2 (en) * | 2012-05-09 | 2015-02-17 | Siemens Medical Solutions Usa, Inc. | Efficient point spread function model for multi-channel collimator in photon imaging including extra-geometric photons |
| CN103164863B (zh) * | 2013-04-02 | 2016-03-02 | 中国科学院高能物理研究所 | 用于重建正电子发射计算机断层成像图像的方法 |
| CN107260194B (zh) * | 2016-04-08 | 2020-08-28 | 山西锦地裕成医疗设备有限公司 | 用于正电子断层成像的方法和图像重建方法及系统 |
| CN107202805B (zh) * | 2017-05-31 | 2020-05-05 | 中国人民解放军信息工程大学 | 基于卷积核的锥束ct散射伪影校正方法 |
| CN108615250B (zh) * | 2018-05-31 | 2022-04-26 | 上海联影医疗科技股份有限公司 | 图像重建方法、装置、系统和计算机可读存储介质 |
-
2018
- 2018-12-29 CN CN201811636380.2A patent/CN109480892B/zh active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN109480892A (zh) | 2019-03-19 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Tong et al. | Image reconstruction for PET/CT scanners: past achievements and future challenges | |
| JP5726910B2 (ja) | 陽子コンピュータ断層撮影のためのシステム及び方法 | |
| US8847166B2 (en) | Imaging device using gamma rays, image signal processor, and image processing method for gamma ray measurement data | |
| CN108615250B (zh) | 图像重建方法、装置、系统和计算机可读存储介质 | |
| CN101681520A (zh) | Pet局部断层摄影 | |
| CN109658390B (zh) | 一种用于正电子检测正弦矩阵图的感兴趣区域提取方法 | |
| CN113112558B (zh) | 一种高清pet图像重建方法 | |
| CN109480892B (zh) | 一种图像的生成方法 | |
| Szlávecz et al. | GPU-based acceleration of the MLEM algorithm for SPECT parallel imaging with attenuation correction and compensation for detector response | |
| EP2360643A1 (en) | Methods and systems for image reconstruction | |
| CN110215223A (zh) | 散射校正方法、系统、可读存储介质和设备 | |
| US9870627B2 (en) | System and method for limited angle positron emission tomography | |
| Xiao et al. | A study on scattering correction for γ-photon 3D imaging test method | |
| Kapadia | Accuracy and patient dose in neutron stimulated emission computed tomography for diagnosis of iron overload: simulations in GEANT4 | |
| CN116577819B (zh) | 一种多头康普顿探测方法及系统 | |
| Brusaferri | Improving quantification in non-TOF 3D PET/MR by incorporating photon energy information | |
| Nasirzadeh et al. | Modeling GE advance PET-scanner using FLUKA simulation code | |
| JPWO2020059135A1 (ja) | データ処理方法、プログラム、データ処理装置および陽電子放出断層撮像装置 | |
| Hirano | Applications to development of PET/SPECT system by use of Geant4 | |
| Francis | 3D Volumetric Reconstruction for Light-Field SPECT | |
| Holmberg | Optimisation of image acquisition and reconstruction of 111In-pentetrotide SPECT | |
| JP2010203807A (ja) | 画像処理装置、画像処理システム、画像処理方法およびプログラム | |
| Torres Espallardo | Image reconstruction and correction methods for madpet-ii based on monte carlo techniques | |
| Zhu | Improved Modeling and Image Generation for Fluorescence Molecular Tomography (FMT) and Positron Emission Tomography (PET) | |
| Curkic | Evaluation of internal dosimetry for 225-Ac using one single measurement based on 111-In imaging |
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 | ||
| 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. |
|
| GR01 | Patent grant | ||
| GR01 | Patent grant |