CN104851080B

CN104851080B - 一种基于tv的三维pet图像重建方法

Info

Publication number: CN104851080B
Application number: CN201510233428.5A
Authority: CN
Inventors: 刘华锋; 陈智
Original assignee: Zhejiang University ZJU
Current assignee: Zhejiang University ZJU
Priority date: 2015-05-08
Filing date: 2015-05-08
Publication date: 2017-11-17
Anticipated expiration: 2035-05-08
Also published as: CN104851080A

Abstract

本发明公开了一种基于TV的三维PET图像重建方法，包括：(1)利用PET的3D采集模式获得3D投影数据；(2)将3D投影数据重组成一系列2D投影数据；(3)从2D投影数据中计算出重建图像频域分布；(4)在直接傅里叶变换模型中加入全变差算子，建立目标函数；(5)采用BOSVS算法求解上述模型，得到2D重建结果；(6)将所有的2D重建结果合并成3D重建图像。本发明在直接傅里叶变换基础上引用全变差作为正则项建立模型，利用BOSVS算法进行迭代求解，获得了更高分辨率和边缘更加清晰的重建图像，在保证直接傅里叶变换快速重建的基础上显著地提高了重建图像的质量，同时在低计数率情况下仍然可以获得很好的结果。

Description

一种基于TV的三维PET图像重建方法

技术领域

本发明属于PET成像技术领域，具体涉及一种基于TV的三维PET图像重建方法。

背景技术

PET(Positron Emission Tomography)即正电子发射型计算机断层成像，作为核医学发展的一项先进技术，已成为目前癌症诊断，心脏功能检测以及脑功能成像非常重要的检查手段。其利用了核物理学和医学影像等多项新技术的优势，能够从分子水平上观察细胞代谢的活动。PET显像技术是在生物活体内注射一种能直接或间接反映生物新陈代谢过程的放射性同位素，通过PET设备中探测器环接收湮灭产生的光子对，进而计算确定正电子湮灭(发射)的位置，最后就可以反应出活体内的生理过程，以达到诊断和分析的目的。

从1976年第一台用于临床的商品化PET面世以来，PET系统不断发展和完善。20世纪80年代，在更多公司像西门子(Siemens)，通用电气公司(GE)等的投入研制下，很多PET新技术被开发了出来。在旧式PET仪器上，只能支持2D模式的数据采集，随着PET仪器的快速发展，3D采集数据越来越普遍。3D采集模式去掉了2D模式中探测器之间金属隔板的限制，考虑了所有探测器环内的符合计数，极大地提高了投影数据的计数率和探测器的灵敏度，但是3D采集的数据中含有更多的散射，进行重建等处理时需要大量的计算消耗和存储需求。

采样方式的不同决定了2D图像重建和3D图像重建的显著差异。在2D PET图像中，数据只采集同一平面或相临平面探测器的符合线(LOR)。这样每个平面都会重建出一个2D图像，所有的2D图像都是3D图像的一个截面，它们可以共同组成最后的重建结果。在3D成像中，符合计数是从所有平面的探测器的响应中产生，这将更好地利用产生的辐射，从而使探测器的灵敏度提高。在给定辐射剂量和成像时间下，3D成像检查到的事件一般是2D成像的5-10倍。随着灵敏度的提高，重建图像中的信噪比也会大大改善。但是，另一方面，3D测量数据将需要更大的数据存储和重建处理时间，同时图像中散射和随机因素将更多。这些不足在早期PET的发展过程中曾经阻碍过3D成像技术的发展，但是随着数据存储，计算速度及误差校正技术的不断更新，3D成像目前已广泛应用于临床医学中。三维PET因其采集的数据量大，重建速度较慢，如何保证重建图像精度的同时提高3D PET的重建速度成为目前的研究热门之一。

发明内容

本发明针对三维PET图像重建中存在的问题，建立了基于TV(total variation，全变差)的三维PET图像重建模型，并在此基础上使用可变迭代步长的伯格曼分裂算子(Bregman Operator Splitting with Variable Step size，简称BOSVS)方法对模型求解，同时从成像质量，参数鲁棒性，噪声特性等几个方面进行了实验方面的探索，初步形成了三维PET图像重建研究的一个新框架。

一种基于TV的三维PET图像重建方法，包括以下步骤：

(1)通过采集获取得到PET的3D投影数据；

(2)将所述的3D投影数据重组成一系列2D投影数据；

(3)根据所述的2D投影数据计算出对应二维PET重建图像u的空间频域分布F_s；

(4)在关于二维PET重建图像u和空间频域分布F_s的矩阵拟合模型中加入关于二维PET重建图像u的全变差算子TV(u)，并进行最优化求解，最终得到二维PET重建图像u；

(5)将所有二维PET重建图像u合并成三维PET图像。

进一步地，所述的步骤(2)中选用傅里叶重组算法(FORE)将3D投影数据重组成一系列2D投影数据；能够较大程度地提高了重建图像质量。

进一步地，所述的步骤(3)中根据中心切片定理对所述的2D投影数据做傅里叶变换后沿各投影角度摆放在空间域后进行双线性插值，从而得到对应二维PET重建图像u的空间频域分布F_s。

进一步地，所述的步骤(4)中在关于二维PET重建图像u和空间频域分布F_s的矩阵拟合模型中加入关于二维PET重建图像u的全变差算子TV(u)，则模型的最优化求解表达式如下：

其中：F()为傅里叶变换算符，μ为模型的调节参数，||||₂表示二范数。

进一步地，所述的步骤(4)中采用BOSVS算法对所述的最优化求解表达式进行求解。

进一步地，根据以下迭代方程对所述的最优化求解表达式进行求解：

λ^k+1＝λ^k-β(w^k+1-D(u^k+1)+1)

其中：D()为差分算符，λ^k+1和λ^k分别为第k+1次和第k次迭代的拉格朗日参数，u^k+1和u^k分别为第k+1次和第k次迭代的二维PET重建图像，w^k+1和w^k分别为第k+1次和第k次迭代的差分参数，Φ()为目标函数，β为目标函数的调节参数，k为迭代次数。

所述的目标函数Φ()的表达式如下：

其中：〈〉为内积算符，G(u)为关于二维PET重建图像u的矩阵拟合函数，w为差分参数，λ为拉格朗日参数。

所述的矩阵拟合函数G(u)的表达式如下：

所述二维PET重建图像u^k+1的迭代算式进一步表示为：

其中：δ为迭代固定步长，▽为梯度算子。

进一步地，利用Barzilai-Borwein步长δ_BB替代所述的迭代固定步长δ，由于Barzilai-Borwein步长δ_BB不再固定，故能够大大加速收敛速度；所述的Barzilai-Borwein步长δ_BB的表达式如下：

其中：为第k次迭代的Barzilai-Borwein步长。

本发明针对3D PET图像重建中存在的问题，建立了基于全变差的3D PET图像重建模型，并在此基础上使用可变迭代步长的伯格曼算子方法对模型求解，可获得更高分辨率和边缘更加清晰的重建图像；同时本发明从成像质量、参数鲁棒性、噪声特性等几个方面进行了实验方面的探索，初步形成了3D PET图像重建研究的一个新框架。故本发明方法在保证直接傅里叶变换快速重建的基础上显著地提高了重建图像的质量，同时在低计数率情况下仍然可以获得很好的结果。

附图说明

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

图2为zubal phantom数字体模的模型示意图。

图3(a)为关于zubal phantom蒙特卡罗仿真数据采用本发明方法重建后的中心切片图像。

图3(b)为关于zubal phantom蒙特卡罗仿真数据采用直接傅里叶重建后的中心切片图像。

具体实施方式

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

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

(1)利用PET的3D采集模式获得3D投影数据。

(2)将3D投影数据重组成一系列2D投影数据。选用傅里叶重组方法来将3D投影数据重组2D投影数据。傅立叶重组是在二维频率空间将倾斜响应线重组到距两探测器环中间平面轴向相位移为某一数值的平面上，是在频率空间依据频率-距离关系来进行的。傅里叶重组方法依然是一种近似方法，但是它较大程度地提高了重建图像质量，是目前现有3DPET中使用最广泛的一种重组算法。

(3)从2D投影数据中计算出重建图像频域分布。本实施方式利用中心切片定理和插值的方法来获取物函数的空间频域图；由中心切片定理，投影函数的一维傅里叶变换等于物函数的二维傅里叶变换沿投影方向的一个切片。根据此原理我们上述φ(x,y)做傅里叶变换后沿个投影角度摆放在空间域从而得到物函数的傅里叶频谱图F_s。F_s中，低频信号代表图像灰度的变化的快慢，而高频信号表示物体的空间分辨率以及图像的保边性。

(4)在直接傅里叶变换模型中加入全变差算子，建立目标函数。其数学模型如下式：

TV(u)＝∫_Ω|▽u(x)|dx

其中：μ是可调参数，||||₂是l₂范数，用来进行数据拟合。TV(u)是差分算子，用来进行平滑图像和保边。F为傅里叶变换算符，F_s是物函数的频域值，u即我们要求的重建图像。

(5)采用BOSVS算法求解上述模型，得到2D重建结果。

对于二维图像，TV可以写成如下表达式：

带入重建模型中得：

令w_i＝D_iu,i＝1,...,N并引入拉格朗日算子，重建模型转化为：

其中：

将变量分解，交替迭代步骤如下：

λ^k+1＝λ^k-β(w^k+1-Du^k+1)

对u的迭代可写为：

对G(u)进行Bregman算子分解(BOS)：

其中：δ为迭代固定步长。

对u的迭代可写为：

利用Barzilai-Borwein步长(δ_BB)代替固定步长δ；δ_BB的表达式如下：

因此，u可以由下式求解：

由于Barzilai-Borwein步长δ_BB不再固定，大大的加速了收敛速度。

(6)最后，将所有二维PET重建图像合并成三维PET图像。

以下通过实验来验证本实施方式的实用性和可靠性；图2为著名的zubalphantom蒙特卡罗仿真数据，最后所成图像像素为64×64；针对该数据，采用本实施方式重建的PET图像与采用直接傅里叶重建的PET图像相比较，从图3可以看出，用本实施方式重建出来的图像分辨率和图像清晰程度都比用直接傅里叶(DF)重建出来的要高得多。为了从统计学的角度分析本发明方法的优点，表1列出了在计数率分别为5×10⁵、1×10⁶、3×10⁶、6×10⁶、9×10⁶的情况下，用本实施方式与用直接傅里叶方法重建得到的图像的偏差和标准差。

表1

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

Claims

1.一种基于TV的三维PET图像重建方法，包括以下步骤：

(1)通过采集获取得到PET的3D投影数据；

(2)将所述的3D投影数据重组成一系列2D投影数据；

(3)根据中心切片定理对所述的2D投影数据做傅里叶变换后沿各投影角度摆放在空间域后进行双线性插值，从而得到对应二维PET重建图像u的空间频域分布F_s；

(4)在关于二维PET重建图像u和空间频域分布F_s的矩阵拟合模型中加入关于二维PET重建图像u的全变差算子TV(u)，则模型的最优化求解表达式如下，进而采用BOSVS算法对该表达式进行求解，最终得到二维PET重建图像u；

其中：F()为傅里叶变换算符，μ为模型的调节参数，|| ||₂表示二范数；

具体根据以下迭代方程对上述表达式进行求解：

<mrow> <msup> <mi>u</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msup> <mo>=</mo> <mi>arg</mi> <munder> <mrow> <mi>m</mi> <mi>i</mi> <mi>n</mi> </mrow> <mi>u</mi> </munder> <mi>&Phi;</mi> <mrow> <mo>(</mo> <mi>u</mi> <mo>,</mo> <msup> <mi>w</mi> <mi>k</mi> </msup> <mo>,</mo> <msup> <mi>&lambda;</mi> <mi>k</mi> </msup> <mo>)</mo> </mrow> </mrow>

<mrow> <msup> <mi>w</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msup> <mo>=</mo> <mi>arg</mi> <munder> <mi>min</mi> <mi>w</mi> </munder> <mrow> <mo>(</mo> <msup> <mi>u</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msup> <mo>,</mo> <mi>w</mi> <mo>,</mo> <msup> <mi>&lambda;</mi> <mi>k</mi> </msup> <mo>)</mo> </mrow> </mrow>

λ^k+1＝λ^k-β(w^k+1-D(u^k+1)+1)

其中：D()为差分算符，λ^k+1和λ^k分别为第k+1次和第k次迭代的拉格朗日参数，u^k+1和u^k分别为第k+1次和第k次迭代的二维PET重建图像，w^k+1和w^k分别为第k+1次和第k次迭代的差分参数，β为目标函数的调节参数，k为迭代次数，Φ()为目标函数且具体表达式如下：

<mrow> <mi>&Phi;</mi> <mrow> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>w</mi> <mo>,</mo> <mi>&lambda;</mi> <mo>)</mo> </mrow> <mo>=</mo> <mo>|</mo> <mo>|</mo> <mi>w</mi> <mo>|</mo> <mo>|</mo> <mo>+</mo> <mi>G</mi> <mrow> <mo>(</mo> <mi>u</mi> <mo>)</mo> </mrow> <mo>-</mo> <mo><</mo> <mi>&lambda;</mi> <mo>,</mo> <mi>w</mi> <mo>-</mo> <mi>D</mi> <mrow> <mo>(</mo> <mi>u</mi> <mo>)</mo> </mrow> <mo>></mo> <mo>+</mo> <mfrac> <mi>&beta;</mi> <mn>2</mn> </mfrac> <mo>|</mo> <mo>|</mo> <mi>w</mi> <mo>-</mo> <mi>D</mi> <mrow> <mo>(</mo> <mi>u</mi> <mo>)</mo> </mrow> <mo>|</mo> <msup> <mo>|</mo> <mn>2</mn> </msup> </mrow>

其中：<>为内积算符，G(u)为关于二维PET重建图像u的矩阵拟合函数，w为差分参数，λ为拉格朗日参数；

所述二维PET重建图像u^k+1的迭代算式进一步表示为：

<mrow> <msup> <mi>u</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msup> <mo>=</mo> <mi>arg</mi> <munder> <mrow> <mi>m</mi> <mi>i</mi> <mi>n</mi> </mrow> <mi>u</mi> </munder> <mo>{</mo> <mfrac> <mn>1</mn> <mrow> <mn>2</mn> <msubsup> <mi>&delta;</mi> <mrow> <mi>B</mi> <mi>B</mi> </mrow> <mi>k</mi> </msubsup> </mrow> </mfrac> <mo>|</mo> <mo>|</mo> <mi>u</mi> <mo>-</mo> <msup> <mi>u</mi> <mi>k</mi> </msup> <mo>+</mo> <msubsup> <mi>&delta;</mi> <mrow> <mi>B</mi> <mi>B</mi> </mrow> <mi>k</mi> </msubsup> <mo>&dtri;</mo> <mi>G</mi> <mrow> <mo>(</mo> <msup> <mi>u</mi> <mi>k</mi> </msup> <mo>)</mo> </mrow> <mo>|</mo> <msup> <mo>|</mo> <mn>2</mn> </msup> <mo>+</mo> <mfrac> <mi>&beta;</mi> <mn>2</mn> </mfrac> <mo>|</mo> <mo>|</mo> <mi>D</mi> <mrow> <mo>(</mo> <mi>u</mi> <mo>)</mo> </mrow> <mo>-</mo> <msup> <mi>w</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msup> <mo>-</mo> <mfrac> <msup> <mi>&lambda;</mi> <mi>k</mi> </msup> <mi>&beta;</mi> </mfrac> <mo>|</mo> <msup> <mo>|</mo> <mn>2</mn> </msup> <mo>}</mo> </mrow>

<mrow> <mfrac> <mn>1</mn> <msubsup> <mi>&delta;</mi> <mrow> <mi>B</mi> <mi>B</mi> </mrow> <mi>k</mi> </msubsup> </mfrac> <mo>=</mo> <msup> <mo>&dtri;</mo> <mn>2</mn> </msup> <mi>G</mi> <mrow> <mo>(</mo> <msup> <mi>u</mi> <mi>k</mi> </msup> <mo>)</mo> </mrow> </mrow>

其中：为第k次迭代的Barzilai-Borwein步长，为梯度算子；

(5)将所有二维PET重建图像u合并成三维PET图像。

2.根据权利要求1所述的三维PET图像重建方法，其特征在于：所述的步骤(2)中选用傅里叶重组算法将3D投影数据重组成一系列2D投影数据。