专利名称：在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先验所带来的负面的过平滑或阶梯状效应均是由于其所依赖的图像局部邻O域内的像素灰度值信息无法有效的区分边缘信息和噪声。为克服局部先验的上述缺点，曾出现 了一些自称为使用图像中全局信息的最大后验方法，例如基于区域的Gibbs先验方法和通过水 平集方法获取图像中全局信息的基于边界的Gibbs先验方法。然而实际上基于区域的Gibbs先 验方法由于受到复杂的区域标识计算或所需的解剖学先验信息的限制，在实际应用中具有一定的局限性。而基于边界的Gibbs先验方法强依赖于水平集的参数化设计，此过程可能会产生 不可预知的结果。发明内容本发明的目的在于提出一种PET成像中引入广义Gibbs先验的最大后验优化图像重建方 法，可以大幅度提高PET重建图像质量。本发明目的可通过以下技术措施来实现，包括步骤如下1. 利用PET成像设备采集成像前的探测器数据，同时获取成像设备中各种数据校正参 数值和系统矩阵2. 根据步骤1获取的成像前校正数据满足的统计特征，构建用于重建图像的数学统计模型；3. 针对步骤2中数学模型的求解，引入广义Gibbs先验，采用最大后验估计方法进行重 建模型转化，得到用于获取PET重建图像的带约^目标函数的优化方程；4. 由步骤3得到的结果，基于对优化方程中全局参数选择的基础上，采用抛物线替换坐 标下降算法进行迭代计算处理，重建出优质的图像。本发明步骤2中用到的数学统计模型为泊松分布或高斯分布，即正电子的探测过程实为 一计数过程，将这些探测值理解为服从独立泊松分布或高斯分布的随机变量。本发明步骤3中得到的图像重建模型实为一带约束的优化问题，其中广义Gibbs先验的 具体设计过程为a、 首先选择一个包含图像中丰富的几何信息的较大方形邻域；同时设计一个相似性测度 用于比较大方形邻域内中像素点和像素点y'处对应的小方形邻域的相似性；b、 随后在选定的大方形邻域的内进行两像素间的灰度值比较的同时，利用两像素间相似 性来获得势能函数中的权值量。本发明步骤3， b中的权值量定义为wf A(H2//^/^ ， 定义为传统Gibbs先验中像素点A:和像素点y'之间的权值，通过图像域内两像点间的欧几里德距离的反比例函 数确定J^和^则设定为以像素点t和像素点乂为中心的小方形邻域；义(W和义(。为此两个邻域中所有像素灰度值数组；Hl代表此两个像素点所在区域的加权欧几里德距离；参数&用于计算像素点间权值的指数函数同邻域相似性测度的反比例衰减关系。上述参数/7的确定过程为首先直接对校正数据采用滤波反投影算法得到用于抛物线替
换坐标下降迭代算法的初始图像；其后对该图像采用金字塔式结构由粗略到精细进行方差分 析得到初始图像中较平滑区域的方差值^ ;最后取A值为方差cr的倍数作为迭代求解时广义 Gibbs先验公式、f中的固定参数。本发明步骤3, a中的相似性测度采用两像素点邻域内所有像素点灰度值的加权欧几里德 距离的反比例函数。本发明步骤4采用抛物线替换坐标下降迭代处理的具体过程如下第一步，首先基于上 一步迭代获得的重建图像为参照图像，获得待重建图像的逐像素点广义Gibbs先验的权值量 wf ，以作为下一步迭代之用；第二步，在第一步获取的权值量基础上，利用抛物线替换坐标下降迭代算法进行迭代重建；第三步，交替进行第一、二步直至收敛，获得最终重建图像。 本发明技术与现有技术的实验比较及结果如下首先应用本发明技术对模拟体模数据进行重建实验，图2是理想体模图像用于说明本发 明的模拟实验对象。其中图2中(a)为体模图像l，表示一个Zubal腹部截面图；(b)为体模图 像2，由一个均匀的背景、 一个方形亮区域、 一个内嵌方形暗区域的模糊圆形亮区域组成；(c) 为体模图像3，表示一个标准的Shepp-Logan体模。本发明实验中设定此三个体模图像的重建 具有相同的重建环境，sinogram数据中均加入了 10%服从泊松分布的随机噪声。转换概率矩阵^，对应于一个平行带状积分几何模型，此几何模型表示一个180。的均匀区域里的具有128个径向取样和128个角采样的系统。由Fessler等人提供的ASPIRE软件系统生成。图3中1), 2), 3)分别显示为采用滤波反投影(FBP)算法和二次Gibbs平滑先验，非二次 Hnber先验及本发明提出的广义Gibbs先验的最大后验方法对图2中三个体模的重建图像。由 图3可看出，视觉上使用本发明提出的广义Gibbs先验的重建结果无论在抑制噪声还是保持 边缘方面均明显优于其他三种重建结果。广义Gibbs先验的引入不仅能够克服二次Gibbs先 验的过平滑效应，而且能够在很大程度上解决非二次先验所导致的阶梯状伪影的问题。.图4中1)， 2), 3)描绘了理想体模图像和分别使用二次Gibbs先验，非二次Huber先验及本 发明提出的广义Gibbs先验的最大后验重建图像的侧面水平轮廓图。由轮廓图所示，相对于 FBP重建算法、二次Gibbs先验，非二次Huber先验的最大后验重建算法所重建的图像，使 用本发明提出的广义Gibbs先验所重建的图像的侧面轮廓图无论在背景区域还是在边缘区域 都更接近于真实图像的侧面轮廓图，从而说明使用本发明提出的广义Gibbs先验的重建算法 能够更好的克服重建中的病态问题，从而能重建出更接近于真实体模的图像。下表给出了所 有以上重建图像相对于其真实图像的信噪比。可以看出使用本发明提出的广义Gibbs先验的
重建图像具有更高的信噪比。FBP重 建二次Gibbs先验 重建非二次Huber先 验重建广义Gibbs先验 重建体模图像l之信 噪比11.3413.1613.5014.78体模图像2之信 噪比11.6713.0713.4914.54体模图像3之信 噪比10.6712.2113.1514.18为进一步说明本发明提出的广义Gibbs先验最大后验重建算法的有效性，对真实的PET 心脏灌注探测数据进行重建实验。图5分别对真实的探测数据使用四种不同重建方法的重建 图像，(a)为FBP重建图像；(b)为二次Gibbs先验重建图像；(c)为非二次Huber先验重建结 果；(d)为本发明提出的广义Gibbs先验重建结果。可以看出，本发明提出的广义Gibbs先验 重建图像在抑制噪声和保持边缘方面均明显优于其他三种重建图像。
图l在使用本发明所述的广义Gibbs先验模型的情况下，目标图(每组图中左图)中的中 心点的邻域权值分布图(每组图中右图)。(a)当中心点位于相对均匀的区域时，权值近似为一 种图像巻积滤波的结果(如高斯滤波)，(b)当中心点位于一条竖直的边缘上时，权值基本上分 布在此竖直边缘线上，(c)当中心点位于一条弯曲的边缘上，原图中位于此弯曲的边缘线上的 点同样具有较大的权值，(d)当中心点位于一个有很多相似结构的背景中时，权值沿着图中相 似的形态结构分布；图2实验中的模拟体模数据。(a)体模图像l， (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,以及全部探测到的随机计数和散射计数/)。根据参数。和r,.进行探测器数据校正， 得到用于广义Gibbs先验最大后验图像重建的数据。3、 构建图像重建模型。采用最大后验方法'，对步骤2得到的校正数据进行数学建模， 完成广义Gibbs先验的设计，得到用于获取PET重建图像的带约束目标函数4)("的优化方程， i = argmax<5W， = -/3f^(义)，其中i(;v，;i)为校正数据y的对数似然能量方程；"goW—SSwf^-义」2为广义Gibbs先验项， 表示像素点y的邻域系统'《。是权值量；p为全局参数。上述步骤3中广义Gibbs先验的具体设计过程首先选择一个较大的方形邻域来包含图 像中更丰富的几何信息，但亦不可过大， 一般取11x11或15x15的大方形邻域。实验中大方 形邻域的取值为11x11的方形邻域；同时设计一个相似性测度用于比较大方形邻域内像素点A和像素点j处小邻域的相似性，小邻域一般取5x5或7x7的方形邻域，相似性测度采用两个像素点邻域内所有像素点灰度值的加权欧几里德距离的反比例函数。实验中此相似性测度计 算小邻域取值为7x7的方形邻域；随后在选定的小邻域内进行两像素间的灰度值比较的同时， 利用两像素间相似性来获得势能函数中的权值量。权值量定义为wf —l卜(W-A^)112/^/ ， 定义为传统Gibbs先验中像素点t和像素点j之间的权值， 通过图像域内两像点间的欧几里德距离的反比函数确定，^和^则设定为以像素点^和像素 点7'为中心的小方形邻域，A(F"和义(。为此两个邻域中所有像素灰度值数组，lll代表此两个像素点所在区域的加权欧几里德距离。参数/7用于计算像素点间权值的指数函数同邻域相似性测度的反比例衰减关系。
为进一步说明权值量的选取，如图1所示，在每组图中，左边的图为原图，右边的图描 绘了左图中心点邻域中各点在广义Gibbs先验模型中权值的取值情况，颜色越亮则相应的权 值越大，此邻域设定为覆盖整个图像的区域。可以看出，权值一般分布于比较相似的结构处， 两个像素点周围的结构越相似，则此两个像素点在先验中的权值就越大，表明广义Gibbs先 验信息能够考虑到图像中的一些更丰富的几何结构形态的信息，从而克服传统Gibbs先验信 息量相关较少的缺点，因此能够对PET图像重建的病态问题提供更为有效的正则化处理。步骤3中广义Gibbs先验公式wf中A的选取过程为首先直接对校正数据采用滤波反投影算法得到用于抛物线替换坐标下降迭代算法的初始图像；其后对该图像采用金字塔式结构 由粗略到精细进行方差分析得到初始图像中较平滑区域的方差值C7;最后取A值为方差cr的 10倍作为广义Gibbs先验公式wf在抛物线替换坐标下降迭代算法迭代求解过程中的固定参数。步骤3中全局参数评估计方法由于PET成像中重建图像的质量很大程度上由探测器探 测到的数据量counts决定。经大量实验验证，本发明中全局参数-可由经验公式-:ibg伊cw她)确定，其中c为一常数，弁co"""表示用于重建的总counts个数。本发明的模拟实验(见图3)中c 取值20， #<:0"她取值3xl05;心脏灌注真实数据实验(见图5)中c取值10， ^做船取值 1.4876x105。4、 完成重建。在步骤2和步骤3相关数据处理和模型设计的基础上，釆用抛物线替换 坐标下降迭代算法对目标函数进行二步式优化求解，具体如下第一步，首先基于上一步迭 代获得的重建图像为参照图像，进行待重建图像的逐像素点广义Gibbs先验的权值量wf计算，以作为下一步迭代之用；第二步，在第一步获取的权值量基础上，第二步中该权值量作 为常数，利用抛物线替换坐标下降迭代算法进行迭代重建。两步交替进行直至收敛，获得最 终重建图像。本发明的实验中，迭代次数为150次时均可得到好的重建图。
权利要求
1、一种在PET成像中引入广义Gibbs先验的最大后验优化图像重建方法，其特征在于包括如下步骤(1)利用PET成像设备采集成像前的探测器数据，同时获取成像设备中各种数据校正参数值和系统矩阵；(2)根据步骤1获取的成像前校正数据满足的统计特征，构建用于重建图像的数学统计模型；(3)针对步骤2中数学模型的求解，引入广义Gibbs先验，采用最大后验估计方法进行重建模型转化，得到用于获取PET重建图像的带约束目标函数的优化方程；(4)由步骤3得到的结果，基于对优化方程中全局参数选择的基础上，采用抛物线替换坐标下降算法进行迭代计算处理，重建出的图像。
2、 根据权利要求1所述的在PET成像中引入广义Gibbs先验的最大后验优 化图像重建方法，其特征在于步骤2中用到的数学统计模型为泊松分布或高斯 分布。
3、 根据权利要求1所述的在PET成像中引入广义Gibbs先验的最大后验优 化图像重建方法，其特征在于步骤3中广义Gibbs先验的具体设计过程为a、 首先选择一个包含图像中丰富的几何信息的较大方形邻域；同时设 计一个相似性测度用于比较大方形邻域内中像素点t和像素点_/处对应的小方 形邻域的相似性；b、 随后在选定的小方形邻域的内进行两像素间的灰度值比较的同时， 利用两像素间相似性来获得势能函数中的权值量。
4、 根据权利要求3所述的在PET成像中一l入广义Gibbs先验的最大后验优 化图像重建方法，其特征在于步骤b中的权值量定义为wf ^xp(-l卜((^-A化)lY/72)/、， 定义为传统Gibbs先验中像素点t和像素点7'之间的权值，通过图像域内两像点间的欧几里德距离的反比例函数确定；K和^ 则设定为以像素点A和像素点_/为中心的小方形邻域；A(KJ和A(。为此两个邻 域中所有像素灰度值数组；11.11代表此两个像素点所在区域的加权欧几里德距离； 参数A用于计算像素点间权值的指数函数同邻域相似性测度的反比例衰减关系。
5、 根据权利要求4所述的在PET成像中引入广义Gibbs先验的最大后验优化图像重建方法，其特征在于上述参数A的确定过程为首先直接对校正数据 采用滤波反投影算法得到用于抛物线替换坐标下降迭代算法的初始图像；其后对该图像采用金字塔式结构由粗略到精细进行方差分析得到初始图像中较平滑区域的方差值c ;最后取值为方差o"的倍数作为迭代求解时广义Gibbs先验公式wf中的固定参数。
6、 根据权利要求3所述的在PET成像中引入广义Gibbs先验的最大后验优 化图像重建方法，其特征在于步骤a中的相似性测度采用两像素点邻域内所有 像素点灰度值的加权欧几里德距离的反比例函数。
7、 根据权利要求1所述的在PET成像中引入广义Gibbs先验的最大后验优 化图像重建方法，其特征在于步骤4采用抛物线替换坐标下降迭代处理的具体过程如下第一步，首先基于上一步迭代获得的重建图像为参照图像，获得待重建图像的每个像素点处广义Gibbs先验的权值量，以作为下一步迭代之用；第二 步，在第一步获取的权值量基础上，利用抛物线替换坐标下降迭代算法进行迭代 重建；第三步，交替进行第一、二步直至收敛，获得最终重建图像。
全文摘要
本发明公开了一种在PET成像中引入广义Gibbs先验的最大后验优化图像重建方法，包括如下步骤(1)利用PET成像设备采集成像前的探测器数据，同时获取成像设备中各种数据校正参数值和系统矩阵；(2)根据步骤1获取的成像前校正数据满足的统计特征，构建用于重建图像的数学统计模型；(3)针对步骤2中数学模型的求解，引入广义Gibbs先验，采用最大后验估计方法进行重建模型转化，得到用于获取PET重建图像的带约束目标函数的优化方程；(4)由步骤3得到的结果，基于对优化方程中全局参数选择的基础上，采用抛物线替换坐标下降算法进行迭代计算处理，重建出图像。本发明可以大幅度提高PET重建图像质量。
文档编号A61B6/00GK101156780SQ200710030079
公开日2008年4月9日 申请日期2007年9月4日 优先权日2007年9月4日
发明者冯前进, 冯衍秋, 阳 陈, 陈武凡, 马建华 申请人:陈武凡