啥是DTAM:DenseTracking&Mapping

啥是DTAM:DenseTracking&Mapping

主要论文出处:

DTAM: Dense Tracking and Mapping in Real-Timewww.robots.ox.ac.uk


0. Overview&Preliminaries

论文提出了使用单个RGB相机进行实时稠密(Dense)定位和跟踪的方法。经过之前两篇文章的解读(啥是KinectFusion啥是DynamicFusion),可以了解到稠密指的是区别基于匹配特征点的方法,而是采用所谓 whole image registration全像素匹配。

首先对后文将会用到的符号作下约定:

c :相机; w :世界空间,则有 :

表示从相机空间到世界空间的变换,其由旋转矩阵和位移构成。则相机空间下的点转换到世界空间即:

对于相机空间下点 x_{c} = (x, y, z) ,定义其透视投影(perspective projection)为:


文中的描述的dense model由一些交叠的关键帧构成,关键帧构成如下:

每一个关键帧 r 对应位姿转换 T_{rw} ,还包括一个反向深度(inverse depth)映射:

对应的RGB图象:

这里的 \Omega 即图象区域;对于像素 u ,可获得其inverse depth:

通过像素和反向深度,经过反投影,获得像素对应的三维坐标:

这里的 K 是相机内参矩阵(可预先通过相机校准获得,在本文算法中作为已知值);而其中的 \dot{u} = (u, v, 1) 即变为齐次坐标。


关于inverse depth,即深度值的倒数

设某点的当前深度值估计为 \rho ,则其inverse depth为 d = 1/\rho

在基于视觉的深度估算中,经常使用inverse depth描述深度,此类算法中多会用到扩展卡尔曼滤波,使用inverse depth可以满足误差高斯分布,有利于扩展卡尔曼滤波发挥效力;此外就是允许描述远点(离相机光心远),举个极端例子如太阳这样的无穷远距离用inverse depth表示的话就是0了。关于inverse dpeth及其参数化描述论文:

Inverse Depth Parametrization for Monocular SLAMwww.doc.ic.ac.uk

1. Dense Mapping

论文构造了一个全局优化能量方程,用来估算关键帧的反向深度映射 ;该能量方程由光测度误差项(photometric error data term )与空间正则项(robust spatial regularisation term)之和组成。

如上图所示,其中

r 表示当前关键帧;

C_{r} 定义的是一个如上图中方格体素区域的光照投影损耗体(projective photometric cost volume)在双摄像机跟踪中也叫视差空间(disparity space) ;

I_{r} 关键帧的RGB图象,其上的像素点 u_{r}

T_{rw} 关键帧变换矩阵;

[\xi_{min}, \xi_{max}] 则表示inverse depth的取值范围;

m\in \mathcal{I}(r) 关键帧附近的连续帧及对应的变换矩阵 T_{mw} ,用于迭代计算;

以下两个关键概念

C_{r}(u) :关键帧RGB图象上每个像素点 对应 中的一纵行,图中红色所示,记为 ,存储了该像素点的光度误差均值 ;

C_{r}(u, d) :如图中红色部分所示,对于同一像素点 ,其inverse depth由最小值取到最大值,每个inverse depth都对应了一个光度误差均值,记为 ;

那么 C_{r}(u, d) 怎么计算呢?

先看给公式吧:

结合之前的定义来解释下公式(2)和(3)即 将cost volume中的点投影到关键帧附近 m 帧上,然后和关键帧 r 上与之对应像素点的值作差,取 L1 norm,将所有结果相加。那么显然(2)(3)的值越小,其对应的 d 就是该像素点要的inverse depth值。

若对光流有所了解,会觉得(3)似曾相识,的确(3)是以类似光流假设中以光照恒定假设为前提的。关于光流可以参看

南瓜番茄:啥是TV-L1 OpticalFlowzhuanlan.zhihu.com图标


(2)和(3)的表现如下图所示

图象中取3个像素点 a, b, c,每条细彩色曲线表示在帧 m \in \mathcal{I}(r) 上该像素点随着inverse depth值在 [\xi_{min}, \xi_{max}] 取值对应的式(3) 即 \rho(u, d) 的变化,粗红线为均值即 C_{r}(u, d) 变化;虚线表示 C_{r}(u, d) 取到的最小值位置。

先总体观察三幅曲线图的相似性,看到单个 \rho(u, d) 即单个曲线会有许多局部最小值,而总体的均值变化 C_{r}(u, d) 会趋向于一个最小值,如图中虚线所示位置。

再看看三个像素点曲线图的不同之处,像素 b 和 c 都有明显的最小值,特别是 b,而像素 a 的整体趋向较为平缓,而局部最小值变化较大说明对 \rho(u, d) 的估算不稳定。

原因呢看回RGB图,像素 a 处于无变化处,梯度趋于0,文中用textureless来描述此类区域;像素 b ,与 a 正好相反,处边界区域,文中称strongly textured区域;像素 c 处于一组交替变化的像素区 linear repeating texture 。可见算法对梯度较大的区域更为友好。


再解释下为啥用多帧均值,而不是单独两个视角估算。看式(3)是两个视角下的估算,即只有两张RGB图( 中取一帧和关键帧 ),若第二张图出现了阻塞(例如有遮挡)的情况明显就无法正确估算视差信息了,而使用多个视角下图象,某种程度上增加了信噪比(Signal-to-noise ratio),阻塞区域被视为离群值,与此同时使用 L1 norm对离群不敏感的特性,可以减少那些被阻塞的帧信息的干扰,而保留未被阻塞的帧中的信息。


像素 u 的inverse depth可以通过 arg\ min_{d}C(u, d) 获得,然而我们已经了解到,对于那些textureless区域的像素点估算并不准确,本该平滑的区域却坑坑洼洼,如下图中机箱底、显示器屏幕...:

左到右本别是2,10,30帧(即 \mathcal{I}(r) 数量)迭代获取的inverse depth结果,可以看到机箱底部,屏幕的marker,都没有被准确估算到。怎么办呢?作者构建了这样一个能量泛函:就如本小节开头所述,它由数据项和正则项之和构成,数据项就是依据光度恒定误差构建的,而正则项用于惩罚本该是空间平滑区域的分离情况(图中坑坑洼洼的情况),保证inverse depth的空间平滑性。增加平滑正则项后,效果对比如下(最后一张图):

这个正则项如果构建呢?


Regularised Cost

使用带有权重的huber norm正则项:

其中 g(u) 是权重, || \nabla \xi(u) ||_{\epsilon} 是对inverse depth map的梯度作huber norm。

这个huber norm,之前在 啥是DynamicFusion 有提到过,其是一个鲁棒方法,实质上混合L1 norm 和 L2 norm,形如本文所示:

其用作正则项与其他norm有啥区别呢,来看图:

http://www.stephanmandt.com/papers/ECML_2016.pdf

结构最简单的鲁棒回归应该就是L1 norm了,但是和 L2 相比显然平滑性太差了,把所有的小的误差(residual/ misfit)都当作大误差给舍弃了,而huber norm兼具了L1鲁棒性和L2的平滑性。显然这里的阈值 \epsilon 的选择比较关键,作者经过实验给出阈值约等于 1.0e-4 时,效果比较理想。

了解了huber norm那么就可以解释 || \nabla \xi(u) ||_{\epsilon} 的用意了,即当像素 u 所对应inverse depth的变化率比较小时,说明可能是空间连续的,只不过是因为其处于textureless区域而出现了梯度波动,那么就使用L2 norm保证其空间平滑;当变化率陡增时,则判断其可能发生了阻塞,使用鲁棒性较高的L1,以保证最终结果少受离群值影响。

再来看权重:

可以看到其从texture角度来考虑,梯度越高的像素区域其权重越小,即原本就是边界区域的像素,其正则项作用越小,以保留其本身的边界(非连续性)。那么最终的能量泛函写作:

撇去三维体积空间的引入,该公式和 啥是TV-L1 OpticalFlow 讨论的带有TV-L1正则项的改进后的Horn-Schunck稠密光流很相似:

那么或许可以借用其TV-L1解法去解式(6) 。事实上作者也确实是这么做的, 将(6)作如下变换,写作:

其中 \alpha 是新引入的一个辅助变量,以配合基于对偶的TV-L1求解法。


Discretisation of the Cost Volume and Solution

前面提到的Cost Volume显然是三个维度的,用 M \times N \times S 表示好了, M \times N 是图象维度(分辨率),而 S 则是在 [\xi_{min}, \xi_{max}] 之间进行线性采样的次数,以表示inverse depth。

下面推导中,分别用 d, a, g 分别表示(7)中的 \xi, \alpha, g 的向量形式。

对于(7)的前半部分,根据 Legendre-Fenchel 变换:

以及 Fenchel 不等式:

推出:

其中 d 左乘的矩阵 A 表示计算其梯度, G = diag(g) 则用于乘以对应的权重值, q 是引入的对偶量。其中 \delta_{q}(q) 某种意义是个指示函数: \delta_{q}(q) \begin{equation} \left\{ \begin{array}{lr} 0, & ||q||_{1} \leq 1\\ \infty, & otherwise. \\ \end{array} \right. \end{equation}

将(8)带入到原始问题中,那么得到的对偶问题变为:

固定辅助变量 \alpha ,当 \partial_{d,q}(E(d,a,q)) = 0 时,得最后解。首先对 q

再对 d

这里依据散度定理(散度定理-维基百科啥是TV-L1 OpticalFlow),有

其中 A^T = - div


首先固定 d ,求辅助变量 \alpha_{u} = \alpha(u) ,对于仍是非凸的问题,使用point-wise search暴力解决咯:


而对于 d 的求解,进行如下迭代:

初值 n = 0;\ q^0 = 0;\ d^0_u=a^0_u = arg\ min_{a_{u}}C(u, a_{u})

1. 固定 \alpha^n ,分别对 q, d\partial_q = 0, \partial_d = 0 按照 半隐性梯度上升和下降迭代:

其中的 \sigma 为步长;再整理如下:

2. 再固定 d^{n+1} ,带到(13)中依然是point-wise exhaustive search

3. 当 \theta^n > \theta_{end} ,按照 \theta^{n + 1} = \theta^{n}(1 - \beta n), n \leftarrow n + 1 ,继续迭代;否则结束迭代。

2. Dense Tracking

最近事太多,之后有时间再捋吧.....


Reference&Acknowledgement

Inverse Depth Parametrization for Monocular SLAM

stephanmandt.com/papers

南瓜番茄:啥是TV-L1 OpticalFlow

Legendre and Legendre-Fenchel transforms

编辑于 2018-09-19

文章被以下专栏收录