三维重建笔记_光束平差法(Bundle Adjustment, BA)

为了求解相机的位姿(R,t)和空间三维点的3D坐标，常采用优化的方式，优化方法中常用的是光束平差法（BA，又叫捆绑调整，束调整等）。BA方法的核心是最小化重投影误差，即有一个图像中的2D点m，算得一个三维点初始的3D坐标P时，将三维点重新投影到图像平面得到m'，最小化原始图像点m与重投影回去的m'之间（所有图像所有对应的图像点）的误差。可以看出，最小化投影误差，是一个非线性最小二乘问题。非线性最小二乘问题求解，常采用的方法是线性化，即使用泰勒展开式展开，忽略高阶项，变成线性最小二乘问题。

对于最小二乘问题，求解方法主要有梯度下降法（一阶梯度法）、牛顿法（二阶梯度法）、高斯牛顿法、LM法，按照这些方法的排序来看，后面一种方法可以说是前面一种方法的改进，每种方法都改进了它前面方法的一些缺点。

一阶，二阶梯度法共有原理

另外一种同样的表示形式：

最速下降法（也叫梯度下降法）--最简单暴力的方法

问题：我们能求出每一个x处的导数J（函数下降速度与方向），怎样求这个函数的最小值？
核心思想：沿着函数变小的方向移动，导数J代表了函数的变化趋势，因此只要顺着函数值变小趋势的方向移动自变量就能迭代得到函数最小值了。
缺点：步长不好确定，往函数最小值方向移动的步长需要人为确定，步长取得太大则容易在最低点附近来回波动，甚至不收敛。取得太小则收敛速度慢
改进思路：利用其他信息来确定步长，想办法确定每一步的步长，使函数技能比较快速的下降，又不会出现不收敛的情况。

牛顿法--计算量大，但是收敛速度快且精确的方法

问题：求F(x0 + dx) 的最小值。注意自变量是dx
核心思想：把F(x0 + dx)泰勒展开成二次函数求最值（注意这里dx是自变量）。
缺点：二阶导数H计算量大
改进思路：避免使用二阶导数H,是否有其他可用信息。
增量正规矩阵：

高斯牛顿法--一种近似型解法，用线性函数近似二次曲线

问题：求F(x0 + dx) 的最小值。注意自变量是dx
核心思想：把F(x0 + dx)泰勒展开成一次函数求最值（注意一次函数必须有限制条件才能求到最值），这样可以避免求海塞矩阵H(使用雅克比矩阵J近似,)。利用限制条件F（x0 + dx）>= 0(针对SLAM问题，因为SLAM中原函数F(x)代表误差大小)
缺点：
- 1、JTJ只有半正定性，若JTJ为奇异矩阵，则会导致算法不收敛。
- 2、步长取太大也可能导致算法不收敛
  - 病态矩阵：求解方程组时如果对数据进行较小的扰动，则得出的结果具有很大波动，这样的矩阵称为病态矩阵。
增量正规矩阵：（相比牛顿法，高斯牛顿法先将残差函数f(x)进行一次泰勒展开，再代入F(x)--即求平方）

LM法--在近似型解法的基础上添加上了可信区间

问题：求F(x0 + dx) 的最小值。注意自变量是dx
核心思想：：把F(x0 + dx)泰勒展开成一次函数求最值，利用函数在每个迭代点处的线性度，来指导更多的使用梯度下降法还是高斯牛顿法
做法：
- 在原来的 H (Hessian Metrix)的对角线上加上一个可以调节的 λ ，当 λ 较小时，重要还是Hessian在起作用，这个时候，该算法就接近于 Gaussian‐Newton 法；当 λ 比较大的时候，就可以忽略H，成为一个梯度下降的方法；即通过调节 λ ，在 Gauss‐Newton method 和 Gradient descent 找到一个平衡。
怎样调节λ？
- 往前执行一步后，看一下，error 是否真的下降了，如果error确实下降了，说明这个 λ 比较合适，则接受这样的迭代，同时减小 λ 一点，让搜索更加精细化一些。如果上一个搜索，使能量函数 error 上升了，说明这个 λ 不合适；这时，就放弃这一步的优化，并 rollback，同时增大 λ ，继续优化。通过这样一个过程，可以得到一个比 Gauss‐Newton method 更好的性能。
缺点：速度有所降低
增量正规矩阵：

1 光束平差法(Bundle Adjustment, BA)

即，使用重投影误差法，怎样做优化（最小化重投影误差）。

bundle: 通过同一个三维点的很多条光线组成的一捆光线集合。（来源：photogrammetry领域，大约1955年，2000年左右引入三维重建领域，然后成为三维重建领域的一个标配。）

1.1 SFM (Structure‐from‐Motion)

1.1.1 定义

从运动相机（可以是网络图片，也可以是 video sequence）中，恢复出场景的三维坐标（点云）。

-------------------------------------------------------------------------

1.1.2 问题描述

输入、输出、目标函数如下所示

-------------------------------------------------------------------------

1.2 重投影误差 (Reprojection Error)

最小化通过相机矩阵投影计算出来的图像特征点坐标 和 检测(detect)出来的图像特征点 之间的误差。

同时希望，计算出来的图像特征点坐标能与其他图像（第二、三...n个相机）计算结果接近。

-------------------------------------------------------------------------

1.2.1 目标函数

是一个非线性的最小二乘法问题

-------------------------------------------------------------------------

1.3 使用 Levenberg‐Marquardt Algorithm 进行优化

1.3.1 问题描述

-------------------------------------------------------------------------

1.3.2 简化目标函数

-------------------------------------------------------------------------

1.3.3 参考 Gauss‐Newton Method

非线性优化都很难，最简单的办法就是泰勒展开，因为泰勒展开，可以把二次项以后的都丢掉，然后就变成线性函数了。“线性函数的东西大家都会，非线性的大家都不会，这是优化里的一个现状。”

所以先把它通过泰勒展开，变成一个线性的，即，给定一个初始值P0，在P0处做泰勒展开。

-------------------------------------------------------------------------

参考：数值优化之高斯-牛顿法（Gauss-Newton）：https://blog.csdn.net/qq_39521554/article/details/79919041

非线性最小二乘问题的求解方法: https://www.cnblogs.com/hustyan/p/11239501.html

1.3.4 使用Gauss-Newton法，优化重投影误差

-------------------------------------------------------------------------

1.3.5 使用 Levenberg‐Marquardt Algorithm，优化重投影误差

Levenberg‐Marquardt Algorithm 是对 Gauss‐Newton method 的一个优化修改。

考虑到 Gauss‐Newton method 收敛得很快，但当初始值选定的不合适时，有可能不收敛；

而另一种方法：梯度下降法（Gradient descent），总是收敛的；

所以 Levenberg‐Marquardt 方法想要达到兼有两者的优点。

-------------------------------------------------------------------------

做法：

在原来的 H (Hessian Metrix)的对角线上加上一个可以调节的 λ ，当 λ 较小时，重要还是Hessian在起作用，这个时候，该算法就接近于 Gaussian‐Newton 法；当 λ 比较大的时候，就可以忽略H，成为一个梯度下降的方法；即通过调节 λ ，在 Gauss‐Newton method 和 Gradient descent 找到一个平衡。

怎样调节λ？

往前执行一步后，看一下，error 是否真的下降了，如果error确实下降了，说明这个 λ 比较合适，则接受这样的迭代，同时减小 λ 一点，让搜索更加精细化一些。

如果上一个搜索，使能量函数 error 上升了，说明这个 λ 不合适；这时，就放弃这一步的优化，并 rollback，同时增大 λ ，继续优化。

通过这样一个过程，可以得到一个比 Gauss‐Newton method 更好的性能。

-------------------------------------------------------------------------

H ( Hessian Matrix) 以及 b的内部结构

根据 J (雅克比矩阵) 行列维数 2 x (3m+6n)，e_ij (第i个三维点在第j个相机的投影误差) 行列维数 1 x 2，可以推出 b_ij 行列维数为 1 x (3m+6n)；

Hij 行列维数 (3m+6n) x (3m+6n); 所以可知，b^T 为一个细长的dense vector； H为一个方阵，元素除对角线外，主要集中在左下角和右上角，下图中，深红色方块表示有不为0的值，浅蓝为值为0的元素。

-------------------------------------------------------------------------

将一个大矩阵的逆的求解变成一个小矩阵的逆的求解，从而大大大大简化了计算。

-------------------------------------------------------------------------

H ( Hessian Matrix ) 的一个样子

下图中黑色表示非零的元素，可以看出它的位置；浅蓝色部分（是 block‐diagonal的）的矩阵的逆，很容易求；从而使计算得到了简化、加速了计算。

-------------------------------------------------------------------------

Levenberg‐Marquardt 法优化的一些工具包

使用方法，基本上都是定义好一个 objective function（这里指 reprojection error），然后丢给下面的算法，这些算法就可以自己进行优化，最后可以得到像素点的三维坐标，和相机的朝向、位姿。

-------------------------------------------------------------------------

1.3.6 参考2：L-M 方法

2 如何参数化相机的旋转矩阵R?

即如何描述三维空间中的旋转问题

2.1 问题引出：雅克比矩阵中的对Rj进行求偏导问题，应该怎么求？

-------------------------------------------------------------------------

2.1.1 方法1：欧拉角（不好的方法）

-------------------------------------------------------------------------

2.1.2 方法2：Axis‐Angle（较常用）

-------------------------------------------------------------------------

2.1.3 方法3：四元数（较常用）

-------------------------------------------------------------------------

参考：四元数和旋转(Quaternion & rotation) https://zhuanlan.zhihu.com/p/78987582

3 Bundle Adjustment 如何初始化

BA是一个非线性优化，所以足够好的初始化是很重要的。非线性优化一般使用Bundle Adjustment来做最小化。

Bundle Adjustment问题的核心在于，重投影误差（reprojection error）的优化；我们根据三维点的坐标、根据相机的姿态，把三维的点投影到图像平面中，然后我们来算它跟图像中，通过特征提取出来的那个点，之间的 distance，我们去minimize这个distance, 这是一个非线性优化。

非线性优化，除了少数的特殊的非线性问题，绝大多数的非线性优化，基本上大家也没有什么好办法，那无非就是梯度下降这些东西，这就非常取决于你的初始化，你的初始值在什么地方，那么那个最后的梯度下降，通过 Gasian-Neton、Levenberg‐Marquardt 这样的方法，最后就收敛到，距离你的初始值不是那么远的一个状态上，所以初始化就变得非常关键。

那么在这里我们要初始化些什么东西呢？如下图所示。

-------------------------------------------------------------------------

3.1 增量式的SfM（Incremental Structure‐from‐Motion）

先从一对图像开始，计算对极几何约束（ epipolar geometry），本质矩阵（Essential Matrix）；然后，可以从 Essential Matrix 中 decompose 两个相机的pose (旋转矩阵R及位移t)，然后就可以使用 triangulation，三角化算出一些三维点（误差累计的来源）；所以可以通过两张图，计算出包含两个相机和一些三维点的初始三维重建，接着再使用 Resection 方法，加第三个相机（第三张图片），即使用三维点（已计算点）和两维点（新加点）的对应关系（correspondence），计算出第三个相机的姿态。逐步增加相机，逐步找到三维点，及对应的图像点。

-------------------------------------------------------------------------

3.1.1 最大问题 -- 误差会不断累计

----------------------------------------------------------------------------------------------

-------------------------------------------------------------------------------------------------

如上图所示，依次加入相机（图像）时，误差会不断累计。

比如加第三张图像时，用的是3D-2D对应点的对应关系，来做 Resection 或者做 PnP。我们之前讲到，在做 Resection时，假设3D点坐标是绝对精确的、没有误差的，误差全部来自于2D点坐标（特征点检测出来的像素坐标值），但增量式SfM在添加相机时，使用的三维点坐标，是前面计算出来的。所以三维点的误差，会导致我们新加进来的第三张图的相机的pose（R和t）的误差，进一步会导致我第四张图有些误差。

这样，使用第三张和第四张图对应相机三角化出来的三维点，也会有误差，这样误差逐步被放大。

怎么办呢？目前没有什么办法，非线性优化呗，即 minimize 重投影误差（reprojection error），所以，通常在做了几次添加相机（图像）后，就需要做一把 Bundle Adjustment。这时候就相对安全了，这时候就可以继续添加相机了。这时就会有一个判断准则，比如添加了20张图像（或者判断 triangulate 出来的相机姿态需要优化了），就需要做 Bundle Adjustment 了。

就这样一直添加图像，直到所有的图像全部添加完以后，还要再最后做一次整体的 Bundle Adjustment 。因为BA是在优化所有的相机，所有的三维点。整套这样的流程就叫做 Incremental Structure‐from‐Motion。如下图所示。

-------------------------------------------------------------------------

另外需要注意的问题，最初始的两张图片很重要！下一张图像怎么选择。

-------------------------------------------------------------------------

3.1.2 缺点

-------------------------------------------------------------------------

3.2 全局式的SfM（Global Structure‐from‐Motion）

先给定一对图像对（pair），同样要计算Essential Matrix，然后同样的decompose出来相机的相对运动( Relative Motion，包含：相对旋转矩阵 R 和相对平移矩阵 t )；使用同样的方法，算出所有的图像对的 E和 Motion，比如有n张图，就有个图像对。

然后根据这些所有的图像对，算出的 Relative Motion，解一个优化问题（称之为：Register all cameras），把所有相机的朝向（Rotation）和中心点坐标（center）一次性的求出来。

最后当然离不开 Bundle Adjustment，因为这才是最后要解的 reprojection error 的 minimize。前面初始化的效果越好，后面Bundle Adjustment 的效果就越好。一般情况下，使用 Global Structure‐from‐Motion，最后收敛的也就越快。

-------------------------------------------------------------------------

3.2.1 Register all cameras

通过 Essential Matrices 和 Decomposition，算出来相对的旋转 R 、相对的平移 t，怎样把他们放到世界坐标系中，如何把他们 register 到一起，这个就是这一步要解决的问题。这步一般可以分成两步 Rotation Averaging 和 Translation Averaging

Rotation Averaging（求解所有相机朝向）

先求解相机的朝向，因为相机的朝向和相机的中心点的位置没有关系。所以可以先把相机的朝向独立的求解出来，有了相机的朝向，再去求解相机的中心点会变得容易。

问题的核心是一对图像对，对应的两个相机的，在世界坐标系下的旋转矩阵 Ri 和 Rj 的表述方式。这就是为什么 Rotation Averaging 这么难，到现在为止（2018）还没有一个好的解法。

旋转矩阵表述方式1：四元数

-------------------------------------------------------------------------

旋转矩阵表述方式2：3x3矩阵表示

当然，这不是一个最优的表述方式，因为旋转矩阵并不是有3x3=9个自由度，这里使用9个参数求出来的结果，肯定不是最优的，特别是，当数据中有噪声、有 outlier 时，很有可能 crash 掉。

-------------------------------------------------------------------------

旋转矩阵表述方式3：nonlinear optimization

-------------------------------------------------------------------------

Translation Averaging（求解所有相机中心点）

使用 Rotation Averaging 的计算结果，来进一步计算 Translation。

待求变量 ci 和 cj ，一对图像对，在世界坐标系下的，相机中心点坐标。

-------------------------------------------------------------------------

虽然上述方程是一个线性方程，但是是严重退化的，不能处理沿着一条指向走的相机，因为Essential Matrix不含平移的尺度信息的。在CVPR 2015年有一篇很理论的文章，论述了用 Essential Matrix，可以把什么样的运动 Translation Averaging，把相机的姿态求出来。他的思想是，将相机的中心点当做是一个节点，如果两两之间，当相机有相对运动时，就把他们连接起来，形成一条边，形成了一个图（称之为 pose graph），当这个图是一个 parallel rigid graph （满足一定条件的一种类型的图），这样的 Translation Averaging 才能够解的出来。

还有一种办法是，引入三维点，因为我一对图像对，算了Essential Matrix之外，还算了一些三维点，我们可以用上这些三维点，使得 Translation Averaging 可以work, 同时推出一个线性的 formulation 出来。

-------------------------------------------------------------------------

列出一个求线段中点的方程如下图：

-------------------------------------------------------------------------

M1 M2的几何含义如下图：

-------------------------------------------------------------------------

任何两个相机，如果可以计算Essential Matrix，我们可以给他们连接一条边，这就形成了一个graph，然后每个三角形（3个相机组成的）可以列出一组线性方程，多个这样的三角形，就可以列出一个关于相机中心的线性方程组，然后就是老套路了，进行SVD，求解所有相机中心点的位姿。这就是 Translation Averaging 比较 general 的一个解法，可以实用到一般的相机运动下，然后，仍然是一个比较 linear 的 formulation。

-------------------------------------------------------------------------

参考

GitCode 开源社区

旨在为数千万中国开发者提供一个无缝且高效的云端环境，以支持学习、使用和贡献开源项目。

更多推荐

[转载]在Windows环境下安装GNU Radio

转自：在Windows环境下安装GNURadio_恐弱智_新浪博客GNU Radio是用Python开发的，大部分开源的工程能够在Linux环境下运行良好，而Windows下却运行的很勉强，而且安装配置都很复杂。GNU Radio算是个例外了，不光提供了Windows的二进制安装，还有比较详细的说明。我是Python小白，所以折腾了好久才弄好，特意记录下来，免得以后再装还折腾。GNU Radio的

GitCode 开源社区

centOS 8 使用dnf安装Docker

DNF是什么？CentOS 8使用YUM软件包管理器版本v4.0.4。现在，该版本使用DNF(已删除YUM)。DNF是软件包管理器。它会在Linux发行版上安装，执行更新并删除软件包。使用DNF安装Docker跳过具有损坏依赖性的程序包一个有效的解决方案是使您的CentOS 8系统使用以下--nobest命令安装最符合条件的版本：sudo dnf install docker...

GitCode 开源社区

定时同步数据库表(mysql+linux+crontab)

sync.sh里面的参数需要改变，ip/username/password/database/tablesync.sh#!/bin/sh# Please change the IP and password of the data source db.# Then change the table name.filename=/home/nington/db/$(date +%Y-%m