DSO windowed optimization 代码 (2)

这里写一点代码与公式的对应关系,解决 Schur Complement 更新逆深度的细节问题。

\[\begin{align} \begin{bmatrix} H_{\rho\rho} & H_{\rho X} \\ H_{X\rho} & H_{XX} \end{bmatrix} \begin{bmatrix} \delta \rho \\ \delta X \end{bmatrix} &= - \begin{bmatrix} J_{\rho}^T r \\ J_X^T r \end{bmatrix} \notag \\
\begin{bmatrix} H_{\rho\rho} & H_{\rho X} \\ 0 & H_{XX} - H_{X\rho} H_{\rho\rho}^{-1} H_{\rho X} \end{bmatrix} \begin{bmatrix} \delta \rho \\ \delta X \end{bmatrix} &= - \begin{bmatrix} J_{\rho}^T r \\ J_X^T r - H_{X\rho} H_{\rho\rho}^{-1} J_{\rho}^T r \end{bmatrix} \notag \end{align}\]

在得到 \(\delta X\) 之后可以计算 \(\delta \rho\),这个计算过程是每一个逆深度分别计算,因为矩阵实在是很大,直接计算无法求逆。

\[\begin{align} H_{\rho\rho}\delta \rho + H_{\rho X} \delta X &= -J_\rho^Tr \notag \H_{\rho\rho}\delta \rho &= - J_\rho^Tr - H_{\rho X} \delta X \notag \\left({\partial r \over \partial \rho}\right)^T {\partial r \over \partial \rho} \delta \rho &= -\left( \left({\partial r \over \partial \rho}\right)^T r + \left({\partial r \over \partial \rho}\right)^T {\partial r \over \partial X}\delta X \right) \notag \end{align}\]

式子中各个矩阵、向量的维度如下:

\(r\) : \(N \times 1\)

\(\rho, \delta \rho\) : \(M \times 1\)

\(\delta X\) : \(68 \times 1\)

\({\partial r \over \partial \rho}\) : \(N \times M\)

\({\partial r \over \partial X}\) : \(N \times 68\)

仅考虑 \(\delta\rho\) 中的一个逆深度更新量。

下面对应代码中的 EnergyFunctional::resubstituteFPt。

由于每一个 residual 只对一个逆深度求导不为 0。 \({\partial r \over \partial \rho}\) 一共有 \(N\) 行,每一行都有且仅有一个不为 0 的数字。这样的结果造成 \(\left({\partial r \over \partial \rho}\right)^T {\partial r \over \partial \rho}\) 是一个对角阵,于是可以分别计算每一个逆深度的更新量。

\(\left({\partial r \over \partial \rho}\right)^T {\partial r \over \partial \rho}\) 对角线上的元素是 \(\sum \left({\partial r_{21} \over \partial \rho_1}\right)^T {\partial r_{21} \over \partial \rho_1}\),这是与该逆深度对应的所有 residual 对它导数的和,对应EFPoint::Hdd_accAF+EFPoint::Hdd_accLF。对角线上元素的倒数是EFPoint::HdiF

对应代码:

p->data->step = -b * p->HdiF;

\(\left({\partial r \over \partial \rho}\right)^T r\)对应EFPoint::bdSum_F,依旧是每个逆深度对应的 residual 相关参数求和。

对应代码:

float b = p->bdSumF;

\[\begin{align} \left({\partial r \over \partial \rho}\right)^T {\partial r \over \partial X}\delta X &= \left({\partial r \over \partial \rho}\right)^T \begin{bmatrix} {\partial r \over \partial C} &{\partial r \over \partial X_F}\end{bmatrix} \begin{bmatrix} \delta C \\ \delta X_F \end{bmatrix} \notag \&= \left({\partial r \over \partial \rho}\right)^T \left( {\partial r \over \partial C} \delta C + {\partial r \over \partial X_F} \delta X_F \right)\notag \&= \left({\partial r \over \partial \rho}\right)^T {\partial r \over \partial C} \delta C + \left({\partial r \over \partial \rho}\right)^T {\partial r \over \partial X_F} \delta X_F \notag \end{align}\]

\(X_F\) 表示所有帧的参数,共 64 个,其中包含 se(3) 也包括 affLight。

\(\left({\partial r \over \partial \rho}\right)^T {\partial r \over \partial C} \delta C\) 对应代码:

b += xc.dot(p->Hcd_accAF + p->Hcd_accLF);

(我觉得Engel博士写了错误的代码,上面的代码是按照我推出的公式写的。)

\(\left({\partial r \over \partial \rho}\right)^T {\partial r \over \partial X_F} \delta X_F\) 对应代码:

b += xAd[r->hostIDX * nFrames + r->targetIDX] * r->JpJdF;

其中r->JpJdF表示\(\left( {r_{21} \over \partial X_{21}} \right)^T {\partial r_{21} \over \partial \rho_1}\),转置一下不是很大关系,因为这个就是一个scalar。

xAd[r->hostIDX * nFrames + r->targetIDX]在 EnergyFunctional::resubstituteF_MT 中构造,对应的公式是

\[\delta X_1^T \left( {\partial X_{21}\over \partial X_1}\right)^T + \delta X_2^T \left( {\partial X_{21}\over \partial X_2}\right)^T\]

所以,最终\(\left({\partial r \over \partial \rho}\right)^T {\partial r \over \partial X_F} \delta X_F\)对于每一个逆深度而言是(做了一次转置)

\[\sum \left( \delta X_1^T \left( {\partial X_{21}\over \partial X_1}\right)^T + \delta X_2^T \left( {\partial X_{21}\over \partial X_2}\right)^T \right)^T \left( {r_{21} \over \partial X_{21}} \right)^T {\partial r_{21} \over \partial \rho_1}\]

回溯到 EnergyFunctional::resubstituteF_MT 函数中,确定一下adHostFadTargetF的所代表的意义。

adHostF[h->idx + nFrames * t->idx]对应\(\left( {\partial X_{21}\over \partial X_1}\right)^T\),按照这个下标可以看做是第t行第h列,也就是adHost[2, 1](用我喜欢的 21 表示 target host)。

adTargetF[h->idx + nFrames * t->idx]对应\(\left( {\partial X_{21}\over \partial X_2}\right)^T\)。

原文地址:https://www.cnblogs.com/JingeTU/p/8395046.html

时间: 2024-10-18 17:31:54

DSO windowed optimization 代码 (2)的相关文章

DSO windowed optimization 代码

这里不想解释怎么 marginalize,什么是 First-Estimate Jacobian.这里只想看看代码,看看Hessian矩阵是怎么构造出来的. 代码中用于存储窗口优化过程Hessian矩阵使用的导数存储在结构体RawResidualJacobian中. https://github.com/JakobEngel/dso/blob/master/src/OptimizationBackend/RawResidualJacobian.h#L32 struct RawResidualJa

DSO 中的Windowed Optimization

该博客内容发表在泡泡机器人公众号上,请尊重泡泡机器人公众号的版权声明 DSO中除了完善直接法估计位姿的误差模型外(加入了仿射亮度变换,光度标定,depth优化),另一个核心就是像okvis一样使用sliding window来优化位姿,Engel也专门用了一节来介绍它.sliding window 就像c++中的队列,队首进来一个新人,队尾出去一个老人,它更像王朝中武将的新老交替,老将解甲归田,新人受window大王的重用,然而安抚老将不得当,会使得SLAM王朝土崩瓦解.对于初次接触slidin

DSO Marginalization

每当有新的关键帧产生时,都会在 FullSystem::makeKeyFrame() 函数的这里调用函数 EnergyFunctional::insertFrame() 为窗口优化添加新的关键帧,在这里将HM和bM resize,为新帧创造空间,新帧对应的位置在末尾. 有新的帧加入,就有旧的帧离开.从窗口中剔除旧的关键帧就是边缘化(marginalization)过程.边缘化过程有两步,一步是点的边缘化,一步是帧的边缘化. 整个流程在 FullSystem::makeKeyFrame() 当中,

Paper Reading: Stereo DSO

开篇第一篇就写一个paper reading吧,用markdown+vim写东西切换中英文挺麻烦的,有些就偷懒都用英文写了. Stereo DSO: Large-Scale Direct Sparse Visual Odometry with Stereo Cameras Abstract Optimization objectives: intrinsic/extrinsic parameters of all keyframes all selected pixels' depth Inte

CSDN日报20170323——《你首先是一个人,然后你才是程序员》

[程序人生]你首先是一个人,然后你才是程序员 作者:左潇龙 程序员这个职业在外面,一直被过于神化,又或者过于丑化.但其实程序员这个职业,与大多数职业一样,并没有什么特别的地方. 唯一不同的是,厨师是靠自己的厨艺吃饭,而程序员是靠自己的技术吃饭,仅此而已罢了. 作为一个厨师,如果想名扬天下,做出一番事业,光靠提高自己的厨艺肯定不行,毕竟做饭好吃的人多了去了,你又算老几? 同样的,程序员也是一样,技术牛逼的人多了去了,你又能在世界上排第几?你又能在中国排第几?你又能在你的城市排第几? 说这个,并不是

DX11学习一-初始化D3D

初始化D3D我们需要以下几步步 一.创建driver type,feature level driver type分为3种类-hardware,warp以及reference,代表GPU,CPU模拟和软件模拟,性能依次下降 feature level代表dx的等级,从11到10依次下调 按driver type以及feature level,从高到低遍历,找到合适的跳出. //Specifying the driver type and feature levels D3D_DRIVER_TYP

OD: GS Bypasing via SEH

通过 SEH 绕过 GS 保护 GS 机制没对 SEH 提供保护,所以可心通过攻击异常来绕过 GS. 实验环境为: VMware : Windows 2000 sp4, 此版本无 SafeSEH 的影响 Visual Studio 2005 Project Properties : Release, Disable Optimization 代码如下: 1 #include <string.h> 2 char shellcode[]= 3 "\xFC\x68\x6A\x0A\x38\

vins-mono的边缘化分析

##marg 基础   摘自贺一家的博客 在我们这个工科领域,它来源于概率论中的边际分布(marginal distribution).如从联合分布p(x,y)去掉y得到p(x),也就是说从一系列随机变量的分布中获得这些变量子集的概率分布.回忆了这个概率论中的概念以后,让我们转到SLAM的Bundle Adjustment上,随着时间的推移,路标特征点(landmark)和相机的位姿pose越来越多,BA的计算量随着变量的增加而增加,即使BA的H矩阵是稀疏的,也吃不消.因此,我们要限制优化变量的

错误和问题解决的成本

问题描写叙述 错误 数据收集 根本原因 版本号   组件:数据修复           在一个实际成本组织中,(平均,先进先出,后进先出) 一个或更 多的下面情况可能发生: 1.导航到物料成本历史表单上的数量信息,与现有量表单的数量不匹配的记录 2. 一些物料前期已计成本的数量与前面的事务处理历史表单的数量不匹配 3. 全部的库存值报表与事务处理值报表不匹配 4. 存货层次成本更新表单的总数量与现有量数量表单不匹配(只在先进先出/后进先出) 5.这些症状的不论什么一个意味着 MMT-CQL不匹配