CN108550172B - 一种基于非局部特性和全变分联合约束的pet图像重建方法 - Google Patents

Abstract

本发明公开了一种基于非局部特性和全变分联合约束的PET图像重建方法，该方法通过建立重建问题的数学模型，加入基于非局部特征的低秩约束以及加入全变分约束联合进行PET图像重建。其中低秩约束根据图像的非局部相似特性保证了图像的结构信息滤除噪声，同时避免了全变分方法存在阶梯效应这个缺点；而全变分作为基于像素点进行保边去噪的约束，也弥补了基于图像块存在的图像边缘棱角化的缺点，同时加快了收敛速度。故本发明有效利用约束联合求解，弥补了两方法的缺点，更进一步改善了计算机在进行PET图像重建的过程中产生的结果低分辨率和噪声干扰的问题；与现有重建方法的实验比较表明，本发明能获得较好的重建效果。

Description

一种基于非局部特性和全变分联合约束的PET图像重建方法

技术领域

本发明属于PET成像技术领域，具体涉及一种基于非局部特性和全变分联合约束的PET图像重建方法。

背景技术

正电子发射断层扫描(Positron EMission Tomography，PET)是当今最高层次的核医学影像技术，是目前惟一可在活体上显示生物分子代谢、受体及神经介质活动的新型影像技术。在临床领域，PET对具有不同代谢特性的组织和器官进行功能性成像，常被用于癌症和肿瘤的诊断；在神经科学领域，可以用于进行脑功能识别和脑功能分析的研究；在医药领域，对药理、药效具有独一无二的窥视效果，可用于辅助新药物的开发。另外，PET所用的示踪药物，是组成有机体组织的基本成分并能参与代谢过程。大多数疾病的生化变化先于解剖学的变化，并且PET对于示踪剂浓度的灵敏度非常高，能高精度地定量地检测出代谢过程的非正常增加，并给出清晰的图像所谓热源成像。因此能提供很多疾病在发展过程中的早期信息，可以进行超前诊断，尤其适合于肿瘤的早期诊断。总之，正电子发射断层成像己成为肿瘤、心、脑疾病诊断与病理研究中的不可缺少的重要方法。

PET显像就是从受到低计数率和噪声影响的投影数据中，重建出临床医生可接受的功能图像。事实上，PET图像的重建问题在理论上是一个病态问题。传统的解析法—滤波反投影法(Filtered Back-Projection，FBP)具有成像速度快的优点，但其重建图像含有大量噪声，图像质量较差。作为另一种重建技术，迭代重建算法对不完全数据的适应性好，能够得到更精确的结果，但由于计算量过大，以前只是被看成一种研究方法。随着计算机运算速度的飞速发展，迭代法已经被越来越广泛地应用于精确重建图像，其中包括著名的MLEM(最大似然期望最大化)、MAP(最大后验)和SAGE(惩罚似然)算法。考虑到临床上对快速成像日益迫切的需求，如何获得更精确、清晰的重建图像是目前研究的热点。

发明内容

鉴于上述，本发明提供了一种基于非局部特性和全变分联合约束的PET图像重建方法，将泊松重建模型、非局部低秩约束、全变分正则项耦合进入同一个目标方程，基于分布式优化的思想对于其中的三个子问题同时求解，最终更高效地获得更精确的PET重建图像。

一种基于非局部特性和全变分联合约束的PET图像重建方法，包括如下步骤：

(1)利用探测器对注入正电子标记示踪药物的生物组织进行探测，采集关于PET的符合事件，从而构建PET的符合计数向量y；

(2)通过对PET测量方程引入Poisson噪声，得到PET的Poisson模型P(x)；

(3)通过对Poisson模型P(x)引入非局部低秩约束和全变分约束，得到基于这两种约束的PET图像重建模型如下：

其中：NL(x,L)为关于PET浓度分布向量x和低秩矩阵L的非局部低秩约束项，TV(x,ω)为关于PET浓度分布向量x和二维差分向量ω的全变分约束项，α和β为权重系数；

(4)对上述PET图像重建模型进行最优化求解，得到PET浓度分布向量x，进而重建获得PET图像。

进一步地，所述步骤(2)中Poisson模型P(x)的表达式如下：

其中：

为对应符合计数向量y的均值向量，G为PET的系统矩阵，r和s对应为随机符合事件和散射符合事件的测量噪声向量，

为均值向量

中的第i个元素值，y_i为符合计数向量y中的第i个元素值，n为探测器中的晶体块个数即符合计数向量y的维度。

进一步地，所述步骤(3)中非局部低秩约束项NL(x,L)的表达式如下：

其中：B_k为块匹配操作矩阵，B_kx为由与PET图像中第k个图像块相似的图像块向量所组成的矩阵，L_k为B_kx对应的低秩矩阵，λ为权重系数，|| ||_F为Frobenius范数，|| ||_*为核范数，K为图像块的数目。

进一步地，所述步骤(3)中全变分约束项TV(x,ω)的表达式如下：

其中：ω_j为PET图像中第j个像素点的二维差分向量，D_j为对应第j个像素点的求差分算子，η为权重系数，|| ||₁为L1范数，|| ||₂为L2范数，m为PET浓度分布向量x的维度。

进一步地，所述步骤(4)中对PET图像重建模型进行最优化求解，即将问题分解成x、L、ω三个子问题进行分布式联合求解，其中x子问题采用EM(ExpectationMaximization，期望值最大)算法进行求解，L子问题采用核范数逼近低秩矩阵的方法进行求解，ω子问题采用软阈值的方法进行求解。

进一步地，对L子问题进行求解过程中关于低秩矩阵L_k则采用奇异值阈值方法求解，计算表达式为L_k＝U(Σ-λ)+V^T，其中U、Σ、V为对B_kx奇异值分解得到的矩阵即UΣV^T＝B_kx，T表示转置，(Σ-λ)₊表示对角矩阵Σ中每一个对角线元素值减去λ后与0进行比较，若大于等于0则保持不变，若小于0则用0替换为对角线元素值，即保证对角矩阵Σ中每一个对角线元素值减去λ后非负。

进一步地，对ω子问题采用软阈值的方法进行求解，即：

其中：v为拉格朗日乘子，D为求差分算子。

进一步地，对x子问题采用EM算法进行求解过程中引入隐藏变量c_qj以求解期望，即：

其中：隐藏变量c_ij表示从第j个体素中发射出被探测器中第i个晶体块探测到的光子数，g_ij为系统矩阵G中第i行第j列元素值，x_j为PET浓度分布向量x中的第j个元素值，r_i和s_i分别为测量噪声向量r和s中的第i个元素值。

本发明通过建立重建问题的数学模型，加入基于非局部特征的低秩约束以及加入全变分约束联合进行PET图像重建。其中低秩约束根据图像的非局部相似特性保证了图像的结构信息滤除噪声，同时避免了全变分方法存在阶梯效应这个缺点；而全变分作为基于像素点进行保边去噪的约束，也弥补了基于图像块存在的图像边缘棱角化的缺点，同时加快了收敛速度。故本发明有效利用约束联合求解，弥补了两方法的缺点，更进一步改善了计算机在进行PET图像重建的过程中产生的结果低分辨率和噪声干扰的问题；与现有重建方法的实验比较表明，本发明能获得较好的重建效果。

附图说明

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

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

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

图2(c)为采用TVAL方法重建脑部体模的PET图像。

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

具体实施方式

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

如图1所示，本发明基于非局部低秩和全变分联合约束的PET图像重建方法，包括如下步骤：

S1.将正电子示踪剂如F-18，注入被检查者体内并随后由探测器接收光子。

采集数据前，首先将正电子示踪剂(如F-18)注入被检查者体内，正电子示踪剂在体内发生衰变不断释放出正电子，正电子在组织内前进1-3mm，然后与体内的负电子发生湮灭反应，同时释放出两个方向互为180°的能量为511千电子伏特的光子，随后由两个探测器接收两个光子，确定发生湮灭的位置。

S2.正电子发射断层扫描仪通过探测湮灭产生的光子对，经过符合和采集系统处理，并以正弦图(sinogram)的方式存为待重建数据。

S3.根据PET探测的原理整理出相应的系统矩阵G，输入正弦图形式的投影矩阵y。

S4.初始化，设置迭代数k＝0，非局部低秩约束的权重α和全变分正则项的权重系数β，用来提取非局部特征的图像块边长p，步长s，全变分方法中的权重系数η；另外采用滤波反投影法x⁰＝FBP(y)进行热启动(warm start)，以加快迭代速度。

S5.求解以下子问题：

5.1 L子问题：根据第k次迭代的重建图像向量x^k，提取对于第i个图像块的矩阵

并根据奇异值阈值算法(SVT)计算低秩矩阵L_i ^k+1。

5.2 ω子问题：根据第k次迭代的重建图像向量x^k，采用软阈值算法更新二维差分向量ω^k+1，同时更新相应的拉格朗日乘子v^k+1←v^k-η(Dx^k-ω^k)。

5.3 x子问题：基于第k次迭代的重建图像向量x^k，按照EM算法求解第k+1次迭代的重建图像向量x^k+1。E步骤中，先根据x^k和正弦图y求解隐藏变量向量c^k+1；在M步骤中，通过求解一元二次方程求解出x^k+1。

S6.判断是否满足迭代停止条件

若不满足该条件则更新迭代次数k←k+1，并执行步骤S5；满足则迭代停止，输出重建图像向量x，实现PET图像重建。

本实验的运行环境为：Windows 10操作系统，8GB内存，Intel i7-6700处理器，3.40GHz。针对重建图像的验证，以Zubal Thorax Phantom的重建为例，我们采用128个投影角度，每个投影角度128条投影线的设置，即有一张重建图像的正弦图大小为128×128＝16384个像素。因此在步骤S3中我们将正弦图数据转化为向量y，维度为16384×1；生成系统矩阵G，为16384×16384的稀疏矩阵。

在步骤S4中，初始化参数，令迭代数k＝0；对于非局部低秩约束，设置权重参数α∈[2²,2⁵]，一般选取α＝8左右；设置图像特征块大小p＝4，提取特征块步长s＝3；对于全变分方法的参数，设置全变分正则项的权重系数β∈[2⁴,2⁶]，一般选取β＝30左右，权重系数η∈[2⁴,2⁶]，一般选取η＝40；滤波反投影预处理采用Fessler教授实验室公布的医学图像重建MATLAB工具包中的FBP函数。

在步骤S5进入迭代。

在求解L子问题时，先根据第k次迭代的重建图像向量x^k，提取对于第i个图像块的矩阵，有

其中

为第i个图像块向量，通过寻找m个图像块向量

使得欧氏距离

最近，构成低秩矩阵X^k+1 _i和操作矩阵

取m＝20，这里有X^k+1 _i∈C^16×20。

在低秩问题中采用核范数进行近似并采用奇异值阈值方法求解：

这里L_i采用奇异值阈值技术(SingularValueThresholding，SVT)进行求解：

L_i ^k+1＝U(∑-λ)₊V^T

其中：U∑V^T为

的奇异值分解，有(x)₊＝max(x,0)。

在求解ω子问题中，有：

基于增广拉格朗日求解方法的思想，得到拉格朗日问题如下：

使用软阈值(Soft-Thresholding)方法进行求解，有：

其中：v为拉格朗日乘子，接下来下来更新：v^k+1←v^k-η(Dx^k-ω^k+1)。

接下来求解x子问题，PET图像重建的模型基于如下方程：

y＝Gx+r+s

其中：重为系统矩阵，y为校正后的符合计数向量，x为PET浓度分布向量，r和s分别表示反射符合事件和散射符合事件的测量噪声矩阵。

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

其中：

表示y服从均值为

的泊松分布，基于独立泊松假设，y的似然方程Pr(y|x)表达式如下：

Poisson模型Pr(y|x)的似然方程取对数并求负后的表达式如下：

其中：

为符合计数向量y的均值向量，

为均值向量

中第i个元素值，i表示第i块探测器晶体块，而n_i表示探测器晶体块的总数，取对数并求负是便于最优化计算。

在EM算法的E步骤中引入隐藏变量c_qj：

其中：隐藏变量c_qj物理意义是从体素j中发射出被探测器晶体块q探测到的光子数，g_qj是系统矩阵G中的第qj项，q表示第q块探测器晶体块，j表示第j个体素。隐藏变量c_qj用期望值

替换：

EM方法的M步骤中求解x：

其中：g_qj是系统矩阵G中的第qj项，q表示第q块探测器晶体块，j表示第j个体素，i表示非局部约束中第i个图像块，x_j是当前估计出的图像x在第j个像素点的值，B_ij表示矩阵提取算子的第ij项，L_ij表示第i个低秩矩阵第j项，D_j表示在第j点的求差分算子，ω_j表示差分向量的第j项。

令Ω(x；x^k)对x_j的偏导为0：

观察得上式等效于关于x_j的一个一元二次方程，故这里直接使用一元二次求根公式求解x：

有：

其中：

解得：

在上述框架中对于三个子问题进行分布式联合最优化求解，得到的PET浓度分布矩阵x，从而获得PET图像。

针对重建图像质量的验证，一般采用高维度的原始采集数据128个投影角度，每个角度下射束为128条，即m＝16384，重建图像大小为128×128，即维度为16384；初始值设定同上。以脑部模拟图像重建为例，图2(a)～图2(d)分别为真值图像、传统ML-EM方法重建的图像、TVAL重建的图像和基于本实施方式重建图像的比对结果，可以直观地看出基于本实施方式重建的图像与ML-EM的结果相比，ML-EM算法重建图像中存在大量噪点且边缘十分不清晰，而TVAL算法重建的图像虽然抑制了噪声的干扰，但由于基于像素点求梯度的特性，区域边界出现阶梯状模糊。对比之下，本实施方式算法重建后的图片各区域有清晰连贯的边界，且图像内部平滑。

对于相同的数据，分别采用本实施方式、传统的ML-EM方法和基于交替方向乘子法的TVAL方法进行比较，如表1所示；应用本实施方式重建结果在与真值的偏差，方差均小于ML-EM和TVAL这两种方法，证明了本发明技术方案在提高精确度和降低噪声方面的可行性。

表1

方法	偏差	方差
			ML-EM	0.18953	0.05445
TVAL	0.15116	0.03165
			本实施方法	0.10027	0.01847

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

Claims (5)

1.一种基于非局部特性和全变分联合约束的PET图像重建方法，包括如下步骤：

(1)利用探测器对注入正电子标记示踪药物的生物组织进行探测，采集关于PET的符合事件，从而构建PET的符合计数向量y；

(2)通过对PET测量方程引入Poisson噪声，得到PET的Poisson模型P(x)表达式如下；

其中：

为对应符合计数向量y的均值向量，G为PET的系统矩阵，r和s对应为随机符合事件和散射符合事件的测量噪声向量，

为均值向量

中的第i个元素值，y_i为符合计数向量y中的第i个元素值，n为探测器中的晶体块个数即符合计数向量y的维度；

(3)通过对Poisson模型P(x)引入非局部低秩约束和全变分约束，得到基于这两种约束的PET图像重建模型如下：

其中：NL(x,L)为关于PET浓度分布向量x和低秩矩阵L的非局部低秩约束项，TV(x,ω)为关于PET浓度分布向量x和二维差分向量ω的全变分约束项，α和β为权重系数，B_k为块匹配操作矩阵，B_kx为由与PET图像中第k个图像块相似的图像块向量所组成的矩阵，L_k为B_kx对应的低秩矩阵，λ为权重系数，|| ||_F为Frobenius范数，|| ||_*为核范数，K为图像块的数目，ω_j为PET图像中第j个像素点的二维差分向量，D_j为对应第j个像素点的求差分算子，η为权重系数，|| ||₁为L1范数，|| ||₂为L2范数，m为PET浓度分布向量x的维度；

(4)对上述PET图像重建模型进行最优化求解，得到PET浓度分布向量x，进而重建获得PET图像。

2.根据权利要求1所述的PET图像重建方法，其特征在于：所述步骤(4)中对PET图像重建模型进行最优化求解，即将问题分解成x、L、ω三个子问题进行分布式联合求解，其中x子问题采用EM算法进行求解，L子问题采用核范数逼近低秩矩阵的方法进行求解，ω子问题采用软阈值的方法进行求解。

3.根据权利要求2所述的PET图像重建方法，其特征在于：对L子问题进行求解过程中关于低秩矩阵L_k则采用奇异值阈值方法求解，计算表达式为L_k＝U(Σ-λ)₊V^T，其中U、Σ、V为对B_kx奇异值分解得到的矩阵即UΣV^T＝B_kx，^T表示转置，(Σ-λ)₊表示对角矩阵Σ中每一个对角线元素值减去λ后与0进行比较，若大于等于0则保持不变，若小于0则用0替换为对角线元素值，即保证对角矩阵Σ中每一个对角线元素值减去λ后非负。

4.根据权利要求2所述的PET图像重建方法，其特征在于：对ω子问题采用软阈值的方法进行求解，即：

其中：v为拉格朗日乘子，D为求差分算子。

5.根据权利要求2所述的PET图像重建方法，其特征在于：对x子问题采用EM算法进行求解过程中引入隐藏变量c_qj以求解期望，即：

其中：隐藏变量c_ij表示从第j个体素中发射出被探测器中第i个晶体块探测到的光子数，g_ij为系统矩阵G中第i行第j列元素值，x_j为PET浓度分布向量x中的第j个元素值，r_i和s_i分别为测量噪声向量r和s中的第i个元素值。

Images