基于亚像素各向异性扩散的低剂量CT重建方法

A low-dose CT reconstruction method using sub-pixel anisotropic diffusion

摘要

目的 提出一种基于亚像素各项异性扩散的低剂量CT重建方法。 方法 通过线性插值技术计算亚像素单元强度值及其二阶差分后，将计算得到的新的梯度信息引入到各项异性扩散过程中，并结合惩罚加权最小二乘模型对低剂量CT投影数据进行滤波，最后使用滤波反投影算法将恢复后的投影数据重建出CT图像。 结果 在Shepp-Logan体模实验中，与FBP、PWLS-Gibbs和PWLS-TV方法相比，新方法滤波后重建的CT图像在结构相似指数上分别提升了28.13%、5.49%和0.91%，在特征相似指数上分别提升了21.08%、1.78%和1.36%，并且在均方根误差上分别降低了69.59%、18.96%和3.90%。在XCAT体模实验中，与FBP、PWLS-Gibbs和PWLS-TV方法相比，新方法在结构相似指数上分别提高了14.24%、1.43%及7.89%，在特征相似指数上分别提高了9.61%、1.78%及5.66%，同时在均方根误差上分别降低了26.88%、9.41%及18.39%。在临床数据实验中，与FBP、PWLS-Gibbs和PWLS-TV方法重建的CT图像相比，新方法在结构相似指数上分别提升了19.24%、15.63%和3.68%，在特征相似指数上分别提升了4.30%、2.92%和0.43%，同时在均方根误差上分别降低了44.60%、36.84%和15.22%，并且在峰值信噪比上提升至28.39。 结论 本文提出的新方法可以有效去除低剂量CT图像的噪声和伪影，并可以保持结构细节信息。

Abstract

Objective We present a new low-dose CT reconstruction method using sub-pixel and anisotropic diffusion. Methods The sub-pixel intensity values and their second-order differences were obtained using linear interpolation techniques, and the new gradient information was then embedded into an anisotropic diffusion process, which was introduced into a penalty-weighted least squares model to reduce the noise in low-dose CT projection data. The high-quality CT image was finally reconstructed using the classical filtered back-projection (FBP) algorithm from the estimated data. Results In the Shepp-Logan phantom experiments, the structural similarity (SSIM) index of the CT image reconstructed by the proposed algorithm, as compared with FBP, PWLS-Gibbs and PWLS-TV algorithms, was increased by 28.13%, 5.49%, and 0.91%, the feature similarity (FSIM) index was increased by 21.08%, 1.78%, and 1.36%, and the root mean square error (RMSE) was reduced by 69.59%, 18.96%, and 3.90%, respectively. In the digital XCAT phantom experiments, the SSIM index of the CT image reconstructed by the proposed algorithm, as compared with FBP, PWLS-Gibbs and PWLS-TV algorithms, was increased by 14.24%, 1.43% and 7.89%, the FSIM index was increased by 9.61%, 1.78% and 5.66%, and the RMSE was reduced by 26.88%, 9.41% and 18.39%, respectively. In clinical experiments, the SSIM index of the image reconstructed using the proposed algorithm was increased by 19.24%, 15.63% and 3.68%, the FSIM index was increased by 4.30%, 2.92% and 0.43%, and the RMSE was reduced by 44.60%, 36.84% and 15.22% in comparison with FBP, PWLS-Gibbs and PWLS-TV algorithms, respectively. Conclusion The proposed method can effectively reduce the noises and artifacts while maintaining the structural details in low-dose CT images.

Graphical abstract

关键词

低剂量CT / 各向异性扩散 / 亚像素 / 图像重建

Key words

low-dose computed tomography / anisotropic diffusion / sub-pixel / image reconstruction

计算机断层成像（CT）因其在时间、空间和密度分辨力的优越表现，已广泛用于临床诊断和治疗^［1］。CT图像的质量与X射线辐射剂量密切相关，剂量越高图像质量越好，然而，过量的X射线照射会引发肿瘤、白血病以及其他遗传性疾病^{［2， 3］}。如何降低CT扫描过程中的剂量，并提供高质量的CT重建图像，是目前低剂量CT在临床应用的困难之一。降低管电流或者管电压是实现低剂量CT（LDCT）最简单且最有效的方法^［4］。但是，降低管电流或者管电压会导致投影数据中光子噪声和电子噪声大幅增加而造成传统滤波反投影（FBP）算法重建图像的质量严重退化。因此，在保证CT图像质量的前提下，最大限度地降低X射线辐射剂量已经成为今后CT发展和应用的一个重要努力方向。

在目前的低剂量CT成像技术中，存在3种主要方法^［5-8］：深度学习方法，统计迭代重建（SIR）方法和基于低剂量CT投影数据的图像恢复方法。近些年来深度学习的迅猛发展，基于深度学习的图像恢复方法获得广泛研究^{［6， 9-11］}。但是深度学习的不可解释性导致它在医学图像处理中的可信度存在很大的质疑，同时深度学习方法需要庞大的数据训练量，这些数据通常是需要标签的，而获得含标签的临床数据CT图像往往需要人工标注，存在较大困难。SIR方法通过建立CT系统模型和优化目标函数来实现图像重建，目标函数通常含有两项，即数据保真项和正则化项^［5］。其中的数据保真项用来度量估计值和测量值之间的差异，正则化项用来吸收先验信息，有助于提升图像平滑度和精度，但SIR方法需要在投影数据和图像之间反复进行投影与反投影，计算时间大幅增加^［12-14］。基于投影数据的图像恢复方法首先对投影数据进行滤波，然后使用经典的滤波反投影算法进行图像重建。由于考虑到投影数据的噪声统计特性，这类方法具有高效且重建图像均匀性好等优点^［7］。Wang等^{［15， 16］}通过实验以及理论分析证明了投影数据经对数变换后近似满足高斯分布，并根据该噪声统计特性提出了多种方法来恢复低剂量CT投影数据正弦图。例如惩罚加权最小二乘方法^［15］以及利用小波变换的多尺度惩罚加权最小二乘方法^［16］。上述研究方法强调了依据噪声统计特性等先验信息对低剂量CT投影数据进行恢复的重要性，说明基于投影数据的图像恢复方法能够有效降低噪声和伪影在低剂量CT图像中的影响，从而实现优质的低剂量CT图像重建。

各向异性扩散利用图像的梯度信息调整像素值，以平滑图像同时保留其边缘和细节。然而，在实际应用中，多次迭代处理可能导致过度平滑现象^［17］。亚像素技术则通过插值分析周围像素的数据，在不增加实际像素的情况下提高图像的清晰度和细节。基于此，为了解决低剂量CT重建图像存在严重的噪声、伪影等问题，同时提升CT重建图像的分辨率，并进一步研究各向异性扩散在低剂量CT处理中的应用前景。本研究将亚像素技术和各向异性扩散结合^［18］，将其作为先验信息引入低剂量CT投影数据的恢复中，提出了一种基于亚像素各向异性扩散的惩罚加权最小二乘低剂量CT重建方法。该方法结合亚像素技术和各向异性扩散的优势，根据投影数据的不同方向特征进行精确的插值，得到新的梯度信息后，使用各向异性扩散对投影数据进行滤波，最后采用FBP算法进行CT图像重建。实验结果表明，该方法能有效去除噪声和伪影，同时提高低剂量CT图像的质量和可视效果。

1 材料和方法

1.1 PWLS模型

低剂量CT扫描下，经对数变换和校准等步骤后的投影数据中的噪声近似满足高斯分布，其噪声统计特性如下^［19］：

\begin{array}{l}{\sigma }_i^{\mathrm{2}}=\frac{\mathrm{1}}{I_{\mathrm{0}}}exp\left({\overline{y}}_i\right)\left(\mathrm{1}+\frac{\mathrm{1}}{I_{\mathrm{0}}}exp\left({\overline{y}}_i\left({\sigma }_e^{\mathrm{2}}-\mathrm{1.25}\right)\right)\right),\\ i=\mathrm{1,2},\cdots ,M\end{array}（1）

其中，I_{\mathrm{0}}表示入射光子总数，{\overline{y}}_i表示第i号探测器上测得的投影数据均值，{\sigma }_e^{\mathrm{2}}表示投影数据中的背景电子噪声方差，M代表探测器个数。

根据上述噪声统计特性，设y={\left(y_{\mathrm{1}},y_{\mathrm{2}},\cdots ,y_M\right)}^T为待恢复的投影数据，q={\left(q_{\mathrm{1}},q_{\mathrm{2}},\cdots ,q_M\right)}^T为投影数据滤波后的理想值，则PWLS投影数据恢复模型可以写为^［20］：

\underset{q\ge \mathrm{0}}{\mathrm{m}\mathrm{i}\mathrm{n}}{\left(y-q\right)}^T{G}^{-\mathrm{1}}\left(y-q\right)+\beta R\left(q\right)（2）

其中，\mathbf{G}=\frac{\mathrm{1}}{\alpha }\mathbf{I}+\mathbf{\Sigma }，\alpha >\mathrm{0}为超参数，\mathrm{\Sigma }=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left\{{\sigma }_{\mathrm{1}}^{\mathrm{2}},{\sigma }_{\mathrm{2}}^{\mathrm{2}},\cdots ,{\sigma }_M^{\mathrm{2}}\right\}是对角矩阵，其对角线上的第i号元素是由（1）式计算得到的投影数据的方差{\sigma }_i^{\mathrm{2}}。\beta >\mathrm{0}为权衡数据保真项和正则化项的超参数。

1.2 各向异性扩散

各项异性扩散在不同方向上的扩散过程，可以描述成能量泛函的形式，并通过寻找能量泛函的最小值来求解。能量泛函的表达式如下^［21］：

\underset{q}{\mathrm{m}\mathrm{i}\mathrm{n}}R(q)={\int }_{\Omega }f(|\nabla q|)d\mathrm{\Omega }（3）

其中，\mathrm{\Omega }表示当前图像空间，\nabla q=({\nabla }_{x}q,{\nabla }_{y}q)是梯度算子，\left|\nabla q\right|=\sqrt[]{{\nabla }_x^{\mathrm{2}}q+{\nabla }_y^{\mathrm{2}}q}。f(\cdot )是扩散系数函数：

f(\omega )=exp(-{\omega }^{\mathrm{2}}/{\epsilon }^{\mathrm{2}})（4）

其中\epsilon是控制扩散强度的一个超参数。能量泛函的最优解可以通过如下的欧拉方程求解^［18］：

\frac{\partial R(q)}{\partial t}=-div[c(\left|\nabla q\right|)\nabla q]=-{\nabla }_x[c({\left|\nabla q\right|}^{\mathrm{2}}){\nabla }_{x}q]-{\nabla }_y[c({\left|\nabla q\right|}^{\mathrm{2}}){\nabla }_{y}q]（5）

其中，t表示当前时间步，div表示散度。c(\cdot )是扩散系数，其表达式如下：

c(\omega )=\frac{f^{\mathrm{\text{'}}}(\sqrt[]{\omega })}{\sqrt[]{\omega }}（6）

其中，符号'代表一阶导数。

1.3 亚像素

假设待恢复的投影数据y大小为m\times n，m为投影数据的采样角数目，n为探测器个数，设定一个亚像素距离h_s，将投影数据中相邻像素之间的像素单元作为亚像素\mathrm{s}，并通过插值的方式获得亚像素的单元强度值。例如，投影数据中有两个相邻的像素单元y_{i,j}和y_{i+\mathrm{1},j}（i和j分别为像素单元在投影数据中的索引，i=\mathrm{0,1},\mathrm{2}\dots m，j=\mathrm{0,1},\mathrm{2}\dots n）（图1）。

两种不同条件下的亚像素单元表示，亚像素距离小于1（图1A），亚像素距离大于1（图1B）。然后采用线性插值方法得到亚像素的强度值，该过程可以表示为^［18］：

s_{i+h_s,j}={\mathrm{y}}_{i,j}-h_s({\mathrm{y}}_{i,j}-{\mathrm{y}}_{i+\mathrm{1},j})（7）

s_{i-h_s,j}={\mathrm{y}}_{i,j}-h_s({\mathrm{y}}_{i,j}-{\mathrm{y}}_{i-\mathrm{1},j})（8）

s_{i,j+h_s}={\mathrm{y}}_{i,j}-h_s({\mathrm{y}}_{i,j}-{\mathrm{y}}_{i,j+\mathrm{1}})（9）

s_{i,j-h_s}={\mathrm{y}}_{i,j}-h_s({\mathrm{y}}_{i,j}-{\mathrm{y}}_{i,j-\mathrm{1}})（10）

其中，y_{i,j}表示投影数据y在(i,j)处的值。y_{i+\mathrm{1},j}，y_{i-\mathrm{1},j}，y_{i,j+\mathrm{1}}和y_{i,j-\mathrm{1}}是y_{i,j}的相邻像素，s_{i+h_s,j}，s_{i-h_s,j}，s_{i,j+h_s}和s_{i,j-h_s}是y_{i,j}的相邻亚像素。特别地，对于边缘处的亚像素，由边界上的投影数据代替，即s_{-h_s,j}=y_{\mathrm{0},j}，s_{m-\mathrm{1}+h_s,j}=y_{m-\mathrm{1},j}，s_{i,-h_s}=y_{i,\mathrm{0}}以及s_{i,n-\mathrm{1}+h_s}=y_{i,n-\mathrm{1}}。

1.4 亚像素各项异性扩散

通过一个离散的拉普拉斯算子{\nabla }_{S}q=({\nabla }_{S_x}q,{\nabla }_{S_y}q)，将亚像素的二阶差分{\nabla }_{S}q作为原先梯度\nabla q的一种泛化形式，过程如下：

{\nabla }_{S_x}{q}_{i,j}=\frac{s_{i+h_s,j}+s_{i-h_s,j}-\mathrm{2}{y}_{i,j}}{h_s^{\mathrm{2}}}（11）

{\nabla }_{S_y}{q}_{i,j}=\frac{s_{i,j+h_s}+s_{i,j-h_s}-\mathrm{2}{y}_{i,j}}{h_s^{\mathrm{2}}}（12）

{\nabla }_S{q}_{i,j}={\nabla }_{S_x}{q}_{i,j}+{\nabla }_{S_y}{q}_{i,j}=\frac{s_{i+h_s,j}+s_{i-h_s,j}+s_{i,j+h_s}+s_{i,j-h_s}-\mathrm{4}{y}_{i,j}}{h_s^{\mathrm{2}}}（13）

其中，{\nabla }_{S_x}{q}_{i,j}和{\nabla }_{S_y}{q}_{i,j}分别是亚像素q_{i,j}在x和y方向的二阶差分。将得到的二阶差分替代式（5）中的\nabla q，进一步，式（5）可写成：

\frac{\partial R(q)}{\partial t}=-{\nabla }_x[c({\left|{\nabla }_{S}q\right|}^{\mathrm{2}}){\nabla }_{S_x}q]-{\nabla }_y[c({\left|{\nabla }_{S}q\right|}^{\mathrm{2}}){\nabla }_{S_{\mathrm{y}}}q]（14）

然后采用迭代策略求解式（14），该过程可以通过如下偏微分方程表示：

\begin{array}{l}{q}^{t+\mathrm{1}}=q^t-\lambda \left({\nabla }_x[c({\left|{\nabla }_S{q}^t\right|}^{\mathrm{2}}){\nabla }_{S_x}{q}^t]+{\nabla }_y[c({\left|{\nabla }_S{q}^t\right|}^{\mathrm{2}}){\nabla }_{S_{\mathrm{y}}}{q}^t]\right),\\ t=\mathrm{0,1},\mathrm{2},\cdots ,\widehat{T}\end{array}（15）

其中，q^{t+\mathrm{1}}表示第t+\mathrm{1}次迭代结果，q^t表示第t次迭代结果，\lambda为时间步长，\widehat{T}是最大迭代次数。

1.5 PWLS-SPAD低剂量CT重建

基于1.1~1.4节，将结合亚像素的各向异性扩散作为正则化项，引入到PWLS模型中，称为基于亚像素各向异性扩散的惩罚加权最小二乘投影数据恢复模型（PWLS-SPAD）。特别地，式（2）中的R\left(q\right)表示引入亚像素各项异性扩散的惩罚项，其具体执行步骤如1.2-1.4节所示。进一步，引入一个新的变量p，采用交替方向优化算法求解式（2）中的PWLS-SPAD投影数据恢复模型：

p^{k+\mathrm{1}}=\mathrm{a}\mathrm{r}\mathrm{g}\underset{p\ge \mathrm{0}}{\mathrm{m}\mathrm{i}\mathrm{n}}{\left(y-p\right)}^T{\mathrm{\Sigma }}^{-\mathrm{1}}\left(y-p\right)+\alpha {‖p-q^k‖}_{\mathrm{2}}^{\mathrm{2}}（16）

q^{k+\mathrm{1}}=\mathrm{a}\mathrm{r}\mathrm{g}\underset{q\ge \mathrm{0}}{\mathrm{m}\mathrm{i}\mathrm{n}}\alpha {‖p^{k+\mathrm{1}}-q‖}_{\mathrm{2}}^{\mathrm{2}}+\beta R\left(q\right)（17）

明显地，第一个子问题式（16）中的目标函数是严格的凸函数，导数设置为0即可得到如下解：

p^{k+\mathrm{1}}=\frac{y+\alpha {\sigma }^{\mathrm{2}}{q}^k}{\mathrm{1}+\alpha {\sigma }^{\mathrm{2}}}（18）

但这个解是由数据保真项求出，一般含有噪声。将这个解作为初值去优化式（17）中的目标函数，可以得到一个更好的结果。第二个子问题式（17）中的目标函数采用梯度下降的方法求解，公式如下：

q^{t+\mathrm{1}}=q^t-\lambda R^{\text{'}}(q^t)（19）

其中，q^{t+\mathrm{1}}和q^t分别是第t+\mathrm{1}和第t次迭代的值，当式（19）达到最大迭代次数\widehat{T}时，此时的q^{t+\mathrm{1}}就是式（17）的解，即q^{k+\mathrm{1}}=q^{t+\mathrm{1}}。\lambda是时间步长，R^{\text{'}}(q^t)表示R(q^t)的一阶导数，结合亚像素各向异性扩散进行计算，公式如下：

R^{\text{'}}(q^t)=\alpha (p^{k+\mathrm{1}}-q^t)+\beta (\nabla R(q^t))（20）

其中，p^{k+\mathrm{1}}为式（18）中得到的解。\nabla R(q^t)是亚像素各向异性扩散的结果，其数值解通过式（14）-（15）计算得到，即\nabla R(q^t)=-{\nabla }_x[c({\left|{\nabla }_{S}q\right|}^{\mathrm{2}}){\nabla }_{S_x}q]-{\nabla }_y[c({\left|{\nabla }_{S}q\right|}^{\mathrm{2}}){\nabla }_{S_{\mathrm{y}}}q]。注意到，交替方向优化算法需要给定初始点q^k=q^{\mathrm{0}}，之后产生的迭代序列为：p^{\mathrm{1}},q^{\mathrm{1}},p^{\mathrm{2}},q^{\mathrm{2}},\cdots ,p^k,q^k,\cdots，直至满足迭代终止条件。

综上所述，基于亚像素各向异性扩散的低剂量CT重建方法执行过程如下所示：

初始化：对参数h_s，\epsilon，\lambda，\widehat{T}，\beta，\alpha以及变量q^{\mathrm{0}}初始化。

步骤1：q^{\mathrm{0}}初始化后代入式（17），由式（19）解得p^{k+\mathrm{1}}；步骤2：通过式（7）~（10）获得亚像素单元强度值；步骤3：通过式（11）~（15）将二阶差分引入到各向异性扩散过程中；步骤4：通过式（20）~式（21）对p^{k+\mathrm{1}}进行滤波，得到q^{k+\mathrm{1}}；步骤5：重复执行步骤1~4直至满足迭代终止条件；步骤6：对步骤5中得到的投影数据q^{k+\mathrm{1}}通过FBP算法得到重建后的CT图像。

PWLS-SPAD方法中有6个参数(h_s,\epsilon ,\lambda ,\widehat{T},\beta ,\alpha )需要初始化。变量q^{\mathrm{0}}的初始值是测量得到的低剂量CT投影数据y。h_s根据文献^［18］中的方法进行设置。\epsilon设置为1.0×10^-2来保留图像边缘。参数\lambda是PDE时间步长，设置为1.0×10^-3确保迭代的收敛。各向异性扩散的最大迭代次数\widehat{T}设置为9。实验结果表明，\alpha对重建结果影响不明显，本文中选取\alpha =\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{3}}。参数\beta是用来平衡数据保真项和惩罚项的影响^［22］，\beta过大会导致图像模糊，过小会导致图像噪声没有很好地去除，在本文中根据噪声水平大小和视觉效果，对不同的实验数据采用试错的方法选取\beta的取值。在Shepp-Logan体模实验中，\beta设置为1×10^-4。在XCAT体模实验中，\beta设置为1.0×10^-3。在临床数据实验中，\beta设置为5×10⁴。在重建过程中，本文使用的迭代终止条件设置为：

\frac{||q^{k+\mathrm{1}}-q^k|{|}_{\mathrm{2}}}{||q^k|{|}_{\mathrm{2}}}\le \mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{3}}（21）

或迭代次数\ge \mathrm{50}。

1.6 实验

1.6.1 实验数据

采用了仿真Shepp-Logan体模（图2A）、XCAT体模（图2B）和真实腹部临床数据进行实验分析。根据文献^［23］中的方法对Shepp-Logan体模进行仿真，CT采样几何为使用弧形探测器的扇形束，其他相关的实验数据如下：X射线源至旋转中心和探测器的距离分别是615.18 mm和1361.2 mm；旋转角在全角度下采样值为1160，每个采样角对应672个探测器，其大小为1.85 mm；（3）模拟的入射光子数目为5.0×10⁴；（4）生成的体模图像大小为512×512，像素大小是1.0 mm×1.0 mm。XCAT体模仿真使用相同的采样几何，其他相关的实验数据如下：X射线源至旋转中心和探测器的距离分别是949 mm和408 mm；旋转角在全角度下采样值为984，每个采样角对应888个探测器，其大小为1.0 mm；模拟的入射光子数目为1.0×10⁵；生成的体模图像大小为512×512，像素大小是1.0 mm×1.0 mm。临床数据由商用CT设备采集，在两种不同的扫描协议下进行CT数据采集。一是在120 kVp球管电压和50 mAs的扫描协议下获取的低剂量CT数据，二是在120 kVp球管电压和400 mAs的扫描协议下获取的正常剂量CT数据，将其重建的图像视作金标准进行比较，其他相关的实验数据如下：X射线源至旋转中心和探测器的距离分别是570 mm和1040 mm；旋转角在全角度下采样值为1160，每个采样角对应672个探测器，其大小为1.40 mm；生成的CT图像大小为512×512，像素大小是0.6 mm×0.6 mm。

1.6.2 对比方法

为了验证PWLS-SPAD方法的有效性，将其与经典的FBP算法，基于吉布斯先验约束的惩罚加权最小二乘（PWLS-Gibbs）^{［7， 24］}方法和基于变分先验约束的惩罚加权最小二乘（PWLS-TV）^{［20， 25］}方法进行比较。

PWLS-Gibbs重建方法的目标函数^［7］为：

\underset{q\ge \mathrm{0}}{\mathrm{m}\mathrm{i}\mathrm{n}}{(y-q)}^T{\mathrm{\Sigma }}^{-\mathrm{1}}(y-q)+{\beta }_{\mathrm{1}}\underset{i}{\Sigma }\underset{m\in N_i}{\Sigma }{w}_{im}{\left(q_i-q_m\right)}^{\mathrm{2}}（22）

其中，{\beta }_{\mathrm{1}}>\mathrm{0}是超参数，N_i表示在投影数据中第i个索引点的一阶邻域信息，w_{im}是在不同方向的权重，一般地，在水平和垂直方向分别为1和0.25。

PWLS-TV方法是基于变分的二阶非线性PDE的应用，本质上也是一种基于能量泛函的去噪方法，通过最小化图像的总变分来去除噪声，其目标函数^［20］为：

\underset{q\ge \mathrm{0}}{\mathrm{m}\mathrm{i}\mathrm{n}}{\left(y-q\right)}^T{\mathbf{G}}^{-\mathrm{1}}\left(y-q\right)+{\beta }_{\mathrm{2}}\mathrm{T}\mathrm{V}\left(q\right)（23）

其中，\mathbf{G}=\frac{\mathrm{1}}{{\beta }_{\mathrm{1}}}\mathbf{I}+\mathbf{\Sigma }，{\beta }_{\mathrm{1}},{\beta }_{\mathrm{2}}>\mathrm{0}是超参数，\mathrm{T}\mathrm{V}\left(q\right)表示q的全变分{‖q‖}_{\mathrm{T}\mathrm{V}}=\sqrt[]{({\nabla }_x{q)}^{\mathrm{2}}+({\nabla }_y{q)}^{\mathrm{2}}}。

1.6.3 评估指标

实验中采用均方根误差（RMSE）^［20］，峰值信噪比（PSNR）^［20］，结构相似性指数（SSIM）^［26］和特征相似度指数（FSIM）^［27］对重建图像进行定量分析。

RMSE的计算公式如下：

\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\left(\frac{{\mathrm{\Sigma }}_{i,j}{\left(\widehat{I}\left(i,j\right)-I\left(i,j\right)\right)}^{\mathrm{2}}}{{\mathrm{\Sigma }}_{i,j}{\left(I\left(i,j\right)\right)}^{\mathrm{2}}}\right)（24）

其中，I(i,j)表示真实体模图像或金标准图像在(i,j)处的强度值，\widehat{I}(i,j)表示重建图像在(i,j)处的强度值。

PSNR的计算公式如下：

\mathrm{P}\mathrm{S}\mathrm{N}\mathrm{R}=\mathrm{10}\mathrm{l}\mathrm{g}\left(\frac{\mathrm{m}\mathrm{a}{\mathrm{x}}^{\mathrm{2}}I(i,j)}{{\sum }_{i,j}(\widehat{I}(i,j)-I{(i,j))}^{\mathrm{2}}/N-\mathrm{1}}\right)（25）

其中，N表示总的像素数目。

SSIM通过亮度、对比度和结构3个方面来度量两幅图像之间的相似性，根据文献［25］进行计算。FSIM考虑了图像的相位一致性（PC）和梯度幅度（GM），可以有效评估图像的特征信息，其计算方法参考文献［27］。

2 结果

2.1 Shepp-Logan体模实验

图3给出了不同方法重建的Shepp-Logan体模图像，低剂量CT投影数据使用FBP算法直接重建得到CT图像（图3A），PWLS-Gibbs方法（图3B）和PWLS-TV滤波后重建的图像（图3C），PWLS-SPAD方法滤波后重建的图像（图3D）。采用FBP直接重建的图像中存在着严重的噪声和伪影，PWLS-Gibbs方法重建的图像去除了部分噪声，但依旧存在较为严重的噪声并且损失了一定的清晰度。PWLS-TV方法重建的图像在图像噪声方面大幅降低，却引入了明显的条纹伪影，在图像四周也可以清楚地观察到阶梯效应引起的环状纹理（图3C中红色箭头）。

与FBP、PWLS-Gibbs和PWLS-TV方法相比，PWLS-SPAD滤波后重建的CT图像在SSIM上分别提升了28.13%、5.49%和0.91%，在FSIM上分别提升了21.08%、1.78%和1.36%，并且在RMSE上分别降低了69.59%、18.96%和3.90%，即PWLS-SPAD方法重建的图像在各项评估指标上有最佳的数值表现（表1）。

不同方法重建结果中感兴趣区域（ROI）的放大图（图2中红色方框），FBP、PWLS-Gibbs、PWLS-TV和PWLS-SPAD方法重建结果中的ROI放大图（图4）。PWLS-Gibbs方法中的噪声没有得到很好地抑制，PWLS-TV方法存在严重的条状伪影，而PWLS-SPAD方法在确保去噪效果的同时表现出了清晰的视觉效果，并且在SSIM和FSIM评估指标上具有最高的数值表现（表2）。

2.2 XCAT体模实验

PWLS-Gibbs方法重建的图像中噪声没有完全去除的同时丢失了清晰度，PWLS-TV方法由于梯度信息的不准确出现了阶梯效应，导致CT图像的两侧在扩散过程中出现了严重的波状纹理（图5C），而PWLS-SPAD方法在有效去除噪声的同时，更好地保留了图像边缘细节，且提供了清晰的CT成像效果。与FBP、PWLS-Gibbs、PWLS-TV方法重建的图像相比，PWLS-SPAD方法在SSIM上分别提高了14.24%、1.43%和7.89%，在FSIM上分别提高了9.61%、1.78%和5.66%，同时在RMSE上分别降低了26.88%、9.41%和18.39%（表3）。

PWLS-Gibbs方法重建图像中的感兴趣区域依然存在一定程度的噪声以及明显的条纹伪影，而PWLS-TV方法重建图像中虽然去除了伪影，但是存在严重的噪声（图6）。PWLS-SPAD方法在去除噪声和伪影以及提高分辨率上实现了更好地均衡。在数值分析上， PWLS-SPAD在结构相似性和特征相似性上都具有更优异的表现（表4）。

2.3 临床数据实验

图7给出了正常剂量下根据临床数据重建的CT图像和低剂量下不同方法重建的CT图像。图7A是正常剂量下得到的CT图像，将其作为金标准进行实验对比。图7B是低剂量下采用FBP算法重建的CT图像，可以看出它存在严重的噪声。采用PWLS-Gibbs（图7C）、PWLS-TV（图7D）和PWLS-SPAD（图7E）方法滤波后重建的CT图像。PWLS-Gibbs方法重建得到的图像严重丢失了清晰度，PWLS-TV方法重建的图像在胸椎周围产生了严重的条纹伪影，同时在图像四周也可以明显地观察到由阶梯效应引起的环状纹理（图7D中红色箭头）。PWLS-SPAD方法滤波后重建的图像较好地抑制了噪声，且提供了均匀、清晰的图像视觉效果。

与FBP、PWLS-Gibbs和PWLS-TV方法重建的CT图像相比，PWLS-SPAD方法在SSIM上分别提升了19.24%、15.63%和3.68%，在FSIM上分别提升了4.30%、2.92%和0.43%，同时在RMSE上分别降低了44.60%、36.84%和15.22%，并且在PSNR上也有明显的数值提升（表5）。

PWLS-Gibbs方法重建的图像损失了清晰度，降低了图像质量，PWLS-TV方法虽然保证了清晰度，但是存在着明显的条状伪影，图像质量依然无法得到有效提升（图8）。PWLS-SPAD方法在去噪和图像清晰度上有更加突出的表现，相似性指数说明PWLS-SPAD方法重建结果中的ROI相似度与金标准最为接近，即PWLS-SPAD方法有效提升了图像质量（表6）。

3 讨论

CT成像质量与X射线辐射剂量息息相关，如何在保证图像质量的前提下最大幅度地降低CT扫描过程中的X射线剂量是当前CT领域的研究热点之一^{［10， 28， 29］}。为提高低剂量CT图像质量，受到亚像素各向异性扩散^［18］的启发，本文提出了一种基于亚像素各向异性扩散的低剂量CT图像重建方法，将其称为PWLS-SPAD，该方法通过线性插值技术和PDE方程将结合了亚像素的各向异性扩散引入到PWLS模型中，并且使用数值仿真Shepp-Logan体模、XCAT体模及真实临床数据进行实验验证提出方法的有效性。

在Shepp-Logan体模实验中，FBP方法重建的图像存在很严重的噪声，PWLS-Gibbs方法和PWLS-TV方法重建的结果存在明显的伪影，并且噪声也没有很好地去除，同时PWLS-TV方法重建图像存在明显的阶梯效应。提出的PWLS-SPAD方法不仅在去除噪声和伪影上有上佳表现，同时遏制了阶梯效应带来的影响。XCAT体模实验中，FBP方法、PWLS-Gibbs方法和PWLS-TV方法都未能很好地滤除噪声，并且PWLS-TV方法重建的CT图像受到了阶梯效应的干扰，导致了波状纹理的产生，严重影响了CT图像质量。PWLS-SPAD方法较对比方法提供了更优质的成像质量。临床数据实验结果中可以看出，相比FBP、PWLS-Gibbs和PWLS-TV方法，PWLS-SPAD方法重建的图像有效地去除了低剂量X射线导致的噪声，并且抑制了条纹伪影的产生和阶梯效应所引起的环状纹理。同时，在SSIM、FSIM、RMSE和PSNR等评估指标上，PWLS-SPAD方法重建的Shepp-Logan体模图像、XCAT体模图像和临床数据的图像都展示了更佳的数值表现。

阶梯效应是由多种原因导致的，迭代次数过多使得图像被过度平滑，或者梯度信息的不准确也会引起阶梯效应的产生^{［24， 30］}。PWLS-SPAD方法根据亚像素技术得到的新的梯度信息，可以更加准确地获取图像梯度信息，避免梯度效应损害低剂量CT图像质量。虽然提出的PWLS-SPAD方法在低剂量CT重建较对比方法有更好地表现，但是PWLS-SPAD方法依然有一定不足。PWLS-SPAD方法为了确保收敛，使用了固定的迭代次数，会对计算时间造成一定的消耗。同时对惩罚加权最小二乘模型中的惩罚参数的选取是根据手动试错及重建图像效果来决定的，目前针对惩罚参数的选取还是一个广泛研究的问题^{［31， 32］}。

综上所述，由于低剂量CT重建的图像存在明显的噪声和伪影，影响了低剂量CT在临床的应用。为了解决低剂量CT重建图像的质量问题，提出的PWLS-SPAD方法将基于亚像素技术计算得到的梯度信息引入到各项异性扩散过程中，并结合PWLS模型对低剂量CT图像进行重建。实验结果表明，PWLS-SPAD方法可以有效地提高低剂量CT图像质量。

参考文献

[1] Vliegenthart R, Fouras A, Jacobs C, et al. Innovations in thoracic imaging: CT, radiomics, AI and X-ray velocimetry[J]. Respirology, 2022, 27(10): 818-33.

[2] Bonney A, Malouf R, Marchal C, et al. Impact of low-dose computed tomography (LDCT) screening on lung cancer-related mortality[J]. Cochrane Database Syst Rev, 2022, 8(8): CD013829.

[3] Gillies RJ, Schabath MB. Radiomics improves cancer screening and early detection[J]. Cancer Epidemiol Biomark Prev, 2020, 29(12): 2556-67.

[4] 牛善洲, 唐诗洲, 黄舒彦, et al. 基于高维PDE投影恢复的低剂量CT重建方法 [J]. 南方医科大学学报, 2024, 44(4): 682-8.

[5] Zhang H, Wang J, Zeng D, et al. Regularization strategies in statistical image reconstruction of low-dose X-ray CT: a review[J]. Med Phys, 2018, 45(10): e886-e907.

[6] 邸江磊, et al. 基于深度学习的稀疏或有限角度CT重建方法研究综述 [J]. 激光与光电子学进展, 2023, 60(8): 42-79.

[7] 牛善洲, 梁礼境, 李 硕, et al. 基于低维流形先验的低剂量CT重建方法[J]. 计算机工程与应用, 2023, 59(18): 242-8.

[8] 李进丹, 陈 龙, 杨聪慧, et al. 低剂量PET/CT及PET/MRI临床研究进展[J]. 中国医学影像技术, 2021, 37(11): 1740-3.

[9] 陈世宣, 曾栋, 边兆英, et al. 基于联邦特征学习的多机型低剂量CT重建算法 [J]. 南方医科大学学报, 2024, 44(2): 333-43.

[10] Yang SH, Pu Q, Lei CT, et al. Low-dose CT denoising with a high-level feature refinement and dynamic convolution network[J]. Med Phys, 2023, 50(6): 3597-611.

[11] Wu DF, Kim K, El Fakhri G, et al. Iterative low-dose CT reconstruction with priors trained by artificial neural network[J]. IEEE Trans Med Imaging, 2017, 36(12): 2479-86.

[12] Zhang H, Ma JH, Wang J, et al. Statistical image reconstruction for low-dose CT using nonlocal means-based regularization. Part II: an adaptive approach[J]. Comput Med Imaging Graph, 2015, 43: 26-35.

[13] Zhang H, Ma JH, Wang J, et al. Statistical image reconstruction for low-dose CT using nonlocal means-based regularization. Part II: an adaptive approach[J]. Comput Med Imaging Graph, 2015, 43: 26-35.

[14] Niu SZ, Liu H, Zhang MZ, et al. Iterative reconstruction for low-dose cerebral perfusion computed tomography using prior image induced diffusion tensor[J]. Phys Med Biol, 2021, 66(11): 1150-24.

[15] Wang J, Li TF, Lu HB, et al. Penalized weighted least-squares approach to sinogram noise reduction and image reconstruction for low-dose X-ray computed tomography[J]. IEEE Trans Med Imaging, 2006, 25(10): 1272-83.

[16] Wang J, Lu HB, Wen JH, et al. Multiscale penalized weighted least-squares sinogram restoration for low-dose X-ray computed tomography[J]. IEEE Trans Biomed Eng, 2008, 55(3): 1022-31.

[17] 张心如, 周先春, 汪志飞, et al. 基于改进各向异性扩散的图像去噪算法 [J]. 电子测量技术, 2022, 45(17): 113-9.

[18] Guo YH, Cheng HD. Image noise removal approach based on subpixel anisotropic diffusion[J]. J Electron Imaging, 2012, 21(3): 033026-1.

[19] Zhang HJ, Zeng D, Lin JH, et al. Iterative reconstruction for dual energy CT with an average image-induced nonlocal means regularization[J]. Phys Med Biol, 2017, 62(13): 5556-74.

[20] Huang J, Zhang YW, Ma JH, et al. Iterative image reconstruction for sparse-view CT using normal-dose image induced total variation prior[J]. PLoS One, 2013, 8(11): e79709.

[21] Zhang HW, Zhang PC, Cheng WT, et al. Learnable PM diffusion coefficients and reformative coordinate attention network for low dose CT denoising[J]. Phys Med Biol, 2023, 68(24): Physicsinmedi-cineandbiologyvol.68, 2410.1088/1361-6560/aced33.11Dec.2023, .

[22] Wang J, Guan HQ, Solberg T. Inverse determination of the penalty parameter in penalized weighted least-squares algorithm for noise reduction of low-dose CBCT[J]. Med Phys, 2011, 38(7): 4066-72.

[23] Liu Y, Ma JH, Fan Y, et al. Adaptive-weighted total variation minimization for sparse data toward low-dose X-ray computed tomography image reconstruction[J]. Phys Med Biol, 2012, 57(23): 7923-56.

[24] 牛善洲, 刘宏, 朱赟, et al. 基于广义惩罚加权最小二乘的低剂量CT重建方法 [J]. 数值计算与计算机应用, 2021, 42(3): 289-302. DOI: 10.12288/szjs.s2020-0671

[25] 牛善洲, 张梦真, 邱 洋, et al. 基于全广义变分约束加权最小二乘的低剂量计算机断层重建方法[J]. 激光与光电子学进展, 2023, 60(4): 3788/LOP212853. DOI: 10.3788/LOP212853

[26] Zhang H, Capaldi D, Zeng D, et al. Prior-image-based CT reconstruction using attenuation-mismatched priors[J]. Phys Med Biol, 2021, 66(6): 064007.

[27] Zhang L, Zhang L, Mou XQ, et al. FSIM: a feature similarity index for image quality assessment[J]. IEEE Trans Image Process, 2011, 20(8): 2378-86.

[28] Tsai MY, Liang HL, Chuo CC, et al. A novel protocol for abdominal low-dose CT scans adapted with a model-based iterative reconstruction method[J]. J Xray Sci Technol, 2023, 31(3): 453-61.

[29] Wilson DO. Is there more to low dose CT scans[J]? Chest, 2022, 161(4): 880-1.

[30] Sun T, Sun NB, Wang J, et al. Iterative CBCT reconstruction using Hessian penalty[J]. Phys Med Biol, 2015, 60(5): 1965-87.

[31] Niu S, Bian Z, Zeng D, et al. Total image constrained diffusion tensor for spectral computed tomography reconstruction[J]. Applied Mathematical Modelling, 2019, 68: 487-508.

[32] Niu SZ, Zhang MZ, Qiu Y, et al. Evaluation of low-dose computed tomography reconstruction using spatial-radon domain total generalized variation regularization[J]. Phys Med Biol, 2024, 69(10): 1050-05.

基金资助

RIGHTS & PERMISSIONS

版权所有©《南方医科大学学报》编辑部2021