《视觉14讲》Eigen库的调用实例。

灰太狼 2022-12-23 09:49 238阅读 0赞

### CMakeLists.txt ###

# 添加Eigen头文件
    include_directories( "/usr/include/eigen3" )

### useEigen.cpp ###

#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> 
                ;
    
        // 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;    // 行列式
        int a = 10;
    
        // 特征值
        // 实对称矩阵可以保证对角化成功
        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;
    }

### useGeometry.cpp ###

#include <iostream>
    #include <cmath>
    
    using namespace std;
    #include <Eigen/Core>
    #include <Eigen/Geometry>
    
    int main(int argc,char** argv)
    {
        
     // 3D 旋转矩阵
      Eigen::Matrix3d rotation_matrix = Eigen::Matrix3d::Identity();
      //旋转向量使用AngleAxis ,它的底层不是Matrix但是可以当做矩阵来运算
      Eigen::AngleAxisd rotation_vector (M_PI/4,Eigen::Vector3d(0,0,1)); //沿着z轴旋转Pi/2
      cout.precision(3);
      cout<<"rotation_matrix = \n"<<rotation_vector.matrix()<<endl; //用 matrix()转化成矩阵
      rotation_matrix = rotation_vector.toRotationMatrix();
      cout<<"rotation_matrix = \n"<<rotation_matrix<<endl;
      
      //用AngleAxis 进行坐标转换
      Eigen::Vector3d v(1,0,0);
      
      //公式 旋转矩阵 乘 V
      Eigen::Vector3d v_rotation = rotation_vector*v;
      
      cout<<"1,0,0 after rotation = \n"<<v_rotation.transpose()<<endl;  // 向量v(1,0,0)也就是三维坐标系中的x轴绕着z轴旋转了pi/2°之后的旋转向量v'
      v_rotation = rotation_matrix  * v; //3D旋转矩阵rotation_vector可以当做旋转矩阵rotation_matrix来用，都是3X3的矩阵
      cout<<"v_rotation = \n"<<v_rotation<<endl;
      
      //将旋转矩阵直接转化成欧拉角
      Eigen::Vector3d euler_angles = rotation_matrix.eulerAngles(2,1,0);
      cout.precision(3);
      cout<<"转化后的欧拉角是\n"<<euler_angles.transpose()<<endl; //也就是yaw pitch roll的顺序
      
      //欧式变换使用 Eigen::Isometry
    Eigen::Isometry3d T = Eigen::Isometry3d::Identity();
    T.rotate(rotation_vector);
    T.pretranslate(Eigen::Vector3d(1,3,4));
    cout<<T.matrix()<<endl;
    
    //用变换矩阵进行坐标转换
    Eigen::Vector3d v_transformed = T*v; //相当于R×v + T  v是（1,0,0） 相当于对向量(1,0,0) 进行旋转与平移操作 具体是 向量(1,0,0)绕着z轴旋转pi/2,然后进行平移，x方向平移1，y方向平移3，z方向平移4
    cout<<" v transformed = "<< v_transformed.transpose()<<endl;
    
    //四元数
    Eigen::Quaterniond q = Eigen::Quaterniond(rotation_vector);
    cout<<"Quaternion is \n"<<q.coeffs()<<endl; //顺序分别是x y z w,w为实部，前三者为虚部
    //除了可以吧AngleAxis 赋给 q （四元数），还可以是矩阵（rotation_matrix）
    q = Eigen::Quaterniond(rotation_matrix);
    cout<<"系数是\n"<<q.coeffs()<<endl;
    
    //用四元数来旋转一个向量，正确的姿势是
    v_rotation = q*v; //数学上实际是 qvq^{-1}
    cout<<"(1,0,0) 被旋转之后是"<<v_rotation.transpose()<<endl;
      return 0;
    }

**创建随机矩阵**

Eigen::Matrix<double,4,5> matrix_45;
    matrix_45 = Eigen::Matrix<double,4,5>::Random();