读论文《A GRADIENT-BASED APPROACH TO FAST AND ACCURATE HEAD MOTION COMPENSATION IN CONE-BEAM CT》

请站在我身后

已于 2024-06-24 14:08:37 修改

阅读量530

点赞数 10

分类专栏：医学图像处理读论文文章标签：人工智能深度学习算法神经网络图像处理

于 2024-06-24 14:08:03 首次发布

本文链接：https://blog.csdn.net/buganything/article/details/139922680

版权

读论文同时被 2 个专栏收录

60 篇文章 2 订阅

订阅专栏

医学图像处理

32 篇文章 0 订阅

订阅专栏

题目：基于梯度的锥束ct头部运动快速精确补偿方法

地址：2401.09283 (arxiv.org)

目的：CBCT运动伪影去除

摘要：

锥形束计算机断层扫描(CBCT)系统，由于其便携性，提供了一个有前途的直接点护理医学成像的途径，特别是在危急情况下，如急性中风评估。然而，将CBCT整合到临床工作流程中面临着挑战，主要与扫描时间长导致患者在扫描过程中运动和导致重建体积中的图像质量下降有关。

本文介绍了一种基于梯度优化算法的CBCT运动估计新方法，该算法利用锥束CT几何形状的反投影算子的广义导数。在此基础上，建立了一个完全可微的目标函数，对当前运动估计在重建空间中的质量进行分级。我们大大加速了运动估计，与现有方法相比，速度提高了19倍。此外，我们研究了用于质量度量回归的网络架构，并提出了预测体素质量地图的建议，倾向于类似自编码器的架构，而不是收缩的架构。这种修改改善了梯度流，导致更准确的运动估计。通过头部解剖实验对该方法进行了验证。在运动补偿后，它实现了从初始平均3mm到0.61mm的重投影误差的减少，并且与现有方法相比始终表现出优越的性能。反投影运算的解析雅可比矩阵是所提出方法的核心，它是公开的。总之，本文提出了一种鲁棒运动估计方法，提高了效率和准确性，解决了时间敏感情况下的关键挑战，从而促进了CBCT与临床工作流程的整合。

Keywords computed tomography · deep learning · differentiable programming · gradient descent · motion compensation

1 Introduction

51%的患者在重建容积中表现出运动伪影，11%的患者由于图像恶化的严重程度而无法临床解释[9]。因此，为了将CBCT系统成功整合到危重患者或中风患者的临床工作流程中，一种防止运动伪影的可靠方法是必不可少的。

提出了一种基于软件的解决方案，该方案从测量数据中估计运动模式并补偿重建过程中的不对准。这些方法共享两步流程，包括(1)运动估计和(2)基于前一步估计的运动模式的运动补偿重建。图1说明了这种两步方法。

刚性运动补偿首先通过求解从测量数据估计运动模式的优化问题来实现。然后，根据估计的运动模式进行运动补偿重建，计算补偿体积

而运动补偿CT重建是相对直接的给定已知的运动模式，估计扫描和患者特定的运动是至关重要的一步。为了以基于图像的方式解决这个问题，必须从2D投影图像中导出3D患者运动。这是一个非平凡的反问题，因为它既取决于测量投影数据，也取决于扫描几何形状。因此，通常将其表述为优化问题，而不是直接估计。

在这项工作中，我们提出了一种新的运动估计方法，用于CBCT中的刚性帧间运动，该方法解决了上述两个缺点。首先，我们观察到所有现有的自动对焦运动补偿方法都使用无梯度优化算法，这导致目标函数评估次数高，计算时间长。相反，我们建议用基于梯度的优化器代替优化器。为此，在我们之前的工作[23,24]的基础上，我们将重构体积的解析导数的推导和实现推广到从扇形梁到锥形梁几何形状的几何参数。这使我们能够为CBCT制定一个完全可区分的自动对焦型目标函数，该函数表示中间重建的质量。由于它是可微的，因此可以得到该目标函数关于运动参数(几何参数的子群)的梯度，并将其用于基于梯度的优化，从而加快运动估计。

此外，我们观察到大多数现有的图像质量度量回归网络使用了收缩架构。空间维度逐渐降低，直到在输出处只回归一个标量值[6,19,22]。然而，在基于梯度的情况下，这种结构并不理想，因为标量输出值相对于体积的梯度倾向于表现不稳定，并且在整个输入端具有不均匀分布的振幅。因此，我们建议回归体素质量地图，从而促进类似自编码器的架构，而不是收缩图像质量度量。我们表明，这些导致更有信息量的梯度信息和改进的整体运动估计。

综上所述，本文的贡献有

•从扇形梁到锥形梁几何的CT重建的解析导数的推广。
•基于梯度优化的运动估计在临床CBCT应用中的应用。
•改善梯度流的体素级质量指标的回归

2 Methods

2.1 Rigid motion compensation for CBCT

在CBCT重建中，从不同方向测量的线积分被转换成空间分辨的x射线衰减系数图。这个过程从根本上依赖于描述x射线源、平板探测器和患者之间三维关系的扫描几何的精确知识。通常，在扫描过程中假定病人是静止的。在这种情况下，CBCT扫描仪的任何可复制的圆形扫描轨迹都具有一组校准的投影矩阵[25]的独特特征，这些矩阵将3D世界坐标系中的一个点映射到特定投影视图的探测器坐标系中的一个点上。给定与投影数据对应的精确投影矩阵，存在许多解析和代数算法来解决重建问题，患者的运动引入了校准投影矩阵和投影数据之间的不匹配。在重建过程中，这种不匹配会导致最终图像中的伪影。

在CBCT上可能出现的运动模式的特征取决于解剖结构。虽然人体的某些区域倾向于非刚性运动，但我们假设头部的刚性运动模式是不可变形的。在这种情况下，运动模式可以由三个旋转和三个平移分量参数化。此外，我们要求运动的频率与单个投影的曝光时间相比要低。这证明了运动只发生在投影之间(帧间)而不是在获取一个投影期间(帧内)的假设是正确的。这是一个有效的假设，因为每次投影的典型曝光时间约为5毫秒，这也符合之前的工作

我们的目标是通过更新投影矩阵来消除投影矩阵与刚体、帧间患者运动所带来的投影数据之间的不匹配。这是任何从观察病人的刚性运动在扫描仪相同可以逆建模的刚性变换扫描仪本身。因此，对于每个投影j，我们寻求估计刚性变换矩阵 $T_j \in R^{4 \times4}$ ,与初始投影矩阵 $P_j \in R^{3 \times 4}$ 相乘，得到更新后的投影矩阵 $\hat P_j \in R^{3\times 4}$ .

$\hat P_j = P_j \cdot T_j(X)$

变换矩阵 $T_j$ 由自由参数向量 $x \in R^N$ 参数化，其中N为自由参数个数(见2.2.1节)。然后将校正后的投影矩阵与原始投影数据一起重建，从而得到运动校正后的重建。该方法的关键部分是变换矩阵 $T_j$ 的估计，变换矩阵 $T_j$ 依赖于投影数据，并且是特定于扫描的。

在这里，我们解决了一个形式为 $x^* =argmin_xf(x) ,x^* \in R^N$ 的无约束优化问题，其中x∗是最优的自由参数集，它理想地补偿了病人的运动，并产生了一致的投影数据集和投影矩阵 $P^*_j =P_j \cdot T_j(x^*)$ 下面，我们将详细介绍我们提出的表达目标函数f(x)和解决优化问题的方法。

2.2 可微分的自动对焦管道

我们提出了一个自动对焦启发的目标函数用于运动估计，它产生描述患者运动的自由参数的信息梯度更新。目标函数仅依赖于投影数据和一组校准的投影矩阵作为扫描几何的初始化。完整的目标函数 $f(x) :R^N \rightarrow R$ ,由图2中突出显示的三个不同的处理步骤组成。

CBCT中刚性运动估计方法综述。三个基本构建块是(1)运动模型，(2)可微的反向投影，实现端到端梯度流，以及(3)经过训练的质量度量，用于回归空间分辨质量映射。

描述运动轨迹的自由参数 $x \in R^N$ 由运动模型 $R^N \rightarrow R^M$ 处理，由初始标定矩阵和x描述的运动模式计算投影矩阵。这里，M是扫描的所有投影矩阵中元素的个数。然后，反投影 $r :R^M \rightarrow R^K$ 使用更新后的投影矩阵将过滤后的投影数据解析反投影到一个有K个元素的体中。请注意，为了简单起见，我们省略了过滤后的投影数据作为反投影函数r的输入。在给定当前自由参数估计的情况下，反投影会导致数据的中间重建。

最后，质量度量 $q:R^K \rightarrow R$ 返回给定中间重构的标量质量值。在我们提出的方法中，我们首先回归与输入体积相同维度的完整质量图，然后将其平均为单个值。完整的目标函数可以写成 $f(x) =q(r(p(x)))$ 。为了计算目标函数对自由参数的梯度，我们应用了链式法则

$\frac {df}{dx} =\frac {dq}{dr} \cdot \frac {dr}{dp} \cdot \frac {dp}{dx}$

这里是dq / dx 是质量度量到中间重构的梯度。dr /dp为反投影函数的雅可比矩阵，dp/dx为运动模型的雅可比矩阵。在下面的2.2.1、2.2.2和2.2.3节中，我们将详细解释这三个步骤，并介绍如何获得相应的导数。此外，在第2.2.3节中，我们将强调为什么在基于梯度的设置中，空间分辨率质量图的回归后跟平均算子是优越的.

2.2.1运动模型

一般来说，任何CBCT扫描的投影几何是由一组投影矩阵 $P_j,j=1, ...,N_p$ 唯一定义的，其中Np是一次扫描中投影的个数。这是描述无约束CBCT轨迹的最灵活的形式。每个投影矩阵只有在比例上是唯一的，并不是每个3 × 4矩阵都是有效的投影矩阵。因此，不建议直接对投影矩阵条目进行优化。相反，对于刚性运动补偿，先验知识将每个视图的自由度限制为6个。在没有进一步约束的情况下，每个投影的刚体运动状态与相邻视图无关。然而，在典型的CBCT扫描中，每秒可以获得多个投影图像。

因此，我们通过分别将Akima样条拟合到六个运动参数 $(t_x, t_y, t_z, r_x, r_y, r_z)$ 中的每一个来对运动模式施加额外的时间平滑约束。样条曲线由 $N_n$ 个节点参数化。假设 $N_n < N_p$ ，这将自由度从 $6N_p$ 降低到 $6N_n$ ，并强制每个运动参数在整个扫描过程中平滑变化。这是一个合理的假设，因为典型的头部运动类型被认为是缓慢的振荡，以及从初始位置均匀增加的偏差[8]。此外，我们的参数化与先前的工作一致，这些工作同样依赖于Akima样条[19]、其他类型的光滑函数[11,28]或目标函数中的正则化项，这些正则化项隐式地惩罚非光滑运动模式

作为运动模型，我们指的是将 $N=6N_n$ 个参数(存储在向量x中并描述每个运动参数的Akima样条节点值)转换为 $M=12N_p$ 个投影矩阵条目的函数。这是通过计算每个投影时间点的样条曲线并根据公式1更新一组初始投影矩阵来实现的。由于这是一个相对便宜的计算，我们在原生PyTorch中实现样条曲线并自动获得相应的导数。如果没有不同的说明，我们使用10个节点来模拟人工运动模式，30个节点用于通过我们的方法估计的运动模式。节点在扫描的时间方向上均匀间隔，一个节点在扫描范围的开始，一个节点在扫描范围的结束。

2.2.2 CBCT反投影算子的解析导数

给出运动模型返回的迭代更新的投影矩阵，实现了解析反投影来执行中间重建。除了投影矩阵之外，反向投影还依赖于投影图像的适当过滤版本。由于我们考虑在全圆轨迹上获得的投影数据，因此我们采用经典的移位不变斜坡滤波器。在反向投影过程中，世界坐标中位置 $p^{world} = (x,y,z)^T$ 的重构值I被计算为每个投影图像中正向投影位置的滤波投影数据值的和

$I(\tilde{P}^{world})= \sum^{N_p}_{j=1} d_j(g(s_j(\tilde{P}^{world})))$

这里，Np是扫描中投影的个数，而 $\tilde{P}^{world} = (x,u,z,1)^T \in P^{3+}$ 是 ${P}^{world}$ 的齐次表示，遵循与Aichert等人[13]中使用的齐次坐标相同的表示法。函数 $s_j :P^{3+} \rightarrow P^{2+}$ 通过与第j个投影矩阵 $s_j(\dot{\tilde{P}}) = P_j\cdot \tilde p$ 相乘，计算第j个投影图像的探测器坐标系中 $\tilde p ^{world}$ 的正投影齐次位置。随后， $g:P^{2+} \rightarrow R^{2}$ 通过将所有元素除以向量的最后一个分量，从齐次坐标恢复欧几里得坐标。最后， $d_j:R^{2} \rightarrow R^{}$ 在正向投影位置双线性插值第j个滤波后的投影图像，并返回对应的值，该值对所有Np个投影求和，从而得到最终的反向投影值I。

为了重建一个完整的体积，公式3需要应用到世界坐标中的所有体素位置，这在顺序评估时非常耗时，并且禁止直接使用自动微分。出于这个原因，我们实现了一个自定义的gpu加速版本的反向投影，并推导并实现了导数的相应解析表达式。这允许我们使用gpu加速实现反向投影和梯度计算，并且仍然以可微的方式将重建图像和投影矩阵连接起来。在我们之前的论文中，我们提出了一个类似的扇束CT反投影的可微公式[23]。然而，它并不适用于锥形梁的几何形状，因此，它的实际适用性是有限的。下面，我们将这个公式推广到锥梁几何。

我们给出了位置 $p^{world} = (x,y,z)^T$ 的反投影值对第γ投影矩阵Pγ的12个分量的梯度表达式。为此，我们定义 $\tilde P_\gamma ^{detector} = s_\gamma (\tilde p^{world}) = (u_\gamma,v_\gamma,w_\gamma)^T \in P^{2+}$ 作为检测器坐标中的齐次正投影位置。梯度 $G_\gamma \in R^{3\times 4}$ 为

为了得到 $g_{u,\gamma} \in R ,G_{v,\gamma} \in R$ ，计算滤波后的投影数据沿探测器坐标系两轴的空间导数。它是数值近似的二阶精确的中心差沿水平和垂直轴的每个过滤投影图像产生两个梯度图像每个视图。这些是在欧几里得检测器位置 $(\frac {u_\gamma} {w_\gamma},\frac {v_\gamma} {w_\gamma})^T \in R^2$ 的双线性插值。γ-th梯度图像相对于第一和第二探测器轴的插值值分别表示为 $g_{u,\gamma} \in R ,G_{v,\gamma} \in R$ 。所有的研究都局限于离散的问题设置，其中投影图像在离散的角度和离散的探测器坐标系中测量。我们参考了之前的工作，对连续推导和相应的离散化版本进行了更深入的比较[23]。

方程4给出了体积中一个位置 $p^{world}$ 的反向投影值对一个投影矩阵Pγ的12个分量的偏导数。它需要对所有K个离散体积位置以及所有Np =M/12投影矩阵进行评估，以填充完整的雅可比矩阵dr dp。gpu加速的反向投影和相应的梯度计算的实现可以在https://github.com/mareikethies/geometry_ gradients_CT上公开获得。所有计算都封装在自定义PyTorch函数中，该函数与库的自动微分功能直接兼容。由于内存的限制，我们从来没有显式地计算完整的雅可比矩阵，而是直接计算雅可比矩阵向量积，给出了下游操作的梯度信号。

2.2.3逐像素质量度量回归

给定从反向投影的中间重建，我们试图根据运动伪影来量化其质量。在推理时，该映射需要是无参考的，这意味着仅在给定受运动影响的体积而没有无运动参考的情况下估计质量。为了解决这个问题，神经网络在模拟的运动影响和相应的无运动扫描的成对数据集上进行训练，以近似给定运动影响体积的基于参考的目标度量。这在之前的工作中也有类似的建议[6,19,22]。与现有方法相比，我们在基于梯度的设置中使用训练好的网络，其中计算网络输出对输入体积的导数以进行运动估计。对于无梯度方法，从输入体积中准确回归质量度量是关键部分。在基于梯度的情况下，通过训练网络反向传播获得的信息量体梯度与前向映射同样重要。为了确保获得信息丰富的体梯度，我们建议训练一个回归完整空间分辨质量图的体到体网络。这与使用契约架构直接回归标量质量度量的现有方法形成对比。为了获得用于优化的标量，我们简单地对体积质量图进行平均，确保每个空间位置对最终值的贡献相等。

我们选择视觉信息保真度(VIF)作为回归目标[30]。这个指标量化了参考图像中存在的信息，以及可以归因于某些失真过程的信息损失，在我们的情况下是运动。它明确地模拟了人类视觉系统，并被认为提供了扭曲图像和参考图像之间的图像相似性度量，这与人类的感知非常一致[6]。VIF的原始公式计算一个标量，其中0为下界，表示参考图像由于失真而丢失了所有信息，1表示失真图像包含与参考图像相同的信息，即它们是相同的[30]。理论上，当考虑的图像具有比参考图像更高的信息含量时，VIF可以大于1。为了获得VIF的空间分辨体积图，我们遵循Shao等人的方法。[31]

给定一对无运动和受运动影响的头部CT体数据集，通过将无运动体设置为参考图像，将受运动影响的体设置为扭曲图像，计算出相应的VIF映射。VIF地图通过体素K的数量进一步缩放，使得平均操作产生的值在[0,1]范围内。作为网络的最终回归目标，我们使用 $VIF^*(I_{dist},I_ref) = 1-K \cdot VIF(I_{dist},I_{ref})$ 。然后，我们训练一个3D U-net架构[32]，从受运动影响的体积中回归调整后的VIF地图VIF∗。U-net由具有三维卷积和ReLU激活函数的基本构建块组成。每层特征映射的数量为8l，其中l = 1，…4的水平。因为我们训练的是一个回归任务，所以最后一层是1×1卷积，没有使用进一步的最终激活函数。该模型使用L1-loss和Adam优化器进行训练，学习率为0.001，批大小为16，输入和输出体积为128 × 128 × 128。训练后，对预测的映射进行简单的平均，得到一个标量值。此外，训练后的网络执行无参考映射，不再需要无运动的体积

2.3基于梯度的优化算法

将预测质量图的平均值最小化以获得扫描特定运动估计x *。重要的是，在优化时，质量度量网络的参数是完全训练和冻结的。大多数现有的运动补偿算法利用无梯度优化器来最小化相对于运动参数的目标。对于无梯度优化算法，只需指定目标函数就足够了，而基于梯度的算法需要梯度的表达式，由于目标值对运动参数的依赖性并不小，因此梯度的表达式很难表述[33]。我们是第一个使用完全可微分的自动对焦目标函数，允许我们应用基于梯度的优化。具体来说，我们利用了一个基本的梯度下降更新方案的形式

$x^{(n+1)} =x^{(n)}-s^{(n)}\cdot \frac {df}{dx}$

更新的估计 $x^{(n)}$ 是从初始估计 $x^{(0)} = 0$ 开始迭代计算的，通过步长s(n)∈R向负梯度方向迈出一步，步长受到指数衰减的影响，使得 $s^{(n)} = s_0 \cdot t^n$ ，其中t∈R是衰减因子，而s0∈R描述初始步长。如果没有特别说明，我们运行算法100次迭代，初始步长为50 = 100，衰减因子为t = 0.97。

3 Data

我们的实验使用的是CQ500数据集[34]。最初，该数据集包括在不同MDCT扫描仪上获得的491个重建头部CT扫描。我们首先对数据集进行过滤，以小切片厚度重建扫描。从剩余的扫描中，我们进一步排除了那些比平均样本更少或更多切片的扫描。这导致了320次扫描，我们在患者水平上依次将其分为训练集(150)、验证集(50)和测试集(120)。滤波后，所有扫描的重建层厚为0.625mm，各向同性层内间距为0.38 ~ 0.58mm(平均0.472mm)。因为我们研究的是锥束几何，所以我们模拟了重建MDCT体中相应的探测器数据。使用PyroNN软件[35]，将每个体块在锥束几何的圆形轨迹上进行正投影，在完整的2π角范围内进行360次投影，源到等心距离为785mm，源到探测器距离为1200mm，探测器形状为500 × 700像素，各向同性像素间距(0.64mm)。斜坡滤波器和余弦滤波器应用于计算投影图像。

为了训练质量度量网络(第2.2.3节)，我们在每次用于训练的样本中动态采样一个新的随机运动扰动，并将其应用于投影矩阵。基于样条的运动模型(章节2.2.1)，每个样条有10个节点，平移参数的最大振幅为10mm，旋转参数的最大振幅为15◦。为了确保在优化过程中可能遇到的所有运动状态都在训练数据中表示出来，我们包括了不同运动参数中振幅不等的运动模式，以及只对数据产生轻微扰动的运动模式。所有样条都是单独的零中心，并使用运动模型转换为每个turbed投影矩阵。过滤后的投影数据从这些扰动矩阵中重建，训练目标VIF∗从扰动和非扰动重建中计算，如第2.2.3节所述。在将重建体的强度馈送到网络之前，用固定的、与样本无关的偏移量和斜率进行归一化，使其近似位于区间[0,1]。

4 Experiments and results

所有的运动补偿实验都是在质量度量模型训练测试集中的30名患者身上进行的。随机运动模式为每个患者采样，平移幅度为5mm，旋转幅度为5◦，在不同的方法和优化算法中保持恒定。运动估计本身运行在大小为128 × 128 × 128、重构信号间距为2mm的网格上，但基于图像的结果是在大小为1mm的256 × 256 × 256体素的更高分辨率信号上计算的，用于最终的运动补偿重建。为了评估，我们将所有运动补偿重建严格注册到各自的三维地面真值重建中。文中所有的箱形图都显示了中位数和四分位数范围以及最小值和最大值。异常值用交叉标记突出显示。运动补偿前的初始度量值用灰框表示。

4.1运动补偿性能

为了评估我们提出的方法的运动补偿性能，我们将其与之前发表的两种方法进行了比较。为了公平比较，所有方法中的数据设置、运动模型和反向投影步骤都是相同的，但只有质量度量被替换。首先，选择3D总变差(TV)作为经典的自动对焦质量度量，通过最小化l1范数来强制稀疏体积梯度。它是一个无参数的质量度量，并已被证明可以为图像对齐提供信息[19,21]。其次，我们重新实现Huang等人提出的质量度量网络，该网络也使用3D架构回归VIF度量[6]。与我们的质量度量相比，它不是在空间分辨率的质量映射上训练的，但作者使用了一个收缩残差卷积架构，然后是一个完全连接的层，用于最终的回归任务。我们在与我们提出的网络相同的数据上训练网络，但直接计算VIF的空间平均值作为训练目标。学习率设置为0.0001，权重通过Adam算法优化500 epoch。

从图3中可以看出，我们的方法在最终运动补偿体积上计算的度量方面始终获得最佳结果。在基于梯度优化的蓝盒图中，只有我们提出的方法在均方根误差(RMSE)、SSIM和VIF方面持续改进了初始运动影响重建(灰盒)，测试集中所有情况的SSIM都超过了0.90。相比之下，其他两种方法都不能在初始状态上稳定地改进。这种观察结果的原因可以在图4中确定。虽然这三种方法都成功地改善了平面外参数的初始运动影响状态，但对于平面内参数却并非如此。使用TV作为质量度量的优化导致特别大的面内误差，通常大于初始状态，同时导致准确拟合面外参数。Huang等人的方法改进了面内和面外参数，但在某些情况下，改进是微不足道的。对于我们的方法，在所有情况下，面外运动几乎完全恢复，绝对误差低于0.5mm和0.5◦，而面内参数tx和ty的中位数绝对误差超过0.5mm。在图5中，我们分析了重投影误差(RPE)，这是通过使用恢复和目标几何形状将一组固定的300个三维点(半径为25mm, 50mm和100mm)在等心周围正投影到探测器平面上来计算的。我们的方法优于两种比较方法，在所有情况下，从初始的中位数RPE约为3mm提高到RPE低于1mm。表1总结了基于梯度的情况下所有指标的平均值。

图6通过两个定性例子验证了我们的发现。重建的切片证实了我们的方法成功地补偿了运动伪影。我们提出的算法的输出在视觉上接近地面真实。Huang等人提出的方法可以减少运动伪影，但仍然存在可见的条纹和模糊的边缘，特别是在例2中。对于电视质量度量，我们可以直观地确认之前的结果，即它导致平面外运动(主要在矢状视图中可见)的可接受结果，但无法恢复平面内运动(参见轴向视图)。在每个切片的放大感兴趣区域(ROI)中可以看到一个小的骨结构。只有我们的方法才能完全恢复这些roi的所有细节，并成功地去除条纹和模糊的边缘。

图7显示了图6中第一个定性实例在扫描范围内的运动模式。该方法估计的运动曲线最接近地面真实运动曲线，但平面内平移tx和ty的误差略高于其他四个参数。参考方法与地面真实值偏差较大，特别是用电视质量度量估计的面内运动误差较大。这些发现与图6中重建图像的特征相吻合。

4.2 Optimizer comparison

为了分析基于梯度优化的影响，我们在基于梯度和无梯度的情况下比较了我们的方法以及两种现有的方法。对于无梯度设置，我们选择了协方差矩阵自适应进化策略(CMAES)算法[36]，因为在相关工作中，CMAES算法经常应用于类似问题[18,22]。它配置了平移和旋转的初始标准偏差分别为0.5mm和0.5◦，并迭代直到目标函数评估的数量达到10000。

图8给出了两种优化设置下所需的平均运行时间。当使用建议的基于梯度的更新时，所有三种研究的方法都需要不到两分钟的时间，而当使用CMAES时，相同的优化需要更长的时间(超过30分钟)。可以观察到，不同的质量指标在评估和梯度计算所需的时间上略有不同，但与两种优化策略之间的差异相比，这种影响可以忽略不计。例如，对于我们提出的质量度量，CMAES优化大约比梯度下降算法慢19倍。除了质量度量和相同的硬件(Nvidia A40 GPU)，所有运行都使用相同的实现来执行所有处理步骤。此外，从图3、图4和图5中，我们可以观察到优化算法的选择并不影响我们提出的方法的运动估计性能。当使用梯度下降时，与使用CMAES进行质量度量相比，所有度量至少是相同的(面向面外运动参数的RPE或MAE)，甚至稍微更好(面向面内运动参数的基于图像的度量或MAE)。我们得出的结论是，基于梯度的优化将所有三种研究方法的速度提高了约19倍。同时，神经网络参数化的质量指标对运动估计的性能没有影响，甚至略有提高。

4.3质量度量选择

在我们提出的方法中，我们使用体素化的VIF地图作为训练质量度量的回归目标。最近的研究表明，VIF是运动补偿的合适度量[6]。然而，任何可以以每体素的方式计算的度量都是我们算法的有效选择。因此，我们通过训练质量度量来回归均方误差(MSE)以及来自运动影响输入的结构相似指数度量(SSIM)来评估目标度量对我们方法性能的影响。对于SSIM，与VIF所描述的策略类似，我们通过SSIM∗(Idist,Iref) = 1−SSIM(Idist,Iref)对调整后的映射进行回归，确保平均值为0对应于理想图像。除了回归目标的计算外，训练设置没有改变。

我们观察到，经过训练分别预测MSE、SSIM和VIF质量地图的所有三种网络都能产生准确的运动估计，并且在RPE方面优于比较方法。VIF的RPE最小，平均为0.61mm，而SSIM和MSE的平均RPE均为0.98mm。因此，VIF在给定的实验设置中工作最准确，但与MSE和SSIM相比只有很小的差距

5 Discussion

实验表明，利用梯度下降代替无梯度CMAES进行优化，运动估计速度提高了19倍。在临床实践中，当运动补偿重建需要在几分钟内可用时，这种运行时间的减少是相关的。例如，在血管造影室的中风成像中，治疗的时间是至关重要的。因此，任何由图像处理引起的延迟都是要避免的。此外，像我们这样的快速算法甚至可以想象在干预设置中。

在运动补偿性能方面，我们的方法明显优于已有的研究方法。总变分是一个简单的非参数函数，适合具有稀疏梯度的体积。然而，它是不可知的内容的体积，在我们的情况下，头部的解剖结构。因此，它可以收敛到表示不合理解剖的解。此外，TV可能更容易受到初始化的影响，并在较大的初始运动模式下发生发散。因此，最近的研究集中在学习质量指标上，这些指标是在目标解剖结构的考试样本上训练的，从而获得对目标解剖结构的隐含理解。Huang等人的工作就是这种方法的一个例子。我们假设他们的方法在我们的实验中不能产生令人满意的结果，因为我们拟合了相当多的自由参数(180个，而他们最初发表的自由参数为72个)。这增加了搜索空间的复杂性，使优化任务更加困难。特别是无梯度优化算法，对于大量的自由参数表现出不稳定和缓慢的收敛性。因此，我们认为基于梯度的优化不仅在速度方面具有优势，而且在增加自由参数数量时也具有优势。

在图4a中，我们观察到所有方法对于平面内平移tx和ty的MAE都高于其他运动参数。这至少部分是由于，相对于扫描轨迹的坐标系的定义。在每个全圆扫描中，两个平面内平移平行于中心射线。重建图像不依赖于平移接近平行于中心射线，因此，所有研究的目标函数对这些参数不像对其他运动参数那样敏感。然而，由于tx和ty的小误差只要与测量射线大致平行，就不会影响最终的重建，因此我们认为这些偏差是可以接受的。

该方法在真实CT重建图像上进行了测试，模拟了锥束图像。现有的关于成像设备中头部运动特征的研究提供了真实的运动模式[8]。在真实的、受运动影响的锥束数据上测试该方法将是一个有价值的实验。这种定量研究的理想数据集由锥束CT测量和相应的运动影响投影数据组成，这些数据可获得地面真实运动模式，例如，在数据采集期间，来自附着在患者头部的光学跟踪系统。我们不知道有公开可用的数据集满足这些要求，但认为它对未来的研究非常有价值。此外，我们认为通过模拟来自重建图像的锥束投影数据不会引入干扰域间隙，因为我们从未在我们的方法中主动使用投影数据。虽然已知模拟的x射线投影与真实的x射线投影具有不同的特征[37]，但这种差异在重建步骤中消失，因为重建算法将模拟过程中假设的理想图像形成反转。由于我们的神经网络是在重建图像上训练的，我们预计在真实CBCT图上的性能不会有实质性的下降。

梯度下降依赖于精确的梯度信息。我们推导并实现了一个解析表达式，将重构体上的梯度转化为几何参数上的梯度。为了获得精确的几何梯度，我们需要信息体梯度，这些体梯度是通过质量度量通过标准反向传播获得的。由于质量度量通常只在前向映射上进行训练，因此不能保证来自反向传播的梯度信息对运动补偿具有信息性。例如，在收缩架构中，学习到的函数可能只依赖于体素的子集，导致那些对输出没有贡献的位置的梯度信息丢失。因此，我们对全质量地图进行回归，从而对体积中的每个体素进行质量度量预测。因此，学习到的映射被强制合并所有体素位置和梯度流回体的所有部分，这确保了更均匀分布的梯度信息。此外，已知U-net架构中的跳过连接可以改善梯度流并缓解梯度消失问题[38]。

用于将梯度从重构体传播到几何参数的解析公式建立在先前发表的工作基础上[23]。这项先前的工作仅限于扇形光束几何形状，与临床应用的相关性有限。在这里，这个想法被扩展到锥束CT扫描仪。此扩展需要仔细实现反向投影和相应的梯度计算，以保持内存需求可管理。虽然明确计算扇形波束几何形状的整个雅可比矩阵dr dp是可行的，但对于大小为256 × 256 × 256的360投影矩阵的体积，这将需要大约290GB。这与现有的消费类gpu不兼容。因此，在锥形梁的情况下，雅可比矩阵从来没有明确计算，但雅可比矩阵向量乘积是使用相对于重建体积的质量度量的梯度实时评估。此外，我们之前的工作并没有纳入运动模型，导致不切实际的快速和暂时不受约束的模拟运动模式，这与有关头部运动特征的现有知识相矛盾。

此外，我们想强调的是，反向投影操作的几何可微实现是公开提供给社区的，作为一个现成的PyTorch函数，可以插入到任何现有的自动对焦启发的方法中，以启用基于梯度的优化。它只需要一个可微分的质量度量。由于大多数最先进的方法使用神经网络进行质量度量回归，因此通过反向传播确保了可微分性。因此，基于梯度的优化不仅适用于这项工作的特定背景，而且有可能加快现有的基于自动对焦的工作。

一些作者认为，由于在图像对齐和运动补偿中遇到的优化问题通常不可能是凸的，因此首选无梯度优化算法[18]。假设基于梯度的算法比无梯度的算法更有可能终止于次优的局部最小值，因此可以得到更精确的解。在我们的实验中，我们不能证实这种行为。我们承认，对于电视质量度量，梯度下降比CMAES表现得更差，但我们提出的质量度量的优化策略与huang等人提出的优化策略之间没有系统差异。如果质量度量和运动模型的结合导致了一个强大的非凸优化景观，我们可以想象一种方法，使用梯度下降进行快速初始对齐，然后使用无梯度算法(如CMAES)对估计的运动模式进行微调，这可能会克服由于其随机行为而导致的一些局部最小值。此外，在深度学习社区中，围绕基于梯度的优化算法进行了研究，该算法在非凸优化环境中具有鲁棒性。例如，这些算法使用附加动量项或自动参数依赖学习率自适应[39]。用一种针对非凸性的基于梯度的优化算法取代普通的梯度下降算法在强非凸场景中也很有前景。这对未来的工作来说是一个有趣的途径。

后续研究的另一个有价值的方向是使用解析反投影型算法将所提出的方法扩展到非圆轨迹[40]。这不仅允许对这些类型的轨迹进行运动补偿，而且可以设想一种类似的方法，用于对无法预先校准的非圆轨迹进行无影、基于图像的校准。