CN111881412A - 一种基于cuda的pet系统矩阵计算方法 - Google Patents

Abstract

本发明公开了一种基于CUDA的PET系统矩阵计算方法。首先，对探测器系统按固定角度进行划分，得到划分后各子集的响应线方程；然后，在CUDA中给每个角度分配grid，每个角度下的每条子响应线分配block，将重建图像的像素坐标数据传至每个block的共享内存，block中的thread快速并行计算各像素的概率权值；各grid依次串行执行，同一个grid中不同block快速并行执行；最后，将概率权值在角度内拼接得到系统矩阵子集，再将系统矩阵子集拼接得到系统矩阵。本发明在保证系统矩阵元素准确的前提下，提高了创建效率，解决了系统矩阵建立速度慢以及GPU显存资源有限而系统矩阵数据量大所导致的难题。

Description

一种基于CUDA的PET系统矩阵计算方法

技术领域

本发明属于射线检测成像领域，特别涉及了一种PET系统矩阵计算方法。

背景技术

正电子发射型计算机断层显像(Posintron Emission Tomography,PET)是临床上应用广泛的一种无创核医学影像诊断技术，它通过对注入活体的放射性示踪剂成像，从而提供活体的新陈代谢等功能信息。目前PET主要应用于医学诊断领域，在工业领域研究较少。而γ光子具有强穿透性和电中性，受材料、结构以及电磁场等因素影响较小，因而将PET技术应用于工业无损检测领域具有重要意义。PET系统的主要工作过程：

(1)探测器得到符合数据：探测器接收湮灭释放的光子对，并判断光子对是否属于同一湮灭事件，从而获得大量符合数据；

(2)数据重组：将符合数据转换成由距离和角度组成的正弦图；

(3)图像重建：根据正弦图和重建算法从而得到切片图。

图像重建过程中，有公式Y＝A*I，其中Y为向量化的投影(向量化的正弦图)，I为向量化的切片图，A为需要计算的系统矩阵，也可以称为概率矩阵，其在重建图像中属于一个已知量或是一个固定参数。系统矩阵描述了探测器对被检测对象的探测过程，将所有的图像像素和响应线(LOR)连接起来，是迭代图像重建的关键，是正弦图数据通向切片图数据的纽带。虽然系统矩阵在图像重建过程极其重要，但是运算量大导致其生成速度很慢，难以满足工业检测速度快的要求。系统矩阵的大小与PET探测器数量和重建切片图大小有关，由于数据量大，其生成速度一直是一个难点，单纯的在CPU上进行改进很难保证速度有大幅度提升。随着GPU的发展，出现了很多运用GPU对计算密集型任务加速的算法。

2006年，NVIDIA公司发布了CUDA(Compute Unified Devices Architecture，统一计算设备架构)体系架构，该构架是一种新的并行编程模型和指令集的通用计算构架，它是基于高速图像处理单元GPU的，在GPU上高速并发执行，极大的提高了程序算法的运行速度，该体系构架为GPU编程提供了一种全新的软硬件构架。在CUDA架构下GPU中程序执行的最小单元是thread(线程)，一个block(线程块)中包含数个thread，同一个block中的thread可进行快速同步，同时还可通过共享内存(share memory)进行数据交换，不同block中的thread之间无法直接进行同步与数据交换。同一个block中thread开辟数量上受限较大，当加速程序需要的thread数量较大时，通常会将block中开辟的thread个数固定，而开辟多个block，多个block的集合称为grid(网格)，不同的thread使用thread ID以及block ID加以区分。利用NVIDIA GPU的并行计算引擎，线程、线程块和网格可以比CPU更高效地解决比较复杂的计算难题。

CUDA编程模型是一个异构模型，需要CPU和GPU协同工作。在异构计算架构中，GPU与CPU通过PCIe总线连接在一起来协同工作，CPU及其内存所在位置称为主机端(host)，而GPU及其内存所在位置称为设备端(device)。CUDA程序中既包含Host程序，又包含Device程序，它们分别在CPU和GPU上运行，CPU负责对逻辑性事务进行处理，以及对串行化运算的控制，而GPU重点处理数据密集型的并行计算程序。

发明内容

为了解决上述背景技术提到的技术问题，本发明提出了一种基于CUDA的PET系统矩阵计算方法，在保证系统矩阵大小和元素不变的情况下，加快系统矩阵建立速度。

为了实现上述技术目的，本发明的技术方案为：

一种基于CUDA的PET系统矩阵计算方法，包括以下步骤：

(1)根据探测器晶体的数量，对探测器系统按固定角度进行划分，得到划分后各个子集的响应线方程；

(2)在CUDA中给每个角度分配一个grid，每个角度下的每条子响应线分配一个block，将重建图像的像素坐标数据传输至每个block的共享内存，block中的thread用来快速并行计算各像素的概率权值；各个grid依次串行执行，同一个grid中不同block快速并行执行；

(3)将步骤(2)所得的概率权值在角度内进行拼接得到对应的系统矩阵子集，再将系统矩阵子集进行拼接，得到最终的系统矩阵。

进一步地，步骤(1)的具体过程如下：

(101)设探测器总数为X，X个探测器以固定的角度φ组成环形结构，构成一个圆形探测视野，并按逆时针方向对探测器依次编号，将连线过探测器圆环视野圆心的两个探测器划分为一个探测器对，则一共划分为N个探测器对，N＝X/2；其中第s号探测器与第s+N号探测器组成第s个探测器对，s＝1,2,…,N；

(102)连接第p号探测器与第k号探测器所形成的圆弧的中点与第p’号探测器与第k’号探测器所形成的圆弧的中点，以此形成一条过探测视野圆心的直线，规定此直线为角度k所在直线，1≤k≤N，共形成N个角度，相邻角度所在直线相差角度φ；其中第p号探测器与第k号探测器为相邻的两个探测器，第p号探测器与第p’号探测器属于一个探测器对，第k号探测器与第k’号探测器属于一个探测器对；

(103)角度k所在直线为第k组所有响应线的中垂线，根据第k组所有响应线距离第k号探测器的距离，依次进行编号，距离最近的为第k组第1条子响应线，距离最远的为第k组第N条子响应线。

进一步地，步骤(2)的具体过程如下：

(201)所述重建图像大小为M×M，将重建图像分割为单位像素为1×1的方格，以M×M单位像素方格组成的正方形中心为原点建立直角坐标系，同时探测器的圆形探测视野内切于此正方形，以角度1所在直线为x轴，角度(N/2+1)所在直线为y轴，各像素的坐标为(x_i,y_j)，i,j＝1,2,...,M，并且按自左向右、自上向下的原则对重建图像进行向量化，按顺序将像素编号为1到m，其中m＝M²，并且将其从CPU内存传输进GPU内存；

(202)对于角度k，给其分配一个grid；给角度k下的N条子响应线分配N个block，子响应线1到子响应线N依次分配给block_k1到block_kN；第k组第n条响应线表达式为Ax+By+C＝0，其中，

B＝1，

即为角度k的值，n＝1,2,…,N；

(203)给每个block分配一定数量的thread，在每个block中申请共享内存用于存放步骤(201)中得到的m个像素的坐标，同一个block中的thread在并行计算像素中心到子响应线的距离时通过共享内存来快速提取坐标数据；利用第k组第n条响应线的block_kn中的thread分别计算m个像素到第n条子响应线的距离

然后计算像素(x_i,y_j)的概率权值a_k,n,i,j：

进一步地，步骤(3)的具体过程如下：

(301)对任一个角度k的权值进行组合，创建一个N×M²的全0二维矩阵A_k，矩阵元素用a_i,j表示；

(302)将第k组第n条响应线的block_kn中的thread计算的概率权值对矩阵A_k中的第n行元素值进行更新，n＝1,2,…,N，得到系统矩阵子集A_k；

(303)创建一个N²×M²全0二维矩阵A，即系统矩阵，利用系统矩阵子集A₁到A_N对矩阵A的第1到第N行元素进行更新，得到最终的系统矩阵。

采用上述技术方案带来的有益效果：

本发明利用系统矩阵子集之间相关度低、计算密集度高的特点，充分发挥CUDA并行计算的优势；首先按角度对响应线进行划分得到各个响应线组，结合CUDA的并行计算能力对各响应线组对应的权值并行计算，最后将各子集有效的拼接，能够实现加速系统矩阵的生成，有效提高PET重建的性能。

本发明提出了使用共享内存的方法，共享内存作为CUDA中block级的内存，比其他内存存取速度更快，几乎无损耗，创新性地提出将重建图像像素的坐标信息提前存入各block的共享内存中，以便thread在并行计算时能够更加快速高效。

本发明针对GPU内存资源有限而系统矩阵数据量又大的问题，提出一种新的利用角度的划分来求取系统矩阵子集的方法，相比于中国专利申请CN109498048A公开的一种加速正电子图像重建的系统矩阵生成与处理方法中利用的角度划分方法，将一对探测器晶体的连线作为一个角度划分的依据，而本方明将两块相邻探测器晶体之间的中点与这两块晶体相对的两块探测器晶体之间的中点连线作为一个角度划分的依据，能够更加准确地划分出N条子响应线，并得到响应线截距大小，从而得到更优的系统矩阵元素。

本发明根据探测器的数量和重建图像大小可以动态分配CUDA的grid数、block数、thread数，从而可以快速生成不同大小的系统矩阵。

本发明有效利用CPU和GPU的优点，在CPU上完成角度划分和矩阵拼接的逻辑串行操作，在GPU上完成计算密集的概率权值并行计算；在GPU上给每一个角度的响应线组分配一个grid，响应线组中每一条子响应线在grid的不同block里进行处理，不同的block之间是并行执行的，进一步利用thread来实现各个像素概率权值的快速并行计算，由此来提高系统矩阵的计算速度。

附图说明

图1为实施例的流程示意图；

图2为实施例中角度2部分响应线示意图；

图3为实施例中CUDA的线程组织结构图；

图4为实施例中响应线穿过像素点以及概率权值计算示意图；

图5为实施例中系统矩阵拼接示意图。

具体实施方式

以下将结合附图，对本发明的技术方案进行详细说明。

以具有12个探测器、重建图像大小为6×6的小型PET系统的系统矩阵为例对本发明进一步说明，不同数量的探测器和不同大小的重建图像的PET系统进行系统矩阵计算时原理与本例相同。

如图1所示，本实施例的基于CUDA的PET系统矩阵计算方法，步骤如下：

步骤1、角度划分：根据12个探测器将PET系统划分为6个角度，即6个响应线组，进一步根据探测器数量，每个响应线组由6条子响应线组成。

步骤2、并行计算：给6个角度的响应线组分别分配一个grid，每个grid分配6个block，分别对应响应线组的6条子响应线，不同grid之间串行执行，同一个grid中不同block之间快速并行执行。再给每个block分配thread，用来计算各像素的概率权值，不同的thread之间并行执行。

步骤3、矩阵拼接：首先将步骤2所得的概率权值在各个角度内进行拼接，得到对应的6个系统矩阵子集，再将6个系统矩阵子集进行拼接，得到最终的系统矩阵。

在本实施例中，上述步骤1采用如下优选方案实现：

步骤101、如图2所示，相邻探测器以固定的角度φ逆时针方向依次排列，形成一个圆形的探测视野，以探测视野的圆心为原点，1号探测器和12号探测器形成的圆弧的中心为点A，OA方向为x轴正方向，3号探测器和4号探测器形成的的圆弧的中心为点B，OB方向为y轴正方向，组成一个坐标系，其中

步骤102、1号探测器和2号探测器形成的圆弧的中心为点C，7号探测器和8号探测器形成的圆弧的中心为点D，角度2对应的直线为L₂，即过原点O的直线CD，对应的角度

直线L₂作为角度2所有子响应线的中垂线，其中1号探测器与2号探测器的连线构成子响应线1(图2中线段Ⅰ)，3号探测器与12号探测器的连线构成子响应线2(图2中线段Ⅱ)，依次排序，7号探测器与8号探测器构成子响应线6，子响应线之间互相平行，且斜率

子响应线1的截距

故角度2子响应线1的方程为

子响应线2的截距

故角度2的子响应线2方程为

以此类推，可以计算出各个角度共6×6条子响应线方程。

在本实施例中，上述步骤2采用如下优选方案实现：

步骤201：如图3所示本发明CUDA的线程结构，Host为CPU端，Device为GPU端。首先将像素坐标信息从CPU内存传输进GPU内存，然后CPU端启动一个核函数，对应在GPU端的grid中处理。本实施例共启动6个核函数，即在GPU端分配6个grid，分别用来处理6个响应线组，核函数之间依次串行执行。在每个grid中分配6个block，分别用来计算6条子响应线，并在每个block中申请共享内存。在本实施例中给每个block分配64个thread，用来计算6×6个单位像素的概率权值，各个thread之间并行执行。

步骤202：如图4所示，一个正方形重建图像由6×6的单位像素构成，圆形探测视野内接于此正方形重建图像。各单位像素由(x_i,y_j)表示，i,j＝1,2,...,36，并且将36个像素的坐标存进每个block的共享内存中，以便thread在计算时能够快速从共享内存中提取坐标数据。AA’，BB’为角度k的两个子响应线，S_k,AA',2,3代表(x₂,y₃)这个像素的中心(-0.5,1.5)到角度k的子响应线AA’的距离，S_k,BB',5,4代表(x₅,y₄)这个像素中心(0.5,-1.5)到角度k的子响应线BB’的距离。根据子响应线的方程Ax_i+By_j+C＝0，计算各像素中心到子响应线的距离

其中n＝1,2,...,6，代表子响应线1到6。角度1中，block₁₁中的thread用来并行计算s_1,1,1,1,s_1,1,1,2,…,s_1,1,6,6，并按下式计算概率权值：

由此得到一组1×6²的数组。在block₁₁运算的同时，block₁₂到block₁₆并行计算像素中心到对应子响应线的距离并判断其大小，最终角度1对应的Grid1将得到一组6×6²的数组，即系统矩阵子集A₁。以此类推，顺序执行Grid2到Grid6都将得到一组6×6²的数组，对应系统矩阵子集A₂到系统矩阵子集A₆。

在本实施例中，上述步骤2采用如下优选方案实现：

通过步骤202计算得到了系统矩阵子集A₁到A₆，创建一个36×36的全0系统矩阵A，如图5所示，系统矩阵子集A₁对系统矩阵A的第1行到第6行进行更新，系统矩阵子集A₂对系统矩阵A的第7行到第12行进行更新，以此类推，系统矩阵子集A₆对系统矩阵A的第31行到第36行进行更新。由此完成了重建图像大小为6×6、探测器个数为12的系统矩阵的建立。

实施例仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明保护范围之内。

Claims (4)

1.一种基于CUDA的PET系统矩阵计算方法，其特征在于，包括以下步骤：

(1)根据探测器晶体的数量，对探测器系统按固定角度进行划分，得到划分后各个子集的响应线方程；

(2)在CUDA中给每个角度分配一个grid，每个角度下的每条子响应线分配一个block，将重建图像的像素坐标数据传输至每个block的共享内存，block中的thread用来快速并行计算各像素的概率权值；各个grid依次串行执行，同一个grid中不同block快速并行执行；

(3)将步骤(2)所得的概率权值在角度内进行拼接得到对应的系统矩阵子集，再将系统矩阵子集进行拼接，得到最终的系统矩阵。

2.根据权利要求1所述基于CUDA的PET系统矩阵计算方法，其特征在于，步骤(1)的具体过程如下：

(101)设探测器总数为X，X个探测器以固定的角度φ组成环形结构，构成一个圆形探测视野，并按逆时针方向对探测器依次编号，将连线过探测器圆环视野圆心的两个探测器划分为一个探测器对，则一共划分为N个探测器对，N＝X/2；其中第s号探测器与第s+N号探测器组成第s个探测器对，s＝1,2,…,N；

(102)连接第p号探测器与第k号探测器所形成的圆弧的中点与第p’号探测器与第k’号探测器所形成的圆弧的中点，以此形成一条过探测视野圆心的直线，规定此直线为角度k所在直线，1≤k≤N，共形成N个角度，相邻角度所在直线相差角度φ；其中第p号探测器与第k号探测器为相邻的两个探测器，第p号探测器与第p’号探测器属于一个探测器对，第k号探测器与第k’号探测器属于一个探测器对；

(103)角度k所在直线为第k组所有响应线的中垂线，根据第k组所有响应线距离第k号探测器的距离，依次进行编号，距离最近的为第k组第1条子响应线，距离最远的为第k组第N条子响应线。

3.根据权利要求2所述基于CUDA的PET系统矩阵计算方法，其特征在于，步骤(2)的具体过程如下：

(201)所述重建图像大小为M×M，将重建图像分割为单位像素为1×1的方格，以M×M单位像素方格组成的正方形中心为原点建立直角坐标系，同时探测器的圆形探测视野内切于此正方形，以角度1所在直线为x轴，角度(N/2+1)所在直线为y轴，各像素的坐标为(x_i,y_j)，i,j＝1,2,...,M，并且按自左向右、自上向下的原则对重建图像进行向量化，按顺序将像素编号为1到m，其中m＝M²，并且将其从CPU内存传输进GPU内存；

(202)对于角度k，给其分配一个grid；给角度k下的N条子响应线分配N个block，子响应线1到子响应线N依次分配给block_k1到block_kN；第k组第n条响应线表达式为Ax+By+C＝0，其中，

B＝1，

即为角度k的值，n＝1,2,…,N；

(203)给每个block分配一定数量的thread，在每个block中申请共享内存用于存放步骤(201)中得到的m个像素的坐标，同一个block中的thread在并行计算像素中心到子响应线的距离时通过共享内存来快速提取坐标数据；利用第k组第n条响应线的block_kn中的thread分别计算m个像素到第n条子响应线的距离

然后计算像素(x_i,y_j)的概率权值a_k,n,i,j：

4.根据权利要求3所述基于CUDA的PET系统矩阵计算方法，其特征在于，步骤(3)的具体过程如下：

(301)对任一个角度k的权值进行组合，创建一个N×M²的全0二维矩阵A_k，矩阵元素用a_i,j表示；

(302)将第k组第n条响应线的block_kn中的thread计算的概率权值对矩阵A_k中的第n行元素值进行更新，n＝1,2,…,N，得到系统矩阵子集A_k；

(303)创建一个N²×M²全0二维矩阵A，即系统矩阵，利用系统矩阵子集A₁到A_N对矩阵A的第1到第N行元素进行更新，得到最终的系统矩阵。

Images