VVBP-Tensor in the FBP Algorithm: Its Properties and Application in Low-Dose CT Reconstruction

FBP算法中的VVBP-Tensor：性质及其在低剂量CT重建中的应用

论文链接：https://ieeexplore.ieee.org/document/8795540

项目链接：https://github.com/xxxxtao/VVBPTensor-net

Abstract

几十年来，商用X射线计算机断层扫描(CT)扫描仪一直使用滤波反投影(FBP)算法进行图像重建。然而，对低辐射剂量的渴望已经将FBP算法推向了极限。以往的研究通过对正弦图进行预处理、修改斜坡滤波器或对重建图像进行后处理来改善FBP的结果。本文主要对FBP算法中逐视图反向投影(VVBP-Tensor)进行分析和处理。我们分析的一个关键挑战在于每个反向投影切片的径向结构。为了克服这个困难，在VVBP-Tensor的z方向(投影视图的方向)引入了一个排序操作。结果表明，经过排序后的张量包含了与物体相似的结构，并且张量的不同切片中的结构具有相关性。然后，我们分析了VVBP-Tensor的性质，包括结构自相似性、张量稀疏性和噪声统计。考虑到这些特性，我们开发了一种使用张量奇异值分解(称为VVBP-tSVD)的算法来对Low-mAs CT成像的VVBP-Tensor进行降噪。实验采用物理phantom和不同mAs水平的临床患者数据进行。结果表明，在正弦图预处理、图像后处理和迭代重建等不同重建方案下，VVBP-tSVD都优于所有竞争方法。我们得出结论，VVBP-Tensor是改善FBP重建质量的适当处理目标，并且提出的VVBP-tSVD是Low-mAs CT成像中有效的噪声降低算法。这一初步工作可能为回顾和反思FBP算法提供一个启发式的视角。

I. INTRODUCTION

X射线计算机断层扫描(CT)已广泛应用于各种临床应用，如筛查、诊断和图像引导干预[1]。在临床环境中，与其他方式相比，CT最重要的特点之一是其快速准确的解剖成像能力。为了保证这一性能，提出了许多先进的技术。一个典型的例子是图像重建，通常使用滤波后的反向投影(FBP)算法，因为它的重建速度快，适用性好[2]。然而，在低剂量CT (LDCT)成像中，如果不进行适当的处理，FBP重建图像会出现严重的噪声伪影[3]。其主要原因是FBP是Radon逆变换的解析近似，当原始数据受到噪声、数据缺失或高衰减材料的影响时，FBP会失败[4]。

一般来说，FBP的实现可分为滤波、逐视图反投影(VVBP)和求和三个步骤，如图1(a)-©所示。在第一步中，通常采用斜坡滤波器来去除纯反向投影中的模糊[5]。在第二步中，通过插值将过滤后的正弦图的每个视图反投影到图像域中。第三步，对所有视图沿z方向的反向投影求和，得到图像。在FBP框架下，许多研究报道了基于FBP的先进算法，包括正弦图预处理[6]-[9]、滤波器修改[10]-[13]、图像后处理[14]-[18]。在这些研究中，直接对每个视图的反向投影求和，得到重建图像。然而，逐视图的反向投影数据很少被讨论。在本研究中，我们重点研究了FBP算法中VVBP数据的分析和处理。

从图1©中可以看出，VVBP数据可以堆叠成一个三阶张量，为简单起见，将其命名为“VVBP-Tensor”。另一方面，在图1(b)中显示的VVBP-Tensor切片包含类似于径向结构，使得随后的分析变得具有挑战性。为了解决这个问题，我们引入了“排序”作为VVBP-Tensor在z方向上的辅助操作。注意，VVBP-Tensor的“z方向”也是投影视图的方向。如图1(d)所示，VVBP-Tensor的每个mode-3向量根据其像素值独立排序。我们发现，经过排序后，VVBP-Tensor中包含与目标相似的结构，并且不同切片中的结构具有相关性，如图1(e)所示。据我们所知，这项工作是第一个在FBP算法中分析和利用VVBP-Tensor的研究。具体而言，本工作的贡献总结如下：

1. 我们分析了VVBP-Tensor的性质，包括结构自相似性、张量稀疏性和噪声统计。这些特性促使我们开发新的算法来处理VVBP-Tensor，以产生高质量的FBP重建。
2. 为了挖掘其在Low-mAs CT成像中的潜力，我们提出了一种称为VVBP-tSVD的VVBP-Tensor去噪算法，该算法充分考虑了VVBP-Tensor的内在特性，结合了三种策略：张量奇异值分解(tSVD)[19]、基于全带patch的块匹配(FPBM)[20]和惩罚加权最小二乘(PWLS)模型[7]。
3. 在不同mAs水平下进行了大量的物理phantom实验和临床患者数据。采用正弦图预处理、图像后处理和迭代重建等策略进行比较。结果表明了VVBP-tSVD算法相对于竞争方法的优越性。

本文组织如下。第二节回顾了一些与FBP相关的CT重建方法。第三节介绍了张量分析的记号和基础。第四节详细介绍了VVBP-Tensor分析。第五节给出了提出的VVBP-tSVD算法。第六节给出了实验设置和结果。最后，第七节给出讨论和结论。

II. FBP-RELATED WORKS

在本节中，我们回顾了一些与FBP相关的CT重建工作，特别是LDCT重建，大致分为三大类：正弦图预处理、滤波器修改和图像后处理。

A. 正弦图预处理

正弦图预处理的目的是使用特定设计的算法对正弦图进行去噪。例如，La Riviere等[6]提出了一种移位泊松模型来描述投影数据中的噪声，并在预对数域中提出了一种高斯模型来对噪声进行惩罚，形成了所谓的惩罚似然(penizize -likelihood, PL)正弦图恢复方法。王等[7]描述了用于恢复正弦图数据的惩罚加权最小二乘(PWLS)准则。Gao等[8]采用自适应正弦图恢复方法和POCS-TV算法，结合Low-mAs和稀疏视图协议重建CT图像。最近，Xie等[9]提出了一种利用正弦图的噪声产生机制，利用二阶总变分(TV)作为正则化项对其进行处理的方法，在已知的正弦图预处理方法中，该方法的处理结果是最先进的。正弦图预处理技术通常利用了正弦图数据的统计特性，而忽略了目标的结构信息，导致重建图像中出现径向模糊伪影。

B. 滤波器修改

除了正弦图预处理之外，FBP还可以通过改进的滤波器进行。一种直接的方法是用窗函数修改斜坡滤波器，如Hann和Hamming窗[4]，以抑制正弦图数据中的高频成分。除了加窗滤波器，Hsieh[10]使用自适应滤波器来减少由光子不足伪影(photon starvation)引起的条纹伪影。kachhelrieß等[11]将自适应滤波器推广到多个维度，用于轴向和螺旋CT重建。Zeng[12]引入了一种逐视图噪声加权(vFBP)的FBP算法来模拟迭代Landweber重建。后来，Zeng等人[13]将view-by-view噪声加权改为ray-by-ray noise weighting (rFBP)，这是自适应滤波器的改进版本。这些方法可以很好地减少条纹伪影；然而，需要额外的后滤波来完全去除图像噪声。

C. 图像后处理

降低CT图像中的噪声和伪影最直观的方法是对重建后的图像进行处理。如Borsdorf等[14]采用基于小波的方法，通过相关分析降低CT图像中的噪声。Supanich等[15]和Leng等[16]采用了HYPR方法对时间分辨CT血管造影和光谱CT图像进行去噪。Ma等[17]提出了一种非局部均值(Nonlocal Means, NLM)框架下正常剂量诱导的低剂量CT恢复方法。关于NLM算法在LDCT成像中的应用的详细综述可以在[21]中找到。Zhang等[18]提出了一种利用从全剂量训练数据库中学习到的特征恢复肺部低剂量CT图像的算法。近年来，基于深度学习的方法在LDCT重建中得到了应用，大多数方法都被用于图像的后处理，并取得了最先进的效果[22]-[27]。图像后处理策略通常只考虑图像的属性，而忽略了原始数据本身的统计属性。

III. NOTATIONS AND PRELIMINARIES

张量是维数大于3的多维数据。一个$N$阶张量可以表示为$\mathcal{X}\in {\mathbb{R}}^{O_1\times O_2\times \cdots \times O_N}$。$\mathbf{X}$的$n$阶向量是在保持其他维度不变的情况下沿第$n^{th}$维$O_n$提取元素得到的向量。例如，图1所示的VVBP-Tensor的模态向量$t_i^3$表示为：${\mathbf{t}}_i^3=\mathcal{T}\left(i_1,i_2,:\right)$，其中$i_1$和$i_2$是$\mathcal{T}\in {\mathbb{R}}^{O_1\times O_2\times O_3}$下$i^{th}$元素相对于矩阵维数$\left(O_1,O_2\right)$的子指标。张量$\mathcal{X}$的$n^{th}$模展开矩阵$X^n$是将张量$\mathcal{X}$沿$n^{th}$维重塑成矩阵得到的，可以表示为：
$$X^n={\mathrm{unfold}}^n(\mathbf{X})\in {\mathbb{R}}^{O_n\times \left(O_1\cdots O_{n-1}{O}_{n+1}\cdots O_N\right)}.\text{\hspace{1em}(1)}$$
低秩和稀疏性是张量分析中考虑的两个常见性质[19]，[20]，[28]-[33]。由于秩优化是一个非凸问题，通常使用核范数来逼近对象的秩[33]。张量的核范数(TNN)可以通过沿张量模的每个方向对展开的矩阵应用奇异值分解(SVD)[31]、基于张量的SVD (tSVD)[19]、Tucker分解[34]等方法得到。TNN最小化问题可以通过阈值化代表系数，如张量奇异值来解决。

在本文中，我们采用tSVD来解决通常采用软阈值法解决的VVBP-Tensor去噪的TNN最小化问题：
KaTeX parse error: Expected 'EOF', got '&' at position 2: &̲D_\tau(\mathcal…
式中$\mathcal{M}\in {\mathbb{R}}^{O_1\times O_2\times O_3}$为待处理的张量；$\mathcal{U}$和$\mathcal{V}$分别是尺寸为$O_1\times O_1\times O_3$ 和 $O_2\times$ ${\mathrm{O}}_2\times {\mathrm{O}}_3$的正交张量；$\mathcal{S}$是一个大小为$O_1 \times O_2 \times O_3 $的矩形f对角张量；$D_t$为软阈值算子；$\sigma$表示$S$中的单个值；“f-diag $(\cdot )$”表示$f$对角张量的变换算子；“$\ast$”表示t积；$\tau$为正阈值；$t_{+}$表示$t$的正部分，即$t_{+}=\max (0,t)$。当$t\ge 0$时，公式(2)符合以下条件[19]，[35]：
$$D_{\tau }(\mathcal{M})=\underset{\mathcal{X}}{\mathrm{argmin}}\left\{\frac{1}{2}\Vert \mathcal{X}-\mathcal{M}{\Vert }_2^2+\tau \Vert \mathcal{X}{\Vert }_{\mathrm{TNN}}\right\}\text{\hspace{1em}(3)}$$
其中“TNN”为张量核范数[31]。

IV. VVBP-TENSOR ANALYSIS

在本节中，我们首先给出了VVBP-Tensor的生成过程，然后给出了一个VVBP-Tensor分析的流，最后分析了它的性质。

A. VVBP-Tensor生成

在不失一般性的前提下，我们以平行光束几何为例。扇形束和锥形束的几何形状可以很容易地推导出来。平行光束几何的FBP可表示为[5]：
$$\mu (x,y)={\int }_0^{\pi }{\int }_{-\infty }^{\infty }p(s,\theta )h(s^{\prime }-s){|}_{s=xcos\theta +ysin\theta }dsd\theta \text{\hspace{1em}(4)}$$
式中，$\mu (x,y)$表示图像域中二维连续坐标$(x,y)$的FBP重构图像；$p(s,\theta )$为detector-bin方向和曝光角$\theta$的一维连续坐标$s$的正弦图；$h(\cdot )$表示滤波器核。将公式(4)离散化表示为：
$$\mu (i,j)=\frac{\pi }{K}\sum _{k=1}^K\left\{{p(\hat{s},k)\ast \hat{h}|}_{\hat{s}=\mathrm{INT}[icos\theta +jsin\theta ]}\right\}\text{\hspace{1em}(5)}$$
式中$(i,j)$为$(x,y)$对应的二维离散坐标；$\hat{s}$为$s$对应的一维离散坐标；$k$为$\theta$对应的投影视图索引；$K$为view总次数；$h$为离散化滤波核；“$\ast$”表示一维卷积；“INT”表示用于反向投影的插值。根据公式(5)，每个视图的反向投影可计算为：
$$\hat{\mathcal{T}}(i,j,k)=p(\hat{s},k)\ast \hat{h}{|}_{\hat{s}=\mathrm{INT}[icos\theta +jsin\theta ]}\text{\hspace{1em}(6)}$$
式中，$\mathcal{T}\in {\mathbb{R}}^{I\times J\times K}$表示重构图像大小为($I\times J$)的VVBP-Tensor(通常为512 × 512)。

B. VVBP-Tensor分析流

图1为本研究的VVBP-Tensor分析流程。随着逐步实现，FBP算法大致可分为滤波、逐视图反投影和沿z方向求和三个步骤(图1中黄色箭头)。从视觉上看，VVBP之后的数据包含径向结构，如图1(b)所示，其mode-3向量为“无序”，如图1©所示。这些症状给VVBP-Tensor的分析带来了困难。

为了便于分析，我们在VVBP-Tensor的z方向上引入了一个排序操作。图3为虚体VVBP-Tensor排序示意图，其中$\mathcal{T}\in {\mathbb{R}}^{I\times J\times K}$为排序后的VVBP-Tensor。排序时，通过对VVBP-Tensor的每个mode-3向量按其值重新排序，独立排序。请注意，排序实际上是像素在z方向上的重新排列，这是线性的和可逆的。排序本身不会改变要重建的图像，可以通过在z方向上对T求和得到：$\mu (i,j)=\frac{\pi }{K}{\sum }_{k=1}^K\mathcal{T}(i,j,k)$。

如图1(e)所示，在排序之后，VVBP-Tensor包含与对象相似的结构，并且不同切片中的结构是相关的。为了进一步验证这一效果，我们收集了几例临床病例，并计算了其对应的分类VVBP-Tensor，如图2所示。这些数据涵盖了从肩膀到臀部的不同身体部位。很明显，所有分类的VVBP-Tensor切片都包含与患者相似的丰富结构。

C. VVBP-Tensor性质

1) 结构自相似性：为了分析VVBP-Tensor的结构相似性特性，我们从phantom VVBP-Tensor的5个相邻切片中裁剪出感兴趣区域(ROI)，并计算它们的相似性图。采用负欧几里得距离指数来描述ROI之间的相似度[17]，表示为：
$$C_{ij,k}=e^{-\Vert T(P_t)-T(P_{ij,k}){\Vert }^2}\text{\hspace{1em}(7)}$$
式中$C_{ij,k}$为第k个切片(i, j)处的相似度指数；$P_t$和$P_{ij,k}$分别表示ROI的目标patch和提取patch；项 T ( P • ) : = { t ( x l ) , x l ∈ P • } T(P_•):= \{t(x_l),x_l \in P_•\} T(P•​):={t(xl​),xl​∈P•​}表示限制在patch $P_{•}$中的邻域像素向量。图4显示了ROI(上一行)及其相似度图(下一行)。通过直观地检查第一行，我们可以在不同的切片中看到相同的结构。此外，第二行所示的相似性图定量地证明了VVBP-Tensor内部切片的结构间和结构内相似性。

2) 张量稀疏性(Tensor Sparsity)：图4所示的切片间隔为5，我们可以看到远离目标patch的切片仍然表现出很高的相似性。这种趋势也可以从图2中观察到，其中不同切片的结构非常相似。这一观察结果表明，VVBP-Tensor中的信息在三维空间上是高度冗余的，并启发我们分析VVBP-Tensor的稀疏性。为了进行稀疏性分析，我们将虚模VVBP-Tensor沿其三个模态展开，并计算每个展开矩阵的奇异值，其表示为：
{ U n , S n , V n } = S V D ( u n f o l d n ( T ) ) , n = 1 , 2 , 3 w h e r e S n = d i a g ( σ i n ) (8) \begin{aligned} \{U^n,S^n,V^n\} &=\mathrm{SVD}(\mathrm{unfold}^n(\mathcal{T})),&n=1,2,3\\ where \quad S^n&=\mathrm{diag}(\sigma_i^n)\end{aligned} \tag{8} {Un,Sn,Vn}whereSn​=SVD(unfoldn(T)),=diag(σin​)​n=1,2,3(8)
式中，$U_n$和$V_n$为两个正交基矩阵；$S_n$是一个对角矩阵，${\sigma }_i^n$表示$\mathcal{T}$的$n$型展开矩阵的第$i$个奇异值；“SVD(·)”表示奇异值分解算子；“diag(·)”表示对角元素。图5(b)绘制了每个矩阵沿其三个模态展开的奇异值。观察到，每条曲线在几个指标后都急剧下降，这表明VVBP-Tensor沿每个模态都位于低秩子空间，因此具有良好的张量稀疏性。

3) 噪声统计：在数据采集过程中，可以假设每个detector bin的噪声是独立的。根据我们之前的工作[36]，正弦图的噪声方差可以计算为：${\sigma }_{p_n}^2=\frac{1}{{\hat{p}}_n}(1+\frac{({\sigma }_{e,n}^2-1.25)}{{\hat{p}}_n})$，其中${\hat{p}}_n=I_0^n{e}^{-{\bar{p}}_n}$，其中$I_0^n$表示入射光子数，${\bar{p}}_n$表示第n个检测器的期望值；${\sigma }_{e,n}^2$为背景电子噪声方差。根据公式(6)，VVBP-Tensor在(i, j, k)处的噪声方差可估计为：
$${\sigma }_{\hat{\mathcal{T}}}^2(i,j,k)={\sigma }_p^2(\hat{s},k)\ast {\hat{h}}^2{|}_{\hat{s}={\mathrm{INT}}^2[icos\theta +jsin\theta ]}.\text{\hspace{1em}(9)}$$
在反向投影的插值过程中，像素$(i,j,k)$的值是其相邻像素的加权平均值。因此，像素$(i,j,k)$处的噪声方差可以估计为相邻像素噪声方差的平方加权平均，表示为式(9)中的“INT2”。由于排序操作是在z方向上对像素进行重排，因此排序后的VVBP-Tensor的噪声方差可以通过对${\sigma }_{\mathcal{T}}^2$按相同的顺序排序来计算。FBP中噪声传播的详细参考文献[37]。

为了验证VVBP-Tensor的噪声估计，我们模拟了100个相同剂量水平(20 mAs)的phantom投影，并计算了它们的VVBP-Tensor和相应的噪声方差图。如图6所示的结果表明，通过公式(9)可以有效地估计VVBP-Tensor的噪声方差。此外，每个像素点的噪声直观地表示为高斯分布。

V. VVBP-TENSOR DENOISING FOR LOW-MAS CT

在第四节分析的基础上，设计了一种VVBP-tSVD算法对Low-mAs CT成像的VVBP-Tensor进行降噪，该算法利用tSVD[19]利用张量稀疏性，利用FPBM[20]利用结构自相似性，利用PWLS[7]纳入噪声统计。

A. VVBP-tSVD去噪模型

由于VVBP-Tensor的噪声具有类高斯分布，因此使用PWLS模型[7]对张量进行去噪是合理的。本文提出的VVBP-tSVD算法的目标函数为：
$$\mathcal{X}=\arg \underset{\mathcal{X}}{\min }\left\{\frac{\lambda }{2}\Vert \mathcal{X}-\mathcal{T}{\Vert }_W^2+\sum _{i=1}^N{\gamma }_i\Vert {\mathcal{X}}_i{\Vert }_{\mathrm{TNN}}\right\}\text{\hspace{1em}(10)}$$
式中$\mathcal{T}$和$\mathcal{X}$分别为低剂量和估计的VVBP-Tensor；$W$是一个加权矩阵，其对角元素为${\sigma }_{\mathcal{T}}^2$的倒数，其他元素为零；$\lambda$是一个正超参数；第二项中的“${\sum }_{i=1}^N$”表示全带patch群的集合；$N$是这些组的总数；${\gamma }_i$为$i^{\text{th }}$全带patch群的正则化参数，表示为${\gamma }_i=\beta {\sigma }_{{\mathbf{T}}_i}^2$，其中$\beta$ 为超参数，${\sigma }_{{\mathcal{T}}_i}^2$为${\mathcal{T}}_i$中各元素噪声方差的平均值，即$\mathcal{T}$中的${\mathcal{X}}_{\mathbf{i}}$。

B. 基于全频带块的块匹配

如IV-C节所分析，VVBP-Tensor包含与目标相似的结构，并且这些结构是相关的。为了充分利用这一特性，我们将FPBM[20]策略纳入了VVBP-Tensor去噪过程。在这里，全带patch被称为子张量，它是通过在每个切片中提取相同位置的2D patch并将其堆叠成一个3D体而获得的。如图7所示，在每次迭代中，首先从${\mathcal{X}}^{l-1}$通过分块匹配和展开得到$N$组${\left\{{\mathcal{X}}_i^{l-1}\right\}}_{i=1}^N$的全带patch。然后使用公式(12)估计过滤后的全带patch组${\left\{{\mathcal{Z}}_i^l\right\}}_{i=1}^N$。通过对所有滤波组进行折叠和聚合，得到去噪的VVBP-Tensor ${\mathcal{Z}}^l$。

C. 优化

由于正则化函数是非线性的，很难直接最小化公式(10)。在这项工作中，我们使用半二次分裂(HQS)算法[38]来放松公式(10)。松弛目标函数为：
$$\frac{\lambda }{2}\Vert \mathcal{X}-\mathcal{T}{\Vert }_W^2+\sum _{i=1}^N\left\{\frac{\rho }{2}{\Vert {\mathcal{Z}}_i-{\mathcal{X}}_i\Vert }_2^2+{\gamma }_i{\Vert {\mathcal{Z}}_i\Vert }_{\mathrm{TNN}}\right\}\text{\hspace{1em}(11)}$$
其中$\mathcal{Z}$ 是一个近似于$\mathcal{X}$的辅助变量，$\rho$是一个超参数。利用HQS算法，公式(11)可以分解为两个子问题：
$$\begin{array}{rl}& {\mathcal{Z}}_i^l=\underset{{\mathcal{Z}}_i}{\mathrm{argmin}}\left\{\frac{{\rho }^l}{2}{\Vert {\mathcal{Z}}_i-{\mathcal{X}}_i^{l-1}\Vert }_2^2+{\gamma }_i{\Vert {\mathcal{Z}}_i\Vert }_{\mathrm{TNN}}\right\}\\ & {\mathcal{X}}^l=\underset{\mathcal{X}}{\mathrm{argmin}}\left\{\lambda \Vert \mathcal{X}-\mathcal{T}{\Vert }_W^2+{\rho }^l\sum _{i=1}^N{\Vert {\mathcal{Z}}_i^l-{\mathcal{X}}_i\Vert }_2^2\right\}\end{array}\text{\hspace{1em}(12)(13)}$$
其中，${\rho }^l$随着迭代的继续而增大，${\rho }^{l+1}=\alpha {\rho }^l$，其中$\alpha$为大于1的用户定义标量；$l=1$, $2,\dots ,L$为迭代步长。更新 ${\rho }^l$ 使得$\mathcal{Z}$成为对$\mathcal{X}$的越来越好的近似，因此使得公式(11)对公式 (10)的逼近越来越好。在实现中，迭代最小化公式(12)和公式(13)求解原目标函数。

利用张量奇异值分解，可将公式(12)更新为：
$$\begin{array}{rl}{z}_i^l& =D_{{\tau }_i}({\mathcal{X}}_i^{l-1})={\mathcal{U}}_i\ast D_{{\tau }_i}({\mathcal{S}}_i)\ast {\mathcal{V}}_i^T,\\ where{D}_{{\tau }_i}({\mathcal{S}}_i)& =\mathrm{f}-\mathrm{diag}(({\sigma }_i-{\tau }_i{)}_{+})\end{array}\text{\hspace{1em}(14)}$$
其中${\tau }_i=\frac{{\gamma }_i}{{\rho }^l}$为正阈值。

公式(13)的闭式解为：
$${\mathcal{X}}^l=\frac{\lambda \mathcal{T}+{\rho }^l{\sigma }_{\mathcal{T}}^2\odot {\sum }_{i=1}^N{\mathcal{Z}}_i^l}{\lambda +{\rho }^l{\sigma }_{\mathcal{T}}^2}\text{\hspace{1em}(15)}$$
其中$\odot$为逐点积。

在实施这些策略时，应注意两点。首先，VVBP-Tensor中不同切片的强度范围可能相差很大。这种变化可以从图1(d)和图3中排序后绘制的向量中观察到。因此，在每个迭代步骤进行块匹配之前，对VVBP-Tensor进行逐片归一化处理。聚集后，将X反规格化回原灰空间。其次，本研究使用的phantom VVBP-Tensor和患者数据分别包含1160个和576个切片。考虑到计算负载，我们对VVBP-Tensor及其在z方向上的噪声方差进行了降采样：
$$\begin{array}{c}\tilde{\mathcal{T}}(:,:,q)=\frac{1}{d}\sum _{p=(q-1)d+1}^{qd}\mathcal{T}(:,:,p)\\ {\sigma }_{\tilde{\mathcal{T}}}^2(:,:,q)=\frac{1}{d^2}\sum _{p=(q-1)d+1}^{qd}{\sigma }_{\mathcal{T}}^2(:,:,p)\end{array}\text{\hspace{1em}(16)}$$
其中$\mathcal{T}\in {\mathbb{R}}^{I\times J\times (K/d)}$为下采样的VVBP-Tensor；${\sigma }_{\mathcal{T}}^2$为$\mathcal{T}$的噪声方差；$d$表示下采样因子；$q=1,2,\dots ,K/d$为降采样后的切片指数。利用仿真数据进行的实验表明，采用降采样策略可以大大缩短处理时间，但结果的质量略有下降。结果见第VI-B节。

算法1给出了所提出的VVBP-tSVD算法的伪代码。

VI. EXPERIMENTS AND RESULTS

A. 实验设置

1) 数据采集：实验使用的数据包括来自低剂量CT大挑战数据集的物理phantom和三名患者[39]。采用临床CT扫描仪(Siemens SOMATOM Sensation 16 CT扫描仪)轴向风扇模式扫描拟人化躯干phantom。设置1级管电压(120 kVp)和3级管电流(17、40和60 mAs)。此外，通过在100 mAs和120 kVp下平均150次重复测量获得的正弦图被用作ground truth。

使用Siemens SOMATOM Definition Flash CT系统在螺旋模式下获取患者数据。正常剂量和模拟四分之一剂量的原始投影由宿主提供。在我们的实验中，我们选择了三个腹部病例(如图2所示的病例3、4和5)。为了测试所提出的方法在较低剂量水平下的能力，使用[40]中提出的方法模拟了1/6剂量和1/8剂量水平下的数据。在本研究中，螺旋数据使用单层重采样方法[41]重新分bin为平行光束几何。

2)对比方法：为了验证算法的性能，选择传统的FBP作为基线对比方法。选择IMAP [9]和BM3D [42]来代表正弦图预处理方法和图像后处理方法，分别缩写为SP-IMAP和IP-BM3D。此外，使用TV[43]和KSVD[44]作为迭代重建的正则化项，分别缩写为IR-TV和IR-KSVD。此外，为了验证tSVD在利用张量稀疏性方面的有效性，我们采用NLM[45]作为正则化项对VVBP-Tensor(简称VVBP-NLM)进行降噪，其表达式为：
$$\arg \underset{\mathcal{X}}{\min }\left\{\frac{\lambda }{2}\Vert \mathcal{X}-\mathcal{T}{\Vert }_W^2+\sum _j\sum _{k\in {\Omega }_j}{w}_{jk}({\mathcal{X}}_j-{\mathcal{X}}_k{)}^2\right\}\text{\hspace{1em}(17)}$$
其中j, k是体素指数；根据非局部斑块相似度计算加权因子$w_{jk}$。注意，VVBP-NLM考虑了VVBP-Tensor的结构自相似性，但忽略了它的稀疏性。另外，利用HQS算法对公式(17)进行优化。为了确保公平的比较，所有竞争方法的参数都是手动调整的，以达到幻象数据的最大峰值信噪比(PSNR)值，并实现患者数据的最佳可视化。

B. 降采样因子的影响

为了充分利用降采样对VVBP-tSVD性能的影响，使用phantom数据进行实验，其中mAs为17，降采样因子变化从1到64。图8显示了一些恢复后的ROI (d = 1,8,16,32)，通过视觉对比，不同降采样因子的结果略有不同。

此外，计算恢复的ROI相对于地真值的PSNR值，并随计算时间绘制在图9中。给出了FBP结果的PSNR。从图9(a)可以看出，VVBP-tSVD得到的ROI恢复后的PSNR值比FBP的要好得多。此外，随着d的增加，由于下采样的信息丢失，PSNR值趋于降低，变得不稳定。从图9(b)可以看出，随着d的减小，计算时间急剧减少，但在降采样因子为20左右后，计算时间收敛。考虑到定量性能和计算时间，本研究的实验将降采样因子设置为8。

C. Phantom研究

图10显示了不同方法在17mAs时的phantom重建结果。可以观察到，用FBP算法重建的图像受到严重的噪声和伪影的污染。另一方面，VVBP-tSVD方法在噪声和伪影抑制以及分辨率保持方面取得了最好的效果。

图11比较了ground truth图像和六种不同方法图像的选定轮廓。轮廓线位于240 ~ 270的像素位置i = 100和j处，如图10(a)中的红线所示。由图可知，在所有方法中，VVBP-tSVD算法的相似度和空间分辨率最高。

所有重建的相对于ground truth的绝对残差图像如图12所示。结果进一步说明，VVBP-tSVD在分辨率上保持比所有其他竞争方法更好的性能。

图13对比了SP-IMAP、IR-TV、IR-KSVD、IP-BM3D、VVBP-NLM、VVBP-tSVD重建的ground truth图像与图像的ROI均值和标准差。选择两个ROI进行比较，如图10(a)中红框所示。与其他方法相比，VVBP-tSVD方法的标准差最小。

表1列出了17、40和60 mAs时phantom结果的PSNR、归一化均方误差(NMSE)和特征相似度(FSIM)[46]。可以看出，在所有方法中，VVBP-tSVD算法产生的数值是最好的。另一方面，VVBP-tSVD的计算时间比其他方法要长得多，这是该方法的缺点。

D. 患者研究

图14为四分之一剂量病例3的结果。很明显，所有的方法都能得到噪声较低的图像。然而，SP-IMAP在图像中产生的颗粒噪声与正常剂量图像不同。IR-TV产生的图像受到分段常数效应的影响。IR-KSVD图像中的结构是扭曲的。IP-BM3D图像中的噪声是不均匀的，并且图像包含额外的伪影。VVBP-NLM图像具有与SP-IMAP图像相似的问题，但在结构保存方面略有增益。此外，与VVBP-tSVD方法相比，所有的竞争方法都倾向于过度平滑图像边缘。视觉上，来自VVBP-tSVD的图像在分辨率保持和噪声特性方面是最接近正常剂量图像的。

图15为正常剂量图像和不同重建方法图像的ROI直方图。ROI取自病例3的左肾部分，如图14(a)中红框所示。可以看出，VVBP-tSVD图像的直方图与正常剂量图像的直方图最接近，这进一步证明了所提出的VVBP-tSVD算法的优越性。

为了测试所提出的算法在较低剂量水平下的能力，我们在1/4、1/6和1/8剂量下重建了所有患者数据。图16为病例3在不同剂量水平下的ROI。可以看出，VVBP-tSVD在降噪和分辨率保持方面产生了很好的增益；然而，在低剂量的图像中，尤其是1/8剂量的图像中，微小的结构可能已经丢失。

类似的结果可以在图17中看到，与其他方法相比，VVBP-tSVD仍然产生最好的结果。

为了评估所有方法的临床表现，邀请了五位放射科医生(来自中国广州南方医院，至少有三年的CT图像阅读经验)，根据噪声和伪影减少，分辨率和结构保存，以及病变的可检测性，对所有患者图像进行评分，从0(最差)到5(最好)。为了确保公平的评估，未压缩的患者图像被随机发送给放射科医生，而没有向他们展示任何有关重建方法的信息。每位放射科医生的图像读取程序保持独立。不同患者评分的平均值和标准差列于表2。可以观察到，VVBP-tSVD方法在所有竞争方法中获得了最高分。患者研究结果表明，该方法在视觉检查和主观评价方面均优于其他方法。

VII. DISCUSSION AND CONCLUSION

本文没有对正弦图进行预处理[6]-[9]，修改斜坡滤波器[10]-[13]，也没有对图像进行后处理[14]-[18]，而是重点对FBP重构过程中的VVBP-Tensor进行分析和处理。我们的关键发现是，VVBP-Tensor包含与物体相似的结构，如图1(e)和2所示。另一方面，VVBP-Tensor的每个像素的噪声方差是可预测的(公式(9)和图6)。因此，处理VVBP-Tensor既可以考虑图像的结构特征，也可以考虑原始数据的噪声统计，而图像后处理或正弦图预处理只能利用其中一个方面。此外，VVBP-Tensor的稀疏特性(图5)可以增强其在Low-mAs CT成像中的降噪适用性。

在分析了VVBP-Tensor特性的基础上，提出了一种Low-mAs CT张量去噪算法(VVBP-tSVD)。VVBP-tSVD算法通过结合三种策略(FBPM[20]、PWLS[7]和tSVD[19])来考虑VVBP-Tensor的特性。不同剂量水平的物理Phantom和临床患者研究结果表明，在不同的重建方案下，包括正弦图预处理、图像后处理和迭代重建，VVBP-tSVD都优于所有对比方法。

注意，排序是VVBP-Tensor分析的重要步骤。在原始的VVBP-Tensor -T中，模式3向量在视觉上是“无序的”(图1©和3)，这使得我们的分析变得困难。为了解决这个问题，我们在z方向上对每个mode-3向量进行排序，因为排序是数据分析的基本方法之一。排序后，VVBP-TensorT中的切片包含与对象相似的结构。本文介绍了分类的基本思想、分类后的VVBP-Tensor的性质、VVBP-tSVD算法及其在低mas CT成像中的应用。建立排序操作的数学模型有助于更深入地理解FBP算法，这是未来的研究课题。

考虑到本文分析的结构自相似性、张量稀疏性和噪声统计，可以利用VVBP-Tensor的更多特性，如mode-3向量之间的相关性，或多维层析成像(如4DCT[47]、频谱CT[48]、灌注CT[32])的高维信息。此外，可以采用其他处理方法来代替所提出的VVBP-tSVD，以利用VVBP-Tensor的特性。例如，可以使用更强大的基于张量的算法(例如KBR[28])来考虑VVBPTensor的稀疏性。我们还注意到，其他层析成像问题(如稀疏视图CT[49]、有限角度CT[50]、PET和SPECT图像重建[51]、[52])也可以从VVBP-Tensor处理中受益。这里讨论的所有要点都是未来工作的潜在主题。

另一点需要注意的是，来自同一物体的High-mAs和Low-mAs CT的VVBP-Tensor是不同的。这种差异的产生是因为FBP中使用的斜坡滤波器在某种程度上是高通滤波器[5]，因此，噪声Low-mAs的正弦图往往比High-mAs的正弦图具有更宽的灰度范围。因此，Low-mAs的VVBP-Tensor的灰度范围比High-mAs的更宽。因此，VVBP-Tensor去噪的目标应该是降低Low-mAs VVBPTensor中的噪声，而不是将其恢复为High-mAs。图18显示了一些典型的ground truth、低剂量和去噪的VVBP-Tensor切片，其中地面真值切片与17 mAs和去噪切片相比处于不同的灰度级。

该方法的一个主要限制是随着数据量的增加而增加的计算负荷，而广泛使用的锥型束CT进一步加剧了这一问题。考虑到VVBP-Tensor中的信息是高度冗余的，本文提出了一种降采样策略来减轻计算负荷。此外，该算法采用基于patch的处理方案，易于并行化。因此，采用多GPU计算等高性能计算技术可以进一步减少计算时间[53]。

目前的工作还有其他局限性。首先，采用试错法对VVBP-tSVD算法的参数进行了调优，耗时长，且这些参数可能因不同类型的数据而有所不同。如何选择最优的参数组合仍然是一个悬而未决的问题。在未来的研究中，将在所提出的算法上测试一些自动参数调整的方法，比如深度学习。其次，本文提出的方法在真实的Low-mAs Phantom数据上进行了测试；然而，其处理实际Low-mAs患者数据的能力尚未得到验证。第三，我们将我们的方法与最先进的正弦图处理方法(IMAP[9])进行了比较，但尚未与最先进的图像域处理/正则化方法(例如，基于深度学习的方法[24]-[27])进行比较。这些比较需要在未来的研究中进行。

综上所述，本文的分析表明，VVBP-Tensor是提高FBP重建质量的合适处理目标；Phantom和患者实验结果表明，所提出的VVBP-tSVD算法在Low-mAs CT成像降噪方面是有效的。本研究为回顾和反思FBP算法提供了一个启发式的视角。