一种基于tv和稀疏约束的动态pet图像和示踪动力学参数同步重建方法

文档序号:10535884阅读:492来源:国知局
一种基于tv和稀疏约束的动态pet图像和示踪动力学参数同步重建方法
【专利摘要】本发明公开了一种基于TV和稀疏约束的动态PET图像和示踪动力学参数同步重建方法,其通过建立了对动态PET图像和动力学参数同时估计的数学模型,加入了基于稀疏约束的字典和全变分算子来进行整个重建的过程;其中加入了全变分算子的动态PET序列的重建子问题采用ADMM算法进行迭代优化求解;结合稀疏约束的字典的PET动力学参数的估计的子问题采用软阈值迭代算法来进行求解。本发明有效解决了同时对动态PET图像序列和动力学参数进行估计的问题,并且加入全变分算子改善了PET图像重建过程中结果分辨率较低和易受噪声干扰的问题,且与其他单独重建动态PET图像或估计动力学参数的算法相比,本发明也都能获得较好的重建结果。
【专利说明】
一种基于TV和稀疏约束的动态PET图像和示踪动力学参数同 步重建方法
技术领域
[0001] 本发明属于PET成像技术领域,具体涉及一种基于TV(total variation,全变分) 和稀疏约束的动态PET图像和示踪动力学参数同步重建方法。
【背景技术】
[0002] 正电子发射断层成像(Positron emission tomography,PET)是核医学成像的一 种,它利用示踪原理来显示生物活体的新陈代谢特性,在医学研究和临床诊断中都有着一 定的作用。相对于只针对静态时间窗内PET的图像重建分析,动态PET还能针对生物体的组 织或器官的真实代谢水平进行定量的分析。从分子水平上观察细胞的代谢活动,为疾病的 早期诊断和预防提供了有效依据;其中,尤以对癌症的发现、诊断等有着重要的应用。
[0003] 传统上,对PET放射性浓度的重建往往采用解析法和迭代统计法的方法,前一类的 方法主要基于雷登变换和中央切片定理对投影数据进行滤波反投影(Foureir Back Pr〇jeCti〇n,FBP)或改进的逆卷积变换从而获得所需的图像。由于其重建的结果中存在伪 影,图像品质精确度并不高。后一类的方法中常见的有ML-EM(最大似然-期望最大化)、MAP (最大后验)等,这些算法一般都只针对同一帧的PET图像信息,忽略了其在时间上提供的生 理信息,一定程度上增加了结果中的噪声。
[0004] 关于示踪动力学参数的重建,一般可将其分为间接重建和直接重建两类。前者需 要先估算出PET放射性浓度分布,再基于此分布对动力学参数进行计算;而后者则是直接从 PET仪器获得到的数据中估计出动力学参数。尽管直接重建的方法实现起来比间接重建困 难,由于减少了二次估计可能造成的误差以及能够对噪声有一个更准确的模拟,往往能获 得更加准确的结果。
[0005] 因此,在直接估计出示踪动力学参数的基础上同时获得准确的动态PET图像序列 是研究PET成像的一个非常关键的问题。

【发明内容】

[0006] 针对现有技术所存在的上述技术问题,本发明提供了一种基于TV和稀疏约束的动 态PET图像和示踪动力学参数同步重建方法,能够同时获得高质量的动态PET重建图像序列 和动力学参数的估计。
[0007] 一种基于TV和稀疏约束的动态PET图像和示踪动力学参数同步重建方法,包括如 下步骤:
[0008] (1)利用探测器对注入有放射性药剂的生物组织进行探测,动态采集得到M组符合 计数向量,进而组建PET的符合计数矩阵Y,M为大于1的自然数;
[0009] (2)通过使动态PET图像序列组合成PET浓度分布矩阵X,根据TOT成像原理,建立 PET测量方程;
[001 0] (3)通过对PET测量方程引入全变分约束,得到基于TV的PET图像重建模型L(X)如 下:
[0012] 其中:a为权重系数,G为系统矩阵,TV(X)为关于X的全变分正则项,I I I |2为L2范 数;
[0013] (4)利用房室模型匹配估计示踪动力学参数,建立关于X和〇的同步重建模型S(X, 〇)如下:
[0015] 其中:y为权重系数,W为字典矩阵,O为示踪动力学参数矩阵,T表示转置,| | | U 为L1范数,|| | |F为F范数;
[0016] (5)将上述两个模型L(X)和S(X,〇 )相结合得到最终的同步重建模型LS(X,〇 )如 下:
[0018]其中:y为权重系数;
[0019] (6)对所述的同步重建模型LS(X,〇)进行最优化求解后即得到PET浓度分布矩阵X 和示踪动力学参数矩阵〇。
[0020] 所述的符合计数矩阵Y即由M组符合计数向量组成,PET浓度分布矩阵X即由M组PET 浓度分布向量组成,所述的PET浓度分布向量即对应某一时刻的一帧PET图像浓度数据。 [0021] 所述PET测量方程的表达式如下:
[0022] Y=GX+R+T
[0023] 其中:R和T分别为关于反射符合事件和散射符合事件的测量噪声矩阵。
[0024] 所述全变分正则项TV(X)的表达式如下:
[0026]其中:为PET浓度分布矩阵X中对应的浓度值,表示第m帧PET图像中的第 j个体素,(Dxm)j为关于的二维差分向量,该向量第一行元素值为与其右边体素的浓 度差值,第二行元素值为1与其下边体素的浓度差值,K为一帧PET图像的体素总个数。 [0027]所述字典矩阵W的表达式如下:
5
[0031 ]其中:(Mt)和CA)分别为t时刻和T时刻放射性药剂在血浆中的浓度值,'和给 分别为关于第m组符合计数向量采集的开始时间和结束时间,0。对应为第c个房室组织指数 函数的系数,m和c均为自然数且1 < c <N,N为大于1的自然数;9i~%的取值为在区 间[0min,0max]内按指数级间隔进行选取,0min和0 max分别为系数的上下限阈值,t和t均表示时 间。
[0032] 优选地,所述的步骤(6)中米用ADMM(Alternating Direction Method of Multipliers,交替方向乘子算法)结合软阈值(Soft-Thresholding)的迭代优化算法对同 步重建模型LS(X,〇)进行最优化求解;其中,ADMM针对PET浓度分布矩阵X进行迭代优化求 解,软阈值针对示踪动力学参数矩阵〇迭代优化求解。
[0033] 本发明通过建立了对动态PET图像和动力学参数同时估计的数学模型,加入了基 于稀疏约束的字典和全变分算子来进行整个重建的过程;其中加入了全变分算子的动态 PET序列的重建子问题采用ADMM算法进行迭代优化求解;结合稀疏约束的字典的PET动力学 参数的估计的子问题采用软阈值迭代算法来进行求解。本发明有效解决了同时对动态PET 图像序列和动力学参数进行估计的问题,并且加入全变分算子改善了在进行PET图像重建 的过程中结果分辨率较低和易受噪声干扰的问题,并且与其他单独重建动态PET图像或估 计动力学参数的算法相比,本发明也都能获得较好的重建结果。
【附图说明】
[0034]图1为脑部phantom模板图。
[0035]图2(a)为采用本发明方法在计数率5 X104下重建蒙特卡洛模拟脑部phantom数据 的第9帧PET图像。
[0036]图2(b)为采用本发明方法在计数率1 X 105下重建蒙特卡洛模拟脑部phantom数据 的第9帧PET图像。
[0037]图2(c)为采用本发明方法在计数率5 X105下重建蒙特卡洛模拟脑部phantom数据 的第9帧PET图像。
[0038]图2(d)为采用本发明方法在计数率1 X 106下重建蒙特卡洛模拟脑部phantom数据 的第9帧PET图像。
[0039]图3(a)为采用ML-EM方法在计数率5 X 104下重建蒙特卡洛模拟脑部phantom数据 的第9帧PET图像。
[0040]图3(b)为采用ML-EM方法在计数率1 X 105下重建蒙特卡洛模拟脑部phantom数据 的第9帧PET图像。
[00411图3(c)为采用ML-EM方法在计数率5 X 105下重建蒙特卡洛模拟脑部phantom数据 的第9帧PET图像。
[0042]图3(d)为采用ML-EM方法在计数率1 X 106下重建蒙特卡洛模拟脑部phantom数据 的第9帧PET图像。
[0043]图4(a)为蒙特卡洛模拟脑部phantom数据的动力学参数VD真值图像。
[0044] 图4(b)为采用本发明方法在计数率5 X104下重建蒙特卡洛模拟脑部phantom数据 的动力学参数Vd结果图像。
[0045] 图4(c)为采用本发明方法在计数率1 X 105下重建蒙特卡洛模拟脑部phantom数据 的动力学参数Vd结果图像。
[0046] 图4(d)为采用本发明方法在计数率5 X105下重建蒙特卡洛模拟脑部phantom数据 的动力学参数Vd结果图像。
[0047]图4(e)为采用本发明方法在计数率1 X 106下重建蒙特卡洛模拟脑部phantom数据 的动力学参数Vd结果图像。
[0048] 图5为本发明重建方法的步骤流程示意图。
【具体实施方式】
[0049] 为了更为具体地描述本发明,下面结合附图及【具体实施方式】对本发明的技术方案 进行详细说明。
[0050] 如图5所示,本发明基于全变分约束和稀疏约束同时重建动态PET图像和示踪动力 学参数的方法,包括如下步骤:
[0051] S1.根据动态PET探测的原理整理测量数据矩阵y和相应的系统矩阵G;
[0052] S2.建立字典矩阵若是模拟数据,则根据模拟的示踪同位核素的半衰期设置指 数参数9和探测器扫描时间间隔建立相应的字典矩阵;若是真实数据,则按照实际操作中使 用的示踪同位核素和实验设置建立相应的字典矩阵。
[0053]根据房室模型可以解出组织中的放射性强度随时间的变化,即时间-放射性函数 (time-activity curve,TAC)是输入组织的TAC和一个5(t)函数以及一组指数函数的卷积, 即是所有组织房室的TAC之和。因此,一个PET示踪房室模型可以用一组一阶线性方程来表 示,每个方程表示一个不同的房室中的示踪剂放射性强度,即:
[0055] 其中:CT(t)表示t时刻组织中的放射性强度。N表示不同房室的总数,也(t)表示第c 种房室对应的放射性强度随时间变化的函数,巾。是M t)相对应的系数。0。是对应房室的指 数函数中的参数。(Mt)表示t时刻输入的放射性强度。
[0056] 重建的每一帧的动态PET图像可看作在这一帧与上一帧的时间间隔内对房室组织 时间-放射性曲线的积分,BP :
[0058]其中,x#是PET图像序列中第m帧的像素j所对应的放射性强度,和分别是第 m帧的开始时间和结束时间,CjG)是方式组织的像素j处对应的TAC。基于此关系,我们将 一组TAC对应的方程组扩展为一个字典矩阵如下:
[0063] S3.初始化,设置_ -的权重系数T,稀疏约束I I巾I li的权重系数4和 软阈值算法中的步长〇。其中:
[0064] 权重系数a在[27'28'5]左右,这个是根据原来已有的算法给出的参考范围进行调 整得到的。权重系数T 一般在[1,5]左右的范围内,这个是基于实验的结果获得的。权重系 数y -般可在[0.01,0. 15]的范围内进行选择得到合适的值(根据参考文献Positron Emission Tomography Compartmental Models&colon;A Basis Pursuit Strategy for Kinetic Modeling中的实验结果),当噪声增加的时候,算法不能准确地从字典中确定出正 确的时间放射性曲线,此时需要适当增大y的取值,增强稀疏约束。
[0065] S4.基于动力学参数巾的第k次迭代后的估计值(^更新动态PET图像序列3+1:提 取出目标参数中仅包括变量x的部分获得以下子优化问题:
[0067]其中:全变分正则项可通过以下定义计算:
[0069] 基于加TV的ADMM算法,上述问题可归纳为:
[0071]以上约束条件中的第一、二个可进一步统一为: H、c =蠢
[0072] , s.t 轟=(:i + ' 為y +
[0073] 此时,动态PET图像x重建的子问题可以转换为与其相对应的增广拉格朗日方程:
[0075] 其中:ujm和入是包含了拉格朗日乘数的向量,在迭代的过程中会随着《 jdPx的更 新而更新,是相应的惩罚参数。接下来将分别对和x进行更新,从而获得此问题的最 优解。
[0076] 4.1对于一个固定的x的估计值以及Ujm和A,《 jm的最优解能够通过仅优化增广拉 格朗日方程中与w jm相关的部分,根据二维收缩公式获得:
[0079] 4.2对于一个固定的的估计值以及〇_和\,可将原公式近似为一个关于x的二 次方程,再使用一个最陡下降步而确定其更新值。最陡下降步的步长是由从Barz i 1 ai Borwein步长开始的逆向跟踪非单调线搜索的方法计算出来。
[0080] 彳^基于更新的&^^以"拉格朗日乘数的向量^^和^^可根据如下公式更 新:
[_ ]螓-p ((P4T)广 )
[0082] Ak+1 = Ak-a(Axk+1-b)
[0083] 4.4判断是否满足迭代停止条件:ADMM算法达到最大迭代数或x的结果满足公差< 1(T4。若不满足则重复开始执行4.1,满足则停止此子迭代算法,执行步骤S5。
[0084] S5.基于S4中对动态PET图像序列的估计值xk+1求解动力学参数的局部最优解 小k+1:提取出目标参数中仅包括变量巾的部分获得以下子优化问题:
[0086]令- (#_)T|^由于Q(x,小)是一个凸函数,因此应用迭代软阈 值算法对上述问题进行优化求解可得:
[0088] 其中:是Q(x,4〇在巾处的导数Aw是软阈值运算子,〇可视为一个步 长。
[0089] S6.判断是否满足迭代停止条件:总算法达到最大迭代数或x和巾的结果同时满足 公差<1(T4。若不满足则执行步骤S4;若满足迭代停止条件则迭代停止,保存PET图像重建结 果x和动力学参数巾,并根据以下公¥
:十算药物分布容积(volumeofdistrib ution,VD) 〇
[0090] 以下我们通过对蒙特卡洛模拟的脑部phantom试验来验证本实施方式重建框架结 果的准确性,图1为实验所用的脑部phantom的模板示意图,将不同的区域分为三个感兴趣 的区域(region of interest,R0I)。实验运行环境为:8G内存,3.40GHz,64位操作系统,CPU 为intel i7-3770。所模拟的PET扫描仪型号为Hamamatsu SHR-22000,设定的放射性核素及 药物为18F-FDG,设置sinogram为64个投影角度在每个角度下64条射束采集到的数据结果。 系统矩阵G为4096X4096维事先计算好的矩阵。为了验证重建系统的鲁棒性,共设置了5 X 104、1 X 105、5X 105、1 X 106四种不同计数率来获得测量数据sinogram。
[0091]本实施方式重建框架中对PET图像序列重建的结果和传统的ML-EM的重建结果做 比较。二者使用相同的测量数据矩阵y和系统矩阵G以控制结果的可比性。图2从左往后分别 是在5 X 104、1 X 105、5 X 105、1 X 106四种不同计数率下,本实施方式重建方法获得的第9帧 的动态PET图像序列的重建结果;图3是相同条件下通过ML-EM方法获得的结果。可以看出, 当计数率更低时,本发明能够使重建图像维持更多的特征,对比度更高,噪声更小;表1为其 进一步的数据量化分析结果。
[0092]弄 1
[0095]本实施方式重建框架中对动力学参数的估计结果用VD来展示,重建结果和一种基 于房室理论的通过数据来估计参数图像的方法(data-driven estimation of parametric images based on compartmental theory,DEPICT)进行比较。图4从左到右依次是真值,5 \104、1\105、5\105、1\106四种不同计数率下7说重建结果,可以看出,相对于真值,计数 率越高,重建的结果越接近真实的情况。即使是最低的计数率,也能一定程度上反映出真值 图像的特性,说明本实施方式重建框架对低计数率情况下的鲁棒性较好。另外表2是针对Vd 的重建数据结果,即为将本发明与DEPICT进行比较的相对均方根误差(relative mean square error,RMSE)结果。从数值上也可以看出本发明不仅能估计出一个相对准确的结 果,当计数率更低时,也能将结果的错误率稳定在一个更合理的范围内,即鲁棒性更好。
[0096]表 2
[0098] 上述的对实施例的描述是为便于本技术领域的普通技术人员能理解和应用本发 明。熟悉本领域技术的人员显然可以容易地对上述实施例做出各种修改,并把在此说明的 一般原理应用到其他实施例中而不必经过创造性的劳动。因此,本发明不限于上述实施例, 本领域技术人员根据本发明的揭示,对于本发明做出的改进和修改都应该在本发明的保护 范围之内。
【主权项】
1. 一种基于TV和稀疏约束的动态PET图像和示踪动力学参数同步重建方法,包括如下 步骤: (1) 利用探测器对注入有放射性药剂的生物组织进行探测,动态采集得到M组符合计数 向量,进而组建PET的符合计数矩阵Y,M为大于1的自然数; (2) 通过使动态PET图像序列组合成PET浓度分布矩阵X,根据PET成像原理,建立PET测 量方程; (3) 通过对PET测量方程引入全变分约束,得到基于TV的PET图像重建模型L(X)如下:其中:α为权重系数,G为系统矩阵,TV(X)为关于X的全变分正则项,I I I |2为L2范数; (4) 利用房室模型匹配估计示踪动力学参数,建立关于X和Φ的同步重建模型S(X,?) 如下:其中:μ为权重系数,ψ为字典矩阵,φ为示踪动力学参数矩阵,τ表示转置,I I I Ii为Li 范数,I I I 为F范数; (5) 将上述两个模型L(X)和S(X,Φ )相结合得到最终的同步重建模型LS(X,Φ )如下:其中:γ为权重系数; (6) 对所述的同步重建模型LS(X,Φ )进行最优化求解后即得到PET浓度分布矩阵X和示 踪动力学参数矩阵Φ。2. 根据权利要求1所述的动态PET图像和示踪动力学参数同步重建方法,其特征在于: 所述的符合计数矩阵Y即由M组符合计数向量组成,PET浓度分布矩阵X即由M组PET浓度分布 向量组成,所述的PET浓度分布向量即对应某一时刻的一帧PET图像浓度数据。3. 根据权利要求1所述的动态PET图像和示踪动力学参数同步重建方法,其特征在于: 所述PET测量方程的表达式如下: Y = GX+R+T 其中:R和T分别为关于反射符合事件和散射符合事件的测量噪声矩阵。4. 根据权利要求2所述的动态PET图像和示踪动力学参数同步重建方法,其特征在于: 所述全变分正则项TV(X)的表达式如下:其中:为PET浓度分布矩阵X中对应ψ的浓度值,ψ表示第m帧PET图像中的第j个 体素,(Dxm)j为关于$^的二维差分向量,该向量第一行元素值为If"与其右边体素的浓度 差值,第二行元素值为'T s与其下边体素的浓度差值,K为一帧PET图像的体素总个数。5. 根据权利要求1所述的动态PET图像和示踪动力学参数同步重建方法,其特征在于: 所述字典矩阵Ψ的表达式如下:其中=CMtWPC1(T)分别为t时刻和τ时刻放射性药剂在血浆中的浓度值,分别 为关于第m组符合计数向量采集的开始时间和结束时间,Θ。对应为第c个房室组织指数函数 的系数,m和c均为自然数且l^m<M,l<c<N,N为大于1的自然数;Θ :~ΘΝ的取值为在区间 [0min,0max]内按指数级间隔进行选取,0 min和0max分别为系数的上下限阈值,t和τ均表示时 间。6. 根据权利要求1所述的动态PET图像和示踪动力学参数同步重建方法,其特征在于: 所述的步骤(6)中采用ADMM结合软阈值的迭代优化算法对同步重建模型LS(X,Φ )进行最优 化求解;其中,ADMM针对PET浓度分布矩阵X进行迭代优化求解,软阈值针对示踪动力学参数 矩阵Φ迭代优化求解。
【文档编号】G06T11/00GK105894550SQ201610196767
【公开日】2016年8月24日
【申请日】2016年3月31日
【发明人】刘华锋, 余海青, 陈舒杭
【申请人】浙江大学
网友询问留言 已有0条留言
  • 还没有人留言评论。精彩留言会获得点赞!
1