光流法和直接法介绍

知也无涯12345

已于 2023-10-11 11:04:44 修改

阅读量768

点赞数 3

分类专栏： slam 文章标签：计算机视觉人工智能

于 2023-10-09 16:29:17 首次发布

本文链接：https://blog.csdn.net/m0_37604894/article/details/133698780

版权

slam 专栏收录该内容

8 篇文章 4 订阅

订阅专栏

1、光流法

光流为追踪源图像某个点在其他图像中的运动，一般分为稀疏光流和稠密光流，其本质是估计像素在不同时刻图像中的运动

稀疏光流以Lucas-Kanade(LK)光流为代表
稠密光流以Horn–Schunck(HS)光流为代表

在这里插入图片描述
本文主要介绍LK光流。

2、光流基础原理

设 t 时刻位于 x,y 处像素点的灰度值为 $I (x, y, t)$ ,在 t+dt 时刻,该像素运动到了 $I (x + d x, y + d y, t + d t)$ ,光流希望计算运动 dx, dy。这里引入灰度不变假设 $I (x + d x, y + d y, t + d t)$ = $I (x, y, t)$ （注：实际当中由于相机曝光、阴影、材质等不同,很可能不成立）。对 t+dt 时刻的灰度进行Taylor展开并保留一阶项:
在这里插入图片描述
由于灰度不变,所以

因此，两边除以 $d t$

希望求解 $d x / d t$ , $d y / d t$ ,但本式是一个二元一次线性方程,欠定,因此，需要引用额外的约束，这里假定一个窗口[w,w]内光度不变，则可以通过超定最小二乘解求得运动 u,v。
在这里插入图片描述

LK光流的结果依赖于图像梯度,但梯度不够平滑,可能剧烈变化,局部的梯度不能用于预测长期图像走向,因此采用图像金字塔，通过coarse to fine的方式进行跟踪。

LK光流问题可以看成最小化像素误差的非线性优化，每次使用了Taylor一阶近似,在离优化点较远时效果不佳,因此往往需要迭代多次，并且对于快速运动来说可以通过IMU或者patch匹配提供先验，在先验的基础上再进行光流跟踪。

3、不同光流方法

按算法分类,有两种分法,一种可以分为迭加法或组合法(additive or compositional),向前或者反向算法( forwards or inverse )
(有人翻译:叠加式(additive)和构造式(compositional),前向(forwards),逆向(inverse))。从而导出四种图像对齐算法( all four image alignment algorithms ),分别是

Forward Additive(FA),
Forward Compositional(FC)
Inverse Compositional(IC) 算法 ,
Inverse Additive(IA)算法。

a、forward-addtive Gauss-Newton 光流

设有图像 1.png,2.png,我们在 1.png 中提取了 GFTT 角点 ,然后希望在 2.png中追踪这些关键点。设两个图分别为 $I_1$ , $I_2$ ,第一张图中提取的点集为 $P = {p_i }$ ,其中 $p_i = [ x_i , y_i ]^T$ 为像素坐标值。考虑第 i 个点,我们希望计算 ∆ $x_i$ , ∆ $y_i$ ,满足:
在这里插入图片描述

for (int x = -half_patch_size; x < half_patch_size; x++) {
	for (int y = -half_patch_size; y < half_patch_size; y++) {
		// TODO START YOUR CODE HERE (~8 lines)
		float p_x = kp.pt.x;
		float p_y = kp.pt.y;
		int k = x + y + 2 * half_patch_size;
		double error = 0;
		Eigen::Vector2d J; // Jacobian
		if (inverse == false) {// Forward Jacobian
			J[0] = GetPixelValue(img2, p_x + x, p_y + y) - GetPixelValue(img2, p_x + x - 1, p_y + y);
			J[1] = GetPixelValue(img2, p_x + x, p_y + y) - GetPixelValue(img2, p_x + x, p_y + y - 1);
		} else {// Inverse Jacobian
			J[0] = GetPixelValue(img1, p_x + x, p_y + y) - GetPixelValue(img1, p_x + x - 1, p_y + y);
			J[1] = GetPixelValue(img1, p_x + x, p_y + y) - GetPixelValue(img1, p_x + x, p_y + y - 1);
			// NOTE this J does not change when dx, dy is updated, so we can store it and only compute error
		}
		// compute H, b and set cost; //H = JT * J
		H += J * J.transpose();
		error = GetPixelValue(img2, p_x + x+ dx, p_y + y + dy) - GetPixelValue(img1, p_x + x, p_y + y);
		b += -error * J;
		cost += error * error / 2;
	}
}

上面只是最简单的一种实现方式，认为窗口内像素只有平移运动，没有旋转变换，实际使用时我们会将cost function模型化为下面这种形式：
在这里插入图片描述
其中 $I$ 表示 $I_2$ 目标图像， $T$ 表示提取特征点的模板图像， $W (x, p)$ 为warp函数，最简单的warp就是：

此时认为窗口内像素块只有平移运动，但如果窗口内的像素存在旋转，那么两个像素块的匹配，只是 $p_1,p_2$ 就无法表达，需要使用下面仿射变换：
在这里插入图片描述
Lucas-Kanade算法假设已知p的当前估计值，然后迭代地求解参数p的增量；也就是说，以下表达式被（近似地）最小化

然后通过如下方式更新：

GN方法求解增量为：

H为：

上面warp对应jacobian为：

b、inverse-addtive Gauss-Newton 光流

在迭代开始时, Gauss-Newton 的计算依赖于 $I_2$ 在 $x_i +∆x i , y_i + ∆y i )$ 处的梯度信息。然而,角点提取算法仅保证了 $I_1 (x_i , y_i )$ 处是角点(可以认为角度点存在明显梯度),但对于 $I_2$ ,我们并没有办法假设 $I_2$ 在 $x_i , y_i$ 处亦有梯度,从而 Gauss-Newton 并不一定成立。反向的光流
法(inverse)则做了一个巧妙的技巧,即用 $I_1 (x_i , y_i)$ 处的梯度,替换掉原本要计算的 $I_2 (x_i + ∆x i , y_i + ∆y i)$ 的梯度。这样做的好处有:

$I_1 (x_i , y_i)$ 是角点,梯度总有意义;
$I_1 (x_i , y_i)$ 处的梯度不随迭代改变,所以只需计算一次,就可以在后续的迭代中一直使用,节省了大量计算时间。

具体代码实现，上面的inverse为true时，即对应inverse的方法。

实验结果：正向法和反向法在计算中体现在计算关键点在 $x_i,y_i$ 处的梯度上，正向法采用当前追踪图像 $I_2$ 的图像像素梯度,反向法使用 $I_1$ 的像素梯度。通常单层图像的反向法比正向法效果要好。

for (int x = -half_patch_size; x < half_patch_size; x++) {
	for (int y = -half_patch_size; y < half_patch_size; y++) {
		// TODO START YOUR CODE HERE (~8 lines)
		float p_x = kp.pt.x;
		float p_y = kp.pt.y;
		int k = x + y + 2 * half_patch_size;
		double error = 0;
		Eigen::Vector2d J; // Jacobian
		inverse = true;
		if (inverse == false) {// Forward Jacobian
			J[0] = GetPixelValue(img2, p_x + x, p_y + y) - GetPixelValue(img2, p_x + x - 1, p_y + y);
			J[1] = GetPixelValue(img2, p_x + x, p_y + y) - GetPixelValue(img2, p_x + x, p_y + y - 1);
		} else {// Inverse Jacobian
			J[0] = GetPixelValue(img1, p_x + x, p_y + y) - GetPixelValue(img1, p_x + x - 1, p_y + y);
			J[1] = GetPixelValue(img1, p_x + x, p_y + y) - GetPixelValue(img1, p_x + x, p_y + y - 1);
			// NOTE this J does not change when dx, dy is updated, so we can store it and only compute error
		}
		// compute H, b and set cost; //H = JT * J
		H += J * J.transpose();
		error = GetPixelValue(img2, p_x + x+ dx, p_y + y + dy) - GetPixelValue(img1, p_x + x, p_y + y);
		b += -error * J;
		cost += error * error / 2;
	}
}

c、forward-compositonal Gauss-Newton 光流

相对于additive的目标函数：
在这里插入图片描述

compositional的目标函数为：
在这里插入图片描述

也就是说，组合方法迭代地求解增量warp W(x; p)，而不是参数p的加法更新。compositional方法和additive方法在一阶p近似中如下等式成立，即：
在这里插入图片描述

展开后为：

compositional操作通常涉及到乘两个矩阵来计算组合warp的参数，就像上式中的仿射warp一样。对于更复杂的warp，两个warp的组合可能更为复杂。因此，我们对warp集合有两个要求：（1）warp集合必须包含identity warp；（2）warp集合必须在组合下是封闭的。因此，变形集合必须形成一个半群。这个要求并不难达到，大多数在计算机视觉中使用的warp，包括单应性和3D旋转，自然就形成了半群。

d、inverse-compositonal Gauss-Newton 光流

cost function为
在这里插入图片描述
更新为：

其中warp的inverse为：

可以看到，向前方法对于输入 $I$ 图像进行参数化(包括warp变换及warp增量). 反向方法则同时参数输入 $I$ 图像和模板图像, 其中输入 $I$ 图像参数化warp变换, 模板图像参数化warp增量. 因此反向方法的计算量显著降低.。由于图像灰度值和运动参数非线性, 整个优化过程为非线性。

4、小结

光流法基于灰度不变假设,在运动较小时基本是合理的。当运动较大时,没法计算对应特征点,光流跟踪不了,我们的区块小于实际运动的位移,无法实现光流的基本假设。光流金字塔就是一种解决方法。此外，也可以使用IMU进行预测，或者基于匀速假设预测，或者通过patch匹配的方法获取初值，结合金字塔方法优化光流估计。
采用金字塔的方法,即窗口固定,将图像生成金字塔,在每一层金字塔上都用同一个大小的窗口来进行光流计算,这样很好的去解决了图像块的问题。这样一来图像块大小并不会带来明显差异。如果没有采用金字塔方法,我们首先要设置一个邻域的窗口,之后计算光流。当窗口较大时,光流计算更鲁棒,当窗口较小时,光流计算更正确。原因在于,当图像中每一个部分的运动都不一致的时候,如果开的窗口过大,很容易违背窗口(邻域)内的所有点光流一致的基本假设,这可能与实际不一致,所以窗口小,包含的像素少,更精确些。这些需要根据实际使用效果进行参数优化。
金字塔层数一般越多效果越好,但是一般图像大于 4~5 层之后都变得太小,特征点像素太过紧密容易出现错误追踪。放大倍率同样的道理,放大倍率小,金字塔的层数可以增加,迭代层数增多,效果自然要好。此时建议使用其他方法提高初值的可靠性。

5、直接法

a、直接法基本思想

光流其实和直接法很相似，直接法或者半直接法更像进化版的光流，主要是光流仅估计了像素间的平移，或者warp变换，但

没有用到相机本身的几何结构、没有考虑相机的旋转和图像缩放，简言之，运动模型过于简单，或者只是近似
对于边界上的点，光流不好追踪
而直接法考虑了上面这两点信息，那么直接法是怎么做的呢

b、基本原理

假设有两个帧,运动未知,但有初始估计 R,t，第1帧上看到了点P,投影为p1，按照初始估计,P在第2帧上投影为p2，投影关系为
在这里插入图片描述

在光流中,我们估计的是每个像素的平移(在 additive 的情况下)，而在直接法当中,我们最小化光流误差,来估计相机的旋转和平移，因此建立如下最小二乘问题，最小化光度误差

这里待估计的量为相机运动，因此需要计算误差相对于相机位姿的导数，jacobian计算就是标准的链式求导：
在这里插入图片描述
为简化表示，令：

三部分连乘：

第一部分：图像梯度
第二部分：像素对投影点的导数
第三部分：投影点对位姿的导数

c、直接法小结

直接法的雅可比项有一个图像梯度因子，因此，在图像梯度不明显的地方,残差对相机运动估计的贡献就小
根据使用的图像信息不同，直接法又分为三种

稀疏直接法:只处理稀疏角点或关键点
稠密直接法:使用所有像素
半稠密直接法:使用部分梯度明显的像素

光流主要还是对角点进行跟踪，而直接法一般可以不提取角点，梯度明显的点甚至整幅图像都可以
直接法的基本思想还是像素灰度主导优化方向，因此为了使优化问题找到正确的解，需要保证从初始估计到最优估计中间的梯度一直下降，这一特点就导致直接法很容易收到图像的非凸性影响，当然工程上通常用金字塔来缓解这个问题
inverse和compositional对直接法没有意义，直接法是通过求解位姿来对应第二帧关键点位置,如果采用 inverse的话,就会使第二帧的信息失效,没有意义。compositional 同理。
直接法的优缺点：