CN100508891C - 在pet成像中最大后验优化图像重建方法 - Google Patents

Abstract

本发明公开了一种在PET成像中引入广义Gibbs先验的最大后验优化图像重建方法，包括如下步骤：(1)利用PET成像设备采集成像前的探测器数据，同时获取成像设备中各种数据校正参数值和系统矩阵；(2)根据步骤1获取的成像前校正数据满足的统计特征，构建用于重建图像的数学统计模型；(3)针对步骤2中数学模型的求解，引入广义Gibbs先验，采用最大后验估计方法进行重建模型转化，得到用于获取PET重建图像的带约束目标函数的优化方程；(4)由步骤3得到的结果，基于对优化方程中全局参数选择的基础上，采用抛物线替换坐标下降算法进行迭代计算处理，重建出图像。本发明可以大幅度提高PET重建图像质量。

Description

在PET成像中最大后验优化图像重建方法

技术领域

本发明涉及一种医学影像成像处理方法，尤其是涉及一种在PET成像中引入广义Gibbs先验的最大后验优化图像重建方法。

背景技术

正电子发射成像(PET)是当前医学影像学的重要检查手段之一，由于受到低计数率和噪声的影响，从探测数据到图像的重建在多数情况下是一个病态的问题。据文献报道，现有多数PET成像算法仍不能得到满意的图像。因此有效地提升PET成像质量在临床上存在巨大需要，也一直是医学PET成像的研究热点同时也是技术难题之一。

作为解决图像重建中的病态问题的有效方法，最大后验方法已被广泛接受。基于Gibbs随机场理论，Gibbs先验能量项可被引入到图像重建中抑制重建图像噪声。根据Gibbs先验能量方程形式的不同可将先验分粗分为两类：二次Gibbs先验和非二次Gibbs先验。常用的二次Gibbs先验形式上简单、理论上易于分析，且不会影响图像重建中整个能量方程的整体凹性；而大多具有边缘保持效果的最大后验方法，包括线位置模型和非连续自适应模型，则能够通过非二次的Gibbs先验能量函数实现保持边缘的图像重建。然而由于待重建图像并不总具备整体平滑特性，图像中存在的边缘和不可避免的附加噪声常带来某些区域一致性的突变。简单的二次Gibbs平滑先验，在消除噪声的同时也将使重建图像中边缘细节模糊，产生过平滑的负面效应。而对于非二次Gibbs先验，如果待重建图像由一些具有明显简单边缘和一些均匀区域组成，非二次Gibbs先验将产生相对较好的结果，但如果待重建图像性质趋于复杂，区域之间往往没有较明显的边界时(这正是PET成像的特征)，非二次Gibbs先验在PET成像时将产生易于被误诊为伪影的阶梯状均匀区域。加之非二次Gibbs先验将引进新的待定参数以及相关数值解法将大大增加整个重建的计算量和复杂度。

简单二次Gibbs先验通过在一个局部区域内像素值的平均效应来平滑含有噪声的重建图像；具有边缘保持作用的非二次Gibbs先验亦需要根据局部邻域内像素值间差量信息来确定目标图像中每个像素点的先验能量的程度和形式。为简化起见，称以上只取决于局部邻域信息的先验为传统Gibbs先验。

传统Gibbs先验所带来的负面的过平滑或阶梯状效应均是由于其所依赖的图像局部邻0域内的像素灰度值信息无法有效的区分边缘信息和噪声。为克服局部先验的上述缺点，曾出现了一些自称为使用图像中全局信息的最大后验方法，例如基于区域的Gibbs先验方法和通过水平集方法获取图像中全局信息的基于边界的Gibbs先验方法。然而实际上基于区域的Gibbs先验方法由于受到复杂的区域标识计算或所需的解剖学先验信息的限制，在实际应用中具有一定的局限性。而基于边界的Gibbs先验方法强依赖于水平集的参数化设计，此过程可能会产生不可预知的结果。

发明内容

本发明的目的在于提出一种PET成像中引入广义Gibbs先验的最大后验优化图像重建方法，可以大幅度提高PET重建图像质量。

本发明目的可通过以下技术措施来实现，包括步骤如下：

1.利用PET成像设备采集成像前的探测器数据，同时获取成像设备中各种数据校正参数值和系统矩阵；

2.根据步骤1获取的成像前校正数据满足的统计特征，构建用于重建图像的数学统计模型；

3.针对步骤2中数学模型的求解，引入广义Gibbs先验，采用最大后验估计方法进行重建模型转化，得到用于获取PET重建图像的带约束目标函数的优化方程；

4.由步骤3得到的结果，基于对优化方程中全局参数选择的基础上，采用抛物线替换坐标下降算法进行迭代计算处理，重建出优质的图像。

本发明步骤2中用到的数学统计模型为泊松分布或高斯分布，即正电子的探测过程实为一计数过程，将这些探测值理解为服从独立泊松分布或高斯分布的随机变量。

本发明步骤3中得到的图像重建模型实为一带约束的优化问题，其中广义Gibbs先验的具体设计过程为：

a、首先选择一个包含图像中丰富的几何信息的较大方形邻域；同时设计一个相似性测度用于比较大方形邻域内中像素点k和像素点j处对应的小方形邻域的相似性；

b、随后在选定的大方形邻域的内进行两像素间的灰度值比较的同时，利用两像素间相似

性来获得势能函数中的权值量。

本发明步骤3，b中的权值量定义为 $w_{\mathrm{kj}}^{\mathrm{GG}}=\exp (-{\left|\right|\mathrm{\λ}\left(V_k\right)-\mathrm{\λ}\left(V_j\right)\left|\right|}^2/h^2)/w_{\mathrm{kj}},$ w_kj定义为传统Gibbs先验中像素点k和像素点j之间的权值，通过图像域内两像点间的欧几里德距离的反比例函数确定；V_k和V_j则设定为以像素点k和像素点j为中心的小方形邻域；λ(V_k)和λ(V_j)为此两个邻域中所有像素灰度值数组；‖‖代表此两个像素点所在区域的加权欧几里德距离；参数h用于计算像素点间权值的指数函数同邻域相似性测度的反比例衰减关系。

上述参数h的确定过程为：首先直接对校正数据采用滤波反投影算法得到用于抛物线替换坐标下降迭代算法的初始图像；其后对该图像采用金字塔式结构由粗略到精细进行方差分析得到初始图像中较平滑区域的方差值σ；最后取h值为方差σ的倍数作为迭代求解时广义Gibbs先验公式

中的固定参数。

本发明步骤3，a中的相似性测度采用两像素点邻域内所有像素点灰度值的加权欧几里德距离的反比例函数。

本发明步骤4采用抛物线替换坐标下降迭代处理的具体过程如下：第一步，首先基于上一步迭代获得的重建图像为参照图像，获得待重建图像的逐像素点广义Gibbs先验的权值量

以作为下一步迭代之用；第二步，在第一步获取的权值量基础上，利用抛物线替换坐标下降迭代算法进行迭代重建；第三步，交替进行第一、二步直至收敛，获得最终重建图像。

本发明技术与现有技术的实验比较及结果如下：

首先应用本发明技术对模拟体模数据进行重建实验，图2是理想体模图像用于说明本发明的模拟实验对象。其中图2中(a)为体模图像1，表示一个Zubal腹部截面图；(b)为体模图像2，由一个均匀的背景、一个方形亮区域、一个内嵌方形暗区域的模糊圆形亮区域组成；(c)为体模图像3，表示一个标准的Shepp-Logan体模。本发明实验中设定此三个体模图像的重建具有相同的重建环境，sinogram数据中均加入了10％服从泊松分布的随机噪声。转换概率矩阵A，对应于一个平行带状积分几何模型，此几何模型表示一个180°的均匀区域里的具有128个径向取样和128个角采样的系统。由Fessler等人提供的ASPIRE软件系统生成。

图3中1)，2)，3)分别显示为采用滤波反投影(FBP)算法和二次Gibbs平滑先验，非二次Huber先验及本发明提出的广义Gibbs先验的最大后验方法对图2中三个体模的重建图像。由图3可看出，视觉上使用本发明提出的广义Gibbs先验的重建结果无论在抑制噪声还是保持边缘方面均明显优于其他三种重建结果。广义Gibbs先验的引入不仅能够克服二次Gibbs先验的过平滑效应，而且能够在很大程度上解决非二次先验所导致的阶梯状伪影的问题。

图4中1)，2)，3)描绘了理想体模图像和分别使用二次Gibbs先验，非二次Huber先验及本发明提出的广义Gibbs先验的最大后验重建图像的侧面水平轮廓图。由轮廓图所示，相对于FBP重建算法、二次Gibbs先验，非二次Huber先验的最大后验重建算法所重建的图像，使用本发明提出的广义Gibbs先验所重建的图像的侧面轮廓图无论在背景区域还是在边缘区域都更接近于真实图像的侧面轮廓图，从而说明使用本发明提出的广义Gibbs先验的重建算法能够更好的克服重建中的病态问题，从而能重建出更接近于真实体模的图像。下表给出了所有以上重建图像相对于其真实图像的信噪比。可以看出使用本发明提出的广义Gibbs先验的重建图像具有更高的信噪比。

为进一步说明本发明提出的广义Gibbs先验最大后验重建算法的有效性，对真实的PET心脏灌注探测数据进行重建实验。图5分别对真实的探测数据使用四种不同重建方法的重建图像，(a)为FBP重建图像；(b)为二次Gibbs先验重建图像；(c)为非二次Huber先验重建结果；(d)为本发明提出的广义Gibbs先验重建结果。可以看出，本发明提出的广义Gibbs先验重建图像在抑制噪声和保持边缘方面均明显优于其他三种重建图像。

附图说明

图1在使用本发明所述的广义Gibbs先验模型的情况下，目标图(每组图中左图)中的中心点的邻域权值分布图(每组图中右图)。(a)当中心点位于相对均匀的区域时，权值近似为一种图像卷积滤波的结果(如高斯滤波)，(b)当中心点位于一条竖直的边缘上时，权值基本上分布在此竖直边缘线上，(c)当中心点位于一条弯曲的边缘上，原图中位于此弯曲的边缘线上的点同样具有较大的权值，(d)当中心点位于一个有很多相似结构的背景中时，权值沿着图中相似的形态结构分布；

图2实验中的模拟体模数据。(a)体模图像1，(b)体模图像2，(c)体模图像3；

图3(1)～(3)分别对图2中三个模拟体模数据的使用四种不同重建方法的重建图像。(a)FBP重建，(b)二次Gibbs先验重建，(c)非二次Huber先验重建，(d)本发明提出的广义Gibbs先验重建；

图4(1)～(3)为对应图3(1)～(3)中最大后验重建图像的水平轮廓图；

图5分别对真实探测数据使用四种不同重建方法的重建图像。(a)FBP重建，(b)二次Gibbs先验重建，(c)非二次Huber先验重建(d)本发明提出的广义Gibbs先验重建。

图6本发明具体实施流程图。

具体实施方式

本发明实施有四个步骤(见图6)，具体如下：

1、利用PET成像系统采集成像前的探测器数据。具体的采集方式可以由用户灵活设定。实验中数据采集方式设计为：在一个180°的角度区间内，取128个径向采样和128个角度采样；系统矩阵A对应于平行束带状积分几何模型。将采样数据存入数组中。

2、数据校正。由系统获取的扫描时间校准系数、探测器的效率、衰减系数和死时间的校正系数c_i以及全部探测到的随机计数和散射计数r_i。根据参数c_i和r_i进行探测器数据校正，得到用于广义Gibbs先验最大后验图像重建的数据。

3、构建图像重建模型。采用最大后验方法，对步骤2得到的校正数据进行数学建模，完成广义Gibbs先验的设计，得到用于获取PET重建图像的带约束目标函数Φ(λ)的优化方程， $\widehat{\mathrm{\λ}}=\underset{\mathrm{\λ}\≥0}{\arg \max }\mathrm{\Φ}\left(\mathrm{\λ}\right),$ Φ(λ)＝L(y，λ)-βU_GG(λ)，其中L(y，λ)为校正数据y的对数似然能量方程； $U_{\mathrm{GG}}\left(\mathrm{\λ}\right)=\frac{1}{2}\underset{j}{\mathrm{\Σ}}\underset{k\∈N_j}{\mathrm{\Σ}}{w}_{\mathrm{kj}}^{\mathrm{GG}}{({\mathrm{\λ}}_k-{\mathrm{\λ}}_j)}^2$ 为广义Gibbs先验项，N_j表示像素点j的邻域系统，

是权值量；β为全局参数。

上述步骤3中广义Gibbs先验的具体设计过程：首先选择一个较大的方形邻域来包含图像中更丰富的几何信息，但亦不可过大，一般取11×11或15×15的大方形邻域。实验中大方形邻域的取值为11×11的方形邻域；同时设计一个相似性测度用于比较大方形邻域内像素点k和像素点j处小邻域的相似性，小邻域一般取5×5或7×7的方形邻域，相似性测度采用两个像素点邻域内所有像素点灰度值的加权欧几里德距离的反比例函数。实验中此相似性测度计算小邻域取值为7×7的方形邻域；随后在选定的小邻域内进行两像素间的灰度值比较的同时，利用两像素间相似性来获得势能函数中的权值量。权值量定义为 $w_{\mathrm{kj}}^{\mathrm{GG}}=\exp (-{\left|\right|\mathrm{\λ}\left(V_k\right)-\mathrm{\λ}\left(V_j\right)\left|\right|}^2/h^2)/w_{\mathrm{kj}},$ w_kj定义为传统Gibbs先验中像素点k和像素点j之间的权值，通过图像域内两像点间的欧几里德距离的反比函数确定，V_k和V_j则设定为以像素点k和像素点j为中心的小方形邻域，λ(V_k)和λ(V_j)为此两个邻域中所有像素灰度值数组，‖‖代表此两个像素点所在区域的加权欧几里德距离。参数h用于计算像素点间权值的指数函数同邻域相似性测度的反比例衰减关系。

为进一步说明权值量的选取，如图1所示，在每组图中，左边的图为原图，右边的图描绘了左图中心点邻域中各点在广义Gibbs先验模型中权值的取值情况，颜色越亮则相应的权值越大，此邻域设定为覆盖整个图像的区域。可以看出，权值一般分布于比较相似的结构处，两个像素点周围的结构越相似，则此两个像素点在先验中的权值就越大，表明广义Gibbs先验信息能够考虑到图像中的一些更丰富的几何结构形态的信息，从而克服传统Gibbs先验信息量相关较少的缺点，因此能够对PET图像重建的病态问题提供更为有效的正则化处理。

步骤3中广义Gibbs先验公式中h的选取过程为：首先直接对校正数据采用滤波反投影算法得到用于抛物线替换坐标下降迭代算法的初始图像；其后对该图像采用金字塔式结构由粗略到精细进行方差分析得到初始图像中较平滑区域的方差值σ；最后取h值为方差σ的10倍作为广义Gibbs先验公式

在抛物线替换坐标下降迭代算法迭代求解过程中的固定参数。

步骤3中全局参数β估计方法：由于PET成像中重建图像的质量很大程度上由探测器探测到的数据量counts决定。经大量实验验证，本发明中全局参数β可由经验公式 $\mathrm{\β}=\frac{1}{c}\log (\#\mathrm{counts})$ 确定，其中c为一常数，#counts表示用于重建的总counts个数。本发明的模拟实验(见图3)中c取值20，#counts取值3×10⁵；心脏灌注真实数据实验(见图5)中c取值10，#counts取值1.4876×10⁵。

4、完成重建。在步骤2和步骤3相关数据处理和模型设计的基础上，采用抛物线替换坐标下降迭代算法对目标函数进行二步式优化求解，具体如下：第一步，首先基于上一步迭代获得的重建图像为参照图像，进行待重建图像的逐像素点广义Gibbs先验的权值量

计算，以作为下一步迭代之用；第二步，在第一步获取的权值量基础上，第二步中该权值量作为常数，利用抛物线替换坐标下降迭代算法进行迭代重建。两步交替进行直至收敛，获得最终重建图像。本发明的实验中，迭代次数为150次时均可得到好的重建图。

Claims (4)

1、一种在PET成像中引入广义Gibbs先验的最大后验优化图像重建方法，其特征在于包括如下步骤：

(1)利用PET成像设备采集成像前的探测器数据，同时获取成像设备中各种数据校正参数值和系统矩阵；

(2)根据步骤1获取的成像前用于校正的数据满足的统计特征，构建用于重建图像的数学统计模型；

(3)针对步骤2中数学统计模型的构建，引入广义Gibbs先验，采用最大后验估计方法进行重建模型转化，得到用于获取PET重建图像的带约束目标函数的优化方程： $\widehat{\mathrm{\λ}}=\underset{\mathrm{\λ}\≥0}{\arg \max }\mathrm{\Φ}\left(\mathrm{\λ}\right),$ Φ(λ)＝L(y，λ)-βU_GG(λ)，其中L(y，λ)为校正数据y的对数似然能量方程； $U_{\mathrm{GG}}\left(\mathrm{\λ}\right)=\frac{1}{2}\underset{j}{\mathrm{\Σ}}\underset{k\∈N_j}{\mathrm{\Σ}}{w}_{\mathrm{kj}}^{\mathrm{GG}}{({\mathrm{\λ}}_k-{\mathrm{\λ}}_j)}^2$ 为广义Gibbs先验项，N_j表示像素点j的邻域系统，

是权值量；β为全局参数；

广义Gibbs先验的具体设计过程为：

a、首先选择一个包含图像中丰富的几何信息的较大方形邻域；同时设计一个相似性测度用于比较所述较大方形邻域内像素点k和像素点j处对应的小方形邻域的相似性；

b、随后在选定的小方形邻域内进行两像素间的灰度值比较的同时，利用两像素间相似性来获得势能函数中的权值量；

步骤a中的相似性测度采用两像素点邻域内所有像素点灰度值的加权欧几里德距离的反比例函数；

步骤b中的权值量定义为 $w_{\mathrm{kj}}^{\mathrm{GG}}=\exp (-{\left|\right|\mathrm{\λ}\left(V_k\right)-\mathrm{\λ}\left(V_j\right)\left|\right|}^2/h^2)/w_{\mathrm{kj}},$ w_kj定义为传统Gibbs先验中像素点k和像素点j之间的权值，通过图像域内两像点间的欧几里德距离的反比例函数确定；V_k和V_j则设定为以像素点k和像素点j为中心的小方形邻域；λ(V_k)和λ(V_j)为此两个邻域中所有像素灰度值数组；‖·‖代表此两个像素点所在区域的加权欧几里德距离；参数h用于计算像素点间权值的指数函数同邻域相似性测度的反比例衰减关系；

参数h的确定过程为：首先直接对成像前用于校正的数据采用滤波反投影算法得到用于抛物线替换坐标下降迭代算法的初始图像；其后对该图像采用金字塔式结构由粗略到精细进行方差分析得到初始图像中较平滑区域的方差值σ；最后取h值为方差σ的倍数作为迭代求解时广义Gibbs先验公式

中的固定参数；

(4)由步骤3得到的结果，基于对优化方程中全局参数选择的基础上，采用抛物线替换坐标下降算法进行迭代计算处理，重建后的图像，具体过程如下：第一步，首先基于上一步迭代获得的重建图像为参照图像，获得待重建图像的每个像素点处广义Gibbs先验的权值量，以作为下一步迭代之用；第二步，在第一步获取的权值量基础上，利用抛物线替换坐标下降迭代算法进行迭代重建；第三步，交替进行第一、二步直至收敛，获得最终重建图像。

2、根据权利要求1所述的在PET成像中引入广义Gibbs先验的最大后验优化图像重建方法，其特征在于：步骤2中用到的数学统计模型为泊松分布或高斯分布。

Images