CN105046744A - 基于gpu加速的pet图像重建方法 - Google Patents

Abstract

本发明公开了一种基于GPU加速的动态PET图像重建方法，包括如下步骤：选取时间帧；采集时间帧内的PET数据；对采集的时间帧内的PET数据进行GPU加速图像重建，得到时间帧的PET图像；判断时间帧的PET图像是否达到预定标准，如果时间帧的PET图像达到预定标准，则停止PET数据采集与图像重建；如果时间帧的PET图像未达到预定标准，则继续上述步骤，直到PET采集后重建的图像达到预定标准为止。本发明的动态PET图像重建方法，以时间帧为单位进行PET数据的采集和PET图像的重建，同时利用GPU并行运算加速对PET图像的重建，能够保证有效图像质量的同时能快速地进行图像更新，能够满足PET图像实时采集与重建的需求。

Description

基于GPU加速的PET图像重建方法

技术领域

本发明涉及核医学PET图像重建领域，尤其涉及一种基于GPU加速的PET图像的采集与重建方法。

背景技术

PET(positronemissiontomography，PET)全称为正电子发射计算机断层扫描，是核医学领域重要的临床检查影像技术。PET技术是利用代谢和受体的功能显像技术，在肿瘤学临床医学影像和癌扩散等方面有着大量重要的应用。

PET的探测原理是将能发射正电子的核素作为特异性示踪剂注入病人体内，示踪剂发射的正电子与人体内的电子湮没产生一对γ光子，这对光子被探测器记录并储存下来，通过图像重建就可以得到示踪剂在病人体内的分布情况，从而实现肿瘤的早期发现和定位。临床现在最为广泛应用的示踪剂有氟代脱氧葡萄糖18F-FDG，18F-FDG为葡萄糖类似物，18F-FDG进入体内后主要会被癌细胞以及对葡萄糖生理性代谢需求高的组织(脑、心脏)所摄取。PET对于18F的衰变产物г光子进行探测，重建后的图像可反映体内器官/组织的葡萄糖代谢水平。所以PET显像与CT显像的结构成像原理不同，PET得到的图像可以反映人体内的新陈代谢情况，而癌症病变的功能改变往往早于结构改变，因而PET较之于CT等结构成像设备对于癌症病变的发现可以更为提前。

PET可以反应机体代谢功能，如果在此基础之上能进一步对于机体代谢进行随时间变化的动态分析，或者可以利用PET图像实时采集和重建进行动态交互，都将对于临床诊断以及各种前沿研究有着重要的价值。但是PET数据采集的低计数率特性限制了PET的动态成像，并且通用的基于CPU的PET重建方法需要在数据采集完成后再进行图像重建并且需要一定的耗时，可以说在时间方面的迟滞性是PET图像实时采集重建的最大瓶颈。

在数据采集方面，早期的PET在探测器环之间加入铅或钨的隔板，仅允许同一个环之间的晶体条的事例进行符合，这种2D采集模式灵敏度低，采集耗时长。随着探测器性能和电子学处理速度的提高，探测器环间的隔板被取消，进而采用3D采集模式，可以大大提高探测灵敏度，缩短采集时间。但是采用3D采集模式，允许不同环层间的事例进行符合，使得3D重建的计算量急剧增多，3D采集在大大地提高了探测效率的同时，对于相应的三维重建过程提出了高强度的计算需求。因而要对于PET进行实时成像，关键的问题在于如何解决图像实时重建的耗时问题，单纯基于CPU的串行解决方案耗时严重，已无法达到快速重建的需求。

发明内容

针对现有技术存在的问题，本发明的目的之一在于提供一种能够进行PET图像快速实时重建的基于GPU加速的PET图像采集与动态重建的方法。

为实现上述目的，本发明的基于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数据采集与图像重建。

进一步，所述时间帧内的PET数据包括PET采集记录下的多个事例，每一个事例对应一响应线(LOR)。

进一步，选取的所述第一时间帧和所述第二时间帧相互连续、相互间隔或部分重叠。

进一步，选取的所述第一时间帧和所述第二时间帧的长度与PET系统的灵敏度相关。

进一步，对采集的时间帧内的PET数据进行GPU加速图像重建包括：

在GPU上采用二维线程网格，将时间帧内的PET数据图形迭代重建的过程在所述二维线程网格上进行分配。

进一步，所述迭代重建的算法的公式为：

$x_j^{(n+1)}=x_j^n\frac{{\Σ}_{i=1}^{n_d}{a}_{ij}\frac{y_i}{{\[{\mathrm{Ax}}^{\left(n\right)}\]}_i}}{{\Σ}_{i=1}^{n_d}{a}_{ij}},\mathrm{Pr}oj={\[{\mathrm{Ax}}^{\left(n\right)}\]}_i,Back\mathrm{Pr}oj=\frac{{\Σ}_{i=1}^{n_d}{a}_{ij}\frac{y_i}{\mathrm{Pr}oj}}{{\Σ}_{i=1}^{n_d}{a}_{ij}}$

其中，n表示迭代次数，为n+1次迭代第j个图像像素的估计值，为n次迭代第j个图像像素的估计值，y_i为表示第i个探测器符合线的测量值，A表示系统矩阵，a_ij表示第j个像素被第i个探测器对探测到的概率，Proj表示正投影算子，BackProj表示反投影算子。

进一步，将时间帧内的PET数据图形迭代重建的过程抽象在所述二维线程网格上进行分配包括：

正投影算子Proj将时间帧内PET数据中的探测器符合线分配给GPU线程网格的一条线程单元，即一共分配i个线程，对应i个探测器符合线对；

反投影算子BackProj将时间帧内PET数据中的每一像素分配给GPU线程网格的一条线程单元，即一共分配j个线程，对应j个像素。

本发明的动态PET图像重建方法，以时间帧为单位进行PET数据的采集和PET图像的重建，同时利用GPU并行运算加速对PET图像的重建，能够在保证图像质量的同时能快速的进行图像更新，能够满足PET图像快速实时重建的需求。

附图说明

图1a和1b为本发明的动态PET图像重建方法的一实施例的流程示意图；

图2a-图2c为本发明的动态PET图像重建方法的时间选取模式的示意图；

图3a和图3b为本发明GPU加速三维迭代重建方法中正投影和反投影的示意图。

具体实施方式

如图1a和图1b所示，本发明的基于GPU加速的动态PET图像重建方法，包括如下步骤：

步骤S11:选取第一时间帧；

步骤S21：采集所述第一时间帧内的PET数据，其中PET采集过程记录下的是一个个的符合事例(event)，每一个事例累计加和到对应的符合响应线上(LOR)。本实例定义的时间帧即为特定时间内的PET数据。

步骤S31：对采集的所述第一时间帧内的PET数据进行GPU加速图像重建，得到第一时间帧的PET图像；

步骤S41、S51：判断所述第一时间帧的PET图像是否达到预定标准，如果所述第一时间帧的PET图像达到所述预定标准，则停止PET数据采集与图像重建；

如果所述第一时间帧的PET图像未达到所述预定标准，则重复上述步骤S11-S41，此处重复上述步骤S11-S41由于仅是在选取、采集、PET图像重建的对象上略有不同，即另外选取第二时间帧、采集第二时间帧内的PET数据以及对采集的第二时间帧内的PET数据进行GPU加速图像重建，其中第二时间帧与第一时间帧是不同的时间帧。

由于选取时间帧的长度的确定与PET系统本身的灵敏度有关，因此为了进一步提高PET图像重建的效率，在此对选取时间帧做特别说明。本实施例中步骤S11选取时间帧的长度与PET系统本身的灵敏度呈反比关系，即PET的灵敏度越高时间帧的选取就可以缩短，本实例按照一般PET的灵敏度选取30s作为时间帧进行使用为例来说明，时间帧最短要大于GPU加速重建时间，对于一般的PET使用高性能的GPU进行加速后，GPU加速一般都能达到10s以内，所以一般时间帧序列在[重建时间，30s]之间，以保证时间帧的更新，以及在下一个时间帧(第二时间帧)采集完成时可以获得上一个时间帧(第一时间帧)的重建图像。然而，选取时间帧的长度并不局限于此，时间帧的长度选取可以根据实际需要另做调整。

本实例使用MLEM迭代重建算法，MLEM为最经典的统计迭代法，基于最大似然函数的数据统计模型。算法流程如下所示：

1)输入初始的估计图像

2)对于估计图像进行正投影得到估计投影

3)真实测量得到的测量投影与估计投影进行比较获得投影比值

4)将投影比值进行反投影得到校正矩阵

5)用校正矩阵更新图像

然后重复以上过程，在一次次的校正迭代过程中，就会得到较为精确的估计图像，直至满足迭代停止条件，则迭代终止得到重建图像。

根据迭代重建算法的一般流程和MLEM算法的似然函数，MLEM算法的迭代过程公式如下所示：

$x_j^{(n+1)}=x_j^n\frac{{\Σ}_{i=1}^{n_d}{a}_{ij}\frac{y_i}{{\[{\mathrm{Ax}}^{\left(n\right)}\]}_i}}{{\Σ}_{i=1}^{n_d}{a}_{ij}},\mathrm{Pr}oj={\[{\mathrm{Ax}}^{\left(n\right)}\]}_i,Back\mathrm{Pr}oj=\frac{{\Σ}_{i=1}^{n_d}{a}_{ij}\frac{y_i}{\mathrm{Pr}oj}}{{\Σ}_{i=1}^{n_d}{a}_{ij}}$

其中，n表示征迭代次数，为n+1次迭代第j个图像像素的估计值，为n次迭代第j个图像像素的估计值，y_i为表示第i个探测器符合线的测量值，A表示整个系统矩阵，a_ij为系统矩阵A的元素，且表示第j个像素被第i个探测器对探测到的概率。

其中Proj为正征投影算子，代表每个探测器对上的估计投影值的计算过程。然后与真实测量值y_i进行对比，获得投影比值。BackProj为反投影算子，代表对于投影比值进行反投影，得到每个像素校正值的过程。

本实例采用了3D的数据采集模式，3D-MLEM重建，对于系统矩阵采用解析方法进行实时计算。对于GPU加速重建过程，主要对于正投影和反投影进行了速度的优化。系统矩阵描述了探测器探测被检测对象的过程，是将LOR和像素联系起来的过程，系统矩阵中各个元素体现了成像过程中特定像素被特定LOR探测到的概率，在迭代重建的过程中至关重要。系统矩阵包括LOR与像素的耦合定位以及LOR与像素的耦合程度。耦合定位反映的是像素发出的湮灭光子是否会对LOR产生贡献，是重建中正投影和反投影的基础，直接影响到重建图像的空间分辨率。而耦合程度则反映的是产生贡献的多少，影响到图像重建中的对比度恢复和信噪比，也是影响图像重建的一个重要因素。

在GPU加速迭代重建的方案中，本实例具体选用解析重建的方法实时计算系统矩阵，不予提前进行计算存储。解析计算法采用比较简单的模型，在准确描述了图像空间和探测器结构的前提下，计算像素被LOR探测的情况，主要反映的是探测器的几何结构。通常采用的计算模型有点线模型、线积分模型、面积分模型和立体角模型。这些模型在3D重建的情况下，由于需要考虑LOR在Z方向上的扩展，在三维空间的射束计算将远比二维形式复杂。

三维迭代重建算法最为耗时的就是正投影(projection)和反投影(backprojection)的计算过程，这两个过程都连接图像空间和投影空间。算法流程中估计图像以及校正矩阵都属于图像空间，而投影图像对应了投影空间的维度。对于GPU加速重建的考虑，正投影过程可以分解为对于某一条特定的LOR(sino投影矩阵的一个元素)可能发生耦合的图像矩阵元素的贡献累加，正投影过程对于每一条LOR而言都是独立的。同样，反投影的过程可以分解为对于某一图像元素其可能发生耦合的LOR的贡献累加，对于反投影的过程，每一个像素之间都是相互独立的。对于正反投影的示意图如图3a和图3b所示：

对于正投影和反投影GPU加速的kernel函数设置，正投影算子Proj将时间帧内PET数据中的探测器符合线分配给GPU线程网格的一条线程单元，即一共分配i个线程，对应i个探测器符合线对，如图3a所示。正投影采用射束遍历(raydriven)的方式，将每一条LOR的处理作为一条线程进行并行运算。系统矩阵的解析计算可以采用ray-driven方式，该方法是由Siddon提出的一种经典算法。在ray-driven方式中，沿着LOR的方向遍历，列举LOR穿过的所有离散的图像像素。在扩展到3D以后，射线仍然是一条理想的直线。但是在分辨率较好的探测器系统中，LOR的宽度一般会大于像素大小，因此，需要对LOR进行插值避免错误。或者使用Tube模型也是行之有效的方式，即在射束计算时，可以将LOR当成tube，即采用TOR模型，而像素则可以当成理想的格点。沿着TOR采样时，被TOR穿过的像素即认为被该TOR探测到，其探测概率可以由晶体条对像素所张立体角决定。

反投影算子BackProj将时间帧内PET数据中的每一像素分配给GPU线程网格的一条线程单元，即一共分配j个线程，对应j个像素，如图3b所示。反投影则采用像素遍历(voxeldriven)的方式，对于每一个像素分配一条线程，将每一条像素(voxel)的处理作为一条线程进行并行运算。Voxel-driven方法是对每一个像素点，求出所有经过该点的LOR，计算像素对LOR的贡献。Voxel-Driven算法流程如下：1.选定一个像素点2.列举LOR所有环差和列差3.求出某一环差和列差情况下LOR的方向向量4.根据方向向量确定一条经过选定像素点的射线5.求出射线与两端探测器的交点，从而确定一条LOR，并计算像素点对LOR的贡献6.用同样的方法计算其他所有环差和列差7.计算出所有可能LOR对该voxel的贡献累计。

本方法对于系统矩阵不进行纹理捆绑，虽然这样会增加总体的计算量，但对于GPU加速方法来说，计算密度高的运算正是其处理强项。另外对于实时计算耦合定位与耦合概率的方式对于所有的PET探测器更为普适。对于系统矩阵的构架以及GPU加速重建方法已有很多研究结果，本实例中列举了一种可行方式，对于具体的GPU加速重建算法不做限定。

对于整个算法流程，在正反重建过程使用GPU的kernel函数进行加速，此外可也以选择性地对于投影比较和图像更新的环节进行GPU加速，这两部分计算耗时较小，但使用GPU加速也可以带来微量的时间上的性能提升。判断第一时间帧的PET图像是否达到预定标准，可以根据第一时间帧的PET图像是否达到药物的代谢变化为依据，如第一时间帧的PET图像已达到药物的代谢变化，则停止PET图像的采集以及重建，否则继续进行穿刺操作，进行下一帧的PET数据采集及PET图像重建，直到判断出PET图像已达到药物的代谢变化为止。如判断出PET图像未达到上述预定标准，则需选取第二时间帧。本实施例中时间帧的选取可以采用连续独占的模式、复用模式、间隔模式以及其他任意定义的模式。对于连续独占的模式即如图2a所示，选取的第一时间帧和第二时间帧之间不存在时间间隔，且彼此之间不存在重叠。对于复用模式即如图2b所示，选取的第一时间帧和第二时间帧之间有部分重叠的部分。因为PET低计数率特性，复用模式可以对于数据进行复用，这样可以提高PET图像质量，适合于灵敏度低的设备。对于间隔模式即如图2c所示，选取的第一时间帧和第二时间帧之间存在时间间隔，间隔模式可以对交互时间内的数据进行去除，这样可以减少数据伪影，使得图像更为准确，适合于灵敏度高的设备，可以保证数据的利用率。然而，本发明中时间帧的选取模式并不局限于此，其可以根据实际需要进行任意调整，例如可以采用连续独占的模式、复用模式、间隔模式中的任意组合。

本发明基于GPU加速的PET图像重建方法与以往技术相比具有以技术效果：

本发明基于GPU加速的动态PET图像重建方法，由于可以对某一特定时间帧进行PET图像重建，因此提高了PET图像重建速度并且可以监测不同时间的PET图像，如此使得PET图像引导下的交互操作成为可能。本发明的动态PET图像重建方法比其他方法而言，能够保证有效图像质量的同时能快速的进行图像更新，以完成PET的实时采集和重建。

Claims (7)

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图像重建方法，其特征在于，所述迭代重建的算法的公式为：

$x_j^{(n+1)}=x_j^n\frac{{\mathrm{\Σ}}_{i=1}^{n_d}{a}_{ij}\frac{y_i}{{\[{\mathrm{Ax}}^{\left(n\right)}\]}_i}}{{\mathrm{\Σ}}_{i=1}^{n_d}{a}_{ij}},$ Proj＝[Ax⁽ⁿ⁾]_i, $BackProj=\frac{{\mathrm{\Σ}}_{i=1}^{n_d}{a}_{ij}\frac{y_i}{\mathrm{Pr}oj}}{{\mathrm{\Σ}}_{i=1}^{n_d}{a}_{ij}}$

其中，n表示迭代次数，为n+1次迭代第j个图像像素的估计值，为n次迭代第j个图像像素的估计值，y_i为表示第i个探测器符合线的测量值，A表示系统矩阵，a_ij表示第j个像素被第i个探测器对探测到的概率，Proj表示正投影算子，BackProj表示反投影算子。

7.如权利要求6所述的动态PET图像重建方法，其特征在于，将时间帧内的PET数据图形迭代重建的过程抽象在所述二维线程网格上进行分配包括：

正投影算子Proj将时间帧内PET数据中的探测器符合线分配给GPU线程网格的一条线程单元，即一共分配i个线程，对应i个探测器符合线对；

反投影算子BackProj将时间帧内PET数据中的每一像素分配给GPU线程网格的一条线程单元，即一共分配j个线程，对应j个像素。