first version

include/Exp_mat.h

@@ -0,0 +1,103 @@
+#ifndef EXP_MAT_H
+#define EXP_MAT_H
+
+#include <math.h>
+#include <Eigen/Core>
+#include <opencv/cv.h>
+// #include <common_lib.h>
+
+#define SKEW_SYM_MATRX(v) 0.0,-v[2],v[1],v[2],0.0,-v[0],-v[1],v[0],0.0
+
+template<typename T>
+Eigen::Matrix<T, 3, 3> Exp(const Eigen::Matrix<T, 3, 1> &&ang)
+{
+    T ang_norm = ang.norm();
+    Eigen::Matrix<T, 3, 3> Eye3 = Eigen::Matrix<T, 3, 3>::Identity();
+    if (ang_norm > 0.0000001)
+    {
+        Eigen::Matrix<T, 3, 1> r_axis = ang / ang_norm;
+        Eigen::Matrix<T, 3, 3> K;
+        K << SKEW_SYM_MATRX(r_axis);
+        /// Roderigous Tranformation
+        return Eye3 + std::sin(ang_norm) * K + (1.0 - std::cos(ang_norm)) * K * K;
+    }
+    else
+    {
+        return Eye3;
+    }
+}
+
+template<typename T, typename Ts>
+Eigen::Matrix<T, 3, 3> Exp(const Eigen::Matrix<T, 3, 1> &ang_vel, const Ts &dt)
+{
+    T ang_vel_norm = ang_vel.norm();
+    Eigen::Matrix<T, 3, 3> Eye3 = Eigen::Matrix<T, 3, 3>::Identity();
+
+    if (ang_vel_norm > 0.0000001)
+    {
+        Eigen::Matrix<T, 3, 1> r_axis = ang_vel / ang_vel_norm;
+        Eigen::Matrix<T, 3, 3> K;
+
+        K << SKEW_SYM_MATRX(r_axis);
+
+        T r_ang = ang_vel_norm * dt;
+
+        /// Roderigous Tranformation
+        return Eye3 + std::sin(r_ang) * K + (1.0 - std::cos(r_ang)) * K * K;
+    }
+    else
+    {
+        return Eye3;
+    }
+}
+
+template<typename T>
+Eigen::Matrix<T, 3, 3> Exp(const T &v1, const T &v2, const T &v3)
+{
+    T &&norm = sqrt(v1 * v1 + v2 * v2 + v3 * v3);
+    Eigen::Matrix<T, 3, 3> Eye3 = Eigen::Matrix<T, 3, 3>::Identity();
+    if (norm > 0.00001)
+    {
+        T r_ang[3] = {v1 / norm, v2 / norm, v3 / norm};
+        Eigen::Matrix<T, 3, 3> K;
+        K << SKEW_SYM_MATRX(r_ang);
+
+        /// Roderigous Tranformation
+        return Eye3 + std::sin(norm) * K + (1.0 - std::cos(norm)) * K * K;
+    }
+    else
+    {
+        return Eye3;
+    }
+}
+
+/* Logrithm of a Rotation Matrix */
+template<typename T>
+Eigen::Matrix<T,3,1> Log(const Eigen::Matrix<T, 3, 3> &R)
+{
+    T &&theta = std::acos(0.5 * (R.trace() - 1));
+    Eigen::Matrix<T,3,1> K(R(2,1) - R(1,2), R(0,2) - R(2,0), R(1,0) - R(0,1));
+    return (std::abs(theta) < 0.001) ? (0.5 * K) : (0.5 * theta / std::sin(theta) * K);
+}
+
+// template<typename T>
+// cv::Mat Exp(const T &v1, const T &v2, const T &v3)
+// {
+    
+//     T norm = sqrt(v1 * v1 + v2 * v2 + v3 * v3);
+//     cv::Mat Eye3 = cv::Mat::eye(3, 3, CV_32F);
+//     if (norm > 0.0000001)
+//     {
+//         T r_ang[3] = {v1 / norm, v2 / norm, v3 / norm};
+//         cv::Mat K = (cv::Mat_<T>(3,3) << SKEW_SYM_MATRX(r_ang));
+
+//         /// Roderigous Tranformation
+//         return Eye3 + std::sin(norm) * K + (1.0 - std::cos(norm)) * K * K;
+//     }
+//     else
+//     {
+//         return Eye3;
+//     }
+// }
+
+#endif

include/common_lib.h

@@ -0,0 +1,140 @@
+#ifndef COMMON_LIB_H
+#define COMMON_LIB_H
+
+#include <Exp_mat.h>
+#include <Eigen/Eigen>
+#include <pcl/point_types.h>
+#include <pcl/point_cloud.h>
+#include <fast_lio/States.h>
+#include <fast_lio/Pose6D.h>
+#include <sensor_msgs/Imu.h>
+#include <nav_msgs/Odometry.h>
+#include <tf/transform_broadcaster.h>
+#include <eigen_conversions/eigen_msg.h>
+
+
+// #define DEBUG_PRINT
+
+#define PI_M (3.14159265358)
+#define G_m_s2 (9.8099)         // Gravaty const in GuangDong/China
+#define DIM_OF_STATES (18)      // Dimension of states (Let Dim(SO(3)) = 3)
+#define DIM_OF_PROC_N (12)      // Dimension of process noise (Let Dim(SO(3)) = 3)
+#define CUBE_LEN  (6.0)
+#define LIDAR_SP_LEN    (2)
+#define INIT_COV  (0.0001)
+
+#define VEC_FROM_ARRAY(v)        v[0],v[1],v[2]
+#define MAT_FROM_ARRAY(v)        v[0],v[1],v[2],v[3],v[4],v[5],v[6],v[7],v[8]
+#define CONSTRAIN(v,min,max)     ((v>min)?((v<max)?v:max):min)
+#define ARRAY_FROM_EIGEN(mat)    mat.data(), mat.data() + mat.rows() * mat.cols()
+#define STD_VEC_FROM_EIGEN(mat)  std::vector<decltype(mat)::Scalar> (mat.data(), mat.data() + mat.rows() * mat.cols())
+
+#define DEBUG_FILE_DIR(name)  (std::string(std::string(ROOT_DIR) + "Log/"+ name))
+
+typedef fast_lio::Pose6D Pose6D;
+typedef pcl::PointXYZINormal PointType;
+typedef pcl::PointCloud<PointType> PointCloudXYZI;
+
+Eigen::Matrix3d Eye3d(Eigen::Matrix3d::Identity());
+Eigen::Matrix3f Eye3f(Eigen::Matrix3f::Identity());
+Eigen::Vector3d Zero3d(0, 0, 0);
+Eigen::Vector3f Zero3f(0, 0, 0);
+// Eigen::Vector3d Lidar_offset_to_IMU(0.05512, 0.02226, 0.0297); // Horizon
+Eigen::Vector3d Lidar_offset_to_IMU(0.04165, 0.02326, -0.0284); // Avia
+
+struct MeasureGroup     // Lidar data and imu dates for the curent process
+{
+    MeasureGroup()
+    {
+        this->lidar.reset(new PointCloudXYZI());
+    };
+    double lidar_beg_time;
+    PointCloudXYZI::Ptr lidar;
+    std::deque<sensor_msgs::Imu::ConstPtr> imu;
+};
+
+struct StatesGroup
+{
+    StatesGroup() {
+		this->rot_end = Eigen::Matrix3d::Identity();
+		this->pos_end = Zero3d;
+        this->vel_end = Zero3d;
+        this->bias_g  = Zero3d;
+        this->bias_a  = Zero3d;
+        this->gravity = Zero3d;
+        this->cov     = Eigen::Matrix<double,DIM_OF_STATES,DIM_OF_STATES>::Identity() * INIT_COV;
+	};
+
+    StatesGroup& operator+=(const Eigen::Matrix<double, DIM_OF_STATES, 1> &state_add)
+	{
+		this->rot_end = this->rot_end * Exp(state_add(0,0), state_add(1,0), state_add(2,0));
+		this->pos_end += state_add.block<3,1>(3,0);
+        this->vel_end += state_add.block<3,1>(6,0);
+        this->bias_g  += state_add.block<3,1>(9,0);
+        this->bias_a  += state_add.block<3,1>(12,0);
+        this->gravity += state_add.block<3,1>(15,0);
+		return *this;
+	};
+
+	Eigen::Matrix3d rot_end;      // the estimated attitude (rotation matrix) at the end lidar point
+    Eigen::Vector3d pos_end;      // the estimated position at the end lidar point (world frame)
+    Eigen::Vector3d vel_end;      // the estimated velocity at the end lidar point (world frame)
+    Eigen::Vector3d bias_g;       // gyroscope bias
+    Eigen::Vector3d bias_a;       // accelerator bias
+    Eigen::Vector3d gravity;      // the estimated gravity acceleration
+    Eigen::Matrix<double, DIM_OF_STATES, DIM_OF_STATES>  cov;     // states covariance
+};
+
+template<typename T>
+T rad2deg(T radians)
+{
+  return radians * 180.0 / PI_M;
+}
+
+template<typename T>
+T deg2rad(T degrees)
+{
+  return degrees * PI_M / 180.0;
+}
+
+template<typename T>
+auto set_pose6d(const double t, const Eigen::Matrix<T, 3, 1> &a, const Eigen::Matrix<T, 3, 1> &g, \
+                const Eigen::Matrix<T, 3, 1> &v, const Eigen::Matrix<T, 3, 1> &p, const Eigen::Matrix<T, 3, 3> &R)
+{
+    Pose6D rot_kp;
+    rot_kp.offset_time = t;
+    for (int i = 0; i < 3; i++)
+    {
+        rot_kp.acc[i] = a(i);
+        rot_kp.gyr[i] = g(i);
+        rot_kp.vel[i] = v(i);
+        rot_kp.pos[i] = p(i);
+        for (int j = 0; j < 3; j++)  rot_kp.rot[i*3+j] = R(i,j);
+    }
+    // Eigen::Map<Eigen::Matrix3d>(rot_kp.rot, 3,3) = R;
+    return std::move(rot_kp);
+}
+
+template<typename T>
+Eigen::Matrix<T, 3, 1> RotMtoEuler(const Eigen::Matrix<T, 3, 3> &rot)
+{
+    T sy = sqrt(rot(0,0)*rot(0,0) + rot(1,0)*rot(1,0));
+    bool singular = sy < 1e-6;
+    T x, y, z;
+    if(!singular)
+    {
+        x = atan2(rot(2, 1), rot(2, 2));
+        y = atan2(-rot(2, 0), sy);   
+        z = atan2(rot(1, 0), rot(0, 0));  
+    }
+    else
+    {    
+        x = atan2(-rot(1, 2), rot(1, 1));    
+        y = atan2(-rot(2, 0), sy);    
+        z = 0;
+    }
+    Eigen::Matrix<T, 3, 1> ang(x, y, z);
+    return ang;
+}
+
+#endif