光流学习（一）：Fast Optical Flow using Dense Inverse Search

光流学习（一）：Fast Optical Flow using Dense Inverse Search

最近读到一篇较好的光流文章：Fast Optical Flow using Dense Inverse Search。opencv中开源了c++、opencl 代码，相关关键词为DISOpticalFlow。个人认为该方法是传统光流的巅峰。而且该方法修改了迭代部分相对较容易做成工程算法并部署到嵌入式平台视频插帧、slam 都要用到光流，这里将这个相关学习做个详细的学习笔记并分享以抛砖引玉，如有错误还望指正。
论文链接：https://arxiv.org/pdf/1603.03590.pdf

一、论文阅读

引论：

1 引论：

翻译：

光流估计一直面临着不断增加质量和速度的压力。这种进步为新的应用提供了可能性。更高的速度使其能够被整合到需要大量后续处理的更大系统中（例如，用于运动分割、跟踪或动作/活动识别的可靠特征），并且能够在计算资源有限的情景下部署（例如嵌入式系统、自主机器人、大规模数据处理）。一个强大的光流算法应该能够处理不连续性（异常值、遮挡、运动不连续性）、外观变化（光照、色度、模糊、变形）和大位移。数十年过去了，Horn 和 Schunck [27] 以及 Lucas 和 Kanade [35] 的开创性研究之后，我们已经找到了针对前两个问题的解决方案 [9, 38]，并且最近的努力在处理大位移方面取得了显著进展 [45, 13, 11, 31, 49, 28, 36, 41, 2, 51, 53, 54, 55, 18]。但这是以通常在计算资源受限的实时应用等情景下通常不可接受的高运行时间为代价的。最近，只有极少数的工作致力于在准确性和运行时间之间取得平衡，以提高效率 [17, 48, 54]，或者利用大规模并行化的专用硬件来实现可接受的运行时间 [4, 40, 18]。相比之下，最近已经注意到对于几个计算机视觉任务 [24, 44, 8, 6]，常常希望权衡强大但复杂的算法，采用简单高效的方法，并依靠高帧率和更小的搜索空间来获得良好的准确性。

本文的重点是在非特定领域场景下提高光流的速度，同时保持接近最先进的光流质量水平。我们提出了两个具有低时间复杂度的新组件，一个利用逆向搜索进行快速补丁对应，另一个基于多尺度聚合实现快速密集光流估计。此外，一个快速的变分细化步骤进一步改善了基于我们密集逆向搜索方法的解决方案。总体而言，在与最先进方法相当的光流质量操作点上，我们实现了1-2个数量级的速度提升（图1）。我们使用普通台式计算机上的单个CPU核心，在1024×436分辨率图像上的运行时间范围为10-600 Hz，达到了人类生物视觉系统的时间分辨率 [8]。据我们所知，这是首次在任何硬件上以如此高的光流质量达到了数百帧每秒的速度。

理解：
其实引论主要是陈述，这里图1还是有点意思的。作者引用了[8]中的经验结果。[8]通过实验描述了一种快速的视频流采集方法，我这里没有详细阅读其中引论中的 10Hz 时间是怎么评估的我这里进行了省略，但至少引论说明一个问题：DIS光流在保证相同效果下真的很快。（[8]论文链接：https://www.researchgate.net/profile/Sio-Hoi-Ieng/publication/258355234_Event-Based_Visual_Flow/links/0deec52bc305b0c69b000000/Event-Based-Visual-Flow.pdf）

引论还提到了HS光流和 LK光流，这里就不赘述。自从HS和LK光流后大家就没有太多很好既兼顾精确度又兼顾快速，言外之意就是本文做到了。

1.1 相关工作：

略

1.2 本文贡献：

翻译：
我们提出了一种基于密集逆向搜索（DIS）的新型光流方法，我们展示其在单个CPU核心上以10-600 Hz的速度提供高质量的流估计。与Sintel[14]和KITTI[21]数据集上的先前结果[51, 49, 54]相比，考虑所有方法在相近的流质量操作点时，该方法的速度快了1-2个数量级。同时，与以相同速度运行的现有方法[17, 35]相比，它的准确性显著更高。这一结果基于两个主要贡献：
快速逆向搜索对应：受到Baker和Matthews的逆合成图像对齐[3, 1]的启发，我们设计了我们的逆向搜索过程（在§2.1中解释），用于快速挖掘两个输入图像之间基于补丁的对应网格。虽然通常不如穷举特征匹配那样稳健，但我们可以在微秒内提取均匀的对应网格。
多尺度推理下的快速光流：许多方法假设稀疏且无异常值的对应，并且严重依赖于变分细化来提取像素级的流[51, 49]。这有助于平滑小误差，并覆盖具有平坦和模糊纹理的区域，而穷举特征匹配通常无法胜任。其他方法直接依赖于像素级的细化[40, 4]。我们选择了一个折中方案，并提出了一个非常快速且稳健的补丁平均方案，仅在提取基于网格的对应后每个尺度执行一次。这一步增强了对异常值对应的鲁棒性，并初始化了一个每个尺度执行一次的像素级变分细化。我们在单个CPU核心上在准确性和速度之间达到了300Hz的最佳折中点，并在牺牲准确性的情况下达到了600Hz的速度，这两个操作点在图1、4、5中标记为(2)和(1)。

与我们方法相关的是[40]。在这里，逆向图像扭曲的想法[1]被用于所有像素，而我们的方法则独立优化补丁。与我们每个尺度执行一次的密集化相反，他们依赖频繁的流插值，这需要高性能的GPU，并且仍然明显慢于我们仅使用CPU的方法。本文结构如下：在§2中，我们介绍了我们的DIS方法。在§3中，我们描述了实验，单独评估了基于补丁的对应搜索，并分析了具有和没有变分细化的完整DIS算法。在§4中，我们总结了本文。

理解：
本文贡献两个：
1 快速逆向搜索对应：
逆向搜索作者说受到[3]的启发（论文链接：https://www.ncorr.com/download/publications/bakerequivalence.pdf)这篇文章对于DIS光流有两个贡献点：
（1）通过连续性， 通过△p可以获得warp 参数
这里逆向搜索过程的数学圆形早在2001年就已经提出，期初是用在图像对齐上的，这里有个原始假设：
我们运动的物体都是连续的，因此可以求解一阶导数，于是有了泰勒展开公式如下：

这是一个很重要的基本假设，如果突然出现场景变换，例如突然闪亮，视频里的场景切换，那么这个假设就不成立，因为泰勒展开就不成立，因此光流就不成立。等式2 到 5 说明了 可以通过每次迭代计算△p 而不是从新计算光流。
这里对公式(3)进行推导：

${\sum }_{i=1}^n[I(W(x;p))+▽I\frac{\delta W}{\delta p}-T(x){]}^2$

=> ${\sum }_{i=1}^n[I(W(x;p))-T(x)+▽I\frac{\delta W}{\delta p}{]}^2$ 为了书写方便， 将 I(W(x;p)) -T(x) 用 g(x)替代 将$\frac{\delta W}{\delta p}$ 用 J(x)表示
=>${\sum }_{i=1}^n[g(x)+▽I\frac{\delta W}{\delta p}{]}^2$为当前最小值 相当于 求导为0
=>$\frac{\delta ||g(x)+J(x)▽I|{|}^2}{\delta p}=0$
=>上面展开二范式$||g(x)+J(x)▽I|{|}^2$得到$g(x{)}^2+J(x)J(x{)}^T▽I^2+2g(x)J(x)▽I$ 注意：▽I 是梯度也就是微分项
=>然后对▽I求导：$2J(x)Jx(x{)}^T▽I+2g(x)J(x)$ 让导数等于0得到最小值
=>$J(x)Jx(x{)}^T▽I+g(x)J(x)=0$
=>$J(x)J(x{)}^T▽I=-J(x)g(x)$ 注意这里 $J(x)J(x{)}^T$ 其实是hessien 阵，但对于图像，x方向和y方向完全是正交的，可以近似用 $J(x)J(x{)}^T$替代，也就是说 $df(x,y)/dy\ast df(x,y)/dx\approx d(df(x,y)/dx)$ 注意这里很重要这个思想也被带到了DIS光流，甚至opencv 也是这么求解hessien阵。自此，我们得到 上图公式(4)和(5)。

同时这篇论文讲了如果是正向求解和反向求解在哪里：

上图左侧是正向计算，右侧则是反向搜索，差别在于正向需要算导数图像的jacob矩阵。为什么差别3、4这三步骤？因为作者指出，如果是正向，那么需要 warped gradient ▽I和 雅可比$\frac{\delta W}{\delta p}$ 都是依赖于P。所以正向求解求的是 △P，反向求解求的是 warp( △P)。
这里图右边其实是Compositional Image Alignment 的inverse过程，2001年那篇文章的作者证明了 additive 和 compositional 是等价的。所以我用来说明正向和反向的区别（注意这里是inverse alignment 和 我们DIS光流中的 inverse search 还是不一样的）。

2 多尺度推理下的快速光流：
多尺度是光流的传统方法，我这里提一嘴：一般光流都会建立由粗到细（coarse-to-fine ）的金字塔，这样可以解决大位移光流估计不准确问题下图是一个网上的金字塔模型case：

2 本文方法：

我这里总结了DIS光流主要步骤，没有完全按照论文死板教条，因为理解是王道 \(^_^)/
2.0 网格化
2.1 反向查找过程
2.2 多尺度推理
2.3 稠密化
2.4 变分优化
DIS光流的整体算法过程如下截图

2.0 ：网格化

DIS光流论文中没有该章节，但是我翻看了opencv源码后这个过程是比较重要的，如果没有这个过程，DIS光流作者说的“The first step requires extraction and bilinear interpolation” 是无法理解的。注意这里因为有overlap 所以后面的很多过程都和这个patch 和overlap有关。

opencv 中的 patch size 为 8×8 ，其实按照正常的假设，这个8×8的patch u 和 v 是一致的，既这个8*8范围内，用我们求得的u 和 v 会有$argmin={\sum }_{i=1}^{8\ast 8}[I(x+u,y+v)-T(x,y){]}^2$

2.1 ：逆向查找对应点

翻译：略，这里理解重要
理解：2.1 小节讲解了怎么样快速地找到两帧之间的对应点，以此来求得光流。本小节的核心计算公式还是基于我上面截图中的公式，这里再列下：
$argmin={\sum }_{i=1}^n[I(W(x;p+△p))-T(x){]}^2$

也就是说本文其实并没有用inverse alignment 而是 inverse search
。在本文有了如下假定：
1 给定的patch：θps × θps
2 将图像分成patch，patch 的中线点坐标：x = （x,y）公式如下：
$argmin={\sum }_x[I_{t+1}(x+△u^{\prime })-T(x){]}^2$

dis光流作者也说明DIS光流这里查找对应点的方法是 inverse L-K光流。不难看出上面公式是个非线性的最小二乘问题。该问题求解时用到梯度下降，也就是我下面的截图（截图来自opencv，因为opencv 中已经实现了 dis光流）
红框内框出的是反向search，每次计每个patch的光流时，需要查找x方向-1 的和 y 方向-1的patch。

所以在这里，和dis光流作者提到的那篇[3](https://www.ncorr.com/download/publications/bakerequivalence.pdf 这篇文章启发了作者)中提到的通过“inverse alignment” 不同的，这里只是反向查找，查找出来还是计算的是 du 和 dv 也就是dx 和 dy。
梯度下降求解current ux 和 current uy

注意这里每次计算完后得到dx 和 dy也就是文中说的△u
其实逼着

2.2 ：多尺度光流推理

翻译：略
理解：
这里其实就是在多尺度上完成求解光流的工作，这里我翻看了openv 源码，opencv 的金字塔层数为 10 - 2 +1 = 9层

每层金字塔按照：
(1)反向查找
(2)稠密化
(3)refine
并且当前层的光流结果带入下一层计算

2.3 ：稠密化

翻译：略
理解：
这里的稠密化是为了将当前金字塔层的稀疏光流UxUy放大尺寸到下一层。所以稠密化只是将光流矩阵长宽个放大2倍，借助的计算还是SSD，只是将SSD的计算diff结果作为惩罚项，因为差的平方的自然量纲放到论文中都不合适只有做惩罚项最合适。
其中有个概念需要声明下，DIS光流的稀疏Ux 和 Uy 是在原尺寸上有overlap，也就是说一个1/4 patch 其实杯4个patch（图中是 橙色、红色、粉色、绿色这四个patch）共享。这是下面公式中λ需要通过双线程差值求得。

论文中给了下图公式：

这个公式有两点信息：
1 λ 是个权重，这个权重通过双线性差值得到：

2 分母是个二范数，既 $d_i(x)=[I_{t+1}(x+u)-T(x){]}^2$，其中max 1 是为了保护系数，而这个二范数说明希望通过SSD来作为惩罚项，既差别过大说明不是特征相近的块，因此要给的权重小一些。opencv 代码和 这里公式高度统一。

2.4 ：变分优化（其实就是变分光流估计模型）

翻译：略
理解：

是本文的光流平滑项其中 $\Psi (x)=\sqrt{x^2+0.001^2}$
这里要说的是，dis光流作者提到了论文[51]、[46]、[59]，其中[51]论文关于“refine“这部分的通篇其实介绍非常详细，[46]纯属水，[59]则是讨论了文章当时截止是常用的 平滑项，并进行了数学分析（比较麻烦笔者没有深究）。因此我们结合[51]进行详细展开[51] 的能参考的论文链接地址其实是这里：https://drive.google.com/file/d/1XCo-Esv8iQru5ujBC4vKjgu5WOyIepbY/view
我们通过其内容截图也可以看到关联度很大，如下图

ok，我们继续理解。
DIS光流中作者用的是$E_I$是亮度惩罚项，$E_G$是梯度惩罚项,$E_S$光滑项惩罚项。注意dis光流作者说了这里：
dis光流作者也说le$E_I$和 $E_G$ 就是一般意义上的光流光流不变性，这就成了HS模型：
$I(x,y,t)=I(x+u,y+v,t+1)$ 这里是理想情况，也就是亮度不变，但实际让这两个近似就好，也就是
$u\frac{\delta I}{\delta x}+v\frac{\delta I}{\delta y}=\frac{\delta I}{\delta t}$接近0最好，越小越好
于是光流的的解变为了：
$E(u,v)=\int \gamma \Psi (|I(x+u,y+v,t+1)-I(x,y,t){|}^2)+\sigma \Psi (|dx{|}^2+|dy{|}^2)+\alpha \Psi (|u{|}^2+|v{|}^2)d(x,y,t)$
于是这成了变分问题，什么样的 u、v、t 让 E最小。
详细见本小节提到的相关文章[12]，链接：High accuracy optical flow estimation based on a theory for warping。这篇文章指出了怎么解这个变分问题，opencv 解法是用了 Red-Black SOR 也是dis光流作者提到的，Red-Black SOR貌似是一种 解法：https://en.wikipedia.org/wiki/Gauss%E2%80%93Seidel_method
完全没看懂。。。
第一阶段分享就到这里

鸣谢