diff --git a/CMakeLists.txt b/CMakeLists.txt
new file mode 100755
index 0000000..979512c
--- /dev/null
+++ b/CMakeLists.txt
@@ -0,0 +1,87 @@
+cmake_minimum_required(VERSION 2.8.3)
+project(point_lio_unilidar)
+
+SET(CMAKE_BUILD_TYPE "Debug")
+
+ADD_COMPILE_OPTIONS(-std=c++14 )
+set(CMAKE_CXX_FLAGS "-std=c++14 -O3" )
+set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -fexceptions" )
+
+set(CMAKE_CXX_STANDARD 14)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+set(CMAKE_CXX_EXTENSIONS OFF)
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++14 -pthread -std=c++0x -std=c++14 -fexceptions")
+
+add_definitions(-DROOT_DIR=\"${CMAKE_CURRENT_SOURCE_DIR}/\")
+message("DROOT_DIR = ${DROOT_DIR}")
+
+message("Current CPU archtecture: ${CMAKE_SYSTEM_PROCESSOR}")
+if(CMAKE_SYSTEM_PROCESSOR MATCHES "(x86)|(X86)|(amd64)|(AMD64)" )
+ include(ProcessorCount)
+ ProcessorCount(N)
+ message("Processer number: ${N}")
+ if(N GREATER 5)
+ add_definitions(-DMP_EN)
+ add_definitions(-DMP_PROC_NUM=4)
+ message("core for MP: 3")
+ elseif(N GREATER 3)
+ math(EXPR PROC_NUM "${N} - 2")
+ add_definitions(-DMP_EN)
+ add_definitions(-DMP_PROC_NUM="${PROC_NUM}")
+ message("core for MP: ${PROC_NUM}")
+ else()
+ add_definitions(-DMP_PROC_NUM=1)
+ endif()
+else()
+ add_definitions(-DMP_PROC_NUM=1)
+endif()
+
+find_package(OpenMP QUIET)
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
+set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")
+
+find_package(PythonLibs REQUIRED)
+find_path(MATPLOTLIB_CPP_INCLUDE_DIRS "matplotlibcpp.h")
+
+find_package(catkin REQUIRED COMPONENTS
+ geometry_msgs
+ nav_msgs
+ sensor_msgs
+ roscpp
+ rospy
+ std_msgs
+ pcl_ros
+ tf
+ # livox_ros_driver
+ message_generation
+ eigen_conversions
+)
+
+find_package(Eigen3 REQUIRED)
+find_package(PCL 1.8 REQUIRED)
+
+message(Eigen: ${EIGEN3_INCLUDE_DIR})
+
+include_directories(
+ ${catkin_INCLUDE_DIRS}
+ ${EIGEN3_INCLUDE_DIR}
+ ${PCL_INCLUDE_DIRS}
+ ${PYTHON_INCLUDE_DIRS}
+ include
+)
+
+catkin_package(
+ CATKIN_DEPENDS geometry_msgs nav_msgs roscpp rospy std_msgs message_runtime
+ DEPENDS EIGEN3 PCL
+ INCLUDE_DIRS
+)
+
+add_executable(pointlio_mapping
+ src/laserMapping.cpp
+ include/ikd-Tree/ikd_Tree.cpp
+ src/parameters.cpp
+ src/preprocess.cpp
+ src/Estimator.cpp
+)
+target_link_libraries(pointlio_mapping ${catkin_LIBRARIES} ${PCL_LIBRARIES} ${PYTHON_LIBRARIES})
+target_include_directories(pointlio_mapping PRIVATE ${PYTHON_INCLUDE_DIRS})
\ No newline at end of file
diff --git a/Log/guide.md b/Log/guide.md
new file mode 100755
index 0000000..883b3d4
--- /dev/null
+++ b/Log/guide.md
@@ -0,0 +1 @@
+Here saved the debug records which can be drew by the ../log/plot.py. The record function can be found frm the MACRO: DEBUG_FILE_DIR(name) in common_lib.h.
diff --git a/Log/imu.txt b/Log/imu.txt
new file mode 100755
index 0000000..e69de29
diff --git a/Log/imu_pbp.txt b/Log/imu_pbp.txt
new file mode 100644
index 0000000..e69de29
diff --git a/Log/mat_out.txt b/Log/mat_out.txt
new file mode 100644
index 0000000..e69de29
diff --git a/Log/plot.py b/Log/plot.py
new file mode 100755
index 0000000..b5d1c1a
--- /dev/null
+++ b/Log/plot.py
@@ -0,0 +1,46 @@
+# import matplotlib
+# matplotlib.use('Agg')
+import numpy as np
+import matplotlib.pyplot as plt
+
+a_out=np.loadtxt('mat_out.txt')
+#######for normal#######
+fig, axs = plt.subplots(3,2)
+lab_out = ['', 'out-x', 'out-y', 'out-z']
+plot_ind = range(7,10)
+time=a_out[:,0]
+axs[0,0].set_title('Attitude')
+axs[1,0].set_title('Translation')
+axs[2,0].set_title('Velocity')
+axs[0,1].set_title('bg')
+axs[1,1].set_title('ba')
+axs[2,1].set_title('Gravity')
+for i in range(1,4):
+ for j in range(6):
+ axs[j%3, j//3].plot(time, a_out[:,i+j*3],'.-', label=lab_out[i])
+ for j in range(6):
+ axs[j%3, j//3].grid()
+ axs[j%3, j//3].legend()
+plt.grid()
+ #######for normal#######
+
+
+#### Draw IMU data
+#fig, axs = plt.subplots(2)
+#imu=np.loadtxt('imu_pbp.txt')
+#time=imu[:,0]
+#axs[0].set_title('Gyroscope')
+#axs[1].set_title('Accelerameter')
+#lab_1 = ['gyr-x', 'gyr-y', 'gyr-z']
+#lab_2 = ['acc-x', 'acc-y', 'acc-z']
+#for i in range(3):
+ #if i==1:
+# axs[0].plot(time, imu[:,i+1],'.-', label=lab_1[i])
+# axs[1].plot(time, imu[:,i+4],'.-', label=lab_2[i])
+#for i in range(2):
+ #axs[i].set_xlim(386,389)
+# axs[i].grid()
+# axs[i].legend()
+#plt.grid()
+
+plt.show()
diff --git a/Log/plot_imu.py b/Log/plot_imu.py
new file mode 100755
index 0000000..0177273
--- /dev/null
+++ b/Log/plot_imu.py
@@ -0,0 +1,125 @@
+# import matplotlib
+# matplotlib.use('Agg')
+import numpy as np
+import matplotlib.pyplot as plt
+
+#### Draw IMU data
+fig, axs = plt.subplots(2)
+imu=np.loadtxt('imu_pbp.txt')
+time=imu[:,0]
+axs[0].set_title('Gyroscope')
+axs[1].set_title('Accelerameter')
+lab_1 = ['gyr-x', 'gyr-y', 'gyr-z']
+lab_2 = ['acc-x', 'acc-y', 'acc-z']
+for i in range(3):
+ #if i==1:
+ axs[0].plot(time, imu[:,i+1],'.-', label=lab_1[i])
+ axs[1].plot(time, imu[:,i+4],'.-', label=lab_2[i])
+for i in range(2):
+ #axs[i].set_xlim(386,389)
+ axs[i].grid()
+ axs[i].legend()
+plt.grid()
+
+#fig, axs = plt.subplots(5)
+#axs[0].set_title('miss')
+#axs[1].set_title('miss')
+#axs[2].set_title('miss')
+#axs[3].set_title('miss')
+#axs[4].set_title('miss')
+#len_time1 = np.arange(0,1977)
+#len_time2 = np.arange(1977, 3954)
+#len_time3 = np.arange(3954,5931)
+#len_time4 = np.arange(5931,7908)
+#len_time5 = np.arange(7908,9885)
+ #if i==1:
+#axs[0].plot(len_time1, time[0:1977],'.-', label='check')
+#axs[1].plot(len_time2, time[1977:3954],'.-', label='check')
+#axs[2].plot(len_time3, time[3954:5931],'.-', label='check')
+#axs[3].plot(len_time4, time[5931:7908],'.-', label='check')
+#axs[4].plot(len_time5, time[7908:9885],'.-', label='check')
+ #axs[i].set_xlim(386,389)
+#axs[0].grid()
+#axs[0].legend()
+#axs[1].grid()
+#axs[1].legend()
+#axs[2].grid()
+#axs[2].legend()
+#axs[3].grid()
+#axs[3].legend()
+#axs[4].grid()
+#axs[4].legend()
+#plt.grid()
+
+#fig, axs = plt.subplots(5)
+#axs[0].set_title('miss')
+#axs[1].set_title('miss')
+#axs[2].set_title('miss')
+#axs[3].set_title('miss')
+#axs[4].set_title('miss')
+#len_time1 = np.arange(9885,9885+1977)
+#len_time2 = np.arange(9885+1977,9885+3954)
+#len_time3 = np.arange(9885+3954,9885+5931)
+#len_time4 = np.arange(9885+5931,9885+7908)
+#len_time5 = np.arange(9885+7908,9885+9885)
+ #if i==1:
+#axs[0].plot(len_time1, time[9885+0:9885+1977],'.-', label='check')
+#axs[1].plot(len_time2, time[9885+1977:9885+3954],'.-', label='check')
+#axs[2].plot(len_time3, time[9885+3954:9885+5931],'.-', label='check')
+#axs[3].plot(len_time4, time[9885+5931:9885+7908],'.-', label='check')
+#axs[4].plot(len_time5, time[9885+7908:9885+9885],'.-', label='check')
+ #axs[i].set_xlim(386,389)
+#axs[0].grid()
+#axs[0].legend()
+#axs[1].grid()
+#axs[1].legend()
+#axs[2].grid()
+#axs[2].legend()
+#axs[3].grid()
+#axs[3].legend()
+#axs[4].grid()
+#axs[4].legend()
+#plt.grid()
+
+# #### Draw time calculation
+# plt.figure(3)
+# fig = plt.figure()
+# font1 = {'family' : 'Times New Roman',
+# 'weight' : 'normal',
+# 'size' : 12,
+# }
+# c="red"
+# a_out1=np.loadtxt('Log/mat_out_time_indoor1.txt')
+# a_out2=np.loadtxt('Log/mat_out_time_indoor2.txt')
+# a_out3=np.loadtxt('Log/mat_out_time_outdoor.txt')
+# # n = a_out[:,1].size
+# # time_mean = a_out[:,1].mean()
+# # time_se = a_out[:,1].std() / np.sqrt(n)
+# # time_err = a_out[:,1] - time_mean
+# # feat_mean = a_out[:,2].mean()
+# # feat_err = a_out[:,2] - feat_mean
+# # feat_se = a_out[:,2].std() / np.sqrt(n)
+# ax1 = fig.add_subplot(111)
+# ax1.set_ylabel('Effective Feature Numbers',font1)
+# ax1.boxplot(a_out1[:,2], showfliers=False, positions=[0.9])
+# ax1.boxplot(a_out2[:,2], showfliers=False, positions=[1.9])
+# ax1.boxplot(a_out3[:,2], showfliers=False, positions=[2.9])
+# ax1.set_ylim([0, 3000])
+
+# ax2 = ax1.twinx()
+# ax2.spines['right'].set_color('red')
+# ax2.set_ylabel('Compute Time (ms)',font1)
+# ax2.yaxis.label.set_color('red')
+# ax2.tick_params(axis='y', colors='red')
+# ax2.boxplot(a_out1[:,1]*1000, showfliers=False, positions=[1.1],boxprops=dict(color=c),capprops=dict(color=c),whiskerprops=dict(color=c))
+# ax2.boxplot(a_out2[:,1]*1000, showfliers=False, positions=[2.1],boxprops=dict(color=c),capprops=dict(color=c),whiskerprops=dict(color=c))
+# ax2.boxplot(a_out3[:,1]*1000, showfliers=False, positions=[3.1],boxprops=dict(color=c),capprops=dict(color=c),whiskerprops=dict(color=c))
+# ax2.set_xlim([0.5, 3.5])
+# ax2.set_ylim([0, 100])
+
+# plt.xticks([1,2,3], ('Outdoor Scene', 'Indoor Scene 1', 'Indoor Scene 2'))
+# # # print(time_se)
+# # # print(a_out3[:,2])
+# plt.grid()
+# plt.savefig("time.pdf", dpi=1200)
+plt.show()
diff --git a/Log/plot_out.py b/Log/plot_out.py
new file mode 100755
index 0000000..8dd80cb
--- /dev/null
+++ b/Log/plot_out.py
@@ -0,0 +1,178 @@
+# import matplotlib
+# matplotlib.use('Agg')
+import numpy as np
+import matplotlib.pyplot as plt
+
+# a_pre=np.loadtxt('mat_pre.txt')
+a_out=np.loadtxt('mat_out.txt')
+if((a_out.shape[1] != 19) & (a_out.shape[1] != 20)):
+ ######for ikfom
+ fig, axs = plt.subplots(4,2)
+ #lab_pre = ['', 'pre-x', 'pre-y', 'pre-z']
+ lab_out = ['', 'out-x', 'out-y', 'out-z']
+ plot_ind = range(7,10)
+ time=a_out[:,0]
+ axs[0,0].set_title('Attitude')
+ axs[1,0].set_title('Translation')
+ axs[2,0].set_title('Velocity')
+ axs[3,0].set_title('Angular velocity')
+ axs[0,1].set_title('Acceleration')
+ axs[1,1].set_title('Gravity')
+ axs[2,1].set_title('bg')
+ axs[3,1].set_title('ba')
+ for i in range(1,4):
+ for j in range(8):
+ #axs[j%4, j//4].plot(time, a_pre[:,i+j*3],'.-', label=lab_pre[i])
+ axs[j%4, j//4].plot(time, a_out[:,i+j*3],'.-', label=lab_out[i])
+ for j in range(8):
+ # axs[j].set_xlim(386,389)
+ axs[j%4, j//4].grid()
+ axs[j%4, j//4].legend()
+ plt.grid()
+ ######for ikfom#######
+else:
+ #######for normal#######
+ fig, axs = plt.subplots(3,2)
+ lab_pre = ['', 'pre-x', 'pre-y', 'pre-z']
+ lab_out = ['', 'out-x', 'out-y', 'out-z']
+ plot_ind = range(7,10)
+ time=a_out[:,0]
+ time1 = a_pre[:,0]
+ axs[0,0].set_title('Attitude')
+ axs[1,0].set_title('Translation')
+ axs[2,0].set_title('Velocity')
+ axs[0,1].set_title('bg')
+ axs[1,1].set_title('ba')
+ axs[2,1].set_title('Gravity')
+ for i in range(1,4):
+ for j in range(6):
+ axs[j%3, j/3].plot(time1, a_pre[:,i+j*3],'.-', label=lab_pre[i])
+ axs[j%3, j/3].plot(time, a_out[:,i+j*3],'.-', label=lab_out[i])
+ for j in range(6):
+ # axs[j].set_xlim(386,389)
+ axs[j%3, j//3].grid()
+ axs[j%3, j//3].legend()
+ plt.grid()
+ #######for normal#######
+
+
+#### Draw IMU data
+fig, axs = plt.subplots(2)
+imu=np.loadtxt('imu_pbp.txt')
+time=imu[:,0]
+axs[0].set_title('Gyroscope')
+axs[1].set_title('Accelerameter')
+lab_1 = ['gyr-x', 'gyr-y', 'gyr-z']
+lab_2 = ['acc-x', 'acc-y', 'acc-z']
+for i in range(3):
+ #if i==1:
+ axs[0].plot(time, imu[:,i+1],'.-', label=lab_1[i])
+ axs[1].plot(time, imu[:,i+4],'.-', label=lab_2[i])
+for i in range(2):
+ #axs[i].set_xlim(386,389)
+ axs[i].grid()
+ axs[i].legend()
+plt.grid()
+
+#fig, axs = plt.subplots(5)
+#axs[0].set_title('miss')
+#axs[1].set_title('miss')
+#axs[2].set_title('miss')
+#axs[3].set_title('miss')
+#axs[4].set_title('miss')
+#len_time1 = np.arange(0,1977)
+#len_time2 = np.arange(1977, 3954)
+#len_time3 = np.arange(3954,5931)
+#len_time4 = np.arange(5931,7908)
+#len_time5 = np.arange(7908,9885)
+ #if i==1:
+#axs[0].plot(len_time1, time[0:1977],'.-', label='check')
+#axs[1].plot(len_time2, time[1977:3954],'.-', label='check')
+#axs[2].plot(len_time3, time[3954:5931],'.-', label='check')
+#axs[3].plot(len_time4, time[5931:7908],'.-', label='check')
+#axs[4].plot(len_time5, time[7908:9885],'.-', label='check')
+ #axs[i].set_xlim(386,389)
+#axs[0].grid()
+#axs[0].legend()
+#axs[1].grid()
+#axs[1].legend()
+#axs[2].grid()
+#axs[2].legend()
+#axs[3].grid()
+#axs[3].legend()
+#axs[4].grid()
+#axs[4].legend()
+#plt.grid()
+
+#fig, axs = plt.subplots(5)
+#axs[0].set_title('miss')
+#axs[1].set_title('miss')
+#axs[2].set_title('miss')
+#axs[3].set_title('miss')
+#axs[4].set_title('miss')
+#len_time1 = np.arange(9885,9885+1977)
+#len_time2 = np.arange(9885+1977,9885+3954)
+#len_time3 = np.arange(9885+3954,9885+5931)
+#len_time4 = np.arange(9885+5931,9885+7908)
+#len_time5 = np.arange(9885+7908,9885+9885)
+ #if i==1:
+#axs[0].plot(len_time1, time[9885+0:9885+1977],'.-', label='check')
+#axs[1].plot(len_time2, time[9885+1977:9885+3954],'.-', label='check')
+#axs[2].plot(len_time3, time[9885+3954:9885+5931],'.-', label='check')
+#axs[3].plot(len_time4, time[9885+5931:9885+7908],'.-', label='check')
+#axs[4].plot(len_time5, time[9885+7908:9885+9885],'.-', label='check')
+ #axs[i].set_xlim(386,389)
+#axs[0].grid()
+#axs[0].legend()
+#axs[1].grid()
+#axs[1].legend()
+#axs[2].grid()
+#axs[2].legend()
+#axs[3].grid()
+#axs[3].legend()
+#axs[4].grid()
+#axs[4].legend()
+#plt.grid()
+
+# #### Draw time calculation
+# plt.figure(3)
+# fig = plt.figure()
+# font1 = {'family' : 'Times New Roman',
+# 'weight' : 'normal',
+# 'size' : 12,
+# }
+# c="red"
+# a_out1=np.loadtxt('Log/mat_out_time_indoor1.txt')
+# a_out2=np.loadtxt('Log/mat_out_time_indoor2.txt')
+# a_out3=np.loadtxt('Log/mat_out_time_outdoor.txt')
+# # n = a_out[:,1].size
+# # time_mean = a_out[:,1].mean()
+# # time_se = a_out[:,1].std() / np.sqrt(n)
+# # time_err = a_out[:,1] - time_mean
+# # feat_mean = a_out[:,2].mean()
+# # feat_err = a_out[:,2] - feat_mean
+# # feat_se = a_out[:,2].std() / np.sqrt(n)
+# ax1 = fig.add_subplot(111)
+# ax1.set_ylabel('Effective Feature Numbers',font1)
+# ax1.boxplot(a_out1[:,2], showfliers=False, positions=[0.9])
+# ax1.boxplot(a_out2[:,2], showfliers=False, positions=[1.9])
+# ax1.boxplot(a_out3[:,2], showfliers=False, positions=[2.9])
+# ax1.set_ylim([0, 3000])
+
+# ax2 = ax1.twinx()
+# ax2.spines['right'].set_color('red')
+# ax2.set_ylabel('Compute Time (ms)',font1)
+# ax2.yaxis.label.set_color('red')
+# ax2.tick_params(axis='y', colors='red')
+# ax2.boxplot(a_out1[:,1]*1000, showfliers=False, positions=[1.1],boxprops=dict(color=c),capprops=dict(color=c),whiskerprops=dict(color=c))
+# ax2.boxplot(a_out2[:,1]*1000, showfliers=False, positions=[2.1],boxprops=dict(color=c),capprops=dict(color=c),whiskerprops=dict(color=c))
+# ax2.boxplot(a_out3[:,1]*1000, showfliers=False, positions=[3.1],boxprops=dict(color=c),capprops=dict(color=c),whiskerprops=dict(color=c))
+# ax2.set_xlim([0.5, 3.5])
+# ax2.set_ylim([0, 100])
+
+# plt.xticks([1,2,3], ('Outdoor Scene', 'Indoor Scene 1', 'Indoor Scene 2'))
+# # # print(time_se)
+# # # print(a_out3[:,2])
+# plt.grid()
+# plt.savefig("time.pdf", dpi=1200)
+plt.show()
diff --git a/Log/pos_log.txt b/Log/pos_log.txt
new file mode 100644
index 0000000..e69de29
diff --git a/README.md b/README.md
index a9271b8..cbe1e09 100644
--- a/README.md
+++ b/README.md
@@ -1,2 +1,145 @@
# point_lio_unilidar
-An open-source robust and high-bandwidth lidar inertial odometry for Unitree L1 LiDAR.
+
+## Introduction
+
+This repository adapts the state-of-the-art lidar inertial odometry algorithm, `Point-LIO`, for use with our lidar, `Unitree Lidar L1`.
+
+`Unitree Lidar L1` is a low-cost large-FOV 3D lidar, which possesses these features:
+- large field of view (360° × 90°)
+- non-repetitive scanning
+- low cost (only $349)
+- suitable for applications in low-speed mobile robots
+
+If you want to learn more about our lidar, you can refer to the official website for details.
+-
+
+
+`Point-LIO` is a robust and high-bandwidth lidar inertial odometry (LIO) with the capability to provide accurate, high-frequency odometry and reliable mapping under severe vibrations and aggressive motions. If you need further information about the `Point-LIO` algorithm, you can refer to their official website and paper:
+-
+- [Point‐LIO: Robust High‐Bandwidth Light Detection and Ranging Inertial Odometry](https://onlinelibrary.wiley.com/doi/epdf/10.1002/aisy.202200459)
+
+![demo](./doc/demo.png)
+
+## Prerequisites
+
+### Ubuntu and ROS
+We tested our code on Ubuntu20.04 with [ROS noetic](http://wiki.ros.org/noetic/Installation/Ubuntu). Ubuntu18.04 and lower versions have problems of environments to support the Point-LIO, try to avoid using Point-LIO in those systems.
+
+You can refer to the official website to install ROS noetic:
+-
+
+Additional ROS package is required:
+```
+sudo apt-get install ros-xxx-pcl-conversions
+```
+
+### Eigen
+Following the official [Eigen installation](eigen.tuxfamily.org/index.php?title=Main_Page), or directly install Eigen by:
+```
+sudo apt-get install libeigen3-dev
+```
+
+### unitree_lidar_ros
+
+You should download and build [unitree_lidar_ros](https://github.com/unitreerobotics/unilidar_sdk/tree/main/unitree_lidar_ros/src/unitree_lidar_ros) follwing these steps:
+
+```
+git clone https://github.com/unitreerobotics/unilidar_sdk.git
+
+cd unilidar_sdk/unitree_lidar_ros
+
+catkin_make
+```
+
+
+## Build
+
+Clone this repository and run `catkin_make`:
+
+```
+mkdir -p catkin_point_lio_unilidar/src
+
+cd catkin_point_lio_unilidar/src
+
+git clone https://github.com/unitreerobotics/point_lio_unilidar.git
+
+cd ../..
+
+catkin_make
+```
+
+
+## Run
+
+### Run with Unilidar
+
+Firstly, you should connect our lidar to your PC serial port, and supply power for the lidar with a 12V charger.
+
+Besides, to ensure proper initialization of the IMU, it is advisable to keep the lidar in a stationary state during the initial few seconds of algorithm execution.
+
+Run `unilidar`:
+```
+cd unilidar_sdk/unitree_lidar_ros
+
+source devel/setup.bash
+
+roslaunch unitree_lidar_ros run.launch
+```
+
+Run `Point-LIO`:
+```
+cd catkin_unilidar_point_lio
+
+source devel/setup.bash
+
+roslaunch point_lio_unilidar mapping_unilidar.launch
+```
+
+
+After completion of the run, all cached pointcloud map will be saved to the following path:
+```
+catkin_point_lio_unilidar/src/point_lio_unilidar/PCD/scans.pcd
+```
+
+You can use the `pcl_viewer` tool to view this pcd file:
+```
+pcl_viewer scans.pcd
+```
+
+### Run with dataset
+
+If you don't have our lidar for now, you can download our dataset recorded with our lidar and run testify this algorithm with it.
+The download address is here:
+- [unilidar-2023-09-22-12-42-04.bag - 百度云](https://pan.baidu.com/s/1PD0e9R5Q9gxHM966nurkDQ) --> 提取码: ge3c
+
+
+Run `Point-LIO`:
+```
+cd catkin_point_lio_unilidar
+
+source devel/setup.bash
+
+roslaunch point_lio_unilidar mapping_unilidar.launch
+```
+
+Play the dataset you downloaded:
+```
+rosbag play unilidar-2023-09-22-12-42-04.bag
+```
+
+
+After completion of the run, all cached pointcloud map will be saved to the following path:
+```
+catkin_point_lio_unilidar/src/point_lio_unilidarPCD/scans.pcd
+```
+
+You can use the `pcl_viewer` tool to view this pcd file:
+```
+pcl_viewer scans.pcd
+```
+
+## Version History
+
+### v1.0.0 (2023.09.22)
+- Adapt `Point-LIO` for `Unitree Lidar L1`
+- Upload codes and dataset
\ No newline at end of file
diff --git a/config/avia.yaml b/config/avia.yaml
new file mode 100755
index 0000000..68d125e
--- /dev/null
+++ b/config/avia.yaml
@@ -0,0 +1,59 @@
+common:
+ lid_topic: "/livox/lidar"
+ imu_topic: "/livox/imu"
+ con_frame: false # true: if you need to combine several LiDAR frames into one
+ con_frame_num: 1 # the number of frames combined
+ cut_frame: false # true: if you need to cut one LiDAR frame into several subframes
+ cut_frame_time_interval: 0.1 # should be integral fraction of 1 / LiDAR frequency
+ time_lag_imu_to_lidar: 0.0 # Time offset between LiDAR and IMU calibrated by other algorithms, e.g., LI-Init (find in Readme)
+ # the timesample of IMU is transferred from the current timeline to LiDAR's timeline by subtracting this value
+
+preprocess:
+ lidar_type: 1
+ scan_line: 6
+ timestamp_unit: 1 # the unit of time/t field in the PointCloud2 rostopic: 0-second, 1-milisecond, 2-microsecond, 3-nanosecond.
+ blind: 1.0 # 剔除半径位于此范围内的点
+
+mapping:
+ imu_en: true
+ start_in_aggressive_motion: false # if true, a preknown gravity should be provided in following gravity_init
+ extrinsic_est_en: false # for aggressive motion, set this variable false
+ imu_time_inte: 0.005 # = 1 / frequency of IMU
+ satu_acc: 3.0 # the saturation value of IMU's acceleration. not related to the units
+ satu_gyro: 17.5 # the saturation value of IMU's angular velocity. not related to the units
+ # default: 35 rad/s
+ # for running `PULSAR.bag`: set it to 17.5 rad/s
+ acc_norm: 1.0 # 1.0 for g as unit, 9.81 for m/s^2 as unit of the IMU's acceleration
+ lidar_meas_cov: 0.001 # 0.001; 0.01
+ acc_cov_output: 500
+ gyr_cov_output: 1000
+ b_acc_cov: 0.0001
+ b_gyr_cov: 0.0001
+ imu_meas_acc_cov: 0.1 #0.1 # 0.1
+ imu_meas_omg_cov: 0.1 #0.01 # 0.1
+ gyr_cov_input: 0.01 # for IMU as input model
+ acc_cov_input: 0.1 # for IMU as input model
+ plane_thr: 0.1 # 0.05, the threshold for plane criteria, the smaller, the flatter a plane
+ match_s: 81
+ fov_degree: 90
+ det_range: 450.0
+ gravity_align: true # true to align the z axis of world frame with the direction of gravity, and the gravity direction should be specified below
+ gravity: [0.0, 0.0, -9.810] # [0.0, 9.810, 0.0] # gravity to be aligned
+ gravity_init: [0.0, 0.0, -9.810] # [0.0, 9.810, 0.0] # # preknown gravity in the first IMU body frame, use when imu_en is false or start from a non-stationary state
+ extrinsic_T: [ 0.04165, 0.02326, -0.0284 ]
+ extrinsic_R: [ 1, 0, 0,
+ 0, 1, 0,
+ 0, 0, 1 ]
+
+odometry:
+ publish_odometry_without_downsample: false # 发布没有降采样的里程计
+
+publish:
+ path_en: true # false: close the path output
+ scan_publish_en: true # false: close all the point cloud output
+ scan_bodyframe_pub_en: false # true: output the point cloud scans in IMU-body-frame
+
+pcd_save:
+ pcd_save_en: false
+ interval: -1 # how many LiDAR frames saved in each pcd file;
+ # -1 : all frames will be saved in ONE pcd file, may lead to memory crash when having too much frames.
\ No newline at end of file
diff --git a/config/horizon.yaml b/config/horizon.yaml
new file mode 100755
index 0000000..070bf90
--- /dev/null
+++ b/config/horizon.yaml
@@ -0,0 +1,57 @@
+common:
+ lid_topic: "/livox/lidar"
+ imu_topic: "/livox/imu"
+ con_frame: false # true: if you need to combine several LiDAR frames into one
+ con_frame_num: 1 # the number of frames combined
+ cut_frame: false # true: if you need to cut one LiDAR frame into several subframes
+ cut_frame_time_interval: 0.1 # should be integral fraction of 1 / LiDAR frequency
+ time_lag_imu_to_lidar: 0.0 # Time offset between LiDAR and IMU calibrated by other algorithms, e.g., LI-Init (find in Readme),
+ # the timesample of IMU is transferred from the current timeline to LiDAR's timeline by subtracting this value
+
+preprocess:
+ lidar_type: 1
+ scan_line: 6
+ timestamp_unit: 1 # the unit of time/t field in the PointCloud2 rostopic: 0-second, 1-milisecond, 2-microsecond, 3-nanosecond.
+ blind: 4.0
+
+mapping:
+ imu_en: true
+ start_in_aggressive_motion: false # if true, a preknown gravity should be provided in following gravity_init
+ extrinsic_est_en: false # for aggressive motion, set this variable false
+ imu_time_inte: 0.005 # = 1 / frequency of IMU
+ satu_acc: 3.0 # the saturation value of IMU's acceleration. not related to the units
+ satu_gyro: 35 # the saturation value of IMU's angular velocity. not related to the units
+ acc_norm: 1.0 # 1.0 for g as unit, 9.81 for m/s^2 as unit of the IMU's acceleration
+ lidar_meas_cov: 0.01 # 0.001
+ acc_cov_output: 500
+ gyr_cov_output: 1000
+ b_acc_cov: 0.0001
+ b_gyr_cov: 0.0001
+ imu_meas_acc_cov: 0.01 #0.1 # 2
+ imu_meas_omg_cov: 0.01 #0.1 # 2
+ gyr_cov_input: 0.01 # for IMU as input model
+ acc_cov_input: 0.1 # for IMU as input model
+ plane_thr: 0.1 # 0.05, the threshold for plane criteria, the smaller, the flatter a plane
+ match_s: 81
+ fov_degree: 100
+ det_range: 260.0
+ gravity_align: true # true to align the z axis of world frame with the direction of gravity, and the gravity direction should be specified below
+ gravity: [0.0, 0.0, -9.810] # [0.0, 9.810, 0.0] # gravity to be aligned
+ gravity_init: [0.0, 0.0, -9.810] # [0.0, 9.810, 0.0] # # preknown gravity in the first IMU body frame, use when imu_en is false or start from a non-stationary state
+ extrinsic_T: [ 0.05512, 0.02226, -0.0297 ]
+ extrinsic_R: [ 1, 0, 0,
+ 0, 1, 0,
+ 0, 0, 1 ]
+
+odometry:
+ publish_odometry_without_downsample: false
+
+publish:
+ path_en: true # false: close the path output
+ scan_publish_en: true # false: close all the point cloud output
+ scan_bodyframe_pub_en: false # true: output the point cloud scans in IMU-body-frame
+
+pcd_save:
+ pcd_save_en: false
+ interval: -1 # how many LiDAR frames saved in each pcd file;
+ # -1 : all frames will be saved in ONE pcd file, may lead to memory crash when having too much frames.
\ No newline at end of file
diff --git a/config/ouster64.yaml b/config/ouster64.yaml
new file mode 100755
index 0000000..2ba04b4
--- /dev/null
+++ b/config/ouster64.yaml
@@ -0,0 +1,57 @@
+common:
+ lid_topic: "/os_cloud_node/points"
+ imu_topic: "/os_cloud_node/imu"
+ con_frame: false # true: if you need to combine several LiDAR frames into one
+ con_frame_num: 1 # the number of frames combined
+ cut_frame: false # true: if you need to cut one LiDAR frame into several subframes
+ cut_frame_time_interval: 0.1 # should be integral fraction of 1 / LiDAR frequency
+ time_lag_imu_to_lidar: 0.0 # Time offset between LiDAR and IMU calibrated by other algorithms, e.g., LI-Init (find in Readme)
+ # the timesample of IMU is transferred from the current timeline to LiDAR's timeline by subtracting this value
+
+preprocess:
+ lidar_type: 3 # 2 #velodyne # 1 Livox Avia LiDAR
+ scan_line: 32 # 32 #velodyne 6 avia
+ timestamp_unit: 3 # the unit of time/t field in the PointCloud2 rostopic: 0-second, 1-milisecond, 2-microsecond, 3-nanosecond.
+ blind: 0.20
+
+mapping:
+ imu_en: true
+ start_in_aggressive_motion: false # if true, a preknown gravity should be provided in following gravity_init
+ extrinsic_est_en: false # for aggressive motion, set this variable false
+ imu_time_inte: 0.01 # = 1 / frequency of IMU
+ satu_acc: 30.0 # the saturation value of IMU's acceleration. not related to the units
+ satu_gyro: 35 # the saturation value of IMU's angular velocity. not related to the units
+ acc_norm: 9.81 # 1.0 for g as unit, 9.81 for m/s^2 as unit of the IMU's acceleration
+ lidar_meas_cov: 0.1 # 0.01
+ acc_cov_output: 500
+ gyr_cov_output: 1000
+ b_acc_cov: 0.0001
+ b_gyr_cov: 0.0001
+ imu_meas_acc_cov: 0.1 #0.1 # 2
+ imu_meas_omg_cov: 0.1 #0.1 # 2
+ gyr_cov_input: 0.01 # for IMU as input model
+ acc_cov_input: 0.1 # for IMU as input model
+ plane_thr: 0.1 # 0.05, the threshold for plane criteria, the smaller, the flatter a plane
+ match_s: 81
+ fov_degree: 180
+ det_range: 150.0
+ gravity_align: true # true to align the z axis of world frame with the direction of gravity, and the gravity direction should be specified below
+ gravity: [0.0, 0.0, -9.810] # [0.0, 9.810, 0.0] # gravity to be aligned
+ gravity_init: [0.0, 0.0, -9.810] # [0.0, 9.810, 0.0] # # preknown gravity in the first IMU body frame, use when imu_en is false or start from a non-stationary state
+ extrinsic_T: [0.0, 0.0, 0.0]
+ extrinsic_R: [ 1, 0, 0,
+ 0, 1, 0,
+ 0, 0, 1 ]
+
+odometry:
+ publish_odometry_without_downsample: false
+
+publish:
+ path_en: true # false: close the path output
+ scan_publish_en: true # false: close all the point cloud output
+ scan_bodyframe_pub_en: false # true: output the point cloud scans in IMU-body-frame
+
+pcd_save:
+ pcd_save_en: false
+ interval: -1 # how many LiDAR frames saved in each pcd file;
+ # -1 : all frames will be saved in ONE pcd file, may lead to memory crash when having too much frames.
\ No newline at end of file
diff --git a/config/unilidar.yaml b/config/unilidar.yaml
new file mode 100755
index 0000000..f3ff469
--- /dev/null
+++ b/config/unilidar.yaml
@@ -0,0 +1,59 @@
+common:
+ lid_topic: "/unilidar/cloud"
+ imu_topic: "/unilidar/imu"
+ con_frame: false # true: if you need to combine several LiDAR frames into one
+ con_frame_num: 1 # the number of frames combined
+ cut_frame: false # true: if you need to cut one LiDAR frame into several subframes
+ cut_frame_time_interval: 0.1 # should be integral fraction of 1 / LiDAR frequency
+ time_lag_imu_to_lidar: 0.0 # Time offset between LiDAR and IMU calibrated by other algorithms, e.g., LI-Init (find in Readme)
+ # the timesample of IMU is transferred from the current timeline to LiDAR's timeline by subtracting this value
+
+preprocess:
+ lidar_type: 5
+ scan_line: 18
+ timestamp_unit: 0 # the unit of time/t field in the PointCloud2 rostopic: 0-second, 1-milisecond, 2-microsecond, 3-nanosecond.
+ blind: 0.5
+
+mapping:
+ imu_en: true
+ start_in_aggressive_motion: false # if true, a preknown gravity should be provided in following gravity_init
+ extrinsic_est_en: false # for aggressive motion, set this variable false
+ imu_time_inte: 0.004 # = 1 / frequency of IMU
+ satu_acc: 30.0 # the saturation value of IMU's acceleration. not related to the units
+ satu_gyro: 35 # the saturation value of IMU's angular velocity. not related to the units
+ acc_norm: 9.81 # 1.0 for g as unit, 9.81 for m/s^2 as unit of the IMU's acceleration
+ lidar_meas_cov: 0.01 # 0.001
+ acc_cov_output: 500
+ gyr_cov_output: 1000
+ b_acc_cov: 0.0001
+ b_gyr_cov: 0.0001
+ imu_meas_acc_cov: 0.1 #0.1 # 2
+ imu_meas_omg_cov: 0.1 #0.1 # 2
+ gyr_cov_input: 0.01 # for IMU as input model
+ acc_cov_input: 0.1 # for IMU as input model
+ plane_thr: 0.1 # 0.05, the threshold for plane criteria, the smaller, the flatter a plane
+ match_s: 81
+ fov_degree: 180
+ det_range: 100.0
+ gravity_align: true # true to align the z axis of world frame with the direction of gravity, and the gravity direction should be specified below
+ gravity: [0.0, 0.0, -9.810] # [0.0, 9.810, 0.0] # gravity to be aligned
+ gravity_init: [0.0, 0.0, -9.810] # [0.0, 9.810, 0.0] # # preknown gravity in the first IMU body frame, use when imu_en is false or start from a non-stationary state
+
+ # transform from imu to lidar
+ extrinsic_T: [ 0.007698, 0.014655, -0.00667] # ulhk # [-0.5, 1.4, 1.5] # utbm
+ extrinsic_R: [ 1, 0, 0,
+ 0, 1, 0,
+ 0, 0, 1 ] # ulhk 4 utbm 3
+
+odometry:
+ publish_odometry_without_downsample: enable
+
+publish:
+ path_en: true # false: close the path output
+ scan_publish_en: true # false: close all the point cloud output
+ scan_bodyframe_pub_en: false # true: output the point cloud scans in IMU-body-frame
+
+pcd_save:
+ pcd_save_en: true # save map to pcd file
+ interval: -1 # how many LiDAR frames saved in each pcd file;
+ # -1 : all frames will be saved in ONE pcd file, may lead to memory crash when having too much frames.
\ No newline at end of file
diff --git a/config/velody16.yaml b/config/velody16.yaml
new file mode 100755
index 0000000..ed59f00
--- /dev/null
+++ b/config/velody16.yaml
@@ -0,0 +1,63 @@
+common:
+ lid_topic: "/velodyne_points"
+ imu_topic: "/imu/data"
+ con_frame: false # true: if you need to combine several LiDAR frames into one
+ con_frame_num: 1 # the number of frames combined
+ cut_frame: false # true: if you need to cut one LiDAR frame into several subframes
+ cut_frame_time_interval: 0.1 # should be integral fraction of 1 / LiDAR frequency
+ time_lag_imu_to_lidar: 0.0 # Time offset between LiDAR and IMU calibrated by other algorithms, e.g., LI-Init (find in Readme)
+ # the timesample of IMU is transferred from the current timeline to LiDAR's timeline by subtracting this value
+
+preprocess:
+ lidar_type: 2
+ scan_line: 32
+ timestamp_unit: 2 # the unit of time/t field in the PointCloud2 rostopic: 0-second, 1-milisecond, 2-microsecond, 3-nanosecond.
+ blind: 2.0
+
+mapping:
+ imu_en: true
+ start_in_aggressive_motion: false # if true, a preknown gravity should be provided in following gravity_init
+ extrinsic_est_en: false # for aggressive motion, set this variable false
+ imu_time_inte: 0.01 # = 1 / frequency of IMU
+ satu_acc: 30.0 # the saturation value of IMU's acceleration. not related to the units
+ satu_gyro: 35 # the saturation value of IMU's angular velocity. not related to the units
+ acc_norm: 9.81 # 1.0 for g as unit, 9.81 for m/s^2 as unit of the IMU's acceleration
+ lidar_meas_cov: 0.01 # 0.001
+ acc_cov_output: 500
+ gyr_cov_output: 1000
+ b_acc_cov: 0.0001
+ b_gyr_cov: 0.0001
+ imu_meas_acc_cov: 0.1 #0.1 # 2
+ imu_meas_omg_cov: 0.1 #0.1 # 2
+ gyr_cov_input: 0.01 # for IMU as input model
+ acc_cov_input: 0.1 # for IMU as input model
+ plane_thr: 0.1 # 0.05, the threshold for plane criteria, the smaller, the flatter a plane
+ match_s: 81
+ fov_degree: 180
+ det_range: 100.0
+ gravity_align: true # true to align the z axis of world frame with the direction of gravity, and the gravity direction should be specified below
+ gravity: [0.0, 0.0, -9.810] # [0.0, 9.810, 0.0] # gravity to be aligned
+ gravity_init: [0.0, 0.0, -9.810] # [0.0, 9.810, 0.0] # # preknown gravity in the first IMU body frame, use when imu_en is false or start from a non-stationary state
+ extrinsic_T: [ 0, 0, 0.28] # ulhk # [-0.5, 1.4, 1.5] # utbm
+ # extrinsic_R: [ 0, 1, 0,
+ # -1, 0, 0,
+ # 0, 0, 1 ] # ulhk 5 6
+ # extrinsic_R: [ 0, -1, 0,
+ # 1, 0, 0,
+ # 0, 0, 1 ] # utbm 1, 2
+ extrinsic_R: [ 1, 0, 0,
+ 0, 1, 0,
+ 0, 0, 1 ] # ulhk 4 utbm 3
+
+odometry:
+ publish_odometry_without_downsample: false
+
+publish:
+ path_en: true # false: close the path output
+ scan_publish_en: true # false: close all the point cloud output
+ scan_bodyframe_pub_en: false # true: output the point cloud scans in IMU-body-frame
+
+pcd_save:
+ pcd_save_en: false
+ interval: -1 # how many LiDAR frames saved in each pcd file;
+ # -1 : all frames will be saved in ONE pcd file, may lead to memory crash when having too much frames.
\ No newline at end of file
diff --git a/doc/demo.png b/doc/demo.png
new file mode 100644
index 0000000..4b65177
Binary files /dev/null and b/doc/demo.png differ
diff --git a/include/FOV_Checker/FOV_Checker.cpp b/include/FOV_Checker/FOV_Checker.cpp
new file mode 100755
index 0000000..69032bf
--- /dev/null
+++ b/include/FOV_Checker/FOV_Checker.cpp
@@ -0,0 +1,472 @@
+#include "FOV_Checker.h"
+
+FOV_Checker::FOV_Checker(){
+ // fp = fopen("/home/ecstasy/catkin_ws/fov_data.csv","w");
+ // fprintf(fp,"cur_pose_x,cur_pose_y,cur_pose_z,axis_x,axis_y,axis_z,theta,depth\n");
+ // fclose(fp);
+}
+
+FOV_Checker::~FOV_Checker(){
+
+}
+
+void FOV_Checker::Set_Env(BoxPointType env_param){
+ env = env_param;
+}
+
+void FOV_Checker::Set_BoxLength(double box_len_param){
+ box_length = box_len_param;
+}
+
+void round_v3d(Eigen::Vector3d &vec, int decimal){
+ double tmp;
+ int t;
+ for (int i = 0; i < 3; i++){
+ t = pow(10,decimal);
+ tmp = round(vec(i)*t);
+ vec(i) = tmp/t;
+ }
+ return;
+}
+
+void FOV_Checker::check_fov(Eigen::Vector3d cur_pose, Eigen::Vector3d axis, double theta, double depth, vector &boxes){
+ round_v3d(cur_pose,4);
+ round_v3d(axis,3);
+ axis = axis/axis.norm();
+ // fp = fopen("/home/ecstasy/catkin_ws/fov_data.csv","a");
+ // fprintf(fp,"%f,%f,%f,%f,%f,%f,%0.4f,%0.1f,",cur_pose(0),cur_pose(1),cur_pose(2),axis(0),axis(1),axis(2),theta,depth);
+ // fclose(fp);
+ // cout << "cur_pose: " << cur_pose.transpose() << endl;
+ // cout<< "axis: " << axis.transpose() << endl;
+ // cout<< "theta: " << theta << " depth: " << depth << endl;
+ // cout<< "env: " << env.vertex_min[0] << " " << env.vertex_max[0] << endl;
+ double axis_angle[6], min_angle, gap, plane_u_min, plane_u_max;
+ Eigen::Vector3d plane_w, plane_u, plane_v, center_point, start_point, box_p;
+ Eigen::Vector3d box_p_min, box_p_max;
+ int i, j, k, index, maxn, start_i, max_uN, max_vN, max_ulogN, u_min, u_max;
+ bool flag = false, box_found = false;
+ boxes.clear();
+ BoxPointType box;
+ axis_angle[0] = acos(axis(0));
+ axis_angle[1] = acos(axis(1));
+ axis_angle[2] = acos(axis(2));
+ axis_angle[3] = acos(-axis(0));
+ axis_angle[4] = acos(-axis(1));
+ axis_angle[5] = acos(-axis(2));
+ index = 1;
+ min_angle = axis_angle[0];
+ for (i=1;i<6;i++){
+ if (axis_angle[i]=0){
+ box_p = box_p_min + plane_u * box_length * (u_min + pow(2,k)) + plane_v * box_length + plane_w * box_length;
+ box.vertex_min[0] = box_p_min(0);
+ box.vertex_min[1] = box_p_min(1);
+ box.vertex_min[2] = box_p_min(2);
+ box.vertex_max[0] = box_p(0);
+ box.vertex_max[1] = box_p(1);
+ box.vertex_max[2] = box_p(2);
+ if (!check_box(cur_pose, axis, theta, depth, box)) u_min = u_min + pow(2,k);
+ k = k-1;
+ }
+ k = max_ulogN;
+ u_max = 0;
+ while (k>=0){
+ box_p = box_p_max - plane_u * box_length * (u_max + pow(2,k)) - plane_v * box_length - plane_w * box_length;
+ box.vertex_min[0] = box_p(0);
+ box.vertex_min[1] = box_p(1);
+ box.vertex_min[2] = box_p(2);
+ box.vertex_max[0] = box_p_max(0);
+ box.vertex_max[1] = box_p_max(1);
+ box.vertex_max[2] = box_p_max(2);
+ if (!check_box(cur_pose, axis, theta, depth, box)) u_max = u_max + pow(2,k);
+
+ k = k-1;
+ }
+ u_max = max(0, max_uN - u_max - 1);
+ box_found = false;
+ //printf("---- u_min -> u_max: %d->%d\n",u_min,u_max);
+ for (k = u_min; k <= u_max; k++){
+ box_p = box_p_min + plane_u * box_length * k;
+ box.vertex_min[0] = box_p(0);
+ box.vertex_min[1] = box_p(1);
+ box.vertex_min[2] = box_p(2);
+ box.vertex_max[0] = box_p(0) + box_length;
+ box.vertex_max[1] = box_p(1) + box_length;
+ box.vertex_max[2] = box_p(2) + box_length;
+ if (check_box_in_env(box)){
+ //printf("---- FOUND: (%0.3f,%0.3f,%0.3f),(%0.3f,%0.3f,%0.3f)\n",box.vertex_min[0],box.vertex_min[1],box.vertex_min[2],box.vertex_max[0],box.vertex_max[1],box.vertex_max[2]);
+ box_found = true;
+ boxes.push_back(box);
+ }
+ }
+ if (box_found) {
+ flag = true;
+ } else {
+ if (j>1) break;
+ }
+ }
+ for (j = 1; j <= max_vN; j++){
+ k = max_ulogN;
+ u_min = 0;
+ box_p_min = start_point.cwiseProduct(plane_w + plane_v) + plane_u * plane_u_min - plane_v * box_length * j;
+ box_p_max = plane_u * plane_u_max + start_point.cwiseProduct(plane_w + plane_v) - plane_v * box_length * (j-1) + plane_w * box_length;
+ //printf("---- DOWNSIDE (%0.3f,%0.3f,%0.3f),(%0.3f,%0.3f,%0.3f)\n",box_p_min[0],box_p_min[1],box_p_min[2],box_p_max[0],box_p_max[1],box_p_max[2]);
+
+ while (k>=0){
+ box_p = box_p_min + plane_u * box_length * (u_min + pow(2,k)) + plane_v * box_length + plane_w * box_length;
+ box.vertex_min[0] = box_p_min(0);
+ box.vertex_min[1] = box_p_min(1);
+ box.vertex_min[2] = box_p_min(2);
+ box.vertex_max[0] = box_p(0);
+ box.vertex_max[1] = box_p(1);
+ box.vertex_max[2] = box_p(2);
+ if (!check_box(cur_pose, axis, theta, depth, box)) u_min = u_min + pow(2,k);
+ k = k-1;
+ }
+ k = max_ulogN;
+ u_max = 0;
+ while (k>=0){
+ box_p = box_p_max - plane_u * box_length * (u_max + pow(2,k)) - plane_v * box_length - plane_w * box_length;
+ box.vertex_min[0] = box_p(0);
+ box.vertex_min[1] = box_p(1);
+ box.vertex_min[2] = box_p(2);
+ box.vertex_max[0] = box_p_max(0);
+ box.vertex_max[1] = box_p_max(1);
+ box.vertex_max[2] = box_p_max(2);
+ if (!check_box(cur_pose, axis, theta, depth, box)) {
+ u_max = u_max + pow(2,k);
+ // printf("-------- Not Included: (%0.3f,%0.3f,%0.3f),(%0.3f,%0.3f,%0.3f)\n",box.vertex_min[0],box.vertex_min[1],box.vertex_min[2],box.vertex_max[0],box.vertex_max[1],box.vertex_max[2]);
+ }
+
+ k = k-1;
+ }
+ u_max = max(0, max_uN - u_max - 1);
+ //printf("---- u_min -> u_max: %d->%d\n",u_min,u_max);
+ box_found = 0;
+ for (k = u_min; k <= u_max; k++){
+ box_p = box_p_min + plane_u * box_length * k;
+ box.vertex_min[0] = box_p(0);
+ box.vertex_min[1] = box_p(1);
+ box.vertex_min[2] = box_p(2);
+ box.vertex_max[0] = box_p(0) + box_length;
+ box.vertex_max[1] = box_p(1) + box_length;
+ box.vertex_max[2] = box_p(2) + box_length;
+ if (check_box_in_env(box)){
+ //printf("---- FOUND: (%0.3f,%0.3f,%0.3f),(%0.3f,%0.3f,%0.3f)\n",box.vertex_min[0],box.vertex_min[1],box.vertex_min[2],box.vertex_max[0],box.vertex_max[1],box.vertex_max[2]);
+ box_found = 1;
+ boxes.push_back(box);
+ }
+ }
+ if (box_found) {
+ flag = true;
+ } else {
+ if (j>1) break;
+ }
+ }
+ if (!flag && i>0) break;
+ }
+}
+
+bool FOV_Checker::check_box(Eigen::Vector3d cur_pose, Eigen::Vector3d axis, double theta, double depth, const BoxPointType box){
+ Eigen::Vector3d vertex[8];
+ bool s;
+ vertex[0] = Eigen::Vector3d(box.vertex_min[0], box.vertex_min[1], box.vertex_min[2]);
+ vertex[1] = Eigen::Vector3d(box.vertex_min[0], box.vertex_min[1], box.vertex_max[2]);
+ vertex[2] = Eigen::Vector3d(box.vertex_min[0], box.vertex_max[1], box.vertex_min[2]);
+ vertex[3] = Eigen::Vector3d(box.vertex_min[0], box.vertex_max[1], box.vertex_max[2]);
+ vertex[4] = Eigen::Vector3d(box.vertex_max[0], box.vertex_min[1], box.vertex_min[2]);
+ vertex[5] = Eigen::Vector3d(box.vertex_max[0], box.vertex_min[1], box.vertex_max[2]);
+ vertex[6] = Eigen::Vector3d(box.vertex_max[0], box.vertex_max[1], box.vertex_min[2]);
+ vertex[7] = Eigen::Vector3d(box.vertex_max[0], box.vertex_max[1], box.vertex_max[2]);
+ for (int i = 0; i < 8; i++){
+ if (check_point(cur_pose, axis, theta, depth, vertex[i])){
+ return true;
+ }
+ }
+ Eigen::Vector3d center_point = (vertex[7]+vertex[0])/2.0;
+ if (check_point(cur_pose, axis, theta, depth, center_point)){
+ return true;
+ }
+ PlaneType plane[6];
+ plane[0].p[0] = vertex[0];
+ plane[0].p[1] = vertex[2];
+ plane[0].p[2] = vertex[1];
+ plane[0].p[3] = vertex[3];
+
+ plane[1].p[0] = vertex[0];
+ plane[1].p[1] = vertex[4];
+ plane[1].p[2] = vertex[2];
+ plane[1].p[3] = vertex[6];
+
+ plane[2].p[0] = vertex[0];
+ plane[2].p[1] = vertex[4];
+ plane[2].p[2] = vertex[1];
+ plane[2].p[3] = vertex[5];
+
+ plane[3].p[0] = vertex[4];
+ plane[3].p[1] = vertex[6];
+ plane[3].p[2] = vertex[5];
+ plane[3].p[3] = vertex[7];
+
+ plane[4].p[0] = vertex[2];
+ plane[4].p[1] = vertex[6];
+ plane[4].p[2] = vertex[3];
+ plane[4].p[3] = vertex[7];
+
+ plane[5].p[0] = vertex[1];
+ plane[5].p[1] = vertex[5];
+ plane[5].p[2] = vertex[3];
+ plane[5].p[3] = vertex[7];
+ if (check_surface(cur_pose, axis, theta, depth, plane[0]) || check_surface(cur_pose, axis, theta, depth, plane[1]) || check_surface(cur_pose, axis, theta, depth, plane[2]) || check_surface(cur_pose, axis, theta, depth, plane[3]) || check_surface(cur_pose, axis, theta, depth, plane[4]) || check_surface(cur_pose, axis, theta, depth, plane[5]))
+ s = 1;
+ else
+ s = 0;
+ return s;
+}
+
+bool FOV_Checker::check_surface(Eigen::Vector3d cur_pose, Eigen::Vector3d axis, double theta, double depth, PlaneType plane){
+ Eigen::Vector3d plane_p, plane_u, plane_v, plane_w, pc, p, vec;
+ bool s;
+ double t, vec_dot_u, vec_dot_v;
+ plane_p = plane.p[0];
+ plane_u = plane.p[1] - plane_p;
+ plane_v = plane.p[2] - plane_p;
+ if (check_line(cur_pose, axis, theta, depth, plane_p, plane_u) || check_line(cur_pose, axis, theta, depth, plane_p, plane_v) || check_line(cur_pose, axis, theta, depth, plane_p + plane_u, plane_v) || check_line(cur_pose, axis, theta, depth, plane_p + plane_v, plane_u)){
+ s = 1;
+ return s;
+ }
+ pc = plane_p + (plane.p[3]-plane.p[0])/2;
+ if (check_point(cur_pose, axis, theta, depth, pc)){
+ s = 1;
+ return s;
+ }
+ plane_w = plane_u.cross(plane_v);
+ p = plane_p - cur_pose;
+ t = (p.dot(plane_w))/(axis.dot(plane_w));
+ vec = cur_pose + t * axis - plane_p;
+ vec_dot_u = vec.dot(plane_u)/plane_u.norm();
+ vec_dot_v = vec.dot(plane_v)/plane_v.norm();
+ if (t>=-eps_value && t<=depth && vec_dot_u>=-eps_value && vec_dot_u<=plane_u.norm() && vec_dot_v>=-eps_value && vec_dot_v <= plane_v.norm())
+ s = 1;
+ else
+ s = 0;
+ return s;
+}
+
+bool FOV_Checker::check_line(Eigen::Vector3d cur_pose, Eigen::Vector3d axis, double theta, double depth, Eigen::Vector3d line_p, Eigen::Vector3d line_vec){
+ Eigen::Vector3d p, vec_1, vec_2;
+ double xl, yl, zl, xn, yn, zn, dot_1, dot_2, ln, pn, pl, l2, p2;
+ double A, B, C, delta, t1, t2;
+ bool s;
+ p = line_p - cur_pose;
+ xl = line_vec(0); yl = line_vec(1); zl = line_vec(2);
+ xn = axis(0); yn = axis(1); zn = axis(2);
+ vec_1 = line_p - cur_pose;
+ vec_2 = line_p + line_vec - cur_pose;
+ dot_1 = vec_1.dot(axis);
+ dot_2 = vec_2.dot(axis);
+ //printf("xl yl zl: %0.4f, %0.4f, %0.4f\n", xl, yl, zl);
+ //printf("xn yn zn: %0.4f, %0.4f, %0.4f\n", xn, yn, zn);
+ //printf("dot_1, dot_2, %0.4f, %0.4f\n",dot_1, dot_2);
+ if ((dot_1<0 && dot_2<0) || (dot_1>depth && dot_2>depth)){
+ s = false;
+ return s;
+ }
+ ln = xl*xn+yl*yn+zl*zn;
+ pn = p(0)*xn+p(1)*yn+p(2)*zn;
+ pl = p(0)*xl+p(1)*yl+p(2)*zl;
+ l2 = xl*xl+yl*yl+zl*zl;
+ p2 = p.norm()*p.norm();
+ //printf("ln: %0.4f\n",ln);
+ //printf("pn:%0.4f\n",pn);
+ //printf("pl:%0.4f\n",pl);
+ //printf("l2:%0.4f\n",l2);
+ //printf("p2:%0.4f\n",p2);
+ //printf("theta, cos(theta):%0.4f %0.4f\n",theta,cos(theta));
+ A = ln * ln - l2 * cos(theta) * cos(theta);
+ B = 2 * pn * ln - 2 * cos(theta) * cos(theta)*pl;
+ C = pn * pn - p2 * cos(theta) * cos(theta);
+ //printf("A:%0.4f, B:%0.4f, C:%0.4f\n", A,B,C);
+ if (!(fabs(A)<=eps_value)){
+ delta = B*B - 4*A*C;
+ //printf("delta: %0.4f\n",delta);
+ if (delta <= eps_value){
+ if (A < -eps_value){
+ s = false;
+ return s;
+ }
+ else{
+ s = true;
+ return s;
+ }
+ } else {
+ double sqrt_delta = sqrt(delta);
+ t1 = (-B - sqrt_delta)/(2*A);
+ t2 = (-B + sqrt_delta)/(2*A);
+ if (t1>t2) swap(t1,t2);
+ //printf("t1,t2: %0.4f,%0.4f\n",t1,t2);
+ // printf("%d\n",check_point(cur_pose, axis, theta, depth, line_p + line_vec * t1));
+ if ((t1>=-eps_value && t1<=1+eps_value) && check_point(cur_pose, axis, theta, depth, line_p + line_vec * t1)){
+ s = true;
+ return s;
+ }
+ // printf("%d\n",check_point(cur_pose, axis, theta, depth, line_p + line_vec * t2));
+ if ((t2>=-eps_value && t2<=1+eps_value) && check_point(cur_pose, axis, theta, depth, line_p + line_vec * t2)){
+ s = true;
+ return s;
+ }
+ if (A>-eps_value && (t21-eps_value)){
+ s = true;
+ return s;
+ }
+ if (A1-eps_value){
+ s = true;
+ return s;
+ }
+ s = false;
+ }
+ } else{
+ if (!(fabs(B)<=eps_value)){
+ s = (B>-eps_value && -C/B<=1+eps_value) || (B=-eps_value);
+ return s;
+ }
+ else {
+ s = C>=-eps_value;
+ return s;
+ }
+ }
+ return false;
+}
+
+bool FOV_Checker::check_point(Eigen::Vector3d cur_pose, Eigen::Vector3d axis, double theta, double depth, Eigen::Vector3d point){
+ Eigen::Vector3d vec;
+ double proj_len;
+ bool s;
+ vec = point-cur_pose;
+ if (vec.transpose()*vec < 0.4 * box_length * box_length){
+ return true;
+ }
+ proj_len = vec.dot(axis);
+ if (proj_len > depth){
+ s = false;
+ return s;
+ }
+ //printf("acos: %0.4f\n",acos(proj_len/vec.norm()));
+ if (fabs(vec.norm()) <= 1e-4 || acos(proj_len/vec.norm()) <= theta + 0.0175)
+ s = true;
+ else
+ s = false;
+ return s;
+}
+
+bool FOV_Checker::check_box_in_env(BoxPointType box){
+ if (box.vertex_min[0] >= env.vertex_min[0]-eps_value && box.vertex_min[1] >= env.vertex_min[1]-eps_value && box.vertex_min[2] >= env.vertex_min[2]-eps_value && box.vertex_max[0]<= env.vertex_max[0]+eps_value && box.vertex_max[1]<= env.vertex_max[1]+eps_value && box.vertex_max[2]<= env.vertex_max[2]+eps_value){
+ return true;
+ } else {
+ return false;
+ }
+}
+
diff --git a/include/FOV_Checker/FOV_Checker.h b/include/FOV_Checker/FOV_Checker.h
new file mode 100755
index 0000000..e7db390
--- /dev/null
+++ b/include/FOV_Checker/FOV_Checker.h
@@ -0,0 +1,33 @@
+// Include Files
+#pragma once
+#include
+#include
+#include "ikd-Tree/ikd_Tree.h"
+#include
+#include
+
+#define eps_value 1e-6
+
+struct PlaneType{
+ Eigen::Vector3d p[4];
+};
+
+class FOV_Checker{
+public:
+ FOV_Checker();
+ ~FOV_Checker();
+ void Set_Env(BoxPointType env_param);
+ void Set_BoxLength(double box_len_param);
+ void check_fov(Eigen::Vector3d cur_pose, Eigen::Vector3d axis, double theta, double depth, vector &boxes);
+ bool check_box(Eigen::Vector3d cur_pose, Eigen::Vector3d axis, double theta, double depth, const BoxPointType box);
+ bool check_line(Eigen::Vector3d cur_pose, Eigen::Vector3d axis, double theta, double depth, Eigen::Vector3d line_p, Eigen::Vector3d line_vec);
+ bool check_surface(Eigen::Vector3d cur_pose, Eigen::Vector3d axis, double theta, double depth, PlaneType plane);
+ bool check_point(Eigen::Vector3d cur_pose, Eigen::Vector3d axis, double theta, double depth, Eigen::Vector3d point);
+ bool check_box_in_env(BoxPointType box);
+private:
+ BoxPointType env;
+ double box_length;
+ FILE *fp;
+
+};
+
diff --git a/include/IKFoM/.gitignore b/include/IKFoM/.gitignore
new file mode 100644
index 0000000..259148f
--- /dev/null
+++ b/include/IKFoM/.gitignore
@@ -0,0 +1,32 @@
+# Prerequisites
+*.d
+
+# Compiled Object files
+*.slo
+*.lo
+*.o
+*.obj
+
+# Precompiled Headers
+*.gch
+*.pch
+
+# Compiled Dynamic libraries
+*.so
+*.dylib
+*.dll
+
+# Fortran module files
+*.mod
+*.smod
+
+# Compiled Static libraries
+*.lai
+*.la
+*.a
+*.lib
+
+# Executables
+*.exe
+*.out
+*.app
diff --git a/include/IKFoM/.vscode/settings.json b/include/IKFoM/.vscode/settings.json
new file mode 100644
index 0000000..d756850
--- /dev/null
+++ b/include/IKFoM/.vscode/settings.json
@@ -0,0 +1,5 @@
+{
+ "files.associations": {
+ "complex": "cpp"
+ }
+}
\ No newline at end of file
diff --git a/include/IKFoM/IKFoM_toolkit/.vscode/settings.json b/include/IKFoM/IKFoM_toolkit/.vscode/settings.json
new file mode 100755
index 0000000..b619d8a
--- /dev/null
+++ b/include/IKFoM/IKFoM_toolkit/.vscode/settings.json
@@ -0,0 +1,62 @@
+{
+ "files.associations": {
+ "iosfwd": "cpp",
+ "core": "cpp",
+ "cctype": "cpp",
+ "clocale": "cpp",
+ "cmath": "cpp",
+ "cstdarg": "cpp",
+ "cstddef": "cpp",
+ "cstdio": "cpp",
+ "cstdlib": "cpp",
+ "cstring": "cpp",
+ "ctime": "cpp",
+ "cwchar": "cpp",
+ "cwctype": "cpp",
+ "array": "cpp",
+ "atomic": "cpp",
+ "strstream": "cpp",
+ "*.tcc": "cpp",
+ "bitset": "cpp",
+ "chrono": "cpp",
+ "complex": "cpp",
+ "cstdint": "cpp",
+ "deque": "cpp",
+ "list": "cpp",
+ "unordered_map": "cpp",
+ "vector": "cpp",
+ "exception": "cpp",
+ "fstream": "cpp",
+ "functional": "cpp",
+ "initializer_list": "cpp",
+ "iomanip": "cpp",
+ "iostream": "cpp",
+ "istream": "cpp",
+ "limits": "cpp",
+ "new": "cpp",
+ "ostream": "cpp",
+ "numeric": "cpp",
+ "ratio": "cpp",
+ "sstream": "cpp",
+ "stdexcept": "cpp",
+ "streambuf": "cpp",
+ "system_error": "cpp",
+ "thread": "cpp",
+ "cfenv": "cpp",
+ "cinttypes": "cpp",
+ "tuple": "cpp",
+ "type_traits": "cpp",
+ "utility": "cpp",
+ "typeinfo": "cpp",
+ "algorithm": "cpp",
+ "iterator": "cpp",
+ "map": "cpp",
+ "memory": "cpp",
+ "memory_resource": "cpp",
+ "optional": "cpp",
+ "random": "cpp",
+ "set": "cpp",
+ "string": "cpp",
+ "string_view": "cpp"
+ }
+}
\ No newline at end of file
diff --git a/include/IKFoM/IKFoM_toolkit/esekfom/.vscode/settings.json b/include/IKFoM/IKFoM_toolkit/esekfom/.vscode/settings.json
new file mode 100755
index 0000000..c807bd5
--- /dev/null
+++ b/include/IKFoM/IKFoM_toolkit/esekfom/.vscode/settings.json
@@ -0,0 +1,5 @@
+{
+ "files.associations": {
+ "core": "cpp"
+ }
+}
\ No newline at end of file
diff --git a/include/IKFoM/IKFoM_toolkit/esekfom/esekfom.hpp b/include/IKFoM/IKFoM_toolkit/esekfom/esekfom.hpp
new file mode 100755
index 0000000..d4524e1
--- /dev/null
+++ b/include/IKFoM/IKFoM_toolkit/esekfom/esekfom.hpp
@@ -0,0 +1,390 @@
+/*
+ * Copyright (c) 2019--2023, The University of Hong Kong
+ * All rights reserved.
+ *
+ * Author: Dongjiao HE
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * * Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * * Redistributions in binary form must reproduce the above
+ * copyright notice, this list of conditions and the following
+ * disclaimer in the documentation and/or other materials provided
+ * with the distribution.
+ * * Neither the name of the Universitaet Bremen nor the names of its
+ * contributors may be used to endorse or promote products derived
+ * from this software without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+ * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+ * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+ * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+ * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
+ * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#ifndef ESEKFOM_EKF_HPP
+#define ESEKFOM_EKF_HPP
+
+
+#include
+#include
+
+#include
+#include
+#include
+#include
+#include
+#include
+
+#include "../mtk/types/vect.hpp"
+#include "../mtk/types/SOn.hpp"
+#include "../mtk/types/S2.hpp"
+#include "../mtk/types/SEn.hpp"
+#include "../mtk/startIdx.hpp"
+#include "../mtk/build_manifold.hpp"
+#include "util.hpp"
+
+namespace esekfom {
+
+using namespace Eigen;
+
+template
+struct dyn_share_modified
+{
+ bool valid;
+ bool converge;
+ T M_Noise;
+ Eigen::Matrix z;
+ Eigen::Matrix h_x;
+ Eigen::Matrix z_IMU;
+ Eigen::Matrix R_IMU;
+ bool satu_check[6];
+};
+
+template
+class esekf{
+
+ typedef esekf self;
+ enum{
+ n = state::DOF, m = state::DIM, l = measurement::DOF
+ };
+
+public:
+
+ typedef typename state::scalar scalar_type;
+ typedef Matrix cov;
+ typedef Matrix cov_;
+ typedef SparseMatrix spMt;
+ typedef Matrix vectorized_state;
+ typedef Matrix flatted_state;
+ typedef flatted_state processModel(state &, const input &);
+ typedef Eigen::Matrix processMatrix1(state &, const input &);
+ typedef Eigen::Matrix processMatrix2(state &, const input &);
+ typedef Eigen::Matrix processnoisecovariance;
+
+ typedef void measurementModel_dyn_share_modified(state &, dyn_share_modified &);
+ typedef Eigen::Matrix measurementMatrix1(state &);
+ typedef Eigen::Matrix measurementMatrix1_dyn(state &);
+ typedef Eigen::Matrix measurementMatrix2(state &);
+ typedef Eigen::Matrix measurementMatrix2_dyn(state &);
+ typedef Eigen::Matrix measurementnoisecovariance;
+ typedef Eigen::Matrix measurementnoisecovariance_dyn;
+
+ esekf(const state &x = state(),
+ const cov &P = cov::Identity()): x_(x), P_(P){};
+
+ void init_dyn_share_modified(processModel f_in, processMatrix1 f_x_in, measurementModel_dyn_share_modified h_dyn_share_in)
+ {
+ f = f_in;
+ f_x = f_x_in;
+ // f_w = f_w_in;
+ h_dyn_share_modified_1 = h_dyn_share_in;
+ maximum_iter = 1;
+ x_.build_S2_state();
+ x_.build_SO3_state();
+ x_.build_vect_state();
+ x_.build_SEN_state();
+ }
+
+ void init_dyn_share_modified_2h(processModel f_in, processMatrix1 f_x_in, measurementModel_dyn_share_modified h_dyn_share_in1, measurementModel_dyn_share_modified h_dyn_share_in2)
+ {
+ f = f_in;
+ f_x = f_x_in;
+ // f_w = f_w_in;
+ h_dyn_share_modified_1 = h_dyn_share_in1;
+ h_dyn_share_modified_2 = h_dyn_share_in2;
+ maximum_iter = 1;
+ x_.build_S2_state();
+ x_.build_SO3_state();
+ x_.build_vect_state();
+ x_.build_SEN_state();
+ }
+
+ // iterated error state EKF propogation
+ void predict(double &dt, processnoisecovariance &Q, const input &i_in, bool predict_state, bool prop_cov){
+ if (predict_state)
+ {
+ flatted_state f_ = f(x_, i_in);
+ x_.oplus(f_, dt);
+ }
+
+ if (prop_cov)
+ {
+ flatted_state f_ = f(x_, i_in);
+ // state x_before = x_;
+
+ cov_ f_x_ = f_x(x_, i_in);
+ cov f_x_final;
+ F_x1 = cov::Identity();
+ for (std::vector, int> >::iterator it = x_.vect_state.begin(); it != x_.vect_state.end(); it++) {
+ int idx = (*it).first.first;
+ int dim = (*it).first.second;
+ int dof = (*it).second;
+ for(int i = 0; i < n; i++){
+ for(int j=0; j res_temp_SO3;
+ MTK::vect<3, scalar_type> seg_SO3;
+ for (std::vector >::iterator it = x_.SO3_state.begin(); it != x_.SO3_state.end(); it++) {
+ int idx = (*it).first;
+ int dim = (*it).second;
+ for(int i = 0; i < 3; i++){
+ seg_SO3(i) = -1 * f_(dim + i) * dt;
+ }
+ MTK::SO3 res;
+ res.w() = MTK::exp(res.vec(), seg_SO3, scalar_type(1/2));
+ F_x1.template block<3, 3>(idx, idx) = res.normalized().toRotationMatrix();
+ res_temp_SO3 = MTK::A_matrix(seg_SO3);
+ for(int i = 0; i < n; i++){
+ f_x_final. template block<3, 1>(idx, i) = res_temp_SO3 * (f_x_. template block<3, 1>(dim, i));
+ }
+ }
+
+ F_x1 += f_x_final * dt;
+ P_ = F_x1 * P_ * (F_x1).transpose() + Q * (dt * dt);
+ }
+ }
+
+ bool update_iterated_dyn_share_modified() {
+ dyn_share_modified dyn_share;
+ state x_propagated = x_;
+ int dof_Measurement;
+ double m_noise;
+ for(int i=0; i z = dyn_share.z;
+ // Matrix R = dyn_share.R;
+ Matrix h_x = dyn_share.h_x;
+ // Matrix h_v = dyn_share.h_v;
+ dof_Measurement = h_x.rows();
+ m_noise = dyn_share.M_Noise;
+ // dof_Measurement_noise = dyn_share.R.rows();
+ // vectorized_state dx, dx_new;
+ // x_.boxminus(dx, x_propagated);
+ // dx_new = dx;
+ // P_ = P_propagated;
+
+ Matrix PHT;
+ Matrix HPHT;
+ Matrix K_;
+ // if(n > dof_Measurement)
+ {
+ PHT = P_. template block(0, 0) * h_x.transpose();
+ HPHT = h_x * PHT.topRows(12);
+ for (int m = 0; m < dof_Measurement; m++)
+ {
+ HPHT(m, m) += m_noise;
+ }
+ K_= PHT*HPHT.inverse();
+ }
+ Matrix dx_ = K_ * z; // - h) + (K_x - Matrix::Identity()) * dx_new;
+ // state x_before = x_;
+
+ x_.boxplus(dx_);
+ dyn_share.converge = true;
+
+ // L_ = P_;
+ // Matrix res_temp_SO3;
+ // MTK::vect<3, scalar_type> seg_SO3;
+ // for(typename std::vector >::iterator it = x_.SO3_state.begin(); it != x_.SO3_state.end(); it++) {
+ // int idx = (*it).first;
+ // for(int i = 0; i < 3; i++){
+ // seg_SO3(i) = dx_(i + idx);
+ // }
+ // res_temp_SO3 = A_matrix(seg_SO3).transpose();
+ // for(int i = 0; i < n; i++){
+ // L_. template block<3, 1>(idx, i) = res_temp_SO3 * (P_. template block<3, 1>(idx, i));
+ // }
+ // {
+ // for(int i = 0; i < dof_Measurement; i++){
+ // K_. template block<3, 1>(idx, i) = res_temp_SO3 * (K_. template block<3, 1>(idx, i));
+ // }
+ // }
+ // for(int i = 0; i < n; i++){
+ // L_. template block<1, 3>(i, idx) = (L_. template block<1, 3>(i, idx)) * res_temp_SO3.transpose();
+ // // P_. template block<1, 3>(i, idx) = (P_. template block<1, 3>(i, idx)) * res_temp_SO3.transpose();
+ // }
+ // for(int i = 0; i < n; i++){
+ // P_. template block<1, 3>(i, idx) = (P_. template block<1, 3>(i, idx)) * res_temp_SO3.transpose();
+ // }
+ // }
+ // Matrix res_temp_S2;
+ // MTK::vect<2, scalar_type> seg_S2;
+ // for(typename std::vector >::iterator it = x_.S2_state.begin(); it != x_.S2_state.end(); it++) {
+ // int idx = (*it).first;
+
+ // for(int i = 0; i < 2; i++){
+ // seg_S2(i) = dx_(i + idx);
+ // }
+
+ // Eigen::Matrix Nx;
+ // Eigen::Matrix Mx;
+ // x_.S2_Nx_yy(Nx, idx);
+ // x_propagated.S2_Mx(Mx, seg_S2, idx);
+ // res_temp_S2 = Nx * Mx;
+
+ // for(int i = 0; i < n; i++){
+ // L_. template block<2, 1>(idx, i) = res_temp_S2 * (P_. template block<2, 1>(idx, i));
+ // }
+
+ // {
+ // for(int i = 0; i < dof_Measurement; i++){
+ // K_. template block<2, 1>(idx, i) = res_temp_S2 * (K_. template block<2, 1>(idx, i));
+ // }
+ // }
+ // for(int i = 0; i < n; i++){
+ // L_. template block<1, 2>(i, idx) = (L_. template block<1, 2>(i, idx)) * res_temp_S2.transpose();
+ // }
+ // for(int i = 0; i < n; i++){
+ // P_. template block<1, 2>(i, idx) = (P_. template block<1, 2>(i, idx)) * res_temp_S2.transpose();
+ // }
+ // }
+ // if(n > dof_Measurement)
+ {
+ P_ = P_ - K_*h_x*P_. template block<12, n>(0, 0);
+ }
+ }
+ return true;
+ }
+
+ void update_iterated_dyn_share_IMU() {
+
+ dyn_share_modified dyn_share;
+ for(int i=0; i z = dyn_share.z_IMU;
+
+ Matrix PHT;
+ Matrix HP;
+ Matrix HPHT;
+ PHT.setZero();
+ HP.setZero();
+ HPHT.setZero();
+ for (int l_ = 0; l_ < 6; l_++)
+ {
+ if (!dyn_share.satu_check[l_])
+ {
+ PHT.col(l_) = P_.col(15+l_) + P_.col(24+l_);
+ HP.row(l_) = P_.row(15+l_) + P_.row(24+l_);
+ }
+ }
+ for (int l_ = 0; l_ < 6; l_++)
+ {
+ if (!dyn_share.satu_check[l_])
+ {
+ HPHT.col(l_) = HP.col(15+l_) + HP.col(24+l_);
+ }
+ HPHT(l_, l_) += dyn_share.R_IMU(l_); //, l);
+ }
+ Eigen::Matrix K = PHT * HPHT.inverse();
+
+ Matrix dx_ = K * z;
+
+ P_ -= K * HP;
+ x_.boxplus(dx_);
+ }
+ return;
+ }
+
+ void change_x(state &input_state)
+ {
+ x_ = input_state;
+
+ if((!x_.vect_state.size())&&(!x_.SO3_state.size())&&(!x_.S2_state.size())&&(!x_.SEN_state.size()))
+ {
+ x_.build_S2_state();
+ x_.build_SO3_state();
+ x_.build_vect_state();
+ x_.build_SEN_state();
+ }
+ }
+
+ void change_P(cov &input_cov)
+ {
+ P_ = input_cov;
+ }
+
+ const state& get_x() const {
+ return x_;
+ }
+ const cov& get_P() const {
+ return P_;
+ }
+ state x_;
+private:
+ measurement m_;
+ cov P_;
+ spMt l_;
+ spMt f_x_1;
+ spMt f_x_2;
+ cov F_x1 = cov::Identity();
+ cov F_x2 = cov::Identity();
+ cov L_ = cov::Identity();
+
+ processModel *f;
+ processMatrix1 *f_x;
+ processMatrix2 *f_w;
+
+ measurementMatrix1 *h_x;
+ measurementMatrix2 *h_v;
+
+ measurementMatrix1_dyn *h_x_dyn;
+ measurementMatrix2_dyn *h_v_dyn;
+
+ measurementModel_dyn_share_modified *h_dyn_share_modified_1;
+
+ measurementModel_dyn_share_modified *h_dyn_share_modified_2;
+
+ int maximum_iter = 0;
+ scalar_type limit[n];
+
+public:
+ EIGEN_MAKE_ALIGNED_OPERATOR_NEW
+};
+
+} // namespace esekfom
+
+#endif // ESEKFOM_EKF_HPP
diff --git a/include/IKFoM/IKFoM_toolkit/esekfom/util.hpp b/include/IKFoM/IKFoM_toolkit/esekfom/util.hpp
new file mode 100755
index 0000000..ab39fc4
--- /dev/null
+++ b/include/IKFoM/IKFoM_toolkit/esekfom/util.hpp
@@ -0,0 +1,82 @@
+/*
+ * Copyright (c) 2019--2023, The University of Hong Kong
+ * All rights reserved.
+ *
+ * Author: Dongjiao HE
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * * Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * * Redistributions in binary form must reproduce the above
+ * copyright notice, this list of conditions and the following
+ * disclaimer in the documentation and/or other materials provided
+ * with the distribution.
+ * * Neither the name of the Universitaet Bremen nor the names of its
+ * contributors may be used to endorse or promote products derived
+ * from this software without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+ * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+ * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+ * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+ * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
+ * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ */
+
+#ifndef __MEKFOM_UTIL_HPP__
+#define __MEKFOM_UTIL_HPP__
+
+#include
+#include "../mtk/src/mtkmath.hpp"
+namespace esekfom {
+
+template
+class is_same {
+public:
+ operator bool() {
+ return false;
+ }
+};
+template
+class is_same {
+public:
+ operator bool() {
+ return true;
+ }
+};
+
+template
+class is_double {
+public:
+ operator bool() {
+ return false;
+ }
+};
+
+template<>
+class is_double {
+public:
+ operator bool() {
+ return true;
+ }
+};
+
+template
+static T
+id(const T &x)
+{
+ return x;
+}
+
+} // namespace esekfom
+
+#endif // __MEKFOM_UTIL_HPP__
diff --git a/include/IKFoM/IKFoM_toolkit/mtk/build_manifold.hpp b/include/IKFoM/IKFoM_toolkit/mtk/build_manifold.hpp
new file mode 100755
index 0000000..115ec7e
--- /dev/null
+++ b/include/IKFoM/IKFoM_toolkit/mtk/build_manifold.hpp
@@ -0,0 +1,248 @@
+// This is an advanced implementation of the algorithm described in the
+// following paper:
+// C. Hertzberg, R. Wagner, U. Frese, and L. Schroder. Integratinggeneric sensor fusion algorithms with sound state representationsthrough encapsulation of manifolds.
+// CoRR, vol. abs/1107.1119, 2011.[Online]. Available: http://arxiv.org/abs/1107.1119
+
+/*
+ * Copyright (c) 2019--2023, The University of Hong Kong
+ * All rights reserved.
+ *
+ * Modifier: Dongjiao HE
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * * Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * * Redistributions in binary form must reproduce the above
+ * copyright notice, this list of conditions and the following
+ * disclaimer in the documentation and/or other materials provided
+ * with the distribution.
+ * * Neither the name of the Universitaet Bremen nor the names of its
+ * contributors may be used to endorse or promote products derived
+ * from this software without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+ * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+ * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+ * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+ * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
+ * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ */
+
+/*
+ * Copyright (c) 2008--2011, Universitaet Bremen
+ * All rights reserved.
+ *
+ * Author: Christoph Hertzberg
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * * Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * * Redistributions in binary form must reproduce the above
+ * copyright notice, this list of conditions and the following
+ * disclaimer in the documentation and/or other materials provided
+ * with the distribution.
+ * * Neither the name of the Universitaet Bremen nor the names of its
+ * contributors may be used to endorse or promote products derived
+ * from this software without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+ * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+ * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+ * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+ * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
+ * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ */
+/**
+ * @file mtk/build_manifold.hpp
+ * @brief Macro to automatically construct compound manifolds.
+ *
+ */
+#ifndef MTK_AUTOCONSTRUCT_HPP_
+#define MTK_AUTOCONSTRUCT_HPP_
+
+#include
+
+#include
+#include
+#include
+
+#include "src/SubManifold.hpp"
+#include "startIdx.hpp"
+
+#ifndef PARSED_BY_DOXYGEN
+//////// internals //////
+
+#define MTK_APPLY_MACRO_ON_TUPLE(r, macro, tuple) macro tuple
+
+#define MTK_TRANSFORM_COMMA(macro, entries) BOOST_PP_SEQ_ENUM(BOOST_PP_SEQ_TRANSFORM_S(1, MTK_APPLY_MACRO_ON_TUPLE, macro, entries))
+
+#define MTK_TRANSFORM(macro, entries) BOOST_PP_SEQ_FOR_EACH_R(1, MTK_APPLY_MACRO_ON_TUPLE, macro, entries)
+
+#define MTK_CONSTRUCTOR_ARG( type, id) const type& id = type()
+#define MTK_CONSTRUCTOR_COPY( type, id) id(id)
+#define MTK_BOXPLUS( type, id) id.boxplus(MTK::subvector(__vec, &self::id), __scale);
+#define MTK_OPLUS( type, id) id.oplus(MTK::subvector_(__vec, &self::id), __scale);
+#define MTK_BOXMINUS( type, id) id.boxminus(MTK::subvector(__res, &self::id), __oth.id);
+#define MTK_HAT( type, id) if(id.IDX == idx){id.hat(vec, res);}
+#define MTK_JACOB_RIGHT_INV( type, id) if(id.IDX == idx){id.Jacob_right_inv(vec, res);}
+#define MTK_JACOB_RIGHT( type, id) if(id.IDX == idx){id.Jacob_right(vec, res);}
+#define MTK_S2_hat( type, id) if(id.IDX == idx){id.S2_hat(res);}
+#define MTK_S2_Nx_yy( type, id) if(id.IDX == idx){id.S2_Nx_yy(res);}
+#define MTK_S2_Mx( type, id) if(id.IDX == idx){id.S2_Mx(res, dx);}
+#define MTK_OSTREAM( type, id) << __var.id << " "
+#define MTK_ISTREAM( type, id) >> __var.id
+#define MTK_S2_state( type, id) if(id.TYP == 1){S2_state.push_back(std::make_pair(id.IDX, id.DIM));}
+#define MTK_SO3_state( type, id) if(id.TYP == 2){(SO3_state).push_back(std::make_pair(id.IDX, id.DIM));}
+#define MTK_vect_state( type, id) if(id.TYP == 0){(vect_state).push_back(std::make_pair(std::make_pair(id.IDX, id.DIM), type::DOF));}
+#define MTK_SEN_state( type, id) if(id.TYP == 4){(SEN_state).push_back(std::make_pair(std::make_pair(id.IDX, id.DIM), type::DOF));}
+
+#define MTK_SUBVARLIST(seq, S2state, SO3state, SENstate) \
+BOOST_PP_FOR_1( \
+ ( \
+ BOOST_PP_SEQ_SIZE(seq), \
+ BOOST_PP_SEQ_HEAD(seq), \
+ BOOST_PP_SEQ_TAIL(seq) (~), \
+ 0,\
+ 0,\
+ S2state,\
+ SO3state,\
+ SENstate ),\
+ MTK_ENTRIES_TEST, MTK_ENTRIES_NEXT, MTK_ENTRIES_OUTPUT)
+
+#define MTK_PUT_TYPE(type, id, dof, dim, S2state, SO3state, SENstate) \
+ MTK::SubManifold id;
+#define MTK_PUT_TYPE_AND_ENUM(type, id, dof, dim, S2state, SO3state, SENstate) \
+ MTK_PUT_TYPE(type, id, dof, dim, S2state, SO3state, SENstate) \
+ enum {DOF = type::DOF + dof}; \
+ enum {DIM = type::DIM+dim}; \
+ typedef type::scalar scalar;
+
+#define MTK_ENTRIES_OUTPUT(r, state) MTK_ENTRIES_OUTPUT_I state
+#define MTK_ENTRIES_OUTPUT_I(s, head, seq, dof, dim, S2state, SO3state, SENstate) \
+ MTK_APPLY_MACRO_ON_TUPLE(~, \
+ BOOST_PP_IF(BOOST_PP_DEC(s), MTK_PUT_TYPE, MTK_PUT_TYPE_AND_ENUM), \
+ ( BOOST_PP_TUPLE_REM_2 head, dof, dim, S2state, SO3state, SENstate))
+
+#define MTK_ENTRIES_TEST(r, state) MTK_TUPLE_ELEM_4_0 state
+
+//! this used to be BOOST_PP_TUPLE_ELEM_4_0:
+#define MTK_TUPLE_ELEM_4_0(a,b,c,d,e,f, g, h) a
+
+#define MTK_ENTRIES_NEXT(r, state) MTK_ENTRIES_NEXT_I state
+#define MTK_ENTRIES_NEXT_I(len, head, seq, dof, dim, S2state, SO3state, SENstate) ( \
+ BOOST_PP_DEC(len), \
+ BOOST_PP_SEQ_HEAD(seq), \
+ BOOST_PP_SEQ_TAIL(seq), \
+ dof + BOOST_PP_TUPLE_ELEM_2_0 head::DOF,\
+ dim + BOOST_PP_TUPLE_ELEM_2_0 head::DIM,\
+ S2state,\
+ SO3state,\
+ SENstate )
+
+#endif /* not PARSED_BY_DOXYGEN */
+
+
+/**
+ * Construct a manifold.
+ * @param name is the class-name of the manifold,
+ * @param entries is the list of sub manifolds
+ *
+ * Entries must be given in a list like this:
+ * @code
+ * typedef MTK::trafo > Pose;
+ * typedef MTK::vect Vec3;
+ * MTK_BUILD_MANIFOLD(imu_state,
+ * ((Pose, pose))
+ * ((Vec3, vel))
+ * ((Vec3, acc_bias))
+ * )
+ * @endcode
+ * Whitespace is optional, but the double parentheses are necessary.
+ * Construction is done entirely in preprocessor.
+ * After construction @a name is also a manifold. Its members can be
+ * accessed by names given in @a entries.
+ *
+ * @note Variable types are not allowed to have commas, thus types like
+ * @c vect need to be typedef'ed ahead.
+ */
+#define MTK_BUILD_MANIFOLD(name, entries) \
+struct name { \
+ typedef name self; \
+ std::vector > S2_state;\
+ std::vector > SO3_state;\
+ std::vector, int> > vect_state;\
+ std::vector, int> > SEN_state;\
+ MTK_SUBVARLIST(entries, S2_state, SO3_state, SEN_state) \
+ name ( \
+ MTK_TRANSFORM_COMMA(MTK_CONSTRUCTOR_ARG, entries) \
+ ) : \
+ MTK_TRANSFORM_COMMA(MTK_CONSTRUCTOR_COPY, entries) {}\
+ int getDOF() const { return DOF; } \
+ void boxplus(const MTK::vectview & __vec, scalar __scale = 1 ) { \
+ MTK_TRANSFORM(MTK_BOXPLUS, entries)\
+ } \
+ void oplus(const MTK::vectview & __vec, scalar __scale = 1 ) { \
+ MTK_TRANSFORM(MTK_OPLUS, entries)\
+ } \
+ void boxminus(MTK::vectview __res, const name& __oth) const { \
+ MTK_TRANSFORM(MTK_BOXMINUS, entries)\
+ } \
+ friend std::ostream& operator<<(std::ostream& __os, const name& __var){ \
+ return __os MTK_TRANSFORM(MTK_OSTREAM, entries); \
+ } \
+ void build_S2_state(){\
+ MTK_TRANSFORM(MTK_S2_state, entries)\
+ }\
+ void build_vect_state(){\
+ MTK_TRANSFORM(MTK_vect_state, entries)\
+ }\
+ void build_SO3_state(){\
+ MTK_TRANSFORM(MTK_SO3_state, entries)\
+ }\
+ void build_SEN_state(){\
+ MTK_TRANSFORM(MTK_SEN_state, entries)\
+ }\
+ void Lie_hat(Eigen::VectorXd &vec, Eigen::MatrixXd &res, int idx) {\
+ MTK_TRANSFORM(MTK_HAT, entries)\
+ }\
+ void Lie_Jacob_Right_Inv(Eigen::VectorXd &vec, Eigen::MatrixXd &res, int idx) {\
+ MTK_TRANSFORM(MTK_JACOB_RIGHT_INV, entries)\
+ }\
+ void Lie_Jacob_Right(Eigen::VectorXd &vec, Eigen::MatrixXd &res, int idx) {\
+ MTK_TRANSFORM(MTK_JACOB_RIGHT, entries)\
+ }\
+ void S2_hat(Eigen::Matrix &res, int idx) {\
+ MTK_TRANSFORM(MTK_S2_hat, entries)\
+ }\
+ void S2_Nx_yy(Eigen::Matrix &res, int idx) {\
+ MTK_TRANSFORM(MTK_S2_Nx_yy, entries)\
+ }\
+ void S2_Mx(Eigen::Matrix &res, Eigen::Matrix dx, int idx) {\
+ MTK_TRANSFORM(MTK_S2_Mx, entries)\
+ }\
+ friend std::istream& operator>>(std::istream& __is, name& __var){ \
+ return __is MTK_TRANSFORM(MTK_ISTREAM, entries); \
+ } \
+};
+
+
+
+#endif /*MTK_AUTOCONSTRUCT_HPP_*/
diff --git a/include/IKFoM/IKFoM_toolkit/mtk/src/SubManifold.hpp b/include/IKFoM/IKFoM_toolkit/mtk/src/SubManifold.hpp
new file mode 100755
index 0000000..a1b13de
--- /dev/null
+++ b/include/IKFoM/IKFoM_toolkit/mtk/src/SubManifold.hpp
@@ -0,0 +1,123 @@
+// This is an advanced implementation of the algorithm described in the
+// following paper:
+// C. Hertzberg, R. Wagner, U. Frese, and L. Schroder. Integratinggeneric sensor fusion algorithms with sound state representationsthrough encapsulation of manifolds.
+// CoRR, vol. abs/1107.1119, 2011.[Online]. Available: http://arxiv.org/abs/1107.1119
+
+/*
+ * Copyright (c) 2019--2023, The University of Hong Kong
+ * All rights reserved.
+ *
+ * Modifier: Dongjiao HE
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * * Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * * Redistributions in binary form must reproduce the above
+ * copyright notice, this list of conditions and the following
+ * disclaimer in the documentation and/or other materials provided
+ * with the distribution.
+ * * Neither the name of the Universitaet Bremen nor the names of its
+ * contributors may be used to endorse or promote products derived
+ * from this software without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+ * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+ * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+ * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+ * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
+ * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ */
+
+/*
+ * Copyright (c) 2008--2011, Universitaet Bremen
+ * All rights reserved.
+ *
+ * Author: Christoph Hertzberg
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * * Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * * Redistributions in binary form must reproduce the above
+ * copyright notice, this list of conditions and the following
+ * disclaimer in the documentation and/or other materials provided
+ * with the distribution.
+ * * Neither the name of the Universitaet Bremen nor the names of its
+ * contributors may be used to endorse or promote products derived
+ * from this software without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+ * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+ * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+ * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+ * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
+ * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ */
+/**
+ * @file mtk/src/SubManifold.hpp
+ * @brief Defines the SubManifold class
+ */
+
+
+#ifndef SUBMANIFOLD_HPP_
+#define SUBMANIFOLD_HPP_
+
+
+#include "vectview.hpp"
+
+
+namespace MTK {
+
+/**
+ * @ingroup SubManifolds
+ * Helper class for compound manifolds.
+ * This class wraps a manifold T and provides an enum IDX refering to the
+ * index of the SubManifold within the compound manifold.
+ *
+ * Memberpointers to a submanifold can be used for @ref SubManifolds "functions accessing submanifolds".
+ *
+ * @tparam T The manifold type of the sub-type
+ * @tparam idx The index of the sub-type within the compound manifold
+ */
+template
+struct SubManifold : public T
+{
+ enum {IDX = idx, DIM = dim /*!< index of the sub-type within the compound manifold */ };
+ //! manifold type
+ typedef T type;
+
+ //! Construct from derived type
+ template
+ explicit
+ SubManifold(const X& t) : T(t) {};
+
+ //! Construct from internal type
+ //explicit
+ SubManifold(const T& t) : T(t) {};
+
+ //! inherit assignment operator
+ using T::operator=;
+
+};
+
+} // namespace MTK
+
+
+#endif /* SUBMANIFOLD_HPP_ */
diff --git a/include/IKFoM/IKFoM_toolkit/mtk/src/mtkmath.hpp b/include/IKFoM/IKFoM_toolkit/mtk/src/mtkmath.hpp
new file mode 100755
index 0000000..a8770e4
--- /dev/null
+++ b/include/IKFoM/IKFoM_toolkit/mtk/src/mtkmath.hpp
@@ -0,0 +1,294 @@
+// This is an advanced implementation of the algorithm described in the
+// following paper:
+// C. Hertzberg, R. Wagner, U. Frese, and L. Schroder. Integratinggeneric sensor fusion algorithms with sound state representationsthrough encapsulation of manifolds.
+// CoRR, vol. abs/1107.1119, 2011.[Online]. Available: http://arxiv.org/abs/1107.1119
+
+/*
+ * Copyright (c) 2019--2023, The University of Hong Kong
+ * All rights reserved.
+ *
+ * Modifier: Dongjiao HE
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * * Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * * Redistributions in binary form must reproduce the above
+ * copyright notice, this list of conditions and the following
+ * disclaimer in the documentation and/or other materials provided
+ * with the distribution.
+ * * Neither the name of the Universitaet Bremen nor the names of its
+ * contributors may be used to endorse or promote products derived
+ * from this software without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+ * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+ * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+ * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+ * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
+ * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ */
+
+/*
+ * Copyright (c) 2008--2011, Universitaet Bremen
+ * All rights reserved.
+ *
+ * Author: Christoph Hertzberg
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * * Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * * Redistributions in binary form must reproduce the above
+ * copyright notice, this list of conditions and the following
+ * disclaimer in the documentation and/or other materials provided
+ * with the distribution.
+ * * Neither the name of the Universitaet Bremen nor the names of its
+ * contributors may be used to endorse or promote products derived
+ * from this software without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+ * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+ * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+ * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+ * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
+ * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ */
+/**
+ * @file mtk/src/mtkmath.hpp
+ * @brief several math utility functions.
+ */
+
+#ifndef MTKMATH_H_
+#define MTKMATH_H_
+
+#include
+
+#include
+
+#include "../types/vect.hpp"
+
+#ifndef M_PI
+#define M_PI 3.1415926535897932384626433832795
+#endif
+
+
+namespace MTK {
+
+namespace internal {
+
+template
+struct traits {
+ typedef typename Manifold::scalar scalar;
+ enum {DOF = Manifold::DOF};
+ typedef vect vectorized_type;
+ typedef Eigen::Matrix matrix_type;
+};
+
+template<>
+struct traits : traits > {};
+template<>
+struct traits : traits > {};
+
+} // namespace internal
+
+/**
+ * \defgroup MTKMath Mathematical helper functions
+ */
+//@{
+
+//! constant @f$ \pi @f$
+const double pi = M_PI;
+
+template inline scalar tolerance();
+
+template<> inline float tolerance() { return 1e-5f; }
+template<> inline double tolerance() { return 1e-11; }
+
+
+/**
+ * normalize @a x to @f$[-bound, bound] @f$.
+ *
+ * result for @f$ x = bound + 2\cdot n\cdot bound @f$ is arbitrary @f$\pm bound @f$.
+ */
+template
+inline scalar normalize(scalar x, scalar bound){ //not used
+ if(std::fabs(x) <= bound) return x;
+ int r = (int)(x *(scalar(1.0)/ bound));
+ return x - ((r + (r>>31) + 1) & ~1)*bound;
+}
+
+/**
+ * Calculate cosine and sinc of sqrt(x2).
+ * @param x2 the squared angle must be non-negative
+ * @return a pair containing cos and sinc of sqrt(x2)
+ */
+template
+std::pair cos_sinc_sqrt(const scalar &x2){
+ using std::sqrt;
+ using std::cos;
+ using std::sin;
+ static scalar const taylor_0_bound = boost::math::tools::epsilon();
+ static scalar const taylor_2_bound = sqrt(taylor_0_bound);
+ static scalar const taylor_n_bound = sqrt(taylor_2_bound);
+
+ assert(x2>=0 && "argument must be non-negative and must not be nan/-nan");
+
+ // FIXME check if bigger bounds are possible
+ if(x2>=taylor_n_bound) {
+ // slow fall-back solution
+ scalar x = sqrt(x2);
+ return std::make_pair(cos(x), sin(x)/x); // x is greater than 0.
+ }
+
+ // FIXME Replace by Horner-Scheme (4 instead of 5 FLOP/term, numerically more stable, theoretically cos and sinc can be calculated in parallel using SSE2 mulpd/addpd)
+ // TODO Find optimal coefficients using Remez algorithm
+ static scalar const inv[] = {1/3., 1/4., 1/5., 1/6., 1/7., 1/8., 1/9.};
+ scalar cosi = 1., sinc=1;
+ scalar term = -1/2. * x2;
+ for(int i=0; i<3; ++i) {
+ cosi += term;
+ term *= inv[2*i];
+ sinc += term;
+ term *= -inv[2*i+1] * x2;
+ }
+
+ return std::make_pair(cosi, sinc);
+
+}
+
+template
+Eigen::Matrix hat(const Base& v) {
+ Eigen::Matrix res;
+ res << 0, -v[2], v[1],
+ v[2], 0, -v[0],
+ -v[1], v[0], 0;
+ return res;
+}
+
+template
+Eigen::Matrix A_inv_trans(const Base& v){
+ Eigen::Matrix res;
+ if(v.norm() > MTK::tolerance())
+ {
+ res = Eigen::Matrix::Identity() + 0.5 * hat(v) + (1 - v.norm() * std::cos(v.norm() / 2) / 2 / std::sin(v.norm() / 2)) * hat(v) * hat(v) / v.squaredNorm();
+
+ }
+ else
+ {
+ res = Eigen::Matrix::Identity();
+ }
+
+ return res;
+}
+
+template
+Eigen::Matrix A_inv(const Base& v){
+ Eigen::Matrix res;
+ if(v.norm() > MTK::tolerance())
+ {
+ res = Eigen::Matrix::Identity() - 0.5 * hat(v) + (1 - v.norm() * std::cos(v.norm() / 2) / 2 / std::sin(v.norm() / 2)) * hat(v) * hat(v) / v.squaredNorm();
+
+ }
+ else
+ {
+ res = Eigen::Matrix::Identity();
+ }
+
+ return res;
+}
+
+template
+Eigen::Matrix S2_w_expw_( Eigen::Matrix v, scalar length)
+ {
+ Eigen::Matrix res;
+ scalar norm = std::sqrt(v[0]*v[0] + v[1]*v[1]);
+ if(norm < MTK::tolerance()){
+ res = Eigen::Matrix::Zero();
+ res(0, 1) = 1;
+ res(1, 2) = 1;
+ res /= length;
+ }
+ else{
+ res << -v[0]*(1/norm-1/std::tan(norm))/std::sin(norm), norm/std::sin(norm), 0,
+ -v[1]*(1/norm-1/std::tan(norm))/std::sin(norm), 0, norm/std::sin(norm);
+ res /= length;
+ }
+ }
+
+template
+Eigen::Matrix A_matrix(const Base & v){
+ Eigen::Matrix res;
+ double squaredNorm = v[0] * v[0] + v[1] * v[1] + v[2] * v[2];
+ double norm = std::sqrt(squaredNorm);
+ if(norm < MTK::tolerance()){
+ res = Eigen::Matrix::Identity();
+ }
+ else{
+ res = Eigen::Matrix::Identity() + (1 - std::cos(norm)) / squaredNorm * hat(v) + (1 - std::sin(norm) / norm) / squaredNorm * hat(v) * hat(v);
+ }
+ return res;
+}
+
+template
+scalar exp(vectview result, vectview vec, const scalar& scale = 1) {
+ scalar norm2 = vec.squaredNorm();
+ std::pair cos_sinc = cos_sinc_sqrt(scale*scale * norm2);
+ scalar mult = cos_sinc.second * scale;
+ result = mult * vec;
+ return cos_sinc.first;
+}
+
+
+/**
+ * Inverse function to @c exp.
+ *
+ * @param result @c vectview to the result
+ * @param w scalar part of input
+ * @param vec vector part of input
+ * @param scale scale result by this value
+ * @param plus_minus_periodicity if true values @f$[w, vec]@f$ and @f$[-w, -vec]@f$ give the same result
+ */
+template
+void log(vectview result,
+ const scalar &w, const vectview vec,
+ const scalar &scale, bool plus_minus_periodicity)
+{
+ // FIXME implement optimized case for vec.squaredNorm() <= tolerance() * (w*w) via Rational Remez approximation ~> only one division
+ scalar nv = vec.norm();
+ if(nv < tolerance()) {
+ if(!plus_minus_periodicity && w < 0) {
+ // find the maximal entry:
+ int i;
+ nv = vec.cwiseAbs().maxCoeff(&i);
+ result = scale * std::atan2(nv, w) * vect::Unit(i);
+ return;
+ }
+ nv = tolerance();
+ }
+ scalar s = scale / nv * (plus_minus_periodicity ? std::atan(nv / w) : std::atan2(nv, w) );
+
+ result = s * vec;
+}
+
+
+} // namespace MTK
+
+
+#endif /* MTKMATH_H_ */
diff --git a/include/IKFoM/IKFoM_toolkit/mtk/src/vectview.hpp b/include/IKFoM/IKFoM_toolkit/mtk/src/vectview.hpp
new file mode 100755
index 0000000..5025071
--- /dev/null
+++ b/include/IKFoM/IKFoM_toolkit/mtk/src/vectview.hpp
@@ -0,0 +1,168 @@
+
+/*
+ * Copyright (c) 2008--2011, Universitaet Bremen
+ * All rights reserved.
+ *
+ * Author: Christoph Hertzberg
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * * Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * * Redistributions in binary form must reproduce the above
+ * copyright notice, this list of conditions and the following
+ * disclaimer in the documentation and/or other materials provided
+ * with the distribution.
+ * * Neither the name of the Universitaet Bremen nor the names of its
+ * contributors may be used to endorse or promote products derived
+ * from this software without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+ * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+ * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+ * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+ * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
+ * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ */
+/**
+ * @file mtk/src/vectview.hpp
+ * @brief Wrapper class around a pointer used as interface for plain vectors.
+ */
+
+#ifndef VECTVIEW_HPP_
+#define VECTVIEW_HPP_
+
+#include
+
+namespace MTK {
+
+/**
+ * A view to a vector.
+ * Essentially, @c vectview is only a pointer to @c scalar but can be used directly in @c Eigen expressions.
+ * The dimension of the vector is given as template parameter and type-checked when used in expressions.
+ * Data has to be modifiable.
+ *
+ * @tparam scalar Scalar type of the vector.
+ * @tparam dim Dimension of the vector.
+ *
+ * @todo @c vectview can be replaced by simple inheritance of @c Eigen::Map, as soon as they get const-correct
+ */
+namespace internal {
+ template
+ struct CovBlock {
+ typedef typename Eigen::Block, T1::DOF, T2::DOF> Type;
+ typedef typename Eigen::Block, T1::DOF, T2::DOF> ConstType;
+ };
+
+ template
+ struct CovBlock_ {
+ typedef typename Eigen::Block, T1::DIM, T2::DIM> Type;
+ typedef typename Eigen::Block, T1::DIM, T2::DIM> ConstType;
+ };
+
+ template
+ struct CrossCovBlock {
+ typedef typename Eigen::Block, T1::DOF, T2::DOF> Type;
+ typedef typename Eigen::Block, T1::DOF, T2::DOF> ConstType;
+ };
+
+ template
+ struct CrossCovBlock_ {
+ typedef typename Eigen::Block, T1::DIM, T2::DIM> Type;
+ typedef typename Eigen::Block, T1::DIM, T2::DIM> ConstType;
+ };
+
+ template
+ struct VectviewBase {
+ typedef Eigen::Matrix matrix_type;
+ typedef typename matrix_type::MapType Type;
+ typedef typename matrix_type::ConstMapType ConstType;
+ };
+
+ template
+ struct UnalignedType {
+ typedef T type;
+ };
+}
+
+template
+class vectview : public internal::VectviewBase::Type {
+ typedef internal::VectviewBase VectviewBase;
+public:
+ //! plain matrix type
+ typedef typename VectviewBase::matrix_type matrix_type;
+ //! base type
+ typedef typename VectviewBase::Type base;
+ //! construct from pointer
+ explicit
+ vectview(scalar* data, int dim_=dim) : base(data, dim_) {}
+ //! construct from plain matrix
+ vectview(matrix_type& m) : base(m.data(), m.size()) {}
+ //! construct from another @c vectview
+ vectview(const vectview &v) : base(v) {}
+ //! construct from Eigen::Block:
+ template
+ vectview(Eigen::VectorBlock block) : base(&block.coeffRef(0), block.size()) {}
+ template
+ vectview(Eigen::Block block) : base(&block.coeffRef(0), block.size()) {}
+
+ //! inherit assignment operator
+ using base::operator=;
+ //! data pointer
+ scalar* data() {return const_cast(base::data());}
+};
+
+/**
+ * @c const version of @c vectview.
+ * Compared to @c Eigen::Map this implementation is const correct, i.e.,
+ * data will not be modifiable using this view.
+ *
+ * @tparam scalar Scalar type of the vector.
+ * @tparam dim Dimension of the vector.
+ *
+ * @sa vectview
+ */
+template
+class vectview : public internal::VectviewBase::ConstType {
+ typedef internal::VectviewBase VectviewBase;
+public:
+ //! plain matrix type
+ typedef typename VectviewBase::matrix_type matrix_type;
+ //! base type
+ typedef typename VectviewBase::ConstType base;
+ //! construct from const pointer
+ explicit
+ vectview(const scalar* data, int dim_ = dim) : base(data, dim_) {}
+ //! construct from column vector
+ template
+ vectview(const Eigen::Matrix& m) : base(m.data()) {}
+ //! construct from row vector
+ template
+ vectview(const Eigen::Matrix& m) : base(m.data()) {}
+ //! construct from another @c vectview
+ vectview(vectview x) : base(x.data()) {}
+ //! construct from base
+ vectview(const base &x) : base(x) {}
+ /**
+ * Construct from Block
+ * @todo adapt this, when Block gets const-correct
+ */
+ template
+ vectview(Eigen::VectorBlock block) : base(&block.coeffRef(0)) {}
+ template
+ vectview(Eigen::Block block) : base(&block.coeffRef(0)) {}
+
+};
+
+
+} // namespace MTK
+
+#endif /* VECTVIEW_HPP_ */
diff --git a/include/IKFoM/IKFoM_toolkit/mtk/startIdx.hpp b/include/IKFoM/IKFoM_toolkit/mtk/startIdx.hpp
new file mode 100755
index 0000000..4dc2958
--- /dev/null
+++ b/include/IKFoM/IKFoM_toolkit/mtk/startIdx.hpp
@@ -0,0 +1,328 @@
+// This is an advanced implementation of the algorithm described in the
+// following paper:
+// C. Hertzberg, R. Wagner, U. Frese, and L. Schroder. Integratinggeneric sensor fusion algorithms with sound state representationsthrough encapsulation of manifolds.
+// CoRR, vol. abs/1107.1119, 2011.[Online]. Available: http://arxiv.org/abs/1107.1119
+
+/*
+ * Copyright (c) 2019--2023, The University of Hong Kong
+ * All rights reserved.
+ *
+ * Modifier: Dongjiao HE
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * * Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * * Redistributions in binary form must reproduce the above
+ * copyright notice, this list of conditions and the following
+ * disclaimer in the documentation and/or other materials provided
+ * with the distribution.
+ * * Neither the name of the Universitaet Bremen nor the names of its
+ * contributors may be used to endorse or promote products derived
+ * from this software without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+ * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+ * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+ * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+ * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
+ * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ */
+
+/*
+ * Copyright (c) 2008--2011, Universitaet Bremen
+ * All rights reserved.
+ *
+ * Author: Christoph Hertzberg
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * * Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * * Redistributions in binary form must reproduce the above
+ * copyright notice, this list of conditions and the following
+ * disclaimer in the documentation and/or other materials provided
+ * with the distribution.
+ * * Neither the name of the Universitaet Bremen nor the names of its
+ * contributors may be used to endorse or promote products derived
+ * from this software without specific prior written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+ * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+ * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
+ * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
+ * COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
+ * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+ * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+ * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
+ * ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ */
+/**
+ * @file mtk/startIdx.hpp
+ * @brief Tools to access sub-elements of compound manifolds.
+ */
+#ifndef GET_START_INDEX_H_
+#define GET_START_INDEX_H_
+
+#include
+
+#include "src/SubManifold.hpp"
+#include "src/vectview.hpp"
+
+namespace MTK {
+
+
+/**
+ * \defgroup SubManifolds Accessing Submanifolds
+ * For compound manifolds constructed using MTK_BUILD_MANIFOLD, member pointers
+ * can be used to get sub-vectors or matrix-blocks of a corresponding big matrix.
+ * E.g. for a type @a pose consisting of @a orient and @a trans the member pointers
+ * @c &pose::orient and @c &pose::trans give all required information and are still
+ * valid if the base type gets extended or the actual types of @a orient and @a trans
+ * change (e.g. from 2D to 3D).
+ *
+ * @todo Maybe require manifolds to typedef MatrixType and VectorType, etc.
+ */
+//@{
+
+/**
+ * Determine the index of a sub-variable within a compound variable.
+ */
+template
+int getStartIdx( MTK::SubManifold Base::*)
+{
+ return idx;
+}
+
+template
+int getStartIdx_( MTK::SubManifold Base::*)
+{
+ return dim;
+}
+
+/**
+ * Determine the degrees of freedom of a sub-variable within a compound variable.
+ */
+template
+int getDof( MTK::SubManifold Base::*)
+{
+ return T::DOF;
+}
+template
+int getDim( MTK::SubManifold Base::*)
+{
+ return T::DIM;
+}
+
+/**
+ * set the diagonal elements of a covariance matrix corresponding to a sub-variable
+ */
+template
+void setDiagonal(Eigen::Matrix &cov,
+ MTK::SubManifold Base::*, const typename Base::scalar &val)
+{
+ cov.diagonal().template segment(idx).setConstant(val);
+}
+
+template
+void setDiagonal_(Eigen::Matrix &cov,
+ MTK::SubManifold Base::*, const typename Base::scalar &val)
+{
+ cov.diagonal().template segment(dim).setConstant(val);
+}
+
+/**
+ * Get the subblock of corresponding to two members, i.e.
+ * \code
+ * Eigen::Matrix m;
+ * MTK::subblock(m, &Pose::orient, &Pose::trans) = some_expression;
+ * MTK::subblock(m, &Pose::trans, &Pose::orient) = some_expression.trans();
+ * \endcode
+ * lets you modify mixed covariance entries in a bigger covariance matrix.
+ */
+template
+typename MTK::internal::CovBlock::Type
+subblock(Eigen::Matrix &cov,
+ MTK::SubManifold Base::*, MTK::SubManifold Base::*)
+{
+ return cov.template block(idx1, idx2);
+}
+
+template
+typename MTK::internal::CovBlock_::Type
+subblock_(Eigen::Matrix &cov,
+ MTK::SubManifold Base::*, MTK::SubManifold Base::*)
+{
+ return cov.template block(dim1, dim2);
+}
+
+template
+typename MTK::internal::CrossCovBlock::Type
+subblock(Eigen::Matrix &cov, MTK::SubManifold Base1::*, MTK::SubManifold Base2::*)
+{
+ return cov.template block(idx1, idx2);
+}
+
+template
+typename MTK::internal::CrossCovBlock_::Type
+subblock_(Eigen::Matrix &cov, MTK::SubManifold Base1::*, MTK::SubManifold Base2::*)
+{
+ return cov.template block(dim1, dim2);
+}
+/**
+ * Get the subblock of corresponding to a member, i.e.
+ * \code
+ * Eigen::Matrix m;
+ * MTK::subblock(m, &Pose::orient) = some_expression;
+ * \endcode
+ * lets you modify covariance entries in a bigger covariance matrix.
+ */
+template
+typename MTK::internal::CovBlock_::Type
+subblock_(Eigen::Matrix &cov,
+ MTK::SubManifold Base::*)
+{
+ return cov.template block(dim, dim);
+}
+
+template
+typename MTK::internal::CovBlock::Type
+subblock(Eigen::Matrix &cov,
+ MTK::SubManifold Base::*)
+{
+ return cov.template block(idx, idx);
+}
+
+template
+class get_cov {
+public:
+ typedef Eigen::Matrix type;
+ typedef const Eigen::Matrix const_type;
+};
+
+template
+class get_cov_ {
+public:
+ typedef Eigen::Matrix type;
+ typedef const Eigen::Matrix const_type;
+};
+
+template
+class get_cross_cov {
+public:
+ typedef Eigen::Matrix type;
+ typedef const type const_type;
+};
+
+template
+class get_cross_cov_ {
+public:
+ typedef Eigen::Matrix type;
+ typedef const type const_type;
+};
+
+
+template