基于加权引导滤波与梯度域卷积稀疏的CT重建

CT Reconstruction Based on Weighted Guided Filtering and Gradient Domain Convolution Sparseness

摘要

对于不完整的医学CT扫描数据， 传统算法无法保证重建图像满足诊断要求。针对这种情况， 提出了一种基于加权引导滤波与梯度域卷积稀疏编码结合的CT重建算法。该算法首先采用惩罚最小二乘法迭代重建初始CT图像； 其次， 利用加权引导滤波获取图像的低频分量， 采用带有梯度约束的卷积稀疏编码处理图像的高频分量； 最后， 将两段分量相结合得到的新图像作为输入继续进行最小二乘逼近， 反复迭代重建， 直到获得更清晰的图像。实验结果表明， 与其他卷积稀疏算法及组稀疏算法相比， 该算法可有效抑制噪声和伪影， 恢复更多图像的结构和边缘细节信息， 获得更优的重建图像。

Abstract

For incomplete medical CT scan data， traditional algorithms cannot ensure that the reconstructed images meet the diagnostic requirements. In order to solve this situation， a CT reconstruction algorithm based on the combination of weighted guided filtering and gradient domain convolutional sparse coding was proposed. Firstly， the penalty least squares method is used to iteratively reconstruct the initial CT image. Secondly， the weighted guided filtering is used to obtain the low-frequency components of the image， and the convolutional sparse coding with gradient constraints is used to process the high-frequency components of the image. Finally， the two components are combined to obtain a new image as input， and the least-squares approximation is continued， and the iterative reconstruction is repeated until a clearer image is obtained. Experimental results show that compared with other convolutional sparse algorithms and group sparse algorithms， the proposed algorithm can effectively suppress noise and artifacts， recover more image structure and edge detail information， and obtain better reconstructed images.

Graphical abstract

关键词

计算机断层成像 / 卷积稀疏编码 / 加权引导滤波 / 稀疏角度 / 图像重建

Key words

computed tomography / convolutional sparse coding / weighted guided filtering / sparse angle / image reconstruction

J. Test Measur. Technol., 2025, 39(5): 558-564.

0 引　言

近几十年来， 计算机断层扫描（Computed Tomography， CT）已广泛应用于临床诊断。然而， X射线可以被人体部分吸收， 并可能导致遗传或癌症疾病， 因此， 人们开始关注降低X射线辐射的潜在风险。减少CT扫描中X射线辐射的有效方法可分为两种： 一是降低X射线源处的光子强度， 但这会增加投影数据中的量子噪声^［1］； 二是减少CT设备获得的投影视图数量， 这将导致投影数据不足， 从而出现少视图或有限角度CT的情况。目前研究重点关注稀疏角度方法， 该方法也可以提高扫描速度。传统CT图像重建算法有两类， 一类是解析法， 如滤波反投影； 另一类是迭代法， 如联合代数重建， 其重建速度相对较慢， 但可以实现稀疏角度下的高精度重建。

压缩感知（Compressed Sensing， CS）理论的出现使得图像处理中的不适定问题可以得到有效解决^［2］。Sidky等^［3］提出将全变分（Total Variation， TV）和凸集投影（Projection Onto Convex Sets， POCS）结合来解决投影数据重建不完整的问题。虽然TV正则项能够很好地恢复图像的边缘信息， 但是过度平滑会导致重建图像中出现梯度伪影。为了解决这个缺点， 提出了几种TV变体， 包括基于梯度方向总变差^［4］、 TV原始对偶算法^［5］和自适应加权高阶TV^［6］。此外， 基于TV正则化的方法被用于CT图像投影领域， 以恢复局部细节信息。

图像在某些变换域（如小波、 离散余弦变换和离散梯度）或字典中通常是稀疏的。近年来， 字典学习已得到广泛应用， 包括超分辨率成像、 图像或视频去噪、 分类和医学图像重建等。然而， 传统的字典学习方法忽略了补丁之间的关系， 并且学习到的字典可能过于冗余。随后， 出现了组稀疏正则化（Group-based Sparse Representation， GSR）^［7］， 其中相似的补丁被聚类成组， 作为基于组稀疏表示的基本单元， 该方法有效地消除了条带伪影并保留了图像细节。然而， 组稀疏方法都是基于补丁的， 并且学习到的特征通常包含相同特征的移位。为了解决这些问题， 提出了卷积稀疏编码（Convolutional Sparse Coding， CSC）^［8］。CSC将滤波器应用于稀疏特征图以自然地重建整个图像， 有效解决了重叠补丁的问题。Bao等^［9］对CSC进一步改进， 提出了一种基于CSC的惩罚加权最小二乘法， 辅以稀疏视图 CT 重建特征图上的梯度正则化（PWLS-CSCGR）。

近年来， 出现了许多基于引导图像滤波引导图像滤波（Guided Image Filtering， GIF）的CT图像重建算法。GIF作为一种边缘保留平滑算子被提出， 与双边滤波相比， 表现出优越的性能和更快的滤波速度。作为局部过滤器， GIF 无法像基于全局优化的过滤器那样保留锐利边缘， 它会将模糊集中在边缘附近并引入光晕。Li等^［10］在GIF中引入了边缘感知权重， 创造了加权引导滤波（Weighted Guided Image Filtering， WGIF）。通过在引导图像中3×3的像素窗口内应用局部方差来计算边缘感知权重， 可以有效保留图像边缘。

传统的CSC往往忽略图像的低频分量， 仅用于处理高频特征， 然而有些强度较大的伪影会出现在图像的低频分量上。WGIF利用边缘像素点与平坦区域像素点的方差差异对引导滤波的归整化因子进行惩罚， 可以获得更好的边缘保持特性， 从而较大程度地去除条纹伪影。梯度域卷积稀疏编码（Convolutional Sparse Coding with Gradient Regularization， CSCGR）算法对特征图施加梯度约束， 可以避免因滤波器不准确产生而新的伪影。受上述研究启发， 本文提出了一种基于WGIF与CSCGR结合的CT重建算法。该算法将WGIF与梯度域卷积稀疏正则化相结合以获取更优的CT重建图像。经过实验验证， 本文算法对重建CT图像的质量有较大提升。

1 理论依据

1.1 CT重建模型

在单能束假设下， X射线CT重建模型可以用离散线性系统表示为

\mathit{P}=\mathit{A}\mathit{X},（1）

式中： 向量\mathit{P}为不同投影角度下测量的投影数据； 向量\mathit{X}为线性衰减系数估计， 即重建图像； \mathit{A}=\left\{a_{ij}\right\}\in {\mathrm{R}}^{I\times Q}为对应于CT图像的特定投影矩阵， I为投影数据总数， Q为图像像素总数， \mathit{A}中的元素a_{ij}表示第j个像素对第i条射线的贡献。但在实际应用中， 由于系统矩阵\mathit{A}很大， 直接对式（1）求逆无法得到重建图像， 且在稀疏角度投影的情况下， 式（1）的解是不适定的逆问题。为解决病态性、 欠定性问题， 一般采用在最小二乘模型中嵌入惩罚项对方程组进行求解。惩罚最小二乘法（Penalty Least Squares， PLS）的CT图像重建模型为

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

式中： 第1项为数据保真度项； 第2项为正则化项， 通常包含重建图像的先验信息； \beta为平衡保真度项和正则化项的参数。

1.2 CSCGR

卷积稀疏CSC假设图像可以表示为一组滤波器及其相应特征图之间的卷积之和。由于CSC对图像的低频分量不敏感， 因此只对图像的高频分量进行计算， 其目标函数为

\underset{\left\{M_i\right\},\left\{f_i\right\}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\frac{\mathrm{1}}{\mathrm{2}}{‖{\displaystyle \sum _{i=\mathrm{1}}^N}{\mathit{f}}_i\mathrm{*}{\mathit{M}}_i-{\mathit{X}}_h‖}_{\mathrm{2}}^{\mathrm{2}}+\lambda {\displaystyle \sum _{i=\mathrm{1}}^N}{‖{\mathit{M}}_i‖}_{\mathrm{1}},（3）

式中： {\mathit{X}}_h为图像的高频分量； \mathrm{*}表示卷积算子； {\left\{{\mathit{f}}_i\right\}}_{i=\mathrm{1,2},\dots ,N}为一组滤波器； {\left\{{\mathit{M}}_i\right\}}_{i=\mathrm{1,2},\dots ,N}为对应滤波器f_i的一组特征图。基于字典学习方法， 存在字典原子不准确的问题， 可能会出现结构丢失或产生新的伪影。同样， 如果存在不准确的过滤器， 基于CSC的方法同样会遇到这个问题。CSCGR是对特征图施加梯度约束， 采用梯度正则化约束来抑制离群点， 从而解决因滤波器不准确而造成的结构丢失或新伪影的问题^［11］。CSCGR在特征图上施加梯度约束， 构建目标函数为

\begin{array}{l}\underset{\left\{{\mathit{M}}_i\right\},\left\{{\mathit{f}}_i\right\}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\frac{\mathrm{1}}{\mathrm{2}}{‖{\displaystyle \sum _{i=\mathrm{1}}^N}{\mathit{f}}_i\mathrm{*}{\mathit{M}}_i-{\mathit{X}}_h‖}_{\mathrm{2}}^{\mathrm{2}}+\\ \lambda \left({\displaystyle \sum _{i=\mathrm{1}}^N}{‖{\mathit{M}}_i‖}_{\mathrm{1}}+\frac{\tau }{\mathrm{2}}{{\displaystyle \sum _{i=\mathrm{1}}^N}‖\sqrt[]{({\mathit{g}}_{\mathrm{1}}\mathrm{*}{\mathit{M}}_i{)}^{\mathrm{2}}+({\mathit{g}}_{\mathrm{2}}\mathrm{*}{\mathit{M}}_i{)}^{\mathrm{2}}}‖}_{\mathrm{2}}^{\mathrm{2}}\right),\end{array}（4）

式中： {\mathit{g}}_{\mathrm{1}}， {\mathit{g}}_{\mathrm{2}}分别为计算{\mathit{M}}_i沿x轴和y轴梯度的差分滤波器； \lambda为稀疏约束项的系数； \tau为超参数。

1.3 WGIF

引导滤波是一种局部线性模型图像滤波器^［12］。设引导图像为\mathit{G}， 滤波器输出\mathit{Z}是引导图像\mathit{G}在以像素k为中心的窗口{\omega }_k中的线性变换，

\mathit{Z}=a_k{\mathit{G}}_i+b_k,\forall i\in {\omega }_k,（5）

式中： 线性系数a_k， b_k为2个待求的参数。这里使用半径为r的方形窗口。通过在窗口中最小化以下问题获得系数，

E(a_k,b_k)=\sum _{i\in {\omega }_k}((a_k{\mathit{G}}_i+b_k-p_i{)}^{\mathrm{2}}+\eta a_k^{\mathrm{2}}),（6）

式中： p为输入图像； \eta为归整化因子， 用来防止系数a_k过大， 也是调节滤波器的重要参数。然而， 规整化因子\eta是固定的， 并未考虑不同窗口内像素之间的纹理差异。因此， 在图像边缘会因为较大的规整化因子而出现模糊， 降低了滤波效果。故将边缘感知加权加入到GIF中以形成WGIF。边缘感知权重是通过使用所有像素的3×3窗口的局部方差来定义的， 如下所示

{\Gamma }_G(k)=\frac{\mathrm{1}}{N}{\displaystyle \sum _{i\mathrm{\text{'}}=\mathrm{1}}^N}\frac{{\sigma }_{G,\mathrm{1}}^{\mathrm{2}}(k)+\epsilon }{{\sigma }_{G,\mathrm{1}}^{\mathrm{2}}(i\mathrm{\text{'}})+\epsilon },（7）

式中： {\sigma }_{G,\mathrm{1}}^{\mathrm{2}}(k)为引导图像\mathit{G}在以k为中心的3×3窗口{\omega }_k内的方差； i\mathrm{\text{'}}为图像\mathit{G}内的像素； \epsilon为一个常数； N为像素总数； 权重{\Gamma }_G(k)衡量整个引导图像中每个像素的重要性。然后， 将权重合并到成本函数E(a_k,b_k)中， 变为

E(a_k,b_k)=\sum _{i\in {\omega }_k}((a_k{G}_i+b_k-p_i{)}^{\mathrm{2}}+\frac{\eta }{{\Gamma }_G(k)}{a}_k^{\mathrm{2}})\mathrm{。}（8）

2 本文算法

基于前述分析， 本文将WGIF与CSCGR相结合用于构建新的CT重建算法， 如式（9）所示。此外， 该算法中使用了预先训练好的滤波器\left\{f_i\right\}。

\underset{\mathit{X},\left\{{\mathit{M}}_i\right\}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\frac{\mathrm{1}}{\mathrm{2}}{‖\mathit{A}\mathit{X}-\mathit{P}‖}_{\mathrm{2}}^{\mathrm{2}}+\beta \left(\frac{\mathrm{1}}{\mathrm{2}}{‖{\displaystyle \sum _{i=\mathrm{1}}^N}{\mathit{f}}_i\mathrm{*}{\mathit{M}}_i-{\mathit{X}}_h‖}_{\mathrm{2}}^{\mathrm{2}}+\right.\left.\lambda {\displaystyle \sum _{i=\mathrm{1}}^N}{‖{\mathit{M}}_i‖}_{\mathrm{1}}+\frac{\tau }{\mathrm{2}}{{\displaystyle \sum _{i=\mathrm{1}}^N}‖\sqrt[]{({\mathit{g}}_{\mathrm{1}}\mathrm{*}{\mathit{M}}_i{)}^{\mathrm{2}}+({\mathit{g}}_{\mathrm{2}}\mathrm{*}{\mathit{M}}_i{)}^{\mathrm{2}}}‖}_{\mathrm{2}}^{\mathrm{2}}\right)\mathrm{。}（9）

式（9）实质上是一个双目标优化问题， 本文采用交替最小化方案来获得其最优解。同时， 在求解过程中采用WGIF提取重建图像\mathit{X}的低频分量{\mathit{X}}_l， 而后将滤波器\left\{{\mathit{f}}_i\right\}与交替更新得到的特征图\left\{{\tilde{\mathit{M}}}_i\right\}卷积求和得到图像的高频分量{\mathit{X}}_h， 与{\mathit{X}}_l相加得到新的重建图像\mathit{X}。求解步骤如下：

1） 固定特征图\left\{{\tilde{\mathit{M}}}_i\right\}获得中间重建图像X

该阶段优化模型为

\underset{\mathit{X}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\frac{\mathrm{1}}{\mathrm{2}}{‖\mathit{A}\mathit{X}-\mathit{P}‖}_{\mathrm{2}}^{\mathrm{2}}+\beta \left(\frac{\mathrm{1}}{\mathrm{2}}{‖{\displaystyle \sum _{i=\mathrm{1}}^N}{\mathit{f}}_i\mathrm{*}{\tilde{\mathit{M}}}_i-\mathit{X}‖}_{\mathrm{2}}^{\mathrm{2}}\right)\mathrm{。}（10）

式（10）中{\mathit{X}}_h用\mathit{X}代替。这是一个二次优化问题。接下来采用HS共轭梯度下降法迭代求解出CT图像\mathit{X}， 具体迭代公式为

{\mathit{X}}^{(t+\mathrm{1})}={\mathit{X}}^{(t)}+{\mu }^{(t+\mathrm{1})}{\mathit{d}}^{(t+\mathrm{1})},（11）

式中： t为迭代次数； {\mu }^{(t+\mathrm{1})}为t+\mathrm{1}次更新后的最优步长； d^{(t+\mathrm{1})}为t+\mathrm{1}次后目标函数负梯度的共轭方向。式（10）的梯度表达式为

{\mathit{g}}^{(t+\mathrm{1})}={\mathit{A}}^{\mathrm{T}}(\mathit{A}{\mathit{X}}^{(t)}-\mathit{P})+\beta ({\mathit{X}}^{(t)}-{\displaystyle \sum _{i=\mathrm{1}}^N}{\mathit{f}}_i\mathrm{*}{\tilde{\mathit{M}}}_i),（12）

则负梯度的共轭方向为

{\mathit{d}}^{(t+\mathrm{1})}=-{\mathit{g}}^{(t+\mathrm{1})}+\frac{({\mathit{g}}^{(t+\mathrm{1})}{)}^{\mathrm{T}}({\mathit{g}}^{(t+\mathrm{1})}-{\mathit{g}}^{(t)})}{({\mathit{d}}^{(t)}{)}^{\mathrm{T}}({\mathit{g}}^{(t+\mathrm{1})}-{\mathit{g}}^{(t)}){\mathit{d}}^{(t)}}\mathrm{。}（13）

式（11）可以写为

{\mathit{X}}^{(t+\mathrm{1})}={\mathit{X}}^{(t)}-

\frac{({\mathit{g}}^{(t+\mathrm{1})}){\mathit{d}}^{(t+\mathrm{1})}}{(\mathit{A}{\mathit{d}}^{(t+\mathrm{1})}{)}^{\mathrm{T}}(\mathit{A}{\mathit{d}}^{(t+\mathrm{1})})+\beta {\mathit{d}}^{(t+\mathrm{1})}}{\mathit{d}}^{(t+\mathrm{1})}\mathrm{。}（14）

2） WGIF提取重建图像的低频分量

该阶段采用WGIF来提取图像的低频分量， 为此需求得使式（8）所示成本函数E(a_k,b_k)最小化的a_k和b_k， 由线性回归可得

a_k=\frac{\frac{\mathrm{1}}{\left|\omega \right|}\sum _{i\in {\omega }_k}{G}_i{p}_i-{\mu }_k{\overline{p}}_k}{{\sigma }_k^{\mathrm{2}}+\frac{\eta }{{\Gamma }_G(k)}},（15）

b_k={\overline{p}}_k-a_k{\mu }_k,（16）

式中： {\mu }_k和{\sigma }_k^{\mathrm{2}}为窗口{\omega }_k中像素的均值和方差； \left|\omega \right|为{\omega }_k中的像素数； {\overline{p}}_k为输入图像\mathit{p}在{\omega }_k中像素的均值。利用WGIF分离出图像{\mathit{X}}^{(t+\mathrm{1})}的低频分量{\mathit{X}}_l， 将输入图像{\mathit{X}}^{(t+\mathrm{1})}作为引导图像， 则低频分量可以表示为

{\mathit{X}}_{li}=a_k{\mathit{X}}_i^{t+\mathrm{1}}+b_k\forall i\in {\omega }_k\mathrm{。}（17）

3） 固定X， 更新特征图\left\{{\mathit{M}}_i\right\}

该阶段首先通过高通滤波获得图像的高频分量{\mathit{X}}_h， 而后结合预先训练好的滤波器\left\{{\mathit{f}}_i\right\}将式（9）模型转化为如下的优化问题，

\underset{\left\{{\mathit{M}}_i\right\}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\frac{\mathrm{1}}{\mathrm{2}}{‖{\displaystyle \sum _{i=\mathrm{1}}^N}{\mathit{f}}_i\mathrm{*}{\mathit{M}}_i-{\mathit{X}}_h‖}_{\mathrm{2}}^{\mathrm{2}}+\lambda {\displaystyle \sum _{i=\mathrm{1}}^N}{‖{\mathit{M}}_i‖}_{\mathrm{1}}+

\frac{\tau }{\mathrm{2}}{{\displaystyle \sum _{i=\mathrm{1}}^N}‖\sqrt[]{({\mathit{g}}_{\mathrm{1}}\mathrm{*}{\mathit{M}}_i{)}^{\mathrm{2}}+({\mathit{g}}_{\mathrm{2}}\mathrm{*}{\mathit{M}}_i{)}^{\mathrm{2}}}‖}_{\mathrm{2}}^{\mathrm{2}}\mathrm{。}（18）

定义线性算子\mathit{F}， {\mathit{G}}_{\mathrm{1}}， {\mathit{G}}_{\mathrm{2}}， 使得{\mathit{F}}_i{\mathit{M}}_i={\mathit{f}}_i\mathrm{*}{\mathit{M}}_i， {\mathit{G}}_l{\mathit{M}}_i={\mathit{g}}_l\mathrm{*}{\mathit{M}}_i(l=\mathrm{1,2})， 式（18）中的梯度项可以写为

\frac{\tau }{\mathrm{2}}{\displaystyle \sum _{i=\mathrm{1}}^N}{‖{\mathit{G}}_{\mathrm{1}}{\mathit{M}}_i‖}_{\mathrm{2}}^{\mathrm{2}}+\frac{\tau }{\mathrm{2}}{\displaystyle \sum _{i=\mathrm{1}}^N}{‖{\mathit{G}}_{\mathrm{2}}{\mathit{M}}_i‖}_{\mathrm{2}}^{\mathrm{2}}\mathrm{。}（19）

令\mathit{F}=({\mathit{F}}_{\mathrm{1}},{\mathit{F}}_{\mathrm{2}},\cdots ,{\mathit{F}}_N)， \mathit{M}=({\mathit{M}}_{\mathrm{1}},{\mathit{M}}_{\mathrm{2}},\cdots，{\mathit{M}}_N{)}^{\mathrm{T}}， {\mathit{\Phi }}_l=\left[\begin{array}{ccc}{\mathit{G}}_l& \mathrm{0}& \cdots \\ \mathrm{0}& {\mathit{G}}_l& \cdots \\ ⋮& ⋮& ⋮\end{array}\right]， 则式（18）可以写为

\underset{\mathit{M}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\frac{\mathrm{1}}{\mathrm{2}}{‖\mathit{F}\mathit{M}-{\mathit{X}}_h‖}_{\mathrm{2}}^{\mathrm{2}}+\lambda {‖\mathit{M}‖}_{\mathrm{1}}+

\frac{\tau }{\mathrm{2}}{‖{\mathit{\Phi }}_{\mathrm{1}}\mathit{M}‖}_{\mathrm{2}}^{\mathrm{2}}+\frac{\tau }{\mathrm{2}}{‖{\mathit{\Phi }}_{\mathrm{2}}M‖}_{\mathrm{2}}^{\mathrm{2}}\mathrm{。}（20）

应用交替方向乘子法（Alternating Direction Method Multipliers， ADMM）^［13］求解式（20）， 引入约束\mathit{B}等于变量\mathit{M}， 导致以下问题

\underset{\mathit{M},\mathit{B}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\frac{\mathrm{1}}{\mathrm{2}}{‖\mathit{F}\mathit{M}-{\mathit{X}}_h‖}_{\mathrm{2}}^{\mathrm{2}}+\frac{\tau }{\mathrm{2}}{‖{\mathit{\Phi }}_{\mathrm{1}}\mathit{M}‖}_{\mathrm{2}}^{\mathrm{2}}+

\frac{\tau }{\mathrm{2}}{‖{\mathit{\Phi }}_{\mathrm{2}}\mathit{M}‖}_{\mathrm{2}}^{\mathrm{2}}+\lambda {‖\mathit{B}‖}_{\mathrm{1}},\mathrm{s}.\mathrm{t}.\mathit{M}=\mathit{B}\mathrm{。}（21）

采用ADMM， 迭代公式为

{\mathit{M}}^{j+\mathrm{1}}=\underset{\mathit{M}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\frac{\mathrm{1}}{\mathrm{2}}{‖\mathit{F}\mathit{M}-{\mathit{X}}_h‖}_{\mathrm{2}}^{\mathrm{2}}+\frac{\tau }{\mathrm{2}}{‖{\mathit{\Phi }}_{\mathrm{1}}\mathit{M}‖}_{\mathrm{2}}^{\mathrm{2}}+\frac{\tau }{\mathrm{2}}{‖{\mathit{\Phi }}_{\mathrm{2}}\mathit{M}‖}_{\mathrm{2}}^{\mathrm{2}}+\frac{\rho }{\mathrm{2}}{‖\mathit{M}-{\mathit{B}}^j+{\mathit{C}}^j‖}_{\mathrm{2}}^{\mathrm{2}},（22）

{\mathit{B}}^{j+\mathrm{1}}=\underset{\mathit{B}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\lambda {‖\mathit{B}‖}_{\mathrm{1}}+\frac{\rho }{\mathrm{2}}{‖{\mathit{M}}^{j+\mathrm{1}}-\mathit{B}+{\mathit{C}}^j‖}_{\mathrm{2}}^{\mathrm{2}},（23）

{\mathit{C}}^{j+\mathrm{1}}={\mathit{C}}^j+{\mathit{M}}^{j+\mathrm{1}}-{\mathit{B}}^{j+\mathrm{1}}\mathrm{。}（24）

将式（22）变换到傅里叶域， 则解可表示为

({\widehat{\mathit{F}}}^H\widehat{\mathit{F}}+\tau {\widehat{\mathit{\Phi }}}_{\mathrm{1}}^H{\widehat{\mathit{\Phi }}}_{\mathrm{1}}+\tau {\widehat{\mathit{\Phi }}}_{\mathrm{2}}^H{\widehat{\mathit{\Phi }}}_{\mathrm{2}}+\rho \mathit{I})\widehat{\mathit{M}}=

{\widehat{\mathit{F}}}^H{\widehat{\mathit{X}}}_{\mathit{h}}+\rho (\widehat{\mathit{B}}-\widehat{\mathit{C}}),（25）

式中： \widehat{\mathit{F}}、 {\widehat{\mathit{\Phi }}}_{\mathrm{1}}、 {\widehat{\mathit{\Phi }}}_{\mathrm{2}}、 \widehat{\mathit{M}}、 {\widehat{\mathit{X}}}_h、 \widehat{\mathit{B}}、 \widehat{\mathit{C}}表示变换到傅里叶域后的表达式； \mathit{I}为单位矩阵。式（23）的闭式解可得

{\mathit{B}}^{j+\mathrm{1}}=S_{\lambda /\rho }({\mathit{M}}^{j+\mathrm{1}}+{\mathit{C}}^j),（26）

式中： S_{\lambda /\rho }(\cdot )表示软阈值函数^［14］， 其计算公式为

S_{\lambda /\rho }(x)=sign(x)\mathrm{m}\mathrm{a}\mathrm{x}(\mathrm{0},\left|x\right|-\lambda /\rho )\mathrm{。}（27）

4） 低频分量{\mathit{X}}_l与高频分量{\mathit{X}}_h结合获得新的重建图像\mathit{X}

该阶段首先将第3）步更新的特征图{\mathit{M}}^{j+\mathrm{1}}与预先训练好的滤波器\left\{{\mathit{f}}_i\right\}卷积求和得到新的高频分量， 令\tilde{\mathit{M}}={\mathit{M}}^{j+\mathrm{1}}， 即高频分量计算公式为

{\mathit{X}}_h={\displaystyle \sum _{i=\mathrm{1}}^N}{\mathit{f}}_i\mathrm{*}{\tilde{\mathit{M}}}_i\mathrm{。}（28）

而后将WGIF提取的低频分量{\mathit{X}}_l与高频分量{\mathit{X}}_h相加， 新的重建图像为

\mathit{X}={\mathit{X}}_l+{\mathit{X}}_h\mathrm{。}（29）

本文算法具体伪代码如下：

3 实验结果

为验证本文方法的性能， 在实验中测试了腹部和胸腔两个代表性的切片。并与PLS^［15］， PLSGSR^［7］， PLSCSC^［8］， PLS-CSCGR^［9］ 等算法进行了对比分析。采用峰值信噪比（PSNR）、 均方误差（RMSE）和特征相似性（FSIM）等评价指标进一步定量评估了各重建方法的性能。测试平台为： 64位Intel Core i5-1035G1处理器、 主机频率1.19 GHz、 8 G运行内存的Windows10操作系统， 数学软件为MATLAB 2024a。为了加速算法， 使用图形处理单元NVIDIA GeForce MX350执行CSC和CSCGR阶段。本实验所用滤波器\left\{{\mathit{f}}_i\right\}是对NIH-AAPM-Mayo Clinic的CT图像数据集的高频分量进行训练。实验中使用的所有滤波器和正则化参数都是相同的， 包括32个滤波器， 尺寸为10×10， β=0.005， λ=0.005， ρ=100λ+1和τ=0.06。

3.1 腹部图像重建结果分析

图 1 为不同方法重建的腹部图像效果比较， 显示窗口为［-150 250］HU， 其中第一行展示了原始腹部图像和不同算法重建的结果， 第二行为相应地放大了红色框区域， 可更好地看到细节信息。由图 1（b） 可看出， PLS出现了严重的条纹伪影， 模糊了细小的结构特征； 由图 1（c）， （d） 可看出， PLS-GSR， PLS-CSC有效地消除了条纹伪影， 但在某些地方看起来很平滑； 由图 1（e） 可看出， PLS-CSCGR方法保留了很多细节， 但仍存在噪声； 由图 1（f） 可看出， 本文算法与其他算法相比， 保留了更多细节， 视觉上纹理清晰， 边缘丰富。所以， 本文算法重建的结果最接近参考图像。

由于通过视觉效果评估算法的重建效果有一定的主观性， 进一步对比分析了各重建算法在迭代过程中PSNR、 MSE、 FISM的变化情况， 如图 2 所示。

从图 2（a）， （c） 可看出， 本文算法重建结果的PSNR和FSIM值相比于其他算法均有一定的提升， 表明该算法重建的图像与原始图像最为接近。在图 2（b） 中， 本文算法的MSE值在迭代过程中低于其他算法， 意味着其重建图像与原始图像相差最小， 具有良好的重建性能， 进一步证明了本文所提算法的有效性。

3.2 胸腔图像重建结果分析

图 3 为不同方法重建的胸腔图像效果比较， 显示窗口为［-1 000 250］HU， 其中， 第一行展示了原始胸腔图像与不同方法重建的结果； 第二行为放大的感兴趣区域（红色框区域）， 可更好地展示可视化细节。PLS重建结果（见图3（b））存在严重条纹伪影和噪声， 丢失了很多细节信息； PLS-GSR重建结果（见图3（c））消除了条状伪影， 是引入了振铃伪影， 这可能是由于字典原子不准确引起的； PLS-CSC与PLS-CSCGR重建结果（见图3（d）， （e））消除了大量伪影， 保留了较完整的细节， 但边缘有些平滑， 结构细节模糊； 本文算法（见图 3（f））边缘比较清晰， 保留了细微结构信息， 边缘保留效果优于比较算法。图 4 展示了在迭代过程中各算法重建胸腔图像的PSNR、 MSE、 FISM值的变化情况。

从图 4（a）~（c） 可清晰地看到在每次迭代过程中， 本文算法具有最高的PSNR和FSIM， 以及最低的MSE。并且在迭代开始前期， 本文算法已达到了较好的性能， 进一步说明了本文算法具有良好的重建优势。

表 1， 表 2 分别给出了不同算法重建腹部图像和胸腔图像的PSNR， MSE以及FSIM的数据。其中， 本文算法与PLS， PLS-GSR， PLS-CSC， PLS-CSCGR相比， 重建腹部图像的PSNR分别提高了7.132 0， 3.207 1， 3.098 9和1.551 6 dB； 重建胸腔图像的PSNR值分别提高了8.625 2， 3.324 8， 3.771和1.810 8 dB。同时， 本文算法重建两幅图像得到的MSE最低， 表明该算法在减少图像误差方面表现最突出。并且本文算法的FSIM最接近1， 表明该算法重建的图像与原始图像有更高的相似性。由此可看出本文算法在重建图像过程中还原了更多结构信息。

4 结　论

本文提出了一种基于WGIF与CSCGR结合的CT重建算法。该算法通过CSC处理高频信息， 保留了重要的图像细节特征， 利用WGIF获得低频信息， 在消除噪声和伪影的同时很好地保留了图像的组织结构和边缘信息。因此， 在CT图像重建实验中， 该算法无论在具体的数值评价指标上还是在重建结果的直观表现上， 相较于对比算法都体现出了更优的算法性能。

参考文献

[1] 周少飞， 柴锐， 秦品乐， 等. 基于Mask RCNN和U-Net结合的三阶段肾脏与肿瘤分割方法［J］. 中北大学学报（自然科学版）， 2022， 43（3）： 236-243.

[2] ZHOU Shaofei， CHAI Rui， QIN Pinle， et al. A three-stage kidney and tumor segmentation method based on the combination of mask RCNN and U-net［J］. Journal of North University of China （Natural Science Edition）， 2022， 43（3）： 236-243. （in Chinese）

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

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

[5] QU Z， ZHAO X， PAN J， et al. Sparse-view CT reconstruction based on gradient directional total variation［J］. Measurement Science and Technology， 2019， 30（5）： 055404.

[6] PAN Z， LIU Z， LUO K， et al. Compressive sensing total-variation primal-dual algorithms for image reconstruction［J］. IEEE Signal Processing Letters， 2024， 31： 1965-1969.

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

[8] BAO P， ZHOU J， 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）， 2018： 5097-5100.

[9] BRISTOW H， ERIKSSON A， LUCEY S. Fast convolutional sparse coding［C］//2013 IEEE Conference on Computer Vision and Pattern Recognition， 2013： 391-398.

[10] BAO P， XIA W， YANG K， et al. Sparse-view ct reconstruction via convolutional sparse coding［C］//2019 IEEE 16th International Symposium on Biomedical Imaging （ISBI 2019）， 2019： 1446-1449.

[11] LI Z， ZHENG J， ZHU Z， et al. Weighted guided image filtering［J］. IEEE Transactions on Image Processing， 2014， 24（1）： 120-129.

[12] 李雨， 史娜， 孔慧华， 等. 基于全变分和梯度域卷积稀疏编码的稀疏角度CT重建算法［J］. 激光与光电子学进展， 2021， 58（12）： 339-348.

[13] LI Yu， SHI Na， KONG Huihua， et al. Sparse angle CT reconstruction algorithm based on total variation and convolutional sparse coding in gradient domain［J］. Laser & Optoelectronics Progress， 2021， 58（12）： 339-348. （in Chinese）

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

[15] BOYD S， PARIKH N， CHU E， et al. Distributed optimization and statistical learning via the alternating direction method of multipliers［J］. Foundations and Trends® in Machine Learning， 2011， 3（1）： 1-122.

[16] BIOUCAS-DIAS J M， FIGUEIREDO M A T. A new twIst： two-step iterative shrinkage/thresholding algorithms for image restoration［J］. IEEE Transactions on Image Processing， 2007， 16（12）： 2992-3004.

[17] ROGGEMANN M C， TYLER D W. Model-based image reconstruction by means of a constrained least-squares solution［J］. Applied Optics， 1997， 36（11）： 2360.

基金资助

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

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

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

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

山西省回国留学人员科研项目(2021-108)