openMVG--GlobalSFM(原理及代码解读)

1、相关概念

1. globalSFM模块，集成较先进的算法（相较于incrementalSFM）
2. 依赖相对定向的质量，鲁棒性较差
3. sfm结果的精度评价基于重投影误差
4. 代码参考论文（阅读笔记链接）：Global Fusion of Relative Motions for Robust, Accurate and Scalable Structure from Motion

2、默认参数设置

3、使用注意：

重要！！！！

4、算法流程

（官方链接）

Require: internal camera calibration (possibly from EXIF data)
Require: pairwise geometry consistent point correspondences
Ensure: 3D point cloud
Ensure: camera poses

compute relative pairwise rotations
detect and remove false relative pairwise rotations

• using composition error of triplet of relative rotations
  compute the global rotation
• using a dense least square and approximated rotations
  compute relative translations
• using triplet of views for stability and colinear motion support
  compute the global translation
• integration of the relative translation directions using a l-∞ method.
  final structure and motion
• link tracks validated per triplets and compute global structure by triangulation,
• refine estimated parameter in a 2 step Bundle Adjustment
  • refine structure and translations
  • refine structure and camera parameters (rotations, translations).

5、代码实现过程

5-0 globalSFM接口函数：

GlobalSfMReconstructionEngine_RelativeMotions::Process()

5-1接口函数内部实现：

5-1-0

几个graph的node-edge变化

1. Geometry graph—view graph: biedge 约束 大部分情况下对pose-graph影响较小

2. View graph—pose graph : 5点法 ransac 进行E估计时，model(E约束)内点数<2.5 minimun_sample(5),该view pair在pose graph中 舍去 大部分情况下对pose-graph影响较小
   实现：

RobustRelativePose函数内部实现：

if (relativePose_info.vec_inliers.size() <
        2.5 * KernelType::Solver::MINIMUM_SAMPLES )
    {
		
      return false; // no sufficient coverage (the model does not support enough samples)
    }
  }

3. Pose(317) –pose final(257)：！！！大部分情况下对pose-graph影响较大
   约束1：利用triplet-error-inference，去除部分posepair
   约束2：biedge 约束

5-1-1 Keep only the largest biedge connected subgraph

--去除与少数视图相关的所有点匹配

实现
1.

graph::CleanGraph_KeepLargestBiEdge_Nodes<Pair_Set, IndexT>(pairs)
先构造biedge connected components –romove edge
// Remove not bi-edge connected edges

根据largest subgraph —remove node（get remain nodes） // Graph is
bi-edge connected, but still many connected components can exist //
Keep only the nodes belonging to the largest Bi-edge component

（等价理解： 先去除弱连接（not
bi-edge）的边，得到含有多个连通分量的graph(集图本身不连通)，随后在保留graph中的最大子图，其为 largest
Bi-edge component ）

KeepOnlyReferencedElement(set_remainingIds, matches_provider_->pairWise_matches_);

3.

Geometry graph—view graph: biedge 约束 大部分情况下对pose-graph影响较小

4. Biedge-connected component graph (重连通分量、重连通子图)

参考链接1 参考链接2 参考链接3

定义：
在一个无向图中，若任意两点间至少存在两条“点不重复”的路径，则说这个图是点双连通的（简称双连通,biconnected）

在一个无向图中，点双连通的极大子图称为点双连通分量（简称双连通分量,Biconnected Component,BCC）
点双连通：删掉一个点之后，图仍联通 边双连通：删掉一条边之后，图仍联通

5-1-2 Compute_Relative_Rotations (2 views matching)

1. 方法步骤

E-R（四种可能解，cheirality 测试） E估计：(基于两视图)

初值：
Pose ：robust—5点法（小孔相机模型）或8点法 + ransac
Structure：利用像点匹配关系以及计算的r\t +intrinsics（已知）基于两视图的空三进行计算

BA调整(intrinsics作为常数，不参与迭代更新)

2. 接口函数：Relative_Pose_Engine::Relative_Pose_Engine::Process() :

（1）compute all the possible relative pose

**step1**

robustRelativePose(cam_I, cam_J,
                   x1, x2, relativePose_info,
                   {cam_I->w(), cam_I->h()},
                   {cam_J->w(), cam_J->h()},
                   256)
                   
robustRelativePose函数内部实现：

1、AContrario adaptor to use the 5 point essential matrix solver

！！！！
View graph---pose graph : 5点法 ransac 进行E估计时，model(E约束)内点数<2.5 minimun_sample(5),该view pair在pose graph中 舍去，大部分情况下对pose-graph影响较小

2、RelativePoseFromEssential(
    						bearing1,
    						bearing2,
    						relativePose_info.essential_matrix,
relativePose_info.vec_inliers, &relative_pose) 

RelativePoseFromEssential函数内部实现：
// Recover plausible relative poses from E.    
// Find which solution is the best:
// count how many triangulated observations are in front of the cameras 【Triangulate2View  ---- cheirality_accumulator】

！！！重要：
bearing vector: 方位向量，指的是像点在相机坐标系（局部三维坐标系）下的空间坐标！！！

**step2**
const bool bRefine_using_BA = true;  // Refine the defined scene  - refine only Structure and Rotations & translations (keep intrinsic constant)
BA实现时需要输入的参数：{initial-iteration-(refinement optical)}

（2）Get_Relative_Poses()

5-1-3 Compute_Global_Rotations(relatives_R, global_rotations)

接口函数： rotation_averaging_solver.Run( eRotation_averaging_method_, eRelativeRotationInferenceMethod, relatives_R, global_rotations)

（1） 计算global rotation时，对relative rotation graph 进行biedge约束，并整理，影响pose入网率（以实验8为例）

（2）Triplet inference (test over the composition error)

采用triplet-error-inference
Global rotation 的计算基于relative rotation, 具有“传递性”，
在relative rotation图中，采用不同的edge，对于同一pose，可能计算得到不同的Global rotation。因此采用triplet进行约束，对于满足条件（一定精度的edge才进行使用，用于后续计算）
条件：基于 Rij*Rjk*Rki vs I （identity matrics）的角度误差
//Rejection triplet that are 'not' identity rotation (error to identity > 5°)
具体实现如下：

getRotationMagnitude()：实现 rotation matrix -- 四元数 ---angle_axis 的转换，进一步计算角度误差

！！！angle_axis： 表示用轴角描述旋转时，对应的角度

此步骤对relative rotation 进行整理更新，影响最终的valid pose数目

Reject前：（以实验1为例）

Reject后：

（3）solve global rotation computation

//Solve the global rotation estimation problem:`
bSuccess = rotation_averaging::l2::L2RotationAveraging(reindexForward.size(),
        											   relativeRotations,
       												   vec_globalR);

//Non linear refinement of the global rotations
  if (bSuccess)
    bSuccess = rotation_averaging::l2::L2RotationAveraging_Refine(relativeRotations,
    															  vec_globalR);
Refine the global rotation:基于角度表示的rotation矩阵计算稀疏最小二乘法，具体由 class ceres::Problem 实现

5-1-4 Compute_Global_Translations(global_rotations, tripletWise_matches)

1. 方法原理

3视图-3焦点张量估计T –P1\P2\P3-P分解得到KRt
在计算T时，使用K && global rotation作为已知量，最后分解得到的t是relative translation（每个triplet单独计算，局部参考坐标系）
以triplets为单位，计算与每条“边”相关的relative motion,即同一边，有多个可能的relative translation

2. 步骤流程
   Poseid—graph –edge—triplet –visibility count per triplets (sum of their 2 view matches)-sort—计算Tensor(triplets)—分解得到t(relative translation\edge : t = 1:n)
3. 接口函数： translation_averaging_solver.Run( eTranslation_averaging_method_, sfm_data_, features_provider_, matches_provider_, global_rotations, tripletWise_matches);

（1） Compute the relative translations and save them to vec_initialRijTijEstimates:

Compute_translations(
    sfm_data,
    features_provider,
    matches_provider,
    map_globalR,
    tripletWise_matches);

Compute_translations（）函数实现:
    	
ComputePutativeTranslation_EdgesCoverage()
//Compute relative translations over the graph of global rotations--Perform a trifocal estimation of the graph contain in vec_triplets with an edge coverage algorithm
PS:Avoid to cover each edge of the graph by using an edge coverage algorithm !!!(An estimated triplets of translation mark three edges as estimated.)


原理：根据global rotation 步骤中得到的有效pose graph 计算每条edge的relative translation时，采用triplets 计算，在此过程中，该triplets所涉及的3条边均已计算（虽然是为了计算某一条边），因此无需重复遍历edge来计算

ComputePutativeTranslation_EdgesCoverage()函数内部：
	Step1   List plausible triplets over the global rotation pose graph Ids
	step2   Try to estimate this triplet of translations
接口函数：Estimate_T_triplet()// Robust estimation and refinement of a triplet of translations 

Estimate_T_triplet函数内部：
//(robust::ACRANSAC(kernel, vec_inliers, ORSA_ITER, &T, dPrecision/min_focal, false)  
//KRt_From_P(T.P1, &K, &R, &vec_tis[0]))

！！！重要：
//dPrecision/min_focal 单位统一：
//将dPrecision的像素单位转换为相机坐标系的mm单位（采用所有像片中focal最小值）
//focal表示相片f(mm)所对应的像素数



	step3   RelativeCameraMotion(RI, ti, RJ, tj, &Rij, &tij)

Ps: 可复用的常用函数、方法与模块：   ApplyTransformationToPoints(x1_, Kinv_, &x1n_);

（2）GLOBAL TRANSLATIONS ESTIMATION from initial triplets t_ij guess

1. 基于eTranslationAveragingMethod,结合Tj=Rij* Ti+λij *tij，计算唯一的global translation

2. 接口函数： Translation_averaging( eTranslationAveragingMethod, sfm_data, map_globalR)

Translation_averaging()内部实现：
step1 // Keep the largest Biedge connected component graph of relative translations
step2 //global translation

 	OSI_CLP_SolverWrapper solverLP(vec_solution.size())
	const bool bFeasible = solverLP.solve();// Solving

5-1-5 Compute_Initial_Structure(tripletWise_matches)

（1）Build tracks from selected triplets (Union of all the validated triplet tracks (_tripletWise_matches))

  TracksBuilder 
  tracksBuilder;

（2） Compute 3D position of the landmark of the structure by triangulation of the observations

接口函数： SfM_Data_Structure_Computation_Blind structure_estimator(true);

SfM_Data_Structure_Computation_Blind structure_estimator函数实现：
structure_estimator.triangulate(sfm_data_);
// Export initial structure  "initial_structure", "ply"

5-1-6 Adjust()

接口函数 bundle_adjustment_obj.Adjust()

3 processs
 
// - refine only Structure and translations    
"structure_00_refine_T_Xi", "ply"

// - refine only Structure and Rotations & translations  
"structure_01_refine_RT_Xi", "ply" 

//- refine all: Structure, motion:{rotations, translations} and optics:{intrinsics}  "structure_02_refine_KRT_Xi", "ply"

5-1-7 Remove outliers (max_angle, residual error)

RemoveOutliers_PixelResidualError(sfm_data_, 4.0)
RemoveOutliers_AngleError(sfm_data_, 2.0)
"structure_03_outlier_removed", "ply"

eraseUnstablePosesAndObservations(sfm_data_, minPointPerPose, minTrackLength)

5-1-8 Final BA. We refine one more time,

 // since some outlier have been removed and so a better solution can be found.
"structure_04_outlier_removed", "ply"

5-1-9 sfmEngine.Get_SfM_Data()

"sfm_data", ".bin"
"cloud_and_poses", ".ply"

6、pose graph 可视化

'neato' 不是内部或外部命令，也不是可运行的程序

不影响整个解算过程，但无法将pose graph 可视化

解决办法：参考链接
下载安装 graphviz后 ，并配置环境变量，命令行运行
dot -version \ dot -c
确保成功安装
官网安装程序执行后， bin文件夹下，无neato.exe？？？
采用其他安装程序，安装旧版本 graphviz