程，是将LOR和像素联系起来的过程，系统矩阵中 各个元素体现了成像过程中特定像素被特定LOR探测到的概率，在迭代重建的过程中至关 重要。系统矩阵包括LOR与像素的耦合定位以及LOR与像素的耦合程度。耦合定位反映的 是像素发出的湮灭光子是否会对LOR产生贡献，是重建中正投影和反投影的基础，直接影 响到重建图像的空间分辨率。而耦合程度则反映的是产生贡献的多少，影响到图像重建中 的对比度恢复和信噪比，也是影响图像重建的一个重要因素。
[0050] 在GPU加速迭代重建的方案中，本实例具体选用解析重建的方法实时计算系统矩 阵，不予提前进行计算存储。解析计算法采用比较简单的模型，在准确描述了图像空间和探 测器结构的前提下，计算像素被LOR探测的情况，主要反映的是探测器的几何结构。通常 采用的计算模型有点线模型、线积分模型、面积分模型和立体角模型。这些模型在3D重建 的情况下，由于需要考虑LOR在Z方向上的扩展，在三维空间的射束计算将远比二维形式复 杂。
[0051] 三维迭代重建算法最为耗时的就是正投影（projection)和反投影 (backprojection)的计算过程，这两个过程都连接图像空间和投影空间。算法流程中估计 图像以及校正矩阵都属于图像空间，而投影图像对应了投影空间的维度。对于GPU加速重 建的考虑，正投影过程可以分解为对于某一条特定的LOR(Sinc)投影矩阵的一个元素）可能 发生耦合的图像矩阵元素的贡献累加，正投影过程对于每一条LOR而言都是独立的。同样， 反投影的过程可以分解为对于某一图像元素其可能发生耦合的LOR的贡献累加，对于反投 影的过程，每一个像素之间都是相互独立的。对于正反投影的示意图如图3a和图3b所示：
[0052] 对于正投影和反投影GPU加速的kernel函数设置，正投影算子Proj将时间帧内 PET数据中的探测器符合线分配给GPU线程网格的一条线程单元，即一共分配i个线程，对 应i个探测器符合线对，如图3a所示。正投影采用射束遍历（ray driven)的方式，将每一 条LOR的处理作为一条线程进行并行运算。系统矩阵的解析计算可以采用ray-driven方 式，该方法是由Siddon提出的一种经典算法。在ray-driven方式中，沿着LOR的方向遍历， 列举LOR穿过的所有离散的图像像素。在扩展到3D以后，射线仍然是一条理想的直线。但 是在分辨率较好的探测器系统中，LOR的宽度一般会大于像素大小，因此，需要对LOR进行 插值避免错误。或者使用Tube模型也是行之有效的方式，即在射束计算时，可以将LOR当 成tube，即采用TOR模型，而像素则可以当成理想的格点。沿着TOR采样时，被TOR穿过的 像素即认为被该TOR探测到，其探测概率可以由晶体条对像素所张立体角决定。
[0053] 反投影算子BackPro j将时间帧内PET数据中的每一像素分配给GPU线程网格的 一条线程单元，即一共分配j个线程，对应j个像素，如图3b所示。反投影则采用像素遍 历（voxel driven)的方式，对于每一个像素分配一条线程，将每一条像素（voxel)的处理 作为一条线程进行并行运算。Voxel-driven方法是对每一个像素点，求出所有经过该点的 L0R，计算像素对LOR的贡献。Voxel-Driven算法流程如下：1.选定一个像素点2.列举LOR 所有环差和列差3.求出某一环差和列差情况下LOR的方向向量4.根据方向向量确定一 条经过选定像素点的射线5.求出射线与两端探测器的交点，从而确定一条L0R，并计算像 素点对LOR的贡献6.用同样的方法计算其他所有环差和列差7.计算出所有可能LOR对该 voxel的贡献累计。
[0054] 本方法对于系统矩阵不进行纹理捆绑，虽然这样会增加总体的计算量，但对于GPU 加速方法来说，计算密度高的运算正是其处理强项。另外对于实时计算耦合定位与耦合概 率的方式对于所有的PET探测器更为普适。对于系统矩阵的构架以及GPU加速重建方法已 有很多研究结果，本实例中列举了一种可行方式，对于具体的GPU加速重建算法不做限定。
[0055] 对于整个算法流程，在正反重建过程使用GPU的kernel函数进行加速，此外可也 以选择性地对于投影比较和图像更新的环节进行GPU加速，这两部分计算耗时较小，但使 用GPU加速也可以带来微量的时间上的性能提升。判断第一时间帧的PET图像是否达到预 定标准，可以根据第一时间帧的PET图像是否达到药物的代谢变化为依据，如第一时间帧 的PET图像已达到药物的代谢变化，则停止PET图像的采集以及重建，否则继续进行穿刺操 作，进行下一帧的PET数据采集及PET图像重建，直到判断出PET图像已达到药物的代谢变 化为止。如判断出PET图像未达到上述预定标准，则需选取第二时间帧。本实施例中时间 帧的选取可以采用连续独占的模式、复用模式、间隔模式以及其他任意定义的模式。对于连 续独占的模式即如图2a所示，选取的第一时间帧和第二时间帧之间不存在时间间隔，且彼 此之间不存在重叠。对于复用模式即如图2b所示，选取的第一时间帧和第二时间帧之间有 部分重叠的部分。因为PET低计数率特性，复用模式可以对于数据进行复用，这样可以提 高PET图像质量，适合于灵敏度低的设备。对于间隔模式即如图2c所示，选取的第一时间 帧和第二时间帧之间存在时间间隔，间隔模式可以对交互时间内的数据进行去除，这样可 以减少数据伪影，使得图像更为准确，适合于灵敏度高的设备，可以保证数据的利用率。然 而，本发明中时间帧的选取模式并不局限于此，其可以根据实际需要进行任意调整，例如可 以采用连续独占的模式、复用模式、间隔模式中的任意组合。
[0056] 本发明基于GPU加速的PET图像重建方法与以往技术相比具有以技术效果：
[0057] 本发明基于GPU加速的动态PET图像重建方法，由于可以对某一特定时间帧进行 PET图像重建，因此提高了 PET图像重建速度并且可以监测不同时间的PET图像，如此使得 PET图像引导下的交互操作成为可能。本发明的动态PET图像重建方法比其他方法而言，能 够保证有效图像质量的同时能快速的进行图像更新，以完成PET的实时采集和重建。
【主权项】
1. 基于GPU加速的动态PET图像重建方法，包括如下步骤： 1) 选取第一时间帧； 2) 采集所述第一时间帧内的PET数据； 3) 对采集的所述第一时间帧内的PET数据进行GPU加速图像重建，得到第一时间帧的 PET图像； 4) 判断所述第一时间帧的PET图像是否达到预定标准，如果所述第一时间帧的PET图 像达到所述预定标准，则停止PET数据采集与图像重建； 5) 如果所述第一时间帧的PET图像未达到所述预定标准，则继续选取第二时间帧； 6) 采集所述第二时间帧内的PET数据； 7) 对采集的所述第二时间帧内的PET数据进行GPU加速图像重建，得到第二时间帧的 PET图像； 8) 判断所述第二时间帧的PET图像是否达到所述预定标准，根据判断结果决定重复步 骤5)-7)直到PET图像达到所述预定标准，停止PET数据采集与图像重建。2. 如权利要求1所述的动态PET图像重建方法，其特征在于，所述时间帧内的PET数据 包括PET采集记录下的多个事例，每一个事例对应一响应线（LOR)。3. 如权利要求1所述的动态PET图像重建方法，其特征在于，选取的所述第一时间帧和 所述第二时间帧相互连续、相互间隔或部分重叠。4. 如权利要求1所述的动态PET图像重建方法，其特征在于，选取的所述第一时间帧和 所述第二时间帧的长度与PET系统的灵敏度相关。5. 如权利要求1所述的动态PET图像重建方法，其特征在于，对采集的时间帧内的PET 数据进行GPU加速图像重建包括： 在GPU上采用二维线程网格，将时间帧内的PET数据图形迭代重建的过程在所述二维 线程网格上进行分配。6. 如权利要求5所述的动态PET图像重建方法，其特征在于，所述迭代重建的算法的公 式为：其中，n表示迭代次数，为n+1次迭代第j个图像像素的估计值，xf为n次迭代 第j个图像像素的估计值，7,为表示第i个探测器符合线的测量值，A表示系统矩阵，au表 示第j个像素被第i个探测器对探测到的概率，Proj表示正投影算子，BackProj表示反投 影算子。7. 如权利要求6所述的动态PET图像重建方法，其特征在于，将时间帧内的PET数据图 形迭代重建的过程抽象在所述二维线程网格上进行分配包括： 正投影算子Proj将时间帧内PET数据中的探测器符合线分配给GPU线程网格的一条 线程单元，即一共分配i个线程，对应i个探测器符合线对； 反投影算子BackProj将时间帧内PET数据中的每一像素分配给GPU线程网格的一条 线程单元，即一共分配j个线程，对应j个像素。
【专利摘要】本发明公开了一种基于GPU加速的动态PET图像重建方法，包括如下步骤：选取时间帧；采集时间帧内的PET数据；对采集的时间帧内的PET数据进行GPU加速图像重建，得到时间帧的PET图像；判断时间帧的PET图像是否达到预定标准，如果时间帧的PET图像达到预定标准，则停止PET数据采集与图像重建；如果时间帧的PET图像未达到预定标准，则继续上述步骤，直到PET采集后重建的图像达到预定标准为止。本发明的动态PET图像重建方法，以时间帧为单位进行PET数据的采集和PET图像的重建，同时利用GPU并行运算加速对PET图像的重建，能够保证有效图像质量的同时能快速地进行图像更新，能够满足PET图像实时采集与重建的需求。
【IPC分类】G06T17/00, A61B6/03
【公开号】CN105046744
【申请号】CN201510401616
【发明人】顾笑悦, 李琳, 魏龙, 尹鹏飞, 刘双全, 贠明凯, 樊馨, 王海鹏, 黄先超, 丰宝桐, 王培林, 卢贞瑞
【申请人】中国科学院高能物理研究所
【公开日】2015年11月11日
【申请日】2015年7月9日