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

摘要

本发明公开了一种基于广义熵与MR先验的PET图像最大后验重建方法，(1)利用PET成像设备采集PET成像前的探测数据，同时获取成像设备中各种数据校正参数值及成像设备的系统矩阵；(2)构建用于重建PET图像的数学统计模型；(3)针对数学统计模型求解，采用最大似然法得到PET初值图像；(4)将预先获取的MR图像与PET初值图像进行配准；(5)通过广义熵与配准后的MR图像引入解剖先验，采用最大后验方法对PET图像的数学统计模型进行重建模型转化，得到带约束目标函数的优化方程；(6)采用迟一步算法进行迭代计算，得到最终的医学影像的重建图像。本发明可提高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}{\Sigma }}{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}{\Sigma }}{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(maximum  likelihood-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}{\Sigma }}\phi \left(\frac{x-f_k}{{\sigma }_x}\right)\phi \left(\frac{y-a_k}{{\sigma }_y}\right).$

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

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

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

图1是本发明实验中的模拟体摸数据，模拟¹⁸F-FDG示踪剂在正常大脑的灰质(gray  matter，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重建，可以有效提高 分辨率，得到更为准确的重建图像。

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