针对不完全投影数据的CT图像重建算法可以分为2类：第1类只适用于特定的扫描结构，首先将不完整的投影数据进行插值补充完整，然后再采用滤波反投影(Filtered Back-Projection，FBP)算法进行图像重建；第2类是迭代类的CT图像重建算法，如代数重建技术(Algebraic Reconstruction Technique，ART)，但传统迭代类算法的计算时间要求高，在早期CT中没有得到应用。近年来，随着大规模并行计算技术的发展以及计算机硬件成本的降低，迭代类算法已经成为相关领域研究人员和CT机生产厂商高度关注的热点。这其中较典型的是Sidky等在2006年提出的TV(Total Variation)算法[7-8]，该算法利用梯度图像稀疏性先验知识提高了稀疏投影条件下的重建图像质量。TV算法具有很高的实用性，但是TV算法的图像平滑效应会降低图像的空间分辨率和对比度。TV算法是一种启发式的最优化计算方法，压缩感知(Compressed Sensing，CS)理论[9]在一定程度上为TV算法良好的应用性能提供理论支撑。随后又出现诸多基于CS理论的CT迭代图像重建算法，较知名的是美国威斯康辛大学Chen博士及其研究组利用动态CT图像序列相关性先验知识提出的先验图像约束压缩感知(Prior Image Constrained Compressed Sensing，PICCS)算法[10]。PICCS算法要求能够从冗余或完整的投影数据集合中选取均匀采样的投影数据并利用解析的FBP算法获得所谓的先验图像，这极大地限制了该算法的应用范围；同时，当先验图像和待重建图像不完全对应时，由该算法得到的重建图像会出现伪影。随后，Debatin等[11-12]将各向异性全变差作为目标函数进行CT迭代图像重建，提出了ATV(Anisotropic Total Variation minimization)算法和GATV(Generalized Anisotropic Total Variation minimization)算法，在稀疏投影数据的情况下，ATV算法和GATV算法较经典的TV算法得到的重建图像边界更清晰、细节更丰富。
综合已有文献并与CT工程技术人员进行交流，将[0, 2π)范围内扇束/锥束扫描不超过20个投影角度称为极稀疏投影，并着力于从极稀疏投影数据足够精确地重建断层图像，从而能够在显著降低CT检查X-射线辐射剂量的前提下，提供充分适宜影像学临床诊断需求的重建图像。为此，提出了投影驱动的系统模型，并设计了一种将CT迭代图像重建与CS理论相结合的重建算法；为检验重建算法性能，分别对圆周扫描扇束和锥束投影得到的极稀疏投影数据进行了图像重建仿真实验；仿真实验结果表明，在极稀疏投影数据的条件下，算法数值精度高，计算复杂度低，内存开销少，有很强的工程实用性。
1 理论与方法 1.1 基于压缩感知的CT图像重建 CS理论指出，如果信号是稀疏可压缩的，那么就可以使用远低于Nyquist频率的抽样频率对原始信号进行抽样，并且能精确恢复出该信号。基于此，稀疏投影角度下的CT迭代图像重建，可以看成是求解如下不适定线性方程组问题[13]。
 | (1) |
式中：g为S维投影数据列矢量; f为T维待重建图像; W为S×T维的系统矩阵(即下文的正投影矩阵)，通常S?T。根据CS理论，如果系统矩阵W满足严格等距同构条件，那么式(1)可以转化为下述有约束的l1范数最优化问题:
 | (2) |
式中：Ψf为待重建图像f的稀疏化表示。
1.2 投影驱动的系统模型 在CT迭代图像重建中，为了确定系统矩阵W元素，必须建立合理的系统模型(即正/反投影运算模型)。系统模型对CT迭代图像重建的数值精度和重建图像质量都有着重要影响。目前，主流的系统模型主要分为3类：像素驱动模型、射线驱动模型和距离驱动模型。像素驱动模型和射线驱动模型分别容易在投影域和图像域引入高频伪影；而距离驱动模型则充分利用存储空间的顺序访问模式，算法复杂度相对较低，能够避免图像域伪影和投影域伪影[14-15]。结合像素驱动和射线驱动模型的优点，并基于距离驱动模型的基本思想，提出了如下投影驱动的系统模型(即投影驱动的正/反投影运算模型)。
1.2.1 投影驱动模型 图 1所示的二维扇束投影几何描述了投影驱动的正投影计算过程。在投影角度一定的情况下，将光源分别与某行像素边界的中点与检测器单元的边界相连接，得到它们各自在公共轴(见图 1中的X轴)上的交点；然后在X轴上计算像素交点与相邻检测器交点的重叠长度，并用归一化的重叠长度表示像素对检测器正投影的贡献；依次处理每一行像素。投影驱动的反投影计算也是如此。至此，投影驱动的正/反投影计算策略可概括为：一次迭代处理一个投影角度下的所有像素。
 |
| 图 1 二维扇束投影示意图与细节图 Fig. 1 Schematic diagram and detailed diagram of two-dimensional fan-beam projection |
| 图选项 |
图 1也详细地展示了检测器边界di、像素边界pi、检测器上投影值dij、像素值pij之间的关系。需要注意，投影驱动的正/反投影中归一化加权系数计算中的分母不同。正投影计算中，为了更精确地模拟线积分，将加权系数除以投影角度β的余弦值，投影驱动的正投影运算表示为
 | (3) |
反投影运算表示为
 | (4) |
1.2.2 正/反投影矩阵 在CT迭代图像重建中，正/反投影矩阵将图像域与投影域联系起来。而且，上述正/反投影运算模型(即系统模型)决定了正/反投影矩阵各元素的值。正投影运算g=Wf将图像从图像域映射到投影域，W为正投影矩阵(即前文所述的系统矩阵)；f和g的含义前文已提及，分别为待重建图像和投影数据列矢量。式(5)所示的反投影运算将投影数据从投影域映射到图像域，M为反投影矩阵。
 | (5) |
正/反投影矩阵W和M的结构一致，但矩阵元素不同。图 2为正/反投影矩阵元素示意图。结合图 2，进一步阐述正/反投影矩阵各元素的含义。
 |
| 图 2 正/反投影矩阵元素示意图 Fig. 2 Schematic diagram of matrix elements of forward/backward projection |
| 图选项 |
正投影矩阵W的结构如图 2所示。图中：ViewNum为投影角度数；R为像素行数；C为像素列数。Wθ为正投影矩阵W在第θ个投影角度的子矩阵；Wθ(r)为在第θ个投影角度下第x行像素对所有检测器单元的相对贡献值，矩阵元素为式(3)中的加权系数。式(6)和式(7)清楚地揭示了以上3个矩阵的关系。
 | (6) |
 | (7) |
式(5)表明反投影矩阵M的转置MT直接参与反投影运算，矩阵MT的结构如图 2所示。与正投影类似，MθT为矩阵MT在第θ个投影角度的子矩阵，Mθ(r)T为在第θ个投影角度下，全部检测器单元对第r行像素的相对贡献值，矩阵元素为式(4)中的加权系数。以上3个矩阵的关系还可以由式(8)和式(9)表示。
 | (8) |
 | (9) |
1.3 压缩感知的投影驱动CT图像重建 在CS理论框架下，综合传统CT迭代图像重建技术的优势，并利用上述投影驱动系统模型，设计了一种CT迭代图像重建算法，称作CS的投影驱动CT图像重建(Compressed Sensing View Driven CT image reconstruction，CSVD)算法，算法由基于投影驱动模型的粗略图像重建和最优化计算2个环节组成。
算法1??CSVD算法。
输入：??初始图像向量，迭代次数N，M，K，投影角度数。
输出：??重建图像。
过程：??
总体迭代???//总体迭代次数N
{
/*基于投影驱动模型的粗略图像重建*/
循环1??//粗略重建次数M
{
循环2??//投影角度数ViewNum
{
??投影驱动正投影
??作差
??投影驱动反投影
??校正
}
}
/*最优化计算*/?//最优化迭代次数K
{
??梯度计算得到梯度图像
??求最速下降一次的图像
}
}
1.3.1 基于投影驱动模型的粗略图像重建 基于投影驱动系统模型，迭代求解不适定的线性方程组g=Wf，以获得众多满足约束项的解。循环2还可以用图 3表示，每次循环依次处理每个投影角度下的投影数据；每个角度的投影数据都进行正投影、作差、反投影和校正操作。图中：Wθ和Mθ分别为上文提及的正/反投影矩阵的子矩阵；f(RR)[n, m, θ]为第n次总体迭代的第m次粗略重建子迭代中在第θ个投影角度下的图像向量；Pθ为第θ个投影角度的实际投影数据；φθ为校正因子。
 |
| 图 3 循环2 Fig. 3 Loop 2 |
| 图选项 |
1.3.2 最优化计算 使用梯度下降法最优化目标函数
 | (10) |
 | (11) |
式中：dA[n]为第n次总体迭代中计算得到的平衡因子; 第n次总体迭代的第k次最优化子迭代中的图像向量表示为f(GRAD)[n, k]；α为步长因子。
2 仿真实验与结果分析 为研究CSVD算法的应用性能，在不同极稀疏投影数据下进行了2组数值仿真实验：圆周扫描扇束投影的二维CT图像重建、圆周扫描锥束投影的三维CT图像重建。上述CT扫描的几何布局如图 4所示，另外，综合考虑参数设置对成像质量等多方面的影响，扫描参数的设置以及仿真模型的信息如表 1所示。需要注意，圆周轨迹的锥束极稀疏投影扫描方式显然不满足Tuy条件[16]，故此算法属于近似的锥束CT图像重建算法。
 |
| 图 4 圆周CT扫描系统结构示意图 Fig. 4 Schematic diagram of circular CT scanning system |
| 图选项 |
表 1 CT扫描参数和仿真模型信息 Table 1 CT scanning parameters and simulation model information
| 仿真参数 | 二维扇束 | 三维锥束 |
| 检测器类型 | 平板检测器 | 平板检测器 |
| 行/列检测器数 | 512 | 256 |
| 行/列检测器中心 | 255.5 | 127.5 |
| 行/列检测器间距/mm | 1.0 | 1.0 |
| 光源到旋转中心距离/mm | 640 | 640 |
| 旋转中心到检测器中心距离/mm | 640 | 640 |
| 仿真模型 | 二维Shepp-Logan | 三维Shepp-Logan模型 |
| 模型尺寸：512×512 | 模型尺寸：256×256×256 | |
| 模型像素尺寸：0.5mm×0.5mm | 体素尺寸：0.5mm×0.5mm×0.5mm |
表选项
为在极稀疏投影条件下获得更多的投影信息，对投影角度θ进行如下设置(ViewgNum数值为偶数)：
 | (12) |
2.1 极稀疏投影数据下的二维扇束图像重建 仿真实验先通过对图 5(a)所示的Shepp-Logan模型扫描获得投影数据，再使用CSVD算法对该模型在不同投影角度数下进行二维CT图像重建，重建结果如图 5(b)所示。为了更精确地比较重建图像与原始图像的像素值，图 5(c)给出两者在水平方向上图像的中心像素值分布。重建图像和像素值分布曲线显示，36和20个投影角度下的重建图像与参考图像几乎相同、像素曲线几乎重合；18个投影角度下的重建图像部分边缘模糊且部分区域有块状伪影，但是各个组织结构都能清晰分辨，像素值分布曲线在像素值跳变部位有冲激，其他部分与参考图像曲线几乎重合。
 |
| 图 5 Shepp-Logan模型和重建图像及其中心行像素值比较 Fig. 5 Shepp-Logan model and reconstructed images, as well as comparison of pixel values of center row between them |
| 图选项 |
为进一步检验CSVD算法在极稀疏投影数据下的重建性能，使用标准均方根误差(Normalized Root Mean Squared Error，NRMSE)、峰值信噪比(Peak Signal to Noise Ratio，PSNR)、全局质量指标(Universal image Quality Index，UQI)[17]3种衡量标准对3种算法(FBP算法、ART-TV算法和CSVD算法)的重建质量进行客观评价，结果如图 6所示。实验结果表明：CSVD算法在各项指标上明显优于其他2种算法，说明CSVD算法在重建图像精确度、噪声抑制和图像恢复程度等方面有着更好的性能。
 |
| 图 6 各个投影角度数下不同算法的重建图像质量对比 Fig. 6 Comparison of reconstructed image quality of different algorithms under various projection angles |
| 图选项 |
最后给出以上3种算法在20个投影角度下的3种衡量指标具体数值以及重建时间，如表 2所示。与FBP算法相比，CSVD算法耗费了较多的时间，但是CSVD算法在不完全投影数据下的重建性能明显优于FBP算法。与ART-TV算法相比，CSVD算法的重建效率提高了2倍多，同时从以上3种衡量指标可以看出CSVD算法的重建质量有明显提升。
表 2 各个投影角度数下不同算法的重建结果 Table 2 Reconstruction results of different algorithms under various projection angles
| 重建算法 | NRMSE | PSNR | UQI | 重建时间/s |
| FBP | 0.9641 | 12.4663 | 0.5222 | 4.094 |
| ART-TV | 0.3387 | 21.5526 | 0.9181 | 69.585 |
| CSVD | 0.1958 | 26.3113 | 0.9734 | 32.997 |
表选项
2.2 极稀疏投影数据下的三维锥束图像重建 与二维扇束图像重建实验类似，首先对自定义的三维Shepp-Logan模型扫描获得投影数据，该模型透视图如图 7(a)所示；其次，使用CSVD算法对该模型分别在20和18个投影角度数下进行三维图像重建，获得中间6层(125层至130层)的横断面切片，选取其中的第125层和128层的重建结果作为展示，如图 7所示，其中图 7(b)、图 7(c)分别为20和18个投影角度下的重建图像，同时还比较了第128层重建图像与参考图像在水平方向上图像中心像素值，得到的像素分布曲线与图 5(c)类似，这里不再列出；最后，为了更客观地评价重建图像的质量，使用NRMSE、PSNR、UQI三种衡量标准对第128层重建图像质量进行分析如表 3所示。
 |
| 图 7 Shepp-Logan模型透视图和重建图像(第125层和128层) Fig. 7 Perspective of Shepp-Logan model and reconstructed images(125th and 128th layers) |
| 图选项 |
表 3 CSVD算法重建图像质量评价 Table 3 Quality evaluation of reconstructed images using CSVD algorithm
| 质量评价指标 | 20个投影角度 | 18个投影角度 |
| NRMSE | 0.2887 | 0.2980 |
| PSNR | 22.9883 | 22.7137 |
| UQI | 0.9406 | 0.9362 |
表选项
重建结果表明：在极稀疏投影数据的条件下，CSVD算法表现出良好的重建性能。在仅有18个投影角度的前提下，除重建图像边缘较模糊外，各个组织结构仍能清晰分辨，像素值分布曲线重合度高，且在NRMSE、PSNR、UQI这3项指标上与20个投影角度相当。
3 结论 研究如何在极稀疏投影数据情况下进行CT迭代图像重建具有重要的临床意义。理论分析与仿真实验结果都表明，CSVD算法能够从极稀疏投影数据足够精确地重建断层图像，从而将能够在显著降低CT检查X-射线辐射剂量的前提下，提供充分适宜影像学临床诊断需求的重建图像。CSVD算法良好的图像重建性能主要体现在：
1) 数值精度高。仿真实验结果表明：极稀疏投影角度数下的重建图像精确地再现了参考图像的图像结构及像素分布，另外CSVD算法的各项图像质量衡量指标明显优于ART-TV算法和FBP算法。
2) 计算复杂度低。由前文的理论分析可知，CSVD算法一次迭代会处理一个投影角度，且在一个投影角度下，只需一次迭代就可以获得一行像素对所有检测器单元的贡献(正投影过程)，减少了大量不必要的遍历运算，使得计算复杂度大幅度降低。
3) 内存开销少。在医学成像应用中，系统矩阵的规模通常都很大，致使基于其他系统模型的CT迭代图像重建一般都需要存储系统矩阵，而基于投影驱动系统模型的图像重建则不需要对系统矩阵进行存储，大大减小了内存的开销，CSVD算法有着很强的工程实用性。
总之，在极稀疏投影数据的条件下，CSVD算法数值精度高，计算复杂度低，内存开销少，有很强的工程实用性。这给相关领域的科研工作者在极稀疏投影数据情况下进行CT迭代图像重建(从而可以降低CT检查的X-射线辐射剂量)提供了一条切实的技术途径。
参考文献
| [1] | SABA L, SURI J S. Multi-detector CT imaging principles, head, neck, and vascular systems[M]. Boca Raton: CRC Press, 2013: 103-122. |
| [2] | OHNESORGE B M, FLOHR T G, BECKER C R, et al. Multi-slice CT technology[M]. Berlin: Springer, 2007: 41-69. |
| [3] | NORELLI L J, COATES A D, KOVASZNAY B M. Cancer risk from diagnostic radiology in a deliberate self-harm patient[J]. Acta Psychiatrica Scandinavica, 2010, 122(5): 427-430. |
| [4] | RAMPINELLI C, DE M P, ORIGGI D, et al. Exposure to low dose computed tomography for lung cancer screening and risk of cancer:Secondary analysis of trial data and risk-benefit analysis[J]. British Medical Journal, 2017, 356(347): 8. |
| [5] | PEARCE M S, SALOTTI J A, LITTLE M P, et al. Radiation exposure from CT scans in childhood and subsequent risk of leukaemia and brain tumors:A retrospective cohort study[J]. Lancet, 2012, 380(9840): 499-505. |
| [6] | MEHYAR L S, ABU-ARJA M H, STANEK J R, et al. The risk of developing secondary central nervous system tumors after diagnostic irradiation from computed tomography in pediatrics:A literature review[J]. Pediatric Neurology, 2019, 98: 18-24. |
| [7] | SIDKY E Y, KAO C M, PAN X. 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. |
| [8] | SIDKY E Y, PAN X. Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization[J]. Physics in Medicine and Biology, 2008, 53(17): 4777-4807. |
| [9] | DONOHO D L. Compressed sensing[J]. IEEE Transactions on Information Theory, 2006, 52(4): 1289-1306. |
| [10] | CHEN G H, TANG J, LENG S. Prior image constrained compressed sensing (PICCS):A method to accurately reconstruct dynamic CT images from highly under-sampled projection data sets[J]. Medical Physics, 2008, 35(2): 660-663. |
| [11] | DEBATIN M, ZYGMANSKI P, STSEPANKOU D, et al.CT reconstruction from few-views by anisotropic total variation minimization[C]//2012 IEEE Nuclear Science Symposium and Medical Imaging Conference Record.Piscataway: IEEE Press, 2012: 2295-2296. |
| [12] | DEBATIN M, HESSER J. Accurate low-dose iterative CT reconstruction from few projections by generalized anisotropic total variation minimization for industrial CT[J]. Journal of X-ray Science and Technology, 2015, 23(6): 701-726. |
| [13] | CANDES E, ROMBERG J, TAO T. Stable signal recovery from incomplete and inaccurate measurements[J]. Communications on Pure and Applied Mathematics, 2006, 59(8): 1207-1223. |
| [14] | MAN B D, BASU S.Distance-driven projection and backprojection[C]//2002 IEEE Nuclear Science Symposium Conference Record.Piscataway: IEEE Press, 2002: 1477-1480. |
| [15] | MAN B D, BASU S. Distance-driven projection and backprojection in three dimensions[J]. Physics in Medicine and Biology, 2004, 49(11): 2463-2475. |
| [16] | TUY H K. An inversion formula for cone-beam reconstruction[J]. SIAM Journal on Applied Mathematics, 1983, 43(3): 546-552. |
| [17] | WANG Z, BOVIK A C. A universal image quality index[J]. IEEE Signal Processing Letters, 2002, 9(3): 81-84. |
