CN102831627A - 一种基于gpu多核并行处理的pet图像重建方法 - Google Patents

Abstract

本发明公开了一种基于GPU多核并行处理的PET图像重建方法，包括：(1)对体素进行粒子采样；(2)多核并行求取粒子权重；(3)对粒子进行重采样；(4)求取粒子浓度真值及其权重真值；(5)计算每一体素的体素值。本发明通过利用粒子滤波算法，将PET中噪声的数据模型定为泊松分布而非高斯分布，更符合PET扫描的真实情况，从而使得重建过程中对噪声的过滤和优化更加有效，得到的重建结果比ML-EM等传统重建方法更接近PET真实情况；同时利用GPU多核并行处理技术对重建过程进行加速，将原本需要依次执行的计算转换成同时并行计算的模式，大大地缩短了计算的时间，使得PET图像重建方法能够适合于临床应用要求。

Description

一种基于GPU多核并行处理的PET图像重建方法

技术领域

本发明属于计算机图像处理技术领域，具体涉及一种基于GPU(图形处理器)多核并行处理的PET图像重建方法。

背景技术

正电子断层成像PET(Positron Emission Tomography)，是一种核医学功能成像技术。与CT成像不同，PET为功能成像，能得到病人体内新陈代谢的情况，从而可以比CT、MRI更早地探测出病灶。

经过同位素标记的药物在病人体内扩散的过程中，同位素不断衰变而产生正电子(e⁺)。正电子与周围的负电子(e^-)发生湮灭后，将发出一对光子。这对光子能够穿透人体组织发射到体外，从而被光电探测器探测到。通过探测这对光子，就可以得到一条符合线(LOR)。将若干个探测器排列在病人周围，可以接收发自病人体内的光子，从而得到一系列符合线。这些符合线和体内同位素标记的药物的位置有一定的关系，依据这些符合线，就可以重建出病人体内药物的浓度分布图，即PET图像。

两个毫不相关的光子如果在同一时间窗内被一对探测器探测到，那探测系统将认为这是由一对正负电子湮灭产生的，将它们的连线记录为符合线，这样的记录是假事件，在PET成像中称为随机符合事件。由于以上原因以及病人体内的干扰，PET采集到的数据包含许多噪声，这些噪声直接影响重建后的图像质量。所以，重建算法的鲁棒性(Robust)决定了PET图像的质量，这也是为什么重建算法成为现今PET领域的研究热点。

目前PET重建的方法大致可以分为解析法、迭代法和状态空间法三种：

解析法是以滤波反投影算法(FBP)为代表的一种方法，该算法是基于傅立叶变换来实现由投影数据得到原始数据。但是该算法假设的理想模型并不十分符合PET扫描的真实情况，加之该算法对噪声的抑制能力差，导致滤波反投影算法得到的图像伪影很重，质量很差。

迭代法中以ML-EM(Maximum Likelihood-Expectation Maximum)和MAP(Maximum a Posteriori)算法较为典型。ML-EM算法根据最大似然估计原理，将病人体素(Voxel)浓度作为待估计值、建立似然函数、求得满足一定条件的极大值。MAP算法根据原始数据和PET的统计特性建立后验函数再求其最大值，从而得到体素浓度分布。迭代算法大部分依赖于数据的统计模型，除了统计特性，决定PET图像的还有生理先验、结构先验等信息，而这些信息在统计迭代算法中并不能有效地被利用，这也是ML-EM和MAP的巨大缺陷。

状态空间法能够将PET成像的过程模型化，使其物理过程得到数学的表达，通过状态空间，能很好地把PET扫描中的统计特性、生理特性和结构特性结合起来，从而达到改进重建质量的目的。目前现有的状态空间求解方法为卡尔曼滤波，该方法假设数据分布符合高斯分布，这与PET数据的统计特性不符，所以重建得到的PET图像效果不是很理想，这也是目前困扰着状态空间方法的问题之一。

由于PET扫面得到的数据量大，重建算法复杂，导致重建所花的时间过长，这与临床需求相违背；而传统利用CPU串行计算也会花费相当长的时间，且需要CPU具有超强的计算性能，计算成本很高，故并不能在商业应用领域普及。

发明内容

针对现有技术所存在的上述技术缺陷，本发明提供了一种基于GPU多核并行处理的PET图像重建方法，在保证PET图像重建质量的前提下，大大提升重建的速度。

一种基于GPU多核并行处理的PET图像重建方法，包括如下步骤：

(1)根据体素值的概率分布区间，利用线性平均采样法对PET图像中的每一体素进行粒子采样，获得体素的每一粒子浓度值；

所述的线性平均采样法是通过对体素值的概率分布区间进行(s-1)等分，得到s个节点值，每个节点值即对应为任一体素的每一粒子浓度值，s为粒子采样个数；

(2)根据采集获取到的符合计数向量，利用FBP法求得体素的先验浓度值；根据体素的先验浓度值和粒子浓度值，利用GPU多核并行处理技术通过粒子权重估计算法求得体素每一粒子浓度值对应的粒子权重值；

(3)利用GPU多核并行处理技术对所述的粒子浓度值及其粒子权重值进行重采样，得到重采样后的粒子浓度值及其粒子权重值；

(4)返回步骤(2)，并使重采样后的粒子浓度值作为步骤(2)中的粒子浓度值，依此循环直至重采样后的粒子浓度值收敛，将收敛时重采样后的粒子浓度值及其粒子权重值作为粒子浓度真值及其粒子权重真值；

(5)根据所述的粒子浓度真值及其粒子权重真值，利用GPU多核并行处理技术求取PET图像中每一体素的体素值。

所述的粒子权重估计算法基于以下方程表达式：

$P(j,i)=\frac{D(i,j)\×\frac{x\left[j\right]}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}x\left[j\right]}}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\{D(i,j)\×\frac{x\left[j\right]}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}x\left[j\right]}\}}$

$w_{\left(k\right)}\left[j\right]=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\left\{y\right[i]\×P(j,i)\}-\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\{D(i,j)\×x_{\left(k\right)}[j\left]\right\}$

其中：x[j]为PET图像中第j体素的先验浓度值，y[i]为符合计数向量中的第i符合计数；D(i，j)为系统矩阵中第i行第j列的元素值；x_(k)[j]为PET图像中第j体素的第k粒子浓度值，w_(k)[j]为x_(k)[j]对应的粒子权重值；i、j和k均为自然数，且1≤i≤m，1≤j≤n，1≤k≤s，m为符合计数向量的维度，n为PET图像的体素个数。

系统矩阵D为m×n维矩阵，其与符合计数向量y以及PET图像浓度分布向量x满足关系式y＝Dx+e，e为测量噪声，PET图像浓度分布向量x为n维向量，其每个元素值对应PET图像每个体素的体素值，系统矩阵D与探测光子的光电探测器的结构和性能参数有关。

所述的步骤(5)中，通过以下方程表达式求取PET图像中每一体素的体素值：

$x^{*}\left[j\right]=\underset{k=1}{\overset{s}{\mathrm{\Σ}}}\left\{{x_{\left(k\right)}}^{*}\right[j]\×{w_{\left(k\right)}}^{*}[j\left]\right\}$

其中：x(k)^*[j]为PET图像中第j体素的第k粒子浓度真值，w_(k) ^*[j]为x_(k) ^*[j]对应的粒子权重真值，x^*[j]为PET图像中第j体素的体素值。

所述的GPU多核并行处理技术的计算处理策略如下：

对于向量与标量的乘法或除法，将该计算量映射到个数与向量维数相等的多个内核中，每个内核完成对应向量元素与标量相乘或相除；

对于矩阵与向量的乘法，将该计算量映射到个数与矩阵维数相等的多个内核中，每个内核完成对应矩阵元素与其相对应的向量元素相乘，然后对所得到的矩阵的各行再求和；

对于矩阵与标量的乘法或除法，将该计算量映射到个数与矩阵维数相等的多个内核中，每个内核完成对应矩阵元素与标量相乘或相除；

对于其余公式的代入计算，则使一个内核对应负责一个公式不同数据的计算。

本发明通过利用粒子滤波算法，将PET中噪声的数据模型定为泊松分布而非高斯分布，更符合PET扫描的真实情况，从而使得重建过程中对噪声的过滤和优化更加有效，得到的重建结果比ML-EM、FBP等传统重建方法更接近PET真实情况；同时利用GPU多核并行处理技术对重建过程进行加速，将原本需要依次执行的计算转换成同时并行计算的模式，大大地缩短了计算的时间，使得PET图像重建方法能够适合于临床应用要求。

附图说明

图1为本发明PET图像重建方法的步骤流程示意图。

图2为Zubal-thorax-phantom数字体模的模型示意图。

图3(a)为关于Zubal-thorax-phantom数字体模采用本发明方法重建后的PET图像。

图3(b)为关于Zubal-thorax-phantom数字体模采用ML-EM方法重建后的PET图像。

具体实施方式

为了更为具体地描述本发明，下面结合附图及具体实施方式对本发明的PET图像重建方法进行详细说明。

如图1所示，一种基于GPU多核并行处理的PET图像重建方法，包括如下步骤：

(1)对体素进行粒子采样。

根据体素值的概率分布区间，利用线性平均采样法对PET图像中的每一体素进行粒子采样，获得体素的每一粒子浓度值；

线性平均采样法是通过对体素值的概率分布区间进行(s-1)等分，得到s个节点值，每个节点值即对应为任一体素的每一粒子浓度值，s为粒子采样个数；本实施方式中，s＝500。

(2)多核并行求取粒子权重。

首先，根据采集获取到的符合计数向量，利用FBP法求得体素的先验浓度值；

其次，根据PET成像原理，故PET图像的状态空间方程可表示如下：

$\left\{\begin{array}{c}y=\mathrm{Dx}+e\\ x_{t+1}={\mathrm{Ax}}_t\end{array}\right.$

其中：D为m×n维的系统矩阵，其与探测光子的光电探测器的结构和性能参数有关；y为采集到的符合计数向量，且为m维向量；x为n维的PET图像浓度分布向量，其每个元素值对应PET图像每个体素的体素值；e为测量噪声；x_t+1为第(t+1)帧PET图像浓度分布向量，x_t为第t帧PET图像浓度分布向量，A为状态矩阵。

根据PET图像的状态空间方程建立粒子权重估计算法的如下方程表达式：

$P(j,i)=\frac{D(i,j)\×\frac{x\left[j\right]}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}x\left[j\right]}}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\{D(i,j)\×\frac{x\left[j\right]}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}x\left[j\right]}\}}$

$w_{\left(k\right)}\left[j\right]=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\left\{y\right[i]\×P(j,i)\}-\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\{D(i,j)\×x_{\left(k\right)}[j\left]\right\}$

其中：x[j]为PET图像中第j体素的先验浓度值，y[i]为符合计数向量中的第i符合计数；D(i，j)为系统矩阵中第i行第j列的元素值；x_(k)[j]为PET图像中第j体素的第k粒子浓度值，w_(k)[j]为x_(k)[j]对应的粒子权重值；i、j和k均为自然数，且1≤i≤m，1≤j≤n，1≤k≤s，m为符合计数向量的维度，n为PET图像的体素个数。

最后，根据体素的先验浓度值和粒子浓度值，利用GPU多核并行处理技术通过粒子权重估计算法求得体素每一粒子浓度值对应的粒子权重值。

该步骤中涉及的很多矩阵向量计算，由于计算量大、并行程度高，可以通过并行计算处理技术来加速；将这些计算量映射到GPU多核中，以提高计算效率；GPU多核的具体计算处理策略如下：

对于向量与标量的乘法或除法，将该计算量映射到个数与向量维数相等的多个内核中，每个内核完成对应向量元素与标量相乘或相除；

对于矩阵与向量的乘法，将该计算量映射到个数与矩阵维数相等的多个内核中，每个内核完成对应矩阵元素与其相对应的向量元素相乘，然后对所得到的矩阵的各行再求和；

对于矩阵与标量的乘法或除法，将该计算量映射到个数与矩阵维数相等的多个内核中，每个内核完成对应矩阵元素与标量相乘或相除；

对于其余公式的代入计算，则使一个内核对应负责一个公式不同数据的计算。

(3)对粒子进行重采样。

利用GPU多核并行处理技术对粒子浓度值及其粒子权重值进行重采样，得到重采样后的粒子浓度值及其粒子权重值；

重采样算法如下(用伪代码表示)：

a.构建一个累计分布函数(CDF)：c₁＝0，

b.从一个归一化的随机分布中取一个数

c.for j＝1：Ns

$u_j=u_1+N_s^{-1}(j-1)$

While u_j＞c_i

*i＝i+1

end while

对样本重新赋值：

end for

重采样过程将权重大的粒子保留下来，权重小的粒子删除，从而使整个粒子的分布更接近真值；经过若干次迭代后，就能得到精度范围内的结果。

由于重采样过程中需要依次用到同一体素中的所有粒子，即各个粒子之间存在时序关系，因此同一体素的粒子重采样不能并行计算；故该步骤的并行计算处理策略是使不同体素的重采样互相并行，即一个内核对应执行一个体素的重采样算法。

(4)求取粒子浓度真值及其权重真值。

返回步骤(2)，并使重采样后的粒子浓度值作为步骤(2)中的粒子浓度值，依此循环直至重采样后的粒子浓度值收敛，将收敛时重采样后的粒子浓度值及其粒子权重值作为粒子浓度真值及其粒子权重真值。

(5)计算每一体素的体素值。

根据粒子浓度真值及其粒子权重真值，利用GPU多核并行处理技术通过以下方程表达式求取PET图像中每一体素的体素值。

$x^{*}\left[j\right]=\underset{k=1}{\overset{s}{\mathrm{\Σ}}}\left\{{x_{\left(k\right)}}^{*}\right[j]\×{w_{\left(k\right)}}^{*}[j\left]\right\}$

其中：x_(k) ^*[j]为PET图像中第j体素的第k粒子浓度真值，w_(k) ^*[j]为x_(k) ^*[j]对应的粒子权重真值，x^*[j]为PET图像中第j体素的体素值。

该步骤的并行计算处理策略也是使不同体素的体素值计算互相并行，即一个内核对应执行一个体素的体素值计算。

以下通过实验来验证本实施方式的实用性和可靠性；图2为著名的Zubal-thorax-phantom数字体模的模型，采样角度为180°，共有128个角度和128个半径对，最后所成图像像素为128×128；针对该模型，采用本实施方式重建的PET图像与采用ML-EM方法重建的PET图像相比较，从图3可以看出，用本实施方式重建出来的图像分辨率和图像清晰程度都比用ML-EM方法重建出来的要高得多。为了从统计学的角度分析本发明方法的优点，表1列出了在原始sinogram(符合计数向量)中分别加入10％、20％、30％、40％噪声的情况下，用本实施方式与用ML-EM方法重建得到的图像的偏差和标准差。

表1

从表1可以看出，本实施方式在偏差和标准差两个统计方面都明显优于现有技术常用的ML-EM方法，即证明了基于粒子滤波的静态PET图像重建方法相较于目前现有的重建方法的可靠性。

接下来将本实施方式与传统CPU串行计算实现的PET重建方法进行比较，如表2所示，比较两者之间在速率上的差异。

表2

sinogram大小	GPU并行计算	CPU串行计算	加速比
				32*32	0.1850s	23.4804s	126.9211
64*64	0.3730s	245.1592s	657.263
				128*128	10.917s	7610.6s	697.133

由表2可知，并行计算的参与大大缩短了重建时间；而并行计算使用的是单精度，从表3中可见，并行计算的单精度重建在精度(偏差和方差)上相对于串行计算并没有太大的损失。

表3

	偏差	方差
			GPU并行计算	0.0216	0.0931
CPU串行计算	0.0056	0.1687

Claims (5)

1.一种基于GPU多核并行处理的PET图像重建方法，包括如下步骤：

(1)根据体素值的概率分布区间，利用线性平均采样法对PET图像中的每一体素进行粒子采样，获得体素的每一粒子浓度值；

(2)根据采集获取到的符合计数向量，利用FBP法求得体素的先验浓度值；根据体素的先验浓度值和粒子浓度值，利用GPU多核并行处理技术通过粒子权重估计算法求得体素每一粒子浓度值对应的粒子权重值；

(3)利用GPU多核并行处理技术对所述的粒子浓度值及其粒子权重值进行重采样，得到重采样后的粒子浓度值及其粒子权重值；

(4)返回步骤(2)，并使重采样后的粒子浓度值作为步骤(2)中的粒子浓度值，依此循环直至重采样后的粒子浓度值收敛，将收敛时重采样后的粒子浓度值及其粒子权重值作为粒子浓度真值及其粒子权重真值；

(5)根据所述的粒子浓度真值及其粒子权重真值，利用GPU多核并行处理技术求取PET图像中每一体素的体素值。

2.根据权利要求1所述的基于GPU多核并行处理的PET图像重建方法，其特征在于：所述的线性平均采样法是通过对体素值的概率分布区间进行s-1等分，得到s个节点值，每个节点值即对应为任一体素的每一粒子浓度值，s为粒子采样个数。

3.根据权利要求1所述的基于GPU多核并行处理的PET图像重建方法，其特征在于：所述的粒子权重估计算法基于以下方程表达式：

$P(j,i)=\frac{D(i,j)\×\frac{x\left[j\right]}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}x\left[j\right]}}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\{D(i,j)\×\frac{x\left[j\right]}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}x\left[j\right]}\}}$

$w_{\left(k\right)}\left[j\right]=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\left\{y\right[i]\×P(j,i)\}-\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\{D(i,j)\×x_{\left(k\right)}[j\left]\right\}$

其中：x[j]为PET图像中第j体素的先验浓度值，y[i]为符合计数向量中的第i符合计数；D(i，j)为系统矩阵中第i行第j列的元素值；x_(k)[j]为PET图像中第j体素的第k粒子浓度值，w_(k)[j]为x_(k)[j]对应的粒子权重值；i、j和k均为自然数，且1≤i≤m，1≤j≤n，1≤k≤s，m为符合计数向量的维度，n为PET图像的体素个数。

4.根据权利要求1所述的基于GPU多核并行处理的PET图像重建方法，其特征在于：所述的步骤(5)中，通过以下方程表达式求取PET图像中每一体素的体素值：

$x^{*}\left[j\right]=\underset{k=1}{\overset{s}{\mathrm{\Σ}}}\left\{{x_{\left(k\right)}}^{*}\right[j]\×{w_{\left(k\right)}}^{*}[j\left]\right\}$

其中：x_(k) ^*[j]为PET图像中第j体素的第k粒子浓度真值，w_(k) ^*[j]为x_(k) ^*[j]对应的粒子权重真值，x^*[j]为PET图像中第j体素的体素值。

5.根据权利要求1所述的基于GPU多核并行处理的PET图像重建方法，其特征在于：所述的GPU多核并行处理技术的计算处理策略如下：

对于向量与标量的乘法或除法，将该计算量映射到个数与向量维数相等的多个内核中，每个内核完成对应向量元素与标量相乘或相除；

对于矩阵与向量的乘法，将该计算量映射到个数与矩阵维数相等的多个内核中，每个内核完成对应矩阵元素与其相对应的向量元素相乘，然后对所得到的矩阵的各行再求和；

对于矩阵与标量的乘法或除法，将该计算量映射到个数与矩阵维数相等的多个内核中，每个内核完成对应矩阵元素与标量相乘或相除；

对于其余公式的代入计算，则使一个内核对应负责一个公式不同数据的计算。

Images