Unsupervised CT Metal Artifact Reduction by Plugging Diffusion Priors in Dual Domains

通过在双域中插入扩散先验来减少无监督 CT 金属伪影

论文链接：https://arxiv.org/abs/2308.16742

项目链接：https://github.com/DeepXuan/DuDoDp-MAR

Abstract

在计算机断层扫描(CT)过程中，患者体内的金属植入物通常会导致重建图像中的破坏性伪影，从而阻碍准确诊断。许多基于监督深度学习的方法被提出用于金属伪影还原(MAR)。然而，这些方法依赖于大量的模拟数据训练，因为在临床环境中获得配对的金属伪影CT和干净CT数据是具有挑战性的。这一限制可能导致在临床实践中应用时性能下降。现有有效的无监督MAR方法，无论是否基于学习，通常都是在图像域或弦图域进行单域处理。在本文中，我们引入了一种基于扩散模型的无监督MAR方法，扩散模型是一种具有强大表示数据分布能力的生成模型。具体来说，我们首先用不含金属伪影的CT图像训练扩散模型。随后，我们在弦图域和图像域迭代利用预训练扩散模型中嵌入的先验，旨在恢复由金属伪影引起的退化部分。双域处理使我们的方法优于现有的无监督MAR方法，包括另一种基于扩散模型的MAR方法，该方法已经在合成数据集上进行了定性和定量验证。此外，我们的方法在临床数据集的监督和非监督方法中都显示出优越的视觉效果。

I. INTRODUCTION

计算机断层扫描在现代医学的诊断和治疗中起着至关重要的作用。然而，患者体内的金属植入物会在重建的 CT 图像中造成严重的伪影 [1]、[2]，这可能会严重影响后续的诊断程序。金属植入物的存在会使部分弦图值出现异常。由这种异常弦图重建的CT图像在非金属区域往往会出现破坏性伪影[3]。因此，人们设计了许多算法来减轻重建CT图像中被金属污染的金属伪影。

最常用的金属伪影减少 (MAR) 传统方法是基于弦图修复，该方法在执行重建过程之前去除并修复弦图中受金属影响的区域。Kalender 等人提出使用线性插值 (LI) 来填充空心弦图 [3]。 Meyer 等人提出了归一化金属伪影减少 (NMAR)，用先前图像的前向投影替换弦图中受金属影响的区域 [4]。 尽管 LI 和 NMAR 都是简单而有效的 MAR 方法，但由于完整弦图边缘可能存在不连续性，重建图像仍可能包含伪影。

随着过去十年深度学习技术的快速发展，许多基于深度学习的MAR方法被提出，与传统方法相比有了显著的改进。大多数基于深度学习的MAR方法都是有监督的，依靠配对的金属伪影CT和干净CT数据来训练神经网络[5]-[9]。与传统方法一样，基于深度学习的MAR方法可以在弦图域、图像域或两者兼而有之。Park等人提出使用卷积神经网络(convolutional neural network, CNN)修复由于线束硬化引起的不一致弦图[2]。Zhang等人提出了CNNMAR，由CNN产生先验图像，然后用先验图像的正投影完成金属影响弦图[10]。此外，后处理深度神经网络试图直接从金属伪像CT图像中学习映射来清洗CT图像。Huang等利用深度残差学习对颈椎CT图像进行MAR[11]。Wang等人使用生成式对抗网络(GAN)[12]对耳部CT图像进行MAR[13]。利用弦图域和图像域的信息，提出了一些双域网络来同时恢复弦图一致性和增强CT图像。Lin等人提出了端到端可训练双域网络(DuDoNet)[5]。Zhou等人提出了一种双域数据一致循环网络(DuDoRNet)[8]。Wang等人提出了一种可解释的双域网络（InDuDoNet）[6]及其改进版本InDuDoNet+[9]。双域网络对弦图数据和图像数据进行连续或循环处理，与单域网络相比，具有更好的性能。

然而，以监督的方式训练MAR神经网络需要大量配对的CT数据，有或没有金属伪影，这在临床中实际上是无法实现的。因此，监督方法通常必须专门针对模拟数据进行训练，这可能导致由于领域差距而导致临床应用中的次优结果。因此，人们提出了一些无监督或弱监督的方法来解决缺乏成对训练数据的问题。Liao等人引入了第一种无监督学习的MAR方法，名为artifact disentanglement network (ADN)，实现了金属影响图像与干净图像之间的翻译[14]。Du等人提出了一种新的基于无监督域自适应(unsupervised domain adaptation, UDA)的MAR方法，称为UDAMAR来解决域间隙问题[15]。此外，一些一般的无监督图像到图像的翻译方法也可以采用forMAR，如CycleGAN [16]， MUNIT[17]。这些流行的无监督MAR方法在图像域内工作。然而，由于金属工件在图像域的严重和非结构化的性质，这些方法的有效性，特别是在视觉质量方面，往往是有限的。Song等人提出了一种新的MAR方法，他们将任务表述为线性逆问题，并使用预训练分数模型[18]作为正则化器[19]。具体来说，他们使用来自分数模型的先验图像迭代地对弦图域中的空心区域进行补全。该方法很有见地，但我们认为它没有充分利用分数模型的强大功能，因为分数模型提供的先验和已知部分的组合只发生在弦图域中。

为了解决现有无监督MAR方法的局限性，我们提出了一种新的MAR方法，该方法在弦图域和图像域使用预训练扩散模型[20]中的嵌入先验。具体来说，我们首先用大量不含金属伪影的CT图像训练扩散模型进行准备。训练良好的算法不仅能够生成视觉逼真的干净CT图像，而且能够用解析概率密度函数表示其分布。随后，我们从弦图中提取金属影响部分，并将MAR模型作为部分测量缺失的重建问题。这个问题是迭代解决的。在扩散生成过程的每个时间步，我们首先将扩散的输出与已知弦图在弦图域中融合。由于扩散模型提供的先验弦图与已知弦图之间的不连续，重构图像可能会在图像域中引入新的伪影。因此，我们在图像域重新引入扩散先验来缓解这一问题。此外，为了保证先验和似然的平衡，我们在图像域精心设计了加权掩膜进行融合，并验证了其有效性。

我们将这种具有扩散先验的双域方法命名为DuDoDp。我们的贡献可以概括如下：

1. 我们提出了一种在双域框架内结合扩散先验的MAR方法，以充分利用来自弦图域和图像域的信息。
2. 提出在每个时间步融合三种类型的图像，分别是扩散模型的先验图像、补全弦图的重构图像和原始金属伪影图像。这些图像的组合校正了补全中单一弦图所造成的伪影。
3. 设计动态权值掩膜来控制扩散先验的加入。根据扩散的迭代生成特性，我们随时间步长动态调整掩模以获得增强的结果。
4. 在合成数据集和临床数据集上验证了该方法，证明了其临床应用潜力。

II. METHOD

A. 初步

1) 扩散模型：扩散模型(Diffusion Model，即去噪扩散概率模型的简称)是一种深度生成模型[18]，[20]。给定一组数据$x_0\sim q(x_0)$，扩散模型定义了一个固定的前向过程，该过程逐渐向原始数据中添加高斯噪声。然后尝试学习一个反向过程来反转正向过程，实现从随机噪声中生成图像。一般将前向过程$q(x_{1:T}|x_0)$定义为马尔可夫链，即每个时间步的状态只与前一个时间步相关。具体来说，有：
$$\begin{array}{rl}q\left(x_{1:T}\mid x_0\right)& :=\prod _{t=1}^{T}q\left(x_t\mid x_{t-1}\right),q\left(x_t\mid x_{t-1}\right)\\ & :=\mathcal{N}\left(x_t;\sqrt{1-{\beta }_t}{x}_{t-1},{\beta }_{t}I\right),\end{array}\text{\hspace{1em}(1)}$$
其中$x_t$为时间步长t的数据，$\mathcal{N}(\cdot ;\mu ,\Sigma )$表示高斯分布，平均值$\mu$，协方差矩阵$\Sigma$， { β 1 , … , β T } \{β_1,…,β_T\} {β1​,…,βT​}为固定方差表。根据公式(1)，可得以下两个条件分布：
$$\begin{array}{r}q\left(x_t\mid x_0\right)=\mathcal{N}\left(x_t;\sqrt{{\alpha }_t}{x}_0,\left(1-{\bar{\alpha }}_t\right)I\right),\\ q\left(x_{t-1}\mid x_t,x_0\right)=\mathcal{N}\left(x_{t-1};{\tilde{\mu }}_t\left(x_t,x_0\right),{\tilde{\beta }}_{t}I\right),\end{array}\text{\hspace{1em}(2)(3)}$$
其中，
$$\begin{array}{rl}& {\bar{\alpha }}_t:=\prod _{i=1}^t(1-{\beta }_i),\\ & {\tilde{\mu }}_t\left(x_t,x_0\right):=\frac{\sqrt{{\bar{\alpha }}_{t-1}}{\beta }_t}{1-{\bar{\alpha }}_t}{x}_0+\frac{\sqrt{{\alpha }_t}\left(1-{\bar{\alpha }}_{t-1}\right)}{1-{\bar{\alpha }}_t}{x}_t,\\ & {\tilde{\beta }}_t:=\frac{(1-{\bar{\alpha }}_{t-1})}{1-{\bar{\alpha }}_t}{\beta }_t.\end{array}\text{\hspace{1em}(4)(5)(6)}$$
然后将反向过程定义为具有可学习参数$\theta$的马尔可夫链，
$$p_{\theta }\left(x_{0:T}\right):=p\left(x_T\right)\prod _{t=1}^T{p}_{\theta }\left(x_{t-1}\mid x_t\right),\text{\hspace{1em}(7)}$$
其中，
$$\begin{array}{rl}& p\left(x_T\right)=\mathcal{N}\left(x_T;0,I\right),\\ & p_{\theta }\left(x_{t-1}\mid x_t\right):=\mathcal{N}\left(x_{t-1};{\mu }_{\theta }\left(x_t,t\right),{\sigma }_t^{2}I\right).\end{array}\text{\hspace{1em}(8)(9)}$$
其中，${\mu }_{\theta }(x_t,t)$由可训练神经网络预测，${\sigma }_t^2$定义为$\frac{\bar{1}-{\bar{\alpha }}_{t-1}}{1-{\bar{\alpha }}_t}{\beta }_t$。训练损失是通过最小化负对数似然$-\log p_{\theta }(x_0)$的变分界得到的[20]。

2) $ϵ$-Prediction和$x_0$-Prediction扩散模型：通常，扩散模型不能直接预测公式(9)中每个时间步长的平均值${\mu }_{\theta }(x_t,t)$。或者，更常见的是预测添加的噪声${ϵ}_t$或从$x_t$中预测初始图像$x_0$，其中$x_t=\sqrt{{\bar{\alpha }}_t}{x}_0+\sqrt{1-{\bar{\alpha }}_t}{ϵ}_t$。预测噪声${ϵ}_{\theta }(x_t,t)$与${\mu }_{\theta }(x_t,t)$的关系为[20]，
$$\begin{array}{c}{\mu }_{\theta }\left({\mathbf{x}}_t,t\right)={\tilde{\mu }}_t\left({\mathbf{x}}_t,\frac{1}{\sqrt{{\bar{\alpha }}_t}}\left(x_t-\sqrt{1-{\bar{\alpha }}_t}{ϵ}_{\theta }\left(x_t\right)\right)\right)\\ =\frac{1}{\sqrt{{\alpha }_t}}\left(x_t-\frac{{\beta }_t}{\sqrt{1-{\bar{\alpha }}_t}}{ϵ}_{\theta }\left(x_t,t\right)\right),\end{array}\text{\hspace{1em}(10)}$$
其中的${\tilde{\mu }}_t$来自公式(5)。从$x_t$预测${ϵ}_{\theta }(x_t,t)$使得扩散模型的训练更加简洁，这等于从噪声图像$x_t$和噪声水平$t$中预测附加噪声${ϵ}_t$[20]。

当然，我们也可以从$x_t$中预测出潜在的干净图像$x_0$。设$x_0$的预测值为$f_{\theta }(x_t,t)$
$$f_{\theta }\left(x_t,t\right)=(x_t-\sqrt{1-{\bar{\alpha }}_t}{ϵ}_{\theta }(x_t,t))/\sqrt{{\bar{\alpha }}_t}.\text{\hspace{1em}(11)}$$
在此设置下，公式(9)中扩散模型的推理过程可重写为[21]：
$$p_{\theta }\left(x_{t-1}\mid x_t\right)=q\left(x_{t-1}\mid x_t,f_{\theta }\left(x_t,t\right)\right),\text{\hspace{1em}(12)}$$
式中$q\left(x_{t-1}\mid x_t,f_{\theta }\left(x_t,t\right)\right)$根据公式(3)计算。

无论预测${\mu }_{\theta }$、${ϵ}_{\theta }$还是$f_{\theta }$，扩散模型的原理及其训练损失的含义都是不变的。唯一的区别在于对神经网络输出的解释。

B. 双域扩散先验的MAR

在本节中，我们介绍了DuDoDp，一种具有双域扩散先验的MAR方法。首先，在II-B.1节通过迭代求解多个优化问题，将扩散先验引入到MAR问题中。这种合并相当于利用扩散模型提供的先验图像在弦图域内进行图像补全。为了减轻因所补全弦图不连续而产生的伪影，我们建议在图像域中进一步整合扩散先验，如第II-B.2节所述。在II-B.3节中设计了图像域融合的权重掩模，考虑到扩散的迭代性，在不同的时间步长采用动态掩模。整体算法流程如图1所示。

1) MAR的扩散先验引入：MAR问题可以被视为从受金属影响的弦图$s_0$重建干净图像$x_0$的任务。 直接应用分析算法（例如滤波反投影 (FBP)）进行重建将得到带有金属伪影的图像 $y_0$，其中 $y=FBP(s_0)$。 通常，我们知道弦图的哪些部分受到金属的影响，用布尔掩膜 $M_s$ 来表示。在这个掩膜中，未受影响的部分标记为 0，而受影响的部分标记为 1。 $s_0$ 和 $x_0$ 之间的关系是，
$$(1-{\mathcal{M}}_s)\odot s_0=(1-{\mathcal{M}}_s)\odot \mathrm{FP}(x_0)\text{\hspace{1em}(13)}$$
因此，MAR问题可表示为：
$$\underset{x_0}{\min }\Vert (1-{\mathcal{M}}_s)\odot (\mathrm{FP}(x_0)-s_0){\Vert }^2+\lambda \mathcal{R}(x_0),\text{\hspace{1em}(14)}$$
式中，FP(·)为CT系统的正投影函数，$\odot$为Hadamard哈达玛积，$\mathcal{R}(\cdot )$为干净图像$x_0$的正则化器，λ为系数。

我们尝试从预训练的扩散模型中得到正则器$\mathcal{R}$。然而，扩散模型并不直接表示数据的概率密度函数，而是迭代地得到它。具体来说，在已知$x_t$的情况下，根据公式(12)，数据在时间步长$t-1$处的概率密度函数为$q(x_{t-1}|x_t,f_{\theta }(x_t,t))$。

为了利用时间步长$t-1$的先验，我们可以根据公式(3)迭代求解与$x_0$有明确关系的$x_{t-1}$，
$$x_{t-1}={\tilde{\mu }}_t\left(x_t,x_0\right)+\sqrt{{\tilde{\beta }}_t}{z}_{t-1},\text{\hspace{1em}(15)}$$
其中${\tilde{\mu }}_t$来自公式(5)，$z_{t-1}\sim \mathcal{N}(0,I)$。同样地，我们可以定义一种方法来得到$s_{t-1}$，
$$s_{t-1}={\tilde{\mu }}_{t-1}\left(\mathrm{FP}(x_t),s_0\right)+\mathrm{FP}({\sqrt{\beta }}_t{z}_{t-1}).\text{\hspace{1em}(16)}$$
这样，我们有，
$$(1-{\mathcal{M}}_s)\odot s_{t-1}=(1-{\mathcal{M}}_s)\odot \mathrm{FP}(x_{t-1})\text{\hspace{1em}(17)}$$
根据公式(17)，在每个时间步长$t-1$，我们可以迭代求解这样的优化问题，
$$\begin{array}{r}\underset{x_{t-1}}{\min }\Vert (1-{\mathcal{M}}_s)\odot (\mathrm{FP}(x_{t-1})-s_{t-1}){\Vert }^2+\lambda q\left(x_{t-1}\mid x_t,f_{\theta }\left(x_t,t\right)\right).\end{array}\text{\hspace{1em}(18)}$$
在公式(18)中，我们通过迭代方法成功地将扩散模型提供的先验纳入到MAR问题的求解中。因为$x_{t-1}$和$s_{t-1}$都是以同样的方式由$x_0$和$s_0$变换而来(见公式15和公式16)，通过消去某些系数，公式(18)变成，
$$\begin{array}{r}\underset{{\tilde{x}}_0^{t-1}}{\min }\Vert (1-{\mathcal{M}}_s)\odot (\mathrm{FP}({\tilde{x}}_0^{t-1})-s_0){\Vert }^2+{\lambda }^{{}^{\prime }}\Vert {\tilde{x}}_0^{t-1}-f_{\theta }(x_t,t){\Vert }^2,\end{array}\text{\hspace{1em}(19)}$$
其中${\tilde{x}}_0^{t-1}$表示时间步 $t-1$ 处的预测初始图像，${\lambda }^{{}^{\prime }}$是不用于求解公式 (19) 的系数。对于公式（19）中的凸优化问题，有许多迭代算法可用于解决该问题。然而，为了缓解与迭代求解相关的计算时间延长的问题，我们直接使用 $f_{\theta }(x_t,t)$ 的前向投影来补全受金属影响的弦图并导出补全结果${\tilde{s}}_0^{t-1}$，它满足：

$${\tilde{s}}_0^{t-1}=\{\begin{array}{llc}{s}_0,& \text{where}& {\mathcal{M}}_s=0,\\ s_0^{prior},& \text{where}& {\mathcal{M}}_s=1,\end{array}\text{\hspace{1em}(20)}$$
其中$s_0^{prior}=\mathrm{FP}(f_{\theta }(x_t,t))$。然后，可通过FBP从${\tilde{s}}_0^{t-1}$近似重构出${\tilde{x}}_0^{t-1}$，
$${\tilde{x}}_0^{t-1}=\mathrm{FBP}({\tilde{s}}_0^{t-1}).\text{\hspace{1em}(21)}$$
${\tilde{x}}_0^{t-1}$可以看作近似于公式(19)。则根据公式(3)的后验分布，可由${\tilde{x}}_0^{t-1}$导出$x_{t-1}$。

通过上述方法，我们将扩散先验整合到MAR框架中，对应于图1(b)中的sinogram inpainting分量。

2) 图像域的融合策略：尽管使用公式(20)和公式(21)中概述的方法求解公式(19)是有效的，但扩散模型提供的先验弦图 ${\mathcal{M}}_s\odot \mathrm{FP}(f_{\theta }(x_t,t))$与已知弦图$(1-{\mathcal{M}}_s)\odot s_0$表现出不连续性。这导致在使用FBP重建时引入新的伪影，如图2所示。

因此，我们提出利用图像域的扩散先验进一步细化重构图像${\tilde{x}}_0^{t-1}$。首先，我们将时间步长为$t-1$时的扩散预测$x_0$，以及$f_{\theta }(x_t,t)$，以一定的权重掩模加入到${\tilde{x}}_0^{t-1}$中，

$${\tilde{x}}_0^{{}^{\prime }t-1}={\mathcal{M}}_f\odot f_{\theta }(x_t,t)+(1-{\mathcal{M}}_f)\odot {\tilde{x}}_0^{t-1},\text{\hspace{1em}(22)}$$
其中$M_f$表示扩散先验$f_{\theta }(x_t,t)$的权重掩模，这有助于纠正${\tilde{x}}_0^{t-1}$的非自然部分。

此外，我们提出进一步将原始的金属伪像$y_0$融合到${\tilde{x}}_0^{{}^{\prime }}$中。这可以提供额外的可能性信息。具体方法包括:
$$\begin{array}{rl}& {\tilde{y}}_0={\mathcal{M}}_y\odot f_{\theta }(x_t,t)+(1-{\mathcal{M}}_y)\odot y_0,\\ & \\ & {\tilde{x}}_0^{{}^{\prime \prime }t-1}=\sqrt{{\bar{\alpha }}_t}{\tilde{x}}_0^{{}^{\prime }t-1}+(1-\sqrt{{\bar{\alpha }}_t}){\tilde{y}}_0.\end{array}\text{\hspace{1em}(23)(24)}$$
公式(23)演示了使用一定权重掩模${\mathcal{M}}_y$将金属伪像图像与先验图像融合。公式(24)说明了使用特定比例与先验图像(如公式(23))和先验增强重建图像(来自公式(22))合并的金属伪影图像的组合。系数由扩散正过程的超参数推导而来。这表明随着t的减小，金属伪像的比例变小。上述图像域融合算法如图1（c）所示。在第IV.A节中，我们验证了将图像域融合策略引入MAR所带来的性能增强。

3) 时间动态权重掩模：选择合适的 ${\mathcal{M}}_f$和 ${\mathcal{M}}_y$具有挑战性。 根据经验，我们建议利用二进制弦图的重建结果来构建权重掩模，如下所示：
$$\mathcal{M}(\delta )=\mathrm{FBP}(\delta {\mathcal{M}}_s),\text{\hspace{1em}(25)}$$
其中$\delta$表示填入金属影响区域的值，而其他区域填入零。通常，随着$\delta$的增大，重构权重掩模的总体值也随之增大，如图3所示。

对于${\mathcal{M}}_y$，使用固定的$\delta$是合理的，因为在公式(23)$y_0$保持不变。因此，我们设它为，
$${\mathcal{M}}_y=\mathcal{M}({\delta }_y).\text{\hspace{1em}(26)}$$
然而，对于${\mathcal{M}}_f$，在不同的时间步长，已知弦图$(1-{\mathcal{M}}_s)\odot s_0$，和先验弦图 ${\mathcal{M}}_s\odot \mathrm{FP}(f_{\theta }(x_t,t))$之间的差异将导致图像域伪影的变化。通常，随着时间步长t的减小，扩散模型的生成变得更加逼真，从而减少了补全弦图的不连续，从而减少了图像域中的伪影。因此，设计一个随时间步长t变化的权重掩模更为合理，
$${\mathcal{M}}_f(t)=\mathcal{M}(\delta (t)),\text{\hspace{1em}(27)}$$
式中$\delta (t)$为随t变化的值。然后，将公式(22)中的${\mathcal{M}}_f$替换为${\mathcal{M}}_f(t)$，可以得到更精确的融合结果。在第IV-A.2节中，我们验证了时间动态Mf对MAR结果的积极影响。

III. IMPLEMENTATION DETAILS

A. 数据集

1) 训练数据集：为了训练我们的扩散模型，我们使用来自DeepLesion数据集[22]的不含金属伪影的CT图像。除去预留用于测试的200张CT图像，总共有927,802张512×512用于训练的CT图像。

2) 合成测试数据集：我们首先在一个合成数据集上评估了我们提出的方法以及比较方法，其中选择了来自DeepLesion数据集的200张干净的CT图像来模拟金属伪像CT图像。每张CT图像中分别植入10个金属掩模，采用扇束投影模拟得到金属影响弦图。总共获得了2000组CT 弦图和受金属伪影影响的图像，用于对MAR方法进行评价。掩膜来自CNNMAR[10]， 10个金属掩膜的尺寸为[2061,890,881,451,254,124,118,112,53,35]像素。扇形波束的几何形状遵循了Y等人先前的工作[23]。在扇束投影过程中，首先将CT图像调整为416x416。然后，在0 ~ 2π的角度范围内均匀地拍摄640个投影。detector有641个bin，因此投影尺寸为640x641。在模拟金属伪影过程中，考虑了多色X射线、部分体积效应、线束硬化和泊松噪声等效应。该合成数据集在前人的工作[5]-[9]，[14]中得到了广泛的应用。

3) 临床T检验数据集：我们进一步在CTPelvis1k数据集上对我们的方法进行检验[24]。CTPelvis1k的临床-金属子数据集包含有金属伪影的术后CT图像。首先将临床CT图像调整为416×416。然后，按照与合成DeepLesion数据相同的投影原理，得到临床数据的弦图。采用2500 HU的阈值对金属植入物进行分割，并通过正向投影得到金属植入物在弦图中的迹线。

B. 网络结构和扩散参数

我们选择patch diffusion[25]进行预训练，这是一种改进版的guided diffusion[26]。Patch diffusion算法首先将图像重构成一个由不重叠的Patch组成的网格，并在通道维度上进行拼接，然后学习重构后图像的生成。这种方法减少了生成高分辨率图像的计算需求。在我们的实验中，我们将512 × 512 × 1的图像分成形状为128 × 128 × 16的小块。

在扩散训练中，我们选择预测$x_0$的方法，其中在每次迭代中，网络输出$f_{\theta }(x_t,t)$来预测初始图像，如公式(11)。在网络架构方面，我们采用了引导扩散的U-net架构。表1列出了一些网络参数，其中Dim表示初始通道维度的数量，Muls ({m[1]，…， m[k]})表示通道乘数，N-Res表示残差块数，ResAttn表示使用自注意残差块的分辨率，Dropout表示丢弃率，ChHead表示每个头的通道数。

我们将批处理大小设置为16，学习率设置为1.0 × 10⁻⁴，并训练模型进行150,000次迭代。训练在单个RTX 3090 GPU上进行。

在采样阶段，我们采用加速技术来减少时间步数，以提高采样效率[26]。具体来说，我们对超参数${\beta }_1,\dots ,{\beta }_T$的正向过程，将它们的计数从1000减少到100。随后，我们在扩散推理中使用了新的超参数调度来实现加速。

C. 评估指标

对于地面真实值可用的模拟数据集，我们采用峰值信噪比(PSNR)和结构相似性指数(SSIM)作为评估指标。我们根据植入金属面具的大小将测试图像从大到小分为五组，并计算每组的平均度量。分组为[(2061)，(890,881)，(451,254)，(124,118,112)，(53,35)]。指标是用HU窗口[- 1000,4208]计算的。同时，我们展示了不同方法在减少不同尺寸金属植入物的金属伪影方面的视觉效果。

在真实临床数据集的情况下，缺乏ground truth限制了我们使用视觉评估作为评估算法性能的标准。

IV. EXPERIMENTS

A. 消融研究

在本节中，我们验证了II-B节中提出的不同模块的有效性，包括sinogram inpainting和image fusion。图像融合模块包括是否引入初始金属伪影(MA)图像。此外，我们还探讨了如何选择动态权重掩膜来进一步提高算法的性能。

1) 模块的验证：我们分别使用sinogram inpainting或image fusion策略，并在DeepLesion合成数据集上同时使用这两种策略，验证了MAR的性能。我们将图像域融合分为两部分进行验证：一部分表示为公式(22)，融合来自扩散模型的先验图像(表2中表示为“+先验图像”);另一个由公式(23)和公式(24)表示，涉及与初始金属伪像(表2中表示为“+MA图像”)的融合。表二给出了量化指标。

从表2可以看出，以sinogram inpainting作为基础可以初步减少金属伪影。图像域融合进一步提高了有效性，特别是在提高PSNR方面。很明显，我们提出的每个模块都有积极的影响。

值得注意的是，我们在表2的实验(b)中测试了仅使用金属伪像进行图像域融合而不进行sinogram inpainting的效果。该过程涉及根据公式(23)进行迭代，而不涉及公式(24)的融合。结果表明性能较差，这可归因于图像域金属伪影的复杂性。这说明了双域处理的必要性。从图4中，我们可以观察到仅采用弦图域处理和图像域处理均可成功恢复金属植入体周围的图像。然而，在sinogram inpainting后加入图像域融合可以改善金属区域周围的恢复。将扩散先验图像与原始金属伪影图像融合，可以获得最佳的视觉效果。

2) 动态权重掩膜的验证：在II-B.3部分引入了时间动态权重掩模${\mathcal{M}}_f$来融合扩散先验图像和补全弦图重建图像。在这里，我们测试如何选择时变的$\delta (t)$来生成动态掩模${\mathcal{M}}_f(t)=\mathcal{M}(\delta (t))$。经验表明，随着时间步长t的减小，图像域内的伪影变小。因此，我们设$\delta (t)$为与时间步长$t$正相关的函数，
$$\delta (t)=(a-1)e^{-n\frac{t}{T}}+1.\text{\hspace{1em}(28)}$$
通过这种设计，对于较大的n值，$\delta (t)$在t = t时接近1，而在t = 0时，$\delta (t)=a$。选择不同的a和n值会产生不同的函数，如图5所示，其中我们选择a =(0.3, 0.4, 0.5)和n =(3,4,5)，并绘制出从这些组合中获得的9个$\delta (t)$函数。

我们在合成DeepLesion数据集上测试了图5中9个$\delta (t)$函数对应的动态掩模。

表III给出了相应的PSNR值，其中a = 1.0表示定权掩模，因为$\delta (t)=1$。

很明显，与使用恒定掩模时的40.80dB的PSNR相比，使用所提出的动态权重掩模后，我们的方法的性能得到了进一步提高。此外，从表III可以看出，动态掩模对参数a和n具有很强的鲁棒性。

B. 比较方法

我们将所提出的方法与各种类型的MAR方法进行了比较，包括经典算法和深度学习算法，无监督方法和有监督方法。

经典的方法包括LI[3]和NMAR[4]，它们涉及在sinogram domain中补全金属影响区域。

无监督深度学习方法包括CycleGAN[16]、ADN[14]和Score-MAR[19]。其中，CycleGAN和ADN是利用未配对的干净图像和金属伪像图像在图像域学习风格迁移的方法。虽然CycleGAN不是专门为MAR设计的，但它可以应用于金属伪影还原。ADN是为MAR量身定制的，因此表现出优越的性能。Score-MAR采用预训练分数模型作为MAR的先验。其解决方案主要涉及将分数先验与已知的似然在弦图域融合。这与我们方法中的sinogram inpainting模块类似，主要区别在于它不是对预测的初始图像进行操作。此外，原始的Score-MAR工作主要侧重于引入求解逆问题的广义框架，因此其MAR的实验条件并不严格，例如仅使用平行光束投影。在Score-MAR的实现中，为了确保公平性，我们使用了与我们的方法相同的预训练模型。

此外，我们还比较了两种监督学习方法。第一种是CNNMAR[10]，这是一种较早的方法，它利用卷积神经网络来预测用于sinogram inpainting的先验图像。第二个是Indudonet+[9]，它展开了实现双域学习的重构过程。它目前在合成数据集上拥有最先进的性能。

C. 合成数据比较

我们首先在合成的DeepLesion数据集上对方法和我们的方法进行了比较。MAR结果的指标如表4所示。根据金属植入物的大小，将测试数据分为五组，计算每组的平均指标。最后，还计算了总体平均指标。

从表4可以看出，我们的方法在无监督算法中总体平均性能最好，超过了早期的监督方法CNNMAR，但与最先进的监督方法InDuDoNet+相比仍有一定差距。在无监督方法中，CycleGAN和ADN仅在图像域操作，在大型金属植入物上表现不佳，但在较小的金属植入物上表现较好。我们的方法得益于双域信息，在大型金属工件上表现更好。Score-MAR的有限性能归因于它只在sinogram inpainting中使用了分数先验。这突出了在双域中使用扩散先验的优点。

图6显示了不同方法MAR结果的视觉效果。

如图6所示，金属植入物会在CT图像中引入严重的条纹伪影。虽然LI和NMAR通过在弦图域中进行补全来显著减轻金属伪影，但由于在弦图中补全部分和已知部分之间的不连续性，它们引入了新的伪影。CNNMAR使用卷积神经网络对NMAR进行了改进，但其结果仍然显示出一些类似NMAR的伪影。InDuDoNet+通过双域监督学习取得了优异的结果，但在其输出图像中可以看到残留的条纹伪影，特别是如图6（b）所示。

CycleGAN和ADN在图像域使用无监督学习将金属伪像图像映射到干净图像。然而，由于金属伪影的复杂性和严重性，如图6所示，他们很难有效地去除条纹伪影，特别是对于较大的金属植入物。

Score-MAR和我们的方法都使用深度生成模型作为先验。他们的结果在视觉上比其他方法更真实。然而，Score-MAR仅使用分数先验进行sinogram inpainting，这导致金属植入物周围存在很大的误差。相比之下，我们的方法，通过迭代双域处理，恢复了金属植入物周围自然和真实的外观。

D. 临床数据比较

图7显示了不同方法在临床数据集上的结果。

骨盆CT片。由于金属植入物的存在，在金属区域周围观察到明显的伪影。然而，对远离金属植入物的区域的影响相对较小。从视觉上看，传统的LI和NMAR方法取得了很好的效果，如浅蓝色感兴趣区域(ROI)所示，它们在金属区域周围实现了良好的修复。同样，CNN增强型NMAR(简称CNNMAR)也表现出类似的效果。然而，这三种方法在远离金属植入物的区域引入了新的伪影，如黄色ROI所示。

纯基于图像的无监督方法，包括CycleGAN和ADN，在整体灰度预测中表现出明显的偏差。这种差异可能是由于学习到的图像域翻译与临床数据不一致造成的。此外，这些方法不能充分利用正弦图中存在的信息。虽然它们不会在黄色ROI区域引入新的伪影，但它们也无法去除金属植入物附近的伪影。

令人惊讶的是，InDuDoNet+在临床数据集上完全失败了。这可能是由于其复杂的架构过度拟合模拟数据，使其不适合泛化到不同的数据集。

Score-MAR在一定程度上减轻了金属假体的影响，但在金属假体周围仍有一个不自然的区域。与上述方法相比，我们的方法不仅成功地减少了金属周围的工件，而且避免了在其他区域引入新的工件。考虑到这些方面，我们的方法在所有比较方法中显示出最有利的视觉结果，突出了其在临床环境中的潜在适用性。

V. DISCUSSION AND CONCLUSION

A. 讨论

对于不同尺寸的金属植入物，采用不同的图像融合策略是进一步提高算法性能的一个潜在途径。从表4可以看出，在涉及小型金属植入物的实验中，我们算法的指标不如纯图像域算法ADN。这一观察结果表明，当处理较小的金属植入物时，金属伪影图像仍然包含大量的干净图像信息。因此，在图像融合阶段增加金属伪影图像的添加可能会导致算法性能的进一步提高。

我们引入了动态权值掩膜，将图像域的扩散先验纳入其中。动态掩膜的设计是为了校正在绘制弦图重建后新引入的伪影，它们与时间步长t相关。尽管动态掩膜提供了增强功能，但其范围仍然有限。因此，我们认为这些权重掩膜也可能是内容依赖的，从而导致不同图像的掩膜不同。这种方法可以提高算法的适应性，是一个值得研究的方面。

在以上两点的提示下，问题出现了:我们如何根据图像内容、金属植入物的形状和时间步长t推导出最合适的重量掩模?在这种情况下，手动设计函数本质上是次优的。因此，我们提出了利用成对数据学习图像融合的想法，将我们的方法转化为基于模型的监督学习方法。虽然这与我们最初的动机有些偏离，但可以预见，这种方法可以显著提高我们算法的性能。这可能会使我们的方法与最先进的监督学习算法(如InDuDoNet+)相媲美。这仍然是一条未来有待探索的道路。

B. 结论

本文提出了一种具有扩散先验的双域方法DuDoDp，用于金属伪影还原。首先，我们用干净的CT图像训练一个去噪扩散概率模型。为了将扩散引入到MAR的求解中，我们在扩散模型的每个时间步建立了一个优化问题，该优化问题的解对应于一个sinogram inpainting模块。在此基础上，我们提出了将重建图像、扩散先验图像和金属伪影图像进行融合，以更好地利用图像域信息。为了更好地控制图像融合的权重，我们设计了动态权重掩模来进一步提高性能。消融研究表明，我们的方法中提出的模块对MAR的结果有积极的贡献。合成和临床数据集的结果表明，我们的方法优于所有其他无监督方法，超过了早期的有监督的MAR方法。我们相信我们的工作对探索深度生成模型在MAR领域的应用具有重要意义。这可能会将金属伪影还原技术的使用扩展到配对训练数据稀缺的场景。

T图像训练一个去噪扩散概率模型。为了将扩散引入到MAR的求解中，我们在扩散模型的每个时间步建立了一个优化问题，该优化问题的解对应于一个sinogram inpainting模块。在此基础上，我们提出了将重建图像、扩散先验图像和金属伪影图像进行融合，以更好地利用图像域信息。为了更好地控制图像融合的权重，我们设计了动态权重掩模来进一步提高性能。消融研究表明，我们的方法中提出的模块对MAR的结果有积极的贡献。合成和临床数据集的结果表明，我们的方法优于所有其他无监督方法，超过了早期的有监督的MAR方法。我们相信我们的工作对探索深度生成模型在MAR领域的应用具有重要意义。这可能会将金属伪影还原技术的使用扩展到配对训练数据稀缺的场景。