CN107346556A

CN107346556A - 一种基于块字典学习和稀疏表达的pet图像重建方法

Info

Publication number: CN107346556A
Application number: CN201710546670.7A
Authority: CN
Inventors: 刘华锋; 朱源原
Original assignee: Zhejiang University ZJU
Current assignee: Zhejiang University ZJU
Priority date: 2017-07-06
Filing date: 2017-07-06
Publication date: 2017-11-14

Abstract

本发明公开了一种基于块字典学习和稀疏表达的PET图像重建方法，其通过块字典学习的方法从解剖图像中学习先验知识，利用稀疏表达建立了基于先验知识的惩罚项，并结合根据泊松分布模型建立的似然函数形成后验函数形式的目标方程；其中，本发明利用解剖图像和块字典学习方法训练的块字典可以代表解剖图像的局部特征，并且稀疏度相同的情况下，表现出比未采用块结构字典更好的性能。因此，本发明PET重建方法有效地克服了ML‑EM方法中出现的棋盘效应等病态解问题，同时相比于其他重建方法，也取得比较好的重建效果。

Description

一种基于块字典学习和稀疏表达的PET图像重建方法

技术领域

本发明属于PET成像技术领域，具体涉及一种基于块字典学习和稀疏表达的PET图像重建方法。

背景技术

正电子发射断层扫描(Positron emission tomography，PET)是放射型断层成像的一种，依靠进入病人体内的放射性示踪剂进行断层成像。PET图像重建主要分为解析方法和迭代方法(均为统计方法)，其中解析方法发展较早，比较著名的是滤波反投影方法(Filteredback-projection，FBP)。FBP方法基于雷登变换(Radon transform)，因为没有考虑到仪器在测量时的噪声，其重建准确率不高。随着计算机性能的提高，迭代方法成为发展热点，但是迭代方法如极大似然-期望值最大化(Maximum likelihood-expectationmaximization，ML-EM)方法在迭代次数高的时候存在病态解等问题。加入惩罚项的极大后验估计方法(Most a posterior，MAP)考虑了图像的先验知识，一定程度上克服了ML-EM方法存在的问题，于是惩罚项的设计成了学者关注的焦点。基于稀疏表达和字典学习的方法因为能够借助解剖图像如CT和MRI的先验信息，因此可以设计成惩罚项。另一方面，块字典学习相比于普通的字典学习方法，能够获得更小的表示误差和更快的汇聚速度。

发明内容

鉴于上述，本发明提出了一种基于块字典学习和稀疏表达的PET图像重建方法，即利用块字典学习方法学习字典和稀疏表达作为惩罚项(先验函数)，利用基于泊松分布假设的似然函数，形成一个最大后验估计问题。

一种基于块字典学习和稀疏表达的PET图像重建方法，包括如下步骤：

(1)利用探测器对注入有放射性示踪剂的生物组织进行探测，根据各个位置探测器对所采集到的符合计数向量，构建得到符合计数矩阵y；

(2)根据PET成像原理建立PET的测量方程如下，通过对所述测量方程引入泊松噪声约束，得到PET基于泊松分布的似然函数L(x)；

y＝Gx+r+v

其中：G为系统矩阵，x为PET浓度分布矩阵，r和v分别为关于反射符合事件和散射符合事件的测量噪声矩阵；

(3)利用解剖图像通过训练学习得到块字典D，并根据其稀疏表达形成PET基于块字典的惩罚函数R(x,α)；

(4)根据似然函数L(x)和惩罚函数R(x,α)建立PET重建目标函数Ψ(x,α)如下，通过对其进行最优化求解重建得到PET浓度分布矩阵x并依此进行PET成像；

Ψ(x,α)＝λL(x)+R(x,α)

其中：λ为权重系数，α为稀疏系数矩阵。

所述似然函数L(x)的表达式如下：

其中：y_i为符合计数矩阵y中对应第i组探测器对所采集到的符合计数向量，为符合计数向量y_i中所有元素的平均值，M为探测器对的数量。

所述惩罚函数R(x,α)的表达式如下：

其中：E_s为分割算子，E_sx为经分割后PET浓度分布矩阵x中的第s个块矩阵，μ为权重系数，稀疏系数矩阵α的大小为m×N，N＝(m-n+1)²，m为PET浓度分布矩阵x的维度，n为预设的块矩阵维度，α_s为稀疏系数矩阵α中第s列稀疏系数向量，||α_s||₀表示稀疏系数向量α_s中非零元素的个数，|| ||₂表示二范数。

所述块字典D由解剖图像通过块字典学习算法训练得到，所述块字典学习算法包括稀疏凝聚算法求解字典块结构、块形式的正交匹配追踪算法求解稀疏系数以及块形式的字典元素更新。

所述块字典学习算法基于以下函数求解得到块字典D：

s.t.||x_s||_0,d≤k且|d_j|≤s_max,j∈d

其中：Y为用于训练的解剖图像且大小为n²×N，块字典D的大小为n²×l，l为设定值即块字典D的列数，d为块字典D列聚类后得到对应列的类别标签数组，X为稀疏系数矩阵且大小为l×N，x_s为稀疏系数矩阵X中的第s列稀疏系数向量，||x_s||_0,d表示稀疏系数向量x_s中连续非0元素组的个数，k为设定的稀疏度，|d_j|为类别标签数组d中类别标签为j的个数，j为自然数且1≤j≤K，K为设定聚类后的类别数，所述连续非0元素组中的元素个数小于等于s_max，s_max为设定的元素个数上限值。

优选地，所述解剖图像Y选用与PET图像对象具有相似性的CT图像。

本发明通过块字典学习的方法从解剖图像中学习先验知识，利用稀疏表达建立了基于先验知识的惩罚项，并结合根据泊松分布模型建立的似然函数形成后验函数形式的目标方程；其中，本发明利用解剖图像和块字典学习方法训练的块字典可以代表解剖图像的局部特征，并且稀疏度相同的情况下，表现出比未采用块结构字典更好的性能。因此，本发明PET重建方法有效地克服了ML-EM方法中出现的棋盘效应等病态解问题，同时相比于其他重建方法，也取得比较好的重建效果。

附图说明

图1(a)为本发明PET图像重建方法的流程示意图。

图1(b)为本发明块字典学习方法的流程示意图。

图2为蒙特卡洛模拟Zubal胸腔数据的模板示意图。

图3(a)为Zubal胸腔数据的真值图像。

图3(b)为采用ML-EM方法对Zubal胸腔数据的重建结果示意图。

图3(c)为采用本发明方法对Zubal胸腔数据的重建结果示意图。

图3(d)为采用SPS-OS方法对Zubal胸腔数据的重建结果示意图。

图4为蒙特卡洛模拟探测目标数据的模板示意图。

图5(a)为探测目标数据的真值图像。

图5(b)为采用ML-EM方法对探测目标数据的重建结果示意图。

图5(c)为采用SPS-OS方法对探测目标数据的重建结果示意图。

图5(d)为采用本发明方法对探测目标数据的重建结果示意图。

图6(a)为对探测目标数据真值图像的聚类结果示意图。

图6(b)为对采用ML-EM方法重建图像的聚类结果示意图。

图6(c)为对采用SPS-OS方法重建图像的聚类结果示意图。

图6(d)为对采用本发明方法重建图像的聚类结果示意图。

具体实施方式

为了更为具体地描述本发明，下面结合附图及具体实施方式对本发明的技术方案进行详细说明。

如图1(a)所示，本发明基于块字典学习和稀疏表达的PET图像重建方法，包括如下步骤：

(1)利用PET成像仪器对注入放射性示踪剂的生物组织进行探测，通过各个位置探测器对收集到符合计数向量，建立符合计数矩阵y；

(2)根据PET成像原理以及PET成像仪器的系统参数，建立PET测量方程如下：

y＝Gx+r+v

其中：G为系统矩阵，由PET测量系统决定，其中的元素g_ij表示的是x中像素j被探测器i探测到的概率，r为因反射符合事件引起的测量噪声，v为因散射符合事件引起的测量噪声。

(3)根据每个探测对上符合计数事件符合泊松分布的假设，建立基于泊松分布的似然函数L(x)如下：

其中：y＝{y_i,i＝1,2…M}，是系统每个探测器对探测到的符合事件的个数和，M是探测器对的个数。每个探测器探测到的光子个数都是独立同分布，且符合伯努利过程，正弦图中每个元素都可以用一个符合泊松分布的随机变量来描述如下：

(4)用块字典学习方法根据解剖图像训练块字典，并根据其稀疏表达，形成惩罚函数R(x，α)如下：

其中：α_s为第s块区域的稀疏表达系数，μ为稀疏项的权重系数。每个字典元素的长度为n²，将图像分成n×n的小块，并且用一个向量表示。用一个分割算子E从图像向量中提取区域块。将图像排列成的向量，并且将其分成N个重叠的区块，则代表每一个小块向量可以表示成p_s＝E_sx，D为由解剖图像通过块字典学习方法训练得到的块字典。

块字典D由解剖图像用块字典学习算法得到，如图1(b)所示，块字典学习算法包括稀疏凝聚算法求解字典块结构，块形式的正交匹配追踪算法求解稀疏系数，块形式的字典元素更新方法；块字典学习算法主要基于以下函数表达：

s.t.||x_i||_0,d≤k,1≤i≤N,|d_j|≤s_max,j∈d

其中：Y为输入的训练信号，X为稀疏系数矩阵，x_i表示稀疏系数矩阵X的第i列，k表示块稀疏度。s_max为字典中块的大小的上限，即块中最多含有的字典元素，|d_j|表示字典中第j块中包含的字典原子个数，即块的大小。

块结构d_j由稀疏凝聚算法(SAC)得到，稀疏系数矩阵X由块正交匹配追踪算法(BOMP)得到，然后以块形式由k次奇异值分解方法更新字典元素，其中稀疏凝聚算法将对字典元素的分组转换为对相应的系数矩阵的行的分组，将具有相似非零元素分布的行归为一组等效于将字典中相似的列归为一块：

其中，w_j(X,d)表示系数矩阵X中块d_j中的字典元素所对应的行是非零值的列集；设置最大块尺寸s,求解以下式子：

其中，|w_j|表示的w_j的大小，B是指字典当前的块数量；在满足|d_j|≤s的前提下，合并相似的块，每一次迭代中，找到一对块满足：

合并和即将w_j1←{w_j1∪w_j2}和重复这一步，直到没有块可以合并。

然后，利用块形式的正交匹配追踪算法更新稀疏系数矩阵，具体地：设置块稀疏度和容忍误差，输入字典D和训练信号矩阵Y，按照以下所示步骤更新系数；块稀疏度和容忍误差都需要根据具体的图像而定，块稀疏度一般设置为2或3，而容忍误差一般为图像中像素均值的0.1左右。

其中：r为信号y和稀疏表达的残差；上述算法得到的稀疏系数矩阵X以及块结构d，结合原字典D，对于字典中的某块d_j，将误差函数写成然后用SVD奇异值分解这里取秩为|d_j|的对角矩阵即得到块d_j更新后的字典元素和相应系数如下；将字典中所有的块更新之后，即得到新的字典D。

(5)将似然函数和惩罚函数相结合，形成新的目标函数如下：

Ψ(x,α)＝λL(x)+R(x,α)

其中：λ为惩罚函数和似然函数之间的关系系数，α为稀疏系数矩阵。

(6)对目标函数Ψ(x,α)最优化求解即得到PET浓度分布矩阵x，目标函数Ψ(x,α)的优化求解表达式如下：

针对上式，首先用期望值最大化的思想，引入一个隐含变量c_ij，表示像素j对于探测器i贡献的符合事件，g_ij为系统矩阵的元素，而x_j表示像素j发出的光子对数，进一步表达如下：

计算出c_ij的期望值再将其带回目标方程得到新的目标方程

将新的目标函数的第二项用可分离的凸替代函数替换：

进一步将[E_sx]_l分解成含x_j的形式：

其中，为第k次迭代中对重建图像第j个元素的估计值，e_s,lj为E_s,l的第j个元素；利用以下不等式，即得到替代的凸函数。

新的目标方程变成：

在以上的目标方程中，对x_j求偏导，即得到的交替计算EM中的期望值和最值直至收敛。

以下我们通过对蒙特卡洛模拟数据进行实验来验证本发明PET重建的准确性。实验运行环境为:8G内存，2.70GHz，64位操作系统，CPU为Intel Core i5；所模拟的PET扫描仪型号为Hamamatsu SHR-22000，设定的放射性核素及药物为18F-FDG，设置正弦图为128个投影角度在每个角度下128条射束采集到的数据结果，系统矩阵G∈R^16384×16384；实验中使1×10⁶计数率下的投影数据作为实验数据。

将本发明PET重建方法和经典的ML-EM方法以及SPS-OS方法进行了比较，为控制结果的可比性，使用了相同的测量数据y和系统矩阵G。图2为实验用的Zubal胸腔数据的模板，除背景外将图像分为三个感兴趣的区域ROI；图3(a)为真值图，图3(b)～图3(d)分别为ML-EM，本发明方法以及SPS-OS方法的重建效果。计算每个ROI区域相当于真值图像的偏差和方差，表1是分析结果：

表1

从结果来看，本发明重建框架不仅比传统ML-EM的方法重建精度更高，更平滑；相比于SPS-OS方法，本发明重建框架的重建效果也有较好的表现。

图4为实验用的6个尺寸不同的探测目标的模板，图5(a)为真值图，图5(b)～图5(d)分别为ML-EM，SPS-OS和本发明方法的重建效果图；图6(a)～图6(d)是相对应的经过k均值聚类后的结果，用Jaccard指数来分析每一个探测目标的探测效果，如表2所示：

表2

从结果来看，本发明重建框架的探测性能与SPS-OS的探测性能相当，均优于传统的ML-EM方法。虽然传统ML-EM方法能探测到最小的目标，但是从重建的图中可以看到，目标周围有许多噪声，不适合实际中使用。

上述对实施例的描述是为便于本技术领域的普通技术人员能理解和应用本发明。熟悉本领域技术的人员显然可以容易地对上述实施例做出各种修改，并把在此说明的一般原理应用到其他实施例中而不必经过创造性的劳动。因此，本发明不限于上述实施例，本领域技术人员根据本发明的揭示，对于本发明做出的改进和修改都应该在本发明的保护范围之内。

Claims

1.一种基于块字典学习和稀疏表达的PET图像重建方法，包括如下步骤：

y＝Gx+r+v

Ψ(x,α)＝λL(x)+R(x,α)

其中：λ为权重系数，α为稀疏系数矩阵。

2.根据权利要求1所述的PET图像重建方法，其特征在于：所述似然函数L(x)的表达式如下：

<mfenced open = "" close = ""> <mtable> <mtr> <mtd> <mrow> <mi>L</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>M</mi> </munderover> <msub> <mover> <mi>y</mi> <mo>&OverBar;</mo> </mover> <mi>i</mi> </msub> <mo>-</mo> <msub> <mi>y</mi> <mi>i</mi> </msub> <mi>l</mi> <mi>o</mi> <mi>g</mi> <mrow> <mo>(</mo> <msub> <mover> <mi>y</mi> <mo>&OverBar;</mo> </mover> <mi>i</mi> </msub> <mo>)</mo> </mrow> </mrow> </mtd> <mtd> <mrow> <mi>s</mi> <mo>.</mo> <mi>t</mi> <mo>.</mo> </mrow> </mtd> <mtd> <mrow> <mi>y</mi> <mo>=</mo> <mi>G</mi> <mi>x</mi> <mo>+</mo> <mi>r</mi> <mo>+</mo> <mi>v</mi> </mrow> </mtd> </mtr> </mtable> </mfenced>

3.根据权利要求1所述的PET图像重建方法，其特征在于：所述惩罚函数R(x,α)的表达式如下：

<mrow> <mi>R</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>&alpha;</mi> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&Sigma;</mo> <mrow> <mi>s</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>N</mi> </munderover> <mo>|</mo> <mo>|</mo> <msub> <mi>E</mi> <mi>s</mi> </msub> <mi>x</mi> <mo>-</mo> <msub> <mi>D&alpha;</mi> <mi>s</mi> </msub> <mo>|</mo> <msubsup> <mo>|</mo> <mn>2</mn> <mn>2</mn> </msubsup> <mo>+</mo> <mi>&mu;</mi> <mo>|</mo> <mo>|</mo> <msub> <mi>&alpha;</mi> <mi>s</mi> </msub> <mo>|</mo> <msub> <mo>|</mo> <mn>0</mn> </msub> </mrow>

其中：E_s为分割算子，E_sx为经分割后PET浓度分布矩阵x中的第s个块矩阵，μ为权重系数，稀疏系数矩阵α的大小为m×N，N＝(m-n+1)²，m为PET浓度分布矩阵x的维度，n为预设的块矩阵维度，α_s为稀疏系数矩阵α中第s列稀疏系数向量，||α_s||₀表示稀疏系数向量α_s中非零元素的个数，||||₂表示二范数。

4.根据权利要求3所述的PET图像重建方法，其特征在于：所述块字典D由解剖图像通过块字典学习算法训练得到，所述块字典学习算法包括稀疏凝聚算法求解字典块结构、块形式的正交匹配追踪算法求解稀疏系数以及块形式的字典元素更新。

5.根据权利要求4所述的PET图像重建方法，其特征在于：所述块字典学习算法基于以下函数求解得到块字典D：

s.t. ||x_s||_0,d≤k且|d_j|≤s_max,j∈d

其中：Y为用于训练的解剖图像，块字典D的大小为n²×l，l为设定值即块字典D的列数，d为块字典D列聚类后得到对应列的类别标签数组，X为稀疏系数矩阵且大小为l×N，x_s为稀疏系数矩阵X中的第s列稀疏系数向量，||x_s||_0,d表示稀疏系数向量x_s中连续非0元素组的个数，k为设定的稀疏度，|d_j|为类别标签数组d中类别标签为j的个数，j为自然数且1≤j≤K，K为设定聚类后的类别数，所述连续非0元素组中的元素个数小于等于s_max，s_max为设定的元素个数上限值。

6.根据权利要求5所述的PET图像重建方法，其特征在于：所述解剖图像Y选用与PET图像对象具有相似性的CT图像。