一、流程
1. 三角化
估计 features 初始深度
2. 优化求解
构造非线性最小二乘问题
3. 边缘化
求解边缘化先验约束
二、三角化
1. 细节
(1). 用所有帧的 pose 三角化,待求解量:启始帧相机坐标系下的齐次坐标
(1)
(2). 通过相机坐标系下的齐次坐标,计算深度
(2)
2. 代码实现
void FeatureManager::triangulate(Vector3d Ps[], Vector3d tic[], Matrix3d ric[])
{
for (auto &it_per_id : feature) {
......
int imu_i = it_per_id.start_frame, imu_j = imu_i - 1; // 起始帧
Eigen::MatrixXd svd_A(2 * it_per_id.feature_per_frame.size(), 4);
int svd_idx = 0;
Eigen::Matrix<double, 3, 4> P0;
Eigen::Vector3d t0 = Ps[imu_i] + Rs[imu_i] * tic[0];
Eigen::Matrix3d R0 = Rs[imu_i] * ric[0];
P0.leftCols<3>() = Eigen::Matrix3d::Identity();
P0.rightCols<1>() = Eigen::Vector3d::Zero();
// 用所有帧的 pose 三角化
for (auto &it_per_frame : it_per_id.feature_per_frame) {
imu_j++;
Eigen::Vector3d t1 = Ps[imu_j] + Rs[imu_j] * tic[0];
Eigen::Matrix3d R1 = Rs[imu_j] * ric[0];
Eigen::Vector3d t = R0.transpose() * (t1 - t0);
Eigen::Matrix3d R = R0.transpose() * R1;
Eigen::Matrix<double, 3, 4> P;
P.leftCols<3>() = R.transpose();
P.rightCols<1>() = -R.transpose() * t;
Eigen::Vector3d f = it_per_frame.point.normalized();
svd_A.row(svd_idx++) = f[0] * P.row(2) - f[2] * P.row(0);
svd_A.row(svd_idx++) = f[1] * P.row(2) - f[2] * P.row(1);
if (imu_i == imu_j)
continue;
}
// SVD 求解
Eigen::Vector4d svd_V = Eigen::JacobiSVD<Eigen::MatrixXd>(svd_A, Eigen::ComputeThinV).matrixV().rightCols<1>();
double svd_method = svd_V[2] / svd_V[3];
// 估计初始深度
it_per_id.estimated_depth = svd_method;
if (it_per_id.estimated_depth < 0.1) {
it_per_id.estimated_depth = INIT_DEPTH;
}
}
}
三、优化
1. marginlization_factor
引入边缘化先验约束
if (last_marginalization_info){
// construct new marginlization_factor
MarginalizationFactor *marginalization_factor = new MarginalizationFactor(last_marginalization_info);
problem.AddResidualBlock(marginalization_factor, NULL,
last_marginalization_parameter_blocks);
}
具体细节 todo
2. imu factor
for (int i = 0; i < WINDOW_SIZE; i++) {
int j = i + 1;
if (pre_integrations[j]->sum_dt > 10.0)
continue;
IMUFactor* imu_factor = new IMUFactor(pre_integrations[j]);
problem.AddResidualBlock(imu_factor, NULL, para_Pose[i], para_SpeedBias[i], para_Pose[j], para_SpeedBias[j]);
}
具体细节参考:Vins-Mono 论文 && Coding 一 4(2). imu factor
3. reprojection factor
int f_m_cnt = 0;
int feature_index = -1;
for (auto &it_per_id : f_manager.feature) {
it_per_id.used_num = it_per_id.feature_per_frame.size();
if (!(it_per_id.used_num >= 2 && it_per_id.start_frame < WINDOW_SIZE - 2))
continue;
++feature_index;
int imu_i = it_per_id.start_frame;
int imu_j = imu_i - 1;
Vector3d pts_i = it_per_id.feature_per_frame[0].point;
for (auto &it_per_frame : it_per_id.feature_per_frame) {
imu_j++;
if (imu_i == imu_j) {
continue;
}
Vector3d pts_j = it_per_frame.point;
if (ESTIMATE_TD) { // 考虑 timeoffset
ProjectionTdFactor *f_td = new ProjectionTdFactor(pts_i, pts_j, it_per_id.feature_per_frame[0].velocity, it_per_frame.velocity, it_per_id.feature_per_frame[0].cur_td, it_per_frame.cur_td, it_per_id.feature_per_frame[0].uv.y(), it_per_frame.uv.y());
problem.AddResidualBlock(f_td, loss_function, para_Pose[imu_i], para_Pose[imu_j], para_Ex_Pose[0], para_Feature[feature_index], para_Td[0]);
} else {
ProjectionFactor *f = new ProjectionFactor(pts_i, pts_j); // 观测:在两帧之间的像素坐标
problem.AddResidualBlock(f, loss_function, para_Pose[imu_i], para_Pose[imu_j], para_Ex_Pose[0], para_Feature[feature_index]); // 影响的状态变量:第 i,j 帧的 pose,外参,特征点位置信息
}
f_m_cnt++; // factor 计数加上一
}
}
4. relocalization factor
if(relocalization_info) {
ceres::LocalParameterization *local_parameterization = new PoseLocalParameterization();
problem.AddParameterBlock(relo_Pose, SIZE_POSE, local_parameterization); // 引入的额外状态变量:double relo_Pose[SIZE_POSE];
int retrive_feature_index = 0;
int feature_index = -1;
for (auto &it_per_id : f_manager.feature) {
it_per_id.used_num = it_per_id.feature_per_frame.size();
if (!(it_per_id.used_num >= 2 && it_per_id.start_frame < WINDOW_SIZE - 2))
continue;
++feature_index;
int start = it_per_id.start_frame;
if(start <= relo_frame_local_index) {
while((int)match_points[retrive_feature_index].z() < it_per_id.feature_id) {
retrive_feature_index++;
}
if((int)match_points[retrive_feature_index].z() == it_per_id.feature_id) {
// 该 feature 同时在 frame_i(窗口内) 和 relo_pose(窗口外) 被观测
Vector3d pts_j = Vector3d(match_points[retrive_feature_index].x(), match_points[retrive_feature_index].y(), 1.0);
Vector3d pts_i = it_per_id.feature_per_frame[0].point;
ProjectionFactor *f = new ProjectionFactor(pts_i, pts_j);
problem.AddResidualBlock(f, loss_function, para_Pose[start], relo_Pose, para_Ex_Pose[0], para_Feature[feature_index]);
retrive_feature_index++;
}
}
}
}
具体细节 todo
三、marginalization
1. marginalization_flag == MARGIN_OLD
<1>. last margin block
ResidualBlockInfo *residual_block_info = new
ResidualBlockInfo(marginalization_factor, NULL,
last_marginalization_parameter_blocks, drop_set);
marginalization_info->addResidualBlockInfo(residual_block_info);
<2>. imu factor block
ResidualBlockInfo *residual_block_info = new ResidualBlockInfo(imu_factor, NULL, vector<double *>{para_Pose[0], para_SpeedBias[0], para_Pose[1], para_SpeedBias[1]}, vector<int>{0, 1});
marginalization_info->addResidualBlockInfo(residual_block_info);
<3>. reprojection error
if (ESTIMATE_TD) {
ProjectionTdFactor *f_td = new ProjectionTdFactor(pts_i, pts_j, it_per_id.feature_per_frame[0].velocity, it_per_frame.velocity, it_per_id.feature_per_frame[0].cur_td, it_per_frame.cur_td, it_per_id.feature_per_frame[0].uv.y(), it_per_frame.uv.y());
ResidualBlockInfo *residual_block_info = new ResidualBlockInfo(f_td, loss_function, vector<double *>{para_Pose[imu_i], para_Pose[imu_j], para_Ex_Pose[0], para_Feature[feature_index], para_Td[0]}, vector<int>{0, 3});
marginalization_info->addResidualBlockInfo(residual_block_info);
} else {
ProjectionFactor *f = new ProjectionFactor(pts_i, pts_j);
ResidualBlockInfo *residual_block_info = new ResidualBlockInfo(f, loss_function, vector<double *>{para_Pose[imu_i], para_Pose[imu_j], para_Ex_Pose[0], para_Feature[feature_index]}, vector<int>{0, 3});
marginalization_info->addResidualBlockInfo(residual_block_info);
}
<4>. marginalization_info->preMarginalize();
<5>. marginalization_info->marginalize();
2. marginalization_flag == MARGIN_SECOND_NEW
<1>. last margin block
ResidualBlockInfo *residual_block_info = new
ResidualBlockInfo(marginalization_factor, NULL,
last_marginalization_parameter_blocks, drop_set);
marginalization_info->addResidualBlockInfo(residual_block_info);
<2>. marginalization_info->preMarginalize();
<3>. marginalization_info->marginalize();
3. 不同 marginalization_flag 的区别
(1). 当次新帧为关键帧时,marginalization_flag 是 MARGIN_OLD,将 marg 掉最老帧
(2). 当次新帧不是关键帧时,marginalization_flag 是 MARGIN_SECOND_NEW, 我们将直接扔掉次新帧及它的视觉约束,而不对次新帧进行 marg