论文阅读《Direct Sparse Odometry》1

摘要

直接稀疏里程计方法(Direct Sparse Odometry, DSO)是一种基于一种新颖、高精度的稀疏直接结构和运动公式的视觉里程计方法。它结合了一个完全直接的概率模型(最小化光度误差)与所有模型参数的一致和联合优化，包括几何(在参考帧中以逆深度表示)和相机运动。这是通过丢弃在其它直接方法中使用的平滑度先验来实现的，本文在整个图像中均匀采样像素。由于我们的方法不依赖于关键点检测器或描述子，它可以自然地从所有具有灰度梯度的图像区域采样像素，包括边缘或平滑灰度变化的基本无特征的墙壁。该模型集成了完整的光度标定，考虑了曝光时间、透镜渐晕和非线性响应函数。我们在包含几个小时视频的三个不同数据集上彻底评估了我们的方法。实验表明，所提出的方法在各种真实世界环境中，无论是在跟踪精度还是鲁棒性方面，都明显优于最先进的直接和间接方法。

1 介绍

同时定位与建图(SLAM)和视觉里程计(VO)是许多新兴技术的基本组成部分——从自动驾驶汽车和无人机到虚拟现实和增强现实。近年来，SLAM和VO的实时方法取得了显著进展。长期以来，该领域以基于特征(间接)的方法为主，近年来，许多不同的方法得到了普及，即直接和稠密的表达。

直接对比间接。所有公式的基础是一个概率模型，它以噪声测量$Y$为输入，并为未知和隐藏的模型参数(3D世界模型和相机运动)计算一个估计器$X$。通常使用最大似然方法，它找到最大概率获得实际测量值的模型参数，即，$X^{\ast }:=argma{x}_X$ P(Y|X)。

间接方法分两步进行。首先，对原始传感器测量数据进行预处理，生成一个中间表示，解决了部分整体问题，如计算图像中对应点的坐标。其次，将计算得到的中间值解释为概率模型中的噪声测量值$Y$，以估计几何和相机运动。注意，第一步通常是通过提取和匹配稀疏的关键点来实现的，但是也存在其它选择，比如以密集的正则化的光流的形式建立对应关系。它还可以包括提取和匹配其它几何基元(如线段或曲线段)的参数表示的方法。

直接方法跳过这个预计算步骤，直接使用实际的传感器值——在特定时间周期内从特定方向接收到的光——作为概率模型中的测量$Y$。

在被动视觉的情况下，直接方法因此优化光度误差，因为传感器提供光度测量。另一方面，间接方法优化了几何误差，因为预先计算的值-位置或流-矢量都是几何量。请注意，对于其它传感器形式，如深度相机或激光扫描仪(直接测量几何量)，直接公式也可以优化几何误差。

稠密对比稀疏。稀疏方法只使用和重构一组选定的独立点(传统的角点)，而密集方法试图使用和重构二维图像域中的所有像素。中间的方法(半密集)避免重建完整的曲面，但仍然旨在使用和重建一个(大部分连接和良好约束)子集。

然而，除了所使用的图像区域的范围之外，一个更根本的也是必然的区别在于增加了一个几何先验。在稀疏表述中，没有邻域的概念，几何参数(关键点位置)是有条件地独立于给定的相机位姿和内参。另一方面，稠密(或半稠密)方法利用所使用的图像区域的连通性来形成一个几何先验，通常有利于平滑。事实上，这样的先验是必要的，以使稠密的世界模型，仅从被动视觉观察。一般来说，这个先验是直接以额外的对数似然能量项的形式表示的。

请注意，稠密和稀疏之间的区别并不等同于直接和间接——事实上，所有四种组合都存在。

• 稀疏+间接：这是最广泛使用的表示，从一组关键点匹配中估计三维几何，从而使用没有几何先验的几何误差。例如Jin等人的研究[12]，monoSLAM[4]， PTAM[16]和ORB-SLAM[20]。
• 稠密+间接：该公式从密集的、正则化的光流场中或与之结合估计三维几何形状，从而将几何误差(偏离流场)与几何先验(流场的平滑度)结合起来，示例包括文献[23]和文献[27]。
• 稠密+直接：该公式采用了光度误差以及几何先验估计稠密或半稠密几何。例如DTAM[21]，其前体[26]，REMODE[22]和LSD-SLAM[5]。
• 稀疏+直接：这就是本文提出的表达。它优化了直接在图像上定义的光度误差，而没有结合几何先验。虽然我们不知道最近有任何使用该公式的工作，但Jin等人已经在2003年[13]提出了一种稀疏和直接的公式。然而，与他们基于扩展卡尔曼滤波器的工作相比，我们的方法使用了非线性优化框架。探索稀疏和直接组合的动机将在下一节中阐述。

此外，还有混合方法，如SVO[9]，使用直接公式进行初始对齐和获取对应关系，然后切换到间接公式进行联合模型优化。

1.1 动机

本文提出的单目视觉里程计的直接和稀疏公式是出于以下考虑。

直接：关键点的主要好处之一是对光照和几何畸变鲁棒。例如自动曝光变化、非线性响应函数(伽马校正/白平衡)、镜头衰减(渐晕)、去马赛克伪影，甚至是由卷帘快门引起的强烈几何扭曲。然而，这种对光照变化的鲁棒性甚至不变性是以丢弃在这些变化中具有潜在价值的信息为代价的。

图1 直接稀疏里程计(DSO)。3D重建和跟踪轨迹，围绕建筑物进行1分40秒的视频骑行(仅限单目视觉里程计)。左下角的插图显示了开始点和结束点的特写，显示了轨迹中的累积漂移。下面一行显示了一些视频帧。  

与此同时，对于介绍中提到的所有用例，数以百万计的设备将(而且已经)装备摄像头，只为计算机视觉算法提供数据，而不是捕捉图像供人类消费。这些摄像头也应该将被设计成提供一个完整的传感器模型，并以最佳服务于算法的方式捕捉数据。例如，自动曝光和伽马校正并不是未知的噪声源，而是能够提供更好的图像数据的特征——这些特征可以被纳入模型中，使获得的数据具有更多的信息。由于直接方法建模完整的图像形成过程到像素灰度，它大大受益于更精确的传感器模型。

直接公式的主要好处之一是，它不需要识别点本身，从而允许更细粒度的几何表示(像素逆深度)。此外，我们可以从所有可用数据(包括边缘和弱灰度变化)中采样，生成更完整的模型，并在纹理稀疏的环境中提供更强的鲁棒性。

稀疏。添加几何先验的主要缺点是引入了几何参数之间的相关性，这使得统计上一致的实时联合优化是不可行的(见图2)。这就是为什么现有的稠密或半稠密方法(a)忽略或粗略近似几何参数之间的相关性(橙色)，和/或几何参数和相机位姿之间的相关性(绿色)，以及(b)对稠密几何部分采用不同的优化方法，如原始对偶公式。

图2 稀疏和稠密海塞矩阵结构。左图为稀疏光束调整的海塞矩阵结构，由于几何-几何块是对角的，因此通过舒尔补可以快速求逆。在右图中，几何先验加入了部分非结构化的几何-几何相关性，这导致矩阵不仅维数大，而且还很难求逆。为了简单起见，我们不显示全局相机的内参。

 

此外，当今先验的表达复杂性是具有局限性的。虽然它们使三维重建更稠密、局部更准确和视觉上更吸引人，但我们发现，先验可能会引入偏差，从而降低而不是提高长期的大规模的场景的精度。请注意，随着从真实数据中学习到的更现实、更公正的先验的引入，这种情况很可能会改变。

1.2 贡献

在本文中，我们提出了一种稀疏和直接的单目视觉里程计方法，一个示例重建如图1所示。据我们所知，这是唯一的完全直接的方法，联合优化所有涉及的模型参数的概率，包括相机位姿，相机内参，和几何参数(逆深度值)。这与混合方法(如SVO[9])形成了对比，后者恢复为联合模型优化的间接公式。

在一个滑动窗口中进行优化，旧相机位姿以及离开相机视野的点被边缘化，灵感来自[17]。与现有方法相比，我们的方法进一步充分利用了光度相机标定，包括镜头衰减、伽马校正和已知曝光时间。这种综合光度标定进一步提高了准确性和鲁棒性。

我们基于CPU的实现在笔记本电脑上实时运行。我们在三个不同的数据集(包括几个小时的视频)上进行了广泛的评估，表明它在鲁棒性和准确性方面优于其他最先进的方法(直接和间接)。通过减少设置(更少的点和活动的关键帧)，它甚至可以以5倍的实时速度运行，同时仍然优于最先进的间接方法。在高且非实时设置(更多的点和活动关键帧)下，它创建的半稠密模型在密度上类似于LSD- SLAM，但更准确。

2 直接稀疏模型

我们的直接稀疏里程计是基于对最近帧的一个窗口的光度误差的持续优化，考虑到图像形成的光度标定模型。与现有的直接方法相比，我们联合优化所有涉及的参数(相机内参、相机外参和逆深度值)，有效地执行了窗口稀疏光束调整的光度等效。我们保留了其他直接方法所采用的几何图形表示方法，例如，三维点在参考坐标系中表示为逆深度(因此有一个自由度)。

记号。在本文中，粗体小写字母($xx$)表示向量，粗体大写字母($HH$)表示矩阵。标量用小写字母($t$)表示，函数(包括图像)用小写字母($I$)表示。相机位姿表示为变换矩阵${TT}_i\in SE(3)$，将一个点从世界坐标系转换到相机坐标系。线性化的姿态增量将被表示为李代数元素$x_i\in \mathfrak{se}(3)$，我们直接将其写成向量$x_i\in {\mathbb{R}}^6$——稍微滥用了符号。我们进一步定义常用的操作符$⊞:\mathfrak{se}(3)\times SE(3)\to SE(3)$，利用左乘表述，即
$${xx}_i⊞{TT}_i:=e^{\widehat{x_i}}\cdot {TT}_i\text{\hspace{1em}(1)}$$

2.1 标定

直接法是对图像建模的综合方法。除了几何相机模型(包含将3D点投射到2D图像的功能)之外，考虑光度相机模型也是有益的，它包含将传感器上的像素接收到的真实能量(辐照度)映射到相应灰度值的函数。请注意，对于间接方法，这是没有什么好处，因此被广泛忽略，因为常见的特征提取器和描述子是不变的(或高度鲁棒)光照变化。

2.1.1 几何相机标定

为了简单起见，我们为众所周知的针孔相机模型制定了我们的方法——在预处理步骤中消除径向失真。而对于广角相机，这确实减少了视野，它允许跨方法的比较，只实现有限的相机模型的选择。在本文中，我们将用${\Pi }_c:{\mathbb{R}}^3\to \Omega$表示投影，用${\Pi }_c^{-1}:\Omega \times \mathbb{R}\to {\mathbb{R}}^3$表示反投影，其中$c$表示相机的内参(对于针孔模型，这是焦距和主点)。请注意，与[2]类似，我们的方法可以扩展到其它(可逆)相机模型，尽管这确实增加了计算需求。

2.1.2 光度相机标定

我们使用[8]中使用的成像模型，它包含非线性响应函数$G:\mathbb{R}\to [0,255]$，以及透镜衰减(渐晕)$V:\Omega \to [0,1]$。图3显示了一个来自TUM monoVO数据集的标定示例。那么，结合的模型可以由下式给出，
$$I_i(xx)=G(t_{i}V(xx)B_i(xx))\text{\hspace{1em}(2)}$$
其中，$B_i$和$I_i$为第$i$帧的辐照度和观察到的像素灰度，$t_i$为曝光时间。该模型的第一步是对每个视频帧进行光度校正，即计算，
$$I_i^{\prime }(xx):=t_i{B}_i(xx)=\frac{G^{-1}(I_i(xx))}{V(xx)}\text{\hspace{1em}(3)}$$
在本文的其余部分中，除另有说明外，$I_i$总是指经过光度校正的图像$I_i^{\prime }$。

图3 光度标定。上图：图1所用相机的逆响应函数G^{-1}和镜头衰减V。下图：包含室内和室外部分的序列的毫秒级曝光t。请注意，它的变化幅度超过500倍，从0.018 ms到10.5 ms。我们不把这些量当作未知的噪声源，而是在光度误差模型中明确地说明它们。

 

2.2 模型表述

我们将参考坐标系$I_i$中$pp\in {\Omega }_i$点在目标坐标系$I_j$中观测到的光度误差定义为像素小邻域上的加权SSD。我们的实验中已经表明，8个像素排列在一个稍微分散的模式(见图4)，在评估所需的计算、对运动模糊的鲁棒性和提供足够信息之间给出了一个很好的权衡。注意，就所包含的信息而言，在如此小的像素邻域上评估SSD类似于为中心像素添加一阶和二阶辐照度导数常数项(除了辐照度常数)。令，
$$E_{p_j}:=\sum _{p\in {\mathcal{N}}_p}{w}_p\Vert (I_j[{pp}^{\prime }]-b_j)-\frac{t_j{e}^{a_j}}{t_i{e}^{a_i}}(I_i[pp]-b_i){\Vert }_{\gamma }\text{\hspace{1em}(4)}$$
其中${\mathcal{N}}_p$是包含在SSD当中的像素集合；$t_i$和$t_j$分别是图像$I_i$和$I_j$的曝光时间；$\Vert \cdot \Vert$是Huber范数。此外，${pp}^{\prime }$表示$(pp,d_p)$的投影点，由下式给出，
$${pp}^{\prime }={\Pi }_c(RR{\Pi }_c^{-1}(pp,d_p)+tt)\text{\hspace{1em}(5)}$$
其中，
$$\left[\begin{array}{cc}RR& tt\\ 0& 1\end{array}\right]:={TT}_j{TT}_i^{-1}\text{\hspace{1em}(6)}$$

图4 残差模式。该模式\mathcal{N}_p用于能量计算。省略右下像素以启用SSE优化处理。注意，由于我们每个点有1个未知值(它的逆深度)，并且不使用正则化，我们需要|N_p|>1，以便在仅优化两帧时，所有模型参数都得到良好约束。图19显示了该模式如何影响跟踪精度的评估。

 

为了使我们的方法能够在没有已知曝光时间的情况下对序列进行操作，我们加入了一个由$e^{-a_i}(I_i-b_i)$给出的附加仿射亮度传递函数。请注意，与大多数以前的公式[6]和[13]相比，标量因子$e^{-a_i}$是对数参数化的。这既防止它变成负的，并避免由乘法(即指数增加)漂移引起的数值问题。

除了使用鲁棒的Huber惩罚外，我们还应用了一个梯度相关的加权$w_p$，
$$w_p:=\frac{c^2}{c^2+||\nabla I_i(pp)|{|}_2^2}\text{\hspace{1em}(7)}$$
降低那些高梯度值像素的权重。该加权函数可以在概率上解释为在投影点位置${pp}^{\prime }$上添加小的独立的几何噪声，并立即边缘化它——近似的小几何误差。总结而言，$E_{p_j}$取决于如下变量：（1）点的逆深度$d_p$；（2）相机内参$c$；（3）相关帧${TT}_i$和${TT}_j$的位姿；（4）它们的亮度转换参数$a_i$、$b_i$、$a_j$和$b_j$。

所有帧和点的全部光度误差由下式给出，
$$E_{photo}:=\sum _{i\in \mathcal{F}}\sum _{p\in P_i}\sum _{j\in obs(p)}{E}_{p_j}\text{\hspace{1em}(8)}$$
其中$i$遍历所有的帧$\mathcal{F}$，$pp$遍历第$i$帧的像素点集合${\mathcal{P}}_i$，$j$遍历能观测到点$pp$的所有帧的集合$obs(pp)$。图5展示了最终的因子图。与经典重投影误差的唯一区别是每个残差附加依赖于主帧的位姿，也就是说，每个残差项依赖于两帧而不是一帧。虽然这增加了非对角条目的位姿块的海塞矩阵，但它不影响稀疏模式后应用舒尔补边缘化的点参数。由此得到的系统可以类比地解为间接公式。请注意，雅可比矩阵相对于两帧的位姿是通过它们的相对位姿的伴随线性相关的。在实践中，当计算海塞矩阵或它的舒尔补时，这个因子可以从求和中提取出来，大大减少了由更多变量依赖引起的额外计算。

图5 直接稀疏模型的因子图。4个关键帧和4个点：KF1中1个点，KF2中2个点，KF3中1个点。每个能量项(公式(4)中定义的)取决于点的主帧(蓝色)、观测到该点的帧(红色)和点的逆深度(黑色)。此外，所有的项取决于全局相机内参向量c，在图中没有显示。

 

如果曝光时间已知，我们进一步添加一个先验使得仿射亮度传递函数为零，
$$E_{prior}:=\sum _{i\in \mathcal{F}}({\lambda }_a{a}_i^2+{\lambda }_b{b}_i^2)\text{\hspace{1em}(9)}$$
如果没有提供光度标定，我们设$t_i=1$和${\lambda }_a={\lambda }_b=0$，在这种情况下，他们需要模拟相机曝光时间的未知变化。作为一个附注，应该提到的是，如果$x_i$和$y_i$都包含噪声测量(参见[7])，那么乘性因子$a^{\ast }=argma{x}_a{\sum }_i(a{x}_i-y_i{)}^2$的ML估计器是有偏差的；导致$a$在无约束情况${\lambda }_a=0$下存在漂移。虽然这通常对估计的位姿影响不大，但如果场景只包含少量的较弱的灰度变化，它可能会引入偏差。

点的维度。在提出的直接模型中，一个点仅由一个参数(参考系中的逆深度)参数化，而在间接模型中有三个未知数。为了理解这种差异的原因，我们首先注意到，在这两种情况下，3D点实际上都是连续的真实世界的3D表面上任意位置的离散样本。不同之处在于曲面上这个2D位置的定义方式。在间接方法中，它隐式定义为在使用的角点响应函数中产生最大值的点(投影到图像中)。这意味着表面和点在表面上的位置都是未知的，需要估计。在我们的直接公式中，点被简单地定义为源像素的光线照射到表面的点，因此只剩下一个未知点。除了减少参数的数量之外，这自然可以实现逆深度参数化，在高斯框架中，它更适合表示基于立体深度估计的不确定性，特别是对于距离较远的点[3]。

一致性。严格地说，提出的直接稀疏模型允许多次使用一些观测值(像素值)，而另一些则完全不使用。这是因为——尽管我们的点选择策略试图通过在空间中均匀分布点来避免这种情况(参见第3.2节)——我们允许点观测重叠，从而依赖于相同的像素值。这种情况尤其发生在纹理很少的场景中，所有的点都必须从纹理图像区域的一个小子集中选择。然而，我们认为这在实践中的影响可以忽略不计，如果需要的话，可以通过删除(或降低权重)使用相同像素值的观测值来避免。

2.3 窗口优化

我们遵循Leutenegger等人[17]的方法，使用高斯-牛顿算法在滑动窗口中优化总误差(公式(8))，这在速度和灵活性之间给出了很好的权衡。为了便于表示，我们将公式(1)中定义的$⊞$操作符扩展到所有优化参数——对于$SE(3)$以外的参数，它表示常规的加法。我们将使用$\zeta \in SE(3{)}^n\times {\mathbb{R}}^m$来表示所有优化的变量，包括相机位姿、仿射亮度参数、逆深度值和相机内参。与[17]中一样，边缘化依赖于$\zeta$中的参数的残差将固定切线空间，在该切线空间中积累有关该参数的任何未来信息(增量更新)。我们将用${\zeta }_0$表示这个切线空间的评估点，用$\zeta =x⊞{\zeta }_0$表示累积的增量更新。图6可视化了不同变量之间的联系。

图6 窗口优化。红色曲线表示参数空间，由$SE(3)$中的非欧几里得相机位姿和其余欧几里得参数组成。蓝色的线对应于\zeta_0周围的切线空间，在这里我们(1)累计x上的二次边缘化先验，(2)计算高斯-牛顿步骤\delta。对于每个参数，一旦该参数成为边缘化项的一部分，切线空间就会被固定。请注意，当我们在我们的表示法中平等地对待所有参数时，对于欧氏参数、切空间和参数空间是重合的。

 

高斯-牛顿优化。我们计算高斯-牛顿系统如下，
$$H=J^{T}WJandb=-J^{T}Wr\text{\hspace{1em}(10)}$$
其中$W\in {\mathbb{R}}^{n\times n}$是包含权重的对角矩阵，$r\in {\mathbb{R}}^n$是堆积的残差向量，$J\in {\mathbb{R}}^{n\times d}$是$r$的雅克比。

注意，对于每个点向能量贡献$|{\mathcal{N}}_p|=8$个残差。为了记号简便，下面我们只考虑一个残差$r_k$，以及雅可比矩阵$J_k$的相关行。在优化期间，以及在边缘化时，残差总是按当前状态估计进行评估，即：
$$\begin{gathered} r_k=r_k(x⊞{\zeta }_0) \\ =(I_j[{pp}^{\prime }(T_i,T_j,d,c)]-b_j)-\frac{t_j{e}^{a_j}}{t_i{e}^{a_i}}(I_i[pp]-b_i)\text{\hspace{1em}(11)} \end{gathered}$$
其中$(T_i,T_j,d,c,a_i,a_j,b_i,b_j):=x⊞{\zeta }_0$是当前残差决定的状态变量。雅可比矩阵$J_k$的值是关于对$x$的附加增量的，
$$J_k=\frac{\partial r_k((\delta +x)⊞{\zeta }_0)}{\partial \delta }\text{\hspace{1em}(12)}$$
它可以分解为，
$$J_k=[\underset{J_I}{\underbrace{\frac{\partial I_j}{\partial {pp}^{\prime }}}}\underset{J_{geo}}{\underbrace{\frac{\partial {pp}^{\prime }((\delta +x)⊞{\zeta }_0)}{\partial {\delta }_{geo}}}},\underset{J_{photo}}{\underbrace{\frac{\partial r_k((\delta +x)⊞{\zeta }_0)}{\partial {\delta }_{photo}}}}]\text{\hspace{1em}(13)}$$
其中${\delta }_{geo}$表示“几何”参数($T_i,T_j,d,c$)，${\delta }_{photo}$表示“光度”参数$(a_i,a_j,b_i,b_j)$。我们采用两种近似方法，如下所述。

首先，$J_{photo}$和$J_{geo}$都是在$x=0$处计算而来的。这种技术被称为“首次估计雅可比矩阵”，用于保持系统的一致性和防止虚假信息的积累。特别是，在能量中存在非线性零空间(在我们的公式中绝对位姿和比例)，在不同的评估点周围添加线性化可以消除这些，从而慢慢地破坏系统。在实践中，这种近似是非常好的，因为$J_{photo}$和$J_{geo}$相比于增量$x$的大小是平滑的。相比之下，$J_I$就不那么平滑了，但是不影响零空间。因此，取$x$的当前值，即与残差$r_k$在同一点。我们利用中心差分计算图像在整数位置的导数，然后用双线性插值。

其次，假设$J_{geo}$对于属于同一点的所有残差都是相同的，并且只对中心像素进行计算。同样，这个近似在实践中是非常好的。虽然它大大减少了所需的计算，但我们没有观察到它对使用的任何数据集的准确性有显著影响。

从产生的线性系统中，增量被计算为$\delta =H^{-1}b$，并添加到当前状态，
$$x^{new}\leftarrow \delta +x\text{\hspace{1em}(14)}$$
注意，由于第一估计雅可比矩阵近似，一个乘式公式(用$\delta ⊞(x⊞{\zeta }_0)$替换公式(12)中的$(\delta +x)⊞{\zeta }_0$)得到完全相同的雅可比矩阵，因此乘式更新步骤$x^{new}\leftarrow log(\delta ⊞e^x)$同样有效。

在每个更新步骤之后，我们使用${\zeta }_0^{new}\leftarrow x⊞{\zeta }_o$和$x\leftarrow 0$更新不属于边缘化项的所有变量的${\zeta }_0$。实际上，这包括所有深度值以及最新关键帧的位姿。每次添加一个新的关键帧时，我们执行多达6次高斯-牛顿实验，如果$\delta$足够小，就会提前中断。我们发现——因为我们从不远离最小值——不需要Levenberg-Marquardt阻尼(减慢收敛速度)。

边缘化。当活动变量集变得太大时，使用舒尔补边缘化旧的变量。与[17]相似，我们去掉所有会影响$H$的稀疏模式的残差项。当边缘化$i$帧时，我们首先边缘化${\mathcal{P}}_i$中的所有点，以及在最后两个关键帧中没有被观察到的点。关键帧$i$中活动点的剩余观测值从系统中删除。

边缘化的进程如下。设$E^{\prime }$表示能量中包含所有依赖于被边缘化的状态变量的残差部分。我们首先在当前状态估计$\zeta =x⊞{\zeta }_0$附近计算一个高斯-牛顿近似$E^{\prime }$。这给出了，
$$\begin{gathered} E^{\prime }(x⊞{\zeta }_0) \\ \approx 2(x-x_0{)}^{T}b+(x-x_0{)}^{T}H(x-x_0)+c \\ =2x^T\underset{=:b^{\prime }}{\underbrace{(b-H{x}_0)}}+x^{T}Hx+\underset{=:c^{\prime }}{\underbrace{(c+x_0^{T}H{x}_0-x_0^{T}b)}}\text{\hspace{1em}(15)} \end{gathered}$$
其中$x_0$表示$x$的当前值($r$的评估值)。常数$c$和$c^{\prime }$可降，$H$和$b$的定义为(10)、(11)、(12)和(13)。这是一个关于$x$的二次函数，我们可以应用舒尔补边缘化变量的子集。写成线性方程组，它变成，
$$\left[\begin{array}{cc}{HH}_{\alpha \alpha }& {HH}_{\alpha \beta }\\ {HH}_{\beta \alpha }& {HH}_{\beta \beta }\end{array}\right]\left[\begin{array}{c}{xx}_{\alpha }\\ {xx}_{\beta }\end{array}\right]=\left[\begin{array}{c}{bb}_{\alpha }^{\prime }\\ {bb}_{\beta }^{\prime }\end{array}\right]\text{\hspace{1em}(16)}$$
其中$\beta$表示我们想要边缘化的变量块，$\alpha$表示我们想要保留的变量块。应用舒尔补，产生$\widehat{H_{\alpha \alpha }}{x}_{\alpha }=\widehat{b_{\alpha }^{\prime }}$，
$$\widehat{H_{\alpha \alpha }}=H_{\alpha \alpha }-H_{\alpha \beta }{H}_{\beta \beta }^{-1}{H}_{\beta \alpha }\text{\hspace{1em}(17)}$$
$$\widehat{b_{\alpha }^{\prime }}=b_{\alpha }^{\prime }-H_{\alpha \beta }{H}_{\beta \beta }^{-1}{b}_{\beta }^{\prime }\text{\hspace{1em}(18)}$$

因此，$x_{\alpha }$上的剩余能量可以写成，
$$E^{\prime }(x_{\alpha }⊞({\zeta }_0{)}_{\alpha })=2x_{\alpha }^T\widehat{b_{\alpha }^{\prime }}+x_{\alpha }^T\widehat{H_{\alpha \alpha }}{x}_{\alpha }\text{\hspace{1em}(19)}$$
这是一个关于$x$的二次函数，可以在后续的所有优化和边缘化操作中简单地添加到全光度误差$E_{photo}$中，取代相应的非线性项。注意，这要求在所有后续优化和边缘化步骤中，对于$E^{\prime }$中出现的所有变量，${\zeta }_0$的切线空间保持不变。