CN105844678A - 基于全广义变分正则化的低剂量x射线ct图像重建方法 - Google Patents

Abstract

一种基于全广义变分正则化的低剂量X射线CT图像重建方法，包括如下步骤：(1)获取CT设备的系统参数和低剂量扫描协议下对数变换后的投影数据；(2)对投影数据进行Anscombe变换，将服从Poission分布的投影数据转化为近似服从方差为1的Gaussian分布数据u；(3)建立基于全广义变分最小化的理想数据恢复模型，使用Chambolle‑Pock算法求解得到恢复后的投影数据f；(4)对步骤(3)得到的恢复后的投影数据f进行Anscombe逆变换，再通过滤波反投影算法得到CT重建图像。该方法可以在投影数据不满足分段常数假设的前提下去除图像的噪声和条形伪影，同时较好地保持图像的空间分辨率。

Description

基于全广义变分正则化的低剂量X射线CT图像重建方法

技术领域

本发明涉及医学影像成像技术领域，特别是涉及一种基于全广义变分正则化的低剂量X射线CT图像重建方法。

背景技术

X射线计算机断层成像(Computed Tomography，CT)因其在时间、空间分辨率上的卓越表现，已经广泛应用于临床诊断和治疗。CT图像的质量与X射线辐射剂量密切相关，剂量越高图像质量越好，然而，过量的X射线照射又会诱发恶性肿瘤、白血病以及其他遗传性疾病。

为了降低X射线的辐射剂量，许多优化的低剂量CT扫描方案以及抑制噪声和伪影的重建算法相继被提出。目前，直接降低管电流、管电压或者减少扫描角向采样数目是实现低剂量CT扫描最简单且最有效的方法。对于这一问题有两种传统的解决方法:一是直接对低剂量CT图像进行滤波,以减少图像的噪声和伪影，属于图像后处理技术；二是根据投影数据满足的统计学规律,完成基于统计的CT图像迭代重建。

第一类方法为后处理技术，直接减少图像的噪声和伪影。图像后处理技术因其简单且易于操作。由于低剂量CT图像中噪声和伪影分布的复杂性，使得高精度的滤波器设计困难极大。

第二种方法是通过CT系统建模，构建图像重建模型，通过优化目标函数实现图像重建。相对于经典的FBP算法，迭代重建算法通过系统建模(系统光学模型和系统统计模型)对CT成像几何、X射线的能谱特性、射束硬化效应、散射和噪声特性进行准确描述，而且易于加入先验信息约束，因此特别适合低剂量CT图像优质重建。统计迭代重建在抑制图像噪声和伪影以及提高空间分辨率等方面都有上佳表现，还可以结合其他硬件相关的低剂量成像技术，进一步降低辐射剂量，提高重建图像的质量。然而，由于统计迭代重建需要反复进行投影与反投影运算，且CT图像数据量庞大，导致CT图像重建速度特别慢，难以满足临床中实时交互的需求。

不同于上述两种方法的另一种策略是，基于低剂量CT投影数据恢复的图像重建方法对CT成像系统进行精确建模(包括X线源、探测器、电子噪声、成像物体等)，可以实现优质的低剂量CT图像重建。其目标函数是根据测量数据的噪声特性来构建的，通常包含两项，即数据保真项和正则化项，前者用于描述投影数据的统计特性，后者用于修正解的特性。然而，许多临床CT图像并不是完全满足分段常量的假设，在剂量特别低时或者投影角度特别少时，TV重建的图像会产生阶梯效应和块状伪影。因此TV正则化在基于投影数据恢复的低剂量CT重建中会失去应有的效力。

因此，针对现有技术不足，有必要提供一种基于全广义变分正则化的低剂量X射线CT图像重建方法，其可以在投影数据不满足分段常数假设的前提下去除图像的噪声和条形伪影，同时较好地保持图像的空间分辨率。

发明内容

本发明的目的在于避免现有技术的不足之处而提供一种基于全广义变分正则化的低剂量X射线CT图像重建方法，该基于全广义变分正则化的低剂量X射线CT图像重建方法，可以在投影数据不满足分段常数假设的前提下去除图像的噪声和条形伪影，同时较好地保持图像的空间分辨率。

本发明的上述目的通过如下技术手段实现。

提供一种基于全广义变分正则化的低剂量X射线CT图像重建方法，包括如下步骤：

(1)获取CT设备的系统参数和低剂量扫描协议下对数变换后的投影数据q；

(2)对步骤(1)获取的投影数据q进行Anscombe变换，将服从复合Poission分布的投影数据q转化为近似服从方差为1的Gaussian分布数据u；

(3)对步骤(2)获取的数据u建立基于全广义变分最小化的理想数据恢复模型，使用Chambolle-Pock算法求解得到恢复后的投影数据f；

(4)对步骤(3)得到的恢复后的投影数据f进行Anscombe逆变换，再通过滤波反投影算法得到CT重建图像。

上述步骤(1)中获取的CT设备的系统参数包括X射线入射光子强度I₀。

上述步骤(2)中对步骤(1)获取的投影数据q进行Anscombe变换，计算式如下：

u＝(u₁,u₂,…,u_N)^T表示投影数据经过Anscombe变化后近似服从方差为1的Gaussian分布数据，其中，T表示转置运算，u₁、u₂、…、u_N是u的元素，N为元素个数；令f＝(f₁,f₂,…,f_N)^T为待估计投影数据的理想值，f₁、f₂、…、f_N是f的分量。

上述步骤(3)中对步骤(2)获取的数据u建立全广义变分最小化模型，具体通过如下方法进行：

(3-1)投影数据恢复的一般模型通过式(2)表示：

其中，第一项为数据保真项，第二项为正则化项，λ＞0是正则化参数；

(3-2)通过全广义变分正则化代替式(2)中的正则化项R(f)，得到基于全广义变分最小化的理想数据恢复模型为：

TGV(f)是f的全广义变分，

$TGV\left(f\right)=\underset{\mathrm{\ω}}{min}{\mathrm{\α}}_1{\∫}_{\mathrm{\Ω}}\left|\▿f-\mathrm{\ω}\right|dx+{\mathrm{\α}}_0{\∫}_{\mathrm{\Ω}}\left|\mathrm{\ϵ}\left(\mathrm{\ω}\right)\right|dx;$

其中α₀,α₁>0是两个正数，是f的梯度，ω是对偶变量， 是ω的梯度。

进一步的，式(3)的离散形式为：

其中，F＝R^N×N,W＝R^2N×N以及微分算子ε，由限差分算子近似得到，R^N×N表示N×N维的实空间；

根据对偶原理，将式(4)转化为式(5)的鞍点问题：

其中，p和q为对偶变量，并且P和Q定义为：

P＝{p∈R^2N×N|||p||_∞≤α₁}……式(6)；

Q＝{q∈R^3N×N|||q||_∞≤α₀}……式(7)；

定义和如下：

是邻近算子，其由下式表示：

其中，τ是步长，τ>0。

本发明的基于全广义变分正则化的低剂量X射线CT图像重建方法，包括如下步骤：(1)获取CT设备的系统参数和低剂量扫描协议下对数变换后的投影数据q；(2)对步骤(1)获取的投影数据q进行Anscombe变换，将服从复合Poission分布的投影数据q转化为近似服从方差为1的Gaussian分布数据u；(3)对步骤(2)获取的数据u建立基于全广义变分最小化的理想数据恢复模型，使用Chambolle-Pock算法求解得到恢复后的投影数据f；(4)对步骤(3)得到的恢复后的投影数据f进行Anscombe逆变换，再通过滤波反投影算法得到CT重建图像。该方法可以在投影数据不满足分段常数假设的前提下去除图像的噪声和条形伪影，同时较好地保持图像的空间分辨率。

附图说明

利用附图对本发明作进一步的说明，但附图中的内容不构成对本发明的任何限制。

图1是本发明一种基于全广义变分正则化的低剂量X射线CT图像重建方法的流程图。

图2是本发明实施例2的数值Clock体模及采用不同方法重建的结果。其中，(a)是真实的Clock体模图像；(b)表示采用斜波滤波反投影方法(FBP算法，Ramp窗)重建的图像；(c)表示采用汉宁滤波反投影方法(FBP算法，Hanning窗)重建的图像；(d)表示采用非单调最小化方法(NTVM方法)重建的图像；(e)表示采用本发明的方法重建的图像。

图3是对应图2数值Clock体模重建结果的局部放大图；每行分别对应一个区域，每行的第一列分别是真实的Clock体模图像，第二列是采用斜波滤波反投影方法(FBP算法，Ramp窗)重建的图像；第三列是采用汉宁滤波反投影方法(FBP算法，Hanning窗)重建的图像；第四列是采用非单调最小化方法(NTVM方法)重建的图像；第五列是采用本发明的方法重建的图像。

图4是对应图2中不同算法重建图像的线剖面图。

图5是本发明实施例2的数值Shepp-Logan体模及采用不同方法重建的结果。其中，(a)是真实的Shepp-Logan体模图像；(b)表示采用斜波滤波反投影方法(FBP算法，Ramp窗)重建的图像；(c)表示采用汉宁滤波反投影方法(FBP算法，Hanning窗)重建的图像；(d)表示采用非单调最小化方法(NTVM方法)重建的图像；(e)表示采用本发明的方法重建的图像。

图6是对应图5数值Shepp-Logan体模重建结果的局部放大图；每行分别对应一个区域，每行的第一列分别是真实的Shepp-Logan体模图像，第二列是采用斜波滤波反投影方法(FBP算法，Ramp窗)重建的图像；第三列是采用汉宁滤波反投影方法(FBP算法，Hanning窗)重建的图像；第四列是采用非单调最小化方法(NTVM方法)重建的图像；第五列是采用本发明的方法重建的图像。

图7是对应图5中不同算法重建图像的线剖面图。

具体实施方式

结合以下实施例对本发明作进一步描述。

实施例1。

提供一种基于全广义变分正则化的低剂量X射线CT图像重建方法，通过以下步骤进行的。

(1)获取CT设备的系统参数和低剂量扫描协议下对数变换后的投影数据q。获取的CT设备的系统参数包括X射线入射光子强度I₀等。

(2)对步骤(1)获取的投影数据q进行Anscombe变换，将服从复合Poission分布的投影数据q转化为近似服从方差为1的Gaussian分布数据u。

对步骤(1)获取的投影数据q进行Anscombe变换，计算式如下：

u＝(u₁,u₂,…,u_N)^T表示投影数据经过Anscombe变化后近似服从方差为1的Gaussian分布数据，其中，T表示转置运算，u₁、u₂、…、u_N是u的分量，N为分量的个数；令f＝(f₁,f₂,…,f_N)^T为待估计投影数据的理想值，f₁、f₂、…、f_N是f的分量。

(3)对步骤(2)获取的数据u建立基于全广义变分最小化的理想数据恢复模型，使用Chambolle-Pock算法求解得到恢复后的投影数据f。

对步骤(2)获取的数据u建立全广义变分最小化模型，具体通过如下方法进行：

(3-1)投影数据恢复的一般模型通过式(2)表示：

其中，第一项为数据保真项，第二项为正则化项，λ＞0是正则化参数。

(3-2)通过全广义变分正则化代替式(2)中的正则化项R(f)，得到基于全广义变分最小化的理想数据恢复模型为：

TGV(f)是f的全广义变分，

$TGV\left(f\right)=\underset{\mathrm{\ω}}{min}{\mathrm{\α}}_1{\∫}_{\mathrm{\Ω}}\left|\▿f-\mathrm{\ω}\right|dx+{\mathrm{\α}}_0{\∫}_{\mathrm{\Ω}}\left|\mathrm{\ϵ}\left(\mathrm{\ω}\right)\right|dx;$

其中α₀,α₁>0是两个正数，是f的梯度，ω是对偶变量， 是ω的梯度。

式(3)的离散形式为：

其中，F＝R^N×N,W＝R^2N×N以及微分算子ε，由限差分算子近似得到，R^N×N表示N×N维的实空间。

根据对偶原理，将式(4)转化为式(5)的鞍点问题：

其中，p和q为对偶变量，并且P和Q定义为：

P＝{p∈R^2N×N|||p||_∞≤α₁}……式(6)；

Q＝{q∈R^3N×N|||q||_∞≤α₀}……式(7)；

定义和如下：

是邻近算子，其由下式表示：

其中，τ是步长，τ>0。

(4)对步骤(3)得到的恢复后的投影数据f进行Anscombe逆变换，再通过滤波反投影算法得到CT重建图像。

本发明的一种基于全广义变分正则化的低剂量X射线CT图像重建方法，构造全广义变分正则化项，可以在投影数据不满足分段常数假设的前提下去除噪声。最后通过Anscombe逆变换和经典的滤波反投影算法对恢复后的投影数据进行解析重建。数值体模实验结果表明，本发明可以有效地抑制低剂量CT图像中的噪声和条形伪影,同时可以很好地保持图像的结构信息和空间分辨率。

实施例2。

为了对本发明的效果做进一步验证，采用图2所示的Clock数字体模图像和图5所示的Shepp-Logan数字体模图像作为本发明的计算机仿真实验对象。

(1)体模图像大小设定位512×512，模拟CT机的X射线源到旋转中心和探测器的距离分别为570mm和1040mm，旋转角在[0,2π]间采样值为1160，每个采样角对应672个探测器单元，探测器单元的大小为1.407mm。通过CT系统仿真生成大小为1160×672的投影数据q，对于Clock体模，入射光子数为5.0×104；相应的Shepp-Logan体模的入射光子数为2.5×10⁵。需要说明的是，在实际的CT数据采集中，投影数据和系统参数即入射光子强度I₀可以直接获取。

(2)对步骤(1)获取的投影数据q进行Anscombe变换，将服从复合Poission分布的投影数据q转化为近似服从方差为1的Gaussian分布数据u。

对步骤(1)获取的投影数据q进行Anscombe变换，计算式如下：

u＝(u₁,u₂,…,u_N)^T表示投影数据经过Anscombe变化后近似服从方差为1的Gaussian分布数据，其中，T表示转置运算，u₁、u₂、…、u_N是u的分量，N为分量的个数；令f＝(f₁,f₂,…,f_N)^T为待估计投影数据的理想值，f₁、f₂、…、f_N是f的分量。

(3)对步骤(2)获取的数据u建立基于全广义变分最小化的理想数据恢复模型，使用Chambolle-Pock算法求解得到恢复后的投影数据f。

对步骤(2)获取的数据u建立全广义变分最小化模型，具体通过如下方法进行：

(3-1)投影数据恢复的一般模型通过式(2)表示：

其中，第一项为数据保真项，第二项为正则化项，λ＞0是正则化参数。

(3-2)通过全广义变分正则化代替式(2)中的正则化项R(f)，得到基于全广义变分最小化的理想数据恢复模型为：

TGV(f)是f的全广义变分，

$TGV\left(f\right)=\underset{\mathrm{\ω}}{min}{\mathrm{\α}}_1{\∫}_{\mathrm{\Ω}}\left|\▿f-\mathrm{\ω}\right|dx+{\mathrm{\α}}_0{\∫}_{\mathrm{\Ω}}\left|\mathrm{\ϵ}\left(\mathrm{\ω}\right)\right|dx;$

其中α₀,α₁>0是两个正数，是f的梯度，ω是对偶变量， 是ω的梯度。

式(3)的离散形式为：

其中，F＝R^N×N,W＝R^2N×N以及微分算子ε，由限差分算子近似得到，R^N×N表示N×N维的实空间。

根据对偶原理，将式(4)转化为式(5)的鞍点问题：

其中，p和q为对偶变量，并且P和Q定义为：

P＝{p∈R^2N×N|||p||_∞≤α₁}……式(6)；

Q＝{q∈R^3N×N|||q||^∞≤α₀}……式(7)；

定义和如下：

是邻近算子，其由下式表示：

其中，τ是步长，τ>0。

(4)对步骤(3)得到的恢复后的投影数据f进行Anscombe逆变换，再通过滤波反投影算法得到CT重建图像。

为了验证全广义变分最小化方法在低剂量CT图像重建中的有效性，通过Clock和Shepp-Logan数值体模进行定量和定性分析。并且与斜波滤波反投影算法、汉宁滤波反投影算法和非单调全变分最小化算法进行比较。CT成像几何采用弧形探测器的扇形束,其中射线源到旋转中心和探测器的距离分别为570mm和1040mm，旋转角在[0,2π]间的采样值为1160，探测器个数为672。体模大小设定为512×512，对于Clock体模，入射光子数为5.0×10⁴；相应的Shepp-Logan体模的入射光子数为2.5×10⁵。

图2是本发明实施例2的数值Clock体模及采用不同方法重建的结果。其中，(a)是真实的Clock体模图像；(b)表示采用斜波滤波反投影方法(FBP算法，Ramp窗)重建的图像；(c)表示采用汉宁滤波反投影方法(FBP算法，Hanning窗)重建的图像；(d)表示采用非单调最小化方法(NTVM方法)重建的图像；(e)表示采用本发明的方法重建的图像。通过图2可以看到，斜波滤波反投影和汉宁滤波反投影算法重建的图像含有大量的噪声和条形伪影，图像质量严重退化。非单调全变分最小化算法重建的图像虽然在一定程度上减少了CT图像的条形伪影和噪声。从视觉效果评价来看，本发明在去除图像的噪声和条形伪影方面具有更佳的表现。为了进一步展现本方法重建效果，由图3中对应图2的局部放大图可以清晰的看出，本发明对噪声和伪影的抑制效果。

通过计算重建图像的信噪比和均方误差可以定量地分析不同方法的差异。信噪比计算公式为：

$SNR=10{\log }_{10}\left(\frac{{\Σ}_i{\left(F_{rec}\right(i)-{\stackrel{\‾}{F}}_{rec})}^2}{{\Σ}_i{\left(F_{rec}\right(i)-F_{ture}(i\left)\right)}^2}\right);$

均方误差计算公式为：

$NMSE=\left(\frac{{\Σ}_i{\left(F_{rec}\right(i)-F_{ture}(i\left)\right)}^2}{{\Σ}_i{\left(F_{ture}\right(i\left)\right)}^2}\right);$

其中，F_rec(i)表示像素点i通过算法得到的重建图像的灰度值，表示重建图像F_rec中所有像素点灰度的平均值；F_ture(i)表示真实图像像素点i处的灰度值。由表1的性噪比和均方误差分析可以看出，相对于滤波反投影和非单调最小化方法，本发明提出的重建图像的方法具有最高的信噪比和最小的均方误差。

表1 图2中重建图像的信噪比和均方误差

为了进一步验证本文提出方法的优越性，分别给出了四种方法重建结果的线剖面图，如图4所示，通过图4可以看出本文提出的重建图像的方法与其他三种方法相比，能够较好的接近真实图像的线剖面。

图5是本发明实施例2的数值Shepp-Logan体模及采用不同方法重建的结果。其中，(a)是真实的Shepp-Logan体模图像；(b)表示采用斜波滤波反投影方法(FBP算法，Ramp窗)重建的图像；(c)表示采用汉宁滤波反投影方法(FBP算法，Hanning窗)重建的图像；(d)表示采用非单调最小化方法(NTVM方法)重建的图像；(e)表示采用本发明的方法重建的图像。通过图5可以看到，相比滤波反投影和非单调全变分最小化方法，本发明提出的全广义变分最小化方法在抑制噪声与条形伪影和保持边缘方面都有良好表现，可以更好的保持图像一致性。为了更清晰地比较四种方法的噪声和伪影的抑制效果，图6给出了低剂量CT重建结果的一个局部的放大图。从放大图可以明显的看出，本发明提出的全广义变分方法重建的图像对噪声和伪影有更好的抑制作用。

由表2的信噪比和均方误差分析可以看出，与其他三种方法相比，本发明提出的全广义变分方法在提高图像重建的信噪比，降低重建图像均方误差方面均有上佳表现。

表2 图5中重建图像的信噪比和均方误差

此外，为了更直观的验证重建图像的质量，画出重建图像的线剖面图，如图7所示，可以看出本发明重建图像的方法与真实图像的线剖面有更高的吻合度，体现出了全广义变分方法的优越性。

本发明在克服低剂量CT重建的全变分(TV)正则化方法的不足基础上提出了一个低剂量CT优质重建的全广义变分(TGV)方法性。本发明构造全广义变分正则化项，可以在投影数据不满足分段常数假设的前提下去除噪声。最后通过Anscombe逆变换和经典的滤波反投影算法对恢复后的投影数据进行解析重建。数值体模实验结果表明，本发明可以有效地抑制低剂量CT图像中的噪声和条形伪影,同时可以很好地保持图像的结构信息和空间分辨率。

最后应当说明的是，以上实施例仅用以说明本发明的技术方案而非对本发明保护范围的限制，尽管参照较佳实施例对本发明作了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的实质和范围。

Claims (5)

1.一种基于全广义变分正则化的低剂量X射线CT图像重建方法，其特征在于，包括如下步骤：

(1)获取CT设备的系统参数和低剂量扫描协议下对数变换后的投影数据q；

(2)对步骤(1)获取的投影数据q进行Anscombe变换，将服从复合Poission分布的投影数据q转化为近似服从方差为1的Gaussian分布数据u；

(3)对步骤(2)获取的数据u建立基于全广义变分最小化的理想数据恢复模型，使用Chambolle-Pock算法求解得到恢复后的投影数据f；

(4)对步骤(3)得到的恢复后的投影数据f进行Anscombe逆变换，再通过滤波反投影算法得到CT重建图像。

2.根据权利要求1所述的基于全广义变分正则化的低剂量X射线CT图像重建方法，其特征在于，所述步骤(1)中获取的CT设备的系统参数包括X射线入射光子强度I₀。

3.根据权利要求2所述的基于全广义变分正则化的低剂量X射线CT图像重建方法，其特征在于，

所述步骤(2)中对步骤(1)获取的投影数据q进行Anscombe变换，计算式如下：

u＝(u₁,u₂,…,u_N)^T表示投影数据经过Anscombe变化后近似服从方差为1的Gaussian分布数据，其中，T表示转置运算，u₁、u₂、…、u_N是u的分量，N为分量的个数；令

f＝(f₁,f₂,…,f_N)^T为待估计投影数据的理想值，f₁、f₂、…、f_N是f的分量。

4.根据权利要求3所述的基于全广义变分正则化的低剂量X射线CT图像重建方法，其特征在于，

所述步骤(3)中对步骤(2)获取的数据u建立全广义变分最小化模型，具体通过如下方法进行：

(3-1)投影数据恢复的一般模型通过式(2)表示：

其中，第一项为数据保真项，第二项为正则化项，λ＞0是正则化参数；

(3-2)通过全广义变分正则化代替式(2)中的正则化项R(f)，得到基于全广义变分最小化的理想数据恢复模型为：

TGV(f)是f的全广义变分，

$TGV\left(f\right)=\underset{\mathrm{\ω}}{min}{\mathrm{\α}}_1{\∫}_{\mathrm{\Ω}}|\▿f-\mathrm{\ω}|dx+{\mathrm{\α}}_0{\∫}_{\mathrm{\Ω}}\left|\mathrm{\ϵ}\left(\mathrm{\ω}\right)\right|dx;$

其中α₀,α₁>0是两个正数，是f的梯度，ω是对偶变量， 是ω的梯度。

5.根据权利要求4所述的基于全广义变分正则化的低剂量X射线CT图像重建方法，其特征在于，

式(3)的离散形式为：

其中，F＝R^N×N,W＝R^2N×N以及微分算子ε，由限差分算子近似得到，R^N×N表示N×N维的实空间；

根据对偶原理，将式(4)转化为式(5)的鞍点问题：

其中，p和q为对偶变量，并且P和Q定义为：

P＝{p∈R^2N×N|||p||_∞≤α₁}……式(6)；

Q＝{q∈R^3N×N|||q||_∞≤α₀}……式(7)；

定义和如下：

是邻近算子，其由下式表示：

其中，其中，τ是步长，τ>0。