CN110811667B - 一种基于gpu加速的高精度pet重建方法及装置 - Google Patents
一种基于gpu加速的高精度pet重建方法及装置 Download PDFInfo
- Publication number
- CN110811667B CN110811667B CN201911286593.1A CN201911286593A CN110811667B CN 110811667 B CN110811667 B CN 110811667B CN 201911286593 A CN201911286593 A CN 201911286593A CN 110811667 B CN110811667 B CN 110811667B
- Authority
- CN
- China
- Prior art keywords
- image
- estimated
- projection
- iteration
- gpu
- 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 75
- 230000001133 acceleration Effects 0.000 title claims abstract description 35
- 238000004364 calculation method Methods 0.000 claims abstract description 29
- 230000035945 sensitivity Effects 0.000 claims abstract description 20
- 230000006870 function Effects 0.000 claims description 29
- 239000011159 matrix material Substances 0.000 claims description 22
- 230000008569 process Effects 0.000 claims description 21
- 238000012937 correction Methods 0.000 claims description 18
- 238000005316 response function Methods 0.000 claims description 13
- 238000007781 pre-processing Methods 0.000 claims description 11
- 238000012545 processing Methods 0.000 claims description 8
- 230000004044 response Effects 0.000 claims description 8
- 230000000704 physical effect Effects 0.000 claims description 6
- 208000032825 Ring chromosome 2 syndrome Diseases 0.000 claims description 5
- 238000013178 mathematical model Methods 0.000 claims description 5
- 230000000694 effects Effects 0.000 claims description 4
- 208000035217 Ring chromosome 1 syndrome Diseases 0.000 claims description 3
- 238000012804 iterative process Methods 0.000 claims description 3
- 238000003384 imaging method Methods 0.000 abstract description 4
- HPTJABJPZMULFH-UHFFFAOYSA-N 12-[(Cyclohexylcarbamoyl)amino]dodecanoic acid Chemical compound OC(=O)CCCCCCCCCCCNC(=O)NC1CCCCC1 HPTJABJPZMULFH-UHFFFAOYSA-N 0.000 description 4
- 238000010586 diagram Methods 0.000 description 4
- 238000004891 communication Methods 0.000 description 3
- 238000007667 floating Methods 0.000 description 3
- 238000013519 translation Methods 0.000 description 3
- 238000007476 Maximum Likelihood Methods 0.000 description 2
- 230000008878 coupling Effects 0.000 description 2
- 238000010168 coupling process Methods 0.000 description 2
- 238000005859 coupling reaction Methods 0.000 description 2
- 238000000265 homogenisation Methods 0.000 description 2
- 238000005457 optimization Methods 0.000 description 2
- 230000006978 adaptation Effects 0.000 description 1
- 239000013078 crystal Substances 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 210000000056 organ Anatomy 0.000 description 1
- 229940121896 radiopharmaceutical Drugs 0.000 description 1
- 239000012217 radiopharmaceutical Substances 0.000 description 1
- 230000002799 radiopharmaceutical effect Effects 0.000 description 1
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/52—Devices using data or image processing specially adapted for radiation diagnosis
- A61B6/5211—Devices using data or image processing specially adapted for radiation diagnosis involving processing of medical diagnostic data
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/15—Correlation function computation including computation of convolution operations
-
- 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
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B2503/00—Evaluating a particular growth phase or type of persons or animals
- A61B2503/40—Animals
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Medical Informatics (AREA)
- Theoretical Computer Science (AREA)
- General Health & Medical Sciences (AREA)
- High Energy & Nuclear Physics (AREA)
- Pathology (AREA)
- Radiology & Medical Imaging (AREA)
- Biomedical Technology (AREA)
- Heart & Thoracic Surgery (AREA)
- Molecular Biology (AREA)
- Surgery (AREA)
- Animal Behavior & Ethology (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Public Health (AREA)
- Veterinary Medicine (AREA)
- Mathematical Physics (AREA)
- Optics & Photonics (AREA)
- Biophysics (AREA)
- Computational Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Algebra (AREA)
- Computing Systems (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Apparatus For Radiation Diagnosis (AREA)
- Nuclear Medicine (AREA)
Abstract
本发明涉及GPU成像领域，具体涉及一种基于GPU加速的高精度PET重建方法及装置。该方法及装置输入初始的估计图像，对估计图像进行正投影操作，获得估计投影；将获得的估计投影与实际测量投影作比较，获得比较结果；对比较结果进行反投影操作，获得校正图像；将校正图像与原始图像对应体素相乘，然后除敏感度图像，完成一次迭代，得到的图像就是下一次迭代的输入估计图像。通过GPU加速策略将PET迭代重建算法改进后进行快速计算，获得高分辨率、高精度的三维PET图像。
Description
技术领域
本发明涉及GPU成像领域，具体而言，涉及一种基于GPU加速的高精度PET重建方法及装置。
背景技术
小动物PET是小型化在临床上显示生理结构的成像机器，使研究人员能够量化多种器官中的生化活性，目前已经在临床前成像研究、生物学和医学基础研究、新的放射性药物开发上广泛应用。迭代重建算法是PET重建中的主流方法之一，相较于解析图像重建，它可以通过对高能光子产生和探测过程进行更精确的物理和统计建模，提升重建图像质量。GPU是为了解决图像处理任务而被设计出来的高度并行化的多线程多核心处理器，浮点计算能力显著高于CPU。越来越多的研究者使用GPU进行PET重建计算，Zhou Jian等用GPU计算正投影和反投影函数，证明能极大节省计算时间。
传统的迭代重建方法使用CPU计算系统矩阵，对估计图像进行正反投影，在迭代中接近真实分布图像，其重建速度依赖于输入sinogram的维度大小和生成的图像维度大小，以及迭代的次数。高精度带有DOI功能的PET能输出光子在晶体中的深度信息，分辨率比一般小动物PET高，因此sinogram的维度非常大，传统算法计算时间需要数天。
发明内容
本发明实施例提供了一种基于GPU加速的高精度PET重建方法及装置，以至少解决现有使用GPU进行PET重建计算时间长的技术问题。
根据本发明的一实施例，提供了一种基于GPU加速的高精度PET重建方法，包括以下步骤：
S100:输入初始的估计图像，对估计图像进行正投影操作，获得估计投影；
S200:将获得的估计投影与实际测量投影作比较，获得比较结果；
S300:对比较结果进行反投影操作，获得校正图像；
S400:将校正图像与原始图像对应体素相乘，然后除敏感度图像，完成一次迭代，得到的图像就是下一次迭代的输入估计图像。
进一步地，方法还包括步骤：
S500:将下一次迭代的输入估计图像返回步骤S100进行迭代处理，最后一次迭代获得的估计图像即是最终PET重建的结果。
进一步地，方法还包括步骤：
S001:对初始的估计图像进行预处理，将初始的估计图像进行统计和合并，生成符合预设要求的四维弦图Sinogram，该四维弦图的每一个元素代表一类响应线Line OfResponse，四维分别是LOR的平面内距中心距离d、平面内角度phi、接收第一个光子的探测器环编号ring1、接收第二个光子的探测器环编号ring2。
进一步地，在迭代过程中估计出图像中每个体素内源的强度值x，构建一个期望是强度值x的泊松分布A，其数学模型为：
x和y的期望的关系是：
其中,k是图像中每个体素内源的出现次数，e是自然常数，y是实际测量的投影值，矩阵p是与探测器几何结构相关的表示在某一位置被探测器探测到的可能性，d是LOR的平面内距中心距离；
在PET重建中找到一个x’的分布，使得P(y|x’)最大，x’是一个变量，其维度与x相同，似然函数为：
利用对数似然函数极值点不变的特性，将公式(2)和公式(3)结合成l(y|x)＝logL(y|x)，其中：
使用最大期望方法求公式(4)的极值点，在不考虑衰减和散射的情况下：
其中N是均一化校正矩阵，i是迭代次数，T表示矩阵转置，d是随机校正矩阵,xi为估计图像，xi+1是本次迭代的结果，也是下一次迭代的输入。
进一步地，进行迭代前先利用GPU计算pTN，把均一化校正系数做反投影，得到敏感度图像sensitivity map。
进一步地，在正投影过程中，使用三维纹理内存对图像矩阵进行绑定；在反投影过程中，用原子操作来确保计算的正确性：原子操作执行时，锁定该地址，其他并行线程不允许对该地址对应的元素进行读取和写入操作，只能等待当前计算结束才能开始下一次访问。
进一步地，方法步骤的计算过程均在GPU中并行处理，将核函数在每个线程中同时并行计算；在线计算LOR穿过体素深度时，根据系统几何结构的对称性，将计算结果重复多次使用。
进一步地，迭代重建过程中，在进行投影之前及反投影得到校正图像后先对估计图像进行卷积操作，具体是将每个估计图像体素值与该体素值附近的等效系统函数相乘后相加，得到新的该体素位置图像值。
进一步地，测量点源在系统中重建图像上的系统函数，用这些物理效应对点源图像的影响来等效；等效系统响应函数H(x,y,z)是一个类三维高斯函数：
其中，r,t,a分别代表坐标轴的径向、切向和轴向，传播参数sigma在拟合中被估计，得到三个参数σr、σt、σa的值后，这三个参数代表系统响应函数的模型，σr、σt、σa分别表示径向、切向和轴向的传播参数。
根据本发明的另一实施例，提供了一种基于GPU加速的高精度PET重建装置，包括：
正投影单元，用于输入初始的估计图像，对估计图像进行正投影操作，获得估计投影；
比较单元，用于将获得的估计投影与实际测量投影作比较，获得比较结果；
反投影单元，用于对比较结果进行反投影操作，获得校正图像；
迭代单元，用于将校正图像与原始图像对应体素相乘，然后除敏感度图像，完成一次迭代，得到的图像就是下一次迭代的输入估计图像。
本发明实施例中的基于GPU加速的高精度PET重建方法及装置，输入初始的估计图像，对估计图像进行正投影操作，获得估计投影；将获得的估计投影与实际测量投影作比较，获得比较结果；对比较结果进行反投影操作，获得校正图像；将校正图像与原始图像对应体素相乘，然后除敏感度图像，完成一次迭代，得到的图像就是下一次迭代的输入估计图像。通过GPU加速策略将PET迭代重建算法改进后进行快速计算，获得高分辨率、高精度的三维PET图像。
附图说明
此处所说明的附图用来提供对本发明的进一步理解，构成本申请的一部分，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：
图1为本发明基于GPU加速的高精度PET重建方法的流程图；
图2为本发明基于GPU加速的高精度PET重建方法的优选流程图；
图3为本发明基于GPU加速的高精度PET重建方法中LOR的轴向平移特性图；
图4为本发明基于GPU加速的高精度PET重建方法中LOR的轴向对称特性图；
图5为本发明基于GPU加速的高精度PET重建方法的整体流程图；
图6为本发明基于GPU加速的高精度PET重建装置的模块图；
图7为本发明基于GPU加速的高精度PET重建装置的优选模块图。
具体实施方式
为了使本技术领域的人员更好地理解本发明方案，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分的实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都应当属于本发明保护的范围。
需要说明的是，本发明的说明书和权利要求书及上述附图中的术语“第一”、“第二”等是用于区别类似的对象，而不必用于描述特定的顺序或先后次序。应该理解这样使用的数据在适当情况下可以互换，以便这里描述的本发明的实施例能够以除了在这里图示或描述的那些以外的顺序实施。此外，术语“包括”和“具有”以及他们的任何变形，意图在于覆盖不排他的包含，例如，包含了一系列步骤或单元的过程、方法、系统、产品或设备不必限于清楚地列出的那些步骤或单元，而是可包括没有清楚地列出的或对于这些过程、方法、产品或设备固有的其它步骤或单元。
实施例1
根据本发明一实施例，提供了一种基于GPU加速的高精度PET重建方法，参见图1，包括以下步骤：
S100:输入初始的估计图像，对估计图像进行正投影操作，获得估计投影；
S200:将获得的估计投影与实际测量投影作比较，获得比较结果；
S300:对比较结果进行反投影操作，获得校正图像；
S400:将校正图像与原始图像对应体素相乘，然后除敏感度图像，完成一次迭代，得到的图像就是下一次迭代的输入估计图像。
本发明实施例中的基于GPU加速的高精度PET重建方法，输入初始的估计图像，对估计图像进行正投影操作，获得估计投影；将获得的估计投影与实际测量投影作比较，获得比较结果；对比较结果进行反投影操作，获得校正图像；将校正图像与原始图像对应体素相乘，然后除敏感度图像，完成一次迭代，得到的图像就是下一次迭代的输入估计图像。通过GPU加速策略将PET迭代重建算法改进后进行快速计算，获得高分辨率、高精度的三维PET图像。
具体地，参见图2，方法还包括步骤：
S500:将下一次迭代的输入估计图像返回步骤S100进行迭代处理，最后一次迭代获得的估计图像即是最终PET重建的结果。
具体地，参见图2，方法还包括步骤：
S001:对初始的估计图像进行预处理，将初始的估计图像进行统计和合并，生成符合预设要求的四维弦图。
下面以具体的实施例，对本发明基于GPU加速的高精度PET重建方法进行详细描述。
本发明提出的基于GPU加速的高精度PET重建方法分为预处理和迭代重建，预处理过程为把探测器得到的所有事件进行统计和合并，生成符合要求的四维弦图(Sinogram)，该四维弦图的每一个元素代表一类响应线(Line Of Response)，四维分别是LOR的平面内距中心距离(d)、平面内角度(phi)、接收第一个光子的探测器环编号(ring1)、接收第二个光子的探测器环编号(ring2)。
迭代重建的原理是PET数据y符合泊松分布的变量，是由系统收集到的所有LOR组成，每一条LOR代表着经过它的源产生的信号的总和。本发明需要估计出图像中每个体素内源的强度x，构建一个期望是x的泊松分布A，其数学模型为：
x和y的期望的关系是：
其中,k是图像中每个体素内源的出现次数，e是自然常数，y是实际测量的投影值，矩阵p是与探测器几何结构相关的表示在某一位置被探测器探测到的可能性，d是LOR的平面内距中心距离。在PET重建中本发明需要找到一个x’的分布，x’是一个变量，其维度与x相同，使得P(y|x’)最大，似然函数为：
利用对数似然函数极值点不变的特性，可以将公式(2)和公式(3)结合成l(y|x)＝logL(y|x)，其中：
使用最大期望方法可以求公式(4)的极值点，通常用于对包含隐变量或缺失数据的概率模型进行参数估计，是通过迭代进行极大似然估计的优化算法，算法的收敛性可以确保迭代至少逼近局部极大值。在不考虑衰减和散射的情况下：
其中N是均一化校正矩阵，i是迭代次数，T表示矩阵转置，d是随机校正矩阵,xi为估计图像，xi+1是本次迭代的结果，也是下一次迭代的输入。
计算机中操作的时候，输入初始估计图像，一般每个图像体素值都为1，在进行迭代前先利用GPU计算pTN，即把均一化校正系数做反投影，得到敏感度图像(sensitivitymap)。按公式先对估计图像xi进行正投影操作，得到的估计投影与实际测量投影作比，得到的结果进行反投影操作，然后把得到的图像与原始图像对应体素相乘，然后除sensitivitymap，完成一次迭代，得到的图像就是下一次迭代的输入估计图像，所有过程均在GPU中计算，GPU并行的原理是将核函数在每个线程中同时并行计算。最后一次迭代获得的估计图像就是最终重建的结果。对于重建中所用的函数，由于不需要线程间进行过多的通信，可以充分发挥GPU浮点计算速度快的特点。
在正投影过程中，如果每种LOR分配一个线程将有大于3*109个线程，如此多个线程进行并行，必然会有同时访问一个图像体素内存情况发生，造成数据竞争的问题，降低运行速度，对此本发明使用三维纹理内存对图像矩阵进行绑定。纹理内存是CUDA全局内存的特殊形式，对其读写操作通过专门的纹理缓存进行，可以让多个线程同时进行访问，且访问速度很快。
在反投影过程中，由于多个线程会同时修改同一个图像体素值，因此本发明使用原子操作来确保计算的正确性。原子操作执行时，锁定该地址，其他并行线程不允许对该地址对应的元素进行读取和写入操作，只能等待当前计算结束才能开始下一次访问，确保了安全。
GPU同时并行的流处理器为1000到3000左右，意味着实际上GPU只能同时计算这么多的线程，因此对于含有DOI信息的高分辨率PET系统，不可能同时并行运算所有分配的线程，线程实际上是一批批的进行计算，而这个过程是串行的。本发明的方法是利用系统的几何结构的对称性，可以减小需要计算的总LOR数量，从而减少并行线程数量，加快计算时间。对称性的考虑是在投影和反投影过程中，即在线计算LOR穿过体素深度时，根据系统的几何结构特性，可以将计算结果重复多次使用。有三类对称可以使用，第一类对称称为LOR的轴向平移特性，是指一条LOR可以在轴向平移，在相应的体素中深度不变，如图3所示，将LOR轴向平移，经过对应的体素深度不变，d1＝d2＝d3,因此可以将这些LOR归为一类进行计算；第二类对称称为LOR的轴向对称特性，是指LOR在轴向上中心对称，因此两条对称的LOR经过对称的像素深度完全相同，如图4所示，LOR轴向中心对称的另一条LOR，经过对应的体素深度不变，d1＝d2,因此可以将这2条LOR归为一类进行计算；第三类对称称为LOR的径向对称特性，因为LOR在径向上根据r和phi中心对称，因此可以将4条LOR归为一类，该对称可以类比第二类对称。
本发明提出的PET重建方法除了对均一化和随机事件进行校正以外，还加入了对系统正电子自由程和光子非线性的物理效应的校正，方法是测量点源在系统中重建图像上的系统函数，用这些物理效应对点源图像的影响来等效。等效系统响应函数H(x,y,z)是一个类三维高斯函数：
其中，r,t,a分别代表坐标轴的径向、切向和轴向，传播参数sigma在拟合中被估计，得到三个参数σr、σt、σa的值后，这三个参数就能代表系统响应函数的模型，σr、σt、σa分别表示径向、切向和轴向的传播参数。迭代重建过程中，在进行投影之前先要对估计图像进行卷积操作，具体是将每个估计图像体素值与该体素值附近的等效系统响应函数相乘后相加，得到新的该体素位置图像值。在反投影得到校正图像后，也要对其类似操作。本发明采集了系统多个不同位置的点源，选择了其中有代表性的图像进行拟合获得了系统响应函数。重建方法整体流程图如图5所示。
实施例2
根据本发明的另一实施例，提供了一种基于GPU加速的高精度PET重建装置，参见图6，包括：
正投影单元201，用于输入初始的估计图像，对估计图像进行正投影操作，获得估计投影；
比较单元202，用于将获得的估计投影与实际测量投影作比较，获得比较结果；
反投影单元203，用于对比较结果进行反投影操作，获得校正图像；
迭代单元204，用于将校正图像与原始图像对应体素相乘，然后除敏感度图像，完成一次迭代，得到的图像就是下一次迭代的输入估计图像。
本发明实施例中的基于GPU加速的高精度PET重建装置，输入初始的估计图像，对估计图像进行正投影操作，获得估计投影；将获得的估计投影与实际测量投影作比较，获得比较结果；对比较结果进行反投影操作，获得校正图像；将校正图像与原始图像对应体素相乘，然后除敏感度图像，完成一次迭代，得到的图像就是下一次迭代的输入估计图像。通过GPU加速策略将PET迭代重建算法改进后进行快速计算，获得高分辨率、高精度的三维PET图像。
具体地，参见图6，装置还包括：
返回迭代单元205，用于将下一次迭代的输入估计图像返回正投影单元201进行迭代处理，最后一次迭代获得的估计图像即是最终PET重建的结果。
具体地，参见图6，装置还包括：
预处理单元200，用于对初始的估计图像进行预处理，将初始的估计图像进行统计和合并，生成符合预设要求的四维弦图。
下面以具体的实施例，对本发明基于GPU加速的高精度PET重建装置进行详细描述。
本发明提出的基于GPU加速的高精度PET重建装置分为预处理和迭代重建，预处理单元200：预处理过程为把探测器得到的所有事件进行统计和合并，生成符合要求的四维弦图(Sinogram)，该四维弦图的每一个元素代表一类响应线(Line Of Response)，四维分别是LOR的平面内距中心距离(d)、平面内角度(phi)、接收第一个光子的探测器环编号(ring1)、接收第二个光子的探测器环编号(ring2)。
迭代重建的原理是PET数据y符合泊松分布的变量，是由系统收集到的所有LOR组成，每一条LOR代表着经过它的源产生的信号的总和。本发明需要估计出图像中每个体素内源的强度x，构建一个期望是x的泊松分布A，其数学模型为：
x和y的期望的关系是：
其中,k是图像中每个体素内源的出现次数，e是自然常数，y是实际测量的投影值，矩阵p是与探测器几何结构相关的表示在某一位置被探测器探测到的可能性，d是LOR的平面内距中心距离。在PET重建中本发明需要找到一个x’的分布，使得P(y|x’)最大，x’是一个变量，其维度与x相同，似然函数为：
利用对数似然函数极值点不变的特性，可以将公式(2)和公式(3)结合成l(y|x)＝logL(y|x)，其中：
使用最大期望方法可以求公式(4)的极值点，通常用于对包含隐变量或缺失数据的概率模型进行参数估计，是通过迭代进行极大似然估计的优化算法，算法的收敛性可以确保迭代至少逼近局部极大值。在不考虑衰减和散射的情况下：
其中N是均一化校正矩阵，i是迭代次数，T表示矩阵转置，d是随机校正矩阵,xi为估计图像，xi+1是本次迭代的结果，也是下一次迭代的输入。
计算机中操作的时候，输入初始估计图像，一般每个图像体素值都为1，在进行迭代前先利用GPU计算pTN，即把均一化校正系数做反投影，得到敏感度图像(sensitivitymap)。按公式先对估计图像xi进行正投影操作，得到的估计投影与实际测量投影作比，得到的结果进行反投影操作，然后把得到的图像与原始图像对应体素相乘，然后除sensitivitymap，完成一次迭代，得到的图像就是下一次迭代的输入估计图像，所有过程均在GPU中计算，GPU并行的原理是将核函数在每个线程中同时并行计算。最后一次迭代获得的估计图像就是最终重建的结果。对于重建中所用的函数，由于不需要线程间进行过多的通信，可以充分发挥GPU浮点计算速度快的特点。
在正投影过程中，如果每种LOR分配一个线程将有大于3*109个线程，如此多个线程进行并行，必然会有同时访问一个图像体素内存情况发生，造成数据竞争的问题，降低运行速度，对此本发明使用三维纹理内存对图像矩阵进行绑定。纹理内存是CUDA全局内存的特殊形式，对其读写操作通过专门的纹理缓存进行，可以让多个线程同时进行访问，且访问速度很快。
在反投影过程中，由于多个线程会同时修改同一个图像体素值，因此本发明使用原子操作来确保计算的正确性。原子操作执行时，锁定该地址，其他并行线程不允许对该地址对应的元素进行读取和写入操作，只能等待当前计算结束才能开始下一次访问，确保了安全。
GPU同时并行的流处理器为1000到3000左右，意味着实际上GPU只能同时计算这么多的线程，因此对于含有DOI信息的高分辨率PET系统，不可能同时并行运算所有分配的线程，线程实际上是一批批的进行计算，而这个过程是串行的。本发明的方法是利用系统的几何结构的对称性，可以减小需要计算的总LOR数量，从而减少并行线程数量，加快计算时间。对称性的考虑是在投影和反投影过程中，即在线计算LOR穿过体素深度时，根据系统的几何结构特性，可以将计算结果重复多次使用。有三类对称可以使用，第一类对称称为LOR的轴向平移特性，是指一条LOR可以在轴向平移，在相应的体素中深度不变，如图3所示，将LOR轴向平移，经过对应的体素深度不变，d1＝d2＝d3,因此可以将这些LOR归为一类进行计算；第二类对称称为LOR的轴向对称特性，是指LOR在轴向上中心对称，因此两条对称的LOR经过对称的像素深度完全相同，如图4所示，LOR轴向中心对称的另一条LOR，经过对应的体素深度不变，d1＝d2,因此可以将这2条LOR归为一类进行计算；第三类对称称为LOR的径向对称特性，因为LOR在径向上根据r和phi中心对称，因此可以将4条LOR归为一类，该对称可以类比第二类对称。
本发明提出的PET重建方法除了对均一化和随机事件进行校正以外，还加入了对系统正电子自由程和光子非线性的物理效应的校正，方法是测量点源在系统中重建图像上的系统函数，用这些物理效应对点源图像的影响来等效。等效系统响应函数H(x,y,z)是一个类三维高斯函数：
其中，r,t,a分别代表坐标轴的径向、切向和轴向，传播参数sigma在拟合中被估计，得到三个参数σr、σt、σa的值后，这三个参数就能代表系统响应函数的模型，σr、σt、σa分别表示径向、切向和轴向的传播参数。迭代重建过程中，在进行投影之前先要对估计图像进行卷积操作，具体是将每个估计图像体素值与该体素值附近的等效系统响应函数相乘后相加，得到新的该体素位置图像值。在反投影得到校正图像后，也要对其类似操作。本发明采集了系统多个不同位置的点源，选择了其中有代表性的图像进行拟合获得了系统响应函数。重建方法整体流程图如图5所示。
本发明将迭代重建算法使用在GPU计算中，并考虑了系统的几何对称性和CUDA C纹理内存加速方法，提出一种速度快、精度高的基于GPU加速的高精度PET重建方法，通过CUDA C并行加速和探测器几何特性将超大规模数据量的重建工作并行处理，进行高精度快速重建，即通过GPU加速策略将PET迭代重建算法改进后进行快速计算，获得高分辨率、高精度的三维PET图像。
在实际应用中，只需将PET系统采集到的数据输入即可获得重建图像。提高重建速度，提高图像分辨率。该方法及装置的保护点在于：基于GPU的利用几何对称性加速高清迭代重建方法；基于实际采集数据的系统响应函数计算方法。
上述本发明实施例序号仅仅为了描述，不代表实施例的优劣。
在本发明的上述实施例中，对各个实施例的描述都各有侧重，某个实施例中没有详述的部分，可以参见其他实施例的相关描述。
在本申请所提供的几个实施例中，应该理解到，所揭露的技术内容，可通过其它的方式实现。其中，以上所描述的系统实施例仅仅是示意性的，例如单元的划分，可以为一种逻辑功能划分，实际实现时可以有另外的划分方式，例如多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另一点，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些接口，单元或模块的间接耦合或通信连接，可以是电性或其它的形式。
作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本实施例方案的目的。
另外，在本发明各个实施例中的各功能单元可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个单元中。上述集成的单元既可以采用硬件的形式实现，也可以采用软件功能单元的形式实现。
集成的单元如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的全部或部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可为个人计算机、服务器或者网络设备等)执行本发明各个实施例方法的全部或部分步骤。而前述的存储介质包括：U盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、移动硬盘、磁碟或者光盘等各种可以存储程序代码的介质。
以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。
Claims (5)
1.一种基于GPU加速的高精度PET重建方法，其特征在于，包括以下步骤：
S001:对初始的估计图像进行预处理，将初始的估计图像进行统计和合并，生成符合预设要求的四维弦图Sinogram，该四维弦图的每一个元素代表一类响应线LOR，四维分别是LOR的平面内距中心距离d、平面内角度phi、接收第一个光子的探测器环编号ring1、接收第二个光子的探测器环编号ring2；
S100:输入初始的估计图像，对估计图像进行正投影操作，获得估计投影；
S200:将获得的估计投影与实际测量投影作比较，获得比较结果；
S300:对比较结果进行反投影操作，获得校正图像；
S400:将校正图像与估计图像对应体素相乘，然后除敏感度图像，完成一次迭代，得到的图像就是下一次迭代的输入估计图像；
所述方法还包括步骤：
S500:将下一次迭代的输入估计图像返回步骤S100进行迭代处理，最后一次迭代获得的估计图像即是最终PET重建的结果；
在迭代过程中估计出估计图像中每个体素内源的强度值x，构建一个期望是强度值x的泊松分布A，其数学模型为：
x和y的期望的关系是：
y＝px+d (2)
其中,k是估计图像中每个体素内源的出现次数，e是自然常数，y是实际测量的投影值，矩阵p是与探测器几何结构相关的表示在某一位置被探测器探测到的可能性，d是LOR的平面内距中心距离；
在PET重建中找到一个x’的分布，使得P(y|x’)最大，x’是一个变量，其维度与x相同，似然函数为：
利用对数似然函数极值点不变的特性，将公式(2)和公式(3)结合成l(y|x)＝logL(y|x)，其中：
l(y|x)＝∑y log y-y-log y！ (4)
使用最大期望方法求公式(4)的极值点，在不考虑衰减和散射的情况下：
其中N是均一化校正矩阵，i是迭代次数，T表示矩阵转置，m是随机校正矩阵,xi为估计图像，xi+1是本次迭代的结果，也是下一次迭代的输入；
进行迭代前先利用GPU计算pTN，把均一化校正系数做反投影，得到敏感度图像。
2.根据权利要求1所述的基于GPU加速的高精度PET重建方法，其特征在于，所述方法步骤的计算过程均在GPU中并行处理，将核函数在每个线程中同时并行计算；在线计算LOR穿过体素深度时，根据系统几何结构的对称性，将计算结果重复多次使用。
3.根据权利要求2所述的基于GPU加速的高精度PET重建方法，其特征在于，迭代重建过程中，在进行投影之前及反投影得到校正图像后先对估计图像进行卷积操作，具体是将每个估计图像体素值与该体素值附近的等效系统函数相乘后相加，得到新的该体素位置图像值。
4.根据权利要求3所述的基于GPU加速的高精度PET重建方法，其特征在于，测量点源在系统中重建图像上的系统函数，用这些物理效应对点源图像的影响来等效；等效系统响应函数H(x,y,z)是一个类三维高斯函数：
其中，r,t,a分别代表坐标轴的径向、切向和轴向，传播参数σ在拟合中被估计，得到三个参数σr、σt、σa的值后，这三个参数代表系统响应函数的模型，σr、σt、σa分别表示径向、切向和轴向的传播参数。
5.一种基于GPU加速的高精度PET重建装置，其特征在于，包括：
预处理单元，对初始的估计图像进行预处理，将初始的估计图像进行统计和合并，生成符合预设要求的四维弦图Sinogram，该四维弦图的每一个元素代表一类响应线LOR，四维分别是LOR的平面内距中心距离d、平面内角度phi、接收第一个光子的探测器环编号ring1、接收第二个光子的探测器环编号ring2；
正投影单元，用于输入初始的估计图像，对估计图像进行正投影操作，获得估计投影；
比较单元，用于将获得的估计投影与实际测量投影作比较，获得比较结果；
反投影单元，用于对比较结果进行反投影操作，获得校正图像；
迭代单元，用于将校正图像与原始图像对应体素相乘，然后除敏感度图像，完成一次迭代，得到的图像就是下一次迭代的输入估计图像；
返回迭代单元，用于将下一次迭代的输入估计图像返回正投影单元进行迭代处理，最后一次迭代获得的估计图像即是最终PET重建的结果；
在迭代过程中估计出估计图像中每个体素内源的强度值x，构建一个期望是强度值x的泊松分布A，其数学模型为：
x和y的期望的关系是：
y＝px+d (2)
其中,k是估计图像中每个体素内源的出现次数，e是自然常数，y是实际测量的投影值，矩阵p是与探测器几何结构相关的表示在某一位置被探测器探测到的可能性，d是LOR的平面内距中心距离；
在PET重建中找到一个x’的分布，使得P(y|x’)最大，x’是一个变量，其维度与x相同，似然函数为：
利用对数似然函数极值点不变的特性，将公式(2)和公式(3)结合成l(y|x)＝logL(y|x)，其中：
l(y|x)＝∑y log y-y-log y！ (4)
使用最大期望方法求公式(4)的极值点，在不考虑衰减和散射的情况下：
其中N是均一化校正矩阵，i是迭代次数，T表示矩阵转置，m是随机校正矩阵,xi为估计图像，xi+1是本次迭代的结果，也是下一次迭代的输入；
进行迭代前先利用GPU计算pTN，把均一化校正系数做反投影，得到敏感度图像。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201911286593.1A CN110811667B (zh) | 2019-12-14 | 2019-12-14 | 一种基于gpu加速的高精度pet重建方法及装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201911286593.1A CN110811667B (zh) | 2019-12-14 | 2019-12-14 | 一种基于gpu加速的高精度pet重建方法及装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110811667A CN110811667A (zh) | 2020-02-21 |
| CN110811667B true CN110811667B (zh) | 2023-08-15 |
Family
ID=69545473
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201911286593.1A Active CN110811667B (zh) | 2019-12-14 | 2019-12-14 | 一种基于gpu加速的高精度pet重建方法及装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110811667B (zh) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111754436B (zh) * | 2020-06-24 | 2024-05-03 | 上海联影医疗科技股份有限公司 | 医学图像伪影校正的加速方法、计算机设备和存储介质 |
| CN111881412A (zh) * | 2020-07-28 | 2020-11-03 | 南京航空航天大学 | 一种基于cuda的pet系统矩阵计算方法 |
| CN112991482B (zh) * | 2021-04-12 | 2023-03-24 | 明峰医疗系统股份有限公司 | 基于gpu的快速重建成像方法、设备及可读存储介质 |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101283913A (zh) * | 2008-05-30 | 2008-10-15 | 首都师范大学 | Ct图像重建的gpu加速方法 |
| CN103784158A (zh) * | 2012-10-29 | 2014-05-14 | 株式会社日立医疗器械 | Ct 装置及ct 图像生成方法 |
| CN105931280A (zh) * | 2016-03-29 | 2016-09-07 | 中北大学 | 基于gpu的快速三维ct迭代重建方法 |
| CN106255994A (zh) * | 2014-02-18 | 2016-12-21 | 皇家飞利浦有限公司 | 针对正电子发射断层摄影(pet)列表模式迭代重建的重建中滤波 |
| CN107220924A (zh) * | 2017-04-11 | 2017-09-29 | 西安电子科技大学 | 一种基于gpu加速pet图像重建的方法 |
| CN109949411A (zh) * | 2019-03-22 | 2019-06-28 | 电子科技大学 | 一种基于三维加权滤波反投影和统计迭代的图像重建方法 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7693318B1 (en) * | 2004-01-12 | 2010-04-06 | Pme Ip Australia Pty Ltd | Method and apparatus for reconstruction of 3D image volumes from projection images |
| US9196063B2 (en) * | 2013-05-02 | 2015-11-24 | Korea Advanced Institute Of Science And Technology | Super-resolution apparatus and method |
-
2019
- 2019-12-14 CN CN201911286593.1A patent/CN110811667B/zh active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101283913A (zh) * | 2008-05-30 | 2008-10-15 | 首都师范大学 | Ct图像重建的gpu加速方法 |
| CN103784158A (zh) * | 2012-10-29 | 2014-05-14 | 株式会社日立医疗器械 | Ct 装置及ct 图像生成方法 |
| CN106255994A (zh) * | 2014-02-18 | 2016-12-21 | 皇家飞利浦有限公司 | 针对正电子发射断层摄影(pet)列表模式迭代重建的重建中滤波 |
| CN105931280A (zh) * | 2016-03-29 | 2016-09-07 | 中北大学 | 基于gpu的快速三维ct迭代重建方法 |
| CN107220924A (zh) * | 2017-04-11 | 2017-09-29 | 西安电子科技大学 | 一种基于gpu加速pet图像重建的方法 |
| CN109949411A (zh) * | 2019-03-22 | 2019-06-28 | 电子科技大学 | 一种基于三维加权滤波反投影和统计迭代的图像重建方法 |
Non-Patent Citations (1)
| Title |
|---|
| 窦少彬.CT图像优化重建算法研究.《中国博士学位论文全文数据库 信息科技辑》.2019,(第8期),第I138-42页. * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110811667A (zh) | 2020-02-21 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Zhou et al. | Fast and efficient fully 3D PET image reconstruction using sparse system matrix factorization with GPU acceleration | |
| Shackleford et al. | On developing B-spline registration algorithms for multi-core processors | |
| CN110811667B (zh) | 一种基于gpu加速的高精度pet重建方法及装置 | |
| Kachelriess et al. | Hyperfast parallel‐beam and cone‐beam backprojection using the cell general purpose hardware | |
| US8620054B2 (en) | Image reconstruction based on accelerated method using polar symmetries | |
| US9111381B2 (en) | Shift-varying line projection using graphics hardware | |
| Kösters et al. | EMRECON: An expectation maximization based image reconstruction framework for emission tomography data | |
| Cui et al. | Distributed MLEM: An iterative tomographic image reconstruction algorithm for distributed memory architectures | |
| EP2875489B1 (en) | Method, computer readable medium and system for tomographic reconstruction | |
| Keck et al. | GPU-accelerated SART reconstruction using the CUDA programming environment | |
| Liu et al. | GPU-based branchless distance-driven projection and backprojection | |
| CN102831627A (zh) | 一种基于gpu多核并行处理的pet图像重建方法 | |
| WO2021114306A1 (zh) | 一种基于gpu加速的高精度pet重建方法及装置 | |
| Nguyen et al. | GPU-accelerated 3D Bayesian image reconstruction from Compton scattered data | |
| Wettenhovi et al. | OMEGA—open-source emission tomography software | |
| Liu et al. | GPU-based acceleration for interior tomography | |
| US20060098857A1 (en) | Method for faster iterative reconstruction for converging collimation spect with depth dependent collimator response modeling | |
| Schellmann et al. | Parallel medical image reconstruction: from graphics processing units (GPU) to grids | |
| Xie et al. | An effective CUDA parallelization of projection in iterative tomography reconstruction | |
| CN113112558B (zh) | 一种高清pet图像重建方法 | |
| Selivanov et al. | Fast PET image reconstruction based on SVD decomposition of the system matrix | |
| KR20120127214A (ko) | Ｇｐｕ를 이용한 양전자 방출 단층 촬영 영상에서의 몬테카를로 시뮬레이션 감마선 산란 추정 방법 및 장치 | |
| Weinlich et al. | Comparison of high-speed ray casting on GPU using CUDA and OpenGL | |
| Pratx et al. | 3-D tomographic image reconstruction from randomly ordered lines with CUDA | |
| Nguyen et al. | GPU-accelerated iterative reconstruction from Compton scattered data using a matched pair of conic projector and backprojector |
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 |