CN105118078A - 欠采样的ct图像重建方法 - Google Patents

Abstract

本案涉及欠采样的CT图像重建方法，它首先通过已经训练构建好的自适应参数选择模型确定迭代算法目标函数中正则化参数的值，然后建立基于字典学习理论的迭代重建模型，在目标函数的正则约束项上添加加权惩罚因子，形成加权字典学习迭代重建模型，该模型的每次迭代过程交替更新加权惩罚因子、自适应字典、图像块的稀疏表示以及重建图像，最终达到迭代停止条件。本案能够很好的适应投影数据不充分的欠采样条件以满足降低辐射剂量的需求，在数据采样率下降到一定程度时依然能够高质量重建CT图像，对低剂量引起的投影噪声也有很好的鲁棒性。

Description

欠采样的CT图像重建方法

技术领域

本发明涉及数字图像处理技术领域，特别涉及一种对欠采样CT扫描获得的原始投影数据进行数据处理的图像重建方法。

背景技术

现有技术中具有代表性的CT图像重建算法是一种基于滤波反投影的解析重建方法。此算法以中心切片定理为依据，即二维图像在某一方向上的投影数据的一维傅里叶变换函数的值与该图像的二维傅里叶变换函数在频域平面上沿同一方向过原点的直线上的值相等。因此，滤波反投影算法通过各个方向上的投影数据得到一系列一维傅里叶变换函数，将这些函数按照对应方向加权叠加并进行滤波变换后得到待重建CT图像的二维傅里叶变换函数，对其做二维傅里叶反变换得到重建结果。

为达到精确重建的目的，上述解析重建方法的采样频率必须大于信号最高频率的2倍，(奈奎斯特抽样定理表明若频带宽度有限，要从抽样信号中无失真地恢复原信号，抽样频率应大于2倍信号最高频率。当抽样频率小于2倍频谱最高频率时，信号的频谱有混叠，解析重建方法从CT图像频率域还原到空间域时，会引入因为频谱混叠产生的伪影，降低重建图像的质量。)这使得投影数据的信息量远远大于重建图像本身，采样数据冗余度过高，高采样率的要求与降低CT辐射剂量的需求相悖。而如果降低CT扫描的采样率，当不满足奈奎斯特采样定律的条件时，重建图像会存在比较明显的伪影，影响CT图像的质量；另一方面，解析重建方法对投影数据的噪声也很敏感。

发明内容

针对现有技术存在的不足，本发明的目的在于通过构建一种全新的数据处理方法，使得可以通过欠采样CT扫描获得的原始数据重建出高质量的CT图像。

解析重建和迭代重建是CT图像重建的两种基本方法。与解析重建算法相比，迭代重建算法具有可在数据不完全和低信噪比(低剂量引起的噪声影响)条件下高质量重建CT图像的优点。针对滤波反投影算法无法应对欠采样投影数据的缺点，本案提出一种结合加权惩罚因子的基于字典学习理论的迭代重建算法。该算法首先通过已经训练构建好的自适应参数选择模型确定迭代算法目标函数中正则化参数的值，然后建立基于字典学习理论的迭代重建模型，在目标函数的正则约束项上添加加权惩罚因子，形成加权字典学习迭代重建模型，该模型的每次迭代过程交替更新加权惩罚因子、自适应字典、图像块的稀疏表示以及重建图像，最终达到迭代停止条件。

为实现上述目的，本发明通过以下技术方案实现：

一种欠采样的CT图像重建方法，包括以下步骤：

步骤1)加载待处理的原始CT扫描数据；

步骤2)基于字典学习的迭代重建算法，假设正则化参数λ为正无穷大，得到此条件下的初步重建图像

步骤3)计算出所述初步重建图像的前向投影数据与所述原始CT扫描数据的相对误差δ_λ→∞；

步骤4)构建一个正则化参数λ的自适应选择模型，该自适应选择模型通过以所述相对误差δ_λ→∞为自变量来获取与所述原始CT扫描数据相匹配的正则化参数λ的值；

步骤5)构建一个基于图像块全变分的加权惩罚因子w_s；

步骤6)将所述加权惩罚因子w_s与字典学习的约束项相结合，构建一个基于加权字典学习迭代重建模型的目标函数一；

步骤7)对所述目标函数一最小化求解，得到图像μ；

步骤8)根据图像μ更新加权惩罚因子w_s；

步骤9)重复执行步骤7)和步骤8)，直至满足迭代终止条件，输出最终的CT重建图像μ。

优选的是，所述的欠采样的CT图像重建方法，其中，步骤2)中的初步重建图像通过以下方法获得：

步骤a)：构建一个基于字典学习迭代重建模型的目标函数二，该目标函数二的表达式为：

$\underset{\mathrm{\μ},\mathrm{\α},D}{min}\underset{i=1}{\overset{I}{\Σ}}\frac{{\mathrm{\ω}}_i}{2}{({\[A\mathrm{\μ}\]}_i-{\hat{l}}_i)}^2+\mathrm{\λ}\left(\underset{s}{\Σ}\right||E_s\mathrm{\μ}-{\mathrm{D\α}}_s|{|}_2^2+\underset{s}{\Σ}{\mathrm{\ν}}_s\left|\right|{\mathrm{\α}}_s\left|{|}_0\right)$

其中，其为CT成像的系统矩阵；其为原始投影数据；其为统计权重；其为待求图像排成的列向量；为数据保真项；其为从图像中提取N₀×N₀维小图像块的算子；其为字典；α_s∈R^K×1，其是以字典D为基底得到的图像块的稀疏表示；为正则约束项；λ和ν_s均为正则化参数；

步骤b)：在字典更新过程中，图像μ的初值取随机矩阵，之后图像μ取上一次图像更新过程产生的更新结果，并假设图像μ在本次字典更新过程中保持恒定，则目标函数二的表达式简化为简式一，具体为：

$\underset{\left(D\right),\mathrm{\α}}{min}\underset{s}{\Σ}\left|\right|E_s\mathrm{\μ}-{\mathrm{D\α}}_s|{|}_2^2+\underset{x}{\Σ}{\mathrm{\ν}}_s|\left|{\mathrm{\α}}_s\right|{|}_0$

其中，字典D的更新采用K-SVD算法获得，稀疏表示α_s的更新采用OMP算法获得；

步骤c)：在图像更新过程中，字典D的初值取离散余弦变换矩阵，稀疏表示α_s的初值以图像μ和字典D的初值为基础通过OMP算法获得，且之后字典D和稀疏表示α_s取上一次字典更新过程产生的更新结果，则目标函数二的表达式简化为简式二，具体为：

$\underset{\mathrm{\μ}}{min}\underset{i=1}{\overset{I}{\Σ}}\frac{{\mathrm{\ω}}_i}{2}{({\[A\mathrm{\μ}\]}_i-{\hat{l}}_i)}^2+\mathrm{\λ}\underset{s}{\Σ}\left|\right|E_s\mathrm{\μ}-{\mathrm{D\α}}_s|{|}_2^2$

以二次可分离代理函数来替代简式二中的数据保真项则图像μ的最终表达式为：

${\mathrm{\μ}}_j^{(t+1)}={\[{\mathrm{\μ}}_j^{\left(t\right)}-\frac{\underset{i=1}{\overset{I}{\Σ}}\left(a_{ij}{\mathrm{\ω}}_i\right({\[{\mathrm{A\μ}}^{\left(t\right)}\]}_i-{\hat{l}}_i\left)\right)+2\mathrm{\λ}\underset{s}{\Σ}\underset{n=1}{\overset{N_0^2}{\Σ}}{e}_{nj}^s({\[E_s{\mathrm{\μ}}^{\left(t\right)}\]}_n-{\[{\mathrm{D\α}}_s\]}_n)}{\underset{i=1}{\overset{I}{\Σ}}\left(a_{ij}{\mathrm{\ω}}_i\underset{k=1}{\overset{N^2}{\Σ}}{a}_{ik}\right)+2\mathrm{\λ}\underset{s}{\Σ}\underset{n=1}{\overset{N_0^2}{\Σ}}{e}_{nj}^s\underset{k=1}{\overset{N^2}{\Σ}}{e}_{nk}^s}\]}_{+}$

j＝1,2,…,N²

其中，t为迭代次数，j表示图像μ的第j个元素，二次可分离代理函数为：

$Q(\mathrm{\μ};{\mathrm{\μ}}^{(t+1)})=\underset{i=1}{\overset{I}{\Σ}}\underset{j=1}{\overset{N\×N}{\Σ}}{\mathrm{\α}}_{ij}\frac{{\mathrm{\ω}}_i}{2}\×{\left(\frac{a_{ij}}{{\mathrm{\α}}_{ij}}\right({\mathrm{\μ}}_j-{\mathrm{\μ}}_j^{\left(t\right)})+{\[{\mathrm{A\μ}}^{\left(t\right)}\]}_i-{\hat{l}}_i)}^2$

其中， ${\mathrm{\α}}_{ij}=\frac{a_{ij}}{\underset{j=1}{\overset{N\×N}{\Σ}}{a}_{ij}};$

步骤d)：在正则化参数λ为正无穷大的条件下，重复步骤b)和c)，直至满足迭代终止条件，获得初步重建图像

优选的是，所述的欠采样的CT图像重建方法，其中，步骤3)中，相对误差δ_λ→∞的表达式为：

${\mathrm{\δ}}_{\mathrm{\λ}\→\mathrm{\∞}}=\frac{\underset{i=1}{\overset{I}{\Σ}}\frac{{\mathrm{\ω}}_i}{2}{({\[{\mathrm{A\μ}}_{\mathrm{\λ}\→\mathrm{\∞}}^{*}\]}_i-l_i)}^2}{\underset{i=1}{\overset{I}{\Σ}}\frac{{\mathrm{\ω}}_i}{2}{l}_i^2}.$

优选的是，所述的欠采样的CT图像重建方法，其中，步骤4)中，所述正则化参数λ的自适应选择模型的表达式为：

$\mathrm{\λ}=f\left({\mathrm{\δ}}_G\right)\frac{\underset{j}{\Σ}\underset{i=1}{\overset{I}{\Σ}}\left(a_{ij}{\mathrm{\ω}}_i\underset{k=1}{\overset{N^2}{\Σ}}{a}_{ik}\right)}{\underset{j}{\Σ}\underset{s}{\Σ}\underset{n=1}{\overset{N_0^2}{\Σ}}{e}_{nj}^s\underset{k=1}{\overset{N^2}{\Σ}}{e}_{nk}^s}$

其中，f(δ_G)的表达式为：

$f\left({\mathrm{\δ}}_G\right)=\left\{\begin{array}{cc}1.74485{\mathrm{\δ}}_G^2+0.58883{\mathrm{\δ}}_G-6.88253& if{\mathrm{\δ}}_G>1.96\\ -0.21545{\mathrm{\δ}}_G^2+1.08602{\mathrm{\δ}}_G-0.32634& if{\mathrm{\δ}}_G\≤1.96\end{array}\right.$

其中，δ_G＝10⁶δ_λ→∞。

优选的是，所述的欠采样的CT图像重建方法，其中，步骤5)中，所述加权惩罚因子w_s的表达式为：

$w_s=\frac{mean\left(\right|\left|{\mathrm{\μ}}_s\right|{|}_{TV})}{\left|\right|{\mathrm{\μ}}_s|{|}_{TV}+\mathrm{\ϵ}}$

其中，μ_s表示从图像μ中抽取出的第s个图像块，且μ_s＝reshape(E_sμ,N₀,N₀)，其表示从图像μ中提取N₀×N₀维图像块的算子；reshape函数表示将E_sμ重新排列为二维图像块；||μ_s||_TV为图像块的全变分，计算公式为G为梯度算子，(Gx)_i,j＝(x_i+1,j-x_i,j,x_i,j+1-x_i,j)，x代表任意矩阵；mean(||μ_s||_TV)表示所有图像块全变分的平均值；ε＞0。

优选的是，所述的欠采样的CT图像重建方法，其中，步骤6)中，所述目标函数一的表达式为：

$\underset{\mathrm{\μ},\mathrm{\α},D}{min}\underset{i=1}{\overset{I}{\Σ}}\frac{{\mathrm{\ω}}_i}{2}{({\[A\mathrm{\μ}\]}_i-{\hat{l}}_i)}^2+\mathrm{\λ}\underset{s}{\Σ}{w}_s\left(\right||E_s\mathrm{\μ}-{\mathrm{D\α}}_s|{|}_2^2+{\mathrm{\ν}}_s\left|\right|{\mathrm{\α}}_s\left|{|}_0\right).$

优选的是，所述的欠采样的CT图像重建方法，其中，步骤7)中，采用与获得初步重建图像相同的方法对目标函数一进行最小化求解，得到的图像μ的表达式为：

${\mathrm{\μ}}_j^{(t+1)}={\[{\mathrm{\μ}}_j^t-\frac{\underset{i=1}{\overset{I}{\Σ}}\left(a_{ij}{\mathrm{\ω}}_i\right({\[{\mathrm{A\μ}}^{\left(t\right)}\]}_i-\hat{l}\left)\right)+2\mathrm{\λ}\underset{s}{\Σ}{w}_s\underset{n=1}{\overset{N_0^2}{\Σ}}{e}_{nj}^s({\[E_s{\mathrm{\μ}}^{\left(t\right)}\]}_n-{\[{\mathrm{D\α}}_s\]}_n)}{\underset{i=1}{\overset{I}{\Σ}}\left(a_{ij}{\mathrm{\ω}}_i\underset{k=1}{\overset{N^2}{\Σ}}{a}_{ik}\right)+2\mathrm{\λ}\underset{s}{\Σ}{w}_s\underset{n=1}{\overset{N_0^2}{\Σ}}{e}_{nj}^s\underset{k=1}{\overset{N^2}{\Σ}}{e}_{nk}^s}\]}_{+}$

其中，加权惩罚因子w_s的初始值取1。

优选的是，所述的欠采样的CT图像重建方法，其中，步骤8)中，所述加权惩罚因子w_s的更新公式为：

${\mathrm{\μ}}_s^{\left(t\right)}=reshape(E_s{\mathrm{\μ}}^{\left(t\right)},N_0,N_0)$

$w_s^{(t+1)}=\frac{mean\left(\right|\left|{\mathrm{\μ}}_s^{\left(t\right)}\right|{|}_{TV})}{\left|\right|{\mathrm{\μ}}_s^{\left(t\right)}|{|}_{TV}+\mathrm{\ϵ}}.$

优选的是，所述的欠采样的CT图像重建方法，其中，在步骤9)中，迭代终止条件为两次迭代的扫描数据相对误差增量小于给定阈值T₁，迭代终止条件的表达式为：

$\frac{\left|\sqrt{\underset{i=1}{\overset{I}{\&Sigma;}}\frac{{\mathrm{\&omega;}}_i}{2}{({\&lsqb;{\mathrm{A\&mu;}}^{(t+1)}\&rsqb;}_i-{\hat{l}}_i)}^2}-\sqrt{\underset{i=1}{\overset{I}{\&Sigma;}}\frac{{\mathrm{\&omega;}}_i}{2}{({\&lsqb;{\mathrm{A\&mu;}}^{\left(t\right)}\&rsqb;}_i-{\hat{l}}_i)}^2}\right|}{\sqrt{\underset{i=1}{\overset{I}{\&Sigma;}}\frac{{\mathrm{\&omega;}}_i}{2}{({\&lsqb;{\mathrm{A\&mu;}}^{(t+1)}\&rsqb;}_i-{\hat{l}}_i)}^2}}<T_1.$

优选的是，所述的欠采样的CT图像重建方法，其中，在步骤d)中，迭代终止条件为两次迭代的扫描数据相对误差增量小于给定阈值T₂，迭代终止条件的表达式为：

$\frac{\left|\sqrt{\underset{i=1}{\overset{I}{\&Sigma;}}\frac{{\mathrm{\&omega;}}_i}{2}{({\&lsqb;{\mathrm{A\&mu;}}^{(t+1)}\&rsqb;}_i-{\hat{l}}_i)}^2}-\sqrt{\underset{i=1}{\overset{I}{\&Sigma;}}\frac{{\mathrm{\&omega;}}_i}{2}{({\&lsqb;{\mathrm{A\&mu;}}^{\left(t\right)}\&rsqb;}_i-{\hat{l}}_i)}^2}\right|}{\sqrt{\underset{i=1}{\overset{I}{\&Sigma;}}\frac{{\mathrm{\&omega;}}_i}{2}{({\&lsqb;{\mathrm{A\&mu;}}^{(t+1)}\&rsqb;}_i-{\hat{l}}_i)}^2}}<T_2.$

本发明的有益效果是：本案能够很好的适应投影数据不充分的欠采样条件以满足降低辐射剂量的需求，在数据采样率下降到一定程度时依然能够高质量重建CT图像，对低剂量引起的投影噪声也有很好的鲁棒性。

附图说明

图1为欠采样的CT图像重建方法的流程示意图。

图2为采用不同方法获得的实物CT图像的比对图。

具体实施方式

下面结合附图对本发明做进一步的详细说明，以令本领域技术人员参照说明书文字能够据以实施。

本案提供一实施例的欠采样的CT图像重建方法，包括以下步骤：

步骤1)加载待处理的原始CT扫描数据；

步骤2)基于字典学习的迭代重建算法，假设正则化参数λ为正无穷大，得到此条件下的初步重建图像

步骤3)计算出初步重建图像的前向投影数据与原始CT扫描数据的相对误差δ_λ→∞；

步骤4)构建一个正则化参数λ的自适应选择模型，该自适应选择模型通过以相对误差δ_λ→∞为自变量来获取与原始CT扫描数据相匹配的正则化参数λ的值；

步骤5)构建一个基于图像块全变分的加权惩罚因子w_s；

步骤6)将该加权惩罚因子w_s与字典学习的约束项相结合，构建一个基于加权字典学习迭代重建模型的目标函数一；

步骤7)对该目标函数一最小化求解，得到图像μ；

步骤8)根据图像μ更新加权惩罚因子w_s；

步骤9)重复执行步骤7)和步骤8)，直至满足迭代终止条件，输出最终的CT重建图像μ。

具体的图像重建方法如下：

步骤2)中的初步重建图像通过以下方法获得：

步骤a)：构建一个基于字典学习迭代重建模型的目标函数二，该目标函数二的表达式为：

$\underset{\mathrm{\&mu;},\mathrm{\&alpha;},D}{min}\underset{i=1}{\overset{I}{\&Sigma;}}\frac{{\mathrm{\&omega;}}_i}{2}{({\&lsqb;A\mathrm{\&mu;}\&rsqb;}_i-{\hat{l}}_i)}^2+\mathrm{\&lambda;}\left(\underset{s}{\&Sigma;}\right||E_s\mathrm{\&mu;}-{\mathrm{D\&alpha;}}_s|{|}_2^2+\underset{s}{\&Sigma;}{\mathrm{\&nu;}}_s\left|\right|{\mathrm{\&alpha;}}_s\left|{|}_0\right)$

其中，其为CT成像的系统矩阵，a_ij表示系统矩阵A中第_i行第j列的元素，代表大小为I行N²列的矩阵集合，I为原始投影的总数目，N²表示图像块的大小；其为排成列向量的原始投影数据，R^I×1代表长度为I的列向量；其为统计权重；其为待求图像排成的列向量，T代表转置；为数据保真项；其为从图像中提取N₀×N₀维小图像块的算子，s代表图像块的编号，即第s个小图像块，nj代表抽取的小图像块以图像中第n行第j列的元素作为小图像块的左上角元素，nj遍历整个图像；其为字典，d代表原子，与图像块的大小相同，K代表原子的数目；α_s∈R^K×1，其是以字典D为基底得到的图像块的稀疏表示；为正则约束项；λ和ν_s均为正则化参数；由于图像块的稀疏表示通过OMP算法(OMP，OrthogonalMatchingPursuit，正交匹配追踪)实现，无需确定ν_s，因此这里阐述的是求解λ这个正则化参数的模型。(需要申明的是，在本案公式中的出现的字母、字符、上标或下标所指代的含义在全文是通用的，因此，本案为了避免造成叙述冗长，未在每个公式后面重复标注每个字母或字符的含义。)

步骤b)：在字典更新过程中，图像μ的初值取随机矩阵，之后图像μ取上一次图像更新过程产生的更新结果，并假设图像μ在本次字典更新过程中保持恒定，则目标函数二的表达式简化为简式一，具体为：

$\underset{\left(D\right),\mathrm{\&alpha;}}{min}\underset{s}{\&Sigma;}\left|\right|E_s\mathrm{\&mu;}-{\mathrm{D\&alpha;}}_s|{|}_2^2+\underset{x}{\&Sigma;}{\mathrm{\&nu;}}_s|\left|{\mathrm{\&alpha;}}_s\right|{|}_0$

其中，字典D的更新采用K-SVD算法获得，稀疏表示α_s的更新采用OMP算法获得；

步骤c)：在图像更新过程中，字典D的初值取离散余弦变换(DCT)矩阵，稀疏表示α_s的初值以图像μ和字典D的初值为基础通过OMP算法获得，且之后字典D和稀疏表示α_s取上一次字典更新过程产生的更新结果，则目标函数二的表达式简化为简式二，具体为：

$\underset{\mathrm{\&mu;}}{min}\underset{i=1}{\overset{I}{\&Sigma;}}\frac{{\mathrm{\&omega;}}_i}{2}{({\&lsqb;A\mathrm{\&mu;}\&rsqb;}_i-{\hat{l}}_i)}^2+\mathrm{\&lambda;}\underset{s}{\&Sigma;}\left|\right|E_s\mathrm{\&mu;}-{\mathrm{D\&alpha;}}_s|{|}_2^2$

以二次可分离代理函数来替代简式二中的数据保真项则图像μ的最终表达式为：

${\mathrm{\&mu;}}_j^{(t+1)}={\&lsqb;{\mathrm{\&mu;}}_j^{\left(t\right)}-\frac{\underset{i=1}{\overset{I}{\&Sigma;}}\left(a_{ij}{\mathrm{\&omega;}}_i\right({\&lsqb;{\mathrm{A\&mu;}}^{\left(t\right)}\&rsqb;}_i-{\hat{l}}_i\left)\right)+2\mathrm{\&lambda;}\underset{s}{\&Sigma;}\underset{n=1}{\overset{N_0^2}{\&Sigma;}}{e}_{nj}^s({\&lsqb;E_s{\mathrm{\&mu;}}^{\left(t\right)}\&rsqb;}_n-{\&lsqb;{\mathrm{D\&alpha;}}_s\&rsqb;}_n)}{\underset{i=1}{\overset{I}{\&Sigma;}}\left(a_{ij}{\mathrm{\&omega;}}_i\underset{k=1}{\overset{N^2}{\&Sigma;}}{a}_{ik}\right)+2\mathrm{\&lambda;}\underset{s}{\&Sigma;}\underset{n=1}{\overset{N_0^2}{\&Sigma;}}{e}_{nj}^s\underset{k=1}{\overset{N^2}{\&Sigma;}}{e}_{nk}^s}\&rsqb;}_{+}$

j＝1,2,…,N²

其中，t为迭代次数，j表示图像μ的第j个元素，二次可分离代理函数为：

$Q(\mathrm{\&mu;};{\mathrm{\&mu;}}^{(t+1)})=\underset{i=1}{\overset{I}{\&Sigma;}}\underset{j=1}{\overset{N\&times;N}{\&Sigma;}}{\mathrm{\&alpha;}}_{ij}\frac{{\mathrm{\&omega;}}_i}{2}\&times;{\left(\frac{a_{ij}}{{\mathrm{\&alpha;}}_{ij}}\right({\mathrm{\&mu;}}_j-{\mathrm{\&mu;}}_j^{\left(t\right)})+{\&lsqb;{\mathrm{A\&mu;}}^{\left(t\right)}\&rsqb;}_i-{\hat{l}}_i)}^2$

其中， ${\mathrm{\&alpha;}}_{ij}=\frac{a_{ij}}{\underset{j=1}{\overset{N\&times;N}{\&Sigma;}}{a}_{ij}};$

步骤d)：在正则化参数λ为正无穷大的条件下，重复步骤b)和c)，直至满足迭代终止条件，获得初步重建图像

步骤3)中，相对误差δ_λ→∞的表达式为：

${\mathrm{\&delta;}}_{\mathrm{\&lambda;}\&RightArrow;\mathrm{\&infin;}}=\frac{\underset{i=1}{\overset{I}{\&Sigma;}}\frac{{\mathrm{\&omega;}}_i}{2}{({\&lsqb;{\mathrm{A\&mu;}}_{\mathrm{\&lambda;}\&RightArrow;\mathrm{\&infin;}}^{*}\&rsqb;}_i-l_i)}^2}{\underset{i=1}{\overset{I}{\&Sigma;}}\frac{{\mathrm{\&omega;}}_i}{2}{l}_i^2}.$

步骤4)中，正则化参数λ的自适应选择模型的表达式为：

$\mathrm{\&lambda;}=f\left({\mathrm{\&delta;}}_G\right)\frac{\underset{j}{\&Sigma;}\underset{i=1}{\overset{I}{\&Sigma;}}\left(a_{ij}{\mathrm{\&omega;}}_i\underset{k=1}{\overset{N^2}{\&Sigma;}}{a}_{ik}\right)}{\underset{j}{\&Sigma;}\underset{s}{\&Sigma;}\underset{n=1}{\overset{N_0^2}{\&Sigma;}}{e}_{nj}^s\underset{k=1}{\overset{N^2}{\&Sigma;}}{e}_{nk}^s}$

其中，f(δ_G)的表达式为：

$f\left({\mathrm{\&delta;}}_G\right)=\left\{\begin{array}{cc}1.74485{\mathrm{\&delta;}}_G^2+0.58883{\mathrm{\&delta;}}_G-6.88253& if{\mathrm{\&delta;}}_G>1.96\\ -0.21545{\mathrm{\&delta;}}_G^2+1.08602{\mathrm{\&delta;}}_G-0.32634& if{\mathrm{\&delta;}}_G\&le;1.96\end{array}\right.$

其中，δ_G＝10⁶δ_λ→∞。

f(δ_G)的获得主要是通过一系列训练样本图像以及对应的经验选取正则化参数λ，随后通过函数拟合的方法求出f(δ_G)，经过分析，该拟合函数还与图像大小有关，当图像大小为256×256维时，用分段二次函数拟合效果最好。

正则化参数λ的传统选取方法是经验选取，需要经过大量的测试来确定合适的正则化参数λ的取值，这需要耗费大量的测试时间，这里提出的正则化参数λ的自适应选取模型节省了这部分测试时间，提高了运算效率。

步骤5)中，加权惩罚因子w_s的表达式为：

$w_s=\frac{mean\left(\right|\left|{\mathrm{\&mu;}}_s\right|{|}_{TV})}{\left|\right|{\mathrm{\&mu;}}_s|{|}_{TV}+\mathrm{\&epsiv;}}$

其中，μ_s表示从图像μ中抽取出的第s个图像块，且μ_s＝reshape(E_sμ,N₀,N₀)，其表示从图像μ中提取N₀×N₀维图像块的算子；reshape函数表示将E_sμ重新排列为二维图像块；||μ_s||_TV为图像块的全变分，计算公式为G为梯度算子，(Gx)_i,j＝(x_i+1,j-x_i,j,x_i,j+1-x_i,j)，x代表任意矩阵；mean(||μ_s||_TV)表示所有图像块全变分的平均值；ε＞0。ε的作用是保证算法的稳定性。ε一般可取0.00001。

相比于未结合加权惩罚因子的算法，加权字典学习迭代重建算法能够保留重建图像更多的细节信息，因为其设计原理是针对细节信息较多的图像块施加较小的加权惩罚因子，而图像块的全变分越大，则包含的细节信息越多，因此设计的加权惩罚因子的表达式与图像块的全变分成反比。加权惩罚因子更好的区分了图像的平滑区域和细节区域，在算法迭代过程中给予不同程度的平滑效果，加权惩罚因子越小，平滑效果越弱，从而在滤除噪声的同时，更好的保留了图像的细节信息。

步骤6)中，目标函数一的表达式为：

$\underset{\mathrm{\&mu;},\mathrm{\&alpha;},D}{min}\underset{i=1}{\overset{I}{\&Sigma;}}\frac{{\mathrm{\&omega;}}_i}{2}{({\&lsqb;A\mathrm{\&mu;}\&rsqb;}_i-{\hat{l}}_i)}^2+\mathrm{\&lambda;}\underset{s}{\&Sigma;}{w}_s\left(\right||E_s\mathrm{\&mu;}-{\mathrm{D\&alpha;}}_s|{|}_2^2+{\mathrm{\&nu;}}_s\left|\right|{\mathrm{\&alpha;}}_s\left|{|}_0\right).$

其中λ的取值已在正则化参数选取步骤中确定，加权字典学习重建模型在原算法基础上加入了基于图像块全变分的加权惩罚因子，其迭代求解过程依然分为两个循环的步骤：字典更新和图像更新，并采用交替最小化策略达到迭代截止条件。

步骤7)中，采用与获得初步重建图像相同的方法对目标函数一进行最小化求解，得到的图像μ的表达式为：

${\mathrm{\&mu;}}_j^{(t+1)}={\&lsqb;{\mathrm{\&mu;}}_j^{\left(t\right)}-\frac{\underset{i=1}{\overset{I}{\&Sigma;}}\left(a_{ij}{\mathrm{\&omega;}}_i\right({\&lsqb;{\mathrm{A\&mu;}}^{\left(t\right)}\&rsqb;}_i-{\hat{l}}_i\left)\right)+2\mathrm{\&lambda;}\underset{s}{\&Sigma;}\underset{n=1}{\overset{N_0^2}{\&Sigma;}}{e}_{nj}^s({\&lsqb;E_s{\mathrm{\&mu;}}^{\left(t\right)}\&rsqb;}_n-{\&lsqb;{\mathrm{D\&alpha;}}_s\&rsqb;}_n)}{\underset{i=1}{\overset{I}{\&Sigma;}}\left(a_{ij}{\mathrm{\&omega;}}_i\underset{k=1}{\overset{N^2}{\&Sigma;}}{a}_{ik}\right)+2\mathrm{\&lambda;}\underset{s}{\&Sigma;}\underset{n=1}{\overset{N_0^2}{\&Sigma;}}{e}_{nj}^s\underset{k=1}{\overset{N^2}{\&Sigma;}}{e}_{nk}^s}\&rsqb;}_{+}$

其中，初始化各迭代参数：μ⁽⁰⁾为随机矩阵，D为DCT矩阵，这两个参数的取值与第一次迭代时相同，加权惩罚因子w_s＝1。

步骤8)中，加权惩罚因子w_s的更新公式为：

${\mathrm{\&mu;}}_s^{\left(t\right)}=reshape(E_s{\mathrm{\&mu;}}^{\left(t\right)},N_0,N_0)$

$w_s^{(t+1)}=\frac{mean\left(\right|\left|{\mathrm{\&mu;}}_s^{\left(t\right)}\right|{|}_{TV})}{\left|\right|{\mathrm{\&mu;}}_s^{\left(t\right)}|{|}_{TV}+\mathrm{\&epsiv;}}.$

在步骤9)中，迭代终止条件为两次迭代的扫描数据相对误差增量小于给定阈值T₁，迭代终止条件的表达式为：

$\frac{\left|\sqrt{\underset{i=1}{\overset{I}{\&Sigma;}}\frac{{\mathrm{\&omega;}}_i}{2}{({\&lsqb;{\mathrm{A\&mu;}}^{(t+1)}\&rsqb;}_i-{\hat{l}}_i)}^2}-\sqrt{\underset{i=1}{\overset{I}{\&Sigma;}}\frac{{\mathrm{\&omega;}}_i}{2}{({\&lsqb;{\mathrm{A\&mu;}}^{\left(t\right)}\&rsqb;}_i-{\hat{l}}_i)}^2}\right|}{\sqrt{\underset{i=1}{\overset{I}{\&Sigma;}}\frac{{\mathrm{\&omega;}}_i}{2}{({\&lsqb;{\mathrm{A\&mu;}}^{(t+1)}\&rsqb;}_i-{\hat{l}}_i)}^2}}<T_1.$

在步骤d)中，迭代终止条件为两次迭代的扫描数据相对误差增量小于给定阈值T₂，迭代终止条件的表达式为：

$\frac{\left|\sqrt{\underset{i=1}{\overset{I}{\&Sigma;}}\frac{{\mathrm{\&omega;}}_i}{2}{({\&lsqb;{\mathrm{A\&mu;}}^{(t+1)}\&rsqb;}_i-{\hat{l}}_i)}^2}-\sqrt{\underset{i=1}{\overset{I}{\&Sigma;}}\frac{{\mathrm{\&omega;}}_i}{2}{({\&lsqb;{\mathrm{A\&mu;}}^{\left(t\right)}\&rsqb;}_i-{\hat{l}}_i)}^2}\right|}{\sqrt{\underset{i=1}{\overset{I}{\&Sigma;}}\frac{{\mathrm{\&omega;}}_i}{2}{({\&lsqb;{\mathrm{A\&mu;}}^{(t+1)}\&rsqb;}_i-{\hat{l}}_i)}^2}}<T_2.$

两次迭代终止条件的阈值T是根据实际情况确定的，两次的阈值T₁和T₂可以相同，也可以不相同。例如，一般地，阈值T可取0.001～0.00001。

图2是采用不同方法获得的实物CT图像的比对图，其中，A为基准图像，它是背景技术里提到的采用滤波反投影的解析重建方法得到的，是目前最成熟的重建方法，此图像从采样次数足够多(即满足奈奎斯特抽样定理)的原始扫描数据重建得到；B和C均为欠采样条件下得到的原始扫面数据重建得到的结果，B为采用本案的图像重建方法获得的结果图像；C为采用原字典学习重建算法(即未引入加权惩罚因子)获得的结果图像；D为B与A的差值图像；E为C与A的差值图像。从图2可以直观地看到，相较于C，B与A的图像更接近，图像噪点少，画面更平滑清晰，而从D与E的对比也能更好地说明这一点。

尽管本发明的实施方案已公开如上，但其并不仅仅限于说明书和实施方式中所列运用，它完全可以被适用于各种适合本发明的领域，对于熟悉本领域的人员而言，可容易地实现另外的修改，因此在不背离权利要求及等同范围所限定的一般概念下，本发明并不限于特定的细节和这里示出与描述的图例。

Claims (10)

1.一种欠采样的CT图像重建方法，其特征在于，包括以下步骤：

步骤1)加载待处理的原始CT扫描数据；

步骤2)基于字典学习的迭代重建算法，假设正则化参数λ为正无穷大，得到此条件下的初步重建图像

步骤3)计算出所述初步重建图像的前向投影数据与所述原始CT扫描数据的相对误差δ_λ→∞；

步骤4)构建一个正则化参数λ的自适应选择模型，该自适应选择模型通过以所述相对误差δ_λ→∞为自变量来获取与所述原始CT扫描数据相匹配的正则化参数λ的值；

步骤5)构建一个基于图像块全变分的加权惩罚因子w_s；

步骤6)将所述加权惩罚因子w_s与字典学习的约束项相结合，构建一个基于加权字典学习迭代重建模型的目标函数一；

步骤7)对所述目标函数一最小化求解，得到图像μ；

步骤8)根据图像μ更新加权惩罚因子w_s；

步骤9)重复执行步骤7)和步骤8)，直至满足迭代终止条件，输出最终的CT重建图像μ。

2.如权利要求1所述的欠采样的CT图像重建方法，其特征在于，步骤2)中的初步重建图像通过以下方法获得：

步骤a)：构建一个基于字典学习迭代重建模型的目标函数二，该目标函数二的表达式为：

$\underset{\mathrm{\&mu;},\mathrm{\&alpha;},D}{min}\underset{i=1}{\overset{I}{\&Sigma;}}\frac{{\mathrm{\&omega;}}_i}{2}{\left({\&lsqb;A\mathrm{\&mu;}\&rsqb;}_i-{\hat{l}}_i\right)}^2+\mathrm{\&lambda;}\left(\underset{s}{\&Sigma;}\right||E_s\mathrm{\&mu;}-{\mathrm{D\&alpha;}}_s|{|}_2^2+\underset{s}{\&Sigma;}{\mathrm{\&nu;}}_s\left|\right|{\mathrm{\&alpha;}}_s\left|{|}_0\right)$

其中，其为CT成像的系统矩阵；其为原始投影数据；其为统计权重；其为待求图像排成的列向量；为数据保真项；其为从图像中提取N₀×N₀维小图像块的算子；其为字典；α_s∈R^K×1，其是以字典D为基底得到的图像块的稀疏表示；为正则约束项；λ和ν_s均为正则化参数；

步骤b)：在字典更新过程中，图像μ的初值取随机矩阵，之后图像μ取上一次图像更新过程产生的更新结果，并假设图像μ在本次字典更新过程中保持恒定，则目标函数二的表达式简化为简式一，具体为：

$\underset{\left(D\right),\mathrm{\&alpha;}}{min}\underset{s}{\&Sigma;}\left|\right|E_s\mathrm{\&mu;}-{\mathrm{D\&alpha;}}_s|{|}_2^2+\underset{x}{\&Sigma;}{\mathrm{\&nu;}}_s|\left|{\mathrm{\&alpha;}}_s\right|{|}_0$

其中，字典D的更新采用K-SVD算法获得，稀疏表示α_s的更新采用OMP算法获得；

步骤c)：在图像更新过程中，字典D的初值取离散余弦变换矩阵，稀疏表示α_s的初值以图像μ和字典D的初值为基础通过OMP算法获得，且之后字典D和稀疏表示α_s取上一次字典更新过程产生的更新结果，则目标函数二的表达式简化为简式二，具体为：

$\underset{\mathrm{\&mu;}}{min}\underset{i=1}{\overset{I}{\&Sigma;}}\frac{{\mathrm{\&omega;}}_i}{2}{\left({\&lsqb;A\mathrm{\&mu;}\&rsqb;}_i-{\hat{l}}_i\right)}^2+\mathrm{\&lambda;}\underset{s}{\&Sigma;}\left|\right|E_s\mathrm{\&mu;}-{\mathrm{D\&alpha;}}_s|{|}_2^2$

以二次可分离代理函数来替代简式二中的数据保真项则图像μ的最终表达式为：

${\mathrm{\&mu;}}_j^{(t+1)}={\&lsqb;{\mathrm{\&mu;}}_j^{\left(t\right)}-\frac{\underset{i=1}{\overset{I}{\&Sigma;}}\left(a_{ij}{\mathrm{\&omega;}}_i\right({\&lsqb;{\mathrm{A\&mu;}}^{\left(t\right)}\&rsqb;}_i-{\hat{l}}_i\left)\right)+2\mathrm{\&lambda;}\underset{s}{\&Sigma;}\underset{n=1}{\overset{N_0^2}{\&Sigma;}}{e}_{nj}^s({\&lsqb;E_s{\mathrm{\&mu;}}^{\left(t\right)}\&rsqb;}_n-{\&lsqb;{\mathrm{D\&alpha;}}_s\&rsqb;}_n)}{\underset{i=1}{\overset{I}{\&Sigma;}}\left(a_{ij}{\mathrm{\&omega;}}_i\underset{k=1}{\overset{N^2}{\&Sigma;}}{a}_{ik}\right)+2\mathrm{\&lambda;}\underset{s}{\&Sigma;}\underset{n=1}{\overset{N_0^2}{\&Sigma;}}{e}_{nj}^s\underset{k=1}{\overset{N^2}{\&Sigma;}}{e}_{nk}^s}\&rsqb;}_{+}$

j＝1,2,...,N²

其中，t为迭代次数，j表示图像μ的第j个元素，二次可分离代理函数为：

$Q\left(\mathrm{\&mu;};{\mathrm{\&mu;}}^{\left(t+1\right)}\right)=\underset{i=1}{\overset{I}{\&Sigma;}}\underset{j=1}{\overset{N\&times;N}{\&Sigma;}}{\mathrm{\&alpha;}}_{ij}\frac{{\mathrm{\&omega;}}_i}{2}\&times;{\left(\frac{a_{ij}}{{\mathrm{\&alpha;}}_{ij}}\left({\mathrm{\&mu;}}_j-{\mathrm{\&mu;}}_j^{\left(t\right)}\right)+{\&lsqb;{\mathrm{A\&mu;}}^{\left(t\right)}\&rsqb;}_i-{\hat{l}}_i\right)}^2$

其中， ${\mathrm{\&alpha;}}_{ij}=\frac{a_{ij}}{\underset{j=1}{\overset{N\&times;N}{\&Sigma;}}{a}_{ij}};$

步骤d)：在正则化参数λ为正无穷大的条件下，重复步骤b)和c)，直至满足迭代终止条件，获得初步重建图像

3.根据权利要求2所述的欠采样的CT图像重建方法，其特征在于，步骤3)中，相对误差δ_λ→∞的表达式为：

${\mathrm{\&delta;}}_{\mathrm{\&lambda;}\&RightArrow;\mathrm{\&infin;}}=\frac{\underset{i=1}{\overset{I}{\&Sigma;}}\frac{{\mathrm{\&omega;}}_i}{2}{({\&lsqb;{\mathrm{A\&mu;}}_{\mathrm{\&lambda;}\&RightArrow;\mathrm{\&infin;}}^{*}\&rsqb;}_i-l_i)}^2}{\underset{i=1}{\overset{I}{\&Sigma;}}\frac{{\mathrm{\&omega;}}_i}{2}{l}_i^2}.$

4.根据权利要求3所述的欠采样的CT图像重建方法，其特征在于，步骤4)中，所述正则化参数λ的自适应选择模型的表达式为：

$\mathrm{\&lambda;}=f\left({\mathrm{\&delta;}}_G\right)\frac{\underset{j}{\&Sigma;}\underset{i=1}{\overset{I}{\&Sigma;}}\left(a_{ij}{\mathrm{\&omega;}}_i\underset{k=1}{\overset{N^2}{\&Sigma;}}{a}_{ik}\right)}{\underset{j}{\&Sigma;}\underset{s}{\&Sigma;}\underset{n=1}{\overset{N_0^2}{\&Sigma;}}{e}_{nj}^s\underset{k=1}{\overset{N^2}{\&Sigma;}}{e}_{nk}^s}$

其中，f(δ_G)的表达式为：

$f\left({\mathrm{\&delta;}}_G\right)=\left\{\begin{array}{ccc}1.74485{\mathrm{\&delta;}}_G^2+0.58883{\mathrm{\&delta;}}_G-6.88253& if& {\mathrm{\&delta;}}_G>1.96\\ -0.21545{\mathrm{\&delta;}}_G^2+1.08602{\mathrm{\&delta;}}_G-0.32634& if& {\mathrm{\&delta;}}_G\&le;1.96\end{array}\right.$

其中，δ_G＝10⁶δ_λ→∞。

5.根据权利要求4所述的欠采样的CT图像重建方法，其特征在于，步骤5)中，所述加权惩罚因子w_s的表达式为：

$w_s=\frac{mean\left(\right|\left|{\mathrm{\&mu;}}_s\right|{|}_{TV})}{\left|\right|{\mathrm{\&mu;}}_s|{|}_{TV}+\mathrm{\&epsiv;}}$

其中，μ_s表示从图像μ中抽取出的第s个图像块，且μ_s＝reshape(E_sμ,N₀,N₀)，其表示从图像μ中提取N₀×N₀维图像块的算子；reshape函数表示将E_sμ重新排列为二维图像块；||μ_s||_TV为图像块的全变分，计算公式为G为梯度算子，(Gx)_i,j＝(x_i+1,j-x_i,j,x_i,j+1-x_i,j)，x代表任意矩阵；mean(||μ_s||_TV)表示所有图像块全变分的平均值；ε＞0。

6.根据权利要求5所述的欠采样的CT图像重建方法，其特征在于，步骤6)中，所述目标函数一的表达式为：

$\underset{\mathrm{\&mu;},\mathrm{\&alpha;},D}{min}\underset{i=1}{\overset{I}{\&Sigma;}}\frac{{\mathrm{\&omega;}}_i}{2}{\left({\&lsqb;A\mathrm{\&mu;}\&rsqb;}_i-{\hat{l}}_i\right)}^2+\mathrm{\&lambda;}\underset{s}{\&Sigma;}{w}_s\left(\right||E_s\mathrm{\&mu;}-{\mathrm{D\&alpha;}}_s|{|}_2^2+{\mathrm{\&nu;}}_s\left|\right|{\mathrm{\&alpha;}}_s\left|{|}_0\right).$

7.根据权利要求6所述的欠采样的CT图像重建方法，其特征在于，步骤7)中，采用与获得初步重建图像相同的方法对目标函数一进行最小化求解，得到的图像μ的表达式为：

${\mathrm{\&mu;}}_j^{(t+1)}={\&lsqb;{\mathrm{\&mu;}}_j^t-\frac{\underset{i=1}{\overset{I}{\&Sigma;}}\left(a_{ij}{\mathrm{\&omega;}}_i\right({\&lsqb;{\mathrm{A\&mu;}}^{\left(t\right)}\&rsqb;}_i-{\hat{l}}_i\left)\right)+2\mathrm{\&lambda;}\underset{s}{\&Sigma;}{w}_s\underset{n=1}{\overset{N_0^2}{\&Sigma;}}{e}_{nj}^s({\&lsqb;E_s{\mathrm{\&mu;}}^{\left(t\right)}\&rsqb;}_n-{\&lsqb;{\mathrm{D\&alpha;}}_s\&rsqb;}_n)}{\underset{i=1}{\overset{I}{\&Sigma;}}\left(a_{ij}{\mathrm{\&omega;}}_i\underset{k=1}{\overset{N^2}{\&Sigma;}}{a}_{ik}\right)+2\mathrm{\&lambda;}\underset{s}{\&Sigma;}{w}_s\underset{n=1}{\overset{N_0^2}{\&Sigma;}}{e}_{nj}^s\underset{k=1}{\overset{N^2}{\&Sigma;}}{e}_{nk}^s}\&rsqb;}_{+}$

其中，加权惩罚因子w_s的初始值取1。

8.根据权利要求7所述的欠采样的CT图像重建方法，其特征在于，步骤8)中，所述加权惩罚因子w_s的更新公式为：

${\mathrm{\&mu;}}_s^{\left(t\right)}=reshape(E_s{\mathrm{\&mu;}}^{\left(t\right)},N_0,N_0)$

$w_s^{\left(t+1\right)}=\frac{mean\left(\right|\left|{\mathrm{\&mu;}}_s^{\left(t\right)}\right|{|}_{TV})}{\left|\right|{\mathrm{\&mu;}}_s^{\left(t\right)}|{|}_{TV}+\mathrm{\&epsiv;}}.$

9.根据权利要求8所述的欠采样的CT图像重建方法，其特征在于，在步骤9)中，迭代终止条件为两次迭代的扫描数据相对误差增量小于给定阈值T₁，迭代终止条件的表达式为：

$\frac{|\sqrt{\underset{i=1}{\overset{I}{\&Sigma;}}\frac{{\mathrm{\&omega;}}_i}{2}{({\&lsqb;{\mathrm{A\&mu;}}^{(t+1)}\&rsqb;}_i-{\hat{l}}_i)}^2}-\sqrt{\underset{i=1}{\overset{I}{\&Sigma;}}\frac{{\mathrm{\&omega;}}_i}{2}{({\&lsqb;{\mathrm{A\&mu;}}^{\left(t\right)}\&rsqb;}_i-{\hat{l}}_i)}^2}|}{\sqrt{\underset{i=1}{\overset{I}{\&Sigma;}}\frac{{\mathrm{\&omega;}}_i}{2}{({\&lsqb;{\mathrm{A\&mu;}}^{(t+1)}\&rsqb;}_i-{\hat{l}}_i)}^2}}<T_1.$

10.根据权利要求2所述的欠采样的CT图像重建方法，其特征在于，在步骤d)中，迭代终止条件为两次迭代的扫描数据相对误差增量小于给定阈值T₂，迭代终止条件的表达式为：

$\frac{|\sqrt{\underset{i=1}{\overset{I}{\&Sigma;}}\frac{{\mathrm{\&omega;}}_i}{2}{({\&lsqb;{\mathrm{A\&mu;}}^{(t+1)}\&rsqb;}_i-{\hat{l}}_i)}^2}-\sqrt{\underset{i=1}{\overset{I}{\&Sigma;}}\frac{{\mathrm{\&omega;}}_i}{2}{({\&lsqb;{\mathrm{A\&mu;}}^{\left(t\right)}\&rsqb;}_i-{\hat{l}}_i)}^2}|}{\sqrt{\underset{i=1}{\overset{I}{\&Sigma;}}\frac{{\mathrm{\&omega;}}_i}{2}{({\&lsqb;{\mathrm{A\&mu;}}^{(t+1)}\&rsqb;}_i-{\hat{l}}_i)}^2}}<T_2.$