最近对CUDA比较感兴趣,看了一下这篇论文.
Abstract
GPUs牛逼.
我们的工作首先回顾了非极大值抑制(non-maxima suppression的问题, 特别是在GPUs上.
然后提出了一个选择局部响应最大的特征检测, 强制了空间特征分布, 同时同步检测特征.
我们的第二个贡献介绍了一个加强的FAST特征检测, 他应用了之前提到的非极大值抑制方法.
我们将我们的方法和其他CPU和GPU版本的比较, 我们的总是比他们牛逼.
1. Introduction
A. Motivation
现成的(off-the-shelf)相机都可以100fps, 但是算法一般都不做到.
B. Related Work
- 特征检测
特征检测没有重大的改变, 一般都用Harris, Shi-Tomasi, FAST. Harris和Shi-Tomasi对于边缘没有那么敏感, 所以被大范围的用作角点检测器. ORB, 作为FAST的扩展也存在在VIO的世界中. 毫无疑问, FAST代表最快的特征检测器.
以我们所知, 最快的CPU应用是KFAST, 比原始的x5.
在GPU端, 我们意识到两个应用在OpenCV和ArrayFire中. 两者都用了查找表来加速决策过程: 后者用了一个64KB的查找表, 前者用了8KB的查找表. 虽然两个解决方案都提供了快速的角点检测, 但是都不能保证空间特征分布.
- 非极大值抑制
它可以被认为是局部最大值搜索在一个候选的Moore neighborhood(摩尔社区). 每个 像素响应的摩尔社区是它的正方形(长度是2n+1).
提出的算法的复杂度是基于比较的数量的, 这个算法需要((2n+1)^2). [14] 提出了按照螺旋顺序的方法. 理论上, 比较大数量没有改变, 但是实际的数量却暴跌, 因为大多数的候选会在一个更小的3x3邻近中被抑制, 而只有小部分保留.
- 特征跟踪
特征跟踪可以被分为三类: 特征匹配, 滤波跟踪, 差分跟踪.
特征匹配需要在每帧上进行特征提取, 然后是匹配. 但是, 特征检测器的重复性会不利的影响其鲁棒性.
[18] EKF使用了滤波, 把特征的位置包含在状态量中, 然后用连续的预测和更新步骤来跟着特征.
第三种差分方法目标是直接使用像素强度, 最小化光度误差的变化. Lucas-Kanade tracker[19, 20, 21]. 因为它直接在像素强度patch上操作, GPU的版本来的更早[24].
C. 贡献
我们的工作介绍了牛逼的非极大值抑制, 也利用了low-level的GPU命令, 由GPU-optimized应用完成.
我们的方法组合了特征检测和非极大值抑制, 保证一致的特征分布.
我们将我们的前段和最先进的VIO后端组合.
2. Methodology
A. Preliminaries on parallelization
NVIDIA GPU是围绕可扩展array of 多线程streaming multiprocessors(SMs)构建的.
引入了阶级的计算单元: 线程, warps, 线程blocks, thread grid. 每个warp有32个threads.
一个warp里的每个线程in a lock-step basis来运行相同的指令. NVIDIA把这个执行模型叫做Single Instruction Multiple Threads(SIMT). 它还需要一个warp的if/else的divergence导致serialized的执行.
随着NVIDIA Kepler GPU microarchitecture, 在一个warp里的线程可以读互相的register with specific instructions.
我们的工作集中在这些warp-level的primitives(原函数), 高效的交流机制for同一个warp里的线程之间的共享数据.
在上一个GPU的时代中, 线程需要趋向使用common memory (通常是shrared memory)来进行分享数据, 这个会导致很大的延迟. 随着Kepler architecture的引入, 可以在warp之间先交互, 然后只在更高抽象的执行使用更慢的memory.
B. Feature Detector Overview
GPU对于特征检测特别适合, 因为可以被认为是线性沟通pattern里的模板运行. 在stencil operation, 每个计算单元需要输入元素(e.g 像素)和它近邻. 所以, 图像可以高效的被分割在CUDA核中.
对于特征提取, 输入图像首先被降采样获得图像金字塔, 对于每个图像精度, 两个函数会在每个像素被衡量:
- a coarse corner response function (CCRF)
- a corner response function (CRF)
CCRF是快速的计算来敏捷的派出绝大部分的候选点, 所以更慢的CRF函数只要接受通过检测的.
图像的均匀的特征分布会提升VIO的稳定性, 为了满足这个要求, [22, 23] 引入了2D grid cell: 图像被分割为固定大小的长方形, 在每个cell里, 只会选择一个特征: 也就是CRF分最高的.
C. Non-maxima suppresion with CUDA
非极大值抑制也可以被认为是一种削减的操作.
我们的方法把角点的响应图分成规则的cell grid, 在grid的金字塔第一层, cell使用32的倍数, i.e. 32w, 因为NVIDIA GPU硬件, 一个warp有32个线程. cell的一个直线可以被分为有32个元素的cell线段.
我们把cell的高度限制在(2^{l-1}), 这里 (l) 是金字塔的层级.
在一个cell线段中, 一个线程被分配来处理一个像素的响应. i.e. 一个warp处理一整个线段(32线程对应32个像素响应). 但是一个warp可以处理多个线, 一个block里的多个warps合作处理一个cell的连续线.
因为角点响应图被使用float型 in a pitched memory layout保存, 水平的cell边界完美的和L1-cache线边界重合, 这个最大化了内存总线的利用.
为了演示的方便, 一个1:1的warp-to-cell映射被用在32x32的cell. 当一个warp读了cell的第一行的时候, warp里的每一个线程处理了一个像素response. 在下一个步骤, 整个warp开始邻域抑制: 基于[14]旋转, 每个线程验证响应是否在Moore领域里有最大值.
一旦邻域验证做完了, 一些线程可能抑制了它的响应. 但是, 目前未知, 没有写的操作, 每个线程保存它的状态(响应分, x-y位置)在寄存器中.
warp会继续下一个线, 然后重复之前的步骤.
一旦处理完cell所有的线之后, 每个线程有它最大值和对应的2D位置. 但是, 因为32个线程也只是处理了独立的列, 最大值只是column-wise的. 所以warp需要做warp-level reduction来获得cell-wise maximum. 他们减少最大值和位置 to the first thread(thread 0)使用warp-level shuffle down reduction[32]. 线程0最终会把结果写在global memory.
为了加速reduction, 多个warps处理一个cell, 所以, 在warp-level reduction, 最大值已经在shared memory里减少了.
在pyramidical特征检测中, 我们只维护一个grid. 在level 0, 使用了上述的算法. 在第一点的金字塔层, 我们虚拟的scale了cell大小, cell就变成了(left(frac{32 cdot w}{k}, frac{2^{l-1} h}{k} ight)) . 在cell宽度小于32时, 一个warp会处理多个线.
回收看算法1, 我们的方法组合了regular neighborhood suppresion (NMS)和cell最大值选择(NMS-C)在一个步骤里.
D. FAST Feature Detector
fast特征的潜在逻辑很简单: 每个像素位置, 我们做一个segment测试, 我们在Bresenham圈上比较像素强度. 这个圈给了我们16个像素位置.
如果每个线程运行SIMT, 这个比较 in if/else会执行不同的代码. 因为所有的线程会运行同一个指令, 有些线程会inactive在if分支, 有些会inactive在else分支. 这个叫做code divergence, 会减少并行的throughput(通量). 但是它可以用一个不同的方法解决: 一个查找表. 如上图.
我们的方法存储了16个比较结果在bit array. 所有可能的16bit向量会被预计算: a bit (b_x) 是 "1"如果在圈上的像素亮度是更亮/更暗, 是"0"如果是类似的. 因为结果是二进制的, 所以结果可以存储在(2^{16}) bits, i.e. 8 KB.
文献区分了三种计算角点分数的办法:
- sum of absolute differences 在 Bresenham circle (SAD-B).
- sum of absolute differences on the continues arc (SAD-A).
- 最大阈值((epsilon)) for 点还被考虑为角点(MT).
我们的方法压缩了每个16-bit到一个bit, 得到一个8KB查找表.
E. Lucas-Kanade Feature Tracker
我们的方案应用了金字塔的近似同步逆成分的LK算法作为特征跟踪. LK[19]算法最小化一个长方形patch在template和新图的光度误差, 通过使用了一个warping function on the image coordinates. 这个逆成分算法是一个扩展来提升每个迭代的计算复杂度, 通过允许预计算Hessian阵和在每个迭代中复用.
这个同时逆成分LK加入了仿射光度变化的估计. 但是因为Hessian变成了外观估计的函数, 它就不能被预计算了, 这就比原版的LK慢. 所以, 我们用了近似的版本, 假设外观的参数不会重大的变化, Hessian可以用初始估计计算.[21]
我们用平移模型(t) with 仿射光度变化估计(lambda). 完成的参数是 (q=[t, lambda]^T = [t_x, t_y, alpha, eta]^T).
(oldsymbol{W}(oldsymbol{x}, oldsymbol{t})=left(egin{array}{l} x+t_{x} \ y+t_{y} end{array} ight))
每个特征的光度误差是:
这里(T(x) 和 I(x))代表原图和现在的图在位置 (x) 光度.
这样, 最小化问题可以写做:
在计算上式的导, 然后设置为0之后, 解就是:
这里有两个GPU问题:
- memory coalescing: memory合并
- warp divergence
因为VIO一般太不需要很多特征(只有50-200个), 这些稀疏特征是散点, 所以在memory也是散的.
这个算法最小化多层金字塔的每个特征的邻近的正方形的光度误差. 结果, 如果一个warp里的线程处理不同的特征, 那么memory权限会被分开, 那么如果一些特征没有收敛/或者在同一层的不同迭代, 那么有些线程就会闲置. 为了解决这些问题, 一个完整的warp会只处理一个特征. 我们也为了长方形patch优化, 使得其可以被一个warp处理: 在高精度的16x16, 在低精度的8x8. 它解决了warp divergence, 因为warp里的每个线程会处理一样的迭代, 直道他们达到一样的金字塔层数.
我们的方法的牛逼之处就是 线程-对-特征的派发.
3. Evaluation
A. 硬件
我们用了NVIDIA Jetson TX2, 和Intel i7-6700HQ和一个NVIDIA 960M显卡.
B. 非极大值抑制
C. 特征检测
红圈是原来的检测器, 黄圈是两边, 蓝圈是false-positives.
D. Feature Tracker
E. Visual Odometry
我们用ICE-BA[33]作为后端.
4. Conclusion
可以达到200Hz, 别的没啥.