diff --git a/.gitignore b/.gitignore
index b5e0a7a67..fb1ed5e75 100644
--- a/.gitignore
+++ b/.gitignore
@@ -7,3 +7,4 @@ doxgen_generated
*.swp
*.swo
doc
+.vscode
diff --git a/ReadMe.md b/ReadMe.md
index 63a0bd4f7..472f356c5 100644
--- a/ReadMe.md
+++ b/ReadMe.md
@@ -17,10 +17,11 @@ details on what the system supports.
* Github project page - https://github.com/rpng/open_vins
* Documentation - https://docs.openvins.com/
* Getting started guide - https://docs.openvins.com/getting-started.html
-* Publication reference - http://udel.edu/~pgeneva/downloads/papers/c10.pdf
+* Publication reference - https://pgeneva.com/downloads/papers/Geneva2020ICRA.pdf
## News / Events
+* **July 14, 2022** - Improved feature extraction logic for >100hz tracking, some bug fixes and updated scripts. See v2.6.1 [PR#259](https://github.com/rpng/open_vins/pull/259).
* **March 14, 2022** - Initial dynamic initialization open sourcing, asynchronous subscription to inertial readings and publishing of odometry, support for lower frequency feature tracking. See v2.6 [PR#232](https://github.com/rpng/open_vins/pull/232) for details.
* **December 13, 2021** - New YAML configuration system, ROS2 support, Docker images, robust static initialization based on disparity, internal logging system to reduce verbosity, image transport publishers, dynamic number of features support, and other small fixes. See
v2.5 [PR#209](https://github.com/rpng/open_vins/pull/209) for details.
diff --git a/config/euroc_mav/estimator_config.yaml b/config/euroc_mav/estimator_config.yaml
index 17d3a1029..98e422462 100644
--- a/config/euroc_mav/estimator_config.yaml
+++ b/config/euroc_mav/estimator_config.yaml
@@ -27,7 +27,7 @@ feat_rep_aruco: "ANCHORED_MSCKF_INVERSE_DEPTH"
# zero velocity update parameters we can use
# we support either IMU-based or disparity detection.
-try_zupt: true
+try_zupt: false
zupt_chi2_multipler: 0 # set to 0 for only disp-based
zupt_max_velocity: 0.1
zupt_noise_multiplier: 50
@@ -39,24 +39,25 @@ zupt_only_at_beginning: true
init_window_time: 2.0 # how many seconds to collect initialization information
init_imu_thresh: 1.5 # threshold for variance of the accelerometer to detect a "jerk" in motion
-init_max_disparity: 1.5 # max disparity to consider the platform stationary (dependent on resolution)
-init_max_features: 25 # how many features to track during initialization (saves on computation)
-
-init_dyn_mle_opt_calib: false
-init_dyn_mle_max_iter: 50
-init_dyn_mle_max_time: 0.05
-init_dyn_mle_max_threads: 6
-init_dyn_num_pose: 6
-init_dyn_min_deg: 10.0
-
-init_dyn_inflation_ori: 10
-init_dyn_inflation_vel: 100
-init_dyn_inflation_bg: 10
-init_dyn_inflation_ba: 100
-init_dyn_min_rec_cond: 1e-12
-
-init_dyn_bias_g: [0.0, 0.0, 0.0]
-init_dyn_bias_a: [0.0, 0.0, 0.0]
+init_max_disparity: 10.0 # max disparity to consider the platform stationary (dependent on resolution)
+init_max_features: 75 # how many features to track during initialization (saves on computation)
+
+init_dyn_use: false # if dynamic initialization should be used
+init_dyn_mle_opt_calib: false # if we should optimize calibration during intialization (not recommended)
+init_dyn_mle_max_iter: 50 # how many iterations the MLE refinement should use (zero to skip the MLE)
+init_dyn_mle_max_time: 0.05 # how many seconds the MLE should be completed in
+init_dyn_mle_max_threads: 6 # how many threads the MLE should use
+init_dyn_num_pose: 6 # number of poses to use within our window time (evenly spaced)
+init_dyn_min_deg: 10.0 # orientation change needed to try to init
+
+init_dyn_inflation_ori: 10 # what to inflate the recovered q_GtoI covariance by
+init_dyn_inflation_vel: 100 # what to inflate the recovered v_IinG covariance by
+init_dyn_inflation_bg: 10 # what to inflate the recovered bias_g covariance by
+init_dyn_inflation_ba: 100 # what to inflate the recovered bias_a covariance by
+init_dyn_min_rec_cond: 1e-12 # reciprocal condition number thresh for info inversion
+
+init_dyn_bias_g: [ 0.0, 0.0, 0.0 ] # initial gyroscope bias guess
+init_dyn_bias_a: [ 0.0, 0.0, 0.0 ] # initial accelerometer bias guess
# ==================================================================
# ==================================================================
@@ -76,14 +77,14 @@ filepath_gt: "/tmp/ov_groundtruth.txt"
# our front-end feature tracking parameters
# we have a KLT and descriptor based (KLT is better implemented...)
-use_klt: true
-num_pts: 200
-fast_threshold: 20
-grid_x: 20
-grid_y: 20
-min_px_dist: 15
-knn_ratio: 0.70
-track_frequency: 21.0
+use_klt: true # if true we will use KLT, otherwise use a ORB descriptor + robust matching
+num_pts: 200 # number of points (per camera) we will extract and try to track
+fast_threshold: 20 # threshold for fast extraction (warning: lower threshs can be expensive)
+grid_x: 5 # extraction sub-grid count for horizontal direction (uniform tracking)
+grid_y: 5 # extraction sub-grid count for vertical direction (uniform tracking)
+min_px_dist: 10 # distance between features (features near each other provide less information)
+knn_ratio: 0.70 # descriptor knn threshold for the top two descriptor matches
+track_frequency: 21.0 # frequency we will perform feature tracking at (in frames per second / hertz)
downsample_cameras: false # will downsample image in half if true
multi_threading: true # if should enable opencv multi threading
histogram_method: "HISTOGRAM" # NONE, HISTOGRAM, CLAHE
diff --git a/config/kaist/estimator_config.yaml b/config/kaist/estimator_config.yaml
index 82dfff5f2..90bc96735 100644
--- a/config/kaist/estimator_config.yaml
+++ b/config/kaist/estimator_config.yaml
@@ -9,17 +9,16 @@ use_rk4int: true # if rk4 integration should be used (overrides imu averaging)
use_stereo: true # if we have more than 1 camera, if we should try to track stereo constraints
max_cameras: 2
-calib_cam_extrinsics: false
+calib_cam_extrinsics: false # degenerate motion
calib_cam_intrinsics: true
-calib_cam_timeoffset: true
+calib_cam_timeoffset: true # degenerate motion
-max_clones: 12
+max_clones: 8
max_slam: 50
max_slam_in_update: 25
max_msckf_in_update: 50
dt_slam_delay: 1
-#gravity_mag: 9.79858
gravity_mag: 9.81
feat_rep_msckf: "ANCHORED_MSCKF_INVERSE_DEPTH"
@@ -29,7 +28,7 @@ feat_rep_aruco: "ANCHORED_MSCKF_INVERSE_DEPTH"
# zero velocity update parameters we can use
# we support either IMU-based or disparity detection.
try_zupt: true
-zupt_chi2_multipler: 1 # set to 0 for only disp-based
+zupt_chi2_multipler: 0.25 # set to 0 for only disp-based
zupt_max_velocity: 0.1
zupt_noise_multiplier: 5
zupt_max_disparity: 0.4 # set to 0 for only imu-based
@@ -39,16 +38,17 @@ zupt_only_at_beginning: false
# ==================================================================
init_window_time: 2.0
-init_imu_thresh: 0.5 #0.5
+init_imu_thresh: 0.5
init_max_disparity: 1.5
-init_max_features: 25
+init_max_features: 75
+init_dyn_use: true
init_dyn_mle_opt_calib: false
init_dyn_mle_max_iter: 50
init_dyn_mle_max_time: 0.05
init_dyn_mle_max_threads: 6
init_dyn_num_pose: 6
-init_dyn_min_deg: 0.0 # traj is mostly straight line motion
+init_dyn_min_deg: 5.0 # traj is mostly straight line motion
init_dyn_inflation_ori: 10
init_dyn_inflation_vel: 100
@@ -78,18 +78,19 @@ filepath_gt: "/tmp/ov_groundtruth.txt"
# we have a KLT and descriptor based (KLT is better implemented...)
use_klt: true
num_pts: 200
-fast_threshold: 50
-grid_x: 20
-grid_y: 15
-min_px_dist: 30
+fast_threshold: 30
+grid_x: 5
+grid_y: 5
+min_px_dist: 20
knn_ratio: 0.65
-track_frequency: 21.0
+track_frequency: 31.0
downsample_cameras: false
multi_threading: true
histogram_method: "HISTOGRAM" # NONE, HISTOGRAM, CLAHE
-fi_max_dist: 500
-fi_max_baseline: 800
+fi_min_dist: 0.25
+fi_max_dist: 150.0
+fi_max_baseline: 200
fi_max_cond_number: 20000
fi_triangulate_1d: false
@@ -103,9 +104,9 @@ downsize_aruco: true
# ==================================================================
# camera noises and chi-squared threshold multipliers
-up_msckf_sigma_px: 1
+up_msckf_sigma_px: 1.5
up_msckf_chi2_multipler: 1
-up_slam_sigma_px: 1
+up_slam_sigma_px: 1.5
up_slam_chi2_multipler: 1
up_aruco_sigma_px: 1
up_aruco_chi2_multipler: 1
diff --git a/config/kaist/kalibr_imu_chain.yaml b/config/kaist/kalibr_imu_chain.yaml
index b3eea522e..df2714f67 100644
--- a/config/kaist/kalibr_imu_chain.yaml
+++ b/config/kaist/kalibr_imu_chain.yaml
@@ -1,22 +1,23 @@
%YAML:1.0 # need to specify the file type at the top!
-# MTI-100 series converted from data sheet, guess on bias random walk
+# MTI-100 series converted from datasheet, guess on bias random walk
# https://www.xsens.com/hubfs/Downloads/usermanual/MTi_usermanual.pdf
+
imu0:
T_i_b:
- [1.0, 0.0, 0.0, 0.0]
- [0.0, 1.0, 0.0, 0.0]
- [0.0, 0.0, 1.0, 0.0]
- [0.0, 0.0, 0.0, 1.0]
- accelerometer_noise_density: 2.0000e-03 # [ m / s^2 / sqrt(Hz) ] ( accel "white noise" )
- accelerometer_random_walk: 3.0000e-03 # [ m / s^3 / sqrt(Hz) ]. ( accel bias diffusion )
- gyroscope_noise_density: 1.6968e-04 # [ rad / s / sqrt(Hz) ] ( gyro "white noise" )
- gyroscope_random_walk: 1.9393e-05 # [ rad / s^2 / sqrt(Hz) ] ( gyro bias diffusion )
- #
- #
- #
- #
+# accelerometer_noise_density: 2.0000e-03 # [ m / s^2 / sqrt(Hz) ] ( accel "white noise" )
+# accelerometer_random_walk: 3.0000e-03 # [ m / s^3 / sqrt(Hz) ]. ( accel bias diffusion )
+# gyroscope_noise_density: 1.6968e-04 # [ rad / s / sqrt(Hz) ] ( gyro "white noise" )
+# gyroscope_random_walk: 1.9393e-05 # [ rad / s^2 / sqrt(Hz) ] ( gyro bias diffusion )
+ accelerometer_noise_density: 5.8860e-03 # [ m / s^2 / sqrt(Hz) ] ( accel "white noise" )
+ accelerometer_random_walk: 1.0000e-04 # [ m / s^3 / sqrt(Hz) ]. ( accel bias diffusion )
+ gyroscope_noise_density: 1.7453e-04 # [ rad / s / sqrt(Hz) ] ( gyro "white noise" )
+ gyroscope_random_walk: 1.0000e-05 # [ rad / s^2 / sqrt(Hz) ] ( gyro bias diffusion )
model: calibrated
rostopic: /imu/data_raw
time_offset: 0.0
diff --git a/config/kaist/kalibr_imucam_chain.yaml b/config/kaist/kalibr_imucam_chain.yaml
index 5f5e715a3..172591ee0 100644
--- a/config/kaist/kalibr_imucam_chain.yaml
+++ b/config/kaist/kalibr_imucam_chain.yaml
@@ -13,6 +13,7 @@ cam0:
intrinsics: [8.1690378992770002e+02,8.1156803828490001e+02,6.0850726281690004e+02,2.6347599764440002e+02]
resolution: [1280, 560]
rostopic: /stereo/left/image_raw
+
cam1:
T_imu_cam: #rotation from camera to IMU R_CtoI, position of camera in IMU p_CinI
- [-0.00768,-0.01509,0.99986,1.73376]
@@ -25,4 +26,4 @@ cam1:
distortion_model: radtan
intrinsics: [8.1378205539589999e+02,8.0852165574269998e+02,6.1386419539320002e+02,2.4941049348650000e+02]
resolution: [1280, 560]
- rostopic: /stereo/right/image_raw
\ No newline at end of file
+ rostopic: /stereo/right/image_raw
diff --git a/config/kaist_vio/estimator_config.yaml b/config/kaist_vio/estimator_config.yaml
index b26db2d65..a23950973 100644
--- a/config/kaist_vio/estimator_config.yaml
+++ b/config/kaist_vio/estimator_config.yaml
@@ -16,18 +16,18 @@ calib_cam_timeoffset: false # disable since this is a degenerate motion
max_clones: 11 # how many clones in the sliding window
max_slam: 50 # number of features in our state vector
max_slam_in_update: 25 # update can be split into sequential updates of batches, how many in a batch
-max_msckf_in_update: 40 # how many MSCKF features to use in the update
+max_msckf_in_update: 50 # how many MSCKF features to use in the update
dt_slam_delay: 1 # delay before initializing (helps with stability from bad initialization...)
gravity_mag: 9.81 # magnitude of gravity in this location
-feat_rep_msckf: "GLOBAL_3D"
+feat_rep_msckf: "ANCHORED_MSCKF_INVERSE_DEPTH"
feat_rep_slam: "ANCHORED_MSCKF_INVERSE_DEPTH"
feat_rep_aruco: "ANCHORED_MSCKF_INVERSE_DEPTH"
# zero velocity update parameters we can use
# we support either IMU-based or disparity detection.
-try_zupt: true
+try_zupt: false
zupt_chi2_multipler: 0 # set to 0 for only disp-based
zupt_max_velocity: 0.02
zupt_noise_multiplier: 10
@@ -38,10 +38,11 @@ zupt_only_at_beginning: false
# ==================================================================
init_window_time: 2.0 # how many seconds to collect initialization information
-init_imu_thresh: 1.5 # threshold for variance of the accelerometer to detect a "jerk" in motion
-init_max_disparity: 0.20 # max disparity to consider the platform stationary (dependent on resolution)
-init_max_features: 25 # how many features to track during initialization (saves on computation)
+init_imu_thresh: 0.60 # threshold for variance of the accelerometer to detect a "jerk" in motion
+init_max_disparity: 5.0 # max disparity to consider the platform stationary (dependent on resolution)
+init_max_features: 75 # how many features to track during initialization (saves on computation)
+init_dyn_use: false
init_dyn_mle_opt_calib: false
init_dyn_mle_max_iter: 50
init_dyn_mle_max_time: 0.05
@@ -79,11 +80,11 @@ filepath_gt: "/tmp/ov_groundtruth.txt"
use_klt: true
num_pts: 200
fast_threshold: 30
-grid_x: 20
-grid_y: 20
+grid_x: 5
+grid_y: 5
min_px_dist: 15
knn_ratio: 0.70
-track_frequency: 10.0
+track_frequency: 31.0
downsample_cameras: false # will downsample image in half if true
multi_threading: true # if should enable opencv multi threading
histogram_method: "HISTOGRAM" # NONE, HISTOGRAM, CLAHE
diff --git a/config/kaist_vio/kalibr_imu_chain.yaml b/config/kaist_vio/kalibr_imu_chain.yaml
index ee06ed6df..b9a3fe5d1 100644
--- a/config/kaist_vio/kalibr_imu_chain.yaml
+++ b/config/kaist_vio/kalibr_imu_chain.yaml
@@ -6,14 +6,14 @@ imu0:
- [0.0, 1.0, 0.0, 0.0]
- [0.0, 0.0, 1.0, 0.0]
- [0.0, 0.0, 0.0, 1.0]
- accelerometer_noise_density: 0.00333388 # [ m / s^2 / sqrt(Hz) ] ( accel "white noise" )
- accelerometer_random_walk: 0.00047402 # [ m / s^3 / sqrt(Hz) ]. ( accel bias diffusion )
- gyroscope_noise_density: 0.00005770 # [ rad / s / sqrt(Hz) ] ( gyro "white noise" )
- gyroscope_random_walk: 0.00001565 # [ rad / s^2 / sqrt(Hz) ] ( gyro bias diffusion )
-# accelerometer_noise_density: 0.07 # [ m / s^2 / sqrt(Hz) ] ( accel "white noise" )
-# accelerometer_random_walk: 0.009 # [ m / s^3 / sqrt(Hz) ]. ( accel bias diffusion )
-# gyroscope_noise_density: 0.001 # [ rad / s / sqrt(Hz) ] ( gyro "white noise" )
-# gyroscope_random_walk: 0.0003 # [ rad / s^2 / sqrt(Hz) ] ( gyro bias diffusion )
+# accelerometer_noise_density: 0.00333388 # [ m / s^2 / sqrt(Hz) ] ( accel "white noise" )
+# accelerometer_random_walk: 0.00047402 # [ m / s^3 / sqrt(Hz) ]. ( accel bias diffusion )
+# gyroscope_noise_density: 0.00005770 # [ rad / s / sqrt(Hz) ] ( gyro "white noise" )
+# gyroscope_random_walk: 0.00001565 # [ rad / s^2 / sqrt(Hz) ] ( gyro bias diffusion )
+ accelerometer_noise_density: 0.07 # [ m / s^2 / sqrt(Hz) ] ( accel "white noise" )
+ accelerometer_random_walk: 0.009 # [ m / s^3 / sqrt(Hz) ]. ( accel bias diffusion )
+ gyroscope_noise_density: 0.001 # [ rad / s / sqrt(Hz) ] ( gyro "white noise" )
+ gyroscope_random_walk: 0.0003 # [ rad / s^2 / sqrt(Hz) ] ( gyro bias diffusion )
model: calibrated
rostopic: /mavros/imu/data
time_offset: 0.0
diff --git a/config/rpng_aruco/estimator_config.yaml b/config/rpng_aruco/estimator_config.yaml
index 1bb04aba1..24a4186f6 100644
--- a/config/rpng_aruco/estimator_config.yaml
+++ b/config/rpng_aruco/estimator_config.yaml
@@ -27,7 +27,7 @@ feat_rep_aruco: "ANCHORED_MSCKF_INVERSE_DEPTH"
# zero velocity update parameters we can use
# we support either IMU-based or disparity detection.
-try_zupt: true
+try_zupt: false
zupt_chi2_multipler: 0 # set to 0 for only disp-based
zupt_max_velocity: 0.1
zupt_noise_multiplier: 50
@@ -39,9 +39,10 @@ zupt_only_at_beginning: true
init_window_time: 2.0
init_imu_thresh: 1.2
-init_max_disparity: 1.5
-init_max_features: 25
+init_max_disparity: 2.0
+init_max_features: 75
+init_dyn_use: false
init_dyn_mle_opt_calib: false
init_dyn_mle_max_iter: 50
init_dyn_mle_max_time: 0.05
@@ -77,9 +78,9 @@ filepath_gt: "/tmp/ov_groundtruth.txt"
use_klt: true
num_pts: 150
fast_threshold: 30
-grid_x: 20
-grid_y: 20
-min_px_dist: 30
+grid_x: 5
+grid_y: 5
+min_px_dist: 20
knn_ratio: 0.85
track_frequency: 21.0
downsample_cameras: false
@@ -96,12 +97,12 @@ downsize_aruco: true
# ==================================================================
# camera noises and chi-squared threshold multipliers
-up_msckf_sigma_px: 1
+up_msckf_sigma_px: 1.5
up_msckf_chi2_multipler: 1
-up_slam_sigma_px: 1
+up_slam_sigma_px: 1.5
up_slam_chi2_multipler: 1
-up_aruco_sigma_px: 1
-up_aruco_chi2_multipler: 8
+up_aruco_sigma_px: 2.0
+up_aruco_chi2_multipler: 10
# masks for our images
use_mask: false
diff --git a/config/rpng_ironsides/estimator_config.yaml b/config/rpng_ironsides/estimator_config.yaml
index dc04c6177..a0c4ee8d3 100644
--- a/config/rpng_ironsides/estimator_config.yaml
+++ b/config/rpng_ironsides/estimator_config.yaml
@@ -9,11 +9,11 @@ use_rk4int: true # if rk4 integration should be used (overrides imu averaging)
use_stereo: true # if we have more than 1 camera, if we should try to track stereo constraints
max_cameras: 2
-calib_cam_extrinsics: false
+calib_cam_extrinsics: false # degenerate motion
calib_cam_intrinsics: true
-calib_cam_timeoffset: true
+calib_cam_timeoffset: false # degenerate motion
-max_clones: 12
+max_clones: 11
max_slam: 50
max_slam_in_update: 25
max_msckf_in_update: 50
@@ -37,12 +37,12 @@ zupt_only_at_beginning: false
# ==================================================================
# ==================================================================
-#init_window_time: 1.0
init_window_time: 2.0
init_imu_thresh: 0.5
init_max_disparity: 1.5
-init_max_features: 25
+init_max_features: 75
+init_dyn_use: false
init_dyn_mle_opt_calib: false
init_dyn_mle_max_iter: 50
init_dyn_mle_max_time: 0.05
@@ -76,21 +76,23 @@ filepath_gt: "/tmp/ov_groundtruth.txt"
# our front-end feature tracking parameters
# we have a KLT and descriptor based (KLT is better implemented...)
use_klt: true
-num_pts: 200
-fast_threshold: 20
-grid_x: 20
-grid_y: 20
-min_px_dist: 20
+num_pts: 100
+fast_threshold: 30
+grid_x: 5
+grid_y: 5
+min_px_dist: 15
knn_ratio: 0.65
-track_frequency: 21.0
+track_frequency: 31.0
downsample_cameras: false
multi_threading: true
-histogram_method: "CLAHE" # NONE, HISTOGRAM, CLAHE
-
-fi_max_dist: 150
-fi_max_baseline: 200
-fi_max_cond_number: 20000
-fi_triangulate_1d: false
+histogram_method: "HISTOGRAM" # NONE, HISTOGRAM, CLAHE
+
+fi_min_dist: 1.0
+fi_max_dist: 500.0
+#fi_max_baseline: 200
+#fi_max_cond_number: 20000
+#fi_max_runs: 10
+#fi_triangulate_1d: false
# aruco tag tracker for the system
# DICT_6X6_1000 from https://chev.me/arucogen/
diff --git a/config/rpng_sim/estimator_config.yaml b/config/rpng_sim/estimator_config.yaml
index 138cf36d7..2ce974f1b 100644
--- a/config/rpng_sim/estimator_config.yaml
+++ b/config/rpng_sim/estimator_config.yaml
@@ -42,6 +42,7 @@ init_imu_thresh: 1.0
init_max_disparity: 1.5
init_max_features: 15
+init_dyn_use: true # if dynamic initialization should be used
init_dyn_mle_opt_calib: false # if we should optimize calibration during intialization (not recommended)
init_dyn_mle_max_iter: 50 # how many iterations the MLE refinement should use (zero to skip the MLE)
init_dyn_mle_max_time: 0.5 # how many seconds the MLE should be completed in
diff --git a/config/tum_vi/estimator_config.yaml b/config/tum_vi/estimator_config.yaml
index def1cb94c..ce5097280 100644
--- a/config/tum_vi/estimator_config.yaml
+++ b/config/tum_vi/estimator_config.yaml
@@ -37,17 +37,18 @@ zupt_only_at_beginning: true
# ==================================================================
# ==================================================================
-init_window_time: 2.0
+init_window_time: 1.5 # make 2sec if using dynamic...
init_imu_thresh: 0.45 # room1-5:0.45, room6:0.25
-init_max_disparity: 5.0
-init_max_features: 25
+init_max_disparity: 15.0
+init_max_features: 75
+init_dyn_use: false
init_dyn_mle_opt_calib: false
init_dyn_mle_max_iter: 50
init_dyn_mle_max_time: 0.05
init_dyn_mle_max_threads: 6
init_dyn_num_pose: 6
-init_dyn_min_deg: 20.0
+init_dyn_min_deg: 10.0
init_dyn_inflation_ori: 10
init_dyn_inflation_vel: 100
@@ -75,10 +76,10 @@ filepath_gt: "/tmp/ov_groundtruth.txt"
# our front-end feature tracking parameters
# we have a KLT and descriptor based (KLT is better implemented...)
use_klt: true
-num_pts: 150
+num_pts: 200
fast_threshold: 20
-grid_x: 20
-grid_y: 20
+grid_x: 5
+grid_y: 5
min_px_dist: 15
knn_ratio: 0.65
track_frequency: 21.0
diff --git a/config/uzhfpv_indoor/estimator_config.yaml b/config/uzhfpv_indoor/estimator_config.yaml
index e423ffb07..9f9d24d1b 100644
--- a/config/uzhfpv_indoor/estimator_config.yaml
+++ b/config/uzhfpv_indoor/estimator_config.yaml
@@ -16,10 +16,10 @@ calib_cam_timeoffset: true
max_clones: 11
max_slam: 50
max_slam_in_update: 25
-max_msckf_in_update: 30
+max_msckf_in_update: 40
dt_slam_delay: 2
-gravity_mag: 9.81
+gravity_mag: 9.8065 # kalibr calibration
feat_rep_msckf: "GLOBAL_3D"
feat_rep_slam: "ANCHORED_MSCKF_INVERSE_DEPTH"
@@ -27,7 +27,7 @@ feat_rep_aruco: "ANCHORED_MSCKF_INVERSE_DEPTH"
# zero velocity update parameters we can use
# we support either IMU-based or disparity detection.
-try_zupt: true
+try_zupt: false
zupt_chi2_multipler: 0 # set to 0 for only disp-based
zupt_max_velocity: 0.5
zupt_noise_multiplier: 20
@@ -38,10 +38,11 @@ zupt_only_at_beginning: false
# ==================================================================
init_window_time: 2.0
-init_imu_thresh: 0.0
-init_max_disparity: 1.5
-init_max_features: 25
+init_imu_thresh: 0.30
+init_max_disparity: 2.0
+init_max_features: 75
+init_dyn_use: false
init_dyn_mle_opt_calib: false
init_dyn_mle_max_iter: 50
init_dyn_mle_max_time: 0.05
@@ -75,10 +76,10 @@ filepath_gt: "/tmp/ov_groundtruth.txt"
# our front-end feature tracking parameters
# we have a KLT and descriptor based (KLT is better implemented...)
use_klt: true
-num_pts: 100
-fast_threshold: 65
-grid_x: 20
-grid_y: 20
+num_pts: 200
+fast_threshold: 50
+grid_x: 5
+grid_y: 5
min_px_dist: 15
knn_ratio: 0.70
track_frequency: 31.0
@@ -96,9 +97,9 @@ downsize_aruco: true
# ==================================================================
# camera noises and chi-squared threshold multipliers
-up_msckf_sigma_px: 2
+up_msckf_sigma_px: 1.5
up_msckf_chi2_multipler: 1
-up_slam_sigma_px: 2
+up_slam_sigma_px: 1.5
up_slam_chi2_multipler: 1
up_aruco_sigma_px: 1
up_aruco_chi2_multipler: 1
diff --git a/config/uzhfpv_indoor_45/estimator_config.yaml b/config/uzhfpv_indoor_45/estimator_config.yaml
index e423ffb07..9f9d24d1b 100644
--- a/config/uzhfpv_indoor_45/estimator_config.yaml
+++ b/config/uzhfpv_indoor_45/estimator_config.yaml
@@ -16,10 +16,10 @@ calib_cam_timeoffset: true
max_clones: 11
max_slam: 50
max_slam_in_update: 25
-max_msckf_in_update: 30
+max_msckf_in_update: 40
dt_slam_delay: 2
-gravity_mag: 9.81
+gravity_mag: 9.8065 # kalibr calibration
feat_rep_msckf: "GLOBAL_3D"
feat_rep_slam: "ANCHORED_MSCKF_INVERSE_DEPTH"
@@ -27,7 +27,7 @@ feat_rep_aruco: "ANCHORED_MSCKF_INVERSE_DEPTH"
# zero velocity update parameters we can use
# we support either IMU-based or disparity detection.
-try_zupt: true
+try_zupt: false
zupt_chi2_multipler: 0 # set to 0 for only disp-based
zupt_max_velocity: 0.5
zupt_noise_multiplier: 20
@@ -38,10 +38,11 @@ zupt_only_at_beginning: false
# ==================================================================
init_window_time: 2.0
-init_imu_thresh: 0.0
-init_max_disparity: 1.5
-init_max_features: 25
+init_imu_thresh: 0.30
+init_max_disparity: 2.0
+init_max_features: 75
+init_dyn_use: false
init_dyn_mle_opt_calib: false
init_dyn_mle_max_iter: 50
init_dyn_mle_max_time: 0.05
@@ -75,10 +76,10 @@ filepath_gt: "/tmp/ov_groundtruth.txt"
# our front-end feature tracking parameters
# we have a KLT and descriptor based (KLT is better implemented...)
use_klt: true
-num_pts: 100
-fast_threshold: 65
-grid_x: 20
-grid_y: 20
+num_pts: 200
+fast_threshold: 50
+grid_x: 5
+grid_y: 5
min_px_dist: 15
knn_ratio: 0.70
track_frequency: 31.0
@@ -96,9 +97,9 @@ downsize_aruco: true
# ==================================================================
# camera noises and chi-squared threshold multipliers
-up_msckf_sigma_px: 2
+up_msckf_sigma_px: 1.5
up_msckf_chi2_multipler: 1
-up_slam_sigma_px: 2
+up_slam_sigma_px: 1.5
up_slam_chi2_multipler: 1
up_aruco_sigma_px: 1
up_aruco_chi2_multipler: 1
diff --git a/config/uzhfpv_outdoor/estimator_config.yaml b/config/uzhfpv_outdoor/estimator_config.yaml
index fe985638a..e34ceaf8a 100644
--- a/config/uzhfpv_outdoor/estimator_config.yaml
+++ b/config/uzhfpv_outdoor/estimator_config.yaml
@@ -16,10 +16,10 @@ calib_cam_timeoffset: true
max_clones: 11
max_slam: 50
max_slam_in_update: 25
-max_msckf_in_update: 30
+max_msckf_in_update: 40
dt_slam_delay: 2
-gravity_mag: 9.81
+gravity_mag: 9.8065 # kalibr calibration
feat_rep_msckf: "GLOBAL_3D"
feat_rep_slam: "ANCHORED_MSCKF_INVERSE_DEPTH"
@@ -27,7 +27,7 @@ feat_rep_aruco: "ANCHORED_MSCKF_INVERSE_DEPTH"
# zero velocity update parameters we can use
# we support either IMU-based or disparity detection.
-try_zupt: true
+try_zupt: false
zupt_chi2_multipler: 0 # set to 0 for only disp-based
zupt_max_velocity: 0.5
zupt_noise_multiplier: 20
@@ -38,10 +38,11 @@ zupt_only_at_beginning: false
# ==================================================================
init_window_time: 2.0
-init_imu_thresh: 0.0
-init_max_disparity: 1.5
-init_max_features: 25
+init_imu_thresh: 0.30
+init_max_disparity: 2.0
+init_max_features: 75
+init_dyn_use: false
init_dyn_mle_opt_calib: false
init_dyn_mle_max_iter: 50
init_dyn_mle_max_time: 0.05
@@ -75,10 +76,10 @@ filepath_gt: "/tmp/ov_groundtruth.txt"
# our front-end feature tracking parameters
# we have a KLT and descriptor based (KLT is better implemented...)
use_klt: true
-num_pts: 100
-fast_threshold: 65
-grid_x: 20
-grid_y: 20
+num_pts: 200
+fast_threshold: 25
+grid_x: 5
+grid_y: 5
min_px_dist: 15
knn_ratio: 0.70
track_frequency: 31.0
@@ -96,9 +97,9 @@ downsize_aruco: true
# ==================================================================
# camera noises and chi-squared threshold multipliers
-up_msckf_sigma_px: 2
+up_msckf_sigma_px: 1.5
up_msckf_chi2_multipler: 1
-up_slam_sigma_px: 2
+up_slam_sigma_px: 1.5
up_slam_chi2_multipler: 1
up_aruco_sigma_px: 1
up_aruco_chi2_multipler: 1
diff --git a/config/uzhfpv_outdoor_45/estimator_config.yaml b/config/uzhfpv_outdoor_45/estimator_config.yaml
index e423ffb07..9f9d24d1b 100644
--- a/config/uzhfpv_outdoor_45/estimator_config.yaml
+++ b/config/uzhfpv_outdoor_45/estimator_config.yaml
@@ -16,10 +16,10 @@ calib_cam_timeoffset: true
max_clones: 11
max_slam: 50
max_slam_in_update: 25
-max_msckf_in_update: 30
+max_msckf_in_update: 40
dt_slam_delay: 2
-gravity_mag: 9.81
+gravity_mag: 9.8065 # kalibr calibration
feat_rep_msckf: "GLOBAL_3D"
feat_rep_slam: "ANCHORED_MSCKF_INVERSE_DEPTH"
@@ -27,7 +27,7 @@ feat_rep_aruco: "ANCHORED_MSCKF_INVERSE_DEPTH"
# zero velocity update parameters we can use
# we support either IMU-based or disparity detection.
-try_zupt: true
+try_zupt: false
zupt_chi2_multipler: 0 # set to 0 for only disp-based
zupt_max_velocity: 0.5
zupt_noise_multiplier: 20
@@ -38,10 +38,11 @@ zupt_only_at_beginning: false
# ==================================================================
init_window_time: 2.0
-init_imu_thresh: 0.0
-init_max_disparity: 1.5
-init_max_features: 25
+init_imu_thresh: 0.30
+init_max_disparity: 2.0
+init_max_features: 75
+init_dyn_use: false
init_dyn_mle_opt_calib: false
init_dyn_mle_max_iter: 50
init_dyn_mle_max_time: 0.05
@@ -75,10 +76,10 @@ filepath_gt: "/tmp/ov_groundtruth.txt"
# our front-end feature tracking parameters
# we have a KLT and descriptor based (KLT is better implemented...)
use_klt: true
-num_pts: 100
-fast_threshold: 65
-grid_x: 20
-grid_y: 20
+num_pts: 200
+fast_threshold: 50
+grid_x: 5
+grid_y: 5
min_px_dist: 15
knn_ratio: 0.70
track_frequency: 31.0
@@ -96,9 +97,9 @@ downsize_aruco: true
# ==================================================================
# camera noises and chi-squared threshold multipliers
-up_msckf_sigma_px: 2
+up_msckf_sigma_px: 1.5
up_msckf_chi2_multipler: 1
-up_slam_sigma_px: 2
+up_slam_sigma_px: 1.5
up_slam_chi2_multipler: 1
up_aruco_sigma_px: 1
up_aruco_chi2_multipler: 1
diff --git a/docs/gs-datasets.dox b/docs/gs-datasets.dox
index 3284ec7de..05bc0634d 100644
--- a/docs/gs-datasets.dox
+++ b/docs/gs-datasets.dox
@@ -204,6 +204,7 @@ Typically we process the datasets at 1.5x rate so we get a ~20 Hz image feed and
| Dataset Name | Length (km) | Dataset Link | Groundtruth Traj. | Example Launch |
|-------------:|--------|--------------|------------------|------------------|
| Urban 28 | 11.47 | [download](https://sites.google.com/view/complex-urban-dataset/download-lidar-stereo?authuser=0) | [link](https://github.com/rpng/open_vins/tree/master/ov_data/kaist) | [config](https://github.com/rpng/open_vins/blob/master/config/kaist) |
+| Urban 32 | 7.30 | [download](https://sites.google.com/view/complex-urban-dataset/download-lidar-stereo?authuser=0) | [link](https://github.com/rpng/open_vins/tree/master/ov_data/kaist) | [config](https://github.com/rpng/open_vins/blob/master/config/kaist) |
| Urban 38 | 11.42 | [download](https://sites.google.com/view/complex-urban-dataset/download-lidar-stereo?authuser=0) | [link](https://github.com/rpng/open_vins/tree/master/ov_data/kaist) | [config](https://github.com/rpng/open_vins/blob/master/config/kaist) |
| Urban 39 | 11.06 | [download](https://sites.google.com/view/complex-urban-dataset/download-lidar-stereo?authuser=0) | [link](https://github.com/rpng/open_vins/tree/master/ov_data/kaist) | [config](https://github.com/rpng/open_vins/blob/master/config/kaist) |
@m_enddiv
@@ -211,7 +212,7 @@ Typically we process the datasets at 1.5x rate so we get a ~20 Hz image feed and
@section gs-data-kaist-vio KAIST VIO Dataset
-The [KAIST VIO dataset](https://github.com/zinuok/kaistviodataset) @cite Jeon2021RAL is a dataset of a MAV in an indoor 3.15 x 3.60 x 2.50 meter environment which undergoes various trajectory motions.
+The [KAIST VIO dataset](https://github.com/url-kaist/kaistviodataset) @cite Jeon2021RAL is a dataset of a MAV in an indoor 3.15 x 3.60 x 2.50 meter environment which undergoes various trajectory motions.
The camera is intel realsense D435i 25 Hz, while the IMU is 100 Hz sensing rate from the pixelhawk 4 unit.
A groundtruth "baseline" trajectory is also provided from a OptiTrack Mocap system at 50 Hz, the bag files have the marker body frame to IMU frame already applied.
This topic has been provided in ov_data for convinces sake.
diff --git a/docs/gs-installing.dox b/docs/gs-installing.dox
index e6bcd1031..bc7bfe427 100644
--- a/docs/gs-installing.dox
+++ b/docs/gs-installing.dox
@@ -145,6 +145,16 @@ cmake ..
make -j4
@endcode
+If you wish to debug and run with assert statements, you can configure your workspace as follows:
+
+
+@code{.shell-session}
+catkin config --cmake-args -DCMAKE_BUILD_TYPE=Debug
+catkin build --cmake-args -DCMAKE_BUILD_TYPE=Debug
+colcon build --cmake-args -DCMAKE_BUILD_TYPE=Debug
+@endcode
+
+
@section gs-install-oveval Additional Evaluation Requirements
@@ -192,7 +202,7 @@ You can disable this with `catkin build -DENABLE_ARUCO_TAGS=OFF` or `cmake -DENA
-@subsection gs-install-ceres Ceres Solver
+@section gs-install-ceres Ceres Solver
Ceres solver @cite ceres-solver is required for dynamic initialization and backend optimization.
Please refer to their [documentation](http://ceres-solver.org/installation.html#linux) for specifics to your platform.
@@ -205,6 +215,8 @@ To install we can perform the following:
@par Ceres Source Installation
Try to first build with your system with `sudo apt-get install libceres-dev`.
Only fall back onto this if it does not allow you to compile, or want a newer version!
+ You will need to build from source if there is an Eigen miss-match:
+ "Failed to find Ceres - Found Eigen dependency, but the version of Eigen found (3.3.7) does not exactly match the version of Eigen Ceres was compiled with (3.3.4)."
@code{.shell-session}
sudo apt-get install -y cmake libgoogle-glog-dev libgflags-dev libatlas-base-dev libeigen3-dev libsuitesparse-dev
diff --git a/ov_core/CMakeLists.txt b/ov_core/CMakeLists.txt
index a4a1bc23d..2308d6d01 100644
--- a/ov_core/CMakeLists.txt
+++ b/ov_core/CMakeLists.txt
@@ -40,6 +40,11 @@ set(CMAKE_CXX_STANDARD 14)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
set(CMAKE_CXX_EXTENSIONS OFF)
+# Default build as release
+if (NOT CMAKE_BUILD_TYPE)
+ set(CMAKE_BUILD_TYPE "Release" CACHE STRING "Choose the type of build, options are: Debug Release RelWithDebInfo MinSizeRel." FORCE)
+endif ()
+
# Enable compile optimizations
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 -fsee -fomit-frame-pointer -fno-signed-zeros -fno-math-errno -funroll-loops")
diff --git a/ov_core/cmake/ROS1.cmake b/ov_core/cmake/ROS1.cmake
index ab2d319e4..181dc3e7d 100644
--- a/ov_core/cmake/ROS1.cmake
+++ b/ov_core/cmake/ROS1.cmake
@@ -34,7 +34,7 @@ list(APPEND thirdparty_libraries
${Boost_LIBRARIES}
${OpenCV_LIBRARIES}
${catkin_LIBRARIES}
-)
+ )
##################################################
# Make the core library
@@ -42,6 +42,8 @@ list(APPEND thirdparty_libraries
list(APPEND LIBRARY_SOURCES
src/dummy.cpp
+ src/cpi/CpiV1.cpp
+ src/cpi/CpiV2.cpp
src/sim/BsplineSE3.cpp
src/track/TrackBase.cpp
src/track/TrackAruco.cpp
@@ -50,9 +52,10 @@ list(APPEND LIBRARY_SOURCES
src/track/TrackSIM.cpp
src/types/Landmark.cpp
src/feat/Feature.cpp
+ src/feat/FeatureDatabase.cpp
src/feat/FeatureInitializer.cpp
src/utils/print.cpp
-)
+ )
file(GLOB_RECURSE LIBRARY_HEADERS "src/*.h")
add_library(ov_core_lib SHARED ${LIBRARY_SOURCES} ${LIBRARY_HEADERS})
target_link_libraries(ov_core_lib ${thirdparty_libraries})
@@ -61,11 +64,11 @@ install(TARGETS ov_core_lib
ARCHIVE DESTINATION ${CATKIN_PACKAGE_LIB_DESTINATION}
LIBRARY DESTINATION ${CATKIN_PACKAGE_LIB_DESTINATION}
RUNTIME DESTINATION ${CATKIN_PACKAGE_BIN_DESTINATION}
-)
+ )
install(DIRECTORY src/
DESTINATION ${CATKIN_GLOBAL_INCLUDE_DESTINATION}
FILES_MATCHING PATTERN "*.h" PATTERN "*.hpp"
-)
+ )
##################################################
# Make binary files!
@@ -79,7 +82,7 @@ if (catkin_FOUND AND ENABLE_ROS)
ARCHIVE DESTINATION ${CATKIN_PACKAGE_LIB_DESTINATION}
LIBRARY DESTINATION ${CATKIN_PACKAGE_LIB_DESTINATION}
RUNTIME DESTINATION ${CATKIN_PACKAGE_BIN_DESTINATION}
- )
+ )
endif ()
@@ -89,7 +92,7 @@ install(TARGETS test_webcam
ARCHIVE DESTINATION ${CATKIN_PACKAGE_LIB_DESTINATION}
LIBRARY DESTINATION ${CATKIN_PACKAGE_LIB_DESTINATION}
RUNTIME DESTINATION ${CATKIN_PACKAGE_BIN_DESTINATION}
-)
+ )
diff --git a/ov_core/cmake/ROS2.cmake b/ov_core/cmake/ROS2.cmake
index 3c8dc116b..43d7f82bd 100644
--- a/ov_core/cmake/ROS2.cmake
+++ b/ov_core/cmake/ROS2.cmake
@@ -23,7 +23,7 @@ include_directories(
list(APPEND thirdparty_libraries
${Boost_LIBRARIES}
${OpenCV_LIBRARIES}
-)
+ )
##################################################
# Make the core library
@@ -31,6 +31,8 @@ list(APPEND thirdparty_libraries
list(APPEND LIBRARY_SOURCES
src/dummy.cpp
+ src/cpi/CpiV1.cpp
+ src/cpi/CpiV2.cpp
src/sim/BsplineSE3.cpp
src/track/TrackBase.cpp
src/track/TrackAruco.cpp
@@ -39,9 +41,10 @@ list(APPEND LIBRARY_SOURCES
src/track/TrackSIM.cpp
src/types/Landmark.cpp
src/feat/Feature.cpp
+ src/feat/FeatureDatabase.cpp
src/feat/FeatureInitializer.cpp
src/utils/print.cpp
-)
+ )
file(GLOB_RECURSE LIBRARY_HEADERS "src/*.h")
add_library(ov_core_lib SHARED ${LIBRARY_SOURCES} ${LIBRARY_HEADERS})
ament_target_dependencies(ov_core_lib rclcpp cv_bridge)
@@ -51,11 +54,11 @@ install(TARGETS ov_core_lib
LIBRARY DESTINATION lib
RUNTIME DESTINATION bin
PUBLIC_HEADER DESTINATION include
-)
+ )
install(DIRECTORY src/
DESTINATION include
FILES_MATCHING PATTERN "*.h" PATTERN "*.hpp"
-)
+ )
ament_export_include_directories(include)
ament_export_libraries(ov_core_lib)
diff --git a/ov_core/src/cam/CamBase.h b/ov_core/src/cam/CamBase.h
index 6bdcd3244..c8de21798 100644
--- a/ov_core/src/cam/CamBase.h
+++ b/ov_core/src/cam/CamBase.h
@@ -49,6 +49,8 @@ class CamBase {
*/
CamBase(int width, int height) : _width(width), _height(height) {}
+ virtual ~CamBase() {}
+
/**
* @brief This will set and update the camera calibration values.
* This should be called on startup for each camera and after update!
diff --git a/ov_core/src/cam/CamEqui.h b/ov_core/src/cam/CamEqui.h
index 1e4eee1af..0709d0f73 100644
--- a/ov_core/src/cam/CamEqui.h
+++ b/ov_core/src/cam/CamEqui.h
@@ -98,6 +98,8 @@ class CamEqui : public CamBase {
*/
CamEqui(int width, int height) : CamBase(width, height) {}
+ ~CamEqui() {}
+
/**
* @brief Given a raw uv point, this will undistort it based on the camera matrices into normalized camera coords.
* @param uv_dist Raw uv coordinate we wish to undistort
diff --git a/ov_core/src/cam/CamRadtan.h b/ov_core/src/cam/CamRadtan.h
index 021a90728..4f8be1564 100644
--- a/ov_core/src/cam/CamRadtan.h
+++ b/ov_core/src/cam/CamRadtan.h
@@ -89,6 +89,8 @@ class CamRadtan : public CamBase {
*/
CamRadtan(int width, int height) : CamBase(width, height) {}
+ ~CamRadtan() {}
+
/**
* @brief Given a raw uv point, this will undistort it based on the camera matrices into normalized camera coords.
* @param uv_dist Raw uv coordinate we wish to undistort
diff --git a/ov_core/src/cpi/CpiBase.h b/ov_core/src/cpi/CpiBase.h
index a18ffab5a..a05be84f5 100644
--- a/ov_core/src/cpi/CpiBase.h
+++ b/ov_core/src/cpi/CpiBase.h
@@ -73,6 +73,8 @@ class CpiBase {
e_3x = skew_x(e_3);
}
+ virtual ~CpiBase() {}
+
/**
* @brief Set linearization points of the integration.
* @param[in] b_w_lin_ gyroscope bias linearization point
diff --git a/ov_core/src/cpi/CpiV1.cpp b/ov_core/src/cpi/CpiV1.cpp
new file mode 100644
index 000000000..7b9224d0b
--- /dev/null
+++ b/ov_core/src/cpi/CpiV1.cpp
@@ -0,0 +1,315 @@
+
+/*
+ * MIT License
+ * Copyright (c) 2018 Kevin Eckenhoff
+ * Copyright (c) 2018 Patrick Geneva
+ * Copyright (c) 2018 Guoquan Huang
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in all
+ * copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+ * SOFTWARE.
+ */
+
+#include "CpiV1.h"
+
+#include "utils/quat_ops.h"
+
+using namespace ov_core;
+
+void CpiV1::feed_IMU(double t_0, double t_1, Eigen::Matrix w_m_0, Eigen::Matrix a_m_0,
+ Eigen::Matrix w_m_1, Eigen::Matrix a_m_1) {
+
+ // Get time difference
+ double delta_t = t_1 - t_0;
+ DT += delta_t;
+
+ // If no time has passed do nothing
+ if (delta_t == 0) {
+ return;
+ }
+
+ // Get estimated imu readings
+ Eigen::Matrix w_hat = w_m_0 - b_w_lin;
+ Eigen::Matrix a_hat = a_m_0 - b_a_lin;
+
+ // If averaging, average
+ if (imu_avg) {
+ w_hat += w_m_1 - b_w_lin;
+ w_hat = 0.5 * w_hat;
+ a_hat += a_m_1 - b_a_lin;
+ a_hat = .5 * a_hat;
+ }
+
+ // Get angle change w*dt
+ Eigen::Matrix w_hatdt = w_hat * delta_t;
+
+ // Get entries of w_hat
+ double w_1 = w_hat(0, 0);
+ double w_2 = w_hat(1, 0);
+ double w_3 = w_hat(2, 0);
+
+ // Get magnitude of w and wdt
+ double mag_w = w_hat.norm();
+ double w_dt = mag_w * delta_t;
+
+ // Threshold to determine if equations will be unstable
+ bool small_w = (mag_w < 0.008726646);
+
+ // Get some of the variables used in the preintegration equations
+ double dt_2 = pow(delta_t, 2);
+ double cos_wt = cos(w_dt);
+ double sin_wt = sin(w_dt);
+
+ Eigen::Matrix w_x = skew_x(w_hat);
+ Eigen::Matrix a_x = skew_x(a_hat);
+ Eigen::Matrix w_tx = skew_x(w_hatdt);
+ Eigen::Matrix w_x_2 = w_x * w_x;
+
+ //==========================================================================
+ // MEASUREMENT MEANS
+ //==========================================================================
+
+ // Get relative rotation
+ Eigen::Matrix R_tau2tau1 = small_w ? eye3 - delta_t * w_x + (pow(delta_t, 2) / 2) * w_x_2
+ : eye3 - (sin_wt / mag_w) * w_x + ((1.0 - cos_wt) / (pow(mag_w, 2.0))) * w_x_2;
+
+ // Updated rotation and its transpose
+ Eigen::Matrix R_k2tau1 = R_tau2tau1 * R_k2tau;
+ Eigen::Matrix R_tau12k = R_k2tau1.transpose();
+
+ // Intermediate variables for evaluating the measurement/bias Jacobian update
+ double f_1;
+ double f_2;
+ double f_3;
+ double f_4;
+
+ if (small_w) {
+ f_1 = -(pow(delta_t, 3) / 3);
+ f_2 = (pow(delta_t, 4) / 8);
+ f_3 = -(pow(delta_t, 2) / 2);
+ f_4 = (pow(delta_t, 3) / 6);
+ } else {
+ f_1 = (w_dt * cos_wt - sin_wt) / (pow(mag_w, 3));
+ f_2 = (pow(w_dt, 2) - 2 * cos_wt - 2 * w_dt * sin_wt + 2) / (2 * pow(mag_w, 4));
+ f_3 = -(1 - cos_wt) / pow(mag_w, 2);
+ f_4 = (w_dt - sin_wt) / pow(mag_w, 3);
+ }
+
+ // Compute the main part of our analytical means
+ Eigen::Matrix alpha_arg = ((dt_2 / 2.0) * eye3 + f_1 * w_x + f_2 * w_x_2);
+ Eigen::Matrix Beta_arg = (delta_t * eye3 + f_3 * w_x + f_4 * w_x_2);
+
+ // Matrices that will multiply the a_hat in the update expressions
+ Eigen::MatrixXd H_al = R_tau12k * alpha_arg;
+ Eigen::MatrixXd H_be = R_tau12k * Beta_arg;
+
+ // Update the measurement means
+ alpha_tau += beta_tau * delta_t + H_al * a_hat;
+ beta_tau += H_be * a_hat;
+
+ //==========================================================================
+ // BIAS JACOBIANS (ANALYTICAL)
+ //==========================================================================
+
+ // Get right Jacobian
+ Eigen::Matrix J_r_tau1 =
+ small_w ? eye3 - .5 * w_tx + (1.0 / 6.0) * w_tx * w_tx
+ : eye3 - ((1 - cos_wt) / (pow((w_dt), 2.0))) * w_tx + ((w_dt - sin_wt) / (pow(w_dt, 3.0))) * w_tx * w_tx;
+
+ // Update orientation in respect to gyro bias Jacobians
+ J_q = R_tau2tau1 * J_q + J_r_tau1 * delta_t;
+
+ // Update alpha and beta in respect to accel bias Jacobians
+ H_a -= H_al;
+ H_a += delta_t * H_b;
+ H_b -= H_be;
+
+ // Derivatives of R_tau12k wrt bias_w entries
+ Eigen::MatrixXd d_R_bw_1 = -R_tau12k * skew_x(J_q * e_1);
+ Eigen::MatrixXd d_R_bw_2 = -R_tau12k * skew_x(J_q * e_2);
+ Eigen::MatrixXd d_R_bw_3 = -R_tau12k * skew_x(J_q * e_3);
+
+ // Now compute the gyro bias Jacobian terms
+ double df_1_dbw_1;
+ double df_1_dbw_2;
+ double df_1_dbw_3;
+
+ double df_2_dbw_1;
+ double df_2_dbw_2;
+ double df_2_dbw_3;
+
+ double df_3_dbw_1;
+ double df_3_dbw_2;
+ double df_3_dbw_3;
+
+ double df_4_dbw_1;
+ double df_4_dbw_2;
+ double df_4_dbw_3;
+
+ if (small_w) {
+ double df_1_dw_mag = -(pow(delta_t, 5) / 15);
+ df_1_dbw_1 = w_1 * df_1_dw_mag;
+ df_1_dbw_2 = w_2 * df_1_dw_mag;
+ df_1_dbw_3 = w_3 * df_1_dw_mag;
+
+ double df_2_dw_mag = (pow(delta_t, 6) / 72);
+ df_2_dbw_1 = w_1 * df_2_dw_mag;
+ df_2_dbw_2 = w_2 * df_2_dw_mag;
+ df_2_dbw_3 = w_3 * df_2_dw_mag;
+
+ double df_3_dw_mag = -(pow(delta_t, 4) / 12);
+ df_3_dbw_1 = w_1 * df_3_dw_mag;
+ df_3_dbw_2 = w_2 * df_3_dw_mag;
+ df_3_dbw_3 = w_3 * df_3_dw_mag;
+
+ double df_4_dw_mag = (pow(delta_t, 5) / 60);
+ df_4_dbw_1 = w_1 * df_4_dw_mag;
+ df_4_dbw_2 = w_2 * df_4_dw_mag;
+ df_4_dbw_3 = w_3 * df_4_dw_mag;
+ } else {
+ double df_1_dw_mag = (pow(w_dt, 2) * sin_wt - 3 * sin_wt + 3 * w_dt * cos_wt) / pow(mag_w, 5);
+ df_1_dbw_1 = w_1 * df_1_dw_mag;
+ df_1_dbw_2 = w_2 * df_1_dw_mag;
+ df_1_dbw_3 = w_3 * df_1_dw_mag;
+
+ double df_2_dw_mag = (pow(w_dt, 2) - 4 * cos_wt - 4 * w_dt * sin_wt + pow(w_dt, 2) * cos_wt + 4) / (pow(mag_w, 6));
+ df_2_dbw_1 = w_1 * df_2_dw_mag;
+ df_2_dbw_2 = w_2 * df_2_dw_mag;
+ df_2_dbw_3 = w_3 * df_2_dw_mag;
+
+ double df_3_dw_mag = (2 * (cos_wt - 1) + w_dt * sin_wt) / (pow(mag_w, 4));
+ df_3_dbw_1 = w_1 * df_3_dw_mag;
+ df_3_dbw_2 = w_2 * df_3_dw_mag;
+ df_3_dbw_3 = w_3 * df_3_dw_mag;
+
+ double df_4_dw_mag = (2 * w_dt + w_dt * cos_wt - 3 * sin_wt) / (pow(mag_w, 5));
+ df_4_dbw_1 = w_1 * df_4_dw_mag;
+ df_4_dbw_2 = w_2 * df_4_dw_mag;
+ df_4_dbw_3 = w_3 * df_4_dw_mag;
+ }
+
+ // Update alpha and beta gyro bias Jacobians
+ J_a += J_b * delta_t;
+ J_a.block(0, 0, 3, 1) +=
+ (d_R_bw_1 * alpha_arg + R_tau12k * (df_1_dbw_1 * w_x - f_1 * e_1x + df_2_dbw_1 * w_x_2 - f_2 * (e_1x * w_x + w_x * e_1x))) * a_hat;
+ J_a.block(0, 1, 3, 1) +=
+ (d_R_bw_2 * alpha_arg + R_tau12k * (df_1_dbw_2 * w_x - f_1 * e_2x + df_2_dbw_2 * w_x_2 - f_2 * (e_2x * w_x + w_x * e_2x))) * a_hat;
+ J_a.block(0, 2, 3, 1) +=
+ (d_R_bw_3 * alpha_arg + R_tau12k * (df_1_dbw_3 * w_x - f_1 * e_3x + df_2_dbw_3 * w_x_2 - f_2 * (e_3x * w_x + w_x * e_3x))) * a_hat;
+ J_b.block(0, 0, 3, 1) +=
+ (d_R_bw_1 * Beta_arg + R_tau12k * (df_3_dbw_1 * w_x - f_3 * e_1x + df_4_dbw_1 * w_x_2 - f_4 * (e_1x * w_x + w_x * e_1x))) * a_hat;
+ J_b.block(0, 1, 3, 1) +=
+ (d_R_bw_2 * Beta_arg + R_tau12k * (df_3_dbw_2 * w_x - f_3 * e_2x + df_4_dbw_2 * w_x_2 - f_4 * (e_2x * w_x + w_x * e_2x))) * a_hat;
+ J_b.block(0, 2, 3, 1) +=
+ (d_R_bw_3 * Beta_arg + R_tau12k * (df_3_dbw_3 * w_x - f_3 * e_3x + df_4_dbw_3 * w_x_2 - f_4 * (e_3x * w_x + w_x * e_3x))) * a_hat;
+
+ //==========================================================================
+ // MEASUREMENT COVARIANCE
+ //==========================================================================
+
+ // Going to need orientation at intermediate time i.e. at .5*dt;
+ Eigen::Matrix R_mid =
+ small_w ? eye3 - .5 * delta_t * w_x + (pow(.5 * delta_t, 2) / 2) * w_x_2
+ : eye3 - (sin(mag_w * .5 * delta_t) / mag_w) * w_x + ((1.0 - cos(mag_w * .5 * delta_t)) / (pow(mag_w, 2.0))) * w_x_2;
+ R_mid = R_mid * R_k2tau;
+
+ // Compute covariance (in this implementation, we use RK4)
+ // k1-------------------------------------------------------------------------------------------------
+
+ // Build state Jacobian
+ Eigen::Matrix F_k1 = Eigen::Matrix::Zero();
+ F_k1.block(0, 0, 3, 3) = -w_x;
+ F_k1.block(0, 3, 3, 3) = -eye3;
+ F_k1.block(6, 0, 3, 3) = -R_k2tau.transpose() * a_x;
+ F_k1.block(6, 9, 3, 3) = -R_k2tau.transpose();
+ F_k1.block(12, 6, 3, 3) = eye3;
+
+ // Build noise Jacobian
+ Eigen::Matrix G_k1 = Eigen::Matrix::Zero();
+ G_k1.block(0, 0, 3, 3) = -eye3;
+ G_k1.block(3, 3, 3, 3) = eye3;
+ G_k1.block(6, 6, 3, 3) = -R_k2tau.transpose();
+ G_k1.block(9, 9, 3, 3) = eye3;
+
+ // Get covariance derivative
+ Eigen::Matrix P_dot_k1 = F_k1 * P_meas + P_meas * F_k1.transpose() + G_k1 * Q_c * G_k1.transpose();
+
+ // k2-------------------------------------------------------------------------------------------------
+
+ // Build state Jacobian
+ Eigen::Matrix F_k2 = Eigen::Matrix::Zero();
+ F_k2.block(0, 0, 3, 3) = -w_x;
+ F_k2.block(0, 3, 3, 3) = -eye3;
+ F_k2.block(6, 0, 3, 3) = -R_mid.transpose() * a_x;
+ F_k2.block(6, 9, 3, 3) = -R_mid.transpose();
+ F_k2.block(12, 6, 3, 3) = eye3;
+
+ // Build noise Jacobian
+ Eigen::Matrix G_k2 = Eigen::Matrix::Zero();
+ G_k2.block(0, 0, 3, 3) = -eye3;
+ G_k2.block(3, 3, 3, 3) = eye3;
+ G_k2.block(6, 6, 3, 3) = -R_mid.transpose();
+ G_k2.block(9, 9, 3, 3) = eye3;
+
+ // Get covariance derivative
+ Eigen::Matrix P_k2 = P_meas + P_dot_k1 * delta_t / 2.0;
+ Eigen::Matrix P_dot_k2 = F_k2 * P_k2 + P_k2 * F_k2.transpose() + G_k2 * Q_c * G_k2.transpose();
+
+ // k3-------------------------------------------------------------------------------------------------
+
+ // Our state and noise Jacobians are the same as k2
+ // Since k2 and k3 correspond to the same estimates for the midpoint
+ Eigen::Matrix F_k3 = F_k2;
+ Eigen::Matrix G_k3 = G_k2;
+
+ // Get covariance derivative
+ Eigen::Matrix P_k3 = P_meas + P_dot_k2 * delta_t / 2.0;
+ Eigen::Matrix P_dot_k3 = F_k3 * P_k3 + P_k3 * F_k3.transpose() + G_k3 * Q_c * G_k3.transpose();
+
+ // k4-------------------------------------------------------------------------------------------------
+
+ // Build state Jacobian
+ Eigen::Matrix F_k4 = Eigen::Matrix::Zero();
+ F_k4.block(0, 0, 3, 3) = -w_x;
+ F_k4.block(0, 3, 3, 3) = -eye3;
+ F_k4.block(6, 0, 3, 3) = -R_k2tau1.transpose() * a_x;
+ F_k4.block(6, 9, 3, 3) = -R_k2tau1.transpose();
+ F_k4.block(12, 6, 3, 3) = eye3;
+
+ // Build noise Jacobian
+ Eigen::Matrix G_k4 = Eigen::Matrix::Zero();
+ G_k4.block(0, 0, 3, 3) = -eye3;
+ G_k4.block(3, 3, 3, 3) = eye3;
+ G_k4.block(6, 6, 3, 3) = -R_k2tau1.transpose();
+ G_k4.block(9, 9, 3, 3) = eye3;
+
+ // Get covariance derivative
+ Eigen::Matrix P_k4 = P_meas + P_dot_k3 * delta_t;
+ Eigen::Matrix P_dot_k4 = F_k4 * P_k4 + P_k4 * F_k4.transpose() + G_k4 * Q_c * G_k4.transpose();
+
+ // done-------------------------------------------------------------------------------------------------
+
+ // Collect covariance solution
+ // Ensure it is positive definite
+ P_meas += (delta_t / 6.0) * (P_dot_k1 + 2.0 * P_dot_k2 + 2.0 * P_dot_k3 + P_dot_k4);
+ P_meas = 0.5 * (P_meas + P_meas.transpose());
+
+ // Update rotation mean
+ // Note we had to wait to do this, since we use the old orientation in our covariance calculation
+ R_k2tau = R_k2tau1;
+ q_k2tau = rot_2_quat(R_k2tau);
+}
diff --git a/ov_core/src/cpi/CpiV1.h b/ov_core/src/cpi/CpiV1.h
index fdf64b311..bb5dbf0bf 100644
--- a/ov_core/src/cpi/CpiV1.h
+++ b/ov_core/src/cpi/CpiV1.h
@@ -27,7 +27,7 @@
*/
#include "CpiBase.h"
-#include "utils/quat_ops.h"
+
#include
namespace ov_core {
@@ -63,6 +63,8 @@ class CpiV1 : public CpiBase {
CpiV1(double sigma_w, double sigma_wb, double sigma_a, double sigma_ab, bool imu_avg_ = false)
: CpiBase(sigma_w, sigma_wb, sigma_a, sigma_ab, imu_avg_) {}
+ virtual ~CpiV1() {}
+
/**
* @brief Our precompound function for Model 1
* @param[in] t_0 first IMU timestamp
@@ -77,288 +79,7 @@ class CpiV1 : public CpiBase {
*/
void feed_IMU(double t_0, double t_1, Eigen::Matrix w_m_0, Eigen::Matrix a_m_0,
Eigen::Matrix w_m_1 = Eigen::Matrix::Zero(),
- Eigen::Matrix a_m_1 = Eigen::Matrix::Zero()) {
-
- // Get time difference
- double delta_t = t_1 - t_0;
- DT += delta_t;
-
- // If no time has passed do nothing
- if (delta_t == 0) {
- return;
- }
-
- // Get estimated imu readings
- Eigen::Matrix w_hat = w_m_0 - b_w_lin;
- Eigen::Matrix a_hat = a_m_0 - b_a_lin;
-
- // If averaging, average
- if (imu_avg) {
- w_hat += w_m_1 - b_w_lin;
- w_hat = 0.5 * w_hat;
- a_hat += a_m_1 - b_a_lin;
- a_hat = .5 * a_hat;
- }
-
- // Get angle change w*dt
- Eigen::Matrix w_hatdt = w_hat * delta_t;
-
- // Get entries of w_hat
- double w_1 = w_hat(0, 0);
- double w_2 = w_hat(1, 0);
- double w_3 = w_hat(2, 0);
-
- // Get magnitude of w and wdt
- double mag_w = w_hat.norm();
- double w_dt = mag_w * delta_t;
-
- // Threshold to determine if equations will be unstable
- bool small_w = (mag_w < 0.008726646);
-
- // Get some of the variables used in the preintegration equations
- double dt_2 = pow(delta_t, 2);
- double cos_wt = cos(w_dt);
- double sin_wt = sin(w_dt);
-
- Eigen::Matrix w_x = skew_x(w_hat);
- Eigen::Matrix a_x = skew_x(a_hat);
- Eigen::Matrix w_tx = skew_x(w_hatdt);
- Eigen::Matrix w_x_2 = w_x * w_x;
-
- //==========================================================================
- // MEASUREMENT MEANS
- //==========================================================================
-
- // Get relative rotation
- Eigen::Matrix R_tau2tau1 = small_w ? eye3 - delta_t * w_x + (pow(delta_t, 2) / 2) * w_x_2
- : eye3 - (sin_wt / mag_w) * w_x + ((1.0 - cos_wt) / (pow(mag_w, 2.0))) * w_x_2;
-
- // Updated rotation and its transpose
- Eigen::Matrix R_k2tau1 = R_tau2tau1 * R_k2tau;
- Eigen::Matrix R_tau12k = R_k2tau1.transpose();
-
- // Intermediate variables for evaluating the measurement/bias Jacobian update
- double f_1;
- double f_2;
- double f_3;
- double f_4;
-
- if (small_w) {
- f_1 = -(pow(delta_t, 3) / 3);
- f_2 = (pow(delta_t, 4) / 8);
- f_3 = -(pow(delta_t, 2) / 2);
- f_4 = (pow(delta_t, 3) / 6);
- } else {
- f_1 = (w_dt * cos_wt - sin_wt) / (pow(mag_w, 3));
- f_2 = (pow(w_dt, 2) - 2 * cos_wt - 2 * w_dt * sin_wt + 2) / (2 * pow(mag_w, 4));
- f_3 = -(1 - cos_wt) / pow(mag_w, 2);
- f_4 = (w_dt - sin_wt) / pow(mag_w, 3);
- }
-
- // Compute the main part of our analytical means
- Eigen::Matrix alpha_arg = ((dt_2 / 2.0) * eye3 + f_1 * w_x + f_2 * w_x_2);
- Eigen::Matrix Beta_arg = (delta_t * eye3 + f_3 * w_x + f_4 * w_x_2);
-
- // Matrices that will multiply the a_hat in the update expressions
- Eigen::MatrixXd H_al = R_tau12k * alpha_arg;
- Eigen::MatrixXd H_be = R_tau12k * Beta_arg;
-
- // Update the measurement means
- alpha_tau += beta_tau * delta_t + H_al * a_hat;
- beta_tau += H_be * a_hat;
-
- //==========================================================================
- // BIAS JACOBIANS (ANALYTICAL)
- //==========================================================================
-
- // Get right Jacobian
- Eigen::Matrix J_r_tau1 =
- small_w ? eye3 - .5 * w_tx + (1.0 / 6.0) * w_tx * w_tx
- : eye3 - ((1 - cos_wt) / (pow((w_dt), 2.0))) * w_tx + ((w_dt - sin_wt) / (pow(w_dt, 3.0))) * w_tx * w_tx;
-
- // Update orientation in respect to gyro bias Jacobians
- J_q = R_tau2tau1 * J_q + J_r_tau1 * delta_t;
-
- // Update alpha and beta in respect to accel bias Jacobians
- H_a -= H_al;
- H_a += delta_t * H_b;
- H_b -= H_be;
-
- // Derivatives of R_tau12k wrt bias_w entries
- Eigen::MatrixXd d_R_bw_1 = -R_tau12k * skew_x(J_q * e_1);
- Eigen::MatrixXd d_R_bw_2 = -R_tau12k * skew_x(J_q * e_2);
- Eigen::MatrixXd d_R_bw_3 = -R_tau12k * skew_x(J_q * e_3);
-
- // Now compute the gyro bias Jacobian terms
- double df_1_dbw_1;
- double df_1_dbw_2;
- double df_1_dbw_3;
-
- double df_2_dbw_1;
- double df_2_dbw_2;
- double df_2_dbw_3;
-
- double df_3_dbw_1;
- double df_3_dbw_2;
- double df_3_dbw_3;
-
- double df_4_dbw_1;
- double df_4_dbw_2;
- double df_4_dbw_3;
-
- if (small_w) {
- double df_1_dw_mag = -(pow(delta_t, 5) / 15);
- df_1_dbw_1 = w_1 * df_1_dw_mag;
- df_1_dbw_2 = w_2 * df_1_dw_mag;
- df_1_dbw_3 = w_3 * df_1_dw_mag;
-
- double df_2_dw_mag = (pow(delta_t, 6) / 72);
- df_2_dbw_1 = w_1 * df_2_dw_mag;
- df_2_dbw_2 = w_2 * df_2_dw_mag;
- df_2_dbw_3 = w_3 * df_2_dw_mag;
-
- double df_3_dw_mag = -(pow(delta_t, 4) / 12);
- df_3_dbw_1 = w_1 * df_3_dw_mag;
- df_3_dbw_2 = w_2 * df_3_dw_mag;
- df_3_dbw_3 = w_3 * df_3_dw_mag;
-
- double df_4_dw_mag = (pow(delta_t, 5) / 60);
- df_4_dbw_1 = w_1 * df_4_dw_mag;
- df_4_dbw_2 = w_2 * df_4_dw_mag;
- df_4_dbw_3 = w_3 * df_4_dw_mag;
- } else {
- double df_1_dw_mag = (pow(w_dt, 2) * sin_wt - 3 * sin_wt + 3 * w_dt * cos_wt) / pow(mag_w, 5);
- df_1_dbw_1 = w_1 * df_1_dw_mag;
- df_1_dbw_2 = w_2 * df_1_dw_mag;
- df_1_dbw_3 = w_3 * df_1_dw_mag;
-
- double df_2_dw_mag = (pow(w_dt, 2) - 4 * cos_wt - 4 * w_dt * sin_wt + pow(w_dt, 2) * cos_wt + 4) / (pow(mag_w, 6));
- df_2_dbw_1 = w_1 * df_2_dw_mag;
- df_2_dbw_2 = w_2 * df_2_dw_mag;
- df_2_dbw_3 = w_3 * df_2_dw_mag;
-
- double df_3_dw_mag = (2 * (cos_wt - 1) + w_dt * sin_wt) / (pow(mag_w, 4));
- df_3_dbw_1 = w_1 * df_3_dw_mag;
- df_3_dbw_2 = w_2 * df_3_dw_mag;
- df_3_dbw_3 = w_3 * df_3_dw_mag;
-
- double df_4_dw_mag = (2 * w_dt + w_dt * cos_wt - 3 * sin_wt) / (pow(mag_w, 5));
- df_4_dbw_1 = w_1 * df_4_dw_mag;
- df_4_dbw_2 = w_2 * df_4_dw_mag;
- df_4_dbw_3 = w_3 * df_4_dw_mag;
- }
-
- // Update alpha and beta gyro bias Jacobians
- J_a += J_b * delta_t;
- J_a.block(0, 0, 3, 1) +=
- (d_R_bw_1 * alpha_arg + R_tau12k * (df_1_dbw_1 * w_x - f_1 * e_1x + df_2_dbw_1 * w_x_2 - f_2 * (e_1x * w_x + w_x * e_1x))) * a_hat;
- J_a.block(0, 1, 3, 1) +=
- (d_R_bw_2 * alpha_arg + R_tau12k * (df_1_dbw_2 * w_x - f_1 * e_2x + df_2_dbw_2 * w_x_2 - f_2 * (e_2x * w_x + w_x * e_2x))) * a_hat;
- J_a.block(0, 2, 3, 1) +=
- (d_R_bw_3 * alpha_arg + R_tau12k * (df_1_dbw_3 * w_x - f_1 * e_3x + df_2_dbw_3 * w_x_2 - f_2 * (e_3x * w_x + w_x * e_3x))) * a_hat;
- J_b.block(0, 0, 3, 1) +=
- (d_R_bw_1 * Beta_arg + R_tau12k * (df_3_dbw_1 * w_x - f_3 * e_1x + df_4_dbw_1 * w_x_2 - f_4 * (e_1x * w_x + w_x * e_1x))) * a_hat;
- J_b.block(0, 1, 3, 1) +=
- (d_R_bw_2 * Beta_arg + R_tau12k * (df_3_dbw_2 * w_x - f_3 * e_2x + df_4_dbw_2 * w_x_2 - f_4 * (e_2x * w_x + w_x * e_2x))) * a_hat;
- J_b.block(0, 2, 3, 1) +=
- (d_R_bw_3 * Beta_arg + R_tau12k * (df_3_dbw_3 * w_x - f_3 * e_3x + df_4_dbw_3 * w_x_2 - f_4 * (e_3x * w_x + w_x * e_3x))) * a_hat;
-
- //==========================================================================
- // MEASUREMENT COVARIANCE
- //==========================================================================
-
- // Going to need orientation at intermediate time i.e. at .5*dt;
- Eigen::Matrix R_mid =
- small_w ? eye3 - .5 * delta_t * w_x + (pow(.5 * delta_t, 2) / 2) * w_x_2
- : eye3 - (sin(mag_w * .5 * delta_t) / mag_w) * w_x + ((1.0 - cos(mag_w * .5 * delta_t)) / (pow(mag_w, 2.0))) * w_x_2;
- R_mid = R_mid * R_k2tau;
-
- // Compute covariance (in this implementation, we use RK4)
- // k1-------------------------------------------------------------------------------------------------
-
- // Build state Jacobian
- Eigen::Matrix F_k1 = Eigen::Matrix::Zero();
- F_k1.block(0, 0, 3, 3) = -w_x;
- F_k1.block(0, 3, 3, 3) = -eye3;
- F_k1.block(6, 0, 3, 3) = -R_k2tau.transpose() * a_x;
- F_k1.block(6, 9, 3, 3) = -R_k2tau.transpose();
- F_k1.block(12, 6, 3, 3) = eye3;
-
- // Build noise Jacobian
- Eigen::Matrix G_k1 = Eigen::Matrix::Zero();
- G_k1.block(0, 0, 3, 3) = -eye3;
- G_k1.block(3, 3, 3, 3) = eye3;
- G_k1.block(6, 6, 3, 3) = -R_k2tau.transpose();
- G_k1.block(9, 9, 3, 3) = eye3;
-
- // Get covariance derivative
- Eigen::Matrix P_dot_k1 = F_k1 * P_meas + P_meas * F_k1.transpose() + G_k1 * Q_c * G_k1.transpose();
-
- // k2-------------------------------------------------------------------------------------------------
-
- // Build state Jacobian
- Eigen::Matrix F_k2 = Eigen::Matrix::Zero();
- F_k2.block(0, 0, 3, 3) = -w_x;
- F_k2.block(0, 3, 3, 3) = -eye3;
- F_k2.block(6, 0, 3, 3) = -R_mid.transpose() * a_x;
- F_k2.block(6, 9, 3, 3) = -R_mid.transpose();
- F_k2.block(12, 6, 3, 3) = eye3;
-
- // Build noise Jacobian
- Eigen::Matrix G_k2 = Eigen::Matrix::Zero();
- G_k2.block(0, 0, 3, 3) = -eye3;
- G_k2.block(3, 3, 3, 3) = eye3;
- G_k2.block(6, 6, 3, 3) = -R_mid.transpose();
- G_k2.block(9, 9, 3, 3) = eye3;
-
- // Get covariance derivative
- Eigen::Matrix P_k2 = P_meas + P_dot_k1 * delta_t / 2.0;
- Eigen::Matrix P_dot_k2 = F_k2 * P_k2 + P_k2 * F_k2.transpose() + G_k2 * Q_c * G_k2.transpose();
-
- // k3-------------------------------------------------------------------------------------------------
-
- // Our state and noise Jacobians are the same as k2
- // Since k2 and k3 correspond to the same estimates for the midpoint
- Eigen::Matrix F_k3 = F_k2;
- Eigen::Matrix G_k3 = G_k2;
-
- // Get covariance derivative
- Eigen::Matrix P_k3 = P_meas + P_dot_k2 * delta_t / 2.0;
- Eigen::Matrix P_dot_k3 = F_k3 * P_k3 + P_k3 * F_k3.transpose() + G_k3 * Q_c * G_k3.transpose();
-
- // k4-------------------------------------------------------------------------------------------------
-
- // Build state Jacobian
- Eigen::Matrix F_k4 = Eigen::Matrix::Zero();
- F_k4.block(0, 0, 3, 3) = -w_x;
- F_k4.block(0, 3, 3, 3) = -eye3;
- F_k4.block(6, 0, 3, 3) = -R_k2tau1.transpose() * a_x;
- F_k4.block(6, 9, 3, 3) = -R_k2tau1.transpose();
- F_k4.block(12, 6, 3, 3) = eye3;
-
- // Build noise Jacobian
- Eigen::Matrix G_k4 = Eigen::Matrix::Zero();
- G_k4.block(0, 0, 3, 3) = -eye3;
- G_k4.block(3, 3, 3, 3) = eye3;
- G_k4.block(6, 6, 3, 3) = -R_k2tau1.transpose();
- G_k4.block(9, 9, 3, 3) = eye3;
-
- // Get covariance derivative
- Eigen::Matrix P_k4 = P_meas + P_dot_k3 * delta_t;
- Eigen::Matrix P_dot_k4 = F_k4 * P_k4 + P_k4 * F_k4.transpose() + G_k4 * Q_c * G_k4.transpose();
-
- // done-------------------------------------------------------------------------------------------------
-
- // Collect covariance solution
- // Ensure it is positive definite
- P_meas += (delta_t / 6.0) * (P_dot_k1 + 2.0 * P_dot_k2 + 2.0 * P_dot_k3 + P_dot_k4);
- P_meas = 0.5 * (P_meas + P_meas.transpose());
-
- // Update rotation mean
- // Note we had to wait to do this, since we use the old orientation in our covariance calculation
- R_k2tau = R_k2tau1;
- q_k2tau = rot_2_quat(R_k2tau);
- }
+ Eigen::Matrix a_m_1 = Eigen::Matrix::Zero());
};
} // namespace ov_core
diff --git a/ov_core/src/cpi/CpiV2.cpp b/ov_core/src/cpi/CpiV2.cpp
new file mode 100644
index 000000000..10e280631
--- /dev/null
+++ b/ov_core/src/cpi/CpiV2.cpp
@@ -0,0 +1,387 @@
+
+/*
+ * MIT License
+ * Copyright (c) 2018 Kevin Eckenhoff
+ * Copyright (c) 2018 Patrick Geneva
+ * Copyright (c) 2018 Guoquan Huang
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in all
+ * copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+ * SOFTWARE.
+ */
+
+#include "CpiV2.h"
+
+#include "utils/quat_ops.h"
+
+using namespace ov_core;
+
+void CpiV2::feed_IMU(double t_0, double t_1, Eigen::Matrix w_m_0, Eigen::Matrix a_m_0,
+ Eigen::Matrix w_m_1, Eigen::Matrix a_m_1) {
+
+ // Get time difference
+ double delta_t = t_1 - t_0;
+ DT += delta_t;
+
+ // If no time has passed do nothing
+ if (delta_t == 0) {
+ return;
+ }
+
+ // Get estimated imu readings
+ Eigen::Matrix w_hat = w_m_0 - b_w_lin;
+ Eigen::Matrix a_hat = a_m_0 - b_a_lin - R_k2tau * quat_2_Rot(q_k_lin) * grav;
+
+ // If averaging, average
+ // Note: we will average the LOCAL acceleration after getting the relative rotation
+ if (imu_avg) {
+ w_hat += w_m_1 - b_w_lin;
+ w_hat = 0.5 * w_hat;
+ }
+
+ // Get angle change w*dt
+ Eigen::Matrix w_hatdt = w_hat * delta_t;
+
+ // Get entries of w_hat
+ double w_1 = w_hat(0, 0);
+ double w_2 = w_hat(1, 0);
+ double w_3 = w_hat(2, 0);
+
+ // Get magnitude of w and wdt
+ double mag_w = w_hat.norm();
+ double w_dt = mag_w * delta_t;
+
+ // Threshold to determine if equations will be unstable
+ bool small_w = (mag_w < 0.008726646);
+
+ // Get some of the variables used in the preintegration equations
+ double dt_2 = pow(delta_t, 2);
+ double cos_wt = cos(w_dt);
+ double sin_wt = sin(w_dt);
+
+ Eigen::Matrix w_x = skew_x(w_hat);
+ Eigen::Matrix w_tx = skew_x(w_hatdt);
+ Eigen::Matrix w_x_2 = w_x * w_x;
+
+ //==========================================================================
+ // MEASUREMENT MEANS
+ //==========================================================================
+
+ // Get relative rotation
+ Eigen::Matrix R_tau2tau1 = small_w ? eye3 - delta_t * w_x + (pow(delta_t, 2) / 2) * w_x_2
+ : eye3 - (sin_wt / mag_w) * w_x + ((1.0 - cos_wt) / (pow(mag_w, 2.0))) * w_x_2;
+
+ // Updated roation and its transpose
+ Eigen::Matrix R_k2tau1 = R_tau2tau1 * R_k2tau;
+ Eigen::Matrix R_tau12k = R_k2tau1.transpose();
+
+ // If averaging, average the LOCAL acceleration
+ if (imu_avg) {
+ a_hat += a_m_1 - b_a_lin - R_k2tau1 * quat_2_Rot(q_k_lin) * grav;
+ a_hat = 0.5 * a_hat;
+ }
+ Eigen::Matrix a_x = skew_x(a_hat);
+
+ // Intermediate variables for evaluating the measurement/bias Jacobian update
+ double f_1;
+ double f_2;
+ double f_3;
+ double f_4;
+
+ if (small_w) {
+ f_1 = -(pow(delta_t, 3) / 3);
+ f_2 = (pow(delta_t, 4) / 8);
+ f_3 = -(pow(delta_t, 2) / 2);
+ f_4 = (pow(delta_t, 3) / 6);
+ } else {
+ f_1 = (w_dt * cos_wt - sin_wt) / (pow(mag_w, 3));
+ f_2 = (pow(w_dt, 2) - 2 * cos_wt - 2 * w_dt * sin_wt + 2) / (2 * pow(mag_w, 4));
+ f_3 = -(1 - cos_wt) / pow(mag_w, 2);
+ f_4 = (w_dt - sin_wt) / pow(mag_w, 3);
+ }
+
+ // Compute the main part of our analytical means
+ Eigen::Matrix alpha_arg = ((dt_2 / 2.0) * eye3 + f_1 * w_x + f_2 * w_x_2);
+ Eigen::Matrix Beta_arg = (delta_t * eye3 + f_3 * w_x + f_4 * w_x_2);
+
+ // Matrices that will multiply the a_hat in the update expressions
+ Eigen::Matrix H_al = R_tau12k * alpha_arg;
+ Eigen::Matrix H_be = R_tau12k * Beta_arg;
+
+ // Update the measurements
+ alpha_tau += beta_tau * delta_t + H_al * a_hat;
+ beta_tau += H_be * a_hat;
+
+ //==========================================================================
+ // BIAS JACOBIANS (ANALYTICAL)
+ //==========================================================================
+
+ // Get right Jacobian
+ Eigen::Matrix J_r_tau1 =
+ small_w ? eye3 - .5 * w_tx + (1.0 / 6.0) * w_tx * w_tx
+ : eye3 - ((1 - cos_wt) / (pow((w_dt), 2.0))) * w_tx + ((w_dt - sin_wt) / (pow(w_dt, 3.0))) * w_tx * w_tx;
+
+ // Update orientation in respect to gyro bias Jacobians
+ Eigen::Matrix J_save = J_q;
+ J_q = R_tau2tau1 * J_q + J_r_tau1 * delta_t;
+
+ // Update alpha and beta in respect to accel bias Jacobian
+ H_a -= H_al;
+ H_a += delta_t * H_b;
+ H_b -= H_be;
+
+ // Update alpha and beta in respect to q_GtoLIN Jacobian
+ Eigen::Matrix g_k = quat_2_Rot(q_k_lin) * grav;
+ O_a += delta_t * O_b;
+ O_a += -H_al * R_k2tau * skew_x(g_k);
+ O_b += -H_be * R_k2tau * skew_x(g_k);
+
+ // Derivatives of R_tau12k wrt bias_w entries
+ Eigen::MatrixXd d_R_bw_1 = -R_tau12k * skew_x(J_q * e_1);
+ Eigen::MatrixXd d_R_bw_2 = -R_tau12k * skew_x(J_q * e_2);
+ Eigen::MatrixXd d_R_bw_3 = -R_tau12k * skew_x(J_q * e_3);
+
+ // Now compute the gyro bias Jacobian terms
+ double df_1_dbw_1;
+ double df_1_dbw_2;
+ double df_1_dbw_3;
+
+ double df_2_dbw_1;
+ double df_2_dbw_2;
+ double df_2_dbw_3;
+
+ double df_3_dbw_1;
+ double df_3_dbw_2;
+ double df_3_dbw_3;
+
+ double df_4_dbw_1;
+ double df_4_dbw_2;
+ double df_4_dbw_3;
+
+ if (small_w) {
+ double df_1_dw_mag = -(pow(delta_t, 5) / 15);
+ df_1_dbw_1 = w_1 * df_1_dw_mag;
+ df_1_dbw_2 = w_2 * df_1_dw_mag;
+ df_1_dbw_3 = w_3 * df_1_dw_mag;
+
+ double df_2_dw_mag = (pow(delta_t, 6) / 72);
+ df_2_dbw_1 = w_1 * df_2_dw_mag;
+ df_2_dbw_2 = w_2 * df_2_dw_mag;
+ df_2_dbw_3 = w_3 * df_2_dw_mag;
+
+ double df_3_dw_mag = -(pow(delta_t, 4) / 12);
+ df_3_dbw_1 = w_1 * df_3_dw_mag;
+ df_3_dbw_2 = w_2 * df_3_dw_mag;
+ df_3_dbw_3 = w_3 * df_3_dw_mag;
+
+ double df_4_dw_mag = (pow(delta_t, 5) / 60);
+ df_4_dbw_1 = w_1 * df_4_dw_mag;
+ df_4_dbw_2 = w_2 * df_4_dw_mag;
+ df_4_dbw_3 = w_3 * df_4_dw_mag;
+ } else {
+ double df_1_dw_mag = (pow(w_dt, 2) * sin_wt - 3 * sin_wt + 3 * w_dt * cos_wt) / pow(mag_w, 5);
+ df_1_dbw_1 = w_1 * df_1_dw_mag;
+ df_1_dbw_2 = w_2 * df_1_dw_mag;
+ df_1_dbw_3 = w_3 * df_1_dw_mag;
+
+ double df_2_dw_mag = (pow(w_dt, 2) - 4 * cos_wt - 4 * w_dt * sin_wt + pow(w_dt, 2) * cos_wt + 4) / (pow(mag_w, 6));
+ df_2_dbw_1 = w_1 * df_2_dw_mag;
+ df_2_dbw_2 = w_2 * df_2_dw_mag;
+ df_2_dbw_3 = w_3 * df_2_dw_mag;
+
+ double df_3_dw_mag = (2 * (cos_wt - 1) + w_dt * sin_wt) / (pow(mag_w, 4));
+ df_3_dbw_1 = w_1 * df_3_dw_mag;
+ df_3_dbw_2 = w_2 * df_3_dw_mag;
+ df_3_dbw_3 = w_3 * df_3_dw_mag;
+
+ double df_4_dw_mag = (2 * w_dt + w_dt * cos_wt - 3 * sin_wt) / (pow(mag_w, 5));
+ df_4_dbw_1 = w_1 * df_4_dw_mag;
+ df_4_dbw_2 = w_2 * df_4_dw_mag;
+ df_4_dbw_3 = w_3 * df_4_dw_mag;
+ }
+
+ // Gravity rotated into the "tau'th" frame (i.e. start of the measurement interval)
+ Eigen::Matrix g_tau = R_k2tau * quat_2_Rot(q_k_lin) * grav;
+
+ // Update gyro bias Jacobians
+ J_a += J_b * delta_t;
+ J_a.block(0, 0, 3, 1) +=
+ (d_R_bw_1 * alpha_arg + R_tau12k * (df_1_dbw_1 * w_x - f_1 * e_1x + df_2_dbw_1 * w_x_2 - f_2 * (e_1x * w_x + w_x * e_1x))) * a_hat -
+ H_al * skew_x((J_save * e_1)) * g_tau;
+ J_a.block(0, 1, 3, 1) +=
+ (d_R_bw_2 * alpha_arg + R_tau12k * (df_1_dbw_2 * w_x - f_1 * e_2x + df_2_dbw_2 * w_x_2 - f_2 * (e_2x * w_x + w_x * e_2x))) * a_hat -
+ H_al * skew_x((J_save * e_2)) * g_tau;
+ J_a.block(0, 2, 3, 1) +=
+ (d_R_bw_3 * alpha_arg + R_tau12k * (df_1_dbw_3 * w_x - f_1 * e_3x + df_2_dbw_3 * w_x_2 - f_2 * (e_3x * w_x + w_x * e_3x))) * a_hat -
+ H_al * skew_x((J_save * e_3)) * g_tau;
+ J_b.block(0, 0, 3, 1) +=
+ (d_R_bw_1 * Beta_arg + R_tau12k * (df_3_dbw_1 * w_x - f_3 * e_1x + df_4_dbw_1 * w_x_2 - f_4 * (e_1x * w_x + w_x * e_1x))) * a_hat -
+ -H_be * skew_x((J_save * e_1)) * g_tau;
+ J_b.block(0, 1, 3, 1) +=
+ (d_R_bw_2 * Beta_arg + R_tau12k * (df_3_dbw_2 * w_x - f_3 * e_2x + df_4_dbw_2 * w_x_2 - f_4 * (e_2x * w_x + w_x * e_2x))) * a_hat -
+ H_be * skew_x((J_save * e_2)) * g_tau;
+ J_b.block(0, 2, 3, 1) +=
+ (d_R_bw_3 * Beta_arg + R_tau12k * (df_3_dbw_3 * w_x - f_3 * e_3x + df_4_dbw_3 * w_x_2 - f_4 * (e_3x * w_x + w_x * e_3x))) * a_hat -
+ H_be * skew_x((J_save * e_3)) * g_tau;
+
+ //==========================================================================
+ // MEASUREMENT COVARIANCE
+ //==========================================================================
+
+ // Going to need orientation at intermediate time i.e. at .5*dt;
+ Eigen::Matrix R_G_to_k = quat_2_Rot(q_k_lin);
+ double dt_mid = delta_t / 2.0;
+ double w_dt_mid = mag_w * dt_mid;
+ Eigen::Matrix R_mid;
+
+ // The middle of this interval (i.e., rotation from k to mid)
+ R_mid = small_w ? eye3 - dt_mid * w_x + (pow(dt_mid, 2) / 2) * w_x_2
+ : eye3 - (sin(w_dt_mid) / mag_w) * w_x + ((1.0 - cos(w_dt_mid)) / (pow(mag_w, 2.0))) * w_x_2;
+ R_mid = R_mid * R_k2tau;
+
+ // Compute covariance (in this implementation, we use RK4)
+ // k1-------------------------------------------------------------------------------------------------
+
+ // Build state Jacobian
+ Eigen::Matrix F_k1 = Eigen::Matrix::Zero();
+ F_k1.block(0, 0, 3, 3) = -w_x;
+ F_k1.block(0, 3, 3, 3) = -eye3;
+ F_k1.block(6, 0, 3, 3) = -R_k2tau.transpose() * a_x;
+ F_k1.block(6, 9, 3, 3) = -R_k2tau.transpose();
+ F_k1.block(6, 15, 3, 3) = -R_k2tau.transpose() * skew_x(R_k2tau * R_G_to_k * grav);
+ F_k1.block(6, 18, 3, 3) = -R_k2tau.transpose() * R_k2tau * skew_x(R_G_to_k * grav);
+ F_k1.block(12, 6, 3, 3) = eye3;
+
+ // Build noise Jacobian
+ Eigen::Matrix G_k1 = Eigen::Matrix::Zero();
+ G_k1.block(0, 0, 3, 3) = -eye3;
+ G_k1.block(3, 3, 3, 3) = eye3;
+ G_k1.block(6, 6, 3, 3) = -R_k2tau.transpose();
+ G_k1.block(9, 9, 3, 3) = eye3;
+
+ // Get state transition and covariance derivative
+ Eigen::Matrix Phi_dot_k1 = F_k1;
+ Eigen::Matrix P_dot_k1 = F_k1 * P_big + P_big * F_k1.transpose() + G_k1 * Q_c * G_k1.transpose();
+
+ // k2-------------------------------------------------------------------------------------------------
+
+ // Build state Jacobian
+ Eigen::Matrix F_k2 = Eigen::Matrix::Zero();
+ F_k2.block(0, 0, 3, 3) = -w_x;
+ F_k2.block(0, 3, 3, 3) = -eye3;
+ F_k2.block(6, 0, 3, 3) = -R_mid.transpose() * a_x;
+ F_k2.block(6, 9, 3, 3) = -R_mid.transpose();
+ F_k2.block(6, 15, 3, 3) = -R_mid.transpose() * skew_x(R_k2tau * R_G_to_k * grav);
+ F_k2.block(6, 18, 3, 3) = -R_mid.transpose() * R_k2tau * skew_x(R_G_to_k * grav);
+ F_k2.block(12, 6, 3, 3) = eye3;
+
+ // Build noise Jacobian
+ Eigen::Matrix G_k2 = Eigen::Matrix::Zero();
+ G_k2.block(0, 0, 3, 3) = -eye3;
+ G_k2.block(3, 3, 3, 3) = eye3;
+ G_k2.block(6, 6, 3, 3) = -R_mid.transpose();
+ G_k2.block(9, 9, 3, 3) = eye3;
+
+ // Get state transition and covariance derivative
+ Eigen::Matrix Phi_k2 = Eigen::Matrix::Identity() + Phi_dot_k1 * dt_mid;
+ Eigen::Matrix P_k2 = P_big + P_dot_k1 * dt_mid;
+ Eigen::Matrix Phi_dot_k2 = F_k2 * Phi_k2;
+ Eigen::Matrix P_dot_k2 = F_k2 * P_k2 + P_k2 * F_k2.transpose() + G_k2 * Q_c * G_k2.transpose();
+
+ // k3-------------------------------------------------------------------------------------------------
+
+ // Our state and noise Jacobians are the same as k2
+ // Since k2 and k3 correspond to the same estimates for the midpoint
+ Eigen::Matrix F_k3 = F_k2;
+ Eigen::Matrix G_k3 = G_k2;
+
+ // Get state transition and covariance derivative
+ Eigen::Matrix Phi_k3 = Eigen::Matrix::Identity() + Phi_dot_k2 * dt_mid;
+ Eigen::Matrix P_k3 = P_big + P_dot_k2 * dt_mid;
+ Eigen::Matrix Phi_dot_k3 = F_k3 * Phi_k3;
+ Eigen::Matrix P_dot_k3 = F_k3 * P_k3 + P_k3 * F_k3.transpose() + G_k3 * Q_c * G_k3.transpose();
+
+ // k4-------------------------------------------------------------------------------------------------
+
+ // Build state Jacobian
+ Eigen::Matrix F_k4 = Eigen::Matrix::Zero();
+ F_k4.block(0, 0, 3, 3) = -w_x;
+ F_k4.block(0, 3, 3, 3) = -eye3;
+ F_k4.block(6, 0, 3, 3) = -R_k2tau1.transpose() * a_x;
+ F_k4.block(6, 9, 3, 3) = -R_k2tau1.transpose();
+ F_k4.block(6, 15, 3, 3) = -R_k2tau1.transpose() * skew_x(R_k2tau * R_G_to_k * grav);
+ F_k4.block(6, 18, 3, 3) = -R_k2tau1.transpose() * R_k2tau * skew_x(R_G_to_k * grav);
+ F_k4.block(12, 6, 3, 3) = eye3;
+
+ // Build noise Jacobian
+ Eigen::Matrix G_k4 = Eigen::Matrix::Zero();
+ G_k4.block(0, 0, 3, 3) = -eye3;
+ G_k4.block(3, 3, 3, 3) = eye3;
+ G_k4.block(6, 6, 3, 3) = -R_k2tau1.transpose();
+ G_k4.block(9, 9, 3, 3) = eye3;
+
+ // Get state transition and covariance derivative
+ Eigen::Matrix Phi_k4 = Eigen::Matrix::Identity() + Phi_dot_k3 * delta_t;
+ Eigen::Matrix P_k4 = P_big + P_dot_k3 * delta_t;
+ Eigen::Matrix Phi_dot_k4 = F_k4 * Phi_k4;
+ Eigen::Matrix P_dot_k4 = F_k4 * P_k4 + P_k4 * F_k4.transpose() + G_k4 * Q_c * G_k4.transpose();
+
+ // done-------------------------------------------------------------------------------------------------
+
+ // Collect covariance solution
+ // Ensure it is positive definite
+ P_big += (delta_t / 6.0) * (P_dot_k1 + 2.0 * P_dot_k2 + 2.0 * P_dot_k3 + P_dot_k4);
+ P_big = 0.5 * (P_big + P_big.transpose());
+
+ // Calculate the state transition from time k to tau
+ Eigen::Matrix Phi =
+ Eigen::Matrix::Identity() + (delta_t / 6.0) * (Phi_dot_k1 + 2.0 * Phi_dot_k2 + 2.0 * Phi_dot_k3 + Phi_dot_k4);
+
+ //==========================================================================
+ // CLONE TO NEW SAMPLE TIME AND MARGINALIZE OLD SAMPLE TIME
+ //==========================================================================
+
+ // Matrix that performs the clone and mariginalization
+ Eigen::Matrix B_k = Eigen::Matrix::Identity();
+ B_k.block(15, 15, 3, 3).setZero();
+ B_k.block(15, 0, 3, 3) = Eigen::Matrix::Identity();
+
+ // Change our measurement covariance and state transition
+ P_big = B_k * P_big * B_k.transpose();
+ P_big = 0.5 * (P_big + P_big.transpose());
+ Discrete_J_b = B_k * Phi * Discrete_J_b;
+
+ // Our measurement covariance is the top 15x15 of our large covariance
+ P_meas = P_big.block(0, 0, 15, 15);
+
+ // If we are using the state transition Jacobian, then we should overwrite the analytical versions
+ // Note: we flip the sign for J_q to match the Model 1 derivation
+ if (state_transition_jacobians) {
+ J_q = -Discrete_J_b.block(0, 3, 3, 3);
+ J_a = Discrete_J_b.block(12, 3, 3, 3);
+ J_b = Discrete_J_b.block(6, 3, 3, 3);
+ H_a = Discrete_J_b.block(12, 9, 3, 3);
+ H_b = Discrete_J_b.block(6, 9, 3, 3);
+ O_a = Discrete_J_b.block(12, 18, 3, 3);
+ O_b = Discrete_J_b.block(6, 18, 3, 3);
+ }
+
+ // Update rotation mean
+ // Note we had to wait to do this, since we use the old orientation in our covariance calculation
+ R_k2tau = R_k2tau1;
+ q_k2tau = rot_2_quat(R_k2tau);
+}
diff --git a/ov_core/src/cpi/CpiV2.h b/ov_core/src/cpi/CpiV2.h
index 8a358501e..44594c124 100644
--- a/ov_core/src/cpi/CpiV2.h
+++ b/ov_core/src/cpi/CpiV2.h
@@ -27,7 +27,7 @@
*/
#include "CpiBase.h"
-#include "utils/quat_ops.h"
+
#include
namespace ov_core {
@@ -79,6 +79,8 @@ class CpiV2 : public CpiBase {
CpiV2(double sigma_w, double sigma_wb, double sigma_a, double sigma_ab, bool imu_avg_ = false)
: CpiBase(sigma_w, sigma_wb, sigma_a, sigma_ab, imu_avg_) {}
+ virtual ~CpiV2() {}
+
/**
* @brief Our precompound function for Model 2
* @param[in] t_0 first IMU timestamp
@@ -94,360 +96,7 @@ class CpiV2 : public CpiBase {
*/
void feed_IMU(double t_0, double t_1, Eigen::Matrix w_m_0, Eigen::Matrix a_m_0,
Eigen::Matrix w_m_1 = Eigen::Matrix::Zero(),
- Eigen::Matrix a_m_1 = Eigen::Matrix::Zero()) {
-
- // Get time difference
- double delta_t = t_1 - t_0;
- DT += delta_t;
-
- // If no time has passed do nothing
- if (delta_t == 0) {
- return;
- }
-
- // Get estimated imu readings
- Eigen::Matrix w_hat = w_m_0 - b_w_lin;
- Eigen::Matrix a_hat = a_m_0 - b_a_lin - R_k2tau * quat_2_Rot(q_k_lin) * grav;
-
- // If averaging, average
- // Note: we will average the LOCAL acceleration after getting the relative rotation
- if (imu_avg) {
- w_hat += w_m_1 - b_w_lin;
- w_hat = 0.5 * w_hat;
- }
-
- // Get angle change w*dt
- Eigen::Matrix w_hatdt = w_hat * delta_t;
-
- // Get entries of w_hat
- double w_1 = w_hat(0, 0);
- double w_2 = w_hat(1, 0);
- double w_3 = w_hat(2, 0);
-
- // Get magnitude of w and wdt
- double mag_w = w_hat.norm();
- double w_dt = mag_w * delta_t;
-
- // Threshold to determine if equations will be unstable
- bool small_w = (mag_w < 0.008726646);
-
- // Get some of the variables used in the preintegration equations
- double dt_2 = pow(delta_t, 2);
- double cos_wt = cos(w_dt);
- double sin_wt = sin(w_dt);
-
- Eigen::Matrix w_x = skew_x(w_hat);
- Eigen::Matrix w_tx = skew_x(w_hatdt);
- Eigen::Matrix w_x_2 = w_x * w_x;
-
- //==========================================================================
- // MEASUREMENT MEANS
- //==========================================================================
-
- // Get relative rotation
- Eigen::Matrix R_tau2tau1 = small_w ? eye3 - delta_t * w_x + (pow(delta_t, 2) / 2) * w_x_2
- : eye3 - (sin_wt / mag_w) * w_x + ((1.0 - cos_wt) / (pow(mag_w, 2.0))) * w_x_2;
-
- // Updated roation and its transpose
- Eigen::Matrix R_k2tau1 = R_tau2tau1 * R_k2tau;
- Eigen::Matrix R_tau12k = R_k2tau1.transpose();
-
- // If averaging, average the LOCAL acceleration
- if (imu_avg) {
- a_hat += a_m_1 - b_a_lin - R_k2tau1 * quat_2_Rot(q_k_lin) * grav;
- a_hat = 0.5 * a_hat;
- }
- Eigen::Matrix a_x = skew_x(a_hat);
-
- // Intermediate variables for evaluating the measurement/bias Jacobian update
- double f_1;
- double f_2;
- double f_3;
- double f_4;
-
- if (small_w) {
- f_1 = -(pow(delta_t, 3) / 3);
- f_2 = (pow(delta_t, 4) / 8);
- f_3 = -(pow(delta_t, 2) / 2);
- f_4 = (pow(delta_t, 3) / 6);
- } else {
- f_1 = (w_dt * cos_wt - sin_wt) / (pow(mag_w, 3));
- f_2 = (pow(w_dt, 2) - 2 * cos_wt - 2 * w_dt * sin_wt + 2) / (2 * pow(mag_w, 4));
- f_3 = -(1 - cos_wt) / pow(mag_w, 2);
- f_4 = (w_dt - sin_wt) / pow(mag_w, 3);
- }
-
- // Compute the main part of our analytical means
- Eigen::Matrix alpha_arg = ((dt_2 / 2.0) * eye3 + f_1 * w_x + f_2 * w_x_2);
- Eigen::Matrix Beta_arg = (delta_t * eye3 + f_3 * w_x + f_4 * w_x_2);
-
- // Matrices that will multiply the a_hat in the update expressions
- Eigen::Matrix H_al = R_tau12k * alpha_arg;
- Eigen::Matrix H_be = R_tau12k * Beta_arg;
-
- // Update the measurements
- alpha_tau += beta_tau * delta_t + H_al * a_hat;
- beta_tau += H_be * a_hat;
-
- //==========================================================================
- // BIAS JACOBIANS (ANALYTICAL)
- //==========================================================================
-
- // Get right Jacobian
- Eigen::Matrix J_r_tau1 =
- small_w ? eye3 - .5 * w_tx + (1.0 / 6.0) * w_tx * w_tx
- : eye3 - ((1 - cos_wt) / (pow((w_dt), 2.0))) * w_tx + ((w_dt - sin_wt) / (pow(w_dt, 3.0))) * w_tx * w_tx;
-
- // Update orientation in respect to gyro bias Jacobians
- Eigen::Matrix J_save = J_q;
- J_q = R_tau2tau1 * J_q + J_r_tau1 * delta_t;
-
- // Update alpha and beta in respect to accel bias Jacobian
- H_a -= H_al;
- H_a += delta_t * H_b;
- H_b -= H_be;
-
- // Update alpha and beta in respect to q_GtoLIN Jacobian
- Eigen::Matrix g_k = quat_2_Rot(q_k_lin) * grav;
- O_a += delta_t * O_b;
- O_a += -H_al * R_k2tau * skew_x(g_k);
- O_b += -H_be * R_k2tau * skew_x(g_k);
-
- // Derivatives of R_tau12k wrt bias_w entries
- Eigen::MatrixXd d_R_bw_1 = -R_tau12k * skew_x(J_q * e_1);
- Eigen::MatrixXd d_R_bw_2 = -R_tau12k * skew_x(J_q * e_2);
- Eigen::MatrixXd d_R_bw_3 = -R_tau12k * skew_x(J_q * e_3);
-
- // Now compute the gyro bias Jacobian terms
- double df_1_dbw_1;
- double df_1_dbw_2;
- double df_1_dbw_3;
-
- double df_2_dbw_1;
- double df_2_dbw_2;
- double df_2_dbw_3;
-
- double df_3_dbw_1;
- double df_3_dbw_2;
- double df_3_dbw_3;
-
- double df_4_dbw_1;
- double df_4_dbw_2;
- double df_4_dbw_3;
-
- if (small_w) {
- double df_1_dw_mag = -(pow(delta_t, 5) / 15);
- df_1_dbw_1 = w_1 * df_1_dw_mag;
- df_1_dbw_2 = w_2 * df_1_dw_mag;
- df_1_dbw_3 = w_3 * df_1_dw_mag;
-
- double df_2_dw_mag = (pow(delta_t, 6) / 72);
- df_2_dbw_1 = w_1 * df_2_dw_mag;
- df_2_dbw_2 = w_2 * df_2_dw_mag;
- df_2_dbw_3 = w_3 * df_2_dw_mag;
-
- double df_3_dw_mag = -(pow(delta_t, 4) / 12);
- df_3_dbw_1 = w_1 * df_3_dw_mag;
- df_3_dbw_2 = w_2 * df_3_dw_mag;
- df_3_dbw_3 = w_3 * df_3_dw_mag;
-
- double df_4_dw_mag = (pow(delta_t, 5) / 60);
- df_4_dbw_1 = w_1 * df_4_dw_mag;
- df_4_dbw_2 = w_2 * df_4_dw_mag;
- df_4_dbw_3 = w_3 * df_4_dw_mag;
- } else {
- double df_1_dw_mag = (pow(w_dt, 2) * sin_wt - 3 * sin_wt + 3 * w_dt * cos_wt) / pow(mag_w, 5);
- df_1_dbw_1 = w_1 * df_1_dw_mag;
- df_1_dbw_2 = w_2 * df_1_dw_mag;
- df_1_dbw_3 = w_3 * df_1_dw_mag;
-
- double df_2_dw_mag = (pow(w_dt, 2) - 4 * cos_wt - 4 * w_dt * sin_wt + pow(w_dt, 2) * cos_wt + 4) / (pow(mag_w, 6));
- df_2_dbw_1 = w_1 * df_2_dw_mag;
- df_2_dbw_2 = w_2 * df_2_dw_mag;
- df_2_dbw_3 = w_3 * df_2_dw_mag;
-
- double df_3_dw_mag = (2 * (cos_wt - 1) + w_dt * sin_wt) / (pow(mag_w, 4));
- df_3_dbw_1 = w_1 * df_3_dw_mag;
- df_3_dbw_2 = w_2 * df_3_dw_mag;
- df_3_dbw_3 = w_3 * df_3_dw_mag;
-
- double df_4_dw_mag = (2 * w_dt + w_dt * cos_wt - 3 * sin_wt) / (pow(mag_w, 5));
- df_4_dbw_1 = w_1 * df_4_dw_mag;
- df_4_dbw_2 = w_2 * df_4_dw_mag;
- df_4_dbw_3 = w_3 * df_4_dw_mag;
- }
-
- // Gravity rotated into the "tau'th" frame (i.e. start of the measurement interval)
- Eigen::Matrix g_tau = R_k2tau * quat_2_Rot(q_k_lin) * grav;
-
- // Update gyro bias Jacobians
- J_a += J_b * delta_t;
- J_a.block(0, 0, 3, 1) +=
- (d_R_bw_1 * alpha_arg + R_tau12k * (df_1_dbw_1 * w_x - f_1 * e_1x + df_2_dbw_1 * w_x_2 - f_2 * (e_1x * w_x + w_x * e_1x))) * a_hat -
- H_al * skew_x((J_save * e_1)) * g_tau;
- J_a.block(0, 1, 3, 1) +=
- (d_R_bw_2 * alpha_arg + R_tau12k * (df_1_dbw_2 * w_x - f_1 * e_2x + df_2_dbw_2 * w_x_2 - f_2 * (e_2x * w_x + w_x * e_2x))) * a_hat -
- H_al * skew_x((J_save * e_2)) * g_tau;
- J_a.block(0, 2, 3, 1) +=
- (d_R_bw_3 * alpha_arg + R_tau12k * (df_1_dbw_3 * w_x - f_1 * e_3x + df_2_dbw_3 * w_x_2 - f_2 * (e_3x * w_x + w_x * e_3x))) * a_hat -
- H_al * skew_x((J_save * e_3)) * g_tau;
- J_b.block(0, 0, 3, 1) +=
- (d_R_bw_1 * Beta_arg + R_tau12k * (df_3_dbw_1 * w_x - f_3 * e_1x + df_4_dbw_1 * w_x_2 - f_4 * (e_1x * w_x + w_x * e_1x))) * a_hat -
- -H_be * skew_x((J_save * e_1)) * g_tau;
- J_b.block(0, 1, 3, 1) +=
- (d_R_bw_2 * Beta_arg + R_tau12k * (df_3_dbw_2 * w_x - f_3 * e_2x + df_4_dbw_2 * w_x_2 - f_4 * (e_2x * w_x + w_x * e_2x))) * a_hat -
- H_be * skew_x((J_save * e_2)) * g_tau;
- J_b.block(0, 2, 3, 1) +=
- (d_R_bw_3 * Beta_arg + R_tau12k * (df_3_dbw_3 * w_x - f_3 * e_3x + df_4_dbw_3 * w_x_2 - f_4 * (e_3x * w_x + w_x * e_3x))) * a_hat -
- H_be * skew_x((J_save * e_3)) * g_tau;
-
- //==========================================================================
- // MEASUREMENT COVARIANCE
- //==========================================================================
-
- // Going to need orientation at intermediate time i.e. at .5*dt;
- Eigen::Matrix R_G_to_k = quat_2_Rot(q_k_lin);
- double dt_mid = delta_t / 2.0;
- double w_dt_mid = mag_w * dt_mid;
- Eigen::Matrix R_mid;
-
- // The middle of this interval (i.e., rotation from k to mid)
- R_mid = small_w ? eye3 - dt_mid * w_x + (pow(dt_mid, 2) / 2) * w_x_2
- : eye3 - (sin(w_dt_mid) / mag_w) * w_x + ((1.0 - cos(w_dt_mid)) / (pow(mag_w, 2.0))) * w_x_2;
- R_mid = R_mid * R_k2tau;
-
- // Compute covariance (in this implementation, we use RK4)
- // k1-------------------------------------------------------------------------------------------------
-
- // Build state Jacobian
- Eigen::Matrix F_k1 = Eigen::Matrix::Zero();
- F_k1.block(0, 0, 3, 3) = -w_x;
- F_k1.block(0, 3, 3, 3) = -eye3;
- F_k1.block(6, 0, 3, 3) = -R_k2tau.transpose() * a_x;
- F_k1.block(6, 9, 3, 3) = -R_k2tau.transpose();
- F_k1.block(6, 15, 3, 3) = -R_k2tau.transpose() * skew_x(R_k2tau * R_G_to_k * grav);
- F_k1.block(6, 18, 3, 3) = -R_k2tau.transpose() * R_k2tau * skew_x(R_G_to_k * grav);
- F_k1.block(12, 6, 3, 3) = eye3;
-
- // Build noise Jacobian
- Eigen::Matrix G_k1 = Eigen::Matrix::Zero();
- G_k1.block(0, 0, 3, 3) = -eye3;
- G_k1.block(3, 3, 3, 3) = eye3;
- G_k1.block(6, 6, 3, 3) = -R_k2tau.transpose();
- G_k1.block(9, 9, 3, 3) = eye3;
-
- // Get state transition and covariance derivative
- Eigen::Matrix Phi_dot_k1 = F_k1;
- Eigen::Matrix P_dot_k1 = F_k1 * P_big + P_big * F_k1.transpose() + G_k1 * Q_c * G_k1.transpose();
-
- // k2-------------------------------------------------------------------------------------------------
-
- // Build state Jacobian
- Eigen::Matrix F_k2 = Eigen::Matrix::Zero();
- F_k2.block(0, 0, 3, 3) = -w_x;
- F_k2.block(0, 3, 3, 3) = -eye3;
- F_k2.block(6, 0, 3, 3) = -R_mid.transpose() * a_x;
- F_k2.block(6, 9, 3, 3) = -R_mid.transpose();
- F_k2.block(6, 15, 3, 3) = -R_mid.transpose() * skew_x(R_k2tau * R_G_to_k * grav);
- F_k2.block(6, 18, 3, 3) = -R_mid.transpose() * R_k2tau * skew_x(R_G_to_k * grav);
- F_k2.block(12, 6, 3, 3) = eye3;
-
- // Build noise Jacobian
- Eigen::Matrix G_k2 = Eigen::Matrix::Zero();
- G_k2.block(0, 0, 3, 3) = -eye3;
- G_k2.block(3, 3, 3, 3) = eye3;
- G_k2.block(6, 6, 3, 3) = -R_mid.transpose();
- G_k2.block(9, 9, 3, 3) = eye3;
-
- // Get state transition and covariance derivative
- Eigen::Matrix Phi_k2 = Eigen::Matrix::Identity() + Phi_dot_k1 * dt_mid;
- Eigen::Matrix P_k2 = P_big + P_dot_k1 * dt_mid;
- Eigen::Matrix Phi_dot_k2 = F_k2 * Phi_k2;
- Eigen::Matrix P_dot_k2 = F_k2 * P_k2 + P_k2 * F_k2.transpose() + G_k2 * Q_c * G_k2.transpose();
-
- // k3-------------------------------------------------------------------------------------------------
-
- // Our state and noise Jacobians are the same as k2
- // Since k2 and k3 correspond to the same estimates for the midpoint
- Eigen::Matrix F_k3 = F_k2;
- Eigen::Matrix G_k3 = G_k2;
-
- // Get state transition and covariance derivative
- Eigen::Matrix Phi_k3 = Eigen::Matrix::Identity() + Phi_dot_k2 * dt_mid;
- Eigen::Matrix P_k3 = P_big + P_dot_k2 * dt_mid;
- Eigen::Matrix Phi_dot_k3 = F_k3 * Phi_k3;
- Eigen::Matrix P_dot_k3 = F_k3 * P_k3 + P_k3 * F_k3.transpose() + G_k3 * Q_c * G_k3.transpose();
-
- // k4-------------------------------------------------------------------------------------------------
-
- // Build state Jacobian
- Eigen::Matrix F_k4 = Eigen::Matrix::Zero();
- F_k4.block(0, 0, 3, 3) = -w_x;
- F_k4.block(0, 3, 3, 3) = -eye3;
- F_k4.block(6, 0, 3, 3) = -R_k2tau1.transpose() * a_x;
- F_k4.block(6, 9, 3, 3) = -R_k2tau1.transpose();
- F_k4.block(6, 15, 3, 3) = -R_k2tau1.transpose() * skew_x(R_k2tau * R_G_to_k * grav);
- F_k4.block(6, 18, 3, 3) = -R_k2tau1.transpose() * R_k2tau * skew_x(R_G_to_k * grav);
- F_k4.block(12, 6, 3, 3) = eye3;
-
- // Build noise Jacobian
- Eigen::Matrix G_k4 = Eigen::Matrix::Zero();
- G_k4.block(0, 0, 3, 3) = -eye3;
- G_k4.block(3, 3, 3, 3) = eye3;
- G_k4.block(6, 6, 3, 3) = -R_k2tau1.transpose();
- G_k4.block(9, 9, 3, 3) = eye3;
-
- // Get state transition and covariance derivative
- Eigen::Matrix Phi_k4 = Eigen::Matrix::Identity() + Phi_dot_k3 * delta_t;
- Eigen::Matrix P_k4 = P_big + P_dot_k3 * delta_t;
- Eigen::Matrix Phi_dot_k4 = F_k4 * Phi_k4;
- Eigen::Matrix P_dot_k4 = F_k4 * P_k4 + P_k4 * F_k4.transpose() + G_k4 * Q_c * G_k4.transpose();
-
- // done-------------------------------------------------------------------------------------------------
-
- // Collect covariance solution
- // Ensure it is positive definite
- P_big += (delta_t / 6.0) * (P_dot_k1 + 2.0 * P_dot_k2 + 2.0 * P_dot_k3 + P_dot_k4);
- P_big = 0.5 * (P_big + P_big.transpose());
-
- // Calculate the state transition from time k to tau
- Eigen::Matrix Phi =
- Eigen::Matrix::Identity() + (delta_t / 6.0) * (Phi_dot_k1 + 2.0 * Phi_dot_k2 + 2.0 * Phi_dot_k3 + Phi_dot_k4);
-
- //==========================================================================
- // CLONE TO NEW SAMPLE TIME AND MARGINALIZE OLD SAMPLE TIME
- //==========================================================================
-
- // Matrix that performs the clone and mariginalization
- Eigen::Matrix B_k = Eigen::Matrix::Identity();
- B_k.block(15, 15, 3, 3).setZero();
- B_k.block(15, 0, 3, 3) = Eigen::Matrix::Identity();
-
- // Change our measurement covariance and state transition
- P_big = B_k * P_big * B_k.transpose();
- P_big = 0.5 * (P_big + P_big.transpose());
- Discrete_J_b = B_k * Phi * Discrete_J_b;
-
- // Our measurement covariance is the top 15x15 of our large covariance
- P_meas = P_big.block(0, 0, 15, 15);
-
- // If we are using the state transition Jacobian, then we should overwrite the analytical versions
- // Note: we flip the sign for J_q to match the Model 1 derivation
- if (state_transition_jacobians) {
- J_q = -Discrete_J_b.block(0, 3, 3, 3);
- J_a = Discrete_J_b.block(12, 3, 3, 3);
- J_b = Discrete_J_b.block(6, 3, 3, 3);
- H_a = Discrete_J_b.block(12, 9, 3, 3);
- H_b = Discrete_J_b.block(6, 9, 3, 3);
- O_a = Discrete_J_b.block(12, 18, 3, 3);
- O_b = Discrete_J_b.block(6, 18, 3, 3);
- }
-
- // Update rotation mean
- // Note we had to wait to do this, since we use the old orientation in our covariance calculation
- R_k2tau = R_k2tau1;
- q_k2tau = rot_2_quat(R_k2tau);
- }
+ Eigen::Matrix a_m_1 = Eigen::Matrix::Zero());
};
} // namespace ov_core
diff --git a/ov_core/src/feat/FeatureDatabase.cpp b/ov_core/src/feat/FeatureDatabase.cpp
new file mode 100644
index 000000000..5cb8de5ce
--- /dev/null
+++ b/ov_core/src/feat/FeatureDatabase.cpp
@@ -0,0 +1,321 @@
+/*
+ * OpenVINS: An Open Platform for Visual-Inertial Research
+ * Copyright (C) 2018-2022 Patrick Geneva
+ * Copyright (C) 2018-2022 Guoquan Huang
+ * Copyright (C) 2018-2022 OpenVINS Contributors
+ * Copyright (C) 2018-2019 Kevin Eckenhoff
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see .
+ */
+
+#include "FeatureDatabase.h"
+
+#include "Feature.h"
+#include "utils/print.h"
+
+using namespace ov_core;
+
+std::shared_ptr FeatureDatabase::get_feature(size_t id, bool remove) {
+ std::lock_guard lck(mtx);
+ if (features_idlookup.find(id) != features_idlookup.end()) {
+ std::shared_ptr temp = features_idlookup.at(id);
+ if (remove)
+ features_idlookup.erase(id);
+ return temp;
+ } else {
+ return nullptr;
+ }
+}
+
+bool FeatureDatabase::get_feature_clone(size_t id, Feature &feat) {
+ std::lock_guard lck(mtx);
+ if (features_idlookup.find(id) == features_idlookup.end())
+ return false;
+ // TODO: should probably have a copy constructor function in feature class
+ std::shared_ptr temp = features_idlookup.at(id);
+ feat.featid = temp->featid;
+ feat.to_delete = temp->to_delete;
+ feat.uvs = temp->uvs;
+ feat.uvs_norm = temp->uvs_norm;
+ feat.timestamps = temp->timestamps;
+ feat.anchor_cam_id = temp->anchor_cam_id;
+ feat.anchor_clone_timestamp = temp->anchor_clone_timestamp;
+ feat.p_FinA = temp->p_FinA;
+ feat.p_FinG = temp->p_FinG;
+ return true;
+}
+
+void FeatureDatabase::update_feature(size_t id, double timestamp, size_t cam_id, float u, float v, float u_n, float v_n) {
+
+ // Find this feature using the ID lookup
+ std::lock_guard lck(mtx);
+ if (features_idlookup.find(id) != features_idlookup.end()) {
+ // Get our feature
+ std::shared_ptr feat = features_idlookup.at(id);
+ // Append this new information to it!
+ feat->uvs[cam_id].push_back(Eigen::Vector2f(u, v));
+ feat->uvs_norm[cam_id].push_back(Eigen::Vector2f(u_n, v_n));
+ feat->timestamps[cam_id].push_back(timestamp);
+ return;
+ }
+
+ // Debug info
+ // PRINT_DEBUG("featdb - adding new feature %d",(int)id);
+
+ // Else we have not found the feature, so lets make it be a new one!
+ std::shared_ptr feat = std::make_shared();
+ feat->featid = id;
+ feat->uvs[cam_id].push_back(Eigen::Vector2f(u, v));
+ feat->uvs_norm[cam_id].push_back(Eigen::Vector2f(u_n, v_n));
+ feat->timestamps[cam_id].push_back(timestamp);
+
+ // Append this new feature into our database
+ features_idlookup[id] = feat;
+}
+
+std::vector> FeatureDatabase::features_not_containing_newer(double timestamp, bool remove, bool skip_deleted) {
+
+ // Our vector of features that do not have measurements after the specified time
+ std::vector> feats_old;
+
+ // Now lets loop through all features, and just make sure they are not old
+ std::lock_guard lck(mtx);
+ for (auto it = features_idlookup.begin(); it != features_idlookup.end();) {
+ // Skip if already deleted
+ if (skip_deleted && (*it).second->to_delete) {
+ it++;
+ continue;
+ }
+ // Loop through each camera
+ // If we have a measurement greater-than or equal to the specified, this measurement is find
+ bool has_newer_measurement = false;
+ for (auto const &pair : (*it).second->timestamps) {
+ has_newer_measurement = (!pair.second.empty() && pair.second.at(pair.second.size() - 1) >= timestamp);
+ if (has_newer_measurement) {
+ break;
+ }
+ }
+ // If it is not being actively tracked, then it is old
+ if (!has_newer_measurement) {
+ feats_old.push_back((*it).second);
+ if (remove)
+ features_idlookup.erase(it++);
+ else
+ it++;
+ } else {
+ it++;
+ }
+ }
+
+ // Debugging
+ // PRINT_DEBUG("feature db size = %u\n", features_idlookup.size())
+
+ // Return the old features
+ return feats_old;
+}
+
+std::vector> FeatureDatabase::features_containing_older(double timestamp, bool remove, bool skip_deleted) {
+
+ // Our vector of old features
+ std::vector> feats_old;
+
+ // Now lets loop through all features, and just make sure they are not old
+ std::lock_guard lck(mtx);
+ for (auto it = features_idlookup.begin(); it != features_idlookup.end();) {
+ // Skip if already deleted
+ if (skip_deleted && (*it).second->to_delete) {
+ it++;
+ continue;
+ }
+ // Loop through each camera
+ // Check if we have at least one time older then the requested
+ bool found_containing_older = false;
+ for (auto const &pair : (*it).second->timestamps) {
+ found_containing_older = (!pair.second.empty() && pair.second.at(0) < timestamp);
+ if (found_containing_older) {
+ break;
+ }
+ }
+ // If it has an older timestamp, then add it
+ if (found_containing_older) {
+ feats_old.push_back((*it).second);
+ if (remove)
+ features_idlookup.erase(it++);
+ else
+ it++;
+ } else {
+ it++;
+ }
+ }
+
+ // Debugging
+ // PRINT_DEBUG("feature db size = %u\n", features_idlookup.size())
+
+ // Return the old features
+ return feats_old;
+}
+
+std::vector> FeatureDatabase::features_containing(double timestamp, bool remove, bool skip_deleted) {
+
+ // Our vector of old features
+ std::vector> feats_has_timestamp;
+
+ // Now lets loop through all features, and just make sure they are not
+ std::lock_guard lck(mtx);
+ for (auto it = features_idlookup.begin(); it != features_idlookup.end();) {
+ // Skip if already deleted
+ if (skip_deleted && (*it).second->to_delete) {
+ it++;
+ continue;
+ }
+ // Boolean if it has the timestamp
+ // Break out if we found a single timestamp that is equal to the specified time
+ bool has_timestamp = false;
+ for (auto const &pair : (*it).second->timestamps) {
+ has_timestamp = (std::find(pair.second.begin(), pair.second.end(), timestamp) != pair.second.end());
+ if (has_timestamp) {
+ break;
+ }
+ }
+ // Remove this feature if it contains the specified timestamp
+ if (has_timestamp) {
+ feats_has_timestamp.push_back((*it).second);
+ if (remove)
+ features_idlookup.erase(it++);
+ else
+ it++;
+ } else {
+ it++;
+ }
+ }
+
+ // Debugging
+ // PRINT_DEBUG("feature db size = %u\n", features_idlookup.size())
+ // PRINT_DEBUG("return vector = %u\n", feats_has_timestamp.size())
+
+ // Return the features
+ return feats_has_timestamp;
+}
+
+void FeatureDatabase::cleanup() {
+ // Loop through all features
+ // int sizebefore = (int)features_idlookup.size();
+ std::lock_guard lck(mtx);
+ for (auto it = features_idlookup.begin(); it != features_idlookup.end();) {
+ // If delete flag is set, then delete it
+ if ((*it).second->to_delete) {
+ features_idlookup.erase(it++);
+ } else {
+ it++;
+ }
+ }
+ // PRINT_DEBUG("feat db = %d -> %d\n", sizebefore, (int)features_idlookup.size() << std::endl;
+}
+
+void FeatureDatabase::cleanup_measurements(double timestamp) {
+ std::lock_guard lck(mtx);
+ for (auto it = features_idlookup.begin(); it != features_idlookup.end();) {
+ // Remove the older measurements
+ (*it).second->clean_older_measurements(timestamp);
+ // Count how many measurements
+ int ct_meas = 0;
+ for (const auto &pair : (*it).second->timestamps) {
+ ct_meas += (int)(pair.second.size());
+ }
+ // If delete flag is set, then delete it
+ if (ct_meas < 1) {
+ features_idlookup.erase(it++);
+ } else {
+ it++;
+ }
+ }
+}
+
+void FeatureDatabase::cleanup_measurements_exact(double timestamp) {
+ std::lock_guard lck(mtx);
+ std::vector timestamps = {timestamp};
+ for (auto it = features_idlookup.begin(); it != features_idlookup.end();) {
+ // Remove the older measurements
+ (*it).second->clean_invalid_measurements(timestamps);
+ // Count how many measurements
+ int ct_meas = 0;
+ for (const auto &pair : (*it).second->timestamps) {
+ ct_meas += (int)(pair.second.size());
+ }
+ // If delete flag is set, then delete it
+ if (ct_meas < 1) {
+ features_idlookup.erase(it++);
+ } else {
+ it++;
+ }
+ }
+}
+
+double FeatureDatabase::get_oldest_timestamp() {
+ std::lock_guard lck(mtx);
+ double oldest_time = -1;
+ for (auto const &feat : features_idlookup) {
+ for (auto const &camtimepair : feat.second->timestamps) {
+ if (!camtimepair.second.empty() && (oldest_time == -1 || oldest_time < camtimepair.second.at(0))) {
+ oldest_time = camtimepair.second.at(0);
+ }
+ }
+ }
+ return oldest_time;
+}
+
+void FeatureDatabase::append_new_measurements(const std::shared_ptr &database) {
+ std::lock_guard lck(mtx);
+
+ // Loop through the other database's internal database
+ // int sizebefore = (int)features_idlookup.size();
+ for (const auto &feat : database->get_internal_data()) {
+ if (features_idlookup.find(feat.first) != features_idlookup.end()) {
+
+ // For this feature, now try to append the new measurement data
+ std::shared_ptr temp = features_idlookup.at(feat.first);
+ for (const auto × : feat.second->timestamps) {
+ // Append the whole camera vector is not seen
+ // Otherwise need to loop through each and append
+ size_t cam_id = times.first;
+ if (temp->timestamps.find(cam_id) == temp->timestamps.end()) {
+ temp->timestamps[cam_id] = feat.second->timestamps.at(cam_id);
+ temp->uvs[cam_id] = feat.second->uvs.at(cam_id);
+ temp->uvs_norm[cam_id] = feat.second->uvs_norm.at(cam_id);
+ } else {
+ auto temp_times = temp->timestamps.at(cam_id);
+ for (size_t i = 0; i < feat.second->timestamps.at(cam_id).size(); i++) {
+ double time_to_find = feat.second->timestamps.at(cam_id).at(i);
+ if (std::find(temp_times.begin(), temp_times.end(), time_to_find) == temp_times.end()) {
+ temp->timestamps.at(cam_id).push_back(feat.second->timestamps.at(cam_id).at(i));
+ temp->uvs.at(cam_id).push_back(feat.second->uvs.at(cam_id).at(i));
+ temp->uvs_norm.at(cam_id).push_back(feat.second->uvs_norm.at(cam_id).at(i));
+ }
+ }
+ }
+ }
+
+ } else {
+
+ // Else we have not found the feature, so lets make it be a new one!
+ std::shared_ptr temp = std::make_shared();
+ temp->featid = feat.second->featid;
+ temp->timestamps = feat.second->timestamps;
+ temp->uvs = feat.second->uvs;
+ temp->uvs_norm = feat.second->uvs_norm;
+ features_idlookup[feat.first] = temp;
+ }
+ }
+ // PRINT_DEBUG("feat db = %d -> %d\n", sizebefore, (int)features_idlookup.size() << std::endl;
+}
diff --git a/ov_core/src/feat/FeatureDatabase.h b/ov_core/src/feat/FeatureDatabase.h
index 46c65fcf0..7eae4de5c 100644
--- a/ov_core/src/feat/FeatureDatabase.h
+++ b/ov_core/src/feat/FeatureDatabase.h
@@ -25,13 +25,13 @@
#include
#include
#include
+#include
#include
-#include "Feature.h"
-#include "utils/print.h"
-
namespace ov_core {
+class Feature;
+
/**
* @brief Database containing features we are currently tracking.
*
@@ -65,17 +65,15 @@ class FeatureDatabase {
* @param remove Set to true if you want to remove the feature from the database (you will need to handle the freeing of memory)
* @return Either a feature object, or null if it is not in the database.
*/
- std::shared_ptr get_feature(size_t id, bool remove = false) {
- std::lock_guard lck(mtx);
- if (features_idlookup.find(id) != features_idlookup.end()) {
- std::shared_ptr temp = features_idlookup.at(id);
- if (remove)
- features_idlookup.erase(id);
- return temp;
- } else {
- return nullptr;
- }
- }
+ std::shared_ptr get_feature(size_t id, bool remove = false);
+
+ /**
+ * @brief Get a specified feature clone (pointer is thread safe)
+ * @param id What feature we want to get
+ * @param feat Feature with data in it
+ * @return True if the feature was found
+ */
+ bool get_feature_clone(size_t id, Feature &feat);
/**
* @brief Update a feature object
@@ -90,33 +88,7 @@ class FeatureDatabase {
* This will update a given feature based on the passed ID it has.
* It will create a new feature, if it is an ID that we have not seen before.
*/
- void update_feature(size_t id, double timestamp, size_t cam_id, float u, float v, float u_n, float v_n) {
-
- // Find this feature using the ID lookup
- std::lock_guard lck(mtx);
- if (features_idlookup.find(id) != features_idlookup.end()) {
- // Get our feature
- std::shared_ptr feat = features_idlookup.at(id);
- // Append this new information to it!
- feat->uvs[cam_id].push_back(Eigen::Vector2f(u, v));
- feat->uvs_norm[cam_id].push_back(Eigen::Vector2f(u_n, v_n));
- feat->timestamps[cam_id].push_back(timestamp);
- return;
- }
-
- // Debug info
- // PRINT_DEBUG("featdb - adding new feature %d",(int)id);
-
- // Else we have not found the feature, so lets make it be a new one!
- std::shared_ptr feat = std::make_shared();
- feat->featid = id;
- feat->uvs[cam_id].push_back(Eigen::Vector2f(u, v));
- feat->uvs_norm[cam_id].push_back(Eigen::Vector2f(u_n, v_n));
- feat->timestamps[cam_id].push_back(timestamp);
-
- // Append this new feature into our database
- features_idlookup[id] = feat;
- }
+ void update_feature(size_t id, double timestamp, size_t cam_id, float u, float v, float u_n, float v_n);
/**
* @brief Get features that do not have newer measurement then the specified time.
@@ -125,46 +97,7 @@ class FeatureDatabase {
* For example this could be used to get features that have not been successfully tracked into the newest frame.
* All features returned will not have any measurements occurring at a time greater then the specified.
*/
- std::vector> features_not_containing_newer(double timestamp, bool remove = false, bool skip_deleted = false) {
-
- // Our vector of features that do not have measurements after the specified time
- std::vector> feats_old;
-
- // Now lets loop through all features, and just make sure they are not old
- std::lock_guard lck(mtx);
- for (auto it = features_idlookup.begin(); it != features_idlookup.end();) {
- // Skip if already deleted
- if (skip_deleted && (*it).second->to_delete) {
- it++;
- continue;
- }
- // Loop through each camera
- // If we have a measurement greater-than or equal to the specified, this measurement is find
- bool has_newer_measurement = false;
- for (auto const &pair : (*it).second->timestamps) {
- has_newer_measurement = (!pair.second.empty() && pair.second.at(pair.second.size() - 1) >= timestamp);
- if (has_newer_measurement) {
- break;
- }
- }
- // If it is not being actively tracked, then it is old
- if (!has_newer_measurement) {
- feats_old.push_back((*it).second);
- if (remove)
- features_idlookup.erase(it++);
- else
- it++;
- } else {
- it++;
- }
- }
-
- // Debugging
- // PRINT_DEBUG("feature db size = %u\n", features_idlookup.size())
-
- // Return the old features
- return feats_old;
- }
+ std::vector> features_not_containing_newer(double timestamp, bool remove = false, bool skip_deleted = false);
/**
* @brief Get features that has measurements older then the specified time.
@@ -172,46 +105,7 @@ class FeatureDatabase {
* This will collect all features that have measurements occurring before the specified timestamp.
* For example, we would want to remove all features older then the last clone/state in our sliding window.
*/
- std::vector> features_containing_older(double timestamp, bool remove = false, bool skip_deleted = false) {
-
- // Our vector of old features
- std::vector