CN102324089A - 基于广义熵与mr先验的pet图像最大后验重建方法 - Google Patents

Abstract

本发明公开了一种基于广义熵与MR先验的PET图像最大后验重建方法，(1)利用PET成像设备采集PET成像前的探测数据，同时获取成像设备中各种数据校正参数值及成像设备的系统矩阵；(2)构建用于重建PET图像的数学统计模型；(3)针对数学统计模型求解，采用最大似然法得到PET初值图像；(4)将预先获取的MR图像与PET初值图像进行配准；(5)通过广义熵与配准后的MR图像引入解剖先验，采用最大后验方法对PET图像的数学统计模型进行重建模型转化，得到带约束目标函数的优化方程；(6)采用迟一步算法进行迭代计算，得到最终的医学影像的重建图像。本发明可提高PET重建图像的视觉效果和量化指标。

Description

基于广义熵与MR先验的PET图像最大后验重建方法

技术领域

本发明涉及一种医学影像的图像重建方法，具体是指一种基于广义熵与MR先验的PET图像最大后验重建方法。

背景技术

正电子发射成像(PET)是一种非介入的定量研究活体功能活度的临床工具。由于低空间分辨率和系统固有噪声，PET重建是一个病态问题。统计图像重建方法，比如最大似然-期望最大法(maximum likelihood-expectation maximization，ML-EM)，能够更好的考虑系统模型的物理效应而且能够针对探测数据和噪声的统计特性建立数学模型，其迭代重建的图像质量要优于传统的以滤波反投影(Filtered Back Projection，FBP)方法为代表的解析重建算法。然而，ML-EM方法在迭代过程中会伴随着图像质量退化而导致的棋盘效应，从而导致非收敛的迭代过程。此病态问题可以通过贝叶斯方法有效的求解。基于贝叶斯理论，先验知识可以对原始的重建进行正则化，所以先验的选择对于最大后验方法极为关键。近二十年来，绝大部分先验以马尔可夫(MRF)先验的形式出现，先验通常反映图像局部邻域的平滑特性。最常用的为二次先验(quadratic prior，QP)，在抑制噪声的同时，使边缘细节模糊。

与简单的使用图像自身先验或马尔可夫形式先验相比，在PET图像重建中引入解剖图像信息已经引起广泛的关注。MR(Magnetic Resonance)图像为磁共振图像，由于来自高分辨率的MR/CT解剖图像能为PET图像重建提供大量的先验信息，所以已有很多在PET重建中引入解剖先验的工作。一般来说，解剖先验的使用可以分为以下两类：第一类方法基于解剖图像的边缘信息，此种方法基于解剖图像与功能图像的边缘相互对应，通过分割或者标记解剖图像，调节先验权值来惩罚功能图像边缘区域的像素灰度变化；第二类方法基于解剖图像的区域信息，此种方法假设在每一个解剖区域或器官中，放射性活度分布是均匀的。显然，这两类方法都需要获取准确的解剖边缘或者对解剖图像进行准确的分割。

近来，在PET图像重建中，信息理论测度(包括互信息与联合熵)引导的解剖信息作为先验已经被深入研究。此种方法不需要对解剖图像进行准确的分割或边缘提取，可以重建得到与解剖图像具有相似灰度分布的PET图像。特别是，Somayajula使用MR图像与PET图像的特征向量的互信息定义先验。之后，Nuyts通过实验说明，在解剖图像与PET图像的整体分布具有差异时，互信息作为先验在图像重建中会引入偏差，而联合熵是一种更具有鲁棒性的测度。总的来说，所有的这些信息理论测度都是基于香农熵(Shannon entropy)或玻尔兹曼吉伯斯(Bolzmann-Gibbs)统计。

发明内容

本发明的目的在于提供一种基于广义熵与MR先验的PET图像最大后验重建方法，该方法可以有效利用MR图像信息，大幅提高PET重建图像的视觉效果和量化指标。

本发明的目的可通过以下的技术措施来实现：一种基于广义熵与MR先验的PET图像最大后验重建方法，其包括如下步骤：

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

(2)根据步骤1获取的探测数据满足的统计特征，即正电子的探测过程是计数过程，将这些探测数据作为服从独立泊松分布的随机变量，构建用于重建PET图像的数学统计模型；

(3)针对步骤2构建的数学统计模型求解，采用最大似然法得到PET初值图像；

(4)将预先获取的MR图像与步骤3得到的PET初值图像进行配准，获得配准后的MR图像；

(5)针对步骤2中的数学统计模型，通过广义熵与步骤4获得的配准后的MR图像引入解剖先验，采用最大后验方法对步骤2构建的PET图像的数学统计模型进行重建模型转化，得到用于获取PET重建图像的带约束目标函数的优化方程；

(6)由步骤5得到的结果，基于对PET重建图像的优化方程中全局参数选择的基础上，采用迟一步算法进行迭代计算，得到最终的医学影像的重建图像。

本发明步骤2中用到的数学统计模型为泊松分布：

$g_i~\mathrm{Poisson}\{c_i\underset{j}{\overset{M}{\mathrm{\Σ}}}{a}_{\mathrm{ij}}{f}_j+r_i\},i=\mathrm{1,2},...,N$

其中，g_i表示第i对探测器探测到的来自发射扫描的光子数据，N表示探测器对的数目，f_j表示第j个重建像素点处的同位素分布，M表示待重建图像像素的总个数，r_i表示在发射扫描中第i对探测器总共探测到的随机计数和散射计数，A＝{a_ij}是成像设备的系统矩阵，a_ij定义为在理想条件下图像像素j被探测器对i探测到的几何概率。c_i表示扫描时间的校准系数，探测器的效率，衰减系数和死时间的校正系数的综合值。

本发明步骤3中采用经典的最大似然-期望最大法得到PET初值图像

最大似然估计：

其中L(g|f)为校正数据g的对数似然能量方程。

本发明步骤4中采用简单的刚性配准即可得到配准后的MR图像。

本发明步骤5中通过广义熵与配准后的MR图像引入解剖先验，采用最大后验方法对步骤2构建的PET图像的数学统计模型进行重建模型转化，得到的图像重建模型为带约束目标函数的优化方程：

L(f)＝L(g|f)-βD(x，y)，其中L(g|f)为校正数据g的对数似然能量方程，β为全局参数，D(X，Y)为基于广义熵与MR图像的先验项。

本发明步骤5中解剖先验的先验方程为：p(f)＝Z^-1×exp(-β×D(X，Y))，其中Z为正火常数，β为全局超参数，X为PET图像的灰度特征，Y为MR图像的灰度特征，D(X，Y)为基于广义熵与MR图像的先验项。

本发明中，所述基于广义熵与MR图像的先验项D(X，Y)的构建过程为：分别选取配准后的MR图像与PET图像的灰度特征作为两个随机变量，基于广义熵定义该两个随机变量的互信息与联合熵作为Gibbs先验的势函数，D(X，Y)中两个随机变量的联合分布与其边缘分布的概率密度函数，采用Parzen窗计算。

上述步骤中，选择高斯窗函数作为Parzen窗计算概率密度函数，使得势函数可微可导。

本发明步骤6中迟一步算法迭代的过程为：第一步，选取步骤3得到的PET初值图像；第二步，采用迟一步算法对步骤5获得的PET重建图像的优化方程进行迭代计算，得到最终的医学影像的重建图像。

与现有技术相比，本发明的基于广义熵与MR先验的PET图像最大后验重建方法不但可以有效利用MR图像信息，与使用二次先验的PET最大后验重建图像相比，还可以大幅提高PET重建图像的视觉效果以及诸如PET重建图像剖面边缘保持效果，偏差与噪声曲线等方面的量化指标。

附图说明

图1为本发明实验中的模拟体摸数据：(a)PET体模图像，(b)MR图像；

图2为本发明体模图像与重建图像：(a)PET体模图像，(b)MI重建图像，(c)-(e)Tsallis-MI重建图像，α分别为0.5，0.9，1.5，(f)QP重建图像，(g)JE重建图像，(h)-(j)Tsallis-JE重建图像，α分别为0.5，0.9，1.5；

图3为本发明MR图像分别与图2(a)-(j)图像的联合分布：MR图像与(a)PET体模图像，(b)MI重建图像，(c)-(e)Tsallis-MI重建图像，α分别为0.5，0.9，1.5，(f)QP重建图像，(g)JE重建图像，(h)-(j)Tsallis-JE重建图像，α分别为0.5，0.9，1.5；

图4A分别对图2中图像预设的三条水平线轮廓图，(a)-(b)水平线L1；

图4B分别对图2中图像预设的三条水平线轮廓图，(c)-(d)水平线L2；

图4C分别对图2中图像预设的三条水平线轮廓图，(e)-(f)水平线L3；

图5为本发明感兴趣区域归一化的均方误差(normalized mean squared error，NMSE)与归一化的标准差(normalized standard deviation，NSD)曲线，(a)-(b)白质区域，(c)-(d)灰质区域；

图6是本发明的框架示意图。

具体实施方式

如图6所示，本发明一种基于广义熵与MR先验的PET图像最大后验重建方法的优化实施例，包括如下步骤：

(1)利用PET成像设备采集PET成像前的探测数据，同时获取成像设备中各种数据校正参数值及成像设备的系统矩阵，具体的采集方式可以根据实际需要设定。本发明的实验中数据采集方式设计为：在一个180°的角度区间内，取128个径向采样和128个角度采样；系统矩阵A对应于平行束带状积分几何模型。将采样数据存入数组中。数据校正由系统获取的扫描时间校准系数、探测器的效率、衰减系数和死时间的校正系数c_i以及全部探测到的随机计数和散射计数r_i。根据c_i和r_i进行探测器数据校正，得到用于重建的数据；

(2)根据步骤1获取的探测数据满足的统计特征，构建用于重建PET图像的数学统计模型，该数学统计模型为泊松分布，即正电子的探测过程是计数过程，将这些探测值理解为服从独立泊松分布的随机变量：

$g_i~\mathrm{Poisson}\{c_i\underset{j}{\overset{M}{\mathrm{\Σ}}}{a}_{\mathrm{ij}}{f}_j+r_i\},i=\mathrm{1,2},...,N$

其中，g_i表示第i对探测器探测到的来自发射扫描的光子数据，N表示探测器对的数目，f_j表示第j个重建像素点处的同位素分布，M表示待重建图像像素的总个数，r_i表示在发射扫描中第i对探测器总共探测到的随机计数和散射计数，A＝{a_ij}是成像设备的系统矩阵，a_ij定义为在理想条件下图像像素j被探测器对i探测到的几何概率。c_i表示扫描时间的校准系数，探测器的效率，衰减系数和死时间的校正系数的综合值；

(3)针对步骤2构建的数学统计模型求解，采用经典10步ML-EM(maximumlikelihood-expectation maximization，ML-EM，最大似然-期望最大法)得到PET初值图像最大似然估计：其中L(g|f)为校正数据g的对数似然能量方程。；

(4)将预先获取的MR图像与步骤3得到的PET初值图像采用简单的刚性配准，获得配准后的MR图像；

(5)针对步骤2中的数学统计模型，通过广义熵与步骤4获得的配准后的MR图像引入解剖先验，采用最大后验方法对步骤2构建的PET图像的数学统计模型进行重建模型转化，得到用于获取PET重建图像的带约束目标函数的优化方程；

(6)由步骤5得到的结果，基于对PET重建图像的优化方程中全局参数选择的基础上，采用迟一步算法进行迭代计算，采用迟一步算法进行迭代计算分两步进行，第一步，选取步骤3得到的PET初值图像；第二步，采用迟一步算法对步骤5获得的PET重建图像的优化方程进行迭代计算，得到最终的医学影像的重建图像。

步骤5中通过广义熵与配准后的MR图像引入解剖先验，采用最大后验方法对步骤2构建的PET图像的数学统计模型进行重建模型转化，得到的图像重建模型为带约束目标函数的优化方程：

L(f)＝L(g|f)-βD(x，y)，其中L(g|f)为校正数据g的对数似然能量方程，β为全局参数，D(X，Y)为基于广义熵与MR图像的先验项。

所述解剖先验的先验方程为：p(f)＝Z^-1×exp(-β×D(X，Y))，其中Z为正火常数，β为全局超参数，X为PET图像的灰度特征，Y为MR图像的灰度特征，D(X，Y)为基于广义熵与MR图像的先验项。

所述基于广义熵与MR图像的先验项D(X，Y)的构建过程为：分别选取配准后的MR图像与PET图像的灰度特征作为两个随机变量，基于广义熵定义该两个随机变量的互信息与联合熵作为Gibbs先验的势函数，D(X，Y)中两个随机变量的联合分布与其边缘分布的概率密度函数，采用Parzen窗计算。

选择高斯窗函数作为Parzen窗计算概率密度函数，使得势函数可微可导，通过Parzen窗计算两个随机变量的联合分布与其边缘分布的概率密度函数。

步骤5中的先验项D(X，Y)的具体设计过程如下：首先由PET图像和MR图像提取N_S个特征向量(包括灰度，梯度等)，分别表示为x_i和y_i，i＝1，2，...，Ns。本发明中，我们只考虑两幅图像的灰度特征。这些特征向量可以看作独立的随机变量X和Y。基于广义熵(本发明中我们考虑Tsallis熵，一种广义香农熵)的互信息可定义为I^α(X，Y)＝H^α(X)+H^α(Y)-(1-α)×H^α(X)H^α(Y)-H^α(X，Y)，H^α(X)＝(1-α)^-1(∑_xp(x)^α-1)，α为可调节参数。基于广义熵的联合熵可定义为H^α(X，Y)＝(1-α)^-1(∑_x，yp(x，y)^α-1)。我们可以通过最大化互信息或最小化联合熵，得到与解剖图像灰度分布相似的图像，则先验项D(X，Y)可以定义如下D(X，Y)＝-I^α(X，Y)或D(X，Y)＝H^α(X，Y)。我们采用非参数方法来估计概率密度函数，设采样点的数目等于图像中像素点的个数N。PET图像中所有像素表示为f₁，f₂，...，

MR图像中所有像素表示为a₁，a₂，...，a_N，那么PET图像的灰度的概率密度p(x)的Parzen窗估计

可以表示为

其中

是宽度为σ的窗函数。本发明中，我们选择均值为0标准差为σ的高斯窗函数作为窗函数。在Parzen窗估计时，将图像灰度范围range划分为M级(通常认为可以区分不同的区域即可，本发明中取M＝128)，则每一个级表示的值Δ＝range/M。σ大小设为5Δ即可。为估计PET图像与MR图像的联合概率密度，我们采用协方差矩阵为

的高斯窗函数。所以，PET图像与MR图像的灰度的联合概率密度可以表示为

$\hat{p}(x,y)=\frac{1}{N_s}\underset{k=1}{\overset{N_s}{\mathrm{\Σ}}}\mathrm{\φ}\left(\frac{x-f_k}{{\mathrm{\σ}}_x}\right)\mathrm{\φ}\left(\frac{y-a_k}{{\mathrm{\σ}}_y}\right).$

上述先验项D(X，Y)中的参数α为可调节参数，

α对1取极限时，上述基于广义熵的互信息与联合熵即为传统的基于香农熵的互信息与联合熵。经我们研究发现，α取值小于1，如α为0.5和0.9，基于广义熵的互信息可取得较传统的基于香农熵的互信息更优的重建结果；而基于广义熵的联合熵可取得与传统的基于香农熵的联合熵相当的重建结果。

本发明步骤2中用到的数学统计模型也可以采用高斯分布，即正电子的探测过程认为是一计数过程，将这些探测数据也服从独立高斯分布的随机变量。

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

图1是本发明实验中的模拟体摸数据，模拟¹⁸F-FDG示踪剂在正常大脑的灰质(graymatter，GM)，白质(white matter，WM)和脑脊液(cerebrospinal fluid，CSF)中产生近似均匀的分布。图1(a)为PET体模图像。采用大小为128×128的Hoffman脑体模切片，在灰质，白质和脑脊液区域的活度比为4∶1∶0。图1(b)为MR图像。采用同一切片，同样大小的Hoffman脑体模切片，在在灰质，白质和脑脊液区域的像素值分别为180，255，0。本发明实验中设定重建环境，sinogram数据中均加入了10％服从泊松分布的随机噪声。转换概率矩阵A，对应于一个平行带状积分几何模型，此几何模型表示一个180°的均匀区域里的具有128个径向取样和128个角采样的系统。由Fessler等人提供的ASPIRE软件系统生成。

图2是本发明PET体模图像与重建图像。(a)为PET体模图像，(b)为MI重建图像，(c)-(e)为Tsallis-MI重建图像，α分别为0.5，0.9，1.5。(f)为QP重建图像，(g)为JE重建图像，(h)-(j)为Tsallis-JE重建图像，α分别为0.5，0.9，1.5。整体来看，使用MR先验得到的重建图像要优于QP重建。由(b)-(e)比较，可以看到，基于Tsallis-MI的最大后验重建在α为0.5，0.9时，较基于MI的最大后验重建在抑制噪声保持边缘上具有更优的效果。由(g)-(j)比较，可以看到基于JE的最大后验重建与基于Tsallis-JE的最大后验重建均有较优的重建结果。

图3是本发明MR图像分别与图2(a)-(j)图像的联合概率分布。MR图像与(a)PET体模图像，(b)MI重建图像，(c)-(e)Tsallis-MI重建图像，α分别为0.5，0.9，1.5。(f)QP重建图像，(g)JE重建图像，(h)-(j)Tsallis-JE重建图像，α分别为0.5，0.9，1.5。可以看到基于Tsallis-MI与Tsallis-JE的最大后验重建在α等于0.5时，最接近于MR与PET体模图像的联合概率分布。以上说明，基于Tsallis熵与MR先验的最大后验重建，在选择合适的参数α下，可以得到较其他方法更优的图像。

图4A至图4C中描绘了PET体模图像和分别使用QP，MI，JE，和本发明提出的Tsallis-MI与Tsallis-JE的重建图像在三条感兴趣水平线的轮廓图。由图所示，基于Tsallis-MI的最大后验重建，在α小于1，即α为0.5，0.9时，要明显优于使用QP与MI的最大后验重建。基于JE的最大后验重建与Tsallis-JE的最大后验重建均有较QP重建更优的结果。

图5中描绘了随参数β变化感兴趣区域归一化的均方误差(normalized mean squared error，NMSE)与归一化的标准差(normalized standard deviation，NSD)曲线。由图所示，QP重建的NMSE与NSD曲线做为参考，可以看到使用Tsallis-MI重建得到的曲线特性要优于QP重建重建，而且当α小于1时，其曲线特性要明显优于MI重建的曲线特性。同时可以看到，JE重建与Tsallis-JE重建的曲线特性不同传统的QP重建，随着参数β增大，其NMSE与NSD同时下降。以上分析表明Tsallis-MI重建与Tsallis-JE重建较传统的QP重建，可以有效提高分辨率，得到更为准确的重建图像。

本发明的实施方式不限于此，在本发明上述基本技术思想前提下，按照本领域的普通技术知识和惯用手段对本发明内容所做出其它多种形式的修改、替换或变更，均落在本发明权利保护范围之内。

Claims (9)

1.一种基于广义熵与MR先验的PET图像最大后验重建方法，其包括如下步骤：

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

(2)根据步骤1获取的探测数据满足的统计特征，即正电子的探测过程是计数过程，将这些探测数据作为服从独立泊松分布的随机变量，构建用于重建PET图像的数学统计模型；

(3)针对步骤2构建的数学统计模型求解，采用最大似然法得到PET初值图像；

(4)将预先获取的MR图像与步骤3得到的PET初值图像进行配准，获得配准后的MR图像；

(5)针对步骤2中的数学统计模型，通过广义熵与步骤4获得的配准后的MR图像引入解剖先验，采用最大后验方法对步骤2构建的PET图像的数学统计模型进行重建模型转化，得到用于获取PET重建图像的带约束目标函数的优化方程；

(6)由步骤5得到的结果，基于对PET重建图像的优化方程中全局参数选择的基础上，采用迟一步算法进行迭代计算，得到最终的医学影像的重建图像。

2.根据权利要求1所述的基于广义熵与MR先验的PET图像最大后验重建方法，其特征在于：所述步骤2中用到的数学统计模型为泊松分布：

$g_i~\mathrm{Poisson}\{c_i\underset{j}{\overset{M}{\mathrm{\Σ}}}{a}_{\mathrm{ij}}{f}_j+r_i\},i=\mathrm{1,2},...,N$

其中，g_i表示第i对探测器探测到的来自发射扫描的光子数据，N表示探测器对的数目，f_j表示第j个重建像素点处的同位素分布，M表示待重建图像像素的总个数，r_i表示在发射扫描中第i对探测器总共探测到的随机计数和散射计数，A＝{a_ij}是成像设备的系统矩阵，a_ij定义为在理想条件下图像像素j被探测器对i探测到的几何概率。c_i表示扫描时间的校准系数，探测器的效率，衰减系数和死时间的校正系数的综合值。

3.根据权利要求1所述的基于广义熵与MR先验的PET图像最大后验重建方法，其特征在于：步骤3中采用经典的最大似然-期望最大法得到PET初值图像最大似然估计：

其中L(g|f)为校正数据g的对数似然能量方程。

4.根据权利要求1所述的基于广义熵与MR先验的PET图像最大后验重建方法，其特征在于：步骤4中采用刚性配准得到配准后的MR图像。

5.根据权利要求1所述的基于广义熵与MR先验的PET图像最大后验重建方法，其特征在于：步骤5中通过广义熵与配准后的MR图像引入解剖先验，采用最大后验方法对步骤2构建的PET图像的数学统计模型进行重建模型转化，得到的图像重建模型为带约束目标函数的优化方程：

L(f)＝L(g|f)-βD(x，y)，其中L(g|f)为校正数据g的对数似然能量方程，β为全局参数，D(X，Y)为基于广义熵与MR图像的先验项。

6.根据权利要求5所述的基于广义熵与MR先验的PET图像最大后验重建方法，其特征在于：所述解剖先验的先验方程为：p(f)＝Z^-1×exp(-β×D(X，Y))，其中Z为正火常数，β为全局超参数，X为PET图像的灰度特征，Y为MR图像的灰度特征，D(X，Y)为基于广义熵与MR图像的先验项。

7.根据权利要求5或6所述的基于广义熵与MR先验的PET图像最大后验重建方法，其特征在于：所述基于广义熵与MR图像的先验项D(X，Y)的构建过程为：分别选取配准后的MR图像与PET图像的灰度特征作为两个随机变量，基于广义熵定义该两个随机变量的互信息与联合熵作为Gibbs先验的势函数，D(X，Y)中两个随机变量的联合分布与其边缘分布的概率密度函数，采用Parzen窗计算。

8.根据权利要求7所述的基于广义熵与MR先验的PET图像最大后验重建方法，其特征在于：选择高斯窗函数作为Parzen窗计算概率密度函数，使得势函数可微可导。

9.根据权利要求1所述的基于广义熵与MR先验的PET图像最大后验重建方法，其特征在于：步骤6中迟一步算法迭代的过程为：第一步，选取步骤3得到的PET初值图像；第二步，采用迟一步算法对步骤5获得的PET重建图像的优化方程进行迭代计算，得到最终的医学影像的重建图像。

Images