0. 内容
T1.
参考SLAM14讲P247直接可写,注意 ξ 1 , ξ 2 \xi_1,\xi_2 ξ1,ξ2之间有约束(关系)。
套用舒尔补公式:
marg掉 ξ 1 \xi_1 ξ1之后,信息被传递到 L 1 和 L 2 L_1和L_2 L1和L2之间了。
T2.
T3.
课上同学的作业分享:
算是勉强看着答案做出来了:
1.
我的维度不对 ,H的分布就是先pose,再路标点,所以先计算完pose的Jacobian才能计算landmark的Jacobian
那个Jacobian看不懂怎么算的。
- 整体代码
//
// Created by hyj on 18-11-11.
//
#include <iostream>
#include <vector>
#include <random>
#include <Eigen/Core>
#include <Eigen/Geometry>
#include <Eigen/Eigenvalues>struct Pose
{Pose(Eigen::Matrix3d R, Eigen::Vector3d t):Rwc(R),qwc(R),twc(t) {};Eigen::Matrix3d Rwc;Eigen::Quaterniond qwc;Eigen::Vector3d twc;
};
int main()
{int featureNums = 20;int poseNums = 10;int diem = poseNums * 6 + featureNums * 3;double fx = 1.;double fy = 1.;Eigen::MatrixXd H(diem,diem);H.setZero();std::vector<Pose> camera_pose;double radius = 8;for(int n = 0; n < poseNums; ++n ) {double theta = n * 2 * M_PI / ( poseNums * 4); // 1/4 圆弧// 绕 z轴 旋转Eigen::Matrix3d R;R = Eigen::AngleAxisd(theta, Eigen::Vector3d::UnitZ());Eigen::Vector3d t = Eigen::Vector3d(radius * cos(theta) - radius, radius * sin(theta), 1 * sin(2 * theta));camera_pose.push_back(Pose(R,t));}// 随机数生成三维特征点std::default_random_engine generator;std::vector<Eigen::Vector3d> points;for(int j = 0; j < featureNums; ++j){std::uniform_real_distribution<double> xy_rand(-4, 4.0);std::uniform_real_distribution<double> z_rand(8., 10.);double tx = xy_rand(generator);double ty = xy_rand(generator);double tz = z_rand(generator);Eigen::Vector3d Pw(tx, ty, tz);points.push_back(Pw);for (int i = 0; i < poseNums; ++i) {Eigen::Matrix3d Rcw = camera_pose[i].Rwc.transpose();Eigen::Vector3d Pc = Rcw * (Pw - camera_pose[i].twc);double x = Pc.x();double y = Pc.y();double z = Pc.z();double z_2 = z * z;Eigen::Matrix<double,2,3> jacobian_uv_Pc;// Jacobian of residual with respect to point in camera coordinatejacobian_uv_Pc<< fx/z, 0 , -x * fx/z_2,0, fy/z, -y * fy/z_2;Eigen::Matrix<double,2,3> jacobian_Pj = jacobian_uv_Pc * Rcw;Eigen::Matrix<double,2,6> jacobian_Ti;jacobian_Ti << -x* y * fx/z_2, (1+ x*x/z_2)*fx, -y/z*fx, fx/z, 0 , -x * fx/z_2,-(1+y*y/z_2)*fy, x*y/z_2 * fy, x/z * fy, 0,fy/z, -y * fy/z_2;H.block(i*6,i*6,6,6) += jacobian_Ti.transpose() * jacobian_Ti;//guan yu pose de Jacobian/// 请补充完整作业信息矩阵块的计算//zhu dui jiao xian 3*3
// H.block(i * 6 + j * 3,i * 6 + j * 3,3,3) += jacobian_Pj.transpose() * jacobian_Pj; //mine
// H.block(i * 6 ,i * 6 + j * 3,6,3) += jacobian_Ti.transpose() * jacobian_Pj;
// H.block(i * 6 + j * 3, i * 6 ,3,6) += jacobian_Pj.transpose() * jacobian_Ti;//zheng que deH.block(poseNums * 6 + j * 3,poseNums * 6 + j * 3,3,3) += jacobian_Pj.transpose() * jacobian_Pj;H.block(i * 6 ,poseNums * 6 + j * 3,6,3) += jacobian_Ti.transpose() * jacobian_Pj;H.block(poseNums * 6 + j * 3, i * 6 ,3,6) += jacobian_Pj.transpose() * jacobian_Ti;}}// std::cout << H << std::endl;
// Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> saes(H);
// std::cout << saes.eigenvalues() <<std::endl;Eigen::JacobiSVD<Eigen::MatrixXd> svd(H, Eigen::ComputeThinU | Eigen::ComputeThinV);std::cout << svd.singularValues() <<std::endl;return 0;
}
最后的SVD分解还不是很熟,矩阵论忘差不多了。