CN106127825A - 一种基于广义惩罚加权最小二乘的x射线ct图像重建方法 - Google Patents

Abstract

一种基于广义惩罚加权最小二乘的X射线CT图像重建方法，包括如下步骤：(1)获取CT设备的系统参数和低剂量扫描协议下的投影数据q_e；(2)对(1)中获取的投影数据q_e进行系统校正和对数变换后得到投影数据y；对y进行逐个数据点上的方差估计，获得每个点上的方差集合(3)根据(2)中得到的方差建立广义惩罚加权最小二乘的低剂量CT重建模型；(4)对(3)中的低剂量CT重建模型进行求解，建立算法的全局收敛性以收敛点作为解，再进行图像重建，得到最终的CT重建图像。本发明可以有效地抑制低剂量CT图像中的噪声和条形伪影,同时可以很好地保持图像的结构信息和空间分辨率。

Description

一种基于广义惩罚加权最小二乘的X射线CT图像重建方法

技术领域

本发明涉及医学影像的计算机断层成像技术领域，特别是涉及一种基于广义惩罚加权最小二乘的X射线CT图像重建方法。

背景技术

X射线计算机断层成像技术是现代医学影像学的杰出代表,已经广泛应用于临床诊断和治疗。然而偏高的X射线照射剂量可诱发癌症、白血病或者其他遗传性疾病。如何最大限度地降低X射线使用剂量已经成为医学CT成像领域研究的关键技术。

当前,临床上采用的低剂量CT扫描方式虽然可以在一定程度上减少X射线照射剂量，但低剂量的投影数据将受到量子噪声和电子噪声污染，而采用滤波反投影算法重建的CT图像质量严重退化，直接影响临床诊断的准确性。

对于这一问题有两种传统的后处理解决方法:一是直接对低剂量CT图像进行滤波，以减少图像的噪声和伪影，属于图像后处理技术；二是根据投影数据满足的统计学规律，完成基于统计的CT图像迭代重建。第一种的处理技术，直接减少图像的噪声和伪影。图像后处理技术因其简单且易于操作已在低剂量CT图像恢复中得到广泛的研究和应用。由于低剂量CT图像中噪声和伪影分布的复杂性，使得高精度的滤波器设计困难极大。第二种通过CT系统建模，构建图像重建模型，通过优化目标函数实现图像重建。相对于经典的FBP算法，迭代重建算法通过系统建模(系统光学模型和系统统计模型)对CT成像几何、X射线的能谱特性、射束硬化效应、散射和噪声特性进行准确描述，而且易于加入先验信息约束，因此特别适合低剂量CT图像优质重建。统计迭代重建在抑制图像噪声和伪影以及提高空间分辨率等方面都有上佳表现，是当前低剂量CT成像领域的研究热点。然而，由于统计迭代重建需要反复进行投影与反投影运算，且CT图像数据量庞大，导致CT图像重建速度特别慢，难以满足临床中实时交互的需求。

不同于上述两种方法的另一种策略是，根据投影数据的噪声统计特性建立数据恢复模型，得到更接近于理想值的投影数据后采用FBP算法重建出CT图像。该类方法不仅可以避免统计迭代重建方法在计算速度上的劣势，而且考虑了投影数据的噪声统计特性，去除噪声的同时可以较好地保持图像的空间分辨率。诸多研究中，一般都是将投影数据的噪声看作高斯或者泊松噪声来处理，大量的临床低剂量CT实验和理论分析，提出投影数据(对数变换后)的噪声近似满足均值与方差呈非线性关系的高斯分布，基于这一重要特点，多种低剂量CT投影数据恢复的惩罚加权最小二乘(penalized weighted least-squares,简记为PWLS)方法相继被提出，例如基于K—L变换的PWLS方法，多尺度PWLS方法等。这种处理方法相比传统的滤波反投影方法能够在一定程度上改善CT重建图像的质量，但重建的图像在对比度和信噪比等方面仍然存在一定的不足，从而影响临床诊断结果。

因此，针对现有技术不足，提供一种可以有效地抑制低剂量CT图像中的噪声和伪影，同时较好地保持图像的空间分辨率的基于广义惩罚加权最小二乘(generalizedpenalized weighted least-squares,简记为GPWLS)的X射线CT图像重建方法以克服现有技术不足甚为必要。

发明内容

本发明的目的在于避免现有技术的不足之处而提供一种基于广义惩罚加权最小二乘的X射线CT图像重建方法，该方法可以有效地抑制低剂量CT图像中的噪声和伪影，同时较好地保持图像的空间分辨率。

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

提供一种基于广义惩罚加权最小二乘的X射线CT图像重建方法，包括如下步骤：

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

(2)对(1)中获取的投影数据q_e进行系统校正和对数变换后得到投影数据y，y＝{y₁,y₂,...,y_M},i＝1,2,L M，y_i是投影数据y的第i个分量，M是正整数，M是投影数据y的分量的个数；对y进行逐个数据点上的方差估计，获得每个点上的方差集合

(3)根据(2)中得到的的方差建立广义惩罚加权最小二乘的低剂量CT重建模型；

(4)对(3)中的低剂量CT重建模型进行求解，建立算法的全局收敛性以收敛点作为解，再进行图像重建，得到最终的CT重建图像。

优选的，上述步骤(1)中获取的CT设备的系统参数包括X射线入射光子强度I₀和系统电子噪声的方差

优选的，上述步骤(2)中的对y进行逐个数据点上的方差估计具体是采用基于小邻域图像的局部方差估计或者基于CT投影数据噪声特性的方差估计。

优选的，上述步骤(3)中建立的低剂量CT重建模型为:

$\underset{q\≥0}{min}{(y-q)}^T{G}^{-1}(y-q)+{\mathrm{\β}}_{2}R\left(q\right)......\left(I\right);$

其中，“T”表示转置运算，I为单位矩阵，β₁、β₂为两个超参数，β₁>0，β₂>0,R(q)是惩罚项，R(q)的表达式为:

$R\left(q\right)=\frac{1}{2}\underset{i}{\Σ}\underset{m\∈N_i}{\Σ}{\mathrm{\ω}}_{im}{(q_i-q_m)}^2......\left(II\right);$

其中，N_i为投影数据中第i个像素的一阶邻域，q_i为投影数据中第i个像素值，q_m为投影数据中第i个像素的邻域的内像素值，ω_im表示q_i的一阶邻域内的q_m所占的权重，并且ω_im在水平方向取值为1，在垂直方向取值为0.25。

优选的，上述步骤(4)具体是通过交替方向优化算法对(3)中的低剂量CT重建模型进行求解。

优选的，上述步骤(4)具体是通过交替方向优化算法对(3)中的低剂量CT重建模型进行求解，具体是：

引入变量p，得到

$\underset{p}{min}{(y-p)}^T{\mathrm{\Σ}}^{-1}(y-p)+{\mathrm{\β}}_1\left|\right|p-q|{|}_2^2={(y-q)}^T{G}^{-1}(y-q)......\left(IV\right);$

其中，

根据式(Ⅳ)，将式(Ⅰ)写成如下形式：

$\underset{p,q}{min}J(p,q)={(y-p)}^T{\mathrm{\Σ}}^{-1}(y-p)+{\mathrm{\β}}_1\left|\right|p-q|{|}_2^2+{\mathrm{\β}}_{2}R\left(q\right)......\left(V\right);$

其中，J(p,q)是J是关于p和q的函数；

(4.1)给定初始点q⁰＝y，令k＝0，k是迭代次数,k+1表示第k+1次迭代；

(4.2)使用交替方向法依次求p^k+1和q^k+1，

$\left\{\begin{array}{cc}{S}_G\left(q^k\right)=p^{k+1}=\underset{p}{\arg min}{(y-p)}^T{\mathrm{\Σ}}^{-1}(y-p)+{\mathrm{\β}}_1\left|\right|p-q^k|{|}_2^2& ......\left(a\right);\\ S_R\left(p^{k+1}\right)=q^{k+1}=\underset{q}{\arg min}{\mathrm{\β}}_1\left|\right|p^{k+1}-q|{|}_2^2+{\mathrm{\β}}_{2}R\left(q\right)& ......\left(b\right);\end{array}\right.\mathrm{......}\left(VI\right);$

其中，S_G(q^k)表示式(Ⅴ)中对p的优化解，S_R(p^k+1)表示式(Ⅴ)中q的优化解；

(4.3)判断是否满足迭代终止条件：如果满足迭代终止条件，则停止迭代，进入步骤(4.5)，否则进入步骤(4.4)；

(4.4)令k＝k+1，进入步骤(4.2)；

(4.5)式Ⅵ(a)中的目标函数为光滑的二次凸函数，对Ⅵ(a)式求导并令其为零，得到Ⅵ(a)的解为

对于式Ⅵ(b)，使用Gauss-Seidel算法求解，其对应的迭代格式为其中 表示像素q_i左上方的领域，表示像素q_i右上方的领域，N_i表示像素的一阶领域。

优选的，上述步骤(4)中，具体通过滤波反投影方法进行图像重建，得到最终的CT重建图像。

本发明的基于广义惩罚加权最小二乘的X射线CT图像重建方法，包括如下步骤：(1)获取CT设备的系统参数和低剂量扫描协议下的投影数据q_e；(2)对(1)中获取的投影数据q_e进行系统校正和对数变换后得到投影数据y，y＝{y₁,y₂,...,y_M},i＝1,2,L M，y_i是投影数据y的第i个元素，M是正整数，M是投影数据y的分量的个数；对y进行逐个数据点上的方差估计，获得每个点上的方差集合(3)根据(2)中得到的的方差建立广义惩罚加权最小二乘的低剂量CT重建模型；(4)对(3)中的低剂量CT重建模型进行求解，建立算法的全局收敛性以收敛点作为解，再进行图像重建，得到最终的CT重建图像。该基于广义惩罚加权最小二乘的X射线CT图像重建方法，构造广义加权最小二乘保真项，并且将惩罚先验信息引入投影数据的恢复过程中，从而达到抑制噪声的目的，最后使用经典的滤波反投影算法对恢复后的投影数据进行解析重建。本发明可以有效地抑制低剂量CT图像中的噪声和条形伪影,同时可以很好地保持图像的结构信息和空间分辨率。

附图说明

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

图1是本发明一种基于广义惩罚加权最小二乘的X射线CT图像重建方法的流程图。

图2为本发明实施例2中用到的数值体模，其中(a)为Clock体模，(b)为Shepp-Logan体模。

图3为数值Clock体模的重建结果，其中(a)为FBP算法(Ramp窗)的重建结果；(b)为FBP算法(Hanning窗)的重建结果；(c)为PWLS方法的重建结果；(d)为本发明方法的重建结果。

图4为数值Shepp-Logan体模的重建结果，其中(a)为FBP算法(Ramp窗)的重建结果；(b)为FBP算法(Hanning窗)的重建结果；(c)为PWLS方法的重建结果；(d)为本发明方法的重建结果。

图5为对应图4中不同方法重建图像的水平剖面图，其中(a)为FBP算法(Ramp窗)重建图像的水平剖面图；(b)为FBP算法(Hanning窗)重建图像的水平剖面图；(c)为PWLS方法重建图像的水平剖面图；(d)为GPWLS方法重建图像的水平剖面图。

图6为对应图4中不同方法重建图像感兴趣区域的局部放大图，其中(a)为FBP算法(Ramp窗)重建图像感兴趣区域的局部放大图；(b)为FBP算法Hanning窗)重建图像感兴趣区域的局部放大图；(c)为PWLS方法重建图像感兴趣区域的局部放大图；(d)为GPWLS方法重建图像感兴趣区域的局部放大图。

具体实施方式

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

实施例1。

一种基于广义惩罚加权最小二乘的X射线CT图像重建方法，如图1所示，包括如下步骤。

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

(2)对(1)中获取的投影数据q_e进行系统校正和对数变换后得到投影数据y，y＝{y₁,y₂,...,y_M},i＝1,2,L M，y_i是投影数据y的第i个分量，M是正整数，M是投影数据y的分量的个数；对y进行逐个数据点上的方差估计，获得每个点上的方差集合

(3)根据(2)中得到的的方差建立广义惩罚加权最小二乘的低剂量CT重建模型，简称构建GPWLS模型；

(4)对(3)中的低剂量CT重建模型进行求解，建立算法的全局收敛性以收敛点作为解，再通过滤波反投影方法进行图像重建，得到最终的CT重建图像。

具体的，上述步骤(1)中获取的CT设备的系统参数包括X射线入射光子强度I₀和系统电子噪声的方差

具体的，上述步骤(2)中的对y进行逐个数据点上的方差估计具体是采用基于小邻域图像的局部方差估计或者基于CT投影数据噪声特性的方差估计。

具体的，上述步骤(3)中建立的低剂量CT重建模型为:

$\underset{q\≥0}{min}{(y-q)}^T{G}^{-1}(y-q)+{\mathrm{\β}}_{2}R\left(q\right)......\left(I\right);$

其中，“T”表示转置运算，I为单位矩阵，β₁、β₂为两个超参数，β₁>0，β₂>0,R(q)是惩罚项，R(q)的表达式为:

$R\left(q\right)=\frac{1}{2}\underset{i}{\Σ}\underset{m\∈N_i}{\Σ}{\mathrm{\ω}}_{im}{(q_i-q_m)}^2......\left(II\right);$

其中，N_i为投影数据中第i个像素的一阶邻域，q_i为投影数据中第i个像素值，q_m为投影数据中第i个像素的邻域的内像素值，ω_im表示q_i的一阶邻域内的q_m所占的权重，并且ω_im在水平方向取值为1，在垂直方向取值为0.25。

具体的，上述步骤(4)具体是通过交替方向优化算法对(3)中的低剂量CT重建模型进行求解。

步骤(4)是通过交替方向优化算法对(3)中的低剂量CT重建模型进行求解，具体是：

引入变量p，得到

$\underset{p}{min}{(y-p)}^T{\mathrm{\Σ}}^{-1}(y-p)+{\mathrm{\β}}_1\left|\right|p-q|{|}_2^2={(y-q)}^T{G}^{-1}(y-q)......\left(IV\right);$

其中，

根据式(Ⅳ)，将式(Ⅰ)写成如下形式：

$\underset{p,q}{min}J(p,q)={(y-p)}^T{\mathrm{\Σ}}^{-1}(y-p)+{\mathrm{\β}}_1\left|\right|p-q|{|}_2^2+{\mathrm{\β}}_{2}R\left(q\right)......\left(V\right);$

其中，J(p,q)是J是关于p和q的函数；

(4.1)给定初始点q⁰＝y，令k＝0，k是迭代次数,k+1表示第k+1次迭代；

(4.2)使用交替方向法依次求p^k+1和q^k+1，

$\left\{\begin{array}{cc}{S}_G\left(q^k\right)=p^{k+1}=\underset{p}{\arg min}{(y-p)}^T{\mathrm{\Σ}}^{-1}(y-p)+{\mathrm{\β}}_1\left|\right|p-q^k|{|}_2^2& ......\left(a\right);\\ S_R\left(p^{k+1}\right)=q^{k+1}=\underset{q}{\arg min}{\mathrm{\β}}_1\left|\right|p^{k+1}-q|{|}_2^2+{\mathrm{\β}}_{2}R\left(q\right)& ......\left(b\right);\end{array}\right.\mathrm{......}\left(VI\right);$

其中，S_G(q^k)表示式(Ⅴ)中对p的优化解，S_R(p^k+1)表示式(Ⅴ)中q的优化解；

(4.3)判断是否满足迭代终止条件：如果满足迭代终止条件，则停止迭代，进入步骤(4.5)，否则进入步骤(4.4)；

(4.4)令k＝k+1，进入步骤(4.2)；

(4.5)式Ⅵ(a)中的目标函数为光滑的二次凸函数，对Ⅵ(a)式求导并令其为零，得到Ⅵ(a)的解为

对于式Ⅵ(b)，使用Gauss-Seidel算法求解，其对应的迭代格式为其中 表示像素q_i左上方的领域，表示像素q_i右上方的领域，N_i表示像素的一阶领域。

下面给出GPWLS方法的全局收敛性分析。根据最优化理论，非扩张渐近正则算子Γ以及Γ的不动点解集非空，则对任意初始点q₀都存在序列{Γ(q^k)}弱收敛到Γ的不动点。因此，要分析GPWLS方法的收敛性，首先应该证明是Γ非扩张渐近正则算子,以及目标函数的解集非空，下面给出几个相关的概念。

定义1如果对任意的x₁,x₂∈R^N都有||P(x₁)-P(x₂)||≤||x₁-x₂||，则P称为非扩张算子。如果存在非扩张算子Α以及κ∈(0,1)使得P＝(1-κ)I+κΑ,则称是κ—平均非扩张算子。

定义2设C是Banach空间X上的一个闭凸集，如果对任意的x∈C,F:C→C满足：

则称F:C→C是渐近正则的。

定义3设函数如果至少存在一点使得并且对任意x∈X都有则称在X上是适定的。如果对X中满足||x_k||→∞的任意子列都则称在X上是强制的。

在给出GPWLS方法的收敛性定理之前，首先给出如下几个引理。

引理1设是凸的下半连续的函数以及β>0，定义如下：对任意x，定义则S为—平均非扩张算子。

引理2公式Γ(·)＝S_R(S_G(·))定义的算子Γ是非扩张的。

证明：由于(Σ^-1+β₁I)是对称正定矩阵且其特征值大于β₁，则对任意的u，v有：

||β₁(Σ^-1+β₁I)^-1(u-v)||₂≤||u-v||₂；

因此，

$\begin{array}{c}\left|\right|\mathrm{\Γ}\left(u\right)-\mathrm{\Γ}\left(v\right)|{|}_2=||S_R\left(S_G\right(u\left)\right)-S_R\left(S_G\right(v\left)\right)|{|}_2\\ \≤\left|\right|S_R\left(u\right)-S_G\left(v\right)|{|}_2\\ \≤\left|\right|{\mathrm{\β}}_1{({\Σ}^{-1}+{\mathrm{\β}}_{1}I)}^{-1}(u-v)|{|}_2\\ \≤\left|\right|u-v|{|}_2\end{array};$

故引理得证。

引理3公式q^k+1＝Γ(q^k)产生的点列{q_k}满足

证明：在点(p_k,q_k)第一个分量p_k+1处对目标函数进行泰勒展开可以得到：

$\begin{array}{c}J(p_k,q_k)=J(p_{k+1},q_k)+{(p_k-p_{k+1})}^T\frac{\∂J}{\∂p}(p_{k+1},p_k)\\ +\frac{1}{2}{(p_k-p_{k+1})}^T\frac{{\∂}^{2}J}{\∂p^2}(p_k-p_{k-1})\end{array};$

由于，

p_k+1＝argmin_p J(p,q_k)，

则，

$\frac{\∂J}{\∂p}(p_{k+1},q_k)=0;$

进一步,

$\frac{{\∂}^{2}J}{\∂p^2}(p_{k+1},q_k)={\mathrm{\Σ}}^{-1}+{\mathrm{\β}}_{1}I;$

由于Σ^-1+β₁I对称正定且其特征值大于β₁,则：

$\begin{array}{c}J(p_k,q_k)-J(p_{k+1},q_k)=\frac{1}{2}{(p_k-p_{k+1})}^T({\mathrm{\Σ}}^{-1}+{\mathrm{\β}}_{1}I)(p_k-p_{k+1})\\ \≥{\mathrm{\β}}_1\left|\right|p_k-p_{k+1}|{|}_2^2\end{array};$

由于S_R(q)为非扩张算子,则：

$\left|\right|p_k-p_{k+1}|{|}_2^2\≥||S_R\left(p_k\right)-S_R\left(p_{k+1}\right)||=||q_k-q_{k+1}|{|}_2^2;$

可以得到：

$J(p_k,q_k)-J(p_{k+1},q_{k+1})\≥{\mathrm{\β}}_1\left|\right|q_k-q_{k+1}|{|}_2^2;$

因此收敛，则

故引理得证。

引理4对任意初始值q₀，设{q_k}是由q^k+1＝Γ(q^k)产生的序列,则Γ为渐近正则算子。

证明：由Γ的定义，则：

$\begin{array}{c}\underset{k\→\mathrm{\∞}}{\lim }\left|\right|{\mathrm{\Γ}}^{k+1}\left(q_0\right)-{\mathrm{\Γ}}^k\left(q_0\right)|{|}_2=\underset{k\→\mathrm{\∞}}{\lim }||{\mathrm{\Γ}}^{k+1}\left(q_1\right)-{\mathrm{\Γ}}^k\left(q_1\right)|{|}_2\\ M\\ =\underset{k\→\mathrm{\∞}}{\lim }\left|\right|{\mathrm{\Γ}}^{k+1}\left(q_k\right)-{\mathrm{\Γ}}^k\left(q_{k-1}\right)|{|}_2\\ =\underset{k\→\mathrm{\∞}}{\lim }\left|\right|q_{k+1}-q_k|{|}_2\\ =0\end{array};$

故引理得证。

引理5由公式(Ⅴ)定义的函数J(p,q)是强制的。

证明：由R(q)的定义，我们可以得到：

$R\left(q\right)=\underset{i}{\Σ}\underset{j}{\Σ}{\mathrm{\ω}}_{ij}{(q_i-q_j)}^2\≥{\mathrm{\ω}}_{min}\left|\right|Dq|{|}_2^2;$

其中，ω_min＝min{ω_ij}，D_h，D_v分别为水平和垂直方向的差分矩阵，注意到：

$\begin{array}{c}J(p,q)={(y-p)}^T{\Σ}^{-1}(y-p)+{\mathrm{\β}}_1\left|\right|p-q|{|}_2^2+{\mathrm{\β}}_{2}R\left(q\right)\\ \≥\mathrm{\γ}\left|\right|p-y|{|}_2^2+{\mathrm{\β}}_1||p-q|{|}_2^2+{\mathrm{\β}}_2{\mathrm{\ω}}_{\min }\left|\right|Dq|{|}_2^2\\ =\left|\right|\left(\begin{array}{cc}\sqrt{\mathrm{\γ}}I& 0\\ \sqrt{{\mathrm{\β}}_1}I& -\sqrt{{\mathrm{\β}}_1}I\end{array}\right)\left(\begin{array}{c}p\\ q\end{array}\right)-\left(\begin{array}{c}y\\ 0\end{array}\right)|{|}_2^2+{\mathrm{\β}}_2{\mathrm{\ω}}_{\min }|\left|\left(\begin{array}{cc}0& D\end{array}\right)\left(\begin{array}{c}p\\ q\end{array}\right)\right|{|}_2^2\end{array};$

其中，进一步，

$\left(\begin{array}{cc}\sqrt{\mathrm{\γ}}I& 0\\ \sqrt{{\mathrm{\β}}_1}I& -\sqrt{{\mathrm{\β}}_1}I\\ 0& D\end{array}\right)=\left(\begin{array}{ccc}I& 0& 0\\ 0& I& 0\\ 0& -\frac{1}{\sqrt{{\mathrm{\β}}_1}}D& I\end{array}\right)\left(\begin{array}{cc}\sqrt{\mathrm{\γ}}I& 0\\ \sqrt{{\mathrm{\β}}_1}I& -\sqrt{{\mathrm{\β}}_1}I\\ D& 0\end{array}\right);$

注意到上述矩阵是满秩，因此，当时，J(p,q)→∞，

故引理得证。

引理6公式Γ(·)＝S_R(S_G(·))定义的Γ的不动点解集非空。

证明：由于J是强制的，因此J的极小值构成的集合非空，设(p^*,q^*)是J(x,u)的一个极小值即，

$\frac{\∂J}{\∂p}(p^{*},q^{*})=0,\frac{\∂J}{\∂q}(p^{*},q^{*})=0.\left(2.32\right);$

因此，我们可以得到

$\left\{\begin{array}{c}{p}^{*}=S_G\left(q^{*}\right)=\arg \min J\left(q^{*}\right)\\ q^{*}=S_R\left(p^{*}\right)=\arg \min J\left(p^{*}\right)\end{array}\right.;$

则q^*＝S_R(S_G(q^*))＝Γ(q^*)，即q^*是Γ的一个不动点。故引理得证。

由于目标函数J是严格凸的，因此，Γ的不动点即是J的全局最优点。根据最优化理论序列{q_k}收敛到J的不动点，即极小值点。下面给出GPWLS算法的全局收敛。

定理1给定初始值q⁰∈Rⁿ，设{q_k}是由q^k+1＝Γ(q^k)产生的序列,则q_k收敛到J的不动点即极小值点。

本发明的基于广义惩罚加权最小二乘的X射线CT图像重建方法，验证了算法的全局收敛性。本发明构造广义加权最小二乘保真项，并且将惩罚先验信息引入投影数据的恢复过程中，从而达到抑制噪声的目的，最后使用经典的滤波反投影算法对恢复后的投影数据进行解析重建。数值体模实验结果表明，本发明可以有效地抑制低剂量CT图像中的噪声和条形伪影,同时可以很好地保持图像的结构信息和空间分辨率。

实施例2。

采用图1所示的Shepp-Logan数字体模图像和Clock数字体模图像作为本发明的计算机仿真实验对象。体模图像大小设定位512×512，模拟CT机的X射线源到旋转中心和探测器的距离分别为570mm和1040mm，旋转角在[0,2π]间采样值为1160，每个采样角对应672个探测器单元，探测器单元的大小为1.407mm。通过CT系统仿真生成大小为1160×672的投影数据，其中X射线的入射光子强度I₀为5.0×104。在实际的CT数据采集中，投影数据和系统参数即入射光子强度I₀和系统电子噪声的方差均可以直接获取。

对步骤1中模拟生成的CT投影数据和系统参数和I₀，通过步骤2系统校正和对数变换后得到的投影数据y，以及方差采用如下公式：

${\mathrm{\σ}}_{y_i}^2=\frac{1}{I_0}\exp \left({\stackrel{\‾}{y}}_i\right)(1+\frac{1}{I_0}\exp ({\stackrel{\‾}{y}}_i\left({\mathrm{\σ}}_e^2-1.25\right)\left)\right)......\left(III\right);$

在本实例中，I₀为第i个数据点的X射线入射光子强度，即I₀为5.0×10⁴。为系统电子噪声的方差。式(Ⅱ)中N_i是投影数据中第i个像素的一阶领域，ω_im在水平方向取值为1，在垂直方向取值为0.25。不同于图像域的二阶邻域相似性,投影数据的相似性仅存在于相邻探测器单元和相邻投影角度这两个维度之间,而且前者的相似概率大于后者,因此水平与垂直方向使用不同的权重。根据步骤3利用本发明所示模型(Ⅰ)对待估计投影数据的理想值q进行求解。

本发明使用(Ⅰ)的广义加权最小二乘代替加权最小二乘。注意到，当β₁充分大时,广义加权最小二乘就等价于加权最小二乘。进一步,引入一个新的变量p，可以得到：

$\underset{p}{min}{(y-p)}^T{\mathrm{\Σ}}^{-1}(y-p)+{\mathrm{\β}}_1\left|\right|p-q|{|}_2^2={(y-q)}^T{G}^{-1}(y-q)......\left(IV\right);$

根据(Ⅳ)式，可以将式(Ⅰ)写成如下形式：

$\underset{p,q}{min}J(p,q):={(y-p)}^T{\mathrm{\Σ}}^{-1}(y-p)+{\mathrm{\β}}_1\left|\right|p-q|{|}_2^2+{\mathrm{\β}}_{2}R\left(q\right)......\left(V\right);$

求解式(Ⅴ)中的目标函数，本发明使用交替方向法依次求p和q，

$\left\{\begin{array}{cc}{S}_G\left(q^k\right):=p^{k+1}=\underset{p}{\arg min}{(y-p)}^T{\mathrm{\Σ}}^{-1}(y-p)+{\mathrm{\β}}_1\left|\right|p-q^k|{|}_2^2& \left(a\right)\\ S_R\left(p^{k+1}\right):=q^{k+1}=\underset{q}{\arg min}{\mathrm{\β}}_1\left|\right|p^{k+1}-q|{|}_2^2+{\mathrm{\β}}_{2}R\left(q\right)& \left(b\right)\end{array}\right.;---\left(VI\right)$

S_G(q^k)表示式子Ⅴ中对p的优化解，S_R(p^k+1)表示式子Ⅴ中q的优化解。

给定初始点p⁰，产生迭代序列：p⁰，p¹，q¹，p²，q²，…，p^k，q^k，…。

满足如下的关系：

q^k+1＝S_R(S_G(q^k))……(Ⅶ)；

可以转化为：

q^k+1＝T(q^k)……(Ⅷ)；

由于Ⅵ(a)式中的目标函数为光滑的二次凸函数，对Ⅵ(a)式求导并令其为零，可以得到Ⅵ(a)的解为对于Ⅵ(b)式，使用Gauss-Seidel算法求解，其对应的迭代格式为其中 表示像素q_i左上方的领域，表示像素q_i右上方的领域，N_i表示像素的一阶领域。

为了验证本发明所示方法的效果，对系统校正和对数变换后的投影数据建立广义惩罚加权最小二乘(GPWLS)模型，然后对估计的投影数据q采用传统的滤波反投影(FBP)算法进行重建，得到重建图像。为了对比本发明技术方案的效果，将本发明方法与Ramp窗FBP算法、Hanning窗FBP以及(PWLS)重建图像的影像进行比较，图3给出了Clock体模不同方法的重建结果。图3(a)为Ramp窗FBP算法重建的结果；图3(b)为Hanning窗FBP算法重建的结果；(c)为PWLS方法重建的结果；(d)为GPWLS方法重建的结果。从视觉效果评价来看，GPWLS方法在抑制图像噪声和条形伪影方面具有更佳的表现。为了进一步定量分析GPWLS方法,计算重建图像的信噪比(SNR)，均方误差(NMSE)和以及重建图像中两个感兴趣区域(如图2(a)中所示)的对比度(CNR)。SNR，NMSE以及CNR的计算公式如下：

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

$NMSE=\frac{{\mathrm{\Σ}}_i{\left(\mathrm{\μ}\right(i)-{\mathrm{\μ}}_{ture}(i\left)\right)}^2}{{\mathrm{\Σ}}_i{\left({\mathrm{\μ}}_{ture}\right(i\left)\right)}^2};$

$CNR=\frac{|{\stackrel{\‾}{\mathrm{\μ}}}_{ROI}-{\stackrel{\‾}{\mathrm{\μ}}}_{BG}|}{\sqrt{{\mathrm{\σ}}_{R0I}^2+{\mathrm{\σ}}_{BG}^2}};$

其中，μ(i)表示图像μ的第i个像素，表示μ的均值，μ_ture(i)表示真实的体模图像μ_ture的第i个像素，表示感兴趣(ROI)区域内图像μ的均值，表示感背景(Background)区域内图像μ的均值，表示ROI区域内图像μ的方差，表示Background区域内图像μ的方差。本发明中选取的ROI区域大小为20×20。由表1的信噪比和均方误差分析可以看出，与其他三种方法相比本发明在提高重建图像的信噪比，降低重建图像均方误差方面均有更佳表现。此外，由表2可以看出，在高对比度区域ROI1，四种方法重建图像的对比度相差不大。然而，在低对比度区域ROI2，GPWLS方法重建的图像具有最高的对比度。

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

表2图3中重建图像的对比度

图4为低剂量的Shepp-Logan投影数据由不同方法重建的结果。图4(a)为Ramp窗FBP算法重建的结果；图4(b)为Hanning窗FBP算法重建的结果；4(c)为PWLS方法重建的结果；4(d)为GPWLS方法重建的结果。从图4可以看出PWLS方法和GPWLS两种方法较好地抑制了图像中的条形伪影。从视觉评价来看，图4(d)中CT图像的噪声水平比图4(c)中CT图像的噪声水平低。

图5是重建结果的水平剖面图，从图中可以看出通过本发明重建的图像更接近真实的体模图像。为了更清晰的比较各种方法对噪声和伪影的抑制效果，图6给出了对应图4中重建结果的感兴趣区域的局部放大图。从图6可以看出，本发明在噪声去除、条形伪影抑制以及结构信息保持方面都比其他方法具有更佳的表现。

本发明的基于广义惩罚加权最小二乘的X射线CT图像重建方法，构造广义加权最小二乘保真项，并且将惩罚先验信息引入投影数据的恢复过程中，从而达到抑制噪声的目的，最后使用经典的滤波反投影算法对恢复后的投影数据进行解析重建。数值体模实验结果表明，本发明可以有效地抑制低剂量CT图像中的噪声和条形伪影,同时可以很好地保持图像的结构信息和空间分辨率。

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

Claims (7)

1.一种基于广义惩罚加权最小二乘的X射线CT图像重建方法，其特征在于，包括如下步骤：

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

(2)对(1)中获取的投影数据q_e进行系统校正和对数变换后得到投影数据y，y＝{y₁,y₂,...,y_M},i＝1,2,L M，y_i是投影数据y的第i个分量，M是正整数，M是投影数据y的分量的个数；对y进行逐个数据点上的方差估计，获得每个点上的方差集合

(3)根据(2)中得到的的方差建立广义惩罚加权最小二乘的低剂量CT重建模型；

(4)对(3)中的低剂量CT重建模型进行求解，建立算法的全局收敛性以收敛点作为解，再进行图像重建，得到最终的CT重建图像。

2.根据权利要求1所述的基于广义惩罚加权最小二乘的X射线CT图像重建方法，其特征在于，步骤(1)中获取的CT设备的系统参数包括X射线入射光子强度I₀和系统电子噪声的方差

3.根据权利要求2所述的基于广义惩罚加权最小二乘的X射线CT图像重建方法，其特征在于，步骤(2)中的对y进行逐个数据点上的方差估计具体是采用基于小邻域图像的局部方差估计或者基于CT投影数据噪声特性的方差估计。

4.根据权利要求1至3任意一项所述的基于广义惩罚加权最小二乘的X射线CT图像重建方法，其特征在于：

所述步骤(3)中建立的低剂量CT重建模型为:

$\underset{q\≥0}{min}{(y-q)}^T{G}^{-1}(y-q)+{\mathrm{\β}}_{2}R\left(q\right)......\left(I\right);$

其中，“T”表示转置运算，I为单位矩阵，β₁、β₂为两个超参数，β₁>0，β₂>0,R(q)是惩罚项，R(q)的表达式为:

$R\left(q\right)=\frac{1}{2}\underset{i}{\Σ}\underset{m\∈N_i}{\Σ}{\mathrm{\ω}}_{im}{(q_i-q_m)}^2......\left(II\right);$

其中，N_i为投影数据中第i个像素的一阶邻域，q_i为投影数据中第i个像素值，q_m为投影数据中第i个像素的邻域的内像素值，ω_im表示q_i的一阶邻域内的q_m所占的权重，并且ω_im在水平方向取值为1，在垂直方向取值为0.25。

5.根据权利要求4所述的基于广义惩罚加权最小二乘的X射线CT图像重建方法，其特征在于，所述步骤(4)具体是通过交替方向优化算法对(3)中的低剂量CT重建模型进行求解。

6.根据权利要求5所述的基于广义惩罚加权最小二乘的X射线CT图像重建方法，其特征在于，所述步骤(4)具体是通过交替方向优化算法对(3)中的低剂量CT重建模型进行求解，具体是：

引入变量p，得到

$\underset{p}{min}{(y-p)}^T{\mathrm{\Σ}}^{-1}(y-p)+{\mathrm{\β}}_1\left|\right|p-q|{|}_2^2={(y-q)}^T{G}^{-1}(y-q)......\left(IV\right);$

其中，

根据式(Ⅳ)，将式(Ⅰ)写成如下形式：

$\underset{p,q}{min}J(p,q)={(y-p)}^T{\mathrm{\Σ}}^{-1}(y-p)+{\mathrm{\β}}_1\left|\right|p-q|{|}_2^2+{\mathrm{\β}}_{2}R\left(q\right)......\left(V\right);$

其中，J(p,q)是J是关于p和q的函数；

(4.1)给定初始点q⁰＝y，令k＝0，k是迭代次数,k+1表示第k+1次迭代；

(4.2)使用交替方向法依次求p^k+1和q^k+1，

$\left\{\begin{array}{cc}{S}_G\left(q^k\right)=p^{k+1}=\underset{p}{\mathrm{argmin}}{(y-p)}^T{\mathrm{\Σ}}^{-1}(y-p)+{\mathrm{\β}}_1\left|\right|p-q^k|{|}_2^2& ......\left(a\right);\\ S_R\left(p^{k+1}\right)=q^{k+1}=\underset{q}{\mathrm{argmin}}{\mathrm{\β}}_1\left|\right|p^{k+1}-q|{|}_2^2+{\mathrm{\β}}_{2}R\left(q\right)& ......\left(b\right);\end{array}\right.......\left(VI\right);$

其中，S_G(q^k)表示式(Ⅴ)中对p的优化解，S_R(p^k+1)表示式(Ⅴ)中q的优化解；

(4.3)判断是否满足迭代终止条件：如果满足迭代终止条件，则停止迭代，进入步骤(4.5)，否则进入步骤(4.4)；

(4.4)令k＝k+1，进入步骤(4.2)；

(4.5)式Ⅵ(a)中的目标函数为光滑的二次凸函数，对Ⅵ(a)式求导并令其为零，得到Ⅵ(a)的解为

对于式Ⅵ(b)，使用Gauss-Seidel算法求解，其对应的迭代格式为其中 表示像素q_i左上方的领域，表示像素q_i右上方的领域，N_i表示像素的一阶领域。

7.根据权利要求6所述的基于广义惩罚加权最小二乘的X射线CT图像重建方法，其特征在于，所述步骤(4)中，具体通过滤波反投影方法进行图像重建，得到最终的CT重建图像。