一种基于Poisson TV的动态PET图像重建方法

摘要

本发明公开了一种基于Poisson TV的动态PET图像重建方法，该方法通过建立重建问题的数学模型，等价问题模型的转化，并基于Low Rank方法重建PET图像；其中，结合Poisson和Low Rank模型对PET图像进行重建的过程中，等价模型问题在求解过程中采用交替最小化算法。故本发明有效利用Low-Rank算法，改善了计算机在进行动态PET图像重建的过程中产生的结果低分辨率和噪声干扰的问题；与现有重建方法的实验比较表明，本发明能获得较好的重建效果。

说明书

技术领域

本发明属于PET成像技术领域，具体涉及一种基于Poisson TV(泊松分布 的全变差)的动态PET图像重建方法。

背景技术

正电子发射断层成像(Positron emission tomography，PET)是一种基于核 物理学、分子生物学的医学影像技术，它能够从分子水平上观察细胞的代谢活 动，为早期疾病的检测和预防提供了有效依据。在进行PET扫描时首先需要将 由放射性同位核素标记的药物注入人体内，通过血液循环系统，这些物质在人 体内各组织器官中将形成一定的分布。由于放射性同位核素不稳定，将发生衰 变，衰变过程中所产生的正电子与组织中的电子发生湮灭反应，产生一对方向 几乎相反的伽玛光子，经由符合采集系统对这些带有放射性药物分布信息的成 对光子进行处理生成投影数据。通过相应的数学方法对投影数据进行反演求解， 可重建出放射性物质的空间浓度分布。

动态正电子发射断层扫描通过连续数据采集，获取多帧生理状态空间分布。 随后，图像通常被分割为多个ROI(感兴趣区域)代谢参数估计，并进一步病 理分析。

传统上，放射性浓度分布重建和感兴趣区域边界分割问题被视为两个顺序 步骤。放射性浓度分布重建往往采用统计迭代方法，由于迭代法基于统计学模 型，对不完全数据适应性好，逐渐成为PET重建算法研究关注点，其中包括著 名的MLEM(最大似然期望最大化)、MAP(最大后验)和SAGE(惩罚似然) 算法，然后再应用以聚类、形变模型或图割(Graph Cut)为基础的方法来分割 PET图像。然而，这些方法由于测量数据的复杂性，所得到的放射性浓度分布 存在低分辨率和噪声干扰的问题，造成后续分割的困难。

发明内容

本发明提供了一种基于Poisson TV的动态PET图像重建方法，能够解决计 算机在进行图像重建的过程中产生的结果低分辨率和噪声干扰的问题。

一种基于Poisson TV的动态PET图像重建方法，包括如下步骤：

(1)利用探测器对注入有放射性物质的生物组织进行探测，动态采集得到 PET的M组符合计数向量，进而构建PET的符合计数矩阵Y；M为大于1的自 然数；

(2)通过使PET图像序列组合成PET浓度分布矩阵X，根据PET成像原 理，建立PET的测量方程；所述的PET浓度分布矩阵X包含M组PET浓度分 布向量，每一组PET浓度分布向量对应一帧PET图像数据；

(3)通过对所述的测量方程引入Poisson噪声约束，得到PET的Poisson 模型Ψ(Y|X)；

(4)通过对Poisson模型Ψ(Y|X)引入Low Rank(低秩)，得到PET的Low  Rank模型如下：

$\underset{L,S}{\min }\{{\left|\right|L\left|\right|}_{*}+\lambda {\left|\right|S\left|\right|}_1+\mathrm{\mu \Psi }\left(Y\right|X)\}\mathrm{}s.t.\mathrm{}X=L+S$

其中：|| ||_*为核范数，|| ||₁为1-范数，L为PET图像序列的背景部分，S为从PET 图像序列分割出活动的生物组织部分，λ和μ均为权重系数；

(5)构建上述LowRank模型的增强型拉格朗日函数，并对其进行最小化 求解，同时获得背景部分L和生物组织部分S，进而根据L与S相加得到的PET 浓度分布矩阵X进行动态PET成像，从而获得M帧连续的PET图像。

所述的测量方程的表达式如下：

Y＝GX+R+T

其中：G为系统矩阵，R和T分别为关于反射符合事件和散射符合事件的测量噪 声矩阵。

所述的Poisson模型Ψ(Y|X)的表达式如下：

$\Psi \left(Y\right|X)=\underset{m=1}{\overset{M}{\Sigma }}\underset{i=1}{\overset{N}{\Sigma }}{\bar{y}}_{\mathrm{im}}-y_{\mathrm{im}}\log \left({\bar{y}}_{\mathrm{im}}\right)$

${\bar{y}}_m={\mathrm{Gx}}_m+r_m+t_m$

其中：G为系统矩阵，x_m为第m组PET浓度分布向量，y_im为第m组符合计数 向量y_m中第i个测量值，r_m和t_m分别为对应第m组关于反射符合事件和散射符 合事件的测量噪声向量，为对应符合计数向量y_m的均值向量，为均值向 量中第i个元素值，N为符合计数向量的维度。

所述的增强型拉格朗日函数的表达式如下：

$\left(\begin{array}{c}{L}_A(L,S,U,Q)=\\ {\left|\right|L\left|\right|}_{*}+\lambda {\left|\right|S\left|\right|}_1-<Z,X-(L+S)>+\frac{\beta }{2}{\left|\right|X-(L+S)\left|\right|}_F^2\\ +\mathrm{\mu \Psi }\left(Y\right|X)+v_L{H}_{\mathrm{vtv}}(▿U)-<Z_L,L-U>+\frac{{\beta }_L}{2}{\left|\right|L-U\left|\right|}_F^2\\ +v_S{H}_{\mathrm{vtv}}(▿Q)-<Z_S,S-Q>+\frac{{\beta }_S}{2}{\left|\right|S-Q\left|\right|}_F^2\end{array}\right)$

其中：L_A(L,S,U,Q)为关于L、S、U和Q的增强型拉格朗日函数，U和Q为辅助 变量矩阵，Z、Z_L、Z_S均为拉格朗日乘子，|| ||_F为F-范数，< >为内积算符，β、 β_L、β_S均为惩罚系数，▽为梯度算子，v_L和v_S均为拉格朗日系数，H_vtv( )为VTV (导数矩阵的最大奇异值的积分)函数。

所述的步骤(5)中通过以下迭代方程组对增强型拉格朗日函数进行最小化 求解：

$S^{k+1}=S_{\lambda /(\beta +{\beta }_S)}\left({\Omega }_S^k\right)$

${\Omega }_S^k=\frac{\beta (X^k-L^k-Z^k/\beta )+{\beta }_S(Q^k+Z_S^k/{\beta }_S)}{\beta +{\beta }_S}$

$L^{k+1}=D_{1/(\beta +{\beta }_L)}\left({\Omega }_L^k\right)$

${\Omega }_L^k=\frac{\beta (X^k-S^k-Z^k/\beta )+{\beta }_L(U^k+Z_L^k/{\beta }_L)}{\beta +{\beta }_L}$

其中：( )表示以λ/(β+β_S)为底的软收缩算子，( )表示以1/(β+β_L)为 底的奇异值阈值算子，L^k和L^k+1分别为第k次和第k+1次迭代的背景部分，S^k和 S^k+1分别为第k次和第k+1次迭代的生物组织部分，X^k为第k次迭代的PET浓 度分布矩阵，Z^k、和分别为对应Z、Z_L、Z_S第k次迭代的拉格朗日乘子，Q^k和U^k均为第k次迭代的辅助变量矩阵，k为迭代次数。

所述的拉格朗日乘子Z^k、Z_L^k和Z_S^k的表达式如下：

Z^k＝Z^k-1-β(X^k-L^k-S^k)

$Z_S^k=Z_S^{k-1}-{\beta }_S(S^k-Q^k)$

$Z_L^k=Z_L^{k-1}-{\beta }_L(L^k-U^k)$

$x_{\mathrm{jm}}^{k+1}=\frac{{-b}_{\mathrm{jm}}+\sqrt{{b_{\mathrm{jm}}}^2-{4a}_{\mathrm{jm}}{c}_{\mathrm{jm}}}}{{2a}_{\mathrm{jm}}}$

$a_{\mathrm{jm}}=\beta b_{\mathrm{jm}}=(\mu \underset{i=1}{\overset{N}{\Sigma }}{g}_{\mathrm{ij}}-\beta {[L^k+S^k+Z^k/\beta ]}_{\mathrm{jm}})$

$C_{\mathrm{jm}}=-\mu \underset{i=1}{\overset{N}{\Sigma }}{\hat{\omega }}_{\mathrm{ij}}^m\mathrm{}{\hat{\omega }}_{\mathrm{ij}}^m=\frac{g_{\mathrm{ij}}{y}_{\mathrm{im}}}{\underset{j=1}{\overset{K}{\Sigma }}{g}_{\mathrm{ij}}{x}_{\mathrm{jm}}^k+r_{\mathrm{im}}+t_{\mathrm{im}}}{x}_{\mathrm{jm}}^k$

其中：Z^k-1、和分别为对应Z、Z_L、Z_S第k-1次迭代的拉格朗日乘子，g_ij为系统矩阵G的第i行第j列元素值，[ ]_jm表示[ ]内矩阵的第j行第m列元素值， y_im为第m组符合计数向量y_m中第i个测量值，为第k次迭代的PET浓度分 布矩阵内第m组PET浓度分布向量中第j个元素值，r_im和t_im分别为对应第m 组关于反射符合事件和散射符合事件的测量噪声向量中第i个元素值，K为PET 浓度分布向量的维度，N为符合计数向量的维度。

所述的辅助变量矩阵Q^k和U^k的表达式如下：

$U^k=F^{-1}\left\{\frac{F\{{\beta }_L{H}^k-{▿}^T{Z}_{\mathrm{vtv}}^k+{\beta }_{\mathrm{vtv}}{▿}^T{E}^k\}}{{\beta }_{L}I+{\beta }_{\mathrm{vtv}}F\{{▿}^T▿\}}\right\}$

$Q^k=F^{-1}\left\{\frac{F\{{\beta }_S{H}^k-{▿}^T{Z}_{\mathrm{vtv}}^k+{\beta }_{\mathrm{vtv}}{▿}^T{E}^k\}}{{\beta }_{S}I+{\beta }_{\mathrm{vtv}}F\{{▿}^T▿\}}\right\}$

$\left(\begin{array}{cc}{H}^k=L^k-Z_U^k/{\beta }_L& Z_U^k=Z_U^{k-1}-{\beta }_U(U^k-L^k)\end{array}\right)$

$Z_{\mathrm{vtv}}^k=Z_{\mathrm{vtv}}^{k-1}-{\beta }_{\mathrm{vtv}}(E^k-▿U^k)$

$E^k=\max \{1-\frac{v_L}{{\beta }_{\mathrm{vtv}}}/\sqrt{\underset{m}{\overset{M}{\Sigma }}{\left|\right|w_{\mathrm{jm}}^k\left|\right|}_2^2},0\}·W^k$

$W^k=▿U^k+Z_{\mathrm{vtv}}^{k-1}/{\beta }_{\mathrm{vtv}}$

其中：Z_U和Z_vtv均为拉格朗日乘子，和分别为对应Z_U第k次迭代和第k-1 次迭代的拉格朗日乘子，和分别为对应Z_vtv第k次迭代和第k-1次迭代的 拉格朗日乘子，β_U和β_vtv均为惩罚系数，F{ }和F^-1{ }分别表示傅里叶变换和傅 里叶反变换，为矩阵W^k的第j行第m列元素值，|| ||₂为2-范数，I为单位矩阵。

本发明通过建立重建问题的数学模型，等价问题模型的转化，并基于Low Rank方法重建PET图像；其中，结合Poisson和Low Rank模型对PET图像进 行重建的过程中，等价模型问题在求解过程中采用交替最小化算法。故本发明 有效利用Low-Rank算法，改善了计算机在进行动态PET图像重建的过程中产 生的结果低分辨率和噪声干扰的问题；与现有重建方法的实验比较表明，本发 明能获得较好的重建效果。

附图说明

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

图2(a)为关于肺部体模的真值图像。

图2(b)为采用EM算法重建肺部体模的PET图像。

图2(c)为采用本发明方法重建肺部体模的PET图像。

图2(d)为采用本发明方法分割出的肺部组织的PET图像。

具体实施方式

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

如图1所示，本发明基于PoissonTV的动态PET图像重建方法，包括如下 步骤：

正电子发射断层扫描仪探测人体内发出的放射性信号，经过符合和采集系 统处理，形成原始投影线，并以正弦图的方式存放于计算机硬盘中；对原始采 集到的sinogram，以sinogram和已知的系统矩阵G为输入项，调用相关模块。

S1.根据PET探测的原理建立重建问题的基本模型；

S2.引入Poisson和Low-Rank来优化问题模型；

S3.初始化，设置初始点U⁰，Q⁰，L⁰和S⁰；设置Z＝Z_L＝Z_S＝0。

S4.设置的初始点开始，按照交替最小化方法计算增强版本的拉格朗日方程 中的噪声r^k+1和重建的图像的Low-Rank部分L^k+1以及Sparse部分S^k+1；

S5.更新Z、Z_L、Z_S；准备将更新后的参数返回步骤S4进行循环；

S6.判断是否满足迭代停止条件，不满足该条件则执行步骤S5，满足则迭 代停止；进而使L+S＝X得到的PET浓度分布向量X实现PET成像。

为了完成PET图像的重建，PET检测过程的基本模型基于如下方程：

Y＝GX+R+S

其中：G为系统矩阵，Y为校正后的符合计数向量，X为PET浓度分布向量，R 和S分别表示反射符合事件和散射符合事件的测量噪声矩阵。

符合计数向量的Poisson模型的表达式如下：

$y~\mathrm{Poisson}\left\{\bar{y}\right\},\bar{y}=\mathrm{Gx}+r+s$

其中：表示y服从均值为的泊松分布，r和s分别表示反射符合事件 和散射符合事件的测量噪声向量，基于独立泊松假设，y的似然方程Pr(y|x)表达 式如下：

$\mathrm{Pr}\left(y\right|x)=\underset{i}{\overset{n_i}{\Pi }}{e}^{-{\bar{y}}_i}\frac{{\bar{y}}_i^{\mathrm{yi}}}{y_i!},\bar{y}=\mathrm{Gx}+r+s$

$\underset{x}{\min }\Psi \left(y\right|x)=-\log \left(\mathrm{Pr}\left(y\right|x)\right)=\underset{i}{\overset{n_i}{\Sigma }}{\bar{y}}_i-y_i\log \left({\bar{y}}_i\right)$

$\bar{y}=\mathrm{Gx}+r+s$

其中：Ψ(y|x)为重建问题的目标函数，y_i为第i个探测器的测量值。进一步的， 可以将上式扩展到图像序列的情况，如下：

$\underset{x}{\min }\Psi \left(Y\right|X)=\underset{m=1}{\overset{n_m}{\Sigma }}\underset{i}{\overset{n_i}{\Sigma }}{\bar{y}}_{\mathrm{im}}-y_{\mathrm{im}}\log \left({\bar{y}}_{\mathrm{im}}\right)$

$\bar{Y}=\mathrm{GX}+R+S$

Y＝[y₁,y₂,…,y_m,…y_nm]

X＝[x₁,x₂,…,x_m,…x_nm]

其中：y_im是第m祯i个探测器测量的数值，x_m是第m祯PET图像。

其Low Rank模型的表达式如下：

$\underset{L,S}{\min }{\left|\right|L\left|\right|}_{*}+\lambda {\left|\right|S\left|\right|}_1+\mathrm{\mu \Psi }\left(Y\right|X)s.t.\mathrm{}X=L+S$

其中：||L||_*为L的核范数，||S||₁为S的L-1范数Y为校正后的符合计数向量，L 为PET图像中的背景，S为PET图像中变化的组织，λ和μ均为权重系数。它对 应的增强型拉格朗日函数的表达式如下：

$\left(\begin{array}{c}L(L,S,U,Q)=\\ {\left|\right|L\left|\right|}_{*}+\lambda {\left|\right|S\left|\right|}_1-<Z,X-(L+S)>+\frac{\beta }{2}{\left|\right|X-(L+S)\left|\right|}_F^2\\ +\mathrm{\mu \Psi }\left(Y\right|X)+v_L{H}_{\mathrm{vtv}}\left(\overline{▿}U\right)-<Z_L,L-U>+\frac{{\beta }_L}{2}{\left|\right|L-U\left|\right|}_F^2\\ +v_S{H}_{\mathrm{vtv}}\left(\overline{▿}Q\right)-<Z_S,S-Q>+\frac{{\beta }_S}{2}{\left|\right|S-Q\left|\right|}_F^2\end{array}\right)$

其中：L(L,S,U,Q)为关于L、S、U和Q的增强型拉格朗日函数，Z，Z_L,Z_S为拉 格朗日乘子；G为系统矩阵，Y为校正后的符合计数向量，L为PET图像中的 背景，S为PET图像中变化的组织，U和Q为辅助变量，|| ||_*为核范数，|| ||₁为 1-范数，|| ||_F为F-范数，为变分算子，<·>为内积算符，λ和μ均为权重系数，β 为惩罚系数；R_vtv( )表示求VTV即R_vtv(u)＝∫_Ωσ₁(Du)dx。

进一步的，通过以下迭代方程组对增强型拉格朗日函数进行最小化求解：

$S^{k+1}=S_{\lambda /(\beta +{\beta }_S)}\left({\Omega }_S\right)$

$s.t.\mathrm{}{\Omega }_S=\frac{\beta (X-L-Z/\beta )+{\beta }_S(v+Z_S/{\beta }_S)}{\beta +{\beta }_S}$

$L^{k+1}={\mathrm{DD}}_{1/(\beta +{\beta }_L)}\left({\Omega }_L\right)V^T$

$s.t.\mathrm{}{\Omega }_L=\frac{\beta (X-L-Z/\beta )+{\beta }_L(u+Z_L/{\beta }_L)}{\beta +{\beta }_L}$

$x_{\mathrm{jm}}^{k+1}=\frac{{-b}_{\mathrm{jm}}+\sqrt{{b_{\mathrm{jm}}}^2-{4a}_{\mathrm{jm}}{c}_{\mathrm{jm}}}}{{2a}_{\mathrm{jm}}}$

$s.t.\mathrm{}{a}_{\mathrm{jm}}=\beta ,b_{\mathrm{jm}}=(\mu \underset{i=1}{\overset{n_i}{\Sigma }}{g}_{\mathrm{ij}}-\beta {[L+S+Z/\beta ]}_{\mathrm{jm}}),$

$C_{\mathrm{jm}}=-\mu \underset{i=1}{\overset{n_i}{\Sigma }}{\hat{\omega }}_{\mathrm{ij}}^m$

${\hat{\omega }}_{\mathrm{ij}}^m=\frac{g_{\mathrm{ij}}{y}_{\mathrm{im}}}{{\Sigma }_{l=1}^{\mathrm{nj}}{g}_{\mathrm{il}}{x}_{\mathrm{lm}}^{\left(k\right)}+r_{\mathrm{im}}+s_{\mathrm{im}}}{x}_{\mathrm{jm}}^k$

其中：L^k和L^k+1分别为第k次和第k+1次迭代PET图像中的背景，S^k和S^k+1分别 为第k次和第k+1次迭代PET图像中变化的组织，Z^k和Z^k+1分别为第k次和第 k+1次迭代的拉格朗日乘子，表示奇异值阈值算子，S_τ/β( )表示软收缩算 子，G为系统矩阵，Y为校正后的符合计数向量，μ为权重系数，β为惩罚系数， 为第k+1次迭代中第m祯第j个像素点的数值。

通过上述迭代方程组进行迭代计算，则使达到最大迭代次数或迭代收敛后 PET图像中的背景和变化的组织相加即作为估计得到的PET浓度分布向量；迭 代收敛条件如下：

$\frac{{\left|\right|L^{k+1}-L^k\left|\right|}_F}{{\left|\right|L^k\left|\right|}_F}<{10}^{-3}$且$\frac{{\left|\right|S^{k+1}-S^k\left|\right|}_F}{{\left|\right|S^k\left|\right|}_F}<{10}^{-3}$

以下我们采用脑部数字体模模型实验来验证本实施方式的有效性，该模型 包含一些高浓度区域。本实验运行环境为：8G内存，3.40GHz，64位操作系统， CPU为intel酷睿双核。

将本实施方式基于Poisson TV的动态PET图像重建方法和传统的EM(最 大期望算法)方法重建结果做比较，二者使用相同的观测值Y和相同的系统矩 阵D以保证结果的可比性，具体参数设置如下：Y为n×n维采集到的sinogram， 令m＝n×n，D为m×m维事先计算好的系统矩阵；这里n＝64，即m＝4096。

针对重建图像质量的验证，采用高维度的原始采集数据64个投影角度，每 个角度下射束为64条，即m＝4096，重建图像大小为64×64，即维度n＝4096； 初始值设定同上。图2(a)～(c)是真值图像、传统EM方法重建的图像和基于本实 施方式重建的图像的比对示意图，可以直观地看出基于本实施方式重建的图像 与EM的结果相比，能够恢复出更多的结构，且能够在重建的同时完成图像分 割，获得活动组织的图像。图2(d)是分割出来的活动的生物组织图像。

对于相同的数据，分别采用本实施方式和传统的EM方法进行比较，如表1 所示；应用本实施方式重建结果在与真值的偏差，方差和均方根误差均小于传 统的EM方法，说明对本发明技术方案在提高精确度和降低噪声方面的可行性。

表1

方法 偏差 方差 均方根误差 EM 0.0602 0.0131 0.1130 本实施方法 0.0312 0.0052 0.0701