基于双向流场引导投影补全的稀疏角度锥束CT重建算法

A sparse-view cone-beam CT reconstruction algorithm based on bidirectional flow field- guided projection completion

摘要

目的 为解决稀疏角度锥束CT成像面临的不适定反问题求解问题，提出基于双向流场引导投影补全的稀疏角度锥束CT重建算法（BBC-Recon）。 方法 BBC-Recon方法包含两个主要模块：投影补全模块和图像恢复模块。投影补全模块基于流场估计的思想，通过设计的双向和多尺度关联体，充分计算投影之间的相关性信息和冗余信息，用以精确指导双向流场和缺失帧的生成，实现对缺失投影的高精度补全，获得伪完备投影；图像恢复模块对获得的伪完备投影进行重建，然后对图像进行细化，去除伪影残留，进一步提升图像质量。 结果 在Mayo和桂林医学院公开数据集上的实验结果表明，对比现有算法，BBC-Recon方法在稀疏4倍角度的情形下，较次优方法：PSNR指标提升1.80%，SSIM指标提升0.29%，RMSE指标降低4.12%；在稀疏8倍角度的情形下，较次优方法：PSNR指标提升1.43%，SSIM指标提升1.49%，RMSE指标降低0.77%。 结论 BBC-Recon充分挖掘了投影之间的相关性信息，不仅能在保持图像结构信息的前提下有效去除条纹伪影，而且在层间一致性的保持上发挥了极大优势。

Abstract

Objective We propose a sparse-view cone-beam CT reconstruction algorithm based on bidirectional flow field guided projection completion (BBC-Recon) to solve the ill-posed inverse problem in sparse-view cone-beam CT imaging. Methods The BBC-Recon method consists of two main modules: the projection completion module and the image restoration module. Based on flow field estimation, the projection completion module, through the designed bidirectional and multi-scale correlators, fully calculates the correlation information and redundant information among projections to precisely guide the generation of bidirectional flow fields and missing frames, thus achieving high-precision completion of missing projections and obtaining pseudo complete projections. The image restoration module reconstructs the obtained pseudo complete projections and then refines the image to remove the residual artifacts and further improve the image quality. Results The experimental results on the public datasets of Mayo Clinic and Guilin Medical University showed that in the case of a 4-fold sparse angle, compared with the suboptimal method, the BBC-Recon method increased the PSNR index by 1.80% and the SSIM index by 0.29%, and reduced the RMSE index by 4.12%; In the case of an 8-fold sparse angle, the BBC-Recon method increased the PSNR index by 1.43% and the SSIM index by 1.49%, and reduced the RMSE index by 0.77%. Conclusion The BBC-Recon algorithm fully exploits the correlation information between projections to allow effective removal of streak artifacts while preserving image structure information, and demonstrates significant advantages in maintaining inter-slice consistency.

Graphical abstract

关键词

稀疏角度锥束CT / 双向多尺度关联性 / 双向流场估计 / 双域训练

Key words

sparse-view cone-beam CT / bidirectional multi-scale correlation / bidirectional flow field estimation / dual-domain training

锥束计算机断层扫描成像作为一种3D成像技术，因其具有的高空间分辨率、快速扫描速度优势，已被广泛用于物体内部结构的检查^［1］。然而，由辐射暴露所导致的公共卫生健康风险越发地受到人们关注^{［2， 3］}，在保证图像质量的前提下尽可能地降低辐射剂量以减小对人体的辐射损失风险成为了目前CT成像研究中颇受关注的话题。稀疏角度扫描通过对常规扫描的角度数进行稀疏，即减少球管在旋转一圈的扫描中曝光的次数，可以有效降低患者所受辐射剂量，进而降低辐射损伤发生的风险。但由于稀疏角度扫描产生的投影数据为不完备数据，使用传统的FDK解析重建算法^［4］会导致在图像中呈现出明显的条纹伪影，易掩盖或改变物体结构，对医学诊断产生误导甚至误诊。为获得高质量的CT图像，研究者们在抑制条纹伪影的同时并尽可能保持物体结构细节不发生改变这一任务上开展了大量研究，方法大致可分为3类：投影域方法、图像域方法、双域方法。

在投影域方法中，投影恢复又称弦图恢复，目的是通过不同手段对缺失角度的投影数据进行恢复以获得伪完备数据。常见的使用TV约束的迭代优化策略^{［5， 6］} 、字典学习策略^［7］、方向插值及其变体来完成缺失投影的恢复^［8-10］，虽然可以获得足量的投影数据来满足FDK重建，但是数据库构建的完备性和多样性差、插值的误差大等问题极易导致结构丢失以及二次伪影的产生。图像域处理方法又称后处理方法，即直接对重建的图像进行噪声和伪影的抑制。早期通过非局部均值算法去伪影^［11］，但是图像噪声统计特性和伪影结构较为复杂，效果不佳；后得益于深度学习^{［12， 13］}（DL）技术的发展和U-Net^［15］的提出，FBPConvNet^［14］、REDCNN^［16］这类经典方法被提出，以及解纠缠策略也被应用稀疏角度中^［17］，这类方法直接从图像入手，虽然推理速度快、去伪影干净，但是图像结构平滑、模糊问题凸显。由于基于投影或图像的单域方法都存在一定缺陷，双域结合成为了必然趋势。双域方法是指将投影域与图像域结合，或是先进行投影的恢复，然后使用图像域网络进行细化，亦或是将投影域与图像域进行联合训练。较为直接的是Sparsier2Sparse方法^［22］，将稀疏投影再次稀疏来获得配对数据；较为复杂的是多阶段训练策略^［18-21］、多域transformer架构^［23］和迭代优化策略^［24］，通过设计的多个模块或者多次迭代来提升图像质量、保持数据一致性，而近年来的研究中，扩散模型类、神经辐射场^［25］（NeRF）类双域方法也被提出^［26-29］，用以解决稀疏角度问题。目前这些双域方法虽不同程度地提高了图像质量，但也具有一定局限性：一方面，现有的投影插值未能或未能较好的利用投影之间的几何关系和冗余信息，导致生成的投影连续性较差或者准确性不够，难以在结构细节信息的保留上发挥作用；另一方面，现有方法多聚焦于2D成像问题，例如平行束或扇形束，难以直接针对锥束进行处理，而扩散模型方法的生成不确定性导致结构丢失、噪点多，NeRF方法面对不同的case都需要重新训练、泛化性差。

在最新的稀疏角度研究中，基于流场估计的方法^［30-32］也崭露头角，它通过估计两帧图像之间像素的单向/双向流动，进行中间帧的合成，或用于图像恢复，或用于视频合成，可以获得非常好的效果。其最新应用是在DSA新视图合成任务中^［30］，用以实现缺失投影的恢复。考虑到现有的投影补全方法多数仅应用简单的插值策略，缺乏对投影之间相关性和冗余信息的有效利用，不可避免地在补全过程中引入较多错误信息，产生二次伪影；而自然图像处理中这种流场估计的手段基于前后两帧图像，可充分挖掘图像之间的相关性，其思想与挖掘两帧投影之间的相关性和冗余信息、提高投影补全精准度的思想非常契合，故而基于此，我们提出了基于双向流场引导投影补全的稀疏角度锥束CT重建算法（BBC-Recon）。BBC-Recon主要由两个模块构成：投影补全模块和图像恢复模块，在投影补全模块中，应用双向流场进行引导，首先，对两帧投影进行计算，获得一个包含了投影相关性和冗余信息的双向多尺度关联体，其次，使用编码器对两帧投影进行编码，获取双向流场和缺失帧的初值，随后，采取尺度查询策略从计算的关联体中查询相关性，对双向流场和缺失帧进行更新，最后通过扭曲得到缺失帧，形成伪完备投影；在图像恢复模块中，对伪完备投影进行重建后，采用一个2D网络进行图像细化，进一步提升图像质量，最终获得高质量的锥束CT图像。

1 资料和方法

1.1 方法概述

BBC-Recon方法的整体框架（图1），本研究的方法主要分为两部分：投影补全模块、图像恢复模块，方法流程是，使用构建的基于双向流场估计的3D投影补全模块对缺失投影进行插值，可获得连续性好、准确性高的伪完备投影，再使用FDK算法进行重建得到初始恢复的CT图像，最后使用一个轻量的图像域网络FBPConvNet^［14］，结合稀疏投影直接重建的伪影图像，进行图像域的细化，得到结构细节高保留、伪影抑制干净的CT图像。

1.2 投影域恢复

区别于单向流场估计，本研究基于全对双向流场估计^［31］，对两帧投影之间的双向流场进行估计，利用投影之间的几何信息，实现对投影间相关性的充分挖掘，提高了插值的准确性和插值投影之间的连续性（图2）。

该模块主要由3部分构成，分别对应3个任务：（1）双向多尺度关联体生成：通过一个特定的关联编码器生成投影之间的相关性，用于指导流场和插值的中间帧的更新；（2）流场与中间帧的初始化：使用内容编码器对输入投影进行编码，获得网络训练的初值——双向流场和缺失的中间帧，并生成已有帧的特征图；（3）关联查询与更新：使用尺度查询和update模块，查询关联性，对双向流场和插值的中间帧进行更新。后续内容介绍均以稀疏4倍角度情形下的恢复为例，即从已有的两帧投影恢复出中间缺失的3帧投影。

1.2.1 双向多尺度关联体生成

投影域恢复模块的设计包含两类编码器（图2）：相关性编码器和内容编码器，此处我们采用前者，给定已知的两帧投影\{{\mathrm{I}}_{\mathrm{0}},{\mathrm{I}}_{\mathrm{1}}\}\in R^{H\times W}，通过相关性编码器进行多次卷积操作获取特征图\{{\mathrm{g}}_{\mathrm{0}},{\mathrm{g}}_{\mathrm{1}}\}\in R^{\frac{H}{\mathrm{8}}\times \frac{W}{\mathrm{8}}}，然后对其应用矩阵乘法，可得到一个特征之间的关联性C_{\mathrm{1}}（图3），具体计算如下：

C_{{\mathrm{1}}_{ijkl}}={\sum }_h{g}_{{\mathrm{0}}_{ijh}\bullet }{g}_{{\mathrm{1}}_{klh}},C_{\mathrm{1}}\in R^{\frac{H}{\mathrm{8}}\times \frac{W}{\mathrm{8}}\times \frac{H}{\mathrm{8}}\times \frac{W}{\mathrm{8}}}，（1）

为了进一步获得多尺度水平的关联性，我们对C_{\mathrm{1}}的后两个维度应用平均池化操作，可获得两个不同尺度水平下的关联性C_{\mathrm{2}},C_{\mathrm{3}}：

C_{\mathrm{2}}=pool\left(C_{\mathrm{1}}\right),C_{\mathrm{2}}\in R^{\frac{H}{\mathrm{8}}\times \frac{W}{\mathrm{8}}\times \frac{H}{\mathrm{16}}\times \frac{W}{\mathrm{16}}}，（2）

C_{\mathrm{3}}=pool\left(C_{\mathrm{2}}\right),C_{\mathrm{3}}\in R^{\frac{H}{\mathrm{8}}\times \frac{W}{\mathrm{8}}\times \frac{H}{\mathrm{32}}\times \frac{W}{\mathrm{32}}}，（3）

最后，我们便得到了一个包含3个尺度水平的单向特征图关联体C=\left\{C_{\mathrm{1}},C_{\mathrm{2}},C_{\mathrm{3}}\right\}，但考虑到C仅反应了从{\mathrm{g}}_{\mathrm{0}}到{\mathrm{g}}_{\mathrm{1}}的多尺度正向对应关系，缺乏从{\mathrm{g}}_{\mathrm{1}}到{\mathrm{g}}_{\mathrm{0}}的反向对应，为充分挖掘帧之间的关联性，需要再次实施矩阵乘法，以获取从{\mathrm{g}}_{\mathrm{1}}到{\mathrm{g}}_{\mathrm{0}}的反向多尺度对应关系。值得注意的是，我们已经计算出了正向关联体C，只需要对其进行转置便可以计算出反向对应关系即关联体C^T=\left\{C_{\mathrm{1}}^T,C_{\mathrm{2}}^T,C_{\mathrm{3}}^T\right\}，从而避免多次计算增加计算量。最终，我们可以得到双向多尺度关联体C=\{C,C^T\}，上式中，H、W分别表示输入投影的高和宽，pool为平均池化操作，(i,j)为{\mathrm{g}}_{\mathrm{0}}的坐标索引，(k,l)为{\mathrm{g}}_{\mathrm{1}}的坐标索引。

1.2.2 流场与中间帧的初始化

流场是指两帧之间像素的移动，如果能将现有帧和缺失帧之间的像素移动刻画出来，便可以利用流场结合现有帧通过扭曲得到缺失帧，所以流场是网络主要学习的目标。图4A显示的是估计的双向流场，使用RGB色彩进行了可视化，中间为5帧投影，左侧为中间帧到0帧的流场，右侧为中间帧到1帧的流场，图4B为颜色图例，不同颜色代表不同移动方向，颜色深浅与移动程度成正比，图4A中对流场的可视化有效刻画了像素的移动，不难看出，随着帧间距离的增加，像素的移动程度也在逐渐增大。

正式学习之前，需要进行初始化，通过内容编码器Context encoder（图5），初始化缺失帧 I_{t_{i=\mathrm{1,2},\mathrm{3}}} 的特征X_{t_i}^{\mathrm{1}}以及现有帧 \{I_{\mathrm{0}},I_{\mathrm{1}}\}与 I_{t_{i=\mathrm{1,2},\mathrm{3}}}之间的流场\{F_{t_i\to \mathrm{0}}^{\mathrm{1}},F_{t_i\to \mathrm{1}}^{\mathrm{1}}\}，获得网络训练的初值，其中t表示缺失帧，下标i表示某一缺失帧的具体位置。除此之外，内容编码器还生成不同水平的特征图\left(X_{\mathrm{0}}^{l=\mathrm{1,2},\mathrm{3}},X_{\mathrm{1}}^{l=\mathrm{1,2},\mathrm{3}}\right)，构成解码器输入的一部分。

1.2.3 关联查询与更新

通过1.2.1和1.2.2操作，我们得到了已有两帧投影之间的双向多尺度关联性以及初始化的流场和待恢复的中间帧，接下来需要从关联体中查询相关性来更新流场和中间帧，最终合成新视图，具体分3步进行：

⑴尺度相关性查询：在获得双向关联体\mathit{ℂ}后，我们需要根据初始化的双边流场\{F_{t_i\to \mathrm{0}}^{\mathrm{1}},F_{t_i\to \mathrm{1}}^{\mathrm{1}}\}从 \mathit{ℂ}中查询相关特征图。这里参考RAFT方法^［31］中的查询策略，但值得注意的是，RAFT估计的流场是\left\{F_{\mathrm{1}\to \mathrm{0}}^{\mathrm{1}},F_{\mathrm{0}\to \mathrm{1}}^{\mathrm{1}}\right\}，即从0帧到1帧的映射，由于其生成的关联体\mathit{ℂ}是关于0帧与1帧的，所以能够直接实施查询操作，而我们估计的流场是0帧与t_i帧或t_i帧与1帧之间的映射，无法直接对\mathit{ℂ}实施查询操作，针对这样的不匹配问题，我们借鉴AMT方法^［33］的思路，对流场进行尺度变换来得到\left\{F_{\mathrm{1}\to \mathrm{0}}^{\mathrm{1}},F_{\mathrm{0}\to \mathrm{1}}^{\mathrm{1}}\right\}：

F_{\mathrm{0}\to \mathrm{1}}^{\mathrm{1}}=\frac{\mathrm{1}}{\mathrm{1}-t_i}{F}_{t_i\to \mathrm{1}}^{\mathrm{1}},F_{\mathrm{1}\to \mathrm{0}}^{\mathrm{1}}=\frac{\mathrm{1}}{t_i}{F}_{t_i\to \mathrm{0}}^{\mathrm{1}},t_i\in \left\{\frac{i}{\lambda }\right\}（4）

式中\lambda为稀疏强度，我们采用\lambda =\mathrm{4,8}。随后，采用RAFT查询策略，将X_{\mathrm{0}}^l中每一个像素x=\left(u,v\right)映射到X_{\mathrm{1}}^l中：x^{\text{'}}=\left(u+F_{\mathrm{1}\to \mathrm{0}}^{\mathrm{1}}\left(u\right),v+F_{\mathrm{0}\to \mathrm{1}}^{\mathrm{1}}\left(v\right)\right)，接着便可以定义一个局部网格：

N{\left(x^{\text{'}}\right)}_r=\left\{x^{\text{'}}+dx|dx\in Z^{\mathrm{2}},{‖dx‖}_{\mathrm{1}}\le r\right\}，（5）

N{\left(x^{\text{'}}\right)}_r作为一个整数偏移的集合，将用于从\mathit{ℂ}中进行索引，并采取双线性插值进行采样。值得注意的是，该查询操作将会在\mathit{ℂ}中的所有尺度水平进行，查询到的双向关联将会被级联成特征图Corr，用于下一步对双向流场和插值的中间特征的更新。

⑵流场与特征更新：使用Update模块（图5右）基于提取的双向关联性Corr，预测双边流场\{F_{t_i\to \mathrm{0}}^{\mathrm{1}},F_{t_i\to \mathrm{1}}^{\mathrm{1}}\}的残差∆F_{t_i\to \mathrm{0}}^{\mathrm{1}},∆F_{t_i\to \mathrm{0}}^{\mathrm{1}}以及插值的中间帧X_{t_i}^{\mathrm{1}}的残差∆X_{t_i}^{\mathrm{1}}。根据估计的残差可以更新双向流场和插值的中间帧：

{\widehat{X}}_{t_i}^{\mathrm{1}}=X_{t_i}^{\mathrm{1}}+∆X_{t_i}^{\mathrm{1}},{\widehat{F}}_{t_i\to \mathrm{0}}^{\mathrm{1}}=F_{t_i\to \mathrm{0}}^{\mathrm{1}}+∆F_{t_i\to \mathrm{0}}^{\mathrm{1}},{\widehat{F}}_{t_i\to \mathrm{1}}^{\mathrm{1}}=F_{t_i\to \mathrm{1}}^{\mathrm{1}}+∆F_{t_i\to \mathrm{1}}^{\mathrm{1}}，（6）

利用更新的双向流场，对输入特征X_{\mathrm{0}}^l,X_{\mathrm{1}}^l进行扭曲得到{\widehat{X}}_{\mathrm{0}}^{\mathrm{1}},{\widehat{X}}_{\mathrm{1}}^{\mathrm{1}}，然后将上述更新的双向流场、中间帧以及扭曲特征输入到解码器D^{\mathrm{1}}中，得到下一阶段的输入：

\left[F_{t_i\to \mathrm{0}}^{\mathrm{2}},F_{t_i\to \mathrm{1}}^{\mathrm{2}},X_{t_i}^{\mathrm{2}}\right]=D^{\mathrm{1}}\left(\left[{\widehat{X}}_{\mathrm{0}}^{\mathrm{1}},{\widehat{X}}_{\mathrm{1}}^{\mathrm{1}},{\widehat{F}}_{t_i\to \mathrm{0}}^{\mathrm{1}},{\widehat{F}}_{t_i\to \mathrm{1}}^{\mathrm{1}},{\widehat{X}}_{t_i}^{\mathrm{1}}\right]\right)（7）

D^{l=\mathrm{1,2},\mathrm{3}}使用的是IFRNet^［32］的解码器模块IFRBlock（图5左下）。

⑶新视图合成：在通过最后一个解码器D^{\mathrm{3}}时，我们将得到双向流场F_{t\to \mathrm{0}},F_{t\to \mathrm{1}}，根据流场生成的掩膜M、以及插值帧的残差R，再结合输入的{\mathrm{I}}_{\mathrm{0}},{\mathrm{I}}_{\mathrm{1}}，便可以插值出所缺失的中间帧：

I_{t_i}=M\odot W\left(I_{\mathrm{0}},F_{t_i\to \mathrm{0}}\right)+\left(\mathrm{1}-M\right)\odot \left(I_{\mathrm{1}},F_{t_i\to \mathrm{1}}\right)+R（8）

\mathit{𝒲}定义为扭曲操作，I_{t_i}表示恢复的所有缺失帧，以4倍稀疏为例的话，则表示为缺失的中间3帧投影。至此，投影恢复模块便完成了缺失投影的恢复，通过生成新视图得到了伪完备投影数据，可直接用于FDK重建得到CT图像。

1.3 图像域细化

由于投影恢复模块通过估计双向流场和双向多尺度关联体，充分考虑了投影之间的相关性信息，利用该信息可以实现较为准确的缺失投影恢复，获得高质量的伪完备投影数据，能够较为彻底的去除条纹伪影。但是，区别于自然图像，投影之间的变化单依靠流场并不能充分刻画，单投影域恢复不可避免的会产生偏差，一定程度上引入二次伪影，影响图像质量。为此，在获得伪完备投影数据进行重建后，我们特地引入了FBPConvNet网络^［14］作为图像域恢复方法，通过将稀疏投影直接重建的图像与恢复的伪完备重建的图像同时作为网络输入，用于网络训练，对图像进行细化处理，可以有效去除二次伪影，这为获得高细节保留、伪影抑制干净的CT图像提供了强有力的保障。

1.4 损失函数

投影恢复网络损失函数由两部分构成：流动蒸馏损失{\mathrm{ℒ}}_{flow}^［32］和Charbonnier损失{\mathrm{ℒ}}_{char}^［34］，其中{\mathrm{ℒ}}_{flow}可以确保预测的流场尽可能的精准，{\mathrm{ℒ}}_{char}作为Charbonnier损失的拓展，构建在投影恢复中，用于确保插值的中间帧投影尽可能的接近真值，具体计算如下：

\begin{array}{l}{\mathrm{ℒ}}_{flow}=\\ \frac{\mathrm{1}}{\mathrm{H}\bullet \mathrm{W}}{\sum }_{i=\mathrm{1}}^{\mathrm{H}}{\sum }_{j=\mathrm{1}}^{\mathrm{W}}{\left({\left({\widehat{F}}_{t_i\to k,i,j}^{\mathrm{1}}-F_{t\to k,i,j}^{gt}\right)}^{\mathrm{2}}+{{\mathrm{10}}^{-\left(\mathrm{10}*{\omega }_k-\mathrm{1}\right)/\mathrm{3}}}^{\mathrm{2}}\right)}^{\frac{{\omega }_k}{\mathrm{2}}}\end{array}（9）

{\mathrm{ℒ}}_{char}=\frac{\mathrm{1}}{\mathrm{H}\bullet \mathrm{W}}{\sum }_{i=\mathrm{1}}^{\mathrm{H}}{\sum }_{j=\mathrm{1}}^{\mathrm{W}}{\left({\left({\widehat{I}}_{t,i,j}^{}-I_{t,i,j}^{gt}\right)}^{\mathrm{2}}+\tau \right)}^{\mathrm{1}/\mathrm{2}}（10）

式（9）中，参数{\omega }_k=\mathrm{e}\mathrm{x}\mathrm{p}\left(-\beta {\left|{\widehat{F}}_{t\to k}-F_{t\to k}^{gt}\right|}_{epe}\right)用于调整多尺度上的损失鲁棒性，其中\beta控制鲁棒性灵敏度，设为0.3，操作{\left|\bullet \right|}_{epe}用于计算像素端点误差，{\widehat{F}}_{t\to k}是估计的流场，F_{t\to k}^{gt}是根据完备投影预先计算的流场真值。式（10）中{\widehat{I}}_{t,i,j}^{}是插值的中间帧投影，I_{t,i,j}^{gt}是真实的中间帧投影，\tau是作为极小常数1e-6，避免计算结果为0。最后得到投影域网络训练的完整损失函数{\mathrm{ℒ}}_{proj}：

{\mathrm{ℒ}}_{proj}={\lambda }_{char}\bullet {\mathrm{ℒ}}_{char}+{\lambda }_{flow}\bullet {\mathrm{ℒ}}_{flow}，（11）

{\lambda }_{char}、{\lambda }_{flow}是对应损失的权重，设为1.0、0.002。

图像域网络（FBPConvNet）损失函数采用均方误差损失：

{\mathrm{ℒ}}_{image}={\mathrm{ℒ}}_{MSE}=\frac{\mathrm{1}}{\mathrm{H}\bullet \mathrm{W}}{\sum }_{i=\mathrm{1}}^H{\sum }_{j=\mathrm{1}}^W{\left({\widehat{I}}_{i,j}^{}-I_{i,j}^{gt}\right)}^{\mathrm{2}}（12）

式中{\widehat{I}}_{i,j}^{}、I_{i,j}^{gt}分别表示网络输出的预测图像、真值图像。

1.5 实验设计

1.5.1 实验数据

本文选用了两个数据集：一个是Mayo公开数据集Low Dose CT（LDCT） Image and Projection^［35］中的数据，该数据是西门子医疗公司的Somatom Definition AS+CT系统在自动曝光模式下扫描获得的，共选取94例胸部和腹部图像，其中90例用于训练，4例用于测试；另一个是桂林医学院公开的口腔CBCT数据集^［36］，该数据是郎视医疗公司的LargeV Smart3D-X设备对真实患者扫描所获得的，均为图像数据，共选取190例用于训练，5例用于测试。对于上述两个数据集，我们均使用全剂量图像作为真值图像进行锥束仿真，每个患者的图像尺寸为512×512×64，图像像素大小为0.56415 mm×0.56415 mm，口腔图像层厚为0.25 mm、胸部图像层厚为1.5 mm、腹部图像层厚为5 mm。探测器阵列规格为64×512，排数为64，通道数为512，探测器单元大小为1.21 mm×1.21 mm。完备投影扫描角度数设为600，对应投影大小为64×512×600，其中源-探测器距离（DSD）为1085.6 mm，源-旋转中心距离（DSO）为595 mm，共仿真了4倍和8倍两类稀疏强度，对应稀疏的投影角度数分别为150、75，投影大小分别为64×512×150、64×512×75，仿真使用ASTRA工具包完成。

1.5.2 实施细节

训练主要分为投影域、图像域两阶段。第1阶段，对完备投影进行稀疏采样，然后将稀疏-完备投影配对数据输入投影域网络进行训练，完成缺失投影的恢复，得到伪完备投影数据；第2阶段，使用ASTRA重建包对上一阶段得到的伪完备投影进行重建，将伪影-真值图像配对数据输入到2D的FBPConvNet网络进行训练，完成对图像的细化。其中，投影域网络的训练epoch数为350，采用AdamW优化算法，学习率为2e-4，每个epoch遍历所有患者，使用一张NVIDIA RTX A6000显卡完成。图像域网络的训练epoch数为100，采用Adam优化算法，学习率为1e-3，每个epoch迭代1000次，实际上50epoch左右便达到最优，使用一张NVIDIA RTX 4090显卡完成。

1.5.3 对比方法

本文共选用了6个对比方法，其中包含2个传统方法：FDK直接重建算法^［4］、基于TV约束的迭代重建方法ADS-POCS^［6］；4个深度学习方法：基于U-Net结构进行恢复的图像域方法FBPConvNet^［14］、基于CNN的投影插值联合图像恢复的双域方法Dual-CNN^［21］、基于Transformer和NeRF架构的结构感知稀疏视图3D重建的双域方法SAX-NeRF^［29］、基于预训练的2D扩散模型和SIRT先验解决3D逆问题的双域方法TOSM^［27］。

1.5.4 评价指标

为定量评价图像质量，采用峰值信噪比（PSNR）、结构相似度（SSIM）、均方根误差（RMSE）指标对稀疏角度锥束CT重建结果进行分析，具体定义如下：

PSNR=\mathrm{10}lo{g}_{\mathrm{10}}\left(\frac{{\left({\mathrm{2}}^n-\mathrm{1}\right)}^{\mathrm{2}}}{MSE}\right)，

（13）

SSIM=\frac{\left(\mathrm{2}{\mu }_x{\mu }_y+c_{\mathrm{1}}\right)\left({\sigma }_{xy}+c_{\mathrm{2}}\right)}{\left({\mu }_x^{\mathrm{2}}+{\mu }_y^{\mathrm{2}}+c_{\mathrm{1}}\right)\left({\sigma }_x^{\mathrm{2}}+{\sigma }_y^{\mathrm{2}}+c_{\mathrm{2}}\right)}，

（14）

RMSE=MS{E}^{\mathrm{1}/\mathrm{2}}（15）

式中n表示图像存储的位数，MSE为均方误差；SSIM指标实际由亮度、对比度、结构比较3部分构成，这里仅展示最终公式（14），式中{\mu }_x、{\mu }_y表示输入图像x、y的均值，{\sigma }_x^{\mathrm{2}}、{\sigma }_y^{\mathrm{2}}为对应的方差，{\sigma }_{xy}是协方差，c_{\mathrm{1}}、c_{\mathrm{2}}为常数，其作用是避免分母为零。

2 结果

为评估BBC-Recon方法的性能，本文与FDK、ADS-POCS、FBPConvNet、Dual-CNN、SAX-NeRF、TOSM方法进行比较，采用4倍和8倍两个稀疏强度，选用两个数据集，包含胸腔、腹腔和口腔3个部位，以及横断、冠状、矢状3个切面，从定量和定性上进行结果对比。最后，本文对各方法的测试用时进行了统计，用以评估方法性能；同时，还开展了3组消融实验，验证所提方法设计的各关键模块的重要性和有效性。值得注意的是，在定量分析中，我们对3个切面的PSNR、SSIM、RMSE指标均进行计算，为排除锥角伪影影响，计算舍弃了前后各10层图像，仅对中间44层图像计算，然后取均值作为最终的评价指标值。

2.1 胸腔部位不同稀疏强度重建效果

图6是对胸腔部位进行锥束CT稀疏4倍角度重建的结果，其中，横断面结果下方为冠状面结果，右侧为矢状面结果。不难看出，FDK算法直接重建存在明显的条纹伪影和噪声（图6B）；ASD-POCS方法结果伪影残留、结构边缘模糊明显（图6C）；FBPConvNet方法也有伪影残留（图6D放大区箭头处）；SAX-NeRF方法结构结构模糊、细节丢失（图6E）；Dual-CNN方法结果偏平滑、引入额外结构（图6F放大区绿色箭头处）；TOSM方法结果噪点较多、有伪影残留（图6G黄色放大区黄色箭头处），同时冠状、矢状面有条纹结构伪影（图6G冠状矢状面箭头处）；BBC-Recon方法结果中，没有伪影残留、没有条纹结构伪影和额外结构的引入（图6H）。表1是定量结果，加粗为最优，下划线为次优，可以看出，BBC-Recon方法结果均为最优，FBPConvNet方法次之。

图7是对胸腔部位进行锥束CT稀疏8倍角度重建的结果，其中，横断面结果下方为冠状面结果，右侧为矢状面结果。可以看到，8倍稀疏强度下，FDK直接重建存在严重的条纹伪影（图7B）；ASD-POCS方法结果的结构严重模糊、伪影残留多（图7C）；FBPConvNet方法结果残留条纹伪影（图7D绿色箭头处和黄色放大区）；SAX-NeRF方法的结果结构严重模糊、细节丢失（图7E）；Dual-CNN方法结果较为平滑、结构异常（图7F红色放大区箭头处）；TOSM方法结果结构清晰，但是噪点多、有伪影残留（图7G黄色放大区箭头处）、冠状和矢状面有条纹结构伪影（图7G冠状矢状面箭头处）；BBC-Recon方法结果有轻微平滑，但没有伪影残留和异常结构（图7H）。表2是定量结果，加粗为最优，下划线为次优，可以看到，BBC-Recon方法在横断面的RMSE指标上表现次优，冠状、矢状面结果指标均为最优，TOSM方法次之。

2. 2 腹腔部位不同稀疏强度重建效果

图8是对腹腔部位进行锥束CT稀疏4倍角度重建结果，其中，横断面结果下方为冠状面结果，右侧为矢状面结果。FDK直接重建存在大量条纹伪影（图8B）；ASD-POCS方法结构模糊、锥体边缘丢失（图8C红色放大区箭头处）；FBPConvNet方法结果存在条纹伪影残留（图8D放大区箭头处）；SAX-NeRF方法结果的结构丢失、模糊（图8E）；Dual-CNN方法结果有伪影残留、锥体边缘轻微畸形（图8F放大区箭头处）；TOSM方法结果中锥体边缘丢失（图8G红色放大区箭头处）、冠状和矢状面中存在明显条纹结构伪影（图8G冠状矢状面箭头处）；BBC-Recon方法结果结构恢复清晰，没有伪影残留、结构丢失和条纹结构伪影出现（图8H）。表3是定量结果，加粗为最优，下划线为次优，可以看出，BBC-Recon方法的所有指标均为最优，FBPConvNet方法次之。

图9是对腹腔部位进行锥束CT稀疏8倍角度重建结果，其中，横断面结果下方为冠状面结果，右侧为矢状面结果。在更高的稀疏强度下，FDK直接重建结果存在严重条纹伪影（图9B）；FBPConvNet方法结果在肝脏区域存在结构丢失、在肝脏和肾脏的交界存在伪影残留（图9D绿色箭头处）；SAX-NeRF方法结果结构丢失、图像模糊（图9E）；Dual-CNN方法结果在肝脏区域存在结构丢失，有条纹伪影残留，同时锥体边缘结构畸变（图9F红色放大区红色箭头处）；TOSM方法结果锥体边缘丢失、锥孔区域存在二次伪影（图9G红色放大区），同时在冠状和矢状面存在条纹结构伪影（图9G冠状矢状面箭头处）；BBC-Recon方法结果有轻微平滑，没有结构畸变、边缘丢失和伪影残留问题（图9H）。表4是定量结果，加粗为最优，下划线为次优，在横断面指标中，BBC-Recon方法表现次优，在冠状、矢状面的峰值信噪比、结构相似度指标上最优，TOSM方法次之。

2.3 口腔数据集不同稀疏强度重建效果

图10是对口腔数据集进行锥束CT稀疏4倍角度重建结果，其中，横断面结果下方为冠状面结果，右侧为矢状面结果。FDK直接重建结果存在条纹伪影（图10B）；ASD-POCS方法结果噪点多、结构模糊（图10C）；FBPConvNet方法结果存在伪影残留（图10D放大区绿色箭头）；SAX-NeRF方法结果平滑、结构恢复不完整（图10E红色放大区红色箭头处）；Dual-CNN方法结构模糊、结构丢失（图10F红色放大区红色箭头处）；TOSM方法结果存在二次伪影（图10G黄色箭头处）；BBC-Recon方法结果没有伪影残留、无二次伪影出现，结构恢复较为完整（图10H）。表5是定量结果，加粗为最优，下划线为次优，除冠状面的结构相似度指标外，其余指标中，BBC-Recon方法结果均为最优，TOSM方法次之。

图11是对口腔数据集进行锥束CT稀疏8倍角度重建结果，其中，横断面结果下方为冠状面结果，右侧为矢状面结果。在该稀疏强度下，FDK直接重建结果存在明显条纹伪影（图11B）；ASD-POCS方法结果中存在伪影残留，同时结构模糊、结构丢失（图11C）；FBPConvNet方法存在异常结构（图11D红色箭头处）、牙髓结构丢失（图D绿色箭头处）；SAX-NeRF方法结果存在结构丢失、结构模糊（图11E放大区域）；Dual-CNN方法结果存在明显条纹伪影（图11F红色箭头处）；TOSM方法结果有二次伪影产生（图11G黄色箭头处）；BBC-Recon方法结构无伪影残留、无二次伪影产生，同时保留了牙髓结构（图11H）。表8是定量结果，加粗为最优，下划线为次优，BBC-Recon方法在除冠状面的结构相似度指标上表现为次优外，其余指标中，BBC-Recon方法结果均为最优，TOSM方法次之。

2.4 不同方法测试用时对比

除TOSM方法外，其余方法均进行10组测试试验，每组输出64层图像，然后对各组用时加权平均，除以图像数64，得到输出一层图像所需的时间；对于TOSM方法，因其测试含大量迭代，1个case推理耗时约19h，所以仅开展5组试验（表7）。为充分客观地展示方法性能，我们完成了在4倍、8倍两个稀疏强度下的测试，所有测试均在一张NVIDIA RTX 4090显卡上完成、batch size均为1。Dual-CNN方法因其网络尺寸小，测试用时最短，一张图像的恢复仅需不到50 ms，FBPConvNet方法次之；所提方法BBC-Recon输出一张图像用时140 ms左右，略逊于前两者，但明显优于其他3个对比方法，尤其值得注意的是，虽然TOSM方法在图像分辨率保持上获得了较好的效果，这实际上得益于测试阶段的多次迭代，这从表7的测试用时可以看出，TOSM输出一张图像需要长达18 min。

2.5 主要模块消融实验

本研究实施了3组消融实验，以稀疏8倍角度的胸腔数据重建为例。基线实验仅用单向和单尺度关联体估计策略，记为E1；第2组实验采用双向和单尺度关联体估计策略，记为E2；第3组实验采用双向和多尺度关联体估计策略，记为E3。每组实验均训练300 epochs，测试时所用的模型均选用训练过程中的最优模型，对3组实验恢复的投影数据均采用FDK重建，并进行量化指标的计算，加粗表示最优，下划线为次优（表8）。基线实验E1在不使用双向和多尺度模块的情形下，对投影间的相关性和冗余信息的挖掘不足，未能充分对缺失投影进行恢复，进而导致在PSNR、SSIM、RMSE 3个指标上获得最差结果；在第2组实验E2中，通过增加双向关联体估计功能，可以有效提高投影补全的准确性，进而提高图像质量，各量化指标均有提升；而在第3组实验E3中，通过在E2的基础上增加多尺度关联体估计策略，发现图像质量得到持续提升，量化指标均获得最优结果，进而证明了多尺度关联体估计策略在提高投影补全精度任务中所发挥的积极作用，与此同时，对比E2和E3发现，E3指标的提升程度大于E2，例如，在横断面的PSNR指标中：E2较E1提升了0.0782，而E3较E2提升了0.1366，其他指标也均有此规律。

3 讨论

针对稀疏角度锥束CT成像面临的不适定反问题求解问题，本文提出了一种基于流场引导投影补全的双域稀疏角度锥束CT成像算法（BBC-Recon），用于在保留结构细节的基础上去除图像的条纹伪影。本文方法从提高投影补全的精度的角度出发，使用流场估计思想^［31-33］，通过设计的双向流场引导投影补全的投影恢复模块，有效实现对两帧投影之间双向的变化进行计算，来指导缺失投影的合成，其中，双向关联体^［33］的设计，更好地理解了像素移动的复杂变化，而多尺度关联体^［31］模块的加入，有效捕捉了从微观到宏观不等的特征，实现了对不同层次的相关性信息的挖掘，双向和多尺度二者结合，极大程度地提高了对投影之间的相关性和冗余信息的挖掘，实现了对补全投影的连续性和准确性的提升，弥补了传统的线性、方向插值或者字典学习等方法对投影补全的精度低的不足。另外，该方法通过引入图像域的细化模块，实现了投影＋图像的双域结合训练，在高精度投影补全的基础上，进行图像的细化，不仅有效实现了对图像条纹伪影的去除，而且获得了传统图像域方法所难以实现的图像细节信息的保留。值得注意的是，本文方法通过构建的基于流场估计的投影恢复模块，实现了对稀疏投影的3D恢复，即一次推理便可获得完整的锥束投影，进而推动对稀疏角度锥束CT成像不适定反问题的求解，这是现有的研究平行束、扇形束等2D不适定反问题求解的双域方法所难以达到的。

实验结果表明，ASD-POCS^［6］作为TV约束的迭代重建方法，计算复杂度较高，而且受限于投影的完备性，其性能在较高的稀疏强度下快速下降；FBPConvNet^［14］直接学习伪影-真值图像映射的方式过于粗糙，伪影图像掩盖了大量有效信息会增加网络学习的难度，难以避免伪影的残留和结构的模糊；双域方法SAX-NeRF^［29］通过学习坐标-衰减值的映射进行3D体素的渲染，实际上其图像的重建直接来源于渲染的体素，图像域损失的梯度回传并未利用投影之间的信息，难以对体素的渲染实现有效更新；Dual-CNN^［21］利用多个独立训练好的网络逐级对投影进行恢复，容易叠加误差，这也是现有的多数双域方法共有的不足；生成式的方法TOSM^［27］，通过在迭代去噪中实现伪影的去除，但其代价是生成的不确定性和长的推理时间；所提方法BBC-Recon，得益于对投影之间的相关性信息和冗余信息进行充分挖掘，联合图像域进行双域训练，对比上述方法，该方法结果不仅没有伪影残留，而且对结构细节的恢复较为完整，同时量化指标也是最优；尤其值得注意的是，得益于该方法是直接基于锥束投影进行的恢复，不需要对其进行平行束或者扇形束的转换，所以其对图像的层间一致性实现了有效的保持，优于其他的投影恢复方法。

但是，本文方法也存在一定的不足，即投影之间的变化并非是自然图像中像素单纯的流动，即使是利用了双向和多尺度的优势，也难以完整无误的挖掘投影之间的相关性和冗余信息，而单独的图像域细化只能从图像上一定程度上缓解这种不足，并不能解决该问题，这可能也是导致图像趋于平滑的原因。理论上，方法应该需要从投影域出发去克服该问题，即利用图像域所具有的先验信息，将图像域与投影域联合训练，实现梯度从图像到投影的回传，完成对投影域恢复误差的校正，这也将是我们下一步开展的工作。

综上，本文提出的BBC-Recon方法，首次将流场估计的思想应用到锥束CT稀疏角度成像研究中，为锥束CT成像研究中的不适定问题求解提供了新思路；该方法利用双向和多尺度关联体设计的优势，实现了对投影间相关性和冗余信息的最大程度挖掘，得以对双向流场和中间帧估计提供精准指导，有效实现了缺失投影的高质量恢复，弥补了现有投影补全方法难以对投影相关性和冗余信息充分利用、进而引入明显插值误差导致二次伪影的不足，与此同时，联合图像域细化，最终有效抑制了条纹伪影、恢复了结构细节信息、获得了高的层间一致性，实现了高质量的锥束CT稀疏角度重建。

参考文献

[1] Scarfe WC, Farman AG, Sukovic P. Clinical applications of cone-beam computed tomography in dental practice[J]. J Can Dent Assoc, 2006, 72(1): 75-80.

[2] Brenner DJ, Hall EJ. Computed tomography: an increasing source of radiation exposure[J]. N Engl J Med, 2007, 357(22): 2277-84.

[3] Berrington de González A, Mahesh M, Kim KP, et al. Projected cancer risks from computed tomographic scans performed in the United States in 2007[J]. Arch Intern Med, 2009, 169(22): 2071-7.

[4] Feldkamp LA, Davis LC, Kress JW. Practical cone-beam algorithm[J]. J Opt Soc Am A, 1984, 1(6): 612.

[5] Andersen AH, Kak AC. Simultaneous algebraic reconstruction technique (SART): a superior implementation of the art algorithm[J]. Ultrason Imaging, 1984, 6(1): 81-94.

[6] Sidky EY, Pan XC. Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization[J]. Phys Med Biol, 2008, 53(17): 4777-807.

[7] Li S, Cao Q, Chen Y, et al. Dictionary learning based sinogram inpainting for CT sparse reconstruction[J]. Optik, 2014, 125(12): 2862-7.

[8] Bertram M, Rose G, Schafer D, et al. Directional interpolation of sparsely sampled cone-beam CT sinogram data[C]//2004 2nd IEEE International Symposium on Biomedical Imaging: Nano to Macro. Arlington, VA, USA. IEEE, 2004: 928-31.

[9] Zhang H, Sonke JJ. Directional sinogram interpolation for sparse angular acquisition in cone-beam computed tomography[J]. J Xray Sci Technol, 2013, 21(4): 481-96.

[10] Muhlich M, Aach T. Analysis of multiple orientations[J]. IEEE Trans Image Process, 2009, 18(7): 1424-37.

[11] Li ZB, Yu LF, Trzasko JD, et al. Adaptive nonlocal means filtering based on local noise level for CT denoising[J]. Med Phys, 2014, 41(1): 011908.

[12] Litjens G, Kooi T, Bejnordi BE, et al. A survey on deep learning in medical image analysis[J]. Med Image Anal, 2017, 42: 60-88.

[13] Razzak MI, Naz S, Zaib A. Deep learning for medical image processing: overview, challenges and future[EB/OL]. 2017: arXiv: 1704.06825.

[14] JinKyong Hwan, McCann MT, Froustey E, et al. Deep convolutional neural network for inverse problems in imaging[J]. IEEE Trans Image Process, 2017, 26(9): 4509-22.

[15] Ronneberger O, Fischer P, Brox T. U-net: Convolutional networks for biomedical image segmentation[C]//Medical image computing and computer-assisted intervention-MICCAI 2015: 18th international conference, Munich, Germany, October 5-9, 2015, proceedings, part III 18. Springer International Publishing, 2015: 234-41.

[16] Chen H, Zhang Y, Kalra MK, et al. Low-dose CT with a residual encoder-decoder convolutional neural network[J]. IEEE Trans Med Imaging, 2017, 36(12): 2524-35.

[17] Zhong WH, Li TL, Hou SK, et al. Unsupervised disentanglement strategy for mitigating artifact in photoacoustic tomography under extremely sparse view[J]. Photoacoustics, 2024, 38: 100613.

[18] Wu WW, Hu DL, Niu C, et al. DRONE: dual-domain residual-based optimization NEtwork for sparse-view CT reconstruction[J]. IEEE Trans Med Imaging, 2021, 40(11): 3002-14.

[19] Wu WW, Guo XD, Chen Y, et al. Deep embedding-attention-refinement for sparse-view CT reconstruction[J]. IEEE Trans Instrum Meas, 2023, 72: 4501111.

[20] Lahiri A, Maliakal G, Klasky ML, et al. Sparse-view cone beam CT reconstruction using data-consistent supervised and adversarial learning from scarce training data[J]. IEEE Trans Comput Imag, 2023, 9: 13-28.

[21] Chao LY, Wang ZW, Zhang HB, et al. Sparse-view cone beam CT reconstruction using dual CNNs in projection domain and image domain[J]. Neurocomputing, 2022, 493: 536-47.

[22] Kim S, Kim B, Lee J, et al. Sparsier2Sparse: self-supervised convolutional neural network-based streak artifacts reduction in sparse-view CT images[J]. Med Phys, 2023, 50(12): 7731-47.

[23] Pan JY, Zhang HY, Wu WF, et al. Multi-domain integrative Swin transformer network for sparse-view tomographic reconstruction[J]. Patterns, 2022, 3(6): 100498.

[24] Genzel M, Gühring I, MacDonald J, et al. Near-exact recovery for tomographic inverse problems via deep learning[J]. Inter Confer Machine Learn, 2022，2022: 7368-81.

[25] Mildenhall B, Srinivasan PP, Tancik M, et al. NeRF: representing scenes as neural radiance fields for view synthesis[C]//European Conference on Computer Vision. Cham: Springer, 2020: 405-421.

[26] Chung H, Ryu D, McCann MT, et al. Solving 3D inverse problems using pre-trained 2D diffusion models[C]//2023 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR). Vancouver, BC, Canada. IEEE, 2023: 22542-51.

[27] Li ZR, Wang YY, Zhang JJ, et al. Two-and-a-half order score-based model for solving 3D ill-posed inverse problems[J]. Comput Biol Med, 2024, 168: 107819.

[28] Liu ZT, Fang Y, Li CJ, et al. Geometry-aware attenuation learning for sparse-view CBCT reconstruction[EB/OL]. 2023: arXiv: 2303. 14739.

[29] Cai YH, Wang JH, Yuille A, et al. Structure-aware sparse-view X-ray 3D reconstruction[C]//2024 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR). Seattle, WA, USA. IEEE, 2024: 11174-83.

[30] Liu ZT, Zhao HX, Qin WH, et al. 3D vessel reconstruction from sparse-view dynamic DSA images via vessel probability guided attenuation learning[EB/OL]. 2024: arXiv: 2405.10705.

[31] Teed Z, Deng J. RAFT: recurrent all-pairs field transforms for optical flow[C]//European Conference on Computer Vision. Cham: Springer, 2020: 402-419.

[32] Kong LT, Jiang BY, Luo DH, et al. IFRNet: intermediate feature refine network for efficient frame interpolation[C]//2022 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR). New Orleans, LA, USA. IEEE, 2022: 1959-68.

[33] Li Z, Zhu ZL, Han LH, et al. AMT: all-pairs multi-field transforms for efficient frame interpolation[C]//2023 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR). Vancouver, BC, Canada. IEEE, 2023: 9801-10.

[34] Charbonnier P, Blanc-Feraud L, Aubert G, et al. Two deterministic half-quadratic regularization algorithms for computed imaging[C]//Proceedings of 1st International Conference on Image Processing. Austin, TX, USA. IEEE, 1994: 168-72.

[35] Moen TR, Chen BY, Holmes DR 3rd, et al. Low-dose CT image and projection dataset[J]. Med Phys, 2021, 48(2): 902-11.

[36] Huang YY, Liu WJ, Yao CQ, et al. A multimodal dental dataset facilitating machine learning research and clinic services[J]. Sci Data, 2024, 11(1): 1291.

基金资助

RIGHTS & PERMISSIONS

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