视觉SLAM十四讲作业练习（2）

花舞の街道

已于 2022-05-03 19:38:45 修改

阅读量622

点赞数 1

分类专栏：视觉SLAM 文章标签：线性代数矩阵几何学

于 2022-05-01 20:24:20 首次发布

本文链接：https://blog.csdn.net/ASUNAchan/article/details/124531937

版权

视觉SLAM 专栏收录该内容

18 篇文章 9 订阅

订阅专栏

2.Eigen 矩阵运算

设线性方程 Ax = b，在 A 为方阵的前提下，请回答以下问题：

在什么条件下，x 有解且唯一？

有解的充分必要条件是：系数矩阵的秩等于增广矩阵的秩，即rank(A)=rank(A, b)

高斯消元法的原理是什么?

逐步消除未知数来将原始线性系统转化为另一个更简单的等价的系统。用初等行变换把增广矩阵转换为行阶梯阵，然后回代求出方程的解。

QR 分解的原理是什么？
Cholesky 分解的原理是什么？
编程实现 A 为 100 × 100 随机矩阵时，用 QR 和 Cholesky 分解求 x 的程序。你可以参考本次课用到的 useEigen 例程。

利用Eigen求解非齐次方程组

1.编写主程序

#include <iostream>
using namespace std;
#include <ctime>
// Eigen 部分
#include <Eigen/Core>
// 稠密矩阵的代数运算（逆，特征值等）
#include <Eigen/Dense>

#define MATRIX_SIZE 100

int main(int argc, char *argv[])
{
    //求解Ax=b 
    //A为100*100随机矩阵，b为1*100的列向量
    //A：matrix_NN
    Eigen::Matrix< double, MATRIX_SIZE, MATRIX_SIZE > matrix_NN;
    matrix_NN = Eigen::MatrixXd::Random( MATRIX_SIZE, MATRIX_SIZE );
    //b:v_Nd
    Eigen::Matrix< double, MATRIX_SIZE,  1> v_Nd;
    v_Nd = Eigen::MatrixXd::Random( MATRIX_SIZE,1 );

    Eigen::Matrix< double, MATRIX_SIZE, 1 > x_inverse;
    Eigen::Matrix< double, MATRIX_SIZE, 1 > x_QR;
    //直接求逆
    x_inverse = matrix_NN.inverse() * v_Nd;
    //QR分解
    x_QR = matrix_NN.colPivHouseholderQr().solve(v_Nd);

    cout << "A矩阵 matrix_NN : " << endl; 
    cout << matrix_NN << endl;
    cout << "b向量 v_Nd : " << endl; 
    cout << v_Nd << endl;
    cout << "通过求逆所得的x x_inverse : " << endl; 
    cout << x_inverse << endl;
    cout << "通过QR分解所得的x x_QR : " << endl; 
    cout << x_QR << endl;

    return 0;
}

2.编写CMakeLists.txt文件

cmake_minimum_required(VERSION 3.0)
project(useEigen)
include_directories(/usr/include/eigen3)
add_executable(useEigen useEigen.cpp)

因为Eigen3仅包含头文件，所以只需添加include_directories(/usr/include/eigen3)

3.例程

#include <iostream>
using namespace std;
#include <ctime>
// Eigen 部分
#include <Eigen/Core>
// 稠密矩阵的代数运算（逆，特征值等）
#include <Eigen/Dense>

#define MATRIX_SIZE 50

/****************************
* 本程序演示了 Eigen 基本类型的使用
****************************/

int main( int argc, char** argv )
{
    // Eigen 中所有向量和矩阵都是Eigen::Matrix，它是一个模板类。它的前三个参数为：数据类型，行，列
    // 声明一个2*3的float矩阵
    Eigen::Matrix<float, 2, 3> matrix_23;

    // 同时，Eigen 通过 typedef 提供了许多内置类型，不过底层仍是Eigen::Matrix
    // 例如 Vector3d 实质上是 Eigen::Matrix<double, 3, 1>，即三维向量
    Eigen::Vector3d v_3d;
	// 这是一样的
    Eigen::Matrix<float,3,1> vd_3d;

    // Matrix3d 实质上是 Eigen::Matrix<double, 3, 3>
    Eigen::Matrix3d matrix_33 = Eigen::Matrix3d::Zero(); //初始化为零
    // 如果不确定矩阵大小，可以使用动态大小的矩阵
    Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > matrix_dynamic;
    // 更简单的
    Eigen::MatrixXd matrix_x;
    // 这种类型还有很多，我们不一一列举

    // 下面是对Eigen阵的操作
    // 输入数据（初始化）
    matrix_23 << 1, 2, 3, 4, 5, 6;
    // 输出
    cout << matrix_23 << endl;

    // 用()访问矩阵中的元素
    for (int i=0; i<2; i++) {
        for (int j=0; j<3; j++)
            cout<<matrix_23(i,j)<<"\t";
        cout<<endl;
    }

    // 矩阵和向量相乘（实际上仍是矩阵和矩阵）
    v_3d << 3, 2, 1;
    vd_3d << 4,5,6;
    // 但是在Eigen里你不能混合两种不同类型的矩阵，像这样是错的
    // Eigen::Matrix<double, 2, 1> result_wrong_type = matrix_23 * v_3d;
    // 应该显式转换
    Eigen::Matrix<double, 2, 1> result = matrix_23.cast<double>() * v_3d;
    cout << result << endl;

    Eigen::Matrix<float, 2, 1> result2 = matrix_23 * vd_3d;
    cout << result2 << endl;

    // 同样你不能搞错矩阵的维度
    // 试着取消下面的注释，看看Eigen会报什么错
    // Eigen::Matrix<double, 2, 3> result_wrong_dimension = matrix_23.cast<double>() * v_3d;

    // 一些矩阵运算
    // 四则运算就不演示了，直接用+-*/即可。
    matrix_33 = Eigen::Matrix3d::Random();      // 随机数矩阵
    cout << matrix_33 << endl << endl;

    cout << matrix_33.transpose() << endl;      // 转置
    cout << matrix_33.sum() << endl;            // 各元素和
    cout << matrix_33.trace() << endl;          // 迹
    cout << 10*matrix_33 << endl;               // 数乘
    cout << matrix_33.inverse() << endl;        // 逆
    cout << matrix_33.determinant() << endl;    // 行列式

    // 特征值
    // 实对称矩阵可以保证对角化成功
    Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> eigen_solver ( matrix_33.transpose()*matrix_33 );
    cout << "Eigen values = \n" << eigen_solver.eigenvalues() << endl;
    cout << "Eigen vectors = \n" << eigen_solver.eigenvectors() << endl;

    // 解方程
    // 我们求解 matrix_NN * x = v_Nd 这个方程
    // N的大小在前边的宏里定义，它由随机数生成
    // 直接求逆自然是最直接的，但是求逆运算量大

    Eigen::Matrix< double, MATRIX_SIZE, MATRIX_SIZE > matrix_NN;
    matrix_NN = Eigen::MatrixXd::Random( MATRIX_SIZE, MATRIX_SIZE );
    Eigen::Matrix< double, MATRIX_SIZE,  1> v_Nd;
    v_Nd = Eigen::MatrixXd::Random( MATRIX_SIZE,1 );

    clock_t time_stt = clock(); // 计时
    // 直接求逆
    Eigen::Matrix<double,MATRIX_SIZE,1> x = matrix_NN.inverse()*v_Nd;
    cout <<"time use in normal inverse is " << 1000* (clock() - time_stt)/(double)CLOCKS_PER_SEC << "ms"<< endl;
    
	// 通常用矩阵分解来求，例如QR分解，速度会快很多
    time_stt = clock();
    x = matrix_NN.colPivHouseholderQr().solve(v_Nd);
    cout <<"time use in Qr decomposition is " <<1000*  (clock() - time_stt)/(double)CLOCKS_PER_SEC <<"ms" << endl;

    return 0;
}

4.几何运算练习

下面我们来练习如何使用 Eigen/Geometry 计算一个具体的例子。

一个机器人上通常会安装许多不同的传感器，而且这些传感器之间还存在固连关系。我们举一个典型的例子。在世界系 W 下，存在一个运动的机器人 R。按照固定的或者某些开发人员或者领导的特殊喜好，R 系定义在机器人脚部的位置。但是机器人在设计的时候，又定义了 B 系（Body 系，或本体系），位于机器人头部的位置。由于沟通不畅，标定人员把一台激光传感器和一台视觉传感器标定在了 B 系下。我们称激光传感器为 L 系，视觉传感器为 C 系。

现在请你完成以下工作：

说明一个激光传感器下的看到的点应该如何计算它的世界坐标。

在雷达坐标下的点p_l，需要左乘T_bl将p_l转换到p_b, 再需要左乘T_rb转换到p_r, 最后再左乘T_wr转换到p_w， p_w = T_wr * T_rb * T_bl * p_l

取：qWR = [0.55, 0.3, 0.2, 0.2], tW R = [0.1, 0.2, 0.3]T ，qRB = [0.99, 0, 0, 0.01], tRB = [0.05, 0, 0.5]T ， qBL = [0.3, 0.5, 0, 20.1], tBL = [0.4, 0, 0.5]T , qBC = [0.8, 0.2, 0.1, 0.1], tBC = [0.5, 0.1, 0.5]T 。现在假设相机传感器观察到自身坐标系下的点 [0.3, 0.2, 1.2]，请计算：

(a)这个点在激光系下的坐标；

(b)这个点在世界系下的坐标。

提示： 1. 本题的数据是随意取的，没有真实意义，产生什么结果都不要奇怪。 2. 四元数在使用前需要归一化。 3. 请注意 Eigen 在使用四元数时的虚部和实部顺序。

useGeometry.cpp

#include <iostream>
#include <Eigen/Core>
// Eigen 几何模块
#include <Eigen/Geometry>

using namespace std;

int main(int argc, char *argv[])
{
    Eigen::Quaterniond q_wr(0.55, 0.3, 0.2, 0.2);
    Eigen::Quaterniond q_rb(0.99, 0.0, 0.0, 0.01);
    Eigen::Quaterniond q_bl(0.3, 0.5, 0.0, 20.1);
    Eigen::Quaterniond q_bc(0.8, 0.2, 0.1, 0.1);

    //四元数归一化
    q_wr.normalize();
    q_rb.normalize();
    q_bl.normalize();
    q_bc.normalize();

    Eigen::Vector3d t_wr(0.1, 0.2, 0.3);
    Eigen::Vector3d t_rb(0.05, 0, 0.5);
    Eigen::Vector3d t_bl(0.4, 0, 0.5);
    Eigen::Vector3d t_bc(0.5, 0.1, 0.5);

    Eigen::Vector3d p_c(0.3, 0.2, 1.2);

    Eigen::Isometry3d T_wr(q_wr);
    T_wr.pretranslate(t_wr);
    Eigen::Isometry3d T_rb(q_rb);
    T_rb.pretranslate(t_rb);
    Eigen::Isometry3d T_bl(q_bl);
    T_bl.pretranslate(t_bl);
    Eigen::Isometry3d T_bc(q_bc);
    T_bc.pretranslate(t_bc);

    Eigen::Vector3d p_l;
    Eigen::Vector3d p_w;

    p_l = T_bl.inverse() * T_bc * p_c;
    p_w = T_wr * T_rb * T_bc * p_c;

    cout << "p_c在相机坐标系下的坐标: " ;
    cout << p_c.transpose() << endl; 

    cout << "p_c在激光坐标系下的坐标: " ;
    cout << p_l.transpose() << endl;

    cout << "p_c在世界坐标系下的坐标: " ;
    cout << p_w.transpose() << endl;
    

    return 0;
}

CMakeLists.txt

cmake_minimum_required(VERSION 3.0)
project(geometry)
set(CMAKE_BUILD_TYPE Release)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++14 -Wall -O3")
include_directories(/usr/include/eigen3)
add_executable(useGeometry useGeometry.cpp)

花舞の街道

关注

1
点赞
踩
5

收藏

觉得还不错? 一键收藏
0
评论
视觉SLAM十四讲作业练习（2）

2.Eigen 矩阵运算设线性方程 Ax = b，在 A 为方阵的前提下，请回答以下问题：在什么条件下，x 有解且唯一？有解的充分必要条件是：系数矩阵的秩等于增广矩阵的秩，即rank(A)=rank(A, b)高斯消元法的原理是什么?逐步消除未知数来将原始线性系统转化为另一个更简单的等价的系统。用初等行变换把增广矩阵转换为行阶梯阵，然后回代求出方程的解。QR 分解的原理是什么？Cholesky 分解的原理是什么？编程实现 A 为 100 × 100 随机矩阵时，用 QR 和 C
复制链接

扫一扫