基于混合正则项的稀疏角度CT图像重建

Sparse Angle CT Image Reconstruction Based on Mixed Regular Terms

摘要

针对传统迭代算法无法保证在数据欠采样情况下重建出高质量医学CT图像的问题， 提出了一种基于表面积最小化和群稀疏表示（GSR）的双正则项惩罚最小二乘的方法来用于稀疏角度CT图像重建。根据网格场中噪声图像比平滑图像表面积要大这一规律， 将表面积项作为重建图像的一个先验约束， 以达到平滑图像噪声的目的； 再根据GSR利用非局部相似性保留图像细微结构这一特性， 将GSR作为另一个先验约束， 建立重建模型。采用交替最小化方案将重建问题解耦为求解中间图像子问题和求解字典以及稀疏表示子问题， 之后利用半二次分裂法和分裂布雷格曼方法分别对两个子问题进一步求解。实验结果表明， 与算法PLS-GSR、 PLS-NLM、 SART-GSR相比， 本文算法的重建图像有良好的视觉效果， 可以较清晰地辨别图像的内部结构和细节， 并且FSIM值、 PSNR值和MSE值均优于对比算法。

Abstract

To solve the problem that the traditional iterative algorithms can not guarantee the reconstruction of high-quality medical CT images in the case of data undersampling， a method of double regular term penalty least squares based on surface area minimization and Group Sparse Representation （GSR） was proposed for sparse angle CT image reconstruction.According to the rule that the surface area of the noisy image is larger than that of the smooth image in the grid field， the surface area term was used as a prior constraint on reconstructing the image to achieve the purpose of smoothing the image noise； based on the feature that GSR preserves the fine structure of the image by using non-local similarity， the reconstruction model was established by using GSR as another prior constraint. An alternate minimization scheme was used to decouple the reconstruction problem into solving the intermediate image subproblem， and solving the dictionary and sparse representation subproblem. Then， the two sub-problems were further solved by the half-quadratic splitting method and the split Bregman method respectively. The experimental results show that compared with the algorithms PLS-GSR， PLS-NLM and SART-GSR， the reconstructed image of the proposed algorithm has good visual effect， which can clearly identify the internal structure and details of the image， and the value of FSIM， PSNR and MSE are superior to the comparison algorithm.

Graphical abstract

关键词

群稀疏表示 / 表面积 / 图像平滑 / 非局部相似性 / 图像重建

Key words

group sparse representation / surface area / image smoothing / non-local similarity / image reconstruction

0 引　言

CT成像技术已经成为病理检测的重要手段， 然而， 使患者置于辐射剂量较高的环境中接受X射线扫描， 会诱发癌症、 白血病或其他遗传性疾病。为了降低CT扫描辐射剂量对患者的影响， 检测中广泛采用了稀疏角CT重建的方法。这种方法通过使用有限数量的投影数据进行CT图像重建， 从而提高了图像重建速度。经典的解析算法如滤波反投影（FBP）凭借其重建速度快、 操作简便易学的优点， 在临床CT成像中得到持续广泛的应用。然而， 由于稀疏角度采样缺乏大量的投影数据， 该算法重建的图像往往无法避免地呈现出过度噪声和条纹伪影， 在低剂量CT成像时性能略显不足。相比而言， 迭代重建算法对数据的完备性要求不高， 可以灵活地引入图像先验约束信息^［1］， 通过先验知识的引导， 迭代重建算法在投影数据不完备的情况下仍然可以得到质量较高的重建图像。

压缩感知（CS）理论的兴起为稀疏角度投影的CT图像重建提供了新的理论依据和思路^［2］。CS理论指出如果一个未知信号在某个变换域中是稀疏的， 那么可用少于Nyquist采样数的测量值恢复原始信号。一种常见的稀疏变换是离散梯度变换， 其系数相加形成全变分（TV）。Sidky等^［3］利用TV最速下降法和凸集投影约束相结合的方法进行了CBCT图像重建， 该算法可在一定程度上消除条纹伪影， 但假设信号是分段平滑的， 会导致其重建图像出现块状阶梯和过平滑效应。Kim等^［4］提出了一类新的全变分， 结合了一阶和二阶导数， 以很好地保持图像的平滑度， 并可以平衡一阶导数引起的阶梯伪影现象。Xi等^［5］利用二阶梯度来构造二阶TV范数， 其注重相邻像素间各项异性的边缘特性， 改善了传统TV算法重建图像块伪影严重的问题。针对TV算法仅考虑水平和垂直方向的梯度而缺乏指向性的问题， Yu等^［6］提出一种基于四向梯度算子的TV模型， 使得重建图像的边缘更加清晰。尽管基于TV的算法在许多情况下是有效的， 但是TV正则项是一个全局约束不能直接反映对象的结构。

相比逐像素更新的方法， 字典学习逐块处理图像， 具有能够更好反映图像局部结构的优点。Gong等^［7］提出了一种用于高光谱图像（HSI）去噪的低秩张量字典学习 （LTDL） 方法， 该方法考虑了一种接近低秩的近似， 更接近真实HSI干净组的潜在低秩结构， 具有优越的去噪性能。 Komolafe等^［8］提出了一种平滑的l_{\mathrm{0}}-范数约束字典学习算法， 该算法在目标函数中引入l_{\mathrm{0}}-范数作为正则项， 在每次迭代过程中， 对中间图像进行平滑处理后作为迭代输出， 避免了传统字典学习在迭代过程中产生伪影的问题。Gui等^［9］突破传统字典学习中所有图像块均由同一个字典表示的限制， 提出将图像块分类为光滑块和细节块后再学习字典的低剂量CT重建方法。

字典学习方法计算复杂度相对较高， 并且只独立考虑每个图像块的内在结构， 忽略了非局部但具有相似特征的图像块之间的关系， 而非局部相似先验可以更好地利用图像块之间的相关性。非局部均值算法（NLM）最初利用非局部相似性实现图像去噪。Zhang等^［10］提出一种将不同的搜索窗口、 块大小和滤波器平滑参数应用于不同区域的自适应NLM方法。刘进等^［11］提出了稀疏张量约束的低剂量CT图像重建方法， 该方法利用非局部相似的先验信息， 将CT图像分割成一系列图像块组； 通过采用张量的多维低秩分解方法， 将这一先验信息引入低剂量CT重建中， 提高了图像重构质量。Zhang等^［12］将具有结构相似性的非局部块组成的群作为稀疏表示的基本单位， 建立了一种群稀疏表示模型， 并设计了一种高效的自适应字典学习方法， 该算法在图像恢复中具有很强的鲁棒性。Komolafe等^［13］将群稀疏应用于双能乳腺CT重建， 把联合代数重建算法（SART）的输出图像作为初始图像， 将其纳入群稀疏的自适应字典学习框架中， 再由稀疏表示和字典重构出精度较高的图像。 姜斌等^［14］在图像去噪模型中引入基于l_{\mathrm{2,1}}范数的群稀疏正则项， 保留了图像边缘特征， 从被混合噪声污染的高光谱图像中恢复出干净的图像。Cai等^［15］将群稀疏表示和低秩群残差学习应用于图像恢复， 噪声和块状伪影得到了有效抑制。

为进一步解决现有稀疏角度CT重建算法噪声平滑不彻底和细节保留不完全的问题， 本文提出一种基于表面积最小化和群稀疏表示的混合正则项CT重建算法， 采用交替最小化方案来优化目标函数。表面积正则项可以直接性地引导重建模型消除噪声； 群稀疏作为传统稀疏表示的延伸， 充分利用整幅图像中相似图像块的结构信息， 以群为单位学习自适应字典， 获取图像更精确的稀疏表示， 在提高图像细节与边缘的清晰度方面表现卓越。通过实验验证， 表面积项与群稀疏相结合的算法在平滑伪影噪声的同时， 还较为清晰地呈现了精细结构。

1 理论依据

1.1 惩罚最小二乘法（PLS）CT重建

基于惩罚最小二乘法的CT图像重建模型可表示为

\underset{x}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{‖\mathit{A}x-p‖}_{\mathrm{2}}^{\mathrm{2}}+\beta R(x),（1）

式中： {‖\mathit{A}x-p‖}_{\mathrm{2}}^{\mathrm{2}}为数据保真项； R(x)为惩罚项； \beta为平衡参数； p为投影数据； x为待重建的图像矢量； \mathit{A}是大小为P\times Q的系统矩阵（P是射线总数， Q是图像x的像素总数）。

1.2 图像稀疏表示（Sparse Representions，SR）

图像内容的有效表示是图像处理领域的一个核心问题。在某些情况下， 图像在某一过完备字典上可以由极少量原子的线性组合表示， 称为稀疏表示（SR）， SR已广泛应用于图像去噪和图像恢复领域， 作为先验信息来抑制接近零的噪声， SR可以建模为

\left\{\mathit{\alpha },D\right\}=\mathrm{a}\mathrm{r}\mathrm{g}\underset{D,\alpha }{\mathrm{m}\mathrm{i}\mathrm{n}}\frac{\mathrm{1}}{\mathrm{2}}{‖x-D\mathit{\alpha }‖}_{\mathrm{2}}^{\mathrm{2}}+\lambda {‖\mathit{\alpha }‖}_{\mathrm{0}},（2）

式中： \sqrt[]{b^{\mathrm{2}}-\mathrm{4}ac}为图像矢量； D为字典； \mathit{\alpha }为稀疏矩阵； {‖\cdot ‖}_{\mathrm{0}}为l_{\mathrm{0}}-范数。字典D的选取是SR逼近真实图像特征的关键， 常见的固定字典有小波字典、 曲波字典等， 此外， 还可以通过MOD、 K-SVD等算法学习自适应字典并得到稀疏表示。

1.3 基于群的图像稀疏表示（GSR）

传统SR方法将图像分割为重叠图像块后逐块处理， 没有考虑到相似图像块之间的有效信息。而GSR充分地利用了图像的非局部自相似性和局部稀疏性， 由具有相似结构的图像块构建为相似组， 从而将所有图像块转换到群域中， 在群域中实现对整个图像的稀疏表示。

1.3.1 群的构建

首先， 以4像素为滑动距离， 把图像x分割成n个大小为\sqrt[]{B}\times \sqrt[]{B}（B为图像块像素数）互相重叠的小块x_k(k=\mathrm{1,2},\cdots ,n)， 采用欧氏距离为相似度度量标准， 在大小为L\times L的搜索窗口内搜集与x_k相似度最高的m个图像块， 形成块集合G_{x_k}（在图 1 中， 紫色方框表示x_k， 绿色方框表示搜索窗口）。其次， 将G_{x_k}中所有块矢量化， 排列成大小为B\times m的矩阵{\mathit{x}}_{G_k}。定义E_{G_k}是从图像中提取{\mathit{x}}_{G_k}的算子， E_{G_k}^T是将x_{G_k}放回重建图像第k个位置的算子， 则{\mathit{x}}_{G_k}和重建图像x可以表示为

{\mathit{x}}_{G_k}=E_{G_k}(x),（3）

x={\displaystyle \sum _{k=\mathrm{1}}^n}{E}_{G_k}^T({\mathit{x}}_{G_k})/{\displaystyle \sum _{k=\mathrm{1}}^n}{E}_{G_k}^T({\mathrm{1}}_{B\times m}),（4）

式中： 运算符/表示两个矩阵对应元素的除法； {\mathrm{1}}_{B\times m}是元素都为1的矩阵。

1.3.2 GSR建模

首先， 用r表示x的近似估计， r_{G_k}表示{\mathit{x}}_{G_k}的近似估计（r和r_{G_k}由下文算法得到）。对r_{G_k}进行奇异值分解（SVD）， 即

r_{G_k}=U_{G_k}{\mathrm{\Sigma }}_{G_k}{V}_{G_k}^T={\displaystyle \sum _{i=\mathrm{1}}^c}{\gamma }_{r_{G_k\times i}}\left(u_{G_k\times i}{v}_{G_k\times i}^T\right)=

{\displaystyle \sum _{i=\mathrm{1}}^c}{\gamma }_{r_{G_k\times i}}{d}_{G_k\times i},（5）

式中： c为字典原子数； {\gamma }_{r_{G_k}}=\{{\gamma }_{r_{G_k\times \mathrm{1}}},{\gamma }_{r_{G_k\times \mathrm{2}}},\cdots ,{\gamma }_{r_{G_k\times c}}\}；

{\mathrm{\Sigma }}_{G_k}=diag({\gamma }_{r_{G_k}})且d_{G_k\times i}=u_{G_k\times i}{v}_{G_k\times i}^T。

由此， 得到自适应字典D_{G_k}， 即

D_{G_k}=\left\{d_{G_k\times \mathrm{1}},d_{G_k\times \mathrm{2}},\cdots ,d_{G_k\times c}\right\}\mathrm{。}（6）

再得到{\mathit{x}}_{G_k}在字典D_{G_k}上的稀疏编码{\alpha }_{G_k}， 则{\mathit{x}}_{G_k}可以表示为

{\mathit{x}}_{G_k}=D_{G_k}{\alpha }_{G_k}\mathrm{。}（7）

整幅图像的x可以由稀疏编码集合\{{\alpha }_{G_k}\}表示为

x=D_G{\alpha }_G={\displaystyle \sum _{k=\mathrm{1}}^n}{E}_{G_k}^T(D_{G_k}{\alpha }_{G_k})/{\displaystyle \sum _{k=\mathrm{1}}^n}{E}_{G_k}^T({\mathrm{1}}_{B\times m}),（8）

式中： D_G表示所有D_{G_k}的级联； {\alpha }_G表示所有{\alpha }_{G_k}的级联。

GSR模型就是寻找D_{G_k}和{\alpha }_{G_k}， 即

\{{\alpha }_G,D_G\}=\underset{{\alpha }_G,D_G}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\frac{\mathrm{1}}{\mathrm{2}}{‖x-D_G{\alpha }_G‖}_{\mathrm{2}}^{\mathrm{2}}+\lambda {‖{\alpha }_G‖}_{\mathrm{0}}\mathrm{。}（9）

1.4 表面积最小化约束

图 2（e）和图2（f）在三维网格化空间中， 分别展示了图 2（a）（清晰图像）和图2（b）（噪声图像）的感兴趣区域图 2（c） 和图 2（d） 中像素点所对应的像素值。

若将其看作曲面， 显然图 2（d） 在网格场中的空间曲面凹凸不平， 而图 2（a） 在网格场中的曲面起伏很少， 图 2（e） 中曲面的表面积显然小于图 2（f）。因此， 表面积项可以作为一种更高效的先验约束以达到保持图像平滑性的目的。

根据微积分学^［16］对曲面面积的定义， 光滑曲面z=z(x,y)的面积可以由下述积分表示为

S={\iint }_{\mathrm{\Omega }}^{}\sqrt[]{\mathrm{1}+{\left(\frac{\partial z}{\partial x}\right)}^{\mathrm{2}}+{\left(\frac{\partial z}{\partial y}\right)}^{\mathrm{2}}}\mathrm{d}x\mathrm{d}y,（10）

式中： \mathrm{\Omega }表示该曲面在Oxy平面上的投影。

图像x在网格场中的曲面面积可表示为

{\rm T}(\nabla x)={\displaystyle \sum _{j=\mathrm{1}}^M}{\displaystyle \sum _{i=\mathrm{1}}^N}\sqrt[]{\mathrm{1}+{\left|\nabla x\right|}_{i,j}^{\mathrm{2}}},（11）

式中： M和N分别表示图像x的宽和高； \left|\nabla x\right|=\sqrt[]{{\left|{\nabla }_{h}x\right|}^{\mathrm{2}}+{\left|{\nabla }_{v}x\right|}^{\mathrm{2}}}， {\nabla }_h为水平梯度算子， {\nabla }_v为垂直梯度算子。

2 本文算法

GSR基于图像的非局部相似性， 在很大程度上保留了图像的细节信息。表面积最小化约束使得图像更加平滑， 因此， 本文将GSR和表面积项作为PLS的惩罚项， 提出一种混合约束最小二乘模型， 即

\underset{x}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{‖\mathit{A}x-p‖}_{\mathrm{2}}^{\mathrm{2}}+\eta {\rm T}(\nabla x)+

\beta (\frac{\mathrm{1}}{\mathrm{2}}{‖x-D_G{\alpha }_G‖}_{\mathrm{2}}^{\mathrm{2}}+\lambda {‖{\alpha }_G‖}_{\mathrm{0}}),（12）

式中： \eta和\beta为正则项的系数； \lambda为GSR中的正则项系数。

本文采用交替最小化方案求解模型（12）， 分为两步， 具体如下：

步骤1： 固定D_G和{\alpha }_G， 更新迭代中间图像x， 即求解如下问题：

\underset{x}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{‖\mathit{A}x-p‖}_{\mathrm{2}}^{\mathrm{2}}+\eta {\rm T}(\nabla x)+

\beta {‖x-D_G{\alpha }_G‖}_{\mathrm{2}}^{\mathrm{2}}\mathrm{。}（13）

运用半二次分裂法优化式（13）， 引入一个辅

助变量U， 使问题（13）近似于

\underset{x,U}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\mathrm{\Phi }(x,U)=\underset{x}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{‖\mathit{A}x-p‖}_{\mathrm{2}}^{\mathrm{2}}+

\theta {‖U-\nabla x‖}_{\mathrm{2}}^{\mathrm{2}}+\eta {\rm T}(U)+\beta {‖x-D_G{\alpha }_G‖}_{\mathrm{2}}^{\mathrm{2}},（14）

式中： \theta为惩罚参数。

然后采用交替最小化方案求解问题（14）。

给定x， \Phi (x,U)关于U的最小化问题为

\underset{U}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\theta {‖U-\nabla x‖}_{\mathrm{2}}^{\mathrm{2}}+\eta {\rm T}(U)\mathrm{。}（15）

这个问题对于所有M\times N个像素点是可分离的， 所以采用并行计算， 求解如下一维优化问题：

\underset{U_{i,j}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}H(U_{i,j})=\theta {‖U_{i,j}-v_{i,j}‖}_{\mathrm{2}}^{\mathrm{2}}+\eta \sqrt[]{\mathrm{1}+{\left|U_{i,j}\right|}^{\mathrm{2}}},（16）

式中： v_{i,j}=\nabla x_{i,j}。

使用牛顿迭代法对问题（16）求解， 即通过式（17）更新U_{i,j}。

U_{i,j}^{(l+\mathrm{1})}=U_{i,j}^{(l)}-\frac{H(U_{i,j}^{(l)})}{H^{\text{'}}(U_{i,j}^{(l)})},（17）

直到达到迭代条件为止。

给定U， 通过求解如下问题得到x。

\underset{x}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{‖\mathit{A}x-p‖}_{\mathrm{2}}^{\mathrm{2}}+\theta {‖U-\nabla x‖}_{\mathrm{2}}^{\mathrm{2}}+

\beta {‖x-D_G{\alpha }_G‖}_{\mathrm{2}}^{\mathrm{2}}\mathrm{。}（18）

使用梯度下降法求解问题（18）， 即

x^{(q+\mathrm{1})}=x^{(q)}-\xi ({\mathit{A}}^{\mathrm{T}}(\mathit{A}{x}^{(q)}-p)+

\theta {\nabla }^{\mathrm{T}}(\nabla x^{(q)}-U)+\beta (x^{(q)}-D_G{\alpha }_G))，（19）

式中： \xi为学习率。

步骤 2： 固定x， 更新字典D_G和稀疏编码{\alpha }_G， 即求解如下问题：

\underset{D_G,{\alpha }_G}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}\frac{\mathrm{1}}{\mathrm{2}}}{‖x-D_G{\alpha }_G‖}_{\mathrm{2}}^{\mathrm{2}}+\lambda {‖{\alpha }_G‖}_{\mathrm{0}}\mathrm{。}（20）

引入变量y， 式（20）等价于

\underset{y,D_G,{\alpha }_G}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\frac{\mathrm{1}}{\mathrm{2}}{‖y-x‖}_{\mathrm{2}}^{\mathrm{2}}+\lambda {‖{\alpha }_G‖}_{\mathrm{0}},y=D_G{\alpha }_G\mathrm{。}（21）

由于D_G可以由y的近似得出， 则式（21）等价于

\underset{y,{\alpha }_G}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\frac{\mathrm{1}}{\mathrm{2}}{‖y-x‖}_{\mathrm{2}}^{\mathrm{2}}+\lambda {‖{\alpha }_G‖}_{\mathrm{0}},y=D_G{\alpha }_G，（22）

式中： y,{\alpha }_G是优化变量， y=D_G{\alpha }_G为等式约束， 其中， \frac{\mathrm{1}}{\mathrm{2}}{‖y-x‖}_{\mathrm{2}}^{\mathrm{2}}是数据保真项， \lambda是正则化参数。采用Split-Bregman方法将式（22）解耦为几个子问题， 进而求解y和{\alpha }_G， 具体为

y^{(t+\mathrm{1})}=\underset{y}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\frac{\mathrm{1}}{\mathrm{2}}{‖y-x‖}_{\mathrm{2}}^{\mathrm{2}}+\frac{\rho }{\mathrm{2}}{‖y-D_G{\alpha }_G^{(t)}-w^{(t)}‖}_{\mathrm{2}}^{\mathrm{2}}，（23）

{{\alpha }_G}^{(t+\mathrm{1})}=\underset{{\alpha }_G}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\lambda ‖{\alpha }_G‖{}_{\mathrm{0}}{}^{}+

\frac{\rho }{\mathrm{2}}{‖y^{(t)}-D_G{\alpha }_G-w^{(t)}‖}_{\mathrm{2}}^{\mathrm{2}}，（24）

w^{(t+\mathrm{1})}=w^{(t)}-y^{(t+\mathrm{1})}+D_G{\alpha }_G^{(t+\mathrm{1})}\mathrm{。}（25）

根据文献［17］可知， 式（23）的解为

y^{(t+\mathrm{1})}=\frac{x+\rho D_G{\alpha }_G^{(t)}+\rho w^{(t)}}{\mathrm{1}+\rho }\mathrm{。}（26）

式（24）等价于

{{\alpha }_G}^{(t+\mathrm{1})}=\underset{{\alpha }_G}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\frac{\lambda }{\rho }‖{\alpha }_G‖{}_{\mathrm{0}}{}^{}+\frac{\mathrm{1}}{\mathrm{2}}{‖r^{(t)}-D_G{\alpha }_G‖}_{\mathrm{2}}^{\mathrm{2}}，（27）

式中： r^{(t)}=y^{(t)}-w^{(t)}为y^{(t)}的近似。由文献［12］可知， 式（27）中稀疏编码{\alpha }_G^{（t+\mathrm{1}）}可以由{\alpha }_k^{（t+\mathrm{1}）}的级联得到， {\alpha }_{G_k}^{(t+\mathrm{1})}可由式（28）计算得出。

{{\alpha }_{G_k}}^{(t+\mathrm{1})}=hard({\gamma }_{r_{G_k}}^{(t)},\sqrt[]{\mathrm{2}\tau })（28）

式中： hard(\cdot )为硬阈值函数； {\gamma }_{r_{G_k}}^{(t)}由式（5）得到； \tau =(\lambda \times B\times m\times n)/(\rho \times Q)。

本文算法的伪代码如下：

Input:A,p,\eta ,\beta ,\theta ,\epsilon ,\lambda ,\rho ,\xi ,m,n,B,L,N^{iter},n^{iter}

Initializtion:U^{\mathrm{0}}=x^{\mathrm{0}}=w^{\mathrm{0}}=\mathrm{0},D_G^{\mathrm{0}}=\mathrm{0},{\alpha }_G^{\mathrm{0}}=\mathrm{0}

for n=1， 2， …， Iter

step1：fix D_G， {\alpha }_G

Update U through Eq.(17)

Update x through Eq.(19)

step2：fix x

for t=1， 2， …， iter

y^{\left(\mathrm{0}\right)}=x

Update D_{G_k}^{(t+\mathrm{1})}， {\alpha }_{G_k}^{(t+\mathrm{1})}through Eq.(6)， (28) respectively

Update D_G^{(t+\mathrm{1})}， {\alpha }_G^{(t+\mathrm{1})} through concatenating all D_{G_k}^{(t+\mathrm{1})}，

{\alpha }_{G_k}^{(t+\mathrm{1})} respectively

end

end

3 实验结果

对所提出的重建算法进行性能验证实验， 测试平台为： 64位Intel Core i7-10700处理器、 主机频率2.90 GHz， 8G运行内存的Windows10操作系统， 数学软件为Matlab R2021a。给出了Shepp-Logan头骨模型、 胸腔和腹部图像在64个投影角度下的实验结果。为了进一步评估本文算法的性能， 将本文算法与PLS-GSR^［18］、 PLS-NLM^［19］、 SART-GSR^［17］进行对比， 并使用结构相似性（FSIM）、 峰值信噪比（PSNR）和均方误差（MSE）作为图像质量的评价指标。

3.1 头骨模型重建结果分析

图 3（a）~图 3（e） 分别展示了头骨模型原始图像和由PLS-GSR、 PLS-NLM、 SART-GSR、 本文算法等4种算法重建出来的全局图像和感兴趣域。PLS-GSR算法重建出的图像边缘较清晰， 但存在大量噪声， 整体不够平滑； PLS-NLM在噪声平滑上有着很好的性能， 但在边缘处过度平滑， 而SART-GSR在平滑噪声和保留边缘上较前二者有较好的效果， 但在边缘处还是不够清晰； 本文算法去除噪声和保持边缘的能力都很突出。

图 4（a）~图 4（c） 分别展示了算法在迭代50次过程中重建图像与原始图像的FSIM值、 PSNR值和MSE值的比较结果及变化趋势， 更好地体现了所提算法的有效性。

3.2 胸腔图像重建结果分析

图 5（a）~图 5（e） 分别展示了胸腔图像原始图像和由PLS-GSR、 PLS-NLM、 SART-GSR、 本文算法等4种算法重建出来的全局图像和感兴趣域。PLS-GSR算法重建出的图像中， 胸腔细节被保留但伪影噪声非常严重， 不能清晰地观察到胸腔内部的血管和组织； PLS-NLM算法重建的图像虽然很好地规避了噪声严重的问题， 但丢失了很多重要的结构细节， 对胸腔内部的真实环境还原度极低； SART-GSR算法则保留了较多图像细节和图像边缘； 本文算法较前者保留了更多的图像细节（如红色箭头所指）， 且边缘保留效果更加显著。

图 6（a）~图 6（c） 分别展示了算法在迭代50次过程中重建图像与原始图像FSIM值、 PSNR值和MSE值的比较结果及变化趋势， 直观地展现出本文所提算法的有效性。

3.3 腹部图像重建结果分析

图 7（a）~图 7（e） 分别展示了胸腔图像原始图像和由PLS-GSR、 PLS-NLM、 SART-GSR、 本文算法等4种算法重建出来的全局图像和感兴趣域。PLS-GSR算法保留了较全面的细节但仍然有条纹伪影； PLS-NLM算法对图像的细节和边缘都过于平滑， 不具备医学诊断参考价值； SART-GSR算法较完整地保留了组织及其边缘； 本文算法不仅平滑了伪影噪声， 还保留了细微组织结构（如红色箭头所指）。图 8（a）~图 8（c） 分别展示了算法在迭代50次过程中重建图像与原始图像的FSIM值、 PSNR值和MSE值的比较结果及变化趋势， 可以看出本文算法在第50次迭代时的3个指标的值都比其余算法高。

表 1~表 3 分别给出各算法基于头骨、 胸腔和腹部的FSIM、PSNR和MSE这3种指标的对比结果， 可以发现所提算法得到的FSIM值较其他算法更接近于1甚至等于1， 这说明本文算法在CT重建中还原了更多的结构特征和图像细节。同时， 本文算法与PLS-GSR、PLS-NLM、SART-GSR算法相比， 重建头骨模型的PSNR值分别提高了13.501 5， 7.589 4， 6.255 9 dB； 重建胸腔图像的PSNR值分别提高了7.420 7，7.507 7，2.981 6 dB； 重建腹部图像的PSNR值分别提高了4.015 9，6.778 9，2.202 6 dB； 另外， MSE值最低也验证了本文算法在平滑噪声和还原图像真实度方面具有良好的性能。

4 结　论

在医学图像重建中， 平滑噪声和保留图像细节都是设计重建模型过程中需要重点考虑的因素。本文提出了一种基于表面积约束和群稀疏正则项的惩罚最小二乘算法。表面积最小约束提升了算法去除噪声的性能， 群稀疏正则项考虑了图像的非局部相似性， 很好地保留了图像的细微结构。各项评价指标和重建结果表明， 该算法在数据欠采样情况下， 既能平滑噪声又能保留较为完整的组织和细节， 有着良好的性能。后续研究工作将对此方法的运行速度和鲁棒性进行优化， 以期设计出兼容性更高的医学图像重建算法。

参考文献

[1] 张莹芳.基于全变分的稀疏角度CT重建保边算法的研究［D］.太原： 中北大学， 2023.

[2] DONOHO D L. Compressed sensing［J］. IEEE Transactions on Information Theory， 2006， 52（4）： 1289-1306．

[3] SIDKY E Y， KAO C M， PAN X C. Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT［J］. Journal of X-ray Science and Technology， 2006， 14（2）： 119-139.

[4] KIM Y， KUDO H. Nonlocal total variation using the first and second order derivatives and its application to CT image reconstruction［J］. Sensors， 2020， 20（12）： 3494.

[5] XI Y R， ZHOU P W， YU H J， et al. Adaptive‐weighted high order TV algorithm for sparse‐view CT reconstruction［J］. Medical Physics， 2023， 50（9）： 5568-5584.

[6] YU Z S， WEN X Y， YANG Y. Reconstruction of Sparse-view X-ray computed tomography based on adaptive total variation minimization［J］. Micromachines， 2023， 14（12）： 2245.

[7] GONG X， CHEN W， CHEN J. A low-rank tensor dictionary learning method for hyperspectral image denoising［J］. IEEE Transactions on Signal Processing，2020，68：1168-1180.

[8] KOMOLAFE T， WANG KANG， DU QIANG， et al. Smoothed L0-constraint dictionary learning for low-dose X-ray CT reconstruction［J］.IEEE Access 2020，8： 116961-116973.

[9] GUI Y J， ZHAO X， BAI Y J， et al. Low-dose CT iterative reconstruction based on image block classification and dictionary learning［J］. Signal， Image and Video Processing， 2023， 17（2）： 407-415.

[10] ZHANG P C， LIU Y， GUI Z G， et al. A region-adaptive non-local denoising algorithm for low-dose computed tomography images［J］. Mathematical Biosciences and Engineering， 2023， 20（2）： 2831-2846.

[11] LIU JIN， KANG Y Q， GU Y B， et al. Low dose computed tomography image reconstruction based on sparse tensor constraint［J］. Acta Optica Sinica， 2019， 39（8）： 0811004.

[12] ZHANG J， ZHAO D B， GAO W. Group-based sparse representation for image restoration［J］. IEEE Transactions on Image Processing， 2014， 23（8）： 3336-3351.

[13] KOMOLAFE T， DU Q， ZHANG Y， et al. Material decomposition for simulated dual-energy breast computed tomography via hybrid optimization method［J］. Journal of X-Ray Science and Technology， 2020， 28（6）： 1037-1054.

[14] 姜斌，叶军.基于群稀疏正则化的高光谱图像去噪［J］.计算机技术与发展，2023，33（12）：171-177.

[15] JIANG Bin， YE Jun. Hyperspectral image denoising based on group sparse regularization［J］. Computer Technology and Development， 2023， 33（12）： 171-177. （in Chinese）

[16] CAI Z Y， XIE X H， DENG J J， et al. Image restoration with group sparse representation and low‐rank group residual learning［J］. IET Image Processing， 2024， 18（3）： 741-760.

[17] 邵文凯， 阮杰昌， 王晓平 等.微积分［M］.重庆： 重庆大学出版社： 2015.

[18] BAO P， ZHOU J L， ZHANG Y. Group Sparsity based Sparse-Sampling CT Reconstruction［C］//2018 40th Annual International Conference of the IEEE Engineering in Medicine and Biology Society （EMBC）. IEEE， 2018： 5097-5100.

[19] GU Y N， LIU Y， LIU W T， et al. Sparse angle CT reconstruction based on group sparse representation［J］. Journal of X-Ray Science and Technology， 2022， 30（6）： 1085-1097.

[20] SALEHJAHROMI M. Iterative Techniques for Spectral Computed Tomography［D］. Massachusetts： University of Massachusetts Lowell， 2019.

基金资助

国家自然科学基金资助项目(61774137)

山西省基础研究计划资助项目(202103021224195)

山西省基础研究计划资助项目(202103021224212)

山西省基础研究计划资助项目(202103021223189)

山西省基础研究计划资助项目(20210302123019)