CN104077763A - 基于压缩感知理论的tof-pet图像重建方法 - Google Patents

Abstract

本发明提供一种基于压缩感知理论的TOF-PET图像重建方法。该方法包括：S1.根据飞行时间定位成像原理，通过TOF重建算法基于TOF-PET探测器得到的数据集以及一初始图像重建出第一图像；S2.根据压缩感知理论，对所述第一图像的一稀疏表达方式以l_p范数最小化为目标求解，并根据求解结果更新所述第一图像，得到第二图像；S3.判断是否达到停止条件：若是，则以所述第二图像为结果图像；若否，则将所述初始图像更新为所述第二图像，并转至步骤S1。通过本发明可以更充分的利用给出的时间信息，减少采集时间，降低药物剂量；同时，不但可以提高图像信噪比，而且对噪声有显著的抑制作用，因此可以获得更优的图像质量。

Description

基于压缩感知理论的TOF-PET图像重建方法

技术领域

本公开涉及核医学影像成像技术领域，尤其涉及一种基于压缩感知理论的TOF-PET图像重建方法。

背景技术

正电子发射断层成像(Positron Emission Tomography，PET)是临床上应用广泛的一种无创核医学影像诊断技术，其通过对注入活体的放射性示踪剂成像，从而提供活体的新陈代谢等功能信息，在临床诊断、疗效评价、基础医学研究以及新药研发中发挥着重要作用。其主要原理为：将标记有正电子核素的放射性示踪药物注射入被检测体内，正电子核素发生衰变发射出正电子，与被检测对象体内的负电子发生湮没反应，产生两个方向相反，能量均为511KeV的γ光子，从而被放置在被检测对象周围的闪烁晶体条探测到。经过电子学符合探测的处理，记录下符合要求的γ光子对，其中一次湮灭反应被称为一个事例。

如图1A以及图1B所示，探测到一个事例的两个伽马光子的晶体条之间的连线称作响应线12(Line of Response，LOR)。其中图1A展示的探测器11为环形探测器，图1B展示的探测器11为平板型探测器，更为具体的，平板型探测器可分为静态型与旋转型两种，但这些探测器成像原理均相同。记录大量这样的事例后，通过图像重建即可获得放射性示踪药物活度分布图。

随着现代PET技术的发展，PET的各部分性能一直在优化，一种基于飞行时间(Time-of-Flight，TOF)技术的PET(即TOF-PET)逐渐成为人们关注的焦点。TOF-PET与传统PET最大的区别在于它能够根据两个γ光子飞到两端晶体条的时间差，来大致确定湮没反应发生在响应线上的位置(即图1A以及图1B中的湮没点13)。由于传统PET无法预知湮没反应发生的位置，因而只能将每条响应线对应的事例等权重的分配到该响应线经过的所有路径，而TOF-PET则能够按照不同的权重(一般来说是高斯分布)对响应线上的计数来进行分配，例如可以如图2中所示。

理论上来说，只要时间信息足够精确，TOF-PET就可以根据γ光子入射到两端晶体条的时间差来完全确定湮没点13所在的位置。然而由于目前探测器的时间分辨率有限，时间测量具有一定的不确定性，仍然需要通过图像重建来获得放射性分布的图像。尽管如此，相比较传统PET而言，TOF-PET仍然具有巨大的优势：由于将反投影限定在了一定范围，TOF-PET可以显著提高图像信噪比以及对比度恢复系数，从而提高小病灶的检出率。同时，TOF-PET可以在保持图像质量的同时，降低药物剂量和减少采集时间，从而降低被探测者及操作者所受的辐射风险。

申请号为CN200780027583的中国专利申请公开了一种“用于改进的TOF-PET重建的方法和***”技术方案，该方案包括以下步骤：基于第一时间分辨率估计第一正电子湮没事件的位置；基于第二时间分辨率估计第二正电子湮没事件的位置；利用所估计的位置来重建正电子湮没事件的图像。

申请号为CN201080030171的中国专利申请公开了一种“利用基于飞行时间信息逐个事件生成的图像内容的飞行时间正电子发射断层摄影重建”技术方案，该方案根据每个事例的飞行时间信息独立定位湮没点以形成生成图像，并将生成图像适当地用作迭代重建的初始图像。

现有技术中，还有其他一些其他基于飞行时间技术的图像重建方法，在此不一一列举。在所有这些算法中，均已经证实相比较于传统PET来说，TOF-PET具有相当大的优势。

然而上述方法仍有改进空间。例如，现有技术中的TOF-PET虽然能显著提高图像信噪比，然而单纯对噪声的改善有限。

发明内容

针对现有技术中的部分或者全部问题，本公开提供一种基于压缩感知理论的TOF-PET图像重建方法，用于更进一步的提升TOF-PET重建图像的质量。

本公开的其他特性和优点将通过下面的详细描述变得显然，或部分地通过本公开的实践而习得。

根据本公开的一个方面，一种基于压缩感知理论的TOF-PET图像重建方法，包括：

S1.根据飞行时间定位成像原理，通过TOF重建算法基于TOF-PET探测器得到的数据集以及一初始图像重建出第一图像；

S2.根据压缩感知理论，对所述第一图像的一稀疏表达方式以l_p范数最小化为目标求解，并根据求解结果更新所述第一图像，得到第二图像；

S3.判断是否达到停止条件：若是，则以所述第二图像为结果图像；若否，则将所述初始图像更新为所述第二图像，并转至步骤S1。

在本公开的一种实例实施方式中，所述TOF-PET探测器得到的数据集的数据格式包括List Mode数据格式或者Sinogram数据格式。

在本公开的一种实例实施方式中，所述TOF重建算法为基于飞行时间的迭代算法；所述基于飞行时间的迭代算法的迭代次数为一次或者多次。。

在本公开的一种实例实施方式中，所述步骤S1中通过TOF重建算法由TOF-PET探测器得到的全部数据集或者全部数据集的子集重建出第一图像。

在本公开的一种实例实施方式中，在第一次执行所述步骤S1时，所述初始图像被设定为某一初始值，或者，所述初始图像由所述TOF-PET探测器得到的数据集进行解析重建获得。

在本公开的一种实例实施方式中，所述稀疏表达方式包括梯度变换、离散余弦变换、Fourier变换、离散小波变换以及基于冗余字典的稀疏变换等。

在本公开的一种实例实施方式中，所述步骤S2中，利用梯度下降法，对所述第一图像的一稀疏表达方式以l_p范数最小化为目标求解。

在本公开的一种实例实施方式中，对所述第一图像的一稀疏表达方式以l_p范数最小化为目标求解为一迭代过程，其迭代次数为一次或者多次。

在本公开的一种实例实施方式中，还包括

设置一松弛因子；

在所述第二图像与第一图像之间的差异程度大于预设大小时，根据所述松弛因子缩减所述第二图像与第一图像之间的差异。

在本公开的一种实例实施方式中，所述停止条件为：

所述步骤S1-S3的运行次数达到一预设次数；或者，

所述步骤S2中，更新前后，所述第一图像之间的差异程度小于一预设阈值。

本公开的示例实施方式所提供的基于压缩感知理论的TOF-PET图像重建方法中，基于压缩感知理论对图像进行还原，从而可以更充分的利用给出的时间信息，减少采集时间，降低药物剂量；同时，通过该方法，不但可以提高图像信噪比，而且对噪声有显著的抑制作用，因此可以获得更优的图像质量。

附图说明

通过参照附图详细描述其示例实施方式，本公开的上述和其它特征及优点将变得更加明显。

图1A为一环型探测器及一个事例的两个γ光子在晶体条间的响应线；

图1B为一平板型探测器及一个事例的两个γ光子在晶体条间的响应线；

图2为TOF-PET的原理示意图；

图3本公开示例实施方式中一种基于压缩感知理论的TOF-PET图像重建方法的流程示意图；

图4为飞行时间定位成像原理示意图；

图5A、5B为梯度变换前后的图像示意图；

图6为GATE软件模拟的TOF-PET探测器示意图；

图7、图8为模拟得到的Phantom对象示意图；

图9A为现有技术中TOF-PET方法的图像重建结果；

图9B为本公开示例实施方式中方法的图像重建结果；

图10为现有技术中TOF-PET方法的图像重建结果和本公开示例实施方式中方法的图像重建结果的对比噪声比对比示意图。

附图标记说明：

11：探测器

12：响应线

13：湮没点

S1-S3：步骤

具体实施方式

现在将参考附图更全面地描述示例实施方式。然而，示例实施方式能够以多种形式实施，且不应被理解为限于在此阐述的实施方式；相反，提供这些实施方式使得本公开将全面和完整，并将示例实施方式的构思全面地传达给本领域的技术人员。在图中相同的附图标记表示相同或类似的结构，因而将省略它们的详细描述。

此外，所描述的特征、结构或特性可以以任何合适的方式结合在一个或更多示例实施方式中。在下面的描述中，提供许多具体细节从而给出对本公开的示例实施方式的充分理解。然而，本领域技术人员将意识到，可以实践本公开的技术方案而没有所述特定细节中的一个或更多，或者可以采用其它的方法、材料、组元等。在其它情况下，不详细示出或描述公知结构或者操作以避免模糊本公开的各方面。

本示例实施方式中首先提供了一种基于压缩感知理论的TOF-PET图像重建方法。如图3中所示，该基于压缩感知理论的TOF-PET图像重建方法中主要包括：

步骤S1.根据飞行时间定位成像原理，通过TOF重建算法基于TOF-PET探测器得到的数据集以及一初始图像重建出第一图像；在第一次执行该步骤时，所述初始图像可以设定为某一初始值(例如，可以为0或1等等)，也可由TOF-PET探测器得到的数据集进行解析重建获得。

其中，飞行时间定位成像原理是指根据两个γ光子入射到两端晶体条的时间差来确定湮没点所在的位置的方法。如图4中粗线所示，其中湮没点13的位置具有一定的不确定范围。在该不确定范围内，湮没点13位置的概率分布应为***矩阵概率模型与标准概率模型的乘积。其中，***矩阵描述了探测器对被检测物的探测过程，它连接了图像空间和投影空间。一般来说，***矩阵反应了两方面的内容：一是像素与响应线之间的耦合定位，也即某一个像素发出的光子是否被某一条响应线探测到；二是像素与响应线之间的耦合程度，也即某一个像素发出的光子被某一条响应线探测到的概率。常用的***矩阵概率模型有点线模型、线积分模型、面积分模型、立体角模型以及点扩展函数***矩阵模型；常用的标准概率模型为高斯模型。

湮没反应的最可几位置可由下式得到：

$\mathrm{\Δx}=\frac{c}{2}\mathrm{\Δt}$

式中，Δx为湮没点偏离响应线中心的距离，Δt是两个γ光子的飞行时间差，c为光速。

所述TOF-PET探测器得到的数据集的数据格式包括List Mode数据格式或者Sinogram数据格式。其中，List Mode数据格式是将检测到的湮没反应事例的信息按照数据流的形式依次记录下来而形成的一种数据存储形式；Sinogram数据格式是将每条响应线上发生的事例进行合并存储而形成的一种数据存储形式。

上述的TOF重建算法为基于飞行时间的迭代算法。在基于飞行时间的迭代算法中，所述步骤S1中通过TOF重建算法由TOF-PET探测器得到的全部数据集或者全部数据集的子集重建出第一图像。此外，所述基于飞行时间的迭代算法可迭代一次，也可迭代多次。

以基于List Mode数据格式的最大似然期望最大(Maximum LikelihoodExpectation Maximization，MLEM)统计迭代重建算法为例，其公式为：

$f^{\left(k\right)}\left(i\right)=\frac{f^{(k-1)}\left(i\right)}{\underset{j\∈s_n}{\mathrm{\Σ}}{p}^{\mathrm{TOF}}(i,j)}\underset{j\∈s_n}{\mathrm{\Σ}}{{\mathrm{\η}}_j}^{\mathrm{att}}{p}^{\mathrm{TOF}}\frac{1/{\mathrm{\η}}^{\mathrm{multi}}}{\underset{i^{\′}}{\mathrm{\Σ}}{{\mathrm{\η}}_j}^{\mathrm{att}}{p}^{\mathrm{TOF}}(i^{\′},f)f^{(k-1)}\left(i^{\′}\right)+S+R}$

其中η^multi＝η^det_geomη^decayη^deadtime，分别表示几何探测效率因子、衰变因子和死时间因子，η^att为衰减校正因子，p^TOF为时间扩展函数约束的***矩阵，S为散射校正数据，R为偶然符合校正数据。

S2.根据压缩感知理论，对所述第一图像的一稀疏表达方式以l_p范数最小化为目标求解，并根据求解结果更新所述第一图像，得到第二图像。

其中，压缩感知理论是一种信号恢复理论，该理论指出假设信号是稀疏的或者可以稀疏表示的，则可以通过较少的采样值精确还原出原始信号。重建图像稀疏表示即原始信号在某个变换域里是稀疏或者近似稀疏的，原始信号通过稀疏变换后，绝大部分变换系数为0或者接近于0。上述的稀疏表达方式只要能满足上述重建图像稀疏表达即可，更为具体的，上述的稀疏表达方式可以为梯度变换、离散余弦变换、Fourier变换、离散小波变换以及基于冗余字典的稀疏变换等，也可以为多种稀疏表达方式的结合。

在医学图像中，将梯度变换作为图像的稀疏变换是非常有效的做法，梯度变换的公式为：

$\left|\stackrel{\→}{\▿}{f}_{s,t}\right|=\sqrt{{(f_{s,t}-f_{s-1,t})}^2+{(f_{s,t}-f_{s,t-1})}^2}$

以图5A中是以Shepp-Logan标准图像为例，对其做一个梯度变换，得到的图像如图5B所示。

l_p范数公式如下：

||x||_p＝(|x₁|^p+|x₂|^p+...+|x_n|^p)^1/p

以l₁范数为例，梯度变换的l₁范数可以近似为全变差(TV)变换：

${\left|\right|\stackrel{\→}{f}\left|\right|}_{\mathrm{TV}}={\mathrm{\Σ}}_{m,n}{({({\stackrel{\→}{f}}_{m+1,n}-{\stackrel{\→}{f}}_{m,n})}^2+{({\stackrel{\→}{f}}_{m,n+1}-{\stackrel{\→}{f}}_{m,n})}^2)}^{1/2}\≈{\left|\right|\▿\stackrel{\→}{f}\left|\right|}_{l_1}$

在该步骤中，可以利用梯度下降法或者其他相关算法，对所述第一图像的一稀疏表达方式以l_p范数最小化为目标求解。以梯度下降法为例，其可以描述为：

$x^{(k+1)}=x^{\left(k\right)}-{\mathrm{\α}}_k\▿\mathrm{\Ψ}\left(x^{\left(k\right)}\right)$

式中，Ψ(x^(k))为目标函数，α_k为搜索步长。应用在本示例实施方式中，即为目标函数，对其求偏导，有：

$\begin{array}{c}{v}_{s,t}=\frac{\∂{\left|\right|\stackrel{\→}{f}\left|\right|}_{\mathrm{TV}}}{{\∂f}_{s,t}}\≈\frac{(f_{s,t}-f_{s-1,t})+(f_{s,t}-f_{s,t-1})}{\sqrt{\mathrm{\ϵ}+{(f_{s,t}-f_{s-1,t})}^2+{(f_{s,t}-f_{s,t-1})}^2}}\\ -\frac{(f_{s+1,t}-f_{s,t})}{\sqrt{\mathrm{\ϵ}+{(f_{s+1,t}-f_{s,t})}^2+{(f_{s+1,t}-f_{s+1,t-1})}^2}}\\ -\frac{(f_{s,t+1}-f_{s,t})}{\sqrt{\mathrm{\ϵ}+{(f_{s,t+1}-f_{s,t})}^2+{(f_{s,t+1}-f_{s-1,t+1})}^2}}\end{array}$

式中，ε为一个极小的正数，搜索步长可有多种确定方式，甚至可以为1。搜索步长的一个示例实施方式为：

α_k＝||f⁽⁰⁾-f^(k)||₂

按照上述示例实施方式，此步骤可描述为：

${\stackrel{\→}{f}}^k={\stackrel{\→}{f}}^{k-1}-a{\mathrm{\α}}_k\frac{{\stackrel{\→}{v}}^{k-1}}{\left|{\stackrel{\→}{v}}^{k-1}\right|}$

其中：

α_k＝||f⁽⁰⁾-f^(k)||₂

${\stackrel{\→}{v}}^k={\frac{\∂{\left|\right|\stackrel{\→}{f}\left|\right|}_{\mathrm{TV}}}{\∂f}|}_{f={\stackrel{\→}{f}}^k}$

a的值可自行设定。由上可知，对所述第一图像的一稀疏表达方式以l_p范数最小化为目标求解为一迭代过程，其迭代次数为一次或者多次。

S3.判断是否达到停止条件：

若是，则以所述步骤S2中得到的第二图像为结果图像；

若否，则将所述初始图像更新为所述第二图像，并转至步骤S1，通过TOF重建算法基于TOF-PET探测器得到的数据集以及更新后的初始图像重建出第一图像，继而执行步骤S2、S3。

在步骤S1-S2的循环迭代过程中，可能出现第二图像与第一图像之间的差异程度过大，而引起非线性迭代过程的发散。因此，本示例实施方式中还包括：设置一松弛因子；在所述第二图像与第一图像之间的差异程度大于预设大小时，根据该松弛因子缩减所述第二图像与第一图像之间的差异，从而避免保证每一步之间循环迭代的结果不至于相差过大。

举例而言，所述停止条件可以是，所述步骤S1-S3的运行次数达到一预设次数，即在循环了预设次数后停止。所述停止条件也可以是，所述步骤S2中，更新前后，所述第一图像之间的差异程度小于一预设阈值，即更新前后，所得的图像并无明显区别，即可停止。

为了对上述示例实施方式进行结果验证，本示例实施方式中还用蒙特卡洛方法模拟了中科院高能物理研究所研发的人体PET***。在此，具体使用的模拟软件为GATE软件，GATE软件是一种专门用于模拟射线成像的蒙特卡洛模拟平台，具有较高的正确性和可靠性，已经成为研究机构和影像公司进行射线成像仪器开发必不可少的一种工具。

GATE软件模拟出来的TOF-PET探测器如图6所示，TOF-PET探测器及采集的具体参数如表1所示。

表1

探测器环内径(直径)	887mm
		module总数	64

block组成	4*1
		Block尺寸(mm2)	39.6*39.6
晶体组成	11*11
		晶体尺寸(mm3)	3.5*3.5*25
闪烁晶体材料	LYSO
		晶体缝隙材料	高反膜
缝隙大小(mm)	0.1
		晶体个数	30976
单环晶体个数	704
		横断面Block间接缝	3.57mm
轴向Block之间接缝	0.4mm
		能量分辨率	0.21
能窗	361keV～661keV
		时间分辨率	500ps
时间窗	6ns

仿照NEMA(National Electrical Manufactures Association，美国电气制造商协会)标准里的图像质量Phantom设计了一个Phantom对象(如图7所示)。此Phantom对象为圆柱形，半径为147mm，长度为214mm，材料填充为水，活度密度为0.14uCi/CC。其中置入六个大小不同的小球(如图8所示)，直径分别为37、28、22、17、13、10mm，其中最大的两个小球为冷源，其余小球为热源，活度密度为背底的8倍。小球用毛细管连接，毛细管和小球的管壁厚度为1mm。另外，还在圆柱中心***直径为50mm，长度为214mm的圆柱冷源，填充肺部材料。

取有效总计数为12M的数据，分为10个子集，分别用现有技术中TOF-PET方法和本发明的方法进行图像重建，总迭代次数为10次，重建图像分别如图9A和图9B所示。

由图9B相比9A图可以看出，本发明提出的方法在现有技术中TOF-PET方法的基础上进一步抑制了图像噪声，对于病灶尤其是小病灶的可识别度大幅度提高。此外，对这两幅图像计算对比噪声比(Contrast-to-noise ratio，CNR)，结果如图10中所示。由结果可知，本发明提出的算法重建出的图像与现有技术中TOF-PET方法重建出的图像相比，对比噪声比有了显著提高。由于对比噪声比是衡量病灶检出率的一个重要指标，因而本方法能够显著提高病灶(尤其是小病灶)的检出率。

本公开已由上述相关示例实施方式加以描述，然而上述示例实施方式仅为实施本公开的范例。必需指出的是，已揭露的示例实施方式并未限制本公开的范围。相反地，在不脱离本公开的精神和范围内所作的更动与润饰，均属本公开的专利保护范围。

Claims (10)

1.一种基于压缩感知理论的TOF-PET图像重建方法，其特征在于，包括：

S1.根据飞行时间定位成像原理，通过TOF重建算法基于TOF-PET探测器得到的数据集以及一初始图像重建出第一图像；

S2.根据压缩感知理论，对所述第一图像的一稀疏表达方式以l_p范数最小化为目标求解，并根据求解结果更新所述第一图像，得到第二图像；

S3.判断是否达到停止条件：若是，则以所述第二图像为结果图像；若否，则将所述初始图像更新为所述第二图像，并转至步骤S1。

2.根据权利要求1所述的方法，其特征在于，所述TOF-PET探测器得到的数据集的数据格式包括List Mode数据格式或者Sinogram数据格式。

3.根据权利要求2所述的方法，其特征在于，所述TOF重建算法为基于飞行时间的迭代算法；所述基于飞行时间的迭代算法的迭代次数为一次或者多次。

4.根据权利要求1所述的方法，其特征在于，在第一次执行所述步骤S1时，所述初始图像被设定为某一初始值，或者，所述初始图像由所述TOF-PET探测器得到的数据集进行解析重建获得。

5.根据权利要求1所述的方法，其特征在于，所述步骤S1中通过TOF重建算法由TOF-PET探测器得到的全部数据集或者全部数据集的子集重建出第一图像。

6.根据权利要求1所述的方法，其特征在于，所述稀疏表达方式包括梯度变换、离散余弦变换、Fourier变换、离散小波变换以及基于冗余字典的稀疏变换。

7.根据权利要求1所述的方法，其特征在于，所述步骤S2中，利用梯度下降法，对所述第一图像的一稀疏表达方式以l_p范数最小化为目标求解。

8.根据权利要求1所述的方法，其特征在于，对所述第一图像的一稀疏表达方式以l_p范数最小化为目标求解为一迭代过程，其迭代次数为一次或者多次。

9.根据权利要求1所述的方法，其特征在于，还包括：

设置一松弛因子；

在所述第二图像与第一图像之间的差异程度大于预设大小时，根据所述松弛因子缩减所述第二图像与第一图像之间的差异。

10.根据权利要求1所述的方法，其特征在于，所述停止条件为：

所述步骤S1-S3的运行次数达到一预设次数；或者，

所述步骤S2中，更新前后，所述第一图像之间的差异程度小于一预设阈值。