基于隐含活动轮廓先验的贝叶斯图像重建方法

专利名称：基于隐含活动轮廓先验的贝叶斯图像重建方法
技术领域：
本发明涉及一种图像重建方法，尤其涉及一种基于隐含活动轮廓先验的贝叶斯图像重建方法。
背景技术：
正电子发射计算机断层显像(Positron emission tomography，PET)重建方法的研究近年来越来越受到人们的重视。现有的图像重建方法主要为解析法和迭代法。解析法主要是以中心切片定理为基础的滤波反投影方法，其特点是速度快，但当测量噪声较大或采样不充分时，这类方法的成像效果不理想。迭代法主要包括代数法，最大似然法等等。迭代法的特点是可以根据具体成像条件引入与空间几何有关的或与测量值大小有关的约束和条件因子，但迭代法收敛速度慢，运算时间长。正则化技术或者最大后验重建(Maximum A Posteriori，MAP)方法的使用是必要的。最大后验重建方法对于提高图像的空间分辨率和噪声特性的效果是非常显著的，但经常会出现过平滑现象，也就是图像在被平滑的同时，边缘也被平滑掉了，因为噪声和边缘都是高频分量。而且使用这种方法的前提是要选择一种合适的先验，不恰当的先验分布会导致完全错误的重建结果。由于先验函数中存在超参数估计问题，有文献表示用其它的高质量成像方法如CT、MRI所提供的人体形态学结构图像来构造一种分割模板先验，以保证先验的可靠性，重建结果采用基于动态模拟来计算，这种方法可以很好解决超参数的选择，但这种方法也必须事先知道CT、MRI重构图像，在通常情况下要病人同时进行PET，CT或MRI扫描会增加病人的痛苦，因此这种方法是不实用的。

发明内容
本发明提供一种重建后的图像即可以保持高的分辨率，又能够降低噪声特性且重建后图像的边缘清晰的基于隐含活动轮廓先验的贝叶斯图像重建方法。
本发明采用如下技术方案一种基于隐含活动轮廓先验的贝叶斯图像重建方法
1)获取投影数据，根据待重建图像的尺寸要求，确定初始图像范围，给定初始的灰度值大于1，并将2维图像变换成1维向量进行计算，2)根据投影数据规模和要求成像图像x的大小，计算系统概率矩阵P，3)将系统概率阵和初始图像x相乘，得到前向投影a，4)将投影数据除以前向投影a，得到对投影数据的校正值c，5)将系统概率矩阵P乘以投影数据的校正值c，得到图像成像迭代过程中的修正值xd，6)将高斯函数对初始图像作卷积计算，然后求这个计算结果的梯度模，再将这个计算结果除以0～1之间的任意一个数β1的平方，得到第一个数e，将第一个数e加1的结果的倒数乘以图像的梯度模|x|，得到一个扩散速度g，将图像的梯度x除以图像的梯度模|x|，并对这个结果求散度，这个散度和气球力F0之和乘以上述扩散速度g，得到能量函数 用差分方法将这个能量函数离散化，7)将系统概率矩阵P对它的每一行求和，把这个对行的和的β倍加上能量函数，得到一个权值，将初始图像乘以修正值xd除以这个权值，得到了一个重建图像，再将此重建图像作为下一次迭代的初始图像，返回到第3步，反复迭代直到重建后的图像收敛。
与现有技术相比，本发明具有如下优点本发明利用一种隐含活动轮廓模型作为一个新的能量函数，并通过使这些能量函数最小这种先验信息来进行Bayesian重建。由于这种方法利用了图像的边缘梯度信息，重建后的图像即可以保持高的分辨率，又能够降低噪声特性，而且重建后图像的边缘得到了有效地保持，不需要病人再进行PET，CT或MRI扫描，可以减轻病人的痛苦。


图1为用来测试成像方法的胸腔模板图像。
图2为用来测试重建方法的投影数据。
图3为用来测试重建方法的含泊松噪声的投影数据。
图4为用现有的贝叶斯(Bayesian)成像方法成像后的结果，此时投影数据不含有泊松噪声。
图5为用现有的贝叶斯(Bayesian)成像方法成像后的结果，此时投影数据含有泊松噪声。
图6为用本发明的能量函数替代现有的贝叶斯(Bayesian)成像方法中的能量函数成像后的结果，此时投影数据不含有泊松噪声。
图7为用本发明的能量函数替代现有的贝叶斯(Bayesian)成像方法中的能量函数成像后的结果，此时投影数据含有泊松噪声。
具体实施例方式
一种基于隐含活动轮廓先验的贝叶斯图像重建方法1)获取投影数据，根据待重建图像的尺寸要求，确定初始图像范围，给定初始的灰度值大于1，并将2维图像变换成1维向量进行计算，2)根据投影数据规模和要求成像图像x的大小，计算系统概率矩阵P，3)将系统概率阵和初始图像x相乘，得到前向投影a，4)将投影数据除以前向投影a，得到对投影数据的校正值c，5)将系统概率矩阵P乘以投影数据的校正值c，得到图像成像迭代过程中的修正值xd，6)高斯函数对初始图像作卷积计算，然后求这个计算结果的梯度模，再将这个计算结果除以0～1之间的任意一个数β1的平方，得到第一个数e，将第一个数e加1的结果的倒数乘以图像的梯度模|x|，得到一个扩散速度g，将图像的梯度x除以图像的梯度模|x|，并对这个结果求散度，这个散度和气球力F0之和乘以上述扩散速度g，得到能量函数 用差分方法将这个能量函数离散化，7)将系统概率矩阵P对它的每一行求和，把这个对行的和的β倍加上能量函数，得到一个权值，将初始图像乘以修正值xd除以这个权值，得到了一个重建图像，再将此重建图像作为下一次迭代的初始图像，返回到第3步，反复迭代直到重建后的图像收敛，上述投影数据的获取是从正电子发射计算机断层显像扫描仪上获取，或者是从仿真模板图像进行雷当(Radon)变换，得到的投影数据。
实施例2本发明是通过对已有PET重建方法进行改进后得到的，具体实施方案的内容如下1.已有的Bayesian重建方法一个正则化的最大似然估计的优化问题为maxx≥0logP(Y=y|x)+φ(x)---(1)]]>
这里φ(x)即为惩罚项，它的意义是指图像的先验概率，因此公式(1)的具体形式可写成maxx≥0logP(Y=y|x)+logP(x)---(2)]]>这里P(x)即表示先验概率，根据Bayesian定理有P(x|Y=y)=P(Y=y|x)·P(x)P(Y=y)---(3)]]>我们注意到由于P(Y＝y)与待估参数x无关，公式(2)等价于最大化后验概率P(x|Y＝y)，这是最大后验估计的含义，亦称Bayesian估计Bayesian重建方法首先要解决的问题就是如何选择一个合适的先验，先验模型主要是指图像的光滑模型，图像的光滑性主要表现为一种局部性质，即邻域内的象素间的相互作用，这种作用越强，相邻近的象素值越均匀，图像越光滑，因此象素取某一值仅与其邻域象素有关，我们可以通过给予邻域象素值较均匀的构型以大的概率来定量反映这种光滑特性，这就构成Markov先验概率模型。一般所采用的先验都是基于某种邻域结构的Markov随机场，但由于Markov随机场没有显式地给出图像的总体概率分布，因此大多用Gibbs随机场来替代，这是因为Gibbs可以显式地给出图像的概率分布，而且在一定的条件下Markov随机场等价于Gibbs随机场。Gibbs随机场给出的图像概率分布为P(x)＝Z-1exp(-βV(x)) (4)这里V(x)=Σj∈SΣi∈NU(xi,xj)---(5)]]>称为系统的能量函数，Z是归一化因子，也称分隔函数Z=Σxexp(-βV(x))---(6)]]>有了Gibbs先验概率模型，就可以实现最大后验估计，结合Poisson数据模型，最大后验估计可以重新写成maxx≥0Σi(yilogΣjpijxj-Σjpijxj-logyi!)-βV(x)---(7)]]>这里能量函数为 常用的结构主要是八邻域结构，对于八邻域结构，一般我们还要对最邻近象素与对角相邻象素通过引入一定的权值加以区别，如果是八邻域结构，能量函数表示为 其中wij为权值，最邻近点被赋值为1，对角邻点则被赋值 函数(x)一般为偶函数。例如(x)＝log(cosh(x)) (10)利用期望最大方法对公式(7)求得的Bayesian重建方法为xj(k+1)=xj(k)(ΣiyipijΣjpijxj(k))/(Σipij+β∂V(x(k))∂xj),j=1,2,...n---(11)]]>其中yi表示第i个探测器所探到的光子数，0≤i≤m，m为探测器总数；xj表示第j个象素处发出的光子数。xj≥0，0≤j≤n，n为象素数；pij表示第j个象素处发出的光子能被第i个探测器检测到的概率。
2.基于隐含活动轮廓先验的Bayesian图像的重建方法活动轮廓模型是一种有效的图像分割方法，基于活动轮廓的图像分割实质上就是用活动轮廓逼近物体的边缘，此过程可以通过能量最小来实现，其形变过程就是活动轮廓在外部能量和内能(内力)的作用下向物体边缘靠近的过程，外力推动活动轮廓向着物体边缘运动，而内力则保持活动轮廓的光滑性和拓扑性。当能量最小时，活动轮廓收敛到所要检测的物体边缘。由于这种方法同时考虑了几何约束条件、图像数据、轮廓形状有关的能量最小等约束条件，所以用这种方法能得到满意的分割效果。
由公式(3)我们可以得到p(x|y)∝p(x)p(y|x) (12)本发明定义Bayesian重建方法中的一个能量函数为V(x)＝V内(x)+V外(x)(13)关于图像x的先验概率分布为 这里Z内为归一化因子，也就是内力分隔函数，对于一个待重建的图像x的观测数据Y的似然为
这里Z外为归一化因子，也就是外力分隔函数，因此 在这个新发明里，我们给出了一种隐含的活动轮廓模型，并利用它来表示能量函数V(x)。
本发明首先定义一个内部的能量函数为V(x)=12∫Ω2λ(∂x∂t)2dx---(17)]]>其中λ为系数，外部的能量函数假设为0。此时∂x∂t=g(|▿Gσ*x|)|▿x|(div▿x|▿x|+F0)---(18)]]>用有限差分方法求解(18)式，并将公式(18)代入公式(17)，并将其离散化，得到新的能量函数，将这个能量函数代入公式(11)得到本发明中的Bayesian重建方法。
利用这种模型构成的能量函数作为Bayesian先验，可使得函数x的水平曲线沿着垂直x方向以g(|Gσ*x|)速度扩散，由此可见，图像在边缘处，也就是梯度最大处的邻域实行弱光滑，对边缘点本身实行较小程度的光滑，而对其它区域实行强光滑。F0为常数，也称气球力，用来平衡内外力的作用，它的取值范围在0～10之间中的任意数，使重建后的图像更清晰，其中Gσ为方差为σ的高斯函数。这里g(x,y)=11+(|▿Gσ(x,y)*I(x,y)|/β1)2---(19)]]>在活动轮廓模型中，g(i，j)为时间独立的停止因子。这种活动轮廓作为Bayesian重建方法中的先验函数的能量函数，重建后的图像质量较好，而且边缘定位较准确。
3.基于隐含活动轮廓先验的Bayesian图像重建方法的实验结果本发明用一个计算机仿真的PET胸腔幻影模板来验证本发明方法的可靠性。图1显示了这个胸腔幻影模板，模板大小为128×128象素矩阵，数据规模为185×180，即180个投影角度，每一个角度上有185条平行线，这里使平行线的间距与图像象素的边长相等，以便系统概率阵的确定。图2和图3分别表示投影数据不包含Poisson噪声和包含Poisson噪声时的情况。
图4和图5分别表示投影数据不包含噪声和包含噪声时用现有Bayesian方法对模板进行重建后的结果。图6和图7分别表示投影数据不包含噪声和包含噪声时用本发明方法对模板进行重建后的结果。
权利要求
1.一种基于隐含活动轮廓先验的贝叶斯图像重建方法，其特征在于采用以下步骤1)获取投影数据，根据待重建图像的尺寸要求，确定初始图像范围，给定初始的灰度值大于1，并将2维图像变换成1维向量进行计算，2)根据投影数据规模和要求成像图像x的大小，计算系统概率矩阵P，3)将系统概率阵和初始图像x相乘，得到前向投影a，4)将投影数据除以前向投影a，得到对投影数据的校正值c，5)将系统概率矩阵P乘以投影数据的校正值c，得到图像成像迭代过程中的修正值xd，6)将高斯函数对初始图像作卷积计算，然后求这个计算结果的梯度模，再将这个计算结果除以0～1之间的任意一个数β1的平方，得到第一个数e，将第一个数e加1的结果的倒数乘以图像的梯度模|x|，得到一个扩散速度g，将图像的梯度x除以图像的梯度模|x|，并对这个结果求散度，这个散度和气球力F0之和乘以上述扩散速度g，得到能量函数 用差分方法将这个能量函数离散化，7)将系统概率矩阵P对它的每一行求和，把这个对行的和的β倍加上能量函数，得到一个权值，将初始图像乘以修正值xd除以这个权值，得到了一个重建图像，再将此重建图像作为下一次迭代的初始图像，返回到第3步，反复迭代直到重建后的图像收敛。
2.根据权利要求1所述的基于隐含活动轮廓先验的贝叶斯图像重建方法，其特征在于投影数据的获取是从正电子发射计算机断层显像扫描仪上获取的。
3.根据权利要求1所述的基于隐含活动轮廓先验的贝叶斯图像重建方法，其特征在于投影数据的获取是从仿真模板图像进行雷当变换，得到的投影数据。
全文摘要
本发明公开了一种基于隐含活动轮廓先验的贝叶斯图像重建方法，首先获取投影数据，确定初始图像范围，计算系统概率矩阵，得到前向投影，再将投影数据除以前向投影，得到对投影数据的校正值，乘以系统概率矩阵，得到图像成像迭代过程中的修正值，然后经计算后得到能量函数，用差分方法将这个能量函数离散化，最后将系统概率矩阵对它的每一行求和乘以β加上能量函数，得到权值，将初始图像乘以修正值除以权值，得到重建图像，作为下一次迭代的初始图像，反复迭代直到重建后的图像收敛，本发明具有保持高分辨率，降低噪声，且保持图像边缘，减轻病人的痛苦等优点。
文档编号A61B8/13GK1640362SQ20051003762
公开日2005年7月20日 申请日期2005年1月6日 优先权日2005年1月6日
发明者朱宏擎, 舒华忠, 周键, 罗立民, 李松毅 申请人:东南大学